kill(all);

assume(G>0,K>0,Cw>0,E>0);

here :
-- lowercase l is the overall, twice clamped beam HALF-LENGTH, 
   i.e. the length of an equivalent clamped-free beam
   i.e. the length of the FE modeled portion, if skew-symm. is imposed at the midspan
-- uppercase L is the overall, twice clamped beam length, 

L = 2l = eta*d

d is the Vlasov cross sectional characteristic length, i.e.

d^2=(E*C_w)/(G*K)

L:eta*d;

l:L/2;

"r" is the average twist rate along the profile, i.e.
-- r*l is the relative rotation between the mispan section, and the two terminals
-- r*L is the relative rotation between the two terminals

the nominal, expected according to the de St. Venant torsion torque is

T0:G*K*r;

cross-section characteristic length with respect to Vlasov theory of restrained warping

Restrained warping effect is negligible if L/d >> 1 .

tentative solution, skew-symm. terms only

psi(z):=C1*sinh(z/d)+C3*z/d;

tentative solution check

fullratsimp(
    -d^2*diff(psi(z),z,4)+diff(psi(z),z,2)=0
);

imposed rotation at the z=l terminal

eq1 : psi(l)=r*l;

due to skew-symmetry, the z=-l terminal counterrotates of an equal angle

imposed zero warping at the z=l terminal (and also at z=-l due to skew-symm.)

diff(psi(z),z,1);
eq2 : ev(%=0,z=l);

solution

linsolve([eq1,eq2],[C1,C3]),globalsolve=true;

psi(z);

torsional moments, De Saint Venant and Vlasov portions

T_dSV :         G*K*diff(psi(z),z,1);
T_Vla : - (d^2)*G*K*diff(psi(z),z,3);

plotme:fullratsimp(
[
    ev(T_Vla / T0,z=  l,infeval),
    ev(T_Vla / T0,z=  0,infeval),
    ev(T_dSV / T0,z=  0,infeval),
    1
]
);
wxplot2d(
    plotme,
    [eta,0.5,100],
    logx,
    [legend,"T_{tot} = T_{Vla}, ext.","T_{Vla}, midspan","T_{dSV}, midspan","1"], 
    [xlabel, "2l / d"],
    [ylabel, "T/T0"],
    [title,"Torsional stiffening and torque components"]  
);

Bimoment B evaluation along the beam

B:fullratsimp( -G*K*d^2*diff(psi(z),z,2) );

normalized bimoments at beam extremities, 
with respect to the T0*d product

plotme:fullratsimp( ev(B,z=l)/(T0*d) );

wxplot2d(
    plotme,
    [eta,1,100],
    logx,logy,
    [legend,"B/(T0 * d)"], 
    [xlabel, "2l / d"],
    [ylabel,"B/(T0 * d)"],
    [title,"Scaling factor for the Vlasov induced stress comps at ends"] 
);

Definition of a stiffening factor due to the Vlasov effect

overall torque, measured e.g. at terminals 

T : ev(T_Vla+T_dSV,z=l);

stiffening factor due to restrained warping at ends

T/T0;
ratsubst(tanh(eta/2)*cosh(eta/2),sinh(eta/2),%)$
RWSF :%;

Newton-Raphson procedure for estracting the 

d/L = eta 

ratio from the Vlasov stiffening factor RWSF 

residual definition; 

1/x = 1/RWSF

such 1/RWSF function is well defined on the 0.001 <= eta <=1000 range

wxplot2d(1/RWSF,[eta,0.001,1000],[logx]);

wxplot2d(1/RWSF,[eta,0.001,1000],[logx],[logy]);

definition of a Newton procedure:
residual function and its derivatives

TARGET is the target value for RWSF(eta), obtained from FE analyses

res : 1/TARGET-1/RWSF;

define(f(TARGET,eta) ,
    expand(res)
);

define(
        df(TARGET,eta),
        expand(diff(f(TARGET,eta),eta))
    );

define(
    NR(TARGET,eta),
    fullratsimp(eta - f(TARGET,eta)/df(TARGET,eta))
);

stiffening effect due to Vlasov

r=0.001 rad/mm
l=1200 mm

RWSF_FEM: [restrained warping stiffness] / [free warping stiffness]

RWSF_FEM: 10;

l_FEM:1200;

a few Newton iterations for the eta_FEM value

1.0;
NR(RWSF_FEM,%),numer;
NR(RWSF_FEM,%),numer;
NR(RWSF_FEM,%),numer;
NR(RWSF_FEM,%),numer;
NR(RWSF_FEM,%),numer;
NR(RWSF_FEM,%),numer;
NR(RWSF_FEM,%),numer;
NR(RWSF_FEM,%),numer;
NR(RWSF_FEM,%),numer;
NR(RWSF_FEM,%),numer;
NR(RWSF_FEM,%),numer;
NR(RWSF_FEM,%),numer;
NR(RWSF_FEM,%),numer;
NR(RWSF_FEM,%),numer;
NR(RWSF_FEM,%),numer;
NR(RWSF_FEM,%),numer;
NR(RWSF_FEM,%),numer;
NR(RWSF_FEM,%),numer;

once (and if) convergence is obtained, let's proceed with the characteristic distance "d" evaluation

eta_FEM:%;

linsolve(
    ev(l,eta=eta_FEM)=l_FEM,
    d
),numer;