Global existence of weak solutions to viscoelastic phase separation: Part I Regular Case

TL;DR

This paper proves the existence of weak solutions for a complex viscoelastic phase separation model in two dimensions, combining Cahn-Hilliard and Peterlin-Navier-Stokes equations, and demonstrates solution behavior via finite element simulations.

Contribution

It establishes the mathematical existence of weak solutions for a coupled viscoelastic phase separation model with polynomial potential in 2D.

Findings

Demonstrates complex solution behavior during spinodal decomposition

Validates the model using Lagrange-Galerkin finite element method

Shows transient polymeric network structures in simulations

Abstract

We prove the existence of weak solutions to a viscoelastic phase separation problem in two space dimensions. The mathematical model consists of a Cahn-Hilliard-type equation for two-phase flows and the Peterlin-Navier-Stokes equations for viscoelastic fluids. We focus on the case of a polynomial-like potential and suitably bounded coefficient functions. Using the Lagrange-Galerkin finite element method complex behavior of solution for spinodal decomposition including transient polymeric network structures is demonstrated.

Equations

$$E_{tot}(\varphi ,q,\mathbf{C},\mathbf{u})$$

$$={\int }_{\Omega }\left(\frac{c_0}{2}{\left|\nabla \varphi \right|}^2+F(\varphi )\right)+{\int }_{\Omega }\frac{1}{2}{q}^2+{\int }_{\Omega }\left(\frac{1}{4}\mathrm{tr}\left(\mathrm{tr}(\mathbf{C})\mathbf{C}-2\ln (\mathbf{C})-\mathbf{I}\right)\right)+{\int }_{\Omega }\frac{1}{2}{\left|\mathbf{u}\right|}^2,$$

$$F(\varphi )=a{\varphi }^2(\varphi -1{)}^2,$$

$$F_{log}(\varphi )=\frac{1}{n_p}\varphi \ln (\varphi )+\frac{1}{n_s}(1-\varphi )\ln (1-\varphi )+\chi \varphi (1-\varphi ),$$

$$|F(s)|\le c_1|s{|}^p+c_2$$

$$F(s)\ge -c_7$$

$$\displaystyle\begin{split}\frac{\partial\phi}{\partial t}+\hskip 10.00002pt\mathbf{u}\cdot\nabla\phi&=\mathrm{div}\big{(}{m(\phi)\nabla\mu}\big{)}-\kappa\mathrm{div}\big{(}{n(\phi)\nabla\big{(}A(\phi)q\big{)}}\big{)}\\ \frac{\partial q}{\partial t}+\hskip 10.00002pt\mathbf{u}\cdot\nabla q&=-\frac{1}{\tau(\phi)}q+A(\phi)\Delta\big{(}A(\phi)q\big{)}-\kappa A(\phi)\mathrm{div}\big{(}{n(\phi)\nabla\mu}\big{)}\\ \frac{\partial\mathbf{u}}{\partial t}+(\mathbf{u}\cdot\nabla)\mathbf{u}&=\mathrm{div}\big{(}{\eta(\phi)\big{(}\nabla\mathbf{u}+(\nabla\mathbf{u})^{\top}\big{)}}\big{)}-\nabla p+\mathrm{div}\,{\mathbf{T}}+\nabla\phi\mu\\ \frac{\partial\mathbf{C}}{\partial t}+(\mathbf{u}\cdot\nabla)\mathbf{C}&=(\nabla\mathbf{u})\mathbf{C}+\mathbf{C}(\nabla\mathbf{u})^{\top}+\mbox{\large$\chi$}(\phi,\mathrm{tr}(\mathbf{C}))\mathbf{I}-\Phi(\phi,\mathrm{tr}(\mathbf{C}))\mathbf{C}+\varepsilon\Delta\mathbf{C}\\ \mathrm{div}\left({\mathbf{u}}\right)&=0\hskip 28.45274pt\mathbf{T}=\mathrm{tr}\left({\mathbf{C}}\right)\mathbf{C}\hskip 28.45274pt\mu=-c_{0}\Delta\phi+f(\phi)\end{split}$$

$${(\varphi ,q,\mathbf{u},\mathbf{C})|}_{t=0}=({\varphi }_0,q_0,{\mathbf{u}}_0,{\mathbf{C}}_0),{{\partial }_n\varphi |}_{\partial \Omega }={{\partial }_n\mu |}_{\partial \Omega }={{\partial }_{n}q|}_{\partial \Omega }=0,{\mathbf{u}|}_{\partial \Omega }=0,{{\partial }_n\mathbf{C}|}_{\partial \Omega }=0.$$

$$\mathbf{T}=\mathrm{tr}\left(\mathbf{C}\right)\mathbf{C}.$$

$$0<{\tau }_1\le \tau (s)\le {\tau }_2,0<A_1\le A(s)\le A_2,0<{\eta }_1\le \eta (s)\le {\eta }_2\text{ for all }s\in \mathbb{R}.$$

$$\mbox{\large$\chi$}(\phi,\mathbf{C})=h(\phi)\mathrm{tr}\left({\mathbf{C}}\right),\hskip 20.00003pt\Phi(\phi,\mathbf{C})=h(\phi)\mathrm{tr}\left({\mathbf{C}}\right)^{2}.$$

$$\mbox{\large$\chi$}(\phi,\mathbf{C})=h(\phi)\mathrm{tr}\left({\mathbf{C}}\right)^{\alpha-1},\hskip 20.00003pt\Phi(\phi,\mathbf{C})=h(\phi)\mathrm{tr}\left({\mathbf{C}}\right)^{\alpha},\;\alpha\in[1,6].$$

$$0<h_1\le h(s)\le h_2,0\le m_1\le m(s)\le m_2\text{ and }0\le n_1<n(s)\le n_2\text{ for all }s\in \mathbb{R}.$$

$$L^{2}_{\text{div}}(\Omega):=\overline{C_{0,\text{div}}^{\infty}(\Omega)}^{\left\|\cdot\right\|_{2}},\hskip 20.00003ptL^{2}_{0}(\Omega):=\left\{u\in L^{2}(\Omega)\mathrel{\bigg{|}}\int_{\Omega}u\,\mathrm{d}x=0\right\},$$

$$V:=H_{0,\text{div}}^1(\Omega {)}^2,H:=L_{\text{div}}^2(\Omega {)}^2,W:=H^1(\Omega {)}^{2\times 2}.$$

$${\Vert u\Vert }_4$$

$${\Vert v\Vert }_4$$

$${\Vert u\Vert }_r$$

$${\Vert v\Vert }_{\infty }$$

$$\int_{D}\big{(}u_{n}v_{n}-uv\big{)}g\longrightarrow 0,n\rightarrow\infty\;\text{ for all }g\in L^{2}(D).$$

$$\left|\int_{D}\big{(}u_{n}v_{n}-uv\big{)}g\right|\leq\left|\int_{D}gv_{n}\big{(}u_{n}-u\big{)}\right|+\left|\int_{D}gu\big{(}v_{n}-v\big{)}\right|$$

$$\left\{f\in L^p(0,T;X),\frac{\mathrm{d}f}{\mathrm{d}t}\in L^q(0,T;Z)\right\}\hookrightarrow L^p(0,T;Y)$$

$$\left\{f\in L^{\infty }(0,T;X),\frac{\mathrm{d}f}{\mathrm{d}t}\in L^p(0,T;Z)\right\}\hookrightarrow C([0,T];Y)$$

$$\left\{f\in L^2(0,T;X),\frac{\mathrm{d}f}{\mathrm{d}t}\in L^2(0,T;Y)\right\}\hookrightarrow C(0,T;[X,Y{]}_{\frac{1}{2}}).$$

$${\Vert \mathrm{tr}\left(\mathbf{D}\right)\Vert }_p^p:={\int }_{\Omega }{\left(\sum _{i=1}^d{\mathbf{D}}_{ii}\right)}^p\mathrm{d}x\le d^{p-1}{\int }_{\Omega }\sum _{i,j=1}^d\left|{\mathbf{D}}_{ij}^p\right|\mathrm{d}x=:{\Vert \mathbf{D}\Vert }_p^p.$$

$$\mathbf{b}(\mathbf{u},\mathbf{v},\mathbf{w})\equiv {\int }_{\Omega }(\mathbf{u}\cdot \nabla )\mathbf{v}\cdot \mathbf{w}\mathrm{d}x\text{ and }\mathbf{B}(\mathbf{u},\mathbf{C},\mathbf{D})\equiv {\int }_{\Omega }(\mathbf{u}\cdot \nabla )\mathbf{C}:\mathbf{D}\mathrm{d}x$$

$$\mathbf{b}(\mathbf{u},\mathbf{u},\mathbf{v})$$

$$\mathbf{B}(\mathbf{u},\mathbf{C},\mathbf{D})$$

$$\left|\mathbf{b}(\mathbf{u},\mathbf{w},\mathbf{v})\right|$$

$$\left|\mathbf{B}(\mathbf{u},\mathbf{C},\mathbf{D})\right|$$

$${\Vert \nabla \mathbf{u}\Vert }_2\le \sqrt{2}{\Vert \mathrm{D}\mathbf{u}\Vert }_2,$$

