1 Introduction

A wide range of publications in rigid body dynamics have been available for centuries; however, nearly all of them consider a strong form of the equations in temporal direction. Accordingly, finite difference schemes are applied leading to classical time-stepping formulations. This framework has been extensively explored for rigid bodies, introducing stabilization techniques (Baumgarte [4]), numerical diffusion in time (Hilbert et al. [16]) or mechanical integrators (Betsch [5]). It is obvious that only a few titles can be cited here, as whole journals are dedicated to this issue.

Discretization techniques of the Lagrangian itself have been taken into account (so-called variational integrators, see Lew et al. [20]), as well as temporal finite elements, see, e.g., Betsch et al. [68] for a series on finite element formulations in time. All of them are constructed in such a way that information from the past only affects the future, as this is the natural way oneself experiences time. Hence, at the end a single or multi-step method in time is derived, avoiding the solution of the whole space-time cylinder at once and the transfer of information backwards in time. However, we will show that this is not the optimal way to consider the temporal direction of constrained mechanical systems.

Discontinuous space-time Galerkin methods have been applied in, among others, Haber [1], who applied them to the whole space-time cylinder. Continuous space-time Bubnov–Galerkin strategies on deformable bodies are presented in [15, 23, 25], see also Hughes et al. [17]. In Langer et al. [19], stabilization techniques are developed for parabolic systems using fundamental consideration of Ladyzhenskaya [18]. For inverse problems, Ströhle et al. [27] demonstrate the superiority of continuous space-time Galerkin methods.

In this contribution, we extend the concept of a continuous space-time Bubnov–Galerkin method towards rigid bodies. To be specific, we derive the formulation directly from the continuous case of nonlinear elastodynamics and constrain the body to obtain in the limit of infinitely large material parameters a set of differential algebraic equations (DAE’s) for the rigid body. To avoid unconvincing complexity, especially in the case of multiple connected rigid bodies, we circumvent the use of rotational variables to describe the orientation of the individual rigid bodies, see Betsch et al. [10] for a more detailed discussion. In this sense, we make use of a rotationless formulation by parametrization of the rotation matrix \(\boldsymbol{R}\), i.e., an element of the Lie group \(SO(n)\) with \(\boldsymbol{R}^{T}\boldsymbol{R} = \boldsymbol{I}\) and \(\operatorname{det}(\boldsymbol{R}) = 1\), denoting all matrices with orthonormal column vectors and positive determinant. This simplifies the design of a space-time formulation dramatically and allows for the design of energy-momentum schemes.

Hence, we introduce a constrained director-based kinematic, see Rubin [22] for details on the theory of Cosserat points. As usual for rigid body dynamics, the spatial integral is executed independently of the temporal integration, and thus we make use of a tensor–product structure of the shape functions in space and time. Analogous to Mortar methods in domain decomposition, a weak form of the constraints is considered, i.e., we interpolate the Lagrange multipliers and their corresponding test or trial functions in temporal direction and implement the time continuous weak form. We consider this approach as the most natural way of dealing with Hamilton’s principle. As we will show, we obtain optimal convergence and allow for an efficient, parallel solution of the space-time problem.

The manuscript is organized as follows. In Sect. 2, we introduce the basic space-time formulation, followed by the discretization in Sect. 3. Section 4 demonstrates the superiority of the chosen approach and the applicability of large systems including large rotations. Conclusions are drawn in Sect. 5.

2 Variational formulation

In this section, we develop basic concepts for a space-time formulation on constrained mechanical systems. The motion of a mechanical system is described by Hamilton’s principle, namely, we make the action integral stationary

$$ \delta \mathfrak{L} = 0,\quad \text{with}\quad \delta \mathfrak{L} = \int _{\mathfrak{B}_{0}} \delta L \,\mathrm{d}W. $$
(1)

Here, \(L := L(\boldsymbol{\varphi },\nabla _{t}\boldsymbol{\varphi })\) is the Lagrangian given as \(L = T-V\), where \(T:=T(\nabla _{t}\boldsymbol{\varphi })\) represents the specific kinetic energy, \(\nabla _{t}(\bullet )\) is the gradient in temporal direction, and \(V:=V(\boldsymbol{\varphi })\) denotes a suitable specific potential energy. Note that we do not require the existence of \(V\), only the variation of V. The Lagrangian is defined within a space-time cylinder \(\mathfrak{B}_{0}:=\Omega _{0}\times \mathcal{T}\subset \mathbb{R}^{d+1}\) with regard to its spatial \(d\)-dimensional reference configuration in the open bounded domain \(\Omega _{0}\subset \mathbb{R}^{d}\) within the time period \(\mathcal{T}:=[0,\,T^{e}]\subset \mathbb{R}\). Accordingly, we can define the mantle of the space-time cylinder \(\Sigma = \partial \Omega _{0}\times \mathcal{T}\) and the start and end configuration \(\Sigma ^{0} = \Omega _{0}\times \{0\}\) and \(\Sigma ^{T^{e}} = \Omega _{0}\times \{T^{e}\}\), respectively. As usual, the mantle can be further decomposed into non-overlapping Dirichlet boundaries \(\Sigma ^{u}\) and Neumann boundaries \(\Sigma ^{\sigma}\). Note that we make use of a space-time deformation map \(\boldsymbol{\varphi }:\mathfrak{B}_{0}\rightarrow \mathfrak{B}\), where \(\mathfrak{B}\) is the actual configuration, and \(\boldsymbol{\varphi }:=(\boldsymbol{\varphi }^{T}(\boldsymbol{X},t),t)^{T}\), where the last component is fixed in classical mechanics.

Following the theory of Cosserat points (see Rubin [22]), the deformation map is defined in terms of the mapping \(\bar{\boldsymbol{\varphi }}\) of a reference material point \(\bar{\boldsymbol{X}}\) along with a time-dependent director coordinate frame \(\{\boldsymbol{d}_{i}(t)\},\,i\in [1,\,\ldots ,\,d]\). A material point \(\boldsymbol{X}\in \Omega _{0}\) of the rigid body is given by

$$ \boldsymbol{X} (\boldsymbol{\theta })= \bar{\boldsymbol{X}} + \sum _{i}\theta ^{i}\boldsymbol{d}_{i}(0), \quad \boldsymbol{\theta }=\left [\theta ^{1} \ldots \theta ^{d}\right ]^{ \text{T}}, $$
(2)

where we have made use of the convective coordinates \(\theta ^{i}\). Accordingly, \(\boldsymbol{\varphi }\in \mathfrak{B}\) at \(t\in \mathcal{T}\) reads

