A Cahn–Hilliard phase field model coupled to an Allen–Cahn model of viscoelasticity at large strains

The study of Cahn–Hilliard phase field models coupled with finite viscoelasticity has gained increasing interest in the recent literature, cf., e.g., [6, 12, 2, 24]. These models describe the phase separation phenomena for multiphase materials in presence of elastic interactions between the materials constituents, and may find applications e.g. in soft matter dynamics [6], tumor growth dynamics [12], neurological and neuromuscular deseases (see the discussion in [2]).

In these works both the phase field and the viscoelasticity governing equations are formulated in the Eulerian reference configuration. Hence, the main state variable for elasticity is the velocity field, and the deformation gradient, entering in the elasticity constitutive assumptions, is determined solving a transport equation in terms of the velocity and the velocity gradient. A similar modeling approach for evolutionary models with finite viscoelasticity was implemented also in [3], which studies the Oldroyd-B model for a dilute polymeric fluid, in [7], which studies the motion of a class of incompressible heat-conducting viscoelastic rate-type fluids with stress-diffusion, in [4, 11], which study problems in magnetoelasticity, in [18], which studies an evolutionary non-isothermal viscoelastic problem, and in [20], which studies the diffusion of a solvent in a saturated hyperelastic porous solid of viscoelastic Kelvin-Voigt type at large strains.

A major challenge in studying the existence of weak solutions for phase field models coupled with finite viscoelasticity is to obtain results in three space dimensions and for generic elastic energy densities which are highly nonlinear in the deformation gradient and which may depend on the phase field variable. Indeed, in [4, 6, 11, 12], employing a second order diffusive regularization in the transport equation for the deformation gradient (i.e. adding a term proportional to the Laplacian of the deformation gradient), the existence of global weak solutions for the model is proved in two space dimensions and for elastic energy densities of Neo–Hookean type, i.e. quadratic in the deformation gradient, and independent of the phase field variable. The diffusive regularization of the transport equation, which is needed to increase the regularity of the deformation gradient and gain compactness of its approximations, breaks the kinematic relationship between the velocity and the deformation gradient variables. In [2], we designed a specific second order diffusive regularization of the transport equation for the deformation gradient, depending on both the phase field variable and on the deformation gradient, which allowed us to prove the existence of global weak solutions for elastic energy densities of Neo–Hookean type which depend on the phase field variable. The degeneracy of the elastic energy density depending on the phase field variable was not allowed in the framework of [2]. Also, existence results were obtained in three space dimensions by adding a further viscous regularization to the Cahn–Hilliard equation. By employing a different regularization approach of the model equations, obtained by adding a dissipative contribution to the Cauchy stress tensor which involves high order nonlinear terms in the small strain rate (i.e. the symmetric part of the velocity gradient), in [18] the author obtained existence and regularity results of a distributional solution for an evolutionary non-isothermal viscoelastic problem with nonlinear and compressible elastic energy. The latter regularization approach preserves the kinematic relationship between the velocity and the deformation gradient variables. Its main drawback is the introduction of high order nonlinear terms involving the velocity gradient in the definition of the Cauchy stress tensor, making numerical approximations of the model not straightforward.

In the present paper, we introduce a new kind of diffusive regularization, of Allen–Cahn type, in the transport equation for the deformation gradient, complemented by the introduction of a regularizing interface term depending on the gradient of the deformation gradient in the free energy density of the system. The designed regularization preserves the dissipative structure of the equations. The idea for the introduction of this kind of regularization comes from an observation by Roubicek [19, Remark 3], which suggested the possibility to consider a Cahn–Hilliard regularization of the transport equation for the deformation gradient in order to obtain stronger existence results with respect to those available in the literature and based on diffusive regularizations of the transport equation for the deformation gradient. The introduction of an interface contribution for the deformation gradient in the free energy enhances the space and time regularity of the deformation gradient, while increasing the degree of nonlinearity of the coupled system. In particular, the Cauchy stress tensor contains second order and quadratic first order terms in the deformation gradient. In our framework, the regularity of the latter terms is obtained by adding an Allen–Cahn type second order diffusive regularization in the transport equation for the deformation gradient. Following the Allen–Cahn literature, we consider a dual mixed Allen–Cahn formulation of the transport equation for the deformation gradient, introducing a dual variable of the deformation gradient which enters also in the expression for the Cauchy stress tensor. This theoretical framework allows us to obtain the existence of global weak solutions in three space dimensions and for generic nonlinear elastic energy densities with polynomial growth, which are coupled to the phase field variable and which possibly degenerate for some values of the phase field variable.