$${\Vert \sqrt{\eta }\nabla \mathbf{u}\Vert }_2\le \sqrt{2}{\Vert \sqrt{\eta }\mathrm{D}\mathbf{u}\Vert }_2$$

$$\mathrm{D}\mathbf{u}=\frac{1}{2}\left(\nabla \mathbf{u}+\nabla {\mathbf{u}}^T\right).$$

$${\int }_M{f}_m\to {\int }_{M}f.$$

Full text

Institute of Mathematics, Johannes Gutenberg-University Mainz

Staudingerweg 9, 55128 Mainz, Germany

abrunk@uni-mainz.de

lukacova@uni-mainz.de

Keywords: two-phase flows, non-Newtonian fluids, Cahn-Hilliard equation, Peterlin viscoelastic model, Navier-Stokes equation, phase separation, Flory-Huggins potential

1 Introduction

Phase separation of binary fluids is a fundamental process in soft matter physics. For Newtonian fluids this phenomenon is well-understood. In this case typically the so-called H model is used, which consists of the conservation of mass and momentum combined with a high-order nonlinear convection-diffusion equation for the phase variable $\phi$. The evolution of the phase variable is mainly governed by a gradient flow for the free energy functional and is deeply connected to the thermodynamics of the involved process. To avoid the formation of discontinuous interfaces the free energy functional is supplied with a penalty term for the gradient of $\phi$. From mathematical and numerical point of view, this allows describing topological changes without tracking interfaces. However, if one of the multiphase fluids happens to be a polymer, the physics becomes much more involved. Thus, one has to consider non-Newtonian rheology together with multiphase effects. In this context Tanaka [51] introduced a new concept of viscoelastic phase separation that governs phase-separation of a dynamically asymmetric mixture, which is composed of fast and slow components. In dynamically asymmetric mixtures the phase separation leads to new interesting structures, such as transient formation of network-like structures of a slow-component-rich phase and its volume shrinking. Unfortunately the Tanaka model was not consistent with the second law of thermodynamics. In [55] Zhou, Zhang, E derived a thermodynamically consistent model for viscoelastic phase separation. The key ingredient of both models is an additional viscoelastic stress tensor which is connected to the velocity difference of both multiphase fluids.

The phase-field dynamics is described by a modified Cahn-Hilliard equation. In literature we can find already a variety of analytical studies for the multiphase flow governed by the Cahn-Hilliard equation. Typically, models with a polynomial-type potential and strictly positive diffusion coefficient are considered, see, e.g. [21, 22, 53]. We will refer to such a situation as the regular case, see Section 2.1 for further details. Other models available in literature deal with constant mobility functions and singular potentials [4, 43] or degenerate mobility and singular potentials [5, 17, 13, 22, 31].

Taking also hydrodynamic effects into account yields the Navier-Stokes-Cahn-Hilliard system which has been studied in, e.g., [4, 9]. Including the viscoelastic effects leads to the Navier-Stokes-Cahn-Hilliard-Oldroyd-B system, for which the well-posedness of strong solutions has been studied in [14]. Recent results focus on the so-called non-local Cahn-Hilliard equations, see, e.g., [26, 27]. The non-local diffusion term has a better structure and allows one to obtain better regularity properties and the maximum principles for degenerate mobility functions and singular potentials [16, 25, 35]. As far as we are aware, it is not clear how to obtain the local Cahn-Hilliard equation from the non-local one, see [42]. In [6, 8, 23, 36] time evolution of compressible polymer mixtures has been studied and the existence of global weak solutions was shown.

Although the model of Zhou, Zhang and E was successfully used for numerical simulations of spinodal decomposition, see, e.g. [55, 50], its mathematical analysis was open. The purpose of this paper is to fill this gap and investigate the existence of global weak solutions of a related viscoelastic phase separation model. More precisely, our model is similar to the viscoelastic phase separation model of Zhou et al., but we consider the diffusive Peterlin model for the time evolution of the elastic conformation tensor instead of the classical Oldroyd-B model. The Peterlin model [38] can be viewed as a nonlinear generalization of the classical Oldroyd-B model. Thus, our existence result also applies to the model of Zhou et al. with additional small diffusive terms in the Oldroyd-B equations. Following Barrett and Süli [1] we note that the diffusive term appearing in the evolution equation for the conformation tensor should not be seen as a regularizing term but rather an outcome of physical modelling allowing a heterogeneous fluid velocity around a dumbbell. Indeed, as investigated by Málek et al. [45] the stress diffusion term can be interpreted either as a consequence of a nonlocal energy storage mechanism or as a consequence of a nonlocal entropy production mechanism, see also [18, 49] for related studies.

The main features of the proposed approach can be summarized as follows.

• •

The existence proof is based on the energy method and compactness arguments in order to pass to limits in nonlinear terms. It should be however pointed out that the viscoelastic phase separation contains cross-diffusion terms that require careful treatment since the energy method fails in such contexts in general.

• •

Extensive numerical simulations, from which a representative choice is presented in Section 9, confirm a good agreement with physical experiments of dynamical asymmetric mixtures, see [51]. In particular, formation of a new phenomena, such as transient network-like polymeric structures and volume shrinking has been observed, too. Mathematically speaking, the cross-diffusion terms in the $(\phi,q)$-subsystem, cf. (2.6), yield a realistic model for structure formation effects, such as transient network-like polymeric structures.

• •

In this paper we consider a regular case with a polynomial-like potential function and suitably bounded mobilities. Based on the current results we will study in our subsequent work [12] the existence of a global weak solution for the case with logarithmic Flory-Huggins potential and degenerate mobilities with possibly different rates. We refer a reader to [22] where similar situation with one mobility function has been studied.

The paper is organized in the following way. In Section 2 we present the mathematical model for viscoelastic phase separation. The weak solution to our viscoelastic phase separation model is introduced in Section 3 and the main result on the existence of global weak solutions is formulated in Section 4. Sections 5-8 are devoted to the proof of our main result. In Section 9 we propose a Lagrange-Galerkin finite element method. Further, we present several numerical simulations of spinodal decomposition. Numerical experiments indicate that our numerical method is energy-stable and mass conservative.

2 Mathematical Model

The viscoelastic phase separation can be described by a coupled system consisting of the Cahn-Hilliard equation for phase field evolution, the Navier-Stokes equation for fluid flow and the time evolution of the viscoelastic conformation tensor. A classical approach to model the evolution of interfaces is the diffusive interface theory.

In the recent paper by Zhou et al. [55] the classical interface theory has been combined with time evolution of a viscoelastic fluid in order to model viscoelastic phase separation. The total energy of the polymer-solvent mixture consists of the mixing energy between the polymer and the solvent, the bulk stress energy, the elastic energy and the kinetic energy.

[TABLE]

where $\phi$ denotes the volume fraction of polymer molecules, $q\mathbf{I}$ the bulk stress arising from polymeric interactions, $\mathbf{C}$ the viscoelastic conformation tensor and $\mathbf{u}$ the volume averaged velocity consisting of solvent and polymer velocity. Furthermore, $c_{0}$ is a positive constant controlling the interface width. In the present work we will work with the Ginzburg-Landau potential,

[TABLE]

which is typically used as a regular approximation of the (more physically relevant) logarithmic Flory-Huggins potential

[TABLE]

where $n_{p},n_{s}$ stand for the molecular weights of the polymer and solvent, respectively. $\chi$ is the so-called Flory interaction parameter which describes the interaction between the two components and is proportional to the inverse temperature.

In the forthcoming analysis it will be enough to assume the potential behaves like a polynomial. More precisely, we assume that the potential $F\in C^{2}(\mathbb{R})$ and there are constants $c_{i}>0$ for $i\in\{1,\ldots,7\}$ and $c_{8}\geq 0$ that the following holds true,

[TABLE]

Further, we will set $F^{\prime}(s)=f(s)$.

Following [55] we can derive a multiphase model for the viscoelastic phase separation, see also Appendix, where a brief derivation in the context of GENERIC [28, 29] and of [45, 46] is presented.

\displaystyle\begin{split}\frac{\partial\phi}{\partial t}+\hskip 10.00002pt\mathbf{u}\cdot\nabla\phi&=\mathrm{div}\big{(}{m(\phi)\nabla\mu}\big{)}-\kappa\mathrm{div}\big{(}{n(\phi)\nabla\big{(}A(\phi)q\big{)}}\big{)}\\ \frac{\partial q}{\partial t}+\hskip 10.00002pt\mathbf{u}\cdot\nabla q&=-\frac{1}{\tau(\phi)}q+A(\phi)\Delta\big{(}A(\phi)q\big{)}-\kappa A(\phi)\mathrm{div}\big{(}{n(\phi)\nabla\mu}\big{)}\\ \frac{\partial\mathbf{u}}{\partial t}+(\mathbf{u}\cdot\nabla)\mathbf{u}&=\mathrm{div}\big{(}{\eta(\phi)\big{(}\nabla\mathbf{u}+(\nabla\mathbf{u})^{\top}\big{)}}\big{)}-\nabla p+\mathrm{div}\,{\mathbf{T}}+\nabla\phi\mu\\ \frac{\partial\mathbf{C}}{\partial t}+(\mathbf{u}\cdot\nabla)\mathbf{C}&=(\nabla\mathbf{u})\mathbf{C}+\mathbf{C}(\nabla\mathbf{u})^{\top}+\mbox{\large\chi}(\phi,\mathrm{tr}(\mathbf{C}))\mathbf{I}-\Phi(\phi,\mathrm{tr}(\mathbf{C}))\mathbf{C}+\varepsilon\Delta\mathbf{C}\\ \mathrm{div}\left({\mathbf{u}}\right)&=0\hskip 28.45274pt\mathbf{T}=\mathrm{tr}\left({\mathbf{C}}\right)\mathbf{C}\hskip 28.45274pt\mu=-c_{0}\Delta\phi+f(\phi)\end{split}

