# Kapitel 3: 
# Eigenfunktionen des beidseitig drehbar
# gelagerten Timoshenko-Balkens (z > zg)
# © 2017  Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
> restart: with(LinearAlgebra): unprotect(Zeta); Digits:=15:
# 
# Für den skizzierten Balken  mit Rechteckquerschnitt der Breite b, der Höhe h und der Länge  sind die Eigenkreisfrequenzen und Eigenfunktionen der zeitfreien Verschiebungen W(x) und Drehwinkel ¨(x) zu berechnen. Die ortsahängige Differenzialgleichung lautet:
# Hinweis: Die in Beispiel 3-2_a hergeleiteten Beziehungen gelten auch für zn > zg: Im Fall z > zg tritt mit sin (»2) = 0 eine zweite Reihe von Eigenwerten auf, die es hier zu berechnen gilt.
> odeW:=diff(W(x),x$4)+kappa[B]^4*iy^2*(1+alpha)*diff(W(x),x$2)-kappa[B]^4*(1-kappa[B]^4*iy^4*alpha)*W(x);
# Diese gewöhnliche Differenzialgleichung ist ein Folge des Separationsansatzes für harmonische Schwingungen der Form  w(x,t) = W(x)Åsin(Ét + Æ) . Sie hat die Lösung:
> dsolve(odeW,W(x)): assign(%); simplify(W(x));
# Zur Verkürzung der Schreibweise führen wir die Variable z ein:
> kappa[B]:=z/iy;
> simplify(map(eval,W(x)));
# Wir wandeln die Exponentialdarstellung noch etwas um und erhalten -jetzt mit neuen Konstanten- für die ortsabhängige Verschiebung:
> W:=unapply(_C1*sin(lambda[1]*x)+_C2*cos(lambda[1]*x)+_C3*sin(lambda[2]*x)+_C4*cos(lambda[2]*x),x);
# Die Argumente »1 und »2 werden weiter unten eingeführt. Die harmonische Schwingung selbst ist dann:
> w:= W(x)*sin(omega*t + phi);
# Mit w(x,t) = W(x) sin(Ét + Æ) ergibt sich die Ortsfunktion des Drehwinkels ¨(x) aus der Differenzialgleichung È´ = ñA/(GAS )  w¨  -  w´´ zu
> PSI:=Zeta[1]*(_C1*cos(lambda[1]*x)-_C2*sin(lambda[1]*x))+Zeta[2]*(_C3*cos(lambda[2]*x)-_C4*sin(lambda[2]*x));
> PSIstr:=simplify(diff(PSI,x));
# Lösung des Randwertproblems: An den Balkenenden müssen jeweils  die Verschiebungen und die Biegemomente verschwinden. Das führt auf  das folgende lineare homogene Gleichungssystem:
> gl1:=eval(subs(x=0,W(x)));
> gl2:=eval(subs(x=0,PSIstr));
> gl3:=eval(subs(x=l,W(x)));
> gl4:=eval(subs(x=l,PSIstr));
> DA, b := GenerateMatrix([gl1,gl2,gl3,gl4],[_C1,_C2,_C3,_C4]);
> EWG:=simplify(Determinant(DA));
> EWG:=sin(lambda[1]*l)*sin(lambda[2]*l);
> lambda[2]:=sqrt(2)/2/iy*z*sqrt((alpha+1)*z^2-sqrt(4+z^4*(alpha-1)^2));
> Zeta[2]:=simplify((alpha*z^4/iy^2 - lambda[2]^2)/lambda[2]);
> EWG2:=sin(lambda[2]*l); _C1:=0; _C2:=0; _C3:=1; _C4:=0;
> gl2:=nu-iy*lambda[2]; loe2:=solve(gl2,z);
# Im obigen Ausdruck wurde folgende Abkürzung eingeführt:
> nu:=n*Pi*iy/l;
# Mit den Werten unseres Beispiels für einen kurzen Stahlbetonbalken ( / h = 5)
> b:=0.2; h:=0.6; A:=b*h; Iyy:=b*h^3/12; iy:=sqrt(Iyy/A); l:=3.; 
> E:=3.05E10; G:=E/2; rho:=2500.; AS:=5/6*A; alpha:=E*A/(G*AS); zg:=1/alpha^(1/4);
# Wir stellen die Eigenwertgleichung als Funktion von z grafisch dar:
> EWG2:=unapply(EWG2,z); 
> plot(EWG2(z),z=zg..1.2,axes=boxed,gridlines);
> Z0:=[]: a:=zg: b:=1.2: dz:=0.01: dmax:=round((b-a)/dz): za:= a: ze:=za+dz:
> for n to dmax do
>   if evalf(EWG2(za)*EWG2(ze)) < 0. then
>    z0:=fsolve(EWG2(z),z=za..ze):
>     Z0:= [op(Z0),z0]:
>   end if:
> za:= ze; ze:=za+dz;
> end do: 
> nmax:=nops(Z0);
# Wir ermitteln die ersten nmax  Eigenkreisfrequenzen É und die Eigenfrequenzen f :
>  ZN:=Vector(nmax): OMEGAN:=Vector(nmax): FRN:=Vector(nmax): 
>  g:=sqrt(E*A/rho/Iyy):znaeh:=g*z1^2;
> for n to nmax do
>   ZN[n]:=Z0[n];
>   OMEGAN[n]:=evalf(g*ZN[n]^2) ; 
>   FRN[n]   :=OMEGAN[n]/2/Pi;
> end do:
> print(`z = `,ZN,`É = `,OMEGAN,`f = `, FRN);
# Wir stellen die Zustandsgrößen zum 1. Eigenwert grafisch dar:
> z:=ZN[1];
> QN:=diff(W(x),x)+PSI: MN:=l*PSIstr: Wstr:=diff(W(x),x):
> para:=x=0..l,gridlines=true,titlefont = ["ARIAL", 15],  labeldirections = ["horizontal", "horizontal"], labelfont = ["HELVETICA", 10],axesfont = ["HELVETICA", "ROMAN", 8],axes=boxed:
> p1:= plot(QN  ,para,title = "\nNormierte Querkraft\n"   ,labels = ["x", "QN(x)"]):
> p2:= plot(MN  ,para,title = "\nNormiertes Biegemoment\n" ,labels = ["x", "MN(x)"]):
> p3:= plot(W(x),para,title = "\nVerschiebung\n",labels = ["x", "W(x)"]):
> p4:= plot(PSI ,para,title = "\nDrehwinkel\n"  ,labels = ["x", "È(x)"]):
> plots[display](Matrix(1,4,[p1,p2,p3,p4]));
# 
