{'fraction_non_alphanumeric': 0.03495096225018505, 'fraction_numerical': 0.010362694300518135, 'mean_word_length': 4.633698201720094, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A new conservative/dissipative time integration scheme for nonlinear mechanical systems

We present a conservative/dissipative time integration scheme for nonlinear mechanical systems. Starting from a weak form, we derive algorithmic forces and velocities that guarantee the desired conservation/dissipation properties. Our approach relies on a collection of linearly constrained quadratic programs defining high order correction terms that modify, in the minimum possible way, the classical midpoint rule so as to guarantee the strict energy conservation/dissipation properties. The solution of these programs provides explicit formulas for the algorithmic forces and velocities which can be easily incorporated into existing implementations. Similarities and differences between our approach and well-established methods are discussed as well. The approach, suitable for reduced-order models, finite element models, or multibody systems, is tested and its capabilities are illustrated by means of several examples.


Introduction
A key feature in the numerical approximations of conservative mechanical systems is their ability to exactly preserve the first integrals of their motion (energy, momenta, symplecticity, . . . ), replicating the properties of the continuous counterparts (see, e.g., [1,2]). This interest in structure preserving integrators is hence justified by the qualitative similarity between the dynamical behaviour of a mechanical system and the discrete dynamics generated by the time integration scheme [3]. In addition, a wealth of evidence supports the fact that this kind of time-stepping methods behaves extremely well for long-term simulations [4,5,6,7,8,9].
It is not easy to formulate numerical schemes that unconditionally preserve one or more invariants of the discrete motion. Generally speaking, this goal is accomplished by ensuring that some of the (abstract) geometric structures that appear in the continuous picture are replicated in the discrete dynamics. Since it is well-known that, in general, all invariants cannot be preserved for a fixed time step size scheme, different families of methods strive for the preservation of specific subsets of the various symmetries of the continuous system. For example, some numerical methods resemble discrete Hamiltonian systems [6], based on discrete gradient operators, and unconditionally preserve the energy and the (at most quadratic) momenta. Other methods emanate from discrete variational principles [10] and obtain the update formula from the stationarity conditions of these principles. In fact, it is possible to formulate methods that preserve energy, momenta, and the symplectic form of the system, if the time step size is added as an unknown to the method's equations [7].
In the context of nonlinear elastodynamics, the first energy and momentum conserving algorithms were developed by Simo and co-workers [4]. This pioneering work showed that for Saint Venant-Kirchhoff materials, such structure preserving methods can be easily obtained by a simple modification of the midpoint rule in which the stress, instead of being evaluated at the midpoint instant, should be taken as the average of the stresses at the endpoints of the time interval. Since the constitutive law is linear in the strain, this turns to be equivalent to compute the algorithmic stress with the average of the strains at the endpoints of the time interval. This simple idea was later applied to the conserving integration of shells [5], rods [11,12], contact mechanics [13], multibody systems [14,15], etc., and generalized to elastic materials of arbitrary type [16,17]. The key idea for such generalization is the definition of a discrete gradient operator, a consistent approximation of the gradient that guarantees the strict conservation of energy in Hamiltonian systems [18,19,16,20]. Alternatively, one might derive conserving methods by defining an average vector field [21,8]. In the context of the continuous Galerkin method, an optimization approach was employed to systematically develop high-order energy conserving schemes [22,23]. Very recently, a new mixed variational framework that takes advantage of the structure of polyconvex stored energy functions was proposed [24], and the properties of several formulas for the discrete gradient that are available in the literature were carefully analyzed in the context of multibody systems [25].
Many Hamiltonian problems are modeled with stiff differential equations for which conserving integration schemes might not be the most robust. For these problems, numerical methods with controllable numerical dissipation in the high-frequency range provide often a practical solution [26,27,28,29,30]. Based on a modification of the discrete gradient operator, Armero and Romero [9,31] developed a family of schemes for nonlinear three-dimensional elastodynamics that exhibits this kind of algorithmic dissipation, while preserving the momenta and providing a strict control of the energy, applicable to elastodynamics, as well as to rods and shells [32,33]. Following an alternative path based on the average vector field, Gebhardt and co-workers have proposed similar conserving/dissipative methods for general solid and structural problems [34,35].
This work considers the conservative/dissipative time integration of the equations of motion that typically arise during the analysis of nonlinear mechanical systems. More specifically, we present a novel approach that renders, by construction, methods with the desired conservation or dissipation properties. These methods discretize the equations of motion and add some perturbations related to the main field variables through a collection of ancillary linearly constrained quadratic programs that guarantee the conservation/dissipation properties. This kind of programs are analytically solvable and therefore, very attractive from the computational point of view. One possible interpretation of the contributions in this article is that it results in conservative/dissipative methods where the geometric arguments typically employed for their design have been replaced by optimality conditions. The perturbations proposed in the new methods are designed to correct some of the unwanted effects coming from the discretization of the governing equations. From a geometric point of view, the idea is to redesign the problem in such a way that the behavior of the system on the discrete constrained sub-manifold remains unaltered, but acts as an attractor for trajectories outside of it. Since the constrained programs can be solved in closed form, corrected formulas for the algorithmic internal forces and generalized velocities can be provided, and thus easily incorporated in existing simulation codes. The similarities and differences of the newly proposed method with respect to existing ones are pointed out and discussed critically.
The remaining of the article is organized as follows: In Section 2, we present the basic framework for nonlinear mechanical systems. In Section 3, we address in a comprehensive manner the new time discretization. In Section 4, we present several examples of increasing complexity for the verification of the method. Finally, conclusions, limitations and future work are presented in Section 5. Additionally, the Appendix introduces the precision quotient, with which the correctness of an implementation can be tested.

Mechanical framework 2.1 Statement
In this work we consider mechanical systems whose configuration is completely defined by a vector q ∈ Q, where Q ⊆ R n . Denoting by t the time, the state of the system at any instant is given by the pair (q, s) ∈ W ≡ T Q, where s =q is the velocity, and in which we have employed the notation( ·) = d(·) dt . The dynamical behavior of this system, for t ∈ [t a , t b ] is governed by the variational equation: where (δq, δs) ∈ T W are admissible variations of the generalized coordinates and velocities, p(q) ∈ T * s S and π(s) ∈ T * s S stand for the generalized-coordinate-based and generalized-velocitybased momenta, respectively, f int ∈ T * q Q is the vector of internal forces, f ext ∈ T * q Q is the vector of external loads that can be of conservative or non-conservative nature, and finally, ·, · represents a suitable pairing. Additionally, we assume the following two conditions: i ) the system possesses a positive-definite symmetric mass matrix M such that and, ii ) both the internal and the external forces derive from potential functions depending only on the configuration q, i.e., and we define the total potential energy of the system as V = V int + V ext . We would like to analyze next the implications that symmetry has on the form of the internal forces and the appropriate notions of linear and angular momentum in the abstract space Q. For that, the relation between the configuration space Q and the ambient space R 3 has to be carefully considered. We start by defining Φ : R 3 × Q → Q to be a smooth action of R 3 on the configuration space such that Φ(a, q) is the configuration of the system after all its points have been translated in space by constant vector a. The infinitesimal generator of this translation at q is the vector τ a (q) ∈ T q Q defined as with ∈ R. Let us now assume internal potential energy is invariant under translations, i.e., Then, choosing a one parameter curve of translations Φ( a, ·) in Eq. (5) and differentiating with respect to , it follows that a translation invariant potential implies that the internal forces satisfy To study the conservation of angular momentum, we must repeat the same argument but considering now a second smooth action Ψ : R 3 × Q → Q such that Ψ(θ, q) is the configuration of the system after all its points have rotated in ambient space by the application of a rotation exp[θ]. Defining, as before, the infinitesimal generator of this action to be the vector ρ θ (q) ∈ T q Q calculated as again with ∈ R. If the potential energy is now rotation invariant, i.e., Then the internal force must satisfy The precise notion of linear and angular momentum for the system defined in this section is provided by the following result: Theorem 2.1. Consider a mechanical system with configuration space Q ⊆ R n and vanishing external forces. Let Φ(a, ·), Ψ(θ, ·) be the translation and rotation actions on the configuration space with infinitesimal generators τ a and ρ θ , respectively, and define the linear momentum l ∈ R 3 and the angular momentum j ∈ R 3 as the two quantities that verify l, a = τ a (q), π , j, θ = ρ θ (q), π .
Then, the linear momentum is conserved if the potential energy is invariant with respect to translations. Similarly, if the potential energy is invariant under rotations, the angular momentum is a constant of the motion. Moreover, the total energy is preserved by the motion, due to its time invariance.

