kill(all);

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

qui :
-- l minuscolo è la SEMIlunghezza della trave
-- L maiuscolo è la     lunghezza della trave, L = 2l

qui :
-- theta minuscolo è la rotazione di un terminale rispetto alla mezzeria
-- Theta maiuscolo è la rotazione relativa tra i due terminali, Theta = 2 theta

termine di lunghezza inversa caratteristica per il warping:  se L*beta >> 1, il warping ristretto non si sente.

beta : sqrt(K*G/Cw/E);

rotazione con warping impedito ad ambo le estremità, prima forma (da testo di riferimento)

Theta = T/Cw/E/beta^3 * (beta*L-2*tanh(beta*L/2));

eq1 : expand(%);

rotazione con warping impedito ad ambo le estremità, seconda forma (utilizzata su wiki)

eq2:Theta = expand( L*T/G/K*(1- 2/beta/L * tanh(beta*L/2)));

equivalenza tra le due forme

eq1-eq2;

f(x) := 1/x*tanh(x);

g(y) := y * tanh(1/y);

diff(f(x),x);

exponentialize(%);

diff(g(y),y);

exponentialize(%);

wxplot2d([f(x)],[x,0,10]);

questa sembra più gestibile da un N-R, se g(t)<1...

qui t=

wxplot2d([g(t),t],[t,0,10],[y,0,1]);

soluzione di tentativo

phi(x):=A*sinh(BETA*x)+B*x;

controllo equazione differenziale

fullratsimp(
    diff(phi(x),x,4)-BETA^2*diff(phi(x),x,2)=0
);

rotazione terminali

eq1 : phi(L/2)=Theta/2;

momento torcente, quota De Saint Venant

T_dSV :              G*K*diff(phi(x),x,1);
T_Vla : - (1/BETA^2)*G*K*diff(phi(x),x,3);

warping impedito ai terminali

0=diff(phi(x),x,1);
eq2 : ev(%,x=L/2);

sols:linsolve([eq1,eq2],[A,B]);

fullratsimp(ev(T_dSV,sols));

plotme:fullratsimp(
[
    ev(T_Vla / (G*K*Theta/L),sols,x=L/2,L=eta/BETA,infeval),
    ev(T_Vla / (G*K*Theta/L),sols,x=  0,L=eta/BETA,infeval),
    ev(T_dSV / (G*K*Theta/L),sols,x=  0,L=eta/BETA,infeval),
    1
]
);
wxplot2d(
    plotme,
    [eta,1,100],
    logx,
    [legend,"T_{tot} = T_{Vla}, est.","T_{Vla}, mezz.","T_{dSV}, mezz.","1"], 
    [xlabel, "2l / d_{Vla}"],
    [ylabel, "(T * 2l)/(G * Kt * [2 theta])"],
    [title,"Rigidezza torsionale trave a torsione"]  
);

diff(-G*K/BETA^2*phi(x)/Theta/G/K,x,2)$
ev(%,sols)$
ev(%,x=L/2)$
ev(%,L=eta/BETA);
wxplot2d(
    %,
    [eta,1,100],
    logx,logy,
    [legend,"B/(G * Kt * [2 theta])"], 
    [xlabel, "2l / d_{Vla}"],
    [ylabel,"B/(G * Kt * [2 theta])"],
    [title,"Fattore di scalatura delle tensioni assiali Vlasov di estremità"] 
);

procedura Newton per la soluzione dell'eq. nonlineare

1/Amp_Vla-1/ev(T_Vla / (G*K*Theta/L),sols,x=L/2,L=eta/BETA,infeval);
ratsubst(
    tanh(eta/2),sinh(eta/2)/cosh(eta/2),% 
);
define(f(Amp_Vla,eta) ,
    expand(%)
);
define(
        df(Amp_Vla,eta),
        expand(diff(f(Amp_Vla,eta),eta))
    );
define(
    NR(Amp_Vla,eta),
    fullratsimp(eta - f(Amp_Vla,eta)/df(Amp_Vla,eta))
);

amp:1.1;
1;
NR(amp,%),numer;
NR(amp,%),numer;
NR(amp,%),numer;
NR(amp,%),numer;
NR(amp,%),numer;
NR(amp,%),numer;
NR(amp,%),numer;
NR(amp,%),numer;
NR(amp,%),numer;
NR(amp,%),numer;
NR(amp,%),numer;
NR(amp,%),numer;
NR(amp,%),numer;
NR(amp,%),numer;
NR(amp,%),numer;
NR(amp,%),numer;
NR(amp,%),numer;
NR(amp,%),numer;