Diagonally implicit Runge–Kutta (DIRK) integration applied to ﬁnite strain crystal plasticity modeling

Diagonally implicit Runge–Kutta methods (DIRK) are evaluated and compared to standard solution procedures for ﬁnite strain crystal plasticity boundary value problems. The structure of the DIRK implementation is similar to that of a conventional implicit backward Euler scheme. It is shown that only very small modiﬁcations are required in order to transform the numerical scheme from one into the other. This similarity permits efﬁcient adaption of the integration procedure to a particular problem. To enforce plastic incompressibility, different projection techniques are evaluated. Rate dependent crystal plasticity, using a single crystal is simulated under various load cases as well as a larger polycrystalline sample. It is shown that the two-stage DIRK scheme combined with a step size control and a time continuous projection technique for the update of the plastic deformation gradient is in general more accurate than the implicit backward Euler and an update based on exponential mapping.


Introduction
A key component in modeling of crystal plasticity is the evolution equations for the plastic slip, which can be formulated as either rate independent or rate dependent, cf. [1,2]. The rate independent assumption requires continuous identification of the active set of slip systems and can lead to an ill-conditioned problem. To avoid non-uniqueness related to the selection of the active slip systems, rate dependent plasticity is frequently used, cf. [3][4][5][6]. One possibility to model rate dependent crystal plasticity models is to model the slip by a power law flow rule where all slip systems are active [7]. The drawback of the rate dependent approach used in the quasi-rate independent region is that the evolution equations become stiff. Another possibility to model rate dependence is to make use of an elastic-viscoplastic model, where the elastic domain is defined by introducing a threshold into the governing equations of the slip [8]. This approach can however lead to stiff behavior and large computational times, see [2]. In order to model polycrystals, a large number of B Sally Issa sally.issa@solid.lth.se 1 Division of Solid Mechanics, Lund University, PO Box 118, 221 00 Lund, Sweden individual grains must be considered at significant computational cost. It is therefore of high relevance to find robust and efficient numerical algorithms for solving any type crystal plasticity model. Different numerical methods used to solve crystal plasticity models are reviewed in [6] where uniaxial loading of a single crystal is considered. It is found that a fourth order Runge-Kutta algorithm used to solve an explicit slip rate formulation is computationally more efficient than an implicit slip rate approach and an implicit deformation gradient approach, where the deformation gradient is taken as the primary unknown. To simulate instability phenomena such as shear band localization in a rate dependent viscoplastic single crystal model, an implicit numerical algorithm is used in [4] where it is shown that the integration of the stiff constitutive equations can be performed using a simple and reliable algorithm. In [9], a second order Runge-Kutta algorithm is utilized to solve a three dimensional rate-dependent crystal plasticity model and it is shown that the Runge-Kutta algorithm provides higher order accuracy than first-order explicit and implicit schemes. Another method, discussed in [10], makes use of a homotopy continuation method combined with a Newton-Raphson algorithm to solve a rate dependent crystal plasticity model. The numerical results show that the algorithm is both stable and accurate. The method of varia-tional constitutive updates, where the internal variables are updated by means of a minimization problem can also be used for crystal plasticity models, cf. [11,12]. It is shown in [12] that the framework of variational constitutive updates allows for analysis of a broad spectrum of engineering applications such as size effects and texture evolution in polycrystals. In [13], it is shown that using decoupling between the grains allows the response of the individual crystal to be updated in parallel which opens up for the possibility to use GPUs.
Another solution procedure is based on the idea of considering both the balance laws and the constitutive evolution equations as a system of differential algebraic equations (DAE). This approach is proposed in [14] where the DIRK scheme is used together with a multilevel Newton algorithm, cf. e.g. [15], to integrate a rate-dependent Chaboche viscoplasticity constitutive model. In [16], the DIRK scheme is used to analyze the behavior of viscoelastic elastomers and in [17] it is applied to a finite strain continuum damage model. In these studies, results are compared to the implicit backward Euler of lower order and it is shown that the embedded DIRK method increases the accuracy of the numerical solution significantly. It is also shown that the DIRK scheme preserves the conventional implicit finite element structure found when the constitutive equations are integrated using the implicit backward Euler scheme. This means that adapting a conventional multilevel finite element code to the DIRK structure involves virtually no additional implementation costs. In [18], comparisons are performed between DIRK and linear implicit Runge-Kutta and half explicit Runge-Kutta in nonlinear finite element analyses. It is concluded that the DIRK approach is the preferred approach due to its efficiency and stability.
Use of either an implicit backward Euler approach or the DIRK method in rate dependent crystal plasticity will not per se guarantee that the plastic incompressibility requirement is fulfilled. This is due to the additive structure of the integration steps. To resolve this, an exponential map is frequently used in the update of the plastic deformation gradient, cf. [19]. As an alternative to the incompressibility preserving exponential map, projection methods can be used. In [20], the inelastic incompressibility in an elastic-viscoplastic constitutive model is considered as a side condition and is solved by a projection method. The projection method can be applied at the end of each time-step or it can be solved continuously in time. Both approaches will be studied in the present work.
In this work, numerical crystal plasticity simulations of fcc crystals are performed. The numerical results obtained using the DIRK scheme are compared to the results obtained when considering an implicit backward Euler method and an exponential map along with the Padé approximation. Different projection methods are evaluated in order to enforce plastic incompressibillity. In addition, to control the time step and to increase the robustness of the solution method, an embedded DIRK scheme is evaluated.
The paper is organized as follows: the governing equations and the constitutive model are described in Sect. 2. The formulation of the numerical method is given in Sect. 3 and in Sect. 4, numerical examples related to both single crystals and a polycrystalline structure are shown. Finally, some concluding remarks close the paper in Sect. 5.

