An elasto-chemical phase-field model for isotropic solids

Accepted Manuscript An elasto-chemical phase-field model for isotropic solids O. Tschukin, D. Schneider, B. Nestler PII:

S0997-7538(17)30378-9

DOI:

10.1016/j.euromechsol.2018.06.014

Reference:

EJMSOL 3632

To appear in:

European Journal of Mechanics / A Solids

Received Date: 16 May 2017 Revised Date:

28 June 2018

Accepted Date: 29 June 2018

Please cite this article as: Tschukin, O., Schneider, D., Nestler, B., An elasto-chemical phasefield model for isotropic solids, European Journal of Mechanics / A Solids (2018), doi: 10.1016/ j.euromechsol.2018.06.014. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

An elasto-chemical phase-field model for isotropic solids O. Tschukina , D. Schneidera,b , B. Nestlera,b a IAM-CMS, Karlsruhe Institute of Technologies, Strasse am Forum 7, 76131 Karlsruhe, Germany Institute of Materials and Processes, Karlsruhe University of Applied Sciences, 76133 Karlsruhe, Germany

RI PT

b

Abstract

SC

A quantitative elasto-chemical phase-field model is presented. The challenge to avoid the falsification of the energetics in the diffuse interface (excess energy) is mastered by the usage of the appropriate thermodynamic potentials. Based on the mechanical jump conditions at the coherent interface between two elastically isotropic materials, the calculation of the required elastic fields is simplified to increase the computational efficiency. The model is verified using an analytical solution for the Eshelby problem.

4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33

TE D

3

The understanding of the mechanisms of microstructure evolution plays a key role for application and for the design of new materials. The presence and the location of the interfaces in the microstructure significantly influence the mechanical, thermal and other engineering properties of the used material. Herein the term ”interface” implies a border between two subdomains, which differ at least in one feature (property). For example, the interface can be the boundary between a crystal and its melt, the common interface between two identical misoriented grains, the fracture surface, etc. Since the nature of such interfaces is ambiguous, the reason for their traction roots on different principles. In order to avoid confusion, we differentiate two fundamentals for the motion of the interfaces. First, we speak about the transitional character, if a phase transformation takes place. The motion of the free boundary between the transitional phases is proportional to the resulting driving force as a jump of the appropriate bulk potentials. As a result the neighbouring phases grow and shrink. On the other hand, a different mechanism governing the notion of the interface is a ”dispositional” one. Here the boundary of the neighbouring ”phases” moves with respect to transposition or deformation. Traditionally, phase-field models have been employed for quantitative description of the solidification processes, e.g. for pure components [1, 2] and alloys [3–6]. Subsequently, it has also been applied for dispositional boundary motion, such as wetting [7], sintering [8], fracturing [9] and others. Regardless of the real process, which is described with the phase-field approach, one is endeavoured to homogenise the bulk properties throughout the diffuse interface. Hence, there are also different procedures for the incorporation of the elastic effects into diffuse interface approach. Thus, for example,

EP

2

1. Introduction

AC C

1

M AN U

Keywords: solid-solid phase transformation, phase-field model, coherent interface, thermomechanical equilibrium PACS: 64.70.K, 81,40Jj

Email address: [email protected] (O. Tschukin) Preprint submitted to European Journal of Mechanics - A/Solids

the Voigt/Taylor (VT) homogenisation scheme implies the direct interpolation of the phase corresponding strain energies and concludes the direct interpolation of stiffness [10, 11] and of eigenstresses [12]. An alternative model presented by the authors in [13] assumes the homogeneity of the stresses in the diffuse interface. As a consequence, the calculation of the effective stiffness tensor follows the harmonic type, which is known as Reuss/Sachs (RS) homogenisation approach [14], but the eigenstrains are directly interpolated inside the diffuse interface. Furthermore, the Khachaturyan’s model [15] is a combination of both mentioned approaches and is widely used in the phase field community, e.g. [15, 16]. Herein, the interfacial stiffness is calculated with respect to the VT homogenisation scheme, but instead the eigenstresses, the eigenstrains follow the arithmetic interpolation type in the diffuse interface. In [17], the authors propose a model by combining the VT and RS interpolation schemes to calculate the elastic driving forces and the stress tensor, in order to avoid additional excess energy in the diffuse interface. But the relevant mismatch between the analytical and numerical results reveals the incompleteness of the model. Using a diffuse interface approach with an artificially widened interface layer, the above mentioned inconsistent interpolations lead to spurious interface energetics, as shown in [17, 18]. Since the error scales with the diffuse interface width W, the excess energy can be controlled by selecting the computational interface width, much smaller as the local curvature κ, so that κW  1. Such a severe restriction on the diffuse interface limits the numerical simulations of the real manufacturing processes, which belies the main goal of the usage of the phase field model. In principle, one should be able to artificially widen the diffuse interface width, while still capturing the real physics. Thus, for the quantitative numerical results but in the presence of the distorted interfacial physics, the use of the simulation parameters for a special physical setup is cumbersome. In our recent work [18], we also investigate the reasons for June 30, 2018

34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69

ACCEPTED MANUSCRIPT

77 78 79 80 81 82 83 84 85 86 87

88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105

106 107 108 109 110 111 112 113

114 115 116 117 118 119 120 121 122

RI PT

76

SC

75

Inspired by the philosophy of the grand chemical potential models in [3–5] we propose in [18] an alternative homogenisation approach for the elastic constants and for the eigenstrain/eigenstress, purely based on the mechanical jump conditions acting on the coherent interface. Since the strain energies are functions of strain tensor the interpolation of elastic energy throughout diffuse interface should be done with phase depending strain energies, evaluated with phase corresponding strains. Consequently, the homogenised strain energy in the diffuse interface should be consistent with the homogeneous variables and to the jumps in the inhomogeneous variables. The validation of the derived model is performed for 1-D test cases as well as for a circular inclusion, [18]. In the analysed scenarios we get an excellent agreement between the theoretical sharp interface predictions and simulation results. Although, the central aim is the formulation of the consistent potential for the interpolation in the diffuse interface and the consistent calculation procedure of the stress tensor.

M AN U

74

123

The evolution equations of the system variables are originally derived by the variational approach of the Lyapunov functional as energetics (entropy) of the system, where the vector of the state variables s incorporates the intensive or extensive thermodynamic quantities such as temperature T , composition c, chemical potential µ, deformation etc. The temporal and spatial change of the total energy in the system E(s) consist of the evolution of the bulk energies, and especially on the evolution of interfacial energies between different phases, according to the morphological change. Herein, the interfacial energy between the neighbouring phases is defined on the real phase boundary with the physical width of some angstroms. The adjustment of the length scales of the bulk and interfacial physics, for the practicable computational performance, is inter alia achieved by the introduction of the phase field φ ∈ G = PN {(φ1 , ...., φN ) : φi ∈ [0, 1], i=1 φi = 1}, which characterises the different bulk phases with uniform material properties. The objective Lyapunov functional takes the form of a GinzburgLandau type [4, 19–21]

E(φ, ∇φ, s) =

TE D

73

EP

72

2. Elasto-chemical phase-field model

the interfacial excess energies in the mentioned models by the underlying different homogenisation assumption and the resulting interpolation schemes of the material constants in the diffuse interface. Furthermore we show, that the elastic models only partially satisfy the mechanical jump conditions at the interface. Although the homogeneous normal stress and tangential strain are consistent with the vanishing traction force and with the Hadamard jump condition, respectively. In contrast to the homogeneous variables, tangential stress and normal strain are inhomogeneous variables with a discontinuity at the common interface. The previous consideration with homogeneous and inhomogeneous strains and stresses is similar to the condition of the same chemical potential for the solid and liquid phases in alloys. Since, at the common interface, the compositions of two neighbouring phases differ from each other, in the quantitative phase-field modelling, the free energy are given as an interpolation of phase corresponding free energies, each evaluated with self-contained concentration, [3].

Z

δE(φ, ∇φ, s) , δφα

K X i=1

The paper is organised as follows. In sec. 2 we present the extended phase-field model, which consists of capillary, chemical and elastic terms. By verifying the model we use the setup, for which the theoretical determination of all required quantities is given in a closed-form solution, sec. 3 and explain the corresponding numerical setup in sec. 4. In sec. 5 we compare the simulation results with analytical prescriptions and discuss the achievements. The last section (sec. 6) contains the concluding part of this work and the outlook.

