# Kapitel 2: 
# Freie Torsionsschwingungen zusammengesetzter Stäbe 
# © 2017  Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
# 
> restart: with(LinearAlgebra): with(plots):
# Beispiel 2-18:
# 
# 
# Für das in der obigen Abbildung skizzierte Stabsystem sind die Eigenwerte, Eigenkreisfrequenzen und Eigenfunktionen zu bestimmen. 
# Die Werte des Beispiels werden in einer Liste von Listen (listlist) geführt. Für jeden Stab wird eine Liste in der Form [ñ, G, Ip, ] mit 
# [ñ: Dichte des Materials, G: Schubmodul, Ip: Polares Flächenträgheitsmoment des Stabquerschnitts, : Stablänge] 
# benötigt. Für unser Beispiel gilt der folgende Eingabedatensatz:
> DATEN:=[[rho,G,  IP1,  a] , 
>         [rho,G,IP1/2,2*a]]:
> NELE:=nops(DATEN);  #Anzahl der Stabelemente
;
# Berechnung der Fortpflanzungsgeschwindigkeiten cT , der Longitudinalwellenzahlen kT und dem Longitudinal-Wellenwiderstand ZT für die einzelnen Stababschnitte:
> cT:=[]: kT:=[]: ZT:=[]: l:=[]:
> for i to NELE do
>   cT:= [op(cT),sqrt(DATEN[i][2]/DATEN[i][1])];
>   kT:= [op(kT),omega/cT[i]];
>   ZT:= [op(ZT),DATEN[i][3]*sqrt(DATEN[i][1]*DATEN[i][2])];
>   l := [op(l) ,DATEN[i][4]];
> end do:
> print(`cT = `,cT);
> print(`kT = `,kT);
> print(`ZT = `,ZT);
> print(`l  = `, l);
# Zur Auswertung der Matrizen Ux und U stellen wir zwei Prozeduren zur Verfügung:
> Ux:= proc(kT,ZT)
> description "Matrix Ux";
> Matrix(2,2,[cos(kT*x),I/ZT*sin(kT*x),I*ZT*sin(kT*x),cos(kT*x)]);
> end proc:
> 
;
> U:= proc(kT,ZT,l)
> description "Matrix U";
> Matrix(2,2,[cos(kT*l),I/ZT*sin(kT*l),I*ZT*sin(kT*l),cos(kT*l)]);
> end proc:
> 
;
# Berechnung der Systemmatrix
> USYS:= IdentityMatrix(2): #Startwert für die Systemmatrix ist die Einheitsmatrix
;
> for i to NELE do
>   USYS:=eval(U(kT[i],ZT[i],l[i]).USYS)  assuming rho::positive, G::positive, a::positive, b::positive;
> end do:
> USYS; #Systemmatrix
;
# Im Fall des beidseitig frei drehbaren Stabsystems steht die Eigenwertgleichung in USYS[2,1]. Wir führen folgende Abkürzung ein:
> omega:=lambda*sqrt(G/rho)/a;
> 
;
> EWG:=simplify(USYS[2,1],trig);
> f1:= cos(lambda)^2-1/3; f2:=sin(lambda);
> EWG:=unapply(f1*f2,lambda);
# Wir stellen die Eigenwertgleichung mit ihren beiden Faktoren grafisch dar:
> plot([f1,f2,EWG(lambda)],lambda=0..3.*Pi,title="\nEigenwertgleichung EWG = f1*f2",titlefont = ["ARIAL", 15],legend=["f1","f2","EWG"],color=[red,blue,black],gridlines=true, linestyle = [dashdot, longdash,solid ], axesfont = ["HELVETICA", "ROMAN", 8], legendstyle = [font = ["HELVETICA", 9], location = bottom],axes=boxed);
# Der kleinste Eigenwert ist (s.h. obige Grafik):
> lambda1:=fsolve(EWG(lambda),lambda=0.8..1.0);
# Und für die kleinste Eigenkreisfrequenz folgt mit den Werten des Beispiels:
> G:=0.81E11; rho:=7850.; a:=2.;
> omega1:=lambda1*sqrt(G/rho)/a;
> 
;
