# Kapitel 3: 
# Balken unter harmonischer Einzelkraftbelastung
# © 2017  Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
# Beispiel 3-9
> 
;
# 
# Der in der oberen Abbildung dargestellte Balken der Breite b und der Höhe h wird durch die harmonische Einzelkraft F(t) = F0 sin ©t belastet. Zum Zeitpunkt t = 0 befindet sich der Balken in Ruhe. Gesucht werden sämtliche Zustandsgrößen. Animieren Sie mit Maple sämtliche Zustandsgrößen. 
# Geg.: b = 20 cm, h = 50 cm,  = 8,0 m, ñ = 2500  kg/m3, E = 3,051010 N/m2, a = 6,0 m, c = 1,0 m, F0 = 4.0 MN, © = 200 s-1.
> restart:with(LinearAlgebra):
# Eigenfunktionen
> Wn:=sin(n*Pi*x/l);
> gl1:={sin(n*Pi)=0};
# Generalisierte Masse
> mgn:=simplify(mu*int(Wn^2,x=0..l),gl1);
# Äußere Belastung Einzelkraft  F0 an der Stelle x = a
> g:=F0*Dirac(x-a); 
> h:=unapply(sin(Omega*t),t);
> gl2:={cos(alpha-beta)-cos(alpha+beta)=2*sin(alpha)*sin(beta)};
> f0n:=simplify(int(g*Wn,x=0..l),gl2)/mgn assuming l > a and a > 0;
# Partikuläre Lösung
> Tn:=f0n/omega[n]*simplify(int(h(tau)*sin(omega[n]*(t-tau)),tau=0..t));
# Verschiebung
> wn:=Wn*Tn;
# Im Resonanzfall (© = Ék ) ist 
> wnEB:=simplify(limit(wn,Omega=omega[n]));
# Geschwindigkeit
> vn:=diff(wn,t);
# Schnittlasten
> Mn:=-E*Iyy*diff(wn,x,x); #Biegemoment
;
> Qn:= diff(Mn,x);         #Querkraft
;
# Werte des Beispiels:
> b:=0.2; h:=0.5; A:=b*h; Iyy:=b*h^3/12; l:=8.;
> E:=3.05E10; rho:=2500.; mu:=rho*A;  K:=sqrt(E*Iyy/rho/A/l^4*Pi^4);
> F0:=4.0E6; Omega:=200.; a:=6.;
> wnS:=0: nmax:=25; eps:=10.^(-Digits+1.);
> for n to nmax do
>   omega[n]:=(n*Pi)^2*K;
>   if abs(omega[n]-Omega) < eps then 
>     wnS:= wnS + wnEB;
>       else
>     wnS:= wnS + wn;
>   end if;
> end do:
> vnS:= diff(wnS,t):
> MnS:=-E*Iyy*diff(wnS,x,x):
> QnS:= diff(MnS,x):
> te:=0.1; #Animationszeit
;
> para:= t=0..te, axes=boxed, frames=100,thickness=3,gridlines=true, labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10], axesfont = ["HELVETICA", "ROMAN", 8]:
# Achtung: Zur besseren Anschauung Verschiebung nach unten positiv gewählt, Vorzeichenwechsel in wnS
> 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]));
> 
;
