# Kapitel 3: 
# Spire of Dublin nach Bernoulli 
# ohne Berücksichtigung der Drehträgheit
# © 2017  Friedrich U. Mathiak, 
# mathiak@mechanik-info.de
# 
> restart: with(LinearAlgebra): with(plots): with(DEtools):
# Beispiel 3-6:
# 
# 
# Der Spire of Dublin, ein bekanntes Denkmal in der irischen Hauptstadt, besteht aus einer 120m hohen Stahlnadel. Der Radius ist linear mit der Höhe veränderlich. Für den vorliegenden dünnen Kreisquerschnitt sind die Eigenfrequenzen und Eigenformen für den Bernoulli-Balken ohne Berücksichtigung der Drehträgheit zu berechnen. Der Einfluss der Normalktraft aus Eigengewicht kann vernachlässigt werden. 
# Differenzialgleichung der zeitfreien Auslenkung W(x):
> ode4:=diff(z^3*diff(W(z),z$2),z$2) - kappa^4*z*W(z);
> odeadvisor(ode4);
> loe:=dsolve(ode4);
> FunctionAdvisor(describe, BesselI);FunctionAdvisor(describe, BesselJ);FunctionAdvisor(describe, BesselK);FunctionAdvisor(describe, BesselY);
# Wir testen die Lösung
> odetest(loe, ode4);
> W:=rhs(loe);
# Die Verschiebung am Stabende x = h ist
> Wh:=simplify(subs(z=1+delta,W));
# Ableitungen der Verschiebungsfunktion W(z) nach z:
> W1:=simplify(diff(W,z$1));
> W2:=simplify(diff(W,z$2));
> W3:=simplify(diff(W,z$3));
# Schnittlasten:
> MS:=simplify(-delta^2/h*z^3*W2);
> QS:=simplify(-delta^3/h*(3*z^2*W2+z^3*W3));
# Lösung des Randwertproblems:
> gl1:=subs(z=1,W);
> gl2:=subs(z=1,W1);
> gl3:=simplify(subs(z=1+delta,MS));
> gl4:=simplify(subs(z=1+delta,QS));
> gl:=[gl1,gl2,gl3,gl4]: unb:=[_C1,_C2,_C3,_C4]:
> A, b := GenerateMatrix(gl, unb):
> DA:=simplify(Determinant(A),symbolic);
# Wir entfernen die unbedeutenden Faktoren aus der Eigenwertgleichung:
> DA:=kappa*(-kappa*(1+delta)*(BesselK(1, 2*kappa)*BesselY(0, 2*kappa)+BesselY(1, 2*kappa)*BesselK(0, 2*kappa))*BesselJ(0, 2*kappa*sqrt(1+delta))+kappa*(1+delta)*(BesselJ(1, 2*kappa)*BesselK(0, 2*kappa)+BesselK(1, 2*kappa)*BesselJ(0, 2*kappa))*BesselY(0, 2*kappa*sqrt(1+delta))-(1/4)*sqrt(1+delta)*((-4+(1+delta)*kappa^2)*(BesselK(1, 2*kappa)*BesselY(0, 2*kappa)+BesselY(1, 2*kappa)*BesselK(0, 2*kappa))*BesselJ(1, 2*kappa*sqrt(1+delta))+(1+delta)*kappa^2*(BesselJ(1, 2*kappa)*BesselY(0, 2*kappa)-BesselY(1, 2*kappa)*BesselJ(0, 2*kappa))*BesselK(1, 2*kappa*sqrt(1+delta))-BesselY(1, 2*kappa*sqrt(1+delta))*(-4+(1+delta)*kappa^2)*(BesselJ(1, 2*kappa)*BesselK(0, 2*kappa)+BesselK(1, 2*kappa)*BesselJ(0, 2*kappa))))*BesselI(0, 2*kappa*sqrt(1+delta))+(1/4)*kappa*(4*kappa*(1+delta)*(BesselY(1, 2*kappa)*BesselI(0, 2*kappa)-BesselY(0, 2*kappa)*BesselI(1, 2*kappa))*BesselK(0, 2*kappa*sqrt(1+delta))+sqrt(1+delta)*((BesselK(1, 2*kappa)*BesselY(0, 2*kappa)+BesselY(1, 2*kappa)*BesselK(0, 2*kappa))*(4+(1+delta)*kappa^2)*BesselI(1, 2*kappa*sqrt(1+delta))+(BesselY(1, 2*kappa)*BesselI(0, 2*kappa)-BesselY(0, 2*kappa)*BesselI(1, 2*kappa))*(4+(1+delta)*kappa^2)*BesselK(1, 2*kappa*sqrt(1+delta))+BesselY(1, 2*kappa*sqrt(1+delta))*kappa^2*(1+delta)*(BesselK(1, 2*kappa)*BesselI(0, 2*kappa)+BesselK(0, 2*kappa)*BesselI(1, 2*kappa))))*BesselJ(0, 2*kappa*sqrt(1+delta))-kappa*(kappa*(1+delta)*(BesselJ(1, 2*kappa)*BesselI(0, 2*kappa)-BesselJ(0, 2*kappa)*BesselI(1, 2*kappa))*BesselY(0, 2*kappa*sqrt(1+delta))+(1/4)*sqrt(1+delta)*((1+delta)*kappa^2*(BesselJ(1, 2*kappa)*BesselY(0, 2*kappa)-BesselY(1, 2*kappa)*BesselJ(0, 2*kappa))*BesselI(1, 2*kappa*sqrt(1+delta))+((-BesselY(1, 2*kappa)*BesselI(0, 2*kappa)+BesselY(0, 2*kappa)*BesselI(1, 2*kappa))*BesselJ(1, 2*kappa*sqrt(1+delta))+BesselY(1, 2*kappa*sqrt(1+delta))*(BesselJ(1, 2*kappa)*BesselI(0, 2*kappa)-BesselJ(0, 2*kappa)*BesselI(1, 2*kappa)))*(-4+(1+delta)*kappa^2)))*BesselK(0, 2*kappa*sqrt(1+delta))-(1/4)*sqrt(1+delta)*kappa*((BesselJ(1, 2*kappa)*BesselK(0, 2*kappa)+BesselK(1, 2*kappa)*BesselJ(0, 2*kappa))*(4+(1+delta)*kappa^2)*BesselI(1, 2*kappa*sqrt(1+delta))+(1+delta)*kappa^2*(BesselK(1, 2*kappa)*BesselI(0, 2*kappa)+BesselK(0, 2*kappa)*BesselI(1, 2*kappa))*BesselJ(1, 2*kappa*sqrt(1+delta))+(BesselJ(1, 2*kappa)*BesselI(0, 2*kappa)-BesselJ(0, 2*kappa)*BesselI(1, 2*kappa))*(4+(1+delta)*kappa^2)*BesselK(1, 2*kappa*sqrt(1+delta)))*BesselY(0, 2*kappa*sqrt(1+delta))+((-BesselK(1, 2*kappa)*BesselY(0, 2*kappa)-BesselY(1, 2*kappa)*BesselK(0, 2*kappa))*BesselJ(1, 2*kappa*sqrt(1+delta))+BesselY(1, 2*kappa*sqrt(1+delta))*(BesselJ(1, 2*kappa)*BesselK(0, 2*kappa)+BesselK(1, 2*kappa)*BesselJ(0, 2*kappa)))*BesselI(1, 2*kappa*sqrt(1+delta))+((-BesselY(1, 2*kappa)*BesselI(0, 2*kappa)+BesselY(0, 2*kappa)*BesselI(1, 2*kappa))*BesselJ(1, 2*kappa*sqrt(1+delta))+BesselY(1, 2*kappa*sqrt(1+delta))*(BesselJ(1, 2*kappa)*BesselI(0, 2*kappa)-BesselJ(0, 2*kappa)*BesselI(1, 2*kappa)))*BesselK(1, 2*kappa*sqrt(1+delta));
> 
;
# Wir beschafffen uns eine Abschätzung für den kleinsten Eigenwert
> simplify(series(DA,kappa,8)); convert(%,polynom); k1:=simplify(solve(%,kappa),symbolic);
# Berechnung der Integrationskonstanten. Da eine Konstante frei gewählt werden kann, setzen wir 
> _C4:=1:
# Das verbleibende Gleichungssystem ist dann
> AR, b := GenerateMatrix([gl1,gl2,gl3],[_C1,_C2,_C3]):
> unb:=LinearSolve(AR,b):
> _C1:=simplify(unb[1]): _C2:=simplify(unb[2]): _C3:=simplify(unb[3]):
# Mit den Werten des Beispiels erhalten wir:
> ro:=0.075: ru:=1.5; delta:=ro/ru-1; t:=0.025; h:=120.;
> A0:=2*Pi*ru*t;I0:=Pi*ru^3*t;iy0:=sqrt(I0/A0);
> E:=2.1E11; rho:=7850.;
> EWG:=unapply(simplify(DA),kappa);
> plot(EWG(kappa),kappa=0..10,-100..100,gridlines,axes=boxed);
# Berechnung der Eigenwerte º:
# In der folgendenSchleife werden dieNullstellen der  Eigenfunktion EWG(z) berechnet. Dazu wird EWG(z)  im Intervall [a,b] mit der Schrittweite dz ausgewertet. Ist EWG(z)*EWG(z+dz)  < 0, dann hat die Eigenfunktion  in diesem Intervall einen Nulldurchgang. In diesem Fall wird die Nullstelle mit  fsolve  nummerisch konkretisiert. 
> Z0:=[]: a:=2.: b:=10: 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: 
> nmax:=nops(Z0): f:=sqrt(E*I0/rho/A0): 
> ZN:=convert(Z0,Vector): OMEGA:=Vector(nmax): FR:=Vector(nmax):
> for n to nmax do
>   OMEGA[n]:=evalf((delta*Z0[n]/h)^2*f);
>   FR[n]   :=OMEGA[n]/2/Pi;
> end do:
> print(`º = `,ZN,`É = `,OMEGA, `f = `, FR);
# Die Abschätzungsformel  für den kleinsten Eigenwert k1 ergibt:
> evalf(k1);
# Wir stellen die Zustandsgrößen zum 1. Eigenwert grafisch dar. 
> kappa:=ZN[1]: z:=1. + delta*x/h:
# Die Konstanten sind:
> _C1;_C2;_C3;_C4;
> Wstr:=delta/h*evalf(W1);
> para:= x=0..h,gridlines=true,titlefont = ["ARIAL", 15],  labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10],axesfont = ["HELVETICA", "ROMAN", 8],axes=boxed:
> p1:= plot(QS  ,para,title = "\nNormierte Querkraft QS\n"    ,labels = ["x", "QS (x)"]):
> p2:= plot(MS  ,para,title = "\nNormiertes Biegemoment MS\n" ,labels = ["x", "MS (x)"]):
> p3:= plot(W,para,title = "\nVerschiebung W\n",labels = ["x", "W (x)"]):
> p4:= plot(Wstr,para,title = "\nTangentenneigung W'\n"  ,labels = ["x", "W'(x)"]):
> plots[display](Matrix(1,4,[p3,p4,p2,p1]));
# In einer weiteren Grafik stellen wir  die ersten drei Verschiebungs-Eigenfunktionen grafisch dar. Dazu werden diese derart normiert, dass jeweils die betragsmäßig größte Verschiebung auf Eins gesetzt wird. 
> WX:=[]:
> for i to 3 do 
>   kappa:=ZN[i];
>   WX   :=[op(WX),W(z)/Wh];
> end do:
> plot(WX,x=0..h,gridlines=true,axes=boxed,title = "\nVerschiebungs-Eigenfunktionen\n", legend = ["Grundschwingung", "1. Oberschwingung","2. Oberschwingung"], titlefont = ["ARIAL", 15],  linestyle = [solid, dot, dash], axesfont = ["HELVETICA", "ROMAN", 8], legendstyle = [font = ["HELVETICA", 9], location =bottom],color=[black,black,black]);
> 
;
