8. Phasefield Modeling of Martensitic Laminates
8.1. Phase Field Approximation of Sharp Interfaces
105
X X
p(X) p(X)
1 1
l l
a) b)
Figure 8.1: Sharp and regularized order parameter distributions p(X) for a one-dimensional two-phasic bar with interface atX= 0. a) Exact but discontinuous description of the phase distribution byp(X) b) Approximate but smooth approximate description of the phase distribution byp(X) using of the artificial regularization length scalel.
phase field. We see from (8.1) that sharp phase boundaries Γ⊂ B lead todiscontinuities [[p(X)]] 6= 0 that are difficult to treat analytically and numerically. The idea is thus to approximate the non-smooth phase field distribution given by (8.1) by use of the smooth function
p(X) = 12
tanh (X/l) + 1
, (8.2)
see also McFadden et al. [107]. This function smears out the interface over the ax-ial domain R of the bar, representing a regularized interface description as depicted in Figure 8.1b. The width of the regularization is governed by the (artificially introduced) length scale parameter l > 0. The discontinuous function (8.1) is recovered in the limit l→0. The regularization (8.2) depicted in Figure 8.1b obviously has the properties
p′(−∞) = 0, p(0) = 12 and p′(+∞) = 0. (8.3) In particular, the curvature of p(X) as given in (8.2) is zero for p ∈ {0,1/2,1}. As a consequence, we find that (8.2) is a solution of the second order differential equation
p′′ = 4
l2 p(1−2p)(1−p) inB, (8.4)
subject to the Dirichlet and Neumann type boundary conditions (8.3). Equation (8.4) is the Euler equation of the variational principle
p= argn
pinf∈WpIB(p)o
with W = p
p(0) = 1/2 , (8.5)
where the functionalIB(p) can easily be constructed by integrating a Galerkin-type weak form of the differential equation (8.4) and reads
IB(p) = Z
B
n
2p2(1−p)2+ l2
2(p′)2o dV =
Z +∞
−∞
n
2p2(1−p)2+l2
2(p′)2o
A dX . (8.6) Note that sinceWdoes not specify Dirichlet boundary conditions forpon∂B, (8.5) implic-itly enforces the “natural” Neumann type boundary conditionsp′ = 0 on ∂B translating top′(−∞) = 0 and p′(+∞) = 0 for the given example. In summary we have
1 2
tanh (X/l) + 1
= argn
pinf∈WpIB(p)o
⇒ inf
p∈WpIB(p) = 13A l , (8.7) where (8.7)2 is obtained by insertingp(X) = 12[tanh(X/l) + 1] intoIB(p). Note that since IB(p) as defined in (8.6) is invariant with respect to the transformation p→ 1−p, (8.5)
8.1 Phase Field Approximation of Sharp Interfaces 107
also admits the alternative solution ˜p(X) = 1−p(X) = 12[1−tanh(X/l)] representing an inverse distribution of the “−” and “+” phases. From (8.7)2 we see that we can introduce a modified functional
ΓlB(p) = 3
lIB(p) = Z
B
n6
lp2(1−p)2+3l
2(p′)2o dV =
Z +∞
−∞
n6
lp2(1−p)2+3l
2(p′)2o A dX
(8.8) instead of (8.6). Clearly, the minimization of this scaled functional also yields the regu-larized interface description (8.2) depicted in Figure 8.1b. However, recalling (8.7), we see that the scaling of IB(p) by the factor 3/l has the consequence that the minimum value of the functional ΓlB(p) now directly yields the interface surface area,
pinf∈WpΓlB(p) = surf(Γ) =A (8.9) Note that the fact that ΓlB(p) returns the interface surface surf(Γ) = A for arbitrary regularization length scales l stems from the fact that we are considering an infinite bar.
This relation is best understood by considering a finite domain ˜B = A×[−L, L]. The functional Γl˜
B(p) for this finite domain then reads ΓlB˜(p) =
Z
B˜
n6
lp2(1−p)2+ 3l
2(p′)2o dV =
Z +L
−L
n6
lp2(1−p)2+ 3l
2(p′)2o
A dX . (8.10) Insertion of p(X) = 12[tanh(X/l) + 1] into (8.10) yields the expression
ΓlB˜(p) =A tanh(L/l)
1 + 12 sech2(L/l)
. (8.11)
from which we see that the exact interface surface area surf(Γ) = A is only recovered in the limit L/l →+∞ which is obtained either for an infinite domain L→ ∞ (see above) or independent of the domain size for a vanishing regularization length scale l →0,
liml→0 inf
p∈WpΓlB˜(p) = surf(Γ) =A . (8.12) Note that since for l → 0 we also recover the exact phase distribution (8.1) from the approximate solutionp(X) = 12[tanh(X/l)+1], we can say that the approximate problem converges to the sharp problem for l→0.
8.1.2. Approximation of Two- and Three-Dimensional Sharp Interfaces
The variational approximation of sharp interface geometries derived in the one-dimensional example above can be extended to multi-dimensional (two- and three-dimensional) solids in a straightforward manner. Let B ⊂ Rd, be the reference configuration of a material body with dimension d ∈ {2,3} in space and ∂B ⊂ Rd−1 its surface with normal N as depicted in Figure 8.2a. The domain B is separated into subdomains B− and B+ by a non-material sharp interface Γt. Again, to describe this distribution of phases more conveniently, we introduce the order parameter p(X, t)∈[0,1], wherep= 0 is associated with the “−” phase andp= 1 is associated with the “+” phase of the material, such that with p(X, t) = 0 for X ∈ B− and p(X, t) = 1 for X ∈ B+. Since our goal will be to to study the referential motion of non-material interfaces Γt in the range T ⊂ R+ of time, we introduce the time-dependent phase field
p:
B ×T →[0,1]
(X, t)7→p(X, t) , (8.13)
B− B−
B+ B+
p= 0 p= 0
p= 1 p= 1
Γt p=12
∂B
∂B
2l
N N
a) b)
Figure 8.2: Sharp and regularized order parameter distributions p(X, t) for a two-phasic bodyB with interface Γt. a) Exact but discontinuous description of the phase distribution byp(X, t). b) Approximate but smooth approximate description of the phase distribution byp(X, t) using of the artificially introduced regularization length scalel.
defined in the reference domainB. Then, a multi-dimensional extension of the regularized interface functional (8.8) can be expressed in the form
ΓlB(p) = Z
B
γl(p,∇Xp)dV , (8.14)
where we have introduced the interface surface density function per unit volume γl(p,∇Xp) = 6
lp2(1−p)2+3l
2k∇Xpk2. (8.15)
This function, locally depending both on the phase field p and its material gradient
∇Xp, will be very useful in the subsequent modeling of interface formation and evolution.
Given asharp interface-surface geometry Γt ⊂ Rd−1 in the solid B at time t as depicted in Figure 8.2a, we obtain its regularized phase field representationp(X, t) onBvisualized in Figure 8.2b in analogy to (8.5) from the minimization principle
p(X, t) = arg
pinf∈WpΓlB(p)
with Wp = p
p= 12 forX ∈Γt . (8.16) The Euler equations of the variational principle (8.16) are obtained from the necessary condition δΓlB= 0 and turn out to be a straight-forward generalization of (8.4),
4p(1−p) (1−2p)−l2∆Xp= 0 inB and ∇Xp·N = 0 on∂B, (8.17) where ∆Xp is the material Laplacian of the phase field and N is the outward normal on ∂B, see Figure 8.2b. Note that similar to (8.5), the pure topology problem (8.16) does not yield a unique solution but is always also minimized by ˜p(X, t) = 1−p(X, t) if p(X, t) is a minimizer due to the invariance of the functional (8.14) with respect to the transformation p → 1−p. Since the one-dimensional bar considered in Section 8.1.1 is a special case of the multi-dimensional theory treated in this section, it is not surprising that the approximation property (8.12) can be generalized to two or three dimensions by
liml→0 inf
p∈WpΓlB(p) = surf(Γt), (8.18) which will be the basis for our modeling of interface energy in Section 8.2. We have seen in the one-dimensional example in Section 8.1.1 that for l → 0 the regularized function
8.1 Phase Field Approximation of Sharp Interfaces 109
converges to the sharp function and the regularized surface converges to the sharp surface and are now stating the same for the multi-dimensional case. Indeed there is a much stronger notion of convergence that rigorously justifies these convergence properties: This is the Γ-convergence of the regularized problem to the related sharp problem, see Dal Maso [40] andBraides [26] for a general treatment of Γ-convergence andModica [119]
and Alberti et al. [6] for specifically related results.
8.1.3. Finite Element Formulation
To numerically verify the approximation property (8.18) and to gain some understanding regarding the approximation accuracy, we carry out a finite element discretization of the variational problem (8.16). To this end, we make use of the framework outlined in Section 4.4 by eliminating ϕt and identifying Iwith p and hence ∇XI with ∇Xp. In analogy to equations (4.67) and (4.68) we then have
γl(p,∇Xp) = γl(C) with C={p,∇Xp} and Ch ={p,∇Xp}h =B(X)d, (8.19) where d∈ RNI is the global nodal state vector (which follows from a reduction of (4.64)) that contains the values pI of the phase fieldp(X, t) at allNI global nodesXI at timet,
d=
p1, .., pI, .., pNI
T
, (8.20)
and where B(X) ∈ R(d+1)×NI is the generalized discrete gradient matrix (which follows from a reduction of (4.69)) given by
B(X) =
N1(X) N1,1(X) . . . N1,d(X) ...
NI(X) NI,1(X) . . . NI,d(X) ...
NNI(X) NNI,1(X) . . . NNI,d(X)
T
, (8.21)
where NI,A(X) stands for ∂XANI(X) =∇XNI(X)·EA, see (4.62). With (8.19), (8.20) and (8.21) at hand, we can now write down the space-discrete version of (8.14) as
ΓlBh(d) = Z
Bh
γl(Bd)dV . (8.22)
The discretization (8.22) allows the reformulation of the continuous variational principle (8.16) as a minimization of a function with respect to an NI-dimensional argument
d= arg
dinf∈WdΓlBh(d)
with Wd= d
Np(XJ)d= 12 forXJ ∈Γht , (8.23) where Np(X)∈ RNI is the global interpolation matrix for p (a reduction of (4.66)),
Np(X) =
N1(X) · · · NI(X) · · · NNI(X)
. (8.24)
As noted in Section 4.4.4, the consistent algorithmic solution of (8.23) follows as dk+1 =dk−AˆT
hAˆ ∂2ddΓlBh(dk) ˆAT i−1h
Aˆ ∂dΓlBh(dk)i
, dk ∈ Wd, (8.25)
Γt
B+
B−
∂B
0.5
0.5 0.5
0.5
0.5
Figure 8.3: Two-dimensional domain with circular interface surface Γtand surf(Γt) =π/2.
l/h infp∈WpΓl Bh/surf(Γt)
1.2 1.0 0.8 0.6 0.4 0.2
0.00 1 2 3 4 5 6
h= 0.01000 h= 0.00500 h= 0.00378 h≤l/2
Figure 8.4: Minimum mesh sizeh of bilinear finite elements for an approximation of an interface regularization with length scalel. A reasonable accuracy of the approximation of the sharp interface geometry infp∈WpΓlBh≈surf(Γt) requires an element size ofh < l/2.
where ˆA is the matrix that allows the extraction of the unconstrained nodal degrees of freedom ˆd from d by ˆd = ˆAd, see (4.82). For completeness, we note that the two derivatives of ΓlBh(d) appearing in (8.25) are obtained from numerical integration of
∂dΓlBh(d) = Z
Bh
BT
∂Cγl(B d)
dV , (8.26)
∂2ddΓlBh(d) = Z
Bh
BT
∂CC2 γl(B d)
BdV , (8.27)
which can be interpreted as integrals over the generalized stress (8.26) and the generalized tangent (8.27). Note that due to the non-convexity of γl(p,∇Xp) in p and due to the invariance with respect to the transformation p → 1− p, the algorithmic solution of (8.25) will generally depend on the choice of the initial guess dk=1.
8.1.4. Numerical Example of Sharp Interface Approximation
The following numerical example demonstrates the regularization of a sharp interface geometry for an elementary model problem. We consider a body of dimensiond= 2 with a sharp interface surface Γt that separates a circular inclusion B+ from the rest of the
8.1 Phase Field Approximation of Sharp Interfaces 111
a) b)
c) d)
0 p 1
Figure 8.5: Regularized interface surfaces described by the phase fieldp(X) for different length scales l. a) l = 0.06 yielding infp∈WpΓlBh = 1.006249 surf(Γ), b) l = 0.03 yielding infp∈WpΓlBh = 1.000996 surf(Γ), c)l = 0.015 yielding infp∈WpΓlBh = 0.998266 surf(Γ) and d)l= 0.0075 yielding infp∈WpΓlBh = 0.991357 surf(Γ) all obtained withh= 0.00378.
body as depicted in Figure 8.3. On the exterior boundary ∂B the Neumann condition (8.17)2 is prescribed. On the interface surface Γtthe Dirichlet conditionp= 12 is enforced for the phase field. We consider the finite element computation of the phase field p in the domain B according to the solution procedure (8.25) for given values of the length scale parameter l. The finite element mesh requires a certain minimum element sizeh in order to resolve this length scale. This is demonstrated in Figure 8.4 for finite element meshes consisting of bilinear quadrilateral elements with different element sizesh. In this elementary study, the ratio infp∈WpΓlBh/surf(Γt) of the finite element approximation of theregularized interface surface infp∈WpΓlBhand thesharp interface surface surf(Γ) =π/2 is plotted over the ratio l/h of thelength scale land themesh size hof the three different underlying meshes. The results show that an element size
h≤l/2 (8.28)
is the minimum resolution needed to approximate surf(Γ) by infp∈WpΓlBh with reasonable accuracy. Clearly however, this fine resolution is only needed in subdomains close to the interface surface. Hence, an h-adaptive finite element solution procedure with a minimum element size close to (8.28) could be applied. Figure 8.5 shows the finite element solution of the phase field p for the given boundary value problem with length scales l = 0.0075,0.015,0.03,0.06 based on a very fine mesh with constant mesh sizeh = 0.00378, consisting of 70,500 four node quadrilaterals. The largest considered length scalel = 0.06 yields the approximated regularized interface surface infp∈WpΓlBh = 1.006249 Γ depicted in Figure 8.5a. The smallest considered length scale l = 0.0075 gives infp∈WpΓlBh = 0.991357 Γ as shown in Figure 8.5d.
B− B−
B+ B+
p= 0 p= 0
p= 1 p= 1
Γt p=12
∂B
∂B
∇Xp yt yt
W W
VM VM
a) b)
Figure 8.6: Relation between the referential motion of the sharp interface Γt and the evolution of the order parameter distributionsp(X, t). a) A pointyton the sharp interface Γtmoves with velocityW b) The interface Γtis recovered from the approximate phase field description as the set of points wherep= 1/2, its normalM is given by−∇Xp/k∇Xpk.
8.1.5. Sharp Interface Kinematics and Phase Field Evolution
Having described a variational concept that links sharp interface geometries Γtto a smooth phase field approximationp(X, t), we now briefly describe how the motion of such a non-material sharp interface Γt with reference velocity W (see also Section 5.1.1.2) is linked to the time evolution ˙p(X, t) of the phase field. To this end, we consider a given phase field approximationp(X, t) of a sharp interface Γt at time t∈T that satisfies (8.16) and note that the sharp interface Γt can be recovered from p(X, t) as the level set
Γt=
X ∈ B
p(X, t) = 1/2 , (8.29)
see e.g. Hou et al. [75]. To understand how a referential motion of Γt is related to the evolution of the phase field p(X, t), recall Section 5.1.1 and consider the motion of a point yt(Θ) on the non-material interface Γt for a fixed Θ ∈ A. Since such a point (by definition) always lies on the interface, we have with (8.29) that
p yt(Θ), t
= 1/2 ∀t ∈T . (8.30)
Taking the total time derivative of (8.30) and making use of (5.8), we obtain the relation d
dtp yt(Θ), t
= 0 ⇔ p˙+∇Xp·W = 0 ∀t ∈T , (8.31) whereW = ˙yt(Θ) is the referential velocity of the non-material interface Γt. Noting that
∇Xp is always perpendicular to the level sets defined in (8.29), we can also recover the interface normalM of Γt from the phase field description p(X, t) as
M =−∇Xp/k∇Xpk ⇔ ∇Xp=−k∇XpkM, (8.32) where the minus in (8.32)1 stems from the fact that we have introducedM as the normal pointing towardsB+ in (5.2) and from the association of B+ with p= 1 in Sections 8.1.1 and 8.1.2, such that at the interface given byp= 1/2, the gradient ∇Xpwill point in the direction of −M. Finally, insertion of (8.32)2 into (8.31)2 yields the relation
˙
p− k∇XpkV = 0 where V =W ·M , (8.33)