3 – Finite Elemente Diskretisierung

 

3.1 Generierung einer Struktur aus Elementen

Wir betrachten ein Raumfachwerk:

raumfachwerk-finite-elemente-methode

Für Raumfachwerke rechnen wir mit dem Gebiet

\Omega = \bigcup\limits_{e = 1}^{{n_e}} {{\Omega ^e}} ,\quad {\Omega ^i} \cap {\Omega ^j} = \emptyset \;\;\forall i \ne j.

Wir definieren weiterhin:

\delta {W_{dyn}} = \sum\limits_{e = 1}^{{n_e}} {\delta W_{dyn}^e} ,\quad \delta {W_{int}} = \sum\limits_{e = 1}^{{n_e}} {\delta W_{int}^e} ,\quad \delta {W_{ext}} = \sum\limits_{e = 1}^{{n_e}} {\delta W_{ext}^e}

\delta {W_{dyn}}+\delta {W_{int}} = \delta {W_{ext}}

Im folgenden wollen wir uns auf eindimensionale Probleme beschränken.

3.2 Approximation von Variablen eindimensionaler Kontinua

Das unbekannte exakte Verschiebungsfeld {u_1} wird durch {\tilde u_1} approximiert. {\tilde u_1} besteht aus einer Ansatzfunktion und diskreten Lösungspunkten. Als Ansatzfunktion haben sich polynomiale Funktionen etabliert. Ebenso werden {\ddot u_1}, \delta {u_1} und {x_1} diskretisiert.

3.2.1 Bedingungen an die Ansatzfunktionen

Die Ansatzfunktionen

  • müssen konform sein (Interaktion und Kompatibilität benachbarter finiter Elemente).
  • müssen mindestens konstante Verzerrungen \varepsilon liefern.
  • dürfen bei Starrkörperbewegungen keine Verzerrungen liefern.

3.2.2 Physikalische und natürliche Koordinaten

Wir betrachten folgendes Stabelement:

physikalische-naturliche-koordinaten-stab

In den physikalischen Koordinaten gilt {x_1} \in \left[ {-\frac{L}{2};\frac{L}{2}} \right]. Es ist oft sinnvoll, die Koordinaten unabhängig von der Elementlänge anzugeben. Hierfür verwenden wir die natürlichen Koordinaten. In diesen gilt dann {\xi _1} \in \left[ {-1;1} \right]. Der Zusammenhang zwischen den verschiedenen Koordinaten lautet {x_1} = \frac{L}{2}{\xi _1}. Für die Transformation von differentiellen Linienelementen d{x_1} und d{\xi _1} brauchen wir die Jakobimatrix (Ableitung nach den natürlichen Koordinaten):

J = {x_{1;1}}: = \frac{{\partial {x_1}}}{{\partial {\xi _1}}} = \frac{L}{2}

Dies ist hier eine 1×1 Matrix.

Lokale polynomiale Approximation

Wir verwenden Lagrange-Polynome des Grundtyps

{u_1}\left( {{\xi _1}} \right) \approx {\tilde u_1}\left( {{\xi _1}} \right) = \sum\limits_{j = 0}^p {{\alpha _j}{{\left( {{\xi _1}} \right)}^j}} = {\alpha _0} + {\alpha _1}\left( {{\xi _1}} \right) + {\alpha _2}{\left( {{\xi _1}} \right)^2} + \ldots + {\alpha _p}{\left( {{\xi _1}} \right)^p}.

Die Koeffizienten {\alpha _j} sind so zu wählen, dass {\tilde u_1} an den durch \xi _1^i gekennzeichneten Elementknoten i den Verschiebungen {\tilde u_1}\left( {\xi _1^i} \right) = u_1^{{e_i}} entspricht.

3.2.3 Konzept der analytischen Lösung (Zielvorstellung)

Vorgehen bei der Lösung:

  • Unterteilung in finite Elemente
  • Definition von Systemknotenverschiebungen u_1^j als Parameter stückweise linearer Verschiebungsapproximationen {\tilde u_1}
  • Definition von Elementen, Elementknoten und Elementknotenverschiebungen u_1^{{e_1}}
  • Aufbau von Verschiebungsapproximationen

Wir zeigen nun das Vorgehen am Beispiel des hier abgebildeten Balkens:

analytische-losung-finite-elemente-vergleich-balken

Unterteilung in Elemente, Systemverschiebungen und Verschiebungsapproximation:

finite-elemente-knoten-systemverschiebung-approximation

Assoziierte Elementknotenverschiebung und Verschiebungsapproximation:

assoziierte-elementknotenverschiebung-approximation

Wir kommen nun zur Verschiebungsapproximation. Für jeden Knoten eines Elements gibt es eine eigene Approximationsfunktion. Als Beispiel betrachten wir das mittlere Element e:

approximationsfunktion-beispiel-finite-elemente

Diese Approximationsfunktionen werden zur Gesamtapproximation überlagert:

gesamtverschiebung-approximation-uberlagerung

Die Approximationsfunktionen \tilde u_1^i bestehen aus Ansatzfunktionen {N^i} (hier linear), die mit einem konstanten Faktor multipliziert werden:

ansatzfunktion-finite-elemente

Formulierung in physikalischen Koordinaten {x_1}:

formulierung-physikalische-koordinaten

Formulierung in natürlichen Koordinaten {\xi _1}:

formulierung-naturliche-koordinaten

Ansatz für {N^1}\left( {{\xi _1}} \right):

{N^1}\left( {{\xi _1}} \right) = {\alpha _0}+{\alpha _1}{\xi _1}

{N^1}\left( {-1} \right)\mathop = \limits^! 1,\quad {N^1}\left( 1 \right)\mathop = \limits^! 0

\Rightarrow \quad {\alpha _0} = \frac{1}{2},\quad {\alpha _1} = -\frac{1}{2}

\Rightarrow \quad {N^1}\left( {{\xi _1}} \right) = \frac{1}{2}\left( {1-{\xi _1}} \right)

Analog:

{N^2}\left( {{\xi _1}} \right) = \frac{1}{2}\left( {1+{\xi _1}} \right)

Dabei sind {N^1}\left( {{\xi _1}} \right) und {N^2}\left( {{\xi _1}} \right) lineare Ansatzfunktionen. Sie erfüllen die Interpolationseigenschaften, d.h. Wert “1″ am zugehörigen Elementknoten und Wert “0″ an anderen Knoten. Die zusammengesetzte lineare Approximation des Verschiebungsfeldes lautet nun:

{\tilde u_1}\left( {{\xi _1}} \right) = \sum\limits_{i = 1}^2 {{N^i}\left( {{\xi _1}} \right)u_1^{{e_i}}} = {N^1}\left( {{\xi _1}} \right)u_1^{{e_1}}+{N^2}\left( {{\xi _1}} \right)u_1^{{e_2}}

Statt der linearen Ansatzfunktion werden oft Lagrange-Interpolationspolynome vom Grad p verwendet. Diese lauten:

{N^i}\left( {{\xi _1}} \right) = \prod\limits_{k = 1,k \ne i}^{p+1} {\frac{{\xi _1^k-{\xi _1}}}{{\xi _1^k-\xi _1^i}}}

Die Interpolationseigenschaft ist auch hier erfüllt:

