Multi-component chemo-mechanics based on transport relations for the chemical potential

TL;DR

This paper develops a chemo-mechanical model for multi-component materials, introducing a semi-analytical inversion of chemical potential to improve numerical stability and applying it to spinodal decomposition with elastic and plastic effects.

Contribution

A novel semi-analytical inversion method for chemical potential in multi-component chemo-mechanical models, enhancing numerical stability and enabling detailed spinodal decomposition simulations.

Findings

Improved numerical conditioning of transport equations.

Successful demonstration of convergence to Cahn-Hilliard solutions.

Insights into elastic and plastic effects on spinodal morphologies.

Abstract

A chemo-mechanical model for a finite-strain elasto-viscoplastic material containing multiple chemical components is formulated and an efficient numerical implementation is developed to solve the resulting transport relations. The numerical solution relies on inverting the constitutive model for the chemical potential. In this work, a semi-analytical inversion for a general family of multi-component regular-solution chemical free energy models is derived. This is based on splitting the chemical free energy into a convex contribution, treated implicitly, and a non-convex contribution, treated explicitly. This results in a reformulation of the system transport equations in terms of the chemical potential rather than the composition as the independent field variable. The numerical conditioning of the reformulated system, discretised by finite elements, is shown to be significantly improved, and convergence to the Cahn-Hilliard solution is demonstrated for the case of binary spinodal decomposition. Chemo-mechanically coupled binary and ternary spinodal decomposition systems are then investigated to illustrate the effect of anisotropic elastic deformation and plastic relaxation of the resulting spinodal morphologies in more complex material systems.

Tables

Table 1. Table 1: Chemo-mechanical material parameters for the binary spinodal system. Ω Ω \Omega is the molar volume, E 1 s ​ o ​ l superscript subscript 𝐸 1 𝑠 𝑜 𝑙 E_{1}^{sol} is the component solution energy, E 11 i ​ n ​ t superscript subscript 𝐸 11 𝑖 𝑛 𝑡 E_{11}^{int} is the component interaction energy, θ 𝜃 \theta is the temperature, κ 1 subscript 𝜅 1 \kappa_{1} is the gradient coefficient, α 1 subscript 𝛼 1 \alpha_{1} is the penalty parameter, M 1 subscript 𝑀 1 M_{1} is the component mobility, γ ˙ 0 subscript ˙ 𝛾 0 \dot{\gamma}_{0} is the reference shear rate, n 𝑛 n is the strain-rate sensitivity exponent, g 0 subscript 𝑔 0 g_{0} and g ∞ subscript 𝑔 g_{\infty} are the initial and asymptotic slip resistances, and h 0 subscript ℎ 0 h_{0} and h α ​ β subscript ℎ 𝛼 𝛽 h_{\alpha\beta} are hardening parameters.

Chemical energy	$\Omega$ (${\text{m}}^3{\text{mol}}^{-1}$)	$E_1^{sol}$ (${\text{J}\text{mol}}^{-1}$)	$E_{11}^{int}$ (${\text{J}\text{mol}}^{-1}$)	$\theta$ (K)
	$1\times {10}^{-5}$	$1.24\times {10}^4$	$-1.24\times {10}^4$	498
	${\kappa }_1$ (${\text{J}\text{m}}^2{\text{mol}}^{-1}$)	${\alpha }_1$ (${\text{J}\text{mol}}^{-1}$)	$M_1$ (${\text{m}}^5{\text{s}}^{-1}{\text{J}}^{-1}$)
	$2.5\times {10}^{-17}$	$5\times {10}^5$	$2.2\times {10}^{-19}$
Crystal plasticity model	${\dot{\gamma }}_0$ ($s^{-1}$)	n	$g_0$ (MPa)	$g_{\infty }$ (MPa)	a
	$1\times {10}^{-3}$	20	31	63	2.25
	$h_0$ (MPa)	Coplanar $h_{\alpha \beta }$ (MPa)	Non-coplanar $h_{\alpha \beta }$ (MPa)
	75	1	1.4
Elastic constants	$C_{11}\left(Pa\right)$	$C_{12}\left(Pa\right)$	$C_{44}\left(Pa\right)$
	$1.06\times {10}^{11}$	$6.0\times {10}^{10}$	$2.8\times {10}^{10}$

Table 2. Table 2: Chemical free energy parameters for the ternary spinodal system. Ω Ω \Omega is the molar volume; E A s ​ o ​ l superscript subscript 𝐸 𝐴 𝑠 𝑜 𝑙 E_{A}^{sol} and E B s ​ o ​ l superscript subscript 𝐸 𝐵 𝑠 𝑜 𝑙 E_{B}^{sol} are the solution energies, κ A subscript 𝜅 𝐴 \kappa_{A} and κ B subscript 𝜅 𝐵 \kappa_{B} are the gradient coefficients, α A subscript 𝛼 𝐴 \alpha_{A} and α B subscript 𝛼 𝐵 \alpha_{B} are the penalty parameters, and M A subscript 𝑀 𝐴 {M}_{A} and M B subscript 𝑀 𝐵 {M}_{B} are the mobilities for component A 𝐴 A and B 𝐵 B respectively; E A ​ A i ​ n ​ t superscript subscript 𝐸 𝐴 𝐴 𝑖 𝑛 𝑡 E_{AA}^{int} and E A ​ B i ​ n ​ t superscript subscript 𝐸 𝐴 𝐵 𝑖 𝑛 𝑡 E_{AB}^{int} are the A − A 𝐴 𝐴 A-A and A − B 𝐴 𝐵 A-B interaction energies; and θ 𝜃 \theta is the temperature.

