Gelöste Aufgaben/FEAD

Aufgabenstellung

Statt Formfunktionen über die ganze Stablänge anzusetzten wie in FEAB gehen wir jetzt nach der Methode der Finiten Elemente vor.

Lageplan

Gesucht ist die Näherungslösung für die statische Auslenkung der Stab-Querschnitte mit dem Prinzip der Virtuellen Verrückungen und dem Ansatz nach der Methode der Finiten Elemente. Dazu wir unterteilen die Struktur in Elemente und setzen Lineare Trialfunctions für die Verformung in den Elemeten an.

Lösung mit Matlab

Header

Die Lösung basiert auf Maxima 16.04.2 - wir interessieren uns also vor allem für die Struktur der Lösung.

📑 toggle code

/*******************************************************/
/* MAXIMA script                                       */
/* version: wxMaxima 16.04.2                           */
/* author: Andreas Baumgart                            */
/* last updated: 2017-10-10                            */
/* ref: FENV step 4 im Prozess: FE-Ansätze             */
/* description: mit dem PvV werden die Bewegungsgl.    */
/*          für einen Stab unter Gewichtskraft erstellt*/
/*******************************************************/

Declarations

Für Maxima brauchen wir einige Deklarationen.

Die Anzahl der Finiten-Elementen I wählen wir zu

$I=3$

Wir wählen Elemente gleicher Länge, also ist die Elementlänge l_i im Element i

$\begin{array}{ll}{\ell }_i& ={\ell }_e\\ & ={\ell }_0/I\end{array}$.

📑 toggle code

/*******************************************************/
/* declare variational variables - see 6.3 Identifiers */
declare("δW", alphabetic);
declare("δA", alphabetic);
declare("δΠ", alphabetic);
declare("δQ", alphabetic);
declare("δu", alphabetic);

/*******************************************************/
/* parameter: Number of Elements */
I : 3; /* change as required */

Formfunctions

Wir können die Trial-Functions als

$u(x)=\sum _{i=0}^I{U}_i\cdot {\varphi }_i(x)$

mit