Mi j ∇

δE(φ, ∇φ, s) , δc j

126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142

143 144 145 146 147 148 149 150 151 152 153 154

155 156 157 158

∀ i, j ∈ [1, K]

2.1. Capillary force The capillary force is the key driving force in the phasefield application. Exemplarily, coarsening can mostly be driven solely by the grain boundary energies. Limiting the consideration to a two-phase system, naming as α- and β-phases and using the volume fraction constraint φα + φβ = 1, there is only one independent variable φ, e.g. φα = φ and φβ = 1 − φ. The 2

125

∀ α ∈ [1, N]

with the mobility M and the sign due to the maximization or minimization type. The Cahn Hilliard equation for conserved variables, like the composition, writes as

∂t ci = ∇ ·

124

1 a(φ, ∇φ) + w(φ) + e(φ, s)dV  V

depending on the state variables s. The order parameter vector φ changes smoothly between two different bulk volumes along the diffuse interface, the width of which is related to the model parameter  and is comparable to the length scale of the bulk processes. Furthermore, the interfacial energy is approximated by the usage of the so-called gradient energy a(φ, ∇φ) and the potential w(φ), varying only in the diffuse interface. Possible choices of the appropriate functions for the surface energy are studied exemplary in [22]. Using the variational derivatives, the time-dependent free boundary problem for the non-conserved phase-field variable φ is given by the Allen-Cahn equation ∂t φα = ±M

In this manuscript we extend our previous phase-field model in [18] and show simplifications that can be achieved by the assumption of the elastic isotropic materials. Moreover we present an expanded phase-field model, which also consists of chemical contributions. Furthermore, we demonstrate the setup for the verification and validate the presented model in more detail. Since the focus of this work lies in the analysis of the elastic model, we use simple capillary and chemical terms.

AC C

70 71

159 160 161 162 163 164 165

ACCEPTED MANUSCRIPT

σαβ 168 169 170 171 172 173 174

The grand chemical potentials are derived straightforwardly using the Legendre transformation to read

total interfacial energy, fie , between the phases with isotropic surface tension σαβ is approximated by the volume integral, Z

ωα,β (µ) = Aα,β 0 −

Z

Γαβ

dnαβ


1 ≈ a(∇φ) + w(φ) dV.  V

(1)

180 181 182 183 184

185

186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207

SC

M AN U

179

or, equivalently, directly as the difference of the Landau potentials

  16 (2) δφ fie = −γ 2∆φ − 2 (1 − 2φ) . π In the thermodynamic equilibrium, the left side vanishes, and by integrating the right hand side along the normal direction, 2 we get W = π4  for the interface width and σαβ = γαβ for the equality between the physical and the model parameter for the surface tension [5]. 2.2. Chemical model

δφ fch = (ωα (µ) − ωβ (µ))h0 (φ),


ω(µ, φ) = ωα (µ)h(φ) + ωβ (µ) 1 − h(φ) ,

In particular, problems of solidification and solid-solid phase transformation are widely studied in order to understand the material behaviour during the manufacturing processes. Wide spectrums of different alloys over a processing temperature regime are studied and often result in thermodynamic databases such as CALPHAD. The equilibrium phase diagrams for multicomponent alloys are constructed with approximated free energies, wherein the Redlich-Kister approach is used. To perform simulations of manufacturing processes with real materials, simplified mathematical expressions are constructed. In [24], Eiken et al. suggested to evaluate the evolution equations with extrapolated quasi-equilibrium thermodynamic data. Using the parabolic extrapolation scheme, the large scale simulations of the directional solidification can be quantitatively performed, [25]. In our further analysis, we limit our consideration to two components and, for simplicity, assume the same molar volume for both transitioning phases. Furthermore, we assume the Helmholtz free energies f α,β (c) of the different phases (α and β) to depend on the independent concentration cA = c (cB = 1 − c) and to be of a parabolic form with appropriate coefficients Aα,β 0 , α,β Aα,β and A 1 2



 1  D∇ · χ(φ)∇µ − cα (µ) − cβ (µ) ∂t h , χ(φ)

212

213 214

215 216

217 218

219 220

221 222 223 224 225 226 227 228 229 230 231 232 233 234

(7)

with the susceptibility χ(φ) =

211

(6)

as outlined in [4, 5]. The mass diffusion can be given as a tempo-spatial differential equation of both, chemical potential µ and composition c. Since both are derived using Fick’s laws, the procedure of the temporal update of the chemical potential can be calculated explicitly in contrast to the implicit step in the calculation of the concentration update. This is because the interfacial composition is a linear combination of the respective phase concentrations, such that the chemical potential is the same for both. In our computational study we focus on the analysis of the elastic model and, for simplicity, use a chemical one side model, and same diffusivities D = Dα = Dβ in both phases. Therefore, we neglect the anti-trapping current [26]. The evolution equation for the chemical potential writes as ∂t µ =

210

(5)

as a consequent result of the grand chemical potential interpolation

TE D

178

the chemical driving force for the phase transformation writes as    δφ fch = f α (cα ) − f β (cβ ) − µ cα − cβ h0 (φ), (4)

EP

177

209

(3)

Therefore, and under the assumption of the homogeneous chemical potential in the diffuse interface ∂ f β (c) ∂ f α (c) = µ= , ∂c cα ∂c cβ

16 φ(1 − φ), φ ∈ [0, 1] π2 whereby the model parameter γαβ in the potential, as well as in the gradient energy, corresponds to the physical surface tension σαβ . Henceforth, the variational derivative of the interfacial energy is given as

208

.


f (c, φ) = f α (cα )h(φ) + f β (cβ ) 1 − h(φ) .

AC C

176

4Aα,β 2

In the phase-field context [3], the free energy in the diffuse interface is interpolated with the smooth interpolation function h(φ), e.g. h(φ) = φ2 (3 − 2φ) as

On the left side, eq. (1) a sharp interface formulation contains as a line integral over the sharp boundary Γαβ and with its normal nαβ and a volume integral of gradient energy density a(∇φ) = γαβ |∇φ|2 and potential 1 w(φ) in the phase-field model. Note, that the thickness of the diffuse interface depends on the choice of the potential, [23]. In the following, we employ a double obstacle potential given as w(φ) = γαβ

175

2 (µ − Aα,β 1 )

RI PT

166 167

235

1  h(φ) 1 − h(φ)  + , 2 Aα2 Aβ2

and with the concentrations cα,β (µ) = −∂ωα,β (µ)/∂µ = (µ−A1 )/2Aα,β 2 . Because the derivation of the elastic driving force is similar to the derivation of the thermodynamically consistent chemical approach, we want to recount known but important findings. α,β

α,β α,β 2 f α,β (c) = Aα,β 0 + A1 c + A2 c .

3

236 237 238 239

ACCEPTED MANUSCRIPT

246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267

268

269 270 271 272 273 274 275 276 277 278 279

2.3. Microelasticity theory

281 282 283 284 285


 1 ∇u + ∇u T . 2

287 288

290

291 292 293 294

295 296 297 298

(12)

Hence, the elastic driving force for the phase transformation is not the difference of the elastic energies, but is corrected by the scalar product h·, ·i of the homogeneous normal stress with the normal strain jump. The challenge in the determination of the phase corresponding strains and stresses in the computational calculations can be redressed by the implicit calculation of these fields with respect to the homogeneous variables, eq. (9) and (11) and in combination with eq. (8). An implicit calculation scheme was also applied in [28], whereby the same jump conditions are assumed at the common interface. In our recent work [18], we proposed an alternative explicit procedure for the calculation of the elastic driving force and for the homogenisation of the stress in the diffuse interface. Using the Voigt notation for the strain and stress, we write, with respect to the strain and stress tensor components

299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314

    εCT = ε1 , ε2 , ε3 , ε4 , ε5 , ε6 = E11 , E22 , E33 , 2E23 , 2E13 , 2E12     σCT = σ1 , σ2 , σ3 , σ4 , σ5 , σ6 = Σ11 , Σ22 , Σ33 , Σ23 , Σ13 , Σ12 The underscript C represents Cartesian coordinate system, in which the Voigt strain and stress are written. Rewriting them in another coordinate system as εB = MCB εC and σB = NCB σC , where the new base is given, for example, by the orthonormal
