Numerical Methods for a System of Coupled Cahn-Hilliard Equations

TL;DR

This paper develops and tests numerical methods for a coupled system of Cahn-Hilliard equations modeling phase separation in polymer blends, combining existing schemes and analyzing their behavior through simulations.

Contribution

It introduces a novel combination of numerical schemes tailored for coupled Cahn-Hilliard equations and demonstrates their effectiveness through detailed simulations.

Findings

Numerical methods accurately capture phase separation dynamics.

Parameter variations significantly affect the numerical solutions.

The approach is verified through comprehensive tests.

Abstract

In this work, we study a system of coupled Cahn-Hilliard equations describing the phase separation of a copolymer and a homopolymer blend. The numerical methods we propose are based on suitable combinations of existing schemes for the single Cahn-Hilliard equation. As a verification for our approach, we present some tests and a detailed description of the numerical solutions' behaviour obtained by varying the values of the parameters.

Tables

Table 1. Table 1. Comparison of the evolution of u 𝑢 u and v 𝑣 v .

Method	$t=0.2$	$t=2$	$t=10$
OD2
LS
EY

Table 2. Table 2. Evolution of u 𝑢 u (left) and v 𝑣 v (right) with parameters ϵ u = ϵ v = 0.05 subscript italic-ϵ 𝑢 subscript italic-ϵ 𝑣 0.05 \epsilon_{u}=\epsilon_{v}=0.05 , τ u = 1 subscript 𝜏 𝑢 1 \tau_{u}=1 , σ = 100 𝜎 100 \sigma=100 , α = 0.04 𝛼 0.04 \alpha=0.04 and β = − 0.9 𝛽 0.9 \beta=-0.9 .

${\tau }_v$	$t=3$	$t=6$	$t=15$
$10$
$250$

Table 3. Table 3. Evolution of u 𝑢 u (left) and v 𝑣 v (right) with parameters ϵ u = 0.05 subscript italic-ϵ 𝑢 0.05 \epsilon_{u}=0.05 , τ u = 1 subscript 𝜏 𝑢 1 \tau_{u}=1 , τ v = 100 subscript 𝜏 𝑣 100 \tau_{v}=100 , σ = 60 𝜎 60 \sigma=60 , α = 0.02 𝛼 0.02 \alpha=0.02 and β = − 0.9 𝛽 0.9 \beta=-0.9 .

${ϵ}_v$	$t=0.9$	$t=6$	$t=15$
$0.01$
$0.03$
$0.05$

Table 4. Table 4. Evolution of u 𝑢 u (left) and v 𝑣 v (right) with parameters ϵ u = ϵ v = 0.05 subscript italic-ϵ 𝑢 subscript italic-ϵ 𝑣 0.05 \epsilon_{u}=\epsilon_{v}=0.05 , τ u = 1 subscript 𝜏 𝑢 1 \tau_{u}=1 , τ v = 100 subscript 𝜏 𝑣 100 \tau_{v}=100 , α = 0.4 𝛼 0.4 \alpha=0.4 and β = − 0.9 𝛽 0.9 \beta=-0.9 .

$\sigma$	$t=3$	$t=6$	$t=15$
$0$
$50$
$150$

Table 5. Table 5. Evolution of u 𝑢 u (left) and v 𝑣 v (right) with parameters ϵ u = ϵ v = 0.05 subscript italic-ϵ 𝑢 subscript italic-ϵ 𝑣 0.05 \epsilon_{u}=\epsilon_{v}=0.05 , τ u = 1 subscript 𝜏 𝑢 1 \tau_{u}=1 , τ v = 100 subscript 𝜏 𝑣 100 \tau_{v}=100 , σ = 100 𝜎 100 \sigma=100 , and α = 0.01 𝛼 0.01 \alpha=0.01 .

$\beta$	$t=3$	$t=6$	$t=15$
$-0.3$
$-0.5$
$-0.9$

