A time dependent Stokes interface problem: well-posedness and space-time finite element discretization

TL;DR

This paper develops a well-posed variational formulation for a time-dependent Stokes interface problem with discontinuous coefficients, and introduces a space-time finite element method with numerical validation for two-phase flow simulations.

Contribution

The paper provides the first well-posed variational formulations for the two-phase Stokes interface problem, including divergence-free and pressure-based variants, and proposes a space-time finite element discretization.

Findings

Well-posed variational formulations are established.

A space-time finite element method is introduced.

Numerical experiments demonstrate the method's effectiveness.

Abstract

In this paper a time dependent Stokes problem that is motivated by a standard sharp interface model for the fluid dynamics of two-phase flows is studied. This Stokes interface problem has discontinuous density and viscosity coefficients and a pressure solution that is discontinuous across an evolving interface. This strongly simplified two-phase Stokes equation is considered to be a good model problem for the development and analysis of finite element discretization methods for two-phase flow problems. In view of the unfitted finite element methods that are often used for two-phase flow simulations, we are particularly interested in a well-posed variational formulation of this Stokes interface problem in a Euclidean setting. Such well-posed weak formulations, which are not known in the literature, are the main results of this paper. Different variants are considered, namely one with suitable spaces of divergence free functions, a discrete-in-time version of it, and variants in which the divergence free constraint in the solution space is treated by a pressure Lagrange multiplier. The discrete-in-time variational formulation involving the pressure variable for the divergence free constraint is a natural starting point for a space-time finite element discretization. Such a method is introduced and results of numerical experiments with this method are presented.

Tables

Table 1. Table 1: Error ‖ 𝐮 − 𝐮 h ‖ L 2 ⊗ H 1 subscript norm 𝐮 subscript 𝐮 ℎ tensor-product superscript 𝐿 2 superscript 𝐻 1 \|\mathbf{u}-\mathbf{u}_{h}\|_{L^{2}\otimes H^{1}} for finite element spaces U h b superscript subscript 𝑈 ℎ 𝑏 U_{h}^{b} and Q h subscript 𝑄 ℎ Q_{h} . The estimated temporal (spatial) order of convergence E ​ O ​ C T 𝐸 𝑂 subscript 𝐶 𝑇 EOC_{T} ( E ​ O ​ C S 𝐸 𝑂 subscript 𝐶 𝑆 EOC_{S} ) is computed using the last row (column).

$N_S$\$N$	4	8	16	32	64	128	$EO{C}_S$
4	0.32167	0.25902	0.25045	0.24888	0.24852	0.24844
8	0.21797	0.13582	0.12703	0.12597	0.12577	0.12573	0.98264
16	0.18799	0.09090	0.08163	0.08095	0.08088	0.08086	0.63669
32	0.17626	0.06511	0.05425	0.05450	0.05466	0.05469	0.56423
$EO{C}_T$		1.43670	0.26328	-0.00674	-0.00421	-0.00068

Table 2. Table 2: Error ‖ 𝐮 − 𝐮 h X ‖ L 2 ⊗ H 1 subscript norm 𝐮 superscript subscript 𝐮 ℎ 𝑋 tensor-product superscript 𝐿 2 superscript 𝐻 1 \|\mathbf{u}-\mathbf{u}_{h}^{X}\|_{L^{2}\otimes H^{1}} for finite element spaces U h b superscript subscript 𝑈 ℎ 𝑏 U_{h}^{b} and Q h X superscript subscript 𝑄 ℎ 𝑋 Q_{h}^{X} . The estimated temporal (spatial) order of convergence E ​ O ​ C T 𝐸 𝑂 subscript 𝐶 𝑇 EOC_{T} ( E ​ O ​ C S 𝐸 𝑂 subscript 𝐶 𝑆 EOC_{S} ) is computed using the last row (column).

$N_S$\$N$	4	8	16	32	64	128	$EO{C}_S$
4	0.29649	0.20900	0.19753	0.19552	0.19510	0.19507
8	0.18572	0.06802	0.04699	0.04390	0.04332	0.04318	2.17546
16	0.17339	0.04718	0.01736	0.01154	0.01064	0.01047	2.04410
32	0.17604	0.04423	0.01326	0.00525	0.00306	0.00267	1.97412
$EO{C}_T$		1.99289	1.73737	1.33761	0.77615	0.20161

Table 3. Table 3: Error ‖ p − p h ‖ L 2 ⊗ L 2 subscript norm 𝑝 subscript 𝑝 ℎ tensor-product superscript 𝐿 2 superscript 𝐿 2 \|p-p_{h}\|_{L^{2}\otimes L^{2}} for finite element spaces U h b superscript subscript 𝑈 ℎ 𝑏 U_{h}^{b} and Q h subscript 𝑄 ℎ Q_{h} . The estimated temporal (spatial) order of convergence E ​ O ​ C T 𝐸 𝑂 subscript 𝐶 𝑇 EOC_{T} ( E ​ O ​ C S 𝐸 𝑂 subscript 𝐶 𝑆 EOC_{S} ) is computed using the last row (column).

$N_S$\$N$	4	8	16	32	64	128	$EO{C}_S$
4	2.30681	2.22456	2.20904	2.20538	2.20451	2.20431
8	1.58330	1.57570	1.57848	1.57931	1.57950	1.57956	0.48080
16	1.18787	1.15924	1.17049	1.17305	1.17346	1.17348	0.42873
32	0.93573	0.83261	0.83562	0.84323	0.84483	0.84507	0.47365
$EO{C}_T$		0.16845	-0.00521	-0.01308	-0.00273	-0.00041

Table 4. Table 4: Error ‖ p − p h X ‖ L 2 ⊗ L 2 subscript norm 𝑝 superscript subscript 𝑝 ℎ 𝑋 tensor-product superscript 𝐿 2 superscript 𝐿 2 \|p-p_{h}^{X}\|_{L^{2}\otimes L^{2}} for finite element spaces U h b superscript subscript 𝑈 ℎ 𝑏 U_{h}^{b} and Q h X superscript subscript 𝑄 ℎ 𝑋 Q_{h}^{X} . The estimated temporal (spatial) order of convergence E ​ O ​ C T 𝐸 𝑂 subscript 𝐶 𝑇 EOC_{T} ( E ​ O ​ C S 𝐸 𝑂 subscript 𝐶 𝑆 EOC_{S} ) is computed using the last row (column).

$N_S$\$N$	4	8	16	32	64	128	$EO{C}_S$
4	1.81460	1.11767	0.99726	0.96497	0.97083	0.98616
8	0.60974	0.26134	0.17501	0.16322	0.16244	0.16315	2.59559
16	0.72751	0.27060	0.10682	0.04883	0.04331	0.04355	1.90560
32	1.79975	0.40515	0.17825	0.07309	0.02646	0.01895	1.20034
$EO{C}_T$		2.15127	1.18458	1.28605	1.46568	0.48186

Table 5. Table 5: Error ‖ 𝐮 − 𝐮 h X ‖ L 2 ⊗ H 1 subscript norm 𝐮 superscript subscript 𝐮 ℎ 𝑋 tensor-product superscript 𝐿 2 superscript 𝐻 1 \|\mathbf{u}-\mathbf{u}_{h}^{X}\|_{L^{2}\otimes H^{1}} for finite element spaces U h b superscript subscript 𝑈 ℎ 𝑏 U_{h}^{b} and Q h X superscript subscript 𝑄 ℎ 𝑋 Q_{h}^{X} . The estimated temporal (spatial) order of convergence E ​ O ​ C T 𝐸 𝑂 subscript 𝐶 𝑇 EOC_{T} ( E ​ O ​ C S 𝐸 𝑂 subscript 𝐶 𝑆 EOC_{S} ) is computed using the last row (column).

$N_S$\$N$	4	8	16	32	64	128	$EO{C}_S$
4	0.13430	0.04973	0.04391	0.04261	0.04213	0.04199
8	0.07766	0.02683	0.01856	0.01771	0.01757	0.01758	1.25641
16	0.08380	0.02728	0.01185	0.00828	0.00787	0.00784	1.16538
32	0.11604	0.03314	0.01404	0.00696	0.00492	0.00473	0.72833
$EO{C}_T$		1.80811	1.23931	1.01213	0.50098	0.05590

Equations