vectors B = n, t, s , we note that the coordinate transformation
tensors NCB , MCB are different, but related to NCB T = (MCB )−1 , [18]. Hence, we only use the tensor Mε = MCB in our further derivations. The tensor Mε does not only transform the Voigt strain into a new base, but it is also formulated in such a manner that we can respectively split the strain εB = (εn , εt )T and stress σB = (σn , σt )T into two parts, with respect to the normal and tangential components:

(8)

Since the displacement vector is spatially continuous, the strain tensor evinces a jump through an interface between two solid phases. Under the assumption of the coherent interface, the Hadamard jump condition reduces the strain jump to the jump in the normal strain (E α −E β )n , 0, so that the tangential strains match E α t = E β t,

286

In analogy to the derivation of the chemical driving force, eq. (4) and by using the constraint in eq. (11), we derived the elastic driving force being the partial derivative of eq. (10), with respect to the phase-field variable [18]

Although we consider deformations, as already mentioned, we neglect the dispositional character of the interface movement in our consideration and only deal with the phase transformation. As key results of our next derivations, we firstly show the formulation of the elastic driving force for the phase transition, and secondly we discuss which computational simplification could be achieved, if the materials are assumed to be elastically isotropic. Following the main principles of continuum mechanics, the strain tensor E is given by the symmetrised displacement gradient E=

280

cf. (3). In the chemical model, the phase-dependent concentrations are calculated due to the same chemical potential. In the elastic model, the normal stress is assumed to be the homogeneous variable on the common interface ∂ felβ (E) ∂ felα (E) α n = Σβ n (11) Σ n= n= ∂E E α ∂E β E


 ∂φ fel = felα (E α ) − felβ (E β ) − hΣn, E α − E β ni h0 (φ).

289

(10)

RI PT

245

SC

244


fel (E, φ) = felα (E α )h(φ) + felβ (E β ) 1 − h(φ) ,

M AN U

243

is given in the diffuse interface as an interpolation of the phasedependent strain energies given by the inherent strains

TE D

242

The equality of both derived chemical driving forces, eq. (4) and (5), shows that thermodynamically consistent modelling does not have a unique form, but can be written in terms of any thermodynamic potential. The usage of more convenient thermodynamical potential, depends on the physical problem. For the pure component solidification, and under the assumption of equal density in solid (s) and liquid (l), the use of Helmholtz free energies, f s (T ) and f l (T ), is well established [27]. The derivative of the directly interpolated free energy f (T ) = f s (T )h s + f l (T )hl , with respect to the temperature, automatically results in the interpolated entropy densities S = S s (T )h s + S l (T )hl , totally in accordance with its extensive thermodynamic nature. The solidification of alloys can also be written in terms of Helmholtz free energies, [3], but the interpolation of phasedependent and extensive compositions is postulated due to the same chemical potential. Furthermore, the compositions differ as arguments of both phase-dependent free energies, which is a fact, that still complicates the analysis but only slightly. On the other hand, using the interpolation of grand potentials, which are written in terms of the intensive variables ω(T, µ), eq. (6), the appropriate interpolations of the extensive variables are consequently given by the corresponding differentiations. Moreover, the derivations of the evolution equations, written with respect to the chosen intensive variables are straightforward and without any additional constraints. Therefore, explicit and less complex numerical methods can be used in the computational calculations.

EP

241

AC C

240

(9)

with n as the normal vector of the interface and t as any vector perpendicular to n. Nevertheless, as a consequence of the strain jump, the elastic energy, as a part of Gibbs free energy, 4

i i i T εin =(Enn , 2Ent , 2Ens )

and

εt =(Ett , E ss , 2Ets )T ,

σn =(Σnn , Σnt , Σns )T

and

σit =(Σitt , Σiss , Σits )T ,

315 316 317 318 319 320 321 322 323 324 325 326

ACCEPTED MANUSCRIPT 327 328 329 330 331

In order to highlight the phase-dependent variables we wrote them with the superscript i ∈ {α, β}. In particular, the homogeneous variables are written without superscripts and are the same for both phases, that are in the diffuse interface. Thus the strain energy of i-phase writes as

 2(1−νi )  1−2νi  σn = Gi  0  0

 2νi   1−2νi 0 0  i i i 1 0 (εn − ε˜ n )+G  0   0 0 1

2νi 1−2νi

0 0

 0
 0 εt − ε˜ it  0

and

334 335 336

339 340 341

1 2



− ∆εn σn + εTt ∆σt − ∆ f˜el h0 (φ)

β

344 345

∂p(φ, εt , σn ) = εαn hα + εαn hβ ∂σn ∂p(φ, εt , σn ) σt (εt , σn ) = = σαt hα + σβt hβ , ∂εt 346 347 348 349 350

351 352 353 354 355 356 357 358

AC C

εn (εt , σn ) = −

  νi  1−2νi νi  1−νi 1−ν  2−2νi 0 0 i    0 1 0 σn −  0 0    0 0 1 0 0 {z } | {z i Tnn

1 σn p (εt , σn ) = 2 εt − ε˜ it i

(15)

feli (εin (σn , εt ), εt )

EP

343

(14)

− of the Legendre transformations p = σn εin (σn , εt ) with respect to the phase-field variable φ, cf. eq. (6). Furthermore, applying the differentiation with respect to the homogeneous variables i

0

 0  i 0 (εt − ε˜ t )  1

Tnti

!T

i Tnn Ttni

! ! Tnti σn − σTn ε˜ in Ttti εt − ε˜ αt

which, in terms of homogeneous variables, writes as

 νi  1−νi  νi σit =  1−ν  i 0

  2 2νi  1−νi 1−ν 0 0 i   i  2νi 2 0 0 σn + G  1−ν   i 1−νi 0 0 0 0 | {z Ttti

(16)

we conclude the interpolation of inhomogeneous variables in the diffuse interface. Since the additive character of the normal strain complies with the extensive volume, the interpolated tangential stress also seems to be extensive to result in a surface stress contribution for the interfacial energy. 2.3.1. Elastically isotropic materials By the assumption of the elastic isotropy in both phases, the stiffness tensors is expressed with two independent parameters. In our further derivations we use the Poisson ratios νi and the shear modulus Gi as constants to relate elastic strain and stress. If both phases are elastically isotropic, the stiffness tensors are independent of the chosen coordinate system, and the normal and tangential stresses are given by 5

362

 0  i 0(εt − ε˜ t )  1 }

ε¯ in =ε˜ αn + Ttn ε˜ it

and

ε¯ n = ε¯ αn hα + ε¯ βn hβ ,

σ ¯ it =Ttti ε˜ it

and

σ ¯t = σ ¯ αt hα + σ ¯ βt hβ .

363 364 365

366 367 368

369

(20)

In the previous equations (18), (19) and (20) we introduced the proportionality tensors between the homogeneous and inhomogeneous variables in the notation of [18]. Since in our recent work the determination of T i is formulated in a general form, its explicit calculation is dependent on the interface normal. It is not the case, if both materials are assumed to be elastically isotropic. They rather only depend on elastic bulk constants, i eqs. (18) and (20). While the tensor Tnn is a compliance, and Ttt has the unity of stiffness, the tensors Tnt = TtnT , with mixed indices, are dimensionless. In contrast to the direct formulation of the proportionality tensors Tkli , the eigenstrains of each phase are related to the interface orientation. Using the coordinate transformation, we introduce the following auxiliary quantities for convenience

(17)

361

(19)

It is also interesting to note that the differentiation of the strain energy and of the interfacial elastic potential, both with respect ∂fi to the tangential strain result in the tangential stress σit = ∂εelt = ∂piel ∂εt ,

360

 0
 0 εt − ε˜ it + ε˜ in (18)  0 }

Knowing the explicit dependencies, eq. (18), we rewrite the new interfacial elastic potentials given as a Legendre transform, feli (εin (σn , εt ), εt ) − σn εin (σn , εt ), as

TE D

p(φ, εt , σn ) = p (εt , σn )h(φ) + p (εt , σn )(1 − h(φ)), 342

1 εin = i G |

