Discrete dynamics of complex bodies with substructural dissipation: variational integrators and convergence

For the linearized setting of the dynamics of complex bodies we construct variational integrators and prove their convergence by making use of BV estimates on the rate fields. We allow for peculiar substructural inertia and internal dissipation, all accounted for by a d'Alembert-Lagrange-type principle.


Introduction
A variational integrator is any recursive rule that allows one to calculate discrete trajectories from initial data and coincides with the discrete Euler-Lagrange equation of some discrete (or discretized) Lagrangian. It is a discrete scheme for evaluating (numerically) the dynamics described by some Hamilton or d'Alembert-Lagrange principle. A detailed description of the essential properties of variational integrators is in [18] (see also [17]).
In particular, when the potential appearing in the Lagrangian functional can be additively subdivided in subsystems, e.g. distinct groups of particles or finite elements in space in case of continuous problems, it is possible to select asynchronous time discretizations. Asynchronous variational integrators (AVIs) proposed in [9] and [10] then follow. They preserve the symplectic structure of the original continuum system. Thus, for non-constant time steps the energy evaluated numerically oscillates around its average value.
Although problems in evaluating individual trajectories are the same as for other traditional alghoritms, AVIs perform better in determining time-averaged quantities. Artificial numerical damping is not introduced. For a desired error value the computational cost is lesser than the one of other methods. AVIs display analogies with subciclyc methods (see [2] and [4]) and other methods in molecular dynamics (appropriate comparisons are in [10]).
AVIs have been introduced in [9] for the dynamics of elastic simple bodies. Their convergence in time has been analyzed variously: A repeated use of Gronwall inequality on the discrete Euler-Lagrange equations at a specific node is the basic ingredient of the technique introduced in [10] for potentials with uniformly bounded second derivative. Γ-convergence has been exploited for the zero-dimensional oscillator in [19], some technical hypotheses removed in [11]. For the linear elastodynamics of simple bodies, the convergence has been also proven in [8] along a path based on the direct exploitation of the action functional.
Here we start the analysis of variational integrators (in the synchronous and asynchronous versions) for the linear dynamics of complex bodies, those bodies for which the material substructure prominently influence the gross behavior. We consider cases in which the material substructures display peculiar inertia (additional to the macroscopic one) and internal dissipation. The dynamics of such bodies is governed by a d'Alembert-Lagrange-type principle.
We do no treat any specific model of complex materials. We develop our analyses in the general model-building framework of the mechanics of complex bodies ( [5], [12]). In this way we do not specify the type of the substructure under analysis. We presume only that its essential geometrical features are described by an element of an abstract finite dimensional manifold M, the so-called manifold of substructural shapes. We construct first the linearized version of the natural non-linear theory. To this aim we need the isometric emebedding of M in a linear space. Although such embedding always exists, it is not unique, so its choice is a constitutive part of the modeling procedure. After establishing the structure of the linear theory, we construct appropriate variational integrators and we prove the convergence of the discrete approximation in time to the continuous counterpart, the convergence in space being assured by standard theorems.
The technique that we use here is new: it is based essentially on BV estimates on the time-rates of the fields involved.
Convergence in the asynchronous case is proven only for substructures admitting quadratic peculiar independent kinetic energy. The expression of the substructural kinetic energy is accepted in fully generality in proving convergence of synchronous variational integrators. Our work prepare the way to a number of numerical experiments in various classes of complex bodies.
Some notations. We use the standard symbols | · | for the euclidean norm, and both · and ·, · for the scalar product in some Euclidean space R j , j ∈ N. We will not explicitly indicate the dependence on the dimension j, for the sake of the conciseness of the notation. The symbol # denotes the cardinality: #T is then the cardinality of the set T . By ∂ z we indicate the partial derivative with respect to z. We use standard notations for Lebesgue spaces, Sobolev spaces, and the space of functions with bounded variation (BV in short). The symbol H d−1 denotes the (d − 1)-dimensional Hausdorff measure in R d . In some computations we find inequalities involving constants which depend on the data and on the space dimensions, but are always independent of the set of discrete instants. Since it is not essential to distinguish from one specific constant to another, we indicate all of them by the same letter c, leaving understood that c may change from one inequality to another.