Time discretization
The interest in the current work is in algorithms to approximate the solution of Eq. (1). To define them, let us start by considering a partition of the interval [t a , t b ] into disjoint subintervals (t n , t n+1 ] with t a = t 0 < t 1 < . . . < t N = t b , and ∆t n = t n+1 − t n . Then, the integration algorithms that we consider are based on the midpoint approximation of Eq. (1) and of the form: where the configuration and rate, respectively, at time t n are approximated by q n , s n , we have defined π n = M s n , and we have used the notation (·) n+1/2 = 1 2 (·) n + 1 2 (·) n+1 . The update depends on the definition of an approximation to the internal force at the midpoint t n+1/2 that we have denoted as f int .
Eq. (12) provides an implicit or explicit update (q n , s n ) → (q n+1 , s n+1 ) that, together with the initial conditions of the configuration and velocity, suffices to generate discrete trajectories. We note that the approximation to the internal force in Eq. (12) is a function of two arguments that, by consistency, must satisfy for all configurations q ∈ Q.
We are interested, in particular, in formulating time integration schemes of the form (12) that preserve (some of) the invariants in the motion of the system (1), while controlling the value of the energy at all times. Let us first consider the update Eq. (12) with variations of the form (δq, δs) = (0, c), where c is an arbitrary but constant vector in T Q. Then, trivially, it follows that q n+1 − q n ∆t n = s n+1/2 .
Next, we would like to explore whether the proposed class of integration schemes preserves momenta for mechanical systems defined in the configuration space Q. The result, for every configuration space, is given next.
Theorem 3.1. Consider the time discretization (12) of a mechanical system with configuration space Q ⊆ R n . Let Φ and Ψ denote, as above, the actions of R 3 on Q representing, respectively, translations and rotations. The integration scheme preserves linear momentum if 0 = τ a (q n+1/2 ), f int (q n , q n+1 ) (15a) τ a (q n+1/2 ), π n+1 − π n = l n+1 − l n , a . (15b) for every a ∈ R 3 . Likewise, the integration algorithm preserves angular momentum if for every θ ∈ R 3 The verification of conditions (15)-(16) depends, first, on the structure of Q. For example, if we consider Q ≡ R 3n , the configuration space of n particles in three-dimensional Euclidean space, conditions (15b)-(16b) are easily verified. Conditions (15a)-(16a) depend not only on Q but also on the form of f int which has been, up to this point, left unspecified. For example, the canonical midpoint rule employs f int (x, y) = f int ((x + y)/2), and preserves both linear and angular momenta, but not energy. In turn, the Energy-Momentum method [4,5] provides an expression for this force that guarantees strict energy conservation in the discrete update map, without upsetting the preservation of momenta. Expanding on this idea, the Energy-Dissipative-Momentum-Conserving method [9,31] adds controllable energy dissipation to the solution, so small that does not upset the accuracy of the solution, yet large enough that can damp out some of the spurious oscillations in the high-frequency part of the solution.

