Control of nonlinear bulk deformation and large shear strain on first-order phase transformation kinetics

Abstract

Phase transformations play a key role in numerous coupled natural processes, and they are important for many industrial applications. However, the kinetics of phase transformations in coupled chemo-mechanical systems undergoing large mechanical deformations still needs to be better quantified. Here, we study the phase transformation kinetics of a two-phase binary mixture using the diffuse interface approach. We couple a Cahn–Hilliard type model with a mechanical model for a compressible viscous flow. The bulk compressibility is a nonlinear function of the pressure, and the shear viscosity is a nonlinear function of the concentration. The mechanical coupling is achieved by employing a pressure-dependent mechanical mixing term in the equation for the Gibbs energy. We derive a dimensionless system of equations which we solve numerically with a pseudo-transient method using conservative finite differences for discretization. We perform numerical simulations in 1D and 2D model setups considering far-field simple shear and pure shear. For a chemo-mechanically coupled system, we show that the velocity of the phase boundary is a linear function of the degree of metastability and, hence, confirm the hypothesis of “normal growth.” A stronger mechanical coupling and a larger volumetric effect of the chemical reaction result in lower phase boundary velocities. The 2D results show a significant impact of the mechanical coupling and the far-field deformation on the orientation and kinetics of the phase transformations. Under far-field simple shear and pure shear in 2D, the phase transformations generate string-like patterns. The orientation of these patterns is controlled by the applied far-field deformation and orientations differ by 45 degrees between simple shear and pure shear.

1 Introduction

Phase transformations play a key role in many natural phenomena, for example, geodynamic processes involving magmatism and volcanism, and they are important for many practical applications, such as in metallurgy or ceramic engineering [1,2,3,4]. A vast number of studies have investigated different aspects and types of phase transformations, and there exists a thorough theoretical framework for the study of phase transformations [3,4,5,6,7,8,9,10,11]. Phase transformation occurs at the boundaries between phases. The kinetics of such phase transformations, i.e., the speed of the phase front propagation, control the way systems behave while out of equilibrium [1, 12]. Understanding kinetics is therefore important, especially when the phase transformations occur simultaneously with other processes, e.g., mechanical deformation. Both experimental and theoretical studies have emphasized the importance of the coupling of chemical reactions and phase transformations with mechanical deformation and the associated pressure field [13,14,15,16,17].

There are two main approaches to modeling the moving phase fronts: sharp interface and diffuse interface. In the sharp interface approach, the phase boundaries are represented as 1D or 2D surfaces in a 2D and 3D continuum, respectively, and are explicitly tracked. To describe the motion of the sharp interface, and, therefore, the kinetics of the phase transformation, it is necessary to invoke the “jump” conditions on the phase front.

The simplest kinetic model of “normal growth” assumes a linear dependence of the velocity of the phase boundary on the overstepping of equilibrium conditions, i.e., the degree of metastability. The phenomenological material parameters required to constrain the front propagation velocity can be constrained using data from laboratory experiments.

On the other hand, the diffuse interface approach involves modeling both the phases and interfaces using either physics-based [18, 19] or ad hoc phase field models [20]. In a diffuse interface approach, the phase boundary is a represented by a layer with a finite thickness.

In this work, we use a Cahn–Hilliard-type diffuse interface model [21,22,23,24]. Such Cahn–Hilliard-type model, which resolves the fine structure of an advancing phase boundary [25], “can derive, rather than postulate, a kinetic relation governing the mobility of the phase boundary and check the validity of the “normal growth” approximation” [23].

There are three characteristic regimes of phase boundary migration occurring during phase transformations: mechanics controlled, diffusion controlled and interface controlled [23, 26]. In the mechanics-controlled regime, the volume change associated with the phase boundary migration is limited by the viscosity, or more generally the mechanical resistance, of the system [26]. Here, we focus on the relation between the diffusion- and interface-controlled regimes. For the diffusion-controlled regime, the growth rate of a stable phase is limited by the supply of a component which is governed by a diffusion type equation [23]. The diffusive transport is assumed to be a rate-limiting process and, hence, the interface kinetic is assumed to be instantaneous. For the interface-controlled regime, one assumes that diffusion is fast so that growth of a stable phase is limited by dissipative processes inside the narrow transition zone. The hypothesis of “normal growth” assumes that there is a linear dependence between the phase boundary velocity and the driving force [23].

The diffusion- and interface-controlled regimes have been studied analytically by Truskinovsky [23]. Here, we expand on the work of Truskinovsky [23] and consider five additional aspects: (1) We consider a phase transformation in a system with two chemical components, (2) we add the contribution of pressure-dependent mechanical mixing to the Gibbs energy of the system [17], (3) we apply a geologically relevant nonlinear equation of state in which the compressibility is a function of pressure [27,28,29], (4) we consider advection and large finite shear strains and (5) we solve the governing equations numerically in 1D and 2D employing a pseudo-transient finite difference method [30,31,32].

The aim of our study is to investigate the kinetics of two-component phase transformations which take place in a deforming system. To this end, we couple a Cahn–Hilliard-type diffusion model with a mechanical model for compressible viscous flow. We quantify the dependence of the velocity of the phase boundary, $v_f$, on various model parameters controlling the chemomechanical coupling, such as the difference in density between the two phases and the bulk compressibility of the phases. Particularly, we aim to quantify the mobility of the phase boundary and to test the hypothesis of the linear dependence between $v_f$ and the degree of metastability. Furthermore, we investigate the impact of deformation on the spatial patterns of phase transformations.

Table 1 Symbols and notation used in the text

2 Theoretical and numerical model

2.1 Mass and momentum conservation equations

We consider the time-dependent diffusion process in a system with two phases, each having two chemical components $\alpha$ and $\beta$, coupled with slow (no inertia) viscous fluid flow (Table  1). In the context of porous media modeling, in our model we resolve the two-phase flow on the pore scale.

The conservation equations for the mass of component (k) and for the total mass are [33]:

$$\begin{aligned}{} & {} \frac{\partial \rho c^k}{\partial t} + \nabla _i\left( \rho c^k v_i + \rho q^k_i\right) = 0, \quad k = \{\alpha , \beta \}, \end{aligned}$$

(1)

$$\begin{aligned}{} & {} \frac{\partial \rho }{\partial t} + \nabla _i\left( \rho v_i\right) = 0, \end{aligned}$$

(2)

where $\rho$ is the density, $c^k$ is the concentration (mass fraction) of component k, t is the time, $v_i$ is the barycentric velocity and $q^k_i$ is the diffusion flux of the component k in the barycentric reference frame. The subscripted index i represents the spatial directions.

For the conservation of linear momentum, we apply the equation [33, 34]:

$$\begin{aligned} \nabla _j\left( -p\delta _{ij} + \tau _{ij}\right) = 0, \end{aligned}$$

(3)