Commentary to the dynamics of complex bodies
In a primitive approach, a body is considered as an abstract set B the elements of which are called material elements. Each element is considered as the smallest piece of matter characterizing the nature of the material constituting the body under examination, that is a patch of matter made of entangled molecules or the characteristic cell of some atomic lattice. The representation of this set is a model of the morphology of the body. It cannot prescind from the selection of the placements of the body in the ambient space R d , placements that are selected by means of one-to-one maps f : B −→ R d . The generic B := f (B) represents the 'gross' shape of the body and is assumed to be a bounded domain with boundary ∂B of finite (d − 1)-dimensional measure, a boundary where the outward unit normal n is defined to within a finite number of corners and edges. When substructural complexity arises, description of the structure inside material elements is necessary. The representation of B is then defined by maps κ : B → M assigning to each material element a descriptor of the peculiar features of its inner morphology (called for this reason a morphological descriptor), with M a finite-dimensional differentiable manifold.
The dynamics of a complex body involves then the description of subsequent changes in gross and substructural shapes.
Convenience suggests to select a reference place B * := f * (B) that can be in principle occupied by the body in R d and to obtain from it new places by means of sufficiently smooth bijective maps (transplacements) such that (i) B = y (B * ) is as regular as B * and (ii) the gradient of deformation is the motion at gross scale. Sufficient smoothness in time is presumed so that the velocityẏ = dy dt (x,t) ∈ T y(x,t) B ≃ R d can be defined as a field over B * . Differentiable maps of the type define fields of morphological descriptors. Then, a substructural 'motion' is then given by Assuming sufficient smoothness in time, the rate of change of the substructural morphology in the referential description isν := ∂ν ∂t (x,t). At each x the spatial derivative of x −→ ν is indicated by N ∈ Hom (T x B * , T ν M).
In a conservative setting the mechanics of a complex body is governed by Hamilton principle of least action with Lagrangian density assumed to be differentiable with respect to all its entries. ρ is the referential mass density (conserved during the motion), e the elastic potential and w the potential of external actions, all per unit mass. The term χ (ν,ν) is a real-valued function vanishing whenν = 0, it is of degree 2 inν and is such that the second derivative ∂ 2 νν χ exists and is positive definite. Moreover, supν (∂νχ ·ν − χ) is attained (as a maximum) at a unique point in the relevant T ν M. It coincides with the substructural kinetic energy when peculiar independent inertia pertains to substructural changes (see remarks in [6], [14]). The derivatives of L with respect to F and N represent contact interactions, namely standard and microscopic stresses, respectively. The derivatives of L with respect to y and ν are standard and substructural bulk actions, the latter being divided additively in external fields acting on the substructure and internal selfactions.
Dissipative effects may occur. For the sake of simplicity it is assumed here that dissipation occurs only within each material element and is accounted for by the local self-actions. A concrete example is given by quasicrystals (see [13], [14]).
A d'Alembert-Lagrange-type type principle then substitutes the Hamiltonian one. One then requires that where δ indicates the first variation, z v = z v (F, ν,ν, N ) ∈ T * ν M represents the dissipative part of the self-actions inside the generic material element and (x, t) −→ φ := φ (x, t) ∈ T ν(x,t) M is an arbitrary smooth field with compact support.
To define variations, for ε ∈ [0, 1] one considers smooth maps The variation δ along (y, ν) has then the meaning of the derivative d dε of L evaluated at ε = 0. It is also assumed that the self-action z v is intrinsically dissipative, that is for any choice ofν. Then z v may be of the type withν ♭ the covector associated withν and η a positive diffusion coefficient.
When the fields x −→ P := −ρ∂ F L and x −→ S := ρ∂ N e are of class C 1 (B * ) ∩ C 0 B * , Euler-Lagrange equations are then P is the Piola-Kirchhoff stress tensor, b := −ρ∂ y w the co-vector of body forces, S the microstress tensor in referential configuration, z := ρ∂ ν e the self-action, β := −ρ∂ ν w external body action over the substructure.
By considering fields over the actual place B, the balance equations above become A requirement of objectivity involving the action of SO (3) on both the ambient space R d and M for the density of the elastic potential e implies that skwP F * = e Ā * z + DĀ * S , whereĀ =Ā (ν) ∈ Hom R d , T ν M , and e is Ricci's alternating symbol (see [5]). Moreover, one may require the invariance of the entire Lagrangian under the action of the group of diffeomorphisms of the ambient space and the group of diffeomorphisms of M into itself. Such invariance requirement opens a path that allows one to prove the covariance of the balance laws (see [14]).
Boundary conditions can be assigned in terms of places and morphological descriptors on parts ∂B * y and ∂B * ν of the boundary ∂B * . Tractions t := P n and ̟ := Sn, with n the normal to the boundary in all places in which it is defined, can be in principle assigned on parts ∂B * t and ∂B * ̟ of the same boundary. Mixed boundary conditions require ∂B * y ∩ ∂B * t = ∅, ∂B * ν ∩ ∂B * ̟ = ∅ and ∂B * y ∪ ∂B * t = ∂B * ν ∪ ∂B * ̟ = ∅.
A loading device able to prescribe at the boundary the substructural contact interactions ̟ is not know in imagination or elsewhere. The sole natural boundary conditions in terms of substructural contact interactions is then ̟ = 0. Different is the case of Dirichlet data in terms of morphological descriptors. Example is the case of quasicrystals in nematic order: surfactants spread along the boundary allow to assign a precise orientation to the rod-like molecules in nematic order.

