# Kapitel 3: 
# Der Kragträger unter konstanter Druckkraft 
# mit Berücksichtigung der Drehträgheit
# © 2017  Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
> restart: with(LinearAlgebra): with(plots): 
# 
# 
> ode4:=diff(W(x),x$4) + K^2*diff(W(x),x$2)-kB^4*W(x);
> dsolve(ode4,W(x));
# Für den in der obigen Abbildung dargestellten prismatischen Kragbalken unter Normalkraftbeanspruchung (Zug oder Druck) sind die Eigenwerte und Eigenkreisfrequenzen unter Berücksichtigung der Drehträgheit zu ermitteln. Mit  È = P/(EIyy) gilt dann die zeitfreie Transversalverschiebung 
> W(x):=_C1*sin(lambda[1]*x)+_C2*cos(lambda[1]*x)+_C3*sinh(lambda[2]*x)+_C4*cosh(lambda[2]*x);
# Wir benötigen die Ableitungen von W(x) 
> W1:=simplify(diff(W(x),x$1));
> W2:=simplify(diff(W(x),x$2));
> W3:=simplify(diff(W(x),x$3));
# Die zeitfreien Schnittlasten sind
# 1.) Zeitfreies Biegemoment MD
> MD:=simplify(-E*Iyy*W2,symbolic);
# 2.) Zeitfreie Querkraft QD
> QD:=-E*Iyy*(K^2*W1 + W3);
# Wir lösen jetzt das Randwertproblem des Kragträgers. Am linken eingespannten Ende müssen W und W` verschwinden. Am freien rechten Ende müssen M und Q verschwinden. Das erfordert
> gl1:=eval(subs(x=0,W(x)));
> gl2:=eval(subs(x=0,W1));
> gl3:=eval(subs(x=l,W2));
# Für die zeitfreie Querkraft am rechten Ende ist zu fordern:
> gl4:=simplify(subs(x=l,K^2*W1 + W3));
> gl:=[gl1,gl2,gl3,gl4]: unb:=[_C1,_C2,_C3,_C4]:
# Wir beschaffen uns die Systemdeterminante des homogenen Gleichungssystems
> AS, b := GenerateMatrix(gl, unb);
> DA:=simplify(Determinant(AS));
> simplify(GaussianElimination(AS));
# Die Gaußelimination zeigt, dass im Zähler des Elementes AS[4,4] gerade ein Vielfaches der Systemdeterminante von AS steht, die zu Null gefordert wird. Damit ist C4 ` 0. Wir setzen:
> _C4:=C;
# Die Berechnung der verbleibenden Konstanten ergibt
> gl:=[gl1,gl2,gl3]: unb:=[_C1,_C2,_C3]:
> ASR, c := (GenerateMatrix(gl, unb));
> UNB:=LinearSolve(ASR, c);
> _C1:=UNB[1];_C2:=UNB[2];_C3:=UNB[3];
# Etwas umsortiert und für die nummerische Bearbeitung günstiger ist die Eigenwertgleichung DA in der Form
> DA:=cos(lambda[1]*l)-k1*sin(lambda[1]*l)*tanh(lambda[2]*l) + k2/cosh(lambda[2]*l);
# Dabei wurden folgende Abkürzungen eingeführt:
> k1:=simplify(2*lambda[1]*(-(1/2)*lambda[1]^2+(1/2)*lambda[2]^2+K^2)*lambda[2]/((2*lambda[2]^2+K^2)*lambda[1]^2-K^2*lambda[2]^2));
> k2:=simplify((lambda[1]^4-K^2*lambda[1]^2+lambda[2]^2*(lambda[2]^2+K^2))/((2*lambda[2]^2+K^2)*lambda[1]^2-K^2*lambda[2]^2));
> lambda[1]:=(1/2)*sqrt(2)*sqrt(sqrt(K^4 + 4*kB^4) + K^2);
> lambda[2]:=(1/2)*sqrt(2)*sqrt(sqrt(K^4 + 4*kB^4) - K^2);
> 
;
# Die folgende Variable  VP  legt das Vorzeichen der Normalkraft N  fest. 
# VP = -1 :  N  ist eine Druckkraft
# VP =  0 :  N  = 0
# VP =  1:   N  ist eine Zugkraft
> VP:=-1:
> if VP = -1 then    #N ist eine Druckkraft
>   K:= sqrt(iy^2*kB^4 + kappa^2);
> elif VP = 1 then   #N ist eine Zugkraft
>   K:= sqrt(iy^2*kB^4 - kappa^2);
> elif VP = 0 then   #N = 0
>   K:=sqrt(iy^2*kB^4);
> end if:
# Mit den Werten des Beispiels:
> A:=52.8*10^(-4); Iyy:=4696.*10^(-8); l:=10.;
> rho:=7850.; E:=2.1E11;  By:=E*Iyy; iy:=sqrt(Iyy/A);
> Pkr:=Pi^2*By/4/l^2; epsilon:=.5; P:=epsilon*Pkr; kappa:=sqrt(P/By);
> plot(DA,kB=0..1,axes=boxed,gridlines);
# Berechnung der Eigenwerte
> EWG:=unapply(DA,kB):
# In der folgenden Schleife 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.: b:=3: dz:=0.1: 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*Iyy/A/rho);
> 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:
# Wir geben die Eigenwerte º[B] , Eigenkreisfrequenzen  É und Frequenzen f aus:
> print(`º[B] = `,ZN,`  É = `,OMEGA, `  f = `, FR);
# Wir stellen die Zustandsgrößen zum 2. Eigenwert grafisch dar:
> n:=2; kB:=ZN[n]; C:=1.;
> kQD:=E*Iyy/l^2: kMD:=E*Iyy/l:
> para:=x=0..l,gridlines=true,titlefont = ["ARIAL", 15],  labeldirections = ["horizontal", "vertical"], labelfont = ["HELVDETICA", 10],axesfont = ["HELVDETICA", "ROMAN", 8],axes=boxed:
> p1:= plot(QD/kQD,para,title = "\nNormierte Querkraft\n"   ,labels = ["x", "QN(x)"]):
> p2:= plot(MD/kMD,para,title = "\nNormiertes Biegemoment\n",labels = ["x", "MN(x)"]):
> p3:= plot(W(x)  ,para,title = "\nNormierte Verschiebung\n",labels = ["x", "W(x)"]):
> p4:= plot(W1    ,para,title = "\nNormierte Tangentenneigung\n",labels = ["x", "W '(x)"]):
> plots[display](Matrix(1,4,[p3,p4,p2,p1]));
> kB:='kB':
# Wie Vergleichsrechnungen zeigen, nehmen die Eigenfrequenzen É mit wachsender Druckkraft P ab. Für P = Pkr wird É1 = 0. Wir untersuchen nun den Fall P  = 1. 2 * Pkr:
> epsilon:=1.2; Pkr:=Pi^2*By/4/l^2; P:=epsilon*Pkr; kappa:=sqrt(P/By);  K:= sqrt(iy^2*kB^4 + kappa^2): DA:=map(eval,DA): EWG:=unapply(DA,kB):
# Wir beschaffen uns einen Überblick über die Lage der ersten Nullstellen der Eigenwertgleichung in der komplexen Ebene
> r:=0.15;
> complexplot3d(EWG(kB),kB = -r - r*I.. r + r*I, view=0..0.5);
# Der obigen Grafik entnehmen wir vier komplexe Nullstellen, die symmetrisch zur reelen und imaginären Achse angeordnet sind. Wir berechnen die Nullstelle im 1. Quadranten:
> kB1:=fsolve(EWG(kB),kB = 0.05*(1+I), complex);
# Wir werten die Teilgrößen für diesen Eigenwert aus:
> kB:=kB1:
> k1;
> k2;
> lambda[1];
> lambda[2];
# Die Integrationskonstanten sind:
> C:=1: evalf(_C1); evalf(_C2); evalf(_C3); evalf(_C4);
# Wir geben die normierte zeitfreie Verschiebung zu diesem Eigenwert grafisch aus:
> W:=map(evalf,W(x)); WL:=subs(x=l,W):
> plot(W/WL ,para,title = "\nNormierte Verschiebung\n",labels = ["x", "W(x)"]);
# Daraus folgt die rein imaginäre Eigenkreisfrequenz für eine Druckkraft zwischen P1 = Pkr  und P2
> omega:=evalf(kB^2*f);
# Hinweis: Mit der komplexen Eigenkreisfrequenz wächst die Zeitfunktion T(t) = A1 cos(Ét) +A2 sin(Ét)  wegen sin(axI) = Isinh(ax) und cos(axI) = cosh(ax) über alle Grenzen. 
