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

Im Vergleich zu linearen Ansatzfunktionen brauchen wir oft Trial-Functions höherer Ordnung. So ist beim Euler-Bernoulli-Balken die virtuelle Formänderungsarbeit

<math>\delta\Pi = \displaystyle \int_\ell E\,I\; w'' \cdot \delta w'' dx</math>.



weimaliges Ableiten der linearen Trial-Funktion von oben führt dazu, dass die virtuelle Formänderungsarbeit komplett verschwindet - diese Ansatzfunktionen passen dann also nicht.

Ein Polynom zweiten Grades würde passen - hat aber nur drei Koeffizienten. Und eine ungerade Anzahl von Koeffizienten kann man nicht auf zwei Knoten gleich verteilen! Wie brauchen ein Polynom mit vier Koeffizienten.

== Arbeitsschritte ==
Wir gehen die Arbeitsschritte zum Finden der Trial-Functions wie beim vorangehenden Schema durch:

<ul>
<li>'''A→B Diskretisierung:'''  Wir teilen die Struktur in ''N Stü''cke ( = Finite Elemente ) der Länge ''l''<sub>E</sub>. An jedem Ende eines Finiten Elements entsteht ein Knoten von insgesamt ''N+1'' Knoten.
</li>
<li>
'''B→C Trial-Funktion je Element:'''Je Element wählen wir ein kubisches Polynom (ein Polynom 3.ten Grades mit vier Koeffizienten) als Trial-Funktion

