# Kapitel 7_e:
# Erzwungene Bewegungen des viskos gedämpften
# Einmassenschwingers
# © 2015  Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
> restart: with(LinearAlgebra):
# 
# Die Anwendung des Netonschen Grundgesetzes auf die freigeschnittene Masse m liefert die Differenzialgleichung (fE(t) = f0 cos(©t))
> ode2:=diff(x(t),t,t)+2*D0*omega*diff(x(t),t)+omega^2 *x(t)-fE(t);
# Wir unterstellen eine harmonisch-periodische Erregung (f0 = F0/m)
> fE(t):=f0*cos(Omega*t);
> loe:=dsolve(ode2,x(t));assign(%);
> x:=unapply(x(t),t);
# und für die Geschwindigkeit folgt
> v:=unapply(diff(x(t),t),t);
> x:='x': v:='v':
# Wir lösen nun das vollständige Anfangswertproblem
> AB:=x(0)=x0,D(x)(0)=v0;
> loe:=dsolve({ode2,AB},x(t)); assign(%);
> x:=unapply(x(t),t,x0,v0);
> v:=unapply(diff(x(t,x0,v0),t),t,x0,v0);
# Im Fall homogener Anfangsbedingungen verbleiben
> x(t,0,0);v(t,0,0);
# Wir stellen die Auslenkung x(t) für homogene Anfangsbedingungen grafisch dar und wählen
> f0:=1; omega:=1/2; D0:=0.1; Omega:=sqrt(3);
> plot(x(t,0,0),t=0..30,title = "\nErzwungene Bewegung\n", titlefont = ["ARIAL", 15], labels = ["Zeit t", "Auslenkung x(t)"], labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10], axesfont = ["HELVETICA", "ROMAN", 8],axes=boxed,gridlines=true);
# Die Vergrößerungsfunktion V(·,D0) mit · = ©/É
> f0:='f0': omega:='omega': D0:='D0': Omega:='Omega':
> V:=unapply(1/sqrt((1-eta^2)^2+(2*D0*eta)^2),eta,D0);
# Wir geben die Vergrößerungsfunktion für verschiedene Dämpfungsgrade D0 in Schritten von 0.2 beginnend bei · = 0 aus:
> WD0:=[0,0.2,0.4,0.6,0.8,1.,1.2]:
> plot([seq(V(eta,WD0[i]),i=1..nops(WD0))],eta=0..3,y=0..3,title = "\nVergrößerungsfunktionen\n", titlefont = ["ARIAL", 15], labels = ["eta", "V(eta)"], labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10], legend = ["D0=0.0","D0=0.2","D0=0.4","D0=0.6","D0=0.8","D0=1.0","D0=1.2"],legendstyle = [font = ["HELVETICA", 9], location = right],axesfont = ["HELVETICA", "ROMAN", 8],axes=boxed,gridlines=true);
# Phasenfrequenzgang Æ(·,D0)
> sin(phi):=2*D0*eta/(V(eta,D0)*((1-eta^2)^2+(2*D0*eta)^2));cos(phi):=(1-eta^2)/(V(eta,D0)*((1-eta^2)^2+(2*D0*eta)^2));
> phi:=unapply(arctan(sin(phi),cos(phi)),eta,D0);
> plot([seq(phi(eta,WD0[i]),i=1..nops(WD0))],eta=0..3,title = "\nPhasenfrequenzgang\n", titlefont = ["ARIAL", 15], labels = ["eta", "phi(eta)"], labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10], legend = ["D0=0.0","D0=0.2","D0=0.4","D0=0.6","D0=0.8","D0=1.0","D0=1.2"],legendstyle = [font = ["HELVETICA", 9], location = right], axesfont = ["HELVETICA", "ROMAN", 8],axes=boxed,gridlines=true);
# Allgemeine Erregerkraftfunktionen
# Wir führen zunächst die Gewichtsfunktion g(t-Ä) ein. 
> delta:=D0*omega; lambda:=omega*sqrt(D0^2-1);
> g:=(exp(-(t-tau)*(delta-lambda))-exp(-(t-tau)*(delta+lambda)))/(2*lambda);
# Abhängig vom Dämpfungsgrad erhalten wir  folgende Gewichtsfunktionen 
# 1. Ungedämpfte Schwingung (D0 = 0)
> g_ung:=simplify(subs(D0=0,g));
# 2. Schwache Dämpfung (D0 < 1)
> omegad:=omega*sqrt(1-D0^2);
> g_schw:=exp(-delta*(t-tau))*sin(omegad*(t-tau))/omegad;
# 3. Der Grenzfall (D0 = 1)
> g_grenz:=simplify(limit(g,D0 = 1));
# Ein Partikularintegral der inhomogenen Differenzialgleichung kann dann in der Form:  xp (t) = (&int;)[tau=t`0`]^(t)f(tau)g(t-tau)&DifferentialD;tau;   dargestellt werden. 
# Wir zeigen das am Beispiel f(Ä) = f0Åcos(©Ä) und wählen t0 = 0:
> f:=unapply(f0*cos(Omega*tau),tau);
> xp:=int(f(tau)*g_schw,tau =0..t);
> D0:=0.2; omega:=1; f0:=1; Omega:= 2; 
> plot(xp,t=0..30,title = "\nErzwungene Bewegung (Partikularintegral)\n", titlefont = ["ARIAL", 15], labels = ["Zeit t", "Auslenkung x(t)"], labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10], axesfont = ["HELVETICA", "ROMAN", 8],axes=boxed,gridlines=true);
# Befindet sich das System zum Zeitpunkt t = t0 (hier t0 = 0) nicht in Ruhe, dann ist dem Partikularintegral ein Integral der homogenen Differenzialgleichung hinzuzufügen. Mit den beiden Konstanten aus der homogenen Lösung lässt sich dann die Gesamtlösung an beliebige inhomogene Anfangsbedingungen anpassen. 
> D0:='D0': omega:='omega': Omega:='Omega': 
# Beispiel 7-3:
# 
# Das System in in der obigen Abbildung wird zum Zeitpunkt t = 0 aus der Ruhe heraus durch eine sprunghafte Belastung der Intensität F0 beansprucht. Solche Belastungen können beispielsweise dann auftreten, wenn sich durch den Ausfall von Unterstützungen in Baukonstruktionen das statische System plötzlich ändert (Erdbeben).
> f:=unapply(F0/m,tau);
# Je nach Dämpfungsgrad erhalten wir folgende partikuläre Lösungen:
# 1. Ungedämpfte Schwingung (D0 = 0) 
> xp_ung:=simplify(int(f(tau)*g_ung,tau =0..t));
# 2. Schwache Dämpfung (D0 < 1)
> xp_schw:=unapply(simplify(int(f(tau)*g_schw,tau =0..t)),D0);
# 3. Starke Dämpfung (D0 > 1)
> xp_stark:=unapply(simplify(int(f(tau)*g,tau =0..t)),D0);
# 4. Der Grenzfall (D0 = 1)
> xp_grenz:=simplify(int(f(tau)*g_grenz,tau =0..t));
> omega:=1; m:=1; F0:=1; Omega:= 2; 
> xp_stark;
> plot([xp_ung,xp_schw(0.2),xp_stark(1.2),xp_grenz],t=0..15,title = "\nErzwungene Bewegungen (Sprungfunktion) \n", titlefont = ["ARIAL", 15], labels = ["t", "xp(t)"], labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10], legend = ["D0=0.0","D0=0.2","D0=1.2","D0=1.0"],legendstyle = [font = ["HELVETICA", 9], location = bottom], axesfont = ["HELVETICA", "ROMAN", 8],axes=boxed,gridlines=true);
# Beispiel 7-4:
# 
# Das System in in der obigen Abbildung wird durch einen Rechteckstoß der Intensität F0 und der Dauer tF belastet. Gesucht ist die dynamische Antwort, wenn das System zum Zeitpunkt t = 0 in Ruhe war. 
# Um eine einheitliche Darstellung der Belastung über den gesamten Wertebereich von t zu erhalten, kann die Sprungfunktion formal mittels der Heaviside-Funktion  
# ausgedrückt werden. Denken wir uns die Erregerkraft durch zwei zeitversetzte Sprungfunktionen zusammengesetzt, dann erhalten wir 
> f:=unapply(f0*(Heaviside(tau)-Heaviside(tau-tF)),tau);
> xp:=int(f(tau)*g_schw,tau=0..t);
> D0:=0.05; tF:=2.; f0:=1;
# Wir stellen das Ergebnis grafisch dar
> plot(xp,t=0..30,title = "\nPartikularintegral Rechteckstoß\n", titlefont = ["ARIAL", 15], labels = ["Zeit t", "Auslenkung x(t)"], labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10], axesfont = ["HELVETICA", "ROMAN", 8],axes=boxed,gridlines=true);
> 
;