{N^i}\left( {\xi _1^k} \right) = \left\{ {\begin{array}{*{20}{c}} 1&{i = k} \\ 0&{i \ne k} \end{array}} \right.

Die allgemeine Approximation des Verschiebungsfeldes vom Polynomgrad p lautet:

{\tilde u_1}\left( {{\xi _1}} \right) = \sum\limits_{i = 1}^{p+1} {{N^i}\left( {{\xi _1}} \right)u_1^{{e_i}}}

approximation-verschiebungsfeld-polynom-grad

3.2.4 Approximation der Ableitung von primären Variablen

Die Verzerrung ist die Ableitung der Verschiebung nach der physikalischen oder natürlichen Koordinate:

\frac{{\partial {u_1}\left( {{x_1}} \right)}}{{\partial {x_1}}} = :{u_{1,1}}\left( {{x_1}} \right) \approx {{\tilde u}_{1,1}}\left( {{x_1}} \right) = \sum\limits_{i = 1}^{p+1} {N_{;1}^i\left( {{x_1}} \right)u_1^{{e_i}}}

\frac{{\partial {u_1}\left( {{\xi _1}} \right)}}{{\partial {\xi _1}}} = :{u_{1;1}}\left( {{\xi _1}} \right) \approx {{\tilde u}_{1;1}}\left( {{\xi _1}} \right) = \sum\limits_{i = 1}^{p+1} {N_{;1}^i\left( {{\xi _1}} \right)u_1^{{e_i}}}

Allgemein gilt:

N_{;1}^i\left( {{\xi _1}} \right) = \sum\limits_{l = 1,l \ne i}^{p+1} {\frac{{-1}}{{\xi _1^l-\xi _1^i}}\prod\limits_{k = 1,k \ne i,k \ne l}^{p+1} {\frac{{\xi _1^k-{\xi _1}}}{{\xi _1^k-\xi _1^i}}} }

Für p = 1 gilt:

N_{;1}^1\left( {{\xi _1}} \right) = -\frac{1}{2},\quad \quad N_{;1}^2\left( {{\xi _1}} \right) = \frac{1}{2}

3.2.5 Approximation des Ortes

Hier nutzen wir den linearen Zusammenhang zwischen {x_1} und {\xi _1}.

approximation-ort-naturlich-physikalisch

{x_1}\left( {{\xi _1}} \right) = x_1^{{e_1}}+\frac{{x_1^{{e_2}}-x_1^{{e_1}}}}{2}\left( {{\xi _1}+1} \right)

\Rightarrow \quad {x_1}\left( {{\xi _1}} \right) = \frac{1}{2}\left( {1-{\xi _1}} \right)x_1^{{e_1}}+\frac{1}{2}\left( {1+{\xi _1}} \right)x_1^{{e_2}}

\Rightarrow \quad {x_1}\left( {{\xi _1}} \right) = {N^1}\left( {{\xi _1}} \right)x_1^{{e_1}} + {N^2}\left( {{\xi _1}} \right)x_1^{{e_2}}

Allgemeine Approximation des Ortes:

{\tilde x_1}\left( {{\xi _1}} \right) = \sum\limits_{i = 1}^{p+1} {{N^i}\left( {{\xi _1}} \right)x_1^{{e_i}}}

Dies ist identisch mit der Approximation des Verschiebungsfeldes. Diese Tatsache nennen wir isoparametrisches Konzept. Aus dem Ansatz folgt ebenso die Jacobi-Matrix (hier 1×1)

\frac{{\partial {x_1}\left( {{\xi _1}} \right)}}{{\partial {\xi _1}}} = \frac{L}{2},

oder invers

\frac{{\partial {\xi _1}\left( {{x_1}} \right)}}{{\partial {x_1}}} = \frac{2}{L}

Umformung:

{x_1}\left( {{\xi _1}} \right) = \frac{{x_1^{{e_1}}+x_1^{{e_2}}}}{2}+\frac{L}{2}{\xi _1}

3.2.6 Allgemeine Approximation weiterer Größen für lineare Elemente

Verschiebungsfeld:

{\tilde u_1}\left( {{\xi _1}} \right) = \sum\limits_{i = 1}^{p+1} {{N^i}\left( {{\xi _1}} \right)u_1^{{e_i}}}

Variation des Verschiebungsfeldes:

\delta {\tilde u_1}\left( {{\xi _1}} \right) = \sum\limits_{i = 1}^{p+1} {{N^i}\left( {{\xi _1}} \right)\delta u_1^{{e_i}}}

Beschleunigungsfeld:

{\ddot \tilde u_1}\left( {{\xi _1}} \right) = \sum\limits_{i = 1}^{p+1} {{N^i}\left( {{\xi _1}} \right)\ddot u_1^{{e_i}}}

Verzerrungen:

{{\tilde \varepsilon }_{11}}\left( {{\xi _1}} \right) = \frac{{\partial {{\tilde u}_1}\left( {{\xi _1}\left( {{x_1}} \right)} \right)}}{{\partial {x_1}}} = \frac{\partial }{{\partial {x_1}}}\sum\limits_{i = 1}^{p+1} {{N^i}\left( {{\xi _1}\left( {{x_1}} \right)} \right)u_1^{{e_i}}}

= \sum\limits_{i = 1}^{p+1} {\frac{{\partial {N^i}\left( {{\xi _1}\left( {{x_1}} \right)} \right)}}{{\partial {x_1}}}u_1^{{e_i}}}

= \sum\limits_{i = 1}^{p+1} {N_{,1}^i\left( {{\xi _1}} \right)u_1^{{e_i}}}

Bemerkung:

\frac{{\partial {N^i}\left( {{\xi _1}} \right)}}{{\partial {\xi _1}}} = N_{;1}^i\left( {{\xi _1}} \right)

\Rightarrow \quad N_{,1}^i\left( {{\xi _1}\left( {{x_1}} \right)} \right) = \frac{{\partial {N^i}\left( {{\xi _1}\left( {{x_1}} \right)} \right)}}{{\partial {\xi _1}}}\frac{{\partial {\xi _1}\left( {{x_1}} \right)}}{{\partial {x_1}}} = N_{;1}^i\left( {{\xi _1}} \right)\frac{2}{L}

Variation der Verzerrungen:

\delta {\tilde \varepsilon _{11}}\left( {{\xi _1}} \right) = \sum\limits_{i = 1}^{p+1} {N_{,1}^i\left( {{\xi _1}} \right)\delta u_1^{{e_i}}}

Für p = 1 gelten folgende Ableitungen:

N_{;1}^1\left( {{\xi _1}} \right) = \frac{\partial }{{\partial {\xi _1}}}\left[ {\frac{1}{2}\left( {1-{\xi _1}} \right)} \right] = -\frac{1}{2},\quad \quad N_{,1}^1\left( {{\xi _1}} \right) = N_{;1}^1\left( {{\xi _1}} \right)\frac{2}{L} = -\frac{1}{L}

N_{;1}^2\left( {{\xi _1}} \right) = \frac{\partial }{{\partial {\xi _1}}}\left[ {\frac{1}{2}\left( {1+{\xi _1}} \right)} \right] = +\frac{1}{2},\quad \quad N_{,1}^2\left( {{\xi _1}} \right) = N_{;1}^2\left( {{\xi _1}} \right)\frac{2}{L} = +\frac{1}{L}

3.3 Approximation der schwachen Form

Wir kommen auf unser Beispiel aus dem letzten Kapitel zurück:

approximation-schwache-form-beispiel-stab

Schwache Form des Gesamtsystems:

\int\limits_0^L {\delta {u_1}\rho {{\ddot u}_1}Ad{x_1}} +\int\limits_0^L {\delta {\varepsilon _{11}}EA{\varepsilon _{11}}d{x_1}} = \int\limits_0^L {\delta {u_1}\rho {b_1}Ad{x_1}}

Schwache Form des Elements e in physikalischen Koordinaten:

\delta W_{dyn}^e+\delta W_{int}^e = \delta W_{ext}^e

\int\limits_{x_1^{{e_1}}}^{x_1^{{e_2}}} {\delta {u_1}\rho {{\ddot u}_1}Ad{x_1}} +\int\limits_{x_1^{{e_1}}}^{x_1^{{e_2}}} {\delta {\varepsilon _{11}}EA{\varepsilon _{11}}d{x_1}} = \int\limits_{x_1^{{e_1}}}^{x_1^{{e_2}}} {\delta {u_1}\rho {b_1}Ad{x_1}}

Schwache Form des Elements e in natürlichen Koordinaten:

d{x_1} = d{\xi _1}\frac{L}{2}

\int\limits_{-1}^1 {\delta {u_1}\rho {{\ddot u}_1}\frac{L}{2}Ad{\xi _1}} +\int\limits_{-1}^1 {\delta {\varepsilon _{11}}EA{\varepsilon _{11}}\frac{L}{2}d{\xi _1}} = \int\limits_{-1}^1 {\delta {u_1}\rho {b_1}\frac{L}{2}Ad{\xi _1}}

Approximation dieser schwachen Form:

\int\limits_{-1}^1 {\delta {{\tilde u}_1}\rho {{\ddot \tilde u}_1}\frac{L}{2}Ad{\xi _1}} +\int\limits_{-1}^1 {\delta {{\tilde \varepsilon }_{11}}EA{{\tilde \varepsilon }_{11}}\frac{L}{2}d{\xi _1}} = \int\limits_{-1}^1 {\delta {{\tilde u}_1}\rho {b_1}\frac{L}{2}Ad{\xi _1}}

\delta \tilde W_{dyn}^e+\delta \tilde W_{int}^e = \delta \tilde W_{ext}^e

3.4 Approximation der inneren virtuellen Arbeit

Wir approximieren die innere virtuelle Arbeit durch

\delta \tilde W_{int}^e = \int\limits_{ - 1}^1 {\delta {{\tilde \varepsilon }_{11}}EA{{\tilde \varepsilon }_{11}}\frac{L}{2}d{\xi _1}}

= \int\limits_{-1}^1 {\left[ {\sum\limits_{i = 1}^{p+1} {N_{,1}^i\left( {{\xi _1}} \right)\delta u_1^{{e_i}}} } \right]EA\left[ {\sum\limits_{j = 1}^{p+1} {N_{,1}^j\left( {{\xi _1}} \right)u_1^{{e_j}}} } \right]\frac{L}{2}d{\xi _1}}

Wir ziehen nun die Summe aus dem Integral heraus:

\delta \tilde W_{int}^e = \sum\limits_{i = 1}^{p+1} {\sum\limits_{j = 1}^{p+1} {\int\limits_{-1}^1 {N_{,1}^i\left( {{\xi _1}} \right)\delta u_1^{{e_i}}EAN_{,1}^j\left( {{\xi _1}} \right)u_1^{{e_j}}\frac{L}{2}d{\xi _1}} } }

= \sum\limits_{i = 1}^{p+1} {\sum\limits_{j = 1}^{p+1} {\delta u_1^{{e_i}}u_1^{{e_j}}\underbrace {\int\limits_{-1}^1 {N_{,1}^i\left( {{\xi _1}} \right)EAN_{,1}^j\left( {{\xi _1}} \right)\frac{L}{2}d{\xi _1}} }_{Elementsteifgkeit\;\;{k^{{e_{ij}}}}}} }

= \sum\limits_{i = 1}^{p+1} {\sum\limits_{j = 1}^{p+1} {\delta u_1^{{e_i}}{k^{{e_{ij}}}}u_1^{{e_j}}} }

Es ergibt sich in Matrixdarstellung:

\delta \tilde W_{int}^e = \underbrace {\left\{ {\begin{array}{*{20}{c}} {\delta u_1^{{e_1}}} \\ {\delta u_1^{{e_2}}} \\ \vdots \\ {\delta u_1^{{e_{p+1}}}} \end{array}} \right\}}_{\delta {{\vec u}^e}} \cdot \underbrace {\left[ {\begin{array}{*{20}{c}} {{k^{{e_{11}}}}}&{{k^{{e_{12}}}}}& \ldots &{{k^{{e_{1\;\left( {p+1} \right)}}}}} \\ {{k^{{e_{21}}}}}& \ddots &{}& \vdots \\ \vdots &{}& \ddots & \vdots \\ {{k^{{e_{\left( {p+1} \right)\;1}}}}}& \ldots & \ldots &{{k^{{e_{\left( {p+1} \right)\left( {p+1} \right)}}}}} \end{array}} \right]}_{{\underline{k} }^e}\underbrace {\left\{ {\begin{array}{*{20}{c}} {u_1^{{e_1}}} \\ {u_1^{{e_2}}} \\ \vdots \\ {u_1^{{e_{p+1}}}} \end{array}} \right\}}_{{{\vec u}^e}}

\Rightarrow \quad \delta \tilde W_{int}^e = \delta {{\vec u}^e} \cdot {\underline{k} }^e{{\vec u}^e}

Beispiel: lineares finites Element

{k^{{e_{11}}}} = \int\limits_{-1}^1 {N_{,1}^1\left( {{\xi _1}} \right)EAN_{,1}^1\left( {{\xi _1}} \right)\frac{L}{2}d{\xi _1}} = \int\limits_{-1}^1 {\left( {-\frac{1}{L}} \right)EA\left( {-\frac{1}{L}} \right)\frac{L}{2}d{\xi _1}} = \frac{{EA}}{L}

{k^{{e_{12}}}} = {k^{{e_{21}}}} = \int\limits_{-1}^1 {N_{,1}^1\left( {{\xi _1}} \right)EAN_{,1}^2\left( {{\xi _1}} \right)\frac{L}{2}d{\xi _1}} = \int\limits_{-1}^1 {\left( {-\frac{1}{L}} \right)EA\frac{1}{L}\frac{L}{2}d{\xi _1} = -\frac{{EA}}{L}}

{k^{{e_{22}}}} = \int\limits_{-1}^1 {N_{,1}^2\left( {{\xi _1}} \right)EAN_{,1}^2\left( {{\xi _1}} \right)\frac{L}{2}d{\xi _1}} = \int\limits_{-1}^1 {\frac{1}{L}EA\frac{1}{L}\frac{L}{2}d{\xi _1}} = \frac{{EA}}{L}

Elementsteifigkeitsmatrix:

\underline{k} ^e = \frac{{EA}}{L}\left[ {\begin{array}{*{20}{c}} 1&{-1} \\ {-1}&1 \end{array}} \right]

3.5 Approximation der dynamischen virtuellen Arbeit

Für die dynamische virtuelle Arbeit gilt:

\delta \tilde W_{dyn}^e = \int\limits_{-1}^1 {\delta {{\tilde u}_1}\rho {{\ddot \tilde u}_1}\frac{L}{2}Ad{\xi _1}}

= \int\limits_{-1}^1 {\left[ {\sum\limits_{i = 1}^{p+1} {{N^i}\left( {{\xi _1}} \right)\delta u_1^{{e_i}}} } \right]\rho A\left[ {\sum\limits_{j = 1}^{p+1} {{N^j}\left( {{\xi _1}} \right)\ddot u_1^{{e_j}}} } \right]\frac{L}{2}d{\xi _1}}

= \sum\limits_{i = 1}^{p+1} {\sum\limits_{j = 1}^{p+1} {\delta u_1^{{e_i}}\ddot u_1^{{e_j}}\underbrace {\int\limits_{-1}^1 {{N^i}\left( {{\xi _1}} \right)\rho A{N^j}\left( {{\xi _1}} \right)\frac{L}{2}d{\xi _1}} }_{Elementmassenmatrix\;\;{\underline{m} }^{{e_{ij}}}}} }

= \sum\limits_{i = 1}^{p+1} {\sum\limits_{j = 1}^{p+1} {\delta u_1^{{e_i}}{m^{{e_{ij}}}}\ddot u_1^{{e_j}}} }

Matrixdarstellung:

\delta \tilde W_{dyn}^e = \underbrace {\left\{ {\begin{array}{*{20}{c}} {\delta u_1^{{e_1}}} \\ {\delta u_1^{{e_2}}} \\ \vdots \\ {\delta u_1^{{e_{p+1}}}} \end{array}} \right\}}_{\delta {{\vec u}^e}} \cdot \underbrace {\left[ {\begin{array}{*{20}{c}} {{m^{{e_{11}}}}}&{{m^{{e_{12}}}}}& \ldots &{{m^{{e_{1\;\left( {p+1} \right)}}}}} \\ {{m^{{e_{21}}}}}& \ddots &{}& \vdots \\ \vdots &{}& \ddots & \vdots \\ {{m^{{e_{\left( {p+1} \right)\;1}}}}}& \ldots & \ldots &{{m^{{e_{\left( {p+1} \right)\left( {p+1} \right)}}}}} \end{array}} \right]}_{{\underline{m} }^e}\underbrace {\left\{ {\begin{array}{*{20}{c}} {\ddot u_1^{{e_1}}} \\ {\ddot u_1^{{e_2}}} \\ \vdots \\ {\ddot u_1^{{e_{p+1}}}} \end{array}} \right\}}_{{{\ddot {\vec u}}^e}}

\Rightarrow \quad \delta \tilde W_{dyn}^e = \delta {{\vec u}^e} \cdot {\underline{m} }^e{{\ddot {\vec u}}^e}

Beispiel für p = 1:

{m^{{e_{11}}}} = \int\limits_{-1}^1 {{N^1}\left( {{\xi _1}} \right)\rho A{N^1}\left( {{\xi _1}} \right)\frac{L}{2}d{\xi _1}}

= \frac{{\rho LA}}{8}\int\limits_{-1}^1 {{{\left( {1-{\xi _1}} \right)}^2}d{\xi _1}} = \frac{{\rho LA}}{8}\frac{8}{3}

{m^{{e_{12}}}} = {m^{{e_{21}}}} = \frac{{\rho LA}}{8}\int\limits_{-1}^1 {\left( {1-\xi _1^2} \right)d{\xi _1}} = \frac{{\rho LA}}{8}\frac{8}{6}

{m^{{e_{22}}}} = \frac{{\rho LA}}{8}\int\limits_{-1}^1 {{{\left( {1+{\xi _1}} \right)}^2}d{\xi _1}} = \frac{{\rho LA}}{8}\frac{8}{3}

\Rightarrow \quad {\underline{m} }^e = \frac{{\rho LA}}{6}\left[ {\begin{array}{*{20}{c}} 2&1 \\ 1&2 \end{array}} \right]

3.6 Approximation der externen virtuellen Arbeit

Die externe virtuelle Arbeit ist definiert als das folgende Integral:

\delta \tilde W_{ext}^e = \int\limits_{-1}^1 {\delta {{\tilde u}_1}\rho {b_1}\left( {{\xi _1}} \right)A\frac{L}{2}d{\xi _1}}

= \int\limits_{-1}^1 {\sum\limits_{i = 1}^{p+1} {{N^i}\left( {{\xi _1}} \right)\delta u_1^{{e_i}}\rho {b_1}\left( {{\xi _1}} \right)A\frac{L}{2}} \;d{\xi _1}}

= \sum\limits_{i = 1}^{p+1} {\delta u_1^{{e_i}}\underbrace {\int\limits_{-1}^1 {\;{N^i}\left( {{\xi _1}} \right)\rho {b_1}\left( {{\xi _1}} \right)A\frac{L}{2}d{\xi _1}} }_{r_1^{{e_i}}}}

= \sum\limits_{i = 1}^{p+1} {\delta u_1^{{e_i}}r_1^{{e_i}}}

Dabei ist r_1^{{e_i}} die zum Elementknoten i assoziierte Elementlast.

\delta \tilde W_{ext}^e = \sum\limits_{i = 1}^{p+1} {\delta u_1^{{e_i}}r_1^{{e_i}}} = \left\{ {\begin{array}{*{20}{c}} {\delta u_1^{{e_1}}} \\ \vdots \\ {\delta u_1^{{e_{p+1}}}} \end{array}} \right\} \cdot \left\{ {\begin{array}{*{20}{c}} {r_1^{{e_1}}} \\ \vdots \\ {r_1^{{e_{p+1}}}} \end{array}} \right\} = \delta {\vec u^e}{\vec r^e}

Hier ist {\vec r^e} der Elementlastvektor. Beispiel für p = 1:

r_1^{{e_1}} = \int\limits_{-1}^1 {{N^1}\left( {{\xi _1}} \right)\rho A{b_1}\left( {{\xi _1}} \right)\frac{L}{2}d{\xi _1}}

r_1^{{e_2}} = \int\limits_{-1}^1 {{N^2}\left( {{\xi _1}} \right)\rho A{b_1}\left( {{\xi _1}} \right)\frac{L}{2}d{\xi _1}}

Die Linienlast {b_1}\left( {{\xi _1}} \right) nehmen wir als konstant an: {b_1}\left( {{\xi _1}} \right) = \operatorname{const}. Für den Elementlastvektor ergibt sich dann:

r_1^{{e_1}} = \int\limits_{-1}^1 {\frac{1}{2}\left( {1-{\xi _1}} \right)\rho A{b_1}\frac{L}{2}d{\xi _1}} = \frac{{\rho {b_1}AL}}{2}

r_1^{{e_2}} = \int\limits_{-1}^1 {\frac{1}{2}\left( {1+{\xi _1}} \right)\rho A{b_1}\frac{L}{2}d{\xi _1}} = \frac{{\rho {b_1}AL}}{2}

Oder in Vektorschreibweise:

{\vec r^e} = \left\{ {\begin{array}{*{20}{c}} {r_1^{{e_1}}} \\ {{r_1}^{{e_2}}} \end{array}} \right\} = \frac{{\rho {b_1}AL}}{2}\left\{ {\begin{array}{*{20}{c}} 1 \\ 1 \end{array}} \right\}

Für eine linear veränderliche Linienlast {b_1}\left( {{\xi _1}} \right) gilt (Beschreibung mit Hilfe der Ansatzfunktionen und Knotenwerten der Linienlast):

{b_1}\left( {{\xi _1}} \right) = {N^1}\left( {{\xi _1}} \right)b_1^{{e_1}}+{N^2}\left( {{\xi _1}} \right)b_1^{{e_2}} = \frac{1}{2}\left( {1-{\xi _1}} \right)b_1^{{e_1}}+\frac{1}{2}\left( {1+{\xi _1}} \right)b_1^{{e_2}}

r_1^{{e_1}} = \frac{{\rho LA}}{4}\int\limits_{-1}^1 {\left( {1-{\xi _1}} \right){b_1}\left( {{\xi _1}} \right)d{\xi _1}} = \ldots = \frac{{\rho LA}}{6}\left( {2b_1^{{e_1}}+b_1^{{e_2}}} \right)

r_1^{{e_2}} = \frac{{\rho LA}}{4}\int\limits_{-1}^1 {\left( {1+{\xi _1}} \right){b_1}\left( {{\xi _1}} \right)d{\xi _1}} = \ldots = \frac{{\rho LA}}{6}\left( {b_1^{{e_1}}+2b_1^{{e_2}}} \right)

Elementlastvektor für p = 1 und lineare Linienkraft:

{\vec r^e} = \frac{{\rho LA}}{6}\left\{ {\begin{array}{*{20}{c}} {2b_1^{{e_1}}+b_1^{{e_2}}} \\ {b_1^{{e_1}}+2b_1^{{e_2}}} \end{array}} \right\}

Zusammenfassung:

Wir ersetzen die verteilte Kraft durch Elementknotenkräfte (Elementlastvektor). Diese Knotenkräfte leisten die gleiche virtuelle Arbeit mit dem Feld der virtuellen Verschiebungen. Zur Umsetzung werten wir die Integrale für den Elementlastvektor aus. Niemals sollte die Last nach Gefühl auf die Knoten verteilt werden.

3.7 Approximierte schwache Form des Elementes

Die approximierte dynamische virtuelle Arbeit plus die approximierte interne virtuelle Arbeit entspricht der approximierten externen virtuellen Arbeit:

\delta \tilde W_{dyn}^e+\delta \tilde W_{int}^e = \delta \tilde W_{ext}^e

\Rightarrow \quad \delta {{\vec u}^e} \cdot {\underline{m} }^e{{\ddot {\vec u}}^e}+\delta {{\vec u}^e} \cdot {\underline{k} }^e{{\vec u}^e} = \delta {{\vec u}^e} \cdot {{\vec r}^e}

3.8 Numerische Lösung des Zugstabes mit konstanter Linienkraft

Wir betrachten folgenden Zugstab:

numerische-losung-zugstab-finite-elemente

Annahmen:

  • p = 1
  • 1 Element e
  • statisches Problem ({\ddot u_1} = 0), daher fällt die dynamische virtuelle Arbeit weg
  • Dirichlet-Randbedingung: {u_1^*}\left( 0 \right) = 0, also auch \delta {u_1^*}\left( 0 \right) = 0
  • Neumann-Randbedingung: N_1^*\left( L \right) = 0

Elementsteifigkeitsmatrix:

\underline{k} ^e = \frac{{EA}}{L}\left[ {\begin{array}{*{20}{c}} 1&{-1} \\ {-1}&1 \end{array}} \right]

Virtuelle Verschiebung:

\delta {\vec u^e} = \left\{ {\begin{array}{*{20}{c}} {\delta u_1^{{e_1}}} \\ {\delta u_1^{{e_2}}} \end{array}} \right\} = \left\{ {\begin{array}{*{20}{c}} 0 \\ {\delta u_1^{{e_2}}} \end{array}} \right\}

Elementlastvektor:

{\vec r^e} = \frac{{{p_1}L}}{2}\left\{ {\begin{array}{*{20}{c}} 1 \\ 1 \end{array}} \right\}

Verschiebungen an den Knoten:

{\vec u^e} = \left\{ {\begin{array}{*{20}{c}} {u_1^{{e_1}}} \\ {u_1^{{e_2}}} \end{array}} \right\} = \left\{ {\begin{array}{*{20}{c}} 0 \\ {u_1^{{e_2}}} \end{array}} \right\}

Alles einsetzen in das Fundamentallemma (ohne den dynamischen Anteil):

\delta {{\vec u}^e} \cdot {\underline{k} }^e{{\vec u}^e} = \delta {{\vec u}^e}{{\vec r}^e}

\Rightarrow \quad \left\{ {\begin{array}{*{20}{c}} 0 \\ {\delta u_1^{{e_2}}} \end{array}} \right\} \cdot \left[ {\begin{array}{*{20}{c}} {\frac{{EA}}{L}}&{-\frac{{EA}}{L}} \\ {-\frac{{EA}}{L}}&{\frac{{EA}}{L}} \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} 0 \\ {u_1^{{e_2}}} \end{array}} \right\} = \left\{ {\begin{array}{*{20}{c}} 0 \\ {\delta u_1^{{e_2}}} \end{array}} \right\} \cdot \left\{ {\begin{array}{*{20}{c}} {\frac{{{p_1}L}}{2}} \\ {\frac{{{p_1}L}}{2}} \end{array}} \right\}

