## Abstract

In contrast to previous approaches that consider dissolution–precipitation creep as a multi-stage process and only simulate its governing subprocess, the present model treats this phenomenon as a single continuous process. The applied strategy uses the minimum principle of the dissipation potential according to which a Lagrangian consisting of elastic power and dissipation is minimized. Here, the elastic part has a standard form while the assumption for dissipation stipulates the driving forces to be proportional to two kinds of velocities: the material-transport velocity and the boundary-motion velocity. A Lagrange term is included to impose mass conservation. Two ways of solution are proposed. The strong form of the problem is solved analytically for a simple case. The weak form of the problem is used for a finite-element implementation and for simulating more complex cases.

## 1. Introduction

The notion ‘creep’ encompasses a large group of deformation processes. The main classification distinguishes between dislocation and diffusion creep dependent on whether the motion of defects or material-transport controls the kinetics of the deformation. A finer classification of the latter is related to the diffusion path of constituents: the motion within the bulk of the crystal is typical of Nabarro–Herring creep [1,2], whereas the motion along the crystal boundary is the main characteristic of Coble-creep [3]. The current contribution is focused on dissolution–precipitation creep (DPC), a deformation mechanism possible in the presence of a liquid in the intercrystalline space [4].

Although sometimes mistaken for Coble-creep, DPC has a significantly different underlying mechanism. DPC is typical of polycrystalline materials containing easily soluble solid components and a fluid-film in the intercrystalline space. It includes three primary phases: the dissolution of the material from the boundary surface of a crystal, the transportation of the material through the intercrystalline space and the precipitation of the material on the boundary surface of the same or of another crystal.

Devoid of a grain-shape restoring mechanism, such as nucleation, DPC is *per se* a self-exhausting deformation mechanism. Furthermore, as for any kind of creep, the geometrical rearrangement of the centres of mass of grains may require some grain boundary sliding [5]. Yet, for DPC, sliding is a potential accommodation process rather than a bulk strain producing process in contrast to its role in high-temperature superplasticity [6]. Nevertheless, grain boundary sliding may in principle affect overall deformation kinetics also in the case of diffusion creep [7,8]. Here, we implicitly assume that grain boundary sliding occurs parallel to the dissolution–precipitation mechanism, produces subordinate strain and bears no effect on overall kinetics. The latter assumption seems well justified considering that the envisioned deformation mechanism, DPC, requires transport along grain boundaries, i.e. the very same micromechanism controlling the boundaries' ‘shear viscosity’ for sliding processes, to use the terminology of Lifshitz [8].

Any of the three sequential steps involved in DPC can be rate-limiting [9–12]. A significant factor influencing creep is the architecture of the fluid-laden contact zone between crystals. For example, the estimates of boundary diffusivity vary by five orders of magnitude depending on whether the fluid appears as a continuous film on crystal boundaries or just as an impurity in the interphase space [13]. Motivated by experimental observations [14–18], some theories [9,11,12] suggest the existence of an absorbed water layer capable of supporting normal stresses and possibly exhibiting properties different from bulk water.