(2.6)

Here we have denoted by $\mu$ the chemical potential. This allows us to write (2.6) in the saddle point form which is more convenient for our further study. System (2.6) is formulated on $(0,T)\times\Omega$, where $\Omega\subset\mathbb{R}^{2}$ is a bounded domain with a Lipschitz-continuous boundary. Our model is equipped with the following initial and boundary conditions

[TABLE]

The viscoelastic stress tensor is given by

[TABLE]

The functions $m(\phi),n(\phi)$ stand for mobility functions, \tau(\phi),\mbox{\large\chi}(\phi,\mathrm{tr}\left({\mathbf{C})}\right)^{-1},\Phi(\phi,\mathrm{tr}\left({\mathbf{C}}\right))^{-1} describe generalized relaxation times, $A(\phi)$ the bulk modulus, $\eta(\phi)$ the viscosity, $\varepsilon$ the viscoelastic diffusion parameter and $\kappa$ a positive constant. These functions will be specified in Section 2.1.

2.1 Assumptions

In what follows, we will assume that all of the coefficient functions are continuous, positive and bounded, i.e.

[TABLE]

Additionally $A(s)$ is a $C^{1}-$function and $\left\|A^{\prime}\right\|_{\infty}\leq A_{2}^{\prime}.$ For the generalized relaxation times $\chi$ and $\Phi$ we choose analogously as in [38]

[TABLE]

Let us point out that it may be possible to generalize our results similarly as in [37] for

[TABLE]

We further assume that the functions $h,m,n$ are continuous, positive and bounded,

[TABLE]

The case $m_{1},n_{1}>0$ is called regular case. In our future work [12] we extend the present study to the so-called degenerate case, where $m_{1}=n_{1}=0.$

2.2 Preliminaries

In this section we introduce suitable notation and required analytical tools. For the standard Lebesgue spaces $L^{p}(\Omega)$ the norm is denoted by $\left\|\cdot\right\|_{p}$. Further, we denote the space of divergence free and mean free functions by the

[TABLE]

respectively. We use the standard notation for the Sobolev spaces and set

[TABLE]

The space $V$ is equipped with the norm $\left\|\cdot\right\|_{V}:=\left\|\nabla\cdot\right\|_{2}$. We denote the norms of the corresponding Bochner space $L^{p}(0,T;L^{q}(\Omega))$ by $\left\|\cdot\right\|_{L^{p}(L^{q})}.$

In what follows, we state several inequalities that will be used later.

Proposition 2.1** ([7],[34]).**

Let $\Omega\subset\mathbb{R}^{2}$ be a bounded smooth domain. Then the following inequalities hold true

[TABLE]

Lemma 2.2**.**

Suppose that there are two sequences $\{u_{n}\},\{v_{n}\}$ and a bounded domain $D$ with the following properties

$u_{n}\rightarrow u$* a.e in $D$ and $\left\|u_{n}\right\|_{\infty}\leq c<\infty$ for all $n$* 2. 2.

$v_{n}\rightharpoonup v$* in $L^{2}(D)$.*

Then the product $u_{n}v_{n}$ converges weakly to $uv$ in $L^{2}(D)$.

Proof.

Weak convergence of $u_{n}v_{n}$ in $L^{2}(D)$ is defined by

[TABLE]

Splitting the integral we obtain

[TABLE]

The first integral can be controlled by $\left\|gu_{n}-gu\right\|_{2}\left\|v_{n}\right\|_{2}$ which goes to zero due to the boundedness of $v_{n}$ in $L^{2}(D)$ and the strong convergence of $gu_{n}$ in $L^{2}(D)$. This strong convergence can be achieved by the generalized dominated convergence theorem. Since $gu\in L^{2}(D)$, the second integral tends to zero due to the weak convergence of $v_{n}$. Note that in the proof we have also shown that $gu_{n}$ converges strongly in $L^{2}(D)$. ∎

We will frequently use the interpolation spaces and their norms.

Lemma 2.3** ([9]).**

Let $X\subset Y\subset Z$ three Hilbert spaces, and suppose that the embedding of $X$ in $Y$ is compact. Then:

1. i)

For any $p,q\in(1,\infty)$ the embedding

[TABLE]

is compact. 2. ii)

For any $p>1$ the embedding

[TABLE]

is compact. 3. iii)

The following continuous embeddings holds

[TABLE]

The space $[X,Y]_{\frac{1}{2}}$ is an interpolation space, see [9], and we will use this result for the spaces $[H^{2},H^{-1}]_{\frac{1}{2}}=H^{1/2}$.

Proposition 2.4**.**

For a matrix valued function $\mathbf{D}\in\mathbb{R}^{d\times d}$ and $p\geq 2$ we have

[TABLE]

For symmetric positive definite matrices both norms are equivalent.

The norm $\left\|\mathbf{D}\right\|_{p}$ is the so-called Frobenius norm. We denote by $\mathbf{C}:\mathbf{D}=\mathrm{tr}\left({\mathbf{C}\mathbf{D}^{\top}}\right)$ the Frobenius inner product.

Proposition 2.5** ([53], [38]).**

For any open set $\Omega\subset\mathbb{R}^{2}$ it holds that the forms

[TABLE]

are continuous and trilinear on $V\times V\times V$ and $V\times W\times W$, respectively. Further the following properties hold

[TABLE]

Lemma 2.6** (Korn inequality, [33]).**

Let $\Omega$ be an open and regular set of $\mathbb{R}^{2}$ and $\mathbf{u}\in H^{1}_{0}(\Omega)$ a vector field on $\Omega$, then we have

[TABLE]

with the equality if $\mathrm{div}\left({\mathbf{u}}\right)=0$. Suppose $\eta$ is a continuous function satisfying (2.7) then the inequality

[TABLE]

holds. Here $\mathrm{D}\mathbf{u}$ denotes the deformation tensor which is the symmetric part of $\nabla\mathbf{u}$, i.e.

[TABLE]

For the limiting process we will need the following lemma which is a consequence of the Vitali theorem.

Lemma 2.7** (Vitali theorem [24]).**

Let $M\subset\mathbb{R}^{d}$ be a measurable and bounded set. Let the sequence $\{f_{m}\}_{m\in\mathbb{N}}$ be uniformly bounded in $L^{q}(M)$ for $q>1$. Finally, let $f_{m}\to f$ a.e. in $M$ for some $f\in L^{q}(M)$. Then

[TABLE]

Lemma 2.8** ([9]).**

If $\phi\in H^{s+2}(\Omega)$, then $\|\phi-(\phi)_{\Omega}\|_{H^{s+2}}\leq C\|\Delta\phi\|_{H^{s}}$. Here $(\phi)_{\Omega}:=\frac{1}{\left|\Omega\right|}\int_{\Omega}\phi\,\mathrm{d}x$ denotes the mean value.

Lemma 2.9** ([8]).**

Let $\mathbf{D}\in H^{2}(\Omega)^{m\times m}\cap C^{1}(\overline{\Omega})^{m\times m},m\in\mathbb{N}$, be a symmetric matrix function, which is uniformly positive definite on $\overline{\Omega}$ and satisfies homogeneous Neumann boundary conditions, then

[TABLE]

3 Weak solution

In this section we will introduce the weak solutions to our viscoelastic phase separation model (2.6) and derive the corresponding energy estimates.

Definition 3.1**.**

Let the initial data $\left(\phi_{0},q_{0},\mathbf{u}_{0},\mathbf{C}_{0}\right)\in H^{1}(\Omega)\times L^{2}(\Omega)\times H\times L^{2}(\Omega)^{2\times 2}$. The quintuple $(\phi,q,\mu,\mathbf{u},\mathbf{C})$ is called a corresponding weak solution of (2.6), if

[TABLE]

[TABLE]

and

[TABLE]

for any test function $(\psi,\zeta,\xi,\mathbf{v},\mathbf{D})\in H^{1}(\Omega)^{3}\times V\times W$ and almost every $t\in(0,T).$ For the initial data it holds that $\left(\phi(0),q(0),\mathbf{u}(0),\mathbf{C}(0)\right)=\left(\phi_{0},q_{0},\mathbf{u}_{0},\mathbf{C}_{0}\right)$.

Theorem 3.2**.**

Let $(\phi,q,\mu,\mathbf{u},p,\mathbf{C})$ be a smooth solution of (2.6) and let the initial datum $\mathbf{C}_{0}$ be a symmetric positive definite matrix. Further assume that $\mathbf{C}$ is uniformly positive definite. Then the free energy given by (2.1) satisfies

[TABLE]

thus the energy is non-increasing in time.

Proof.

We will differentiate every term of (2.6) with respect to time and show the energy dissipation

[TABLE]

Summing up the above equations we obtain

[TABLE]

