# Kapitel 3: 
# Eigenfunktionen des links eingespannten und 
# rechts drehbar gelagerten Timoshenko-Balkens (z < zg)
# © 2017  Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
> restart: with(LinearAlgebra): with(plots): unprotect(Zeta); 
# Beispiel 3-1
# 
# 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:
> 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, wobei iy den Trägheitsradius des Rechteckquerschnitts bezeichnet:
> kappa[B]:=z/iy;
> simplify(map(eval,W(x)));
# Wir wandeln die Exponentialdarstellung noch etwas um und erhalten mit neuen Konstanten für die ortsabhängige Verschiebung:
> W:=unapply(_C1*sin(lambda[1]*x)+_C2*cos(lambda[1]*x)+_C3*sinh(lambda[2]*x)+_C4*cosh(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*cosh(lambda[2]*x)+_C4*sinh(lambda[2]*x));
> PSIstr:=simplify(diff(PSI,x));
# Lösung des Randwertproblems: Am linken Balkenende müssen die Durchsenkung sowie der Drehwinkel und am rechten Balkenende die Durchsenkung und das Biegemoment verschwinden. Das führt auf das lineare homogene Gleichungssystem:
> gl1:=eval(subs(x=0,W(x)));
> gl2:=eval(subs(x=0,PSI));
> gl3:=eval(subs(x=l,W(x)));
> gl4:=eval(subs(x=l,PSIstr));
> DA, b := GenerateMatrix([gl1,gl2,gl3,gl4],[_C1,_C2,_C3,_C4]);
> GDA:=GaussianElimination(DA);
> EWG:=simplify(Determinant(DA));
# oder kürzer und für die nummerische Abarbeitung besser geeignet
> EWG:=sin(lambda[1]*l)+cos(lambda[1]*l)*tanh(lambda[2]*l)*Zeta[1]/Zeta[2];
# Berechnung der Integrationskonstanten. Da eine Konstante frei gewählt werden kann, setzen wir 
> _C4:=1;
# Das verbleibende Gleichungssystem ist
> AR, b := GenerateMatrix([gl1,gl2,gl3],[_C1,_C2,_C3]);
> unb:=LinearSolve(AR,b);
# Damit sind 
> _C1:=unb[1]; _C2:=unb[2]; _C3:=unb[3];
# Nach dem Ersatz der Exponentialfunktionen durch trigonometrische und hyperbolische Funktionen, sind  »1 und »2 wie folgt zu setzen:
> lambda[1]:=sqrt(2)/2/iy*z*sqrt(sqrt(4+z^4*(alpha-1)^2)+(alpha+1)*z^2);
> lambda[2]:=sqrt(2)/2/iy*z*sqrt(sqrt(4+z^4*(alpha-1)^2)-(alpha+1)*z^2);
# und  ¶1 und ¶2 ergeben sich zu
> Zeta[1]:=simplify((alpha*z^4/iy^2 - lambda[1]^2)/lambda[1]);
> Zeta[2]:=simplify((alpha*z^4/iy^2 + lambda[2]^2)/lambda[2]);
# Beispiel 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);
> lambda[1];lambda[2];
> Zeta[1]; Zeta[2];
> zg:=1/alpha^(1/4); 
> EWG:=unapply(EWG,z); 
> plot(EWG(z),z=0..zg,title = "\nEigenwertgleichung EWG(z)\n",titlefont = ["ARIAL", 15],axes=boxed,gridlines=true);
# In der folgendenSchleife werden dieNullstellen der  Eigenfunktion EWG(z) berechnet. Dazu wird EWG(z)  im Intervall [a,b] mit der Schrittweite dz ausgewertet. Ist EWG(z)*EWG(z+dz)  < 0, dann hat die Eigenfunktion  in diesem Intervall einen Nulldurchgang. In diesem Fall wird die Nullstelle mit  fsolve  nummerisch konkretisiert. 
> Z0:=[]: a:=0.1: b:=zg: dz:=0.01: dmax:=round((b-a)/dz): za:= a: ze:=za+dz:
> for n to dmax do
>   if evalf(EWG(za)*EWG(ze)) < 0. then
>     Z0:=[op(Z0),fsolve(EWG(z),z=za..ze)]:
>   end if:
>   za:= ze; ze:=za+dz;
> end do: 
> nmax:=nops(Z0): f:=sqrt(E*A/rho/Iyy):
> ZN:=convert(Z0,Vector): OMEGA:=Vector(nmax): FR:=Vector(nmax):
> for n to nmax do
>   OMEGA[n]:=evalf(Z0[n]^2*f);
>   FR[n]   :=OMEGA[n]/2/Pi;
> end do:
> print(`z = `,ZN,`É = `,OMEGA, `f = `, FR);
# 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", "vertical"], 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 W(x)\n",labels = ["x", "W(x)"]):
> p4:= plot(PSI ,para,title = "\nDrehwinkel ¨(x)\n"  ,labels = ["x", "È(x)"]):
> plots[display](Matrix(1,4,[p1,p2,p3,p4]));
# In einer weiteren Grafik geben wir abschließend die ersten drei Durchsenkungs-Eigenfunktionen grafisch aus:
> WX:=[]:PSIX:=[]:
> for i to 3 do 
>   z   :=ZN[i];
>   WX  :=[op(WX),W(x)];
>   PSIX:=[op(PSIX),PSI];
> end do:
> plot(WX,x=0..l,gridlines=true,axes=boxed,title = "\nDurchsenkungs-Eigenfunktionen\n", legend = ["Grundschwingung", "1. Oberschwingung","2. Oberschwingung"], titlefont = ["ARIAL", 15],  linestyle = [solid, dot, dash], axesfont = ["HELVETICA", "ROMAN", 8], legendstyle = [font = ["HELVETICA", 9], location =bottom],color=[black,black,black]);
> plot(PSIX,x=0..l,gridlines=true,axes=boxed,title = "\nDrehwinkel-Eigenfunktionen\n", legend = ["Grundschwingung", "1. Oberschwingung","2. Oberschwingung"], titlefont = ["ARIAL", 15],  linestyle = [solid, dot, dash], axesfont = ["HELVETICA", "ROMAN", 8], legendstyle = [font = ["HELVETICA", 9], location =bottom],color=[black,black,black]);
> 
;