The resulting PDE system for the considered model is the following:

valid in $\Omega_{T}$, endowed with the boundary conditions

on $\partial\Omega\times[0,T]$, and with proper initial conditions. Here, $\mathbf{v}$ is the velocity field, $s$ is the pressure, $\mathbf{F}$ is the elastic deformation gradient, $\mathbf{M}$ is its dual variable, $\phi$ is the phase field variable and $\mu$ is the chemical potential. Moreover, $\nu$ is a physical parameter, representing the viscosity of the material, while $\gamma$ and $\lambda$ are positive regularization parameters. The function $b(\phi)$ represents a positive phase-dependent mobility, $\psi(\phi)$ represents the bulk potential of the Cahn–Hilliard equation, and $w(\phi,\mathbf{F})$ is the elastic free energy density.

By means of an iterative argument based on elliptic regularity bootstrap steps applied to (1)₄, we find the maximum allowed polynomial growths of the Cahn–Hilliard potential and the elastic energy density which guarantee existence of a solutions in three space dimensions.

We observe that the dual mixed formulation of System (1) is particularly suitable to design standard and efficient numerical approximations. Hence, we also propose two kinds of unconditionally energy stable finite element approximations of the model, based on convex splitting ideas and on the use of a scalar auxiliary variable, proving the existence and stability of discrete solutions. We finally show numerical tests for different test cases with shape memory alloy type free energy, proving the interplay between phase separation and finite elasticity in determining the topology of stationary states with pure phases characterized by different elastic properties.

Hence, the novelties of the present work with respect to previous studies in literature are the following:

• •

  Proof of existence of global weak solutions for a Cahn–Hilliard phase field model coupled to viscoelasticity at large strains in three space dimensions, based on a new Allen–Cahn type regularization of finite viscoelasticity, for generic finite elastic energy densities with polynomial growth, which are coupled to the phase field variable and which possibly degenerate for some values of the phase field variable;

• •

  Design and implementation of standard and efficient well posed and unconditionally energy stable finite element approximations of the model.

The paper is organized as follows. In Section $2$ we introduce some notation regarding the employed tensor calculus and functional analysis. In Section $3$ we develop the model derivation. In Section $4$ we report the study of existence for a global weak solution of System (1). In Section $5$ we report the design of finite element approximations of the model and the study of existence and stability of discrete solutions. In Section $6$ we show results of numerical simulations. Finally, Section $7$ contains concluding remarks and future perspectives.