Parameter	Value	Parameter	Value
$\Omega$	$1\times {10}^{-5}$ (${\text{m}}^3{\text{mol}}^{-1}$)	${\kappa }_A$	$1\times {10}^{-18}$ (${\text{J}\text{m}}^2{\text{mol}}^{-1}$)
$E_A^{sol}$	$1.24\times {10}^4$ (${\text{J}\text{mol}}^{-1}$)	${\kappa }_B$	$1\times {10}^{-18}$ (${\text{J}\text{m}}^2{\text{mol}}^{-1}$)
$E_B^{sol}$	$1.24\times {10}^4$ (${\text{J}\text{mol}}^{-1}$)	${\alpha }_A$	$1\times {10}^5$ (${\text{J}\text{mol}}^{-1}$)
$E_{AA}^{int}$	$-1.24\times {10}^4$ (${\text{J}\text{mol}}^{-1}$)	${\alpha }_B$	$1\times {10}^5$ (${\text{J}\text{mol}}^{-1}$)
$E_{AB}^{int}$	$-1.24\times {10}^4$ (${\text{J}\text{mol}}^{-1}$)	$M_A$	$2.2\times {10}^{-19}$ (${\text{m}}^5{\text{s}}^{-1}{\text{J}}^{-1}$)
$\theta$	498(K)	$M_B$	$2.2\times {10}^{-19}$ (${\text{m}}^5{\text{s}}^{-1}{\text{J}}^{-1}$)

Equations

$${\dot{c}}_m=-\mathrm{Div}{\mathbf{j}}_m+{\sigma }_m,\text{for}m=1,\dots M,$$

$$1=\sum _{m=1}^M{c}_m,0=\sum _{m=1}^M{\mathbf{j}}_m,0=\sum _{m=1}^M{\sigma }_m,$$

$$\mathbf{F}={\mathbf{F}}_{\text{e}}{\mathbf{F}}_{\text{c}}{\mathbf{F}}_{\text{p}},$$

$$0=\mathrm{Div}\mathbf{P},{\mathbf{P}}^{\text{T}}\mathbf{F}=\mathbf{P}{\mathbf{F}}^{\text{T}},\dot{\epsilon }=\mathbf{P}\cdot \dot{\mathbf{F}},$$

$$\dot{\eta }={\theta }^{-1}\delta -{\theta }^{-1}\mathrm{Div}\sum _{m=1}^M{\mu }_m{\mathbf{j}}_m+{\theta }^{-1}\sum _{m=1}^M{\mu }_m{\sigma }_m,$$

$$\delta =\mathbf{P}\cdot \dot{\mathbf{F}}+\sum _{m=1}^{M-1}{\tilde{\mu }}_m{\dot{c}}_m-\dot{\psi }-\sum _{m=1}^{M-1}{\mathbf{j}}_m\cdot \nabla {\tilde{\mu }}_m$$

$$\psi \left(\mathbf{F},{\mathbf{F}}_{\text{p}},{\mathbf{F}}_{\text{c}},c_m,{\tilde{c}}_m,\nabla {\tilde{c}}_m\right)={\psi }_{\mathrm{e}}\left(\mathbf{F},{\mathbf{F}}_{\text{p}},{\mathbf{F}}_{\text{c}}\right)+{\psi }_{\mathrm{c}}\left(c_m,{\tilde{c}}_m,\nabla {\tilde{c}}_m\right),$$

$${\tilde{\mu }}_m={\partial }_{c_m}\psi +{\partial }_{{\mathbf{F}}_{\text{c}}}\psi \cdot {\partial }_{c_m}{\mathbf{F}}_{\text{c}},\mathbf{P}={\partial }_{\mathbf{F}}\psi .$$

$$\delta =-{\partial }_{{\mathbf{F}}_{\text{p}}}\psi \cdot {\dot{\mathbf{F}}}_{\text{p}}-\sum _{m=1}^{M-1}{\mathbf{j}}_m\cdot \nabla {\tilde{\mu }}_m$$

$${\psi }_{\mathrm{e}}(\mathbf{F},{\mathbf{F}}_{\text{p}},{\mathbf{F}}_{\text{c}})=\frac{1}{2}{\mathbf{E}}_{\text{e}}\cdot \mathbb{C}{\mathbf{E}}_{\text{e}},$$

$${\mathbf{E}}_{\text{e}}=\frac{1}{2}{{\mathbf{F}}_{\text{c}}}^{\text{T}}\left({{\mathbf{F}}_{\text{e}}}^{\text{T}}{\mathbf{F}}_{\text{e}}-\mathbf{I}\right){\mathbf{F}}_{\text{c}},$$

$$\mathbf{S}=\mathbb{C}{\mathbf{E}}_{\text{e}}.$$

$$\mathbf{P}={\mathbf{F}}_{\text{e}}{\mathbf{F}}_{\text{c}}\mathbf{S}{{\mathbf{F}}_{\text{p}}}^{\text{-T}}.$$

$${\dot{\mathbf{F}}}_{\text{p}}={\mathbf{L}}_{\text{p}}{\mathbf{F}}_{\text{p}},$$

$${\mathbf{M}}_{\text{p}}=-{\partial }_{{\mathbf{F}}_{\text{p}}}\psi {{\mathbf{F}}_{\text{p}}}^{\text{T}}={({\mathbf{F}}_{\text{e}}{\mathbf{F}}_{\text{c}})}^{\text{T}}{\mathbf{F}}_{\text{e}}{\mathbf{F}}_{\text{c}}\mathbf{S}\approx {{\mathbf{F}}_{\text{c}}}^{\text{T}}{\mathbf{F}}_{\text{c}}\mathbf{S},$$

$${\mathbf{L}}_{\text{p}}={\dot{\mathbf{F}}}_{\text{p}}{{\mathbf{F}}_{\text{p}}}^{-1}=\sum _{\alpha }{\dot{\gamma }}^{\alpha }{\mathbf{s}}^{\alpha }\otimes {\mathbf{n}}^{\alpha },$$

$${\dot{\gamma }}^{\alpha }={\dot{\gamma }}_0{\left|\frac{{\tau }^{\alpha }}{g^{\alpha }}\right|}^n\mathrm{sgn}\left({\tau }^{\alpha }\right),$$

$${\dot{g}}^{\alpha }$$

$$-{\partial }_{{\mathbf{F}}_{\text{p}}}\psi \cdot {\dot{\mathbf{F}}}_{\text{p}}=\sum _{\alpha }{\tau }^{\alpha }{\dot{\gamma }}^{\alpha },\text{where}{\tau }^{\alpha }={\mathbf{M}}_{\text{p}}\cdot ({\mathbf{s}}^{\alpha }\otimes {\mathbf{n}}^{\alpha }).$$