where p and $\tau _{ij}$ are the pressure (mean stress) of the fluid and the components of the deviatoric stress tensor, respectively, and $\delta _{ij}$ is the Kronecker delta. Here, we do not consider inertial forces, body forces arising due to gravity and additional body forces that can arise due to chemical concentration gradients [24].

2.2 Thermodynamically admissible model

To derive a closed system of governing equations, constitutive relations for diffusive fluxes $q^k_i$ and deviatoric stresses $\tau _{ij}$ are required. We apply here the concepts of classic irreversible thermodynamics (CIT, [33, 34]). Assuming local thermodynamic equilibrium in an open system, the variation in the internal energy per unit mass (in kg), u, is:

$$\begin{aligned} \textrm{d}u = T\textrm{d}s-p\textrm{d}(1/\rho ) + \mu ^\alpha ~ \textrm{d}c^\alpha + \mu ^\beta ~ \textrm{d}c^\beta , \end{aligned}$$

(4)

where T is the temperature, s is the entropy, $\mu ^\alpha$ and $\mu ^\beta$ are the chemical potentials of components $\alpha$ and $\beta$, respectively, and $c^\alpha$ and $c^\beta$ are the concentrations of components $\alpha$ and $\beta$, respectively. Since $c^\alpha + c^\beta = 1$, we use c and $1-c$ instead of $c^\alpha$ and $c^\beta$, respectively.

For the considered system, an entropy production term can be derived following the concepts of CIT by substituting the corresponding mass and momentum conservation equations as well as entropy balance laws into Eq. (4) for the internal energy [33, 34]. The entropy production term is [17, 35, 36]:

$$\begin{aligned} TQ_S = -q^\alpha _i \nabla _i\mu ^\alpha - q^\beta _i \nabla _i\mu ^\beta + \tau _{ij}e^\textrm{d}_{ij}, \end{aligned}$$

(5)

where $e^\textrm{d}_{ij}$ are the components of the viscous deviatoric strain rate tensor, $e^\textrm{d}_{ij} = 1/2(\nabla _i v_j + \nabla _j v_i) - 1/3\delta _{ij}\nabla _k v_k$. The term $\tau _{ij}e^\textrm{d}_{ij}$ represents the amount of dissipative mechanical work which is converted into heat. For the derivation of the entropy production term $TQ_S$, we refer the reader to chapter III, §2 in [36].

In a two-component system, the simplest definition for the diffusion flux $q^k_i$ that guarantees non-negativity of the entropy production as given by Eq. (5) while ensuring the Fickian limit for the ideal solution model is [24]:

$$\begin{aligned} q^k_i = -D c^k\nabla _i\mu ^k, \quad k = \{\alpha , \beta \}, \end{aligned}$$

(6)

where D is a concentration-dependent mobility coefficient [33].

Similar to the diffusion flux, we define the deviatoric stress tensor components $\tau _{ij}$ to be proportional to the deviatoric strain rate $e^\textrm{d}_{ij}$ to ensure non-negative entropy production $Q_S$ [35]:

$$\begin{aligned} \tau _{ij} = 2\eta e^\textrm{d}_{ij}, \end{aligned}$$

(7)

where $\eta$ is the effective viscosity of the fluid which is a positive number.

2.3 Constitutive relations and equation of state

We assume that the effective viscosity $\eta$ is a function of the concentration c and apply a simple equation of the form:

$$\begin{aligned} \log _{10}(\eta ) = \log _{10}(\eta _0) + \nu \cdot c, \end{aligned}$$

(8)

where $\eta _0$ and $\nu$ are a reference viscosity (for $c=0$) and a model parameter controlling the dependence of $\eta$ on c, respectively. Parameter $\nu$ is a number of orders of magnitude in a viscosity contrast between pure endmember substances.

The total density of the system is:

$$\begin{aligned} \rho = \rho ^\alpha c + \rho ^\beta (1-c), \end{aligned}$$

(9)