:: <math>\bar{w}_i(x) = \left\{\begin{array}{l}a_{i,3}\cdot x^3 + a_{i,2}\cdot x^2 + a_{i,1}\cdot x + a_{i,0} \text{ im Element } i\\0 \text{ sonst}\end{array}\right.</math>

so dass wir als Approximation nun

::<math>\tilde{w}(x) = \displaystyle \sum_{4N} a_{n} \cdot \phi_n(x)</math>

haben, mit den {''a<sub>i,3</sub>, ''a<sub>i,2</sub>, ''a''''<sub>i,1</sub>'' ,''a<sub>i,0</sub>''} als gesuchte Wichtungsfaktoren.
</li>
<li>
'''C→D Anpassen der Trial-Funktion an Übergangsbedingungen:'''  Die ''4N'' Koeffizienten {''a<sub>i,3</sub>, ''a<sub>i,2</sub>, ''a''''<sub>i,1</sub>'' ,''a<sub>i,0</sub>''} können wir wiederum nicht anschaulich interpretieren. Mit der Koordinatentransformation auf die lokale Element-Koordinate ''ξ<sub>i</sub>'' ist:

::<math>\bar{w}_i(\xi) = b_{i,3}\cdot\xi_i^3 + b_{i,2}\cdot\xi_i^2 + 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</math>
''ξ''<sub>i</sub> ist nun die lokale, dimensionslose Ortskoordinate im Element. Für das kubische Polynom müssen wir nun noch sicherstellen, dass die Übergangsbedingungen zwischen den Elementen erfüllt sind. Das ist etwas komplizierter als bei der linearen Ansatzfunktion. Sie sehen, dass bei einem Balken nicht nur die Auslenkung der Elemente an den Übergangsknoten gleich sein muss - sondern auch der Kipp-Winkel ''Φ'':

[[Datei:FEM-cubicTrialFcts-koordinaten.png|ohne|mini|Koordinaten des Stabes.]]
Die Übergangsbedingungen sind also

::<math>\bar{w}_i(1) = \bar{w}_{i+1}(0) \text{ und } \underbrace{\bar{w}'_i(1) = \bar{w}'_{i+1}(0)}_{\displaystyle \text{mit } w' = \Phi}</math>

Das heißt: neben den Gleichgewichtsbedingungen müssen wir nun auch ''2∙(N-1)'' algebraische Nebenbedingungen zwischen den ''b<sub>i,3</sub>'',  ''b<sub>i,2</sub>'',  ''b<sub>i,1</sub>'', ''b<sub>i,0</sub>'' erfüllen, damit an den Knoten die geometrischen Übergangsbedingungen erfüllt sind. Wieder ersetzen wir die ''b<sub>ij</sub>'' durch die Knoten-Verschiebungen und Kipp-Winkel ''Φ'' des Balkens in den Knoten:

::<math>\begin{array}{ll}w_i(0) &= W_{i-1}\\w'_i(0) &= \Phi_{i-1}\\w_i(1) &= W_{i}\\w'_i(1) &= \Phi_{i}\end{array}</math>

{{Template:MyCodeBlock
|title=Trial Functions ϕ<sub>i</sub>
|text=Durch Einsetzten unseres Polynoms in diese Randbedingungen können wir die vier generischen Koeffizienten bi,3,  bi,2,  bi,1, bi,0  gegen die Knoten-Koordinaten Wi-1, Φi-1, Wi, Φi tauschen und damit die Übergangsbedingungen implizit erfüllen! Wir finden

::<math>\begin{array}{ll}{{b}_{i,3}}&=\left({{\Phi}_{i-1}}+{{\Phi}_{i}}\right) \cdot {{l}_{E}}-2\cdot {{W}_{i}}+2\cdot {{W}_{i-1}}\\ {{b}_{i,2}}&=\left( -2\cdot {{\Phi}_{i-1}}-{{\Phi}_{i}}\right) \cdot {{l}_{E}}+3\cdot {{W}_{i}}-3\cdot {{W}_{i-1}}\\ {{b}_{i,1}}&={{\Phi}_{i-1}}\cdot l_E\\ {{b}_{i,0}}&={{W}_{i-1}}\end{array}</math>

und damit

::<math>\begin{array}{ll}\bar{w}_i(x) &= {{\xi}^{3}}\cdot \left( 2\cdot {{W}_{i-1}}-2\cdot {{W}_{i}}+\left( {{\Phi}_{i}}+{{\Phi}_{i-1}}\right) \cdot {{l}_{E}}\right)\\ &+{{\xi}^{2}}\cdot \left( -3\cdot {{W}_{i-1}}+3\cdot {{W}_{i}}+\left( -{{\Phi}_{i}}-2\cdot {{\Phi}_{i-1}}\right) \cdot {{l}_{E}}\right) \\&+\xi\cdot {{\Phi}_{i-1}}\cdot {{l}_{E}}\\&+{{W}_{i-1}}\end{array}</math>

mit den 2(N+1) Koordinaten der Knoten-Verschiebungen Wn,Φn  als gesuchten Größen.
|code=
<syntaxhighlight lang="notmuch" line='line' style="border:1px solid blue">
/* Maxima */
w(xi) := b[i,3]*xi^3+b[i,2]*xi^2+b[i,1]*xi+b[i,0];

/* boundary conditions */
BC: [w(0)=W[i-1],
     subst(0,xi,diff(w(xi),xi))/l[i] = Phi[i-1],
     w(1)=W[i],
     subst(1,xi,diff(w(xi),xi))/l[i] = Phi[i]];

/* solve for polynom coeff. */
sol: solve(BC,[b[i,3],b[i,2],b[i,1],b[i,0]])[1];
print('w(xi),"=",subst(sol,w(xi)))$

/* new coordinates = nodal displacements */
Q : [W[i-1],Phi[i-1],W[i],Phi[i]];
print('w(xi),"=",facsum(subst(sol,w(xi)),Q))$

/* plot trial fcts */
trials : expand(subst(sol,w(xi)));
trials: makelist(coeff(trials,Q[j]), j,1,4);
plot2d(subst([l[i]=5],trials),[xi,0,1], [xlabel, "ξ→"],
       [legend, "ϕ1", "ϕ2", "ϕ3", "ϕ4"],
       [style, [lines,3]] );
</syntaxhighlight>

}}
</li>
<li>
'''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  
 
::<math>\begin{array}{crc}\bar{w}_i(x) =&\;\; W_{i-1}\cdot &\underbrace{\left( {{\left( \xi-1\right) }^{2}}\cdot \left( 1+2\cdot \xi\right)  \right)}_{\displaystyle =:\phi_{1}}\\& + \Phi_{i-1}\cdot &\underbrace{\left( {{\left( \xi-1\right) }^{2}}\cdot \xi\cdot {{l}_{E}} \right)}_{\displaystyle =:\phi_{2}}\\ &+ W_{i}\cdot &\underbrace{\left( {{\xi}^{2}}\cdot \left( 3-2\cdot \xi\right) \right)}_{\displaystyle =:\phi_{3}}\\& + \Phi_{i}\cdot &\underbrace{\left( \left( \xi-1\right) \cdot {{\xi}^{2}}\cdot {{l}_{E}} \right)}_{\displaystyle =:\phi_{4}}\end{array}</math>,

so dass wir uns den Verschiebungsverlauf nun als Linearkombination der Funktionen ''ϕ<sub>1</sub>, ϕ<sub>2</sub>, ϕ<sub>3</sub>'' und ''ϕ<sub>4</sub>'' denken können. Die Form-Funktionen der Koordinaten sind hier für ''l<sub>E</sub>=5'' aufgetragen.
[[Datei:FEM-cubicTrialFcts-trialFunctions.png|mini|Trial-Functions.|alternativtext=|ohne]]

== Die virtuelle Formänderungs-Energie je Euler-Bernoulli-Element für kubische Ansatz-Polynome ==
Die gesamte Formänderungsenergie setzen wir jetzt aus aus den Formänderungsenergien je Element zusammen, also

<math>\delta\Pi = \displaystyle \sum_I \delta \Pi_i</math>

also für einen Eul-Bernoulli-Balken

<math>\delta\Pi_i = \displaystyle \int_{\displaystyle \ell_i} E\,I\; w_i'' \cdot \delta w_i'' dx_i</math>


{{Template:MyCodeBlock|title=Virtual Strain Energy
|text=
Wir setzten also die Ansatzfunktionen von oben in dieses Integral ein und erhalten

::<math>{{\mathit{δΠ}}_{i}}=\left({{\mathit{δW}}_{i-1}},{{\mathit{δΦ}}_{i-1}},{{\mathit{δW}}_{i}},{{\mathit{δΦ}}_{i}}\right)\cdot \displaystyle \frac{\mathit{EI}}{{{l}_{i}^{3}}} \cdot
\begin{pmatrix}
12 & 6\cdot {{l}_{i}} & -12 & 6\cdot {{l}_{i}}\\
 6\cdot {{l}_{i}} & 4\cdot {{l}_{i}^{2}} & -6\cdot {{l}_{i}} & 2\cdot {{l}_{i}^{2}}\\
 -12 & -6\cdot {{l}_{i}} & 12 & -6\cdot {{l}_{i}}\\
 6\cdot {{l}_{i}} & 2\cdot {{l}_{i}^{2}} & -6\cdot {{l}_{i}} & 4\cdot {{l}_{i}^{2}}
\end{pmatrix}\cdot \begin{pmatrix}
{{W}_{i-1}}\\
 {{\Phi}_{i-1}}\\
 {{W}_{i}}\\
 {{\Phi}_{i}}
\end{pmatrix}</math>

Hier die Maxima-Kopiervorlage:

<syntaxhighlight lang="notmuch" line='line' style="border:1px solid blue">
phi : [ 2*xi^3-3*xi^2+1,
       (xi^3-2*xi^2+xi)*l[i],
        3*xi^2-2*xi^3,
       (xi^3-xi^2)*l[i]];
Q[i] : [ W[i-1], Phi[i-1], W[i],Phi[i]];
K[i] : EI[i]/l[i]^3 * matrix(
           [    12,   6*l[i],     -12,   6*l[i]],
           [6*l[i], 4*l[i]^2, -6*l[i], 2*l[i]^2],
           [   -12,  -6*l[i],      12,  -6*l[i]],
           [6*l[i], 2*l[i]^2, -6*l[i], 4*l[i]^2]);
</syntaxhighlight>
|code=
<syntaxhighlight lang="notmuch" line='line' style="border:1px solid blue">
/* Maxima ( continued) */
declare("δW", alphabetic);
declare("δΠ", alphabetic);
declare("δw", alphabetic);
declare("δW", alphabetic);
declare("δW", alphabetic);
declare("δΦ", alphabetic);

Q : [[ W[i-1], Phi[i-1], W[i],Phi[i]],
     [δW[i-1],  δΦ[i-1],δW[i], δΦ[i]]];

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

PvV : δΠ[i] = 'integrate(EI*'diff(w,x[i],2)*'diff(δw,x[i],2),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,4),j,1,4));

print(δΠ[i],"=", Q[2],"*",EI/l[i]^3,"*",K/(EI/l[i]^3),"*",transpose(Q[1]))$
</syntaxhighlight>
}}



{{Template:MyCodeBlock|title=Virtual Work of d'Alembert Forces
|text=
::<math></math>
Hier die Maxima-Kopiervorlage:

<syntaxhighlight lang="notmuch" line='line' style="border:1px solid blue">

</syntaxhighlight>
|code=
<syntaxhighlight lang="notmuch" line='line' style="border:1px solid blue">

</syntaxhighlight>
}}