# Kapitel 2: 
# Freie Stabschwingungen zusammengesetzter Stäbe 
# © 2017  Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
# 
> restart: with(LinearAlgebra): with(plots):
# Beispiel 2-10:
# 
# 
# 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 [ñ, E, A, ] mit 
# [ñ: Dichte des Materials, E: Elastizitätsmodul, A: Flächeninhalt des Stabquerschnitts, : Stablänge] benötigt. Für unser Beispiel gilt 
# der folgende Eingabedatensatz:
> DATEN:=[[7850.,2.1*10^(11),7.6*10^(-4),3.00] ,
>         [2500.,0.2*10^(11),5.2*10^(-4),2.00]]:
>  NELE:=nops(DATEN);  #Anzahl der Stabelemente
;
# Berechnung der Fortpflanzungsgeschwindigkeiten cL , der Longitudinalwellenzahlen kL und dem Longitudinal-Wellenwiderstand ZL für die einzelnen Stababschnitte:
> cL:=[]: kL:=[]: ZL:=[]: l:=[]:
> for i to NELE do
>   cL:= [op(cL),sqrt(DATEN[i][2]/DATEN[i][1])];
>   kL:= [op(kL),omega/cL[i]];
>   ZL:= [op(ZL),DATEN[i][3]*sqrt(DATEN[i][1]*DATEN[i][2])];
>   l := [op(l) ,DATEN[i][4]];
> end do:
> print(`cL = `,cL,`kL = `,kL,`ZL = `,ZL,`l = `,l);
# Zur Auswertung der Matrizen Ux und U stellen wir zwei Prozeduren zur Verfügung:
> Ux:= proc(kL,ZL)
> description "Matrix Ux";
> Matrix(2,2,[cos(kL*x),I/ZL*sin(kL*x),I*ZL*sin(kL*x),cos(kL*x)]);
> end proc:
> 
;
> U:= proc(kL,ZL,l)
> description "Matrix U";
> Matrix(2,2,[cos(kL*l),I/ZL*sin(kL*l),I*ZL*sin(kL*l),cos(kL*l)]);
> end proc:
> 
;
# Berechnung der Systemmatrix
> USYS:= IdentityMatrix(2): #Startwert für die Systemmatrix
;
> for i to NELE do
>   USYS:=evalc(U(kL[i],ZL[i],l[i]).USYS);
> end do:
> 
;
# Im Fall des links eingespannten und rechts freien Stabsystems steht die Eigenwertgleichung in USYS[2,2]. 
> EWG:=unapply(USYS[2,2],omega):
# Grafische Ausgabe der Eigenwertgleichung
> plot(EWG(omega),omega=0..12000,title = "\nEigenwertgleichung\n",labels = ["É", "EWG(É)"],titlefont = ["ARIAL", 15],  labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 12],axesfont = ["HELVETICA", "ROMAN", 10],gridlines,axes=boxed);
# In der folgenden Prozedur F wird zur Berechnung der Nullstellen die Eigenwertgleichung im Intervall [a,b] mit der Schrittweite dx ausgewertet. Ist EWG(x)*EWG(x+dx)  < 0, dann hat die Eigenwertgleichung in diesem Intervall einen Nulldurchgang. In diesem Fall wird die Nullstelle mit  fsolve  nummerisch berechnet und in der Liste Omega gespeichert.
> F:= proc(omega)
> description "Auswertung der Eigenwertgleichung";
> EWG(omega);
> end proc:
> 
;
> Proc_Calc_03:=proc(a,b,dx)
> local nmax, x, n;
> global Omega,EF; 
> description "Berechnung der Eigenwerte É";
> nmax:=round((b-a)/dx);
> x:= a; Omega:= [];
> for n to nmax do
>   if EWG(x)*EWG(x+dx) < 0 then
>     Omega:= [op(Omega),fsolve(F(z),z=x..x+dx)];
>   end if;
> x:= x + dx;
> end do;
> EF:= Omega/(2.*Pi);
> print(`Eigenkreisfrequenzen É =`,Omega);
> print(`Eigenfrequenzen f =`,EF);
> end proc:
# Wir rufen die Prozedur zur Nullstellenberechnung im Intervall [a = 1, b = 12000] mit der Schrittweite dx = 200 auf:
> Proc_Calc_03(1,12000,200);
> 
;
# Sind die Eigenfrequenzen bekannt, dann können elementweise die Eigenfunktionen für die Verschiebungen und die Normalkräfte  berechnet werden. Dazu werden folgende Startwerte benötigt:
> z01:=Vector(2,[0,1/USYS[1,2]]): vs:=IdentityMatrix(2).z01: ZM:=Matrix(2,NELE):
# Wir berechnen elementweise die Verschiebungs-Eigenfunktionen für den kleinsten Eigenwert É1 (j = 1)
> j:=1:
> for i to NELE do
>   omega:=Omega[j];
>   zx:=Ux(kL[i],ZL[i]).vs:
>   zl:= U(kL[i],ZL[i],l[i]).vs:
>   Ul:= U(kL[i],ZL[i],l[i]);
>   vs:= Ul.vs;
>   ZM[1,i]:=evalc(zx[1]); ZM[2,i]:=evalc(zx[2]); 
> end do:
> ZM;
> 
;
# Grafische Ausgabe der Verschiebungs-Eigenfunktionen
> xqa:=0: pl:={}: 
> for i to NELE do
> Z:=subs(x=xq-xqa, ZM[1,i]);
> pl :=pl union {plot(Z,xq = xqa..xqa + l[i],y=-1..1)};
> xqa:=xqa + l[i];
> end do:
> display(pl,title = "\nVerschiebungs-Eigenfunktion\n",labels = ["x", ""],titlefont = ["ARIAL", 15],  labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 12],axesfont = ["HELVETICA", "ROMAN", 10],gridlines,axes=boxed);
# Grafische Ausgabe der Normalkraft-Eigenfunktionen
> xqa:=0: pl:={}:
> for i to NELE do
> Z  :=subs(x=xq-xqa,evalc(I*ZM[2,i]));
> pl :=pl union {plot(Z,xq = xqa..xqa + l[i])};
> xqa:=xqa + l[i];
> end do:
> display(pl,title = "\nNormalkraft-Eigenfunktion\n",labels = ["x", ""],titlefont = ["ARIAL", 15],  labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 12],axesfont = ["HELVETICA", "ROMAN", 10],gridlines,axes=boxed);
> 
;