Governing equations
For the sake of completeness, the key components of the adopted crystal plasticity model are briefly presented below along with the establishment of the mechanical balance laws. For further details on the specific crystal plasticity model, the reader is referred to [19].

Kinematics and thermodynamics
A material particle at the position X in the reference configuration Ω 0 is mapped to the position x in the current configuration Ω at time t by the unique mapping x = ϕ(X, t). The deformation gradient is given by F = ∂ X ϕ, and the right Cauchy-Green deformation tensor is defined as C = F T F. A multiplicative split of the deformation gradient allows F to be separated into an elastic and a plastic part as where superscripts e and p are introduced to identify the elastic and the plastic components, respectively. The volume change associated with the deformation is described by J = det(F), where det(·) denotes the determinant. Since metal plasticity is considered, it is assumed that the volume change is purely elastic, i.e. plastic incompressibility is enforced, whereby it holds that J = det(F e ).
The spatial velocity gradient is defined as where a superposed dot denotes the material time derivative. Making use of (1) together with (2) results in that L can be additively decomposed as L = L e + F e L p F e−1 , where the elastic and plastic velocity gradients are defined as Assuming isothermal conditions, the dissipation inequality can be expressed as where ψ is the Helmholtz free energy and τ is the Kirchhoff stress tensor, sym(·) is the symmetric part of a tensor and (·) : (·) is the contraction of two second order tensors. The mass density in the reference configuration is denoted ρ 0 . The Helmholtz free energy is assumed to be decoupled into an elastic part and a plastic part according to where C e i = (J ) −2/3 F eT F e is the isochoric part of the elastic right Cauchy-Green deformation tensor and where g α are internal variables associated with the hardening, as further discussed in Sect. 2.2.
Using the Helmholtz free energy in (5) and the assumption that the elastic deformation is a non-dissipative process, the Kirchhoff stress tensor can be identified as which allows the dissipation inequality in (4) to be reduced to where the Mandel stress tensor Σ is given as and G α are thermodynamic forces associated with the hardening variables g α defined as The crystal plasticity model is outlined in the following subsections, assuming isothermal conditions and isotropic hardening. For the fcc structure considered here, the slip occurs on three different directions on four closed packed planes, providing the {111} 110 systems, which results in a total number of n = 12 slip systems. As the plastic deformation is accommodated by concurrent slip on the individual slip systems, the evolution of L p is given by whereγ α is the slip rate, M α is the slip direction and N α the normal to the slip plane, both defined in the intermediate configuration. It is assumed that the intermediate configuration is isoclinic, meaning that the crystal lattice in the intermediate configuration is undistorted. Based on the above, the dissipation inequality in (7) is given by (11) where the resolved shear stress, τ α , is defined as

