# Kapitel 3: 
# Kragbalken unter impulsartiger Einzelkraftbelastung
# © 2017  Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
> 
;
# Beispiel 3-10:
# Der in der unteren Abbildung dargestellte Zweig einer Buche wird an der Stelle x = a zum Zeitpunkt t = t0 durch die Einzelkraft der Intensität G0  belastet. Für Zeitpunkte t < t0 befindet sich der Balken in Ruhe. Unmittelbar nach dem Stoß setzt ein Schwingungsvorgang ein, den es zu berechnen gilt. Der Buchenzweig (ñ: Dichte, E: Elastizitäsmodul) besitzt einen kreisförmigen Querschnitt (r: Radius) und befindet sich vor der Belastung in gestreckter Lage. Animieren Sie sämtliche Zustandsgrößen. 
# Geg.: r = 1,0 cm,  = 1,50 m, ñ = 470 kg/m3, E = 3,05108 N/m2, G0 = 0,55 N.
> 
;
# 
> restart: with(LinearAlgebra):
> alias(H = Heaviside);
# Die Eigenfunktionen des Kragträgers sind
> Wn:=sin(lambda*x/l)-sinh(lambda*x/l)+(cosh(lambda*x/l)-cos(lambda*x/l))*Wn1;
> f0n:=G0*int(Dirac(x-a)*Wn,x=0..l)/mgn assuming a > 0, l > a;
# Für den zeitabhängigen Teil der Lösung erhalten wir:
> TDn:=f0n/omega*int((H(tau-t0)-H(tau-t1))*sin(omega*(t-tau)),tau=0..t) assuming t0 > 0,t1 > 0;
> limit(TDn,t1=t0,right);
# Der Faktor bei der Eigenfunktion ist:
> Wn1:=(sinh(lambda)+sin(lambda))/(cosh(lambda)+cos(lambda)):
# Generalisierte Masse
> mgn:=mu*int(Wn^2,x=0..l);
# Ein Reihenglied der Bewegungsgleichung ist dann
> wn:=simplify(Wn*TDn,trig);
# Die Eigenwerte für den Kragträger folgen aus der Eigenwertgleichung
> EWG:=unapply(cos(lambda)+1/cosh(lambda),lambda);
> plot(EWG(lambda),lambda=0..30,gridlines);
# Die obige Grafik zeigt uns im Bereich 0 d » d 30 genau 10 Nullstellen.
> Z0:=[]; a:=0.; b:=30.: dz:=0.1; dmax:=round((b-a)/dz); za:= a; ze:=za+dz;
> for n to dmax do
>   if evalf(EWG(za)*EWG(ze)) < 0. then
>     Z0:=[op(Z0),fsolve(EWG(z),z=za..ze)]:
>   end if:
>   za:= ze; ze:=za+dz;
> end do: 
# Mit den Werten des Beispiels:
> r:=0.005;  A:=r^2*Pi; Iyy:=Pi/4*r^4; l:=1.;
> rho:=470.; E:=3.05E8;  mu:=rho*A;   By:=E*Iyy;  K:=sqrt(E*Iyy/mu/l^4);
> G0:=.40;    a:=0.3*l;   t0:= 0.1;    t1:= 0.6;  te:= 1.0;
# Berechnung der Eigenkreisfrequenzen und der Eigenfrequenzen
> nmax:=nops(Z0): 
> ZN:=convert(Z0,Vector): OMEGA:=Vector(nmax): FR:=Vector(nmax):
> for n to nmax do
>   OMEGA[n]:=evalf(Z0[n]^2*K);
>   FR[n]   :=OMEGA[n]/2/Pi;
> end do:
# Wir geben die Eigenwerte » , Eigenkreisfrequenzen  É  und Frequenzen f aus:
> print(`Lambda = `,ZN,`  Omega = `,OMEGA, `  F = `, FR);
> wnS:=0:
> for n to nmax do
>   lambda:= ZN[n];
>    omega:= OMEGA[n];
>      wnS:= wnS + evalf(wn);
> end do:
> wnS1:=diff(wnS,t$1):
# Geschwindigkeit
> vnS:=diff(wnS,t):
# Schnittlasten
> MnS:=-By*diff(wnS,x,x): #Biegemoment
;
> QnS:=      diff(MnS,x): #Querkraft
;
> para:= t=0..te, axes=boxed, frames=101,thickness=3,gridlines=true, labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10], axesfont = ["HELVETICA", "ROMAN", 8]:
> p1:=plots[animate](plot,[-wnS, x=0..l], para, labels = ["x", "Verschiebung in [m]"]):
> p2:=plots[animate](plot,[ vnS, x=0..l], para, labels = ["x", "Geschwindigkeit in [m/s]"]):
> p3:=plots[animate](plot,[ MnS, x=0..l], para, labels = ["x", "Biegemoment in [Nm]"]):
> p4:=plots[animate](plot,[ QnS, x=0..l], para, labels = ["x", "Querkraft in [N]"]):
> plots[display](Matrix(1,4,[p1,p2,p3,p4]));
> 
;