The geomechanical approaches treat the process mentioned either as quasi-static or by using principles of non-equilibrium thermodynamics [19,20]. However, two aspects are commonly found independent of the chosen framework: first, all models are linear in driving force. Second, strain rate is either inversely proportional to grain size *d* when the material transport is the governing step, or proportional to 1/*d*^{3} when the process is limited by boundary diffusion. The identification of the rate-controlling step from the experimental point of view is a rather problematic task. Many tests have been performed for this purpose using mono-mineralic aggregates of NaCl and KCl [10,21,22]. Further investigations clearly show that the deformation is faster if dissimilar materials participate in the process, than if a single material is considered [23].

In this paper, we present a novel model for DPC with a variational structure where no governing process is chosen *a priori*. Models with a variational structure have already been established as a standard within the mechanical community [24–27]. One of the reasons for the attractiveness of this strategy is that it eases the formulation of the coupling of different processes, another one is its suitability for a combination with the finite-element method [28,29]. The paper furthermore aims at resolving a problem associated with the formulations found in the literature. Previous models assume continuity of the stress component normal to the interface at triple and quadruple points to obtain a mathematically well-posed problem [30–34]. This assumption, however, results in a physically unrealistic hydrostatic stress state at those points in an aggregate subjected to a macroscopically deviatoric stress. Contrary to such approaches, we propose an alternative definition for driving forces using the Eshelby tensor as basis. The well-posedness is then achieved by introducing the notion of the smooth approximation of the normal component of this stress tensor.

The paper is structured as follows. Section 2 explains the mechanical model based on the minimization of a Lagrangian consisting of the elastic power and a dissipation term. While the elastic part of the model has the standard form, the dissipation term is an original contribution derived for the particular process. The process of minimization yields Euler's equations whose physical interpretation is illustrated by a simple example in §3. Subsequently, §4 considers the weak formulation necessary for the FE implementation, and the structure of the program code. Finally, some characteristic examples are shown in §5. The paper closes with conclusions and an outlook.

## 2. Continuum-mechanical model

To formulate an appropriate model for DPC, we first focus on the behaviour of an individual crystal *Ω*_{i} within its natural polycrystalline environment *Ω* (figure 1).

As the process yields significant deformation for geologically relevant, long time-scales, the model is set into the framework of the theory of finite deformations. It is furthermore assumed that the total deformation ** ϕ** is split into an elastic and an inelastic part which yields a multiplicative decomposition of the deformation gradient

**F**

Here, our approach differs from the formulation common in standard plasticity theory, where only a decomposition of the deformation gradient exists, but not of the deformation itself. According to the standard theory, the instantaneous unloading of an inelastically deformed body yields a so-called intermediate configuration that is stress free [35,36]. The intermediate configuration is not necessarily compatible, such that, in general, the elastic and inelastic deformation gradient at a material point, **F**^{E} and **F**^{I}, respectively, cannot be expressed as gradients of position vectors. An exception arises for cases for which the deformation causes a compatible intermediate configuration. This actually fits the situation characteristic of DPC where an individual grain changes its shape not by mechanical deformation, but due to motion of its surface. Thus ** ϕ**,

**u**,

For the simulation of the inelastic processes, the approach proposed by Carstensen *et al.* [37], Hackl *et al.* [38], Mielke [39,40] and Ortiz & Stainier [41] is followed. According to this concept, known as the minimum principle of the dissipation potential [42,43], the actual state of fluxes minimizes the Lagrangian consisting of elastic power *D*
*E* expressed in terms of the Helmholtz energy *Ψ*(**F**^{E}) as

The dissipation potential is assumed in a form strictly adapted to the DPC process. The reasoning is based on the physical description of the process according to which the material transport causes the motion of the boundary of the crystals. Thus, we relate the dissipation term to the material-transfer velocity **Q**_{i} and the velocity of the boundary motion **Q**_{i} constitute internal variables in the commonly used sense that they correspond to the internal structure, hidden to an external observer. While they can be in principle measured—albeit not continuously but only at certain points in time—they are not directly controllable.

The assumed specific form of the dissipation potential stipulates the velocities to be proportional to the respective driving forces
*γ* can be understood as a mobility controlled by a diffusion coefficient while *κ* is related to parameters quantifying reaction kinetics.

The mathematical background and the choice of two dissipation mechanisms in (2.6) bear important physical consequences. The majority of previous approaches assumed transport in the intercrystalline space to be the rate-limiting process, i.e. setting *κ*=0. First of all, mathematically this assumption leads to an ill-posed variational problem, i.e. the functional in (2.2) has no minimum, as arbitrarily small values can be achieved by assuming arbitrarily small boundary layers. As a physical consequence, the driving force is the normal component of the Eshelby tensor (approximately the normal stress) on the interfaces, as actually assumed in most theories. However, this driving force had to be continuous at triple or quadruple points, strictly requiring hydrostatic stress states for all of these points in a polycrystal, a requirement that appears physically unrealistic. Thus, the case *κ*=0 leads to physically as well as mathematically unacceptable consequences. One of the main goals of this work is to remedy this situation.

With the expressions (2.2)–(2.6), the complete Lagrangian takes the form
*α*_{i}, the Lagrange multiplier whose physical meaning is discussed in the following section.

## 3. Strong form

### (a) Derivation of Euler's equations

To derive the strong form of the problem, first the total Lagrangian *L*^{tot} is defined as a sum of the Lagrangian *L* and the power of the external forces *L*^{ext}
**B** represents the volume forces and **T** the tractions. The latter two quantities are related to the undeformed configuration. The same holds for the gradient operator ∇=∂/∂**X**, where **X** represents the position vector in the undeformed configuration. The first equation is valid for each point within the representative volume element *Ω* consisting of individual crystals *Ω*_{i}. The second equation is related to the boundary ∂*Ω* of the entire sample *Ω*. Furthermore, the minimization of equation (3.1) with respect to **Q**_{i} and *α*_{i} yields three equations
**b**_{i} represents the Eshelby stress tensor also known as the elastic energy-momentum tensor [36]
**F**≈**I**, *ψ*≪∥**P**∥ and **P**≈** σ**. Thus, it follows

*β*

_{i}≈

**n**

_{i}⋅

**⋅**

*σ***n**

_{i}, the normal component of the Cauchy-stress.

The reformulation of (3.3b) and (3.3c) starting with constraint (2.5) representing mass conservation leads to the final expressions for the inelastic processes on the crystal boundaries
*α*_{i} is a smoothing of *β*_{i}. For example, at triple points, the Lagrange parameter *α*_{i} is continuous but *β*_{i} is discontinuous (figure 3). Equation (3.7), on the other hand, shows that inelastic deformation occurs only in the zones exhibiting a difference between *β*_{i} and its smoothing.

### (b) Single crystal uniaxially loaded

The physical meaning of the governing equations is illustrated by considering the idealized case of a rectangular crystal uniaxially loaded. Such a situation is met if two crystals are pressed against each other. However, the symmetry of the problem allows restricting to the simulation of only one of these crystals. The example is concerned with the instant in time when the normal component of the Eshelby tensor has the following properties: it is constant over the loaded side (*β*=*c*=const.), it is equal to zero at the vertical sides and it has a jump from *c* to zero at corner points (figure 4). For such a one-dimensional case, equation (3.6) turns into an ordinary inhomogeneous differential equation
*k*_{1} and *k*_{2} are integration constants dependent on the boundary conditions. The index *i* is omitted as an individual crystal is considered. Let us specify a particular case where the following values are assumed. The normal component *β*=*c*=0.1 kN mm^{−2} on the loaded edge, *β*=0 at the endpoints and *κ*/*γ*=0.02 mm^{2}. The length of the loaded side is *l*=2 mm. The solution of the differential equation then turns into
*α* is measured in kN mm^{−2} and *x* in millimetre. The diagram of this function clearly shows that the Lagrange parameter *α* has the same value as *β* everywhere except in the vicinity of the corners (figure 4). In these zones, *α* smoothly changes from *c* to zero. The slope of the curve depends on the chosen ratio *κ*/*γ*. The obtained solution (3.8) corresponds to a small time step only as *β* is time dependent. Namely, the evolution of inelastic deformations causes the change of total deformations and thus of the Eshelby tensor and of its normal component.

According to the second governing equation (3.7), the difference between *α* and *β* is the driving force for the material motion. This supposition has already been proved by experimental observations showing that a constant stress does not initiate DPC except at the very corners of the sample [44].

## 4. Weak form and finite-element method implementation

### (a) The weak form of the problem

As illustrated by the above example, analytical solutions for the strong form of the problem can be found in simple cases of crystal geometry and load configuration. However, modelling a polycrystal affected by DPC exceeds the limitations of the analytical methods, which motivates the application of numerical procedures. The finite-element method (FEM) is used here for this purpose.

The starting point is the total potential (3.1) that, when using (2.7), (2.5) and (3.3c) takes the form
**u**, inelastic displacements *α*_{i}, and with respect to these, its variation is calculated. The variation with respect to the total displacements yields the classical principle of virtual work

Conditions (4.2)–(4.4) serve as a basis for the development of a time-integration scheme appropriate for simulation of DPC. The scheme involves a split into two stages. At the first stage, the total deformations are evaluated on the assumption that the inelastic deformations are known and constant. At the second stage, the situation is opposite: it is supposed that the total deformations are known in order to evaluate inelastic deformations.

### (b) Evaluation of total deformations

The evaluation of total deformation is based on the solution of equation (4.2) and requires the free Helmholtz energy in equation (2.3) to be specified. Here, we assume the three-field formulation as Simo *et al.* [45] adapted to a combination of elastic and inelastic deformations (appendix Ab,c)
**u**, this model depends on two new variables: hydrostatic pressure *p* and volume change *Θ*. The latter ones are introduced to decouple the volumetric and deviatoric parts of the potential, rendering the model applicable to nearly incompressible materials. With symbol **C**^{E}_{dev}, we denote the right Cauchy–Green deformation tensor corresponding to a pure elastic distortion
**S** is the second Piola Kirchhoff stress tensor and **S** and *δL*^{res} are presented in appendix Ab,c. In the particular case considered here, the total deformations are continuous everywhere, but their elastic and inelastic parts are discontinuous at the boundaries of crystals, resulting in the discontinuity of stresses and of tangent operator, an important advantage of our model. Equation (4.7) is the basis for the development of the modified *Q*_{1}–*P*_{0} element, suitable for the case when the elastic and inelastic deformations are coupled.

### (c) Evaluation of inelastic deformations

The second phase of evaluations is concerned with the determination of inelastic deformation on the basis of equations (4.3) and (4.4). Obviously, these equations are only defined on the boundaries of crystals. They represent a well-defined problem if the total deformations are known. Total deformations (equation (4.7)), by contrast, are defined over the complete volume of a polycrystal and their evaluation requires knowledge of the inelastic deformation at each point and not only on the boundaries of the crystal. Therefore, systems (4.3) and (4.4) must be extended by introducing additional conditions governing inelastic deformation inside the crystals. The non-uniqueness of the decomposition of the total deformation gradient (2.1) provides a venue for this task.

The decomposition (2.1) is non-unique because the stored energy of a system (2.3) only depends on the elastic deformations and thus different constellations of total and inelastic deformations can be associated with the same elastic state. In particular applications, the decomposition is made unique by introducing additional specifications dictated by the nature of the problem. This step can be done in different ways, for example, one can impose

The approach we employ consists of two steps. In the first step, inelastic deformations on the crystal boundaries are calculated as the solution of equations (4.3) and (4.4). Subsequently, the inelastic deformations inside the crystals are calculated by using equation (A 12) corresponding to the *Q*_{1}–*P*_{0} element in its standard form (appendix Aa). Here, the values obtained in the first step are treated as the boundary conditions. The chosen approach yields the inelastic deformations that are compatible inside the crystal and incompatible at the crystal boundaries. Moreover, the incompressibility of the inelastic deformations is guaranteed. Formally, equation (A 12) is adapted for use in the case of inelastic deformation by introducing superscript ‘I’ for all quantities
*δ*_{u}*L*^{I,ext} is omitted. The inelastic deformations are discontinuous, which requires the potential and its variations to be expressed as a sum of the contributions of single crystals.

The discrepancy between equation (4.8) written in terms of displacements, and equations (4.3) and (4.4) depending on velocities, is removed by writing equations (4.3) and (4.4) in an incremental form. Therefore, equation (4.3) is first multiplied by the time step Δ*t*
*n* is used to denote the newly updated value for the inelastic deformations while the remaining quantities have to be considered as the values known from the previous time step or from the elastic step. Beside the Lagrange parameter *α*_{i}, index *n* is introduced in order to emphasize that it has to be calculated again in each time step, but contrary to *n* rather than the time derivative **u**^{I}_{i} which is compatible with the expression (4.7).

### (d) Finite-element approximation

As the FE formulation of equation (4.7) has already been discussed in §4b and appendix A, our attention will be focused on the reformulation of the expressions (4.10) and (4.11). We consider an element *e* of the crystal *i* and particularly the boundary part of this element coinciding with the boundary part of the crystal *i* which will be denoted by
**N** represents the matrix with terms dependent on the shape functions, the superscripts *x*, **u** and *α* indicate the approximated function, and the hat symbol is introduced to denote the nodal values. In illustrative examples, we will subsequently consider two-dimensional cases with quadrilateral elements and bilinear shape functions but the derivations shown here have a general character and they are valid for other types of approximation, too. As the inelastically deformed configuration represents the reference configuration for the Lagrange parameter *α*_{i}, matrix **B** includes derivatives evaluated in the usual manner.

Using notation (4.12) and expression (4.15), the FE approximation of the first term of (4.10) turns into
*J*_{ξ} is introduced in order to relate the surface *S* in the physical coordinates with the surface *S*' in the natural coordinates. Repeating the procedure for the remaining terms of (4.10) and (4.11), the following approximations are obtained
*e*
**R**_{1}=**R**_{11}+**R**_{12} is used and the symmetry of the stiffness matrix can be easily shown as

### (e) Program code

Equations (4.2)–(4.7) and their FE-counterparts (4.25) and (A 35)–(A 39) are used as the basis for formulating the DPC element (figure 5). Here, inelastic deformations in the first time step are assumed to be zero, and the total deformations are calculated for such a state. In the further steps, crystal geometry is checked to calculate the normal vector on the inelastically deformed surface *Q*_{1}–*P*_{0} element) is applied for calculating inelastic deformations inside the crystals and subsequently the total deformations. Please note that the total deformations and the Lagrange parameters must be compatible everywhere while the inelastic deformations are compatible inside the crystals but not on their boundaries.

## 5. Numerical examples

The numerical model is applied to three kinds of microstructure. In all cases, it is assumed that the stored energy density *Ψ* (4.5) corresponds to a Neo-Hooke material, i.e.
*μ* denotes the shear modulus and *K* the bulk modulus.

The actual values of inelastic parameters *κ* and *γ* are difficult to estimate as they can vary over a large range. However, we can use the scale-invariance of the model at hand to make our results independent on the absolute values of these parameters. Only the ratio *κ*/*γ* will enter our calculations. If the parameters and variables *t*_{0}, then so will *t*_{0} for any scaling factor λ>0. All the following examples will be valid for parameter values *κ*=λ kNs mm^{−3}, *γ*=50λ kNs mm^{−5} and time-increments of the size Δ*t*=1/500λ s.