Constitutive model
The elastic part of the Helmholtz free energy is chosen to be of Neo-Hookean type and is formulated as where K and G are the bulk and shear moduli in the limit of small strains. The trace of a tensorial quantity is denoted by tr(·). The internal variable g α associated with hardening on slip system α enters the plastic part of the Helmholtz free energy as where Q is a material parameter and h αβ is a hardening matrix chosen as h αβ = δ αβ + q(1 − δ αβ ). The parameter q governs the ratio between the hardening of the current slip system and the hardening due to interaction with neighboring slip systems. The slip resistance G α that enters the dissipation inequality (11) can be identified as To complete the model description, the evolution law for the slip rateγ α , and the evolution law for the local sliṗ g α have to be postulated. The slip rate is modeled as ratedependent by using the power law formaṫ whereγ 0 and m are material parameters which define the reference shear rate and the rate sensitivity, respectively. The material parameter G 0 describes the lattice friction. For a high value of the rate sensitivity parameter m, a rateindependent response is approached which results in that the resolved shear stress on slip system α has to exceed the slip resistance, i.e. G 0 +G α , in order for the slip system to become active. The evolution law for the hardening variable g α is given bẏ where B is a material parameter defining the saturation of the hardening. Insertion of (16) and (17) into the dissipation inequality (11) allows the dissipation inequality to be reduced to

Balance of linear momentum
The balance of linear momentum in the reference configuration Ω 0 can be stated as which is fulfilled for all admissible δu such that {δu|δu = 0 on ∂Ω u } where ∂Ω u denotes the part of the boundary where the displacement is prescribed and ∂Ω 0T is the part of the boundary where the traction vector T is prescribed. The virtual strain, δ E, is defined as δ E = 1 2 (δ F T F + F T δ F) and the second Piola-Kirchhoff stress tensor S is obtained by the following pull-back operation on the Kirchhoff stress tensor Referring to (6), it can be concluded that since τ depends implicitly on the displacement field u and the internal variable F p it also holds that S = S(∇u, F p ).

Numerical solution procedure
The governing equations comprise the balance law (19) and the constitutive model defined by (13)- (17). Introducing the auxiliary variable x α = ln(1 − Bg α ) permits the evolution law in (17) to be rewritten aṡ Making use of the notationż = Ḟ pẋ α T , the evolution equations can be conveniently written on the compact forṁ The system of governing equations can now be written as The differential algebraic system of equations (DAE) in (23) will be spatially discretized by the finite element method, i.e. the displacement field u is interpolated as u = N a where N and a contain the shape functions and the nodal displacements, respectively. From the Galerkin method it follows that δu = Nδa and exploiting the fact that δa can be chosen arbitrarily yields the residual equation The domain Ω 0 is discretized by N elm elements and the integrals over each element are evaluated using N gp Gauss quadrature points. The internal variables z are evaluated in all N gp · N elm points in the domain Ω 0 . Letting z ξ denote the internal variables at ξ = 1 . . . N gp · N elm , the spatially discretized system can be stated as The superscript ξ is dropped in the following subsections.