$$\{\begin{array}{rl}{\rho }_i(\frac{\partial \mathbf{u}}{\partial t}+(\mathbf{u}\cdot \nabla )\mathbf{u})& =\mathrm{div}{\mathbf{\sigma }}_i+{\mathbf{g}}_i\\ \mathrm{div}\mathbf{u}& =0\end{array}$$

$$[\mathbf{\sigma }{\mathbf{n}}_{\Gamma }]=-\tau \kappa {\mathbf{n}}_{\Gamma }\text{on}\Gamma (t),$$

$$[\mathbf{u}]=0\text{on}\Gamma (t),$$

$$V_{\Gamma }=\mathbf{u}\cdot {\mathbf{n}}_{\Gamma }\text{on}\Gamma (t).$$

$$\dot{v}:=\frac{\partial v}{\partial t}+\mathbf{w}\cdot \nabla v.$$

$$\{\begin{array}{rl}\rho \dot{\mathbf{u}}-\mathrm{div}(\mu D(\mathbf{u}))+\nabla p& ={\mathbf{g}}_i\\ \mathrm{div}\mathbf{u}& =0\end{array}$$

$$[(-p\mathbf{I}+\mu D(\mathbf{u})){\mathbf{n}}_{\Gamma }]=-\tau \kappa {\mathbf{n}}_{\Gamma }\text{on}\Gamma (t),$$

$$[\mathbf{u}]=0\text{on}\Gamma (t),$$

$$\mathbf{u}(x,t)=0\quad\text{for}~{}~{}(x,t)\in\big{(}\partial\Omega\times I\big{)}\cup\big{(}\Omega\times\{0\}\big{)}.$$

$$\mathcal{V}:=\{v\in H_0^1(\Omega {)}^d|\mathrm{div}v=0\},X:=L^2(I;\mathcal{V}).$$

$$(\rho \dot{\mathbf{u}},\mathbf{v}{)}_{L^2}+(\mu D(\mathbf{u}),D(\mathbf{v}){)}_{L^2}=(\rho \mathbf{g},\mathbf{v}{)}_{L^2}-\tau {\int }_0^T{\int }_{\Gamma (t)}\kappa {\mathbf{n}}_{\Gamma }\cdot \mathbf{v}dsdt,$$

$$\{\begin{array}{l}\frac{\partial \Phi }{\partial t}(y,t)=\mathbf{w}(\Phi (y,t),t),t\in I,\\ \Phi (y,0)=y.\end{array}$$

$$(\mathbf{u},\mathbf{v})_{X}:=\int_{0}^{T}\big{(}\mathbf{u}(t),\mathbf{v}(t)\big{)}_{1,\Omega}\,dt=\sum_{i=1}^{d}\int_{0}^{T}\int_{\Omega}\nabla u_{i}(t)\cdot\nabla v_{i}(t)\,dx\,dt=(\nabla\mathbf{u},\nabla\mathbf{v})_{L^{2}}.$$

$${\mathcal{D}}_0:=\{\sum _{i=1}^n{g}_i{\mathbf{\varphi }}_i|n\in \mathbb{N},g_i\in C_0^{\infty }(I),{\mathbf{\varphi }}_i\in C_0^1(\Omega {)}^d\cap \mathcal{V}\}\subset C_0^1(Q{)}^d$$

$$⟨\rho \dot{\mathbf{v}},\mathbf{\varphi }⟩=-(\rho \mathbf{v},\dot{\mathbf{\varphi }}{)}_{L^2}\text{for}\mathbf{\varphi }\in {\mathcal{D}}_0.$$

$$W=\{\mathbf{v}\in X|\rho \dot{\mathbf{v}}\in X^{\prime }\},\Vert \mathbf{v}{\Vert }_W^2=\Vert \mathbf{v}{\Vert }_X^2+\Vert \rho \dot{\mathbf{v}}{\Vert }_{X^{\prime }}^2.$$

$$(\rho \mathbf{v},\dot{\mathbf{\varphi }}{)}_{L^2}=(\rho \mathbf{v},\frac{\partial \mathbf{\varphi }}{\partial t}{)}_{L^2}+(\rho \mathbf{v},\mathbf{w}\cdot \nabla \mathbf{\varphi }{)}_{L^2}.$$

$${\int }_0^T{\int }_{\Omega }\rho \dot{\mathbf{v}}\cdot \mathbf{\varphi }+\rho \mathbf{v}\cdot \dot{\mathbf{\varphi }}={\int }_0^T{\int }_{\Omega }\rho \dot{\overbrace{\mathbf{v}\cdot \mathbf{\varphi }}}$$

$$\displaystyle=\big{(}\rho(\cdot,T)\mathbf{v}(\cdot,T),\boldsymbol{\phi}(\cdot,T)\big{)}_{L^{2}(\Omega)}-\big{(}\rho(\cdot,0)\mathbf{v}(\cdot,0),\boldsymbol{\phi}(\cdot,0)\big{)}_{L^{2}(\Omega)}.$$

$$⟨\rho \dot{\mathbf{v}},\mathbf{\varphi }⟩=(\rho \dot{\mathbf{v}},\mathbf{\varphi }{)}_{L^2}\text{for}\mathbf{v}\in C^1(\bar{Q}{)}^d\cap X,\mathbf{\varphi }\in {\mathcal{D}}_0,$$

$$⟨\rho \dot{\mathbf{u}},\mathbf{v}⟩+(\mu D(\mathbf{u}),D(\mathbf{v}){)}_{L^2}=F(\mathbf{v})\text{for all}\mathbf{v}\in X.$$

$$U$$

$$V$$

$$U\to V\to W.$$

$${\overline{\mathcal{D}(C_0^1(\Omega {)}^d\cap \mathcal{V}{)}_{|I}}}^{\Vert \cdot {\Vert }_W}=V.$$

$$\underset{0\le t\le T}{\sup }\Vert \mathbf{u}(t){\Vert }_{L^2(\Omega )}\le c\Vert \mathbf{u}{\Vert }_V\text{for all}\mathbf{u}\in V,$$

$$\left\langle\rho\dot{\mathbf{u}},\mathbf{v}\right\rangle+\left\langle\rho\dot{\mathbf{v}},\mathbf{u}\right\rangle=\big{(}\rho(T)\mathbf{v}(T),\mathbf{u}(T)\big{)}_{L^{2}(\Omega)}-\big{(}\rho(0)\mathbf{v}(0),\mathbf{u}(0)\big{)}_{L^{2}(\Omega)}.$$

$$\Vert \rho (t_e{)}^{\frac{1}{2}}\mathbf{u}(t_e){\Vert }_{L^2(\Omega )}^2-\Vert \rho (t{)}^{\frac{1}{2}}\mathbf{u}(t){\Vert }_{L^2(\Omega )}^2=2(\rho \dot{\mathbf{u}},\mathbf{u}{)}_{L^2(\tilde{Q})}.$$

$$\Vert \mathbf{u}(t){\Vert }_{L^2(\Omega )}^2$$

$$\displaystyle\leq 2\rho_{\min}^{-1}\big{(}|(\rho\sigma^{\prime}\sigma\mathbf{u},\mathbf{u})_{L^{2}(\tilde{Q})}|+|(\rho\dot{\mathbf{u}},\sigma^{2}\mathbf{u})_{L^{2}(\tilde{Q})}|\big{)}.$$

$$|(\rho \dot{\mathbf{u}},{\sigma }^2\mathbf{u}{)}_{L^2(\tilde{Q})}|$$

$$\le c\underset{\mathbf{v}\in X}{\sup }\frac{(\rho \dot{\mathbf{u}},\mathbf{v}{)}_{L^2(Q)}}{\Vert \mathbf{v}{\Vert }_X}\Vert \mathbf{u}{\Vert }_X\le c\Vert \rho \dot{\mathbf{u}}{\Vert }_{X^{\prime }}\Vert \mathbf{u}{\Vert }_X\le c\Vert \mathbf{u}{\Vert }_V^2.$$

$$\Vert \mathbf{u}(t){\Vert }_{L^2(\Omega )}\le c\Vert \mathbf{u}{\Vert }_V,$$

$$⟨\rho \dot{\mathbf{u}},\mathbf{v}⟩+⟨\rho \dot{\mathbf{v}},\mathbf{u}⟩=(\rho \dot{\mathbf{u}},\mathbf{v}{)}_{L^2}+(\rho \dot{\mathbf{v}},\mathbf{u}{)}_{L^2}.$$

$$(\rho \dot{\mathbf{u}},\mathbf{v}{)}_{L^2}+{\int }_0^{T}a(\mathbf{u}(t),\mathbf{v}(t))dt={\int }_0^T(\mathbf{f}(t),\mathbf{v}(t){)}_{L^2(\Omega )}dt\text{for all}\mathbf{v}\in X.$$

Full text

1 Introduction

Let $\Omega\subset\mathbb{R}^{d}$ be an open bounded connected domain and $I:=(0,T)$ a time interval. On the space-time cylinder $\Omega\times I$ we consider the following standard sharp interface model (in strong formulation) for the fluid dynamics of a two-phase incompressible flow, cf. [2, 35, 22]:

[TABLE]

Here $\Gamma(t)=\overline{\Omega_{1}(t)}\cap\overline{\Omega_{2}(t)}$ denotes the (sharp) interface, \boldsymbol{\sigma}_{i}=-p\mathbf{I}+\mu_{i}\big{(}\nabla\mathbf{u}+(\nabla\mathbf{u})^{T}\big{)} the Newtonian stress tensor and $V_{\Gamma}$ is the normal velocity of the interface. The density and viscosity, $\rho_{i}$ and $\mu_{i},i=1,2,$ are assumed to be constant in each phase. The constant $\tau\geq 0$ is the surface tension coefficient and $\kappa$ is the mean curvature of $\Gamma$, i.e., $\kappa(x)=\mathop{\rm div}\mathbf{n}_{\Gamma}(x)$ for $x\in\Gamma$. Unknowns are the velocity $\mathbf{u}=\mathbf{u}(x,t)$, the pressure $p=p(x,t)$ and the (evolving) interface $\Gamma(t)$. To make the problem well-posed one needs suitable initial and boundary conditions for $\mathbf{u}$ and $\Gamma$. Due to the coupling of the interface dynamics and the fluid dynamics in the two bulk phases, this is a highly nonlinear problem. There is extensive literature on existence of solutions and well-posedness of different formulations of this problem. Most publications on these topics study quite regular solutions (in Hölder spaces) and deal with well-posedness locally in time or global existence of solutions close to equilibrium states (e.g., [14, 39, 15, 33, 34]). Often simplifying assumptions are used, for example, $\mathbf{g}_{i}=0$, $\tau=0$ or constant density ($\rho_{1}=\rho_{2}$). In other studies weaker solution concepts are used, for example, in [31, 32] the notion of renormalized solutions [16] of transport equations is used to derive an existence result (for $\tau=0$) and in [1] existence of so-called measure-valued varifold solutions is shown (for constant density). Here we do not give an overview of the extensive literature in this research field; for this we refer the interested reader to the literature discussion in the recent book [35].

We are interested in the development and analysis of finite element discretization methods for the two-phase flow problem given above. Finite element methods for this problem class can be found in, e.g., [22, 5, 13, 7]. We are not aware of any literature in which rigorous error analysis of such finite element methods is presented. Only very few partial results, e.g. on discrete stability as in [5], are known. This lack of analysis is clearly related to the strong nonlinearity of the problem (1)-(4). We also note that approaches and results available in the mathematical literature on existence of solutions and well-posedness of this problem turn out not to be very useful for the analysis of finite element discretization methods. In view of this, we introduce and analyze a much simpler (linear) Stokes interface problem which, however, is motivated by and closely related to the two-phase flow problem given above. We now derive this Stokes interface problem. In almost all numerical simulation methods for (1)-(4) one uses an iterative decoupling technique in which the interface evolution is decoupled from the flow problems in the subdomains. For the interface representation and numerical propagation one can use, for example, the level set method and given an approximation of $\Gamma(t)$ for $t$ in a (small) time interval one then discretizes the coupled Navier-Stokes equations in the subdomains. These Navier-Stokes equations are usually linearized by inserting a known approximation of the velocity in the first argument of the quadratic term $(\mathbf{u}\cdot\nabla)\mathbf{u}$. These two subproblems (interface propagation and solution of flow problem in the subdomains) can be coupled by several different iterative methods. This decoupling and linearization procedure motivates the following simplifying assumptions. Firstly, we assume a given sufficiently smooth (specified below) flow field $\mathbf{w}=\mathbf{w}(x,t)\in\mathbb{R}^{d}$ on $Q:=\Omega\times I$, with $\mathop{\rm div}\mathbf{w}=0$ on $Q$, which transports the interface. Instead of the interface dynamics condition $V_{\Gamma}=\mathbf{u}\cdot\mathbf{n}_{\Gamma}$ we impose $V_{\Gamma}=\mathbf{w}\cdot\mathbf{n}_{\Gamma}$. This implies that the interface evolution is completely determined by $\mathbf{w}$ and the Navier-Stokes flow problem in the two subdomains $\Omega_{i}(t)$ is decoupled from the interface dynamics. Secondly, we use a linearization of the Navier-Stokes equation in which $(\mathbf{u}\cdot\nabla)\mathbf{u}$ is replaced by $(\mathbf{w}\cdot\nabla)\mathbf{u}$. Thus we obtain a time dependent (generalized) Stokes problem (also called Oseen problem) in each of the subdomains, with coupling conditions as in (2)-(3). We introduce the usual notation for the material derivative along the flow field $\mathbf{w}$:

[TABLE]

We also introduce the piecewise constant functions $\rho,\mu$ with $\rho(x,t):=\rho_{i}$, $\mu(x,t):=\mu_{i}$ in $\Omega_{i}(t)$ and the deformation tensor $D(\mathbf{u}):=\nabla\mathbf{u}+(\nabla\mathbf{u})^{T}$. Thus we obtain the following much simpler linear problem: determine $\mathbf{u}$ and $p$ such that

[TABLE]

combined with suitable initial and boundary conditions for $\mathbf{u}$. We restrict to homogeneous Dirichlet boundary and initial conditions for $\mathbf{u}$:

[TABLE]

Both for the analysis and numerical simulations it is very convenient to reformulate this simplified model in a one-fluid Stokes interface model that combines the flow equations in the subdomains (5) and the interface conditions (6)-(7). We consider this Stokes interface problem to be an interesting and relevant subproblem for the numerical simulation of the full two-phase flow problem (1)-(4). For example, a finite element method that is stable and accurate for this Stokes interface problem can be expected to be an efficient discretization for the Navier-Stokes flow equations (with small Reynold’s numbers) in the full two-phase flow problem. The main contribution of this paper is the derivation of a well-posed space-time variational formulation of this Stokes interface model and, based on this, a (Galerkin) space-time finite element discretization.

We mention a few relevant properties of the interface Stokes problem (5)-(7). The discontinuity of the coefficients $\rho_{i},\mu_{i}$ across the interface and the interface force induced by the surface tension in (6) lead, even if the data is otherwise smooth, to a discontinuity in the pressure $p$ and to a discontinuity in the derivative of the velocity $\mathbf{u}$ on the space-time interface $\mathcal{S}$, cf. [35]. Hence, we have to deal with moving discontinuities. Typically the interface is not constant in time and thus we do not have a tensor product structure. These properties make this interface Stokes problem significantly more difficult to solve numerically than a standard time-dependent Stokes equation. Even for this strongly simplified problem we are not aware of any rigorous (sharp) error bounds for finite element discretization methods.

As a first step towards such an error analysis we need a suitable well-posed variational formulation. Concerning this we distinguish two different approaches. Firstly, the formulation and corresponding analysis is based on Lagrangian techniques, in which a suitable (coordinate) transformation is used to transform the given problem into one with a tensor product structure (i.e., a stationary interface). Such an approach is used in e.g. [35] (various parabolic two-phase problems) or [36] (free boundary Stokes problem). Such Lagrangian formulations are useful in the context of ALE (arbitrary Lagrangian Eulerian) discretizations and fitted finite elements. For a class of parabolic interface problems error bounds for fitted finite element methods have been derived in the literature, e.g. [11]. Alternatively, one can consider a formulation and analysis in an Eulerian setting (no coordinate transformations). Such formulations, which are standard for one-phase (Navier-)Stokes equations ([17, 41, 43]) are better suited for unfitted finite element techniques. Finite element methods for fluid-stucture interaction have analysed in this setting in [30, 37]. If in the original two-phase flow problem an interface capturing method such as the very popular level set method is used, this very often leads to the application of unfitted finite element discretization methods for the flow problem (meaning that the triangulations are not fitted to the evolving interface). This then requires special finite element spaces, for example an XFEM [20, 28], unfitted FEM with a Nitsche penalty term [24, 6] or a CutFEM [9, 8, 10]. In this paper we restrict to the Eulerian approach.

Hence, for the time dependent Stokes interface problem described above we are interested in a well-posed variational formulation in an Euclidean setting, similar to those for one-phase (Navier-)Stokes equations known in the literature ([17, 41, 43]). For one-phase (Navier-)Stokes equations new well-posed space-time formulations have been developed in the recent papers [23, 38]. These formulations do not cover the Stokes interface model described above, due to the lack of a tensor product structure. It turns out that in particular the discontinuity in the mass density $\rho$ across the interface causes significant difficulties considering the analysis of well-posedness, as explained in Remark 2.3. As a main contribution of this paper we develop an analysis resulting in a well-posed space-time variational formulation. Main results on well-posedness are given in Corollary 8, Theorem 9 and Theorem 10.

Our analysis is rather different from the analyses used in the derivation of a well-posed variational one-phase Stokes problems [17, 41, 43].

Based on this space-time variational formulation we propose an space-time unfitted finite element method. The method combines standard Discontinuous Galerkin time discretization [3, 42, 40] with an XFEM or CutFEM approach [20, 28, 8, 10] to account for the jump in pressure across the space-time interface $\mathcal{S}.$ We present results of numerical experiments with this method. An error analysis of this method is a topic of current research and not considered in this paper.

The remainder of the paper is organized as follows. In Section 2 we introduce a variational formulation of the Stokes interface problem (5)-(7) in an obvious space-time Sobolev space of divergence free functions. In Remark 2.3 we explain why the analysis of well-posedness of this formulation is problematic. This motivates the introduction of other (related) spaces, for which well-posedness of a variational formulation can be proved. This analysis is presented in Section 3. In Section 4 a standard discontinuous Galerkin approach is applied to derive a well-posed space-time variational formulation that allows a time stepping procedure. In Section 5 we study a space-time variational problem involving the pressure variable to satisfy the divergence free constraint. Based on this variational formulation we introduce an unfitted space-time finite element method in Section 6 and give results of numerical experiments with this method. We finally give a summary and outlook in Section 7.

2 Space-time variational formulation

We start with an assumption concerning the required smoothness of the space-time interface $\mathcal{S}:=\bigcup_{t\in I}\Gamma(t)\times\{t\}$ and the given velocity field $\mathbf{w}$.

Assumption 2.1**.**

Throughout the paper we assume that $\mathcal{S}$ is a connected Lipschitz hypersurface in $\mathbb{R}^{d+1}$ and that the given velocity field $\mathbf{w}$ is divergence free and $\mathbf{w}\in C(\bar{I};L^{2}(\Omega)^{d})$. The latter guarantees that the material derivative $\dot{v}=\frac{\partial v}{\partial t}+\mathbf{w}\cdot\nabla v$ is well-defined in a weak sense as in [16]. The piecewise constant density $\rho$ and the velocity field $\mathbf{w}$ are assumed to satisfy the compatibility condition $\dot{\rho}=0$.

Furthermore, we make the assumption $\mathbf{w}\in L^{\infty}(Q)^{d}$. This condition can be replaced by another (more natural) one which depends on the dimension $d$, cf. Remark 3.2.

As is usually done in the analysis of (Navier-)Stokes equations, we restrict to suitable subspaces of divergence free velocity fields and thus eliminate the pressure. We derive well-posedness of a suitable variational formulation in these subspaces. Therefore we introduce the spaces

[TABLE]

Assume that the strong formulation (5)-(7) has a sufficiently smooth solution $\mathbf{u}$. Multiplication by test function $\mathbf{v}\in X$ and partial integration then implies:

[TABLE]

where $(\cdot,\cdot)_{L^{2}}$ denotes the (vector) $L^{2}$ scalar product over the space-time cylinder $Q$.

Remark 2.1**.**

Note that the second term in the right hand-side of (9) corresponds to a force that acts only on the space-time interface $\mathcal{S}$. This induces a discontinuity in the pressure Lagrange multiplier. Below, instead of the specific right hand-side in (9) we consider a generic $F\in X^{\prime}$. If $\mathbf{g}\in L^{2}(Q)^{d}$ and $\kappa\mathbf{n}_{\Gamma}\in L^{2}(\mathcal{S})^{d}$ then the right hand-side satisfies $F\in X^{\prime}$. In order to have the normal $\mathbf{n}_{\Gamma}$ and the curvature $\kappa$ in the classical (strong) sense, we need additional ($C^{2}$) smoothness of $\Gamma(t)$. The regularity of $\Gamma(t)$ depends on the regularity of $\Gamma(0)$ and $\mathbf{w}$ in the following way. The advection field $\mathbf{w}$ defines a Lagrangian flow $\Phi$ (see [12]): for a given $y\in\Omega$ the function $t\mapsto\Phi(y,t)$ is defined by the ODE system

[TABLE]

For $\mathbf{w}\in C^{1}(\bar{I};C^{2}(\bar{\Omega}))^{d}$ this Lagrangian flow is uniquely defined and $\Phi\in C^{2}(\bar{Q})^{d}$. This is known from classical Cauchy-Lipschitz theory, see [12, §1.3]. Since $\Gamma(t)=\Phi(\Gamma(0),t)$, we can conclude that $\Gamma(t)$ is $C^{2}$ if $\mathbf{w}\in C^{1}(\bar{I};C^{2}(\bar{\Omega})^{d})$ and $\Gamma(0)$ is $C^{2}$.

Weaker notions of curvature have been developed for cases with less smoothness. This issue, however, is not relevant for the well-posedness results in the remainder of the paper.

A suitable weak material derivative can be defined in the standard distributional sense. For this we first introduce further notation. Elements $\mathbf{u}\in X$ have values $\mathbf{u}(x,t):=\mathbf{u}(t)(x)\in\mathbb{R}^{d}$, $(x,t)\in Q$. Due to the zero boundary values on $\partial\Omega$, the norm $\|v\|_{1}=(\|v\|_{L^{2}(\Omega)}^{2}+\|\nabla v\|_{L^{2}(\Omega)}^{2})^{\frac{1}{2}}$ on $H_{0}^{1}(\Omega)$ is equivalent to $|v|_{1}:=\|\nabla v\|_{L^{2}(\Omega)}$. In the remainder we use the latter norm, with corresponding scalar product denoted by $(\cdot,\cdot)_{1,\Omega}$ on $\mathcal{V}$. The scalar product on $X$ is denoted by

[TABLE]

Recall that $C_{0}^{1}(\Omega)^{d}\cap\mathcal{V}$ is dense in $\mathcal{V}$ and using the tensor product structure of $X$ we get that

[TABLE]

is dense in $X$, i.e., $\overline{\mathcal{D}_{0}}^{\|\cdot\|_{X}}=X$. For the case of an evolving interface and with the material derivative in (5) it is natural to introduce the following weak material derivative for functions from $X$. For $\mathbf{v}\in X$ we define the functional $\rho\dot{\mathbf{v}}$ by

[TABLE]

Note that in the $L^{2}$ scalar product we use a weighting with the strictly positive piecewise constant function $\rho$. We introduce the following analogon of the space $\{\mathbf{v}\in X~{}|~{}\frac{\partial\mathbf{v}}{\partial t}\in X^{\prime}\}$:

[TABLE]

An important difference between $\{\mathbf{v}\in X~{}|~{}\frac{\partial\mathbf{v}}{\partial t}\in X^{\prime}\}$ and $W$ is that, if $\rho$ varies with $t$ (i.e., $\rho_{1}\neq\rho_{2}$ and the interface is not stationary), the latter does not have a tensor product structure.

Remark 2.2**.**

Inserting the definition of the material derivative we get

[TABLE]

There is a constant $c$, which depends on $\|\mathbf{w}\|_{L^{\infty}(Q)}$, such that $|(\rho\mathbf{v},\mathbf{w}\cdot\nabla\boldsymbol{\phi})_{L^{2}}|\leq c\|\mathbf{v}\|_{X}\|\boldsymbol{\phi}\|_{X}$ for all $\mathbf{v}\in X$, $\boldsymbol{\phi}\in\mathcal{D}_{0}$. This implies that $\rho\dot{\mathbf{v}}\in X^{\prime}$ iff $\frac{\partial(\rho\mathbf{v})}{\partial t}\in X^{\prime}$, and $\|\rho\dot{\mathbf{v}}-\frac{\partial(\rho\mathbf{v})}{\partial t}\|_{X^{\prime}}\leq c\|\mathbf{v}\|_{X}$. Therefore the norms $\|\mathbf{v}\|_{X}+\|\frac{\partial(\rho\mathbf{v})}{\partial t}\|_{X^{\prime}}$ and $\|\mathbf{v}\|_{W}$ are equivalent.

For smooth functions $\mathbf{v},\boldsymbol{\phi}\in C^{1}(\bar{Q})^{d}\cap X$ we obtain, using [16, Theorem II.6] (applied to $\rho$ and $\mathbf{v}\cdot\boldsymbol{\phi}$) and $\mathop{\rm div}\mathbf{w}=0$, the following partial integration identity:

[TABLE]

For $\boldsymbol{\phi}\in\mathcal{D}_{0}$ the boundary terms vanish, and thus we get

[TABLE]

which means that the weak material derivative $\rho\dot{\mathbf{v}}$ can be identified with the function $\rho\dot{\mathbf{v}}$. By a continuity argument it follows that the result in (14) also holds for all $\boldsymbol{\phi}\in X$.

A natural weak formulation of (5)-(7) is as follows, cf. (9). Given $F\in X^{\prime}$, determine $\mathbf{u}\in W$ with $\mathbf{u}(0)=0$ and

[TABLE]

Remark 2.3**.**

As noted above, the spaces $X$ and $W$ are very natural ones. We are, however, not able to prove well-posedness of this formulation. The key difficulty is to show that smooth functions are dense in $W$. For the case that the mass density $\rho$ is constant or $\mathcal{S}$ does not depend on $t$ (stationary interface), density of smooth functions can be proved using mollification procedures in Bochner spaces as in e.g., [43, Chapter 25]. For the general case, however, we do not have a tensor product structure and these techniques fail. We tried to develop a mollification technique in the full space-time cylinder $\mathbb{R}^{d+1}$. Such a mollification needs to satisfy a commutation property between mollification and distributional differentation (12) (which involves the discontinuous function $\rho$) and furthermore must respect the divergence free property and the homogeneous Dirichlet boundary condition. We were not able to develop such a mollification technique. If we would have a density of smooth functions property of $W$, it can be shown that there is a bounded trace operator $W\to L^{2}(\Omega)^{d}$, $\mathbf{u}\to\mathbf{u}(\cdot,t)$, which ensures that $\mathbf{u}(0)$ is well-defined, and partial integration rules can be derived. Well-posedness of (15) can then be derived using fairly standard arguments as in e.g. [43, Chapter 26]. The density of smooth functions property, however, is an open problem.

Remark 2.4**.**

In (15) we consider a variational formulation in which a weak material derivative $\dot{\mathbf{u}}$ is scaled with $\rho$, as in (5). The scaling with $\rho$ (which does not have tensor product structure) causes significant difficultities in the theoretical analysis (Remark 2.3). One might consider a rescaling of the momentum equation in (5) that eliminates the $\rho$ term in front of the material derivative $\dot{\mathbf{u}}$. The two obvious possibilities are to introduce $\tilde{p}:=\rho^{-1}p$ or $\tilde{\mathbf{u}}:=\rho\mathbf{u}$. In both cases we rescale $\mu$, using $\tilde{\mu}:=\rho^{-1}\mu$. If we use $\tilde{p}$, then partial integration of the momentum equation (multiplied by a test function $\mathbf{v}$) over the domain $\Omega=\Omega_{1}(t)\cup\Omega_{2}(t)$ results is an interface term of the form $\int_{\Gamma(t)}[(-\tilde{p}\mathbf{I}+\tilde{\mu}D(\mathbf{u}))\mathbf{n}_{\Gamma}]v\,ds$. This term can not be treated as a natural interface term, because in the interface condition (6) we have the quantities $p,\,\ \mu$ instead of $\tilde{p},\,\tilde{\mu}$. If we use $\tilde{\mathbf{u}}$, then the term $\int_{\Gamma(t)}[(-p\mathbf{I}+\tilde{\mu}D(\tilde{\mathbf{u}}))\mathbf{n}_{\Gamma}]v\,ds$ occurs, which can be handled as a natural interface condition, due to $\tilde{\mu}D(\tilde{\mathbf{u}})=\mu D(\mathbf{u})$. However, from (7) we now obtain the interface condtion $[\rho^{-1}\tilde{\mathbf{u}}]=0$, which implies that one has to use the space $H_{0}^{1}(\Omega_{1}(t)\cup\Omega_{2}(t))^{d}$ for $\tilde{\mathbf{u}}$ instead of the (much nicer) space $H_{0}^{1}(\Omega)^{d}$ for $\mathbf{u}$. Using these rescalings we are not able to derive a simpler analysis for well-posedness and therefore we keep the original formulation (5)-(7), which is closer to physics.

3 Space-time variational formulation in modified spaces

As explained in Remark 2.3, we encounter difficulties in the analysis of well-posedness of the variational formulation (15) using the space $W$. In this section we study a variational formulation as in (15), but with $W$ replaced by a (possibly) smaller space $V$ (introduced below). The structure of the analysis is as follows. In Section 3.1 we introduce further spaces $U$ and $V$ and derive properties of these spaces. In Section 3.2 we study an intermediate weak formulation and derive a well-posedness result. In Section 3.3 we introduce and analyze the final desired weak formulation analogon of (15), with $W$ replaced by $V$. The main results are given in Corollary 8 and Theorem 9.

3.1 Spaces $U\subset V\subset W$

For $\mathbf{v}\in X=L^{2}(I;\mathcal{V})$ its weak derivative $\frac{\partial\mathbf{v}}{\partial t}\in\mathcal{D}^{\prime}(I;\mathcal{V}^{\prime})$ is defined in the usual distributional sense [43]. We define the spaces

[TABLE]

These are Hilbert spaces with continuous embeddings

[TABLE]

The norm $\|\cdot\|_{U}$ is equivalent to $\|\cdot\|_{H^{1}(Q)^{d}}$. The space $U$ has a tensor product structure and we can use standard arguments to show that smooth functions are dense in $U$. More precisely, let $\mathcal{D}(\mathcal{V})$ be the space of all functions $\mathbf{v}:\mathbb{R}\to\mathcal{V}$ which are infinitely differentiable and have a compact support. Then, cf. Lemma 25.1 in [43], $\mathcal{D}(\mathcal{V})_{|I}$ is dense in $U$. Using the density of $C_{0}^{1}(\Omega)^{d}\cap\mathcal{V}$ in $\mathcal{V}$ we obtain that the space of smooth functions $\mathcal{D}(C_{0}^{1}(\Omega)^{d}\cap\mathcal{V})_{|I}$ is dense in $U$. From the density of $U$ in $V$ we thus get the density of smooth functions in $V$:

[TABLE]

Remark 3.1**.**

It seems reasonable (based on analogous results for the tensor product case) to claim that $\overline{\mathcal{D}(C_{0}^{1}(\Omega)^{d}\cap\mathcal{V})_{|I}}^{\|\cdot\|_{W}}=W$ holds, i.e., $V=W$. We are, however, not able to prove this claim, cf. Remark 2.3. Note that the well-posedness result derived for $V$ in Corollary 8 below implies that either $V\neq W$ or the well-posedness result holds for $W$.

Using the density result (18) important properties of $V$ are derived in the following lemma.

Lemma 1**.**

(i) For a.e. $t\in[0,T]$ the trace operator $\mathbf{u}\to\mathbf{u}(\cdot,t)=\mathbf{u}(t)$ can be extended to a bounded linear operator from $V$ into $L^{2}(\Omega)^{d}$. Moreover, the inequality

[TABLE]

holds with a constant $c$ independent of $\mathbf{u}$.

(ii) For all $\mathbf{u},\mathbf{v}\in V$, the following integration by parts identity holds:

[TABLE]

Proof.

Take $t\in[0,\frac{1}{2}T]$ (the case $t\in[\frac{1}{2},T]$ can be treated with very similar arguments). Define $t_{e}:=t+\frac{1}{4}T$, $\tilde{I}:=(t,t_{e})$, $\tilde{Q}:=\Omega\times\tilde{I}\subset Q$. It suffices to prove the result in (19) for the dense subspace $\mathcal{D}:=\mathcal{D}(C_{0}^{1}(\Omega)^{d}\cap\mathcal{V})_{|I}$ of smooth functions. Take $\mathbf{u}\in\mathcal{D}$. The partial integration identity (13) on $\tilde{Q}$ yields

[TABLE]

Let $\sigma$ be a smooth decreasing scalar function with compact support and $\sigma(t)=1$, $\sigma(t_{e})=0$. Note that $\sigma\mathbf{u}\in\mathcal{D}$ holds. If in (21) we use $\sigma\mathbf{u}$, instead of $\mathbf{u}$, we get, with $\rho_{\min}:=\min\{\rho_{1},\rho_{2}\}$:

[TABLE]

Note that $|(\rho\sigma^{\prime}\sigma\mathbf{u},\mathbf{u})_{L^{2}(\tilde{Q})}|\leq c\|\mathbf{u}\|_{L^{2}(\tilde{Q})}^{2}\leq c\|\mathbf{u}\|_{X}^{2}\leq c\|\mathbf{u}\|_{V}^{2}$ holds. Furthermore, with $\tilde{X}:=L^{2}(\tilde{I};\mathcal{V})$, and extending $\mathbf{v}\in\tilde{X}$ by zero outside $\tilde{I}$, we have:

[TABLE]

Thus we get

[TABLE]

with a constant (depending on $T$) that is independent of $\mathbf{u}\in\mathcal{D}$. Due to density of $\mathcal{D}$ this proves the result in (19), and thus (i).

We consider (ii). Due to density and the continuity result in (19) it suffices to prove (20) for $\mathbf{u},\mathbf{v}\in\mathcal{D}$. The identity in (14) holds for $\boldsymbol{\phi}\in X$ and thus for $\mathbf{u},\mathbf{v}\in\mathcal{D}$ it follows from (14) that

[TABLE]

From this and the partial integration identity (13) the result (20) follows. ∎

3.2 Well-posed space-time variational formulation in $U$

We define $U_{0}:=\{\,\mathbf{u}\in U~{}|~{}\mathbf{u}(0)=0\,\}$, where $\mathbf{u}(0)$ is well-defined (in $L^{2}(\Omega)$ sense) due to (19). In the following theorem we treat a variational problem with a sufficiently smooth right hand-side $\mathbf{f}$ and a bilinear form $a(\cdot,\cdot)$ on $\mathcal{V}\times\mathcal{V}$ that is independent of $t$. These assumptions are such that we can apply a standard Galerkin procedure to show existence of a unique solution in $U_{0}$. This intermediate problem will be used in the next section to derive well-posedness of a weak formution as in (15), with $W$ replaced by the space $V\subset W$.

Theorem 2**.**

Take $\mathbf{f}\in C(\overline{I};L^{2}(\Omega)^{d})$ and let $a(\cdot,\cdot)$ be a continuous elliptic bilinear form on $\mathcal{V}\times\mathcal{V}$ (with norm $|\cdot|_{1}$) that does not depend on $t$. Then there exists a unique $\mathbf{u}\in U_{0}$ such that

[TABLE]

Furthermore

[TABLE]

holds, with a constant $c$ independent of $\mathbf{f}$.

Proof.

The proof is based on a standard Galerkin technique known in the literature, e.g. [18]. Let $(\mathbf{v}_{k})_{k\geq 1}$ be a total orthonormal set in $\mathcal{V}$ and define $\mathcal{V}_{m}:={\rm span}\{\mathbf{v}_{1},\ldots,\mathbf{v}_{m}\}$, $X_{m}:=L^{2}(I;\mathcal{V}_{m})$. We consider the following problem: determine $\mathbf{u}_{m}\in X_{m}$ with $\mathbf{u}_{m}(0)=0$ and such that:

[TABLE]

Using the representation $\mathbf{u}_{m}(t)=\sum_{j=1}^{m}g_{j}(t)\mathbf{v}_{j}$ and with $\mathbf{g}_{m}(t):=(g_{1}(t),\ldots,g_{m}(t))^{T}$ this problem can be reformulated as a system of ODEs:

[TABLE]

with a symmetric positive definite matrix $M_{m}\in C(\bar{I};\mathbb{R}^{m\times m})$, $(M_{m}(t))_{i,j}=(\rho(t)\mathbf{v}_{j},\mathbf{v}_{i})_{L^{2}(\Omega)}$ and $B_{m}\in C(\bar{I};\mathbb{R}^{m\times m})$, $(B_{m}(t))_{i,j}=(\mathbf{w}(\cdot,t)\cdot\nabla\mathbf{v}_{j},\mathbf{v}_{i})_{L^{2}(\Omega)}+a(\mathbf{v}_{j},\mathbf{v}_{i})$ and $F_{m}\in C(\bar{I};\mathbb{R}^{m})$, $(F_{m}(t))_{i}=(\mathbf{f}(t),\mathbf{v}_{i})_{L^{2}(\Omega)}$. Standard theory for ODEs implies that (25) has a unique solution $\mathbf{g}_{m}\in C^{1}(\bar{I})^{m}$, and thus (24) has a unique solution $\mathbf{u}_{m}$. We take $\mathbf{v}=\mathbf{u}_{m}$ in (24):

[TABLE]

The ellipticity of $a(\cdot,\cdot)$ on $\mathcal{V}$ implies that $\int_{0}^{T}a(\mathbf{u}_{m}(t),\mathbf{u}_{m}(t))\,dt\geq\gamma\|\mathbf{u}_{m}\|_{X}^{2}$ for a $\gamma>0$ independent of $\mathbf{u}_{m}$. Combining this with partial integration, a Cauchy inequality and $\mathbf{u}_{m}(0)=0$ yields

[TABLE]

which implies a uniform bound $\|\mathbf{u}_{m}\|_{X}\leq\gamma^{-1}\|\mathbf{f}\|_{L^{2}}$.

We take $\mathbf{v}=\frac{\partial\mathbf{u}_{m}}{\partial t}=\sum_{j=1}^{m}g_{j}^{\prime}(t)\mathbf{v}_{j}\in X_{m}$ in (24), and thus get:

[TABLE]

From $a(\mathbf{u}_{m}(t),\mathbf{u}_{m}(t))=\mathbf{g}_{m}(t)^{T}A\mathbf{g}_{m}(t)$, with $A_{i,j}=a(\mathbf{v}_{i},\mathbf{v}_{j})$ and $\mathbf{g}_{m}(0)=0$ it follows that $a(\mathbf{u}_{m}(0),\mathbf{u}_{m}(0))=0$. Using this we get

[TABLE]

Using Cauchy inequalities and the uniform bound $\|\mathbf{u}_{m}\|_{X}\leq\gamma^{-1}\|\mathbf{f}\|_{L^{2}}$ we obtain $\|\frac{\partial\mathbf{u}_{m}}{\partial t}\|_{L^{2}}\leq c\|\mathbf{f}\|_{L^{2}}$ with a constant $c$ which only depends on $\rho,\|\mathbf{w}\|_{\infty}$ and $\gamma$. Hence we have a uniform boundedness result

[TABLE]

Hence there is a subsequence, which we also denote by $(\mathbf{u}_{m})_{m\geq 0}$, that weakly converges $\mathbf{u}_{m}\rightharpoonup\mathbf{u}\in U$, which implies $\mathbf{u}_{m}\rightharpoonup\mathbf{u}$ in $X$ and $\frac{\partial\mathbf{u}_{m}}{\partial t}\rightharpoonup\frac{\partial\mathbf{u}}{\partial t}$ in $L^{2}(Q)$. Passing to the limit and using continuity arguments we conclude that $\mathbf{u}\in U$ satisfies (22). We now show that $\mathbf{u}(0)=0$ holds, i.e., $\mathbf{u}\in U_{0}$. Take an arbitrary $\mathbf{v}\in C^{1}(\bar{I};\mathcal{V}_{N})\subset X_{N}$ with $\mathbf{v}(T)=0$. From (22) and partial integration we obtain

[TABLE]

We also get from (24), for $m\geq N$, and using $\mathbf{u}_{m}(0)=0$:

[TABLE]

Comparing (28), (29) and using $\mathbf{u}_{m}\rightharpoonup\mathbf{u}$ in $U$ it follows that $((\rho(0)\mathbf{u}(0),\mathbf{v}(0))_{L^{2}(\Omega)}=0$ holds. This implies $\mathbf{u}(0)=0$ in $L^{2}(\Omega)$. To show the uniqueness of $\mathbf{u}$ we take $\mathbf{f}=0$ and $\mathbf{v}=\mathbf{u}$ in (22):

[TABLE]

Using $(\rho\dot{\mathbf{u}},\mathbf{u})_{L^{2}}=\frac{1}{2}\|\rho^{\frac{1}{2}}(T)\mathbf{u}(T)\|_{L^{2}(\Omega)}^{2}$ and the ellipticity of $a(\cdot,\cdot)$ it follows that $\|\mathbf{u}\|_{X}=0$, hence we have uniqueness. The bound in (23) follows from (27). ∎

The assumption that the bilinear form $a(\cdot,\cdot)$ is independent of $t$ is used for the derivation of a bound on $\|\frac{\partial\mathbf{u}}{\partial t}\|_{L^{2}}$, for which the estimate in (26) is a key ingredient. Using a standard Gronwall argument in (26), one can derive a similar result if the bilinear form $a(t;\cdot,\cdot)$ is time dependent and differentiable with respect to time. Such arguments, however, fail when $a(t;\cdot,\cdot)$ is not differentiable, which is the case we consider. The extension to a time-dependent bilinear form $a(t;\cdot,\cdot)$, with a possibly nonsmooth dependence on $t$, is treated in Theorem 7.

Remark 3.2**.**

In the proof above we used the assumption $\mathbf{w}\in L^{\infty}(Q)^{d}$. This assumption can be replaced by a different (more natural) assumption by using alternative estimates for the trilinear form $(\rho\mathbf{w}\cdot\nabla\mathbf{u},\mathbf{v})$ which depend on the dimension $d$, see [41, Section 2.3]. The assumption $\mathbf{w}\in L^{\infty}(Q)$ can be replaced by $\mathbf{w}\in V$ for $d=2$. For $d=3$ we additionally need $\mathbf{w}\in L^{4}(I;H^{1}(\Omega)^{3})$.

Corollary 3**.**

Using that $C(\overline{I};L^{2}(\Omega)^{d})$ is dense in $L^{2}(I;L^{2}(\Omega)^{d})$ one can now derive the following well-posedness result: for each $\mathbf{f}\in L^{2}(I;L^{2}(\Omega)^{d})$ there exists a unique $\mathbf{u}\in U_{0}$ such that (22) and (23) hold.

3.3 Well-posed space-time variational formulation in $V$

We define $V^{0}:=\{\,\mathbf{v}\in V~{}|~{}\mathbf{v}(0)=0\,\}$, where the trace is well-defined due to (19). As an easy consequence of the result obtained in Theorem 2 we obtain the following.

Corollary 4**.**

Let $a(\cdot,\cdot)$ be a continuous elliptic bilinear form on $\mathcal{V}\times\mathcal{V}$ that does not depend on $t$. For every $F\in X^{\prime}$ there exists a unique $\mathbf{u}\in V^{0}$ such that

[TABLE]

Furthermore

[TABLE]

holds, with a constant $c$ independent of $F$.

Proof.

Take $F\in X^{\prime}$. Due to the density of $C(\bar{I};L^{2}(\Omega)^{d})$ in $X^{\prime}$ we can take a sequence $\mathbf{f}_{n}\in C(\bar{I};L^{2}(\Omega)^{d})$, $n\in\mathbb{N}$, with $\lim_{n\to\infty}\mathbf{f}_{n}=F$ in $X^{\prime}$. Let $\mathbf{u}_{n}\in U_{0}$ be the unique solution of (22). As test function we take $\mathbf{v}=\mathbf{u}_{n}$ in (22). Using partial integration, $\mathbf{u}_{n}(0)=0$ and ellipticity of $a(\cdot,\cdot)$ we get $\gamma\|\mathbf{u}_{n}\|_{X}^{2}\leq\|\mathbf{f}_{n}\|_{X^{\prime}}\|\mathbf{u}_{n}\|_{X}$, with ellipticity constant $\gamma>0$, and thus $\|\mathbf{u}_{n}\|_{X}\leq\gamma^{-1}\|\mathbf{f}_{n}\|_{X^{\prime}}$. This implies that $(\mathbf{u}_{n})_{n\in\mathbb{N}}$ is a Cauchy sequence in $X$. Take $\mathbf{u}\in X$ such that $\lim_{n\to\infty}\mathbf{u}_{n}=\mathbf{u}$ in $X$. Note that

[TABLE]

Hence, $\|\rho\dot{\mathbf{u}}_{n}\|_{X^{\prime}}\leq c(\|\mathbf{u}_{n}\|_{X}+\|\mathbf{f}_{n}\|_{X^{\prime}})$. This implies that $(\rho\dot{\mathbf{u}}_{n})_{n\in\mathbb{N}}$ is a Cauchy sequence in $X^{\prime}$. Therefore $(\mathbf{u}_{n})_{n\in\mathbb{N}}$ is a Cauchy sequence in $V$ and $\lim_{n\to\infty}\rho\mathbf{u}_{n}=\rho\dot{\mathbf{u}}$ in $X^{\prime}$ holds. Thus we get $\lim_{n\to\infty}\mathbf{u}_{n}=\mathbf{u}$ in $V$. From this and the trace inequality (19) we get $\mathbf{u}(0)=0$, hence $\mathbf{u}\in V^{0}$. If in (32) we take $n\to\infty$ it follows that $\mathbf{u}$ satisfies (30). Uniqueness of $\mathbf{u}$ follows by taking $F=0$ and $\mathbf{v}=\mathbf{u}$ in (30), partial integration identity (20) and elliptcity of $a(\cdot,\cdot)$. From the estimates above we get $\|\mathbf{u}_{n}\|_{X}+\|\rho\dot{\mathbf{u}}_{n}\|_{X^{\prime}}\leq c\|\mathbf{f}_{n}\|_{X^{\prime}}$. Taking $n\to\infty$ we obtain the result in (31). ∎

If the (diffusion) coefficient $\mu$ in (15) would be constant, i.e., $\mu_{1}=\mu_{2}$ the result in Corollary 4 yields a well-posed weak formulation. In view of our applications, however, the case $\mu_{1}\neq\mu_{2}$ is highly relevant. Therefore, in the remainder of this section we present an analysis that can handle the latter case. In that analysis the result derived in Corollary 4 will play an important role.

For $F\in X^{\prime}$ we consider the following generalization of the problem in (30). Determine $\mathbf{u}\in V^{0}$ such that

[TABLE]

In the remainder of this section we assume that the (possibly) $t$-dependent bilinear form $a(t;\cdot,\cdot)$ has the following properties:

[TABLE]

In the remainder we prove well-posedness of the variational problem in (33). For this we first use the framework of the BNB-conditions, cf. [17], to prove well-posedness under the additional assumption that the bilinear form $a(t;\cdot,\cdot)$ is symmetric. We then extend the well-posedness result to $a(t;\cdot,\cdot)$ that may be nonsymmetric.

From $|\int_{0}^{T}a(t;\mathbf{u}(t),\mathbf{v}(t))\,dt|\leq\Gamma\int_{0}^{T}|\mathbf{u}(t)|_{1,\Omega}|\mathbf{v}(t)|_{1,\Omega}\,dt\leq\Gamma\|\mathbf{u}\|_{X}\|\mathbf{v}\|_{X}$ for all $\mathbf{u},\mathbf{v}\in X$ it follows that

[TABLE]

Hence $b(\cdot,\cdot)$ is continuous on $V^{0}\times X$.

Lemma 5**.**

The inf-sup inequality

[TABLE]

holds with $c_{s}=\frac{\sqrt{2}\,\gamma}{2(1+\Gamma^{2})}$.

Proof.

Take $\mathbf{u}\in V^{0}$. From the uniform ellipticity of $a(t;\cdot,\cdot)$ and the partial integration result (20), combined with $\mathbf{u}(0)=0$, we get

[TABLE]

This establishes the control of $\|\mathbf{u}\|_{X}$. We also need control of $\|\rho\dot{\mathbf{u}}\|_{X^{\prime}}$ to bound the full norm $\|\mathbf{u}\|_{V}$. This is achieved by using a duality argument between the Hilbert spaces $X$ and $X^{\prime}$. By Riesz’ representation theorem, there is a unique $\mathbf{z}\in X$ such that $\left\langle\rho\dot{\mathbf{u}},\mathbf{v}\right\rangle=(\mathbf{z},\mathbf{v})_{X}$ for all $\mathbf{v}\in X$, and $\|\mathbf{z}\|_{X}=\|\rho\dot{\mathbf{u}}\|_{X^{\prime}}$ holds. Thus we obtain

[TABLE]

Therefore, using the uniform continuity of $a(t;\cdot,\cdot)$, we get

[TABLE]

This establishes control of $\|\rho\dot{\mathbf{u}}\|_{X^{\prime}}$ at the expense of the $X$-norm, which is controlled in (37). Therefore, we make the ansatz $\mathbf{v}=\mathbf{z}+\delta\mathbf{u}\in X$ for some sufficiently large parameter $\delta\geq 1$. We have the estimate

[TABLE]

From (37) and (38) we conclude

[TABLE]

Taking $\delta:=\frac{1}{2\gamma}(1+\Gamma^{2})\geq 1$, we get

[TABLE]

This completes the proof. ∎

Lemma 6**.**

Assume that for all $t\in I$ the bilinear form $a(t;\cdot,\cdot)$ is symmetric on $X$. If $b(\mathbf{u},\mathbf{v})=0$ holds for all $\mathbf{u}\in V^{0}$, then $\mathbf{v}=0$.

Proof.

Take $\mathbf{v}\in X$ such that

[TABLE]

From Corollary 4 with $F(\mathbf{w}):=\int_{0}^{T}\Gamma(\mathbf{v}(t),\mathbf{w}(t))_{1,\Omega}\,dt$, $\mathbf{w}\in X$, it follows that there exists a unique $\mathbf{z}\in V^{0}$ such that

[TABLE]

We take $\mathbf{w}=\mathbf{z}$ in (41), and use (20), $\mathbf{z}(0)=0$. We get

[TABLE]

Hence, $\|\mathbf{z}\|_{X}\leq\|\mathbf{v}\|_{X}$ holds. Using (40) and taking $\mathbf{w}=\mathbf{v}$ in (41) we obtain:

[TABLE]

We define

[TABLE]

If $S$ has measure 0, then (42) shows that $\mathbf{v}=0$. Thus it suffices to prove that $|S|>0$ leads to a contradiction. Assume that $|S|>0$ holds. We apply, for $t\in S$, the Cauchy-Schwarz inequality to the symmetric positive semi-definite bilinear form $\Gamma(\cdot,\cdot)_{1,\Omega}-a(t,\cdot,\cdot)$ and use the ellipticity property (34):

[TABLE]

which results in a contradiction. Hence $\mathbf{v}=0$ must hold. ∎

As a direct consequence of the preceding two lemmas and the continuity of $b(\cdot,\cdot)$ on $V^{0}\times X$ we obtain the following main well-posedness result.

Theorem 7**.**

Assume that $a(t;\cdot,\cdot)$ satisfies (34)-(35) and is symmetric. For any $F\in X^{\prime}$, the problem (33) has a unique solution $\mathbf{u}\in V^{0}$. This solution satisfies the a-priori estimate

[TABLE]

We can apply this result to the time dependent bilinear form used in the weak formulation of our original problem, cf. (9). Hence, we obtain the following result, which shows well-posedness of the problem (15) with $W$ replaced by the (possibly) smaller subspace $V$.

Corollary 8**.**

For $F\in X^{\prime}$ there exists a unique $\mathbf{u}\in V^{0}$ such that

[TABLE]

Furthermore $\|\mathbf{u}\|_{V}\leq c\|F\|_{X^{\prime}}$ holds with a constant $c$ independent of $F$.

We derive a generalization of Theorem 7 in which the condition that $a(t;\cdot,\cdot)$ is symmetric is not needed.

Theorem 9**.**

Assume that $a(t;\cdot,\cdot)$ satisfies (34)-(35). For any $F\in X^{\prime}$, the problem (33) has a unique solution $\mathbf{u}\in V^{0}$. This solution satisfies the a-priori estimate

[TABLE]

Proof.

Recall the Neumann series result, that if $A\in\mathcal{L}(X,X)$ for some Banach space $X$ and $\|A\|_{\mathcal{L}(X,X)}<1$, then $I+A$ is an isomorphism on $X$ and $(I+A)^{-1}\in\mathcal{L}(X,X)$ (see [4, §5.7]). We introduce some notation. Define the anti-symmetric part of $a(t;\cdot,\cdot)$:

[TABLE]

We split the problem into a problem that we have treated in Theorem 7: $B\mathbf{u}=b(\mathbf{u},\cdot)-\int_{0}^{T}c(t;\mathbf{u},\cdot)\in X^{\prime}$ and a anti-symmetric part $C\mathbf{u}=\int_{0}^{T}c(t;\mathbf{u},\cdot)\in X^{\prime}$, hence (33) has the operator representation $(B+C)\mathbf{u}=F$. For $k\in\mathbb{N}$ we set $C_{k}:=\frac{1}{k}C$. Take $N\in\mathbb{N}$ sufficiently large such that $\|C_{N}\|_{\mathcal{L}(X,X^{\prime})}\leq\frac{\gamma}{2}$ holds. We prove the following statement by induction: for $k\in\mathbb{N}$ the operator $B+kC_{N}\in\mathcal{L}(V^{0},X^{\prime})$ is an isomorphism and $\|(B+kC_{N})^{-1}\|_{\mathcal{L}(X^{\prime},X)}\leq\frac{1}{\gamma}$ holds.

For $k=0$ we can apply Theorem 7, because the symmetric part of $a(t;\cdot,\cdot)$ also satisfies (34)-(35). Hence $B\in\mathcal{L}(V^{0},X^{\prime})$ is an isomorphism. The estimate $\|B^{-1}\|_{\mathcal{L}(X^{\prime},X)}\leq\frac{1}{\gamma}$ follows from (37). We now treat the induction step. Assume that for given $k$ the statement holds. This implies

[TABLE]

and thus by the Neumann series result we get that $I+C_{N}(B+kC_{N})^{-1}\in\mathcal{L}(X^{\prime},X^{\prime})$ is an isomorphism on $X^{\prime}$. Using this, the induction hypothesis and the relation

[TABLE]

it follows that $B+(k+1)C_{N}\in\mathcal{L}(V^{0},X^{\prime})$ is an isomorphism. Using the antisymmetry property of $C$, i.e., $\left\langle C_{N}\mathbf{u},\mathbf{u}\right\rangle=0$ and the ellipticity of $B$, cf. (37), we get for arbitrary $\mathbf{u}\in V^{0}$:

[TABLE]

hence, $\|(B+(k+1)C_{N})^{-1}\|_{\mathcal{L}(X^{\prime},X)}\leq\frac{1}{\gamma}$, which completes the induction. Taking $k=N$ we obtain that $B+C\in\mathcal{L}(V^{0},X^{\prime})$ is an isomorphism. From (36) and $b(\mathbf{u},\cdot)=F$ we get

[TABLE]

which completes the proof. ∎

4 Space-time variational formulation in a broken space

In view of the fact that we want to use a DG method in time, we will now study a time-discontinuous weak formulation. Let $N\in\mathbb{N}$, let $0=t_{0}<\dots<t_{N}=T$ and let $I_{n}=(t_{n-1},t_{n})$ for $n=1,\dots,N$. For $\mathbf{v}\in X$ we define $\mathbf{v}_{n}:=\mathbf{v}|_{I_{n}}\in X_{n}:=L^{2}(I_{n};\mathcal{V})\subset X$, $1\leq n\leq N$. Furthermore

[TABLE]

We define jumps at $t_{n}$ in the usual way. For $\mathbf{u}\in V^{b}$:

[TABLE]

Note that the superscript $n$ denotes an evaluation at $t=t_{n}$, whereas $\mathbf{v}_{n}$ denotes the restriction of $\mathbf{v}$ to $I_{n}$. Note that

[TABLE]

For $\mathbf{u}_{n}\in V_{n}$ we define

[TABLE]

Hence $\rho\dot{\mathbf{u}}_{n}\in X_{n}^{\prime}$.

Remark 4.1**.**

On $Q_{n}:=I_{n}\times\Omega$ we can define a set of smooth functions analogous to (11) by

[TABLE]

which is dense in $X_{n}$. Thus we get

[TABLE]

Hence $\rho\dot{\mathbf{u}}_{n}$ is the same weak material derivative as in Section 2, with $I$ replaced by $I_{n}$. Thus we have analogous results, e.g. as in (13). In particular, for $\mathbf{u}_{n}\in C^{1}(\bar{Q}_{n})^{d}\cap X$ we have

[TABLE]

We also have

[TABLE]

Using $X^{\prime}=L^{2}(I;\mathcal{V}^{\prime})=\oplus_{n=1}^{N}L^{2}(I_{n};\mathcal{V}^{\prime})$ we get

[TABLE]

A broken weak time derivative is defined in the canonical way:

[TABLE]

Hence,

[TABLE]

For controlling the jumps at the interval end points we introduce the usual discontinuous Galerkin bilinear form

[TABLE]

with $\mathbf{z}=(\mathbf{z}^{0},\ldots,\mathbf{z}^{N-1})\in(L^{2}(\Omega)^{d})^{N}$.

As test space in the weak formulation below we use $Y:=X\times H^{N}=\oplus_{n=1}^{N}(X_{n}\times H)$, where $H:=\overline{\mathcal{V}}^{L^{2}}$. We consider the following weak formulation: given $F\in X^{\prime},G\in H^{\prime}$ determine $\mathbf{u}\in V^{b}$ such that

[TABLE]

Note that with $b(\cdot,\cdot)$ as in (33) we have

[TABLE]

In the next theorem we derive equivalence results between different variational formulations.

Theorem 10**.**

*Let the assumptions as in Theorem 9 be satisfied. For $F\in X^{\prime}$ let $\mathbf{u}\in V^{0}$ be the unique solution of (33). Then $\mathbf{u}$ is also the unique solution of each of the following variational problems:

1. The problem (51) with $G=0$.

2. Determine $\mathbf{u}\in V^{b}$ such that*

[TABLE]

with $\rho\mathbf{v}_{+}:=(\rho(t_{0})\mathbf{v}_{+}^{0},\ldots,\rho(t_{N-1})\mathbf{v}_{+}^{N-1})$.

Proof.

Let $\mathbf{u}\in V^{0}$ be the unique solution of (33). Then $d(\mathbf{u},\mathbf{z})=0$ for all $\mathbf{z}\in H^{N}$ and, cf. (49), $\left\langle\rho\dot{\mathbf{u}},\mathbf{v}\right\rangle_{b}=\left\langle\rho\dot{\mathbf{u}},\mathbf{v}\right\rangle$. Hence, $\mathbf{u}\in V^{0}\subset V^{b}$ solves (51) with $G=0$. Let $\mathbf{u}\in V^{b}$ be a solution of (51) with $G=0$. Taking $\mathbf{v}=0$ we get $d(\mathbf{u},\mathbf{z})=0$ for all $\mathbf{z}\in H^{N}$. This implies $[\mathbf{u}]^{n}=0,~{}0\leq n\leq N-1,$ and thus, cf. (45), $\mathbf{u}\in V^{0}$. Take $\mathbf{z}=0$ and using $\left\langle\rho\dot{\mathbf{u}},\mathbf{v}\right\rangle_{b}=\left\langle\rho\dot{\mathbf{u}},\mathbf{v}\right\rangle$ we conclude that $\mathbf{u}$ solves (33). Hence, the unique solution $\mathbf{u}\in V^{0}$ of (33) is also the unique solution of (51) with $G=0$.

Let $\mathbf{u}\in V^{0}$ be the unique solution of (33), which is also the unique solution of (51) with $G=0$. Taking arbitrary $\mathbf{v}\in V^{b}\subset X$ and $\mathbf{z}\in H^{N}$ such that $d(\cdot,\mathbf{z})=d(\cdot,\rho\mathbf{v}_{+})\in(H^{N})^{\prime}$ in (51) it follows that $\mathbf{u}$ is a solution of (53). Let $\mathbf{u}\in V^{b}$ be a solution of (53). The space $\{\,(\mathbf{v},d(\cdot,\rho\mathbf{v}_{+}))~{}|~{}\mathbf{v}\in V^{b}\,\}$ is dense in $X\times(H^{N})^{\prime}$. Note that $\mathbf{v}\to\left\langle\rho\dot{\mathbf{u}},\mathbf{v}\right\rangle_{b}$, $\mathbf{v}\to\int_{0}^{T}a(t;\mathbf{u}(t),\mathbf{v}(t))\,dt$ are continuous functionals on $X$ and $\mathbf{z}\to d(\mathbf{u},\mathbf{z})$ is continuous on $H^{N}$. Using a density argument it follows that $\mathbf{u}$ solves (51) with $G=0$. Hence, the unique solution $\mathbf{u}\in V^{0}$ of (33) is also the unique solution of (53). ∎

The factor $\rho$ in the coupling term $d(\cdot,\cdot)$ in (53) is not essential. It is introduced to obtain a natural scaling, namely one that corresponds to the scaling with $\rho$ in the weak time derivative. Note that in (53) the initial condition $\mathbf{u}(0)=0$ is treated in a weak sense (applies also to $\mathbf{u}(0)=\mathbf{u}^{0}\neq 0$).

Remark 4.2**.**

In Theorem 10 we (only) show that the problem (51) with $G=0$ has a unique solution. For the variational problem (51) a more general well-posedness result can be derived, namely that the bilinear form $B(\cdot,\cdot)$ defines a homeomorphism $V^{b}\to Y^{\prime}$, with norms

[TABLE]

Note that $(V^{b},\|\cdot\|_{V^{b}})$ and $(Y,\|\cdot\|_{Y})$ are Banach spaces. Continuity of the bilinear form $B(\cdot,\cdot)$ on $V^{b}\times Y$ is easy to show. Furthermore, provided $a(t;\cdot,\cdot)$ satisfies (34)-(35), it can be shown that the BNB infsup conditions are satisfied. We do not include a proof in this paper. Given these results one obtains that under the above assumptions on $a(t;\cdot,\cdot)$, for any $F\in X$, $G\in(H^{N})^{\prime}$ the problem (51) has a unique solution $\mathbf{u}\in V^{b}$ and the estimate

[TABLE]

holds with a constant $c$ depending only on $\gamma,\,\Gamma$ from (34)-(35).

Remark 4.3**.**

From the results above it follows that if the assumptions as in Theorem 9 are satisfied, then the weak formulation (53) is a well-posed variational formulation of the original Stokes problem (5). This variational formulation, in which the same trial and test space $V^{b}$ is used, can be reformulated using a time stepping procedure. The unique solution $\mathbf{u}\in V^{b}$ of (53) can be decomposed as $\mathbf{u}=(\mathbf{u}_{1},\ldots,\mathbf{u}_{N})$, with $\mathbf{u}_{n}\in V_{n}$, and the solution of (53) is also the unique solution of the problem: for $n=1,\ldots,N$, determine $\mathbf{u}_{n}\in V_{n}$ such that

[TABLE]

This is the usual form of a discontinuous Galerkin method for parabolic PDEs, cf. [42]. If $\mathbf{u}_{n}$ has sufficient smoothness, e.g. $\mathbf{u}_{n}\in C^{1}(\bar{Q}_{n})\cap V_{n}$, the weak material derivative reduces to the usual strong one: $\left\langle\rho\dot{\mathbf{u}}_{n},\mathbf{v}_{n}\right\rangle_{n}=\int_{I_{n}}(\rho\dot{\mathbf{u}}_{n}(t),\mathbf{v}_{n}(t))_{L^{2}}\,dt$. This formulation is a reasonable starting point for a Galerkin finite element discretization in which the space $V^{b}$ is replaced by a (space-time) finite element subspace. This, however, requires exactly divergence free finite element functions. Recently, such divergence free finite element methods have been further developed using techniques from finite element exterior calculus, e.g. [19]. Most finite element methods, however, treat the divergence constraint by means of a pressure Lagrange multiplier, see [26]. Therefore, in Section 5 we introduce a variant of the weak formulation (53) that involves the pressure Lagrange multiplier to satisfy the divergence free constraint.

5 Existence of a pressure Lagrange multiplier in $L^{2}(Q)$

In this section we reconsider the problem (33), for which a well-posedness result is given in Theorem 9. In the variational problem (33), both in the solution space $V^{0}$ and test space $X$ we restrict to functions $\mathbf{v}$ which satisfy $\mathop{\rm div}\mathbf{v}=0$ on $\Omega$. In this section we derive a formulation in which we eliminate this condition from the trial and test space and instead introduce the pressure Lagrange multiplier for satisfying the divergence free constraint. For this one typically needs additional regularity properties of the solution $\mathbf{u}$ of (33), cf. Section 6.2.1 in [17]. The regularity property that we require in Theorem 11 below will be discussed in Remark 5.1. We use an analysis along the same lines as given for a time dependent Stokes problem with constant coefficients (density and viscosity) in [17].

We first introduce a space-time variant of de Rham’s theorem. Let $\nabla:\,L_{0}^{2}(\Omega)\to H^{-1}(\Omega)^{d}$ be the weak gradient. A standard application of de Rham’s theorem, e.g., Corollary 2.4. in [22], yields:

[TABLE]

where $L^{2}_{0}(\Omega)=\{\,p\in L^{2}(\Omega)~{}|~{}\int_{\Omega}p=0\,\}$. We define $\nabla_{\otimes}={\rm id}\otimes\nabla:\,L^{2}(I;L_{0}^{2}(\Omega))=L^{2}(I)\otimes L_{0}^{2}(\Omega)\rightarrow L^{2}(I;H^{-1}(\Omega)^{d})=L^{2}(I)\otimes H^{-1}(\Omega)^{d}$ in the usual way, i.e., for $g\in L^{2}(I;L_{0}^{2}(\Omega))$, $g(t)=\sum_{i=0}^{\infty}\alpha_{i}(t)\phi_{i}$ with $\alpha_{i}\in L^{2}(I)$, $\phi_{i}\in L_{0}^{2}(\Omega)$ we define $(\nabla_{\otimes}g)(t):=\sum_{i=1}^{\infty}\alpha_{i}(t)\nabla\phi_{i}\in H^{-1}(\Omega)^{d}$. From (55) it follows that

[TABLE]

Furthermore, for $g\in L^{2}(I;L^{2}_{0}(\Omega)),\,\mathbf{v}\in L^{2}(I;H_{0}^{1}(\Omega)^{d})$ we have

[TABLE]

We introduce notation for spaces. Recall $U=\{\,\mathbf{v}\in X~{}|~{}\frac{\partial\mathbf{v}}{\partial t}\in L^{2}(I;L^{2}(\Omega)^{d})\,\}$, cf. (16). We define

[TABLE]

Hence, $U=\{\,\mathbf{v}\in\tilde{U}~{}|~{}\mathop{\rm div}\mathbf{v}(t)=0\quad\text{a.e. for}~{}t\in I\,\}$. Clearly, opposite to $U$ and $V=\overline{U}^{\|\cdot\|_{W}}$, the space $\tilde{U}$ does not involve the divergence free constraint. Below we use this space as trial space and $L^{2}(I;H_{0}^{1}(\Omega)^{d})$ (instead of $X$) as test space for the velocity. In order to do this we assume that the bilinear from $a(t;\cdot,\cdot)$ is not only defined on $\mathcal{V}\times\mathcal{V}$ but on $\mathcal{V}\times H^{1}_{0}(\Omega)^{d}$ and satisfies

[TABLE]

for a positive constant $\tilde{\Gamma}$, independent of $\mathbf{v},\tilde{\mathbf{v}}$.

Theorem 11**.**

Let the assumptions of Theorem 9 hold and assume that, for given $F\in L^{2}(I;H^{-1}(\Omega)^{d})\subset X^{\prime}$, the unique solution $\mathbf{u}$ of (33) has smoothness $\frac{\partial\mathbf{u}}{\partial t}\in L^{2}(I;L^{2}(\Omega)^{d})$, i.e., $\mathbf{u}\in U$. Assume that the bilinear form $a(t;\cdot,\cdot)$ is defined on $\mathcal{V}\times H^{1}_{0}(\Omega)^{d}$ and satisfies (58). Consider the following problem: determine $\mathbf{u}\in\tilde{U}$, $p\in L^{2}(I;L_{0}^{2}(\Omega))$ such that

[TABLE]

for all $\mathbf{v}\in L^{2}(I;H_{0}^{1}(\Omega)^{d})$, $q\in L^{2}(I;L^{2}_{0}(\Omega))$. This problem has a unique solution $(\mathbf{u},p)$ and $\mathbf{u}$ equals the unique solution of (33).

Proof.

Let $\mathbf{u}$ be the unique solution of (33), which by assumption has smoothness $\frac{\partial\mathbf{u}}{\partial t}\in L^{2}(I;L^{2}(\Omega)^{d})$. Hence, $\left\langle\rho\dot{\mathbf{u}},\mathbf{v}\right\rangle=(\rho\dot{\mathbf{u}},\mathbf{v})_{L^{2}}$ for all $\mathbf{v}\in X$ and (59) holds for all $\mathbf{v}\in X$. Define

[TABLE]

Then $l\in L^{2}(I;\mathcal{V}^{0})$. From (56) it follows that there exists a unique $p\in L^{2}(I;L_{0}^{2}(\Omega))$ such that

[TABLE]

Combining this with (57) we conclude that $(\mathbf{u},p)$ satisfies (59) for all $\mathbf{v}\in L^{2}(I;H_{0}^{1}(\Omega)^{d})$. Furthermore, $\mathbf{u}$ trivially satisfies (60), due to $\mathop{\rm div}\mathbf{u}(t)=0$. Hence, the unique solution $\mathbf{u}$ of (33) and the corresponding unique $p\in L^{2}(I;L_{0}^{2}(\Omega))$ solve (59)-(60) for all $\mathbf{v}\in L^{2}(I;H_{0}^{1}(\Omega)^{d})$, $q\in L^{2}(I;L^{2}_{0}(\Omega))$.

We now consider the other direction. Let $(\mathbf{u},p)\in\tilde{U}\times L^{2}(I;L_{0}^{2}(\Omega))$