$${\mathbf{F}}_{\text{c}}=\sum _{m=1}^{M-1}\left(1+{\nu }_m{c}_m\right)\mathbf{I},$$

$$\Omega {\psi }_{\mathrm{c}}(c_m,{\tilde{c}}_m,\nabla {\tilde{c}}_m)=$$

$$+\frac{1}{2}\sum _{m=1}^{M-1}{\alpha }_m{\left(c_m-{\tilde{c}}_m\right)}^2+\frac{1}{2}\sum _{m=1}^{M-1}{\kappa }_m|\nabla {\tilde{c}}_m{|}^2,$$

$$0$$

$$={\alpha }_m\left({\tilde{c}}_m-c_m\right)-\mathrm{Div}{\kappa }_m\nabla {\tilde{c}}_m+{\alpha }_m\left(c_m-{\tilde{c}}_m\right).$$

$${\mathbf{j}}_m=-M_m\nabla {\tilde{\mu }}_m,$$

$$\Omega {\tilde{\mu }}_m=E_m^{sol}+\sum _{n=1}^{M-1}{E}_{mn}^{int}{c}_n+R\theta \ln \frac{c_m}{c_M}+{\alpha }_m\left(c_m-{\tilde{c}}_m\right)-\Omega {\nu }_m\left({{\mathbf{F}}_{\text{c}}}^{\text{T}}\mathbf{S}\right)\cdot \mathbf{I}.$$

$$d=\frac{1}{n+1}\sum _{\alpha }{g}^{\alpha }{\dot{\gamma }}_0{\left|\frac{{\tau }^{\alpha }}{g^{\alpha }}\right|}^{n+1}+\frac{1}{2}\sum _{m=1}^{M-1}\nabla {\tilde{\mu }}_m\cdot M_m\nabla {\tilde{\mu }}_m,$$

$$\zeta :=\sum _{m=1}^{M-1}\left({\tilde{\mu }}_m{\dot{c}}_m+{\partial }_{{\tilde{c}}_m}\psi {\dot{\tilde{c}}}_m\right)+{\partial }_{\nabla \mathbf{\chi }}\psi \cdot \nabla \dot{\mathbf{\chi }}-\sum _{\alpha }{\tau }_{\alpha }{\dot{\gamma }}_{\alpha },p_{\mathrm{s}}=\sum _{m=1}^{M-1}{\tilde{\mu }}_m{\sigma }_m.$$

$$P={\int }_{{\mathcal{B}}_0}{p}_{\mathrm{v}}\text{d}\mathbf{x}+{\int }_{\partial {\mathcal{B}}_0^{\mathbf{f}}}{p}_{\mathrm{f}}\text{d}\mathbf{s},\text{where}{p}_{\mathrm{v}}=\zeta +d+p_{\mathrm{s}}.$$

$$\delta P=$$

$$+$$

Full text

††journal: Computer Methods in Applied Mechanics and Engineering

1 Introduction

The main goal of modern materials science is the theory-guided tailoring of materials, including chemical composition and microstructure, in order to obtain improved properties for a sustainable technological development. While there has been a tremendous growth over recent years in the use of modelling and simulation tools towards this goal (Shanthraj and Zikry, 2012; Wu et al., 2014; Liu et al., 2018; Roters et al., 2019), the realistic prediction of the thermo-chemo-mechanical interactions relevant to industrial processes is still a key development required to enable technological advances in material design, manufacturing and product development for harsh-service environments. Among several full-field simulation approaches, the phase-field method is particularly well-suited to model interface kinetics (Boettinger et al., 2002; Vaithyanathan and Chen, 2002; chen2002phase; Moelans et al., 2008; Emmerich, 2008; steinbach2009phase). It has been successfully applied to describe many thermo-chemo-mechanical processes including solidification (Nestler and Wheeler, 2000), precipitation (Zhou et al., 2010; Jokisaari et al., 2017; Rudraraju2016), fracture (spatschek2011phase; Schneider et al., 2017; Shanthraj et al., 2016, 2017) and dislocation motion (wang2001nanoscale; shen2003phase; Mianroodi et al., 2019). Further methodological developments including chemo-mechanical interface modelling and homogenisation have been treated recently in Svendsen2018.

The use of diffuse-interface models to describe interfacial phenomena dates back to Cahn and Hilliard (1958). The original Cahn-Hilliard (CH) equation was used to model spinodal decomposition in binary alloys, but has since been extended to multi-component systems and coupled with microelasticity (hu2001phase; fischer2003kinetics; garcke2005cahn; steinbach2006multi; Moelans et al., 2008). A critical challenge in simulating the thermodynamics of multi-component chemo-mechanical systems is the numerical approximation of a generally non-convex chemical free energy. Numerical methods have been developed to solve such systems using finite difference (Furihata, 2001), mixed finite element with composition and chemical potential treated as independent fields (Barrett et al., 1999; Gomez and Hughes, 2011), isogeometric (Gomez et al., 2008) and spectral (Zhu et al., 1999) spatial discretisations. In the context of numerical time-integration, the stability, robustness and efficiency of the resulting solution algorithm is sensitive to the non-convexity of the chemical free energy. A successful approach in this regard is the splitting of the chemical free energy into concave and convex components, which are then approximated separately (Elliott and Stuart, 1993; Eyre, 1998). In particular, splitting methods have recently been applied to the multi-component, multi-phase CH equation with a generalized non-convex Landau energy (Boyer and Minjeaud, 2011; Tavakoli, 2016). In the context of a mixed finite-element-based weak formulation of the classic CH relation, Gomez and Hughes (2011) employed a chemical energy splitting approach based on the sign of the fourth-order derivatives of the underlying energy functional. A recent mixed weak formulation of chemo-mechanics for two-phase, two-component finite-deformation gradient elastic solids based on unconditionally stable, second-order accurate time-integration and Taylor expansion of the non-convex energy has been given by Sagiyama et al. (2016).