Discrete derivative
As already mentioned, the direct evaluation of the internal forces at the midpoint configuration does not guarantee, in general, the preservation of energy. There exist however, consistent approximations of these forces that strictly enforce this property of conservative equations.
To introduce the form of this "conservative" approximation of the internal energy let us assume as in Section 2 that the internal forces derive from a smooth potential V : Q → R. Be aware that from now on, we remove the superindex int, since no external force is longer considered along the derivation presented next. The type of approximations we search for are referred in the literature as "discrete derivatives" [6] and are functions f : Q×Q → R that satisfy, for every x, y ∈ Q, two properties, namely: ii. Consistency: where D denotes the standard derivative operator.
When Q ⊂ R, there only exists one discrete derivative [8] and its closed form expression is given by with the well-defined limit lim In higher dimensions, there are actually an infinite number of discrete derivatives [20,8] since only the component of f along the direction of y − x needs to have a precise value in order to guarantee energy conservation, and its orthogonal complement is free to vary, as long as consistency of the approximation is preserved. This statement is formalized next: Any discrete derivative can be rewritten as with g(x, y) a vector-valued function such that and lim where P ⊥ (y−x) is the projection on the component perpendicular to y − x.
Proof. Let g(x, y) be defined as which tends to zero as y → x.

Conservative algorithmic force
We explore next a type of discrete derivative that is different to the one usually employed in nonlinear mechanics [4,16]. For that, we construct first a convex combination of the (exact) derivative at two configurations, i.e., (26) or in a more compact form In this expression, the scalar α cons has to be determined in order to guarantee directionality, f a is the averaged force and∆f is one half of the force jump between the configurations at times t n and t n+1 , i.e., Notice that this conservative approximation satisfies, by construction, the consistency condition (18). The satisfaction of directionality depends, as advanced, on the choice of the parameter α cons . To enforce it, we select α cons by means of an optimality condition [20], namely, as the scalar that minimizes Assuming f (x) = f (y), this optimization problem is a linearly constrained quadratic program that can be solved in closed form. Moreover, the only requirement for the optimization problem to be convex is that f (y)−f (x) 2 G > 0. Its solution can be interpreted as the discrete derivative that is closest to f m , the continuous force at the midpoint f m , namely, Here, G is a metric tensor. The Lagrangian of the optimization problem is where λ cons is a Lagrange multiplier that enforces directionality. To find the stationarity condition, the variation of L is calculated as: Now for the sake of brevity, let us introduce a discrete function defined as From now on, we refer to this function as a conservation function, which allows the preservation energy in the discrete setting for a given fixed time step ∆t. This conservation function is not unique and depends, in principle, on the shape of the approximated discrete form. The stationarity condition for the associated Lagrangian function can be reformulated as the linear system that is explicitly given by with The solution of this program is and Notice that α cons does not depend on the chosen metric. Then, we can claim that the adopted construction affords a unique definition. This feature represents a main innovation of the current work. However and up to this point, it is not clear to which extent the current formula approaches the commonly used formulas like the one due to Gonzalez [6] or the one due to Harten et al. [21]; a comparison of that second method with the current one is beyond the scope of this work. In contrast, the Lagrange multiplier λ cons depends on the chosen metric. Finally, the conservative part of the discrete force takes the following explicit form: In the context of nonlinear elastodynamics, the first term of the formula is equivalent to the definition proposed by Simo and Tarnow [4] that was derived in the context of Saint Venant-Kirchhoff materials. The second term can be interpreted as a correction for the most general hyperelastic case. The formula proposed by Gonzalez [6] cannot be algebraically reduced to the proposed expression, because the former is basically a correction for f m and the latter, for f a .
The conserving force given by Eq. (42) can be rewritten in the form of Eq. (21) and therefore, is a discrete derivative with where P y−x is the projection parallel to y − x, and

