kill(   
    all   
);

assume(a>0,r>0,t>0,t<1);

isotropic material with customary poisson ratio

nu : 3/10      $ 
G  : E/2/(1+nu) $

circular, closed thin walled cross section

Kt:2*J*eta;

I force the usage of the two adimensional factors from the beginning

b:r*a;
h:t*a;

l : a-h;

I declare two auxiliary functions for linear interpolation along a segment

note: the evalution of the RHS of := is suspended (not performed) at define/assigment time; the suspended evaluation is performed in addition to the usual ones at call time

N1(xi) := 1-xi;

the define function, instead, evaluates as usual its second argument, then it assigns it to the first (which is NEVER evaluated, since object of further assigment, this is an exception of the general rule)

define(
    N2(goofy),
    goofy
);

define(
    mylinearinterpolationonBC(s),
    factor( vBp*N1(s/l)+vCm*N2(s/l) )
);

I initialize the strain energy as a null term

U:0$

Components of internal action along the AB segment, out-of-plane components alone

M2Ap :  fullratsimp( 0*H + 0*ZA + 1*ThetaA + 0*PhiA + 0*I );
MtAp :  fullratsimp( 0*H + 0*ZA + 0*ThetaA + 1*PhiA + 0*I );

M2Bm :  fullratsimp( 0*H + h*ZA + 1*ThetaA + 0*PhiA + 0*I );
MtBm :  fullratsimp( 0*H + 0*ZA + 0*ThetaA + 1*PhiA + 0*I );

define(
    M2AB(s),
    expand(M2Ap*N1(s/h)+M2Bm*N2(s/h))
);
define(
    MtAB(s),
    expand(MtAp*N1(s/h)+MtBm*N2(s/h))
);

U : U + integrate( 
    M2AB(s)^2/2/E/J + MtAB(s)^2/2/G/Kt, s , 0 , h 
);

BC segment

M2Bp : fullratsimp( 0*H + h*ZA + 1*ThetaA + 0*PhiA + 0*I);
MtBp : fullratsimp( 0*H + 0*ZA + 0*ThetaA + 1*PhiA + 0*I);

M2Cm : fullratsimp( 0*H + a*ZA + 1*ThetaA + 0*PhiA + 0*I);
MtCm : fullratsimp( 0*H + 0*ZA + 0*ThetaA + 1*PhiA + 0*I);

define(
    M2BC(s),
    expand(M2Bp*N1(s/l)+M2Cm*N2(s/l))
);
define(
    MtBC(s),
    expand(MtBp*N1(s/l)+MtCm*N2(s/l))
);

U : U + integrate( 
    M2BC(s)^2/2/E/J + MtBC(s)^2/2/G/Kt, s , 0 , l 
);

fullratsimp(U);

CD segment

M2Cp : fullratsimp( 0*H + 0*ZA + 0*ThetaA - 1*PhiA + 0*I);
MtCp : fullratsimp( 0*H + a*ZA + 1*ThetaA + 0*PhiA + 0*I);

M2Dm : fullratsimp(-b*H + b*ZA + 0*ThetaA - 1*PhiA + 0*I);
MtDm : fullratsimp( 0*H + a*ZA + 1*ThetaA + 0*PhiA + 0*I);

define(
    M2CD(s),
    expand(M2Cp*N1(s/b)+M2Dm*N2(s/b))
);
define(
    MtCD(s),
    expand(MtCp*N1(s/b)+MtDm*N2(s/b))
);

U : U + integrate( 
    M2CD(s)^2/2/E/J + MtCD(s)^2/2/G/Kt, s , 0 , b 
);

DE segment

M2Dp : fullratsimp( 0*H - a*ZA - 1*ThetaA + 0*PhiA + 0*I);
MtDp : fullratsimp(-b*H + b*ZA + 0*ThetaA - 1*PhiA + 0*I);

M2Em : fullratsimp(-a*H + 0*ZA - 1*ThetaA + 0*PhiA - a*I);
MtEm : fullratsimp(-b*H + b*ZA + 0*ThetaA - 1*PhiA + 0*I);

define(
    M2DE(s),
    expand(M2Dp*N1(s/a)+M2Em*N2(s/a))
);
define(
    MtDE(s),
    expand(MtDp*N1(s/a)+MtEm*N2(s/a))
);

U : U + integrate( 
    M2DE(s)^2/2/E/J + MtDE(s)^2/2/G/Kt, s , 0 , a 
);

fullratsimp(U);

I proceed evaluating the partial derivatives of the strain energy with respect to the applied gen. loads in order to obtain the associated gen. displacements

diff(U,H)$
ratsimp(%)$
wtC:%$

diff(U,I)$
ratsimp(%)$
wtD:%$

diff(U,ZA)$
ratsimp(%)$
wA:%$

diff(U,ThetaA)$
ratsimp(%)$
thetaA:%$

diff(U,PhiA)$
ratsimp(%)$
phiA:%$

After its usage within the Castigliano theorem application, I may be finally declared as "fictitious"

I:0;

I now enforce the constraints at A by imposing null deflections there.

eqns : [ wA=0 , thetaA=0, phiA=0 ]$ 

unks : [ ZA, ThetaA, PhiA ];

implicit_ev( linsolve ( eqns , unks ) , globalsolve=true )

linsolve ( eqns , unks ) , globalsolve=true$

variant without "globalsolve"
sols:linsolve ( eqns , unks );
[XA,YA,PsiA]:ev([XA,YA,PsiA],sols);

it's time to substitute the obtained values for the constraint reactions at A within the sampled displacements at B and C.

ev([wtC,wtD])$
[wtC,wtD]:factor(%);

## DO NOT RUN RIGHT NOW

Evaluation of a reference deflection obtained for a cantilever beam structure equivalent to a single upright frame member. 

bending moment for the ref. structure

Mfref:H*s$

strain energy for the ref. structure

Uref:integrate(Mfref^2/2/E/J,s,0,a);

I evaluate the deflection by resorting to the II Castigliano theorem

dref : diff(Uref,H);

I compare the frame stiffness (k=P/utC) with the stiffness of the one of its upright members alone, thus obtaining a "stiffness ratio"

oopkratio : (H/wtC) / (H/dref);

I plot this normalized stiffness ratio

wxplot2d( 
    [
        ev(oopkratio,eta=1),
        ev(oopkratio,eta=0.1),
        ev(oopkratio,eta=0.001)
    ], 
    [ r , 0.01 , 100 ] , 
    [logx] , 
    [legend, 
        "oopkratio, Kt/(2J)=1",
        "oopkratio, Kt/(2J)=0.1",
        "oopkratio, Kt/(2J)=0.01" 
    ],
    [xlabel, "r=b/a ratio"],
    [ylabel, "oopkratio"]
);

oopkratio, r=1,eta=1,numer;