Verschiebungsmethode

Grundidee und physikalisches Prinzip

Die Verschiebungsmethode ist die Standardformulierung der Finite-Elemente-Methode (FEM) in der Festkörpermechanik. Ihre primären Unbekannten sind die Verschiebungen der Körperpunkte. Diese beschreiben Translation, Rotation und mögliche Verformung eines Festkörpers. Die Methode eignet sich sowohl für lineare Aufgaben wie Eigenschwingungen als auch für stark nichtlineare Vorgänge wie Crashtests. Üblich sind isoparametrische Elemente, bei denen Geometrie und Verschiebungen mit denselben Formfunktionen angenähert werden, sowie die Galerkin-Methode.

Grundlage ist das Prinzip von d’Alembert in Lagrange’scher Fassung: ∫_V δu⃗·ρ₀ü⃗ dV + ∫V T̃:δE dV = ∫{A^σ} δu⃗·t⃗₀ dA + ∫_V δu⃗·ρ₀k⃗₀ dV für alle zulässigen virtuellen Verschiebungen δu⃗.

Virtuelle Verschiebungen sind gedachte, mit den Randbedingungen verträgliche Verschiebungen. Der erste Term beschreibt ihre virtuelle Arbeit an der Impulsänderung ρ₀ü⃗, der zweite die virtuelle Deformationsarbeit der Spannungen T̃ an den virtuellen Verzerrungen δE. Rechts steht die Arbeit äußerer Oberflächen- und Volumenkräfte. Auf Bereichen mit vorgeschriebenen Verschiebungen müssen die virtuellen Verschiebungen verschwinden. Ist die Gleichung für alle zulässigen virtuellen Verschiebungen erfüllt, stimmen Verschiebungen, Verzerrungen und Spannungen mit der Impulsbilanz überein; die vorausgesetzte Symmetrie des Spannungstensors erfüllt zusätzlich die Drehimpulsbilanz.

Diskretisierung und globales Gleichungssystem

Der Körper wird lückenlos und überschneidungsfrei in kompatible finite Elemente zerlegt. An ihren Knoten werden globale Koordinaten und Knotenverschiebungen gespeichert. Innerhalb eines Elements werden die Verschiebungen durch u⃗ = N û angenähert. N ist eine ortsabhängige 3×m-Formfunktionsmatrix, û ein zeitabhängiger Vektor aus m Knotenverschiebungen. Entsprechend gelten ü⃗ = N ü̂ und in der Galerkin-Methode δu⃗ = N δû. Die Beschränkung auf endlich viele Formfunktionen verursacht den Diskretisierungsfehler.

Spannungen und Verzerrungen werden in Voigt’scher Notation als Vektoren geschrieben: T̂ = (T̃xx, T̃yy, T̃zz, T̃xy, T̃yz, T̃xz)ᵀ, δÊ = (δExx, δEyy, δEzz, 2δExy, 2δEyz, 2δExz)ᵀ. Die Faktoren 2 sorgen dafür, dass T̃:δE = δÊᵀT̂ gilt; die doppelten Schubverzerrungen entsprechen den Gleitungen. Mit der Verzerrungsverschiebungsmatrix B gilt δÊ = Bδû.

Daraus entsteht für ein Element δûᵀ[Mü̂ + r̂ − f̂] = 0 für alle zulässigen δû. Dabei sind M = ∫_V Nᵀρ₀N dV die konstante Massenmatrix, r̂ = ∫V BᵀT̂ dV der Vektor der inneren Knotenreaktionen und f̂ = ∫{A^σ} Nᵀt⃗₀ dA + ∫_V Nᵀρ₀k⃗₀ dV der äußere Knotenkraftvektor. Durch Assemblierung, also das Aufsummieren zusammengehöriger Elementbeiträge, entsteht das globale System.

Für Verschiebungsrandbedingungen wird û in unbekannte Komponenten ûu und vorgegebene Komponenten ûb partitioniert. An vorgegebenen Verschiebungen verschwinden die zugehörigen virtuellen Verschiebungen. Für den unbekannten Teil folgt Muu ü̂u = f̂u − r̂u − Mub ü̂b. Nach der Lösung können aus der ausgeschiedenen unteren Gleichungszeile Lagerreaktionen bestimmt werden. Der weitere Text nimmt statische Festlager mit verschwindenden Verschiebungen und Beschleunigungen an und vernachlässigt Geschwindigkeitsabhängigkeiten. Dann lautet die zentrale Bewegungsgleichung einfach Mü̂ = f̂ − r̂.

Zeitintegration und lineare Berechnung

Obwohl in der Bewegungsgleichung Beschleunigungen auftreten, bleiben die Verschiebungen die primären Unbekannten. Beide Größen sind über die zweite Zeitableitung verbunden. Verbreitet sind Einschrittverfahren, die Werte zum Zeitpunkt tⁿ⁺¹ aus bekannten Werten bei tⁿ bestimmen.