Linearization
In the standard theory of elasticity, the space of transplacements can be considered as an infinite dimensional manifold. Linearization is developed over the tangent space of such a manifold. In such a procedure it is useful to define the displacement field The initial condition is then given by (x, t 0 ) → u 0 (x, t 0 ), and the boundary condition is given along ∂B * u = ∂B * y .
The condition |∇u| << 1 defines the infinitesimal deformation regime in which, essentially, σ ≈ P , and the stress may depend linearly on the infinitesimal strain ε := sym∇u.
To linearize the multifield setting described in previous section, an additional difficulty arises because M does not coincide in general with a linear space isomorphic to R k for some k. For this reason, the construction of a Sobolev space of maps taking values on M requires the embedding of the manifold of substructural shapes itself in a linear space. We assume that M is endowed with a C 1 Riemannian structure and the relevant Levi-Civita parallel transport. By Nash theorem an isometric embedding of M in a linear space is always available but it is not unique. The selection of an embedding becomes then a constitutive ingredient of each special model. Once the embedding of M has been chosen, one may consider the space C of pairs of maps (y, ν) can be considered as an infinite dimensional differentiable manifold.
We maintain the assumption of infinitesimal deformation setting, a regime in which we can 'confuse' B * with B, in the sense that if ϕ indicates a vector field tangent to y, then ϕ ≈ u at any x in B * . For this reason, from now on we write B instead of B * for the sake of conciseness. Moreover, we also confuse ν with ν e = ν • y −1 for consequent obvious reasons.
We consider differentiable fields such that ϕ and φ vanish where Dirichlet data are imposed. We continue to write ν and M as before the embedding. The notation L (A) (ū, ϕ) (ν, φ) indicates the linearization of A about the pair of maps (ū,ν). A superposed bar denotes maps calculated at (ū,ν). At each x and t, the linearizations of the maps x −→ F and ν −→ N about (ū,ν) are then given by Previous formulas make sense because of pointwise parallel transport of F over curves on R d , and of N over M. The next step is to consider maps where ς = (F, ν, N ), and to linearize them about (ū,ν). The infinitesimal deformation setting implies P ≈ σ, z ≈ z e , S ≈ S e and we continue to write P , z, S in this sense. By presuming sufficient smoothness, we then get where the L (·) 's are linear forms of their entries. For the sake of simplicity we may assume that (ū,ν) are associated with a (so-called) natural state in whichP = 0, z = 0,S = 0 (that isσ = 0,z e = 0,S e = 0), so that In this case, since the measures of interaction listed above are the partial derivatives of e with respect to F , ν and N , respectively, the elastic energy is then a quadratic form in (∇ϕ,φ,∇φ). By some abuse of notation, by taking the common notation of standard linear elasticity, we then write leaving understood the meaning of the entries, with U = F −Id and Id the identity. For the sake of simplicity and without reducing drastically the generality of the physical cases covered by our treatment, the assumptions listed below apply.
(A1) The density of the elastic potential e : R d×d × R k × R k×d → [0, +∞) is a continuous and coercive quadratic form. There exist constants λ, Λ > 0 such that for every (U, ν, N ) (A2) The potential of external actions w ∈ C 1 (R d × R k , [0, +∞)) has quadratic growth, i.e. there exists Ξ 1 > 0 such that for every (u, ν) Moreover, a homogeneous boundary condition is prescribed for the displacement u on ∂B u := ∂B \ ∂B t , and for the morphological descriptor ν on a closed subset ∂B ν of ∂B. Such assumptions cover all special cases of complex bodies that we know in such linearized setting (see [15] for a list of examples).