$$ \boldsymbol{\varphi } (\boldsymbol{\theta },t) = \bar{\boldsymbol{\varphi }}(\bar{\boldsymbol{X}},t) + \sum _{i} \theta ^{i}\boldsymbol{d}_{i}(t). $$
(3)

A symmetric tensor-valued kinematic constraint \(\boldsymbol{\Phi }^{rb}(t)\) must be established to prevent deformation of the director triad with coefficients

$$ \Phi ^{rb}_{ij}(t):=\frac{1}{2}\left (\boldsymbol{d}_{i}(t)\cdot \boldsymbol{d}_{j}(t) - \boldsymbol{d}_{i}(0)\cdot \boldsymbol{d}_{j}(0) \right ) = 0;\quad i,j\in [1,\,\ldots ,\,d]. $$
(4)

The satisfaction of the kinematic constraint (4) leads to the introduction of a symmetric tensor-valued Lagrange multiplier field \(\boldsymbol{\lambda }_{rb}(t)\). Alternatively, \(\boldsymbol{\lambda }_{rb}(t)\) can be understood as a reactive second Piola–Kirchhoff stress tensor introduced to ensure that the metric tensor defined by the director triad remains constant throughout the deformation. The specific kinetic energy is usually written in classical mechanics as follows:

$$ T = \frac{1}{2}\rho _{0}\nabla _{t}\boldsymbol{\varphi }\cdot \nabla _{t} \boldsymbol{\varphi },\quad \text{and}\quad \boldsymbol{\pi } = \frac{\partial T(\nabla _{t}\boldsymbol{\varphi })}{\partial \nabla _{t}\boldsymbol{\varphi }}= \rho _{0}\nabla _{t}\boldsymbol{\varphi }, $$
(5)

where \(\rho _{0}\in \mathbb{R}\) is the reference density and \(\boldsymbol{\pi }\) is the linear momentum. Inserting the rigid body kinematic in (3) and executing spatial integration independent of the temporal integration, we obtain for the virtual work of the inertial forces

$$ \int _{\mathfrak{B}_{0}} \delta T \,\mathrm{d}W = \int _{\mathcal{T}} \begin{bmatrix} \nabla _{t}\delta \bar{\boldsymbol{\varphi }} \\ \nabla _{t}\delta \boldsymbol{d}^{T}_{i} \end{bmatrix} \cdot \underbrace{\mathbb{M}^{ij} \begin{bmatrix} \nabla _{t}\bar{\boldsymbol{\varphi }}\\ \nabla _{t}\boldsymbol{d}_{j} \end{bmatrix}}_{=[\bar{\boldsymbol{\pi }},\,\boldsymbol{\pi }_{i}^{d}]^{T}} \,\mathrm{d}t, $$
(6)

where

$$ \begin{aligned} \mathbb{M}^{ij} &= \begin{bmatrix} M_{\varphi}\boldsymbol{I} & S^{j}\boldsymbol{I} \\ S^{i}\boldsymbol{I} & E^{ij}\boldsymbol{I} \end{bmatrix} ,\quad M_{\varphi }= \int _{\Omega _{0}} \rho _{0} \,\mathrm{d}V, \\ S^{i} &= \int _{\Omega _{0}} \rho _{0}\theta ^{i} \,\mathrm{d}V, \quad \text{and} \quad E^{ij} = \int _{\Omega _{0}} \rho _{0}\theta ^{i} \theta ^{j} \,\mathrm{d}V \end{aligned} $$
(7)

introduce a symmetric positive-definite mass matrix \(\mathbb{M}\in \mathbb{R}^{(d+1)d\times (d+1)d}\) and the unity matrix \(\boldsymbol{I}\in \mathbb{R}^{d\times d}\).

Remark 2.1

Using a coordinate system placed in the center of mass of \(\Omega _{0}\) yields obviously \(S^{i} = 0\). Additionally, the coordinate axes are often adjusted along the principal axis of the tensor \(\boldsymbol{E} = E^{ij}\boldsymbol{d}_{i}\otimes \boldsymbol{d}_{j}\). Moreover, the Euler tensor can be linked to the common inertia tensor via the relationship \(\boldsymbol{J} = \operatorname{tr}(\boldsymbol{E}) \boldsymbol{I}-\boldsymbol{E}\).

Assuming that the kinetic energy \(T(\nabla _{t}\boldsymbol{\varphi })\) is convex and differentiable, we can introduce the conjugate function (see Appendix A for details)

$$ T^{*}(\boldsymbol{\pi }) = \underbrace{\operatorname{sup}}_{\nabla _{t} \boldsymbol{\varphi }}(\boldsymbol{\pi }\cdot \nabla _{t} \boldsymbol{\varphi }-T(\nabla _{t}\boldsymbol{\varphi })) $$
(8)

and obtain \(T^{*}(\boldsymbol{\pi }) = \frac{1}{2}\rho _{0}^{-1}\| \boldsymbol{\pi }\|^{2}\). Adapting the rigid body kinematic, we conclude that

$$ \int _{\mathfrak{B}_{0}} \delta T^{*} \,\mathrm{d}W = \frac{1}{2}\int _{\mathcal{T}} \begin{bmatrix} \delta \bar{\boldsymbol{\pi }} \\ \delta \boldsymbol{\pi }^{d}_{j} \end{bmatrix} \cdot \left [\mathbb{M}^{-1}\right ]^{ij} \begin{bmatrix} \bar{\boldsymbol{\pi }} \\ \boldsymbol{\pi }^{d}_{j} \end{bmatrix} \,\mathrm{d}t $$
(9)

and obtain