More recently, starting from a grand potential functional (Plapp, 2011; Choudhury and Nestler, 2012), the transport relations have been formulated for multi-component and multi-phase systems (Aagesen2018). In such an approach, the thermodynamics of the system is set in a grand canonical ensemble and formulated in terms of the chemical potential. As a result, boundary and interface conditions based on the chemical potential (Kim et al., 1999), needed for a realistic representation of macroscopic systems where the total number of particles scannot be fixed, are naturally incorporated. However, these approaches are based on a Legendre transformation, which does not exist for more general non-convex forms of the chemical free energy. In this work, a convex splitting of the chemical free energy is instead used to invert the multi-component chemical potential–composition relation. This results in an expression for component transport in terms of the chemical potential rather than the chemical composition, analogous to the grand potential approach, but applicable to more general forms of the chemical free energy. This paper is organised as follows: the basic model formulation for a single phase elasto-viscoplastic solid is presented in Section 2, followed by an outline of its numerical implementation in Section 3. In particular, since a dissipation potential exists for the flux-force relations of the current model, the corresponding initial boundary-value problems (IBVP) can be formulated with the help of rate variational methods (Svendsen, 2004; Miehe, 2011, 2014). In Section 4 representative examples are used to benchmark and compare the proposed method with conventional transport relations. A summary is provided in Section 5 along with perspectives for future applications.

2 Model formulation

2.1 Basic relations

Let $\mathcal{B}_{0}\subset\mathbb{R}^{3}$ be a microstructural domain of interest, with boundary $\partial\mathcal{B}_{0}$. Attention is restricted in this work to the case of a single-phase elasto-viscoplastic solid mixture of $M$ diffusing chemical components. Component diffusion and reaction in the mixture is modeled as usual by the corresponding component mass or number balance relations

[TABLE]

in the mixture in terms of component (mass or number) concentration $0\leq c_{m}\leq 1$, the corresponding flux density $\mathbf{j}_{m}$, and the corresponding specific supply-rate, $\sigma_{m}$. Assuming the mixture is closed with respect to constituent mass or number, the corresponding sum relations

[TABLE]

and (1) imply $\sum_{m=1}^{M}\dot{c}_{m}=0$.

Besides this, the deformation resulting from an applied loading defines a field $\bm{\chi}(\mathbf{x}):\mathbf{x}\in\mathcal{B}_{0}\to\mathbf{y}\in\mathcal{B}$ mapping points $\mathbf{x}$ in the undeformed configuration $\mathcal{B}_{0}$ to points $\mathbf{y}$ in the deformed configuration $\mathcal{B}$. Following Shanthraj et al. (2017), the total deformation gradient, given by $\mathbf{F}=\partial\bm{\chi}/\partial\mathbf{x}=\nabla{}\bm{\chi}$, is multiplicatively decomposed as

[TABLE]

where $\mathbf{F}_{\text{p}}$ is a lattice preserving isochoric mapping due to plastic deformation, $\mathbf{F}_{\text{c}}$ represents the local deformation due to solute misfit, and $\mathbf{F}_{\text{e}}$ is a mapping to the deformed configuration. In the current approach the stress relaxation due to the component diffusion process is captured through the stress-free deformation gradient, $\mathbf{F}_{\text{c}}$.

Restricting attention to isothermal and quasi-static processes with no external supplies of momentum or energy, the balance relations for linear momentum, angular momentum, and internal energy, are given by

[TABLE]

where, $\mathbf{P}$ is the first Piola-Kirchhoff stress, and $\varepsilon$ is the referential internal energy density. As well the balance relation,

[TABLE]

for the entropy holds (e.g., de Groot and Mazur, 1962, Chapter 2), with dissipation-rate density, $\delta$, absolute temperature, $\theta$, and the chemical potential, $\mu_{m}$, of component $m$.

Combination of (1)-(5) yields the form

[TABLE]

for the dissipation-rate density in terms of $\tilde{\mu}_{a}:=\mu_{a}-\mu_{M}$, and the free energy density of the mixture, $\psi:=\varepsilon-\theta\eta$.

The current model formulation is based on the basic constitutive form

[TABLE]

for $\psi$ in terms of elastic, $\psi_{\mathrm{e}}$, and chemical, $\psi_{\mathrm{c}}$, parts. A non-local field, $\tilde{c}_{m}$ for independent components, $m=1,\ldots M-1$, is introduced to weakly enforce any dependence of the energy on chemical inhomogeneity through its gradient $\nabla{}\tilde{c}_{m}$ (Ubachs et al., 2004). Modelling $\mathbf{P}$ and $\tilde{\mu}_{m}$ as purely energetic, i.e. work conjugate to $\mathbf{F}$ and $c_{m}$ respectively, results in the constitutive relations

[TABLE]

Together, (6)-(8) result in the residual form

[TABLE]

for the dissipation-rate density.

2.2 Energetic constitutive relations

The elastic energy, $\psi_{\mathrm{e}}$, is modelled here relative to the intermediate configuration by the form

[TABLE]

in terms of the Green-Lagrange strain measure

[TABLE]

and anisotropic elastic stiffness $\mathbb{C}$. The work conjugate second Piola-Kirchhoff stress is given by,

[TABLE]

The 1st Piola-Kirchhoff stress tensor, $\mathbf{P}$, is then related to $\mathbf{S}$ through

[TABLE]

Further details are provided in Roters et al. (2019).

2.2.1 Crystal plasticity

The plastic deformation gradient is given in terms of the plastic velocity gradient, $\mathbf{L}_{\text{p}}$, by the flow rule

[TABLE]

where $\mathbf{L}_{\text{p}}$ is work conjugate with the Mandel stress in the plastic configuration,

[TABLE]

assuming small elastic strains. A crystal plasticity model is used, where the plastic velocity gradient $\mathbf{L}_{\text{p}}$ is composed of the slip rates $\dot{\gamma}^{\alpha}$ on crystallographic slip systems, which are indexed by $\alpha$

[TABLE]

where $\mathbf{s}^{\alpha}$ and $\mathbf{n}^{\alpha}$ are unit vectors along the slip direction and slip plane normal, respectively (Roters et al., 2010; Shanthraj and Zikry, 2011). The slip rates are given by the phenomenological description of Peirce et al. (1983),

[TABLE]

in terms of the reference shear rate $\dot{\gamma}_{0}$, and strain-rate sensitivity exponent $n$. The slip resistances on each slip system, $g^{\alpha}$, evolve asymptotically towards $g_{\infty}$ with shear $\gamma^{\beta}$ ($\beta=1,\ldots,12$) according to the relationship