Asynchronous Variational Integrators
The asynchronous variational integrator scheme necessary for the d'Alembert-Lagrange-type principle (2.1) can be constructed easily by following the guidelines of [9] and [10]. However, the analysis of the convergence of AVIs in presence of an abstract form of the substructural kinetic co-energy is hard. For this reason we restrict here the attention to the case in which the kinetic co-energy χ has the form where ρ is a material parameter. The general case is treated later for synchronous variational integrators.
Of course the absence of dependence of χ on ν itself implies that M is a linear space per se and that an isometric parallelism is selected over M. The Levi-Civita connection resting on a Riemaniann structure mentioned earlier is the case. Microcraked soft bodies and incommensurate intergrowth compounds fall within this scheme (see [13] and [15] for relevant examples).
In this case, we then consider the linearized dynamics of a complex body with zeroth-order dissipation governed by the d'Alembert-Lagrange-type principle which is valid for arbitrary choices of test functions ϕ, φ, defined earlier. In (4.1), A is defined by and D is the dissipation given by η is a positive parameter.
4.1. Spatial discretization. Discretization in space is obtained by means of standard finite elements given by a tessellation T of B, for instance a regular triangulation chosen to be consistent with the partition of the boundary ∂B into ∂B t and ∂B u (see, e.g., [3], [7]). By adopting a triangulation we assume that B is a polyhedral domain. The integration nodes, the generic one being indicated by a, are selected as the vertexes of the elements K of the triangulation T . We consider the space P A (T ) of continuos functions which are linear polynomials on each K ∈ T . With fixed a bounded open interval I = (t 0 , t f ) ⊆ R, we are interested in the vector subspace V u of P A (T ) ⊗ H 1 I, R d which includes all the displacement mappings satisfying u | ∂Bu ≡ 0. Any map u ∈ V u has then the form where N a ∈ P A (T ) is the nodal shape function corresponding to node a and u a (t) is the value of the displacement at a at time t. Of course we assume u a (t) ≡ 0 if a ∈ ∂B u . The null condition is not a restriction to our treatment.
With a slight abuse of notation, we denote by u (t) the vector in R D , with D the number of degrees of freedom of all nodal placements at time t, relative to the map u (·, t). Specifically, D = d#T , with #T the total number of nodes in T .
It is possible to check that there exists a constant c T > 0, depending only on the tessellation T , such that When restricted to a generic finite element K, the map u(x, t) in (4.5) is indicated by u K (x, t). Analogously u K (t) indicates a vector in R d×(d+1) having as components the displacements of all nodes of the generic element K.
We perform an analogous spatial discretization for the field ν. In this case we have ν a (t) ∈ R k for every a ∈ T , t ∈ I as a consequence of the embedding of the manifold of substructural shapes in R k . We denote by V ν the subspace of P A (T ) ⊗ H 1 I, R k which consists of all mappings ν satisfying ν | ∂Bν = 0. Furthermore, we set M = k#T .
The spatial discretization A T of A is then obtained by restricting to V T := V u × V ν the domain of A. More precisely, we consider the semi-discrete action Notice that the nodal pair (u(·), ν(·)) belongs to H 1 I, R D×M for every map (u, ν) ∈ V T . Hence, the expansion in (4.5) allows one to regard A T as a functional defined over H 1 I, R D×M . We agree to do that in the sequel in order to simplify the notation. A further simplification concerns the referential mass density ρ. In the sequel we simply presume that ρ is constant along the motion. Moreover, we assume that the mass matrix associated with the spatial discretization can be expressed in diagonal form in some basis. Then, a straightforward computation in that basis leads to The diagonalization of the mass matrix leaves invariant the Lagrangian: it can be considered, in a certain sense, as a change in observer.
We call stationary point for A T any element (u, ν) of The semi-discretized dissipation D T that appears above is defined by Moreover, in case a ∈ ∂B u \∂B ν only equation (4.10) has to be satisfied since u a (t) ≡ 0, while if a ∈ ∂B ν \∂B u we have ν a (t) ≡ 0 and the sole (4.9) has to be satisified.
Results about finite elements for complex bodies imply the theorem below (see examples in [15], [16]).
Each K ∈ T is endowed with an elemental time set which is an ordered subset Θ K of Θ. By relabeling the elements we write We assume Θ = ∪ K∈T Θ K and, for the sake of simplicity, that Θ K ∩ Θ K ′ = {t 0 , t f } for any K, K ′ ∈ T with K = K ′ . We denote by T Θ the maximum of the elemental time sizes, namely . The circumstance that each finite element can be endowed with a different time set is the basic characteristic of AVIs as already mentioned: if one imposes appropriate choices of elemental time sets, one may prove conservation of energy in discrete time. Take note that in principle such choices could not exist (see related comments in [9]).
For a node a in T , elemental time sets define also the nodal time set Θ a by As a measure of the asynchronicity of Θ, we consider the ratio and notice that 1 In the sequel we assume that each node a ∈ T follows a linear trajectory within each time interval with end-points that are consecutive instants in Θ a . Such a choice characterizes the class of AVIs that we analyze here. Then we denote by Y Θ the subspace of functions in H 1 I, R D×M which are continuous and with piecewise constant time rates in the intervals in Θ a , and such that u a ≡ 0 if a ∈ ∂B u and ν a ≡ 0 if a ∈ ∂B ν . Thus, for each (u, ν) ∈ Y Θ , a ∈ T , t i a ∈ Θ a and t ∈ t i a , t i+1 a we set Inequality (4.12) then follows by definition of τ Θ (see (4.11)).
By following [9], the fully discrete action sum in time is defined for (u, ν) ∈ Y Θ by where for the generic finite element K, A j (K; u K , ν K ) is defined by with V K given by (4.7). The choice of the approximation gives rise to explicit integrators of central-difference type and is only one of the possible schemes that can be used. It is also convenient to define all action sums on the same function space to avoid to link the function space itself to the choice of Θ. For this reason we extend A T ,Θ to the whole Y := H 1 I, R D×M . With a little abuse of notation we set The discrete stationary points of A T ,Θ are couples (u, ν) ∈ Y Θ satisfying the fully discrete analogous of the variational principle (4.8). More precisely, for all test functions (ϕ, φ) ∈ Y Θ ∩ H 1 0 (I, R D×M ) there holds (4.14) δA where the fully discretized dissipation is given by Alternatively, we will refer to the discrete stationary points also as to discrete solutions. For ν and φ as above, by a direct computation we infer The discretized form of the dissipation has been chosen to be (4.15) for the sake of simplicity. Furthermore, notice that by (4.15) and (4.16) we have In Proposition 2 below we will show that given initial conditions (u(t 0 ), ν(t 0 )) and (u(t 0 ),ν(t 0 )), the discrete d'Alembert-Lagrange-type variational principle (4.14) defines inductively a unique trajectory (u, ν) ∈ Y Θ .