Durch Streichen der ersten Zeile (Nullzeile) erhalten wir:

\delta u_1^{{e_2}}\frac{{EA}}{L}e_1^{{e_2}} = \delta u_1^{{e_2}}\frac{{{p_1}L}}{2}

Fundamentallemma der Variationsrechnung (\delta u_1^{{e_2}} beliebig):

\frac{{EA}}{L}u_1^{{e_2}} = \frac{{{p_1}L}}{2}\quad \Rightarrow \quad u_1^{{e_2}} = \frac{{{p_1}{L^2}}}{{2EA}}

Vergleich mit der analytischen Lösung

{u_1}\left( {{x_1}} \right) = \frac{{{p_1}}}{{2EA}}\left[ {2L{x_1}-x_1^2} \right]\quad \Rightarrow \quad {u_1}\left( {{x_1} = L} \right) = \frac{{{p_1}{L^2}}}{{2EA}}

Nachlaufrechnung:

{{\tilde u}_1}\left( {{\xi _1}} \right) = \underbrace {{N^1}\left( {{\xi _1}} \right)u_1^{{e_1}}}_{ = 0}+{N^2}\left( {{\xi _1}} \right)u_1^{{e_2}} = \frac{1}{2}\left( {1+{\xi _1}} \right)u_1^{{e_2}},\quad {\xi _1} = \frac{2}{L}{x_1}-1