[TABLE]

with parameters $h_{0}$ and $a$. The interaction between different slip systems is captured by the hardening matrix $h_{\alpha\beta}$.

The plastic dissipation can then be reduced to the simpler slip system based conjugate pair

[TABLE]

On this basis, non-negativity of $\dot{\gamma}_{0}$, $g_{\infty}$, $h_{0}$ and $h_{\alpha\beta}$ is sufficient to ensure non-negative plastic dissipation.

2.2.2 Multi-component chemo-mechanics

Following Hüter et al. (2018), the local deformation arising from volumetric mismatch between solute components is given by

[TABLE]

where the volumetric mismatch coefficient, $\nu_{m}$, is related to the Vegard’s coefficient, i.e. $\nu_{m}=V_{m}/a_{M}$, where $V_{m}$ is the Vegard’s coefficient for component $m$ dissolved as a solid solution in a matrix of component $M$, and $a_{M}$ is the lattice parameter of the matrix.

The chemical energy, $\psi_{\mathrm{c}}$, is modelled here using the following regular solution form

[TABLE]

where $\Omega$ is the molar volume; $E_{m}^{sol}$, $\kappa_{m}$ and $\alpha_{m}$ are respectively the solution energy, gradient coefficient, and penalty parameter to weakly enforce $c_{m}=\tilde{c}_{m}$, for component $m$; $E_{mn}^{int}$ is the interaction energy between component $m$ and $n$; and $R$ is the universal gas constant. The non-local field, $\tilde{c}_{m}$, is obtained through the equilibrium relation,

[TABLE]

Following Onsager1931, a linear flux-force form is assumed for the component diffusion,

[TABLE]

with $M_{m}$ the component mobility in the lattice-fixed frame of reference, and chemical potential, $\tilde{\mu}_{m}$, given by

[TABLE]

On this basis, non-negativity of $M_{m}$ is sufficient to ensure non-negative dissipation due to component diffusion.

3 Numerical methods

The constitutive model introduced in Section 2 is implemented in the freeware material simulation kit, DAMASK (Roters et al., 2019), and a large-scale parallel finite element (FE) code using the PETSc numerical library (Balay2015) is developed to handle the discretisation and numerical solution of the coupled field equations.

3.1 Rate variational formulation of the initial boundary-value problems

The following represents a special case of the variational treatment of combined Cahn-Hilliard and Allen-Cahn modelling of finite-deformation gradient elastic solids in Gladkov and Svendsen (2015).

As detailed in Svendsen (2004), a rate-variational-based formulation of the IBVPs is contingent in particular on the existence of a dissipation potential for the model in question. In the current case, the force-based form

[TABLE]

of this potential, $d$, determines the fluxes $\dot{\gamma}^{\alpha}=\partial_{\tau_{\alpha}}d$ and $\mathbf{j}_{m}=-\partial_{\nabla{}\tilde{\mu}_{m}}d$, consistent with (17) and (23), respectively. The convexity of $d$ in the forces, and its non-negativity $d\geqslant 0$, together imply $\delta\geqslant 0$ in the context of (9). Besides this potential, the rate-variational formulation is based on the energy storage-rate density, $\zeta$, and the supply rate density, $p_{\mathrm{s}}$. For the current model, this is given by

[TABLE]

Together, these determine the volumetric part of the rate functional

[TABLE]

The flux boundary conditions on $\partial\mathcal{B}_{0}^{\mathbf{f}}$ is given by $p_{\mathrm{f}}$. The first variation of $P$ with respect to $\tilde{\mu}_{m}$, $\dot{\bm{\chi}}$, $\dot{\tilde{c}}_{m}$, and $\tau_{\alpha}$ takes the form

[TABLE]

via integration by parts and the divergence theorem. Necessary for stationarity, $\delta P=0$, of (28) are the weak forms for the field relations

[TABLE]

and constitutive relation $0=-\dot{\gamma}^{\alpha}+\partial_{\tau_{\alpha}}d$.

3.2 Finite element implementation

The rate-variational formulation yields the weak form of the field relations required for their finite-element (FE) implementation. In particular, equation 29 yields directly the weak momentum balance relation

[TABLE]

where $p_{\mathrm{f}}=0$ is assumed for simplicity, the weak multi-component transport relation

[TABLE]

and the weak non-local relation

[TABLE]

where $\delta\dot{}\bm{\chi}$, $\delta\tilde{\mu}_{m}$ and $\delta\dot{\tilde{c}}_{m}$ are the virtual deformation rate, chemical potential and non-local concentration fields respectively. No-flux boundary conditions are assumed for the sake of simplicity.

The deformation field, $\bm{\chi}\,(\mathbf{x})$, chemical potential, $\tilde{\mu}_{m}\,(\mathbf{x})$, and non-local concentration field, $\tilde{c}_{m}\,(\mathbf{x})$, in addition to their virtual counterparts are discretised using a FE basis of shape functions, $N_{i}^{\bm{\chi}}$, $N_{i}^{\mu}$, $N_{i}^{\tilde{c}}$, $N_{i}^{\delta\dot{}\bm{\chi}}$, $N_{i}^{\delta\mu}$, and $N_{i}^{\delta\tilde{c}}$, where $[\bm{\chi}]_{i}$, $[\tilde{\mu}_{m}]_{i}$, $[\tilde{c}_{m}]_{i}$, ${[\delta\dot{}\bm{\chi}]_{i}}$, ${[\delta\tilde{\mu}_{m}]_{i}}$, and $[\delta\tilde{c}_{m}]_{i}$ are the respective degrees of freedom. The corresponding discrete differential operator matrices are $\mathbf{B}_{i}^{\bm{\chi}}$, $\mathbf{B}_{i}^{\mu}$, $\mathbf{B}_{i}^{\tilde{c}}$, $\mathbf{B}_{i}^{\delta\dot{}\bm{\chi}}$, $\mathbf{B}_{i}^{\delta\mu}$ and $\mathbf{B}_{i}^{\delta\tilde{c}}$. Under these approximations, the weak forms equations 30, 31 and 32 can be rewritten as