Convergence in Time.
The following theorem clarifies the stage.
To prove the theorem above, preliminary results are necessary. As a first step by Lemma 1 we establish an explicit form of the discrete variational principle (4.14). From such a form, it is easy to infer existence and uniqueness of discrete solutions with given initial conditions. Then we show BV estimates for the velocities of discrete solutions in Proposition 1.
By making variations that fix all the nodes except one at a given time, one can find easily the following equivalent form of the d'Alembert-Lagrange-type principle (4.14). Lemma 1. Let (u, ν) ∈ Y Θ , then (u, ν) solves (4.14) if and only if given any node a ∈ T and any nodal time t i a ∈ (t 0 , t f ], denoting by K the sole element for which t i a ∈ Θ K and t j K = t i a , if a ∈ T \ (∂B u ∪ ∂B ν ), the following balance equations hold: Moreover, if a ∈ ∂B u \∂B ν , only the first equation has to be satisfied, while, if a ∈ ∂B ν \∂B u , only the second equation has to be satisfied.

Clearly, equations (4.19) and (4.20) define inductively a unique discrete trajectory in Y Θ once initial conditions have been specified.
Estimates on the L ∞ norm on the time interval I of the velocity of stationary points can be derived by exploiting directly the discrete stationarity conditions and the growth conditions of the potential energy densities as suggested in [11] (see also [8]). Here, we provide also an estimate on the pointwise variation of the velocities of stationary points which has been overlooked in previous works. In turn this estimate let us gain strong compactness properties and permits us to pass to the limit directly into equations (4.19) and (4.20), or equivalently in (4.14). This method is alternative to the ones already developed to study the convergence properties of variational integrators (see [10], [19], [11], [8]). Furthermore, the BV estimate wil be instrumental in Section 5 (see the related comments there).