Dissipative algorithmic force
To account for dissipation, let us assume the existence of a dissipative part of the algorithmic internal force that is proportional to∆f , this is where α diss is a scalar whose precise definition is still open. This construction is supported by the analysis done by Romero [20], which showed that other choices may destroy the accuracy of the approximation. We will see later that this expression is very attractive since it provides an unifying treatment of both conservative and dissipative parts of the algorithmic internal force.
To find the value of α diss , we define a discrete dissipation functionD f (x, y), which must be positive semi-definite, at least second order in y − x to avoid spoiling the accuracy of the algorithm, and tend to 0 as x tends to y. Then α diss can be obtained as the scalar that minimizes This is also a linearly constrained quadratic program. The solution of this optimization problem can be interpreted as the smallest perturbation force that satisfies the dissipation relation y − x, f diss =D f (x, y). Once again, G is a given metric tensor and the associated Lagrangian is simply where λ diss is a Lagrange multiplier that enforces the dissipation constraint. To formulate the stationarity condition, the variation of the associated Lagrangian has to be computed. This procedure yields Noting that δf diss = δα diss∆ f , the stationarity condition of the Lagrangian can be written explicitly as The solution of this linearly constrained quadratic program is and As in the case of the conservative part of the algorithmic force, the parameter α diss does not depend on the chosen metric and therefore, it is unique. However and up to this point, it is not clear to which extent the current formula approaches already established formulas, especially the one due to Armero and Romero [9,31]. As in the case of the conservative part of the approximation, the multiplier λ diss depends on the chosen metric. Finally, the formula for the dissipative part of the algorithmic internal force takes the following explicit form: Notice that this formula has the same structure as the second term of the conservative part of the algorithmic force. However, instead of the conservation function, the dissipation function appears in the numerator. This fact suggests that a unifying formula containing both conservative and dissipative parts is possible, which makes the approach very attractive.

