Numerical aspects associated with the CPFE method

8. Numerical aspects associated with the CPFE method

8.1. General remarks

As far as the FE method is concerned, CPFE approaches can be regarded as a class of constitutive material models. Therefore, they can be implemented directly into an FE code when it is available in source form. In the case of commercial FE codes, CPFE constitutive laws are implemented in the form of a user subroutine, e.g. HYPELA2 in MSC.Marc [388] or UMAT/VUMAT in Abaqus [389]. Depending on whether the FE code is implicit or explicit, the purpose of a material model is one-or twofold: 1. Calculate the stress σ required to reach the ﬁnal deformation gradient (implicit and explicit); 2. Calculate the Jacoby matrix J = dσ/dE (implicit only, E is the symmetric strain tensor).

The stress calculation is usually implemented using a predictor–corrector scheme. Figure 23 visualizes the set-up of the clockwise loop of calculations to be performed. In principle, one could start predicting any of the quantities involved, follow the circle, and compare the resulting quantity with the predicted one. Subsequently the prediction would be updated using for instance a Newton–Raphson scheme. Various implementations were suggested using either the elastic deformation gradient F_e [390], the plastic deformation gardient F_p [391], the second Piola–Kirchhoﬀ stress S[109], or the shear rates γ˙^α[84] as a starting point. While they certainly all should lead to the same results, there are two numerical aspects to consider: First, the inversion of the Jacoby matrix occurring in the Newton–Raphson algorithm; second, the character of the equations to evaluate.

Figure 23: Clockwise loop of calculations during stress determination. (S second Piola{ Kirchhoff stress, ·γα shear rate, Lp plastic velocity gradient, mα slip direction, nα slip plane normal, Fp plastic deformation gradient, Fe elastic deformation gradient, C elasticity tensor, I identity matrix) — **Figure 23:** Clockwise loop of calculations during stress determination. (S second Piola{ Kirchhoff stress, ^·γ^α shear rate, L_p plastic velocity gradient, m^α slip direction, n^α slip plane normal, F_p plastic deformation gradient, F_e elastic deformation gradient, C elasticity tensor, I identity matrix)

Max-Planck-Institut für Eisenforschung GmbH

**Figure 23:** Clockwise loop of calculations during stress determination. (S second Piola{ Kirchhoff stress, ^·γ^α shear rate, L_p plastic velocity gradient, m^α slip direction, n^α slip plane normal, F_p plastic deformation gradient, F_e elastic deformation gradient, C elasticity tensor, I identity matrix)

Max-Planck-Institut für Eisenforschung GmbH

Regarding the ﬁrst point, one has to realize that the dimension of the Jacoby matrix is equal to the number of independent variables of the quantity that is used as predictor. These are nine variables for F_e, eight for F_p(due to volume conservation) and six for S(due to the symmetry of the stress tensor). However, if the ˙γ^αare chosen, there are at least 12 variables (slip systems in fcc crystals), up to 48 (slip systems in bcc crystals), or even more in the case of additional twinning. Inverting such large matrices (i.e. 48 × 48) is numerically quite demanding, which is the reason why such implementations require some eﬀort in reducing the number of (active) slip systems [392].

The second point concerns the numerical convergence behavior of the overall system. When starting an iteration from any other quantity than γ˙^α, the procedure involves calculating the slip rates from the stress. This is usually done using a power or exponential law. The slope of these functions is rapidly increasing, i.e. small variations in stress lead to increasingly larger deviations in the strain rate. Therefore, for large deformations, where convergence becomes a main issue, the iteration behavior of the stress loop becomes worse. However, when starting from ˙γ^α, the inverse tendency applies, i.e. stress variations with varying shear rates get smaller and smaller. This is why the second approach promises a better numerical stability at large strains, however, at the cost of dealing with a large Jacoby matrix.

8.2. Explicit versus implicit integration methods

When discussing explicit versus implicit integration schemes one has to distinguish between two aspects. Firstly, the FE solver can follow an explicit or implicit approach, and secondly the material model can be iterated using explicit or implicit integration schemes. Concerning the ﬁrst point, Harewood and McHugh [182] recently compared the efficiency of both methods when applying crystal plasticity models to forming problems. As could be expected, to some extent the outcome of this comparison is problemdependent. As a rule the explicit scheme generally seemed favorable when contact is involved.

Regarding the material model itself, e.g. the material subroutine in the case of commercial FE solvers, anything from explicit to fully implicit integration methods are possible. The task of the material model is twofold. Firstly, the stress necessary to achieve the prescribed deformation has to be determined. Secondly, the material state has to be updated. In most codes ﬁrst the stress is determined implicitly for a ﬁxed state of the material and in a subsequent step the material state is updated. In the case of a fully implicit implementation the stress then has to be determined again until convergence is achieved, while in a semi-implicit code the calculation is stopped after the state was updated. An advantage of the fully implicit scheme is, that it truly converges to the correct solution (if it converges at all), whereas the explicit solution converges generally but not necessarily to the correct solution. Since explicit schemes typically use very small time steps, semi-implicit integration schemes should generally work satisfactorily with respect to precision, while they are at the same time faster than a fully implicit scheme.

8.3. Element types

CPFE constitutive models as introduced in section 4 are formulated in a tensorial way to account for material anisotropy. Therefore, they are based on a three-dimensional stress tensor. In terms of ﬁnite element design this means, that crystal plasticity works best for 3D models and, when used for 2D models, is restricted to plane strain boundary conditions. However, it does not work for plane stress boundary conditions.

Most CPFE simulations use linear elements, i.e. elements using linear interpolation functions for the displacements. Therefore, these elements cannot describe strain gradients within one element. When the resolution of the FE mesh is reasonably ﬁne, this can be tolerated for single phase materials. However, when strong strain gradients occur, either due to boundary conditions or due to the presence of multiple phases, linear elements are usually not sufficient to correctly capture these strong gradients. In such cases higher order elements should be used.

In cases with advanced CPFE material models, such as introduced in section 4.3.2, that include strain gradients the situation becomes more complicated. The standard element formulations are only continuous in the displacements (C₀-continuous). This implies that strains can be calculated as displacement gradients, but strain gradients might be undeﬁned. To overcome this problem one has to use enhanced element formulations as in Evers et al. [26], Arsenlis et al. [9]. However, the deﬁnition of boundary conditions becomes rather complicated in the case of complex loadings for such element formulations. Therefore, many authors still use standard elements for such simulations and derive the necessary gradients from multi-element patches as described e.g. in Han et al. [393].