Unterschied zwischen FDM, FVM und FVM

Posted on by

Hinweis: Dieser Artikel ist unsauber und unvollständig und enthält Fehler. Ich beschäftige mich aktuell nicht mit dem Thema und habe keine Ambition ihn zu perfektionieren. Wenn du Detailfragen oder Anregungen hast, dann hinterlasse bitte trotzdem ein Kommentar, dann werde ich den Artikel an entsprechender Stelle nachbessern.

1. Einleitung

Wenn man im Internet den Unterschied zwischen den Methoden der Finiten Differenzen, der Finiten Volumina und der Finiten Elemente sucht, bekommt man viele entweder sehr komplizierte oder sehr unvollständige (falsche) Erklärungen. Ich hatte – bevor ich dieses Thema schwerpunktmäßig studiert habe – nach einfachen Erklärungen gesucht gehabt diese aber nicht gefunden. Ein grundlegendes Verständnis der Verfahren ist aber notwendig, da man sonst schwierig die tiefer gehenden Konzepte verstehen kann, die notwendig sind, damit die Erkenntnisfindung durch eine Simulation tatsächlich mehr als nur”viele bunte Bilder” ist. Dieser Artikel versucht daher ohne allzu tief mathematisch werden die konzeptionellen und praktischen Unterschiede zu erläutern.

2. Grundlagen

Hier werden wir die FDM, FVM und FEM anhand der quasistatischen Kontinuumsmechanik für kleine Verformungen erläutert. Die FEM als Konzept ist aber nicht auf Festkörpermechanik limitiert! Genauso wenig wie die FVM oder die FDM auf Fluide oder Wärmeübertragung beschränkt sind. Man kann alle 3 Methoden, FEM, FVM, FDM für Probleme der Festkörpermechanik, Strömungslehre und Wärmeübertragung, eigentlich für jedes System von Gleichungen benutzen! Ein anschauliches Beispiel hierfür ist die Software COMSOL, in welcher man unterschiedliche FE-Methoden auf alles werfen kann, darunter vordefinierte Modelle wie Strömungen, Festkörpermechanik, chemikalische Reaktionen als auch  komplett eigene Modelle (was nicht bedeutet, dass es sinnvolle Ergebnisse liefert, wenn man ungeeignetes FE-Verfahren verwendet). Für die Strömungsmechanik beispielsweise funktioniert die einfachste FEM, die Galerkin-FEM (weiter unten erklärt) nicht und man muss auf alternative Formulierungen zurückgreifen (die im Falle von COMSOL schon implementiert sind).

In diesem Artikel orientieren wir uns an der Festkörpermechanik als Spezialbeispiel.

2.1 Impulserhaltung

Für die lokale Impulserhaltung gilt für Mechanik im Allgemeinen (Festkörpermechanik und Strömungslehre)

\textrm{(1)} \quad\quad \rho \ddot{u} = \rho \vec{b} + \textrm{div}(\sigma),

dabei ist \rho = \rho(\vec{x}, t) die Dichte, \ddot{\vec{u}} = \frac{\text{D}^2\vec{u}}{\text{D}t^2} die substantielle Zeitableitung der Euler’schen Geschwindigkeit, \vec{b} die Volumenkraftdichte und \sigma der Spannungstensor. Im Falle eines Festkörpers definieren wir üblicherweise ein Verzerrungsmaß. Sofern ein Verzerrungmaß existiert enthält dieses implizit die Massebilanz (Beweis hier nicht ausgeführt), was dazu führt dass wir bei Festkörpern die Impulsbilanz nicht auswerten müssen (im Gegensatz zur Strömungsmechanik).
Ohne Massenkraftdichten, für den Sonderfall der Statik und näherungsweise für die Quasistatik ergibt sich:

\textrm{(2)} \quad\quad \vec{0} = \textrm{div}(\sigma)

Für den hyperelastischen Fall, dass der Steifigkeitstensor hauptsymmetrisch ist C^\intercal = C, das heißt C_{ijkl} = C_{klij}. Wir nehmen das Verzerrungsmaß für kleine Deformationen und kleine Rotationen \varepsilon = \frac{1}{2} \left(grad(\vec{u}) + grad(\vec{u})^T\right)

\textrm{(3)} \quad\quad \vec{0} = \textrm{div}(C:\textrm{sym}(\textrm{grad}(\vec{u}))),

wobei der Operator : als lineare Abbildung im Sinne von C_{ijkl} \sigma_{kl}~ zu verstehen ist.
Für Dirichlet-RB gilt: u = u_{\Gamma} für x \in \Gamma_\text{D}