Proposition 1.
There exist a constant T 0 = T 0 (η, ρ, λ, Λ, T ) > 0 such that given initial conditions (u(t 0 ), ν(t 0 )) and (u(t 0 ),ν(t 0 )) for every entire time set Θ with T Θ ∈ (0, T 0 ), the solution (u, ν) ∈ Y Θ of equations (4.19) and (4.20) satisfies for some constant κ > 0 depending on the initial conditions themselves and on the data of the problem. Moreover, the functionsu,ν have (pointwise) bounded variation pV on I, with Proof. Fix a node a in the tessellation T and a nodal time t i a ∈ Θ a . By assumption there exists a unique K ∈ T such that t i a ∈ Θ K , with t i a = t j K . The definition of nodal time set yields t i a = t i(α) α ∈ Θ α for all the other nodes α ∈ K. Moreover, recalling that Θ = {t l } l=1,...,NΘ , we assume t i a = t s+1 ∈ Θ for some s. We first estimate the velocity of u at time t s+1 in terms of the velocities of u itself and ν at previous times. In doing that we use the discrete d'Alembert-Lagrangetype equation (4.19), assumptions (A1)-(A3), and (4.6) which entail By taking into account that the maps u α , ν α are piecewise affine in time, we get By recalling that t i a = t s+1 , and since t i(α)−1 α ≤ t s for every α ∈ K, the latter inequality yields (4.24) u L ∞ ((t0,ts+1),R D ) ≤ cT Θ 1 + ν L ∞ ((t0,ts),R M ) + (1 + cT Θ ) u L ∞ ((t0,ts),R D ) .
We prove the analogous estimate forν, the only difference with the previous argument being that we have to take into account, in addition to the previous case, for the presence of the discretized dissipation. Thus from (4.20), assumptions (A1)-(A3) and (4.6) we get (4.25) which in turn yields Thus, there exists a constant T 0 = T 0 (η, ρ, λ, Λ, T ) > 0 such that for every T Θ ∈ (0, T 0 ) we have (1 − cT Θ ) ≥ 1/2. Finally, we deduce (4.26) In particular, by setting

