kill(all);

load("newton1") ;

assume(G>0,Kt>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
. "d" is the characteristic length of the cross section with respect 
  to the restrained warping Vlasov torsion theory

2*l = eta*d

"eta" is the normalized profile span length, 
measured in number of "d" characteristic lengths;

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

l:eta*d/2;

here:
. "r   " is the average twist rate along the profile
. "r· l" is the relative rotation between the mispan section, and the two terminals
. "r·2l" is the relative rotation between the two terminals

tentative solution, in terms of the twist angle psi(z), and of the C1,C3 constants, to be determined
based on the b.c.; the -l <= z <= l axial coordinate is zero at the profile mispan.

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

the symmetric terms have been removed from the tentative solution (the problem is skew-symmetric).

tentative solution check, with respect to the ODE

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

Retrieval of the A,B constants

imposed counterotation at the two terminals

eq1 : psi(l)-psi(-l)=2*r*l;

imposed zero warping at terminals

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

solution

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

torsional moment, De Saint Venant and Vlasov contributions

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

overall applied torsional moment, measured as T_Vla at the ends (TdSV=0 there)

T:ev(T_Vla,z=l);

graphical representation of the results, normalized with respect to the
torsional moment predicted by the unrestrained torsion theory

T_ref:(G*Kt*r);

plotme:fullratsimp(
[
    ev(T     / T_ref      ,infeval), /* overall applied torsional moment*/
    ev(T_Vla / T_ref,z=  0,infeval), /* T_Vla contribution, measured at midspan */
    ev(T_dSV / T_ref,z=  0,infeval), /* T_dSV contribution, measured at midspan */
    1
]
);
wxplot2d(
    plotme,
    [eta,0.1,100],[y,0.1,2000],
    logx,logy,
    [legend,"T_{tot} = T_{Vla}, ext.","T_{Vla}, midspan","T_{dSV}, midspan","1"], 
    [xlabel, "eta = 2l / d"],
    [ylabel, "T/T_{ref}"],
    [title,"Torsional stiffness"]  
);

Bimoment generalized stress component for the Vlasov beam

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

axial stresses at both ends vary proportionally with such a quantity

ev(B/(r*2*l)/G/Kt,sols,z=l);
wxplot2d(
    %,
    [eta,1,100],
    [logx],[logy],
    [xlabel, "η = 2l / d"],
    [ylabel,"B/(G * Kt * [2 theta])"],
    [title,"Scaling factor for the Vlasov induced stress comps at the end"]
);

Definition of the "S" stiffening factor due to the restrained warping, Vlasof effect

S:  T / T_ref;

reformulation based on a single hyperbolic function

S : ratsubst(tanh(eta/2)*cosh(eta/2),sinh(eta/2),S);

stiffening effect due to Vlasov: results extracted from FE simulations
S: [restrained warping stiffness] / [free warping stiffness]
        ctw         otw
l=40x01 6.209E-02   5.618E02
l=40x08 4.192E-03   1.154E02
l=40x16 2.099E-03   3.389E01
l=40x32 1.048E-03   8.988E00
l=40x64 5.242E-04   2.413E00

ctw01:[myl=40*01, myS=1+6.209e-02]$
ctw08:[myl=40*08, myS=1+4.192E-03]$
ctw16:[myl=40*16, myS=1+2.099E-03]$
ctw32:[myl=40*32, myS=1+1.048e-03]$
ctw64:[myl=40*64, myS=1+5.242e-04]$

otw01:[myl=40*01, myS=1+5.618e02]$
otw08:[myl=40*08, myS=1+1.154e02]$
otw16:[myl=40*16, myS=1+3.389e01]$
otw32:[myl=40*32, myS=1+8.988e00]$
otw64:[myl=40*64, myS=1+2.413e00]$

characteristic lengths vs. cross section size, ctw case

csize:sqrt(80^2+40^2),numer$ /*diagonal of the hex20profile cross section*/

2*myl/newton(myS/S-1,eta,1,1e-6),ctw08$[%,%/csize];
2*myl/newton(myS/S-1,eta,1,1e-6),ctw16$[%,%/csize];
2*myl/newton(myS/S-1,eta,1,1e-6),ctw32$[%,%/csize];
2*myl/newton(myS/S-1,eta,1,1e-6),ctw64$[%,%/csize];

2*myl/newton(myS/S-1,eta,1,1e-6),otw08$[%,%/csize];
2*myl/newton(myS/S-1,eta,1,1e-6),otw16$[%,%/csize];
2*myl/newton(myS/S-1,eta,1,1e-6),otw32$[%,%/csize];
2*myl/newton(myS/S-1,eta,1,1e-6),otw64$[%,%/csize];

H profile control FE model
"H_profile_vlasov.mud"

analytical predictions

_t:1$_b:30$_h:40$_E:1000$_G:_E/2.6$
_ECw:_E*_h^2*(_b^3*_t/12)/2$
_GKt:_G*(2*_b+_h)*_t^3/3$
_csize:sqrt(_b^2+_h^2),numer;

_d:sqrt(_ECw/_GKt)$[%,%/_csize];

FE model stiffening factor

Hprof:[myl=80, myS=1.214e-1/1.977e-3]$

_d_FE : 2*myl/newton(myS/S-1,eta,1,1e-6),Hprof$[%,%/_csize];

relative mismatch between _d ad _d_FE

(2*(_d-_d_FE)/(_d+_d_FE));