The first example is concerned with a rectangular crystal whose upper horizontal boundary is loaded. The vertical displacements on the lower horizontal boundary are constrained. A uniaxial load is chosen to be *q*=0.1 kN mm^{−2} for the sample of unit thickness *d*=1 mm. The material parameters *E*=23 kN mm^{−2} and *ν*=0.16 are assumed, corresponding to the characteristics of halite crystal [44]. Discretization of the crystal is realized by a mesh with 40×20 elements (figure 6*b*).

Some of the numerical results for the chosen example can be compared with the analytical solution discussed in §3a. The smoothing obeys *α*≈*q*=0.1 kN mm^{−2} along the loaded boundary (figure 6*c*). On the vertical sides, *α*=0 as these sides are not loaded. In the edge zones there is a smooth transition of *α* between 0.1 kN mm^{−2} and 0. These results agree with the analytical ones (see §3a and figure 6*a*). Inelastic displacements **u**^{I} (figure 6*d*,*e*) are concentrated at the edges, in agreement with the physical interpretation of the second governing equation (3.7). Finally, figure 6*f*,*g* show total displacements after 50 time steps. Subscripts in **u**^{I} and *α* can be omitted as a single crystal is treated. Results of this and subsequent examples are shown on the initial undeformed configuration by using different colour shades to indicate their intensity. Results are not normalized but they quantitatively depend on the choice of material parameters.

