A Hybrid High-Order method for the convective Cahn-Hilliard problem in mixed form

TL;DR

This paper introduces a new Hybrid High-Order numerical method for solving the convective Cahn-Hilliard problem, capable of handling complex meshes and demonstrating robustness across different flow regimes.

Contribution

It develops a flexible, high-order hybrid method applicable in 2D and 3D with general meshes, advancing numerical solutions for convection-influenced phase separation.

Findings

Method is valid in 2D and 3D

Supports arbitrary approximation orders

Shows robustness with respect to Péclet number

Abstract

We propose a novel Hybrid High-Order method for the Cahn-Hilliard problem with convection. The proposed method is valid in two and three space dimensions, and it supports arbitrary approximation orders on general meshes containing polyhedral elements and nonmatching interfaces. An extensive numerical validation is presented, which shows robustness with respect to the P\'eclet number.

Equations

$${\mathrm{d}}_{t}c-\frac{1}{\mathrm{Pe}}\Delta w+\nabla \cdot (\vec{u}c)$$

$$w$$

$$c(0)$$

$${\partial }_{\vec{n}}c={\partial }_{\vec{n}}w$$

$${\delta }_t{\phi }^n:=\frac{{\phi }^n-{\phi }^{n-1}}{\tau }\in V.$$

$$\underline{U}_{h}^{k}:=\Bigg{(}\times_{T\in\mathcal{T}_{h}}\mathbb{P}^{k+1}(T)\Bigg{)}\times\Bigg{(}\times_{F\in\mathcal{F}_{h}}\mathbb{P}^{k}(F)\Bigg{)}$$

$$(\nabla {\mathrm{p}}_T^{k+1}{\underline{v}}_T,\nabla z{)}_T$$

$$a_T({\underline{u}}_T,{\underline{v}}_T):=(\nabla {\mathrm{p}}_T^{k+1}{\underline{u}}_T,\nabla {\mathrm{p}}_T^{k+1}{\underline{v}}_T{)}_T+s_T({\underline{u}}_T,{\underline{v}}_T),$$

$$s_T({\underline{u}}_T,{\underline{v}}_T):=\sum _{F\in {\mathcal{F}}_T}{h}_F^{-1}({\pi }_F^k(u_F-u_T),{\pi }_F^k(v_F-v_T){)}_F,$$

$$({\mathrm{G}}_{\vec{u},T}^{k+1}{\underline{v}}_T,w{)}_T$$

$$b_{\vec{u},T}({\underline{u}}_T,{\underline{v}}_T):=({\mathrm{G}}_{\vec{u},T}^{k+1}{\underline{u}}_T,v_T{)}_T+s_{\vec{u},T}({\underline{u}}_T,{\underline{v}}_T).$$

$$s_{\vec{u},T}({\underline{u}}_T,{\underline{v}}_T):=\sum _{F\in {\mathcal{F}}_T}(\frac{|\vec{u}\cdot {\vec{n}}_{TF}|-\vec{u}\cdot {\vec{n}}_{TF}}{2}(u_F-u_T),v_F-v_T{)}_F.$$

$$a_h({\underline{u}}_h,{\underline{v}}_h):=\sum _{T\in {\mathcal{T}}_h}{a}_T({\underline{u}}_T,{\underline{v}}_T),b_{\vec{u},h}({\underline{u}}_h,{\underline{v}}_h):=\sum _{T\in {\mathcal{T}}_h}{b}_{\vec{u},T}({\underline{u}}_T,{\underline{v}}_T).$$

$$({\delta }_t{c}_h^n,{\phi }_h)+\frac{1}{\mathrm{Pe}}{a}_h({\underline{w}}_h^n,{\underline{\phi }}_h)+b_{\vec{u},h}({\underline{c}}_h^n,{\underline{\phi }}_h)=0$$

$$(w_h^n,{\psi }_h)=({\Phi }^{\prime }(c_h^n),{\psi }_h)+{\gamma }^2{a}_h({\underline{c}}_h^n,{\underline{\psi }}_h)$$

$$c_0(\vec{x}):=\tanh (\frac{2x_1-1}{2\sqrt{2}{\gamma }^2}),\vec{u}(\vec{x}):=20\cdot \left(\begin{array}{c}{x}_1(x_1-1)(2x_2-1)\\ -x_2(x_2-1)(2x_1-1)\end{array}\right),\forall \vec{x}\in \Omega .$$