$$ \begin{aligned} \int \limits _{\mathcal{T}} \begin{bmatrix} \delta \bar{\boldsymbol{\pi }} \\ \delta \boldsymbol{\pi }^{d}_{i} \end{bmatrix} \cdot \left ( \begin{bmatrix} \nabla _{t}\bar{\boldsymbol{\varphi }} \\ \nabla _{t}\boldsymbol{d}_{i} \end{bmatrix} - \left [\mathbb{M}^{-1}\right ]^{ij} \begin{bmatrix} \bar{\boldsymbol{\pi }} \\ \boldsymbol{\pi }^{d}_{j} \end{bmatrix} \right )\,\mathrm{d}t =&\, 0, \\ \int \limits _{\mathcal{T}} - \begin{bmatrix} \nabla _{t}\delta \bar{\boldsymbol{\varphi }} \\ \nabla _{t}\delta \boldsymbol{d}^{T}_{i} \end{bmatrix} \cdot \begin{bmatrix} \bar{\boldsymbol{\pi }} \\ \boldsymbol{\pi }^{d}_{i} \end{bmatrix} + \begin{bmatrix} \delta \bar{\boldsymbol{\varphi }} \\ \delta \boldsymbol{d}^{T}_{i} \end{bmatrix} \cdot \bar{\lambda}^{kl}_{rb} \begin{bmatrix} \boldsymbol{0} \\ \frac{\partial \Phi ^{rb}_{kl}}{\partial \boldsymbol{d}_{i}} \end{bmatrix} \,\mathrm{d}t =& \\ \int \limits _{\mathfrak{B}_{0}}\delta \boldsymbol{\varphi }\cdot \boldsymbol{B}\,\mathrm{d}W + \int \limits _{\Sigma ^{ \sigma}}\delta \boldsymbol{\varphi }&\cdot \boldsymbol{T} \,\mathrm{d}A\,\mathrm{d}t - \begin{bmatrix} \delta \bar{\boldsymbol{\varphi }} \\ \delta \boldsymbol{d}^{T}_{i} \end{bmatrix} \cdot \begin{bmatrix} \bar{\boldsymbol{\pi }} \\ \boldsymbol{\pi }^{d}_{i} \end{bmatrix} \biggl|_{0}^{T^{e}}, \\ \int \limits _{\mathcal{T}}\delta \bar{\boldsymbol{\lambda }}_{rb}: \boldsymbol{\Phi }^{rb}\,\mathrm{d}t =& \,0, \end{aligned} $$
(10)

where we have multiplied the second equation by −1 for the ease of presentation.

The last term in the second equation refers to the incoming and outgoing linear momentum of the time interval under consideration, called by Hamilton the “law of varying action”. Often \(\delta \boldsymbol{\varphi }|_{0}^{T}\) is set to zero so that the initial and end terms vanish. It is obvious that for classical engineering tasks this is no useful assumption as we do not know the end position in advance; however, for a classical time-stepping scheme, this term does not matter at all. Applying a space-time concept, this term compares to the Neumann boundary conditions in space, either predefined or set to zero. In temporal direction, the incoming and outgoing linear momentum may neither be predefined nor zero, hence this term does not vanish, see Appendix B and also Bailey [2, 3] for further details.

Eventually, we set \(\bar{\boldsymbol{\lambda }}_{rb} = \boldsymbol{\lambda }_{rb}V\), scaled by the volume \(V\) of \(\Omega _{0}\), \(\boldsymbol{B}\) denotes a body force per unit of mass and \(\boldsymbol{T}\) is the traction vector acting on the mantle \(\Sigma ^{\sigma}\) of the space-time cylinder. Application of integration by parts in temporal direction of the first term in (10) yields

$$ \begin{aligned} \int \limits _{\mathcal{T}} \begin{bmatrix} \delta \bar{\boldsymbol{\pi }} \\ \delta \boldsymbol{\pi }^{d}_{i} \end{bmatrix} \cdot \left ( \begin{bmatrix} \nabla _{t}\bar{\boldsymbol{\varphi }} \\ \nabla _{t}\boldsymbol{d}_{i} \end{bmatrix} - \left [\mathbb{M}^{-1}\right ]^{ij} \begin{bmatrix} \bar{\boldsymbol{\pi }} \\ \boldsymbol{\pi }^{d}_{j} \end{bmatrix} \right )\,\mathrm{d}t =& 0, \\ \int \limits _{\mathcal{T}} \begin{bmatrix} \delta \bar{\boldsymbol{\varphi }} \\ \delta \boldsymbol{d}_{i} \end{bmatrix} \cdot \begin{bmatrix} \nabla _{t}\bar{\boldsymbol{\pi }} \\ \nabla _{t}\boldsymbol{\pi }^{d}_{i} \end{bmatrix} + \begin{bmatrix} \delta \bar{\boldsymbol{\varphi }} \\ \delta \boldsymbol{d}_{i} \end{bmatrix} \cdot \bar{\lambda}^{kl}_{rb} \begin{bmatrix} \boldsymbol{0} \\ \frac{\partial \Phi ^{rb}_{kl}}{\partial \boldsymbol{d}_{i}} \end{bmatrix} \,\mathrm{d}t =& \int \limits _{\mathfrak{B}_{0}}\delta \boldsymbol{\varphi }\cdot \boldsymbol{B}\,\mathrm{d}W + \int \limits _{\Sigma ^{\sigma}}\delta \boldsymbol{\varphi }\cdot \boldsymbol{T}\,\mathrm{d}A\,\mathrm{d}t, \\ \int \limits _{\mathcal{T}}\delta \bar{\boldsymbol{\lambda }}_{rb}: \boldsymbol{\Phi }^{rb}\,\mathrm{d}t =& 0. \end{aligned} $$
(11)

The spaces of virtual or admissible test and solution functions for (11) are

$$\begin{aligned} \mathcal{V}^{\varphi }&= \{\boldsymbol{\varphi }\in H^{1}(\mathcal{T}) | \boldsymbol{\varphi } = \bar{\boldsymbol{\varphi }}( \bar{\boldsymbol{X}},t) + \sum _{i}\theta ^{i}\boldsymbol{d}_{i}(t) \;\text{and}\; \boldsymbol{\Phi }^{rb} = 0\}, \end{aligned}$$
(12)
$$\begin{aligned} \mathcal{V}^{\delta \varphi} &= \{\delta \boldsymbol{\varphi }\in L^{2}( \mathcal{T}) | \delta \boldsymbol{\varphi } = \boldsymbol{0}\; \text{on}\;\Sigma ^{\varphi }\cup \Sigma ^{0} \}, \end{aligned}$$
(13)
$$\begin{aligned} \mathcal{V}^{\pi }&= \{\boldsymbol{\pi }\in H^{1}(\mathcal{T})\}, \end{aligned}$$
(14)
$$\begin{aligned} \mathcal{V}^{\delta \pi} &= \{\delta \boldsymbol{\pi }\in L^{2}( \mathcal{T}) | \delta \boldsymbol{\varphi } = \boldsymbol{0}\; \text{on}\;\Sigma ^{\pi }\cup \Sigma ^{0}\}. \end{aligned}$$
(15)