Using assumptions (2.7), (2.9), Lemma 2.9 and the fact that for a symmetric positive definite matrix, cf. [44]

[TABLE]

we have

[TABLE]

This inequality can also be obtained by testing the weak formulation (3.1) with

[TABLE]

respectively. ∎

In addition to the physical energy (2.6) also fulfills an additional energy estimate, which does not rely on the positive definiteness of the conformation tensor.

Lemma 3.3**.**

Each smooth solution of (2.6) satisfies the following energy inequality

[TABLE]

The integrated version reads

[TABLE]

Proof.

We mainly repeat the calculation from the proof of Theorem 3.2 and use the assumptions (2.7), (2.9). The inequality (3.2) can be obtained by testing (3.1) with $(\mu,q,\frac{\partial\phi}{\partial t},\mathbf{u},\frac{1}{2}\mathbf{C})$, respectively. The integrated version is obtained by integrating over the time from [math] to $t$. We only present the calculations for the viscoelastic contributions.

[TABLE]

∎

4 Main Results

In this section we formulate the main result on the existence of global weak solutions to the viscoelastic phase separation system (2.6).

Theorem 4.1**.**

Under assumptions (2.5),(2.7) - (2.9) with $m_{1}>0$ and $\kappa\in[0,1)$ there exists a weak solution of the viscoelastic phase separation model (2.6) in the sense of Definition 3.1 that satisfies the energy inequality in the form (3.3).

Remark 4.2**.**

• •

The results hold also if we substitute \mathrm{div}\big{(}{\varepsilon(\phi)\nabla\mathbf{C}}\big{)} instead of $\varepsilon\Delta\mathbf{C}$, where $\varepsilon(\phi)$ is a continuous function such that $0<\varepsilon_{1}\leq\varepsilon(\phi)\leq\varepsilon_{2}$.

• •

Following [9] it is also possible to show that $\phi\in L^{2}(0,T;H^{3}(\Omega))$.

• •

A similar version of the proof can be applied in three space dimensions, if the viscoelastic effects are not considered, i.e. $(\ref{eq:full_model})_{4}$ is neglected.

Remark 4.3**.**

The proof will be realized in Sections 5 - 8 and consists of the following steps.

We introduce a Galerkin approximation of the system in suitable finite-dimensional subspaces. 2. 2)

We derive the energy inequality for the approximate solutions and by means of the Gronwall inequality we obtain a priori estimates which are independent of the local existence time. Hence we prove existence of approximate solutions for an arbitrarily fixed $T>0$. 3. 3)

In order to pass to the limit in nonlinear terms, we apply the Lions-Aubin lemma to obtain the required compact embeddings. 4. 4)

We pass to the limit in the Galerkin approximation of (3.1). 5. 5)

We pass the limit in the energy inequality to obtain (3.3).

5 Galerkin Approximation

Let $\psi_{j},\mathbf{v}_{j},\mathbf{D}_{j},j=1,\ldots,\infty,$ be smooth basis functions of

[TABLE]

respectively. Where $\mathbf{D}_{j}$ is subjected to homogeneous Neumann boundary conditions and $\mathbf{v}_{j}$ is divergence-free and subjected to homogeneous Dirichlet boundary conditions. Furthermore, $\psi_{j}$ are eigenfunctions of $-\Delta$ with the Neumann boundary conditions. We define the $m$-th Galerkin approximation via

[TABLE]

We use the following notation

[TABLE]

The Galerkin approximation of (3.1) then reads

[TABLE]

The system (2.6) is the coupled nonlinear system of ordinary differential equations for $\phi_{m}(t),$ $q_{m}(t),\mathbf{u}_{m}(t),\mathbf{C}_{m}(t)$. Since all parameters, $\tau,m,n,A,h$ as well as $f$ are continuous functions of $\phi_{m}$, we can show via the Peano theorem that we have a solution of (5.1) till a time $t_{m}>0$. This local existence times depends on the data and $m$. In the next section we will derive a priori estimates and show global existence of the Galerkin approximation.

6 A priori Estimates

6.1 Energy inequality

By multiplying system (5.1) with $(\theta_{jm},\zeta_{jm},g_{jm},G_{jm},\lambda_{jm}^{\prime})$ and summing over all $j$ leads to

[TABLE]

Applying the assumptions (2.7), (2.9) we find

[TABLE]

The term on the right hand side can be bounded by the Gronwall inequality since $h_{2}$ is a constant. Consequently, we derive the following estimates

[TABLE]

which directly implies that we can continuously prolongate the solution of the system of ordinary differential equations (5.1) up to $t_{m}=T$. Furthermore, we have that $\mathbf{T}_{m}\in L^{2}(0,T;L^{2}(\Omega))$,

[TABLE]

We observe that testing $(\ref{eq:disc_weak})_{1}$ with $\psi=1$ yields $\frac{\mathrm{d}}{\mathrm{dt}}\int_{\Omega}\phi_{m}(t)\,\mathrm{d}x=0$ . This implies that the mean value of $\phi_{m}$ is constant in time. Using the Poincaré inequality we find $\phi_{m}\in L^{\infty}(0,T;H^{1}(\Omega))$, because

[TABLE]

The same can be done for $\left(\mu_{m}\right)_{\Omega}$ by testing $(\ref{eq:disc_weak})_{5}$ with $\psi=1/\left|\Omega\right|$. Using assumptions (2.5) we find

[TABLE]

Thus we obtain $\mu_{m}\in L^{2}(0,T;H^{1}(\Omega))$ via the Poincaré inequality.

We want to show that $\phi_{m}\in L^{2}(0,T;H^{2}(\Omega))$. By testing $(\ref{eq:disc_weak})_{5}$ with $-\Delta\phi_{m}(t)$ and using (2.5) we find

[TABLE]

Integrating the last inequality over time we obtain that $\Delta\phi_{m}$ has the same regularity as $\nabla\mu_{m}$, i.e. $L^{2}(0,T;L^{2}(\Omega))$, since $\phi_{m}\in L^{\infty}(0,T;H^{1}(\Omega))$. Lemma 2.8 together with the Agmon inequality (2.13) implies that $\phi_{m}\in L^{2}(0,T;L^{\infty}(\Omega))$.

We need to derive an estimate on $\nabla q_{m}$ using the estimate of \nabla\big{(}A(\phi_{m})q_{m}\big{)}.

[TABLE]

Due to the estimate (6.3) and (2.7) we can conclude that $\nabla q_{m}\in L^{2}(0,T;L^{2}(\Omega))$.

7 Compact Embeddings

In order to pass to the limit in the nonlinear terms we need to derive further estimates via compact embeddings to obtain a strong convergence. To this end, we define firstly suitable operators and rewrite the weak formulation as operator equations.

7.1 Cahn-Hilliard part

For the Cahn-Hilliard part $(\ref{eq:full_model})_{1}$ we have the following operators

[TABLE]

The corresponding operator equation for $(\ref{eq:weak_sol_sys})_{1}$ reads

[TABLE]

where

[TABLE]

We obtain from (7.1) - (7.5) that $\phi^{\prime}\in L^{2}(0,T;H^{-1}(\Omega))$. Due to the Lions-Aubin lemma we have the strong convergence of $\phi_{m}\rightarrow\phi$ in $L^{2}(0,T;H^{1}(\Omega))$. Clearly, the above estimates hold for any $\phi_{m}$ satisfying (7.1), thus $\left\|\phi_{m}^{\prime}\right\|_{L^{2}(H^{-1})}\leq c$.

7.2 Bulk stress

For the bulk stress equation $(\ref{eq:full_model})_{2}$ we consider

[TABLE]

and rewrite the weak formulation for $q,\;(\ref{eq:weak_sol_sys})_{2},$ as the operator equation

[TABLE]

First, let us recall that

[TABLE]

where $r=ak,s=a^{\star}k/2$ with $\frac{1}{a}+\frac{1}{a^{\star}}=1$ with $a,a^{\star}\geq 1$. We can choose $a=3/2,a^{\star}=3$ to obtain $k=4/3$.

[TABLE]

Further, we have

[TABLE]

where $a=3/2,a^{\star}=3$ and $k=4/3$.

[TABLE]

The estimates (7.6) - (7.12) imply that $q^{\prime}\in L^{4/3}(0,T;H^{-1}(\Omega))$. Analogously as in (7.1) due to the Lions-Aubin lemma we obtain the strong convergence of $q_{m}\rightarrow q$ in $L^{2}(0,T;L^{p}(\Omega))$ for every $p\in(1,\infty)$.

7.3 The Navier-Stokes part

We continue with the fluid equations $(\ref{eq:full_model})_{3}$ and define the following operators

[TABLE]

We can rewrite the velocity equation $(\ref{eq:weak_sol_sys})_{3}$ in the following operator form

[TABLE]

where

[TABLE]

Finally, (7.13) - (7.17) imply $\mathbf{u}^{\prime}\in L^{2}(0,T;V^{*})$ and by means of the Lions-Aubin lemma we get the compact embedding into $L^{2}(0,T;L^{4}_{\text{div}}(\Omega))$.

7.4 Conformation tensor

Finally, we define the operators for the viscoelastic part $(\ref{eq:weak_sol_sys})_{4}$ of our model

[TABLE]

[TABLE]

where

[TABLE]