where $\rho ^\alpha$ and $\rho ^\beta$ are the densities of pure endmember components $\alpha$ and $\beta$, respectively. We apply a nonlinear equation of state to relate the density to the pressure using a nonlinear compressibility, $\frac{1}{K_0+K'p}$, which depends on p [27,28,29, 37]:

$$\begin{aligned} \frac{\textrm{d}\rho ^k}{\rho ^k} = \frac{\textrm{d} p}{K_0+K'p}, \quad k = \{\alpha , \beta \}, \end{aligned}$$

(10)

where $K_0$ and $K'$ are material parameters. Equation (10) can be integrated to yield a closed-form expression for the density:

$$\begin{aligned} \rho ^k = \rho ^k_0\left( \frac{K'}{K_0}p + 1\right) ^{1/K'}, \quad k = \{\alpha , \beta \}. \end{aligned}$$

(11)

A pressure-dependent compressibility is important for many geological and geodynamic applications which involve large pressures and often large spatial dimensions causing large pressure variations within the considered system.

2.4 Relation between chemical potential and concentration

2.4.1 Relation between Gibbs energy, chemical potential and concentration

To derive a closed system of equations, a constitutive relation between $\mu$ and c is required. Usually, this relation is formulated in terms of molar concentrations. In the presence of pressure gradients, it is important to consider chemical potentials in mass units to derive a thermodynamically admissible formulation for diffusion fluxes [17].

For the transparency of our derivation, we start with the common definition of chemical potentials and Gibbs energy in molar quantities. We introduce the molar fraction $\hat{c}$ that is related to the mass fraction c as follows:

$$\begin{aligned} \hat{c} = \frac{c}{M^\alpha }\Bigg /\left( \frac{c}{M^\alpha } + \frac{1-c}{M^\beta }\right) , \end{aligned}$$

(12)

where $M^\alpha$ and $M^\beta$ are the molar masses of the components $\alpha$ and $\beta$, respectively.

The generalized molar chemical potentials for components $\alpha$ and $\beta$ are [5, 24]:

$$\begin{aligned} \hat{\mu }^\alpha&= \hat{g} + (1-\hat{c})\frac{\delta \hat{g}}{\delta \hat{c}}, \end{aligned}$$

(13)

$$\begin{aligned} \hat{\mu }^\beta&= \hat{g} -\hat{c} \frac{\delta \hat{g}}{\delta \hat{c}}, \end{aligned}$$

(14)

where $\hat{g}$ is the molar Gibbs energy and $\delta \hat{g}/\delta \hat{c}$ is the variational partial derivative defined as [24]:

$$\begin{aligned} \frac{\delta \hat{g}}{\delta \hat{c}} = \frac{\partial \hat{g} }{\partial \hat{c}} - \nabla _i\frac{\partial \hat{g}}{\partial \nabla _i \hat{c}}. \end{aligned}$$

(15)

The variational derivative includes derivatives with respect to both concentration $\hat{c}$ and the spatial derivative of concentration $\nabla _i \hat{c}$ [21, 24].

Next, we introduce the Gibbs energy of the mixture, chemical potential and the diffusion flux of the first component per unit mass:

$$\begin{aligned} g = \frac{\hat{g}}{\hat{c}M^\alpha + (1-\hat{c})M^\beta }, \quad \mu = \frac{\hat{\mu }^\alpha }{M^\alpha },\quad q_i = q^\alpha _i. \end{aligned}$$

(16)

We use the Gibbs energy and chemical potential per unit mass to obtain a consistent coupling between “chemical” and “mechanical” quantities in the final governing system of equations.

2.4.2 Effective diffusivity

To understand the fundamental characteristics of the effective diffusivity, which is related to $q_i$ in the mass conservation Eq. (1), we assume for the moment that the Gibbs energy depends only on c and not on $\nabla _i c$, so that $\delta g/\delta c = \partial g/\partial c$. Assuming that $\mu$ depends on c, i.e., $\mu (c)$, and applying the chain rule of differentiation, we get:

$$\begin{aligned} q_i = -D c \nabla _i\mu = -D c \frac{\partial \mu }{\partial c}\nabla _i c = -D_\textrm{eff}\nabla _i c, \end{aligned}$$

(17)

where the effective diffusivity is

$$\begin{aligned} D_\textrm{eff} = D c \frac{\partial \mu }{\partial c} = Dc(1-c)\frac{\partial ^2 g}{\partial c^2}. \end{aligned}$$

(18)

The effective diffusivity $D_\textrm{eff}$, hence, depends on the concentration and on the second derivative w.r.t. the concentration c (the “curvature” in g–c phase space) of the Gibbs energy. Assuming that D is always a positive number, the value of $\partial ^2 g/\partial c^2$ must be positive in order to guarantee a stable diffusion of c which is controlled by the mass conservation Eq. (1). It is evident from substituting definitions of the diffusive flux and effective diffusivity, given by (17) and (18), respectively, into Eq. (1), which in the simplest case of $v_i = 0$ reduces to the parabolic diffusion equation:

$$\begin{aligned} \frac{\partial c}{\partial t} + \nabla _i\left( D_\textrm{eff}\nabla _i c\right) = 0. \end{aligned}$$

(19)

If $D_\textrm{eff}$ is negative, Eq. (19) becomes unstable. In the regions where $D_\textrm{eff} < 0$, it becomes necessary to stabilize the modeling of the “uphill” diffusion process by adding a gradient term into the governing equations [24].

It can be also useful to consider a characteristic diffusivity defined as:

$$\begin{aligned} D_\textrm{ch} = \frac{D RT}{M^\alpha }, \end{aligned}$$

(20)

where R is the universal gas constant. Such characteristic diffusivities have the same physical units as typical diffusivities, namely m$^2/$s in SI units. We use the characteristic diffusivities $D_\textrm{ch}$ in the dimensional analysis of the governing equations (see Sect. 2.5).

2.4.3 The Gibbs energy

The Gibbs energy per unit mass, g, of the simplest two-component system is composed of minimum required three energy contributions:

$$\begin{aligned} g&= g_\textrm{mech} + g_\textrm{mix} + g_\textrm{surf}, \end{aligned}$$

(21)

$$\begin{aligned} g_\textrm{mech}&= \int _{p_0}^p\left( \frac{c}{\rho ^\alpha } + \frac{1-c}{\rho ^\beta }\right) \textrm{d}p, \end{aligned}$$

(22)

$$\begin{aligned} g_\textrm{mix}&= \frac{RT}{M}\left( \hat{c}\ln (\hat{c}) + (1-\hat{c})\ln (1-\hat{c}) + \chi \hat{c}(1-\hat{c})\right) , \end{aligned}$$

(23)

$$\begin{aligned} g_\textrm{surf}&= \frac{RT}{M}\frac{\gamma ^2}{2}\left( \nabla c\right) ^2, \end{aligned}$$

(24)

where $\chi$ is an interaction parameter for the Flory–Huggins mixing model [38, 39] describing the deviation of the real material from the ideal mixing, and $\gamma$ is a gradient energy parameter [24]. Here, $g_\textrm{mech}$ is the mechanical term related to the Gibbs energy of the endmembers [17]; $g_\textrm{mix}$ is the free energy of mixing based on the Flory–Huggins model [24]; $g_\textrm{surf}$ is the non-local contribution [24]. The pressure dependence in $g_\textrm{mech}$ is in agreement with a recently published discussion on the thermodynamic admissibility of this formulation [17].

Our choice of the Gibbs energy of mixing $g_\textrm{mix}$ as a non-convex function of concentration $\hat{c}$ is different from the model developed by Truskinovsky [23], where a similar non-convex function was used to describe the dependence of the specific stored energy on deformation in a pure substance. Our formulation enables modeling of phase transformations in a reacting multi-component mixture.

In the following, we assume for simplicity that the molar masses of components $\alpha$ and $\beta$ are equal to each other, i.e., $M^\alpha = M^\beta$, and therefore, $\hat{c}M^\alpha + (1-\hat{c})M^\beta = M^\alpha$. The development of a coupled chemo-mechanical model for the non-equimolar case is an important direction for future research, but beyond the scope of our study.

2.5 Dimensionless form of governing equations and link between chemistry and mechanics

We transform the governing equations into a dimensionless form by using three independent characteristic scales for length, $l^*$, time, $t^*$, and pressure, $p^*$ (Table 2).

Table 2 Applied characteristic scales and dimensionless numbers

In the following, we use the same notation for the non-dimensional parameters as for the dimensional parameters. The dimensionless governing system of equations is:

$$\begin{aligned}{} & {} \frac{\partial \rho c}{\partial t} + \nabla _i\left( \rho [v_ic + q_i]\right) = 0, \end{aligned}$$

(25)

$$\begin{aligned}{} & {} \frac{\partial \rho }{\partial t} + \nabla _i\left( \rho v_i\right) = 0, \end{aligned}$$

(26)

$$\begin{aligned}{} & {} \nabla _j\left( -p\delta _{ij} + \tau _{ij}\right) = 0, \end{aligned}$$

(27)

$$\begin{aligned}{} & {} \tau _{ij} = \eta \left( \nabla _i v_j + \nabla _j v_i - \frac{2}{3}\delta _{ij}\nabla _k v_k\right) , \end{aligned}$$

(28)

$$\begin{aligned}{} & {} \log _{10}(\eta ) = \log _{10}(\eta _0) + \nu \cdot c, \end{aligned}$$

(29)

$$\begin{aligned}{} & {} q_i = - c\nabla _i(\mu ), \end{aligned}$$

(30)

$$\begin{aligned}{} & {} \nabla _i(\mu ) = \frac{\partial \mu }{\partial c}\nabla _i(c) + \frac{\partial \mu }{\partial p}\nabla _i(p) + \nabla _i(\gamma '^2\nabla ^2 c), \end{aligned}$$

(31)

$$\begin{aligned}{} & {} \frac{\partial \mu }{\partial c} = \frac{1}{c} - 2(1-c)\chi , \end{aligned}$$

(32)

$$\begin{aligned}{} & {} \frac{\partial \mu }{\partial p} = -\frac{\textrm{Ba}\,\delta _v}{c(\textrm{De}\,p + 1)^{1/K'}}, \end{aligned}$$

(33)

$$\begin{aligned}{} & {} \rho = (1 + \delta _v\,c)(\textrm{De}\,p + 1)^{1/K'}. \end{aligned}$$

(34)

The non-dimensionalization resulted in five dimensionless numbers, which we term: $\textrm{De}$, $\delta _v$, $\textrm{Ba}$, $\gamma '$ and $\textrm{Pe}$. The definitions of these numbers are given in Table 2. $\textrm{De}$ stands for Deborah number, $\textrm{Ba}$ is a dimensionless number that quantifies the dependence of chemical potential on pressure, and $\textrm{Pe}$ stands for Péclet number. All dimensionless numbers appear in the system of equations, except $\textrm{Pe}$ which is considered by the boundary conditions for the velocities which are specified at the outer model boundary.

The system of equations shows the coupling between the considered chemical and mechanical processes. For example, the chemical potential difference, $\mu$, depends on both c and p. Furthermore, the density, $\rho$, also depends on both c and p. Moreover, in the conservation equation for the mass of component $\alpha$ there is an advection term and the velocity, $v_i$, is controlled by the viscous flow.

2.6 Numerical method

We study the kinetics of first-order phase transformations with a series of numerical simulations on a 1D axisymmetric domain with radius $r\in [R_\textrm{min};R_\textrm{max}]$ and on a 2D axisymmetric domain $(r,\phi )\in [R_\textrm{min}; R_\textrm{max}] \times [0; 2\pi ]$.

To closer match the numerical discretization and the axisymmetrical domain, we use the polar coordinate system. In 2D polar coordinates, Eqs. (25)–(31) take the following form:

$$\begin{aligned}{} & {} \frac{\partial \rho c}{\partial t} + \frac{1}{r}\left[ \frac{\partial }{\partial r}\left( r\rho [v_r c + q_r]\right) + \frac{\partial }{\partial \phi }\left( \rho [v_\phi c + q_\phi ]\right) \right] = 0, \end{aligned}$$

(35)

$$\begin{aligned}{} & {} \frac{\partial \rho }{\partial t} + \frac{1}{r}\left[ \frac{\partial }{\partial r}\left( r\rho v_r\right) + \frac{\partial }{\partial \phi }\left( \rho v_\phi \right) \right] = 0, \end{aligned}$$

(36)

$$\begin{aligned}{} & {} -\frac{\partial p}{\partial r} + \frac{\partial \tau _{rr}}{\partial r} + \frac{1}{r}\frac{\partial \tau _{r\phi }}{\partial \phi } + \frac{\tau _{rr} - \tau _{\phi \phi }}{r} = 0, \end{aligned}$$

(37)

$$\begin{aligned}{} & {} -\frac{1}{r}\frac{\partial p}{\partial \phi } + \frac{1}{r}\frac{\partial \tau _{\phi \phi }}{\partial \phi } + \frac{\partial \tau _{r\phi }}{\partial r} + 2\frac{\tau _{r\phi }}{r} = 0, \end{aligned}$$

(38)

$$\begin{aligned}{} & {} \tau _{rr} = 2\eta \left( \frac{\partial v_r}{\partial r} - \frac{1}{3}\left[ \frac{1}{r}\frac{\partial }{\partial r}\left( r v_r\right) + \frac{1}{r}\frac{\partial u_\phi }{\partial \phi }\right] \right) , \end{aligned}$$

(39)

$$\begin{aligned}{} & {} \tau _{\phi \phi } = 2\eta \left( \frac{1}{r}\left[ \frac{\partial v_\phi }{\partial \phi } + v_r\right] - \frac{1}{3}\left[ \frac{1}{r}\frac{\partial }{\partial r}\left( r v_r\right) + \frac{1}{r}\frac{\partial v_\phi }{\partial \phi }\right] \right) , \end{aligned}$$

(40)

$$\begin{aligned}{} & {} \tau _{r\phi } = \eta \left( \frac{1}{r}\frac{\partial v_r}{\partial \phi } + \frac{\partial v_\phi }{\partial r} - \frac{v_\phi }{r} \right) , \end{aligned}$$

(41)

$$\begin{aligned}{} & {} q_r = - c\left( \frac{\partial \mu }{\partial c}\frac{\partial c}{\partial r} + \frac{\partial \mu }{\partial p}\frac{\partial p}{\partial r} - \gamma ^2\frac{\partial }{\partial r} \left[ \frac{1}{r}\frac{\partial }{\partial r}\left( r\frac{\partial c}{\partial r}\right) + \frac{1}{r^2}\frac{\partial ^2 c}{\partial \phi ^2}\right] \right) , \end{aligned}$$

(42)

$$\begin{aligned}{} & {} q_\phi = - \frac{c}{r}\left( \frac{\partial \mu }{\partial c}\frac{\partial c}{\partial \phi } + \frac{\partial \mu }{\partial p}\frac{\partial p}{\partial \phi } - \gamma ^2\frac{\partial }{\partial \phi } \left[ \frac{1}{r}\frac{\partial }{\partial r}\left( r\frac{\partial c}{\partial r}\right) + \frac{1}{r^2}\frac{\partial ^2 c}{\partial \phi ^2}\right] \right) . \end{aligned}$$

(43)

For 1D models, the angular derivatives $\partial /\partial \phi$ are ignored.

To solve the coupled system (35)–(43), we employ the accelerated pseudo-transient (PT) method described in details in [30, 32]. In accordance with this method, the steady-state solution to the mechanical problem is obtained as a limit of a transient physical process. For the details of the iterative PT numerical method and reference implementations for the mechanical and diffusion problems, we refer the reader to the work by Räss et al. [32].

We augment the right-hand sides of Eqs. (35)–(38), with the pseudo-time derivatives:

$$\begin{aligned}{} & {} \frac{\partial \rho c}{\partial t} + \frac{1}{r}\left[ \frac{\partial }{\partial r}\left( r\rho [v_r c + q_r]\right) + \frac{\partial }{\partial \phi }\left( \rho [v_\phi c + q_\phi ]\right) \right] = \widetilde{\rho }_c\frac{\partial c}{\partial \tau }, \end{aligned}$$

(44)

$$\begin{aligned}{} & {} \frac{\partial \rho }{\partial t} + \frac{1}{r}\left[ \frac{\partial }{\partial r}\left( r\rho v_r\right) + \frac{\partial }{\partial \phi }\left( \rho v_\phi \right) \right] = \frac{\rho }{\widetilde{K}}\frac{\partial p}{\partial \tau }, \end{aligned}$$

(45)

$$\begin{aligned}{} & {} -\frac{\partial p}{\partial r} + \frac{\partial \tau _{rr}}{\partial r} + \frac{1}{r}\frac{\partial \tau _{r\phi }}{\partial \phi } + \frac{\tau _{rr} - \tau _{\phi \phi }}{r} = \widetilde{\rho }_v\frac{\partial v_r}{\partial \tau }, \end{aligned}$$

(46)

$$\begin{aligned}{} & {} -\frac{1}{r}\frac{\partial p}{\partial \phi } + \frac{1}{r}\frac{\partial \tau _{\phi \phi }}{\partial \phi } + \frac{\partial \tau _{r\phi }}{\partial r} + 2\frac{\tau _{r\phi }}{r} = \widetilde{\rho }_v\frac{\partial v_\phi }{\partial \tau }. \end{aligned}$$

(47)

Here, $\partial /\partial \tau$ is the pseudo-time derivative, $\widetilde{\rho }_c$ and $\widetilde{\rho }_v$ are the numerical pseudo-densities for the diffusion and mechanical problems, respectively, and $\widetilde{K}$ is the numerical pseudo-bulk modulus.

The pseudo-transient method defined by Eqs. (44)–(47) leads to the coupled system of equations which could be integrated explicitly in pseudo-time until the pseudo-time derivatives become smaller than a predefined tolerance. The values of numerical parameters are chosen to satisfy the stability conditions for the explicit time integration scheme. Equations (44), (43) and (42) describe a nonlinear Cahn–Hilliard type diffusion, and in combination these equations involve the fourth order spatial derivative of c. Theoretically, the fourth-order derivative leads to the maximum pseudo-time step to be proportional to the fourth power of grid spacing $\Delta r^4$. However, the coefficient $\gamma$ is comparable to the $\Delta r$ which effectively leads to only quadratic dependence on grid spacing, i.e., $\Delta r^2$. Therefore, the number of iterations, i.e., the number of pseudo-time steps required for convergence, grows quadratically with increasing resolution, which limits the practical application of the PT method. We accelerate the convergence of the pseudo-transient method by introducing additional pseudo-time derivatives into the constitutive equations for stress components $\tau _{ij}$ defined by Eqs. (39)–(41), thus converting the viscous rheology into a Maxwell viscoelastic rheology in pseudo-time:

$$\begin{aligned} \frac{1}{2\widetilde{G}}\frac{\tau _{rr}}{\partial \tau } + \frac{\tau _{rr}}{2\eta }= & {} \frac{\partial v_r}{\partial r} - \frac{1}{3}\left[ \frac{1}{r}\frac{\partial }{\partial r}\left( r v_r\right) + \frac{1}{r}\frac{\partial u_\phi }{\partial \phi }\right] , \end{aligned}$$

(48)

$$\begin{aligned} \frac{1}{2\widetilde{G}}\frac{\tau _{\phi \phi }}{\partial \tau } + \frac{\tau _{\phi \phi }}{2\eta }= & {} \frac{1}{r}\left[ \frac{\partial v_\phi }{\partial \phi } + v_r\right] - \frac{1}{3}\left[ \frac{1}{r}\frac{\partial }{\partial r}\left( r v_r\right) + \frac{1}{r}\frac{\partial v_\phi }{\partial \phi }\right] , \end{aligned}$$

(49)

$$\begin{aligned} \frac{1}{2\widetilde{G}}\frac{\tau _{r\phi }}{\partial \tau } + \frac{\tau _{r\phi }}{2\eta }= & {} \frac{1}{2}\left( \frac{1}{r}\frac{\partial v_r}{\partial \phi } + \frac{\partial v_\phi }{\partial r} - \frac{v_\phi }{r}\right) . \end{aligned}$$

(50)

Here, $\widetilde{G}$ is a numerical parameter which can be interpreted as numerical pseudo-elastic shear modulus that controls pseudo-relaxation time.

Similar to the definitions of stresses, we convert the flux definitions (42) and (43) to the following transient equations:

$$\begin{aligned} \widetilde{\theta }\frac{\partial q_r}{\partial \tau } + q_r&= - c\frac{\partial \mu }{\partial r}, \end{aligned}$$

(51)

$$\begin{aligned} \widetilde{\theta }\frac{\partial q_\phi }{\partial \tau } + q_\phi&= - \frac{c}{r}\frac{\partial \mu }{\partial \phi }, \end{aligned}$$

(52)

where $\widetilde{\theta }$ is a relaxation pseudo-time for diffusion fluxes. The values of parameters $\widetilde{\theta }$ and $\widetilde{G}$ can be freely chosen to achieve fastest convergence. For linear equations, the optimal values for $\widetilde{G}$ and $\widetilde{\theta }$ can be determined from a dispersion analysis of the governing equations as reported in [32]. In the linear case, the number of iterations required for convergence grows linearly with increasing grid resolution. For significantly nonlinear coupled systems of equations used in this work, values of these parameters are calibrated empirically.

Fig. 1

Spatial arrangement of unknowns and related quantities on a staggered grid in polar coordinates

We discretize the computational domain on a polar grid with $n_r$ and $n_r \times n_\phi$ cells for the 1D and 2D cases, respectively. In all results presented in this work, we set $n_r = 500$ and $n_\phi = 1570$. We use the staggered grid to store field variables (Fig. 1). We approximate the spatial derivatives using conservative finite differences.

We implement the described PT method in a numerical code written in Julia language [40]. We use the ParallelStencil.jl package [41] to automatically generate the numerical code for different hardware architectures, including multi-core CPUs and graphics processing units (GPUs). The source code of the implementation is available on Github [42] and in a persistent data repository [43].

2.7 Model configuration

We perform 1D and 2D numerical simulations. The 1D numerical model employs an axisymmetric formulation and the 2D model is based on polar coordinates. Both 1D and 2D models apply to a circular domain. The internal boundary of the model is located at the radial position $R_\textrm{min}$ and the outer boundary at $R_\textrm{max}$ (Fig. 2a).

Fig. 2

Initial conditions for tree model setups and common thermodynamic data for all three model setups. Panel a initial radial concentration profiled for three different model setups; panel b relation between c and Gibbs energy g. g is the integral of $\mu$ (panel c) along c; panel c: corresponding constitutive relation between $\mu$ and c

We perform 1D and 2D simulations with different initial and boundary conditions as well as with different model parameters. For all simulations, the initial pressure, p, is always zero. For all 1D simulations, the velocities at the interior and outer model boundaries $r = R_\textrm{min}$ and $r = R_\textrm{max}$ are always zero. For all 2D simulations, the velocities at the interior model boundary $r = R_\textrm{min}$ are always zero and the velocities at the outer boundary can be different depending on the applied far-field deformation. We consider two modes of far-field deformation: pure shear and simple shear. In pure shear, we specify compression in the horizontal direction and tension in the vertical direction:

$$\begin{aligned} v_r = -\textrm{Pe}\,r\,\cos (2\phi ), \quad v_\phi = \textrm{Pe}\,r\,\sin (2\phi )/2 \quad \textrm{at} \quad r = R_\textrm{max}. \end{aligned}$$

(53)

Here, $\textrm{Pe}$ is the Peclet number which characterizes the ratio between advective and diffusive fluxes (Table 2).

In simple shear, we specify the horizontally oriented velocity profile:

$$\begin{aligned} v_r = \textrm{Pe}\,r\,\sin (2\phi )/2, \quad v_\phi = -\textrm{Pe}\,r\,\sin (2\phi )/2 \quad \textrm{at} \quad r = R_\textrm{max}. \end{aligned}$$

(54)

The particular additional initial conditions and model parameters for the performed simulations will be given in the results section with the description of the simulation results.

3 Results

3.1 1D results without mechanical coupling

For the first series of 1D simulations, the mechanical coupling is not considered. The mechanical coupling is removed by setting chemo-mechanical coupling parameters $\textrm{Ba}$ and $\nu$ to 0. The initial profile of c is shown in Fig. 2a (solid black line). The step in the initial c-profile indicates the initial location of the phase boundary.

Fig. 3

Results of 1D simulations with zero chemo-mechanical coupling parameters $\textrm{Ba}$ = 0 and $\nu = 0$. Panel a Mobility as function of Flory–Huggins parameter,$\chi$, and energy gradient parameter,$\gamma$. $[\mu ]$ is the difference of $\mu$ at the right model side minus $\mu$ at the left model side. The mobility is the ratio of $v_f / [\mu ]$. Panel b Data collapse showing all calculated values of $v_f / [\mu ]$ plotted versus the ratio of $\gamma / \chi ^3$. The values of other dimensionless parameters used in the study: $\delta _v$ = 0.05, $\textrm{De} = 0.0001$, $\textrm{Pe} = 0.01$, $K' = 0$

In the internal model domain (the internal phase), the initial value of c is 0.9 and in the outer domain (the outer phase) the initial value of c is 0.2 (Fig. 2a). Figure 2b shows the corresponding initial profile of g versus the initial profile of c, and Fig. 2c shows the initial profile of $\mu$ versus c. The value of c for the internal phase corresponds to the equilibrium value since it corresponds to the value for which g is minimal. The value of c for the outer phase corresponds to a value in the metastable domain.

We perform a series of 1D simulations with different values for $\chi$ and $\gamma$. From the numerical simulations, we determine the values of the velocity of the moving phase boundary, $v_f$. We calculate the value of $[\mu ]$ which is the absolute value of the difference between the values of $\mu$ at the internal and outer model boundary.

The ratio of $v_f / [\mu ]$ represents the mobility of the phase boundary, and this mobility varies with varying values of $\chi$ and $\gamma$ (Fig. 3a). Plotting the values of $v_f / [\mu ]$ for all 1D simulations versus the ratio of $\gamma / \chi ^3$ shows that all values of $v_f / [\mu ]$ lie within a narrow region (Fig. 3b). For larger values of $\gamma / \chi ^3$ values of $v_f / [\mu ]$ are essentially constant while for smaller values of $\gamma / \chi ^3$, the values of $v_f / [\mu ]$ increase with increasing values of $\gamma / \chi ^3$. For larger values of $\gamma / \chi ^3$, the mobility is essentially constant and the phase transformation is controlled by the diffusion-controlled regime, since the mobility is independent of $\gamma$.

3.2 1D results with mechanical coupling

We perform four 1D simulations with different chemo-mechanical coupling. The coupling is controlled by the two parameters $\textrm{Ba}$ and $\nu$, which control the coupling with the pressure, p, and the shear viscosity, $\eta$. Figure 4 shows the time evolution of the c-profiles for four simulations involving different combinations of $\textrm{Ba}$ and $\nu$. Generally, the evolution of the four c-profiles is similar. However, the simulations show that the phase boundary moves slower if the value of $\textrm{Ba}$ is not zero (here $\textrm{Ba}$ = 40), that is, if the chemo-mechanical coupling via the pressure is activated. The slowest movement of the phase boundary is observed if both couplings via p and $\eta$ are activated, here for $\textrm{Ba} = 40$ and $\nu = 1$ (Fig. 4d).

Fig. 4

Time evolution of the c profiles for four combinations of model parameters. Panel a $\eta _0=0$, $\textrm{Ba}=0$; panel b $\nu =0$, $\textrm{Ba}=40$; panel c $\eta _0=1$, $\textrm{Ba}=0$; panel d: $\nu =1$, $\textrm{Ba}=40$

Figure 5 shows the evolution of $v_f$ versus $[\mu ]$ with progressive time for two simulations. During the initial stages of the simulations, values of $v_f$ and $[\mu ]$ are largest. After some time, values of $[\mu ]$ start to decrease toward the equilibrium value of $[\mu ]$ = 0. The values of $v_f$ decrease linearly with decreasing values of $[\mu ]$. During this stage of decreasing values of $[\mu ]$, the value of the mobility, $v_f / [\mu ]$, is approaching a constant value.

Fig. 5

Representative plots of $v_f$ and $[\mu ]$ and their ratio versus time for two simulations with different chemo-mechanical coupling parameter $\textrm{Ba}$. With progressive time, the values of $[\mu ]$ and $v_f$ decrease toward equilibrium (i.e., $v_f =[\mu ] = 0$). After some time, $v_f$ becomes a linear function of $[\mu ]$. The values of other dimensionless parameters used for these two simulations are: $\delta _v$ = 0.05, $\nu = 0$, $\textrm{De} = 0.01$, $\textrm{Pe} = 0$, $\chi = 2.6$, $\gamma ' = 2\times 10^{-2}$, $K' = 4$

We performed a systematic series of numerical simulations with different values of key parameters $\delta _v$ and $\textrm{Ba}$ characterizing chemo-mechanical coupling to assess their impact on the mobility $v_f/[\mu ]$ (Fig. 6), including the two simulations shown in Fig. 5. We plot the values of the mobility as a function of $\delta _v$ (Fig. 6a), color coding marks the values of $\textrm{Ba}$. Scatter of the values is significantly reduced when plotting the mobility, $v_f/[\mu ]$, as a function of $\delta _v\cdot \textrm{Ba}$ (Fig. 6b). These results quantify the phenomena of “mechanical closure” [26], i.e., the slowing down of the transformation rate in systems with large volumetric effect under mechanical resistance to the volume change.

Fig. 6

Mobility coefficient as a function of volumetric effect of reaction $\delta _v$ and chemo-mechanical coupling parameter $\textrm{Ba}$. The values of other dimensionless parameters used in the study: $\nu = 0$, $\textrm{De} = 0.01$, $\textrm{Pe} = 0$, $\chi = 2.6$, $\gamma ' = 2\times 10^{-2}$, $K' = 4$. Panel a data of the 16 simulations; panel b data collapse to the master curve of mobility coefficient as a function of product of $\delta _v$ and $\textrm{Ba}$. Color coding marks the values of $\textrm{Ba}$ parameter. Note the decrease in the mobility coefficient with increase in values of the coupling parameter and the volumetric effect of reaction

3.3 2D results

To investigate the characteristic patterns of phase separation, we first perform a 2D simulation without mechanical coupling, i.e., only the diffusion problem is solved. We specify as initial condition a random distribution of concentration c, given by:

$$\begin{aligned} c(r,\phi ) = c_\textrm{sys} + c_\textrm{a}\textrm{rand}(s, r\cos (\phi ),r\sin (\phi ))\quad \textrm{at} \quad t = 0. \end{aligned}$$

(55)

where $c_\textrm{sys}$ is the average initial concentration, $c_\textrm{a}$ is the perturbation magnitude, s is the characteristic size of an area with a constant value and the function $\textrm{rand}$ returns a random field with values in the range $[-1; 1)$. In all simulations, we set $s = 0.004$ and $c_a = 0.01$.

Fig. 7

Representative distribution of concentration field c for spinodal (panel a), nucleation (panel b), and “normal growth” (panel c) model setups at early stages and later (panels d, e and f) stages of phase transformation

We consider two random initial distributions of c: for the first distribution, $c_\textrm{sys}$ is set to 0.5. For this distribution, the system is in an unstable state and phase transformations occur simultaneously in the entire domain (Figs. 2a, 7a and d). We term this distribution the spinodal mode distribution.

For the second distribution, $c_\textrm{sys}$ is set to 0.25. For this distribution, the average concentration is shifted toward the equilibrium value for the first phase (Figs. 2a and 7b,e). We term this distribution the nucleation-mode distribution.

In addition, we perform one 2D simulation with initial distribution of c similar to the one in the 1D model setup, i.e., a step-like distribution of c described in Sect. 3.1 (Figs. 2a, 7c and f).

Table 3 Values of dimensionless parameters used in the 2D numerical simulations

We then perform a series of 2D simulations featuring chemo-mechanical coupling. The values of the dimensionless parameters applied in the simulations are summarized in Table 3.

Fig. 8

Snapshots of three 2D simulations showing color plots of pressure, p, and concentration, c (black and white colors) at $t = 0.15$. Panel a simulation without any applied far-field shear deformation; panel b simulation with an applied far-field pure shear deformation; panel c simulation with an applied far-field simple shear deformation

We first run 2D simulations that have along their radial directions an initial step-like c-profile similar to the 1D simulations (Fig. 8). Therefore, in these 2D simulations the initial distribution of c is independent of the angular direction, specified by the angle $\phi$. We performed three simulations: without any applied far-field shear deformation (Fig. 8a), with an applied far-field pure shear deformation (Fig. 8b) and with an applied far-field simple shear deformation (Fig. 8c). The resulting distributions of p and c are different for the three simulations showing that the applied far-field deformation has a significant impact on the evolving phase transformation.

Fig. 9

Time evolution of density field for a 2D simulation with far-field simple shear and an initial random perturbation of c, corresponding to the nucleation mode (see text). Panels a (initial situation) to d show four stages of the model evolution at non-dimensional times $t = 10^{-4}$, $10^{-3}$, $2\times 10^{-3}$, and $3\times 10^{-3}$, respectively. Color plot shows the dimensionless density

In Fig. 9, we show four snapshots at different times which show the distribution of the dimensionless density field for a simulation with far-field simple shear deformation and an initial nucleation-mode distribution of c. The value of the density is controlled by both the concentration and the pressure and is, hence, a quantity which is affected by both chemical and mechanical processes. The results show the formation and evolution of string-like patterns. For this simulation, the distributions of viscosity, $\eta$, pressure, p, and concentration, c, show similar patterns (Fig. 10a, c and d). Higher values of c correspond to higher values of $\eta$ and p.

Fig. 10

Representative results for a 2D simulation with far-field simple shear and an initial random perturbation of c corresponding to the nucleation mode at non-dimensional time $t = 3\times 10^{-3}$ (see text). The four panels show color plots of viscosity (panel a), the absolute magnitude of velocity, $(v_1^2+v_2^2)^{0.5}$ (panel b), pressure (panel c) and concentration (panel d). Black arrows in panel c indicate the velocity field

Figure 10b shows the difference between the absolute magnitude of the velocity and the velocity corresponding to homogeneous simple shear that is for simple shear not affected by the rigid internal boundary (visible due to the white quarter-circle). The velocity field deviates strongest from the applied far-field simple shear around the rigid internal model boundary. The string-like regions are oriented with an angle of approximately 45 degrees to the horizontal direction. The orientation of the strings corresponds to the direction of the maximal compressive stress for the applied far-field simple shear for which the upper regions are moving horizontally toward the left side.

Fig. 11

Representative results for a 2D simulation with far-field pure shear and an initial random perturbation of c corresponding to the nucleation mode at non-dimensional time $t = 3\times 10^{-3}$ (see text). The four panels show color plots of viscosity (panel a), the absolute magnitude of velocity, $(v_1^2+v_2^2)^{0.5}$ (panel b), pressure (panel c) and concentration (panel d). Black arrows in panel c indicate the velocity field

We also perform this simulation for a far-field pure shear deformation with horizontal shortening and vertical extension (Fig. 11). Generally, the correlations between the distributions of $\eta$, p and c are the same as for the simulation with simple shear. However, the orientation of the string-like regions with high values of c has a different orientation compared to the ones for a simulation with far-field simple shear. For pure shear, the string-like regions are approximately parallel to the horizontal shortening direction which corresponds to the direction of maximal compressive stress.

Fig. 12

Black and white plots of concentration c for four simulations with far-field simple shear deformation and an initial random distribution of c at $t=3\times 10^{-3}$. Panel a spinodal mode (see text) initial c distribution without mechanical coupling, only advection is activated in mass conservation equation. Panel b spinodal mode initial c distribution with mechanical coupling. Panel c nucleation mode (see text) initial c distribution without mechanical coupling, only advection is activated in mass conservation equation. Panel d nucleation mode initial c distribution with mechanical coupling

We further performed three additional simulations for far-field simple shear with different initial distributions of c and with or without chemo-mechanical coupling (Fig. 12). For the simulations without chemo-mechanical coupling, only the advection caused by the far-field simple shear is considered, but couplings due to pressure and viscosity are not activated. For a spinodal-mode initial c distribution, the chemo-mechanical coupling has a considerable impact on the evolving distribution of c (Fig. 12a and b). Without chemo-mechanical coupling, the c distribution shows a distribution typical for a spinodal decomposition without any preference for a spatial direction (Fig. 12a). In contrast, with chemo-mechanical coupling the c distribution exhibits string-like patterns whereby the orientation of the strings is oriented approximately 45 degrees to the horizontal direction (Fig. 12b). For a nucleation-mode initial c distribution, the differences in the developed c distributions are even more significant between the simulation with (Fig. 12d) and without chemo-mechanical coupling (Fig. 12c).

Fig. 13

Black and white plots of concentration c for four simulations with far-field pure shear deformation and an initial random distribution of c at $t=3\times 10^{-3}$. Panel a spinodal mode (see text) initial c distribution without mechanical coupling, only advection is activated in mass conservation equation. Panel b spinodal mode initial c distribution with mechanical coupling. Panel c nucleation mode (see text) initial c distribution without mechanical coupling, only advection is activated in mass conservation equation. Panel d nucleation mode initial c distribution with mechanical coupling

Finally, we performed three additional simulations for far-field pure shear with different initial distributions of c and with or without chemo-mechanical coupling (Fig. 13). Generally, the impact of the initial c distribution and the chemo-mechanical coupling is similar compared to results of the simulations with far-field simple shear. The main difference is that the string-like patterns have a different orientation and they are approximately parallel to the horizontal shortening direction.

4 Discussion

For a binary system, the equilibrium state requires that the spatial gradients of the two chemical potentials are zero everywhere in the system. The chemical potential depends on pressure; therefore, in the presence of pressure gradients it is possible that the concentration, or mass fraction, profiles in the phases at equilibrium could also exhibit a spatial variation [17]. The kinetics of a phase transformation can be described by the diffusive redistribution of the concentrations of the two components. The two concentrations are not independent because their sum must be one. At the phase boundary, thermodynamic equilibrium requires the equality of both chemical potentials. A diffusive redistribution, constrained by the conservation of mass, can generate a concentration profile so that the chemical potential of the first component is equal at the phase boundary. However, the chemical potentials of the second component may still be different at the boundary. To equate also the chemical potentials of the second component, the phase boundary must move in order to change the absolute value of the concentration while conserving mass. Phenomenological proportionality constants are needed to quantify the velocity and “mobility” of the phase boundary and the intensity of mass exchange between phases. Here, we consider as mobility the ratio of the phase boundary velocity to the degree of metastability, which depends on the spatial variation in the chemical potential difference across the system.

Our simulations show that chemo-mechanical coupling has a considerable impact on phase transformation kinetics. The dimensionless number $\textrm{Ba}$ controls the intensity of the mechanical coupling and the larger the value of $\textrm{Ba}$ the stronger the mechanical coupling. The results presented in Fig. 6b show that an increase in $\textrm{Ba}$ causes a decrease in the mobility of the phase boundary and, hence, a slower phase transformation kinetics.

Furthermore, the applied far-field deformation and the orientation of the corresponding stress field impact the orientation of the evolving phase boundaries. The orientation of the maximal compressive stress for pure shear, with horizontal shortening, and simple shear, with horizontal shearing, varies by an angle of 45 degrees. The string-like regions with high values of c and p show such orientation difference in simulations with far-field simple shear and pure shear. Therefore, our results show that chemo-mechanical coupling can significantly impact the kinetics of phase transformation and the orientation of the evolving phase boundaries.

Our results are general and can be applied to study phase transformations kinetics in deforming systems involving various materials. Particular examples for potential applications are metamorphic reactions and phase transformations occurring in deforming rock [31, 44,45,46]. Although many observations of natural rock suggest the importance of chemo-mechanical coupling in deforming rock, the chemical and mechanical processes are frequently studied independently. Our results show the importance of chemo-mechanical coupling for the kinetics of phase transformations in deforming systems. Therefore, our results have important implications for the interpretation of phase transformations observed in deformed or stressed rocks. Moreover, our mathematical model can be elaborated in the future to consider additionally porous flow related to dehydration or melting reactions [31, 47].

Our 1D results for the time evolution of a phase boundary shown in Fig. 4 are conceptually similar to 1D numerical results of the interaction of acoustic wave propagating through a phase boundary presented in [11], where a diffuse interface approach with a non-convex double-well potential function for the free energy as a function of density. In [11], the case of inviscid pure substance, i.e., involving only one chemical component, is considered. Truskinovsky [23] studied the “normal growth” scenario for a propagating phase boundary in the same pure substance, but in the elastostatic limit.

We follow the work by Truskinovsky [23], but we consider the inertial-less viscous fluid described by Stokes equations, and we extend Truskinovsky’s analysis to a binary system, using a non-convex double-well dependence on the concentration of the first component in a solution model.

We use the principles of classic irreversible thermodynamics (CIT) to derive the expressions for diffusive fluxes $q^k_i$ and deviatoric stresses $\tau _{ij}$ and to constrain the form of the constitutive relations (Sect. 2.2). Numerous mixture theories have been developed on the basis of rational thermodynamics [48,49,50], and subsequently extended to multiphase systems [8,9,10]. Our model is in agreement with the mixture theory of multi-phase systems presented in [10, 49], yielding equivalent formulation. We use a diffuse interface approach, and thus, we don’t use a separate description for jump conditions on the interface obtained by Dell’Isola [10]. Validation of these jump conditions is outside the scope of this study.

5 Conclusion

We present a mathematical model for the kinetics of first-order phase transformations in a system with two chemical components which is coupled to compressible viscous flow. The bulk compressibility is a nonlinear function of the pressure, and the shear viscosity is a nonlinear function of the concentration.

The performed 1D simulations show that the velocity of the phase boundary, $v_f$, increases linearly with increasing degree of metastability, which is quantified by the difference (between outer and inner model boundaries) of the chemical potential difference of the two components, $[\mu ]$. Our results, hence, confirm the hypothesis of “normal growth” of phase boundaries for the investigated coupled chemo-mechanical system. For systems without chemo-mechanical coupling, the mobility of the phase boundary, quantified by the ratio of $v_f / [\mu ]$, becomes essentially constant for larger values of the surface energy parameter $\gamma$. For such large values of $\gamma$, the phase transformation occurs in the diffusion-controlled regime. For systems with chemo-mechanical coupling, the mobility of the phase boundary decreases with increasing intensity of chemo-mechanical coupling, and phase transformations shift toward the mechanics-controlled regime.

The performed 2D simulations, with a random initial distribution in concentration, c, show that the applied far-field deformation, either simple shear or pure shear, exerts a strong control on the evolving phase boundaries. Mechanical coupling and far-field deformation cause the formation of string-like patterns in the phase distribution. The orientation of the string-like patterns is controlled by the applied far-field deformation and corresponding stress field.