# Lektion 10
# Lösung von Differenzialgleichungen
# © 2017  Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
> restart: with(plots):
# Wir beschränken uns auf gewöhnliche Differenzialgleichungen n-ter Ordnung, das sind Gleichungen vom Typ  F[x, y(x), y'(x), ... ,y(n)(x)], in denen Ableitungen einer unbekannten Funktion y = y(x) bis zur n-ten Ordnung auftreten. Als Beispiel betrachten wir die Differenzialgleichung 1-ter Ordnung  N'(x) = -n(x) eines Dehnstabes. Darin bezeichnen N(x) die unbekannte Normalkraft und n(x) die in der Stabachse wirkende Linienkraftschüttung. Gesucht wird also eine Funktion N(x), deren Ableitung -n(x) entspricht. Die Lösung dieser Differenzialgleichung erfolgt durch Quadratur. Die äußere Normalkraftschüttung sei beispielsweise mit  n(x) = n + (nr - n) x/  linear verteilt. 
> ode1:=D(N)(x)=-n(x);
> n:=unapply(nl+(nr-nl)*x/l,x);
# Die DGL wird mittels des Befehls dsolve gelöst:
> dsolve(ode1);
# Die allgemeine Lösung enthält noch eine freie Konstante (hier _C1), die aus einer Randbedingung (RB) bestimmt wird. Ist beispielsweise der Stab bei x = 0 eingespannt und am rechten Ende bei x =   frei gelagert, dann muss dort die Normalkraft verschwinden. Die Einspannbedingung am linken Ende (dort ist die Verschiebung gleich Null)  kann als RB nicht genutzt werden. Wir fordern deshalb
> RB:=N(l)=0;
# Die DGL wird nun abermals unter Vorgabe der Randbedingung gelöst.
> loe_1:=dsolve({ode1,RB});
# Mit dem Befehl odetest  können wir das Ergebnis unter der angegebenen Randbedingung überprüfen. Wird bei diesem Test eine Null zurückgegeben, dann ist die Lösung  gefunden.
> odetest({loe_1,RB},ode1);
> N:=unapply(rhs(loe_1),x);
# Faktorisieren wir diesen Ausdruck, dann erkennen wir mit x =   die Randbedingung am rechten Ende:
> factor(N(x));
# Am linken eingespannten Ende ist
> N0:=factor(N(0));
# Für den weiteren Rechengang wählen wir folgende Systemparameter:
> nl:=1: nr:=2: l:=1:
# und stellen das Ergebnis grafisch dar. 
> par:=x=0..l,titlefont=[HELVETICA,14],color="NavyBlue",thickness=2,filled=[color="Blue",transparency =.5],gridlines=true,axes=framed:
> p1:=plot(n(x),par,title="\nLinienkraft n(x)\n"):
> p2:=plot(N(x),par,title="\nNormalkraft N(x)\n"):
> plotmatrix:=Matrix(1,2,[p1,p2]):
> display(plotmatrix);
# Die größte Normalkraft tritt an der Stelle x = 0 auf
> N(0);
> n:='n':N:='N':l:='l':nl:='nl':nr:='nr':
# Zwischen der Normalkraft N(x) und der Änderung der Stabachsverschiebung u(x)  besteht der Zusammenhang:  N(x) = EA u'(x), wobei das Produkt  EA  Dehnsteifigkeit genannt wird. 
> ode_N:=D(N)(x)=-n(x); ode_u:=D(u)(x)=1/EA*N(x);
# Wir haben also mit ode_N und ode_u ein System von zwei Differenzialgleichungen 
> sys:={ode_N,ode_u};
# unter  den Randbedingungen des links eingespannten (u = 0) und rechts frei gelagerten Stabes (N = 0) zu lösen.
> RB := {u(0)=0,N(l)=0};
# Wir wählen wieder eine linear veränderliche Linienkraftbelastung. 
> n:=unapply(nl+(nr-nl)*x/l,x);
# und lösen das Differenzialgleichungssystem:
> loe:=dsolve(`union`(sys,RB),{u(x),N(x)}); assign(%);
# Mit dem Befehl assign wurden den Variablen N(x) und u(x) die zuletzt in der Lösung rechts stehenden Ausdrücke als Werte zugewiesen.
> N:=unapply(simplify(N(x)),x);u:=unapply(simplify(u(x)),x);
# Die größte Normalkraft tritt wieder bei x = 0 auf, und die größte Verschiebung ergibt sich am freien Ende bei x =  
> simplify(N(0));simplify(u(l));
# Wir wenden uns nun der Schnittlastenermittlung eines statisch bestimmt  gelagerten Balkens zu. Die Differenzialgleichungen für die Querkraft und das Biegemoment (Schnittlastendifferenzialgleichungen) lauten: Q'(x) = -q(x) und M'(x) = Q(x).  Hierbei bezeichen Q(x) die Querkraft, M(x) das Biegemoment und q(x) die Linienkraft senkrecht zur Balkenachse, die hier linear verteilt angenommen wird:
> odeQ:=D(Q)(x)=-q(x); odeM:=D(M)(x)=Q(x);
> q:=unapply(ql+(qr-ql)*x/l,x);
> sys:={odeQ,odeM};
# Wir formulieren einige Randbedingungen für den Träger auf zwei Stützen:
# 1.) Randbedingungen für den bei x = 0  eingespannten und bei  x =   frei gelagerten Balken  (Kragträger)
> RB_01 := {Q(l)=0,M(l)=0};
# 2.) Randbedingungen für den bei x = 0  frei gelagerten  und bei  x =   eingespannten Balken  (Kragträger)
> RB_02 := {Q(0)=0,M(0)=0};
# 3.) Randbedingungen für den beidseits drehbar gelagerten Balken (Träger auf zwei Stützen)
> RB_03 := {M(0)=0,M(l)=0};
# und lösen das System unter den Randbedingungen RB_03  (Träger auf zwei Stützen)
> loe_03:=dsolve(`union`(sys,RB_03),{Q(x),M(x)});assign(%);
> Q:=unapply(simplify(Q(x)),x); M:=unapply(factor(M(x)),x);
# Für den weiteren Rechengang wählen wir 
> ql:=0:qr:=1:l:=1:
# und geben die Schnittlasten für den Träger auf zwei Stützen  mittels der folgenden Prozedur grafisch aus. 
> grafik:=proc(f,x,xmin,xmax,TEXT::string)
> description "Plot einer Zustandsgröße für den Balken";
> plot(f,x=xmin..xmax,title=TEXT,titlefont=[HELVETICA,14],color="NavyBlue",thickness=2,filled=[color="Blue",transparency =.5],gridlines=true,axes=boxed);
> end proc:
> grafik(Q(x),x,0,1,"\nQuerkraft Q(x)\n");
# Die Nullstelle der Querkraft liegt  im Intervall x = (0.5...0.6). 
> xnull:=fsolve(Q(x),x=0.5..0.6);
# An dieser Stelle wird das Biegemoment extremal
> maxM:=M(xnull);
> grafik(M(x),x,0,1,"\nBiegemoment M(x)\n");
# Auch unstetige Belastungen können von Maple verarbeitet werden, etwa partiell wirkende Linienkraftbelastungen mit Unstetigkeiten im Lösungsgebiet.
> l:='l':M:='M':Q:='Q':
> q:=unapply(q0*(Heaviside(x-a)-Heaviside(x-b)),x);
> q(x):=subs(a=1/4,b=1/2,q0=1,q(x));
> grafik(q(x),x,0,1,"\nBelastung q(x)\n");
> sys:={odeQ,odeM};
# Wir wählen wieder die Randbedingungen für den beidseits drehbar gelagerten Balken (RB_03)
> loe_02:=dsolve(`union`(sys,RB_03),{Q(x),M(x)}):assign(%):
> Q:=unapply(simplify(Q(x)),x);
> M:=unapply(simplify(M(x)),x);
# Wir stellen die Ergebnisse grafisch dar und wählen dazu folgende Systemwerte:
> q0:=1: l:=1: a:=1/4: b:=1/2:
> grafik(Q(x),x,0,l,"\nQuerkraft Q(x)\n");
# Der Nulldurchgang der Querkraft liegt im Interval x = (0.3 ... 0.5)
> xnull:=fsolve(Q(x),x=0.3..0.5);
# An dieser Stelle wird das Biegemoment extremal
> maxM:=M(xnull);
> grafik(M(x),x,0,1,"\nBiegemoment M(x)\n");
> l:='l':ql:='ql':qr:='qr':
# Wirkt auf den Balken eine Linienkraft q(x), die im Intervall  [x, xr]  linear von q  auf  qr anwächst, dann gilt:
> q:=unapply(Heaviside(x-xl)*Heaviside(xr-x)*(ql+(qr-ql)/(xr-xl)*(x-xl)),x);
> M:='M':Q:='Q':sys:={odeQ,odeM};
# Beispielrechnung für den links eingespannten Balken  (Kragträger)
> loe_01:=dsolve(`union`(sys,RB_01),{Q(x),M(x)}):assign(%):
> Q:=unapply(Q(x),x);
> M:=unapply(M(x),x);
# Wir wählen folgende Systemparameter
> l:=1: xl:=0.2: xr:=0.9: ql:=1: qr:=1.5:
# und geben die Ergebnisse auf dem Bildschirm aus.
> grafik(q(x),x,0,l,"\nLinienkraft q(x)\n");
> grafik(Q(x),x,0,l,"\nQuerkraft Q(x)\n");
> grafik(M(x),x,0,l,"\nBiegemoment M(x)\n"); #Das kann etwas dauern
;
# Wirkt auf den Träger an der Stelle  x = xF  eine Einzelkraft F, dann kann diese singuläre Belastung mittels der Dirac-Funktion dargestellt werden. Die Dirac-Funktion ist keine Funktion im üblichen Sinne der Analysis, vielmehr handelt es sich um eine verallgemeinerte Funktion (Distribution), mit der Maple rechnen kann.
> l:='l':M:='M':Q:='Q':
> q:=unapply(F*Dirac(x-xF),x);
> sys:={odeQ,odeM};
# Lösung für den Träger auf zwei Stützen mit Einzelkraft an der Stelle x = xF:
> loe_03:=dsolve(`union`(sys,RB_03),{Q(x),M(x)}):assign(%):
> Q:=unapply(Q(x),x);
> M:=unapply(M(x),x);
# Wir stellen die Ergebnisse grafisch dar und wählen
> F:=1: l:=1: xF:=l/4:
> grafik(Q(x),x,0,l,"\nQuerkraft Q(x)\n");
> grafik(M(x),x,0,l,"\nBiegemoment M(x)\n");
# Sind wir an der Verformung des Balkens interessiert, also an seiner Durchbiegung w(x) und Verdrehung Õ(x) = w'(x), dann greifen wir im statisch bestimmten Fall auf die Differenzialgleichung der Balkenbiegung w''(x) = -M(x)/B zurück. Darin bezeichnen M(x) das Biegemoment und B = EIyy die konstante Biegesteifigkeit als Produkt aus dem Elastizitätsmodul E und dem Flächenträgheitsmoment Iyy. 
# Wir wollen im Folgenden die Biegelinie und deren Tangente für den links eingespannten und rechts frei gelagerten Balken (Kragträger) unter einer konstanten Linienkraftbelastung q0 berechnen. Dazu beschaffen wir uns zunächst die Schnittlasten Q(x) und M(x) durch Integration der Schnittlastendifferenzialgleichungen: 
> l:='l': M:='M':Q:='Q':q0:='q0':
> q:=unapply(q0,x);
> sys:={odeQ,odeM};
> loe_w:=dsolve(`union`(sys,RB_01),{Q(x),M(x)});assign(%):
> Q:=unapply(Q(x),x);
> M:=unapply(M(x),x);
# Die Differenzialgleichung zur Berechnung der Biegelinie  lautet
> ode2:=(D@@2)(w)(x)=-M(x)/B;
# Für den bei x = 0 eingespannten Kragträger sind folgende Randbedingungen (Einspannbedingungen) einzuhalten:
> RB_04:=w(0)=0,D(w)(0)=0;
> loe_05:=dsolve({ode2,RB_04},w(x)):assign(%):
> w:=unapply(w(x),x);
# Die Neigung der Tangente an die Biegelinie ist
> phi:=unapply(D(w)(x),x);
# Wir geben die Ergebnisse grafisch aus und wählen
> q0:=1: l:=1: B:=1:
> grafik(w(x),x,0,l,"\nDurchbiegung w(x)\n");
> grafik(phi(x),x,0,l,"\nTangentenneigung phi(x)\n");
# Im Folgenden nehmen wir an, dass der Balken einen Rechteckquerschnitt der Breite b besitzt, dessen Höhe  h(x)  linear von h[`&ell;`];  am linken Ende auf h[r]; am rechten Ende anwächst, also h(x) = h[`&ell;`][1+(h[r]/h[`&ell;`]-1)*x/`&ell;`] and h[`&ell;`][1+(h[r]/h[`&ell;`]-1)*x/`&ell;`] = (h*`&ell;`*h)*n; . Im Bauwesen wird in diesem Fall von einem gevouteten Träger gesprochen. Setzen wir zur Abkürzung eta = h[r]/h[`&ell;`];,  xi=(x)/(`&ell;`),  I[yy,0]=(b (h[`&ell;`])^(3))/(12)   und   W[y,0]=(b (h[`&ell;`])^(2))/(6),;dann erhalten wir die Flächenmomente bezüglich einer Achse durch den Flächenmittelpunkt:
#  I[yy]  =I[yy,0] [1+(eta-1)xi]^(3)   =;  I[yy, 0]*I[yyn];(¾)         (Flächenmoment zweiten Grades, Flächenträgheitsmoment),
#  W[y] = W[y, 0]*[1+(eta-1)*xi]^2;  =  W[y, 0]*W[n](xi);       (Widerstandsmoment).
# Die normierten Größen sind
> hn:=unapply((1+(eta-1)*xi),xi);
> Iyyn:=unapply(hn(xi)^3,xi);
> Wn:=unapply(hn(xi)^2,xi,eta);
# Die Schnittlasten des statisch bestimmt gelagerten Balkens unter Linienkraftbelastung  q0  ergeben sich aus reinen Gleichgewichtsbedingungen und sind demzufolge unabhängig von der Biegesteifigkeitsverteilung längs der Balkenachse. Für den weiteren Rechengang normieren wir das Biegemoment 
# M(x) = ((1/2)*q[0]*`&ell;`^2*(x/`&ell;`))*(1-x/`&ell;`) and ((1/2)*q[0]*`&ell;`^2*(x/`&ell;`))*(1-x/`&ell;`) = M[0]*M[n]*mit*M[0] and M[0]*M[n]*mit*M[0] = q[0]*`&ell;`^2, M[n] = (1/2)*xi(1-xi), xi = x/`&ell;`;
# Die Differenzialgleichung lautet 
# w''(¾) = - M(xi)*`&ell;`^2/(E*I[yy](xi)) = -M[0]*`&ell;`^2/(E*I[yy, 0]);M[n](xi)/I[yyn](xi); = -KM[n](xi)/I[yyn](xi);   mit   K = M[0]*`&ell;`^2/(E*I[yy, 0]); = q[0]*`&ell;`^4/(E*I[yy, 0]);
# 
# Wir führen noch die normierte Verschiebung w[n](xi) = w(xi)/K; ein und erhalten dann die normierte Differenzialgleichung: wn''(¾) = -;M[n](xi)/I[yyn](xi);.;
> Mn:=unapply(1/2*xi*(1-xi),xi);
> ode2:=(D@@2)(wn)(xi)=-Mn(xi)/Iyyn(xi);
# Beim Träger auf zwei Stützen müssen an den Balkenenden  die Verschiebungen verschwinden:
> RB_05:=wn(0)=0,wn(1)=0;
> loe_06:=dsolve({ode2,RB_05},wn(xi)): assign(%);
> wn:=unapply(wn(xi),xi,eta);
# Die Neigung der Tangente an die Biegelinie ist
> ws:=D[1](wn);
# Für die wahren Größen folgt unter Berücksichtigung der Konstanten K 
> w:=unapply(K*wn(xi,eta),xi,eta);
> l:='l':
> phi:=unapply(K/l*ws(xi,eta),xi,eta);
# Im Fall · = 1 (Balken  mit konstanter Höhe h) versagen zunächst die obigen Lösungen. Wir ersetzen diese durch ihre  Grenzwerte
> wk  := unapply(factor(limit(w(xi,eta)  ,eta=1)),xi);  wk(1/2);
> phik:= unapply(factor(limit(phi(xi,eta),eta=1)),xi);
# Wir konkretisieren unsere Aufgabe und wählen als Beispiel einen Stahlbetonbalken (Einheiten m und N). 
> l:=10.; hl:=0.30: hr:=0.50: eta:=hr/hl; b:=0.3; 
> q0:=2000.: E:=2.0*10^(10): B0:=E*b*hl^3/12.; K:=q0*l^4/B0;
# Durchbiegung
> grafik(w(xi,eta),xi,0,1,"\nDurchbiegung w(xi) in [m]\n");
> grafik(phi(xi,eta),xi,0,1,"\nTangentenneigung phi(xi)\n");
# Die größte Durchbiegung tritt beim gevouteten beidseitig drehbar gelagerten Balken nicht in Feldmitte auf, sondern bei
> xinull:=fsolve(phi(xi,eta),xi=0.4..0.5);
# Maximale Durchbiegung
> maxw:=evalf(w(xinull,eta));
# Die Biegerandspannung am oberen Rand errechnet sich aus der Beziehung  sigma; =M/W[y];=M[0]/W[y, 0];M[n]/W[n]; = C M[n]/W[n];,  C  = M[0]/W[y0];  = (6 q[0] `&ell;`^(2))/(bh[`&ell;`]2);
> eta:='eta':
> sigma:=unapply(C*Mn(xi)/Wn(xi,eta),xi,eta);
# Die Änderung von  Ã  mit  ¾  ist
> sigmastr:=D[1](sigma);
# Dort, wo  sigmastr verschwindet, ist der Ort der größten Randspannung, der nicht identisch ist mit dem Ort der größten Durchbiegung 
> maxs:=solve(sigmastr(xi,eta),xi);
# Wir konkretisieren unser Problem
> C:=6*q0*l^2/(b*hl^2); eta:=hr/hl;
> maxs;sigma(maxs,eta);
> grafik(10^(-6)*sigma(xi,eta),xi,0,1,"\nBiegerandspannung  sigma(xi)  in  [N/mm^2]\n");
> a:='a':hl:='hl':hr:='hr':l:='l':x:='x':eta:='eta':wn:='wn':
# Wir betrachten einen links eingespannten Kragträger unter einer Einzelkraft F  bei x =  (rechtes Balkenende). Der Balken besitzt die konstante Breite b mit einer am linken Rand (x = 0) beginnenden zur x-Achse symmetrischen Voute, die sich jetzt aber nicht über die gesamte Trägerlänge erstrecken soll, sondern nur noch bis zur Stelle x = a and a = `&ell;`*alpha*mit*alpha and `&ell;`*alpha*mit*alpha = a/`&ell;`;. 
# Damit ändert sich die Balkenhöhe längs der x-Achse in der Form
> hn:=unapply(1+(eta-1)*(xi/alpha+(1-xi/alpha)*Heaviside(xi-alpha)),xi);
> alpha:='alpha':
> Iyyn  :=unapply(hn(xi)^3,xi);
> Wn    :=unapply(hn(xi)^2,xi);
> Mn    :=unapply(-(1-xi),xi);
# Wir wählen folgende Zahlenwerte
> alpha :=1/3;hl:=0.6;hr:=0.3;eta:=hr/hl;
> par:= xi=0..1,thickness=2,titlefont=[HELVETICA,14],color="NavyBlue",thickness=2,filled=[color="Blue",transparency =.5],gridlines=true,axes=boxed:
> p1:=plot([hn(xi),-hn(xi)],xi=0..1,thickness=2,color=[black,black],title="\nBezogene Querschnittshöhe\n",titlefont=[HELVETICA,14]):
> p2:=plot(Iyyn(xi),par,title="\nBezogenes Flächenträgheitsmoment\n"):
> p3:=plot(Wn(xi)  ,par,title="\nBezogenes Widerstandsmoment\n"):
> p4:=plot(Mn(xi)  ,par,title="\nBezogenes Biegemoment\n"):
> display(Matrix(1,4,[p1,p2,p3,p4]));
> ode2  :=(D@@2)(wn)(xi)=-Mn(xi)/Iyyn(xi);
> RB_07 :=wn(0)=0,D(wn)(0)=0;
> loe_07:=dsolve({ode2,RB_07},wn(xi)): assign(%);
> wn:=unapply(wn(xi),xi);
> b:=0.3;l:=10.:F:=10^3;E:=2.0*10^(10);
> B0:=E*b*hl^3/12.; K:=F*l^3/B0;
> grafik(K*wn(xi),xi,0,1,"\nDurchbiegung w(xi) in [m]\n");
# Die größte Durchbiegung tritt am rechten Balkenende bei ¾ = 1 auf. 
> maxw:=evalf(K*wn(1));
# Da  die verallgemeinerten Ableitungen der Heaviside-Funktionen keine Rolle spielen, bleiben im Folgenden alle Dirac-Terme unberücksichtigt. Die Neigung der Tangente an die Biegelinie erhalten wir dann mit 
> wstr:=subs(Dirac=0,D(wn));
> phi:=unapply(K/l*wstr(xi,eta),xi);
> grafik(phi(xi),xi,0,1,"\nTangentenneigung phi(xi)\n");
# Die Randspannung am oberen Rand (unten gilt der negative Wert, Druckspannung) ergibt sich mit der dimensionsbehafteten Konstanten C zu 
> C:=6*F*l/(b*hl^2); sigma:=unapply(-C*Mn(xi)/Wn(xi,eta),xi);
> grafik(10^(-6)*sigma(xi),xi,0,1,"\nBiegerandspannung sigma(xi) in [N/mm^2]\n");
# Die größte Randspannung tritt bei  ¾ = ±  auf. Da die Funktion sigma in diesem Punkt zwar stetig aber nicht stetig differenzierbar ist, bilden wir den links- bzw. rechtseitigen Grenzwert:
> sigmax:=10^(-6)*limit(sigma(xi),xi=alpha,left);
> sigmax:=10^(-6)*limit(sigma(xi),xi=alpha,right);
# Als Beispiel für ein Anfangswertproblem, das auf eine Differenzialgleichung 1. Ordnung führt, soll der Schwingungsversuch an einem Kelvin-Modell untersucht werden. Die theoretischen Untersuchungen können dem Buch : F.U. Mathiak (2010), Strukturdynamik diskreter Systeme, Oldenbourg Wissenschaftsverlag GmbH  entnommen werden. 
> restart: with(plots):
# 
# Ein Kelvin-Modell (s.h. obige Abbildung) besteht aus der Parallelschaltung von Feder und Dämpfer. Ist l0 die Elementlänge zum Zeitpunkt t = t0, dann lautet die zugehörige Differenzialgleichung:
> ode1:=diff(l(t),t)+k/c*l(t)-1/c*(F(t)+k*l0);
# Wir führen noch dl(t) = l(t) - l0 ein, also
> l(t):=dl(t)+l0;
# dann ist
> ode1;
# Die äußere Belastung sei harmonisch (A: Amplitude, ©: Erregerkraftfrequenz)
> F(t):=A*cos(Omega*t);
# Zum Zeitpunkt t = t0 sei das Element entspannt. 
> AB:=dl(t0)=0;
# Das führt auf folgende Lösung für die Längenänderung dl(t)
> loe:=dsolve({ode1,AB},dl(t));
# Ist die Feder zum Zeitpunkt t = t0 = 0 entspannt, dann verbleibt
> t0:=0;dl:=unapply(simplify(rhs(loe)),t);
# Für den weiteren Rechengang wählen wir folgende Systemwerte
> A:=1; Omega:=1; k:=1; c:=3/4;
# Das führt auf die Längenänderung
> dl(t);
# Damit ergeben sich folgende Kräfte in Feder und Dämpfer
> FF:=unapply(k*dl(t),t);   #Federkraft
;
> FD:=unapply(c*D(dl)(t),t); #Dämpferkraft
;
# Wir stellen die Ergebnisse grafisch dar. Zunächst werden die Kräfte im Kelvin-Modell als Funktion der Zeit t geplottet. 
> plot([FF,FD,FF+FD],0..20,gridlines,title ="\nKräfte im Kelvin-Modell\n",titlefont=["ROMAN",15],legend =["FF","FD","FD+FD"]);
# In einem zweiten Schritt erfolgt die Ausgabe der Kräfte als Funktion den Elementlängenänderung dl.
> tE:=2*Pi/Omega;
> plot({[dl(t),FF(t),t=0..tE],[dl(t),FD(t),t=0..tE],[dl(t),FF(t)+FD(t),t=0..tE]},title ="\nKräfte im Kelvin-Modell\n",titlefont=["ROMAN",15],scaling=constrained,axes = boxed,gridlines,color=[black,red,blue]);
# Nach hinreichend großer Zeit verschwindet der Exponentialanteil in den Lösungen und es verbleibt der Partikularanteil
> dlp:=unapply(eval(subs(exp=0,dl(t))),t); #Längenänderung
;
> FFp:=unapply(k*dlp(t),t);                #Kraft in der Feder
;
> FDp:=unapply(c*D(dlp)(t),t);             #Kraft im Dämpfer
;
> plot([FFp,FDp,FFp+FDp],0..20,gridlines,title ="\nKräfte im Kelvin-Modell\n",titlefont=["ROMAN",15],legend =["FFp","FDp","FDp+FDp"]);
> plot({[dlp(t),FFp(t),t=0..tE],[dlp(t),FDp(t),t=0..tE],[dlp(t),FFp(t)+FDp(t),t=0..tE]},scaling=constrained,axes = boxed,gridlines,color=[black,red,blue]);
# Abschließend untersuchen wir noch den Kriechversuch am Kelvin-Element. Dazu belasten wir das Kelvin-Modell zum Zeitpunkt t = t0 mit einer konstante Kraft F = F0, die zum Zeitpunkt t = t1 wieder entfernt wird. 
> restart: with(plots):
# 
# 
> F(t):=piecewise(t < t1 and t > t0,F0,0);
> ode1:=diff(l(t),t)+k/c*l(t)-1/c*(F(t)+k*l0);
> l(t):=dl(t)+l0;
> ode1;
# Wir wählen folgende Systemwerte
> F0:=0.5; k:=0.5; c:=2; t0:=5; t1:=30; tE:= 50;
# Die Retardationszeit ist
> tRet:=c/k;
# Die Anfangsbedingung lautet
> AB:=dl(t0)=0;
# Damit ergibt sich folgende Lösung
> loe:=dsolve({ode1,AB},dl(t)): dl:=unapply(rhs(loe),t);
# Wir geben die Element-Längenänderung als Funtion der Zeit aus:
> plot(dl,0..tE,gridlines,title = "\nKriechversuch am Kelvin-Modell\n",titlefont = ["ARIAL", 15],labels = ["Zeit t", "Längenänderung dl"], labeldirections = ["horizontal", "vertical"],axes=boxed,filled = [color = "Blue", transparency = .5]);
# Beim Kelvin-Modell tritt nach der Entlastung keine bleibende Verformung auf. Das Modell verhält sich unmittelbar nach Lastaufbringung viskos, langfristig dagegen elastisch. Dieser Effekt wird verzögerte Elastizität, elastische Nachwirkung oder auch Anelastizität genannt. Wird die Last entfernt, dann tritt Rückkriechen (Kriecherholung) des Modells auf die Länge l0 auf (dl = 0). Wird die Länge l konstant gehalten, dann resultiert daraus eine Kraft F = k dl.
> 
;