We need to find $k,a,a^{\star}$ such that $2ka\leq 4,ka^{\star}\leq 4,\frac{1}{a}+\frac{1}{a^{\star}}=1$ and $k>1$. Choosing $a=3/2,a^{\star}=3$ we find that $k\leq 4/3$. Combining (7.18) - (7.24) we obtain that $\mathbf{C}^{\prime}\in L^{4/3}(0,T;W^{\star})$. Due to the Lions-Aubin we obtain the compact embedding of $\mathbf{C}_{m}$ into $L^{2}(0,T;L^{p}(\Omega))$ for every $p\in(1,\infty)$.

7.5 Convergence of subsequences

Due to the uniform boundedness and the compact embeddings derived above we get the following convergence results for suitable subsequences

[TABLE]

Applying the weak continuity results from [52] or the Lions-Aubin-Simons Lemma 2.3 we obtain also continuity in time, i.e.

[TABLE]

8 Limit passing

Our aim now is to pass to the limit as $m\to\infty$ in the Galerkin approximation (5.1). To this end we multiply (5.1) by a time test function $\varphi(t)\in L^{\infty}(0,T)$ and integrate over $(0,T)$. Then we pass to the limit in each term and show that the limits identified in (7.25) satisfy the weak formulation (3.1). We will concentrate only on nonlinear terms since the linear ones follow directly from the weak convergences in (7.25).

8.1 The Cahn-Hilliard part

We first consider the terms from the Cahn-Hilliard part $(\ref{eq:disc_weak})_{1}$. We start with the main operator of the Cahn-Hilliard part

[TABLE]

If we know that $m(\phi_{m})\nabla\mu_{m}$ converges weakly to $m(\phi)\nabla\mu$ in $L^{2}(0,T;L^{2}(\Omega))$ then it is easy to see that the integral (8.1) vanishes for $m\to\infty$. First it follows from (7.25) that $\phi_{m}$ converges strongly to $\phi$ in $L^{2}(0;T;L^{2}(\Omega))$ and therefore it converges almost everywhere in $(0,T)\times\Omega$.

Combining this with the continuity of $m$ we obtain convergence almost everywhere in $(0,T)\times\Omega$ of $m(\phi_{m})$ to $m(\phi)$. Due to assumption (2.9) $m(\phi_{m})$ is bounded in $L^{\infty}((0,T)\times\Omega)$. By applying Lemma 2.2 we obtain the weak convergence of $m(\phi_{m})\nabla\mu_{m}$ to $m(\phi)\nabla\mu$ in $L^{2}(0,T;L^{2}(\Omega))$.

Now we consider the convective term

[TABLE]

The first term goes to zero due to the strong convergence of $\mathbf{u}_{m}$ in $L^{2}(0,T;L_{\text{div}}^{4}(\Omega))$. The second term converges due to the strong convergence of $\nabla\phi_{m}$ in $L^{2}(0,T;L^{2}(\Omega))$.

Now we consider the nonlinear coupling term with the bulk stress.

[TABLE]

Analogously to (8.1) we can treat this term by using Lemma 2.2. Here we explain why the weak limit of \nabla\big{(}A(\phi_{m})q_{m}\big{)} is \nabla\big{(}A(\phi)q\big{)}. To this end, we calculate the gradient

[TABLE]

The first term converges weakly to \big{(}A(\phi)\nabla q\big{)} by analogous arguments as in (8.1). For the second term we recall that $\nabla\phi_{m}q_{m}$ is bounded in $L^{2}(0,T;L^{2}(\Omega))$ and therefore weakly converging to some limit $f\in L^{2}(0,T;L^{2}(\Omega))$. By calculations one can show that $\nabla\phi_{m}q_{m}$ converges strongly in $L^{1}(0,T;L^{1}(\Omega))$ and due to the uniqueness of the limit we can identify $f=\nabla\phi_{m}q_{m}$. Therefore the second term of (8.3) converges weakly to $A^{\prime}(\phi)\nabla\phi q$ in $L^{2}(0,T;L^{2}(\Omega))$.

All terms arising from $\mathbf{b}(\mathbf{u},\mathbf{v},\mathbf{w})$ or $\mathbf{B}(\mathbf{u},\mathbf{C},\mathbf{D})$ can be treated in a similar way as (8.2).

8.2 Bulk Stress

The convergence in the relaxation term

[TABLE]

can be shown analogously as in (8.1).

Next we consider the nonlinear diffusion term of the bulk stress equation

[TABLE]

The second integral is straightforward due to the weak convergence of \nabla\big{(}A(\phi_{m})q_{m}\big{)} in $L^{2}(0,T;L^{2}(\Omega))$. Convergence of the first integral term follows directly, if we know that \nabla\big{(}A(\phi_{m})\zeta\big{)} converges strongly in $L^{2}(0,T;L^{2}(\Omega))$. Due to the Sobolev embedding $\nabla\phi_{m}$ converges strongly in $L^{2}(0,T;L^{3}(\Omega))$. Calculating

[TABLE]

we see that both terms converge strongly in $L^{2}(0,T;L^{2}(\Omega))$ using the generalized dominated convergence theorem. Combining this with the weak convergence of \nabla\big{(}A(\phi_{m})q_{m}\big{)} in $L^{2}(0,T;L^{2}(\Omega))$ the first integral vanishes as $m\to\infty$.

Furthermore, we control the coupling term in the following way

[TABLE]

The first integral goes to zero due to the strong convergence of \nabla\big{(}A(\phi_{m})\zeta\big{)} in $L^{2}(0,T;L^{2}(\Omega))$. The second integral term vanishes due to the weak convergence of $n(\phi_{m})\nabla\mu_{m}$ in $L^{2}(0,T;L^{2}(\Omega))$.

8.3 Chemical Potential

The only nonlinear term in the chemical potential $\mu$, cf. (2.6), is

[TABLE]

Applying the mean value theorem, we have for $z\in[\phi_{m}(t),\phi(t)]$ or $z\in[\phi(t),\phi_{m}(t)],$ and $r,q,s>1,$ $\frac{1}{r}+\frac{1}{q}+\frac{1}{s}=1$ that (8.5)

[TABLE]

Taking into account (2.5) we know that

[TABLE]

therefor we can estimate (8.6) as follows

[TABLE]

This goes to zero due to the strong convergence of $\phi_{m}$ in $L^{2}(0,T;L^{s}(\Omega))$ for every $s\in(1,\infty)$ and the boundedness of $\phi_{m},\phi\in L^{\infty}(0,T;H^{1}(\Omega))$.

8.4 The Navier-Stokes part

We focus on the terms from the Navier-Stokes part $(\ref{eq:disc_weak})_{4}$. First, we see that

[TABLE]

goes to zero due to the weak convergence of $\nabla\mathbf{u}_{m}$ in $L^{2}(0,T;V)$ and the strong convergence of $\eta(\phi_{m})\nabla\mathbf{v}\varphi(t)$ in $L^{2}(0,T;L^{2}(\Omega))$, analogously as in (8.1). Therefore, we can also pass to the limit for the formulation with $\mathrm{D}\mathbf{u}:\mathrm{D}\mathbf{v}$.

We turn now our attention to the stress term in the Navier-Stokes part arising from the coupling to the Cahn-Hilliard part

[TABLE]

The first integral term goes to zero due to the weak convergence of $\nabla\mu_{m}\in L^{2}(0,T;L^{2}(\Omega))$, the strong convergence of $\phi_{m}\in L^{2}(0,T;H^{1}(\Omega))$ and the fact that $\mathbf{v}\varphi(t)\in L^{\infty}(0,T;H^{1}(\Omega))$. The second integral term vanishes due to the strong convergence of $\phi_{m}\in L^{2}(0,T;L^{p}(\Omega))$, for any $p\in(1,\infty)$.

The last term in the Navier-Stokes equation is the stress term arising from the viscoelastic contribution

[TABLE]

where we have used the strong convergence of $\mathbf{C}_{m},\mathrm{tr}\left({\mathbf{C}_{m}}\right)\in L^{2}(0,T;L^{4}(\Omega))$.

8.5 Conformation tensor

Finally, we consider the terms of the viscoelastic part $(\ref{eq:disc_weak})_{5}$. We present only one part of the upper convective term and the nonlinear term, all other terms are treated analogously. We can write

[TABLE]

and find that

[TABLE]

due to the strong convergence of $\mathbf{C}_{m},\mathrm{tr}\left({\mathbf{C}_{m}}\right)\in L^{2}(0,T;L^{4}(\Omega))$.

Further, we rewrite the nonlinear term in the following way

[TABLE]

This yields

[TABLE]

Now, due to the continuity of $h(\phi)$, it is clear that $h(\phi_{m})\mathrm{tr}\left({\mathbf{C}_{m}(t)}\right)^{2}\mathbf{C}_{m}(t)\rightarrow h(\phi)\mathrm{tr}\left({\mathbf{C}(t)}\right)^{2}\mathbf{C}(t)$. Consequently, we can pass to the limit in every term of the system and obtain the existence of a weak solution. This concludes the proof of existence.

8.6 Limiting process in the energy inequality