Here, \(H^{1}\) denotes the Sobolev functional space of square-integrable functions and derivatives. Note that we applied initial (essential) boundary conditions on \(\mathcal{V}^{\delta \varphi}\) and \(\mathcal{V}^{\delta \pi}\) for the initial position and initial linear momentum. However, we can also apply boundary conditions on the endpoint (e.g., the crash has happened and one wants to calculate where the car comes from) by simply applying the Dirichlet-like conditions to \(\Sigma ^{T}\) instead of \(\Sigma ^{0}\).

3 Space-time discretization

Next, we introduce the space-time discretization, following the ideas developed in Hesch et al. [13]. Although the general methodology developed in Schuss et al. [25] allows for unstructured meshes in the space-time cylinder, we make use of a tensor product structure to integrate the spatial domain independent of the temporal. The reference configuration \(\Omega _{0}\) is split into a subdivision of non-overlapping, spatial elements \(e \in \mathbb{E}\), leading to an approximate domain \(\Omega _{0}^{h}\), which satisfies

$$ \Omega _{0}^{h} = \bigcup _{e\in \mathbb{E}}\Omega _{0,e}, \quad \forall e \in \mathbb{E}\quad \text{and}\quad \mathcal{T}^{h} = \bigcup _{e_{t}\in \mathbb{T}}\mathcal{T}_{e_{t}}, \quad \forall e_{t} \in \mathbb{T}, $$
(16)

using non-overlapping elements \(e_{t}\in \mathbb{T}\) for the 1D subdivision of the temporal domain \(\mathcal{T}\). A finite-dimensional polynomial approximation of the test and trial functions for the unknown event field in the space-time domain is defined as follows:

$$ \begin{aligned} \boldsymbol{\varphi }^{h}(\boldsymbol{X},t) = \sum \limits _{A\in \omega} \left (M_{t}(t)N_{\varphi}(\boldsymbol{X})\right )^{A} \boldsymbol{q}_{A}(t),\quad \delta \boldsymbol{\varphi }^{h} ( \boldsymbol{X},t)= \sum \limits _{B\in \omega}\left (M_{t}(t)N_{ \varphi}(\boldsymbol{X})\right )^{B}\delta \boldsymbol{q}_{B}(t), \\ \boldsymbol{\pi }^{h}(\boldsymbol{X},t) = \sum \limits _{A\in \omega} \left (M_{t}(t)N_{\varphi}(\boldsymbol{X})\right )^{A} \boldsymbol{\pi }_{A}(t),\quad \delta \boldsymbol{\pi }^{h} ( \boldsymbol{X},t)= \sum \limits _{B\in \omega}\left (M_{t}(t)N_{ \varphi}(\boldsymbol{X})\right )^{B}\delta \boldsymbol{\pi }_{B}(t), \end{aligned} $$
(17)

where \(\left (M_{t}(t)N_{\varphi}(\boldsymbol{X})\right )^{A}:\mathfrak{B}_{0} \rightarrow \mathbb{R}\) are temporal and spatial shape functions associated with nodes where \(\omega = \omega _{s}\otimes \omega _{t}\), the tensor product of the set of spatial nodes \(\omega _{s}\) and the set of temporal nodes \(\omega _{t}\). Note that we require the partition of unity for the spatial and temporal shape functions. As a rigid body kinematic allows to separately solve the integral over space and time, we separate the multi-index \(A\) into \(A_{s}\) and \(A_{t}\) for the spatial and temporal shape functions, such that \(\left (M_{t}(t)N_{\varphi }(\boldsymbol{X})\right )^{A} = M_{t}^{A_{t}}(t)N_{ \varphi }^{A_{s}}(\boldsymbol{X})\). The nodal values \(\boldsymbol{q}_{A}\), \(\boldsymbol{\pi }_{A}\) and \(\delta \boldsymbol{q}_{B}\), \(\delta \boldsymbol{\pi }_{B}\) admit the following rigid body representation (refer to equation (3)):

$$ \begin{aligned} \boldsymbol{q}_{A}(\boldsymbol{\theta },t) = \bar{\boldsymbol{\varphi }}_{A_{t}} + \sum \limits _{i}\theta ^{i}_{A_{s}} \boldsymbol{d}_{i,A_{t}},\quad \delta \boldsymbol{q}_{B}( \boldsymbol{\theta },t) =\delta \bar{\boldsymbol{\varphi }}_{B_{t}} + \sum \limits _{i}\theta ^{i}_{B_{s}}\delta \boldsymbol{d}_{i,B_{t}}, \\ \boldsymbol{\pi }_{A}(\boldsymbol{\theta },t) = \bar{\boldsymbol{\pi }}_{A_{t}} + \sum \limits _{i}\theta ^{i}_{A_{s}} \boldsymbol{\pi }^{d}_{i,A_{t}},\quad \delta \boldsymbol{\pi }_{B}( \boldsymbol{\theta },t) =\delta \bar{\boldsymbol{\varphi }}_{B_{t}} + \sum \limits _{i}\theta ^{i}_{B_{s}}\delta \boldsymbol{\pi }^{d}_{i,B_{t}}, \end{aligned} $$
(18)

where \(\theta ^{i}_{A_{s}}\) is the \(i\)th convective coordinate of the node \(A_{s}\) measured with respect to the reference material point \(\bar{\boldsymbol{X}}\). Insertion into the virtual work expressions (6) or (9), respectively, allows us to rewrite for the components of the discretized mass matrix \(\mathbb{M}_{h}\) given in (7) in terms of a finite element discretization of the spatial reference configuration as follows:

$$ \begin{aligned} M_{\bar{\varphi }} &= \sum \limits _{e\in \mathbb{E}}\left (\,\int \limits _{\Omega ^{h}_{0,e}}\rho _{0}\,\mathrm{d}V \right ),\quad S^{i} = \sum \limits _{e\in \mathbb{E}}\left (\,\int \limits _{\Omega ^{h}_{0,e}}\rho _{0} \sum \limits _{A_{s}\in \omega _{s}} \left (N^{A_{s}}_{\varphi }\theta ^{i}_{A_{s}}\right ) \,\mathrm{d}V\right ), \\ E^{ij} &= \sum \limits _{e\in \mathbb{E}}\left (\,\int \limits _{ \Omega ^{h}_{0,e}}\rho _{0} \sum \limits _{A_{s},B_{s}\in \omega _{s}} \left ( N^{A_{s}}_{\varphi }N^{B_{s}}_{\varphi }\theta ^{i}_{A_{s}} \theta ^{j}_{B_{s}}\right )\,\mathrm{d}V\right ). \end{aligned} $$
(19)

