> restart:
# HF, Übung 7:  Schwingung viskoelastischer Strukturen, Materialdämpfung
# 
# Zunächst die Stoffgleichung
# 
# Zur Bezeichnung: Die Zeitableitung von z.B.  edev   ist vorerst mit  
# dedev   bezeichnet
# 
> glkugel:=trt-3*K*tre;
> gldev:=tdev-2*(G[0]+Sum(G[j],j=1..n))*edev
> -2*eta*dedev+2*Sum(G[j]*a[j],j=1..n);
> glevol:=da[j]-beta[j]*(edev-a[j]);
> with(linalg):
# 
# Spezialisierung auf den einachsigen Spannungszustand
# 
# Alle auftretenden rotationssymmetrischen Deviatoren sind durch Angabe
# ihrer 1,1-Komponente
# vollständig bestimmt:
# 
> T:=diag(sigma,0,0);
> 
> e:=diag(epsilon,epsilon0,epsilon0);
> trt:=trace(T);
> tre:=trace(e);
> id:=diag(1,1,1);
> tdev:=evalm(T-trt/3*id);
> edev:=evalm(e-tre/3*id);
> dedev:=subs(epsilon=depsilon,epsilon0=depsilon0,eval(edev));
> a[j]:=diag(a1[j],-a1[j]/2,-a1[j]/2);
> da[j]:=subs(a1=da1,eval(a[j]));
> glkugel;
> gldeviator:=map(value,evalm(gldev));
> simplify(gldeviator[1,1]+2*gldeviator[2,2]);
> simplify(gldeviator[1,1]+2*gldeviator[3,3]);
> glevolution:=evalm(glevol);
> simplify(glevolution[1,1]+2*glevolution[2,2]);
> simplify(glevolution[1,1]+2*glevolution[3,3]);
# 
# Es verbleiben also nur drei Gleichungen, d.h. die folgenden drei
# Ausdrücke müssen gleich 0 sein:
# 
> g1:=glkugel;
> g2:=gldeviator[1,1];
> g3:=glevolution[1,1];
# 
# 
# Physikalische Bedeutung besitzen nur die Realteile der Dehnungen und
# Spannungen
# 
> epsilon:=Ep*exp(I*omega*t);
> epsilon0:=-nu*epsilon;#nu ist komplex
> a1[j]:=A[j]*exp(I*omega*t);
> depsilon:=diff(epsilon,t);
> depsilon0:=diff(epsilon0,t);
> da1[j]:=diff(a1[j],t);
> sigma:=S*exp(I*omega*t);
# 
# Aus den drei Gleichungen wird:
# 
> g1;
> g2;
> g3;
# Diese Zeitfunktionen sind genau dann für alle t gleich Null, wenn sie
# für t=0 verschwinden, also:
> fa1:=eval(subs(t=0,g1));
> fa2:=eval(subs(t=0,g2));
> fa3:=eval(subs(t=0,g3));
# 
# Lösung der Gleichungen:
# 
> loesung3:=solve({fa3},{A[j]});
> assign(loesung3);
> loesung12:=solve({fa1,fa2},{S,nu});
> assign(loesung12);
# 
# Definition des komplexen E-Moduls
# Er ist eine Funktion von omega und verknüpft die komplexe
# Spannung und Dehnung gemäß    sigma = E*epsilon
# 
> E:=S/Ep;
# 
# Physikalische Deutung der komplexen Größen:
# 
> eps:=(Epr+I*Epim)*exp(I*(omr+I*omim)*t);
> dehnung:=evalc(Re(eps));
> sig:=(Er+I*Eim)*eps;
> spannung:=evalc(Re(sig));
# 
# Einführung von Betrag und Phasenwinkel (psi) des komplexen E-Moduls
# 
> Er:=Ebetrag*cos(psi); Eim:=Ebetrag*sin(psi);
> sp:=combine(spannung,trig):
> factor(sp);
> factor(dehnung);
# Man erkennt:
# Wenn der E-Modul E komplex ist, eilt die Spannung der Dehnung voraus
# um den Winkel psi
# 
# Zahlenbeispiele dazu:
# 
# a) Stationäre Schwingung
# 
> Epr:=1: Epim:=0:Ebetrag:=1:omr:=1:
> plot([subs(omim=0,psi=Pi/10,dehnung),subs(omim=0,psi=Pi/10,spannung),
> t=0..2*Pi],labels=['epsilon','sigma']);
# 
# Die umschlossene Fläche (Hysterese-Schleife) gibt die dissipierte
# mechanische Arbeit in einem Zyklus
# 
# 
# b)Abklingende Schwingung
# 
> plot([subs(omim=0.2,psi=Pi/5,dehnung),subs(omim=0.2,psi=Pi/5,spannung)
> ,
> t=0..10*Pi],labels=['epsilon','sigma']);
# 
# Noch einmal der komplexe E-Modul:
# 
> E;
# 
# Sonderfälle:
# 
# 
# Elastisches Material
# 
> hooke:=value(subs(n=0,eta=0,[E,nu]));
> Eelast:=hooke[1];
# 
# Viskoelastische Materialien, Verhalten bei erzwungener Schwingung, 
# also reellem omega
# 
# 
# Kelvin-Material
# 
> kelvin:=value(subs(n=0,[E,nu]));
# 
# Wir legen fest, dass die Langzeit-Querkontraktion (also nu für omega
# -- > 0) den Wert 
# 1/3 besitzen soll. Das liefert die folgende Aussage über den
# Kompressionsmodul:
# 
> komprk:=solve({subs(omega=0,kelvin[2])=1/3},{K});
# Auf den elastischen Modul bezogener komplexer E-Modul, ausgedrückt
# durch die
# dimensionslos gemachte Kreisfrequenz  o:
# 
> Emod:=simplify(subs(omega=o/eta*G[0],komprk,kelvin[1]/Eelast));
> 
> plots[complexplot](Emod,o=0..200);
> Er:=evalc(Re(Emod));
> Eim:=evalc(Im(Emod));
> Ebetrag:=evalc(abs(Emod));
> Ewinkel:=evalc(argument(Emod));
> plot([Er,Eim],o=0..100,color=[red,blue]);
> plot([Ebetrag,tan(Ewinkel)],o=0..50,color=[red,blue]);
# 
# Maxwell-Material
# 
> maxwell:=value(subs(n=1,G[0]=0,eta=0,[E,nu]));
# Wir legen fest, dass die Kurzzeit-Querkontraktion (also nu für omega
# -- > Unendlich) den Wert 
# 1/3 besitzen soll. Das liefert die folgende Aussage über den
# Kompressionsmodul:
# 
> komprm:=solve({limit(maxwell[2],omega=infinity)=1/3},{K});
# 
# Auf den (jetzt mit G[1] statt G[0] berechneten)elastischen Modul
# bezogener komplexer E-Modul, ausgedrückt durch die dimensionslos
# gemachte Kreisfrequenz  o:

> Emod:=simplify(subs(omega=o*beta[1],komprm,maxwell[1]/
> subs(G[0]=G[1],Eelast)));
> plots[complexplot](Emod,o=0..20);
> Er:=evalc(Re(Emod));
> Eim:=evalc(Im(Emod));
> Ebetrag:=evalc(abs(Emod));
> Ewinkel:=evalc(argument(Emod));
> plot([Er,Eim],o=0..10,color=[red,blue]);
> plot([Ebetrag,tan(Ewinkel)],o=0..10,0..10,color=[red,blue]);
# 
# Maxwell und Kelvin parallel
# 
> parallel:=value(subs(n=1,[E,nu]));
> su:=omega=o*beta[1],beta[1]=0.2*G[0]/eta,G[1]=5.*G[0],komprk;
> Emod:=simplify(subs(su,parallel[1]/Eelast));
> plots[complexplot](Emod,o=0..800,numpoints=200);
> Er:=evalc(Re(Emod));
> Eim:=evalc(Im(Emod));
> Ebetrag:=evalc(abs(Emod));
> Ewinkel:=evalc(argument(Emod));
> plot([Er,Eim],o=0..200,color=[red,blue]);
> plot([Ebetrag,tan(Ewinkel)],o=0..200,color=[red,blue]);
# 
# Allgemeineres Material (3 Maxwell parallel neben elastischem Element,
# also 3 Relaxationszeiten)
> allgemein:=value(subs(n=3,eta=0,[E,nu]));
> su:=beta[1]=0.1*beta[2],beta[3]=10.*beta[2],omega=o*beta[2],G[1]=G[0],
> G[2]=G[0],G[3]=G[0],komprk;
> Emod:=simplify(subs(su,allgemein[1]/Eelast));
> plots[complexplot](Emod,o=0..80,numpoints=200);
> Er:=evalc(Re(Emod));
> Eim:=evalc(Im(Emod));
> Ebetrag:=evalc(abs(Emod));
> Ewinkel:=evalc(argument(Emod));
> plot([Er,Eim],o=0..80,color=[red,blue]);
> plot([Er,Eim],o=0..10,color=[red,blue]);
> plot([Ebetrag,tan(Ewinkel)],o=0..80,color=[red,blue]);
> plot([Ebetrag,tan(Ewinkel)],o=0..10,color=[red,blue]);
# 
# Biegeschwingungen
# 
# Balken auf zwei Stützen
# 
# Biegedifferentialgleichung
# 
> dgl:=Em*Iy*diff(w(x,t),x,x,x,x)+rho*Area*diff(w(x,t),t,t)=q;
# 
# Komplexer Ansatz, der die Randbedingungen  w(0)=0,w(l)=0,
# Mb(0)=0,Mb(l)=0  erfüllt:
# 
> w:=(x,t)->W*sin(m*Pi*x/l)*exp(I*omega*t);
> q:=Q*sin(m*Pi*x/l)*exp(I*omega*t);
> dgl;
> loes:=solve(dgl,{Q});
> assign(loes);
# 
# Erster Sonderfall: Freie Schwingung, keine Anregung (q=0),
# die komplexe Eigenkreisfrequenz omega ist gesucht:
> ewgleichung:=solve({Q=0},{omega});
# Das ist die Bestimmungsgleichung für die Eigenkreisfrequenz omega. 
# 
# Bei elastischem Material ist der E-Modul Em reell und konstant und
# damit omega reell (rein harmonische Eigenschwingung, omega wächst mit
# dem Quadrat der Nummer m der Oberschwingung).
# Bei viskoelastischem Material ist der E-Modul Em eine komplexe
# Funktion von omega und damit omega ebenfalls komplex (gedämpfte
# harmonische Eigenschwingung).
# 
# Zahlenbeispiele
# 
> ewg:=evalf(subs(Em=Eelast*Emod,omega=o*beta[2],komprk,
> ewgleichung[1]));
> ewgneu:=(rhs(op(ewg))/beta[2])^2-(lhs(op(ewg))/beta[2])^2;
> G[0]:=lambda*rho*Area*l^4*beta[2]^2/Iy;
> ewgneu;
> ewgneu1:=simplify(ewgneu);
> ewgneu2:=collect(op(2,ewgneu1),o);#Das ist die Klammer im Zähler
# 
# Ermittlung der komplexen Eigenfrequenzen
# 
# Man erkennt im Folgenden:
# 
# Es gibt drei imaginäre Werte. Sie beschreiben Kriechmoden. Die
# zugehörigen Kehrwerte
# (die Retardationszeiten) verhalten sich etwa wie die
# Relaxationszeiten.
# 
# Es gibt zwei komplexe Werte, deren Realteile sich im Vorzeichen
# unterscheiden.
# Der Imaginärteil ist positiv, d.h. die Schwingung ist gedämpft
# 
> lambda:=0.0000001: # Zahlenbeispiel
> loes:=fsolve(subs(m=1,ewgneu2),o);
> loes:=fsolve(subs(m=2,ewgneu2),o);
> loes:=fsolve(subs(m=10,ewgneu2),o);
> loes:=fsolve(subs(m=100,ewgneu2),o);
> lambda:='lambda':
# Zweiter Sonderfall: Erzwungene Schwingung (omega reell, vorgegeben)
# 
# Benötigte Amplitude der erregenden Linienkraft:
# 
> Q;
# 
# Vergrößerungsfunktion:
# 
> V:=abs(W/Q*rho*Area*beta[2]^2);
> V1:=simplify(subs(Em=Eelast*Emod,omega=o*beta[2],komprk,V));
# 
# Dieser Betrag eines komplexen Ausdrucks ist natürlich rein reell,
# seine Auswertung geschieht so:
# 
> V2:=evalc(V1);
> lambda:=0.0000001: # Wieder das Zahlenbeispiel
> plot(subs(m=1,V2),o=0..0.02);
> plot(subs(m=2,V2),o=0..0.1);
> plot(subs(m=10,V2),o=0..3);
> plot(subs(m=100,V2),o=0..150);
> 
>