Für Neumann-RB (also eine Randbedingung auf das betrachtete Feld an sich) gilt: \textrm{sym}(\textrm{grad}(\vec{u})) = C^{-1}:\sigma_{\Gamma} für x \in \Gamma_\text{N}

Mit \Gamma_\text{D} \cup \Gamma_\text{N} \subset \delta V wobei \delta V der Rand des zu untersuchenden Gebietes V ist. Es gilt weiterhin, dass \Gamma_\text{D} und \Gamma_\text{N} disjunkt sein müssen, Cauchy-Randbedingunen sind also nicht möglich.

2.2 Konstitutivgleichung (Materialgesetz)

2.2.1 Festkörpermechanik: Hooke’sches Gesetz

Das einfachste Gesetz ist das isotrope Hooke’schen Gesetz \sigma = 2\mu \varepsilon + \lambda sp(\varepsilon) mit den zwei Skalaren Schubmodul \mu und Kompressionsmodul \lambda.

2.2.2 Fluidmechanik: Newton’sche Konstitutivgleichung

Für linear viskoses Fließverhalten gilt \sigma = -p + 2 \mu^* \left( \frac{1}{2}\left[\text{grad}\,\dot{u} +(\text{grad}\,\dot{u})^\text{T}\right] \right) mit Druck p, der Viskosität \mu^* und materielle Geschwindigkeit v = \dot{u}

2.2.3 Wärmeübertragung: Fourrier’sches Gesetz

Für die Wärmeübetragung ist das Fourrier’sche Gesetz die Konstitutivgleichung.

Wärmestrom pro Zeit = Wärmeleitfähigkeitstensor \cdot \textrm{grad}\left(T\right).

3. FDM – Finite Differenzen

Die Finite Differenzen Methode ist die am einfachsten zu implementierende Methode und lässt sich leicht per Hand nachrechnen. Die Diskretisierung, das heißt die Aufteilung des zu untersuchenden Raumes gestaltet sich hier am einfachsten. Man teilt den Raum in eine diskrete Anzahl von Punkten ein. Um die numerische Genauigkeit zu wahren gibt es hierbei zu beachten, dass sich der Abstand zwischen zwei Punkten von Punkt zu Punkt um nicht mehr als 15% verändert. Ansonsten ist man recht frei.

Nun wertet man direkt die Impulsgleichung (3) aus. Die Ableitungen werden hierbei numerisch dargestellt über das Euler-Verfahren, die Mittelpunktregel oder gar Terme höherer Ordnung.

So gibt sich für die Kontinuumsmechanik zum Beispiel ein FDM-Verfahren in 1D mit dem isotropen Hooke nach der Mittelpunktsregel zu

\textrm{(4)} \quad\quad 0 = E \frac{\partial^2 u}{\partial x^2} \\ 0 = E \frac{ u_{i+1} - 2 u_{i} + u_{i-1} }{ (\Delta x)^2 } ,

mit u_i = u(x_i)~ und der Schrittweite \Delta x.