[TABLE]

[TABLE]

[TABLE]

which defines a non-linear system of equations for the unknowns $[\bm{\chi}]_{i}$, $[\tilde{\mu}_{m}]_{i}$, and $[\tilde{c}_{m}]_{i}$. A time-discrete system of equations is obtained by using a backward Euler approximation

[TABLE]

of the rate $\dot{c}_{m}$ in equation 34.

The solution approach followed in this work involves solving the coupled system of equations 33, 34 and 35 within a staggered iterative loop until a self consistent solution is achieved for a time increment.

3.3 Chemical potential solution

The solution of the coupled system of equations 33, 34 and 35, for the deformation field, $\bm{\chi}$, chemical potential, $\tilde{\mu}_{m}$, and non-local concentration field, $\tilde{c}_{m}$, requires the inversion of equation 24 in order to express $c_{m}:=c_{m}({\tilde{\mu}_{n}})$ for $m,n=1,\ldots,M-1$. This is achieved algorithmically in the current work through a semi-implicit splitting of the chemical potential relation,

[TABLE]

into a convex

[TABLE]

and non-convex, i.e.,

[TABLE]

contribution. Equation 38 can then be inverted to express $c_{m}(t_{n})$ in terms of equations 37 and 39, yielding,

[TABLE]

where

[TABLE]

Equation 40 is an implicit system of equations to be numerically solved for $c_{m}(t_{n})$ for a given set ${\tilde{\mu}_{n}}(t_{n})$. A fixed point iteration is employed, which is unconditionally convergent for $\alpha_{m}\geq 0$ as the the fixed point operator, i.e. the RHS in equation 40, is guaranteed to be a contractive mapping.

4 Results and discussion

In this section, the developed chemo-mechanical model for multi-component finite-strain elasto-viscoplastic materials is validated, benchmarked and showcased through illustrative examples.

4.1 Validation and benchmarking of the numerical scheme

Diffusion simulations are first performed to study the convergence, accuracy and performance of the proposed numerical scheme for the CH model.

4.1.1 Convergence behaviour

A non-local concentration field, $\tilde{c}$, is introduced in the current model to account for the gradient energy contributions, thereby reducing the strong fourth-order CH PDE to two weakly non-local second-order PDEs. The conditions which should be fulfilled for the proposed method to approach the results by the strong CH solution are discussed in this section. To investigate the effect of the penalty parameter, $\alpha$, and the gradient coefficient, $\kappa$, on the equilibrium numerical solution, we study one-dimensional (1-D) diffusion simulations without mechanical deformation. Parameters for the chemical free energy density function in equation 21 are taken as $\Omega=1\times 10^{-5}\text{m}^{3}\text{mol}^{-1},E_{1}^{sol}=1.24\times 10^{4}\text{J}\text{mol}^{-1},E_{11}^{int}=-1.24\times 10^{4}\text{J}\text{mol}^{-1},\theta=498\text{K},\kappa_{1}=1\times 10^{-16}\text{J}\text{m}^{2}\text{mol}^{-1},M_{1}=2.2\times 10^{-19}\text{m}^{5}\text{s}^{-1}\text{J}^{-1}$ and the penalty parameter $\alpha_{1}$ is varied from $4\times 10^{3}$ to $2.5\times 10^{6}\text{J}\text{mol}^{-1}$. The resulting spinodal compositions are $c_{\alpha}=0.07$ and $c_{\beta}=0.93$, respectively. The chemical free energy without mechanical deformation of the studied system is represented by the black curve in figure 1. The minima in the double-well free energy curve represent the spinodal compositions resulting from the decomposition of an initial homogeneous mixture.

The 1-D domain, $\mathcal{B}_{0}=[0,2L]$, is initially into two regions, $\mathcal{B}_{\alpha}=[0,L)$ and $\mathcal{B}_{\beta}=(L,1.0]$. The initial concentrations for these two regions are taken as $c_{\alpha}=0.49$ and $c_{\beta}=0.51$, respectively. The system size, $2L$, is $3\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and meshed with $300$ regular hexahedral elements.

Figure 2(a) shows the maximum of the pointwise difference between concentration and non-local concentration, $\text{max}_{x\in[0,2L]}|c(x)-\tilde{c}(x)|$, at steady state, with increasing values of the penalty parameter. The difference decreases substantially from $3.4\times 10^{-1}$ to $6.8\times 10^{-4}$ with increasing the penalty parameter from $4.0\times 10^{3}$ to $2.5\times 10^{6}\text{J}\text{mol}^{-1}$. Convergence of the concentration, $c$, to the non-local concentration, $\tilde{c}$, is observed with an exponential reduction in the difference on increasing the value of the penalty parameter, $\alpha$. Figure 2(b) presents the variation of the interface width, $d$, with increasing the penalty parameter. The interface width converges to a constant value of $4.3\times 10^{-7}m$ with minimal differences when the penalty parameter is larger than $1.0\times 10^{5}\text{J}\text{mol}^{-1}$. The convergence of the interface width in such studies is a good indicator for choosing an appropriate penalty parameter.

An alternative way to characterize the convergence of the proposed numerical algorithm is to study the relationship between the gradient coefficient, $\kappa_{m}$, and the interface width, $d$. Based on the classical CH theory (Cahn and Hilliard, 1958), the gradient coefficient should have quadratic relationship with the interface width. Figure 3 shows the simulated interface width squared as a function of the gradient coefficient at steady state. The gradient coefficient $\kappa_{1}$ varies from $1\times 10^{-17}$ to $1\times 10^{-15}\text{J}\text{m}^{2}\text{mol}^{-1}$. Note that for each simulation, the penalty parameter is chosen large enough to guarantee that the interface width converges to the constant value as discussed above. Figure 3 verifies the linear relationship between the interface width and the gradient coefficient, and shows that the penalty parameter does not affect the simulation results when it is sufficiently large. The above results indicate that the numerical solutions of the proposed approach converges to the conventional CH model for a sufficiently large penalty parameter, and guidelines for choosing the penalty parameter are provided.

4.1.2 Numerical performance