Let $\Omega\subset\mathbb{R}^{3}$ be an open bounded domain in $\mathbb{R}^{3}$. Let $T>0$ denote some final time, and set $\Omega_{T}:=\Omega\times(0,T)$. We start by introducing the notation for vectorial and tensorial calculus. Given $\mathbf{a},\mathbf{b}\in\mathbb{R}^{3}$, we denote by $\mathbf{a}\cdot\mathbf{b}\in\mathbb{R}$ their canonical scalar product in $\mathbb{R}^{3}$, with associated norm $|\mathbf{a}|:=(\mathbf{a}\cdot\mathbf{a})^{\frac{1}{2}}$, and by $\mathbf{a}\otimes\mathbf{b}\in\mathbb{R}^{3\times 3}$ their tensorial product. We indicate by $\{\mathbf{e}_{i}\}_{i=1,2,3}$ the canonical basis of $\mathbb{R}^{3}$. Given two second order tensors $\mathbf{A},\mathbf{B}\in\mathbb{R}^{3\times 3}$, we denote by $\mathbf{A}\colon\mathbf{B}\in\mathbb{R}$ their Frobenius scalar product in $\mathbb{R}^{3\times 3}$, i.e. by components $\mathbf{A}\colon\mathbf{B}:=\sum_{i,j=1}^{3}A_{ij}B_{ij}$, with associated norm $|\mathbf{A}|:=(\mathbf{A}\colon\mathbf{A})^{\frac{1}{2}}$. We indicate by $\{\mathbf{e}_{ij}:=\mathbf{e}_{i}\otimes\mathbf{e}_{j}\}_{i,j=1,2,3}$ the canonical basis of $\mathbb{R}^{3\times 3}$. Given a second order tensor $\mathbf{A}\in\mathbb{R}^{3\times 3}$, we indicate with the notation $\mathbf{A}^{[i]}\in\mathbb{R}^{3}$ the $i$-th column vector of $\mathbf{A}$. Given two third order tensors $\mathbf{C},\mathbf{D}\in\mathbb{R}^{3\times 3\times 3}$, we denote by $\mathbf{C}\mathbin{\setstackgap{S}{0.4ex}\mathrel{\Shortstack{{.}{.}{.}}}}\mathbf{D}\in\mathbb{R}$ their scalar product in $\mathbb{R}^{3\times 3\times 3}$, i.e. by components

We also introduce the operation, defined by components,

which contracts respectively a third order tensor $\mathbf{C}\in\mathbb{R}^{3\times 3\times 3}$ and a second order tensor $\mathbf{B}\in\mathbb{R}^{3\times 3}$ to a vector $\mathbf{C}\odot\mathbf{B}\in\mathbb{R}^{3}$ and two third order tensors $\mathbf{C},\mathbf{D}\in\mathbb{R}^{3\times 3\times 3}$ to a second order tensor $\mathbf{C}\odot\mathbf{D}\in\mathbb{R}^{3\times 3}$.

The phase field model coupled with finite viscoelasticity which we analyze in this paper is a particular case of a class of phase field models derived in our previous contribution [2], obtained by considering a binary, saturated, closed and non reactive mixture of a solid elastic component and a liquid component, whose dynamics is driven by the microscopic interactions between its constituents and by their macroscopic visco-elastic behaviour. The mass and momentum balance of the mixture were derived using a generalized form of the principles of virtual powers, giving constitutive assumptions satisfying the first and second law of thermodynamics in isothermal situations.

In particular, we consider the following form for the free energy $E$ of the system in Eulerian coordinates, associated to an arbitrary subregion of the mixture $R(t)\subset\Omega$ moving with the mixture:

where $\phi$ is the volume concentration of the elastic phase, $\mathbf{F}$ is the deformation gradient associated to the motion of the elastic phase, $w(\phi,\mathbf{F})$ is the hyperelastic free energy density and $\psi(\phi)$ represents a bulk energy due to the mechanical interactions of the micro–components. A term proportional to $\frac{1}{2}|\nabla\phi|^{2}+\psi(\phi)$ represents a surface energy contribution of the interface between the phases, expressed through a diffuse interface approach, while the term proportional to $|\nabla\mathbf{F}|^{2}$ is associated to elastic energy contributions from interfaces in the elastic material. Here, $\lambda>0$ is a regularization parameter, while the interface thickness parameter associated to the surface energy between the phases is taken to be equal to $1$ for ease of notation. We also consider the incompressibility constraint for the partial volume of the elastic phase. The particular model is obtained from [2, System (37)], with the regularization parameters $\theta=\delta=0$ and employing the following simplifying assumptions:

• •

  The momentum transfer in the mixture due to shear stresses in the liquid is negligible with respect to the momentum transfer between the solid and the liquid components;

• •

  The momentum exchange between the phases is induced by a Darcy-like flow of the liquid phase through the porous-permeable solid matrix associated to the solid phase, with constant permeability.

With these assumptions, the PDE system takes the form

valid in $\Omega_{T}$, endowed with the boundary conditions

