Phase-inherent linear visco-elasticity model for infinitesimal deformations in the multiphase-field context

A linear visco-elasticity ansatz for the multiphase-field method is introduced in the form of a Maxwell-Wiechert model. The implementation follows the idea of solving the mechanical jump conditions in the diffuse interface regions, hence the continuous traction condition and Hadamard’s compatibility condition, respectively. This makes strains and stresses available in their phase-inherent form (e.g. \(\varepsilon ^{\alpha }_{ij}\), \(\varepsilon ^{\beta }_{ij}\)), which conveniently allows to model material behaviour for each phase separately on the basis of these quantities. In the case of the Maxwell-Wiechert model this means the introduction of phase-inherent viscous strains. After giving details about the implementation, the results of the model presented are compared to a conventional Voigt/Taylor approach for the linear visco-elasticity model and both are evaluated against analytical and sharp-interface solutions in different simulation setups.

The phase-field method (PFM) is an established simulation technique in science and engineering due to its capability to predict free boundary movement numerically in an efficient way. The key is to replace free boundaries between (physically) separable regions—or phases—by a transition region—or diffuse interface—which allows to track the position and orientation of the free boundary implicitly. In the context of microstructure evolution this reflects e.g. in the movement of grain boundaries driven by the influence of different thermodynamical fields. Typical applications are solidification, solid-state phase transformations, precipitate growth and coarsening, martensitic transformations and grain growth [1,2,3,4]. Further fields of application include topology optimisation [5] and crack propagation [6]. Apart from the capability of capturing free boundary movements, another asset of the PFM is the ability to represent complex structures and geometries without requiring a complex mesh [7]. Within the scope of the PFM the geometry is described by the phases and their diffuse interfaces, and not by a conforming mesh.

By extending the scope of materials and processes examined, different material behaviours may be observed. One of the four basic material behaviours, besides elasticity, plasticity and visco-plasticity, is visco-elasticity [8], which may describe the behaviour of amorphous and semi-crystalline polymers or metals at high temperatures, amongst others. In the context of solid mechanics on the basis of the PFM the predominantly used physically non-linear material behaviours are plasticity (e.g. Guo et al. [9], Ammar et al. [10], Schneider et al. [11] and Herrmann et al. [12], amongst others) and visco-plasticity (e.g. de Rancourt et al. [13], amongst others). Visco-elasticity is nearly exclusively applied to crack propagation simulations in polymeric and rubbery materials (e.g. Liu et al. [14], Shen et al. [15] and Yin et al. [16], amongst others).

In the aforementioned cases of crack propagation, the material model implementation is restricted to a single phase, as common phase-field crack models only consist of two phase formulations: a crack and a solid phase. Therefore these implementations do not have to deal with phase-field-specific questions as how to model the visco-elastic material behaviour at diffuse solid-solid interfaces.

A common approach to treat the material behaviour at diffuse interfaces is to use classical homogenisation schemes, namely those introduced by Voigt/Taylor [17, 18] (VT) or Reuss/Sachs [19, 20] (RS). As they assume locally uniform strain or stress in multiphase domains, respectively, they lead only in (quasi) 1D cases to physically exact behaviour on a continuum scale. In all other scenarios they only render the upper and lower bound, respectively, of the actual behaviour, because they do not consider shape and orientation of the phases. However, due to their simple implementation, and the argument in the limiting case of zero interface thickness they lead to a sharp-interface (SI) behaviour, they are applied regularly. A slightly different procedure is proposed by Khachaturyan [21], who suggests to use the VT scheme for elastic strains, and to treat inelastic strains by the RS scheme. More recently, the Hashin-Shtrikman bounds [22] (HS, which are second-order bounds while VT and RS are first-order bounds) have been applied by Chen and Shu [23] to simulate martensitic phase transformations. Durga et al. [24] examined how VT and RS schemes lead to excess stress, strain and elastic energy in diffuse interfaces. Their results mainly confirm that VT and RS homogenisations are only free of excess energy in special 1D cases, i.e. a straight interface does not introduce excess quantities if it is purely uniaxially loaded in directions tangential or normal to the interface, respectively. A more wholistic approach on the continuum scale is to consider mechanical jump conditions within the diffuse interface region to ensure kinematic compatibility and continuous traction. A model for two phases and finite deformations was introduced by Mosler et al. [25]. Schneider et al. [26] formulated a model for infinitesimal deformations and two phases, which assembles a locally compatible stiffness matrix based on mechanical jump conditions. Furthermore, by using their homogenisation scheme, they introduced a driving force resembling configurational forces (see, e.g. Gurtin [27]), which could match an analytical solution of the Gibbs-Thomson effect for a spherical inclusion. At the same time, but independently, Durga et al. [28] used their previous work to derive a similar model for the two-component, two-phase chemo-mechanical case, but without examination of a driving force. The work of Mosler et al. is continued in Kiefer et al. [29] and Bartels and Mosler [30] by extension to configurational and micro forces, accompanied by comparisons to VT and RS homogenisations and convergence studies. The basis remains two phase systems and finite deformations. Concurrently Schneider et al. [31] proposed a procedure for multiphase domains and finite deformations, which uses a locally implicit scheme to solve for both mechanical jump conditions. Again, the work also covers the derivation and evaluation of a consistent driving force. Schneider et al. [32] furthermore generalised their ansatz from Schneider et al. [26] to the multiphase-field method by using a generalised interface normal vector in multiphase interface regions. Svendsen et al. [33] extended the finite deformation ansatz to model multicomponent, multiphase systems. Herrmann et al. [12] used the procedure introduced by Schneider et al. [31] to implement small strain elasto-plasticity. They showed that within this framework it is possible to completely separate the material behaviour in diffuse interfaces and therefore to assign each phase its own distinct material model. This is possible because these advanced approaches are based on calculating stresses and strains phase-inherently. Lately, Henning et al. [34] used an approach based on Schneider et al. [26] and Durga et al. [28] to study material interfaces and point out that convergence and numerical stability may be improved compared to pure sharp-interface simulations.

To further show the versatile applicability of this ansatz of separable material behaviour, linear visco-elasticity in the form of the Maxwell-Wiechert model is considered. Details of the derivation of the theory and simulations are provided in Schwab [35]. The visco-elastic model is implemented on a multiphase-field basis and for infinitesimal deformations and is evaluated in a set of scenarios with a varying number of phases. The corresponding driving force may be calculated following the approach of Schneider et al. [26, 31], but this is not in the scope of this work.

In the following, the PFM is described and the approach modelling linear visco-elasticity is introduced. Then the PFM-based approach is applied to various simulation setups and the results are compared and discussed before conclusions are drawn.

The basis for this work is the multiphase-field model as formulated by Nestler et al. [36], and the continuum mechanics models considering jump conditions as introduced by Schneider et al. [31] and Herrmann et al. [12]. Hence, only relevant parts of these underlying works are recalled and in the following the focus lies on the new models implemented. As the mechanics are modelled in an infinitesimal deformation framework, no distinction is made between reference and current configuration.