Table 6. Table 6. Evolution of u 𝑢 u (left) and v 𝑣 v (right) with parameters ϵ u = ϵ v = 0.05 subscript italic-ϵ 𝑢 subscript italic-ϵ 𝑣 0.05 \epsilon_{u}=\epsilon_{v}=0.05 , τ u = 1 subscript 𝜏 𝑢 1 \tau_{u}=1 , τ v = 100 subscript 𝜏 𝑣 100 \tau_{v}=100 , σ = 100 𝜎 100 \sigma=100 , and β = − 0.3 𝛽 0.3 \beta=-0.3 .

$\alpha$	$t=3$	$t=6$	$t=15$
$-0.4$
$-0.08$
$0.01$
$0.1$
$0.4$

Table 7. Table 7. Evolution of v 𝑣 v with parameters ϵ u = ϵ v = 0.05 subscript italic-ϵ 𝑢 subscript italic-ϵ 𝑣 0.05 \epsilon_{u}=\epsilon_{v}=0.05 , τ u = 1 subscript 𝜏 𝑢 1 \tau_{u}=1 , τ v = 100 subscript 𝜏 𝑣 100 \tau_{v}=100 , σ = 100 𝜎 100 \sigma=100 , and α = 0.01 𝛼 0.01 \alpha=0.01 , β = − 0.9 𝛽 0.9 \beta=-0.9 . Notice that ∫ Ω v ​ ( t ) = 0.01146 subscript Ω 𝑣 𝑡 0.01146 \int_{\Omega}v(t)=0.01146 .

$\overline{v}$	$t=1$	$t=7.5$	$t=15$
$0$	${\int }_{\Omega }v\left(t\right)=0.0042$	${\int }_{\Omega }v\left(t\right)\approx {10}^{-6}$	${\int }_{\Omega }v\left(t\right)\approx {10}^{-7}$
${\int }_{\Omega }{v}_0$	${\int }_{\Omega }v\left(t\right)=0.01139$	${\int }_{\Omega }v\left(t\right)=0.01146$	${\int }_{\Omega }v\left(t\right)=0.01146$
$0.6$	${\int }_{\Omega }v\left(t\right)=0.38179$	${\int }_{\Omega }v\left(t\right)=0.59967$	${\int }_{\Omega }v\left(t\right)=0.59999$

Equations