Generalization and preservation of momenta
Assuming an additive composition of the algorithmic approximation of the force, that is, we can write, using Eqs. (42) and (53), This formula is very compact and a simple inspection confirms that when the dissipation is zero, the directionality condition is exactly verified.
To accommodate the preservation of linear and angular momenta as discussed in Eqs. (15) and (16), Eq. (55) must be modified as indicated next: Let G be a Lie group with algebra g and coalgebra g * , which acts on the configuration space Q ⊆ R 3n by means of the action χ : G × Q → Q. For every ξ ∈ g, let ξ Q : Q → T Q denote the infinitesimal generator of the action. Following again Gonzalez [6], we can define G-equivariant derivatives. If V : Q → R is a G-invariant function, its G-invariant discrete derivative is a smooth map f G : Q × Q → R that satisfies the requirements of discrete derivatives and, moreover, the equivariance and orthogonality condition, namely, for all x, y ∈ Q, g ∈ G, and To construct a G-equivariant discrete derivative, consider invariant functions under the symmetry action denoted by π i , i = 1, 2, ...., q where q is the dimension of the quotient space Q/G. Let Π = (π 1 , π 2 , ..., π q ). If V : Q → R is G-invariant, a reduced functionṼ can be defined by the relation V =Ṽ • Π. If each of the invariants is at most of degree two, then a G-equivariant discrete derivative for V can be constructed as where, as before x, y ∈ Q. In particular, the formulation of a G-equivariant discrete derivative that preserves linear and angular momenta (cf. Eqs. (15) and (16)) is straightforward. The expression of the G-equivariant force in the current context is given by with z = (x + y)/2 and

Interpretation of the conservative algorithmic force
There exist infinite second order accurate approximations of the midpoint force that that lead to energy and momentum conserving discretizations [20]. A general expression for the algorithmic force that is in agreement with the definition of the discrete derivative is given by where G is the matrix representation of a suitable metric tensor and its associated conservation function isĈ As shown in [20], the choice of the metric tensor in Eq. (61) is crucial, since a wrong choice can destroy the accuracy of the solution even when the directionality and consistency properties are verified. The general expression can be reduced to the original formula proposed in [6] just by adopting the standard Euclidean metric tensor, i.e., Until now, this original formula has been regarded as the optimal one from the implementation point of view. By visual inspection of Eq. (42), it is also possible to claim that the formula derived in this work is as easy to implement as the original one given in Eq. (63). Moreover, the new formula requires only evaluations at the endpoints of the time interval and not at the midpoint.
Next, we would like to analyze the formula (42) in terms of the general expression provided by Eq. (61). First, defining z to be the average z = (x + y)/2, we make use of Taylor's theorem to compute and The averaged force can be expressed as and The force jump can be written as and Now putting everything together, we can rewrite Eq. (42) as Taking a look at this expression, it is apparent that the discrete force (42) is a second order perturbation of the midpoint approximation and that the metric employed for the definition of the conserving correction is just We conclude that the an integration scheme based on Eq. (27) would behave locally in a very similar manner to a method based on Eq. (61) with metric (71). Their global behavior can, in general, differ.

Dissipative algorithmic velocity
As proposed in [9,31], the generalized velocity can also be expressed as the linear combination of a conservative and a dissipative component. If the mass matrix is configuration-independent, the midpoint rule provides precisely the conservative part of the velocity. Following the ideas adopted for the formulation of the dissipative part of the algorithmic force, we find next the smallest perturbation of the midpoint velocity that guarantees dissipation according to a given dissipation functionD s (u, v), which must be non-negative, at least second order accurate in v − u , and tend to 0 as u tends to v. According to the midpoint rule, the conservative part of the algorithmic velocity can be expressed as where the velocity u corresponds to time instant t n , the velocity v corresponds to time instant t n+1 and w is the averaged velocity. This particular choice preserves linear and angular momenta. The preservation of energy is guaranteed when the relation is satisfied, with T being the kinetic energy. Equivalently, this expression can be obtained from δs n+1/2 , π(s cons ) = T (v) − T (u) when the variation is of the form δs n+1/2 = v − u. It is apparent that the conservative part of the algorithmic velocity adopted fulfills this condition without need for further corrections. Now we follow an idea that is slightly different to the one previously used to derive the dissipative part of the algorithmic force. Henceforth, let us assume the existence of a dissipative part of the algorithmic velocity proportional to s cons , this is where β diss is a scalar to be found that must guaranteed the dissipation of energy according to a given dissipation functionD s (u, v). This scalar can be obtained as the minimizer of where M is the mass matrix. Eq. (75) defines a quadratic program with linear constraints. Its solution can be interpreted as the smallest non-conservative velocity perturbation that satisfies for a discrete variation of the form δs n+1/2 = v − u, an energy dissipation according to the adopted rule. The associated Lagrangian function of the optimization problem is where µ diss is a Lagrange multiplier that enforces the dissipation constraint. To formulate the stationarity condition, the variation of the associated Lagrangian function has to be computed. This procedure yields Noting that δs diss = δβ diss s cons , the stationarity condition of the Lagrangian can be written explicitly as with The solution of this linearly constrained quadratic program is and Finally, the formula for the dissipative part of the algorithmic velocity takes the following explicit form: Assuming an additive composition of the algorithmic approximation of the velocity, that is, we can write This formula is identical to the formula proposed in [9,31] that was derived by employing only geometric arguments. This time, it can be clearly interpreted as an optimal approximation.