${\varphi }_i(x)=\{\begin{array}{cl}{\xi }_{i-1}& \text{ wobei }x={\ell }_e\cdot ((i-1)+{\xi }_{i-1}),0<{\xi }_{i-1}<1\\ (1-{\xi }_i)& \text{ wobei }x={\ell }_e\cdot (i+{\xi }_i),0<{\xi }_i<1\end{array}$.

Für i=2 sieht die Trial-Function dann so aus:

Trial-Function für i=2.

Für die Methode der Finiten Elemente ist allerdings die Sichtweise je Elemement anschaulicher, also

$u_i(x_i)=\{\begin{array}{c}{U}_{i-1}\cdot (1-\left(\frac{x_i}{{\ell }_i}\right))+U_i\cdot \left(\frac{x_i}{{\ell }_i}\right)\\ 0\text{ sonst }\end{array}$,

Und für Element i=2 wird der Verlauf der Querschnitts-Verschiebung also durch

$u_2(x_2)=\{\begin{array}{c}{U}_1\cdot (1-\left(\frac{x_2}{{\ell }_2}\right))+U_2\cdot \left(\frac{x_2}{{\ell }_2}\right)\\ 0\text{ sonst }\end{array}$.

entsprechend

$\delta u_2(x_2)=\{\begin{array}{c}\delta U_1\cdot (1-\left(\frac{x_2}{{\ell }_2}\right))+\delta U_2\cdot \left(\frac{x_2}{{\ell }_2}\right)\\ 0\text{ sonst }\end{array}$

beschreiben.

Bei 3 Finiten Elementen sind die System-Koordinaten dann

$\underset{\_}{Q}=\left(\begin{array}{l}{U}_0\\ U_1\\ U_2\\ U_3\end{array}\right)$

📑 toggle code

/* coordinates and their variations */
 Q : makelist( U[i],i,0,I);
δQ : makelist(δU[i],i,0,I);

/* trial functions */
Phi : [1-(x/l[i]),(x/l[i])];

/* Ansatz for element i */
 u[i]: sum( U[i+j-2]*Phi[j],j,1,2);
δu[i]: sum(δU[i+j-2]*Phi[j],j,1,2);

Element-wise contributions to δW

Die gesamte virtuelle Arbeit am System können wir nun als Summe der virtuellen Einzelarbeiten je Element hinschreiben, also

$\begin{array}{ll}\delta W& =\delta W^a-\delta \mathrm{\Pi }\\ & =\sum _{i=1}^I(\delta W_i^a-\delta {\mathrm{\Pi }}_i)\end{array}$.

Für jedes Element erhalten wir die virtuelle Formänderungsenergie

$\begin{array}{ll}\delta {\mathrm{\Pi }}_i& =\underset{0}{\overset{{\ell }_i}{\int }}EA\cdot {u_i}^{\prime }\cdot \delta {u_i}^{\prime }d{x}_i\\ {\displaystyle }& =(\delta U_{i-1},\delta U_i)\cdot \frac{EA}{{\ell }_i}\cdot \left(\begin{array}{rr}1& -1\\ -1& 1\end{array}\right)\cdot \left(\begin{array}{l}{U}_{i-1}\\ U_i\end{array}\right)\end{array}$

und die virtuelle Arbeit der äußeren, eingeprägten Gewichtskraft

$\begin{array}{ll}\delta W_i^a& =\underset{0}{\overset{{\ell }_i}{\int }}\varrho \cdot A\cdot g\cdot \delta u_{i}d{x}_i\\ {\displaystyle }& =\frac{1}{2}\cdot (\delta U_{i-1},\delta U_i)\cdot \varrho \cdot A\cdot g\cdot {\ell }_i\cdot \left(\begin{array}{l}1\\ 1\end{array}\right)\end{array}$.

📑 toggle code

/* Equilibrium */
δΠ[i] : E*A*integrate(diff(u[i],x)*diff(δu[i],x),x,0,l[i]);
δA[i] : integrate(rho*g*A*δu[i], x,0,l[i]);

K[i] : funmake('matrix, makelist(makelist(coeff(coeff(expand( δΠ[i]), δU[i+j]),U[i+k]), j,-1,0), k,-1,0));
P[i] : funmake('matrix,          makelist(     [coeff(expand( δA[i]), δU[i+j])]       , j,-1,0)         );

Equlilibrium Conditions

Für jedes Element müssen wir nun die virtuellen Arbeiten zum Gesamt-Gleichungssystem zusammenaddieren.

Die Gleichgewichtsbedingung

$\begin{array}{ll}\delta W& =\delta W^a-\delta \mathrm{\Pi }\\ & \stackrel{!}{=}0\end{array}$

liefern ein Gleichungssystem, in das wir nun die kinematische Randbedingung U₀=0, δU₀=0 durch Streichen der ersten Zeile des Gleichungssystems und der ersten Spalte der Gesamt-Steifigkeitsmatrix einarbeiten.

Wir wählen nun noch für jedes Element die gleiche Element-Länge ℓ_i = ℓ_e, mit ℓ_e = ℓ₀/3 und es bleibt

$\frac{EA}{{\ell }_i}\left(\begin{array}{ccc}2& -1& 0\\ -1& 2& -1\\ 0& -1& 1\end{array}\right)\cdot \left(\begin{array}{l}{U}_1\\ U_2\\ U_3\end{array}\right)=\varrho A{\ell }_{i}g\cdot \left(\begin{array}{l}1\\ 1\\ \frac{1}{2}\end{array}\right)$.

📑 toggle code

/* Assembling the System Equations */

K[0] : zeromatrix(length(Q),length(Q))$
P[0] : zeromatrix(length(Q),    1    )$

for j:1 thru I do
   (for row:1 thru 2 do
       (P[0][j-1+row][  1    ]: P[0][j-1+row][  1    ]+P[i][row][ 1 ],
        for col:1 thru 2 do
           K[0][j-1+row][j-1+col]: K[0][j-1+row][j-1+col]+K[i][row][col]));

/* Incorporate Boundary Condition U[0]=0 */
K[0] : submatrix(1,K[0],1);
P[0] : submatrix(1,P[0]);

Solving

Die Gleichgewichtsbedingungen lösen wir nun und erhalten

$U_0=0,U_1=\frac{5g{\ell }_i^2\varrho }{2E},U_2=\frac{4g{\ell }_i^2\varrho }{E},U_3=\frac{9g{\ell }_i^2\varrho }{2E}$

📑 toggle code

/* solution */
sol[0]: linsolve_by_lu(K[0], P[0]);
sol[1]: append([U[0]=0],makelist(U[i]=sol[0][1][i][1],i,1,I))

Post-processing

Wir tragen die Ergebnisse für die numerische Näherungslösung gegen die exakte Lösung auf. Dabei setzen wir l_i=l₀/3. Die Verschiebungen sind elementweise für die Elemente e = 1,2,3 aufgetragen:

Verschiebung u(x) - analytische und FE-Lösung

📑 toggle code

/* plot results*/
/* define piecewise, normalized functions u[0]     */
u[0]: subst(sol[1],makelist(subst([i=j],subst([x=xi*l[i]],u[i])),j,1,I))$
/* normalize */
u[0]: expand(subst([l[i]=l[0]/I],u[0])/(rho*g*l[0]^2/E));

/* analytic solution */
analytic : subst([xi=xi/I],(xi)*(1-(xi)/2))$
toPlot: append(        [[parametric,     t, analytic, [t, 0, I]]], 
               makelist([parametric, i-1+t,  u[0][i], [t, 0, 1]],i,1,I))$
/* use t as independent variable for parametric plot */
toPlot: subst([xi=t], toPlot)$

/* plot */
legende : append([legend, "analytic"],makelist(simplode(["e = ",i]),i,1,I))$
plot2d(toPlot, legende,
               [title, sconcat("number of elements I = ",I)],
               [xlabel, "x/l[e] ->"], [ylabel, "u(x)/(rho*g*l/E)"],
               [gnuplot_preamble,"set key bottom"])$

✔ Spannungen im Stab:

Tragen Sie auch die Spannungen im Stab über die Stablänge an! Berechnen Sie die Spannungen auf Basis der Dehnung  $\epsilon =\frac{du}{dx}.$

Links

• Gelöste Aufgaben/FEAB
• Gelöste Aufgaben/FEAC

Literature

• Strang 2008