Hello, I need to study problems related to fluid free surfaces based on the level set method.
My ideas are as follows:
we consider the following two PDE problems at the same time:
The unsteady Navier Stokes equations
Find (u,p):[0,T] \to (H_{0,D}^1)^d \times L^2, s.t.
\begin{align}
\int_{\Omega} \partial_t u \cdot v + \int_{\Omega} \nu(c) \nabla u \nabla v + u \cdot \nabla u v - \int_{\Omega} \operatorname{div}(v) p &= \int f v && \forall v \in (H_{0,D}^1)^d, \\
- \int_{\Omega} \operatorname{div}(u) q &= 0 && \forall q \in L^2, \\
\quad u(t=0) & = u_0
\end{align}
with \nu(c) depends on the solution c of the unsteady scalar linear transport problem, which is the level set equation of free surface
Find c: [0,T] \to V_D := \{ c \in L^2(\Omega), u \cdot \nabla c \in L^2(\Omega), c|_{\Gamma_{in}} = c_D\}, s.t.
\int_{\Omega} \partial_t c d + u \cdot \nabla c d = \int_{\Omega} f d \qquad \forall d \in V_0 = \{ c \in L^2(\Omega), u \cdot \nabla c \in L^2(\Omega), c|_{\Gamma_{in}} = 0\}
We use a Taylor-Hood discretization of degree k. The finite element space for the (vectorial) velocity \mathbf{u} and the pressure p are:
\begin{align}
\mathbf{u} \in \mathbf{V} &= \{ v \in H^1(\Omega) | v|_T \in \mathcal{P}^k(T) \}^2 \\
p \in Q &= \{ q \in H^1(\Omega) | q|_T \in \mathcal{P}^{k-1}(T) \}
\end{align}
Stokes problem for initial values
Find \mathbf{u} \in \mathbf{V}, p \in Q so that
\begin{align}
\int_{\Omega} \nu(c) \nabla \mathbf{u} : \nabla \mathbf{v} - \int_{\Omega} \operatorname{div}(\mathbf{v}) p &= \int \mathbf{f} \cdot \mathbf{v} && \forall \mathbf{v} \in \mathbf{V}, \\- \int_{\Omega} \operatorname{div}(\mathbf{u}) q &= 0 && \forall q \in Q,
\end{align}
with boundary conditions:
- the inflow boundary data (‘inlet’) with mean value 1
- “do-nothing” boundary conditions on the ouflow (‘outlet’) and
- homogeneous Dirichlet conditions on all other boundaries (‘wall’).
To solve the Stokes (and later the Navier-Stokes) problem, we introduce the bilinear form:
a((u,p),(v,q)) := \int_{\Omega} \nu(c) \nabla u : \nabla v - \operatorname{div}(v) p - \operatorname{div}(u) q
IMEX time discretization for the unsteady Navier Stokes equations
For the time integration we consider an semi-implicit Euler method (IMEX) where the convection is treated only explicitly and the remaining part implicitly:
Find (\mathbf{u}^{n+1},p^{n+1}) \in X = \mathbf{V} \times Q, s.t. for all (\mathbf{v},q) \in X = \mathbf{V} \times Q
\begin{align*}
\underbrace{m(\mathbf{u}^{n+1},\mathbf{v}) + \Delta t ~\cdot~a((\mathbf{u}^{n+1},p^{n+1}),(\mathbf{v},q))}_{ \to M^\ast} ~=~m(\mathbf{u}^{n},\mathbf{v}) - \Delta t ~\cdot~c(\mathbf{u}^{n}; \mathbf{u}^{n},\mathbf{v})
\end{align*}
with
\begin{align*}
m(\mathbf{u},\mathbf{v}) = \int_{\Omega} \mathbf{u} \cdot \mathbf{v}
\end{align*}
and
\begin{align*}
c(\mathbf{w},\mathbf{u},\mathbf{v}) = \int_{\Omega} \mathbf{w} \cdot \nabla \mathbf{u} \cdot \mathbf{v}
\end{align*}
We prefer the incremental form (as it homogenizes the boundary conditions):
\begin{align*}
m(\delta \mathbf{u}^{n+1},\mathbf{v}) + \Delta t ~\cdot~a((\delta \mathbf{u}^{n+1},\delta p^{n+1}),(\mathbf{v},q)) ~=~ \Delta t (-a((\mathbf{u}^{n},p^n),(\mathbf{v},q)) -c(\mathbf{u}^{n}; \mathbf{u}^{n},\mathbf{v}))
\end{align*}
Explicit time stepping with a DG formulation for the unsteady scalar linear transport problem
We consider an Upwind DG discretization (in space):
Find c: [0,T] \to V_h := \bigoplus_{T\in\mathcal{T}_h} \mathcal{P}^k(T) so that
\sum_{T} \int_T \partial_t c d - u \cdot \nabla d c + \int_{\partial T} u_n \hat{c} d = \int_{\Omega} f d, \quad \forall d \in V_h.
Here \hat{c} is the Upwind flux, i.e. \hat{c} = c on the outflow boundary part \partial T_{out} = \{ x\in \partial T \mid u(x) \cdot n_T(x) \geq 0 \} of T and \hat{c} = c^{other} else, with c^{other} the value from the neighboring element.
Upwind-in-space combined with an explicit Euler scheme in time yields:
\sum_{T} \int_T c^{n+1} d = \sum_{T} \int_T c^{n} d - \Delta t \sum_{T} \left\{ - \int_T u \cdot \nabla d c
+ \int_{\partial T} u_n \hat{c} d \right\} + \Delta t \int_{\Omega} f d, \quad \forall d \in V_h,
M c^{n+1} = M c^{n} - \Delta t C(c^n) + \Delta t f
In our first example we set f = 0.
\begin{align*}
{m({c}^{n+1},{d})} ~=~m({c}^{n},{d}) - \Delta t ~\cdot~c(\mathbf{u}^{n}; {c}^{n},{d})
\end{align*}
with
\begin{align*}
m({c},{d}) = \int_{T} {c} \cdot {d}
\end{align*}
and
\begin{align*}
c(\mathbf{e}; {c},{d}) = - \int_{T} \mathbf{e} \cdot \nabla {d} \cdot {c} + \int_{\partial T} e_n \hat{c} d
\end{align*}
We prefer the incremental form (as it homogenizes the boundary conditions):
\begin{align*}
{m({\delta c}^{n+1},{d})} ~=~ - \Delta t ~\cdot~c({c}^{n}; {c}^{n},{d})
\end{align*}
combined form
m((\delta \mathbf{u}^{n+1},\delta c^{n+1}),(\mathbf{v},d)) + \Delta t ~\cdot~a((\delta \mathbf{u}^{n+1},\delta p^{n+1}),(\mathbf{v},q)) ~=~ \Delta t (-a((\mathbf{u}^{n},p^n),(\mathbf{v},q)) -c(\mathbf{u}^{n}; (\mathbf{u}^{n},c^{n}),(\mathbf{v},d)))
with :
m((\mathbf{u},c),(\mathbf{v},d)) = \int_{\Omega} \mathbf{u} \cdot \mathbf{v} + \int_{T} {c} \cdot {d}
a((u,p),(v,q)) := \int_{\Omega} \nu \nabla u : \nabla v - \operatorname{div}(v) p - \operatorname{div}(u) q
c(\mathbf{u}; (\mathbf{u},c),(\mathbf{v}.d))= \int_{\Omega} \mathbf{u} \cdot \nabla \mathbf{u} \cdot \mathbf{v} - \int_{T} \mathbf{u} \cdot \nabla {d} \cdot {c} + \int_{\partial T} u_n \hat{c} d
but when i use the following code, i find that the mass matrix can not be calculated in the above way shown in the formula, can somebody give me some advice?
k = 2
# V = VectorH1(mesh,order=k, dirichlet="Walls|In")
V = VectorH1(mesh,order=k, dirichlet="Walls|In",dirichlety="Bottoms")
Q = H1(mesh,order=k-1)
D = L2(mesh,order=k)
X = V*Q*D
(u,p,c), (v,q,d) = X.TnT()
a = BilinearForm(X, nonassemble = True)
nu = IfPos(c,500,0.02)
stokes = (nu*InnerProduct(grad(u),grad(v))-div(u)*q-div(v)*p-1e-8*p*q)*dx
a += stokes
dt = 1e-5
mstar = BilinearForm(X)
# D.Mass(1) since use the DG
mstar += InnerProduct(u,v)*dx + dt*stokes + D.Mass(1)
references:
[1] Operator applications for Discontinuous Galerkin discretizations - linear transport
[2] Incompressible Navier-Stokes equations