Final equations
The combination of all ingredients discussed here yields the full discrete formulation of the dynamic equilibrium for nonlinear mechanical systems. These consist of two residuals, one for the generalized velocities and another for the generalized coordinates, namely, where both residuals have to be minimized at every time step. This task is accomplished by means of a Newton-Raphson algorithm.
The generalized-velocity-based momentum term in its algorithmic form is The generalized-coordinate-based momentum term in its discrete version becomes The generalized-coordinate-based momentum rate term in its algorithmic form iṡ π(s n , s n+1 ) = M s n+1 − s n ∆t n .
Finally, the generalized internal force becomes In the case of accommodating the preservation of linear and angular momenta, we require the G-equivariant version given by and Π n = Π(q n ). The discrete conservation function is given bỹ and its G-equivariant version given bỹ The most basic discrete dissipation function at the level of the generalized internal force that can be chosen isD or its G-equivariant counterpart expressed as whereq n could correspond to an intermediate configuration, and D is constant, symmetric and positive semi-definite. The dissipation function for the velocity is of the form or its G-equivariant version given bỹ wheres n could correspond to an intermediate configuration. With this setting, unconditional stability in the nonlinear sense can be achieved. The chosen dissipation functions correspond to those proposed in [9,31] for the EDMC-1/2.

Numerical results
In this section, we present four numerical examples which were chosen to show the potentialities of the proposed approach. With these, we do not pretend to test the new approach exhaustively, but at rather provide some insight on its properties. For this purpose, we study first two examples involving two-mass systems with potential functions that can arise in the context of reduced-order models, and then two examples of nonlinear elastic shell structures employing a neo-Hookean material. Additionally, we briefly discuss the dissipation properties of the proposed scheme in the high-frequency range.

Reduced-order models
The first example is a mechanical system with two degrees of freedom whose potential function possesses polynomial complexity. The second one considers another mechanical system with also two degrees of freedom, but whose potential function shows non-polynomial complexity. Here, we adopt the most basic discrete dissipation functions at the level of the generalized internal force and at the level of the generalized velocity that are given bỹ in which D is a constant symmetric semi-positive definite matrix, and where χ f and χ s in R ≥0 are merely dissipation parameters. For both examples, four cases are considered: i ) fully conservative, i.e., χ f = χ s = 0; ii ) dissipative at the level of the generalized internal force, i.e., χ f = 0 and χ s = 0; iii ) dissipative at the level of the generalized velocities, i.e., χ f = 0 and χ s = 0; and, iv ) fully dissipative, i.e., χ f = 0 and χ s = 0. Additionally, to numerically gain some insight about the accuracy of the method, we provide for these two examples and all cases the second quotient of precision Q II (t) computed on the basis of the corresponding states ξ ∈ Q × S, namely ξ = (q, s). The definition of Q II (t) is presented in the Appendix.

