Newmark-beta-Verfahren

Newmark-beta-Verfahren sind Methoden zur impliziten numerischen Integration von Differentialgleichungen. Die Verfahren gehören zu den Einschrittverfahren, da zur Berechnung der Werte zur Zeit t n + 1 {\displaystyle t_{n+1}} nur die Werte des vorangegangenen Zeitschritts zur Zeit t n {\displaystyle t_{n}} benötigt werden. Dabei werden zwei Parameter γ {\displaystyle \gamma } und β {\displaystyle \beta } eingeführt, mit denen die Stabilität und die Genauigkeit des Verfahrens gesteuert werden. Die Verfahrensklasse ist in der numerischen Analyse der Dynamik von Festkörpern wie in der Finite-Elemente-Methode weit verbreitet. Benannt ist sie nach Nathan M. Newmark, der sie 1959 für die Anwendung in der Strukturdynamik entwickelte.[1]

Herleitung

Annahme linearer oder konstanter Beschleunigung

Im Zeitintervall [ t 0 , t e ] {\displaystyle [t_{0,}t_{e}]} , in dem eine Lösung x ( t ) {\displaystyle x(t)} einer Differentialgleichung zweiter Ordnung in der Zeit gesucht wird, sei eine streng monoton steigende Folge von Zeitpunkten { t 0 , t 1 , , t e } {\displaystyle \{t_{0},t_{1},\dotsc ,t_{e}\}} vorgegeben, zu denen die Lösung x ( t ) {\displaystyle x(t)} berechnet werden soll. Der Wert der Variable x n {\displaystyle x_{n}} , ihre Rate x ˙ n {\displaystyle {\dot {x}}_{n}} und Beschleunigung x ¨ n {\displaystyle {\ddot {x}}_{n}} seien zur Zeit t n {\displaystyle t_{n}} bekannt. Die Beschleunigung wird im Intervall t [ t n , t n + 1 ] {\displaystyle t\in [t_{n},t_{n+1}]} linear interpoliert, siehe Bild:

(I)     x ¨ h ( t ) = t n + 1 t t n + 1 t n x ¨ n + t t n t n + 1 t n x ¨ n + 1 {\displaystyle {\ddot {x}}^{h}(t)={\frac {t_{n+1}-t}{t_{n+1}-t_{n}}}{\ddot {x}}_{n}+{\frac {t-t_{n}}{t_{n+1}-t_{n}}}{\ddot {x}}_{n+1}}

worin x h {\displaystyle x^{h}} eine Näherungslösung der gesuchten Funktion x {\displaystyle x} bezeichnet. Integration über die Zeit liefert mit Δ t = t t n {\displaystyle \Delta t=t-t_{n}} :

(II)     x ˙ h ( t ) = x ˙ n + t n t x ¨ h ( τ ) d τ = x ˙ n + Δ t [ ( 1 γ ) x ¨ n + γ x ¨ h ( t ) ] {\displaystyle {\dot {x}}^{h}(t)={\dot {x}}_{n}+\int _{t_{n}}^{t}{\ddot {x}}^{h}(\tau )\mathrm {d} \tau ={\dot {x}}_{n}+\Delta t[(1-\gamma ){\ddot {x}}_{n}+\gamma {\ddot {x}}^{h}(t)]}

(III)     x h ( t ) = x n + t n t x ˙ h ( τ ) d τ = x n + Δ t x ˙ n + Δ t 2 [ ( 1 2 β ) x ¨ n + β x ¨ h ( t ) ] {\displaystyle x^{h}(t)=x_{n}+\int _{t_{n}}^{t}{\dot {x}}^{h}(\tau )\mathrm {d} \tau =x_{n}+\Delta t{\dot {x}}_{n}+\Delta t^{2}\left[\left({\frac {1}{2}}-\beta \right){\ddot {x}}_{n}+\beta {\ddot {x}}^{h}(t)\right]}

Mit

γ = 1 2 {\displaystyle \gamma ={\frac {1}{2}}}    und    β = 1 6 {\displaystyle \beta ={\frac {1}{6}}}

sind diese Formeln für lineare Systeme exakt und liefern das lineare Beschleunigungsverfahren. Die von Newmark ursprünglich angegebenen Werte