The second example simulates the behaviour of a polycrystal consisting of nine crystals (figure 7). Each of the crystals has four nodes and is discretized by a mesh with 14×14 elements. A test similar to the one before is modelled. A uniaxial load *q*=0.1 kN mm^{−2} is applied on the upper boundary, and vertical displacements on the lower boundary are suppressed. An important observation, not apparent in the case of a single crystal, is that here total displacements **u** and Lagrange parameter *α*_{i} are continuous everywhere, but inelastic displacements *c*,*d*) thus initiating the creep of material along the boundaries manifesting as the total deformation (figure 7*e*,*f*). In figures, a gradual change of colour shades corresponds to a compatible solution, whereas a sudden change of colour shades indicates discontinuous jumps. The compatibility of total deformations at higher than the investigated strains might require a coupling of DPC with grain boundary sliding [5,8].

The last example studies the behaviour of a polycrystal consisting of crystals with an irregular shape (figure 8). The prescribed conditions are analogous to those from the second example and the achieved results confirm the features already described: the inelastic deformations are concentrated in corners of the crystal and material flows along the crystal boundaries. This third example gives insight into the behaviour of triple and quadruple points. They represent sources of inelastic deformations, owing to the prevailing discontinuity in stress *β*_{i}, in contrast to conventional models where no specific kinetics takes place at triple and quadruple points. In the present model, the grain boundaries reorient in the vicinity of triple points trying to reduce the jumps in the normal components of the stress. The possibility of modelling crystals with an arbitrary shape is an important precondition for the simulation of more realistic polycrystalline microstructure as well as for the applicability of the model for different kinds of materials with various characteristic crystal shapes.