In the following, the notation uses Einstein’s summation convention [37]. A fixed Cartesian coordinate system exists, and hence only subscripts appear in the notation (spatial indices are \(i = \left\{ 1, 2, 3 \right\} \)). Therefore, there is no summation over superscripts, especially phases \(\alpha \), if not explicitly mentioned.

Multiphase-field model

The multiphase-field model by Nestler et al. [36] introduces an N-tuple of order parameters \(\phi ^{\alpha }\) with \(\alpha = \left\{ 1, \dots , \text {N} \right\} \), \(\text {N} \in \mathbb {N}\). These order parameters denote phase volume fractions, and hence the individual phases themselves. This allows different physical regions within a domain to be parametrised in space \(x_i\) and time t by such order parameters \(\phi ^{\alpha } = \phi ^{\alpha }\left( x_i, t\right) \). In their function as a volume fraction, the phases furthermore fulfil the condition \(0 \le \phi ^{\alpha } \le 1\) and the local constraint \(\sum _{\alpha } \phi ^{\alpha }\left( x_i, t\right) = 1\) is enforced. When two neighbouring physical regions meet, the phases form a diffuse interface, which indicates that the involved phases smoothly transit into each other. Accordingly the same happens at triple- or multi-lines, where several physical regions coexist. Here, the transition between phases is governed by a multi-obstacle potential, cf. Nestler et al. [36], which results in a sine-shaped profile across a two-phase interface region. In this context, for a phase \(\phi ^{\alpha }\), this is expressed by

where \(x_{\mathrm{n}}\) is the distance perpendicular to the interface’s centre line (the sharp-interface position) and \(\epsilon \) stretches or compresses the sine-profile in this direction. Hence \(\epsilon \) is a measure for the width of the diffuse interface. To identify the location of diffuse interfaces, the gradient of the involved phases may be calculated, denoted by \(\phi ^{\alpha }_{,i}\). Alike the orientation or normal vector of selfsame diffuse two-phase \(\alpha \!\beta \)-interface may be retrieved (see, e.g. Beckermann et al. [38]):

where \(\left| \phi ^{\alpha }_{,j} \right| = \sqrt{\phi ^{\alpha }_{,j} \phi ^{\alpha }_{,j}}\). For regions where more than two phases meet, generalised normal vectors may be constructed. For example, following Moelans et al. [39], Schneider et al. [31] use a generalised normal vector of the form

in their mechanics model.

The energy densities associated with the diffuse interfaces, \(\gamma \), and with the bulk regions, \(\psi \), are collected in a free energy functional f over the volume of interest v, which generally reads as

which may, indicated by the ellipses, incorporate further dependencies on physical quantities apart from the domain parametrisation by the phase-fields \(\phi ^{\alpha }\).

By taking the variational derivate of the free energy functional f an evolution equation for the phases \(\phi ^{\alpha }\), and hence the diffuse interfaces, may be constructed. A common approach is the Allen-Cahn equation [40, 41] in the form given by Steinbach and Pezzolla [42],

which allows for binary interactions of the locally involved phases. Here, \(M^{\alpha \!\beta }\) is a positive mobility parameter, which may be individual for each \(\alpha \!\beta \)-interface, \(\tilde{N}\) is the number of locally existing phases and \(\epsilon \) controls the width of the diffuse interfaces as specified by Eq. (1). The operator \({\delta }\left( \bullet \right) \) describes a variational derivative.

In this section the authors intentionally did not introduce and discuss a specific energy functional, because phase changes of any kind—as shortly noted in the introduction section—are not considered in this work. By only vaguely defining an energy functional the authors furthermore want to point out that their method may not only be applied to phase-field models (see, e.g. [12, 31]) or continuum mechanics based approaches (see, e.g. [35]) but to diffuse interface approaches in general.

Multiphase continuum mechanics with infinitesimal deformations and phase-inherent quantities

By assuming the interfaces to be singular surfaces as well as the per-volume bulk free energy to be a volume average, \(\psi =\sum _{\alpha } \phi ^{\alpha } \psi ^{\alpha }\), and following Schneider et al. [31] and Herrmann et al. [12], the mechanical quantities have to fulfil the following conditions:

(1):

Total quantities are a volume average of their phase-inherent counterparts, i.e.

(a):

the local total stress is \(\sigma _{ij} =\sum _{\alpha } \phi ^{\alpha } \sigma _{ij}^{\alpha }\).

(b):

the local total strain is \(\varepsilon _{ij} =\sum _{\alpha } \phi ^{\alpha } \varepsilon _{ij}^{\alpha }\).

(2):

In each diffuse \(\alpha \!\beta \)-interface region it must hold

(a):

Hadamard’s compatibility condition (no slip): \(\llbracket u_{i,j} \rrbracket ^{\alpha \!\beta }= a_{i}^{\alpha \!\beta }n_{j}^{\alpha \!\beta }\).

(b):

the continuous traction condition: \(\llbracket \sigma _{ij} \rrbracket ^{\alpha \!\beta }n_{j}^{\alpha \!\beta } = 0\).

In the above lines, \(\sigma _{ij}\) is a Cauchy stress tensor, \(\varepsilon _{ij}\) is an infinitesimal strain tensor and \(a_{i}^{\alpha \!\beta }\) is the jump vector at a specific \(\alpha \!\beta \)-interface. Furthermore, \(\llbracket \bullet \rrbracket ^{\alpha \!\beta }=\left( \bullet \right) ^{\alpha } - \left( \bullet \right) ^{\beta }\) are the jump brackets. When considering non-linearities, all four conditions may be brought in line by using a local Newton-Raphson scheme. The basic idea here is to solve a system of equations which accounts for all locally existing traction conditions (the residual), where each line is of the form

The involved phase-inherent stresses, \(\sigma _{ij}^{\alpha }\), are calculated via Hooke’s law,

where \(C_{ijkl}^{\alpha }\) is the stiffness tensor and \(\varepsilon _{kl}^{\alpha ,\mathrm{el}}\) is the elastic part of the phase-inherent infinitesimal total strain tensor. The phase-inherent infinitesimal (total) strain tensor \(\varepsilon _{ij}^{\alpha }\), in turn, is calculated from combining Hadamard’s compatibility condition and the volume average condition for the local total strain \(\varepsilon _{ij}\).

While based on a physically sound motivation, it unfortunately turns out that, in regions with more than two phases, the system above is usually overdetermined and it is generally not possible to satisfy all equations above simultaneously. In fact, considering for simplicity the purely elastic case, i.e. \(\varepsilon _{ij}^{\alpha ,\mathrm{el}} = \varepsilon _{ij}^{\alpha }\), the phase-inherent (PI) approach outlined above introduces \(\tilde{N}\) phase-inherent displacement gradients \(u_{i,j}^{\alpha }\) as well as \(\tilde{N} \left( \tilde{N} - 1 \right) \) jump vectors \(a_{i}^{\alpha \!\beta }\) as additional degrees of freedom. Using the condition on the volumetric average of the displacement gradient (cf. Item 1b), the phase-inherent quantity \(u_{i,j}^{\alpha }\) may be retrieved in terms of the volume-averaged displacement gradients and the jumps \(\llbracket u_{i,j} \rrbracket ^{\beta \!\alpha }\)