\Rightarrow \quad {{\tilde u}_1}\left( {{x_1}} \right) = \frac{{{x_1}}}{L}u_1^{{e_2}} = \frac{{{p_1}L}}{{2EA}}{x_1}

\Rightarrow \quad {{\tilde u}_1}\left( {{x_1} = \frac{L}{2}} \right) = \frac{{{p_1}{L^2}}}{{4EA}}

An dieser Stelle ergibt sich aber analytisch:

{u_1}\left( {{x_1} = \frac{L}{2}} \right) = \frac{{3{p_1}{L^2}}}{{8EA}} \ne \frac{{{p_1}{L^2}}}{{4EA}}

Die lokale Form der DGL wird also verletzt, die Integrale Form stimmt aber. Skizze:

finite-elemente-approximation-analytisch-vergleich

Numerisch:

{{\tilde \varepsilon }_{11}}\left( {{x _1}} \right) = \frac{1}{2}\frac{2}{L}u_1^{{e_2}}\quad \Rightarrow \quad {{\tilde \varepsilon }_{11}}\left( {{x_1}} \right) = \frac{1}{L}u_1^{{e_2}} = \frac{{{p_1}L}}{{2EA}}

{{\tilde \sigma }_{11}}\left( {{x_1}} \right) = E{{\tilde \varepsilon }_{11}}\left( {{x_1}} \right) = \frac{{{p_1}L}}{{2A}}