## 6. Conclusion

The present model for DPC uses the minimum principle of the dissipation potential according to which the Lagrangian is assumed as a sum of elastic power and dissipation. For the particular problem considered here, the elastic part has the standard form and depends on the Helmholtz energy while the dissipation potential is expressed in terms of the velocity of the material transport and boundary motion. The mentioned velocities are coupled to each other through the condition of mass conservation.

The model is first used as a base for the derivation of the strong form of the problem giving a new aspect to the simulation of DPC. Namely, the derived equations show that not the Eshelby tensor itself but the difference of its normal component *β*_{i} and its smoothing *α*_{i} is the driving force of the process. The triggered material transport eliminates the difference between *α*_{i} and *β*_{i} with time. Thus the inelastic deformation is intrinsically transient. Prevailing deformation requires changing external triggers, e.g. reorientation or fluctuation of the external stresses. The physical meaning of this consequence is elucidated on the example of a single crystal and endorsed by a comparison with the experimental results from the literature.

The second solution method presented in this contribution is concerned with the weak form of the problem further used for the FE implementation and development of the program code. This numerical procedure is first applied to the simulation of a single crystal where a comparison with the analytical solution shows an excellent agreement. The other two examples are concerned with the polycrystal behaviour widely overcoming the possibilities of the analytical solution. The last example is especially important as it shows that the model can be applied in the case of arbitrary crystal geometry characteristic for natural polycrystals.