Beim impliziten Newmark-beta-Verfahren gilt ü̂ⁿ⁺¹ = p ûⁿ⁺¹ + p̂ⁿ, wobei p aus dem Integrationsalgorithmus stammt und p̂ⁿ bekannt ist. Bei expliziter Zeitintegration gilt dagegen ûⁿ⁺¹ = q ü̂ⁿ + q̂ⁿ, mit dem Algorithmusparameter q und dem bekannten Vektor q̂ⁿ.

Im linearen Fall sind Verzerrungen und Spannungen linear von den Knotenverschiebungen abhängig: ÊL = BL û, T̂ = T̂⁰ + CL BL û. T̂⁰ enthält Eigenspannungen, CL ist die von den Verzerrungen unabhängige Stoffmatrix. Damit wird r̂ = r̂L⁰ + KL û, mit r̂L⁰ = ∫_V BLᵀT̂⁰ dV und der konstanten linearen Steifigkeitsmatrix KL = ∫_V BLᵀCLBL dV. Die Bewegungsgleichung lautet Mü̂ⁿ⁺¹ + KLûⁿ⁺¹ = f̂ⁿ⁺¹ − r̂L⁰. Mit dem Newmark-beta-Zusammenhang ergibt sich (pM + KL)ûⁿ⁺¹ = f̂ⁿ⁺¹ − r̂L⁰ − Mp̂ⁿ. Für ein statisches Gleichgewicht mit ü̂ = 0 bleibt KLûⁿ⁺¹ = f̂ⁿ⁺¹ − r̂L⁰. Lineare Systeme können außerdem in den Modalraum übertragen werden, etwa zur Untersuchung freier Schwingungen.

Nichtlineare Lösungsverfahren

Nichtlinearitäten können aus dem Materialverhalten, etwa Plastizität oder nichtlinearer Elastizität, aus verschiebungsabhängigen Randbedingungen wie Kontakt und verformungsabhängigen Kräften oder aus großen Drehungen und Verformungen, Knicken und Beulen entstehen.

Bei der impliziten Lösung wird gewöhnlich das Newton-Verfahren verwendet. Die innere Reaktion wird an einer Näherung ûⁱ linearisiert: r̂(ûⁱ + Δû) ≈ r̂ⁱ + (Gⁱ + Kⁱ)Δû. Gⁱ ist die geometrische Steifigkeitsmatrix. Sie ist nur bei geometrischer Nichtlinearität erforderlich, weil nur dann B von den Verschiebungen abhängt. Die Materialsteifigkeit ist Kⁱ = ∫_V BⁱᵀCBⁱ dV, wobei C = dT̂/dÊ der konsistente Tangentenoperator ist. Bei verschiebungsabhängigen äußeren Kräften gilt zusätzlich f̂(ûⁱ + Δû) ≈ f̂ⁱ + FⁱΔû. Mit KNLⁱ = Gⁱ + Kⁱ − Fⁱ und Δü̂ = pΔû entsteht dynamisch (pM + KNLⁱ)Δû = f̂ⁱ − r̂ⁱ − Mü̂ⁱ. Im statischen Fall gilt KNLⁱΔû = f̂ⁱ − r̂ⁱ.

Die Iteration startet mit einer Näherung für ûⁿ⁺¹. Daraus werden Matrix und Residuum aufgebaut und Δû berechnet. Sind geeignete Normen von Δû und Residuum kleiner als vorgegebene Grenzen, wird die Lösung akzeptiert und das Materialmodell aktualisiert. Andernfalls setzt man ûⁿ⁺¹,ⁱ⁺¹ = ûⁿ⁺¹,ⁱ + Δû und wiederholt die Rechnung.

Bei expliziter Zeitintegration ist keine Newton-Linearisierung nötig: Mü̂ⁿ⁺¹ = f̂ⁿ − r̂ⁿ. Eine diagonalisierte Massenmatrix, die „lumped mass matrix“, erlaubt die schnelle Berechnung ü̂ⁿ⁺¹ = M⁻¹(f̂ⁿ − r̂ⁿ). Das Verfahren ist jedoch nur unterhalb einer kritischen Zeitschrittweite Δtc stabil. Nach der Courant-Friedrichs-Lewy-Bedingung darf sich ein Signal innerhalb eines Schrittes höchstens ungefähr bis zum nächsten Knoten ausbreiten. Für lc = 10 mm, Stahl mit E = 200.000 MPa und ρ = 7870 kg/m³ gilt Δtc ≤ lc/√(E/ρ) ≈ 2·10⁻⁶ s. Daher sind für Bewegungen von Zehntelsekunden oft zehntausende Schritte nötig. Weil Nichtlinearitäten ohne Linearisierung behandelt werden und der Aufwand für die Beschleunigungsberechnung nur linear mit der Zahl der Unbekannten wächst, eignet sich das Verfahren besonders für große, kurzzeitige dynamische Probleme wie Crashtests und auch für sehr große quasistatische Aufgaben.