Assuming Hadamard’s compatibility condition (cf. Item 2a) to hold between all phases, this allows expressing the \(\tilde{N}\) phase-inherent displacement gradients as

and thus as a function of the average displacement gradient and the jump vector alone.

A priori, this leaves \(3 \tilde{N} \left( \tilde{N} - 1 \right) \) additional degrees of freedom for satisfying the same number of continuous traction conditions in Item 2b. As \(n_{j}^{\alpha \!\beta } = - n_{j}^{\beta \!\alpha }\) and \(\llbracket u_{i,j} \rrbracket ^{\alpha \!\beta }=a_{i}^{\alpha \!\beta }n_{j}^{\alpha \!\beta } = - \llbracket u_{i,j} \rrbracket ^{\beta \!\alpha } = -a_{i}^{\beta \!\alpha }n_{j}^{\beta \! \alpha } = a_{i}^{\beta \!\alpha }n_{j}^{\alpha \!\beta }\), one necessarily has \(a_{i}^{\beta \!\alpha } = a_{i}^{\alpha \!\beta }\), which allows reducing the number of unknown jump vectors to \(\nicefrac {\tilde{N} \left( \tilde{N} - 1 \right) }{2}\), given by \(a_{i}^{\alpha \!\beta }\) for \(1 \le \alpha < \beta \le \tilde{N}\). This does not impose any additional difficulties with regards to the solvability since, at the same time, with \(\llbracket \sigma _{ij} \rrbracket ^{\alpha \!\beta }n_{j}^{\alpha \!\beta } = - \llbracket \sigma _{ij} \rrbracket ^{\beta \!\alpha } n_{j}^{\alpha \!\beta } = \llbracket \sigma _{ij} \rrbracket ^{\beta \!\alpha } n_{j}^{\beta \!\alpha }\), one may simultaneously reduce the number of continuous traction conditions to be satisfied to those for \(1 \le \alpha < \beta \le \tilde{N}\). In regions with more than two phases though, the algebraic equalities

for \(\delta \ne \alpha , \beta \) impose an additional set of compatibility conditions on the jump vectors. This may significantly reduce the number of remaining degrees of freedom and in general one may therefore not satisfy all the continuous traction conditions simultaneously.

One way of partially circumventing this difficulty is to choose to enforce Hadamard’s compatibility and the continuous traction condition with respect to one reference phase R only. As Eq. (9) remains valid for the reference phase R, and so does the relation \(u_{i,j}^{\alpha } =u_{i,j}^{\mathrm{R}} + a_{i}^{\alpha \!\mathrm{R}} n_{j}^{\alpha \!\mathrm{R}}\) (as well as Eq. (8)) for all \(\alpha \ne \text {R}\), the phase-inherent strains are still completely fixed in terms of the \(\tilde{N} - 1\) jump vectors \(a_{i}^{\alpha \!\mathrm{R}}\), \(\alpha \ne R\), which may be determined by the corresponding number of conditions on the normal components of the stresses \(\llbracket \sigma _{ij} \rrbracket ^{\alpha \!\text {R}} n_{j}^{\alpha \!\mathrm{R}} = 0\). When solved numerically, the updates yield an adjusted, volume-averaged local total stress \(\sigma _{ij}\), which may be used to solve the static linear momentum balance, here without body forces, \(\sigma _{ij,j} = 0\). For more details the reader is referred to the aforementioned work of Schneider et al. [31] and Herrmann et al. [12].

Note that, since this way of restoring solvability is based on dropping some of the equations, the jump in the displacement gradient (cf. Eq. (10)) between two phases \(\alpha , \beta \ne \text {R}\) is of the form

and thus in general consists of a rank-two tensor which is incompatible with Hadamard’s compatibility condition. Similarly, in general, there is no reason for the continuous traction condition to hold for any such phase pairings. Nevertheless, for pure two-phase regions the overall approach is definite and for multiphase regions still sufficient, as the same problem would exist for multilines or -junctions in a sharp-interface context.

Phase-inherent linear visco-elasticity using the Maxwell-Wiechert model

In contrast to elasticity, visco-elastic material behaviour shows a rate-dependence, but both have no equilibrium hysteresis (see, e.g. Haupt [8]). Hence, Eq. (7), which formulates rate-independent, linear elastic behaviour, is reformulated and written in a more general form for a phase \(\alpha \) as

The equilibrium elastic behaviour is retrieved if \(\dot{\varepsilon }_{mn}^{\alpha } = 0\). A general visco-elasticity law covers creep and relaxation behaviour. The former describes a behaviour when strain increases towards its equilibrium under a constant stress load, and the latter happens under a constant strain load which causes the stress to decrease towards its equilibrium.

Such a general model capturing the behaviour outlined by Eq. (12) is the Maxwell-Wiechert model. As depicted in Fig. 1, the model splits the rate-dependent part of the equation, \(h_{ij}^{\alpha }\), into \(\text {P}^{\alpha }\) functions of the same type, so-called Maxwell elements. Each element describes an individual creep and relaxation process, often linked to a specific time scale. As the model allows to connect an arbitrary number of Maxwell elements in parallel, the model may be adjusted to any visco-elastic behaviour. Therefore it is also called Generalised Maxwell model.

Maxwell-Wiechert model: A rheological representation with \(m = \left\{ 1, \dots , \text {P}^{\alpha } \right\} \), \(\text {P}^{\alpha } \in \mathbb {N}\), Maxwell branches, where \(\sigma ^{\mathrm{ext}}\) is an outer stress load, \(C_{ijkl}\) is a stiffness tensor, \(V_{ijkl}\) is a viscosity tensor, \(\varepsilon _{kl}^{\mathrm{el}}\) is an elastic strain and \(\varepsilon _{kl}^{\mathrm{v}}\) is a viscous strain. A similar representation may be found in e.g. Haupt [8]

Full size image

The Maxwell-Wiechert model, cf. Fig. 1, is composed of springs, an idealisation of the time-independent elastic behaviour (stiffness tensor \(C_{ijkl}\), related quantities superscripted by “el” for elastic), and dampers, which render the time-dependent viscous behaviour (viscosity tensor \(V_{ijkl}\), related quantities superscripted by “v” for viscous). As all branches of the model are connected in parallel, the response of a material \(\alpha \) to a load with a total stress and strain is

Here, \(m = \left\{ 1, \dots , \text {P}^{\alpha } \right\} \), \(\text {P}^{\alpha } \in \mathbb {N}\), indicates a Maxwell element, \(\varepsilon _{ij}^{\alpha ,\mathrm{M},m}\) its strain and \(\sigma _{ij}^{\alpha ,\mathrm{M},m}\) its stress, respectively. Quantities representing the overall response of a Maxwell branch are also superscripted by a capital M. The strain of the zeroth branch, which consists only of a spring, is \(\varepsilon _{ij}^{\alpha ,0} =\varepsilon _{ij}^{\alpha ,\mathrm{el},0}\), and its stress is \(\sigma _{ij}^{\alpha ,0}\). Within a single Maxwell element spring and damper are connected in series, giving the relations

