Gelöste Aufgaben/FEB1

Aufgabenstellung

Auch wenn es nicht so aussieht: für das rotierende Rotorblatt suchen wir eine statische Lösung - das Problem heißt "quasistatisch".

Lageplan

Ein Hubschrauber-Rotor dreht mir der konstanten Winknelgeschwindigkeit Ω. Das Rotor-Blatt ist aus Aluminium. Gesucht ist die FEM-Lösung für der Verschiebung der Querschnitte und die Dehnung der Querschnitte.

Lösung mit Maxima

Declarations

Für die Komposition der Bewegungsgleichungen brauchen wir die Element-Steifigkeitsmatrix und die Element-Lastmatrix.

Wir gehen von gleichlangen Elementen aus, hier von

$I=4$

Elementen, also

${\ell }_i=\frac{{\ell }_0}{I}$.

Wie im Abschnitt "Finite Elemente Methode" verwenden wir die linearen Ansatzfunktionen

${\varphi }_1=1-\xi \text{ und }{\varphi }_1=\xi$,

also

$u_i({\xi }_i)=U_{i-1}\cdot {\varphi }_1({\xi }_1)+U_i\cdot {\varphi }_2({\xi }_i)$.

Die Element-Steifigkeitsmatrix ergibt sich für den Dehnstab zu

${\underset{\_}{\underset{\_}{K}}}_i=\frac{EA}{{\ell }_i}\cdot \left(\begin{array}{rr}1& -1\\ -1& 1\end{array}\right)$.

Die Element-Lastmatrix ist etwas schwieriger. Für ein einzelnes Element i ist

$\delta W^a=\underset{r_i-1}{\overset{r_i}{\int }}\varrho Ar{\mathrm{\Omega }}^2\cdot \delta u(r)dr$.

Mit

$r=r_{i-1}+{\ell }_i{\xi }_i\text{ und }{r}_{i-1}=(i-1)\cdot {\ell }_i$

schreiben die Beziehung um zu

$\delta W^a=\varrho A{\mathrm{\Omega }}^2{\ell }_i^2\cdot \underset{0}{\overset{1}{\int }}(i+\xi )\cdot \delta u(\xi )d\xi$,

wir erhalten

${\underset{\_}{P}}_i=\varrho A{\mathrm{\Omega }}^2{\ell }_i^2\left(\begin{array}{c}\frac{3i+1}{6}\\ \frac{3i+2}{6}\end{array}\right)$.

📑 toggle code

/*******************************************************/
/* start: FEM solution */
/* Trial-Fucntions */
phi : [(1-xi), xi];

/* declare System Matrices */
K[i] : funmake('matrix, E*A*l[i]*
          makelist(
             makelist(
	        integrate(diff(phi[j],xi)/l[i] *
                          diff(phi[k],xi)/l[i], xi,0,1),
		          j,1,2),k,1,2));
P[i] : funmake('matrix, rho*A*Omega^2*l[i]^2*
             makelist(
	        integrate([(n+xi)*phi[j]], xi,0,1),
		          j,1,2));

/* number of elements */
I : 4;

Equlibrium Conditions

Die Gleichgewichtsbeziehungen des Gesamtsystems erhalten wir durch aus der Addition aller virtuellen Arbeiten des Systems, praktisch durch das Hinzuaddieren der Anteile je Element in die System-Matrizen für K und P:

$\underset{\_}{\underset{\_}{K}}=\left(\begin{array}{ccccc}\frac{AE}{{\ell }_i}& -\frac{AE}{{\ell }_i}& 0& 0& 0\\ -\frac{AE}{{\ell }_i}& \frac{2AE}{{\ell }_i}& -\frac{AE}{{\ell }_i}& 0& 0\\ 0& -\frac{AE}{{\ell }_i}& \frac{2AE}{{\ell }_i}& -\frac{AE}{{\ell }_i}& 0\\ 0& 0& -\frac{AE}{{\ell }_i}& \frac{2AE}{{\ell }_i}& -\frac{AE}{{\ell }_i}\\ 0& 0& 0& -\frac{AE}{{\ell }_i}& \frac{AE}{{\ell }_i}\end{array}\right)$

und