An equivalent formula we find also in [11] and similar to the derivation of the chemical driving force with a grand potential formalism, we derive the same elastic driving force by differentiating the interpolation α

2νi 1−2νii 2(1−ν ) 1−νi

respectively. Thus, the normal phase dependent strains can be written directly as



1 − εαn − εβn T σn + εTt σαt − σβt 2  
− ε˜ αn − ε˜ βn T σn + (ε˜ αt )T σαt − (ε˜ βt )T σβt h0 (φ).

∂φ fel = 338

ε˜ in

or, written in a shorthand notation, it leads to

 2(1−νi )   1−2νi 0 0  i i  2νi i 0 0 (εn − ε˜ n )+G  1−ν i   0 0 0

 2νi  1−2νi  2νi σit = Gi  1−2ν i  0

where ε˜ denotes the eigenstrain and are its tangential and normal parts, respectively. Thus, by the substitution of the eq. (13) into eq. (12), the elastic driving force for the phase-field evolution equation results in

∂φ fel =

337

(13)

RI PT

333

ε˜ it ,

i

359

M AN U

332

 1 = (εin − ε˜ in )T σn + (εt − ε˜ it )T σit 2

SC

feli

370 371 372 373 374 375 376 377 378 379 380 381 382 383

ACCEPTED MANUSCRIPT

387 388 389 390

−1 −1 −1 σn =Tnn εn + Tnn Tnt εt − Tnn ε¯ n

different mobilities in both perpendicular and longitudinal directions, given with harmonic and direct interpolations types, respectively. Similarly, it follows that the diffuse interface between two elastically isotropic phases is elastically anisotropic. Consequently, the effective stiffness tensor, can be written as a sparse matrix  1−ν  ν  0   0  C¯ B = λ⊥  1    1  0

(21)


−1 −1 −1 σt =Ttn Tnn εn + Ttt + Ttn Tnn Tnt εt − σ ¯ t − Ttn Tnn ε¯ n (22)

393 394 395

396 397

399 400 401 402 403 404

 2(1−ν) ! hβ −1  1−2ν hα  0 + =  Gα Gβ | {z } 0 G⊥

AC C

 2    1−ν 2ν Ttt = Gα hα + Gβ hβ  1−ν | {z } 0 k G

405

407 408 409 410 411 412 413 414 415 416 417 418 419

ε˜ t = Ttt−1 σ ¯t

 0 0  1 0 ,  0 1

⊥

2ν 1−ν 2 1−ν

0

1 0 0

1 0 0

0 0 0 0

1−2ν 2ν

0 0 0

R+ν2 ν(1−ν) R+ν 1−ν

R+ν 1−ν R+ν2 ν(1−ν)

0

0

0 0 0 0 0 R 2ν(1−2ν)

2G⊥ ν 1−2ν .

But the stiffness tensor is dense in Cartesian and

421 422 423 424 425 426

427 428 429 430

CC = MεT CB Mε .

 ν  1−ν  ν α β ε¯ n = ε˜ hα + ε˜ hβ +  1−ν  0

431 432 433

 0 0    0 0 ε˜ αt hα + ε˜ βt hβ − ε¯ t  0 0

with

434

ε¯ t =

 1  α α β β G ε ˜ h + G ε ˜ h α β t t Gk

3. Validation of the elasto-chemical phase field model

 0  0 .  1

420

       ,   

Likewise, the expression for the effective tensor can be given explicitly in a shorthand notation, the effective eigenstrain in the diffuse interface interpolates as

(24)

2.3.2. Solid phases with the same Poisson ratios The effective stiffness tensor can be further simplified if both elastically isotropic phases have the same Poisson ratios να = ν = νβ . If this is the case, the dimensionless blocks Tntα = Tntβ = Tnt match, and the interpolations in Tnn and Ttt result in the interpolations of the shear moduli. Thus, we write the following for the interpolated proportionality tensors

Cnn

406

eter λ = writes as ⊥

EP

398

and

0 0

in the base B with phase-field dependent, dimensionless paramk eter R = G (φ)(1−2ν) and harmonic interpolated first Lam´e paramG⊥ (φ)

−1 −1 −1 given as Cnn = Tnn , Cnt = Tnn Tnt and Ctt = Ttt + Ttn Tnn Tnt ; and the effective eigenstrain results in

ε˜ n = ε¯ n − Tnt Ttt−1 σ ¯t

1−2ν 2ν

M AN U

392

where we also used the interpolation of the appropriate tensors Tkl = Tklα hα + Tklβ hβ . Alternatively we can write the stress tensor as a matrix vector product with the effective stiffness tensor C and the effective eigenstrain ε˜ ! ! C Cnt εn − ε˜ n σB = nn (23) Ctn Ctt εt − ε˜ t

TE D

391

0

RI PT

386

In the next step, we interpolate inhomogeneous variables, the normal strain εn (σn , εt ) and the tangential stress σt (σn , εt ), respectively in accordance with eq. (16) and (17). By the inversion of the expression of normal strain ε−1 n (σn , εt ), we derive the normal stress depending on strain σ(εn , εt ). Then, we substitute the normal stress in the tangential stress and can write the following for the stress components

SC

384 385

435

In the previous sections we derived thermodynamically consistent driving forces for the phase field evolution, eqs. (2), (5), and (14). Subsequently, the phase-field equation writes as

436 437 438

k

The different interpolation types G and G of the shear moduli for the different directions are remarkable. The harmonic interpolation type of the shear moduli is used in the formulation of the stiffness block, which relates the homogeneous normal stress to the inhomogeneous normal strain. But for the tangential fields we get direct interpolation of the shear modulus. A relative form, which is based on the similar decomposition in normal and tangential directions can be found in the formulation of the mobilities in the diffusion/heat conduction equations, wherein scalar and thus isotropic bulk mobilities in the diffuse interface receive a tensorial form, see cf. [26, 29, 30]. Therefore, the mass/heat fluxes in the bulk, which are given as a scalar scaled conservative vector field, can be interpreted as spatially isotropic. But in the diffuse interface, the fluxes are spatially anisotropic, because the gradient of the potential is scaled with

  τ∂t φ = − δφ fie + δφ fch + δφ fel ,

(25)

With the equation for the strain calculation, eq. (8) together with the Hook’s law, which is written in the Voigt notations as σ = MεT CB (Mε ε− ε˜ B ), the momentum balance equation follows as

439 440 441 442

ρu¨ = ∇ · Σ. whereby the relations between the tensorial fields (E, Σ) and the corresponding fields in the Voigt notation (ε, σ) are selfexplanatory. For a model without the chemical part (δφ fch = 0), the match between analytical and numerical solutions in 1-D is presented in [18]. We performed 2-D simulations of a circular inclusion 6

443 444 445 446 447 448

ACCEPTED MANUSCRIPT

453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476

µ =const, 0 =σαβ κ

480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499

αβ

+ ∆ω (µ) +

∆pαβ el .

G H(x, y)εˆ 1−ν ε M (x, y) =C−1 σ M (x, y)

σ M (x, y) =

(28)

Eq. (26) means that in the end state there is no mass flux ∇µ = 0 → ∂t c = 0) and, furthermore, the chemical driving force ∆ωαβ (µ) is constant along the interface. The mechanical equilibrium, eq. (27) implies steady distribution of the elastic fields with vanishing traction vector Σα nαβ = Σβ nαβ on the common interface and the Gibbs-Thomson equation (28) balances the forces acting on the interface, namely (1.) capillary force, written in terms of the isotropic surface tension σαβ and the interfacial curvature καβ , (2.) constant chemical force eq. (5) and αβ αβ (3.) elastic driving force ∆pαβ el = ∆ fel − σn ∆εn , [31]. If the force balance is achieved, the interface is stationary. Whereby the anisotropy of the interfacial energy due to the inhomogeneous elastic field distribution is assumed to be well-posed. By the comparison of eqs. (25) and (28), the similarity of both equations is obvious, and thus the diffuse interface formulation should resolve the sharp interface condition [32]. In order to validate the presented model, we investigate a setup, for which all required quantities and, in particular, the elastic fields are known. Pursuing our goals, we found the Eshelby inclusion as an ideal setup for the validation. For an elliptical inclusion or inhomogeneity in the embedded infinite matrix, which is without any residual strains, the elastic fields are given in a closed-form solution. Herein, we differ between