$$\vec{u}(\vec{x}):=\frac{1}{2}\left(1+\tanh (80-200\Vert (x_1-0.5,x_2-0.5){\Vert }_2)\right)\cdot \left(\begin{array}{c}2x_2-1\\ 1-2x_1\end{array}\right),\forall \vec{x}\in \Omega .$$

$$\vec{u}(\vec{x}):=\left(\begin{array}{c}\sin (\pi x_1)\cos (\pi x_2)\\ -\cos (\pi x_1)\sin (\pi x_2)\end{array}\right),\forall \vec{x}\in \Omega .$$

Full text

MSC (2010): 65N08, 65N30 65N12

1 Cahn–Hilliard equation

Let $\Omega\subset\mathbb{R}^{d}$, $d\in\{2,3\}$, denote a bounded connected convex polyhedral domain with Lipschitz boundary $\partial\Omega$ and outward normal $\vec{n}$, and let $t_{F}>0$. The convective Cahn–Hilliard problem consists in finding the order-parameter $c:\Omega\times\left(0,t_{F}\right]\to\mathbb{R}$ and the chemical potential $w:\Omega\times\left(0,t_{F}\right]\to\mathbb{R}$ such that

[TABLE]

where $\gamma>0$ is the interface parameter (usually taking small values), Pe$>0$ is the Péclet number, $\vec{u}$ the velocity field such that $\nabla\cdot\vec{u}=0$ in $\Omega$ and $\Phi$ the free-energy such that $\Phi(c):=\frac{1}{4}(1-c^{2})^{2}$. This formulation is an extension of the Cahn–Hilliard model originally introduced in Cahn.Hilliard:58 ; Cahn:61 and a first step towards coupling with the Navier–Stokes equations.

In this work we extend the HHO method of Chave.Di-Pietro.ea:16 to incorporate the convective term in (1a). Therein, a full stability and convergence analysis was carried out for the non-convective case, leading to optimal estimates in $(h^{k+1}+\tau)$ (with $h$ denoting the meshsize and $\tau$ the time step) for the the $C^{0}(H^{1})$-error on the order-parameter and $L^{2}(H^{1})$-error on the chemical potential. The convective term is treated in the spirit of Di-Pietro.Droniou.ea:14 , where a HHO method fully robust with respect to the Péclet number was presented for a locally degenerate diffusion-advection-reaction problem.

The proposed method offers various assets: (i) fairly general meshes are supported including polyhedral elements and nonmatching interfaces; (ii) arbitrary polynomial orders, including the case $k=0$, can be considered; (iii) when using a first-order (Newton-like) algorithm to solve the resulting system of nonlinear algebraic equations, element-based unknowns can be statically condensed at each iteration.

The rest of this paper is organized as follows: in Section 2, we recall discrete setting including notations and assumptions on meshes, define localy discrete operators and state the discrete formulation of (1). In Section 3, we provide an extensive numerical validation.

2 The Hybrid High-Order method

In this section we recall some assumptions on the mesh, introduce the notation, and state the HHO discretization.

2.1 Discrete setting

We consider sequences of refined meshes that are regular in the sense of (Di-Pietro.Ern:12, , Chapter 1). Each mesh $\mathcal{T}_{h}$ in the sequence is a finite collection $\{T\}$ of nonempty, disjoint, polyhedral elements such that $\overline{\Omega}=\bigcup_{T\in\mathcal{T}_{h}}\overline{T}$ and $h=\max_{T\in\mathcal{T}_{h}}h_{T}$ (with $h_{T}$ the diameter of $T$). For all $T\in\mathcal{T}_{h}$, the boundary of $T$ is decomposed into planar faces collected in the set $\mathcal{F}_{T}$. For admissible mesh sequences, card($\mathcal{F}_{T}$) is bounded uniformly in $h$. Interfaces are collected in the set $\mathcal{F}_{h}^{\text{i}}$, boundary faces in $\mathcal{F}_{h}^{\text{b}}$ and we define $\mathcal{F}_{h}:=\mathcal{F}_{h}^{\text{i}}\cup\mathcal{F}_{h}^{\text{b}}$. For all $T\in\mathcal{T}_{h}$ and all $F\in\mathcal{F}_{T}$, the diameter of $F$ is denoted by $h_{F}$ and the unit normal to $F$ pointing out of $T$ is denoted by $\vec{n}_{TF}$.

