Hyperbolische partielle Differentialgleichung

12.5. Hyperbolische partielle Differentialgleichung#

In der Struktur Dynamik und auch anderen Anwendungen trifft man oft auf Systeme von der semidiskreten Form

(12.10)#\[M \ddot{d} + C \dot{d} + K d = F,\]

wobei \(M\) die Massenmatrix, \(C\) eine viskose Dämpfung, \(K\) die Steifigkeitsmatrix und \(F\) ein Vektor gegeben durch die angewandten Kräfte ist. Wir nehmen an, dass \(M, C, K\) symmetrisch sind. Die Massenmatrix \(M\) sei positiv definit und \(C, K\) positiv-semidefinit. Für die vollständige Beschreibung des Anfangswertproblems (12.10) benötigen wir Anfangswerte für

\[d(0) = d_0\quad \text{und}\quad \dot{d}(0) = v_0.\]

12.5.1. Newmark Verfahren#

Die wohl am meisten benutzten Verfahren gehöhren in die Familie der Newmark Verfahren:

(12.11)#\[\begin{split}\begin{split} M\,a_{n+1} + C\,v_{n+1} + K\,d_{n+1} & = F_{n+1}\\ d_{n+1} & = d_n + \Delta t\,v_n + \frac{\Delta t^2}{2} ((1-2\beta) a_n + 2\beta a_{n+1})\\ v_{n+1} & = v_n + \Delta t\,((1-\gamma)\,a_n + \gamma\,a_{n+1}), \end{split}\end{split}\]

wobei \(d_n, v_n, a_n\) Approximationen von \(d(t_n), \dot{d}(t_n), \ddot{d}(t_n)\) sind. Die Parameter \(\beta\) und \(\gamma\) bestimmen die Stabilitäts- und Genauigkeits-Charakteristik.

Es gibt verschiedene Formen für die Implementierung. Oft wird die sogenannte \(a\)-Form benutzt, welche auf der Beschleunigung \(a\) aufbaut. Analoges mit der Deformation \(d\) führt zu äquivalenten Implementationen. Wir betrachten hier die a-Form und definieren folgende Predictors:

(12.12)#\[\begin{split}\begin{split} \tilde{d}_{n+1} & = d_n + \Delta t\,v_n + \frac{\Delta t^2}{2} (1-2\beta) a_n\\ \tilde{v}_{n+1} & = v_n + (1-\gamma) \Delta t\, a_n \end{split}\end{split}\]

und für \(d_{n+1}\) und \(v_{n+1}\) aus (12.11) folgend die Correctors:

(12.13)#\[\begin{split}\begin{split} d_{n+1} & = \tilde{d}_{n+1} + \beta\,\Delta t^2\, a_{n+1}\\ v_{n+1} & = \tilde{v}_{n+1} + \gamma\, \Delta t\, a_{n+1}.\end{split}\end{split}\]

Die initiale Beschleunigung \(a_0\) kann aus

\[M a_0 = F - C v_0 - K d_0\]

berechnet werden. Die noch unbekannte Beschleunigung zum neuen Zeitschritt \(a_{n+1}\) ist nun gegeben durch die Lösung des Gleichungssystems

\[(M + \gamma\, \Delta t\, C + \beta\, \Delta t^2\, K)\,a_{n+1} = F_{n+1} - C\,\tilde{v}_{n+1} - K\,\tilde{d}_{n+1}.\]

Mit den Gleichungen (12.13) folgt \(d_{n+1}\) und \(v_{n+1}\), die Schätzwerte \(\tilde{d}_{n+1}, \tilde{v}_{n+1}\) werden korrigiert. Entsprechend nennt man die Gleichungen (12.13) Correctors.

Remark 12.3

Für \(\beta = \frac{1}{4}\) und \(\gamma = \frac{1}{2}\) wird quasi die mittlere Beschleunigung benutzt. Wir erhalten in dem Fall die Trapezmethode, welche bedingungslos stabil ist.