In order to obtain the differential form of the energy inequality (3.2) we need the following a priori estimates: $q^{\prime}\in L^{2}(0,T;H^{-1}(\Omega)),\mathbf{C}^{\prime}\in L^{2}(0,T;W^{*}(\Omega))$ and $\nabla\phi^{\prime}\in L^{2}(0,T;H^{-1}(\Omega))$. Since we have no control on these terms we can only deduce the integrated inequality (3.3). Let us denote

[TABLE]

Since $\nabla\phi_{m},q_{m},\mathbf{u}_{m},\mathbf{C}_{m}$ converge strongly in $L^{2}(0,T;L^{2}(\Omega))$ they also converge strongly for almost every $t\in(0,T)$ in $L^{2}(\Omega)$. Taking the difference between $E_{m}(t)$ and $E(t)$ we obtain two types of terms. The first one is of the form

[TABLE]

We use (8.11) with $g$ being $\nabla\phi,q,\mathbf{u}$ or $\mathbf{C}$. If $\left\|F^{\prime}(z)\right\|_{2}$ is bounded the second term is

[TABLE]

Here we have used the mean value theorem, where $z=\lambda\phi_{m}(t)+(1-\lambda)\phi(t)$ denotes a convex combination of $\phi_{m}(t)$ and $\phi(t)$, $\lambda\in(0,1)$. Applying analogous calculations as in (8.5), (8.7) now for $F^{\prime},|F^{\prime}(z)|=|f(z)|\leq c_{3}|z|^{p-1}+c_{4}$ we find $\left\|F^{\prime}(z)\right\|_{2}\leq c$.

This implies that we can pass to the limit in $E_{m}(t)$ strongly in $L^{1}(\Omega)$ for almost every $t\in(0,T)$.

The convergence of $E_{m}(0)$ follows from the fact that the initial conditions $\nabla\phi_{m0},q_{m0},\mathbf{u}_{m0},\mathbf{C}_{m0}$ converges strongly in the appropriate spaces. In detail $q_{m0},\mathbf{u}_{m0},\mathbf{C}_{m0}$ are $L^{2}$-projections with $H^{1}$ basis functions so their $L^{2}$-norms converge strongly. However, $\phi_{m0}$ converges strongly in $H^{1}(\Omega)$, due to the projection using the eigenfunctions of the Laplace operator. Therefore, $\nabla\phi_{0m}$ converges strongly to $\nabla\phi_{0}$, see [47]. For the convergence of $\int_{\Omega}F(\phi_{0m})$ to $\int_{\Omega}F(\phi_{0})$ we can employ an analogous calculation as in (8.12).

Integrating (6.1) in time from [math] to $t$ we obtain

[TABLE]

We also recall that due to the Fatou lemma for a weakly/weakly-⋆ convergent sequence $\{g_{m}\}$ we have

[TABLE]

Here we can choose $g$ as \sqrt{m(\phi)}\nabla\mu,\sqrt{\eta(\phi)}\mathrm{D}\mathbf{u},\nabla\mathbf{C},\big{(}\tau(\phi)\big{)}^{-1/2}q,\nabla\big{(}A(\phi)q\big{)} or $\sqrt{h(\phi)}\mathrm{tr}\left({\mathbf{C}}\right)^{2}$. Similarly as in (8.1) we can employ Lemma 2.2, it implies that these terms converge at least weakly in $L^{2}(0,T;L^{2}(\Omega))$, due to the fact that the square root of a continuous positive function is continuous and the $L^{\infty}$-bound is preserved. This allows us to pass to the limit in the dissipative terms. We continue with the term

[TABLE]

which will be treated in a different way. Indeed,

[TABLE]

The first integral term goes to zero due to (8.11) and the second integral term vanishes due to the generalized dominated convergence theorem. Thus, $h(\phi_{m})\left|\mathrm{tr}\left({\mathbf{C}}\right)\right|^{2}$ converges strongly in $L^{2}(0,T;L^{2}(\Omega))$.

Consequently, we can pass to the limit in each term and obtain (3.3), i.e.

[TABLE]

9 Numerical Simulations

In order to illustrate properties of our model for viscoelastic phase separation, we present results of some numerical experiments.

9.1 Lagrange-Galerkin finite element method

We start by describing a Lagrange-Galerkin finite element method for solving system (2.6) based on its weak formulation (3.1).

We employ piecewise-linear finite elements for $(\phi,q,\mu,\mathbf{u},p,\mathbf{C})$ in space and the characteristic scheme in time. Since the Cahn-Hilliard equation contains a fourth order operator and we have a saddle point problem in the Navier-Stokes part we have to use either stable spaces or apply suitable stabilization techniques, the so-called pressure stabilization, cf. [39].

Let $\mathcal{T}_{h}$ be a regular and quasi-uniform triangulation of $\bar{\Omega}$, with $h_{K}$ the diameter of $K\in\mathcal{T}_{h}$, $h:=\max(h_{K})$, $\rho_{K}:=\sup\{\textrm{diam}(S):S\text{ is a ball contained in }K\}$. Thus, we have [15]

[TABLE]

We define the discrete spaces

[TABLE]

Here we will shortly describe the characteristic scheme that has been employed in [39, 48] for similar fluid flow problems. Let $\Delta t$ be a time increment, $N_{T}=T/\Delta t$ the total number of time steps and $t^{n}=n\Delta T$. We consider the first order characteristic approximation of the material derivative

[TABLE]

where $X_{1}^{n+1}:\Omega\to\mathbb{R}^{2}$ is a mapping defined by

[TABLE]

We set $\phi^{n+1/2}:=(\phi^{n+1}+\phi^{n})/2,\,q^{n+1/2}:=(q^{n+1}+q^{n})/2$.

In order to solve numerically our viscoelastic phase separation model (2.6) we apply the operator splitting and solve in the first step $(\phi,q,\mu)$ system $(\ref{eq:full_model})_{1,2,5}$ and in the second step the Navier-Stokes-Peterlin equations $(\ref{eq:full_model})_{3,4}$.

**Step 1:

**For given $(\phi_{h}^{n},q_{h}^{n},\mu^{n}_{h},\mathbf{u}_{h}^{n},p_{h}^{n},\mathbf{C}_{h}^{n})\in(\Psi_{h})^{3}\times V_{h}\times Q_{h}\times W_{h}$ find $(\phi_{h}^{n+1},q_{h}^{n+1},\mu_{h}^{n+\frac{1}{2}})\in(\Psi_{h})^{3}$ such that:

$\displaystyle\left(\frac{\phi_{h}^{n+1}-\phi_{h}^{n}\circ X_{1}^{n+1}}{\Delta t},\psi_{h}\right)+\left(m(\phi_{h}^{n})\nabla\mu_{h}^{n+1/2},\nabla\psi_{h}\right)+\Delta t\left((\phi^{n}_{h})^{2}\nabla\mu_{h}^{n+1/2},\nabla\psi_{h}\right)$

\displaystyle\hskip 227.62204pt-\left(n(\phi_{h}^{n})\nabla\big{(}A(\phi_{h}^{n})q^{n+1/2}\big{)},\nabla\psi_{h}\right)=0

\displaystyle\left(\frac{q_{h}^{n+1}-q_{h}^{n}\circ X_{1}^{n+1}}{\Delta t},\psi_{h}\right)+\left(\nabla\big{(}A(\phi_{h}^{n})q_{h}^{n+1/2}\big{)},\nabla\big{(}A(\phi_{h}^{n})\psi_{h}\big{)}\right)

\displaystyle\hskip 128.0374pt+\left(\left(\tau(\phi_{h}^{n})\right)^{-1}q_{h}^{n+1/2},\psi_{h}\right)-\left(n(\phi_{h}^{n})\nabla\mu_{h}^{n+1/2},\nabla\big{(}A(\phi_{h}^{n})\psi_{h}\big{)}\right)=0

$\displaystyle\left(\mu_{h}^{n+1/2},\psi_{h}\right)-\left(\nabla\phi^{n+1/2},\nabla\psi_{h}\right)-\left(f(\phi_{h}^{n},\phi_{h}^{n+1}),\psi_{h}\right)=0$

Let us point out that step 1 is based on a mixed finite element method for $\phi$ and $\mu$, see [19], and the Taylor approximation $f(\phi_{h}^{n+1})\approx f(\phi_{h}^{n},\phi_{h}^{n+1}):=f(\phi_{h}^{n})+\frac{1}{2}\left(\phi_{h}^{n+1}-\phi_{h}^{n}\right)f^{\prime}(\phi_{h}^{n})$, cf. [50]. Moreover, we have applied the first order correction of the convective velocity in the Cahn-Hilliard equation by introducing $\mathbf{u}^{\star}=\mathbf{u}^{n}-\Delta t\phi_{h}^{n}\nabla\mu_{h}^{n+1/2}$ to stabilize the operator splitting, see [50].

**Step 2:

**For given $(\phi_{h}^{n+1},q_{h}^{n+1},\mu_{h}^{n+\frac{1}{2}},\mathbf{u}_{h}^{n},p_{h}^{n},\mathbf{C}_{h}^{n})\in(\Psi_{h})^{3}\times V_{h}\times Q_{h}\times W_{h}$ find $(\mathbf{u}_{h}^{n+1},p_{h}^{n+1},\mathbf{C}_{h}^{n+1})\in V_{h}\times Q_{h}\times W_{h}$ such that:

$\displaystyle\left(\frac{\mathbf{u}^{n+1}_{h}-\mathbf{u}^{n}_{h}\circ X_{1}^{n+1}}{\Delta t},\mathbf{v}_{h}\right)$ $\displaystyle+\left(\eta(\phi_{h}^{n+1/2})\mathrm{D}\mathbf{u}^{n+1}_{h},\mathrm{D}\mathbf{v}_{h}\right)-\left(\mathrm{div}\,{\mathbf{v}}_{h},p^{n+1}_{h}\right)$

$\displaystyle\hskip 56.9055pt-\left(\mathrm{div}\,{\mathbf{u}}^{n+1}_{h},\psi_{h}\right)-\delta_{0}\sum_{K\in\mathcal{T}_{h}}h_{K}^{2}\left(\nabla p_{h}^{n+1},\nabla\psi_{h}\right)_{K}$

$\displaystyle\hskip 56.9055pt+\left(\mathrm{tr}\left({\mathbf{C}^{n+1}_{h}}\right)\,\widetilde{\mathbf{C}^{n+1}_{h}},\nabla\mathbf{v}_{h}\right)+\left(\phi_{h}^{n}\nabla\mu_{h}^{n+1/2},\mathbf{v}_{h}\right)=0$

$\displaystyle\left(\frac{\mathbf{C}^{n+1}_{h}-\mathbf{C}^{n}_{h}\circ X_{1}^{n+1}}{\Delta t},\mathbf{D}_{h}\right)$ $\displaystyle+\varepsilon\left(\nabla\mathbf{C}^{n+1}_{h},\nabla\mathbf{D}_{h}\right)-2\left((\nabla\mathbf{u}^{n+1}_{h})\widetilde{\mathbf{C}^{n+1}_{h}},\mathbf{D}_{h}\right)+$

$\displaystyle+\left(h(\phi_{h}^{n+1/2})\widetilde{\mathrm{tr}\left({\mathbf{C}^{n+1}_{h}}\right)^{2}}\mathbf{C}^{n+1}_{h},\mathbf{D}_{h}\right)-\left(h(\phi_{h}^{n+1/2})\widetilde{\mathrm{tr}\left({\mathbf{C}^{n+1}_{h}}\right)}\,\mathbf{I},\mathbf{D}_{h}\right)=0$

Note that in step 2 we have applied the Brezzi-Pitkaränta pressure stabilization [10] for the Navier-Stokes equation. Further, we have an inner fixed-point iteration to iterate nonlinear terms with respect to $\mathbf{C}$. We denote by $\widetilde{g}$ the function which will be replaced by the old inner iteration and we iterate equations in step 2 until the relative norm differences of the iterations is less than a given tolerance.

It has been proven in [50] that the above splitting scheme indeed satisfies the discrete dissipation at each time step for the classical Oldroyd-B model instead of $(\ref{eq:full_model})_{4}$. In our case it would mean

[TABLE]

as soon as

[TABLE]

A rigorous proof will be an interesting question for future study. Nevertheless, our numerical experiments indicate that the discrete energy dissipation property (9.2) indeed holds also for (2.6) with the Peterlin model.

9.2 Numerical Experiments

The aim of this section is to illustrate complex behaviour of our viscoelastic phase separation model (2.6) in the context of spinodal decomposition. Instead of the classical Ginzburg-Landau potential (2.2) we work with a modified Ginzburg-Landau potential (9.2) for the following reasons. In the spinodal decomposition the potential usually has two minima which corresponds to the equilibrium volume fractions of the two phases. Physically relevant Flory-Huggins potential (2.3) attains both minima $a,b$ inside $[0,1]$ and not at the boundary of $[0,1]$. Therefore, we have opted for the modification

[TABLE]

with $a=0.134791$ and $b=1-a$. These values are consistent with Flory-Huggins potential (2.3) for $n_{p}=n_{s}=1\text{ and }\chi=28/11$. In general, for $a=1,b=0$ it may happen that $\phi$ attains values outside of $[0,1]$ due to the roundoff errors.

We triangulate a computational domain $\Omega=[0,128]^{2}$ into $128^{2}$ isometric triangles. The following functions and parameters will be used

[TABLE]

Experiment 1: In this test we set the initial data to $\phi_{0}(x)=0.4+\xi(x),q_{0}=0,\mathbf{u}_{0}=\mathbf{0},\mathbf{C}_{0}=\sqrt{2}\mathbf{I}$, where $\xi(x)$ is a random perturbation from $[-10^{-3},10^{-3}]$.

Figures 1-3 present the time evolution of numerical solutions for $\phi,q,p,\mu$ and $\left|\mathbf{u}\right|_{2}$. We can clearly recognize the frozen phase $(t\approx 600)$, the elastic regime with solvent-rich droplets $(t\approx 1200)$, volume shrinking $(t\approx 1800)$ and network pattern $(t\approx 2400)$ and later the relaxation regime $(t\approx 4000)$.

For longer simulations we will also recover the hydrodynamic regime with a typical droplet pattern [50]. These results are qualitatively in a good agreement with physical experiments presented in [51], see Figure 8, and with [55].

In Figure 4 we have plotted the time evolution of partial contributions of the total free energy. We can also realize that the numerical scheme is practically energy dissipative and mass conservative, see Figure 8 and 11.

In Figure 2 we observe that the variables $\mu$ and $q$ have a quite similar appearance despite of the difference in magnitude, which results from the fact that the evolution of $q$ is induced by $\mu$. Further, the pressure $p$ also exhibits a network-like structure, although to velocity looks quite unstructured. In all experiments the conformation tensors are almost constant and we omit their plots.

Experiment 2: The next simulation differs only in the potential and the initial condition for the conformation tensor. More precisely, we have now used the Flory-Huggins potential with $n_{p}=n_{s}=1\text{ and }\chi=28/11$ and $\mathbf{C}_{0}=\mathbf{I}$. We observe similar phase evolution as in Experiment 1 but on a different timescale. We can notice that simulations with the Flory Huggins potential are roughly twice faster than with the Ginzburg-Landau potential.

Experiment 3: In this test we want to study the influence of non-zero initial velocity on the dynamics of the system with the Ginzburg-Landau potential. We denote by $\mathds{1}_{B}$ the characteristic function of a ball with the radius $50$ and midpoint $(64,64)$. The initial conditions are

[TABLE]

All other initial data are the same as in the previous experiments.

In Figure 9 we can recognize typical phase separation patterns as in the previous experiments. However, their evolution is much faster than in Experiment 1 which is due to the initial velocity field. In fact the phases occur now at similar times as in Experiment 2. This effect can be related to the so-called shear-induced phase separation. In Figure 10 we can observe the slow down of the vector field up to $t\approx 750$ and then the dispersion of the initial rotation. At this point the phase separation effects start to dominate the rotation and therefore we observe similar patterns as in Experiment 1 and 2.

10 Conclusion and Outlook

In this paper we have analyzed mathematical properties of the viscoelastic phase separation model (2.6) that has been proposed in [55]. Instead of the Oldroyd-B equations for the time evolution of viscoelastic stress tensor we have used a more general Peterlin model as proposed in [39].

A regular potential for the phase separation, i.e. the Ginzburg-Landau potential, allowed us to prove the existence of global weak solutions to (2.6). To this end we had to modify slightly the original model of [55] by introducing a parameter $\kappa\in[0,1)$ in front of the viscous coupling terms in the equations for volume fraction $\phi$ and bulk stress $q$, see $(\ref{eq:full_model})_{1},(\ref{eq:full_model})_{2}$, respectively. We have derived the necessary a priori estimates and passed to the limit in the Galerkin approximations. Applicability of the model (2.6) to describe complex dynamics of spinodal decomposition has been documented by numerical experiments presented in Section 9, see, also [50]. Our numerical experiments demonstrate that the choice $\kappa=1$ makes no problem for the existence of a weak solution in practice. Therefore, it will be interesting to investigate in the future whether analogous existence results hold also for $\kappa\rightarrow 1$.

Building on the present study our next goal is to investigate a degenerate model with the Flory-Huggins potential and degenerate mobilities, see [12].

Acknowledgement

This research was supported by the German Science Foundation (DFG) under the Collaborative Research Center TRR 146 Multiscale Simulation Methods for Soft Matters (Project C3) and partially by the VEGA grant 1/0684/17 and the Mainz Institute of Multiscale Modeling (M3odel). We would like to thank B. Dünweg, D. Spiller and J. Kat’uchová for fruitful discussions on the topic.

11 Appendix

11.1 Model derivation

The aim of this section is to discuss a derivation of the viscoelastic phase separation model (2.6). Systematic model derivation and further discussions on various model variants can be found in our recent works [3, 11]. A complete rigorous derivation from a microscopic description is however not yet available in the literature.

To keep our paper self-consistent we present important steps of the model derivation. First, as it is standardly used, cf. [28, 29, 55], the derivation of a thermodynamically consistent model is divided into the description of reversible and irreversible dynamics. The reversible part can be obtained by the virtual work principle [20], while the irreversible part will be treated in a dissipative setting as in [30]. Our further discussion is structured in the following way.

Discussion on the free energy and introduction of a conservative model. 2. 2.

Consideration of non-equilibrium dynamics using the relative velocity $\mathbf{w}$. 3. 3.

Consideration of dissipative effects.