Man erhält dann aus Gleichung 4 ein lineares Gleichungssystem in dem man für jeden freien Knoten eine Gleichung erhält. Zusammen mit den Randbedingungen (z.B. eine Dirichlet-RB u_0 = 1 oder einer Neumann-RB \sigma_1 = \sigma_{RB} \Rightarrow \frac{u_1 - u_0}{\Delta x} = \frac{\sigma_{RB}}{E}.

Anmerkung 1: Für die Strömungslehre ist es wichtig die Bewegungsrichtung des Fluids zu kennen um ein Aufwärtsverfahren zu machen, das heißt, dass man die Differenz immer vom aktuellen Punkt zum Punkt von wo die Strömung herkommt berechnet.

Anmerkung 2: Gerade in der Strömungslehre betrachtet man i. d. R. keine quasistatischen Strömungen. Dementsprechend muss die Gleichung angepasst werden.

Vorteile

  • Die FDM ist am einfachsten zu implementieren

Nachteile

  • Die Ergebnisse einer FDM können Erhaltungssätze der klassischen Physik (Energie und Masse) verletzen, das heißt es kann zum Beispiel Masse oder Impuls aus dem nichts hinzukommen oder verschwinden.

5. FVM – Finite Volumina

Grundlage der Methode der Finiten Volumina (FVM) ist wieder die Impulsgleichung (1) in welches die Newton’sche Konstitutivgleichung eingesetzt wird sowie die Kontinuitätsgleichung. Hierbei wird nun jedoch über die Erhaltungsgleichung integriert und der Satz von Gauß angewendet.

\textrm{(5)} \quad\quad \int_V \rho \ddot{u} ~\textrm{d}V = \int_V \rho \vec{b} ~\textrm{d}V + \int_{\delta V} \sigma ~\textrm{d}A .

Der Vorteil der FVM gegenüber der FDM: die Impulsgleichung wird für jedes Kontrollvolumen erfüllt und damit – sofern die Kontrollvolumen korrekt definiert sind, d.h. ohne Hohlräume, Überschneidungen etc, ist die Impulsbilanz, Massenbilanz und Energiebilanz für das gesamte Gebiet erfüllt.

Für die FVM muss das Gebiet wieder diskretisiert werden. Dabei wird die Funktionswerte – im linearen Ansatz – immer für die Zellmittelpunkte ausgewertet und abgespeichert. Anschließend muss mit 3 Approximationen gerechnet werden:

  1. Interpolation der Ströme (für Fluidmechanik), hier zum Beispiel muss man entscheiden wie der Strom der nach Osten aus einem Volumen abfließt berechnet wird
  2. Diskretisierung der Differentialoperatoren, analog zur FDM
  3. Approximation der Oberflächenintegrale, d.h. am einfachsten eine Multiplikation einer konstanten mittleren Größe mit der entsprechenden Oberfläche,oder auch kompliziertere Vorgehen
  4. Approximation der Volumenintegrale, analog zur Approximation der Oberflächenintegrale

Vorteile

  • Erhaltungssätze der klassischen Physik sind erfüllt
  • Immer noch vergleichsweise intuitiv zu verstehen

Nachteile

  • weniger potent was das Lösungsspektrum angeht. Konkret: es können nur Lösungsfunktionen gefunden werden sofern die Lösungsfunktion zweimal differenzierbar ist.
  • aufwändige Implementierung

6. FEM – Finite Elemente

Grundlage für die FEM ist ebenso die Impulsgleichung (1). Diese Gleichung wird auch als starke Form der Impulsgleichung bezeichnet. Die FEM basiert auf der schwachen Form. Diese erhält man, indem man Gleichung (1) mit einer vektoriellen Funktion \vec{v} multipliziert und anschließend partiell integriert. Wir betrachten hier in unserem Falle der Festkörpermechanik das ganze mit der Verschiebungs-FEM, das heißt wir lösen nach der Verschiebung. Diese Funktion \vec{v} heißt Testfunktion. Wir nennen das Elementgebiet nun \Omega umd Verwechslungen mit der Geschwindigkeit und der Testfunktion zu vermeiden.

Der entscheidende Vorteil der FEM gegenüber der FVM ist folgendes:

  1.  Die Ansprüche an die Verschiebungsfunktion werden reduziert. Sie muss nur einmal differenzierbar sein und nicht zweimal (in der Gleichung (3) kommt eine zweite Ableitung nach dem Ort vor). Dadurch lässt sich ein größeres Spektrum an Problemen bearbeiten
  2. Die Elemente können fast beliebige Formen annehmen. Durch isoparametrische Elemente können fast beliebige Elementformen dargestellt werden. Während man in der FVM und FDM mehr oder minder an stark strukturierte Gitter gebunden ist und dadurch zwischen hohem Formdiskretisierungsfehler oder hoher Anzahl an Elementen abwägen muss.

Dann ergibt sich die schwache Form unter Beachtung, dass \text{div} \left( C \text{grad} u \right) = \text{div}(C) \text{grad}(u) + C \text{grad}\left(\text{grad}(u)\right) , der Annahme dass \text{div}(C) = 0 (das ist nur für ein homogenes Material erfüllt!) nach einigem umformen zu

\textrm{(6)} \quad\quad \vec{0} = \int_{\Omega}{{C:\textrm{grad}(\vec{u})\vec{v}}\cdot ~\text{d}v- \int_{\Omega}{\textrm{grad}(\vec{u})\cdot C: \textrm{grad}(\vec{v})}~\text{d}v}+\int_{\Omega}{\rho\vec{b}\cdot\vec{v}~\text{d}v} .

Anschließend legt man fest, dass die Verschiebung u in jedem Element als Summe von Ansatzfunktionen darzustellen ist. Diese Ansatzfunktionen wählt man für gewöhnlich aus dem selben Raum wie die Testfunktion und erhält damit die Galerkin-FEM. Die Galerkin-FEM funktioniert aus mathematischen Überlegeungen die mir gerade zu hoch sind nicht bei der Strömungsmechanik.

Es gilt also im eindimensionalen Fall für eine Interpolation N-ter Ordnung

\textrm{(7)}\quad\quad u(x) = \sum_{i=1}^N {u_i p_i(x) } ,

wobei u_i = u_i(x_i) diskrete Werte sind und p_i das Lagrange’sche Interpolationspolynom ist, das am i-ten Knoten 1 ist (und an den anderen 0) und p'_i seine Ableitung.
Man wählt nun die Testfunktion v(x) = p_j(x) und es ergibt sich im 1D-Fall für die Steifigkeit c

für die rechte Seite für jedes Element

\textrm{(8)}\quad\quad \sum_{i=1}^N u_i \int_{\Omega} p'_i(x) p'_j(x) \textrm{d}x + \sum_{i=1}^N { u_i \int_{\Omega} c p_i(x) p_j(x) \textrm{d}x } ,

und für die linke Seite analog

\textrm{(9)}\quad\quad \sum_{i=1}^N \int_{\Omega} \rho b_i p_j(x) \textrm{d}x ,

Es ergibt sich dann für alle Elemente zusammen ein lineares Gleichungssystem bestehend aus Steifigkeitsmatrix, Verschiebungsvektor und rechter Seite:

\textrm{(10)}\quad\quad K u = f \\ K_ij = \int_{\Omega} p'_i(x) p'_j(x) \textrm{d}x + \int_{\Omega} {c p_i(x) p_j(x) \textrm{d}x } \\ f_i = \int_{\Omega} f(x) \textrm{d}x ,

Vorteile

  • Mächtigstes Lösungsverfahren, Lösungen müssen nur einmal differenzierbar sein
  • Erhaltungssätze der klassischen Physik sind immer erfüllt

Nachteile

  • aufwändiger zu verstehen als FEM und FDM
  • aufwändiger zu implementieren als FDM

Gibt es noch andere Verfahren?

Zunächst sei hier angemerkt, dass allein obige Verfahren schon tausendfache Möglichkeiten bieten. Aber ja es gibt Verfahren die in keines der oben genannten einzuordnen sind, dazu in Kürze mehr.

Letzte Artikeländerung:

  • 2021-01-31: Zersörtes Latex repariert
  • 2021-01-Mitte: Vorteil-Nachteil Tabelle eingefügt
  • Korrektur von Formel (6), sowie ergänzende Erklärung. Außerdem Inhaltsverzeichnis für den Artikel.

 

Adams-Bashforth-Coefficients

Posted on by

Das Adam-Bashforth Verfahren ist ein mehrstufiges Zeitintegrationsverfahren mit

$latex y(t_{n+1}) = y(t_n) + \Delta t \sum_{j=0}^m b_j f(t_{n-j},~y(t_{n-j}))$

Dabei sind $latex b_j$ die AB-Koeffizienten. Und $latex f(t, y(t)) = \frac{\partial y}{\partial t}(t,y(t))$. The AB-Coefficients can be calculated in the following way:

$latex b_j = \frac{(-1)^j}{j!(m-j)!} \int_0^1 \prod_{i=0,i\neq j}^{m} (v+i) \text{d}v$

 Table of Adams-Bashforth-Coefficients up to order 10

The following table contains the resulting coefficients calculated with the formula mentioned above.

$latex b_0$ $latex b_1$ $latex b_2$ $latex b_3$ $latex b_4$ $latex b_5$ $latex b_6$ $latex b_7$ $latex b_8$ $latex b_9$ $latex b_{10}$
0. order 1
1. order 1.5 -0.5
2. order 1.9167 -1.3333 0.4167
3. order 2.2917 -2.4583 1.5417 -0.375
4. order 2.6403 -3.8528 3.6333 -1.7694 0.3486
5. order 2.9701 -5.5021 6.9319 -5.0681 1.9979 -0.3299
6. order 3.2857 -7.3956 11.6658 -11.3799 6.7318 -2.2234 0.3156
7. order 3.59 -9.5252 18.0545 -22.0278 17.3797 -8.6121 2.4452 -0.3042
8. order 3.8848 -11.8842 26.3108 -38.5404 38.0204 -25.1247 10.7015 -2.6632 0.2949
9. order 4.1718 -14.4669 36.642 -62.6463 74.1793 -61.2836 34.8074 -12.9943 2.8776 -0.287
10. order 4.452 -17.2688 49.2505 -96.2691 133.0191 -131.8914 93.6472 -46.617 15.4862 -3.0889 0.2802

The table can be downloaded as an excel file: AB-Coefficients.