In this section, the accuracy and performance of the proposed numerical approach is compared to a conventional CH solution scheme. A two-dimensional spinodal decomposition example in the absence of mechanical deformation is solved to compare the two solution schemes. The model parameters used are as follows: $\Omega=1\times 10^{-5}\text{m}^{3}\text{mol}^{-1},E_{1}^{sol}=2\times 10^{4}\text{J}\text{mol}^{-1},E_{11}^{int}=-5\times 10^{4}\text{J}\text{mol}^{-1},R\theta=0.1\times 10^{4}\text{J}\text{mol}^{-1},\kappa_{1}=1\times 10^{-16}\text{J}\text{m}^{2}\text{mol}^{-1},M_{m}=1.0\times 10^{-19}\text{m}^{5}\text{s}^{-1}\text{J}^{-1}$, $\alpha_{1}=5\times 10^{5}\text{J}\text{mol}^{-1}$. A square domain, $\mathcal{B}_{0}=[0,L]\times[0,L]$, of size $L=$2.56\text{,}\mathrm{\SIUnitSymbolMicro m}$$ is meshed with $256\times 256$ regular hexahedral elements. The initial concentration is taken as 0.5 with a random perturbation of amplitude 0.1.

Figure 4 shows the temporal evolution of the concentration field during spinodal decomposition, starting from a homogeneous mixture with an average concentration of 0.5 and random perturbation with an amplitude of 0.1, into two spinodal compositions of $1.0\times 10^{-4}$ and $9.9\times 10^{-1}$, determined by the minima of the free energy, $\psi$. Negligible differences in the concentration fields between the two solution approaches are observed. In order to quantitatively compare the simulation results, figure 5 shows the evolution of the maximum and minimum values of the concentration during the spinodal decomposition process, with spinodal decomposition initiating at 500s. It can be seen that the simulated results for these two different approaches completely overlap during the entire spinodal decomposition process, thus validating the accuracy of the proposed numerical approach.

In order to compare the numerical performance of the two numerical approaches, the same Newton-Raphson scheme is used in the above simulations. The number of iterations required to obtain the solution for the nonlinear Newton-Raphson solver at each time step is used to quantify the performance and efficiency of the numerical scheme. The number of Newton iterations for both approaches during the calculations is shown in figure 6. The number of Newton iterations for solving the inverted transport relations (2-4 iterations) is significantly lower than that for solving the conventional transport relations (12-14 iterations). These comparisons demonstrate the performance and the efficiency of the proposed numerical scheme.

4.2 Chemo-mechanical coupling

In this section, the chemical diffusion model is coupled to a finite-strain mechanical model to study the role of mechanical deformation on the spinodal decomposition process. Both elastic and plastic deformation accompanying diffusion are considered in the simulations. The material parameters for the chemical energy and crystal plasticity model are listed in table 1. To illustrate the effect of mechanical deformation, through the deformation arising from volumetric mismatch between solute components, $\mathbf{F}_{\text{c}}$ in equation 20, on the free energy of the system, figure 1 shows the free energy of a 1-D system with various values of the solute volumetric mismatch coefficient. Note that in this plot, it is assumed that the local deformation arising from volumetric mismatch between solute components can only be accommodated by elastic deformation. It can be seen that in terms of kinetics, the driving force for spinodal decomposition decreases with increasing levels of volumetric mismatch. From an energetic point of view, the miscibility gap is also modified, with the solute-poor spinodal being more energetically favoured.

A square domain, $\mathcal{B}_{0}=[0,L]\times[0,L]$, of size $L=$2.56\text{,}\mathrm{\SIUnitSymbolMicro m}$$ is meshed with $256\times 256$ regular hexahedral elements. Periodic boundary conditions on both concentrations and displacements were applied. The temporal evolution of the concentration field during the spinodal decomposition process, in the absence of mechanical deformation (i.e. setting $\nu_{1}=0.00$), is shown in figure 7 (a). At the early stage of the spinodal decomposition, the initial homogeneous mixture decomposes into two regions with different compositions, following their corresponding composition branch on the generalized muti-well Landau energy landscape figure 7 ($\text{a}_{\text{1}}$). The decomposition stage is driven by the minimization of the chemical free energy. Following decomposition, coarsening can be observed from figure 7 ($\text{a}_{\text{2}}$) to ($\text{a}_{\text{4}}$), which is driven by the minimization of interface energy, that occurs at a longer time scale compared to the initial decomposition stage.

The effect of the volumetric mismatch between solute components on spinodal decomposition and coarsening is investigated by first considering the case of an elastically deforming material. The volumetric mismatch coefficient, $\nu_{1}$, is varied between $0.01\ \text{and}\ 0.05$, and the material parameters used are listed in table 1. The temporal evolution of the concentration field during spinodal decomposition is shown in figure 7 (b) to (d). It can be seen that the decomposition is minimally affected when the solute induced deformation is relatively small ($\textit{i.e.},\nu_{1}=0.01$), comparing figure 7 (a) and (b). However, as the volumetric mismatch coefficient increases, the spinodal decomposition kinetics is significantly reduced, as shown in figure 7 (c) and (d). This is due to the increasing elastic energy contribution to the driving force, which has the effect of suppressing the spinodal decomposition (figure 1). The spinodal compositions are also affected by the chemo-mechanical coupling, with the solute-rich spinodal point decreasing from 0.93 to 0.82 when the volumetric mismatch coefficient increases up to 0.05. Furthermore, one can observe that the morphology of the solute-rich region is also affected by the deformation arising from volumetric mismatch between solute components, comparing figure 7 (a) and (d). In the absence of mechanical coupling, the solute-rich regions exhibit a spherical morphology due to the isotropic nature of the interface energy (figure 7 ($\text{a}_{\text{4}}$)), but with increasing volumetric mismatch between solute components, the solute-rich regions are observed to order themselves along the softer directions of the cubically anisotropic elastic stiffness tensor used here (figure 7 ($\text{d}_{\text{4}}$)).

To our knowledge, the numerical investigations of spinodal decomposition, presented in the literature, are limited to the coupling of elastic deformation with diffusion. The capability of the developed chemo-mechanical model allows us to explore the more challenging influence of elasto-plasticity on the spinodal decomposition and coarsening, which is observed in materials applications such as rafting in superalloys (Kontis et al., 2018) and hydride formation (Korbmacher et al., 2018). Parameters used for the crystal plasticity model are listed in table 1.