In comparison to the previous models for dissolution–precipitation creep, the current model shows a clear advantage: it supports the experimental observation that the material transport is initiated at points with discontinuous stresses. This special feature of the model is due to the new expression for the driving forces. Determination of the inelastic kinetics parameters *κ* and *γ* certainly is one of the tasks to be considered in the future. These kinetic parameters involve many different aspects. While *γ* includes effects typical of the diffusion caused by the gradient of the chemical potential, the parameter *κ* is related to dissolution from the boundary and the precipitation on it. The testing of the model over a longer period of time as well as for a more illustrative polycrystalline microstructure are possible topics for the future, too. The latter one leads to the application of homogenization theory to model true macroscopic processes.

## Authors' contributions

All authors actively participated in the conceptualization and realization of the manuscript. S.K. developed the mathematical model together with K.H. and performed the numerical simulations. J.R. was predominately responsible for interpreting the physical and geological aspects of the model.

## Competing interests

We declare we have no competing interest.

## Funding

This work was supported by the German Science Foundation within the Collaborative Research Center ‘Rheology of the earth’ (SFB 526).

## Appendix A. The mixed variational principle and the derivation of the *Q*_{1}–*P*_{0} element

**(a) Purely elastic material behaviour**

In the case of incompressible or nearly incompressible materials, the standard variational formulation shows the effect of volume locking manifested by obtaining results corresponding to the behaviour of an apparently stiffer than modelled material. To improve the results, a few mixed methods are developed, with the three-field variational description proposed by Simo *et al.* [45]. In this model, the potential is defined as
**u** denotes displacements, *p* is the hydrostatic pressure, *Θ* the volume change and **C**_{dev}=*J*^{−2/3}**C** the distortion part of the right Cauchy–Green deformation tensor.

According to the standard procedure, the actual state of the variables minimizes the potential, which yields three conditions
*p* and *Θ* are constant over an element, transformation of (A 4) and (A 5) leads to the following expressions
*δΠ*^{res} represents the variation of the internal potential evaluated for the current displacements
*e* are expressed dependent on the shape functions placed in the matrix **N**
**G** is defined so that the following condition is satisfied
**B**, the matrix **G** consists of the derivative of the shape functions but differently ordered. Expressions (A 18) and (A 20) depend on the integral over a volume of an element *V*_{e}. Except in the case of incompressible materials, the principle of decoupling of deformation into a volumetric and a deviatoric part is often used in plastic deformation modes as they are exclusively deviatoric. In the scope of the FEM, the described derivation corresponds to the so-called *Q*_{1}–*P*_{0} element.

**(b) Coupling of elastic and inelastic material behaviour**

In this section, only the differences in comparison with a purely elastic case will be highlighted. These differences go back to two suppositions: first, the modelling of a combination of the elastic and inelastic material behaviour is based on the assumption for the multiplicative decomposition of deformation quantities

**(c) Influence of the discontinuity of the inelastic and elastic parts of deformations**

Finally, it should be taken into consideration that the inelastic and elastic part of deformations are discontinuous over the boundaries of crystals while the total deformations are continuous everywhere. Bearing this in mind, equation (A 31) turns into

- Received December 27, 2014.
- Accepted July 9, 2015.

- © 2015 The Author(s)

Published by the Royal Society. All rights reserved.