# Kapitel 2: 
# Ein einfaches Diskretisierungsverfahren für den Dehnstab
# © 2017  Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
> interface(displayprecision = 5); restart: with(LinearAlgebra): with(plots): with(plottools): Digits:=15:
> Proc_Calc_03:=proc(Daten,T)
> #----------------------------------------------------------------
> #Eingabe:
> # Daten[1,i]: Elementdaten für das Element i
> #               Daten[1,i,1]: rho: Dichte
> #               Daten[1,i,2]:   E: E-modul
> #               Daten[1,i,3]:   A: Flächeninhalt
> #               Daten[1,i,4]:   l: Elementlänge
> # Daten[2,j]: Daten für den (k-m-c)-Schwinger j 
> #               Daten[2,j,1]:   j: Systemknotennummer
> #               Daten[2,j,2]:  kF: Federkonstante 
> #               Daten[2,j,3]:   m: Masse
> #               Daten[2,j,4]:   c: Dämpfungskonstante 
> #          T: Simulationszeit
> #Ausgabe:
> #             M: Massenmatrix
> #             K: Steifigkeitsmatrix
> #             C: Dämpfungsmatrix (wenn flag = 1)
> #            ©2: Massennormalisierte Steifigkeitsmatrix
> #            EW: Vektor der Eigenwerte
> #           Phi: Matrix der Eigenvektoren
> #          u(t): Vektor der Knotenverschiebungen
> #          v(t): Vektor der Knotengeschwindigkeiten
> #          F(t): Vektor der Stabkräfte
> #-------------------------------------------------------------------
> local n,MT,i,j,E,N,flag,eps,EVi,z0,u0,u0p,ej,Tej,caj,saj,p1,p2,p3,Z,EXZ,cd,zdvt,uv,para;
> global nsta,nkno,nkmc,M,ksta,F,WMhp,WMhm,K,C,Delta,AM,Omegaq,EW,PHI,A,EX,u,v,al;
> description "Ein einfaches Diskretisierungsverfahren für den  Dehnstab";
> nsta:=nops(Daten[1]);  #Anzahl der Stababschnitte
> nkno:=nsta+1;          #Anzahl der Knoten
> nkmc:=nops(Daten[2]);  #Anzahl der (k-m-c)-Schwinger
> u0 :=convert(Daten[3,1],Vector); #Anfangsauslenkung
> u0p:=convert(Daten[3,2],Vector); #Anfangsheschwindigkeit
> M:=Matrix(nkno,nkno);  WMhp:=Matrix(nkno,nkno); WMhm:=Matrix(nkno,nkno);
> K:=Matrix(nkno,nkno,shape = symmetric); ksta:=Vector(nsta);
> C:=Matrix(nkno,nkno);  E:=IdentityMatrix(nkno); N:=Matrix(nkno,nkno,shape = zero);
> AM:=Vector(nkno); al:=Vector(nkno); F:=Vector(nsta);
> for n to nsta do
>   MT:=0.5*Daten[1,n,1]*Daten[1,n,3]*Daten[1,n,4];
>   M[n  ,  n]:= M[n  ,  n] + MT ;
>   M[n+1,n+1]:= M[n+1,n+1] + MT ;
>   ksta[n]   := Daten[1,n,2]*Daten[1,n,3]/Daten[1,n,4];
> end do;
> #Aufbau der Steifigkeitsmatrix 
> for n from 2 to nkno-1 do
>   K[n,  n]:= ksta[n-1]+ksta[n];
>   K[n,n+1]:=-ksta[n];
> end do;
> K[1,1]:= ksta[1]; K[1,2]:=-ksta[1]; K[nkno,nkno]:= ksta[nsta];
> for n to nkmc do      #Anteile der (k-m-c)-Schwinger an der Steifigkeits-und Massenmatrix
>   i:= Daten[2,n,1];
>   K[i,i]:= K[i,i] + Daten[2,n,2];
>   M[i,i]:= M[i,i] + Daten[2,n,3];
> end do;
> print(`Massenmatrix M = `,M);
> print(`Steifigkeitsmatrix K = `,K);
> #Berechnung der Eigenwerte und Eigenvektoren
> for n to nkno do 
>   WMhm[n,n]:=1/sqrt(M[n,n]);
>   WMhp[n,n]:=  sqrt(M[n,n]);
> end do;
> Omegaq:=WMhm.K.WMhm; 
> print(`Massennormalisierte Steifigkeitsmatrix ©2  = `,Omegaq);
> flag:=0; eps:=10.^(-Digits+2);
> for n to nkmc do 
>   if abs(Daten[2,n,4]) <> 0 then flag:=1; end if;
> end do;
> if flag = 0 then
> EW,PHI:=LinearAlgebra[Eigenvectors](K,M);   #Das allgemeine Eigenwertproblem
> EW:= LinearAlgebra[Map](sqrt,evalc(Re(EW))): PHI:=evalc(Re(PHI));
> #Normierung der Eigenvektoren auf die Länge 1
> for i to nkno do EVi:=LinearAlgebra[Normalize](LinearAlgebra[Column](PHI,[i]),Euclidean);
> for j to nkno do PHI[j,i]:=EVi[j]; end do; end do;
> print(`Eigenkreisfrequenzen É und Eigenvektoren ¦ = `,EW,PHI);
> #Lösung des Anfangswertproblems
> u:=0;
> for j to nkno do
>   ej   :=LinearAlgebra[Column](PHI,[j]);
>   Tej  :=LinearAlgebra[Transpose](ej);
>   AM[j]:=sqrt((Tej.M.u0)^2 + (Tej.M.u0p)^2/EW[j]^2)/(Tej.M.ej);
>   caj  :=Tej.M.u0/(Tej.M.ej); saj:=Tej.M.u0p/(Tej.M.ej)/EW[j];
>   al[j]:=arctan(saj,caj);
>   u    := u + AM[j]*ej*cos(EW[j]*t - al[j]);
> end do;
>   v    := map(diff,u,t); 
> else      #Aufbau der Dämpfungsmatrix, hier nur aus (k-m-c)-Schwingern
>   for n to nkmc do
>     i:= Daten[2,n,1];
>     C[i,i]:= Daten[2,n,4];
>   end do;
> Delta:=0.5*WMhm.C.WMhm; 
> print(`Massennormalisierte Dämpfungsmatrix Delta = `,Delta);
> A:= Matrix(2*nkno,2*nkno,[[N,E],[-Omegaq,-2*Delta]]);
> EW,PHI:=LinearAlgebra[Eigenvectors](A); #Das spezielle Eigenwertproblem
> #Normierung der Eigenvektoren auf die Länge 1
> for i to nkno do EVi:=LinearAlgebra[Normalize](LinearAlgebra[Column](PHI,[i]),Euclidean);
> for j to nkno do PHI[j,i]:=EVi[j]; end do; end do;
> print(`Eigenwerte É = `,EW);
> print(`Eigenvektoren ¦ = `,PHI);
> #Lösung des Anfangswertproblems
> Z   :=Matrix(EW, shape = diagonal);
> EXZ :=LinearAlgebra[MatrixExponential](Z,t);
> u0  :=Vector([Daten[3,1],Daten[3,2]]);
> z0  :=Matrix([[WMhp,N],[N,WMhp]]).u0;
> cd  :=MatrixInverse(PHI).z0;
> zdvt:=PHI.EXZ.cd:
> uv  :=evalc(Re(Matrix([[WMhm,N],[N,WMhm]]).zdvt)):
> u   :=LinearAlgebra[SubVector](uv,[1..nkno]);
> v   :=LinearAlgebra[SubVector](uv,[nkno+1..2*nkno]);
> end if;
> #Berechnung der Stabkräfte
> for i to nsta do F[i]:=ksta[i]*(u[i+1]-u[i]); end do;
> #Grafische Ausgabe der Zustandsgrößen
> para:=title = "\nDer Stab mit diskret verteilter Masse\n",titlefont = ["ARIAL", 15],labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10],gridlines=true,axes=boxed:
> p1:= plot(u,t=0..T,para,labels = ["t", "Auslenkung u(t)"]): 
> p2:= plot(v,t=0..T,para,labels = ["t", "Geschwindigkeit v(t)"]):
> p3:= plot(F,t=0..T,para,labels = ["t", "Stabkräfte F(t)"]):
> plots[display](Matrix(1,3,[p1,p2,p3]));
> end proc:
> 
;
> Proc_Calc_04:=proc(j)
> local PHIN,li,n,s,i,p1,p2;
> global PHI,nkno;
> description "Berechnung des Eigenvektors zum j-ten Eigenwert eines Stabes mit diskret verteilter Masse";
> PHIN:=LinearAlgebra[Column](PHI,[j]);
> li:=Vector(nkno):
> for n from 2 to nkno do li[n]:=li[n-1] + Daten[1,n-1,4]: end do:
> s:=[0,PHIN[1]]:
> for i from 2 to nkno do s:=s,[li[i],PHIN[i]]: end do:
> p1:=plot([s],thickness=1,color=blue,linestyle=dash):
> p2:=pointplot([s],color=red,symbol=circle,symbolsize=20):
> plots[display](p1,p2,titlefont=["ARIAL", 15],title="\nNormierter Eigenvektor zum kleinsten Eigenwert É1\n",labeldirections = ["horizontal", "horizontal"],labelfont = ["HELVETICA", 10],labels = ["x","u"],gridlines=true,axes=boxed);
> end proc:
> 
;
> Proc_Calc_01:=proc(Daten,T,alpha)
> #----------------------------------------------------
> #Eingabe: 
> #         Daten: Elementdaten (global vereinbart)
> #             T: Simulationszeit
> #         alpha: Überhöhungsfaktor der Auslenkung
> #Ausgabe:
> #            An: Animation des Bewegungsvorganges
> #----------------------------------------------------
> local nkno,li,n,i;
> global p3,p4;
> description "Animation der Verschiebungen eines Stabes mit diskret verteilter Masse";
> nkno:= nops(Daten[1])+1;          #Anzahl der Knoten
> li:= Vector(nkno):
> for n from 2 to nkno do li[n]:=li[n-1] + Daten[1,n-1,4]; end do;
> p3:={};
> for i from 1 to nkno do
>   p3:= p3 union {plots[animate](pointplot,[[li[i]+alpha*u[i],0],symbol=solidcircle,symbolsize=20,color=blue],
>        t=0..T,frames = 150)};
> end do:
> p4:=plots[animate](polygonplot,[[[0,0],[li[nkno]+alpha*u[nkno],0]],thickness = 3,color = blue],t=0..T,frames = 150):
> display(p3,p4,title = "\nDer Stab mit diskret verteilter Masse\n", titlefont = ["ARIAL", 15], labels = ["\nVerschiebung u(t) (mit dem Faktor ± überhöht)\n", ""], labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10], axesfont = ["HELVETICA", "ROMAN", 8],gridlines=true,axes=boxed);
> end proc:
# Beispiel 2-13:
# 
# Beispiel für einen ungedämpften Schwinger mit Anfangsauslenkungen der Knoten 1 bis 5
# Zur Erfüllung allgemeiner Randbedingungen an den Stabenden wird dort jeweils ein (k-m-c)-Schwinger erwartet. Beispielsweise wird an einem eingespannten Rand die dortige Federsteifigkeit mit  kF = 10.^(Digits-2)  des Schwingers sehr hoch angesetzt.
#           stabi  =[           ñ ,               E,                    A,                  ]
> stab1:=[7850., 2.1E11,  7.60E-4,    1.50]:
> stab2:=[7850., 2.1E11,  7.60E-4,    1.50]: 
> stab3:=[2500., 0.2E11,  5.20E-4,    1.00]: 
> stab4:=[2500., 0.2E11,  5.20E-4,    1.00]:
> kF:=10.^(Digits-2):               #Sehr steife Feder, simuliert die Einspannung am linken Rand
;
#        kmci = [ i,           k,      m,   c]
> kmc1:=[1, kF, 0.,0.];   
> u0:=[0,0.003,0.006,0.008,0.01];   #Anfangswerte Auslenkung
;
> up0:=[0,0,0,0,0];                 #Anfangswerte Geschwindigkeit
;
> T:=0.02;                          #Simulationsdauer in [s]
;
> Daten:=[[stab1,stab2,stab3,stab4],[kmc1],[u0,up0]]:
> Proc_Calc_03(Daten,T);
# Grafische Darstellung des normierten Eigenvektors zum kleinsten Eigenwert É1
> Proc_Calc_04(min[index](EW));
# Animation des Bewegungsvorgangs mit dem Überhöhungsfaktor ±
> alpha:=50:
> Proc_Calc_01(Daten,T,alpha);
# Näherungslösung für Beispiel 2-6: Stab mit Endmasse (µ = 0.5)
> stab1:=[7850., 2.1E11,  7.60E-4,    2.5]:
> stab2:=[7850., 2.1E11,  7.60E-4,    2.5]:
> stab3:=[7850., 2.1E11,  7.60E-4,    2.5]:
> stab4:=[7850., 2.1E11,  7.60E-4,    2.5]:
> kmc1:=[1,  kF,    0.,0.]:
> kmc5:=[5,  0., 29.83,0.]:
> Daten:=[[stab1,stab2,stab3,stab4],[kmc1,kmc5],[u0,up0]]:
> Proc_Calc_03(Daten,T);
> alpha:=100;
> Proc_Calc_01(Daten,T,alpha);
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
# 