Diagonally implicit Runge-Kutta method
The system of DAEs defined in (25) will be discretized in time using the diagonally implicit Runge-Kutta method (DIRK). For an ordinary differential equationẏ = f ( y, t), the time interval t ∈ [t 0 , t N ] is divided into n time steps, Δt n , and at each time step the solution y(t n+1 ) is approximated using s stages as where c i and b i are quadrature points and weights at stage i, respectively. At time step n and stage i the stage values y(t n + c i Δt n ) are required and they are approximated in a similar manner to provide where a i j is another set of weight factors and where the time is denoted by T ni = t n + c i Δt n . Making use of (27), the approximation in (26) can be written as The DIRK method is based on the weight factors having the property a si = b i and a i j = 0 for i > j, cf. the Butcher tableau in Table 1. This choice of weight factors implies that the solution y n+1 coincides with the solution at the last stage, i.e. with Y ns . Hence, the approximation in (27) can be rewritten as The sought approximation Y ni in (29) can be rewritten on a more compact form as where Y S,ni , with the subscript S, denotes the start values which are known at time step n and stage i and which are obtained by evaluating Using the definition in (31), it can be noted that the only unknown variable in (30) is Y ni .
Returning to the elasto-plastic boundary value problem defined by (25) and by identifying the state variables [a z] T as y, the DIRK stage values associated with (25) are given by where A ni and Z ni represent the displacements and the internal variables at time step n and stage i, respectively. The nonlinear system in (32) is solved using a multilevel Newton-Raphson scheme in each time step n and at each stage i. In [14], different DIRK schemes are investigated and it is observed that a three-stage DIRK scheme is not able to provide third-order accuracy due to lack of smoothness of the underlying equations. The two-stage DIRK method proposed in [14] is compared to high-order DIRK methods in [21] and it is proven that the two-stage DIRK method is the most efficient of the methods in the study. Therefore the two-stage DIRK, as proposed in [14], is utilized in the present work. The Butcher tableau for the two-stage DIRK scheme is provided in Table 1. Applying the two-stage DIRK scheme on the DAE system in (32) yields the system where L p ni = n α=1 γ α ni M α ⊗ N α according to (10). The DIRK scheme provides an error control that can be used to determine the time step adaptively. Denoting the DIRK solution of order q at time step n + 1 by y n+1 and the solution of order q − 1 at the same time step byŷ n+1 , the solutionŷ n+1 can be estimated using an embedded error control [14]. The embedded DIRK scheme makes use of the same weight factors a i j and c i but different b i , denoted bŷ b i . Thus, the solution of order q − 1 at time step n + 1 is obtained aŝ Based on (28) and (36), the local error can then be estimated as Comparing the solution of order q in (28) with the solution of order q − 1 in (36) it is evident that the solutions are both obtained from the stage derivatives f (Y ni , T ni ), thus the local error is provided at virtually no extra cost using the embedded DIRK method. The new time step will be based on an error measure of the displacements and internal variables by adopting the approach in [14] which provides where N dof is the number of degrees of freedom, u l is the displacement at degree of freedom l, z ξ are the internal variables at Gauss point ξ , r is a relative tolerance and u and z are absolute tolerances. Based on the error estimate, the new time step can be estimated as where f min and f max are parameters chosen such that oscillations of the adaptive time step are prevented.

Implicit backward Euler method
Application of the implicit backward Euler (BE) method on the DAE system in (25) provides Making use of the Butcher tableau in Table 1c reveals that the BE scheme can be seen as a special case of the DIRK scheme in which only one stage is employed at each time step. Moreover a comparison between (33)-(35) and (40)-(42) reveals that the structure of the BE update and the DIRK stage update are identical.

Exponential update
A frequently used alternative numerical scheme is an exponential update cf. [22]. Following this approach, the system in (25) is solved by applying the implicit backward Euler method on the balance equation and the evolution equation for the auxiliary variable x α , cf. (21). An exponential temporal discretization is used to update the plastic deformation gradient cf. [4,19]. Applying the exponential update scheme on the system in (25) provides To calculate the exponential matrix in (43), a first order Padé approximation can be used to provide

Plastic incompressibility
For small strains, plastic incompressibility is easily preserved due to the additive structure of the small strain tensor. For finite strains, however, plastic incompressibility must be handled with care since the DIRK and the BE schemes do not, in general, preserve plastic incompressibility. Since both the DIRK and the BE updates for the plastic deformation gradient are based on (34), it can not be guaranteed that holds. As a consequence, the condition det(F p ni ) = 1 must be enforced separately. In [20], the plastic incompressibility is enforced as a side-condition by the construction The construction of the projection method is based on solving the optimization problem whereF p ni is the solution obtained from (34) and F p ni is the projected deformation gradient which fulfills the plastic incompressibility. Applying the Lagrangian multiplier method leads to the following coupled non-linear system To reduce the computational cost, the following simplification