Two-mass system with a polynomial potential
Here, we consider a nonlinear oscillatory mechanical system with potential function This kind of systems naturally arises in the context of reduced-order models, see for instance [36,37]. To perform our computations, we adopt a model with two degrees of freedom used as a demonstrator in [36].  Fig. 1 shows the idealized mechanical system under consideration and Fig. 2 presents a plot of the potential function, which is clearly convex within the region where the dynamics of the system takes place. Fig. 3 shows different plots for the solution of the fully conservative case. We can observe the very complex and nonlinear oscillatory behavior, which is also in excellent agreement with those results presented in [36].
On the left of Figs. 4, 5, 6 and 7, the evolution of the kinetic, potential, and total energies is shown. On the right of these figure, we show the second precision quotient also as a function of time. Fig. 4 evidently corresponds to the fully conservative case. Fig. 5 shows the dissipative case at the level of the generalized internal forces. Fig. 6 corresponds to the dissipative case at the level of the generalized velocities. Finally, Fig. 7 corresponds to the fully dissipative case. In the latter, the energy decay is larger than the two previous cases.
For all cases the second quotient of precision is almost constant and its value is approximately 4. Therefore, as expected, the numerical method is second-order accurate. According to Eq. (115), a method of a given order is unable to produce solutions with higher quotients of precision. In [38], it is stated that even if the method is correctly implemented, it is not trivial to find the right set of parameters in order to numerically obtain precision quotients of a high quality like the one presented herein.

Two-mass system with a non-polynomial potential
Next, we consider a nonlinear oscillatory mechanical system whose potential function is Figure 1: Two masses connected by linear springs between two walls. The first mass is also connected to the left wall through a nonlinear spring.        in which the conditions V II ab = V II ba , V N ab = V N ba and V D ab = V D ba are used for the sake of simplicity. This kind of systems could arise in the context of Euler-Bernoulli beams [39,40] or for structures with softening behavior [41]. To

Finite elasticity models
Here we analyze two finite elasticity models. The first one is a tumbling cylinder and the second one is a free-flying, single-layer, shell. In both cases, the spatial discretizations are based on a four-node shell element [42,34], i.e., an extensible-director-based solid-degenerate shell model, in which the shear locking and the artificial thickness strains are controlled by means of the assumed natural strain method. Also, the enhancement of the strain field in the thickness direction and the cure of the membrane locking are achieved by means of the enhanced assumed strain method. Such element allows to consider unmodified three-dimensional constitutives laws.
For the current study, we adopt the neo-Hookean hyperelastic material model, whose strain energy density is given byW with C, the right Cauchy-Green deformation tensor, J = det(C), I 1 = trace(C), and λ and µ are the first and second Lamé parameters, respectively.

Tumbling cylinder
This structure is a cylindrical shell subject to body loads with a prescribed time variation and was already investigated, for instance, in [5,43,42] and in many other works. The geometrical and material properties are the following: mean radius 7.5 m, height 3.0 m, thickness 0.02 m, first Lamé parameter 80 MPa, second Lamé parameter 80 MPa and mass density per volume unit 1.0 Kg/m 3 . The cylinder is discretized with 48 elements, in which 16 elements are located along the circumference and 3 elements along the height. The total number of nodes is 60. Moreover, no kinematic boundary conditions are enforced. For the dissipative case we set χ = 0.25. Fig.  15 shows the finite element model of the tumbling cylinder. Additionally, the line segments A, B, C and D, to which the spatial loads are applied, are indicated in magenta. Tab. 1 presents the values for the loads that are applied to the structure. The loads are then multiplied with a function that describes the variation of the applied force over the time, which is defined in Eq. (105). Then f ext in which the last expression is the applied load.
for 0 ≤ t < 0.5 5 − 10t for 0.5 ≤ t < 1 0 for t ≥ 1 (105) Fig. 16 shows a motion sequence for the conservative case, where the original configuration is located at the upper-left corner of the plot, and some deformed configurations are sequentially Figure 15: Tumbling cylinder -finite element representation.  Fig. 17 shows the time history of momenta and energy for the conservative case. It can be observed that the linear momentum, angular momentum and total energy vary during the time in which the external load is active, i.e., the first 1 s. After the external loads vanish, these three quantities are identically preserved through the time. These results confirm that the newly proposed integration scheme preserves momenta and energy. Although the total energy remains constant, the potential and kinematic energies vary in time, complementing each other in such a way that the total energy is perfectly constant. Fig. 18 shows the time history of momenta and energy for the dissipative case. Clearly, the momenta is identically preserved and energy is dissipated.

Free-flying single-layer plate
The structure considered in this last example is a rectangular flat plate, which consisting of a single material layer, subject to spatial loads with a prescribed time variation and was considered, for example, in [26,44,42] and in many other works. The geometrical and material properties are the following: length 0.3 m, width 0.06 m, thickness 0.002 m, first Lamé parameter 0.0 Pa, second Lamé parameter 103.0 GPa and mass density per volume unit 7.3 × 10 3 Kg/m 3 . The plate is then discretized with 120 elements, 30 elements being located along the largest dimension and 4 elements along the smallest dimension. The total amount of nodes is 155 and for the dissipative case we set χ = 0.5. Fig. 19 depicts the finite element discretization of this structure. The loads      Fig. 21 for the conservative case, proving that after the removal of the force, they all remain constant. Fig. 22 shows the values in time of momenta and energy for the dissipative case. Once again, momenta is perfectly preserved and energy is artificially dissipated.