11.2 Free energy and a conservative model

First, we review the energy of the original viscoelastic phase separation model from [55] rewritten using the conformation tensor $\mathbf{C}.$ The total energy reads

[TABLE]

where $F(\phi)$ denotes an appropriate mixing potential such as (2.2) or (2.3). The first integral represents the mixing energy, while the second stands for the total kinetic energy. We note that the last integral represents the free energy modeled by the viscoelastic conformation tensor. Following [51] the second last term models the energy contribution due to the bulk stress $q\mathbf{I}.$ The latter originates from (microscopic) intermolecular attractive interactions and can be seen as additional viscoelastic pressure. We refer a reader to our recent work [11] for a more detailed and physically motivated discussion on the connection between the conformation stress tensor $\mathbf{C}$ and $q\mathbf{I}$. Indeed, the conformation tensor $\mathbf{C}$ is connected to the viscoelastic effects of the volume-averaged velocity $\mathbf{u}$ and $q\mathbf{I}$ is connected to the relative velocity $\mathbf{w}=\mathbf{u}_{p}-\mathbf{u}_{s}.$ Here $\mathbf{u}_{p},\mathbf{u}_{s}$ are the velocities of polymer and solvent phases, respectively.

We proceed with formulating the total energy of our desired model by

[TABLE]

We note that the quadratic term in the last integral results from a nonlinear dependence of the viscoelastic stress tensor on the conformation tensor and the term $-\mathrm{tr}({\ln(\mathbf{C})})$ represents the free energy of non-interacting macro-molecules, see [54] and [2].
This choice of the elastic energy yields a phase separation model, where viscoelastic effects are described by the Peterlin model. The latter was extensively studied in our recent works [38, 39, 44] both from the analytical and numerical point of view. We note that this viscoelastic model has been also derived as a rigorous closure of a kinetic model in [32]. The derivation of a full viscoelastic phase separation model presented below is general and applicable also to other choices of the elastic energy, see, e.g., [41] for the FENE-P model or [40] for a generalized Peterlin model. Using such generalizations within a full phase separation model is an interesting topic for future research.

The conservative dynamics of viscoelastic phase separation can be derived applying the Poisson bracket formalism, as presented in [3]. In general, one should derive a compressible model via the standard Poisson brackets and afterwards apply a reduction to an incompressible model, similar to a derivation presented in [3]. Such an approach yields the following basic equations, see [11],

[TABLE]

It can be shown that formally the model conserves the free energy given by (11.1). Clearly, model (11.2) only describes the reversible part of the dynamics.

11.3 Closure relation for the relative velocity

In what follows we introduce the relative velocity $\mathbf{w}$. The conservative model (11.2) can be understand as a reduced model with $\mathbf{w}=0$. Following our derivation in [11] $q\mathbf{I}$ is associated with the divergence of $\mathbf{w}$, while the viscoelastic stress tensor $\mathbf{T}$ arising from $\mathbf{C}$, i.e. $\mathbf{T}=\mathrm{tr}(\mathbf{C})\mathbf{C}-\mathbf{I}$, is associated to the deformation gradient of the velocity $\mathbf{u}$. Finally, due to the divergence-freedom of the velocity $\mathbf{u}$ there is no relevant stress associated with $\mathrm{div}\left({\mathbf{u}}\right)$. This will be include via the linear term $A_{1}(\phi)\mathrm{div}\left({\mathbf{w}}\right)$ into $q$. A short calculation shows that for the upper convected derivative $\frac{\mathrm{\overset{\triangledown}{D}}\mathbf{I}}{\mathrm{Dt}}=2\mathrm{D}\mathbf{u}$, hence an evolution equation for the viscoelastic stress tensor $\mathbf{T}$ already includes a dependence on the deformation gradient.

Using the relative velocity we obtain the following equations from (11.2)

[TABLE]

We already know that the model is conservative for $\mathbf{w}=0$. In order to close the system we will choose a suitable closure relation for $\mathbf{w}$, which is based on the second law of thermodynamics. Computing the change of the free energy (11.1) yields

[TABLE]

Hence, the above expression is non-negative for $\mathbf{w}=-[\phi(1-\phi)\nabla\frac{\delta E}{\delta\phi}-\nabla(A(\phi)q)]$. We again refer to [11] where we have showed that such a term arises even in more general processes based on the relative velocity $\mathbf{w}$. Note that the model reduction by choosing a closure relation for $\mathbf{w}$ introduces dissipation into the system.

11.4 Dissipative part

In this part, we will introduce additional physically motivated dissipation terms. First, we add a diffusive term to the equations of the conformation tensor $\mathbf{C}$. As mentioned in the introduction the existence of such a diffusive term in thermodynamically consistent viscoelastic models has been investigated and verified in [1, 45]. Furthermore, we add Navier-Stokes-type viscosity terms to the momentum equations.

Since both $q\mathbf{I},\mathbf{C}$ are modeling the viscoelastic effects, we expect them to decay to an equilibrium in the absence of all forces, which will be modeled by the relaxation terms. The generic form for such a relaxation is given by

[TABLE]

We note in passing that in order to preserve frame invariant structure, the choice of $g(\mathbf{C})$ is restricted. In our work we have opted for the Peterlin model, i.e. we choose $g(\mathbf{C})=\mathrm{tr}(\mathbf{C})\mathbf{C}$. For the bulk stress we choose the simplest available mechanism, i.e. $g(q)=1$. Altogether we obtain the following resulting system

[TABLE]

11.5 Reformulation and remarks

In this subsection, we will suitably reformulate (11.5) to obtain our model (2.6).

First, the Leslie-Erickson stress tensor $-c_{0}\mathrm{div}\left({\nabla\phi\otimes\nabla\phi}\right)$ can be rewritten according to [50]

[TABLE]

where the second term can be absorbed into the pressure $p$. Reordering and expanding terms with the choices

[TABLE]

yields the desired system governing viscoelastic phase separation.

\displaystyle\begin{split}\frac{\partial\phi}{\partial t}+\mathbf{u}\cdot\nabla\phi&=\mathrm{div}\left({m(\phi)\nabla\mu}\right)-\mathrm{div}\left({n(\phi)\nabla\big{(}A(\phi)q\big{)}}\right)\\ \frac{\partial q}{\partial t}+\mathbf{u}\cdot\nabla q&=-\frac{1}{\tau_{b}(\phi)}q+A(\phi)\mathrm{div}\left({\nabla\big{(}A(\phi)q\big{)}}\right)-A(\phi)\mathrm{div}\left({n(\phi)\nabla\mu}\right)\\ \frac{\partial\mathbf{u}}{\partial t}+(\mathbf{u}\cdot\nabla)\mathbf{u}&=\mathrm{div}\left({\eta(\phi)\Big{(}\nabla\mathbf{u}+(\nabla\mathbf{u})^{T}\Big{)}}\right)-\nabla p+\mathrm{div}\,{(}\mathbf{T})+\mu\nabla\phi\\ \frac{\partial\mathbf{C}}{\partial t}+(\mathbf{u}\cdot\nabla)\mathbf{C}&=(\nabla\mathbf{u})\mathbf{C}+\mathbf{C}(\nabla\mathbf{u})^{T}+h(\phi)\mathrm{tr}(\mathbf{C})\mathbf{I}-h(\phi)\mathrm{tr}(\mathbf{C})^{2}\mathbf{C}+\varepsilon\Delta\mathbf{C}\\ \mathrm{div}\,{(}\mathbf{u})&=0\hskip 28.45274pt\mathbf{T}=\mathrm{tr}\left({\mathbf{C}}\right)\mathbf{C}\hskip 28.45274pt\mu=-c_{0}\Delta\phi+f(\phi)\end{split}

(11.6)

A formal proof of thermodynamic consistency of (11.6) is given in Theorem 3.2 with $\kappa=1$.

Remark 11.1**.**

We note that the above derivation of model (11.6) is consistent with the GENERIC formalism [29, 28, 11, 3]. GENERIC (general equation for the non-equilibrium reversible–irreversible coupling) is a well-accepted framework for the derivation of thermodynamically admissible evolution equations for non-equilibrium systems. Analogously, one can consider the approach proposed recently by Málek et al. [45, 46] and derive a thermodynamically consistent model by prescribing the state variables, here $(\phi,q,\mathbf{u},\mathbf{C})$, a free energy and a dissipation process modeling the energy dissipation. The free energy is given by

[TABLE]

and dissipation process by

[TABLE]

As pointed out above the conservative part of the model, i.e. (11.2), can be derived by means of the Poisson bracket formalism using only the state variables and the free energy, i.e. the dissipation process is zero. In the second step, the dissipative contributions are incorporated using a non-zero dissipation process.

As one can quickly recognize (11.6) differs only slightly from the studied model (2.6). The energy dissipation (11.4) yields the estimate

[TABLE]

However, we do not have independent a priori estimates for $\nabla\mu$ and $\nabla(A(\phi)q)$. This is due to the cross-diffusion structure in the $(\phi,q)$ subsystem. A typical technique to deal with such problems is to dominate the cross-diffusion by a diagonal diffusion. To cure this problem we have introduced in (11.6) a parameter $\kappa\in[0,1)$, i.e. $n(\phi):=\kappa n(\phi)$ which leads to our final viscoelastic phase separation model (2.6).