Elementmatrizen und Verzerrungen

Da die Elementintegrale meist nicht exakt lösbar sind, werden sie beispielsweise durch Gauß-Quadratur angenähert: Das Integral wird als gewichtete Summe der Integranden an Integrationspunkten berechnet.

Bei einem dreidimensionalen Element interpolieren k Formfunktionen Nⁱ die globalen Koordinaten aus lokalen Koordinaten Θ = (ξ, η, ζ)ᵀ ∈ [−1,1]³: X⃗ = Σᵢ Nⁱ(Θ)X⃗ⁱ = N(Θ)X̂. Die Jacobi-Matrix J vermittelt zwischen lokalen und globalen Ableitungen. Aus den analytisch berechenbaren Ableitungen nach ξ, η und ζ folgen über Jᵀ⁻¹ die Ableitungen nach X, Y und Z. Für das Volumenelement gilt dV = det(J) dξ dη dζ.

In isoparametrischen Elementen werden Verschiebungen ebenso interpoliert: u⃗ = Nû. Aus ihren globalen Ableitungen entsteht der Verschiebungsgradient H = GRAD(u⃗). Der Green-Lagrange-Verzerrungstensor lautet E = ½(H + Hᵀ + HᵀH). Seine sechs unabhängigen Komponenten bilden Ê = (Exx, Eyy, Ezz, 2Exy, 2Eyz, 2Exz)ᵀ. Die Matrix B = dÊ/dû ∈ ℝ^(6×3k) verknüpft virtuelle Verzerrungen und virtuelle Knotenverschiebungen durch δÊ = Bδû.

Im geometrisch linearen Fall wird der quadratische Anteil vernachlässigt: EL = ½(H + Hᵀ), ÊL = BLû. BL hängt nur von den Ableitungen der Formfunktionen ab. Im geometrisch nichtlinearen Fall kommt ENL = ½HᵀH hinzu; dadurch gilt B = BL + BNL, und BNL sowie B hängen von den aktuellen Knotenverschiebungen ab.

Eine Materialroutine berechnet die aktuellen Spannungen aus den Spannungen des letzten Zeitschritts, dem Verzerrungsinkrement und möglichen inneren Variablen. Für das Newton-Verfahren liefert ihre Ableitung den konsistenten Tangentenoperator C = dT̂ⁱ/dÊ, eine 6×6-Matrix. „Konsistent“ bedeutet, dass C aus der Ableitung der numerischen Materialroutine und nicht nur aus dem zugrunde liegenden analytischen Materialmodell stammt. Bei linearer Elastizität ist C = CL. Die geometrische Steifigkeitsmatrix G entsteht aus der verschiebungsbedingten Änderung von BᵀT̂ bei konstant gehaltener Spannung und besitzt eine aus den aktuellen Spannungen und Formfunktionsableitungen aufgebaute Blockstruktur.

Beispiel eines Zugstabs

Als Beispiel dient ein einseitig eingespannter Zugstab mit linear abnehmender Querschnittsfläche A = a − bX und einer Einzelkraft F am freien Ende. Ein eindimensionales zweiknotiges Stabelement mit ξ ∈ [−1,1] verwendet die linearen Formfunktionen N = ½(1−ξ, 1+ξ). Mit der Elementlänge L gilt J = L/2 und N,X = (−1, 1)/L. Damit sind im geometrisch linearen Bereich die konstante Dehnung EL = BLÛ und BL = (−1, 1)/L. Bei T = CEL ergibt die Integration über den veränderlichen Querschnitt die Elementsteifigkeitsmatrix KL = C·Am/L · ((1, −1), (−1, 1)), wobei Am die Querschnittsfläche in der Elementmitte ist. Die Endlast wird als f̂ = (0, F)ᵀ angesetzt.

Zwei Elemente halber Länge werden am gemeinsamen Knoten assembliert. Mit U¹ = 0 und der Kraft F am dritten Knoten liefert das im Artikel angegebene Zahlenbeispiel U² = 0,032 und U³ = 0,109. Verwendet werden L = 1000 mm, a = 100 mm², b = 0,09 mm, C = 200.000 MPa und F = 1000 N.

Zum Vergleich besitzt das Problem die analytische Lösung u(X) = F/(Cb) · ln(a/(a−bX)), mit u(0) = 0. Bei zunehmender Netzverfeinerung konvergiert die FEM-Lösung gegen einen Grenzwert. Das lineare Stabelement bildet innerhalb jedes Elements nur konstante Dehnung und Spannung ab, während die tatsächliche Spannung wegen der abnehmenden Querschnittsfläche kontinuierlich steigt. Kleinere Elemente nähern diesen glatten Verlauf durch feinere Stufen immer genauer an.