502 503 504 505 506

507 508 509 510 511 512

513 514

515 516 517 518

(31) (32)

with H(x, y), depending on the coordinate (x, y) outside the inclusion. Note the equality of both shear modules G = G M = G Inc .

σ2 (x)

501

(30)

For completeness all used tensors in the calculations are listed in Appendix A and are simplified formulated specified for the plain strain case due to the derivations in [34]. The elastic fields outside the inclusion are given as

EP

479


σInc = C(εInc − ε) ˆ = C S(t, ν) − 1(4×4) ε. ˆ

(27)

500

(29)

The corresponding constant inner stress field is given by Hook’s law as

519 520 521

σ4 (x)

y

σ1 (y)

AC C

477 478

εInc = S(t, ν)ε, ˆ

(26)

0 =∇ · Σ, αβ

3.1. Elliptical inclusion in the infinite matrix For an elliptical inclusion with a ratio of semi-axes of t = b/a, see Figure 1, the Eshelby tensor S depends on the geometry of the inclusion t as well as on the Poisson ratio ν = ν M = νInc . For the plane strain eigenstrain of the inclusion, we use the notation εˆ = (εˆ 1 , εˆ 2 , 0, εˆ 4 )T and the constant inner strain reads

RI PT

452

SC

451

both, in the sense of Eshelby [33], whereby the inclusion has a non-vanishing eigenstrain but the same stiffness as the environmental matrix (CInc = C M ), and the inhomogeneity is also characterised by the different stiffness tensor (CInh , C M ) and the eigenstrain [34]. We are considering a 2-D setup and for simplicity we limit our analysis and following simulations on the plane strain case.

in a surrounded matrix, which was influenced by elastic transformational and capillary forces. In the thermodynamical equilibrium (∂t φ = 0) and in the mechanical equilibrium (∇ · Σ = 0), the capillary and configurational forces equates δφ fie = −δφ fel , which is also known as a Gibbs-Thomson equation in the sharp interface description σαβ κ = ∆ fel −hΣnαβ , ∆F nαβ i, [31]. Hence, by matching the equilibrium properties for the sharp interface, and by applying them in the diffuse interface simulations, we can assign the accordance between the modelling and the analysis. The disadvantage in this procedure is the metastability of the quasi-equilibrium state. Thus the equality of capillary and elastic forces is not stable and any small perturbation of the ideally equilibrated system results in a growth or shrinkage of the circular inclusion. Therefore, we performed a series of simulations with same initial circular shape and elastic field contributions, but with a different surface tension parameter around the analytical value. The assignment of the numerical parameter to the analytical value was done by the evolution change from a shrinking (∂t φ < 0) to a growing (∂t φ > 0) behaviour and consequently from balancing the capillary and elastic driving forces. In this work we incorporate the chemical force into the phase field equation to overcome the quasi-equilibrium state between the elastic transformational and the capillary forces in the thermodynamic equilibrium. Hence we solve the mass diffusion equation. Furthermore the thermodynamical and mechanical equilibrium are given by the conditions in the sharp interface formulation

M AN U

450

TE D

449

−σ4 (y)

b

σ1 (y) a

x

σ4 (y)

σ2 (x) −σ4 (x)

Figure 1: Simulation setup of an elliptical inclusion in a matrix phase.

Since all elastic fields are determined by the eigenstrain ε, ˆ we write the following for the elastic potential jump depending solely on εˆ 7

522 523 524

ACCEPTED MANUSCRIPT ∆pel (ε) ˆ = ∆ fel (ε) ˆ − σn (ε)∆ε ˆ ˆ n (ε) 525

(33)



1 ˆ ∆pel (ε, ˆ CInh ) = ∆pel (ε) ˆ + σInh CInh −1 C M S − 14×4 ε, 2

In particular, the elastic driving force at both vertices is

changes with varying stiffness constants in the inhomogeneity CInh . For this reason, the chemical force should be adopted for every inhomogeneity setup in order to satisfy the GibbsThomson equation for the chosen capillary force. In this manner, we evaluate the elastic driving force at the vertex ∆pel (a, 0), add it to the capillary force, and thus obtain the corresponding jump in the chemical force, the chemical potential, and consequently the compositions for both phases. The foregoing procedure allows to perform a more detailed and accurate validation of the elastic driving force, because different scenarios can be resolved with a minimal computational effort.

G (2 + t)εˆ 21 + 2tεˆ 1 εˆ 2 + (t2 − t − 1)εˆ 22 ∆p(a, 0) = , 1−ν (1 + t)2 and

RI PT

526

G (1 − t − t2 )εˆ 21 + 2tεˆ 1 εˆ 2 + t(1 + 2t)εˆ 22 , ∆p(0, b) = 1−ν (1 + t)2

531 532 533 534 535 536 537

538

539 540 541 542 543 544 545 546 547 548

3.2. Equivalent inclusion method

In the case, in which the stiffness tensor of the embedded elliptical precipitate differs from its counterpart in the matrix CInh , C M , we speak of inhomogeneity. Furthermore, denote ε˜ Inh for the eigenstrain in the inhomogeneity. To calculate the corresponding elastic fields in this setup, a spatial equivalent system, but with an inclusion, is considered. The challenge herein is the determination of the eigenstrain εˆ for this equivalent inclusion. By using the equivalent internal stress σInh = σInc and the strains εInh = εInc in both the inhomogeneity and the inclusion, the following equation should be valid 
CInh − C M S + C M εˆ = CInh ε˜ Inh .

(35)

551

Thus, for a given residual strain of the inhomogeneity ε˜ there is a unique equivalent inclusion with a corresponding eigenstrain.

552

3.3. Analytical solution

550

553 554 555 556 557 558 559 560

SC

Consequently, for the desired ellipse with the semi-axes a and b, and with the surface tension σI M , the components of the residual strain in the inclusion are not arbitrary, but should satisfy the foregoing equality. Hence, the previous equation prescribes the match between the elastic constants of the equilibrated inclusion on the left and the corresponding capillary values on the right.

 549

4. Simulation setup

562 563 564 565 566 567 568 569 570 571 572

573

As a first step, we perform theoretical calculations of the elastic fields for an ellipse with semi-axes a = 60 and b = 90 corresponding to the final equilibrium state. The dimensionless shear modulus and the Poisson ratio in the matrix are G M = 840 and ν M = 0.25, respectively. With the surface tension σI M = 0.0942214, the components of the residual strain are chosen as εˆ = (0.00205, 0.00165, 0, 0), with respect to eq. (34). The analytical expressions for the inner and outer elastic fields (29)-(32) are derived for an infinite plane strain plate. To replicate such a scenario we select a rectangular simulation domain with large dimension L x = Ly = 429 and use ∆x = ∆y = 1 for the spatial step. Thus, to adopt the analytical and computational results, we calculate the stress fields along the boundary of our quadratic simulation domain, approximate them with the suitable approach see Figure 2 and impress them as a stress boundary condition, see Figure 1. To reduce the influence of the boundary load, the distance of the inclusion to the boundary is several times longer than its semi-axes. For the chemical model we employ a vanishing grand potential for the inhomogeneity ωI = 0 and the following functions

M AN U

530

TE D

529

EP

528

respectively. Moreover, the curvatures at the vertexes are known, κ(a, 0) = ba2 and κ(0, b) = ab2 . Then, writing the GibbsThomson equation (28) for both critical points of the ellipse, and taking the difference of both, results in !
G b a εˆ 2xx − εˆ 2yy = σI M 2 − 2 . (34) 1−ν a b

AC C

527

561

574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594

ω M (µ) =0.0703(µ + 1.18108)(1 + 0.083215(µ + 1.18108)),

Inh

cI (µ) =0.4682 + 0.0722(µ + 1.18108),

χI = 0.0722,

c M (µ) =0.3979 + 0.0605(µ + 1.18108),

χ M = 0.0605.

In both phases, the dimensionless diffusion coefficient is the same D M,I = 50. On the boundary of the simulation domain we apply the homogeneous Neumann boundary condition for the mass diffusion, so that the total composition inside the simulation domain remains constant. This fact allows a more accurate adjustment of the composition in the domain. In general, the final shape of the inhomogeneity should be the same in all simulations, only the compositions of different phases should be in accordance with the constant chemical potential µ, ¯ satisfying the Gibbs-Thomson equation (28).

