37. Parameter Dependent Problems

Many problems from mechanics are of the form

\[ A(u,v) = f(v) \quad \forall \, v \]

with the bilinear-form depending on a small parameter \(\varepsilon\):

\[ A(u,v) = a(u,v) + \frac{1}{\varepsilon} c(T u, T v) \]

where \(a(.,.)\) and \(c(.,.)\) are non-negative bilinear-forms on Hilbert spaces \(V\) and \(Q\), and

\[ T : V \rightarrow Q \]

is a continuous, linear operator.

We introduce a new variable

\[ p := \frac{1}{\varepsilon} T u \]

Using this new variable in the first row, and writing down the definition of \(p\) in weak form we get the system

\[\begin{split} \begin{array}{ccccll} a(u,v) & + & c(p, T v) & = & f(v) & \forall \, v \in V \\ c(Tu, q) & - & \varepsilon c(p,q) & = & 0 & \forall \, q \in Q \end{array} \end{split}\]

We rename the off-diagonal forms as

\[ b(u,q) := c(Tu, q) \]

This structure is called a mixed formulation with penalty term.

In the mixed formulation we can set the parameter \(\varepsilon = 0\).

We are interested in the convergence

\[ p^{(\varepsilon)} = \frac{1}{\varepsilon} T u^{(\varepsilon)} \rightarrow p^{(0)} \]

37.1. Example: Dirichlet boundary condition by penalty

\[ \int \nabla u \nabla v + \frac{1}{\varepsilon} \int_{\Gamma_D} u v = \int f v \]

Here,

\[ \frac{1}{\varepsilon} u^{(\varepsilon)} \rightarrow p^{(0)} = \frac{\partial u}{\partial n} \]

37.2. Example: Nearly incompressible materials

\[ \int 2 \mu \varepsilon (u) : \varepsilon(v) + \int \lambda \operatorname{div} u \operatorname{div} v = \int f v \]

with \(\varepsilon(u) = \tfrac{1}{2} (\nabla u + (\nabla u)^T)\), and \(\lambda \gg \mu\). Introduce \(p := \lambda \operatorname{div} u\)

from ngsolve import *
from ngsolve.webgui import Draw
mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))

fes = VectorH1(mesh,order=2, dirichlet="bottom|top")
u,v = fes.TnT()
mu,lam = 1, 1e5
a = BilinearForm(2*mu*InnerProduct (Sym(Grad(u)), Sym(Grad(v))) * dx \
                 + lam*div(u)*div(v)*dx).Assemble()
f = LinearForm(v[0]*ds("left")).Assemble()

gfu = GridFunction(fes)
gfu.vec.data = a.mat.Inverse(fes.FreeDofs())*f.vec

The deformation looks reasonable, however the pressure looks strange:

Draw (gfu, deformation=True);

Draw (lam*div(gfu), mesh);

Now, we replace the term

\[ \int \lambda \operatorname{div} u \, \operatorname{div} v \]

by introducing an \(L_2\)-projection into elements-wise constants:

\[ \int \lambda \, \Pi_{L_2}^{P^0} (\operatorname{div} u) \, \operatorname{div} v \]

fesP = L2(mesh, order=0)

a = BilinearForm(2*mu*InnerProduct (Sym(Grad(u)), Sym(Grad(v))) * dx \
                 + lam*Interpolate(div(u),fesP)*div(v)*dx).Assemble()
f = LinearForm(v[0]*ds("left")).Assemble()

gfu = GridFunction(fes)
gfu.vec.data = a.mat.Inverse(fes.FreeDofs())*f.vec

Now, also the pressure looks good (as good as it can be in \(P^0\)):

Draw (gfu, deformation=True);
Draw (lam*Interpolate(div(gfu), fesP), mesh);

How to analyze that ?

Instead of solving the parameter dependent problem for \(u\), we define the pressure \(p\):

\[ \int 2 \mu \varepsilon (u) : \varepsilon(v) + \int \underbrace {\lambda \operatorname{div} u}_{p :=} \operatorname{div} v = \int f v \]

Inserting \(p\), and using the defining equation for \(p\) we obtain the saddle-point problem:

\[\begin{split} \begin{array}{ccccll} \int 2 \mu \varepsilon (u) : \varepsilon(v) & + & \int p \operatorname{div} v & = & f(v) & \forall \, v \in V \\ \int \operatorname{div} v \, q & - & \lambda^{-1} \int p q & = & 0 & \forall \, q \in Q \end{array} \end{split}\]

Now, the large parameter \(\lambda\) got a small paramter \(\lambda^{-1}\). The limit problem is a Stokes equation.

Discretize it with a stable method for Stokes leads to the linear system

\[\begin{align*} A u + B^T p &= f \\ B u - C p &= 0 \end{align*}\]

The meaning of the second discrete equation is

\[ p_h = \lambda \Pi_{L_2}^{P^0} \operatorname{div} u_h \]

If we use discontinuous finite elements for the pressure, the matrix \(C\) is block-diagonal, and thus cheaply invertible. This allows to eliminate \(p\) as

\[ p = C^{-1} B u, \]

and insert it into the first equation to obtain

\[ \big( A + B^T C^{-1} B \big) u = f \]

We obtain exactly the same system by inserting the projector into the bilinear-form.