ENERGY-CONSERVING ALGORITHMS FOR A COROTATIONAL FORMULATION∗

Standard nonlinear schemes for the the simulation of elastodynamic problems have several shortcomings when considering high-speed rotations. To tackle these problems, we use a corotational framework and a corresponding specific linearization to design new energy-conserving numerical schemes. In the two-dimensional case, an algorithm preserving also angular momentum is presented. The existence of a solution for the fully discrete setting of this algorithm is established. Numerical results illustrate the flexibility and efficiency of the proposed algorithms.


Introduction.
The time discretization of elastodynamic problems involving large rotations and translations is a well-known challenge in the simulation of continuum mechanics. Usually, geometrically nonlinear strain measures such as Green-Lagrange strains are applied to tackle this kind of problem. The resulting nonlinear equation requires iterative solvers which can be very costly, for example when considering high-speed rotation in combination with large time steps. Moreover, such methods show unsatisfying numerical results in these cases as the global motion is not reproduced correctly and spurious oscillations can be observed.
In order to accelerate the computations, one can use the framework of linear elasticity, i.e., linearize the evolution equation with respect to parameters associated with the displacements of the body. However, this kind of linearization does not show satisfying results for rotating bodies. In particular, rigid body rotations wrongly contribute to the elastic energy as the small deformation assumption is not fulfilled anymore.
As a remedy, the motion can be split into a rigid body motion and a purely elastic effective deformation. This idea was initially introduced in [10] to trace the global motion of elastic bodies in the context of aircraft and spacecraft dynamics. Then, the approach was put in an abstract framework by Fraejis de Veubeke [20], where different criteria were presented to guarantee the uniqueness of the decomposition. Nowadays, this idea is often used by the engineering community where it is referred to as "corotational formulation." In this paper, we use this splitting to carry out the linearization with respect to the effective displacements. To this end, we make use of the Fraejis de Veubeke's criterion and choose the corotational decomposition that minimizes the L 2 -norm of the effective displacements. The well-posedness of the associated system has been proved recently in [23,24].
The problem of handling large rotations in the Lagrangian kinematic description of continuum mechanics has spawned several publications, e.g., [1,2,3,25]. The application of the splitting of the deformation in combination with finite elements was introduced by Wempner [37] and Belytschko and Hsieh [7]. In the engineering literature, the term "corotational finite elements" was established for this method starting with [5,6] where the splitting was used for each element, resulting in additional degrees of freedom per element. Some developments of this approach were presented in [17,18,34,36]. The method was improved by defining a corotated frame per element [8,9]; see also the overview [27]. Crisfield [12,13] developed the concept of consistent linearization in the corotational formulation. Nour-Omid and Rankin [26], Rankin and Brognan [32], and Rankin and Nour-Omid [33] finally invented the element-independent corotational formulation, where projectors were used to separate the total deformation. Recently, the corotational formulation has been used for control and regulation [4,28], sensitivity analysis [29,30,31], and modern solution techniques for large-scale simulations [15,16]. For a more exhaustive overview of the development of the corotational formulation, we refer the reader to [19]. However, the design of suitable time integration methods for the corotational framework is still an open question; see topic E of section 7 in [19].
An important property of a suitable time discretization is the conservation of energy and angular momentum. For standard schemes, one can proceed using the discrete energy momentum method by Simo and Tarnow; see [35] or similar variants, e.g., [22]. This method has also been successfully adapted to the corotational framework for special finite element models; see [14,21]. However, conservation of the angular momentum is not always guaranteed and the effect of linearizations with respect to effective deformations has not been considered so far.
In this paper, we show how the linearization induced by a corotational formulation enables us to design a new energy-conserving time discretization of elastodynamic problems in two and three dimensions. In two dimensions, one of our algorithms also preserves the angular momentum. The existence of a solution for the fully discrete setting is established in this case. In the other cases, the angular momentum is asymptotically preserved and a priori estimates are obtained. These new algorithms are compared to standard algorithms used in the simulation of elastodynamics involving large rigid body motion.
Our new schemes are quite attractive from the computational point of view and show high accuracy in the case of large rigid body motions. Moreover, the approach can be easily extended to nonlinear strain measures and nonlinear material laws. Numerical simulations indicate that the resulting algorithms are almost as efficient as the linearized version.
Though most of the concepts presented in this paper hold for both two and three dimensions, we focus on the two-dimensional (2D) case for simplicity. However, the necessary adaptations to the three-dimensional (3D) case are presented together with a 3D version of one of our algorithms.
The rest of the paper is organized as follows. In section 2, we derive the 2D model and present the continuous linearized system of equations describing the evolution of the system. We state two time-discretized energy conservative algorithms in section 3. In section 4, the space discretization on a constrained space and the existence of a solution is discussed. In section 5, we extend one of our algorithms to the 3D case. Finally, in section 6, we show the efficiency and flexibility of our approach by means of various numerical tests.