Using the general behaviour of a spring, \(\sigma _{ij} =C_{ijkl}\varepsilon _{kl}\), and a damper, \(\sigma _{ij} =V_{ijkl}\dot{\varepsilon }_{kl}\), the total stress is retrieved as

In the same way, by combining the stresses in a single Maxwell element, the evolution equation for the viscous strain of each Maxwell element is given by

A common simplification is to assume that viscous strains are incompressible, i.e. only the deviatoric part effects the behaviour \(\varepsilon _{st}^{\mathrm{v},\mathrm{dev}} =\varepsilon _{st}^{\mathrm{v}} -\varepsilon _{st}^{\mathrm{v}, \mathrm{sph}}\), with the spherical part being \(\varepsilon _{st}^{\mathrm{v}, \mathrm{sph}} = \nicefrac {1}{3} \varepsilon _{kk} \delta _{st}\), and that the springs and dampers act isotropically, resulting in the stiffness and the viscosity being describable by the bulk modulus K, the shear modulus G and the scalar viscosity \(\eta \). With this, Eqs. (15) and (16) change to

and

Of course, with these assumptions, the rheological model in Fig. 1 only represents the behaviour due to the deviatoric part of the strain. In such cases, the spherical part is now connected in series, as the above equations suggest.

Now, to change the material behaviour of a single phase within the continuum mechanics model derived in the previous subsection, Eq. (7) has to be changed to the stress-strain relation of the Maxwell-Wiechert model in either of its forms.

Implementation of the phase-inherent visco-elastic model

The introduction of phase-inherent strains and viscous strains in the form above leads to a model containing a relatively large number of additional unknowns as compared to more traditional models based on phase-averaged quantities only. After an implicit discretisation in time, these consist, in addition to the average displacement field \(u_i\), of the \(\tilde{N}-1\) jump vectors \(^{n+1}{}a_i^{\mathrm{R}\!\alpha }\) allowing for the determination of the \(\tilde{N}\) phase-inherent total strains \(^{n+1}{}\varepsilon _{ij}^{\alpha ,\mathrm{tot}}\) as well as the \(\sum _{\alpha =1}^{\tilde{N}} \text {P}^{\alpha }\) different phase-inherent viscous strains \(^{n+1}{} \varepsilon _{ij}^{\alpha ,\mathrm{v},m}\). These are accompanied by an equal number of additional conditions to be satisfied, namely the update rule (from Eq. (16) by temporal discretisation)

with

for the viscous strains and the continuous traction conditions

with respect to the reference phase R in the sense of Eq. (6). Inserting the expression for the new viscous strains and, for convenience, dropping the superscript \(\left( \bullet \right) ^\text {tot}\) for the total strain of a phase \(\alpha \) as well as dropping the superscript \(^{n+1}{}(\bullet )\) for the quantities at the new time step, the residual in Eq. (21) becomes

with \(E_{irst}^{\alpha ,m} = C_{iruv}^{\alpha , m} \left( \delta _{us} \delta _{vt} - D_{uvst}^{\alpha , m}\right) \). Unlike \(D_{uvst}^{\alpha ,m}\), \(E_{irst}^{\alpha , m}\) is again a symmetric tensor since the term \(C_{iruv}^{\alpha , m} D_{uvst}^{\alpha , m}\) is given by

Rewriting \(E_{irst}^{\alpha , m}\) based on this expression, and using \(C_{opst}^{\alpha ,m} = \frac{1}{\Delta t} V_{opst}^{\alpha ,m} +C_{opst}^{\alpha ,m} - \frac{1}{\Delta t} V_{opst}^{\alpha ,m}\), gives

From the first line it follows that \(E_{irst}^{\alpha , m}\) is in addition positive (semi-)definite provided \(V_{ijkl}^{\alpha ,m}\) and \(C_{ijkl}^{\alpha ,m}\) are also positive (semi-)definite. The examination of the second line for \(\Delta t \rightarrow \infty \) and \(\Delta t \rightarrow 0\) reveals that \(E_{irst}^{\alpha , m}\) tends towards zero and towards \(C_{irst}^{\alpha , m}\), respectively, which corresponds to the long-term and short-term behaviour of the linear visco-elastic material model.

Combined with the expressions

for the phase-inherent displacement gradients in terms of the jump vectors \(a_i^{\alpha \!\mathrm{R}}\) based on Eq. (9), Hadamard’s compatibility condition (cf. Item 2a) with respect to the reference phase and making use of the subsymmetries of the tensors used, a short calculation allows to transform Eq. (22) into the linear system

for the jump vectors, where \(\varepsilon _{st}\) is the volume-averaged strain tensor and

is the algorithmically consistent phase-inherent effective stiffness tensor (in terms of \(\Delta t\)), which, by the discussion of Eqs. (23) and (24), is symmetric and positive (semi-)definite.

A similar calculation may also be performed for the isotropic model depending on the deviatoric viscous strains described by Eqs. (17) and (18). With respect to the properties of \(\hat{C}_{irst}^{\alpha }\) it is important to stress that the use of strictly deviatoric viscous strains will only lead to a symmetric effective stiffness if the projector onto the deviatoric tensors commutes with \(E^{\alpha , m}_{irst}\). This is obviously the case if all materials are isotropic, and may also be fulfilled for other relatively simple anisotropies, but it is not generally true. As, for a complex anisotropy, the use of purely deviatoric viscous strains does not allow for any major simplifications in the algorithm above anyway, it may be better to avoid such an assumption in such cases.

Solving the system in Eq. (26) leads to the jump vectors as a function of the (given) phase-inherent viscous strains \(^{n}{}\varepsilon _{ij}^{\alpha ,\mathrm{v},m}\) from the previous time-step and the volume-averaged strain \(\varepsilon _{st}\). From these, one may first calculate the phase-inherent total strains \(\varepsilon _{st}^{\alpha }\), then the updated viscous strains \({\varepsilon }_{ij}^{\alpha ,\mathrm{v},m}\), and finally the volume-averaged stress \(\sigma _{ij}\) using the material law in Eq. (15) for the phase-inherent total stresses combined with the volume-averaging condition in Item 1a.

Note that the procedure above essentially just corresponds to a Schur complement approach based on an implicit block-elimination of all unknown (local) phase-inherent quantities \(\varepsilon _{ij}^{\alpha }\), \(a_{i}^{\mathrm{R}\!\alpha }\) and \(\varepsilon _{st}^{\alpha ,\mathrm{v},m}\) in terms of the volume-averaged strain and the known phase-inherent viscous strains \(^{n}{}\varepsilon _{ij}^{\alpha ,\mathrm{v},m}\) from the previous time step (see e.g. Vassilevski [43] for a detailed discussion of the Schur complements associated with block-eliminations). In terms of the global problem given by the static linear momentum balance, \(\sigma _{ij,j} = 0\), regarding the determination of the displacement field \(u_{i}\), the compatibility with any standard Krylov solution algorithm should be mentioned. This is due to the fact that, in addition to the initial residual calculation, they only require the ability to evaluate the action of the system based on increment-type auxiliary vectors \(\delta u_i\). The viscous strains \(^{n}{}\varepsilon _{ij}^{\alpha ,\mathrm{v},m}\), being fixed, need to be dropped from the update rules in Eqs. (19) and (26) in order to obtain the corresponding increments in the viscous strains and the jump vectors.