γ = 1 2 {\displaystyle \gamma ={\frac {1}{2}}}    und    β = 1 4 {\displaystyle \beta ={\frac {1}{4}}}

entsprechen dem #Konstante Durchschnittsbeschleunigungsverfahren mit

x ¨ h ( t ) = 1 2 ( x ¨ n + x ¨ n + 1 ) {\displaystyle {\ddot {x}}^{h}(t)={\frac {1}{2}}({\ddot {x}}_{n}+{\ddot {x}}_{n+1})} .

Unter der Voraussetzung, dass die Extremwerte der Beschleunigung im Intervall [ t n , t n + 1 ] {\displaystyle [t_{n},t_{n+1}]} an den Grenzen des Intervalls auftreten, stellen die Integrale in Gleichungen (II) und (III) eine abgebrochene Taylorreihe mit Restglied dar, wobei mit 0 < γ < 1 {\displaystyle 0<\gamma <1} und 0 < β < 1 {\displaystyle 0<\beta <1} andere Approximationen x h {\displaystyle x^{h}} gefunden werden. So können auch andere Werte für die Konstanten γ {\displaystyle \gamma } und β {\displaystyle \beta } motiviert werden.

Start der Berechnung

Der Newmark-Algorithmus startet zur Zeit t = t 0 {\displaystyle t=t_{0}} mit n = 0 {\displaystyle n=0} . Zumeist wird angenommen, dass für t < t 0 {\displaystyle t<t_{0}} die Beschleunigungen verschwinden. Mit dieser Annahme ist der Algorithmus unter Vorgabe der Anfangswerte x 0 {\displaystyle x_{0}} und Anfangsgeschwindigkeit x ˙ 0 {\displaystyle {\dot {x}}_{0}} selbststartend, d. h. die Anfangsbeschleunigungen brauchen nicht in einem ersten Schritt berechnet zu werden.

Aktualisierung der Variablen

Mit dem Newmark-Algorithmus werden aus gegebenen Werten x n , x ˙ n {\displaystyle x_{n},{\dot {x}}_{n}} und x ¨ n {\displaystyle {\ddot {x}}_{n}} zur Zeit t n {\displaystyle t_{n}} die entsprechenden Werte zur Zeit t n + 1 {\displaystyle t_{n+1}} berechnet. Die im Intervall [ t n , t n + 1 ] {\displaystyle [t_{n},t_{n+1}]} liegenden Werte können mit den Gleichungen (I) bis (III) interpoliert werden. Mit t = t n + 1 {\displaystyle t=t_{n+1}} und x ¨ h ( t n + 1 ) = x ¨ n + 1 {\displaystyle {\ddot {x}}^{h}(t_{n+1})={\ddot {x}}_{n+1}} bekommt man aus Gleichungen (II) und (III):

(IV)     x ˙ h ( t n + 1 ) = x ˙ n + 1 = x ˙ n + Δ t [ ( 1 γ ) x ¨ n + γ x ¨ n + 1 ] {\displaystyle {\dot {x}}^{h}(t_{n+1})={\dot {x}}_{n+1}={\dot {x}}_{n}+\Delta t[(1-\gamma ){\ddot {x}}_{n}+\gamma {\ddot {x}}_{n+1}]} ,

(V)     x h ( t n + 1 ) = x n + 1 = x n + Δ t x ˙ n + Δ t 2 [ ( 1 2 β ) x ¨ n + β x ¨ n + 1 ] {\displaystyle x^{h}(t_{n+1})=x_{n+1}=x_{n}+\Delta t{\dot {x}}_{n}+\Delta t^{2}\left[\left({\frac {1}{2}}-\beta \right){\ddot {x}}_{n}+\beta {\ddot {x}}_{n+1}\right]} .

Die beiden Gleichungen (IV) und (V) enthalten drei Unbekannte x n + 1 , x ˙ n + 1 {\displaystyle x_{n+1},{\dot {x}}_{n+1}} und x ¨ n + 1 {\displaystyle {\ddot {x}}_{n+1}} . Die dritte zum Abschluss benötigte Gleichung liefert die zu lösende Differentialgleichung.

Bei β 0 {\displaystyle \beta \neq 0} kann auch x n + 1 {\displaystyle x_{n+1}} als primäre Unbekannte gewählt werden:

x ¨ n + 1 = 1 β Δ t 2 ( x n + 1 x n ) 1 β Δ t x ˙ n 1 2 β β x ¨ n {\displaystyle {\ddot {x}}_{n+1}={\frac {1}{\beta \Delta t^{2}}}(x_{n+1}-x_{n})-{\frac {1}{\beta \Delta t}}{\dot {x}}_{n}-{\frac {{\frac {1}{2}}-\beta }{\beta }}{\ddot {x}}_{n}}
x ˙ n + 1 = γ β Δ t ( x n + 1 x n ) + ( 1 γ β ) x ˙ n + Δ t β 1 2 γ β x ¨ n {\displaystyle {\dot {x}}_{n+1}={\frac {\gamma }{\beta \Delta t}}(x_{n+1}-x_{n})+\left(1-{\frac {\gamma }{\beta }}\right){\dot {x}}_{n}+\Delta t{\frac {\beta -{\frac {1}{2}}\gamma }{\beta }}{\ddot {x}}_{n}} .

Sind einmal die Werte x n + 1 , x ˙ n + 1 {\displaystyle x_{n+1},{\dot {x}}_{n+1}} und x ¨ n + 1 {\displaystyle {\ddot {x}}_{n+1}} berechnet, wird der Zähler n {\displaystyle n} inkrementiert und die Berechnung fortgesetzt, bis das Ende des interessierenden Zeitintervalls erreicht ist.

Spezialfälle

Konstante Durchschnittsbeschleunigungsverfahren

Die ursprüngliche Form des Newmark-Verfahrens entspricht einer konstanten mittleren Beschleunigung

x ¨ h ( t ) = 1 2 ( x ¨ n + 1 + x ¨ n ) {\displaystyle {\ddot {x}}^{h}(t)={\frac {1}{2}}({\ddot {x}}_{n+1}+{\ddot {x}}_{n})}

mit der man in den obigen Formeln (IV) und (V)

γ = 1 2 {\displaystyle \gamma ={\frac {1}{2}}} und β = 1 4 {\displaystyle \beta ={\frac {1}{4}}}

bekommt.

Gleichung Folgerung
x ˙ h ( t ) = x ˙ n + t n t x ¨ h ( τ ) d τ = x ˙ n + Δ t 2 ( x ¨ n + x ¨ n + 1 ) {\displaystyle {\dot {x}}^{h}(t)={\dot {x}}_{n}+\int _{t_{n}}^{t}{\ddot {x}}^{h}(\tau )\mathrm {d} \tau ={\dot {x}}_{n}+{\frac {\Delta t}{2}}({\ddot {x}}_{n}+{\ddot {x}}_{n+1})} γ = 1 2 {\displaystyle \gamma ={\frac {1}{2}}}
x n + t n t x ˙ h d τ = x n + t n t ( x ˙ n + τ t n 2 ( x ¨ n + x ¨ n + 1 ) ) d τ = x n + Δ t x ˙ n + Δ t 2 4 ( x ¨ n + x ¨ n + 1 ) {\displaystyle x_{n}+\int _{t_{n}}^{t}{\dot {x}}^{h}\mathrm {d} \tau =x_{n}+\int _{t_{n}}^{t}\left({\dot {x}}_{n}+{\frac {\tau -t_{n}}{2}}({\ddot {x}}_{n}+{\ddot {x}}_{n+1})\right)\mathrm {d} \tau =x_{n}+\Delta t{\dot {x}}_{n}+{\frac {\Delta t^{2}}{4}}({\ddot {x}}_{n}+{\ddot {x}}_{n+1})} β = 1 4 {\displaystyle \beta ={\frac {1}{4}}}

Zentrale Differenzenquotienten

Hauptartikel: Differenzenquotient

Die zentralen Differenzenquotienten

(VI)    x ˙ n = 1 2 Δ t ( x n + 1 x n 1 ) {\displaystyle {\dot {x}}_{n}={\frac {1}{2\Delta t}}(x_{n+1}-x_{n-1})}

(VII)     x ¨ n = 1 Δ t 2 ( x n + 1 2 x n + x n 1 ) {\displaystyle {\ddot {x}}_{n}={\frac {1}{\Delta t^{2}}}(x_{n+1}-2x_{n}+x_{n-1})}

entsprechen den obigen Formeln (IV) und (V) mit