Eventually, the spatial discretized external contributions can be written as follows:

$$ \begin{aligned} &\int \limits _{\mathfrak{B}^{h}_{0}}\delta \boldsymbol{\varphi }^{h} \cdot \boldsymbol{B}\,\mathrm{d}W + \int \limits _{ \Sigma ^{\sigma ,h}_{0}}\delta \boldsymbol{\varphi }^{h}\cdot \boldsymbol{T}\,\mathrm{d}A\,\mathrm{d}t = \\ &\sum \limits _{e_{t}\in \mathbb{T}}\left (\,\int \limits _{ \mathcal{T}_{e_{t}}^{h}}\sum \limits _{A_{t}\in \omega _{t}} \delta \bar{\boldsymbol{\varphi }}_{A_{t}}\cdot M^{A_{t}}_{t} \sum \limits _{e \in \mathbb{E}}\left (\,\int \limits _{\Omega _{0,e}^{h}} \sum \limits _{B_{s}\in \omega _{s}}\left (N^{B_{s}}_{\varphi }\right ) \boldsymbol{B}\,\mathrm{d}V +\right.\right. \\ &\left.\left. \int \limits _{\partial \Omega _{0,e}} \sum \limits _{B_{s}\in \omega _{s}}\left (N^{B_{s}}_{ \varphi }\right )\boldsymbol{T}\,\mathrm{d}A\right ) \,\mathrm{d}t\right ) + \\ &\sum \limits _{e_{t}\in \mathbb{T}}\left (\,\int \limits _{ \mathcal{T}_{e_{t}}^{h}}\sum \limits _{A_{t}\in \omega _{t}}\delta \boldsymbol{d}_{i,A_{t}}\cdot M_{t}^{A_{t}} \sum \limits _{e\in \mathbb{E}}\left (\,\int \limits _{\Omega _{0,e}^{h}} \sum \limits _{B_{s} \in \omega _{s}}\left (N^{B_{s}}_{\varphi }\theta ^{i}_{B_{s}}\right ) \boldsymbol{B}\,\mathrm{d}V + \right.\right. \\ &\left.\left. \int \limits _{\partial \Omega _{0,e}}\sum \limits _{B_{s}\in \omega _{s}} \left (N^{B_{s}}_{ \varphi }\theta ^{i}_{B_{s}}\right )\boldsymbol{T} \,\mathrm{d}A\right )\,\mathrm{d}t \right ). \end{aligned} $$
(20)

Note that we obtain for an external torque vector \(\boldsymbol{m}\) the relation

$$ \int _{\Sigma ^{\sigma ,h}_{0}}\delta \boldsymbol{\varphi }^{h} \cdot \boldsymbol{T}\,\mathrm{d}A \,\mathrm{d}t = \sum _{e_{t}\in \mathbb{T}} \left (\int _{\mathcal{T}_{e_{t}}^{h}} \sum _{A_{t} \in \omega _{t}}\delta \boldsymbol{d}_{i,A_{t}}\cdot M_{t}^{A_{t}} \frac{1}{2} \left (\boldsymbol{m}\times \sum _{B_{t}\in \omega _{t}}M_{t}^{B_{t}}\boldsymbol{d}^{i}_{B_{t}}\right ) \,\mathrm{d}t\right ), $$
(21)

where \(\boldsymbol{d}^{i} = ([\boldsymbol{d}_{i}\cdot \boldsymbol{d}_{j}])^{-1} \boldsymbol{d}_{j}\) is a contravariant director. We refer to [12] for details on director-based formulations written in skew coordinates and to [24], Sect. 2.2.3 for the application on external torque. Moreover, the Lagrange multipliers to enforce the orthonormality of the director triad are approximated via

$$ \bar{\lambda }^{kl}_{rb} = \sum _{A_{t}\in \omega _{t}}M_{t}^{A_{t}} \bar{\lambda }^{kl}_{A_{t}}, \quad \delta \bar{\lambda }^{kl}_{rb} = \sum _{B_{t}\in \omega _{t}}M_{t}^{B_{t}}\delta \bar{\lambda }^{kl}_{B_{t}}. $$
(22)

For conservation of angular momentum, see Appendix C.

Remark 3.1

Using the same shape functions for the Lagrange multipliers (and the associated trial functions) as for the solution and trial functions of the rigid body is possible, although other choices can be applied as well, see Brivadis et al. [11] for a detailed discussion on this issue in the spatial domain for NURBS. For the ease of construction, we make use of the same functions, simplifying the whole concept and the implementation of finite elements for rigid bodies in space-time. As we will show in the following section, the convergence remains optimal.

Finally, we obtain the following algebraic system of equations:

$$ \begin{aligned}[c] [\delta \bar{\boldsymbol{\pi }}_{A_{t}},\, \delta \boldsymbol{\pi }^{d}_{i,A_{t}}] \cdot \int \limits _{\mathcal{T}} M^{A_{t}}_{t} \left ( \begin{bmatrix} \nabla _{t}(M^{B_{t}}_{t}) \bar{\boldsymbol{\varphi }}_{B_{t}} \\ \nabla _{t}(M^{B_{t}}_{t})\boldsymbol{d}_{i,B_{t}} \end{bmatrix} - \left [\mathbb{M}_{h}^{-1}\right ]^{ij} \begin{bmatrix} M^{B_{t}}_{t}\bar{\boldsymbol{\pi }}_{B_{t}} \\ M^{B_{t}}_{t}\boldsymbol{\pi }^{d}_{j,B_{t}} \end{bmatrix} \right )\,\mathrm{d}t &= 0, \\ {}[\delta \bar{\boldsymbol{\varphi }}_{A_{t}},\, \delta \boldsymbol{d}_{i,A_{t}}]\cdot \int \limits _{\mathcal{T}} M^{A_{t}}_{t} \left ( \begin{bmatrix} \nabla _{t}(M^{B_{t}}_{t})\bar{\boldsymbol{\pi }}_{B_{t}} \\ \nabla _{t}(M^{B_{t}}_{t})\boldsymbol{\pi }^{d}_{i,B_{t}} \end{bmatrix} + M^{B_{t}}_{t}\bar{\lambda }^{kl}_{B_{t}} \begin{bmatrix} \boldsymbol{0} \\ \left (\frac{\partial \Phi ^{rb}_{kl}}{\partial \boldsymbol{d}_{i}} \right )^{h} \end{bmatrix} \right ) \,\mathrm{d}t & \\ -\int \limits _{\mathfrak{B}_{0}}\delta \boldsymbol{\varphi }^{h} \cdot \boldsymbol{B}\,\mathrm{d}W - \int \limits _{ \Sigma ^{\sigma }}\delta \boldsymbol{\varphi }^{h}\cdot \boldsymbol{T}\,\mathrm{d}A \,\mathrm{d}t &= 0, \\ \delta \bar{\boldsymbol{\lambda }}_{A_{t}}: \int \limits _{ \mathcal{T}}M_{t}^{A_{t}}\boldsymbol{\Phi }^{h} \,\mathrm{d}t &= 0, \end{aligned} $$
(23)