Alternatively, if the solution of the linear momentum balance requires more than a few iterations or more generally for preconditioning purposes, it is often computationally preferable to form the resulting Schur complement system explicitly. In fact, even though based on the solution of small local subsystems only, repeatedly applying the update procedure above may become expensive due to the potentially large number of phase-inherent quantities. As the phase-inherent quantities are, except for the finally converged values, only auxiliary quantities for evaluating the effective volume-averaged stresses in terms of the volume-averaged strains, their action may simply be replaced by the use of an average algorithmically consistent stiffness tensor \(\hat{C}_{ijkl}\). Abbreviating the matrix corresponding to the linear system in Eq. (26) as \(A_{ij}^{\alpha \!\beta }\),

the increments \(\Delta a_{i}^{\alpha \!\mathrm{R}}\) of the jump vectors corresponding to increment \(\Delta \varepsilon _{ij}\) of the volume-averaged strain are given by

Using

together with the definition of the volume-averaged stress in Item 1a and the phase-inherent effective stiffness tensors \(\hat{C}^{\alpha }_{ijkl}\) from Eq. (27), yields

With \(\sum _{\alpha } \phi ^{\alpha } = 1\) and exchanging the summation variable in the last term, this may be rewritten more compactly as

and thus finally, combined with the expression for the jump increments in Eq. (29), leads to the average effective stiffness

based on the phase-inherent stiffnesses \(\hat{C}^{\alpha }_{ijkl}\) in Eq. (27).

Despite its cumbersome explicit expression, this effective stiffness only needs to be precalculated and stored once at the beginning of each time step. After the calculation of the initial residual, all subsequent iterations may then be performed completely in terms of the volume-averaged strains and stresses (without any reference to phase-inherent variables) simply using the precalculated effective stiffness in Eq. (33). It is only after convergence on the global problem has been achieved that the jump vectors and the phase-inherent viscous strains need to be updated based on the converged strain \(\varepsilon _{ij}\).

Finally, it should be noted that - despite the incremental character in time of the viscous strains - both approaches introduced may in principle be implemented using only a single field for the storage of the phase-inherent viscous strains as the updates may always be generated locally based on the update rule, Eq. (19), and only need to be stored at the end of each time step. This reduction may be quite important for problems involving a larger number of phases and Maxwell elements as storing all the strains involved during such calculations would require a relatively large amount of memory.

Conventional approach on the basis of a VT homogenisation

For the sake of comparison, the previously introduced Maxwell-Wiechert model is also briefly noted using a VT approach. Here, the basic assumption is that there is no phase-inherent deformation (cf. Item 1b), hence \(\varepsilon _{ij} = \varepsilon _{ij}^{\alpha }\) and the same applies for the viscous strains. Thus, the material parameters are averaged over all phases (compared to averaging strains and stresses), and, for an isotropic material, this yields

With this, the total stress response in a local point (cf. Item 1a) is directly the stress-strain relation of the Maxwell-Wiechert model in its modified form,

and analogously the evolution equation of each single Maxwell element is averaged over all locally existing phases,

The linear visco-elasticity models presented in the previous section are implemented in the multiphysics and multiphase-field software environment Pace3d [44]. In the following, the mechanical model is solved implicitly on an equidistant finite element grid with linear elements and reduced integration. The spatial discretisation is \(\Delta x = \Delta y = \Delta z =1.0\,\upmu \hbox {m}\) in all simulations. Furthermore, only isotropic materials are considered and viscous behaviour under shear is assumed. The 2D simulations shown were conducted as 3D simulations with a thickness of \(1.0\,\upmu \hbox {m}\) discretised by one cell, and the boundary conditions in z-direction mimic a plane strain state. For visualisation, the software Visit [45, 46] is partly used.

Behaviour in a single material point

Initially, the basic implementation of the visco-elasticity equations is validated. This is done for a single material point with only one phase existent, hence whether the model is implemented phase-inherently or conventionally does not play any role. The usage of a standard linear solid material behaviour (see Fig. 2) allows the calculation of analytical solutions of the temporal behaviour (cf. Haupt [8], Kaliske and Rothert [47] or Careglio et al. [48]).

Standard linear solid: Maxwell-Wiechert model with one Maxwell element

Full size image

These are used in the following to compare the numerical solutions of the model under relaxation and creep conditions.

Relaxation

Using Fig. 2, the analytical solution for the relaxation of stress under a constant shear strain load reads as

For the comparison, a constant shear strain \(\varepsilon _{12} = 0.1\) is considered, and a time period of \(t = 20\,\hbox {s}\) is simulated by using \(\Delta t = 0.01\,\hbox {s}\). The simulation is performed with five different viscosities \(\eta ^{\alpha ,1}\), a complete list of material properties is found in Table 1.

As depicted in Fig. 3, the numerical implementation produces values which are nearly identical to the analytical solution.

Relaxation test, one phase: Comparison of numerical and analytical solution

Full size image

Creep

Similarly to the relaxation test, the analytical solution for the creep test under a constant shear stress load may be constructed from the rheological model in Fig. 2. This yields the equation

In this case a constant shear stress \(\sigma _{12} =50\,\hbox {MPa}\) is chosen. Again the simulation covers \(t = 20\,\hbox {s}\) with a temporal discretisation of \(\Delta t = 0.01\,\hbox {s}\), and the used material properties are composed in Table 2.

The creep behaviour simulated by the numerical implementation agrees well with the analytical solutions (see Fig. 4).

Creep test, one phase: Comparison of numerical and analytical solution

Full size image

Behaviour in two-phase interface regions

After verifying the material behaviour in a single material point, the two model implementations are used in three different scenarios of two-phase interface regions. These are a domain with two phases, connected by a straight interface, a circular inclusion phase within a matrix phase, hence connected by a curved interface, and a domain divided by an inclined interface, whose two subdomains, thus created, contain a circular inclusion each. Not all phases are modelled linear visco-elastically in these scenarios, but also purely linear elastic phases occur to better show the capability of the phase-inherent modelling approach. Sharp-interface solutions, computed by Pace3d, are given to better classify the results. If non-horizontal and non-vertical interfaces are part of a simulation domain, the sharp-interface simulation relies on Abaqus [49]. In Abaqus, for better comparability, linear elements with reduced integration are used as well. The element size of Abaqus simulations is leaned on the cell size in the Pace3d equidistant grid, but due to the complex interface geometries the element size is non-constant. In such cases, the element size of Abaqus simulations tends to be smaller than in Pace3d, hence small deviations in the solutions may be attributed to the different resolutions.

Two phases, straight interface