$\underset{\_}{P}=\varrho A{\ell }_i^2{\mathrm{\Omega }}^2\cdot \left(\begin{array}{c}\frac{1}{6}\\ 1\\ 2\\ 3\\ \frac{11}{6}\end{array}\right)$

📑 toggle code

/* initiate system matrices */
K[0] : zeromatrix(I+1,I+1);
P[0] : zeromatrix(I+1, 1 );

/* compose system matrices */
for e : 1 thru I do
   for row : 1 thru 2 do
     (P[0][e-1+row][   1   ] : P[0][e-1+row][   1   ]+subst([n=e-1],P[i][row][ 1 ]),
      for col : 1 thru 2 do
         K[0][e-1+row][e-1+col] : K[0][e-1+row][e-1+col]+              K[i][row][col]);

Boundary Conditions

Die geometrische Randbedingung U₀ = 0 arbeiten wir ein, indem wir die erste Zeile des Gleichungssystem streichen sowie die erste Spalte der Steifigketsmatrix.

📑 toggle code

/* incorporate geometric boundary conditions */
K[0] : submatrix( 1, K[0], 1);
P[0] : submatrix( 1, P[0]   );

Solving

Die Lösung des linearen Gleichungssystems

$\underset{\_}{\underset{\_}{K}}\cdot \underset{\_}{U}=\underset{\_}{P}$

ist

$\underset{\_}{U}=\left(\begin{array}{c}\frac{47l_0^3{\mathrm{\Omega }}^2\rho }{384E}\\ \frac{11l_0^3{\mathrm{\Omega }}^2\rho }{48E}\\ \frac{39l_0^3{\mathrm{\Omega }}^2\rho }{128E}\\ \frac{l_0^3{\mathrm{\Omega }}^2\rho }{3E}\end{array}\right)$.

Oder - in dimensionsloser Form

$\frac{\underset{\_}{U}}{u_s}=\left(\begin{array}{c}\frac{47}{128}\\ \frac{11}{16}\\ \frac{117}{128}\\ 1\end{array}\right)$.

📑 toggle code

/* solve */
print(K[0],"*Q=",P[0])$
sol[2] : subst([l[i]=l[0]/I],linsolve_by_lu(K[0],P[0])[1]);

Post-Processing

Auslenkung u(x).

Die Ergebnisse der FE-Rechnung und der analytischen Lösung können wir jetzt übereinander auftragen:

Die Dehnungen (und Spannungen) im Bauteil sind

${\epsilon }_{rr}=\frac{du}{dr}$.

Auch dieses Ergebnis können wir auftragen:

Dehnung ε(x).

✔ Konstante Dehnung je Element:

Was ausschaut wie ein Fehler - nämlich die "Treppenfunktion" für die Dehnung im FE-Modell - ist in Wirklichkeit die Folge unserer linearen Ansatzfunktionen. Diese abgeleitet liefern eine konstante Dehnung je Element.

📑 toggle code

/* post-processing */
U : append([0],args(transpose(sol[2])[1]/subst(params,u[s])));
X : makelist(i,i,0,I)/I;

plot2d([ratsimp(subst([r=xi*l[0]],subst(disp,u(r))/subst(params,u[s]))),
       [discrete, X, U],[discrete, X, U]], [xi,0,1],
       [legend, "analytic", "FEM", ""],
       [style, [lines,2,1], [lines,2,2], [points,2,2]],
       [point_type, asterisk],
       [xlabel, "ξ →"], [ylabel, "u →"])$

eps: subst(params,u[s])*makelist((U[i+1]-U[i]),i,1,I)/(l[0]/I);
toPlot : append([expand(subst([r=xi*l[0]],diff(subst(disp,u(r)),r))/((l[0]^2*Omega^2*rho)/E))],
                 makelist([discrete,[i-1,i]/I, [eps[i],eps[i]]/((l[0]^2*Omega^2*rho)/E)],i,1,I));
legends : append([legend, "analytic"],makelist(simplode(["Elem. ",i]),i,1,I));

plot2d(toPlot, [xi,0,1], legends, [y,0,0.6],
       [style, [lines,2]],
       [xlabel, "ξ →"], [ylabel, "ε/((l^2 ω^2 ρ/E) →"])$

Links

• ...

Literature

• ...