γ = 1 2 {\displaystyle \gamma ={\frac {1}{2}}} und β = 0 {\displaystyle \beta =0} .
Gleichung Folgerung
x ¨ n = 1 Δ t 2 ( x n + 1 2 x n + x n 1 ) x n 1 = Δ t 2 x ¨ n x n + 1 + 2 x n x ˙ n = 1 2 Δ t ( x n + 1 x n 1 ) Δ t x ˙ n = 1 2 ( x n + 1 x n 1 ) = x n + 1 Δ t 2 2 x ¨ n x n x n + 1 = x n + Δ t x ˙ n + Δ t 2 2 x ¨ n {\displaystyle {\begin{array}{lcl}{\ddot {x}}_{n}&=&{\frac {1}{\Delta t^{2}}}(x_{n+1}-2x_{n}+x_{n-1})\\[1ex]\rightarrow x_{n-1}&=&\Delta t^{2}{\ddot {x}}_{n}-x_{n+1}+2x_{n}\\[1ex]{\dot {x}}_{n}&=&{\frac {1}{2\Delta t}}(x_{n+1}-x_{n-1})\\[1ex]\rightarrow \Delta t{\dot {x}}_{n}&=&{\frac {1}{2}}(x_{n+1}-x_{n-1})=x_{n+1}-{\frac {\Delta t^{2}}{2}}{\ddot {x}}_{n}-x_{n}\\[1ex]\rightarrow x_{n+1}&=&x_{n}+\Delta t{\dot {x}}_{n}+{\frac {\Delta t^{2}}{2}}{\ddot {x}}_{n}\end{array}}} β = 0 {\displaystyle \beta =0}
Δ t 2 x ¨ n = 1 2 Δ t ( x n + 1 2 x n + x n 1 ) Δ t 2 x ¨ n + 1 = 1 2 Δ t ( x n + 2 2 x n + 1 + x n ) x ˙ n = 1 2 Δ t ( x n + 1 x n 1 ) x ˙ n + 1 = 1 2 Δ t ( x n + 2 x n ) Δ t 2 ( x ¨ n + x ¨ n + 1 ) = 1 2 Δ t ( x n + 2 x n x n + 1 + x n 1 ) = x ˙ n + 1 x ˙ n x ˙ n + 1 = x ˙ n + Δ t 2 ( x ¨ n + x ¨ n + 1 ) {\displaystyle {\begin{array}{lcl}{\frac {\Delta t}{2}}{\ddot {x}}_{n}&=&{\frac {1}{2\Delta t}}(x_{n+1}-2x_{n}+x_{n-1})\\[1ex]{\frac {\Delta t}{2}}{\ddot {x}}_{n+1}&=&{\frac {1}{2\Delta t}}(x_{n+2}-2x_{n+1}+x_{n})\\[2ex]{\dot {x}}_{n}&=&{\frac {1}{2\Delta t}}(x_{n+1}-x_{n-1})\\[1ex]{\dot {x}}_{n+1}&=&{\frac {1}{2\Delta t}}(x_{n+2}-x_{n})\\[1ex]\rightarrow {\frac {\Delta t}{2}}({\ddot {x}}_{n}+{\ddot {x}}_{n+1})&=&{\frac {1}{2\Delta t}}(x_{n+2}-x_{n}-x_{n+1}+x_{n-1})={\dot {x}}_{n+1}-{\dot {x}}_{n}\\[1ex]\rightarrow {\dot {x}}_{n+1}&=&{\dot {x}}_{n}+{\frac {\Delta t}{2}}({\ddot {x}}_{n}+{\ddot {x}}_{n+1})\end{array}}} γ = 1 2 {\displaystyle \gamma ={\frac {1}{2}}}

Explizite Zeitintegration

Hauptartikel: Explizites Euler-Verfahren

Das explizite Zeitintegrationsverfahren gehört nicht zur Familie der (impliziten!) Newmark-beta Algorithmen und wird hier nur zu Vergleichszwecken angegeben. Die obigen Formeln (VI) und (VII) für die zentralen Differenzen sind äquivalent zu