The basic setup for a \(100\,\upmu \hbox {m} \times 100\,\upmu \hbox {m}\) plate with two phases and an x-centred, straight interface is given in Fig. 5. The domain is loaded by outward and inward pointed orthogonal and constant displacement of \(\bar{u} = 0.05\,\upmu \hbox {m}\) for the x- and y-boundaries, respectively. This means the remaining displacement components on the boundaries are not set and hence free to evolve.

Setup: Two phases, purely linear elastic \(\alpha \) and linear visco-elastic \(\beta \), are connected by a diffuse or sharp interface

Full size image

The diffuse interface exhibits twelve points within the range \(0 \le \phi ^{\alpha } \le 1\), here giving a width of \(l=12\,\upmu \hbox {m}\). The system is relaxed over a time of \(t = 20\,\hbox {s}\) and the temporal discretisation is \(\Delta t = 0.01\,\hbox {s}\). In Table 3 the material properties for the purely linear elastic phase \(\alpha \) and the linear visco-elastic phase \(\beta \) are given.

The results are plotted along x for a constant \(y=49.5\,\upmu \hbox {m}\) and at times \(t = \left\{ 0.5, 1.5, 20.0\right\} \,\hbox {s}\). In Fig. 6, the conventional/VT and the PI approaches are compared to a sharp-interface solution.

Results for a \(\varepsilon _{11}\) and b \(\sigma _{11}\): Comparison of the VT and the PI approach in mechanics to the sharp-interface solution. Three different times (\(t = \left\{ 0.5, 1.5, 20.0\right\} \,\hbox {s}\)) are plotted along x with \(y=49.5\,\upmu \hbox {m}\)

Full size image

The phase-inherent ansatz is able to match the sharp-interface solution at all times, obviously for \(\varepsilon _{11}\) having its jump stretched over the diffuse interface region (Fig. 6a). Furthermore, the centre of the jump coincides with the centre of the interface. For the conventional VT ansatz the \(\varepsilon _{11}\) values are shifted towards the visco-elastic phase, and it is unable to match the bulk values. For \(\sigma _{11}\) (Fig. 6b) both approaches have a continuous stress field, but with the VT approach not matching the sharp-interface solution. In Fig. 7 the basic behaviour of the PI approach is illustrated.

Phase-inherent quantities in the interface region: Development of \(\varepsilon _{11}^{\alpha }\) and \(\varepsilon _{11}^{\beta }\) and the volume averaged quantity \(\varepsilon _{11}\) over x with \(y = 49.5\,\upmu \hbox {m}\)

Full size image

In this example it is shown that indeed two separate strains are existing within the binary diffuse interface region, one for each phase. As a sharp-interface \(\varepsilon _{11}\) exhibits a jump in the direction normal to the interface, the phase-inherent quantities \(\varepsilon _{11}^{\alpha }\) and \(\varepsilon _{11}^{\beta }\) do not meet but stay discontinuous. The volume-averaged \(\varepsilon _{11}\) forms a continuous field, which is nothing else than the jump distributed across the interface width.

For this setup the influence of the interface width, the grid resolution and material behaviour have been briefly examined as well.

The results for the interface widths \(l = \left\{ 6, 12, 22\right\} \,\upmu \hbox {m}\) are depicted in Fig. 8.

Influence of different interfaces widths: Results for \(\varepsilon _{11}\) and \(\sigma _{11}\) compared to a sharp-interface solution at \(t=20\,\hbox {s}\). The results of the VT approach are depicted in a and c, while the PI approach is shown in b and d. Three interface width of \(l=\left\{ 6, 12, 22\right\} \,\upmu \hbox {m}\) are plotted along x with \(y=49.5\,\upmu \hbox {m}\)

Full size image

The PI approach in Fig. 8b and d matches the sharp-interface solution for all examined interface widths. In case of \(\varepsilon _{11}\) it is stressed, that the (diffuse) jump of the quantity remains symmetric to the centre of the interface for all l. The VT approach (see Fig. 8a, c) shows a clear trend towards growing deviations in values with a growing interface width. Especially quantities which exhibit a jump across the interface, here \(\varepsilon _{11}\), also show a growing deviation in shape, i.e. the jump is not symmetric to the interface centre but shifted. This implicates if \(l \rightarrow 0\,\upmu \hbox {m}\) the deviations will fully disappear, but since the PFM method needs a certain number of grid points in the interface to work properly, the deviations of the VT approach will de facto never fully vanish.

In Fig. 9 the grid spacing \(\Delta x\) has been varied, while the number of grid points in the interface is kept constant.

Influence of different resolutions: Results for \(\varepsilon _{11}\) compared to a sharp-interface solution at \(t=20\,\hbox {s}\). Three resolutions of \(\Delta x =\left\{ 1.00, 0.50, 0.25\right\} \,\upmu \hbox {m}\) are plotted along a x with \(y = 49.5\,\upmu \hbox {m}\) and b the grid points in the interface (i.e. the interface is entered from the left at interface grid point 0 and left towards the right at interface grid point 12)

Full size image

As visible in Fig. 9a, the PI approach yields results in the bulk regions identical to the sharp-interface solution for all grid spacings examined. For finer resolutions, the physical interface width becomes smaller and hence the (diffuse) jump of \(\varepsilon _{11}\) is performed over a shorter distance. The VT scheme slowly approaches the sharp-interface solution in the bulk regions the finer the grid is resolved, and the jump of the quantity approaches a sharp jump but stays slightly shifted compared to the symmetric shape resulting from the PI approach. Examining the interface region not on the basis of the physical width, but from the number of grid points within the region (see Fig. 9b), it becomes obvious, that besides the gain in accuracy in the bulk regions, the VT approach stays nearly equivalently shifted and distorted for all grid spacings. Hence, even by using physically narrow interfaces, the driving force calculation and therefore the movement of the interface might stay equally defective (e.g. distortion of the interface, excess energy) as with physically or computationally broad interfaces.

The last test with this simulation setup varies the material properties and hence behaviour of the purely elastic phase \(\alpha \). Compared to the values listed in Table 3 the values for \(K^\alpha \) and \(G^\alpha \) have once been multiplied by 0.2 and once by 5.0. The simulation results are shown in Fig. 10.

Influence of differently stiff material behaviour: Results for a, b \(\varepsilon _{11}\) and c, d \(\sigma _{11}\) compared to a sharp-interface solution at \(t=20\,\hbox {s}\). The results of a less stiff elastic phase \(\alpha \) are depicted in a and c, while stiffer material behaviour is shown in b and d. The quantities are plotted along x with \(y = 49.5\,\upmu \hbox {m}\)

Full size image

As previously, the PI approach stays very close to the sharp-interface solution and renders jumps symmetrically independently from changes in material properties. For the VT approach, the simulation results nearly match the sharp-interface solution when the material properties are not too different (see Fig. 10a, c). In such cases the values in the bulk regions are quite identical and the jump of the quantity is only marginally distorted. The more different the material properties and behaviours are, the greater are the deviations in values and shape of the jump (see Fig. 10b, d).

Two phases, curved interface

