Randwertprobleme/Methoden zur Lösung von Randwertproblemen/Finite Elemente Methode/FEM: Trial-Functions für lineare Ansatz-Polynome

Auch wenn es so aussieht, als wäre alles zur FEM gesagt, sind die Trial-Funktionen der FEM eine besondere Delikatesse!

Wie bei allen wirklichen Delikatessen steht vor dem Genuss noch etwas Arbeit, denn:

Es gibt einen Trick, wie diese Trial-Funktionen aussehen, damit wir Sie am Rechner einfach zusammenbauen können.

Dafür ist dieses Beispiel (aus FEAB):

Stabmodell im Erdschwerefeld.

Ein dehnbarer Stab ist an der "Decke" im Erd-Schwerefeld aufgehängt. Der Stab hat die Dehnsteifigkeit EA, die Dichte ρ und die Gesamtlänge ''ℓ''.

Gesucht ist eine Approximation der Auslenkung ''u(x)''.

Arbeitsschritte

Die Einführung der Trial-Funktionen der Finiten Elemente ist hier schematisch aufgezeichnet:

STEP A	B	C	D	E
Einführung der Trial-Functions

Dies sind die Arbeitsschritte zur Umsetzung:

• A→B Diskretisierung: Wir teilen die Struktur in N=6 Stücke ( = Finite Elemente ) gleicher Länge ℓ_E=ℓ/N. An jedem Ende eines Finiten Elements entsteht ein Knoten, die nummerieren wir durch von 0 ...6.
• B→C Trial-Funktion je Element:Je Element wählen wir eine Geradenfunktion (ein Polynom 1sten Grades) als Trial-Funktion

  ${\overline{u}}_i(x)=\{\begin{array}{l}{a}_{i,1}\cdot x+a_{i,0}\text{ im Element }i\\ 0\text{ sonst}\end{array}$

  so dass wir als Approximation nun

  $\tilde{u}(x)=\sum _{2N}{A}_n\cdot {\varphi }_n(x)$

  haben, mit den {a_i,1, a_i,0} als gesuchte Wichtungsfaktoren. Hinzu kommen die Übergangsbedingungen: an N-1 Knoten müssen wir für die Geradengleichungen fordern, dass sie stetig aneinander anschließen, also

  $u_{i-1}(x){|}_{\text{unten}}=u_i(x){|}_{\text{oben}}$

• C→D Anpassen der Trial-Funktion an Übergangsbedingungen: Die 2N Koeffizienten a_i,1, a_i,0 können wir nicht anschaulich interpretieren - und das ist bei komplexen Zusammenhängen wie hier immer ein Nachteil. Das ändert sich, wenn wir eine Koordinatentransformation auf eine lokale Element-Koordinate ξ_i durchführen:

  ${\overline{u}}_i(\xi )=b_{i,1}\cdot {\xi }_i+b_{i,1}\text{ mit }x=x_i+{\ell }_E\cdot {\xi }_i,x_i=(i-1)\cdot {\ell }_E$.

  ξ_i ist nun die lokale, dimensionslose Ortskoordinate im Element. Für die Geradenfunktion müssen wir nun noch sicherstellen, dass die Übergangsbedingungen zwischen den Elementen erfüllt sind, also

  $\begin{array}{rcl}{\overline{u}}_i(1)& =& {\overline{u}}_{i+1}(0)\\ \text{ also}\\ b_{i,1}+b_{i,0}& =& b_{i+1,0}\end{array}$

  Das heißt: neben den Gleichgewichtsbedingungen müssen wir nun auch N-1 algebraische Nebenbedingungen zwischen den b_i,1, b_i,0 erfüllen, damit an den Knoten die geometrischen Übergangsbedingungen erfüllt sind (damit zwischen den Elementen im Verschiebungszustand keine Lücke klafft). Das ist extrem aufwändig wenn Sie bedenken, dass FE-Programme Millionen von Freiheitsgraden = Koordinaten haben können. Also ist das so nicht gewollt. Elegant wird es, wenn wir statt der b_ij die Knoten-Verschiebungen als Koordinaten nutzen: ξ_i = 0: der obere Rand und ξ_i = 1: der untere Rand also

  ${\overline{u}}_i(x)=\underset{=b_{i,1}}{\underbrace{(U_i-U_{i-1})}}\cdot \xi +\underset{=b_{i,0}}{\underbrace{U_{i-1}}}$

  Denn nun erhalten wir

  $\tilde{u}(x)=\sum _{N+1}{U}_n\cdot {\varphi }_n(x)$

  mit den N+1 Koordinaten der Knoten-Verschiebungen U_n als gesuchten Größen.
• D→E Einführung der Trial-Funktionen: Jetzt fehlt nur noch eine Kleinigkeit: wir sortieren die Gleichung für die Element-Verschiebung nach den Knoten-Variablen um und erhalten

  ${\overline{u}}_i(x)=U_{i-1}\cdot \underset{=:{\varphi }_1}{\underbrace{(1-{\xi }_i)}}+U_i\cdot \underset{=:{\varphi }_2(\xi )}{\underbrace{{\xi }_i}}$

  so dass wir uns den Verschiebungsverlauf nun als Linearkombination der Funktionen ϕ₁ und ϕ₂ denken können. Das zeigt das folgende Bild:

  

  Überlagerung der Trial-Functions.

  Die resultierenden Form-Funktionen der Koordinaten U_n der Verschiebung sind in E aufgetragen.

STEP A: Lageplan des Systems	Step B: Diskretisierung in Elemente und Knoten	Step C: Lineare Trial-Funktion im Element 4	Step D: Alle Trial-Funktionen auf einen Blick	Step E: Beitrag der einzelnen Knoten und Ihrer Form-Funktionen.

Numerische Ergebnisse finden Sie in Aufgabe FEAD.

Die virtuelle Formänderungs-Energie je Dehnstab-Element für lineare Ansatz-Polynome

Die gesamte Formänderungsenergie setzen wir jetzt aus aus den Formänderungsenergien je Element zusammen, also

$\delta \mathrm{\Pi }=\sum _I\delta {\mathrm{\Pi }}_i$

also für eine Dehnstab-Element

$\delta {\mathrm{\Pi }}_i=\underset{{\ell }_i}{\int }EA{u_i}^{\prime }\cdot \delta {u_i}^{\prime }d{x}_i$

Virtual Strain Energy

Wir setzten also die Ansatzfunktionen von oben in dieses Integral ein und erhalten

$\delta {\mathrm{\Pi }}_i=(\delta U_i,\delta U_{i-1})\cdot \frac{EA}{{\ell }_i}\left(\begin{array}{cc}+1& -1\\ -1& +1\end{array}\right)\cdot \left(\begin{array}{c}{U}_{i-1}\\ U_i\end{array}\right)$

Hier die Maxima-Kopiervorlage:

phi : [ 1-xi,
          xi];
Q[i] : [ U[i-1], U[i]];
K[i] : EA[i]/l[i] * matrix(
           [+1,-1],
           [-1,+1]);

/* Maxima ( continued) */
declare("δW", alphabetic);
declare("δΠ", alphabetic);
declare("δu", alphabetic);
declare("δU", alphabetic);

phi : [ 1-xi, xi];
Q : [[ U[i-1],  U[i]],
     [δU[i-1], δU[i]]];

trials: [u = Q[1].phi, δu = Q[2].phi];

PvV : δΠ[i] = 'integrate(EA*'diff(u,x[i])*'diff(δu,x[i]),x[i],0,l[i]);
PvV : subst([xi=x[i]/l[i]],subst(trials,PvV));
PvV : expand(ev(PvV,nouns));
K : funmake('matrix,
        makelist(
              makelist(
                   coeff(coeff(rhs(PvV),Q[2][i]),Q[1][j]),
                                      i,1,2),j,1,2));

print(δΠ[i],"=", Q[2],"*",EA/l[i],"*",K/(EA/l[i]),"*",transpose(Q[1]))$

Die virtuelle Arbeit der d'Alembert'schen Trägheitskräfte je Dehnstab-Element für lineare Ansatz-Polynome

Die gesamte virtuelle Arbeit der d'Alembert'sche Trägheitskräfte setzen wir jetzt aus aus den virtuellen Arbeiten je Element zusammen, also

$\delta W^a=\sum _I\delta W_i^a$

also für eine Dehnstab-Element

$\delta W_i^a=\underset{{\ell }_i}{\int }\varrho A(-{\ddot{u}}_i)\cdot \delta u_{i}d{x}_i$

Virtual Work of d'Alembert Forces

Wir setzten also die Ansatzfunktionen von oben in dieses Integral ein und erhalten

$\delta W_i^a=-(\delta U_i,\delta U_{i-1})\cdot \frac{\varrho A{\ell }_i}{3}\left(\begin{array}{cc}1& \frac{1}{2}\\ \frac{1}{2}& 1\end{array}\right)\cdot \left(\begin{array}{c}{\ddot{U}}_{i-1}\\ {\ddot{U}}_i\end{array}\right)$

Hier die Maxima-Kopiervorlage:

phi : [ 1-xi,
          xi];
Q[i] : [ U[i-1](t), U[i](t)];
M[i] : rho*A*l[i]/3 * matrix(
           [ 1 ,1/2],
           [1/2, 1 ]);

/* Maxima ( continued) */
declare("δW", alphabetic);
declare("δΠ", alphabetic);
declare("δu", alphabetic);
declare("δU", alphabetic);

phi : [ 1-xi, xi];
Q : [[ U[i-1](t),  U[i](t)],
     [δU[i-1]   , δU[i]   ]];

trials: [u(x[i],t) = Q[1].phi, δu(x[i]) = Q[2].phi];

PvV : δW^a[i] = -'integrate(rho*A*'diff(u(x[i],t),t,2) * δu(x[i]), x[i],0,l[i]);
PvV : subst(subst([xi=x[i]/l[i]],trials),PvV);
PvV : expand(ev(PvV,nouns));
M : funmake('matrix,
        makelist(
              makelist(
                   coeff(coeff(rhs(-PvV),Q[2][i]),'diff(Q[1][j],t,2)),
                                      i,1,2),j,1,2));

print(δW^a[i],"=", Q[2],"*",rho*A*l[i]/3,"*",M/(rho*A*l[i]/3),"*",transpose(diff(Q[1],t,2)))$