x ˙ n + 1 / 2 = 1 Δ t ( x n + 1 x n ) {\displaystyle {\dot {x}}_{n+1/2}={\frac {1}{\Delta t}}(x_{n+1}-x_{n})}
x ¨ n = 1 Δ t ( x ˙ n + 1 / 2 x ˙ n 1 / 2 ) {\displaystyle {\ddot {x}}_{n}={\frac {1}{\Delta t}}({\dot {x}}_{n+1/2}-{\dot {x}}_{n-1/2})} .

Hier fällt auf, dass die Geschwindigkeiten immer in der Mitte der Zeitintervalle berechnet werden. Mit der Annahme

x ˙ n + 1 x ˙ n + 1 / 2 = x ˙ n 1 / 2 + Δ t x ¨ n x n + 1 = x n + Δ t x ˙ n + 1 / 2 = x n + Δ t ( x ˙ n 1 / 2 + Δ t x ¨ n ) {\displaystyle {\begin{array}{lcl}{\dot {x}}_{n+1}&\approx &{\dot {x}}_{n+1/2}={\dot {x}}_{n-1/2}+\Delta t{\ddot {x}}_{n}\\x_{n+1}&=&x_{n}+\Delta t{\dot {x}}_{n+1/2}=x_{n}+\Delta t({\dot {x}}_{n-1/2}+\Delta t{\ddot {x}}_{n})\end{array}}}

können die Werte x n + 1 {\displaystyle x_{n+1}} und die Geschwindigkeiten x ˙ n + 1 {\displaystyle {\dot {x}}_{n+1}} zum Zeitpunkt t n + 1 {\displaystyle t_{n+1}} auf bereits bekannte Ergebnisse x n , x ˙ n 1 / 2 , x ¨ n {\displaystyle x_{n},{\dot {x}}_{n-1/2},{\ddot {x}}_{n}} zurückgeführt werden und die Differentialgleichung liefert die Bestimmungsgleichung für die nunmehr einzige Unbekannte x ¨ n + 1 {\displaystyle {\ddot {x}}_{n+1}} .

Beispiel

Zeitintegration mit Algorithmen der Newmark Familie

Eine Schwingung gehorche in Abwesenheit einer Erregung der homogenen Differentialgleichung

x ¨ ( t ) + x ( t ) = 0 {\displaystyle {\ddot {x}}(t)+x(t)=0} .

Mit den Anfangsbedingungen

x ( t = 0 ) = x 0 = 0 {\displaystyle x(t=0)=x_{0}=0}
x ˙ ( t = 0 ) = x ˙ 0 = 1 {\displaystyle {\dot {x}}(t=0)={\dot {x}}_{0}=1}

hat die Differentialgleichung die analytische Lösung

x ( t ) = sin ( t ) {\displaystyle x(t)=\sin(t)}

zu der die Anfangsbeschleunigung

x ¨ ( t = 0 ) = sin ( 0 ) = 0 {\displaystyle {\ddot {x}}(t=0)=-\sin(0)=0}

gehört. Die Differentialgleichung liefert die Gleichung für die primäre Unbekannte x ¨ n + 1 {\displaystyle {\ddot {x}}_{n+1}}  :

x ¨ n + 1 + x n + 1 = 0 x ¨ n + 1 = x n + 1 {\displaystyle {\ddot {x}}_{n+1}+x_{n+1}=0\rightarrow {\ddot {x}}_{n+1}=-x_{n+1}}

Die Zeitintegration mit dem Newmark-Verfahren ergibt die Gleichungen für die Werte x n + 1 {\displaystyle x_{n+1}} und Raten x ˙ n + 1 {\displaystyle {\dot {x}}_{n+1}} aus der Tabelle