Initial state Uniaxial tension
Simple shear x y z (c) (b) (a) Fig. 1 Illustration of a initial state, b uniaxial tension and c simple shear of a single crystal, discretized by a single 8-node brick element can be employed which requires only a single non-linear equation in λ ni to be solved In [20], the nonlinear equation in (51) is to be solved at the last stage in each time step. Hence, inserting the sidecondition (46) into (51) permits the side-condition to be written as Solving the nonlinear Eq. (52) by, for example, a Newton-Raphson method, the projected plastic deformation gradient can be obtained as The method provided by (46)-(53) will in general introduce small inconsistencies in the solution due to the fact that the plastic deformation gradient is updated after determining the internal variables g α ni associated with hardening. For the set of evolution laws (16) and (17), the inconsistency caused by using (46)-(53) results in that the projected plastic deformation gradient F p ni and the hardening variables g α ni may no longer satisfy (32), and due to the exponential character of the evolution laws these small inconsistencies can cause large errors. To overcome this problem, a different approach is used in this work. To arrive at the method adopted here, (34) is first rewritten as The system of equations in (33)-(35) together with (56) can now be solved simultaneously. It can be noted that the adopted approach is similar to the time continuous projection technique discussed in [20]. As mentioned previously, for a model in which small strains are assumed, neither the exponential update nor the projection method is needed and the plastic incompressibility is fulfilled without any extra computational costs.

Numerical examples
In order to compare the different numerical solution schemes, a single crystal is represented by a single 8-node brick element. The crystal is subjected to uniaxial tension and simple shear, respectively, as illustrated in Fig. 1. In addition, a polycrystalline plate with a central hole subjected to uniaxial tension is also considered. The material parameters are taken from [19] to mimic the response of an austenitic steel, see Table 2. In all examples a reference simulation is performed using the Exponential Update (EU) scheme with the very small time step Δt = 1 × 10 −6 s. In order to compare the methods for the same computational cost, the tolerances in (38) are adjusted such that the same number of implicit steps are taken, i.e. the average time step for the two-stage DIRK scheme is set to twice that used in the EU and BE schemes, see Table 3.
In all simulations, it is assumed that all grains in the same Gauss point are subjected to the same deformation gradient. Thus, in the single crystal simulations, the von Mises effective stress is the same in all Gauss points in the element. For the polycrystalline case, the Cauchy stress tensor in one Gauss point is determined by averaging the Cauchy stress tensor over all grains at that point by evaluating where N is the number of grains in each Gauss point. Results will be shown for the internal variables by finding the maximum hardening variable g α in all slip systems. The maximum component of F p will also be studied. To simplify the notation, the following definitions will be used subsequently The relative errors of the average von Mises effective stress as well as the maximum of the internal variables are evaluated as

Uniaxial tensile deformation of a one-element single crystal
In the first example, uniaxial tensile deformation of a single crystal is modeled using a single three dimensional 8-node element of dimensions 1 × 1 × 1 mm 3 , see Fig. 1a. The simulation is performed for an initial orientation of the crystal defined by the Bunge Euler angle set (φ 1 , Φ, φ 2 ) = (0, 0, 0). This choice of angle set mimics uniaxial loading along the crystal's [100] axis. The simulations are performed with a constant displacement rate ofu x = 0.1 mm/s, applied as shown in Fig. 1b [20], is used together with the BE and DIRK schemes and a comparison where the projection method presented in (56) is used together with the BE and DIRK schemes.
The von Mises effective stress is plotted in Fig. 2 where a maximum effective stress slightly above 800 MPa is obtained at a displacement of u x = 3 mm. In Fig. 3, the relative errors of the von Mises effective stress, g max and F p max are shown for the different methods. Both the BE and the DIRK schemes simulated with the projection method in (46)-(53) result in significant errors in the von Mises effective stress and in the relative error of F p max . The relative errors in the stresses and in F p max obtained from BE are the highest whereas the DIRK scheme performs slightly better. Considering the relative error of g max shown in Fig. 3b, the BE scheme generates the largest relative error.
Next, the proposed projection method in (56) is evaluated. For this projection method, the accuracy of all integration schemes is slightly increased, see Fig. 3. The relative error of g max obtained when using the DIRK scheme is small for both projection methods considered here, cf. Fig. 3b.
It can be concluded that the BE scheme is systematically less accurate compared to the DIRK and EU schemes and that use of a time continuous projection technique will improve   the results in terms of a lower relative error, at least when considering the effective stress and the plastic deformation gradient. Based on this conclusion, only results from using the time continuous projection technique will be considered from here on. The relative error of the quantity g max is plotted as function of the time step in Fig. 4. It can be concluded that the order of the EU scheme and the DIRK scheme is higher than the order of the BE scheme. It can also be seen that the relative error of the EU and DIRK schemes are lower than the relative error of the BE scheme.