In our analysis, we invert the presented approach and fix the eigenstrain εˆ of the inclusion with respect to the matching condition (34). By the variation of the elastic constants in the inhomogeneity CInh , we adopt the appropriate eigenstrain ε˜ Inh with respect to eq. (35). Doing so, the elastic fields are given with respect to eqs. (29)-(32) and thus is calculated once. In spite of the fact that the elastic field contribution remains the same for all setups, the elastic driving force 8

595 596 597 598 599 600 601 602 603 604

ACCEPTED MANUSCRIPT Σ xy H y L

Σ xx H y L -100

100

200

0.15

630

5.1. Matches in the form

631

0.10

-0.05

0.05

-0.10 -200

-0.15

-100

100

200

y

-0.05

-0.20

-0.10

-0.25

-0.15

(a) σ1 on right/left boundaries. -100

In the first considered simulation scenario (inclusion), the numerical shapes, compared to the theoretical form, are shown in Figure 3.

(b) σ4 on right boundary.

Σ yy H x L

-200

order to investigate its influence on the driving forces.

Σ xy H x L 100

200

100

0.10

-0.05

0.05 -200

-100

100

200

60

x

-0.05

-0.15

100

80

-0.10

40

-0.10

20

-0.20

(d) σ4 on top boundary.

0

0

20

40

60

Wκ ≈ 0.25 Figure 2: Black solid lines and green dashed lines are the theoretical profiles and their respective approximation of the corresponding stress components on the respective boundary as external load. G Inh 840 420 1260 840 840

νInh 0.25 0.25 0.25 0.1 0.3

ε˜ Inh xx [%] 0.205 0.2476 0. 1908 0.228848 0.197051

ε˜ Inh xx [%] 0.165 0.243067 0.138978 0.184592 0.158469

Inh [%] ε˜ zz 0.0 0.0 0.0 0.02896 -0.00965

µ¯ [-] -1.10995 -1.09204 -1.11593 -1.10179 -1.11267

609 610 611 612 613 614 615 616 617

618

619 620 621 622 623 624 625 626 627 628 629

TE D

608

EP

607

In summary, we test different scenarios with inhomogeneities related to different stiffness parameter G Inh and νInh . The corresponding eigenstrain ε˜ Inh and the appropriate chemical potential are listed in Table 1. The first scenario is the case of the inclusion. In the second and third test cases, the inhomogeneity is one softer and one harder, but the Poisson ratio of the inhomogeneity equates to the Poisson ratio in the matrix. In the last two scenarios, we use the same shear module as in the matrix, but vary the Poisson ratio. In this way, the eigenstrain component ε˜ zz is not zero. Furthermore, in order to investigate the interfacial excess energy, in the case, where the elastic driving force is modelled incorrectly, we vary the diffuse interface width for each test case.

60

40

40

20

20

0

20

40

60

Wκ ≈ 0.5

80

0

0

20

40

60

80

Wκ ≈ 0.75

The images from left to right depict the final formation for different diffuse interface width. Because of the symmetry, only a representing quarter of the ellipse is illustrated. Since the solid line corresponds to the theoretical shape, the dashed lines sign the contour lines in the diffuse interface with φ = 0.1, 0.5, 0.9. With increasing interface width, the locations of the locus φ = 0.5, as a representation of the sharp contour and the sharp interface position show an insignificant mismatch, Figure 3. We get similar final forms in all listed scenarios Fig 4, except for the second setup in Table 1 as demonstrated in Figure 5. In this case, the shear module of the inlay of G Inh = 0.5G M = 420 is a half of the shear module of the matrix. Therefore, the inhomogeneity is much softer than its environment. Herein, we observe the strong dependency of the inhomogeneity form on the diffuse interface width, Figure 5. Thus, with increasing diffuse interface width, the inhomogeneity becomes rounder and in the next section we will try to explain this discrepancy, but first we compare the chemical potentials. 5.2. Matches in the chemical potential

AC C

605

60

Figure 3: Comparison of the theoretical (red) and numerical (black) shapes.

Table 1: Parameter for the analytical and numerical calculations

606

80

0

633

100

80

M AN U

sc. 1. 2. 3. 4. 5.

80

632

νInh = 0.25

SC

(c) σ2 on top/bottom boundaries.

1. G Inh = 840,

x

RI PT

-200

y

The values of the final chemical potentials in all simulations are plotted in Figure 6. Using the definition of the characteristic width on the abscissa, the ordinate corresponds to the sharp interface case (W = 0). In this plot, we see an excellent agreement between the theoretical predictions as lines and the numerical results as dots in all scenarios, see Table 1. This is also the case in the second scenario, in spite of the fact that the form deviates significantly. Herein, the highest relative deviation of the chemical potential is at around 0.5‰. Consequently, the chemical driving force will be resolved accurately for all scenarios, independent of the diffuse interface width and elastic constants. Thus, if the chemical potential, and subsequently the chemical driving force, matches well with the theoretical predictions, the discrepancy in the capillary force, as a result of incorrect shape Figure 5, should be implied by the elastic contribution.

5. Results and discussion

Since in the thermodynamic equilibrium there are three contributions, namely capillary, chemical, and elastic driving forces, eq. (28), we analyse the corresponding final numerical shapes, the chemical potentials and the elastic field distribution as representative quantities for all balancing driving forces. Doing so, we introduce the characteristic width Wκ as the dimensionless quantity, given by the product of the diffuse interface 2 width W ≈ π4  with the maximal curvature of the ellipse κ = 0.025. We report on simulations with Wκ ≈ 0.25, 0.5, 0.75, and also at around 10, 20 and 30 grid cells in the diffuse interface for the accurate resolving of the capillary force and in 9

634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652

653

654 655 656 657 658 659 660 661 662 663 664 665 666 667 668

ACCEPTED MANUSCRIPT νInh = 0.25

80

80

80

60

60

60

40

40

40

20

20

20

0

0

0

20

40

60

80

0

Wκ ≈ 0.25

20

40

0

80

100

60

60

60

40

40

40

20

20

20

0

0

0

40

60

80

0

20

40

60

εI =(0.001624, 0.00086933, 0, 0)T

80

60

80

0

20

40

60

5. G Inh = 840, 100

80

νInh = 0.3

σ1 (x)

−1,4

100

80

−1,5

80 60

40

3

40

20

10

20

30

40

50

60

70

σ1 (y)

1

20

−1 −2

0

0

20

40

60

0

80

0

Wκ ≈ 0.25

20

40

60

0

80

0

0

Wκ ≈ 0.5

20

40

60

Wκ ≈ 0.75

G Inh = 420, 80

60

60

40

40

20

20 0

20

40

60

80

80 60 40 20

0

20

40

60

80

Wκ ≈ 0.5

0

0

20

40

60

80

Wκ ≈ 0.75

Figure 5: Comparison of the theoretical (red) and numerical (black) shapes

−1,1

−1,11 −1,12

0

AC C

−1,09

µ [−]

0,2

0,4

0,6

0,8

672

682

σ3 (x) −0,5

40

60

80

100

0

−1

10

20

30

−1 40

50

60

70

σ2 (y)

−1,5

0

10

20

30

40

50

60

70

σ3 (y)

0

−0,5 −2 −1 −2,5

0

20

40

60

80

100

0

20

40

60

80

100

1

5.3. Matches in the elastic fields The distribution of the elastic fields in the matrix is effected by both, the boundary load and the residual strain in the inhomogeneity. Because of the boundary load, due to the analyti-

6. Conclusion and Outlook 669

Summarizing our work, we can denote following findings • The elastic driving force formulation is based on the mechanical jump conditions and is derived by the variational approach.

10

683 684 685 686 687 688 689 690 691 692 693

i

si = components are given in percent in Table 2. Except the second test case we get a remarkable agreement in all scenarios. As long as the computational shape correlates to the theoretical form, Figure 3 and Figure 4 the elastic fields match the theoretical profile. Moreover, in the last column of Table 2 we also give the maximal local relative deviation of the strain energy in the inhomogeneity, and though the characteristic width Wκ ≈ 0.75 is of order one, the strain energy deviates less than 5% from the theoretical value.