To discretize in time, we consider for sake of simplicity a uniform partition $(t^{n})_{0\leq n\leq N}$ of the time interval $[0,t_{F}]$ with $t^{0}=0$, $t^{N}=t_{F}$ and $t^{n}-t^{n-1}=\tau$ for all $1\leq n\leq N$. For any sufficiently regular function of time $\varphi$ taking values in a vector space $V$, we denote by $\varphi^{n}\in V$ its value at discrete time $t^{n}$, and we introduce the backward differencing operator $\delta_{t}$ such that, for all $1\leq n\leq N$,

[TABLE]

2.2 Local space of degrees of freedom

For any integer $l\geq 0$ and $X$ a mesh element or face, we denote by $\mathbb{P}^{l}(X)$ the space spanned by the restrictions to $X$ of $d$-variate polynomials of order $l$. Let

[TABLE]

be the global degrees of freedoms (DOFs) space with single-valued interface unknowns. We denote by $\underline{v}_{h}=((v_{T})_{T\in\mathcal{T}_{h}},(v_{F})_{F\in\mathcal{F}_{h}})$ a generic element of $\underline{U}_{h}^{k}$ and by $v_{h}$ the piecewise polynomial function such that $v_{h}|_{T}=v_{T}$ for all $T\in\mathcal{T}_{h}$. For any $T\in\mathcal{T}_{h}$, we denote by $\underline{U}_{T}^{k}$ and $\underline{v}_{T}=(v_{T},(v_{F})_{F\in\mathcal{F}_{T}})$ the restrictions to $T$ of $\underline{U}_{h}^{k}$ and $\underline{v}_{h}$, respectively.

2.3 Local diffusive contribution

Consider a mesh element $T\in\mathcal{T}_{h}$. We define the local potential reconstruction $\mathrm{p}_{T}^{k+1}:\underline{U}_{T}^{k}\rightarrow\mathbb{P}^{k+1}(T)$ such that, for all $\underline{v}_{T}:=(v_{T},(v_{F})_{F\in\mathcal{F}_{T}})\in\underline{U}_{T}^{k}$ and all $z\in\mathbb{P}_{T}^{k+1}$,

[TABLE]

with closure condition $\int_{T}(\mathrm{p}_{T}^{k+1}\underline{v}_{T}-v_{T})=0$. We introduce the local diffusive bilinear form $a_{T}$ on $\underline{U}_{T}^{k}\times\underline{U}_{T}^{k}$ such that, for all $(\underline{u}_{T},\underline{v}_{T})\in\underline{U}_{T}^{k}\times\underline{U}_{T}^{k}$

[TABLE]

with stabilization bilinear form $s_{T}:\underline{U}_{T}^{k}\times\underline{U}_{T}^{k}\rightarrow\mathbb{R}$ such that

[TABLE]

where, for all $F\in\mathcal{F}_{h}$, $\pi_{F}^{k}:L^{1}(F)\rightarrow\mathbb{P}^{k}(F)$ denotes the $L^{2}$-orthogonal projector onto $\mathbb{P}^{k}(F)$.

2.4 Local convective contribution

For any mesh element $T\in\mathcal{T}_{h}$, we define the local convective derivative reconstruction $\mathrm{G}_{\vec{u},T}^{k+1}:\underline{U}_{T}^{k}\rightarrow\mathbb{P}^{k+1}(T)$ such that, for all $\underline{v}_{T}:=(v_{T},(v_{F})_{F\in\mathcal{F}_{T}})\in\underline{U}_{T}^{k}$ and all $w\in\mathbb{P}^{k+1}(T)$,

[TABLE]

The local convective contribution $b_{\vec{u},T}$ on $\underline{U}_{T}^{k}\times\underline{U}_{T}^{k}$ is such that, for all $(\underline{u}_{T},\underline{v}_{T})\in\underline{U}_{T}^{k}\times\underline{U}_{T}^{k}$

