Analysis of a time-discrete scheme for the Navier-Stokes/Allen-Cahn model

TL;DR

This paper develops and analyzes a stable, implicit time-discrete scheme for simulating two-fluid flows with different densities and viscosities using the Navier-Stokes/Allen-Cahn model, including convergence and numerical validation.

Contribution

It introduces a novel stabilized fixed point algorithm for the coupled Navier-Stokes/Allen-Cahn system, ensuring existence, uniqueness, and maximum principle in a semi-discrete setting.

Findings

Proved well-posedness and convergence of the scheme.

Established maximum principle for the numerical solution.

Numerical experiments demonstrate the effectiveness of the iterative process.

Abstract

This paper address the approximation of the dynamic of two fluids with non matching densities and viscosities modeled by the Allen-Cahn equation coupled with the time dependent Navier-Stokes equations. Existence, uniqueness and a maximum principle are obtained for a totally implicit semi-discrete in time formulation. These results are based on an original stabilized fixed point algorithm for which well posedness and convergence is analyzed. Numerical experiments are performed to show the influence of the iterative process.

Tables

Table 1. Table 1: Precision ( L 2 superscript 𝐿 2 L^{2} error) and computational cost (total number of algebraic system solved) with respect to the three methods (fully implicit Newton method (FIN), fully implicit Picard method (FIP) and strongly coupled explicit method (SCE)).

Method	$\beta$	$\Delta t$	N	# of solves	$E_u$	$E_{\phi }$
FIN	0	1/1 300	10	40	6.74601e-07	9.89655e-11
FIN	9/8	1/1 300	10	90	6.74601e-07	3.02437e-10
FIP	0	1/1 300	10	97	6.74601e-07	2.47127e-10
FIP	2	1/1 300	10	102	6.74601e-07	3.32047e-10
SCE	0	1/1 300	10	31	0.00176404	0.031862
SCE	0	1/13 000	100	295	0.00020879	0.003257

Equations