{{\tilde N}_1}\left( {{x_1}} \right) = {{\tilde \sigma }_{11}}\left( {{x_1}} \right)A = \frac{{{p_1}L}}{2}

Analytisch:

{\varepsilon _{11}}\left( {{x_1}} \right) = \frac{{{p_1}}}{{EA}}\left[ {L-{x_1}} \right]

{\sigma _{11}}\left( {{x_1}} \right) = \frac{{{p_1}}}{A}\left[ {L-{x_1}} \right]

{N_1}\left( {{x_1}} \right) = {p_1}\left[ {L-{x_1}} \right]

Die analytischen Funktionen sind also linear, die numerischen hingegen konstant. Die Integrale über den Bereich sind aber gleich:

verzerrungsfeld-integral-mittel-fehler-finite-elemente-analytisch

Probe:

\int\limits_0^L {\frac{{\delta u_1^2}}{L}\frac{{Eu_1^2}}{L}d{x_1}} -\int\limits_0^L {\frac{{du_1^2}}{L}{x_1}{p_1}d{x_1}} = \delta u_1^2\left[ {\frac{{{p_1}L}}{2}-\frac{{{p_1}L}}{2}} \right] = 0

Möglichkeiten zur Verbesserung der Approximationsqualität:

  • kleinere Teilgebiete, d.h. mehr finite Elemente mit p = 1 (h-Methode)
  • höherwertige Ansätze, z.B. p = 2 (p-Methode)