on $\partial\Omega\times[0,T]$, and with proper initial conditions. Here, $\nu>0$ is a viscosity parameter, $s$ is the solid pressure and $b(\phi)$ is a positive mobility function. System (4) is analogous to [2, System (39)], with (3) as the free energy of the system.

Using the relation

renaming $s\leftarrow s+e$ and adding a proper diffusive regularization, with regularization parameter $\gamma>0$, in the transport equation (4)₃, we get

We introduce the auxiliary variable

and rewrite system (8) as (cf. also System (1) in the Introduction)

valid in $\Omega_{T}$, endowed with the boundary conditions

on $\partial\Omega\times[0,T]$, and with initial conditions $\mathbf{F}(\mathbf{x},0)=\mathbf{F}_{0}(\mathbf{x})$, $\phi(\mathbf{x},0)=\phi_{0}(\mathbf{x})$ for $\mathbf{x}\in\Omega$.

In order to proceed, we make the following assumptions:

• A0

  $\Omega\subset\mathbb{R}^{3}$ is a bounded domain and the boundary $\partial\Omega$ is of class $C^{2}$;

• A1

  $b\in C^{0}(\mathbb{R})$ and there exist $b_{0},b_{1}>0$ such that $b_{0}\leq b(r)\leq b_{1},$ $\forall r\in\mathbb{R}$;

• A2

  $w\in C^{1}(\mathbb{R}\times\mathbb{R}^{3\times 3};\mathbb{R})$, and there exist $d_{1}\geq 0,d_{2}>0$ such that $-d_{1}\leq w(r,\mathbf{T})\leq d_{2}\left(1+|\mathbf{T}|^{p}\right)$ for all $r\in\mathbb{R},\mathbf{T}\in\mathbb{R}^{3\times 3}$,with $p\in[0,6)$. Moreover, there exist $d_{3},d_{4}>0$ such that $\left|\partial_{\mathbf{T}}w(r,\mathbf{T})\right|\leq d_{3}\left(1+|\mathbf{T}|^{p-1}\right)$ and $\left|\partial_{r}w(r,\mathbf{T})\right|\leq d_{4}\left(1+|\mathbf{T}|^{p}\right)$ for all $r\in\mathbb{R},\mathbf{T}\in\mathbb{R}^{3\times 3}$;

• A3

  $\psi\in C^{1}(\mathbb{R})$ and there exist $c_{1}>0,c_{2}\geq 0$ such that $|\psi^{\prime}(r)|\leq c_{1}\left(|r|^{l}+1\right)$, $\psi(r)\geq-c_{2}$, $\forall r\in\mathbb{R}$, with $l\in\left[0,6+\frac{8}{\max(2,2p-6)}\right)$. Moreover, there exists a convex decomposition of $\psi=\psi_{+}+\psi_{-}$, where $\psi_{+}$ is convex and $\psi_{-}$ is concave, such that $|\psi_{-}^{\prime\prime}(r)|\leq c_{1}\left(|r|^{q}+1\right)$, $\forall r\in\mathbb{R}$, with $q\in[0,4)$;

• A4

  The initial data have the regularity $\mathbf{F}_{0}\in H^{1}(\Omega;\mathbb{R}^{3\times 3})$, $\phi_{0}\in H^{1}(\Omega)$.

We state here the main theorem of the present work concerning the existence of a global weak solution to (9) in three space dimensions, which will be proved in the forthcoming sections.

In the following, we will introduce a proper truncation of the growth behavior of $w(\phi,\mathbf{F})$, depending on a truncation parameter $R>1$, and we will define a proper Faedo–Galerkin approximation of a truncated version of (9), proving the existence of a discrete solution and studying its convergence to a continuous weak solution, as the discretization parameter tends to zero, in three space dimensions. We will then study the limit problem, at the continuous level, as the truncation parameter $R\to\infty$, thus removing the truncation operation and recovering an existence result for system (9) associated to the full growth laws assumed in Assumption A2.

Before proceeding, we state some preliminary results which will be used in the analysis.

We first need to obtain an existence result for a truncated system.