A plate with a circular inclusion exhibiting constant normal stress loads of \(100\,\hbox {MPa}\) is depicted in Fig. 11. The domain measures \(400\,\upmu \hbox {m} \times 400\,\upmu \hbox {m}\) and the centred inclusion of phase \(\beta \), embedded in phase \(\alpha \), has a radius of \(r = 50\,\upmu \hbox {m}\). For the simulation, only a quarter of the plate is considered, as pointed out by the dash-dotted symmetry lines.

Setup: Two linear visco-elastic phases, where an \(\alpha \)-phase surrounds a \(\beta \)-phase inclusion

Full size image

As previously, the diffuse interface width is \(l=12\,\upmu \hbox {m}\). A physical time of \(t=2\,\hbox {s}\) is simulated while the temporal discretisation is again \(\Delta t = 0.01\,\hbox {s}\). This time, both phases have a linear visco-elastic material behaviour with parameters as listed in Table 4.

In Fig. 12 the results are plotted along the \(\xi \) coordinate, which goes from the centre of the inclusion to the upper right corner (cf. Fig. 11).

Results for a \(\varepsilon _{\mathrm{tt}}\) and b \(\sigma _{\mathrm{tt}}\) along \(\xi \): Comparison of quantities tangential to the interface of the VT and the PI approach to the sharp-interface solution computed with Abaqus. Three different times (\(t=\left\{ 0.05, 1.00, 2.00\right\} \,\hbox {s}\)) are plotted along \(\xi \)

Full size image

In Fig. 12a the continuous quantity \(\varepsilon _{\mathrm{tt}}\) (component orthogonal to the direction of \(\xi \), hence tangential to the interface) is depicted. The solution of the PI approach follows the sharp-interface solution simulated with Abaqus closely, except for small non-smooth parts within the diffuse interface regions of course. The VT approach on the other hand differs from the sharp-interface solution in each time step, but with locally varying intensity. E.g., while at the beginning the difference in the results is mostly within the inclusion, towards the end of the simulation the deviation in the matrix region increases. The discontinuous quantity \(\sigma _{\mathrm{tt}}\), which is oriented orthogonally to the \(\xi \)-direction, shows for all cases a distinct jump, which is smoothed in case of the phase-field simulations. The overall picture is basically the same as before: The PI approach matches the sharp-interface simulation results, while the VT approach shows deviations in the inclusion and matrix regions.

Four phases, straight and curved interfaces

The last example in this section of binary interfaces is a two-dimensional plate of dimensions \(200\,\upmu \hbox {m} \times 200\,\upmu \hbox {m}\). It consists of two regions \(\alpha \) and \(\beta \), which are separated by an interface inclined by \(45^{\circ }\) equipartitioning the domain. Within both regions there is a circular inclusion embedded each. The positions and further dimensions may be taken from Fig. 13. Furthermore a variable displacement load (Fig. 13, right) is applied to the \(u_{x}\)-components on the x-boundaries. All other boundaries are orthogonally fixed.

Setup: 2D plate consisting of two triangular parts with an inclusion each. The x-boundaries are exposed to a variable load (see graph on the right), while the y-boundaries are fixed in y-direction

Full size image

The width of each diffuse interface is once more \(l=12\,\upmu \hbox {m}\). The variable load is applied over a time of \(t = 10\,\hbox {s}\), where during the first 5 s there is a linear increase from \(\bar{u}_x =0\,\upmu \hbox {m}\) to \(\bar{u}_x = 1\,\upmu \hbox {m}\). Thereafter the load is kept constant at \(\bar{u}_x =1\,\upmu \hbox {m}\). The temporal discretisation used is \(\Delta t = 0.1\,\hbox {s}\). The material properties of the three linear visco-elastic phases (\(\alpha \), \(\beta \), \(\gamma \)) and the purely linear elastic phase \(\delta \) are given in Table 5.

Firstly, in Fig. 14, the results for sharp-interface (simulated with Abaqus), VT and PI simulation approaches are plotted at times \(t = \left\{ 1.0, 5.0, 10.0\right\} \,\hbox {s}\). The plots go along \(\xi \) as depicted in Fig. 13. For all approaches the strain component \(\varepsilon _{\mathrm{nn}}\) (see Fig. 14a), which is oriented in the same direction as \(\xi \), has noticeable jumps in the interface regions. Especially at later times it is obvious that the VT approach fails to follow the sharp-interface solution when material properties and/or material behaviour differ too much (rate-independent elasticity versus rate-dependent visco-elasticity). In the respective bulk regions the behaviour visibly deviates. Furthermore, the jump in \(\varepsilon _{\mathrm{nn}}\) is shifted from the centre of the diffuse interface region. Compared to this, the PI approach manages to stay close to the sharp-interface solution at all times in the bulk regions and to have the jump centred in the interfacial regions.

Results for a \(\varepsilon _{\mathrm{nn}}\) and b \(\sigma _{\mathrm{nn}}\) along \(\xi \) (whole diagonal): comparison of nn-components of the VT and the PI approach to the sharp-interface solution. Three different times (\(t = \left\{ 1.0, 5.0, 10.0\right\} \,\hbox {s}\)) are plotted along \(\xi \)

Full size image

A similar picture may be observed in Fig. 14b, where the continuous stress component \(\sigma _{\mathrm{nn}}\) (component in \(\xi \)-direction) is depicted. Here again the VT approach differs in values throughout the domain, while the PI approach nearly resembles the sharp-interface solution. At \(t=10\,\hbox {s}\) there is also a deviation from the sharp-interface result seen in the phase-inherent solution. The values in the purely linear elastic phase \(\delta \) and the linear visco-elastic phase \(\gamma \) differ slightly from the reference solution, but several times less than those of the VT approach.

In Fig. 15 two points are picked along the previous plot, namely points close to the centre of the purely elastic phase \(\delta \) and the centre of the visco-elastic phase \(\gamma \).

Results for a, b \(\varepsilon _{\mathrm{nn}}\) and c, d \(\sigma _{\mathrm{nn}}\) along the simulation time: The data is plotted for \(x = y = \left\{ 49.5, 149.5\right\} \,\upmu \hbox {m}\). Compared to the VT approach, the PI approach stays for stresses and strains closer to the sharp-interface solution in both, a purely linear elastic (\(x=y=49.5\,\upmu \hbox {m}\)) and a linear visco-elastic (\(x=y=149.5\,\upmu \hbox {m}\)) material

Full size image

In Fig. 15b, d, and slightly less in a,c, the influence of the treatment of the mechanical quantities in the diffuse interface region becomes obvious. While the PI approach stays close or at least near the sharp-interface solution, the VT approach largely fails to predict the correct material behaviour. A further observation is that the deviations visibly grow over time when using the VT approach (compared to the PI approach). This happens during times when the displacement load is increased as well as in the period when the load is kept constant. As a conclusion, the mixing of material properties as done in the VT (or similar) approach leads not only to a mixture of material behaviour in respective interfaces (e.g. purely linear elastic mixed with linear visco-elastic), but also to a different evolution of the viscous strains, which introduces a second source of deviation to the material behaviour. In contrast, the PI approach allows for a clear separation of material behaviours (here rate-dependent and -independent behaviour) which minimises the effects of having a diffuse interface region present.

Behaviour in multiphase interface regions