Problem setting.
Throughout this paper, we use the notation Given a matrix A, we denote by A its transposed and by tr(A) its trace. The identity matrix of R 2 is denoted by I and R θ represents the rotation of θ, R θ := ( cos θ − sin θ sin θ cos θ ). The matrix Π is defined by Π := R π 2 = ( 0 −1 1 0 ) and A : B := tr(A B) is the usual matrix scalar product. We note that using this notation the derivative of the rotation matrix reads d dθ R θ = ΠR θ .

Nonlinear evolution equation.
Let Ω ⊂ R 2 be a bounded domain with a piecewise sufficiently smooth boundary Γ := ∂Ω describing an elastic body with mass density ρ(x). The coordinate system is chosen such that the center of gravity has the coordinates (0, 0), i.e., Ω ρx = 0. The body is subject to some volume forces f (x, t) and boundary tractions g(x, t). Let ϕ : Ω × [0, T ] → R 2 be the deformation mapping. In the fully nonlinear setting, the evolution of ϕ is described by the weak form of the equilibrium condition [11], where η is a function in V := [H 1 (Ω)] 2 , d := ϕ − x, and F := ∇ϕ = I + ∇d is the deformation gradient. The so-called Green-Lagrange strain is given by We Energy-and momentum-conserving time integration schemes for this setting are wellknown; see, e.g., [35]. However, when considering high-speed rotation in combination with large timesteps, these standard time integration schemes show poor accuracy: In this case, the rotation speed is not reproduced correctly which results in a wrong solution; see section 6.1.

Corotational decomposition of the motion.
As a remedy, the corotational formulation proposes to work with effective displacements in smaller spaces by removing the large rigid body component of the displacements involved in the elastic term of (2.1). The deformation is decomposed as follows [24]: In this equation, τ (t) ∈ R 2 describes the rigid body translation of the body and R θ its rigid body rotation. After this change of variable, the effective displacements are described by u. An important property of the decomposition (2.3) is that which shows that the nonlinear model correctly cancels the rigid body motion appearing in ϕ. In order to define a unique triple (τ, θ, u) and to work with small effective displacements, we consider the decomposition that minimizes the weighted L 2 -norm of u, Ω ρu 2 ; see [20]. A straightforward computation shows that if Ω ρ(x + u) · x = 0, the minimization condition for u is equivalent to (see [23,24]) We note that the first condition results from the translation and the second condition from the rotation term. Thus, we are looking for the effective displacements u in the linear space X := X trans ∩ X rot defined by Let us now consider the effects of this change of variable on (2.1). We start with the elastic term of (2.1). DefiningF := I + ∇u, we have, thanks to (2.4), The inertial term of (2.1) can also be rewritten in this new setting. To this end, we define the relative velocity by The acceleration term can then be expressed by (2.8) .. .
ϕ and reducing the test space to X, (2.1) now reads where we have set v = R θ η and We note that the term Ω ..
To compensate this reduction of the test space, additional constraints have to be imposed which are then used to determine the correct values of θ and τ . Thus, in addition to (2.9), we require that the linear and angular momentum are preserved. To obtain the corresponding conditions, we formally test (2.9) with v 1 trans = R θ (0, 1) , v 2 trans = R θ (1, 0) , and v rot = Π(x + u), respectively, to get the linear and angular momentum equations. Using u ∈ X trans , this gives the conservation of the linear momentum, and the conservation of the angular momentum, (2.12) .
where M and J are defined by