The temporal evolution of the concentration field during spinodal decomposition on an elastic-plastically deforming material is shown in figure 8. Comparing figure 7 and figure 8 reveals that visco-plastic relaxation significantly affects the decomposition process. Spinodal decomposition is observed at 6000s in the plastic material (figure 8 ($\text{c}_{\text{3}}$)) compared to 10000s in the elastic material (figure 7 $\text{a}_{\text{4}}$), for the largest volumetric mismatch coefficient considered. In addition, the spinodal compositions are close to those obtained from the simulations in the absence of mechanical deformation, since the deformation arising from volumetric mismatch between solute components can be effectively relaxed by the plastic deformation.

Figure 9 presents the evolution of the hydrostatic stress and plastic strain for the cases with elastic-only and elasto-plastic deformation at a volumetric mismatch coefficient of 0.03. The plastic deformation relaxes the hydrostatic stress generated due to solute agglomeration during spinodal decomposition, and thus dissipates a significant portion of the stored elastic energy. The plastic strain localization surrounding the solute-rich regions is observed, with a maximum value of 0.6. More interestingly, even though the spinodal compositions are similar to the case without mechanical deformation, the morphology of the solute-rich region differs significantly. When considering elasto-plastic deformation, the solute-rich regions exhibit an elliptic morphology at the intermediate level of volumetric misfit, as shown in figure 8 ($\text{b}_{\text{4}}$). At the high levels of volumetric misfit, the solute-rich regions merge together and multiple lamellar regions can be observed, as shown in figure 8 ($\text{c}_{\text{4}}$).

4.3 Ternary chemo-mechanical spinodal decomposition

In order to illustrate the applicability of the modelling approach to more complex material systems, a large three-dimensional (3D) simulation of a chemo-mechanically coupled ternary spinodal decomposition and coarsening process is considered next. The cubic domain, $\mathcal{B}_{0}=[0,L]\times[0,L]\times[0,L]$, of size $L=$1.28\text{,}\mathrm{\SIUnitSymbolMicro m}$$ is meshed with $128\times 128\times 128$ regular hexahedral elements. The initial concentrations for solute A and B are 0.3 and 0.3, respectively, with a random perturbation of amplitude 0.01. The chemical free energy parameters are listed in table 2, and the crystal plasticity material parameters are listed in table 1. The chemical free energy parameters used will result in the formation of three decomposition regimes: (A-rich, B-poor), (A-poor, B-rich), and (A-poor, B-poor).

Figure 10 shows the temporal evolution of the morphologies of different decomposition regimes, the accompanying hydrostatic stress, and the plastic strain. In this fully coupled 3D ternary spinodal-elastic-viscoplastic mechanical decomposition simulation, the separated regions exhibit inter-connected morphologies instead of the isolated islands like those in the 2D binary cases, as the three components have relatively the same concentrations (0.3, 0.3, 0.4), close to the off-critical point. The hydrostatic stress distribution is significantly heterogeneous, even inside the same decomposition region, due to the accompanying plastic deformation, as shown in figure 10(b). It can be seen that the magnitude of the hydrostatic stress only changes slightly during the phase coarsening stage, while the total plastic strain increases gradually (figure 10(c)).

The corresponding 3D simulations of the ternary spinodal decomposition in the absence of mechanical deformation and only considering elastic deformation are presented in figure 11. Comparing figure 10(a) and figure 11, the separated regions produced by these three simulations exhibit qualitatively similar inter-connected morphologies. The approach of (Gameiro et al., 2005) could be used to quantify and distinguish the geometrical differences between these complicated patterns, however, a full discussion is beyond the scope of the current work and represents work in progress to be reported in the future.

5 Conclusions and Outlook

In this work, a numerical approach has been introduced, using a semi-implicit convex splitting to invert the multi-component chemical potential relations, and resulting in a numerically advantageous expression for their transport in terms of the chemical potential. The approach is validated using a generalized spinodal Landau-type system as a benchmark, exposed to different chemo-mechanical starting and boundary conditions. A convergence of the proposed weakly non-local approach to the strong CH form is demonstrated with increasing values of the penalty parameter. However, for practical purposes, a penalty parameter of $1\times 10^{5}\text{J}\text{mol}^{-1}$ was found to give sufficiently converged results. A significant reduction of the numerical cost and stability of the proposed approach is also demonstrated, allowing for the use of larger time steps in long-term diffusion simulations.

The influence of chemo-mechanical coupling, including both elastic and anisotropic crystalline plastic deformation, on the spinodal decomposition and coarsening process is also investigated. It is found that the stored elastic energy can narrow the miscibility gap and alter the spinodal compositions. The decomposition kinetics is significantly reduced with increasing volumetric mismatch between solute components Furthermore, the morphology of the solute-rich phase changes from spherical shape distributions to an ordered cubic distribution due to the anisotropy of the elastic stiffness tensor. When taking in to account additional visco-plastic relaxation of the volumetric mismatch, as in the case of realistic materials applications, it is found that the plastic relaxation accounts for a significant amount of the dissipated total energy. The spinodal compositions are close to those in the absence of mechanical deformation, but the solute-rich regions tend to coalesce into elongated lamellar morphologies.

Extension of the proposed approach to multi-phase systems is straightforward. This will be interesting as, in the context of multiphase systems, working directly with the chemical potential lends itself naturally to the Kim-Kim-Suzuki (KKS) description of interfaces (Kim et al., 1999), since, algorithmically, no partitioning of an average solute composition into the individual phases is required. The proposed approach is also particularly amenable to CALPHAD-based descriptions of the chemical energy (Lukas et al., 2007), which is essential in the drive towards predictive simulations.

6 Acknowledgements

PS CL and BS gratefully acknowledge the financial support by the DFG through the Priority Programme SPP 1713: Strong Coupling of Thermo-chemical and Thermo-mechanical States in Applied Materials. PS is also grateful to the EPSRC for financial support through the associated programme grant LightFORM (EP/R001715/1) and the Airbus–University of Manchester Centre for Metallurgical Excellence for supporting aspects of this research.

7 References