Figure 6: Lines correspond to the theoretical chemical potentials from Table 1. The points around the solid lines correspond to the appropriate scenario in dependence on the diffuse interface width on the abscissa.

671

20

−3

num σana i −σi σana i

Wκ [−]

670

681

κ W ≈0.75

In this figure we exemplary plot the final stress distribution along the axes for the first test case. The black lines correspond to the analytical profiles and the coloured lines represent the stress profiles for the appropriate characteristic widths. Independent on the interface width, we observe excellent agreement between the numerical and analytical solutions. The major discrepancies are observed in the inhomogeneity. But nevertheless, as long as the computational form is comparable with the analytical shape, Figure 4, the elastic fields contributions well correspond with the analytical prescriptions. The maximal relεana −εnum ative discrepancy of the elastic strain ei = iεana −iε˜ i and stress

EP

Wκ ≈ 0.25

0

100

680

Figure 7: Comparison between the computational and analytical stress profiles along the symmetry axes x and y.

TE D

80

0

νInh = 0.25

100

679

80

Figure 4: Comparison of the theoretical (red) and numerical (black) shapes.

100

678

0

σ2 (x)

0

0

20

κ W ≈0.5

−2

0

2

40

1

M AN U

60

677

−1

−1,7 −1,8

60

κ W ≈0.25

analytic

Wκ ≈ 0.75

−1,6

80

676

σ =(−1.72928, −2.32512, −1.0136, 0)T

−1,3

100

675

and

In the sharp interface consideration, the stress and strain components exhibit jumps throughout the common phase boundary. Whereas in the diffuse interface, all elastic fields vary smoothly. The respective profiles are given in Figure. 7 for the scenario of the inclusion.

Wκ ≈ 0.5

674

I

100

80

20

40

νInh = 0.1

80

Wκ ≈ 0.25

20

Wκ ≈ 0.75

80

0

0

Wκ ≈ 0.5 4. G Inh = 840,

100

60

cal solution, the numerical and theoretical elastic fields distributions correlate well in the matrix. Certain deviations occur solely in the immediate vicinity of the diffuse interface. Following the Eshelby theory, the strain eq. (29) and stress eq. (30) components in the inclusion/inhomogeneity are constant and of

673

100

RI PT

100

SC

3. G Inh = 1260,

100

694 695 696 697 698 699 700 701 702

703

704

705 706 707

ACCEPTED MANUSCRIPT

1

3

4

5

2

Wκ 0.25 0.5 0.75 0.25 0.5 0.75 0.25 0.5 0.75 0.25 0.5 0.75 0.25 0.5 0.75

e1 -2.8 -4.1 -5.7 -1.3 -2.0 -2.7 -1.9 -3.3 -4.6 -3.0 -4.4 -5.8 -9.4 -31.2 -42.7

e2 1.7 2.5 3.4 1.0 1.6 2.2 1.2 2.1 2.9 1.9 2.7 3.5 5.4 17.8 24.5

s1 -1.9 -2.8 -3.9 -1.4 -1.8 -2.5 -1.4 -2.4 -3.4 -1.9 -3.0 -4.1 -3.3 -10.9 -14.9

s2 1.2 1.7 2.3 1.1 1.8 2.5 1.0 1.6 2.3 1.2 1.7 2.3 1.8 5.7 7.9

s3 -0.12 -0.24 -0.36 0.06 0.22 0.44 -0.02 -0.05 -0.06 -0.19 -0.37 -0.54 -0.39 -1.41 -1.84

∆ fel -1.44 -2.17 -3.0 -0.65 -0.91 -1.21 -0.98 -1.69 -2.31 -1.52 -2.33 -3.22 -4.10 -12.3 -15.5

744

714 715 716

717 718 719 720

721 722 723 724

725 726 727 728 729 730

731 732 733 734 735 736 737 738 739 740 741 742 743

M AN U

713

Based on the numerical results we conclude the following. In the scenarios in which the Young or the shear modules of the transitional phases are comparable we observe:

Acknowledgements

• The location of the diffuse interface overlaps the sharp contour line in such a manner that the isoline φ = 0.5, as a representation of the sharp interface, insignificantly deviate from the theoretical formation.

TE D

712

• The determinations of the homogenised stiffness tensor and eigenstrain significantly simplifies for isotropic materials

• The incorporation of the chemical contribution into the phase-field model allows to stabilize the equilibrium shape. The numerical deviation in the chemical potential is less than 0.5‰.

EP

711

• The correspondence between the sharp interface approach and the phase-field model was presented for the equilibrium state.

• As long as the numerical and analytical shapes match well, the elastic field distribution distinguished agree with the theoretical predictions. Furthermore, although the characteristic widths are chosen extremely large Wκ ≈ 0.25, 0.5, 0.75, the capillary, chemical and elastic fields show a negligible dependency on the interface width.

Appendix A. Relevant tensors




2(1+t)(1−ν)+1 t 2(1+t)2 (1−ν) 2(1+t)ν−t 2(1+t)2 (1−ν)

2(1+t)ν−1 t 2(1+t)2 (1−ν) 2(1+t)(1−ν)+t 2(1+t)2 (1−ν)

νt (1+t)(1−ν) ν (1+t)(1−ν)

0 0

0 0

0 0

0 0 0

x˜2 − 1 y˜ 2a − t2 R(x, y) = a + + 2 2

1−

t (1+t)2 (1−ν)

s x˜a2 − 1 y˜ 2a − t2 − 2 2

11

r

748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767

a 1 + aR(x, y)

769 770 771 772

774 775 776

777 778 779 780

!2 + x˜a2 y˜ 2a

and further write ρa (x, y) =

747

      

and is used to calculate the fields inside the inclusion in eqs. (29) and (30). The following auxiliary quantities are required for the formulation of the tensor H(x, y). Utilizing the short hand notation of x˜a = x/a and y˜ a = y/a we introduce

Albeit the noted excellent agreements, the numerical results for the second scenario indicate the necessity of a further extension of the presented model. In the particular setting, the inhomogeneity is much softer than the matrix and, with increasing interface width, the form of the inhomogeneity becomes rounder. Although the chemical potential in the final state excellently agrees with the theoretical value, the elastic driving force is falsified to some extent, so that the formation disagrees from the theoretical shape. Because of the inhomogeneity shape, which is not appropriate to the analytical form, the stress and strain fields strongly deviate from the analytical distributions. Since the calculation of the homogenised stress tensor can be presumed as consistent, we expect the evaluation of the elastic

746

773

Some following auxiliary quantities used for the validation are written in a simplified manner of the derivation presented in [34]. Thus the Eshelby tensor for the elliptical inclusion reads


745

768

The authors thank the German research foundation for the funding of the project NE 822/16-1. We are grateful to the computational time provided by the High Performance Computing Center in Stuttgart.

    S =   

AC C

710

SC

Table 2: Maximal relative deviation of relevant elastic quantities in percent for different scenarios and varying interface width

709

driving force to be improved. This issue will be considered in a forthcoming work. Basically, the elastic driving force was derived by the variational approach and we considered solely the variational derivative of the interpolated strain energy with respect to the phase field variable. But because of the mechanical jump conditions through the interface, the inhomogeneity in the phase own strains and stresses depending on the interface orientation, given by the gradient, ∇φ. Consequently, the variational derivative with respect to ∇φ should be taken into account to ensure the complete variation of the objective functional. Thus, also for other models [17, 28] describing elastic inhomogeneous materials in the phase-field context and with missing variational derivative of the interfacial elastic potential jump, with respect to the phase-field gradient, the deviation in the numerical results are to be expected. Likewise, the presented formalism impedes the explicit calculation of the missed term. In our next works 708we aim to overcome this gap and to extend the elasto-chemical phase field model for multi-phases. Nevertheless, if the solid-solid phase transformation takes place and the elastic driving force plays a role, then, in the most cases, both stiffness tensors are slightly different. These cases are well captured by the present model wherein neglecting the additional term, like was shown in our presented study.

RI PT

sc.

781

and

ρb (x, y) = t

r

a , t + aR(x, y)

ACCEPTED MANUSCRIPT

784

x˜a m x (x, y) = 1 + R(x, y) 785

786

y˜ a my (x, y) = 2 . t + R(x, y)

and