3.9 Die h-Methode am Beispiel des Zugstabes mit konstanter Linienlast

Wir betrachten wieder folgenden Zugstab:

h-methode-finite-elemente-verfeinerung-beispiel-zugstab

Unterteilung in zwei finite Elemente:

h-methode-aufteilung-finite-elemente

Für das Element {e_1} gilt mit {L^{{e_1}}} = L/2:

\left\{ {\begin{array}{*{20}{c}} {\delta u_1^{{e_1},1}} \\ {\delta u_1^{{e_1},2}} \end{array}} \right\} \cdot \left[ {\begin{array}{*{20}{c}} {\frac{{EA}}{{{L^{{e_1}}}}}}&{-\frac{{EA}}{{{L^{{e_1}}}}}} \\ {-\frac{{EA}}{{{L^{{e_1}}}}}}&{\frac{{EA}}{{{L^{{e_1}}}}}} \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} {u_1^{{e_1},1}} \\ {u_1^{{e_1},2}} \end{array}} \right\} = \left\{ {\begin{array}{*{20}{c}} {\delta u_1^{{e_1},1}} \\ {\delta u_1^{{e_1},2}} \end{array}} \right\} \cdot \left\{ {\begin{array}{*{20}{c}} {\frac{{{p_1}{L^{{e_1}}}}}{2}} \\ {\frac{{{p_1}{L^{{e_1}}}}}{2}} \end{array}} \right\}

Für das Element {e_2} gilt mit {L^{{e_2}}} = L/2:

\left\{ {\begin{array}{*{20}{c}} {\delta u_1^{{e_2},1}} \\ {\delta u_1^{{e_2},2}} \end{array}} \right\} \cdot \left[ {\begin{array}{*{20}{c}} {\frac{{EA}}{{{L^{{e_2}}}}}}&{-\frac{{EA}}{{{L^{{e_2}}}}}} \\ {-\frac{{EA}}{{{L^{{e_2}}}}}}&{\frac{{EA}}{{{L^{{e_2}}}}}} \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} {u_1^{{e_2},1}} \\ {u_1^{{e_2},2}} \end{array}} \right\} = \left\{ {\begin{array}{*{20}{c}} {\delta u_1^{{e_2},1}} \\ {\delta u_1^{{e_2},2}} \end{array}} \right\} \cdot \left\{ {\begin{array}{*{20}{c}} {\frac{{{p_1}{L^{{e_2}}}}}{2}} \\ {\frac{{{p_1}{L^{{e_2}}}}}{2}} \end{array}} \right\}

Ensemblierung der Struktur

Wir haben nun diskrete Strukturverschiebungen u_1^1, u_1^2, u_1^3 und virtuelle diskrete Strukturverschiebungen \delta u_1^1, \delta u_1^2, \delta u_1^3. Der Zusammenhang zwischen System- und Elementknotenverschiebungen ist:

\delta u_1^1 = \delta u_1^{{e_1},1},\quad \delta u_1^2 = \delta u_1^{{e_1},2} = \delta u_1^{{e_2},1},\quad \delta u_1^3 = \delta u_1^{{e_2},2}

u_1^1 = u_1^{{e_1},1},\quad u_1^2 = u_1^{{e_1},2} = u_1^{{e_2},1},\quad u_1^3 = u_1^{{e_2},2}

Wir generieren nun die Systemsteifigkeitsmatrix durch Ensemblierung der Elementsteifigkeitsmatrizen:

\delta {\tilde W_{int}} = \underbrace {\left\{ {\begin{array}{*{20}{c}} {\delta u_1^2} \\ {\delta u_1^2} \\ {\delta u_1^3} \end{array}} \right\}}_{ = :\delta \vec u} \cdot \underbrace {\left[ {\begin{array}{*{20}{c}} {{k^{{e_1},11}}}&{{k^{{e_1},12}}}&0 \\ {{k^{{e_1},21}}}&{{k^{{e_1},22}}+{k^{{e_2},11}}}&{{k^{{e_2},12}}} \\ 0&{{k^{{e_2},21}}}&{{k^{{e_2},22}}} \end{array}} \right]}_{ = :\underline{K} }\underbrace {\left\{ {\begin{array}{*{20}{c}} {u_1^2} \\ {u_1^2} \\ {u_1^3} \end{array}} \right\}}_{ = :\vec u}

Hierfür kann wie in diesem Artikel beschrieben eine Indextafel verwendet werden.

Aufgrund der Randbedingung u_1^1 = 0,\;\delta u_1^1 = 0 können wir die erste Zeile streichen:

ensemblierung-kondensation-finite-elemente-vereinfachung-gleichungssystem

Es bleibt:

\delta {{\tilde W}_{int}} = \left\{ {\begin{array}{*{20}{c}} {\delta u_1^2} \\ {\delta u_1^3} \end{array}} \right\} \cdot \left[ {\begin{array}{*{20}{c}} {{k^{{e_1},22}}+{k^{{e_2},11}}}&{{k^{{e_2},12}}} \\ {{k^{{e_2},21}}}&{{k^{{e_2},22}}} \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} {u_1^2} \\ {u_1^3} \end{array}} \right\}

= \left\{ {\begin{array}{*{20}{c}} {\delta u_1^2} \\ {\delta u_1^3} \end{array}} \right\} \cdot \left[ {\begin{array}{*{20}{c}} {\frac{{4EA}}{L}}&{-\frac{{2EA}}{L}} \\ {-\frac{{2EA}}{L}}&{\frac{{2EA}}{L}} \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} {u_1^2} \\ {u_1^3} \end{array}} \right\}

Generierung des Systemlastvektors:

\vec r = \left\{ {\begin{array}{*{20}{c}} {r_1^1} \\ {r_1^2} \\ {r_1^3} \end{array}} \right\} = \left\{ {\begin{array}{*{20}{c}} {\frac{{{p_1}L}}{4}} \\ {\frac{{{p_1}L}}{4}+\frac{{{p_1}L}}{4}} \\ {\frac{{{p_1}L}}{4}} \end{array}} \right\} = \left\{ {\begin{array}{*{20}{c}} {\frac{{{p_1}L}}{4}} \\ {\frac{{{p_1}L}}{2}} \\ {\frac{{{p_1}L}}{4}} \end{array}} \right\}

Das Fundamentallemma (ohne dynamischen Anteil) führt auf folgendes lineares Gleichungssystem:

\underline{K} \vec u = \vec r

Ausgeschrieben ohne die erste Zeile, die wir wie oben beschrieben streichen konnten:

\left[ {\begin{array}{*{20}{c}} {\frac{{4EA}}{L}}&{-\frac{{2EA}}{L}} \\ {-\frac{{2EA}}{L}}&{\frac{{2EA}}{L}} \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} {u_1^2} \\ {u_1^3} \end{array}} \right\} = \left\{ {\begin{array}{*{20}{c}} {\frac{{{p_1}L}}{2}} \\ {\frac{{{p_1}L}}{4}} \end{array}} \right\}

Durch Lösen des LGS erhalten wir:

u_1^2 = \frac{3}{8}\frac{{{p_1}{L^2}}}{{EA}},\quad u_1^3 = \frac{1}{2}\frac{{{p_1}{L^2}}}{{EA}}