Parameter Aktualisierungsvorschrift
γ = 1 2 , β = 1 6 {\displaystyle \gamma ={\frac {1}{2}},\;\beta ={\frac {1}{6}}} x n + 1 = x n + Δ t x ˙ n + Δ t 2 6 ( 2 x ¨ n x n + 1 ) x n + 1 = 6 x n + 6 Δ t x ˙ n + 2 Δ t 2 x ¨ n 6 + Δ t 2 x ˙ n + 1 = x ˙ n + Δ t 2 ( x ¨ n x n + 1 ) {\displaystyle {\begin{array}{lcl}x_{n+1}&=&x_{n}+\Delta t{\dot {x}}_{n}+{\frac {\Delta t^{2}}{6}}(2{\ddot {x}}_{n}-x_{n+1})\rightarrow x_{n+1}={\frac {6x_{n}+6\Delta t{\dot {x}}_{n}+2\Delta t^{2}{\ddot {x}}_{n}}{6+\Delta t^{2}}}\\[1ex]{\dot {x}}_{n+1}&=&{\dot {x}}_{n}+{\frac {\Delta t}{2}}({\ddot {x}}_{n}-x_{n+1})\end{array}}}
γ = 1 2 , β = 1 4 {\displaystyle \gamma ={\frac {1}{2}},\;\beta ={\frac {1}{4}}} x n + 1 = x n + Δ t x ˙ n + Δ t 2 4 ( x ¨ n x n + 1 ) x n + 1 = 4 x n + 4 Δ t x ˙ n + Δ t 2 x ¨ n 4 + Δ t 2 x ˙ n + 1 = x ˙ n + Δ t 2 ( x ¨ n x n + 1 ) {\displaystyle {\begin{array}{lcl}x_{n+1}&=&x_{n}+\Delta t{\dot {x}}_{n}+{\frac {\Delta t^{2}}{4}}({\ddot {x}}_{n}-x_{n+1})\rightarrow x_{n+1}={\frac {4x_{n}+4\Delta t{\dot {x}}_{n}+\Delta t^{2}{\ddot {x}}_{n}}{4+\Delta t^{2}}}\\[1ex]{\dot {x}}_{n+1}&=&{\dot {x}}_{n}+{\frac {\Delta t}{2}}({\ddot {x}}_{n}-x_{n+1})\end{array}}}
γ = 1 2 , β = 0 {\displaystyle \gamma ={\frac {1}{2}},\;\beta =0} x n + 1 = x n + Δ t x ˙ n + Δ t 2 2 x ¨ n x ˙ n + 1 = x ˙ n + Δ t 2 ( x ¨ n x n + 1 ) {\displaystyle {\begin{array}{lcl}x_{n+1}&=&x_{n}+\Delta t{\dot {x}}_{n}+{\frac {\Delta t^{2}}{2}}{\ddot {x}}_{n}\\[1ex]{\dot {x}}_{n+1}&=&{\dot {x}}_{n}+{\frac {\Delta t}{2}}({\ddot {x}}_{n}-x_{n+1})\end{array}}}
explizit x ˙ n + 1 = x ˙ n + 1 / 2 = x ˙ n 1 / 2 + Δ t x ¨ n x n + 1 = x n + Δ t x ˙ n 1 / 2 + Δ t 2 x ¨ n {\displaystyle {\begin{array}{lcl}{\dot {x}}_{n+1}&=&{\dot {x}}_{n+1/2}={\dot {x}}_{n-1/2}+\Delta t{\ddot {x}}_{n}\\[1ex]x_{n+1}&=&x_{n}+\Delta t{\dot {x}}_{n-1/2}+\Delta t^{2}{\ddot {x}}_{n}\end{array}}}

Die Lösungen im Intervall t [ 0 , 10 π ] {\displaystyle t\in [0,10\pi ]} und Δ t = 0 , 2 {\displaystyle \Delta t=0,2} haben den Verlauf im Bild. Die mittlere Abweichung

e = n = 1 159 ( x n sin ( t n ) ) 2 159 {\displaystyle e={\sqrt {\frac {\sum _{n=1}^{159}{(x_{n}-\sin(t_{n}))}^{2}}{159}}}}

gibt die Tabelle:

Verfahren Mittlere Abweichung e {\displaystyle e}
Lineares Beschleunigungsverfahren 0.021778594324355638
Zentrales Differenzenverfahren 0.022202937295615111
Konstante mittlere Beschleunigung 0.043283257071468406
Explizites Verfahren 0.022202937295615576

Literatur

  • Robert Gasch, Klaus Knothe, Robert Liebich: Strukturdynamik, Springer Verlag 2012, ISBN 978-3-540-88977-9
  • T. Belytschko, T.J.R. Hughes (Hrsg.): Computational methods for transient analysis. North-Holland 1986. ISBN 9780444864796

Einzelnachweise

  1. Newmark, Nathan M.: A method of computation for structural dynamics. In: Journal of Engineering Mechanics. 85 (EM3). ASCE, 1959, S. 67–94.