Normalization gives my (x, y) m x (x, y) n x (x, y) = and ny (x, y) = . |m(x, y)| |m(x, y)| p with |m(x, y)| = m x (x, y)2 + my (x, y)2 . Using L(x, y) = ρa (x, y) ny (x, y) + ρb (x, y) n x (x, y) , 2

h12 =h44 = h13 h23

ρa ρ2b


ρa t+ A A

− L + 4n2x n2y (L − 1),

and finally write

793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812

h13 h23 0 h34

 h14   h24  . h34  h44

EP

792

h12 h22 h23 h24

[1] R. Kobayashi, Modeling and numerical simulations of dendritic crystal growth, Physica D: Nonlinear Phenomena 63 (3) (1993) 410 - 423. [2] M. Berghoff, B. Nestler, Scale-bridging phase-field simulations of microstructure responses on nucleation in metals and colloids, The European Physical Journal Special Topics 223 (3) (2014) 409–419. [3] S. G. Kim, W. T. Kim, T. Suzuki, Phase-field model for binary alloys, Phys. Rev. E 60 (1999) 7186–7197. [4] M. Plapp, Unified derivation of phase-field models for alloy solidification from a grand-potential functional, Phys. Rev. E 84 (2011) 031601. [5] A. Choudhury, B. Nestler, Grand-potential formulation for multicomponent phase transformations combined with thin-interface asymptotics of the double-obstacle potential, Phys. Rev. E 85 (2012) 021602. [6] P. Steinmetz, M. Kellner, J. H¨otzer, A. Dennstedt, B. Nestler, Phase-field study of the pattern formation in al–ag–cu under the influence of the melt concentration, Comp. Mat. Science 121 (2016) 6 – 13. [7] M. B. Said, M. Selzer, B. Nestler, D. Braun, C. Greiner, H. Garcke, A phase field approach for wetting phenomena of multiphase droplets on solid surfaces, Langmuir 30 (14) (2014) 4033–4039. [8] J. H¨otzer, V. Rehn, W. Rheinheimer, M. J. Hoffmann, B. Nestler, Phasefield study of pore-grain boundary interaction, Journal of the Ceramic Society of Japan 124 (4) (2016) 329–339. [9] D. Schneider, E. Schoof, Y. Huang, M. Selzer, B. Nestler, Phase-field modeling of crack propagation in multiphase systems, Computer Methods in Applied Mechanics and Engineering 312 (2016) 186-195.

AC C

791

h34 = −2νn x ny


tρa − n2x ), h14 = n x ny 1 − 2ρ2a n2y − L + 4n2x (L − 1) , =2ν( A
ρb =2ν( − n2y ), h24 = n x ny 1 − 2ρ2b n2x − L + 4n2y (L − 1) A

 h11 Gρa ρb h12  H(x, y) = 1 − ν h13 h14 790

SC

tρa ρb
= 1 + ρ2a (1 + ) − A A
− n2x 2(1 − (2 − ρ2a )n2y ) + ρ2a + (3 − 4n2x )L , tρa
ρb )− = 1 + ρ2b (1 + A A
− n2y 2(1 − (2 − ρ2b )n2x + ρ2b + (3 − 4n2y )L ,

M AN U

h22

789

2

we express hi j (x, y) as

h11

788

2

TE D

787

2

¨ [10] W. Voigt, Uber die Beziehung zwischen den beiden Elastizit¨atskonstanten isotroper K¨orper, Annalen der Physik 274 (12) (1889) 573–587. [11] R. Spatschek, C. M¨uller-Gugenberger, E. Brener, B. Nestler, Phase field modeling of fracture and stress-induced phase transitions, Physical Review E 75 (6) (2007) 1–14. [12] C. Mennerich, F. Wendler, M. Jainta, B. Nestler, A phase-field model for the magnetic shape memory effect, Arch. Mech. 63 (2011) 549–571. [13] I. Steinbach, M. Apel, Multi phase field model for solid state transformation with elastic strain, Physica D 217 (2006) 153–160. [14] A. Reuss, Berechnung der Flie¨ssgrenze von Mischkristallen auf Grund der Plastizit¨atsbedingung f¨ur Einkristalle., Z. Angew. Math. Mech. 9 (1929) 49–58. [15] A. Khachaturyan, Theory of Structural Transformation in Solids, John Wiley & Sons Inc, 1983. [16] L.-Q. Chen, Phase-field models for microstructure evolution, Annual Review of Materials Research 32 (1) (2002) 113–140. [17] A. Durga, P. Wollants, N. Moelans, Evaluation of interfacial excess contributions in different phase-field models for elastically inhomogeneous systems, Modelling and Simulation in Materials Science and Engineering 21 (5) (2013) 55018. [18] D. Schneider, O. Tschukin, A. Choudhury, M. Selzer, T. B¨ohlke, B. Nestler, Phase-field elasticity model based on mechanical jump conditions, Computational Mechanics 55 (5) (2015) 887–901. [19] E. M. Lifschitz, L. P. Pitaewskii, Statistical Physics, Part 2: Theory of the Condensed State Vol. 9 (1st ed.), Butterworth-Heinemann, 1980. [20] J. W. Cahn, J. E. Hilliard, Free Energy of a Nonuniform System. I. Interfacial Free Energy, Journal of Chemical Physics 28 (1958) 258–267. [21] B. Nestler, H. Garcke, B. Stinner, Multicomponent alloy solidification: Phase-field modeling and simulations, Phys. Rev. E 71 (2005) 041609. [22] H. Garcke, B. Nestler, B. Stoth, On anisotropic order parameter models for multi-phase systems and their sharp interface limits, Physica D: Nonlinear Phenomena 115 (1) (1998) 87 – 108. [23] J. E. Taylor, J. W. Cahn, Diffuse interfaces with sharp corners and facets: Phase field models with strongly anisotropic surfaces, Physica D: Nonlinear Phenomena 112 (3) (1998) 381 – 411. [24] J. Eiken, B. B¨ottger, I. Steinbach, Multiphase-field approach for multicomponent alloys with extrapolation scheme for numerical application, Phys. Rev. E 73 (2006) 066122. [25] J. H¨otzer, M. Jainta, P. Steinmetz, B. Nestler, A. Dennstedt, A. Genau, M. Bauer, H. K¨ostler, U. R¨ude, Large scale phase-field simulations of directional ternary eutectic solidification, Acta Materialia 93 (2015) 194– 204. [26] M. Ohno, T. Takaki, Y. Shibuta, Variational formulation and numerical accuracy of a quantitative phase-field model for binary alloy solidification with two-sided diffusion, Phys. Rev. E 93 (2016) 012802. [27] A. Karma, Phase-field formulation for quantitative modeling of alloy solidification, Phys. Rev. Lett. 87 (2001) 115701. [28] J. Mosler, O. Shchyglo, H. Montazer Hojjat, A novel homogenization method for phase field approaches based on partial rank-one relaxation, Journal of the Mechanics and Physics of Solids 68 (2014) 251–266. [29] M. Nicoli, M. Plapp, H. Henry, Tensorial mobilities for accurate solution of transport problems in models with diffuse interfaces, Phys. Rev. E 84 (2011) 046707. [30] J. Ettrich, A. Choudhury, O. Tschukin, E. Schoof, A. August, B. Nestler, Modelling of transient heat conduction with diffuse interface methods, Modelling and Simulation in Materials Science and Engineering 22 (2014) 085006. [31] W. C. Johnson, J. I. D. Alexander, Interfacial conditions for thermomechanical equilibrium in two-phase crystals, J. Appl. Phys. 59 (1986) 2735–2746. [32] G. Caginalp, W. Xie, Phase-field and sharp-interface alloy models, Phys. Rev. E 48 (1993) 1897. [33] J. D. Eshelby, The elastic field outside an ellipsoidal inclusion, Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 252 (1271) (1959) 561–569. [34] X. Jin, Z. Wang, Q. Zhou, L. M. Keer, Q. Wang, On the solution of an elliptical inhomogeneity in plane elasticity by the equivalent inclusion method, Journal of Elasticity 114 (2014) 1–18.

782

RI PT

783

to define A(x, y) = tρa (x, y) + ρb (x, y). Thus, for the imagi-813 x˜a2 y˜ 2a nary ellipse 1+R(x,y) + t2 +R(x,y) = 1 the normal components are expressed as

12

814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879