12.5.2. Anwendung#

Wir betrachten als Anwendung die Wellengleichung

\[\begin{split}\begin{split}\ddot{u}(t,x) - \Delta u(t,x) & = f(t)\quad\text{für}\ x\in\Omega = [0,1]^2\\ u(t, x) & = 0\quad \text{für}\ x\in\partial \Omega\\ u(0) & = 0\quad \text{für}\ x\in \Omega\\ \dot{u}(0) & = 0\quad \text{für}\ x\in \Omega\\ \end{split}\end{split}\]
Hide code cell content 
from netgen.geom2d import unit_square
from ngsolve import *
from ngsolve.webgui import Draw

import numpy as np
import matplotlib.pyplot as plt

from netgen.geom2d import CSG2d, Circle
circle = Circle((0,0),2,bc='wall')
geo = CSG2d()
geo.Add(circle)

#mesh = Mesh(unit_square.GenerateMesh(maxh=0.05))
mesh = Mesh(geo.GenerateMesh(maxh=0.3))# etwas grob 0.15
V = H1(mesh,order = 4, dirichlet='wall')#'bottom|right|top|left')
u = V.TrialFunction()
v = V.TestFunction()
gfu = GridFunction(V)
gfv = GridFunction(V)# 1. Zeitableitung
gfa = GridFunction(V)# 2. Zeitableitung

V.ndof

2305

K = BilinearForm(V)
K += 0.1*grad(u)*grad(v)*dx
K.Assemble()

M = BilinearForm(V)
M += u*v*dx
M.Assemble();

Als Quelle benutzen wir einen Puls von der Zeitlänge \(\Delta t\) mit einer beschränkten örtlichen Ausdehnung:

\[f(t) = e^{-((x-x_0)^2+(y-y_0)^2)/\tau^2}\]
x0=0
y0=0
f = 5000*exp(-((x-x0)**2+(y-y0)**2)/0.01)

F = LinearForm(V)
F += f*v*dx
F.Assemble();

Draw(f,mesh,'f')

BaseWebGuiScene

Sei \(M^* = M + \beta\, \Delta t^2\, K\). Wir wählen die Trapezmethode:

beta = 1/4
gamma = 1/2

dt = 1/125.0
mstar = K.mat.CreateMatrix()
mstar.AsVector().data = M.mat.AsVector() + beta*dt**2 * K.mat.AsVector()

Wir berechnen nun die zeitabhängige Lösung.

scene = Draw(gfu,mesh,'u', min=-2,max=2,autoscale=False)

freq=1/(10*dt)*2

t = 0
gfu.vec[:] = 0
gfv.vec[:] = 0
gfa.vec[:] = 0
T = 10

# predictors
utilde = gfu.vec.CreateVector()
vtilde = gfu.vec.CreateVector()
b = gfu.vec.CreateVector()
mip = mesh(0,0)
ui = [gfu(mip)]
ti = [t]
with TaskManager():
    while t < T:
        #F.Assemble() # muss bei zeitabhängiger Linearform aufgerufen werden.
        # predictors berechnne
        utilde.data = gfu.vec + dt * gfv.vec + dt**2/2*(1-2*beta)*gfa.vec
        vtilde.data = gfv.vec + (1-gamma)*dt*gfa.vec
        # a neu berechnen
        if t < 1/(2*freq): # ein Impuls
            b.data = F.vec
            b.data *= np.sin(2*np.pi*freq*t)**2
            b.data -= K.mat*utilde
        else:
            b.data = - K.mat*utilde
        gfa.vec.data = mstar.Inverse(freedofs=V.FreeDofs(),inverse='sparsecholesky')*b
        # corrector berechnen
        gfu.vec.data = utilde + beta*dt**2*gfa.vec
        gfv.vec.data = vtilde + gamma*dt*gfa.vec
        # Redraw
        scene.Redraw()
        ui.append(gfu(mip))
        ti.append(t)
        t+=dt
Inhalt