The last simulation example focuses on a multiphase region. Therefore a plate of size \(100\,\upmu \hbox {m} \times 100\,\upmu \hbox {m}\) is divided in four equally sized regions measuring \(50\,\upmu \hbox {m} \times 50\,\upmu \hbox {m}\). The setup is illustrated in Fig. 16. At each x- and y-boundary an outward pointing orthogonal and variable displacement load is applied. The temporal development of the load is depicted on the right side of Fig. 16.

Setup: Plate consisting of four regions with the same sizes and meeting at a multi-point in the centre. The x- and y-boundaries are exposed to a variable load (see graph on the right). The coordinate \(\zeta \) starts at \((0.0\,\upmu \hbox {m},74.5\,\upmu \hbox {m})\) and \(\eta \) at \((74.5\,\upmu \hbox {m},0.0\,\upmu \hbox {m})\)

Full size image

The diffuse interface has a width of \(l = 8\,\upmu \hbox {m}\). The development and the magnitude of the variable load, the simulation time span and the temporal discretisation are the same as in the previous example. The phases and their material properties are listed in Table 6.

In Fig. 17 strain components \(\varepsilon _{\mathrm{tt}}\) and \(\varepsilon _{11}\) are plotted along \(\xi \) and \(\zeta \), \(\eta \), respectively (cf. Fig. 16). As previously, the tt-component is orthogonal to the direction of \(\xi \).

Results for a \(\varepsilon _{\mathrm{tt}}\) and b, c \(\varepsilon _{11}\) along \(\xi \) and \(\zeta \), \(\eta \), respectively: Comparison of the VT and the PI approach to the sharp-interface solution after \(t = 10.0\,\hbox {s}\)

Full size image

The PI approach again stays for all three cuts close to the sharp-interface solution. In Fig. 17a, where the multiphase region is cut, the PI approach slightly deviates from the SI solution already before the coordinate enters the multiphase region. This may indicate, that treating multiphase regions as superposed jumps of binary interfaces is not fully consistent. But compared to the VT approach, which also mostly - as before - deviates from the SI solution in the bulk areas, the PI solution is still quite close to the SI results. In both approaches, the singularity in the centre is smeared out, but in a different way. In Fig. 18, the stress components \(\sigma _{\mathrm{tt}}\) and \(\sigma _{11}\) are plotted along the same coordinates as the strain components.

Results for a \(\sigma _{\mathrm{tt}}\) and b, c \(\sigma _{11}\) along \(\xi \) and \(\zeta \), \(\eta \), respectively: Comparison of the VT and the PI approach to the sharp-interface solution after \(t = 10.0\,\hbox {s}\)

Full size image

For \(\sigma _{\mathrm{tt}}\), depicted in Fig. 18a, the picture is similar to \(\varepsilon _{\mathrm{tt}}\) previously. The PI approach matches the SI solution in the bulk regions, and also manages to stay close to the SI solution when already being in the multiphase region. But the closer \(\sigma _{\mathrm{tt}}\) is plotted to the singularity, the more deviation is seen, which further foster the thought, that additional measures have to be taken in multiphase regions. Yet the PI approach renders better results compared to the VT approach. This manifests in the purely binary interface regions shown in Fig. 18b, c.

In this work a model for the treatment of mechanical quantities on the basis of jump conditions is used and adapted to incorporate a linear visco-elastic material behaviour in a multiphase-field model. The linear visco-elastic material behaviour is represented by a Maxwell-Wiechert model and makes use of phase-inherent strains and stresses. In the simulative application this PI approach shows a close agreement to sharp-interface solutions in binary interface regions, while a conventional VT approach cannot reach the same accuracy and with rate-dependence leads to noticeable deviations in values and shape of the material behaviour. This is associated to the mixture of material properties and mechanical quantities in two- or multiphase regions in the case of the VT approach, which means a direct coupling between different material behaviours (here, e.g. linear elasticity and linear visco-elasticity, but also in the sense of differences in material properties or rate- or similar dependencies). Through the PI approach the model presented is capable to strictly separate material behaviours leading to results in close proximity to the desired sharp-interface solution on this length scale. Naturally, through this separation of mechanical quantities any other material behaviour may be easily incorporated in the basic model (e.g. plasticity or visco-plasticity, or other non-linearities). Furthermore, this indicates that the application of more than two different material behaviours in one simulation might show a similar good approximation of the sharp-interface behaviour.

Despite the good agreement in binary interfaces, the presented approach may show deviations in values in multiphase regions. Still, in such regions, the PI approach is closer to the sharp-interface solution than the VT approach. As briefly discussed, these occurring deviations for the phase-inherent model in multiphase regions may be remedied by using more sophisticated theories in such regions. Improvements may be expected from using gradient continua, as e.g. suggested by Svendsen et al. [33], or adjusted formulations for the force balance and the kinematic compatibility condition in the junction region (e.g. outlined by Simha and Bhattacharya [50]).

Furthermore, the phase-inherent modelling (compared to conventional VT or similar approaches) manages to map the sharp-interface behaviour reliably onto the diffuse interface, no matter of the diffuse interface width, the grid resolution or the difference in material properties or behaviour. It is again stressed, that in the sense of driving forces this is unconditionally necessary, especially as correct interface movement needs a certain number of grid points within.

No datasets were generated or analysed during this study. All details necessary to reproduce the simulations shown are given in the respective sections.

ABQ:

Abaqus

HS:

Hashin–Shtrikman

MWM:

Maxwell–Wiechert model

PFM:

Phase-field method

PI:

phase-inherent

RS:

Reuss/Sachs

SI:

Sharp interface

VT:

Voigt/Taylor

The research work has mainly been funded by the German Research Foundation (DFG). The model development and implementation have been funded by the DFG projects Gottfried-Wilhelm-Leibniz prize NE822/31-1, the cluster of excellence POLiS (EXC_2154) and the International Research and Training Group 2078 (IRTG 2078). Data processing is supported by the initiative “Future Fields” of the excellence strategy of the KIT. Some theoretical background was contributed through the programme “Energy Efficiency, Materials and Resources” of the Helmholtz Association.

Open Access funding enabled and organized by Projekt DEAL.

Authors and Affiliations

1. Institute for Applied Materials (IAM-CMS), Karlsruhe Institute for Technology (KIT), Kaiserstraße 12, 76131, Karlsruhe, Germany

   Felix K. Schwab, Andreas Reiter, Christoph Herrmann, Daniel Schneider & Britta Nestler

2. Insitute for Digital Materials Science (IDM), University of Applied Sciences Karlsruhe, Moltkestraße 30, 76133, Karlsruhe, Germany

   Andreas Reiter, Daniel Schneider & Britta Nestler

Contributions

FS adapted the linear visco-elasticity theory and the Maxwell-Wiechert model to the phase-inherent approach. FS, CH, DS and BN worked on the theoretical background of the phase-inherent simulation approach within the multiphase-field theory. FS and AR implemented the PI and VT approach in the simulation code. FS conducted the numerical simulations. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Felix K. Schwab.

Competing interests

The authors declare that they have no competing interests.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and permissions