where

$$ \boldsymbol{\Phi }^{h} = \frac{1}{2}\left (M_{t}^{B_{t}} \boldsymbol{d}_{i,B_{t}} (t)\cdot M_{t}^{B_{t}}\boldsymbol{d}_{j,B_{t}}(t) - M_{t}^{B_{t}} \boldsymbol{d}_{i,B_{t}}(0) \cdot M_{t}^{B_{t}} \boldsymbol{d}_{j,B_{t}}(0)\right ), $$
(24)

and \(\left (\frac{\partial \Phi ^{rb}_{kl}}{\partial \boldsymbol{d}_{i}} \right )^{h}\) follows in a straightforward manner from the last equation.

Remark 3.2

To solve the nonlinear system of equations (23), we apply a Newton–Raphson iteration. In particular, we obtain a residual vector \(\boldsymbol{R}_{A_{t}}\), derive the tangent matrix \(K_{A_{t},B_{t}}\), and solve

$$ \begin{aligned} K_{A_{t},B_{t}}^{n} \begin{bmatrix} \Delta \bar{\boldsymbol{\varphi }}_{B_{t}} \\ \Delta \boldsymbol{d}_{i,B_{t}} \\ \Delta \bar{\boldsymbol{\pi }}_{B_{t}} \\ \Delta \boldsymbol{\pi }^{d}_{i,B_{t}} \\ \Delta \bar{\boldsymbol{\lambda }}_{B_{t}} \end{bmatrix} &= -\boldsymbol{R}_{A_{t}}^{n},\quad \text{and}\\ \begin{bmatrix} \bar{\boldsymbol{\varphi }}_{A_{t}} \\ \boldsymbol{d}_{i,A_{t}} \\ \bar{\boldsymbol{\pi }}_{A_{t}} \\ \boldsymbol{\pi }^{d}_{i,A_{t}} \\ \bar{\boldsymbol{\lambda }}_{A_{t}} \end{bmatrix} ^{n+1}& = \begin{bmatrix} \bar{\boldsymbol{\varphi }}_{A_{t}} \\ \boldsymbol{d}_{i,A_{t}} \\ \bar{\boldsymbol{\pi }}_{A_{t}} \\ \boldsymbol{\pi }^{d}_{i,A_{t}} \\ \bar{\boldsymbol{\lambda }}_{A_{t}} \end{bmatrix} ^{n} + \begin{bmatrix} \Delta \bar{\boldsymbol{\varphi }}_{A_{t}} \\ \Delta \boldsymbol{d}_{i,A_{t}} \\ \Delta \bar{\boldsymbol{\pi }}_{A_{t}} \\ \Delta \boldsymbol{\pi }^{d}_{i,A_{t}} \\ \Delta \bar{\boldsymbol{\lambda }}_{A_{t}} \end{bmatrix} , \quad \forall A_{t},B_{t}\in \omega _{t}, \end{aligned} $$
(25)

until a predefined stop criterion \(\sum _{A_{t}\in \omega _{t}}\|\boldsymbol{R}_{A_{t}}^{n+1}\| < \epsilon \) is obtained.

4 Results

In this section, we provide a series of numerical examples to demonstrate the accuracy and superiority of the space-time system compared to classical time-stepping schemes. Therefore, we assume that the reader is familiar with time-stepping schemes like mid- and endpoint rules (i.e., implicit Euler scheme), see [9] for details on time-stepping schemes applied to director-based rigid body formulations.

4.1 Single rigid body

We start with a prototypical 3D example to conduct suitable convergence tests. Therefore, a rectangular rigid body with dimensions \([-1,\, 1]\), \([-0.5,\, 0.5]\), and \([-0.5,\, 0.5]\) is introduced with reference density of \(\rho _{0} = 1.5\). The initial velocity of the center of mass is \([5,\,5,\,5]\); the (rotational) initial velocity of the directors is zero. Two different tests are carried out: First, a simple test with constant external momentum of \(\boldsymbol{m} = [1,\,1,\,1]\) throughout the whole time interval \(\mathcal{T}^{\operatorname{rot}} = [0,\,1]\) is applied. Note that we make use of equation (21) to apply the external momentum \(\boldsymbol{m}\) on the director formulation, i.e., \(\boldsymbol{m}\) is written in the global coordinate system, mapped to the director frame via (21). Figure 1, left displays the motion of the body, and although we call this example “simple”, the body performs a nearly \(90^{\circ}\) degree rotation within the time interval, hence large rotations are involved.

Fig. 1
figure 1

Motion of a single body. Left: Constant load. Right: Rectangular load

Figure 2 shows the results of a convergence test using the approach as presented in (11). For this demonstration, we made use of an overkill solution to obtain \(\boldsymbol{\varphi }^{ref}(T)\) with 1 million time steps using a mid-point rule and compared the norm of the final position and the directors minus the overkill solution

$$ Err = \|\bar{\boldsymbol{\varphi }}(\bar{\boldsymbol{X}},T) - \bar{\boldsymbol{\varphi }}^{ref}(\bar{\boldsymbol{X}},T)\| + \sum _{i = 1}^{3} \|\boldsymbol{d}_{i}(T) - \boldsymbol{d}^{ref}_{i}(T) \|. $$
(26)

As can be seen, the implicit Euler “Timestepping endpoint”) converges linearly, whereas the midpoint rule (“Timestepping midpoint”) converges quadratically. Linear finite elements in space-time (“Space-time linear”) converge quadratic, which was expected; using quadratic Lagrange shape functions (“Space-time quad”) converges cubic as expected up to \(10^{-10}\), which was the predefined stop criterion of the Newton–Raphson iteration. Hence, the space-time formulation converges near machine precision towards the results of the overkill time-stepping approach.

Fig. 2
figure 2

Convergence test for constant load. \(Err\) is plotted over the element size h