Linearization of the elastic term.
In order to accelerate the computation or when considering numerical methods to solve (2.9), some additional linearization can be introduced. For small effective displacements u, the strain relationship can be linearized. Using the following approximation of (2.2): andF ≈ I, the equilibrium equation (2.9) becomes linear with respect to u, where σ(u) := 2με(u) + λ tr ε(u)I. We will consider the nonlinear models (2.1) and (2.9) as a reference to measure the efficiency of our algorithms.
Let us now introduce the space Y defined by In a recent paper [23] (see also [24]), Grandmont, Maday, and Metier have obtained the following existence result concerning the system (2.11)-(2.13).

Energy conservative algorithms.
We now present our energy conservative algorithms which are based on two different time discretizations of (2.14). The first one guarantees the conservation of the energy, but not of the angular momentum. Next, an alternative scheme is presented that preserves both quantities. We briefly outline the strategies we have followed to design our algorithms. Denoting by d dt the time derivation operator, we find in terms of (2.7) and (2.8), The first algorithm is now obtained by discretizing each term of (3.2). In the second algorithm, we first discretize τ ) via (3.1) (see the first term of (3.15)).

Midpoint time discretization and discrete setting.
Let us denote the timestep by Δt and the time index by n with t n := nΔt. We use a midpoint time discretization which has the property where we have used the following notations for a generic time-discretized quantity : Let us also define the total energy of the linearized system at time t n by The angular momentum is given by Here s n is a time discretization of (2.7), depending on the algorithm under consideration.
Note that without linearization of the strain term, the energy corresponding to (2.9) is defined by which coincides with (3.5) when neglecting terms of second order with respect to ∇u n . Finally, we discretize B θ at time t n+1/2 by