Such inequality yields by iteration
Eventually, since T Θ #Θ ≤ (t f −t 0 )τ Θ #T , with τ Θ defined in (4.11), from inequality (4.27) we deduce and then (4.21) follows. In order to prove the BV estimate (4.22) we note first that sinceu a ,ν a are piecewise constant, their (pointwise) variation over an interval [t 1 , t 2 ] ⊆ I is given by the sum of the jumps in [t 1 , t 2 ]. Hence, by using (4.23) and (4.25), and taking advantage of (4.21), we obtain To compute the pointwise variation ofu a ,ν a on [t 1 , t 2 ] we sum (4.28) over the set of indices Σ t1,t2

Moreover, it holds
Proof. The convergence property stated in (4.29) follows easily from the strong convergence of ((v h , z h )) h∈N in H 1 and the regularity conditions of the integral densities of V assumed in (A1)-(A3).
. Let (δ j ) be a positive vanishing sequence. By taking into account (4.29), for every j ∈ N we get Hence, with fixed any subsequence (h i ), by a diagonal argument we find a further subsequence (h ij ) such that (4.32) Notice that by definition of first variation Moreover, for every h the action A T ,Θ h is a C 1 functional over Y Θ h by the regularity assumptions (A1)-(A3) on the densities of the potential V . Thus, for some ε j ∈ (0, δ j ) we have Actually, the same regularity assumptions imply that Thus, by (4.32)-(4.35) we infer (4.30) for the subsequence (h ij ). Urysohn property justifies (4.30) for the whole sequence.
To prove (4.31), note that the semi-discretized dissipation D T is continuous along sequences strongly converging in H 1 , so that To obtain (4.31) we need only to prove The limit above is a consequence of (4.17) (since T Θ h → 0 as h → +∞) and the strong convergences in H 1 of (z h ) h∈N and (φ h ) h∈N to z and φ, respectively. Proposition 1 and Lemma 2 are necessary tools for the proof of Theorem 2.
Proof of Theorem 2. According to (4.18) we may fix τ 0 > 0 such that for every h Denote by (u h , ν h ) the solution to equations (4.19) and (4.20), associated with Θ h .
Eventually, we show that (u, ν) solves (4.8) by passing to the limit in (4.14). Take test functions (ϕ, φ) ∈ H 1 0 I, R D×M , and let (ϕ h , φ h ) ∈ Y Θ h ∩ H 1 0 (I, R D×M ) be a sequence of piecewise affine test functions converging to (ϕ, φ) strongly in H 1 I, R D×M (see for instance [19,Lemma 4.3]). Then the discretized d'Alembert-Lagrange-type principle (4.14) is satisfied, namely The conclusion then follows by Lemma 2 and the strong convergence in H 1 of all the relevant sequences.

Remark 1.
A standard diagonalization argument implies that we can allow for intervals I = (t 0 , +∞) provided that all the statements in Theorem 2 are intended in a local sense.
Remark 2. The methods developed to prove Theorem 2 can be refined to deal with kinetic co-energies χ of the form where Ω is a positive definite symmetric element of Hom(T ν M, T * ν M). The idea is to choose for the fields u, ν different shape functions: for u the family {N a } a∈T as above, and for ν a family {N ′ a } a∈T in such a way that the mass matrix associated with the spatial discretization of χ is in diagonal form.