A direct comparison of the runtime is difficult as modern computers scale the frequency of the CPU unpredictably depending on the temperature, peak values of the power consumption, the implementation, and many more factors. Just to give an impression: On a six-core MacBook (2.6 GHz Intel Core i7), the space-time simulation using 1000 linear finite elements and a direct solver (UMFPACK) with METIS reordering took about 12.5 seconds. The midpoint time-stepping scheme used 92.7 seconds. On our 24 core Workstation, 1000 linear finite elements took even longer as the overhead for this small number of elements is too large. However, for larger systems up to 100,000 elements, the multi-core workstations scale up as expected without any problems, resulting in a similar or even stronger advantage of one digit in the runtime. Hence, the overall energy consumption is dramatically reduced.

Remark 4.1

Interestingly, using Lagrange shape functions, only (11) worked correctly, using (10), the Newton–Raphson iteration diverges. Using B-Splines, (10) converges without any problems .Footnote 1 As we focus on the preferable formulation in (11), we do not discuss this further here. Using B-Splines or NURBS shape functions is undoubtedly possible, and we leave it to the reader to make use of the preferred framework. We refer to our preliminary work in [25] for details on the use of B-Splines in space-time.

Next, we investigate the condition number of the finite element formulation. For classical time-stepping methods using constrained dynamics (e.g., quarterions or director formulations), the condition number of the iteration matrix of the differential-algebraic equations (DAEs) can be estimated as of \(O(\Delta t^{-3})\), see Betsch [5], Appendix A. For the constrained finite element formulation at hand, we obtain the condition data shown in Fig. 3. Note that this is the condition number of the matrix for all time steps at once.

Fig. 3
figure 3

Condition number plotted over the element size h

As a second example, a more complex system is introduced using the same external momentum of \(\boldsymbol{m}\), however, applied during the time interval \(\mathcal{T}^{\operatorname{rot}} = [0.5,\,1.5]\) to initiate large rotations in the temporal domain \(\mathcal{T} = [0,\,4]\). Remark that this setup can be considered as a nonsmooth system with respect to the external contributions. Figure 1, right, shows the motions throughout the time interval. It is obvious that we consider large, three-dimensional rotations of the free-flying body. Note that a rectangular load, zero for \(t = [0,\,0.5]\) and \(t = [1.5,\,4]\) and the prescribed external momentum else, is challenging for finite elements that rely on a predefined smoothness; however, this is a typical setup for rigid body problems.

Figure 4 shows the results of the rectangular load. The implicit Euler yields linear convergence, the linear space-time a quadratic convergence, and the quadratic Lagrange shape functions a cubic convergence. Note that the results cannot be as perfect as for the first example. If the time step or the element size matches exactly with the nonsmooth changes in the external load, i.e., the load initiates or terminates exactly at the element or time step boundary, the total error is better than otherwise. For the overkill solution, time-stepping schemes cannot be used any more since the condition number for the index-3 problem raises dramatically with \(O(\Delta t^{-3})\), cf. Betsch [5], and thus the stop criterion of the Newton–Raphson iteration must be lowered dramatically. Thus, we make use of a blockwise calculation of the space-time system using quadratic Lagrange shape functions with 40,000 elements within each of the 40 blocks to obtain an element size of \(2.5\times 10^{-6}\) for the overkill solution. The final position of the center of mass after 4 seconds using an initial velocity of \([5,\,5,\,5]\) is obviously \([20,\,20,\,20]\), the directors are given at \(\boldsymbol{d}_{1} = [0.392,\,-0.271,\,0.879]\), \(\boldsymbol{d}_{2} = [-0.233,\,-0.954,\,-0.190]\), and \(\boldsymbol{d}_{3} = [0.890,\,-0.131,\,-0.437]\). Note that the midpoint rule shows a quadratic convergence only below a time step size \(10^{-4}\).

Fig. 4
figure 4

Convergence test for rectangular load. \(Err\) is plotted over the element size h

Finally, we plotted the convergence of the Lagrange multipliers in Fig. 5. As shown in [11] for the case of Mortar-type boundary conditions, the convergence of a \(p/p\) pairing is restricted.

Fig. 5
figure 5

Convergence test for rectangular load. \(Err\) of the Lagrange multiplier is plotted over the element size h

4.2 Double spherical pendulum

Here, we aim at several additional investigations following the example of a double spherical pendulum presented in Betsch [5], Sect. 4.2. First, we apply the space-time formulation on a two-body system, connected via two additional, spherical constraints using Lagrange multipliers. Second, we apply the space-time formulation block-wise, i.e., we make use of time blocks, calculated using a standard Newton–Raphson iteration for the space-time problem. Therefore, the position \(\bar{\boldsymbol{\varphi }}\) and directors \(\boldsymbol{d}_{i}\) as well as their momentum \(\bar{\boldsymbol{\pi }}\) and \(\boldsymbol{\pi }^{d}_{i}\) at the end of the time block are transferred to the next time block as initial conditions. This allows for arbitrary long temporal domains without jeopardizing the Newton–Raphson convergence as the solutions may be far away from the initial guess.

Again, we conduct two different tests. First, we apply the initial configuration as presented in Fig. 6 using two cubic geometries with length, width, and height of \([1,\,1,\, 1]\) and \([0.8,\,0.8,\, 0.8]\). The initial position of the center of mass is \(\bar{\boldsymbol{\varphi }}_{1} = 4/\sqrt{2}\times [1,\, 0,\, 1]\) for the first and \(\bar{\boldsymbol{\varphi }}_{1} = 8/\sqrt{2}\times [1,\, 0,\, 1]\) for the second body. The directors form a unity matrix for both bodies and the initial velocities are zero. A vector of gravitational acceleration is \(\boldsymbol{g}= [0,\, 0,\, -9.81]\in \mathbb{R}^{3}\), such that the corresponding energy function reads \(V(\bar{\boldsymbol{\varphi }}_{1},\bar{\boldsymbol{\varphi }}_{2}) = M_{ \varphi}^{1}\bar{\boldsymbol{\varphi }}_{1}\cdot \boldsymbol{g}+M_{ \varphi}^{2}\bar{\boldsymbol{\varphi }}_{2}\cdot \boldsymbol{g}\). The rigidity of the connecting rods gives rise to constraint functions written as

$$ \Phi _{1} = \frac{1}{2}\left (\|\bar{\boldsymbol{\varphi }}_{1}\|^{2} - L_{1}^{2}\right ),\quad \text{and}\quad \Phi _{2} = \frac{1}{2}\left ( \|(\bar{\boldsymbol{\varphi }}_{2}-\bar{\boldsymbol{\varphi }}_{1})\|^{2} - L_{2}^{2}\right ), $$
(27)