Dies ist identisch zur analytischen Lösung. Vergleich mit der Lösung mit einem finiten Element:

h-methode-losung-vergleich-analytisch-finite-elemente

Durch die Aufteilung in zwei Elemente ist die approximierte Lösung also deutlich genauer geworden.

Die approximierte Verschiebungsfunktion lautet im Bereich 1:

{\tilde u_1}\left( {{\xi _1}} \right) = \underbrace {{N^1}u_1^1}_0+{N^2}u_1^2 = \frac{1}{2}\left( {1+{\xi _1}} \right)u_1^2

Die normierte Variable ist dabei:

{\xi _1} = \frac{2}{{{L^{e1}}}}{x_1}-1,\quad \quad {L^{e1}} = \frac{L}{2}

Damit lautet die approximierte Lösung in physikalischen Koordinaten:

{\tilde u_1}\left( {{x_1}} \right) = \frac{1}{2}\frac{4}{L}{x_1}u_1^2 = \frac{2}{L}u_1^2{x_1},\quad \quad {x_1} \in \left[ {0,\frac{L}{2}} \right]

Bereich 2:

{{\tilde u}_1}\left( {{\xi _1}} \right) = {N^1}u_1^2+{N^2}u_1^3 = \frac{1}{2}\left( {1-{\xi _1}} \right)u_1^2+\frac{1}{2}\left( {1+{\xi _1}} \right)u_1^3

\Rightarrow \quad {{\tilde u}_1}\left( {{x_1}} \right) = \left( {1-\frac{2}{L}{x_1}} \right)u_1^2+\frac{2}{L}{x_1}u_1^3,\quad \quad {x_1} \in \left[ {0,\frac{L}{2}} \right]

Wir führen nun die Nachlaufrechnung für das Verzerrungsfeld durch. Für den Bereich 1 gilt:

{\tilde \varepsilon _{11}}\left( {{x_1}} \right) = \frac{2}{L}u_1^2 = \frac{3}{4}\frac{{{p_1}L}}{{EA}}

Dies ist ein konstantes Verzerrungsfeld. Für den Bereich 2 gilt:

{\tilde \varepsilon _{11}}\left( {{x_1}} \right) = -\frac{2}{L}u_1^2+\frac{2}{L}u_1^3 = \frac{1}{4}\frac{{{p_1}L}}{{EA}}

Diese Ergebnisse können wir mit der analytischen Lösung vergleichen.

verzerrungsfeld-vergleich-analytisch-finite-elemente

Fazit: Die Lösung ist im Integralmittel erfüllt, die lokale Form wird aber verletzt.

3.10 Die p-Methode am Beispiel des Zugstabes mit konstanter Linienlast

Wir betrachten wieder folgenden Zugstab:

numerische-losung-zugstab-finite-elemente

Wie auch bei der h-Methode fügen wir einen dritten Knoten in der Mitte ein. Dieses Mal benutzen wir aber alle drei Knoten für ein einziges finites Element. Die Knoten liegen bei:

x_1^1 = -\frac{L}{2},\quad x_1^2 = 0,\quad x_1^3 = \frac{L}{2}

\xi _1^1 = -1,\quad \xi _1^2 = 0,\quad \xi _1^3 = 1

Quadratische Ansatzfunktionen (p = 2):

{N^i}\left( {{\xi _1}} \right) = \prod\limits_{k = 1,k \ne i}^3 {\frac{{\xi _1^k-{\xi _1}}}{{\xi _1^k-\xi _1^i}}}

Also in diesem Fall:

{N^1}\left( {{\xi _1}} \right) = \frac{{\xi _1^2-{\xi _1}}}{{\xi _1^2-\xi _1^1}} \cdot \frac{{\xi _1^3-{\xi _1}}}{{\xi _1^3-\xi _1^1}} = -\frac{{{\xi _1}}}{1} \cdot \frac{{1-{\xi _1}}}{2} = \frac{1}{2}\left( {{\xi _1}-1} \right){\xi _1}

{N^2}\left( {{\xi _1}} \right) = \frac{{\xi _1^1-{\xi _1}}}{{\xi _1^1-\xi _1^2}} \cdot \frac{{\xi _1^3-{\xi _1}}}{{\xi _1^3-\xi _1^2}} = \frac{{-1-{\xi _1}}}{{-1}} \cdot \frac{{1-{\xi _1}}}{1} = 1-{\left( {{\xi _1}} \right)^2}

{N^3}\left( {{\xi _1}} \right) = \frac{{\xi _1^1-{\xi _1}}}{{\xi _1^1-\xi _1^3}} \cdot \frac{{\xi _1^2-{\xi _1}}}{{\xi _1^2-\xi _1^3}} = \frac{{-1-{\xi _1}}}{{-1-1}} \cdot \frac{{-{\xi _1}}}{{-1}} = \frac{1}{2}\left( {{\xi _1}+1} \right){\xi _1}

\vec N\left( {{\xi _1}} \right) = \left\{ {\begin{array}{*{20}{c}} {{N^1}\left( {{\xi _1}} \right)} \\ {{N^2}\left( {{\xi _1}} \right)} \\ {{N^3}\left( {{\xi _1}} \right)} \end{array}} \right\}

Approximation der physikalischen Koordinate:

{x_1}\left( {{\xi _1}} \right) \approx {{\tilde x}_1}\left( {{\xi _1}} \right) = \vec N{\left( {{\xi _1}} \right)^T}{{\vec x}^e} = \sum\limits_{i = 1}^3 {{N^i}\left( {{\xi _1}} \right)x_1^{{e_i}}}

= \frac{1}{2}\left( {{\xi _1}-1} \right){\xi _1}\left( {-\frac{L}{2}} \right)+\left( {1-{{\left( {{\xi _1}} \right)}^2}} \right) \cdot 0+\frac{1}{2}\left( {1+{\xi _1}} \right){\xi _1}\frac{L}{2}

= \frac{L}{2}{\xi _1}

Wir erhalten also für die Jacobi-Matrix:

J = \frac{{\partial {{\tilde x}_1}}}{{\partial {\xi _1}}} = \frac{L}{2}

Approximation der primären Variablen:

{{\vec u}^e} = \left\{ {\begin{array}{*{20}{c}} {u_1^{{e_1}}} \\ {u_1^{{e_2}}} \\ {u_1^{{e_3}}} \end{array}} \right\},\quad \quad \delta {{\vec u}^e} = \left\{ {\begin{array}{*{20}{c}} {\delta u_1^{{e_1}}} \\ {\delta u_1^{{e_2}}} \\ {\delta u_1^{{e_3}}} \end{array}} \right\}

{u_1}\left( {{\xi _1}} \right) \approx {{\tilde u}_1}\left( {{\xi _1}} \right) = \vec N{\left( {{\xi _1}} \right)^T}{{\vec u}^e}

\delta {u_1}\left( {{\xi _1}} \right) \approx \delta {{\tilde u}_1}\left( {{\xi _1}} \right) = \vec N{\left( {{\xi _1}} \right)^T}\delta {{\vec u}^e}

Approximation der Verzerrungen:

{\varepsilon _{11}}\left( {{\xi _1}} \right) \approx {{\tilde \varepsilon }_{11}}\left( {{\xi _1}} \right) = \frac{2}{L}{{\vec N}_{;1}}{\left( {{\xi _1}} \right)^T}{{\vec u}^e}

N_{;1}^1 = {\xi _1}-\frac{1}{2},\quad N_{;1}^2 = -2{\xi _1},\quad N_{;1}^3 = {\xi _1}+\frac{1}{2}

Virtuelle innere Arbeit:

\delta W_{int}^e = \sum\limits_{i = 1}^3 {\sum\limits_{j = 1}^3 {\delta u_1^{{e_i}}{k^{{e_{ij}}}}u_1^{{e_j}}} }

Die Elemente der Elementsteifigkeitsmatrix sind dabei:

{k^{{e_{ij}}}} = \int\limits_{-1}^1 {N_{,1}^i\left( {{\xi _1}} \right)EAN_{,1}^j\left( {{\xi _1}} \right)\frac{L}{2}d{\xi _1}}

Wir setzen ein und erhalten:

{k^{{e_{11}}}} = \int\limits_{-1}^1 {N_{,1}^1\left( {{\xi _1}} \right)EAN_{,1}^1\left( {{\xi _1}} \right)\frac{L}{2}d{\xi _1}}

= \int\limits_{-1}^1 {\frac{2}{L}\left( {{\xi _1}-\frac{1}{2}} \right)EA\frac{2}{L}\left( {{\xi _1}-\frac{1}{2}} \right)\frac{L}{2}d{\xi _1}}

= \frac{{2EA}}{L}\int\limits_{-1}^1 {{{\left( {{\xi _1}-\frac{1}{2}} \right)}^2}d{\xi _1}}

= \frac{7}{3}\frac{{EA}}{L}

Für die anderen Elemente erhalten wir analog:

{k^{{e_{12}}}} = -\frac{8}{3}\frac{{EA}}{L},\quad {k^{{e_{13}}}} = \frac{1}{3}\frac{{EA}}{L},\quad {k^{{e_{22}}}} = \frac{{16}}{3}\frac{{EA}}{L},\quad {k^{{e_{23}}}} = -\frac{8}{3}\frac{{EA}}{L},\quad {k^{{e_{33}}}} = \frac{7}{3}\frac{{EA}}{L}

Es ergibt sich die symmetrische Elementsteifigkeitsmatrix:

\underline{k} ^e = \frac{{EA}}{{3L}}\left[ {\begin{array}{*{20}{c}} 7&{-8}&1 \\ {-8}&{16}&{-8} \\ 1&{-8}&7 \end{array}} \right]

Elementlastvektor für p = 2 und {p_1} = \operatorname{const}:

r_1^{{e_i}} = \int\limits_{-1}^1 {{N^i}\left( {{\xi _1}} \right){p_1}\frac{L}{2}d{\xi _1}}

Durch Einsetzen erhalten wir:

r_1^{{e_1}} = \frac{{{p_1}L}}{4}\int\limits_{-1}^1 {\left( {{{\left( {{\xi _1}} \right)}^2}-{\xi _1}} \right)d{\xi _1}} = \frac{{{p_1}L}}{6}

Analog erhalten wir für die anderen Elemente:

r_1^{{e_2}} = \frac{{{p_1}L}}{2}\int\limits_{-1}^1 {\left( {1-{{\left( {{\xi _1}} \right)}^2}} \right)d{\xi _1}} = \frac{{4{p_1}L}}{6}

r_1^{{e_3}} = \frac{{{p_1}L}}{4}\int\limits_{-1}^1 {\left( {{\xi _1}+{{\left( {{\xi _1}} \right)}^2}} \right)d{\xi _1}} = \frac{{{p_1}L}}{6}

Als Vektor:

\vec r_1^e = \frac{{{p_1}L}}{6}\left\{ {\begin{array}{*{20}{c}} 1 \\ 4 \\ 1 \end{array}} \right\}

Aufstellen des linearen Gleichungssystems mit Hilfe des Fundamentallemmas der Variationsrechnung (hier Spezialfall: Element LGS entspricht Struktur LGS):

\left\{ {\begin{array}{*{20}{c}} {\delta u_1^{{e_1}}} \\ {\delta u_1^{{e_2}}} \\ {\delta u_1^{{e_3}}} \end{array}} \right\} \cdot \frac{{EA}}{{3L}}\left[ {\begin{array}{*{20}{c}} 7&{-8}&1 \\ {-8}&{16}&{-8} \\ 1&{-8}&7 \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} {u_1^{{e_1}}} \\ {u_1^{{e_2}}} \\ {u_1^{{e_3}}} \end{array}} \right\} = \left\{ {\begin{array}{*{20}{c}} {\delta u_1^{{e_1}}} \\ {\delta u_1^{{e_2}}} \\ {\delta u_1^{{e_3}}} \end{array}} \right\} \cdot \frac{{{p_1}L}}{6}\left\{ {\begin{array}{*{20}{c}} 1 \\ 4 \\ 1 \end{array}} \right\}

Da die erste Verschiebung gleich Null ist, können wir die erste Zeile und die erste Spalte der Matrix streichen und erhalten so ein Gleichungssystem mit zwei Unbekannten:

\frac{{EA}}{{3L}}\left[ {\begin{array}{*{20}{c}} {16}&{-8} \\ {-8}&7 \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} {u_1^2} \\ {u_1^3} \end{array}} \right\} = \frac{{{p_1}L}}{6}\left\{ {\begin{array}{*{20}{c}} 4 \\ 1 \end{array}} \right\}

\Rightarrow \quad u_1^3 = \frac{{{p_1}{L^2}}}{{2EA}},\quad u_1^2 = \frac{3}{8}\frac{{{p_1}{L^2}}}{{EA}}

Diese Werte sind identisch mit der analytischen Lösung. Das konnten allerdings die linearen Ansatzfunktionen auch schon. Den Unterschied sehen wir erst bei der nun folgenden Nachlaufrechnung.

Verschiebungsfeld:

{u_1}\left( {{\xi _1}} \right) \approx {{\tilde u}_1}\left( {{\xi _1}} \right) = \underbrace {{N^1}\left( {{\xi _1}} \right)u_1^1}_{ = 0}+{N^2}\left( {{\xi _1}} \right)u_1^2+{N^3}\left( {{\xi _1}} \right)u_1^3

{x_1}\left( {{\xi _1}} \right) \approx {{\tilde x}_1}\left( {{\xi _1}} \right) = \frac{L}{2}{\xi _1}\quad \Rightarrow \quad {\xi _1} = \frac{2}{L}{x_1}

\Rightarrow \quad {{\tilde u}_1}\left( {{x_1}} \right) = \left( {1 - \frac{{4{{\left( {{x_1}} \right)}^2}}}{{{L^2}}}} \right)\frac{3}{8}\frac{{{p_1}{L^2}}}{{EA}} + \frac{1}{2}\left( {1 + \frac{2}{L}{x_1}} \right)\frac{2}{L}{x_1}\frac{{{p_1}{L^2}}}{{2EA}}

= \frac{{{p_1}{L^2}}}{{EA}}\left[ {\frac{3}{8}+\frac{{{x_1}}}{{2L}}-\frac{{{{\left( {{x_1}} \right)}^2}}}{{2{L^2}}}} \right],\quad {x_1} \in \left[ {-\frac{L}{2},\frac{L}{2}} \right]

Für die Verzerrungen ergibt sich:

{\tilde \varepsilon _{11}}\left( {{x_1}} \right) = \frac{{{p_1}{L^2}}}{{EA}}\left[ {\frac{1}{{2L}}-\frac{{{x_1}}}{{{L^2}}}} \right],\quad {x_1} \in \left[ {-\frac{L}{2},\frac{L}{2}} \right]

Die Verschiebung ist nun quadratisch, die Verzerrung linear. Vergleich mit analytischer Lösung:

verschiebung-verzerrung-feld-p-methode-analytisch-finite-elemente-vergleich

Approximierte und analytische Lösung stimmen genau überein. Es war also in diesem Fall besser, den Grad der Ansatzfunktion zu erhöhen, als die Anzahl der Elemente.