Algorithm I.
A time discretization of (2.14) can be done in terms of (3.2). Special care has to be taken about the energy conservation property. To this end, we choose suitable approximations, for example ( Algorithm I then reads: Let τ 0 , θ 0 , u 0 , as well as the external forces (f n+1/2 ) n∈N and (g n+1/2 ) n∈N be given, and suppose that τ n , . τ n , θ n , . θ n and u n , Equations (3.8)-(3.10) correspond to the three equations of (2.14). We refer the reader to section 4 for details about how to guarantee that u n+1 ∈ X.
We note that using the discretized relative velocity and the fact that u n , . u n ∈ X, the angular momentum J n can be written as This quantity is not preserved exactly by Algorithm I because of the third term of the left-hand side of (3.9). It is, however, possible to specify the local order with respect to Δt. Indeed, we have Although the momentum is asymptotically preserved if Δt tends to zero, we observe that local oscillations occur; see Figure 6.11. This motivates the introduction of Algorithm II.

Algorithm II.
We now describe a time discretization of (2.14) that preserves both the energy and the angular momentum.
Remark. We note that the time discretization of the auxiliary variable s is done at the timestep n (see (3.11)) in Algorithm I and at the staggered timestep (n + 1/2) (see (3.16)) in Algorithm II. Moreover, given a parameter and denoting by I and II its time discretizations in Algorithms I and II, respectively, we have that

Energy conservation.
Let us state the main property of our algorithms.
Here, E n is the total energy defined by (3.5) where s is defined by (3.11) when considering Algorithm I and by (3.16) when considering Algorithm II. The time discretization of . ϕ is done by Proof. We give only the proof for the case of Algorithm II, since the one corresponding to Algorithm I can be obtained in a similar way.

Multiplying (3.18) by Δt
. θ n+1/2 and (3.17) by Δt and adding the two resulting equations, we get by means of (3.16), and thus Using (3.13), we end up with which yields the energy conservation in combination with (3.19).
In contrast to Algorithm I(see (3.12)), Algorithm II preserves the angular momentum; thanks to (3.14), we have

Space discretization.
In this section, we present a suitable space discretization for the time-discretized systems. We focus on (3.15), although the same method can be used for the discretization of (3.10). The space X motivates the use of a Galerkin method based on the spectral decomposition of the operator −div(σ(·)). Indeed, we note that X is spanned by the eigenvectors with nonzero eigenvalue and is orthogonal (in the sense of the L 2 (Ω) scalar product) on the kernel of this operator.
Unfortunately, finite elements do not fit into this setting. We observe that the conditions (2.5) are not satisfied by nodal finite element basis functions. The orthogonal projection of these basis functions onto the space X leads to a nonlocal support of the finite element basis and additional computational costs. We therefore use a Lagrange multiplier approach to guarantee (2.5) and state the problem in a saddle point framework. For this fully discretized setting, we show the existence of a solution of (3.13)-(3.15).

Existence of a solution for the fully discretized algorithm.
This section is devoted to the existence of a solution of the fully discretized system associated with Algorithm II.
Theorem 4.1. Given E max > 0, B max > 0, and ν > 0, suppose that at time t n we have E n ≤ E max , B n+1/2 ≤ B max , where · is a given norm on V h and Then there exists Δt * > 0 depending on E max , B max , and ν such that, given Δt ≤ Δt * , τ n , . τ n , θ n , u n , and s n , as well as the external forces f n+1/2 and g n+1/2 , there exists θ n+1 ∈ R satisfying (4.5) (with u n+1/2 the corresponding solution of (4.4)). The inversion of (4.4) provides By means of (4.6), we get for δ = 0, There exists Δt 0 depending on E max , B max , and ν such that We have then proved the existence of two roots of δ → q(δ, Δt) in L for all Δt ≤ Δt 0 . By decreasing Δt 0 , these two roots can be guaranteed to be unique: indeed, suppose that there exists a sequence (Δt n ) n∈N of [0, Δt 0 ] converging towards 0 and (r n ) n∈N a sequence of roots of h(·, Δt n ) such that ∀n ∈ N, r n = δ + (Δt n ), r n = δ − (Δt n ).
Remark. Note that the hypothesis (4.6) is the discretized version of the one of Theorem 2.1 and relies on the fact that the model is valid only for motions close to the rigid body motion, whereas the hypothesis on the parameters B max and E max stem from the discretization.

Extension to the 3D case.
This section is devoted to the extension of Algorithm II to the 3D case. Here the rotation is parameterized by a vector. As a consequence, Algorithm I cannot be extended to three dimensions since it strongly relies on the scalar nature of the angular parameter; see, e.g., (3.9).
The subset Ω is now a domain in R 3 . We denote by × the cross product defined by ⎛

Parameterization of 3D rotations.
The major problem met when extending Algorithm II to three dimensions is to parameterize the rotation corresponding to the previous R θ . We choose to describe it by using the Euler angles and refer the reader to Appendix A of [19] for a review on other parameterizations. We remark that the representation of a rotation with respect to Euler angles is not unique. Given θ ∈ R, let us introduce the rotation matrices R 1 θ , R 2 θ , and R 3 θ given by We also define for i = 1, 2, 3 the matrices Π i : . Consider now θ = (θ 1 , θ 2 , θ 3 ) ∈ R 3 . We can express a general 3D rotation matrix R θ as follows: For θ depending on the time variable t, the derivative of R θ reads where Θ θ (t) is the skew-symmetric matrix defined by Note that Θ θ (t) depends on both θ(t) and . θ(t).

Corotational model.
Before considering the system of evolution equations, we introduce the space X = X trans ∩ X rot , defined by Let us now look at the time derivatives of the deformation ϕ. As for the 2D case, we decompose ϕ into with τ (t) ∈ R 3 and u ∈ X.
Introducing the relative velocity s defined by we obtain

s(x, t) + Θ θ (t)s(x, t) .
The angular momentum is given by Note that, in contrast to the 2D case, we have to keep the rotation R θ in this definition, since Ω ρs×(x+u) cannot be collinear to the axis of the rotation. This fact motivates the introduction of the relative angular momentum m defined by In terms of m, the time derivative of J reads The dynamics of our system can be described by the following three equations: and Ω ρ(

3D Algorithm.
We are now able to extend Algorithm II to the 3D case. We keep the notation corresponding to the midpoint time discretization; see (3.4).

Angular momentum and energy conservation.
The energy conservation of the 3D Algorithm can be established as in the 2D case. However, regarding the angular momentum, we only get a weaker conservation property, namely that the norm of the angular momentum vector is preserved.
The proof of Theorem 5.1 can be obtained by following the one of Theorem 3.1. Note that we do not have exact angular momentum conservation since . J cannot be written as an explicit variation of a time discretization of J between two timesteps. As for Algorithm I in two dimensions, it is however possible to specify the local order with respect to Δt. Defining the time-discretized angular momentum J by where θ n , u n , and s n are computed by the 3D Algorithm, it can be shown that there exists C n+1/2 , depending only on where · is a norm in R 3 . Moreover, the norm of the angular momentum is preserved in the following sense:

Numerical examples.
In this section, we present some numerical examples which show the flexibility and efficiency of our method. We start with several 2D examples based on the linear Saint Venant-Kirchhoff material model, which are numerically studied using the energy-and momentum-conserving algorithm (Algorithm II). Then, we apply our corotational scheme to a nonlinear material law. At the end, we present a 3D example.
In our 2D examples, we compare the two corotational formulations (with and without the linearization of strains introduced in section 2.3) with the standard nonlinear formulation, where the linearization is not possible due to the rotation. In what follows, we call these three formulations "corotational linear," "corotational nonlinear," and "(standard) nonlinear," where "nonlinear" is always referred to using the nonlinearized Green-Lagrange strains (2.2). We remark that in the corotational nonlinear setting, two nonlinearities occur: the nonlinearity to determine the correct rotation angle and the nonlinearity in the stiffness term. For the nonlinear formulations, we use an energy-conserving time discretization as introduced, e.g., in [35].
To compare the results, we compute from the solution ϕ n of the standard nonlinear setting in a postprocessing step the translation τ n and the rotation angle θ n such that the minimization condition (2.5) is fulfilled, i.e., with d n := ϕ n − x, whereas the angular momentum is given by For the corotational formulations, the energies are computed by (3.5) and (3.7), respectively, and the angular momentum by (3.6).
As initial deformations, we use either u 0 = 0, i.e., ϕ 0 = x, or u 0 = u rot ( . θ 0 ) and ϕ 0 = ϕ rot ( . θ 0 ), respectively, where u rot , ϕ rot are the effective displacements and deformations of a body rotating with a given constant angular speed . θ 0 . We obtain u rot ( . θ 0 ) for the linearized setting by solving the system We note that multiplying (6.1) by Πx gives Mu rot · Πx = 0 and thus u rot ∈ X h . For the nonlinearized strains, ϕ rot ( . θ 0 ) is given by the solution of Here, S(·) denotes the discretized nonlinear elasticity operator. The nonlinear equation (4.5) is solved for our examples using a Newton scheme. We note that the derivative d dθ u n+1/2 can be computed by solving the system (4.4) with a modified right-hand side. For the corotational nonlinear setting, we do not use a full Newton scheme to resolve both nonlinearities but a nested iteration scheme leading here to an inner iteration where the tangential matrix can be computed in the same way as in the standard nonlinear setting. Restricting the number of inner iterations (also called inexact solving), we get an efficient solver that is easy to implement. 6.1. Examples using the standard nonlinear setting. We start with two examples demonstrating the effect of large timesteps when using the standard nonlinear setting to simulate high-speed rotating bodies; see Figure 6.1. First, we simulate a rigid body motion with a constant rotation speed of 20; for details see Example 1 later. Here we see that the larger the stepwidth is, the worse the rotational speed is approximated. Second (see Example 3), a rigid body motion with high-speed elastic oscillations is depicted. These high-speed oscillations cannot be resolved for large timesteps. Yet, it turns out that the standard nonlinear setting with energy-and momentum-conserving time integration scheme not only reproduces a wrong mean value of the rotational speed but also artificially increases the oscillation.
We assume no boundary and no volume forces, i.e., g = 0 and f = 0. As initial rotation speed, we use    With these initial conditions, it is easy to verify that the solution of (3.13)-(3.15) is given by θ n = nΔt τ 0 , u n = u 0 , s n = s 0 . The solution at various timesteps is given in Figure 6.3.
For the standard nonlinear setting, we set ϕ 0 = ϕ rot ( . θ 0 ) and θ 0 Πϕ 0 . As timesteps we use Δt = 0.05, 0.01, 0.001. We note that in our corotational formulations the solution is independent of the choice of Δt. Moreover, only one Newton iteration is needed if the iteration is started with θ n+1 such that θ n which is the correct solution. This holds for both corotational formulations with and without linearized strain. For the standard nonlinear setting, we observe an increasing number of Newton steps with increasing Δt; see Figure 6.4. Moreover, it turns out that for large Δt the nonlinear approach is not able to handle the rotation correctly resulting in a wrong rotational speed.   and the material parameters E = 500, ν = 0.2, and ρ = 1.0. We assume no boundary and volume forces and use the timestep Δt = 0.001. As initial conditions for the corotational settings, we use . θ 0 = 20 and assume no initial displacements, therefore u 0 = 0. The initial speed is given by s 0 = . θ 0 Πx and the initial rotation angle by θ 0 = 0. For the standard nonlinear setting, we have ϕ 0 = x, θ 0 Πx. We assume zero translation, i.e., τ = (0, 0) , . τ = (0, 0) . In Figure 6.5, the solution of the corotational linear setting at various times is shown. The calculated angular speed for the corotational and standard nonlinear approaches is depicted in Figure 6.6. We observe that the corotational linear approach fails to resolve the frequency of the solution. This is due to the linearized strain which leads to a wrong approximation of the elastic term for large deformations. However, the corotational nonlinear setting perfectly reproduces the results obtained by the standard nonlinear setting. When comparing the energies and angular momenta for the three formulations ( Figure 6.7), we see that all settings yield conservation of these quantities.
In Figure 6.8, the necessary Newton steps per timestep for the standard nonlinear and the two corotational settings are shown. For Δt = 0.01 we see that the standard nonlinear setting even performs somewhat better than the corotational formulations. In a second test, we increase the stepwidth to Δt = 0.1. Here, a comparable  performance is obtained for all three schemes.

Example 3: Rotation of a hard disc with no initial deformation.
Now, we consider a rotating hard disc given by a circle with radius r = 1.0 and the material parameters E = 1.62 · 10 5 , ν = 0.2, and ρ = 1.0. We use the same initial and boundary conditions as in the previous example and the timesteps Δt = 0.001, 0.005, 0.01, 0.02.
In Figure 6.9, the necessary Newton steps per timestep for the three formulations are shown. One can see that for small timesteps all settings perform quite well. For larger timesteps, however, the number of iterations using the standard nonlinear approach increases while the corotational formulations stay at low iteration numbers. When comparing the corotational linear with the corotational nonlinear setting, we can clearly see that in this example the linearization enhances the performance, which is expectable because the very stiff body is close to a rigid body motion.
The angular speed . θ obtained by the three settings is depicted in Figure 6.10. We see that for large timesteps the time integration scheme is no longer able to handle the oscillations correctly. Moreover, the standard nonlinear approach gives a wrong average rotation speed. This points out the fact that the standard time integration scheme with a large timestep is not able to handle the rotation correctly even when using the nonlinear strain relationship. Here, the corotational approach yields a good approximation of the angular speed for both strain-displacement relationships.
In Figure 6.11, we compare the angular momentum for the two different algorithms presented in this paper. While Algorithm II preserves the angular momentum, we see a small oscillation when using Algorithm I.

Example 4: Rotation of a beam under surface traction.
We consider a rotating beam under surface traction. For such rotational problems, boundary forces normal or perpendicular to the rotated configuration are quite common. This leads to a boundary force g = R θ g rel , where g rel is the boundary force relative to the rotated configuration and θ is the actual rotation angle that can be defined, for example, by (2.5). This setting causes extra difficulties when solving with standard time integration schemes as the right-hand side depends on the actual solution. However, in the corotational framework, this type of boundary condition can be modeled in a natural way and moreover gives a right-hand side that is independent of the actual rotation.
We use the same setting and initial conditions as in Example 1 and a timestep of Δt = 0.01. Here, the surface traction is given by g rel = (0, −20) on the right boundary and g rel = (0, 20) on the left boundary. On the top and bottom, we assume no surface traction. The rotation angle and angular speed using the corotational linear framework are depicted on the left of Figure 6  momentum are shown on the right. In Figure 6.13, the effective displacement u is plotted at various timesteps. We note that for this example the number of necessary Newton iterations is two for all timesteps; see Figure 6.14. Computing the same setting with the corotational nonlinear scheme leads to iteration numbers between 4 and 5 without improving the results. 6.6. Example 5: Nonlinear material law. In this example, we apply our corotational formulation to a nonlinear material law. We use the Mooney-Rivlin material model with E = 500, ν = 0.2, and c m = 1 11 . The computational domain and initial conditions are the same as in Examples 2 and 3, i.e., a circle with initial speed 20 and no initial displacements. Again, we use a stepwidth of Δt = 0.001. To achieve energy conservation, we use the method of Gonzalez [22]. The results for the rotation angle and angular speed comparing the standard nonlinear scheme with the corotational nonlinear schemes are given in Figure 6  Saint Venant-Kirchhoff material model used in Example 2, we see that the body is artificially hardened by the linearization in the material law. Good agreement of the angular speed is achieved for the two schemes. This can also be observed in the energy plots. However, in contrast to our new corotational scheme, the standard nonlinear scheme does not preserve the angular momentum; see Figure 6.15. The number of iterations per timestep is compared in Figure 6.16 for stepwidths Δt = 0.001 and Δt = 0.01. Here, both schemes yield similar iteration numbers with the standard nonlinear scheme performing somewhat better. 6.7. Example 6: Rotation of a 3D cube. In our last example, we apply our algorithm to a 3D problem. We model a rotating cube with constant rotation axis (1, 1, 1) . We note that for this kind of rotation, R can be represented in terms of just one rotation angle. Moreover, the angular momentum is exactly preserved in this example as the algorithm preserves the norm of the angular momentum. The domain is given by [0, 1] × [0, 1] × [0, 1] and the material parameters by E = 1.62 · 10 4 , ν = 0.2, and ρ = 1.0. We use Δt = 0.001, an initial angular speed of 20, and zero initial effective displacement. The deformed mesh at various timesteps is presented in Figure 6.17. As the cube is not rotationally symmetric with respect to the rotation