where \(L_{1}\) and \(L_{2}\) are the initial lengths. No external forces are applied so that the system is expected to behave like a two-dimensional double pendulum, equipped with two mass points. The initial velocity of both bodies is zero, the time interval is \(\mathcal{T} = [0,\,80]\), 40 blocks are used with 100 temporal elements each. Figure 7 displays the x-position of both bodies over time, whereas Fig. 8 displays the z-position. Note that for both bodies, the y-position remains below \(10^{-30}\), i.e., we obtain in fact a two-dimensional motion.

Fig. 6
figure 6

Double spherical pendulum with two rigid, cubic-shaped bodies

Fig. 7
figure 7

Double pendulum, x-position of both bodies over time

Fig. 8
figure 8

Double pendulum, z-position of both bodies over time

This changes dramatically if we constrain the second body not in the center of mass. Therefore, the reference point \(\bar{\varphi }_{2}\) is moved out of the center of mass at \([-0.1,\,0,\,0]\). Thus, we obtain for the mass \(M_{\bar{\varphi }_{2}} = 0.768\), which is obviously unchanged, for \(\boldsymbol{S}_{\bar{\varphi }_{2}} = [0.0768,\,0,\,0]\) and \(\boldsymbol{E}_{\bar{\varphi }_{2}} = \operatorname{diag}([0.0486,\,0.0410,\,0.0410])\). The trajectory is shown in Fig. 11; as expected the trajectory is not two-dimensional any more. For this specific example, we make use of 200 time blocks with 200 finite elements within each block. The coordinates at the end point are \(\bar{\varphi }_{1} = [0.0585,\,-0.0002,\,-3.9996]\) and \(\bar{\varphi }_{2} = [-2.3571,\,0.0002,\,-71{,}878]\).

A direct evaluation of the energy is not reasonable as this is a strong evaluation of an integral formulation. As we have shown in Hesch et al. [14], we can conduct integration by parts and compare the initial and the final values. Here, we conducted a different evaluation: For each element, we plotted the action integral in Fig. 9. As can be seen, we obtain for each block a constant function for the action integral; since the total energy is the derivative of the action integral with respect to time, the total energy is a conserved value over time. Note that we observe of an extremely small disturbance at the beginning of each block. We enforce the constraints integral in the temporal domain, i.e., they are not exactly fulfilled at a single point. Hence, if we transfer the end values of one block as initial conditions to the next block, the spherical constraints may be violated. Consequently, we do not observe this phenomenon using a single block and this instability vanishes with smaller time elements, as can be seen in Fig. 9.

Fig. 9
figure 9

Double spherical pendulum, action integral over time for different resolutions using \([500,\,1000,\,2000]\) elements per block. Note that the energy is conserved below the Newton–Raphson stop criterion within each block, as the action integral remains constant near machine precision in each block

In Fig. 10 all three components of the angular momentum, integrated over time for each element (see Appendix C) and added up over both bodies, are plotted. As expected, only one component is active, the other two are constants below the Newton–Raphson stop criterion over time, hence both components of the angular momentum are conserved.

Fig. 10
figure 10

Components of angular momentum, integrated over time using 1000 elements per block. Results are plotted over time. Note that the two components L1 and L3 overlap

Fig. 11
figure 11

Double spherical pendulum, trajectory of both bodies using an out of the center of mass connection of the second body. Blue: Trajectory of \(\bar{\boldsymbol{\varphi }}_{1}\). Red: Trajectory of \(\bar{\boldsymbol{\varphi }}_{2}\)

For the final test, we make use of the same initial configuration of the double pendulum; however, the initial velocity is predefined for both bodies using \(\nabla _{t}\bar{\boldsymbol{\varphi }} = [0,\,1,\,0]\) and \(\nabla _{t} \boldsymbol{d}_{3} = [0,\,1,\,0]\), the initial velocity for both other directors \(\boldsymbol{d}_{1}\) and \(\boldsymbol{d}_{2}\) is set to zero. Figure 12 shows the trajectory of both centers of mass \(\bar{\boldsymbol{\varphi }}_{1}\) and \(\bar{\boldsymbol{\varphi }}_{2}\). Additionally, we obtain large, three-dimensional rotations of both rigid bodies throughout the whole time interval. As shown in Fig. 13, the action integral remains constant in each block. Except for the boundaries of the blocks, the change of the action integral element by element is below the Newton tolerance.

Fig. 12
figure 12

Double spherical pendulum, trajectory of both bodies. Blue: Trajectory of \(\bar{\boldsymbol{\varphi }}_{1}\). Red: Trajectory of \(\bar{\boldsymbol{\varphi }}_{2}\)

Fig. 13
figure 13

Double spherical pendulum, action integral of the 3D motion plotted over time for different resolutions using \([500,\,1000,\,2000]\) elements per block

Finally, in Fig. 14 all three components of the angular momentum, again integrated separately for each element over time and added up over both bodies, are displayed for this three-dimensional motion. As can be seen, the \(L_{3}\) component remains constant, i.e., the third angular momentum component is zero.

Fig. 14
figure 14

Angular momentum, integrated over time using 1000 elements per block. Results are plotted over time

5 Conclusions

In this contribution, we have shown how to apply space-time concepts to rigid multibody dynamics. In particular, constrained rigid body formulations based on a director formulation are fitted within a space-time cylinder. As the spatial components of the kinetic energy like mass matrix as well as first and second moment of inertia can be evaluated independently from the temporal components as usual in rigid body dynamics, we do not have a wave-equation like structure. Hence, we do not have to stabilize the formulation as generally necessary for mechanical space-time problems, cf. the preliminary contribution of the authors in Schuss et al. [25]. Again, we could apply an interpretation of Livens’ theorem; this allows us to simplify the application of initial and/or end conditions dramatically.

Moreover, we could show in Appendix B that the initial and end boundary conditions arise in a most natural way within the variational formulation. This resolves in our own way a decade-old discussion on the formulation of the Lagrangian itself. In the numerical examples, we could not only demonstrate the applicability of the chosen approach and its superiority against classical time-stepping formulation. Following current trends in computer hardware, we can now make use of the rising core numbers in modern architectures for rigid body simulation. Sophisticated parallelization techniques, as developed for decades, can now be applied in a straightforward manner. As we have shown, this reduces the energy consumption of multibody simulation dramatically, which is a major goal in research today.