[TABLE]

with local upwind stabilization bilinear form $s_{\vec{u},T}:\underline{U}_{T}^{k}\times\underline{U}_{T}^{k}\rightarrow\mathbb{R}$ such that

[TABLE]

Notice that the actual computation of $\mathrm{G}_{\vec{u},T}^{k+1}$ is not required, as one can simply use its definition to expand the cell-based term in the bilinear form $b_{\vec{u},T}$.

2.5 Discrete problem

Denote by $\underline{U}_{h,0}^{k}:=\{\underline{v}_{h}=((v_{T})_{T\in\mathcal{T}_{h}},(v_{F})_{F\in\mathcal{F}_{h}})\in\underline{U}_{h}^{k}|\int_{\Omega}v_{h}=0\}$ the zero-average DOFs subspace of $\underline{U}_{h}^{k}$. We define the global bilinear forms $a_{h}$ and $b_{\vec{u},h}$ on $\underline{U}_{h}^{k}\times\underline{U}_{h}^{k}$ such that, for all $(\underline{u}_{h},\underline{v}_{h})\in\underline{U}_{h}^{k}\times\underline{U}_{h}^{k}$

[TABLE]

The discrete problem reads: For all $1\leq n\leq N$, find $(\underline{c}_{h}^{n},\underline{w}_{h}^{n})\in\underline{U}_{h,0}^{k}\times\underline{U}_{h}^{k}$ such that

[TABLE]

where $\underline{c}_{h}^{0}\in\underline{U}_{h,0}^{k}$ solves $a_{h}(\underline{c}_{h}^{0},\underline{\varphi}_{h})=-(\Delta c_{0},\varphi_{h})$ for all $\underline{\varphi}_{h}\in\underline{U}_{h}^{k}$.

3 Numerical test cases

In this section, we numerically validate the HHO method.

3.1 Disturbance of the steady solution

For the first test case, we use a piecewise constant approximation ($k=0$), discretize the domain $\Omega=(0,1)^{2}$ by a triangular mesh ($h=1.92\cdot 10^{-3}$) with $\gamma=5\cdot 10^{-2}$, $\tau=\gamma^{2}$ and $\mathrm{Pe}=1$. The initial condition for the order-parameter and the velocity field are given by

[TABLE]

The result is depicted in Figure 2 and shows that the method is well-suited to capture the interface dynamics subject to a strong velocity fields.

3.2 Thin interface between phases

For the second example, we also use a piecewise constant approximation ($k=0$) with a Cartesian discretization of the domain $\Omega=(0,1)^{2}$, where $h=1.95\cdot 10^{-3}$. The interface parameter is taken to be very small $\gamma=5\cdot 10^{-3}$, the time step is $\tau=1\cdot 10^{-5}$ and $\mathrm{Pe}=50$. The initial condition for the order-parameter is taken to be a random value between $-1$ and $1$ inside a circular partition of the Cartesian mesh and $-1$ outside. The velocity field is given by

[TABLE]

See Figure 3 for the numerical result. The method is robust with respect to $\gamma$ and is also well-suited to approach the thin high-gradient area of the order-parameter.

3.3 Effect of the Péclet number

The Péclet number is the ratio of the contributions to mass transport by convection to those by diffusion: when ${\rm Pe}$ is greater than one, the effects of convection exceed those of diffusion in determining the overall mass flux. In the last test case, we compare several time evolutions obtained with different values of the Péclet number (${\rm Pe}\in\{1,50,200\}$), starting from the same initial condition. We use a Voronoi discretization of the domain $\Omega=(0,1)^{2}$, where $h=9.09\cdot 10^{-3}$, and use piecewise linear approximation ($k=1$). We choose $\gamma=1\cdot 10^{-2}$, $\tau=1\cdot 10^{-4}$ and $t_{F}=1$. The initial condition is given by a random value between $-1$ and $1$ inside a circular domain of the Voronoi mesh and $-1$ outside. The convective term is given by

[TABLE]

Snapshots of the order parameter at several times are shown on Figure 4 for each value of the Péclet number. For each case, the method takes into account the value of $\mathrm{Pe}$ and appropriately models the evolution of the order parameter by prevailing advection to diffusion when ${\rm Pe}\gg 1$.