Fully Implicit RUNGE‐KUTTA Schemes for Ideal Elastoplasticity

The present paper derives a multifield formulation of ideal elastoplasticity to enable the consistent application of higher order accurate RUNGE‐KUTTA schemes. With the help of an axisymmetric model problem the appearing time discretization errors of distinct schemes are determined and compared.


The Principle of JOURDAIN for Dynamic Ideal Elastoplasticity
The principle of JOURDAIN for dissipative continua is based on the stationarity of the balance of energy statṗ u,εp supp σ infp γ≥0 P (u,ü,ε p , ε p , σ, γ) with P =K +Ė + P * + D. (1) Equation (1) consists of the rate of kinetic energy K, the rate of internal energy E, the power of external loads P * and the pseudopotential D characterizing plastic effects, cf. [2,4]. It furthermore depends on the stress tensor σ, the plastic multiplicator γ, the displacement vector u, the plastic strain tensor ε p as well as on corresponding first and second time derivatives. Evaluating Equation (1), performing a spatial discretization, and applying a semi-smooth NEWTON method yields the linearized semidiscrete form . Therein the vector of unknowns w = [u, ε p , σ, γ], the generalized mass matrix M, the generalized damping matrix D, the generalized stiffness matrices K 1 , K 2 , the generalized residual load vectors R 1 , R 2 , the external load vector R * 1 and the inactive as well as active set I, A are included. Moreover, r i,k+1 2 represents the nodal evaluation i of the VON MISES yield function at iteration k + 1 and d ≥ 0 an arbitrary constant, cf. [3,4].

RUNGE-KUTTA Schemes in Multifield Plasticity
For a temporal solution of Equation (2) higher order stiffly accurate fully implicit RUNGE-KUTTA integrators can be applied. Therefore, the time interval of interest [0, T ] is split into time steps ∆t = t n+1 − t n and the stages t ni = t n + c i t ni with c i ∈ [0, 1] and i = 1, .., s are introduced. At the latter points in time Equation (2) is solved simultaneously. Therefore, the link a 11 · · · a 11 a 12 · · · a 12 · · · a 1s · · · a 1s a 21 · · · a 21 a 22 · · · a 22 · · · a 2s · · · a 2s . . .
a s1 · · · a s1 a s2 · · · a s2 · · · a ss · · · a ss is established, whereby the matrix A t represents an extension of the RUNGE-KUTTA coefficients a ij , cf. [1]. Furthermore, the stage values are aggregated in the vector w t = [w n1 , ..., w ns ], the corresponding time derivatives inẇ t and the vector w t,n consists of s copies of w n . Inserting the linearization of Relation (3) into Equation (2), yields the linear system of equations where the incorporated matrices are extended properly. For a final solution a solver for linear systems of equations is applied.

Time Discretization Error Analysis for an Axisymmetric Model Problem
With the preceding procedure at hand, a time discretization error analysis for the axisymmetric benchmark problem of Fig. 1, where a time dependent sinusoidal load is applied to a cylinder with LAMÉ constants λ, µ, density ρ 0 and yield stress σ y , is performed, cf. [4]. Hence, at the end of each time interval the global h-error of a quantity X is calculated by subtracting the numerically determined values X ns from those obtained with a very small time step size X ∆t/16 ns , while identical initial values are taken into account. The local h-error estimator, however, is obtained, by performing two calculations with two distinct time step sizes simultaneously and comparing the results after each time step. Additionally, the results from the calculation using the bigger time step size are passed to the other one as initial condition for the next time step.
ns − X ns ||, q loc/glob = mean(linear fit(log(∆t), log(e loc/glob ns ))). (5) If these errors are determined for distinct time step sizes and collocated in vectors, expression (5) 3 can be used to determine the order of convergence of the time discretization scheme. A more detailed explanation is given in [4,6] and references therein. For distinct RUNGE-KUTTA schemes the orders of convergence for different field variables, obtained from the local and the mation of an Axisymmetric Steel Shaft  Moreover, Figure 6.13 (a) does not show an increase in the stress if the number of elements is extended in g Z -direction, as it is the case in Figure 6.13 (a). A variation of the number of elements in g R -direction has results in only small changes in the stress-strain diagram as well as in the plastic strain-time diagram, see Figure 6.14 (c)-(d). Hence, with the investigated biquadratic meshes a good spatial discretization seems to be found. For a further analysis the horizontal section a-a and the vertical section b-b of Figure 6.6 (a) are considered for a variety of biquadratic meshes. In Figure 6.15 and Figure 6.16 the courses of the normed axial displacement as well as the normed radial displacement are depicted in these sections at two arbitrarily chosen points in time t = 120∆t and t = 240∆t. Figure 6.15 (a) as well as Figure 6.16 (a) demonstrate the expected result that the equally applied inhomogeneous boundary condition leads to a constant axial displacement in section a-a. The radial displacement in this section follows a linear behavior at both points in time, see Figure 6.15 (c) and Figure 6.16 (c). Within the section b-b, however, nonlinearities are perceptible which can be attributed to the prevalent inertia effects, see Figure 6.15 (b)/(d) and Figure 6.16 (b)/(d). Nevertheless, the influence of the spatial discretization on the results obtained in sections a-a and b-b is negligible. Together with the results depicted in Figure 6.13 and Figure 6.14 it can be seen that a biquadratic mesh containing NE R = 2 as well as NE Z = 15 elements leads to a spatially convergent solution.
Hence, further on only this spatial configuration will be investigated. In order to emphasize the differences between the quasi-static and the dynamic analysis a series of contour plots for 16 selected points in time is displayed. These points are marked by black circles in Figure 6.17, where the applied inhomogeneous boundary condition and its time derivatives can be observed. Thus, it can be recognized that t 1 − t 3 and t 9 − t 11 coincide with a compression of the steel shaft, whereas t 4 − t 8 embody a tensional deformation.    (2) Table 1. Considering the local estimates, it becomes apparent that the Lobatto IIIC(2) method is quite far away from its theoretical order of two. In contrast, the Lobatto IIIC(3) scheme almost reaches its theoretical order of four for all considered field variables. A similar behavior can be recognized within the Radau IIA schemes. Analyzing the corresponding global measurements, however, yields a deviating picture. Apart from the Lobatto IIIC(2) scheme, all methods are estimated to reach orders well above two but not more.

Conclusion and Outlook
In the present paper the applicability of higher order accurate fully implicit RUNGE-KUTTA schemes to ideal elastoplasticity using a multifield approach is demonstrated. With the help of a model problem time discretization errors and the linked orders of convergence for distinct time integrators are compared. Orders of convergence of about five are obtained if the local h-error estimator is considered. But for the global h-error estimator an order of convergence greater than three could not be determined. An explanation for this order reduction phenomenon is still missing and hence an open task for further investigations.