$$\phi (t,\mathbf{x})=\{\begin{array}{ll}-1& \mathbf{x}\text{ in fluid a}\\ 1& \mathbf{x}\text{ in fluid b}\end{array}$$

$$\rho (\phi )=\bar{\rho }+\frac{\delta \rho }{2}\phi ,\mu (\phi )=\bar{\mu }+\frac{\delta \mu }{2}\phi$$

$$\begin{array}[]{c}\displaystyle\bar{\rho}=\frac{\rho_{a}+\rho_{b}}{2},\ \ \delta\rho=\rho_{a}-\rho_{b},\ \ \rho_{I}=\min\{\rho_{a},\rho_{b}\},\ \ \rho_{S}=\max\{\rho_{a},\rho_{b}\}.\vspace*{8pt}\\ \displaystyle\bar{\mu}=\frac{\mu_{a}+\mu_{b}}{2},\ \ \delta\mu=\mu_{a}-\mu_{b},\ \ \mu_{I}=\min\{\mu_{a},\mu_{b}\},\ \ \mu_{S}=\max\{\mu_{a},\mu_{b}\}.\end{array}$$

$${\partial }_t\phi +\left(\mathbf{u}\cdot \nabla \right)\phi =\gamma \left(\Delta \phi -\frac{d}{d\phi }F(\phi )\right)$$

$$F(\phi )=\frac{1}{4{\eta }^2}({\phi }^2-1{)}^2\frac{dF(\phi )}{d\phi }=f(\phi )=\frac{1}{{\eta }^2}\phi \left({\phi }^2-1\right)$$

$$\sqrt{\rho }{\partial }_t\left(\sqrt{\rho }\mathbf{u}\right)+\rho \left(\mathbf{u}\cdot \nabla \right)\mathbf{u}+\frac{1}{2}\nabla \cdot \left(\rho \mathbf{u}\right)\mathbf{u}-\nabla \cdot \left(\mu \mathbf{D}(\mathbf{u})\right)+\nabla p$$

$$=-\sigma \nabla \cdot \mathbf{T}(\phi )+\mathbf{G}(\phi )$$

$$\mathbf{D}(\mathbf{u})=\frac{1}{2}(\nabla \mathbf{u}+\nabla {\mathbf{u}}^t).$$

$$\mathbf{T}(\phi )=\sigma \nabla \phi \otimes \nabla \phi ,$$

$$\nabla \cdot \mathbf{T}(\phi )=\sigma \nabla \cdot \left(\nabla \phi \otimes \nabla \phi \right)$$

$$=\frac{\sigma }{\gamma }\left(\partial {\phi }_t+\mathbf{u}\cdot \nabla \phi \right)\nabla \phi -\nabla \hat{p},$$

$$\{\begin{array}{rl}& \sqrt{\rho }{\partial }_t\left(\sqrt{\rho }\mathbf{u}\right)+\rho \left(\mathbf{u}\cdot \nabla \right)\mathbf{u}+\frac{1}{2}\nabla \cdot \left(\rho \mathbf{u}\right)\mathbf{u}-\nabla \cdot \left(\mu \mathbf{D}(\mathbf{u})\right)+\nabla p\\ & =-\frac{\sigma }{\gamma }\left({\partial }_t\phi +\mathbf{u}\cdot \nabla \phi \right)\nabla \phi +\mathbf{G}(\phi ),\\ & \nabla \cdot \mathbf{u}=0,\\ & {\partial }_t\phi +\left(\mathbf{u}\cdot \nabla \right)\phi =\gamma \left(\Delta \phi -\frac{1}{{\eta }^2}\phi \left({\phi }^2-1\right)\right).\end{array}$$

$$L_0^2(\Omega )=\left\{q\in L^2(\Omega );{\int }_{\Omega }q(\mathbf{x})d\mathbf{x}=0\right\}.$$

$$W^{m,p}(\Omega )=\left\{u\in L^p(\Omega ):D^{\alpha }u\in L^p(\Omega )\forall |\alpha |\le m\right\}$$

$$L^p(0,T;E)=\left\{u:(0,T)\mapsto E:{\left({\int }_0^T\Vert u(\tau ){\Vert }_E^{p}d\tau \right)}^{1/p}<\infty \right\}$$

$$r\in \left[q,\frac{qd}{d-q}\right],\text{mbox}ifq\in [1,d),\text{mbox}andr\in [q,+\infty )\text{mbox}ifq\ge d.$$

$$\Vert u{\Vert }_{L^r}\le C\Vert u{\Vert }_{L^q}^{1-s}\Vert u{\Vert }_{1,q}^s\text{mbox}wheres=\frac{d(r-q)}{rq}.$$

$$-\underset{q\in L_0^2(\Omega )}{\inf }\underset{\mathbf{v}\in {\mathbf{H}}_0^1(\Omega )}{\sup }\frac{\left(q,\nabla \cdot \mathbf{v}\right)}{\Vert q\Vert \Vert \mathbf{v}{\Vert }_1}\ge \kappa .$$

$$\begin{array}{rl}& {\int }_{\Omega }(\mu ({\phi }_0)\mathbf{D}({\mathbf{w}}_0)-\mu ({\phi }_1)\mathbf{D}({\mathbf{w}}_1)):\mathbf{D}\left({\mathbf{w}}_0-{\mathbf{w}}_1\right)=\\ & {\int }_{\Omega }\mu ({\phi }_0){\left|\mathbf{D}\left({\mathbf{w}}_0-{\mathbf{w}}_1\right)\right|}^2+\frac{\delta \mu }{2}{\int }_{\Omega }\left({\phi }_0-{\phi }_1\right)\mathbf{D}({\mathbf{w}}_1):\mathbf{D}\left({\mathbf{w}}_0-{\mathbf{w}}_1\right).\end{array}$$

$$\begin{array}{r}{\int }_{\Omega }\rho ({\phi }_0)({\mathbf{w}}_0\cdot \nabla ){\mathbf{w}}_1\cdot {\mathbf{w}}_1=-\frac{1}{2}{\int }_{\Omega }\nabla \cdot (\rho ({\phi }_0){\mathbf{w}}_0){\mathbf{w}}_1\cdot {\mathbf{w}}_1\end{array}$$

$$\begin{array}{r}{\int }_{\Omega }({\mathbf{w}}_0\cdot \nabla {\phi }_0){\phi }_1=-{\int }_{\Omega }({\mathbf{w}}_0\cdot \nabla {\phi }_1){\phi }_0,{\int }_{\Omega }({\mathbf{w}}_0\cdot \nabla {\phi }_0){\phi }_0=0\end{array}$$

$$\begin{array}{rl}{\int }_{\Omega }({\mathbf{w}}_0\cdot \nabla {\phi }_0-{\mathbf{w}}_1\cdot \nabla {\phi }_1)\left({\phi }_0-{\phi }_1\right)& ={\int }_{\Omega }\left({\mathbf{w}}_0-{\mathbf{w}}_1\right)\cdot \nabla {\phi }_1\left({\phi }_0-{\phi }_1\right)\\ & =-{\int }_{\Omega }\left({\mathbf{w}}_0-{\mathbf{w}}_1\right)\cdot \nabla \left({\phi }_0-{\phi }_1\right){\phi }_1.\end{array}$$

$$a^{2}b-c^{2}d$$

$$ab-cd$$

$$\underset{{\mu }_I\to \infty }{\lim }\mathcal{K}=0$$

$$\mathcal{K}=\mathcal{C}+\mathcal{K}.$$

$$\left\{\begin{array}[]{rcl}-a\Delta\varphi+{\boldsymbol{b}}\cdot\nabla\varphi+c\varphi&=&f\qquad\mbox{in}\qquad\Omega\\ \partial_{\boldsymbol{n}}\varphi&=&0\qquad\mbox{on}\qquad\partial\Omega\end{array}\right.$$

$$\Vert \phi {\Vert }_{2+r}$$

$$\left(\mathbf{u},{\partial }_{\mathbf{n}}\phi \right)$$

$${\left(\mathbf{u},\phi \right)}_{t=0}$$

$$\mathbf{G}(0)=0,\Vert {\mathbf{G}}^{\prime }{\Vert }_{L^{\infty }}\le 1\text{mbox}and\forall s\in \mathbb{R},|\mathbf{G}(s)|\le s.$$

$${\mathbf{u}}^0\in {\mathbf{H}}_0^1(\Omega )\cap {\mathbf{H}}^2(\Omega )\text{mbox}and{\phi }^0\in H^3(\Omega )$$

Full text

1 Introduction

In the broad spectrum of approaches to model multi-fluid (multi-phase) flow and the capture of its interfacial behaviour, the phase-field approach is quite attractive as it is a physically motivated model based on the competition between the different species/phases. We refer the readers to [1, 2, 3] for a review of some of the most common approaches (both theoretically and numerically).

This work is part of an effort to produce a new numerical algorithm for the approximation of the complex behaviour of a binary mixture of fluids. More precisely, we are interested in the use of the Allen-Cahn model for the description of the dynamic of the mixture of two non miscible fluids coupled to the unsteady Navier-Stokes equations describing fluids displacements. This leads to what is known as an unsteady Navier-Stokes-Allen-Cahn (NS-AC) model.

We consider an unsteady flow of two immiscible incompressible fluids of nonnegative constants densities $\rho_{a},\rho_{b}>0$ and viscosities $\nu_{a},\nu_{b}>0$. Introducing the scalar function $\varphi$ (called a phase field)

[TABLE]

representing the difference of volume fraction of both fluids, we define the density and viscosity of the fluids mixture as

[TABLE]

with

[TABLE]

The behavior of $\varphi$ is characterized using the Allen-Cahn model which describes the separation of a mixture composed of two species (see [4, 3, 2] concerning the physical details of the model). Although this model raises questions (mass conservation is violated), we neglect modifications to the system in order to recuperate the mass conservation as discussed in [5, 6] for instance. Then the phase field satisfy

[TABLE]

where the nonnegative constant $\gamma>0$ describes the mobility coefficient, $F$ is the free energy density of the mixture, (a double well potential), defined as

[TABLE]

and $0<\eta$ is a parameter related to the thickness of the transition zone between both fluids. Assuming this thickness is small with respect to the other dimensions, in fact we assume that $\eta^{2}\ll\gamma$, the interface between both fluids can be described as the zero level of $\varphi$.

The incompressible nature of the fluids imposes a divergence free velocity ${\boldsymbol{u}}$, and following [7] (for example) the momentum equation of the mixture is

[TABLE]

where $p$ is the fluids pressure and $\boldsymbol{D}({\boldsymbol{u}})$ is the shear rate defined by

[TABLE]

As for the forces involved, ${\boldsymbol{G}}$ represents external volumetric force (such as gravity), and the extraneous elastic stress induced by the surface tension is defined as

[TABLE]

where $\sigma$ represents the density of the mixing energy. Before going any further, based on the potential of the free energy (4) and (3) we can replace the surface tension term

[TABLE]

where $\hat{p}=\displaystyle\sigma(\frac{1}{2}|\nabla\varphi|^{2}+F(\varphi))$ is of the same nature as a pressure. As in [7], we introduce a new pressure, $p+\hat{p}$, still denoted $p$. Then we can rewrite the momentum (5), and the system describing the behaviour of the mixture is

[TABLE]

As this system is nonlinear and strongly coupled, the numerical approximation of solutions is challenging. Existence of solutions for this system is known (see [8] and the references therein), and semi-explicit numerical strategies are documented (see, for example, the approach proposed in [7]). Such explicit schemes, if they are generally perceived as efficient (at least from a computational standpoint), frequently lacks in accuracy and require important constraint on the time step to recover the loss in precision.

From empirical observations, for a given precision, implicit scheme allows larger time step when compared to explicit methods. This makes implicit methods attractive as they frequently lead to less computational effort for a fixed time horizon (Table 1 for example). Aiming for a more accurate while flexible approach, this paper addresses some theoretical aspects related to the numerical approximation of (6) using implicit time discretization.

Based on an original fixed point (the iterates are defined by (23)–(25)), this paper presents the existence and uniqueness of the solution for a semi-discrete formulation of (6). The result is based on the study of the well-posedness (Theorem 4.1), uniform boundness (Theorems 4.2) and the strong convergence (Theorem 4.4) of the fixed point problem. Finally, as the fixed point problem produces bounded phase fields, a maximum principle is obtained for the limit problem (19).

2 Preliminaries

2.1 Spaces

In what follows, $\Omega$ is a bounded open convex of class $C^{2,1}$ of $\mathbb{R}^{d},~{}d=$ 2 or 3, and $\partial\Omega$ denotes the boundary of $\Omega$. $L^{p}(\Omega)$ represent the usual set of $p-th$ power measurable functions, and $\left(L^{p}(\Omega)\right)^{d}={\boldsymbol{L}}^{p}(\Omega)$. The scalar product defined on $L^{2}(\Omega)$ or ${\boldsymbol{L}}^{2}(\Omega)$ is denoted (without distinction) by $\left(\cdot,\cdot\right)$ and its norm $\|\cdot\|$. The space $L^{2}_{0}(\Omega)$ is defined as

[TABLE]

The Sobolev spaces, denoted $W^{m,p}(\Omega)$ (we denote ${\boldsymbol{W}}^{m,p}(\Omega)=\left(W^{m,p}(\Omega)\right)^{d}$), with $p\in[1,+\infty)$ and $m$ integer, is defined as

[TABLE]

where $\alpha$ is a multi-index in $\mathbb{N}^{d}$. These spaces are equipped with the norm $\|\cdot\|_{m,p}$ and semi-norm $|\cdot|_{m,p}$.

The spaces $W^{s,2}(\Omega)$ and ${\boldsymbol{W}}^{s,2}(\Omega)$, $s\in\mathbb{R}$, are denoted $H^{s}(\Omega)$ and ${\boldsymbol{H}}^{s}(\Omega)$ respectively. Their norm are denoted $\|\cdot\|_{s}$ and semi-norm $|\cdot|_{s}$. Without distinction for the dimension, we denote the duality pairing between $H^{1}_{0}(\Omega)$ and its dual $H^{-1}(\Omega)$ (or between ${\boldsymbol{H}}^{1}_{0}(\Omega)$ and ${\boldsymbol{H}}^{-1}(\Omega)$)) by $\left<\cdot,\cdot\right>$. More generally, for a space $V$ and $V^{\prime}$ its dual, we denote the duality pairing by $\left<\cdot,\cdot\right>_{V^{\prime},V}$.

For a fixed positive real number $T$ (representing the final time) and a separable Banach space $E$ equipped with the norm $\|\cdot\|_{E}$, we denote by $\mathcal{C}^{0}(0,T;E)$ the space of continuous functions from $[0,T]$ with values in $E$. For a positive integer $p$, we introduce the Bochner spaces,

[TABLE]

We refer to [9], [10, Chapter 2] for details concerning those spaces. Let us remind the following interpolation inequality for Sobolev spaces (see [11, Exercise II.3.12] for instance):

Lemma 2.1

Assume $\Omega$ locally Lipschitz and let

[TABLE]

The following inequality holds for all $u\in W^{1,q}(\Omega)$

[TABLE]

The assumption on $\Omega$ are sufficient to verify the inf–-sup or Ladyzhenskaya-–Babus̆ka–-Brezzi (LBB) condition (see [12, 13, 10]), that is, there exists a positive constant $\kappa>0$ depending only on $\Omega$, such that

[TABLE]

Most proofs will rely on combination of Sobolev embedding, Poincaré, Korn, Poincaré-Wirtinger, Cauchy-Schwarz, Holder and Young inequalities (see [9, 14] for details on these concept). We will denote by SPKCSHY–sequence such combination.

2.2 Equalities and constants

In order to simplify the proofs, we introduce the following equalities, for arbitrary fields ${\boldsymbol{w}}_{0},{\boldsymbol{w}}_{1}\in\mathbb{R}^{d}$ and $\varphi_{0},\varphi_{1}\in\mathbb{R}$

[TABLE]

If $\nabla\cdot{\boldsymbol{w}}_{0}=\nabla\cdot{\boldsymbol{w}}_{1}=0$ then

[TABLE]

[TABLE]

[TABLE]

Some algebraic identities frequently used in this paper, $\forall\,a,b,c,d\in\mathbb{R}$

[TABLE]

In this work, the dependency of constants with respect to physical parameters plays an important role, in particular the behaviour of constants in relation with the parameters $\mu_{I}$ and $\mu_{S}$ defining $\mu$ the viscosity. $C$ will denote generic constants depending at most on $\Omega$ and its boundary (for example, the Poincaré and the Korn constant). $K$ will denote generic constants depending on any physical data. $\mathscr{C}$ will denote generic constants that could depend on the physical parameter and state variables with the exception of $\mu$, $\mu_{I}$ and $\mu_{S}$. $\mathscr{K}$ will denote constant depending on all the physical values with the addition that

[TABLE]

typically $\mathscr{K}=\mathscr{C}\mu_{I}^{-\alpha}$ with $\alpha>0$. Finally $\mathcal{K}$ will denote special constant of the form

[TABLE]

Lastly, when there is no possible confusion, we will make no distinction between the various values of $C$, $\mathscr{C}$, $\mathscr{K}$, $\mathcal{K}$ and $K$ when manipulating expressions.

2.3 Existence and regularity results

Essential element of this work, we recall a lemma for elliptic problem with Neumann boundary conditions (for instance [15, Chap. 3] or [16, 17, 18]). Once again we underline that the assumptions previously made on $\Omega$ and its boundary makes it possible to apply such regularity lemma.

Lemma 2.2

Let $a(\cdot),{\boldsymbol{b}}(\cdot)$ and $c(\cdot)$ be three functions such that $a,c\in L^{\infty}(\Omega)$ and ${\boldsymbol{b}}\in\mathbf{H^{1}_{0}}(\Omega)$. Assume $a({\boldsymbol{x}})\geq a_{0}>0,\,c({\boldsymbol{x}})\geq c_{0}>0$, $\nabla\cdot{\boldsymbol{b}}=0$ and $f\in H^{r}(\Omega),\,r\geq 0$ . Then the following problem

[TABLE]

admits a unique solution which satisfies

[TABLE]

where $C(\Omega,a_{0},c_{0},{\boldsymbol{b}},r)$ is a constant depending on $\Omega,a_{0},\,{\boldsymbol{b}},c_{0}$ and $r$.

Of course, (6) is a well defined system for $({\boldsymbol{u}},p,\varphi)$ once it is completed with boundary and initial conditions. For the sake of simplicity, we consider homogeneous Dirichlet boundary conditions for ${\boldsymbol{u}}$ and homogeneous Neumann condition on $\varphi$

[TABLE]

The usual hypothesis leading to a Boussinesq-type approximation (such as having a ratio of densities near 1) have no effects on the results presented here. Therefore a more general case is considered with the following assumptions.

Assumption 2.1

The function ${\boldsymbol{G}}:\;\mathbb{R}\rightarrow\mathbb{R}^{d}$  is a $C^{1}(\mathbb{R})$ function such that

[TABLE]

The existence of a solution to (6), (16)–(17) is proved provided basic compatibility conditions are satisfied, namely

Assumption 2.2

The data $\left({\boldsymbol{u}}^{0},\varphi^{0}\right)$ satisfies the regularity conditions

[TABLE]

and the compatibility condition

[TABLE]

for some $p_{0}\in H^{1}(\Omega)$.

Finally, we have an existence and uniqueness result for the continuous problem thanks to [19, 20, 8]:

Theorem 2.1

Under Assumptions 2.1–2.2 there exists $T^{*}>0$ and a unique solution $\left({\boldsymbol{u}},p,\varphi\right)$ to the problem (6), (16)–(17) such that

[TABLE]

for some $q,\;d<q<\displaystyle\frac{2d}{d-2}$.

The rest of this work relies on the regularity of $\varphi$ and ${\boldsymbol{u}}$ given by Theorem 2.1. We emphasize that the regularity imposed by Assumption 2.2 is not excessive. The proofs of most of the results presented here hinges on these assumptions (and consequently on (18)). Nevertheless, we do not exclude that results similar to those presented here could be obtained with less regular assumptions on the initial values.

3 A fully implicit semi-discrete formulation

For the sake of simplicity and clarity, a first order approximation of the time derivative was chosen for both the Navier-Stokes and advection-diffusion equations. Different treatment of (6) could be considered, from a backward Euler (BDF1) to a forward Euler scheme. The resulting system contains two important difficulties: it is non linear and strongly coupled in ${\boldsymbol{u}}$ and $\varphi$.

Excluding the time derivative, the terms in (6) can be treated: implicitly (leading to various fixed-point algorithm), semi-explicitly (using known values of the state variables , see [7] for example) or totally explicitly (in which case no system has to be solved). Of course it is only in the implicit case that we have a proper backward Euler scheme with, a priori, all its properties.

We emphasizes that for all these strategies (implicit, semi-explicit and explicit), the resulting system will still be a coupled system. Different approach can be used to deal with the coupling of ${\boldsymbol{u}}$ and $\varphi$, a strongly coupled approach would consist in solving both equation as a system (possibly non linear). A weakly coupled approach would consist in solving each equation once at each time step, expressing the coupling terms in some explicit way (for example by using a Richardson extrapolation).

Any weakly coupled strategy reduces the work load at each time step. However, as for explicit or semi-explicit time scheme, this leads inevitably to conditional stability and certainly imposes conditions on the time step length offsetting any numerical advantages. A strongly coupled approach avoiding these inconvenient is preferable. Such strategies imply to solve a system in $({\boldsymbol{u}},\varphi)$ using a fixed-point approach. An implicit approach for the non linear terms has been retained as well, leading to a fully implicit approach.

Using a uniform time-step $\Delta t>0$ and denoting

[TABLE]

the time discretization of (6) result in a sequence of nonlinear coupled (strongly) problems of the form

[TABLE]

completed with the boundary conditions (16) and $\left({\boldsymbol{u}}^{0},p^{0},\varphi^{0}\right)=\left({\boldsymbol{u}}_{0},p_{0},\varphi_{0}\right).$

Remark 3.1

Through simple modifications of the terms involved in the system (19) different schemes could be considered. Moreover, limiting the fixed point loop needed to solve (19) to one iteration at each time step would produce a weakly coupled scheme.

To finish this section, let us present the main result of this work: the well posed character of (19) and the respect of a maximum principle. More precisely, in the next section we intend to demonstrate the following,

Theorem 3.1

Under the assumptions of Theorem 2.1, with ${\boldsymbol{u}}^{0}$ divergence free. For $\mu_{I}$ sufficiently large there exist a $\tau\in]0,1[$ depending on $\mu_{I}$ such that for all $\Delta t\leq\tau\eta^{2}/\gamma$ the system (19), (16)–(17) admits a unique solution $({\boldsymbol{u}}^{n},p^{n},\varphi^{n})$ in $\boldsymbol{H^{1}_{0}}(\Omega)\cap\boldsymbol{H^{2}}(\Omega)\times L^{2}_{0}(\Omega)\cap H^{1}(\Omega)\times H^{3}(\Omega)$ with

[TABLE]

at each time step $t^{n}=n\Delta t$.

4 Existence, uniqueness and a maximum principle.

In order to prove Theorem 3.1, we introduce, at each time step, a linearized coupled fixed point sequence. We will then show the convergence of this sequence and, through identification of the limit as a solution of (19), we will get the result announced here.

4.1 A coupled iterative scheme

As we consider a fixed point loop, at each time step we introduce ${\boldsymbol{u}}^{n+1}_{k}$ $p^{n+1}_{k}$ and $\varphi^{n+1}_{k}$ the state variables at the $k-th$ iteration of the fixed point loop at time $t^{n+1}$. The index denoting the iteration number of the fixed point loop and the superscript the time step number.

At time $t^{n+1}$, knowing $\left({\boldsymbol{u}}^{n},p^{n},\varphi^{n}\right)$, we consider the following initialization

[TABLE]

In what follows, to simplify the notation, when there is no ambiguity, we will neglect the time step superscript on ${\boldsymbol{u}}$, $p$ and $\varphi$. For solving the nonlinear phase field equation

[TABLE]

two obvious choices are a Picard fixed point or a Newton-type method. Considering the robustness of the method and fact that $f$ is a third degree polynomial in $\varphi^{n+1}$, the Newton approach is an appropriate choice. However, from numerical experiments the Picard fixed point seems to have comparable performance to the Newton like fixed point (see the Numerical tests in section 5).

Replacing $f(\varphi_{k+1})$ by its first order development:

[TABLE]

we get for (20)

[TABLE]

The convergence analysis will rely heavily on the fact that the phase function satisfy, at each fixed point iteration

[TABLE]

To insure such condition, we introduce the supplementary term

[TABLE]

where $\beta$ is an arbitrary nonnegative real constant. Obviously, at convergence this term will be zero. We will show that, provided $\beta$ is above a specific lower bound, $\varphi_{k}$ will satisfy (21) without any additional condition.

Using a simple linearization of the convective term in the momentum equation in (19), we get the following algorithm: at each time step $t^{n+1}$,

1. Initialization:   $\left({\boldsymbol{u}}_{0},p_{0},\varphi_{0}\right)=\left({\boldsymbol{u}}^{n},p^{n},\varphi^{n}\right).$

2. Until convergence, knowing $\varphi_{k}$ and ${\boldsymbol{u}}_{k}$ compute:

[TABLE]

[TABLE]

[TABLE]

4.2 Analysis of the iterative scheme

The convergence of the fixed point (23)–(25) to a solution of the time discrete system (19), will be demonstrated in four steps: existence of the solution, regularity, convergence of the fixed point and finally identification of the limit as a solution of (19).

Although an existence result for (23)–(25) is our main goal in this first step, as $\varphi^{n}$ must be a phase field, the bound (21) must be established.

4.2.1 Well-posedness

Theorem 4.1

Assuming the hypothesis of Theorem 2.1 with ${\boldsymbol{u}}^{0}$ divergence free. For all step length $\Delta t>0$ and at each time step $n\geq 0$ and iteration $k\geq 0$ there exist a solution $\left({\boldsymbol{u}}_{k+1},p_{k+1},\varphi_{k+1}\right)$ of (23)–(25) in $\boldsymbol{H^{1}_{0}}(\Omega)\cap\boldsymbol{H^{2}}(\Omega)\times L^{2}_{0}(\Omega)\cap H^{1}(\Omega)\times H^{3}(\Omega)$. Moreover

[TABLE]

Proof 1

First we prove that the result holds for $n=0$ and $k=0$ which is the first step in the iterative process for $t^{1}=\Delta t$. Introducing

[TABLE]

we have to find $\varphi_{1}$ solution of

[TABLE]

where ${\boldsymbol{u}}_{0}={\boldsymbol{u}}^{0}$ and $\varphi_{0}=\varphi^{0}$. Existence and uniqueness of a solution in $H^{3}(\Omega)$ results from Lemma 2.2. However the coercivity needed induces a condition on $\Delta t$, namely

[TABLE]

From the definition of $\alpha$ we get

[TABLE]

If $\beta\geq 1$ there is no condition on $\Delta t$ on the other hand if $\Delta t<\eta^{2}/\gamma$ then $\beta$ can be an arbitrary positive constant. Therefore existence and uniqueness is obtained for all positive values of $\Delta t$.

Assume $\varphi_{1}$ achieves its maximum and minimum in $\Omega$ at ${\boldsymbol{x}}^{*},{\boldsymbol{x}}_{*}$ respectively

[TABLE]

From (27) at ${\boldsymbol{x}}_{*}$,${\boldsymbol{x}}^{*}$, using $|\varphi_{0}({\boldsymbol{x}})|\leq 1$ and $\beta\geq 0$, we have

[TABLE]

The role of the parameter $\beta$ is clear, it will induce a bound on $\varphi_{1}$ without the need for any condition on $\Delta t$. Denoting $\theta=\varphi_{0}({\boldsymbol{x}})=\varphi^{0}({\boldsymbol{x}})$, we want establish the values of $\beta$ giving simultaneously

[TABLE]

*Note that the coerciveness of $a(x)$ gives a strictly positive denominator. We introduce *

[TABLE]

[TABLE]

These polynomials have complex roots if $\beta\geq 9/8$. Under this condition there is only one real root for each of those polynomial and we get

[TABLE]

which gives us (28). For $\beta<9/8$ there is no general conclusion regarding these inequalities therefore $\beta\geq 9/8$ is only a sufficient condition for the uniform bound for $\varphi_{k+1}$.

Concerning the Navier-Stokes system (25), notice that at this stage of the loop $\varphi_{k+1}\in H^{3}(\Omega)$, $\rho_{k+1}$ and $\mu_{k+1}$ are known. Therefore the existence, uniqueness (based on the large viscosity assumption) and regularity of the solution $({\boldsymbol{u}}_{1},p_{1})$ comes form classical results (see [21, 22, 23] for instance). Proceeding by induction on $k$ for fixed $n=0$, we get the existence of $\left({\boldsymbol{u}}_{k+1},p_{k+1},\varphi_{k+1}\right)$ solution of (23)–(25) with $\varphi_{k+1}$ satisfying $\|\varphi_{k+1}\|_{L^{\infty}(\Omega)}\leq 1$.

Finally, having proved the assumptions for $n=0$, using again an induction argument over the time iteration $n$ with the same technique as presented here, we get the results for $\left({\boldsymbol{u}}_{k+1},p_{k+1},\varphi_{k+1}\right)$ solution of (23)–(25) at time $(n+1)\Delta t$.

Remark 4.1

As noted in the proof, the condition on $\beta$ is not a necessity, and we could have a uniform $L^{\infty}$ bound for smaller values of $\beta$. This kind of ”maximum principle” for $\varphi_{k+1}$ and $\varphi^{n}$ is fundamental as it defines a phase field on $\Omega$.

4.2.2 Regularity estimates

In what follows we rely on uniform bounds on $\varphi_{k}$ and ${\boldsymbol{u}}_{k}$ to get the convergence of the sequence $({\boldsymbol{u}}_{k},\varphi_{k})_{k\geq 0}$ in $\boldsymbol{H^{2}}(\Omega)\times H^{3}(\Omega)$. Applying Lemma 2.2 to problem (23) is not totally satisfactory since we get a $H^{2}$ bound on $\varphi_{k+1}$ depending (implicitly) on ${\boldsymbol{u}}_{k}$. However, with some additional work, this lemma can be used to get uniform bound on $\varphi_{k+1}$ and ${\boldsymbol{u}}_{k+1}$.

Lemma 4.1** (Bound on $(\varphi_{k})_{k\geq 0}$)**

At each time step $t^{n+1}=(n+1)\Delta t$, for all nonnegative integer $k$ and for all real number $\beta\geq 9/8$, the iterative solution $\varphi_{k+1}$ of (23) satisfies the following uniform bound

[TABLE]

where $\mathscr{C}^{\scriptscriptstyle{\nabla\varphi}}$ is a constant depending on $\beta,\eta,\gamma,|\Omega|$ and $\Delta t$.

Proof 2

The first bound is a direct consequence of $\|\varphi_{k+1}\|_{\infty}\leq 1$. The second one can be obtained by multiplying (23) by $\varphi_{k+1}$ and integrating over $\Omega$:

[TABLE]

Lemma 4.2** (Bound on $({\boldsymbol{u}}_{k})_{k\geq 0}$)**

With the hypothesis of Lemma 4.1 the iterative solution ${\boldsymbol{u}}_{k+1}$of (25) satisfies the following uniform bound

[TABLE]

Proof 3

We multiply (25) by ${\boldsymbol{u}}_{k+1}$ and we integrate over $\Omega$, we find

[TABLE]

By (9), the second line vanishes and by (10), the third one equal to

[TABLE]

Using (2) and Assumption 2.1, yields

[TABLE]

Based on Lemma 4.1 and Korn inequality, we obtain the desired estimates.

These last lemmas, although giving uniform bounds on $\varphi_{k}$, ${\boldsymbol{u}}_{k}$ and the tensor ${\boldsymbol{D}}({\boldsymbol{u}}_{k})$ do not provide sufficient regularity for the solution of (23)–(25) to get the convergence results sought for. Based on Lemmas 2.2, 4.1 and estimate (32) of Lemma 4.2, we can get uniform bounds in spaces more suited for the rest of the analysis.

Theorem 4.2

With the hypothesis of Lemma 4.1 the iterative solutions of schemes (23)–(25) satisfy the following regularity results

[TABLE]

*Where $\mathcal{K}^{\scriptscriptstyle{\varphi}}_{2},\mathcal{K}^{p},\mathcal{K}^{u},\mathscr{K}^{u}$ and $\mathscr{K}^{f}$ are non-negative constants following the convention in section 2.2. *

Proof 4

To bound $\varphi_{k+1}$ in $H^{2}(\Omega)$, we rewrite (23) as follow

[TABLE]

with

[TABLE]

Using Lemma 2.2 with $r=0$ and (29) we get

[TABLE]

In order to bound ${\boldsymbol{u}}_{k}\cdot\nabla\varphi_{k+1}$ in $L^{2}(\Omega)$, using a SPKCSHY–sequence

[TABLE]

Then (33) comes from (37) and (38). We are now in position to bound $p_{k}$. To do so, we first rewrite (25) as a Stokes system with variable viscosity and no-slip condition on $\partial\Omega$

[TABLE]

where

[TABLE]

Then, we have

[TABLE]

Since $p_{k+1}$ belongs to $L^{2}_{0}(\Omega)$, Poincaré inequality (see for instance [10] ) gives

[TABLE]

Based on Lemma 4.2, estimation (33) and Cauchy-Schwartz inequality, the first line in (40) can be bounded as follow

[TABLE]

As for the second and third line in (40), it can be bounded by

[TABLE]

Finally

[TABLE]

To bound $\|{\boldsymbol{u}}_{k}\|_{2}$, we first rearrange (39) using $\hat{p}_{k+1}=p_{k+1}/\mu_{k+1}$

[TABLE]

where ${\boldsymbol{f}}^{k,n}_{\boldsymbol{u}}$ is defined by (40). Next, following the regularity result for Stokes problem (see for instane [10]), we have

[TABLE]

The first two terms on the right hande side can be bounded (using a SPKCSHY–sequence) as follow

[TABLE]

Then

[TABLE]

And

[TABLE]

Then

[TABLE]

Combining (46), (45) with (44) there is two constants $\mathcal{K}^{u}$ and $\mathscr{K}^{u}$ such that,

[TABLE]

Going back to (40), each term in ${\boldsymbol{f}}^{k,n}_{\boldsymbol{u}}$ can be bounded using (31), (33) and applying a SPKCSHY–sequence. For $\alpha>0$

[TABLE]

For the second term (as ${\boldsymbol{u}}_{k}$ is divergence free)

[TABLE]

*For the third term *

[TABLE]

*Finally *

[TABLE]

Then, combining the inequalities of each terms, there is a $\alpha$ and a constant $\mathscr{K}^{f}$ such that

[TABLE]

and from (47)–(48),

[TABLE]

*which completes the proof. *

These uniform bounds suffice to conclude to the existence of a converging subsequence with very few constraints on the physical data. However, we are interested in establishing the convergence of the algorithm composed of the loop (23)–(25).

4.2.3 Convergence analysis and passage to the limit

For the next proofs, we introduce the following quantities

[TABLE]

It is clear that $1\leq\alpha_{0}$ for $\beta\geq 9/8$.

Theorem 4.3

Let $({\boldsymbol{u}}_{k},\varphi_{k})_{k\geq 0}$ be the sequence of solutions of (23)–(25) at time $t^{n+1}$. Assuming $\beta\geq 9/8$ then for all nonnegative integers $m,k$:

[TABLE]

*Where $\mathscr{C}_{\varphi 0}$ depends only on $\gamma,\eta,\beta$ and $\Delta t$. *

Proof 5

To simplify the presentation we introduce the notation

[TABLE]

*The proof will be split in two parts, beginning by the analysis of the sequence $\varphi_{k}$.

Analysis of the sequence $(\varphi_{k})_{k\geq 0}$.* Taking (23) at two different iterations, $k$ and $m$, multiplying each equalities by $\delta\varphi$ and subtracting them, we get*

[TABLE]

From the bound on $(\varphi_{k})_{k}$, we have

[TABLE]

and using (11) for the term in ${\boldsymbol{u}}$ we have

[TABLE]

From which we get

[TABLE]

and (50) follows by choosing $\mathscr{C}_{\varphi 0}=\displaystyle\max\left(\frac{\alpha_{2}}{\alpha_{0}},\sqrt{\frac{\Delta t}{\alpha_{0}\gamma}}\right).$

In same manner we get a bound on the gradient of the difference and a $H^{1}$ bound

[TABLE]

with

[TABLE]

Analysis of the sequence $({\boldsymbol{u}}_{k})_{k\geq 0}$.* We multiply (25) at $(k+1)$ and $(m+1)$ steps by $\delta{\boldsymbol{u}}$ and we integrate over $\Omega$, yielding a sum of seven terms*

[TABLE]

To establish (51) we will take advantage of the fact that the first three terms of this sum contains squared norms. Labeling the terms $A_{0}$ through $A_{6}$ we have

[TABLE]

Adding and subtracting $\rho_{m+1}{\boldsymbol{u}}_{k+1}$ in $A_{0}$ and using (24),

[TABLE]

Using (8) we have

[TABLE]

We can rewrite $A_{2}$ as follow:

[TABLE]

From these equality we get

[TABLE]

We will now establish an appropriate bound for the right hand side of this last inequality. Using SPKCSHY–sequences,

[TABLE]

For $A_{3}$, using (10)

[TABLE]

For $A_{4}+A_{5}$, we use (12), moreover since all velocities are divergence free we get parts of both terms canceling each others, leaving

[TABLE]

For those four terms we get (through SPKCSHY–sequences)

[TABLE]

[TABLE]

[TABLE]

[TABLE]

From (59)–(62), we have

[TABLE]

From Assumption 2.1 and $\displaystyle\left|\sqrt{\rho_{m+1}}-\sqrt{\rho_{k+1}}\right|\leq\frac{|\delta\rho|}{4\sqrt{\rho_{I}}}\left|\delta\varphi\right|,\,\forall(k,m)\in\mathbb{N}^{2}$,

[TABLE]

Combining (57)–(58),(63)–(64) in (56) and introducing the notation

[TABLE]

we get

[TABLE]

Combining (50), (66), and (54), we get (51) with

[TABLE]

The rest of this analysis is based on the fact that, using (50) and (51), we can clearly provide conditions on $\mathscr{C}_{\varphi 0}$ and $K_{u0}$ such that the sequence $\left({\boldsymbol{u_{k}}},\varphi_{k}\right)$ converges strongly in $\boldsymbol{H^{1}}(\Omega)\times L^{2}(\Omega)$. This relies on our ability to control those constants.

For $\mathscr{C}_{\varphi 0}$ as it depends on two arbitrary parameters, $\beta$ and $\Delta t$ its control is relatively simple. As for $K_{u0}$, using $\mu_{I}$ seems to be the simplest route. However this means to impose a lower bound on the minimal value of the viscosity. Let us underline that adding such condition is not specific to the approach proposed here and can be found in other papers such as [24, 25, 26]. Furthermore, imposing a lower bound on the viscosity, as a condition for existence and uniqueness of solution, is relatively frequent even in the stationnary case [10, Chap. 3].

This leads us to study the behaviour of $K_{u0}$ with respect to $\mu_{I}$. From the convention defining the constants (section 2.2), $\mathscr{C}_{\varphi 1}$ do not depend on $\mu_{I}$. As for $\mathscr{K}_{uu}$ it is a linear expression with respect to $1/\sqrt{\mu_{I}}$. Which leaves $K_{u\varphi}$ to examine. From (42), (49), and (65), we get, again based on the convention for the constants,

[TABLE]

From this we can conclude that $K_{u0}$ goes to zero when $\mu_{I}$ goes to infinity if

[TABLE]

Remark 4.2

The condition (67) is quite general, a simpler, less generic, condition can be considered. To have

[TABLE]

we need to assume that $\mu_{I}$ is sufficiently large but also that $\mu_{S}/\mu_{I}$ is small enough. Therefore it is sufficient to assume that there is a $C_{\mu}>0$ such that

[TABLE]

to make sure that we can control $K_{u0}$ through $\mu_{I}$.

With (67) or (68), we get, for $\mu_{I}$ sufficiently large, the convergence of the whole sequence $\left({\boldsymbol{u_{k}}},\varphi_{k}\right)$ in $\boldsymbol{H^{1}}(\Omega)\times L^{2}(\Omega)$ (and in $\boldsymbol{H^{1}}(\Omega)\times H^{1}(\Omega)$). This will be the subject of Corollary 4.1.

Corollary 4.1

*Under the assumptions of Theorem 2.1, with ${\boldsymbol{u}}^{0}$ divergence free. Assuming $\mu_{I}$ is sufficiently large and the following conditions are satisfied *

[TABLE]

then, at each time step, the sequence $({\boldsymbol{u}}_{k},\varphi_{k})_{k\geq 0}$ of solutions of (23)–(25) converge strongly in $\boldsymbol{H^{1}}(\Omega)\times H^{1}(\Omega)$ and for all integers $k$,$m$

[TABLE]

where

[TABLE]

Proof 6

*Since $\beta\geq 9/8$ Theorem 4.3 is valid, adding (50), (51), we have (70). It is relatively easy to show that $K<1$ and $({\boldsymbol{u}}_{k},\varphi_{k})$ converges if the conditions on $\beta$ and $\Delta t$ are satisfied. *

Remark 4.3

*Since we assume that $\eta^{2}/\gamma\ll 1$ (i.e. $\eta$ is nearly (or goes) to zero) in most case the condition $\beta\geq 9/8$ is sufficient. For a fixed time step, as $\mathscr{C}_{\varphi 0}>0$, $\beta$ can be seen as an ”acceleration coefficient”, as it permits to lower the convergence rate of the sequence. Alternatively, if $K_{u0}$ approches $1$ then $\beta$ should be large to enforce the bound on the phase field and the convergence of the sequence. *

Theorem 4.4

Under the hypothesis of Corollary 4.1, at each time step $t^{n}=n\Delta t$, the sequence $({\boldsymbol{u}}_{k},p_{k},\varphi_{k})_{k\geq 0}$ of solutions of (23)–(25) converge strongly in $\boldsymbol{H^{1}_{0}}(\Omega)\cap\boldsymbol{H^{2}}(\Omega)\times L^{2}_{0}(\Omega)\cap H^{1}(\Omega)\times H^{3}(\Omega)$ to a solution $({\boldsymbol{u}}^{n},p^{n},\varphi^{n})$ of (19), (16)–(17). Moreover

[TABLE]

Proof 7

Assuming the sequence of the solutions of (23)–(25) converge strongly in $\boldsymbol{H^{1}_{0}}(\Omega)\cap\boldsymbol{H^{2}}(\Omega)\times L^{2}_{0}(\Omega)\cap H^{1}(\Omega)\times H^{3}(\Omega)$, it is obvious that the limit of this sequence $({\boldsymbol{u}}^{n},p^{n},\varphi^{n})$ is a solution of (19) with (16)–(17). Moreover, under this assumption, using Theorem 4.1 we have $\|\varphi^{n}\|_{{\boldsymbol{L^{\infty}}}(\Omega)}\leq 1.$

To complete the proof all we need is to demonstrate the strong convergence of the solutions of (23)–(25) in $\boldsymbol{H^{1}_{0}}(\Omega)\cap\boldsymbol{H^{2}}(\Omega)\times L^{2}_{0}(\Omega)\cap H^{1}(\Omega)\times H^{3}(\Omega)$. Subtracting (23) for two integers $k\neq m$ then using the test function

[TABLE]

(from Theorem 4.1 we have $\varphi_{k}\in H^{3}(\Omega)$) and using (12) once again

[TABLE]

which gives

[TABLE]

From Corollary 4.1 we get the convergence of $\Delta\varphi_{k}$ in $L^{2}(\Omega)$ hence, recalling Theorem 4.2, the convergence of $\varphi_{k}$ in $H^{2}(\Omega)$.

Concerning the convergence of $\varphi_{k}$ in $H^{3}(\Omega)$, we use the same technique. From Theorem 4.1, we can apply the $\nabla$ operator on (23) at $k$ and $m$. Subtracting both equations and using the test function

[TABLE]

we get the convergence of $\Delta\nabla\varphi_{k}$ in $L^{2}(\Omega)$ from which we conclude the convergence of $\varphi_{k}$ in $H^{3}$. For the convergence of the sequence ${\boldsymbol{u}}_{k}$ in $\boldsymbol{H^{2}}(\Omega)$, we introduce

[TABLE]

*Based on (25) (Corollary 4.1 gives the convergence in $\boldsymbol{H^{1}}(\Omega)$ for $({\boldsymbol{u}}_{k})_{k\geq 0}$ and we prooved the convergence in $H^{3}(\Omega)$ for $(\varphi_{k})_{k\geq 0}$), subtracting (25) at $m\neq k$, we conclude that ${\boldsymbol{v}}_{k,m}$ and ${\boldsymbol{v}}_{k,m}^{R}$ converge to zero. Then, from *

[TABLE]

*we get $({\boldsymbol{u}}_{k})_{k\geq 0}$ converges in $\boldsymbol{H^{2}}(\Omega)$. Finally, for the convergence of the pressure sequence, the usual argument based on the inf–sup condition is used (see [10] for example). *

This also completes the demonstration of Theorem 3.1 which is merely a corollary of Theorem 4.4.

5 Numerical results

The goal of this section is to illustrate the effectiveness of a totally implicit scheme for the Navier-Sokes/Allen-Cahn model. With that in mind, a single numerical test is sufficient: an analytical test based on a two dimensional manufactured solution. The fixed point strategy constructed to establish Theorem 3.1, offers a first numerical recipe to approximate the solution of (6).

5.1 Methodology

Obviously various strategies, efficient and well adapted to the physical context, could be used to discretise and solve (23) and (25). Here, the finite element method has been retained for the spatial discretisation.

We used proper degrees of interpolation for each variables. For the fluid a Taylor-Hood interpolation was retained, for the velocity a quadratic interpolation (P2) and for the pressure a linear (P1) interpolation, insuring the respect of the inf–sup condition (see [27, 13]). As for the phase field we chose a P2 interpolation. Since efficiency of the fully discrete algorithm is not the purpose of this work, we chose to leave the system as is: using a $LU$ solver on a pressure penalised algebraic system for $\eqref{AlgNSDivu}$, with a value of $10^{-8}$ for the pressure penalisation.

For these tests, a uniform triangular mesh of sufficiently small size, $100\times 100$, is used, this insure a negligible spatial error for the chosen range of time steps. The algorithm is implemented using FreeFem++ a finite element software, see [28]. As for the fixed point (Step 5 in Algorithm 1), a tolerance of $10^{-9}$ on the sum of the norm of the variations of $\varphi$ and ${\boldsymbol{u}}$ is used as a stopping criteria.

5.2 Algorithms

Applying the spatial discretization to (23)–(25), a basic strongly coupled implicit algorithm, called fully implicit Newton method (FIN) is obtained

Since we are using a fixed point approach, we also consider the following linearization

[TABLE]

in (20). We call this modification a Picard fixed point approach. In that case $\beta\geq~{}2$ suffice to insure uniform boundedness of the phase field and Theorem 4.1 is valid with the modified bound on $\beta$. This simple linearisation seems a natural choice, we propose to use it as a comparative. Replacing (23) by

[TABLE]

we still get a strongly coupled fully implicit semi-discretization. Using the finite element method on (72), (24)–(25) we get a new algorithm, called the fully implicit Picard method (FIP).

As a second comparative approach, we built an explicit formulation by putting $\beta=0$ and expliciting the non linear term. Replacing (23) by

[TABLE]

we get a strongly coupled explicit semi-discretization approach. Using the finite element method on (73), (24)–(25) we get a new algorithm, called the strongly coupled explicit method (SCE).

5.3 Numerical test

For the numerical experiments a problem having an analytical solution is used. It is inspired by the the finite element tests proposed in [7, Accuracy test 4.2]

[TABLE]

Note that $\varphi(t,x,y)\in[-1,1]$ is a polynomial of degree two in space and linear in time, therefore the approximation, using degree two in space and BDF1 in time, should be exact. Suitable boundary conditions on ${\boldsymbol{u}}$ and $\varphi$ and forcing terms depending on the analytical solution are added.

Table 1 gives the $L^{2}$ error on ${\boldsymbol{u}}$ and $\varphi$ at time $10/1\,300$. As a measure of precision we define

[TABLE]

The computational cost is quantified using the total number of times the algebraic system corresponding to the Navier-Stokes equations (25) is solved. We limit ourselves to a few time step, as it suffice to convincingly demonstrate our point. For the fully implicit Newton method, the time step $\Delta t=\eta^{2}/(13\gamma)$ and parameter $\beta=9/8$ are fixed following Corollary 4.1. For the FIP method, we use the same time step but with $\beta=2$, a value insuring uniform boundedness of $\varphi$ for this method. As the lower bound on $\beta$ is only a sufficient condition, two tests using $\beta=0$ for the FIN and FIP methods are also presented.

As usual in numerical methods the context of use is important. This example does not exclude the relevance of simpler approaches such as explicit or weakly coupled methods. It simply underline the effectiveness of implicit and strongly coupled methods. Concerning the use of a weakly coupled scheme, this relatively simple example does not seem to land itself to such methods. Weakly coupled strategies were excluded from this analysis as in our very few experiment they all exhibited poor precision even when compared to the SCE method. Let summarize some observations resulting from Table 1 and Figure 1

For the FIN and FIP method the level of precision for $\varphi$ is such that we can consider the error on the velocity ${\boldsymbol{u}}$ as produced by the finite element approximation of the Navier-Stokes equation.

• 2.

There is a clear gain in precision for the implicit methods, FIN and FIP, when compared to the explicit approach. At equal precision (or at equal time step), the choice of an implicit scheme is obvious. Based on the two experiments with the SCE method and the gap in precision, argument in favor of the explicit method, based on efficiency, cannot be made.

• 3.

Generally speaking, considering the loss of information in the Picard method, a Newton-type approach seems preferable. However, in this specific case, from Figure 1, these two methods are relatively similar (in efficiency and precision) when using $\beta\in[0.75,4]$.

• 4.

The parameter $\beta$ is basically a ”theoretical trick” producing Theorem 3.1. In the FIN method, $\beta>0$ can be interpreted as perturbing the Newton method. For the FIN method, from a numerical perspective, the use of $\beta=0$, is optimal in term of efficiency (Tableau 1 and Figure 1).

6 Conclusion

This work propose a proof of the well posed character of a fully implicit strongly coupled semi-discrete NS-AC system (Theorem 3.1). We consider this theoretical result as the basis needed for the analysis of various numerical scheme (not restricted to finite element approaches)) aiming at the approximation of the solution of the NS-AC system. The demonstration offered a simple and relatively efficient time scheme which allowed us to produce a finite element approximation. The numerical tests have established clearly the relevance of implicit and strongly coupled schemes for the NS-AS system, therefore reinforcing the significance of Theorem 3.1.