Simple shear deformation of a single crystal
Next, simple shear deformation of a single crystal is modeled, again using a single 3D 8-node element of dimensions 1 × 1×1 mm 3 , cf. Fig. 1c. The initial orientation of the crystal is defined by the Euler angle set (φ 1 , Φ, φ 2 ) = (0, 0, 0) and the BE and EU schemes are used with the time step Δt = 0.4 s. The simulations are performed with a constant displacement As discussed in [12,23], the crystal orientation will, in this deformation mode, continuously change. Figure 5 shows that the von Mises stress is continuously changing due to the evolution of the crystal orientation and as a result the relative errors in the stresses shown in Fig. 6a are significantly higher compared to the uniaxial load case. The peaks in the errors can be explained by additional slip systems being activated as deformation progresses.
The relative error of g max and F p max are plotted versus the displacement in Fig. 6b, c, respectively. The relative errors resulting from using the DIRK and EU schemes are of the same order, and significantly lower than the error generated by the BE scheme. The error estimate e from the embedded DIRK, cf. (38), is plotted together with the time step versus the displacement in Fig. 7. The fluctuations in the error and the time step are due to activation of different slip systems.
As a conclusion from the single crystal examples, the DIRK scheme with the time continuous projection technique and the EU scheme are the two integration methods which generate the smallest relative errors compared to the other methods.

Uniaxial tension of a plate with hole
As a final example, a polycrystalline plate with a centrally located hole is considered, cf. Fig. 8. The plate has the dimensions 100 × 100 × 10 mm 3 and due to symmetry only one eighth of the structure is modeled. The finite element mesh consists of 370 brick elements and the plate is subjected to tensile deformation with a constant displacement rateu x = 5 mm/s for a total time of 1 s. Each Gauss point comprises N = 100 randomly oriented grains. The rate sensitivity parameter m, cf. (16), is chosen as m = 20, which is lower than the value in Table 2, to decrease the computational time. The time step used in the BE and EU schemes is Δt = 0.004 s.
In Fig. 9, the average effective von Mises stress, cf. (57), is plotted for element 1 marked in Fig. 8. The structure is stretched 10 mm which results in a maximum stress slightly  Fig. 10a, and it can be concluded that the relative error obtained from the DIRK method is lower than the relative error obtained from using the EU and BE schemes. Using the DIRK scheme results also in higher accuracy when considering the relative errors of g max and F p max shown in Fig. 10b, c.
In Fig. 11, the error e from the embedded DIRK method, cf. (38), is plotted together with the time step Δt. The error is fluctuating around or slightly below 1 and thus the time step changes are negligible.
Summarizing the results from the single and polycrystal examples, the DIRK method provides higher accuracy when combined with the time continuous projection technique. In addition, the embedded DIRK method offers a time step control at no additional computational costs. It can be concluded that DIRK integration method should be the method of choice when solving crystal plasticity boundary value problems.

Conclusion
The implications of using a diagonally implicit Runge-Kutta (DIRK) scheme in crystal plasticity simulations have been scrutinized. The DIRK approach is compared to a classical implicit backward Euler (BE) scheme and an exponential update procedure. Different projection techniques have been evaluated to enforce plastic incompressibility. For the DIRK scheme, an embedded error estimate is used to control the step size. Both single crystal and polycrystal simulation examples are considered The relative error in terms of stresses, hardening variables and plastic deformation are investigated. It is shown that the choice of projection method to ensure plastic incompressibility is crucial for the DIRK and the implicit backward Euler schemes. A time-continuous projection approach is shown to be able to enforce plastic incompressibility without reducing the accuracy.
The single crystal simulations show that the exponential update procedure generates relative errors of the same order as the relative error from the implicit backward Euler approach, while the DIRK scheme provides higher accuracy. It is also shown that in polycrystalline simulations, the DIRK approach provides significantly higher accuracy than the methods used for comparison in the present work.
In conclusion, since the DIRK scheme results in higher accuracy with negligible additional cost, the DIRK scheme appears to be the preferable time integration scheme for crystal plasticity, when combined with a time-continuous projection approach.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.