Synchronous Variational Integrators
Although the choices 1 2 ρν ♭ ·ν and 1 2 Ων ·ν as expressions of the substructural kinetic energy cover a wide set of special cases, there are circumstances in which more complicated expressions arise (see for example the dynamics of liquids with a dense distribution of bubbles [5]). For this reason the analysis of the dicrete schemes in which the substructural kinetic co-energy is left unspecified has its own importance. In this case we prove the convergence of the discrete dynamics only for synchronous variational integrators, i.e. when Θ K coincide with Θ for every K ∈ T .
(A5) χ(ν, ·) is uniformly convex with respect to ν, i.e. there exists a constant γ > 0 such that for all ν, ζ, ξ ∈ R k The spatial discretization follows the same lines indicated in Subsection 4.1. In the sequel we use the same notations taken there, the only difference being that the family of shape functions {N a } a∈T is chosen to form an orthonormal system in L 2 (B) (see Remark 3). A discrete version in space of the d'Alembert-Lagrange-type principle (4.1) arises. We then get where the test functions satisfy ϕ a ∈ H 1 0 (I, R d ) and ϕ a ≡ 0 if a ∈ ∂B u , and φ a ∈ H 1 0 (I, R k ) and φ a ≡ 0 if a ∈ ∂B ν . As regards the time discretization, we fix a discrete time set Θ = {t i } i=1,...,NΘ for I, with time size h Θ = max Θ (t i+1 − t i ), and consider the subspace Y Θ of H 1 (I, R D×M ) of continuous functions which are affine on each interval (t i , t i+1 ) and such that u a ≡ 0 if a ∈ ∂B u and ν a ≡ 0 if a ∈ ∂B ν . We then define the discrete action as and the potential V is defined in (4.3). Notice that in this case Θ K ≡ Θ for every K ∈ T , and thus T Θ = h Θ , τ Θ = 1. We denote by τ ′ Θ the ratio Note that we have chosen the approximation In this setting, the d'Alembert-Lagrange-type discrete principle is then given by where (ϕ, φ) ∈ Y Θ ∩ H 1 0 (I, R D×M ) and D T ,Θ is defined as in (4.15), namely The main result of the section reads as follows.
The proof follows the same strategy developed for the proof of Theorem 2, once all the preparatory steps have been accomplished. Thus we will not give the details of the proof of Theorem 3. Instead, we will focus our attention to show existence and uniqueness of discrete solutions to (5.4) and the BV estimates on their velocities. In doing that the main ingredient that we will exploit is the uniform convexity of χ(ν, ·) assumed in (A5).
In the following lemma we collect all the properties of χ which will be useful in the sequel. The proof of the lemma is standard and then omitted.

Convergence in time.
As a first step we write an equivalent form of (5.4) which can be obtained by varying only one node at each time t i ∈ Θ.
If for a fixed time t i ∈ Θ ∩ [t 0 , t f ) a solution to (5.6) and (5.7) exists for every node in T , a unique solution at time t i+1 is determined. An induction argument provides the conclusion once the initial conditions are specified.
It is clear that in this setting the variational approach based on Γ-convergence proposed in [19] cannot be pursued. In the conservative case a more conventional approach based on the repeated use of Gronwall inequality has been developed in [10].
The following proposition is the analogue of Proposition 1. Let us point out that the BV estimates proven in the sequel are instrumental. Indeed, they enable us to pass to the limit directly into the discrete equations (5.6) and (5.7) and to prove that their solutions have cluster points solving (5.1) and (5.2). Proposition 3. There exists a constant T 1 = T 1 (γ, Ξ, η, λ, Λ, T ) > 0 such that given initial conditions (u(t 0 ), ν(t 0 )) and (u(t 0 ),ν(t 0 )) for every entire time set Θ with T Θ ∈ (0, T 1 ), the solution (u, ν) ∈ Y Θ to equations (5.6) and (5.7) satisfies (5.11) u L ∞ (I,R D ) + ν L ∞ (I,R M ) ≤ κ exp (κτ ′ Θ ) , for some constant κ > 0 depending on the data of the problem and on the initial conditions themselves.
Moreover, the functionsu,ν have (pointwise) bounded variation on I with for every interval [t 1 , t 2 ] ⊆ I.
Proof. Fix a time t i ∈ Θ. To estimate the velocity of u at time t i in terms of the velocities of u itself and ν at previous times we follows the path indicated in Proposition 1. We use equation (5.6) and, for every node a ∈ T , we get (5.13) |u a (t i ) −u a (t i−1 )| ≤ cT Θ 1 + u L ∞ ((t0,ti−1),R D ) + ν L ∞ ((t0,ti−1),R M ) .
By iterating this inequality and arguing as in the proof of Proposition 1, one can deduce the L ∞ estimate (5.11). Eventually, in order to prove the BV estimate (5.12) we can reason as in the proof of Proposition 1 by taking into account (5.11), (5.13) and (5.19).
Remark 3. The orthogonality condition of the family of shape functions {N a } a∈T has been imposed only for convenience. Indeed, the same model of Section 4, in which the potential is decomposed element-by-element, can be developed for the synchronous setting, too.