# Kapitel 7_i:
# Freie viskos gedämpfte Schwingungen
# mit allgemein n Freiheitsgraden
# © 2015  Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
> interface(displayprecision =4): 
> restart: with(LinearAlgebra):
# Wir stellen eine Prozedur zur Verfügung, mit der die freien viskos gedämpften Bewegungen eines Schwingers mit allgemein n Freiheitgraden automatisiert ermittelt werden können.
> 
;
> Proc_Calc_31:=proc(M::Matrix,C::Matrix,K::Matrix,AB::listlist,tE)
> #----------------------------------------------------------------------
> # Eingabe: 
> #            M: Massenmatrix        [n x n]
> #            C: Dämpfungsmatrix     [n x n]
> #            K: Steifigkeitsmatrix  [n x n]
> #           AB: Anfangsbedingungen: [[x10,v10],[x20,v20],...,[xn0,vn0]]
> #           tE: Ende der Simulationszeit
> # Ausgabe:      (global vereinbart)
> #           EW: Eigenwerte
> #          Phi: Eigenvektoren
> #            x: Auslenkungen
> #            v: Geschwindigkeiten
> #            a: Beschleunigungen
> #----------------------------------------------------------------------          
> local n,DMpeh, DMmeh, EW1, EV1, EVI, DM, i,Delta,OMQ, DSM, zeta,Poly,k,A,G,gl,Unb,loe,m,qd0,qdp0,pd0,pdp0,W1,W1S,Z1,Z0,CV,CVT,pd,para,p1,p2,p3,PM,Mhp,Mhm;
> global EW,Phi,x,v,a;
> description "Freie viskos gedämpfte Schwingungen mit endlich vielen Freiheitsgraden";
> n      :=LinearAlgebra[RowDimension](M);
> DMpeh  :=Matrix(n,n); DMmeh  :=Matrix(n,n); Phi:=Matrix(n,2*n);
> EW1,EV1:=Eigenvectors(M);
> EW1    :=Re(EW1); EV1:=Re(EV1);
> EVI    :=MatrixInverse(EV1);
> DM     :=EVI.M.EV1;
> for i to n do
>   DMpeh[i,i]:= sqrt(DM[i,i]);
>   DMmeh[i,i]:= 1/DMpeh[i,i];
> end do;
> Mhp:=EV1.DMpeh.EVI;
> Mhm:=EV1.DMmeh.EVI;
> Delta:=0.5*Mhm.C.Mhm;  OMQ:=Mhm.K.Mhm;
> #Berechnung der Eigenwerte
> DSM :=zeta^2*IdentityMatrix(n) + 2.*zeta*Delta + OMQ;
> Poly:=factor(LinearAlgebra[Determinant](DSM));
> EW  :=solve(Poly,zeta);
> for k to 2*n do
>   A  :=Vector(n,{seq((i)=c[i,k],i=1..n)}): A[n]:=1:
>   G  :=subs(zeta=EW[k],DSM).A;
>   gl :=[seq(G[k],k=1..n-1)];
>   Unb:=[seq(A[i],i=1..n-1)];
>   loe:=solve(gl,Unb); assign(%);
>   A  :=map(eval,LinearAlgebra[Normalize](A,Euclidean));
>   for m to n do
>     Phi[m,k]:=A[m]; 
>   end do;                        
> end do;
> # Transformation der Anfangswerte
> qd0 :=Vector([seq(AB[k][1],k=1..n)]);
> qdp0:=Vector([seq(AB[k][2],k=1..n)]);
> pd0 :=Mhp.qd0;
> pdp0:=Mhp.qdp0;
> # Berechnung der Matrizen Z0 und Z1
> W1:=Matrix(n): Z1:=Matrix(n):
> for i to 2*n do
>   CV := LinearAlgebra[Column](Phi,i);
>   CVT:= LinearAlgebra[Transpose](CV);
>   W1S:= CV.CVT/(CVT.OMQ.CV-EW[i]*EW[i]*(CVT.CV))*exp(EW[i]*t);
>   W1 := W1 + W1S;
>   Z1 := Z1 - EW[i]*W1S;
> end do:
> Z0:=W1.OMQ:
> # Berechnung von pd
> pd:=evalc(Z0.pd0 + Z1.pdp0);
> x:=Mhm.pd;         #Auslenkung
> v:=map(diff,x,t);  #Geschwindigkeit
> a:=map(diff,v,t);  #Beschleunigung
> # Grafische Darstellung der Ergebnisse
> para:=axes=boxed,gridlines, titlefont = ["ARIAL", 15], labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10],axesfont=["HELVETICA", "ROMAN", 8], color=[black,blue,gold]:
> p1:=plot(x,t=0..tE,para,title = "\nAuslenkungen\n"     , labels = ["Zeit t", "Auslenkung"]):
> p2:=plot(v,t=0..tE,para,title = "\nGeschwindigkeiten\n", labels = ["Zeit t", "Geschwindigkeit"]):
> p3:=plot(a,t=0..tE,para,title = "\nBeschleunigungen\n" , labels = ["Zeit t", "Beschleunigung"]):
> PM:=Array(1..3): PM[1]:=p1: PM[2]:=p2: PM[3]:=p3:
> plots[display](PM);
> end proc:
> 
;
# Beispiel
# 
# Für den obigen Schwinger mit n = 3 Freiehitsgraden sind sämtliche Zustandsgrößen zu berechnen
# 
> M:=Matrix([[m1,0,0],[0,m2,0],[0,0,m3]]);               #Massenmatrix
;
> C:=Matrix([[c1+c2,-c2,0],[-c2,c2+c3,-c3],[0,-c3,c3]]); #Dämpfungsmatrix
;
> K:=Matrix([[k1+k2,-k2,0],[-k2,k2+k3,-k3],[0,-k3,k3]]); #Steifigkeitsmatrix
;
# Systemwerte:
> m1:=20.; c1:=10.; k1:=1600.;
> m2:= 5.; c2:= 3.; k2:=1000.;
> m3:=10.; c3:= 5.; k3:=1200.;
> M; C; K;
# Anfangsbedingungen:
> AB:=[[x10,v10],[x20,v20],[x30,v30]]; tE:=3;
> x10:=0.1;  v10:=0.;
> x20:=0.;  v20:=0.;
> x30:=0.;  v30:=0.;
> Proc_Calc_31(M,C,K,AB,tE);
> 
;