On the dissipation properties
Dissipative schemes would be of little interest if the energy dissipation did not take place mostly in the high-frequency range. It is well-known that the dissipation of the high frequencies results in enhanced stability for the integration of stiff differential equations. Therefore, we could conclude that the only useful dissipative schemes are those that can annihilate the high-frequency content of the response without radically affecting the low frequency content of the response. In a nonlinear mechanical context, and to the best of our knowledge, there exists no dissipation function that only eliminates the high-frequency content and leaves untouched the low-frequency content. Dissipation always takes place along the whole frequency range. Moreover, there is no formal proof that the dissipation can be split in that sense for nonlinear problems and thus, we can only claim that some dissipation functions seem to be effective to address the high-frequency problem, fact that is mainly justified by experience. derivation of a new structure preserving schema that is enriched with the inclusion of numerical dissipation. The choice of a particular dissipation function is left to the structural analyst based on the special demands of the problem to be solved. Keeping these limitations in mind, the free-flying single-layer plate turns to be a suitable example to show the good dissipation properties of the new proposed scheme. Fig. 23, to the left, presents the amplitude spectrum based on the fast Fourier transform of the potential energy for both, the conservative and dissipative cases within the time range 0.06 − 0.1 s such that the direct influence of the initial transient is avoided. The subsequent analysis corresponds to the frequency range 100 − 2000 Hz and to the energy amplitude range 0 − 20 J. Fig. 23, to the right, presents the same information, but for the energy amplitude range 0 − 2 J. Clearly, the dissipative algorithm works very effectively beyond 600 Hz. For the kinetic energy, Fig. 24, to the left and to the right, shows almost identical dissipative properties. Up to 600 Hz, even if slightly different due to some dissipation within 200 − 210 Hz, the behavior for the conservative  and dissipative cases looks similar. Thus, we can claim that the proposed scheme seems to have very interesting dissipation properties. Figure 23: Effectiveness of the dissipative scheme at the potential energy level. Figure 24: Effectiveness of the dissipative scheme at the kinetic energy level.

Concluding remarks
We considered the conservative/dissipative time integration in the context of nonlinear mechanical systems. A systematic approach to derive algorithmic internal forces and generalized velocities that ensure the preservation or the controlled dissipation of energy was presented. As a main concrete result, we proposed a new second-order formula for the algorithmic internal forces. Moreover, this formula was investigated from a geometric point of view and also interpreted for the fully conservative case in terms of a general approach available in the literature. In contrast with conservative/dissipative methods available in the literature and based on the midpoint rule, the proposed formulas are perturbations of averaged evaluations, and thus not equivalent to existing ones.
The proposed methods are able to preserve the total energy of conservative equations, or add artificial dissipation in a controllable fashion, while preserving, in both cases, the linear and angular momenta of the system. Numerical tests verify all the previous assertions.
The proposed methods could be extended to integrate differential-algebraic equations or to include consistently dissipation functions involving derivatives with fractional orders, among others. The reformulation in the context of polyconvex large strain elasticity as well as of Lie Groups may yield interesting results. Beyond that, rigorous mathematical proofs on the robustness are still necessary.
where the numerator is computed as ξ(t, h, 1) − ξ(t, h, 2) = h Notice that this concept removes intrinsically the need for the exact solution of the considered initial value problem. Once again, it is possible to show that for sufficiently small time steps, the second precision quotient can be approximated by 2 n as well as in the case of the first precision quotient, namely For the integration scheme considered in this work (an energy-conservative/dissipative method), accuracy of second order can be guaranteed, meaning that log 2 [Q I (t)] ≈ 2 and log 2 [Q II (t)] ≈ 2. Let us note that for the calculation of precision quotients, h has to be chosen small enough, and the choice may vary from case to case. In addition, if ψ n (t) is very small, both tests may fail even if the implementation is right. For this reason it is sometime necessary to experiment with several initial conditions and time step sizes in order to achieve correct pictures. As a general rule, the quotients of accuracy show better performance when the trajectories are periodic or quasi-periodic.