$$\{\begin{array}{ll}{\tau }_u{u}_t=-\Delta \{{ϵ}_u^2\Delta u+(1-u)(1+u)u-\alpha v-\beta v^2\}& \text{ in }\Omega \times [0,T],\\ {\tau }_v{v}_t=-\Delta \{{ϵ}_v^2\Delta v+(1-v)(1+v)v-\alpha u-2\beta uv\}-\sigma (v-\overline{v})& \text{ in }\Omega \times [0,T],\\ u(\cdot ,0)=u_0& \text{ in }\Omega ,\\ v(\cdot ,0)=v_0& \text{ in }\Omega ,\\ {\partial }_{\mathbf{n}}u={\partial }_{\mathbf{n}}\{{ϵ}_u^2\Delta u+(1-u)(1+u)u-\alpha v-\beta v^2\}=0& \text{ on }\partial \Omega \times [0,T],\\ {\partial }_{\mathbf{n}}v={\partial }_{\mathbf{n}}\{{ϵ}_v^2\Delta v+(1-v)(1+v)v-\alpha u-2\beta uv\}=0& \text{ on }\partial \Omega \times [0,T].\end{array}$$

$$F\equiv F_{\epsilon_{u},\epsilon_{v},\sigma}(u,v)=\int_{\Omega}\biggl{\{}\frac{\epsilon_{u}^{2}}{2}|\nabla u|^{2}+\frac{\epsilon_{v}^{2}}{2}|\nabla v|^{2}+W(u,v)+\frac{\sigma}{2}|(-\Delta)^{-1/2}(v-\overline{v})|^{2}\biggr{\}},$$

$$W(u,v)=\frac{(u^2-1{)}^2}{4}+\frac{(v^2-1{)}^2}{4}+\alpha uv+\beta u{v}^2.$$

$$\{\begin{array}{ll}{\tau }_u{u}_t=-\Delta w_u& \text{ in }\Omega \times [0,T],\\ w_u={ϵ}_u^2\Delta u+\varphi (u)-\alpha v-\beta v^2& \text{ in }\Omega \times [0,T],\\ {\tau }_v{v}_t=-\Delta w_v-\sigma (v-\overline{v})& \text{ in }\Omega \times [0,T],\\ w_v={ϵ}_v^2\Delta v+\varphi (v)-\alpha u-2\beta uv& \text{ in }\Omega \times [0,T],\\ u(\cdot ,0)=u_0& \text{ in }\Omega ,\\ v(\cdot ,0)=v_0& \text{ in }\Omega ,\\ {\partial }_{\mathbf{n}}u={\partial }_{\mathbf{n}}v={\partial }_{\mathbf{n}}{w}_u={\partial }_{\mathbf{n}}{w}_v=0& \text{ on }\partial \Omega \times [0,T].\end{array}$$

$$u_t\approx \frac{u^{(k+1)}-u^{(k)}}{\Delta t},$$

$$\{\begin{array}{l}{\int }_{\Omega }{\tau }_u\frac{u^{(k+1)}-u^{(k)}}{\Delta t}\phi ={\int }_{\Omega }\nabla w_u^{(k+1)}\cdot \nabla \phi ,\\ {\int }_{\Omega }{w}_u^{(k+1)}\eta ={\int }_{\Omega }\{-{ϵ}_u^2\nabla u^{(k+1)}\cdot \nabla \eta +\varphi (u^{(k+1)})\eta -\alpha v^{(k+1)}\eta -\beta (v^{(k+1)}{)}^2\eta \},\\ {\int }_{\Omega }{\tau }_v\frac{v^{(k+1)}-v^{(k)}}{\Delta t}\psi ={\int }_{\Omega }\{\nabla w_v^{(k+1)}\cdot \nabla \psi -\sigma (v^{(k+1)}-\overline{v})\psi \},\\ {\int }_{\Omega }{w}_v^{(k+1)}\zeta ={\int }_{\Omega }\{-{ϵ}_v^2\nabla v^{(k+1)}\cdot \nabla \zeta +\varphi (v^{(k+1)})\zeta -\alpha u^{(k+1)}\zeta -2\beta u^{(k+1)}{v}^{(k+1)}\zeta \},\\ u(\cdot ,0)=u_0,\\ v(\cdot ,0)=v_0.\end{array}$$

$$\{\begin{array}{l}{\int }_{\Omega }{\tau }_u\frac{u_h^{(k+1)}-u_h^{(k)}}{\Delta t}{\phi }_h={\int }_{\Omega }\nabla w_{u,h}^{(k+1)}\cdot \nabla {\phi }_h,\\ {\int }_{\Omega }{w}_{u,h}^{(k+1)}{\eta }_h={\int }_{\Omega }\{-{ϵ}_u^2\nabla u_h^{(k+1)}\cdot \nabla {\eta }_h+\varphi (u_h^{(k+1)}){\eta }_h-\alpha v_h^{(k+1)}{\eta }_h-\beta (v_h^{(k+1)}{)}^2{\eta }_h\},\\ {\int }_{\Omega }{\tau }_v\frac{v_h^{(k+1)}-v_h^{(k)}}{\Delta t}{\psi }_h={\int }_{\Omega }\{\nabla w_{v,h}^{(k+1)}\cdot \nabla {\psi }_h-\sigma (v_h^{(k+1)}-\overline{v}){\psi }_h\},\\ {\int }_{\Omega }{w}_{v,h}^{(k+1)}{\zeta }_h={\int }_{\Omega }\{-{ϵ}_v^2\nabla v_h^{(k+1)}\cdot \nabla {\zeta }_h+\varphi (v_h^{(k+1)}){\zeta }_h-\alpha u_h^{(k+1)}{\zeta }_h-2\beta u_h^{(k+1)}{v}_h^{(k+1)}{\zeta }_h\},\\ u_h(\cdot ,0)=u_{0,h},\\ v_h(\cdot ,0)=v_{0,h}.\end{array}$$

$$P_r=\{v_h\in C^0(\bar{\Omega }):v_h{|}_K\in {\mathbb{P}}_r(K),\forall K\in {\mathcal{T}}_h(\Omega )\}r=1,2,\dots$$

$$\varphi (u^{(k+1)})=-\frac{3}{2}(u^{(k)}{)}^2{u}^{(k+1)}+\frac{1}{2}(u^{(k)}{)}^3+\frac{u^{(k+1)}+u^{(k)}}{2},$$

$${\int }_a^{b}g(x)dx=(b-a)g(a)+\frac{1}{2}(b-a{)}^2{g}^{\prime }(a)+C(b-a{)}^3{g}^{\prime \prime }(\xi ).$$

$$\varphi (u^{(k+1)})=\{\begin{array}{ll}2u^{(k+1)}+u^{(k+1)}-u^{(k)}-2& \text{ if }{u}^{(k)}<-1,\\ & \\ \begin{array}{rl}& 2u^{(k+1)}+u^{(k+1)}-u^{(k)}+3u^{(k)}-(u^{(k)}{)}^3\\ & +\frac{1}{2}(u^{(k+1)}-u^{(k)})(3-3(u^{(k)}{)}^2)\end{array}& \text{ if }{u}^{(k)}\in [-1,1],\\ & \\ 2u^{(k+1)}+u^{(k+1)}-u^{(k)}+2& \text{ if }{u}^{(k)}>1.\end{array}$$

$$F(u)=\begin{dcases}\biggl{(}u^{2}+\frac{1}{4}\biggr{)}-\biggl{(}-2u-\frac{3}{4}\biggr{)}\quad&\text{ if }u<-1,\\ \biggl{(}u^{2}+\frac{1}{4}\bigr{)}-\biggl{(}\frac{3}{2}u^{2}-\frac{1}{4}u^{4}\biggr{)}\quad&\text{ if }u\in[-1,1],\\ \biggl{(}u^{2}+\frac{1}{4}\biggr{)}-\biggl{(}2u-\frac{3}{4}\biggr{)}&\text{ if }u>1.\\ \end{dcases}$$

$${\int }_{\Omega }{w}_u^{(k+1)}\psi$$

$$\displaystyle+\biggl{[}\int_{\Omega}\{-\frac{\epsilon^{2}}{2}\nabla u^{(k)}\cdot\nabla\psi-\alpha\nabla(u^{(k+1)}-u^{(k)})\cdot\nabla\psi\}\biggr{]},$$

$$\varphi (u^{(k+1)})=-2u^{(k+1)}-(u^{(k)}{)}^3+3u^{(k)},$$

$$\varphi (u^{(k+1)})=(1-u^{(k)})(1+u^{(k)})u^{(k+1)}$$

Full text

1. 1 Introduction

2. 2 The mathematical model

3. 3 Numerical approximations

4. 3.1 Semi-discretization in time

5. 3.2 Full discretization

6. 3.3 Linearization

7. 3.4 OD2, WVV, LS and EY methods: a comparison

8. 4 Fully coupled system

9. 5 Conclusions

1. Introduction

The Cahn-Hilliard equation is a fourth order nonlinear PDE introduced in [Cah61, EW59, EW58] to model the phase separation in binary systems. This equation and its variants present many practical applications, see for example [Liu+13, Tre03, Wis+08]. In the last decades, many methods and numerical techniques to treat this kind of equations have been developed. See [J+14] and [Mir17] for a review of the numerical and theoretical results respectively.

The system we study describes the spontaneous separation that occurs in a copolymer and a homopolymer blend. A copolymer consists of two monomers, while a homopolymer is made of a single component. The diffuse interface model we consider avoids the explicit treatment of the sharp interface condition between the copolymer and the homopolymer. This model has been accurately investigated in a three-dimensional setting in [E+16], using the numerical methods developed in [Eyr98] and [TN10].

In this work we study the two-dimensional, using a different numerical technique, namely Finite Element Method (FEM). In particular, we adapt some basic techniques used for the single Cahn-Hilliard equation and suitably couple them in order to deal with the fully coupled system in an efficient and reliable way. Finally we describe the behavior of the numerical solutions obtained by varying the characteristic parameters of the system. This is done in order to check that the results we obtain reflect the physical properties of the system and performing a stability analysis.

The paper is structured as follows: in section 2 we introduce the mathematical model used to describe the physical phenomenon; in section 3 we present the numerical methods for the single equation that will constitute the building blocks for our algorithms. In section 4 we describe the behavior of the numerical solutions obtained by varying the parameters of the system. At the end, in section 5, we draw some conclusions and we discuss some possible further developments.

2. The mathematical model

In this section we give a brief description of the model whose detailed derivation can be found in [E+]. The system we consider consists of a copolymer and a homopolymer blend. In particular we have two phase separations: macrophase separation and microphase separation. The first one takes place between the copolymer and the homopolymer and it is represented by the order parameter $u$, which takes values in $[-1,1]$. The ending points of this interval correspond to a homopolymer rich domain $(-1)$ and a copolymer rich domain $(+1)$. The second one is described by the order parameter $v$ and takes place inside the copolymer, between its two component, say A and B. The function $v$ also takes values in [-1,1] with ending points corresponding to A and B rich domains, respectively.

We focus on the two-dimensional case, mainly due to computational reasons and we study the problem with homogeneous Neumann boundary conditions. Notice that this is a different situation from [E+16], where only periodic boundary conditions are considered.

Let $\Omega\subset\mathbb{R}^{2}$ be a bounded Lipschitz domain and let $\tau_{u},\tau_{v},\alpha,\beta,\sigma,\epsilon_{u},\epsilon_{v},\bar{v},T\in\mathbb{R}$ with $\epsilon_{u},\epsilon_{v},\tau_{u},\tau_{v},$ $T,\sigma>0$. The model is given by the following system of coupled Cahn-Hilliard equations:

[TABLE]

where $\mathbf{n}$ is the outward normal of $\Omega$ and $u_{0},v_{0}:\Omega\to\mathbb{R}$ are given initial conditions.

Notice that the second equation of the system is a Cahn-Hilliard-Ono equation due to the presence of $\sigma(v-\bar{v})$. The time constants $\tau_{u},\tau_{v}$ control the speed of the evolution of $u$ and $v$ and, specifically, to smaller coefficients correspond a faster evolution. The parameters $\epsilon_{u},\epsilon_{v}$ are proportional to the thickness of the propagation fronts of each component and they control the size of the interface between the macrophases and the microphases. The parameter $\sigma$ is related to the bonding between A and B in the copolymer and it controls the nonlocal interactions in the second equation. More details on the nonlocal effects in phase separation problems can be found in [OK86]. From a practical point of view, values of $\sigma$ different from zero prevent the copolymer from forming large macroscopic domains. As $\sigma$ grows, finer structures and different morphologies arise. The average value of $v$ in $\Omega$ is represented by $\bar{v}$ and $T$ is the time horizon. The first coupling parameter $\alpha$ controls the interaction between the confined copolymer and the confining surface. If $\alpha=0$, the configuration is symmetric and thus $u$ has the same preference for any value of $v$. If $\alpha\neq 0$, the configuration symmetry is broken, meaning that $u$ has preference for one specific block of the copolymer. The coupling parameter $\beta$ affects the free energy depending of the value of $u$, as $v^{2}>0$.

As explained in [E+16], the state dynamics of these two coupled systems evolves as a gradient flow, up to the reaction term. More precisely, the state variables minimize the value of the following energy functional:

[TABLE]

where

[TABLE]

3. Numerical approximations

Before presenting our numerical methods for the approximate solution of (2.1), we revise some known techniques for the single Cahn-Hilliard equation. To this aim, let us set:

• •

$w_{u}:=\epsilon_{u}^{2}\Delta u+(1-u)(1+u)u-\alpha v-\beta v^{2}$

• •

$w_{v}:=\epsilon_{v}^{2}\Delta v+(1-v)(1+v)v-\alpha u-2\beta uv$

and $\phi(x)=(1-x^{2})x$. Notice that $w_{u}$ and $w_{v}$ are the chemical potentials of the Cahn-Hilliard equations.

Then (2.1) becomes

[TABLE]

3.1. Semi-discretization in time

In this section we introduce the semi-discretization in time. Let us fix $n>0$, $\Delta t=T/n$ and consider a discretization of $[0,T]$ given by $\{0=t_{0},t_{1},\dots,t_{n},t_{n+1}=T\}$, where $t_{k}=k\Delta t$ and let $u(t_{k})\cong u^{(k)}$ and the same for $v,w_{u},w_{v}$. Consider the following approximation for the time derivative of $u$

[TABLE]

and the same for $v$. By multiplying (3.1) by the test functions $\varphi,\eta,\psi,\zeta$ and integrating in space over $\Omega$, we obtain: for each $k=0,\dots,n$ find $(u^{(k)},v^{(k)},w_{u}^{(k)},w_{v}^{(k)})\in[H^{1}(\Omega)]^{4}$ s.t.

[TABLE]

for all $\varphi,\eta,\psi,\zeta\in H^{1}(\Omega)$.

3.2. Full discretization

To introduce the space discretization of (3.2) based on FEM, we project the semi-discretized problem on a finite dimensional Hilbert space $V_{h}\subset H^{1}(\Omega)$. Then the Galerkin approximation of (3.2), is: for each $k=0,\dots,n$ find $(u_{h}^{(k)},v_{h}^{(k)},w_{u,h}^{(k)},w_{v,h}^{(k)})\in V_{h}^{4}$ s.t.

[TABLE]

for all $\varphi_{h},\eta_{h},\psi_{h},\zeta_{h}\in V_{h}$, where $u_{h}^{(k)},v_{h}^{(k)},w_{u,h}^{(k)},w_{v,h}^{(k)},u_{0,h},v_{0,h}$ are the projection of $u^{(k)},v^{(k)},w_{u}^{(k)},w_{v}^{(k)},u_{0},v_{0}$ on $V_{h}$.

In order to construct such a space, we consider a triangular conforming mesh $\mathcal{T}_{h}(\Omega)$, where $h=\underset{K\in\mathcal{T}_{h}}{\max}\operatorname{diam}(K)$. Then a family of possible choices for $V_{h}$ is

[TABLE]

where $\mathbb{P}_{r}(K)$ is the space of $r$ degree polynomials defined over $K$. For a more detailed theoretical framework see, for example, Chapter 4 in [Qua17].

3.3. Linearization

In this section, we introduce suitable linear approximations for the nonlinear terms $\phi(u^{(k+1)})$ and $\phi(v^{(k+1)})$ and for the coupling terms $(v^{(k+1)})^{2}$ and $u^{(k+1)}v^{(k+1)}$.  

For the last two terms we use the following approximations:

• •

$(v^{(k+1)})^{2}\approx v^{(k)}v^{(k+1)}$,

• •

$u^{(k+1)}v^{(k+1)}\approx u^{(k)}v^{(k+1)}$.

The first one is the simplest possible choice and seems to be natural. The second choice is motivated by the fact that we expect $u$ to evolve faster then $v$, as pointed out in [E+16]. Therefore $u$ will reach a stationary condition before $v$. Then we can reasonably think that the difference $|u^{(k+1)}-u^{(k)}|$ will be smaller than $|v^{(k+1)}-v^{(k)}|$ after the initial macrophase separation; thus, with this strategy, we will achieve a better numerical approximation of the physical phenomenon.

 

 

Then the key point is to find suitable approximations for the potentials $\phi(u^{(k+1)})$ and $\phi(v^{(k+1)})$. To this aim, we propose the following strategies:

1. (1)

Optimal dissipation method (OD2) both for $u$ and $v$, i.e.

[TABLE]

cf. [FG13]. Here the authors develop the so-called optimal dissipation approach. In particular a second order in time linear approximation of the potential term $F(u)=1/4(\phi^{2}-1)^{2}$ is derived by the following Hermite quadrature formula

[TABLE]

Using this approximation in a single Cahn-Hilliard equation we get a second order linear scheme. 2. (2)

Wu-Van Zwieten-Van der Zee’s method (WVV) both for $u$ and $v$:

[TABLE]

which is discussed in [WVV14]. Here the authors consider a diffuse-interface tumor-growth system consisting of a reactive Cahn-Hilliard equation and a reaction-diffusion equation. The schemes are of the Crank-Nicolson type with a convex-concave splitting of the free-energy. The potential term is approximated employing an implicit Taylor expansion of the convex part and an explicit one of the non-convex part. The splitting considered for the potential term is

[TABLE]

Moreover, this method stabilizes the system by modifying the second and the fourth equations of (3.2) as follows:

[TABLE]

where $\alpha>0$ is a stabilization parameter. The expression in squared brackets is completely new and in addition there is a weighted contribution of $\nabla u^{(k+1)}$ which is multiplied by $\frac{1}{2}$. Clearly the same applies to $v$. According to [FG13], this method is second order accurate both in space and time in the single Cahn-Hilliard equation case. 3. (3)

Eyre’s method (EY), both for $u$ and $v$:

[TABLE]

The key point in Eyre’s work [Eyr98a] is to approximate the potential term by introducing a positive phobic numerical dissipation in the discrete energy law, to ensure the unconditional energy-stability of the scheme. 4. (4)

Linear splitting method (LS) both for $u$ and $v$; i.e.

[TABLE]

suggested in [E+16].To treat nonlinearities, the cubic term is split into the product of two terms; the first one is quadratic and related to the state of the system at the present time step and the second one is linear and related to the state of the system at the next time step.

3.4. OD2, WVV, LS and EY methods: a comparison

We first verify that, when the coupling parameters $\alpha$ and $\beta$ are set equal to [math], the computed solutions (for $u$ and $v$ taken separately) coincide with the ones of the single equation case. After that, we set nonzero values to the coupling parameters and we analyze the results obtained by the different numerical methods detailed in subsection 3.3. In the following test, the time step of each method is adjusted in order to have stable numerical solutions, in particular for OD2 $\Delta t=5\cdot 10^{-3}$, for LS $\Delta t=10^{-4}$ and for EY $\Delta t=10^{-4}$. Notice that the time step used for the three methods are different; in particular, the one for OD2 is larger than the others. The setting used for the simulations is the following:

$\Omega=(0,1)\times(0,1)$, $u_{0}=\sin(xy)$, $\tau_{u}=1$, $v_{0}=\cos(10(x-y))xy$, $\tau_{v}=1$, $\sigma=0.3$, $\overline{v}=\int_{\Omega}v_{0}\approx 0.0114559$, $\epsilon_{u}=0.05$, $\epsilon_{v}=0.05$, $\alpha=0.5$, $\beta=0.8$, $T=10$. In all the following simulations the chosen FEM space is $P_{1}$ and the mesh $\mathcal{T}_{h}$ is a $20\times 20$ grid as in Figure 1.

Notice that this parameters have been chosen to verify the methods. In Figure 2 the snapshots of the initial conditions are reported. As we can see in Table 1, $u$ and $v$ are very similar, both in the shape of the solution and in the numerical values attained. Then we choose to use the OD2 method to perform tests on parameters. Since this method achieves reasonably accurate numerical results even with large temporal steps (in some cases $0.01$), it is easier from a computational point of view to use it for tests on parameters.

We conclude that

1. (1)

OD2 method gives encouraging results: the numerical solutions seem to be physically relevant, even with "large" time steps (order $10^{-2}$). We analyze in further details the output of this method in the forthcoming section; 2. (2)

Method WVV does not give reasonable outputs: even with time steps of $10^{-6}$ and high number of nodes, the numerical solutions tend immediately to reach very high values (order $10^{6}$) and to oscillate very fast. We conclude that this approximation is not suitable for the coupled case, as this behavior is not physically acceptable; 3. (3)

Methods LS and EY give accurate results only for sufficiently small time steps (order $10^{-4}$).

4. Fully coupled system

In this section, we verify that the numerical solutions we obtained by varying the parameters respect the physical properties of the system. To achieve this purpose, we present some tests, performed with OD2 linearization method. We remark that this method allows larger time steps than EY and LS, speeding up the computations.

For each test the general setting is the following: $\Omega=(0,1)\times(0,1)$, $u_{0}=\sin(10xy)$, $v_{0}=\cos(10(x-y))xy$, $\overline{v}=\int_{\Omega}v_{0}\approx 0.0114559$, $\Delta t=0.005$, $T=15$ and the FEM space used is $P1$ with a $20\times 20$ grid. For each parameter we report in a figure the behavior of the solution at different times and for different values of the parameters. The plots of the initial conditions for all test are reported in Figure 3 and the values chosen for our tests are inspired by [E+]. In particular we vary them in the ranges considered by the authors.

**Test on $\tau_{v}$

**In this test we vary the value of $\tau_{v}$, in order to verify that its growth implies a slower evolution of the copolymer part of the system. As we can see in Table 2 the system reaches the same final configuration for both values of the parameter, while in the case $\tau_{v}=250$ it evolves slower: for $\tau_{v}=10$ the part of the copolymer in the upper-right corner has already disappered at time $t=3$, while for $\tau_{v}=250$, it is still present at time t = 6. The same happens for $\tau_{u}$.

Test on $\epsilon_{v}$

This parameter controls the separation interface thickness between the two components of the copolymer. For bigger values of $\epsilon_{v}$ the interface between the two copolymers must be thicker so that it prevents the generation of fine patterns. As we can see in Table 3, for small values of $\epsilon_{v}$, we have more complicated patterns which become coarser, according to the growth of $\epsilon_{v}$. The same happens for $\epsilon_{u}$.

**Test on $\sigma$

**The parameter $\sigma$ is related to the connectivity between the two components of the copolymer. As we can see in Table 4, as the value of $\sigma$ increases, the components of the copolymer tend no longer to form big macroareas, but they separate in more and smaller parts.

**Test on $\beta$

**This parameter controls the energetic interaction between $u$ and $v$. As we can see in Table 5, if $\beta$ decreases there is a more neat separation of the two macrophases.

Test on $\alpha$

This value controls the interaction between the confined copolymer and the confining surface, causing symmetry-breaking between microphase separated domains. If $\alpha=0$, both microphases in $v$ can reach $u$ with the same probability, but if $\alpha\neq 0$ we are changing the preference of $u$ for the positive or negative values of $v$ and thus the confining surface will change according to this preference. Indeed, as we can see in Table 6, for negative values of $\alpha$, at the interface with $u$ we have prevalence of $v=-1$, instead, with $\alpha$ positive, at the interface with $u$ we have prevalence of $v=1$.

Mass conservation

Finally, we verify that if $\overline{v}=\int_{\Omega}v_{0}$, the quantity $\int_{\Omega}v(t)$ remains constant during the evolution. On the contrary, if we set an arbitrary value of $\overline{v}$ we can see that there is no mass conservation but $\int_{\Omega}v(t)\to\overline{v}$ as $t\to+\infty$. We report the following test in Table 7 with setting: $\Omega=(0,1)\times(0,1)$, $u_{0}=\sin(10xy)$, $v_{0}=\cos(10(x-y))xy$, $\Delta t=0.005$, $T=15$ and the FEM space we use is still $P1$ on a $20\times 20$ grid.

5. Conclusions

We have proposed and tested different strategies to solve numerically a system of fully coupled Cahn-Hilliard equations. The numerical solutions we have obtained seem to give a fair description of the physical process. From the sensitivity analysis on the robustness with respect to the physical parameters, we can see that the evolution of the copolymer is confined where $u=1$. This is very important, since the solution behaves coherently with the phenomenon without imposing conditions at the interface.

Our results in the two-dimensional framework are in accordance with the ones in [E+16], even if we consider Neumann boundary conditions instead of periodic. From a numerical point of view, we notice that the linearization methods used for the single Cahn-Hilliard equation can be used also in the context of the coupled system, with the exception of the Wu-Van Zwieten-Van der Zee’s method. Moreover, the results are good even with a coarse grid in space ($20\times 20$) and with a low order (i.e $1$) FEM approximations.

This work represents a first step in developing efficient schemes for fully coupled systems of Cahn-Hilliard equations. It may be interesting to further carry out some experimental validations and to fully address the theoretical analysis. Specifically, a further numerical development could be the theoretical analysis in the coupled case of the methods we applied. From an analytical point of view, the study of the well posedness, regularity and longtime behavior of the solutions is also a topic of further research.