David González

Elías Cueto
email: ecueto@unizar.es

Francisco Chinesta

Real-time direct integration of reduced solid dynamics equations

Keywords: Solid dynamics, PGD, POD, model order reduction, real-time simulation

A new method is developed here for the real-time integration of the equations of solid dynamics based on the use of POD-PGD approaches and direct time integration. The method is based upon the formulation of solid dynamics equations as a parametric problem, depending on their initial conditions. A sort of black-box integrator is obtained that takes the resulting displacement field of the current time step as input and (via POD) provides the result for the subsequent time step at feedback rates on the order of 1 kHz. To avoid the so-called curse of dimensionality produced by the large amount of parameters in the formulation (one per degree of freedom of the full model) a combined POD-PGD strategy is implemented. Examples are included that show the promising results of this technique.

Introduction

Despite the well-known (and very often mentioned) growth in computer capabilities, some engineering problems still constitute today an intensive task. Structural dynamics is one of these problems, where the set of differential equations describing the system is usually discretized in space by finite elements (leading to an ODE in time) and then explicit or implicit methods in time are used to further discretize this ODE. This results in a system of very large number of degrees of freedom if fine enough meshes are employed in space, and (specially for explicit approaches) small-enough time steps are employed to avoid instability issues.

For this reason, structural dynamics has been one of the first fields of engineering practice where model order reduction has been employed. Classical modal analysis can be seen as a sort of model order reduction while, for instance, [START_REF] Nickell | Nonlinear dynamics by mode superposition[END_REF] seems to have been the first in employing eigenmodes of the structure with the purpose of reducing the number of degrees of freedom in non-linear applications. In essence, (projection-based) model order reduction consists in changing the original degrees of freedom of the model (i.e., finite element mesh nodal displacements) by an as much as possible reduced number of generalized variables that describe adequately the large-scale structures dominating the solution of the problem [START_REF] Meyer | Efficient model reduction in non-linear dynamics using the Karhunen-Loève expansion and dual-weighted-residual methods[END_REF]. Recently, a number of works have pursued this same idea in the field of solid dynamics, see, to name a few, [START_REF] Noor | Reduced basis technique for nonlinear analysis of structures[END_REF] [START_REF] Idelsohn | A reduction method for nonlinear structural dynamics analysis[END_REF] [START_REF] Krysl | Dimensional model reduction in non-linear finite element dynamics of solids and structures[END_REF] [25] [START_REF] Taylor | A reduced order finite element algorithm for surgical simulation[END_REF]. In these last works it is shown that model order reduction allows not only for a significant reduction of the number of spatial degrees of freedom, but also for larger time steps verifying stability constraints. However, in these references, an implicit simulation of 1.5 seconds of real time consumes 8540 seconds of computing time. Even if savings on the order of ten times have been reported for a judicious choice of the number of degrees of freedom, this is still far from real-time constraints. Some applications are specially demanding for feedback responses of structural dynamics simulations. For instance, it is well known that haptic (force) peripherals require feedback responses on the order of 500 Hz to 1 kHz. This is equivalent to solve non-linear solid dynamics each millisecond, or one thousand times per second. To the best of our knowledge, only linear elastic models undergoing large deformations (Saint Venant-Kirchhoff models) have been solved under such an astringent requirement by employing either GPUs [START_REF] Taylor | High-speed nonlinear finite element analysis for surgical simulation using graphics processing units[END_REF] or explicit algorithms combined with model reduction [START_REF] Taylor | A reduced order explicit dynamic finite element algorithm for surgical simulation[END_REF] [7] [START_REF] Barbič | Time-critical distributed contact for 6-DoF haptic rendering of adaptively sampled reduced deformable models[END_REF]. The extension of these techniques to account for the stability provided by implicit time integrators and general non-linear models seems to be difficult, even for reduced order models.

In this work a completely different (yet still based on some form of model order reduction) philosophy is employed. A computational vademecum [START_REF] Chinesta | PGD-Based Computational Vademecum for Efficient Design, Optimization and Control[END_REF] is developed that acts as a time integrator. This sort of black box is fed with the results of the previous time step, t, (in the form of the vector of nodal displacements of the model) and provides, under the severe feedback restrictions commented before, the solution at time t+ t, still providing the stability features of your preferred time integrator, implicit or not.

To achieve this kind of "universal" time integrator, the approach here presented first poses the problem in a parametric framework. This means that the vector of displacements at time t is considered as parameters of the problem (in addition to being the initial conditions of the subsequent time step), the displacements at time t + t being the unknowns. Of course, the number of degrees of freedom in a nowadays problem of industrial interest is huge, so the number of parameters could be simply unaffordable. To alleviate this problem, a reduced order basis, obtained via Proper Orthogonal Decomposition (POD) [START_REF] Karhunen | Uber lineare methoden in der wahrscheinlichkeitsrechnung[END_REF] [27] [START_REF] Park | The use of the Karhunen-Loève decomposition for the modeling of distributed parameter systems[END_REF] [42] [START_REF] Radermacher | Proper orthogonal decomposition-based model reduction for nonlinear biomechanical analysis[END_REF] is employed, thus guaranteeing the parametrization of the space of initial conditions with a minimum of degrees of freedom. Even with this efficient approach for the initial conditions, the number of parameters (initial conditions) usually reaches some tens.

The method is then structured in an off-line/on-line strategy. The parametric problem, depending on the initial conditions in the interval (t, t + t) is then solved off-line, and once for life, for any value of the parameters within a prescribed interval. This parametric problem can be solved by using your favorite time integrator, implicit or explicit, and the time step for these integrators need not to be the same t (can be smaller if, for some applications, the feedback response requirement is bigger than the stability constraints of the particular integrator used). Once this off-line phase of the method is computed (as already mentioned, only once for life), the on-line part can run under severe real-time constraints by simply feeding the "vademecum" with the results at the previous converged time step.

The only conceptual difficulty in this approach is to solve the parametric problem. Note that this problem considers every degree of freedom of the (reduced) model as a parameter. Thus, it is frequent to consider some tens of parameters for large (full) problems comprising thousands of finite element degrees of freedom. It is well-known that it is not possible to solve this kind of parametric problems by simply employing any mesh-based technique (finite elements, finite differences), since the number of degrees of freedom increases exponentially with the number of state space dimensions (equal to the number of spatial dimensions plus the number of parameters) [START_REF] Laughlin | The theory of everything[END_REF].

To alleviate this curse of dimensionality, the authors, among others, have developed a technique coined as Proper Generalized Decomposition (PGD) [5] [4] [10] [START_REF] Chinesta | A short review on model order reduction based on proper generalized decomposition[END_REF] [3] [START_REF] González | Recent advances on the use of separated representations[END_REF][15] that approximates the solution as a finite sum of separable functions. These functions are determined in the off-line part of the method by invoking a greedy algorithm (whose solution involves a non-linear problem, usually solved by simple fixed-point, alternate directions algorithm). The method, whose roots can be traced back to the space-time separation in the so-called radial loading technique within the LATIN method [START_REF] Ladeveze | Nonlinear Computational Structural Mechanics[END_REF] [23] [START_REF] Ladeveze | On the verification of model reduction methods based on the proper generalized decomposition[END_REF] [39] [START_REF] Ch | A rational strategy for the resolution of parametrized problems in the {PGD} framework[END_REF], allows for an efficient solution of problems defined in high dimensional parameter spaces. In the PGD framework, the complexity of the problem scales linearly with the number of state space dimensions (parameters, actually), instead of exponentially.

With all these ingredients, that will be detailed in the sequel, the method here proposed provides with a sort of meta-model or surface response (computational vademecum, in words of [START_REF] Chinesta | PGD-Based Computational Vademecum for Efficient Design, Optimization and Control[END_REF]) of the system subjected to any type of force and initial conditions (taking values in a given interval).

The outline of the paper is as follows. In Section 2 we describe the basics of the formulation of the problem of solid dynamics here addressed. In Section 3 the parametric form of the problem is introduced as a pre-requisite for the establishment of the off-line phase of the method. Section 4 develops the PGD-based formulation of the resulting multi-dimensional problem while Section 5 explains how to obtain an efficient parametrization of the space of initial conditions. Finally, Section 6 includes some numerical examples showing the potentiality of the proposed method. The paper is completed with some conclusions and comments on future works.

Solid dynamics

The strong form of the problem of solid dynamics is well-known, but it is reviewed here for completeness. We consider first the linear elastic case. Non linear strain measures will be analyzed in Section 6.2. For the former case the problem can be stated as: given f , g, h, u 0 and u0 , find

u : Ω×]0, T ] → R 3 such that ρü = ∇ • σ + f in Ω×]0, T ] (1a) u = g on Γ u ×]0, T ] (1b) σn = h on Γ t ×]0, T ] (1c) u(x, 0) = u 0 (x) x ∈ Ω (1d) u(x, 0) = u0 (x) x ∈ Ω. (1e)
The corresponding weak form is: given f , g, h, u 0 and u0 find u(t)

∈ S t = {u|u(x, t) = g(x, t), x ∈ Γ u , u ∈ H 1 (Ω)}, t ∈ [0, T ], such that for all w ∈ V{u|u(x, t) = 0, x ∈ Γ u , u ∈ H 1 (Ω)}, (w, ρü) + a(w, u) = (w, f ) + (w, h) Γ (2a) (w, ρu(0)) = (w, ρu 0 ) (2b) (w, ρ u(0)) = (w, ρ u0 ), (2c)
where we have used the so-called direct notation:

a(w, u) = Ω ∇ s w : C : ∇ s u dΩ, (w, f ) = Ω wf dΩ, (w, h) Γ = Γt whdΓ.
This constitutes the main ingredient for a finite element semi-discretization of the problem, leading to an ODE to be discretized in time by your favorite time integrator, implicit or explicit. The method here proposed proceeds by stating the problem (2) in a parametric form. In this case, the parameters will be the initial conditions.

Parametric, multidimensional, form of the problem

The main modification in the standard approach to the problem given by Eqs. ( 1)-( 2) is to reinterpret the parametric dependence of the problem on its initial conditions u 0 and u0 (assumed for the time being constant; the general case will be addressed later on) as being actually multi-dimensional:

u : Ω×]0, T ] × I × J → R 3 ,
where I = [u - 0 , u + 0 ] and J = [ u-0 , u+ 0 ] represent the considered intervals of variation of initial boundary conditions, u 0 and u0 . This makes it necessary to define a new (triply-) weak form where:

a(w, u) = I J Ω ∇ s w : C : ∇ s u dΩd u0 du 0 , (w, f ) = I J Ω wf dΩd u0 du 0 , (w, h) Γ = I J Γt wh dΓd u0 du 0 .
With this modification in mind, one can proceed as in standard textbooks of finite elements and semi-discretize, with the help of appropriate finite-dimensional approximations to S t and V, S h t and V h , respectively, the weak form so as to obtain the following problem: given f , g, h, u 0 and u0 find u h (t) = v h + g h ∈ S h t (note that g(x, t) = u(x, t) on Γ u ) such that for every

w h ∈ V h , (w h , ρv h ) + a(w h , v) = (w h , f ) + (w h , h) Γ -(w h , ρg h ) -a(w h , g h ), (5a)
(w h , ρv h (0)) = (w h , ρu 0 ) -(w h , ρg h (0)), (5b)
(w h , ρ vh (0)) = (w h , ρ u0 ) -(w h , ρ ġh (0)). (5c)
This results in an absolutely general solution (thus the name computational vademecum coined in [START_REF] Chinesta | PGD-Based Computational Vademecum for Efficient Design, Optimization and Control[END_REF]) to the problem (2) for any initial conditions. A similar philosophy, but for simpler problems was implemented in [START_REF] Gonzalez | Proper generalized decomposition based dynamic data driven inverse identification[END_REF]. It can also be viewed as a sort of response surface or meta-model obtained in an entirely computational way and in one execution of the algorithm. The method here developed computes off-line and once for life this general solution so as to be able to simply post-process during the on-line phase of the algorithm.

However, there remain some difficulties to be addressed, that will soon become apparent. The most important one is the high-dimensional nature of the problem [START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids[END_REF]. Once semi-discretized, Eqs. ( 5) involve u 0 and u0 as parameters. When these initial conditions evolve in the domain they could be expressed by means of their respective nodal values at the finite element mesh. This means that the number of actual parameters is six times the number of nodes in the mesh (three for u 0 and three more for u0 at each nodal location). This is simply unaffordable for any known method to date (several tens or hundreds of thousands parameters for industrial applications).

In order to divide and conquer, we address first discrete problems in Section 4 and then move to continuous problems discretized with the help of finite elements in Section 6.

PGD approximation to parametric solid dynamics

High-dimensional approximation

To efficiently solve solid dynamics Eqs. (2), one possibility lies in considering space-time approximations to the problem, as in [START_REF] Ladeveze | Nonlinear Computational Structural Mechanics[END_REF] [START_REF] Boucinha | Space-time proper generalized decompositions for the resolution of transient elastodynamic models[END_REF]. Here, we follow a more standard, although multidimensional, PGD approximation [10] [12] to the problem given by Eqs. [START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids[END_REF], by considering direct time integration † , as:

v h (x, t, u 0 , u0 ) = N i=1 F i (x) • G i (u 0 ) • H i ( u0 ) • d(t), (6)
where the nodal coefficients d(t) carry out all the time-dependency of the solution and the symbol "•" stands for the entry-wise Hadamard or Schur multiplication of vectors. Functions F , G and H will be expressed in terms of low (here, three-) dimensional finite element basis functions. The computation of these basis is one of the more salient features of PGD, since they are obtained on the fly (although in the off-line phase of the method), a priori, without the need to prior computer experiments, typical of POD-based techniques. They are computed by means of a greedy algorithm in which one sum is computed at a time, while one product is computed in a fixed point, alternated directions algorithm (here the most classic PGD algorithm is employed very satisfactorily, see [START_REF] Nouy | A priori model reduction through Proper Generalized Decomposition for solving time-dependent partial differential equations[END_REF] for alternative formulations). Thus, having an approximation to v h converged at iteration n, the (n + 1)-th term is obtained as

v n+1 (x, t, u 0 , u0 ) = n i=1 F i (x) • G i (u 0 ) • H i ( u0 ) + R(x) • S(u 0 ) • T ( u0 ) • d(t).
Remark 1. Note that here n and n + 1 refer to the number of terms in the PGD separated representation, not to the time step.

In this framework, an admissible variation of v h is computed as

w h = R * • S • T + R • S * • T + R • S • T * ,
where functional dependencies have been omitted for clarity.

If we inject the approximations to v h and w h into the weak form of the problem, Eqs. [START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids[END_REF], we arrive at a semi-discrete problem that can be approximated in time using your favourite time integrator.

Still a more general expression than that of Eq. ( 6) can be envisaged so as to include the dependence of the response of the system with respect to any load history. Thus, we consider an approximation of the type:

v h (x, t, u 0 , u0 , h 0 ) = N i=1 F i (x) • G i (u 0 ) • H i ( u0 ) • J i (h 0 ) • d(t),
where at each time instant t it is assumed that the modulus of the load h(t) (assumed constant through that time step and equal to its value at the beginning, h(0)), can take any value within a prescribed interval h 0 ∈ [h -, h + ]. This renders a completely general, multidimensional, expression for v h that includes its dependence on initial conditions as well as on the value of the applied load at each time interval. Knowing in advance the analytical expression of the load and particularizing it at each time instant allows to obtain the response of the structure for any load taking values within

[h -, h + ].
The very last detail in the implementation is to search for an approximation not for the whole time interval of the problem, ]0, T ], but for ]0, t]:

v : Ω×]0, t] × I × J × [h -, h + ] → R 3 ,
where t represents the necessary time to response prescribed by the particular envisaged application. For instance, for haptic feedback it has been already mentioned that a physical sensation of touch needs for some 500 Hz to 1 kHz feedback rate. This means that t = 0.001 seconds. This value t is not the necessary time step to achieve stability in the time integration chosen (that can be smaller if needed), although it can be coincident (and will be for all the examples shown hereafter).

In this way we obtain a completely general expression of v that can be particularized for any initial displacement and velocity fields. So to speak, a sort of black-box that, conveniently fed with the last time step results, provides with the displacement field at time t + t.

In what follows we show how the proposed technique behaves when applied to discrete problems. These problems posses analytical solutions that are employed to measure the influence of model reduction in the accuracy of the method.

Example 1: Frame structure

We consider here a simple example of validation with two degrees of freedom. It represents the frame structure depicted in Fig. 1. Slabs are assumed perfectly rigid, of masses m 2 = 2m 1 = 4000 kg, while the stiffness of the columns is such that 12EI Consider now that the structure is subjected to a step load h(t) = 10000 N, applied for t > 0. Note that the frame has been modeled with the aid of one single beam representing both columns, discretized into two finite elements. Only the degrees of freedom relative to the displacements u 1 and u 2 are free. In Fig. 2 the response of the structure is shown by comparing the PGD results with those of the analytical solution. A forward-Euler explicit time integration algorithm has been employed. As can be noticed, results are almost indistinguishable. Eight modes (i.e., n = 8 in Eq. (4.1)) have been considered in this case for almost a perfect agreement with the analytical solution, as noticed in Fig. 2. In general, some error estimation must be available to determine the number of modes necessary for a prescribed error tolerance. In the framework of PGD some previous works have addressed this issue [START_REF] Ammar | An error estimator for separated representations of highly multidimensional models[END_REF] [24] [START_REF] Moitinho De Almeida | A basis for bounding the errors of proper generalised decomposition solutions in solid mechanics[END_REF]. Functions F i , G i , H i and J i , i = 1, . . . , 8 for this problem are depicted in Fig. 3. The choice of the assumed limits of variation for each magnitude (intervals I, J ) is heuristic. If, during the simulation, any parameter escapes from its assumed interval of validity, a new simulation is started with bigger intervals. It has been noticed during this work, however, that uniformly-meshed intervals should be preferred in favor of a better convergence of the PGD algorithm.

L 3 = 10 5 N/m. h(t) u 1 u 2 m 1 m 2 EI, L EI, L

Example 2: Vibration of a clamped beam

We consider now the problem of a doubly-clamped beam subjected to a force in its middle point as sketched in Fig. 4. In IS units, the beam has an area of 2850 × 10 -6 , inertia moment 1.94 × 10 -5 , Young's modulus of 2.2 × 10 7 , length L = 1, and mass density ρ = 325.

As mentioned before, one only ejecution of the program provides with the response of the system under virtually any load (at least those taking values within a prescribed interval). Thus, it could be evaluated in real-time for a load h(t) = 20 sin(2πt), rendering the response in Fig. 5. In it, the sensibility of the explicit (forward-Euler) algorithm to different time steps, from 0.05 to 0.005 seconds, is shown. It must be emphasized that in this approach load is considered constant within every time increment, but variable in the whole simulation.

The proposed algorithm is by no means restricted to explicit time integration. Any time integrator can be used. For instance, we compare in Fig. 6 the response of the beam computed by an explicit, forward-Euler algorithm and an Energy-Momentum conserving integrator proposed in [START_REF] Bathe | Conserving energy and momentum in nonlinear dynamics: A simple implicit time integration scheme[END_REF]. The increase in accuracy can be noticed. Thus, energy conservation shows to be a question not only of the number , where they have been assumed to lie in the interval [-3, 3]m/s. The same applies, but for degree of freedom u2, with (d) and (e). Finally, (f) represents functions J i(h0). Even if the load is constant, it has been assumed to vary between 1000 and 15000N . All amplitudes are expressed in IS units.

of modes retained in the model (for which an appropriate error estimation should be available) but, notably, the time integrator used. At the limit, the just developed PGD approach converges to the full finite element solution is sufficient number of modes are considered, as proved in many previous references [START_REF] Chinesta | A short review on model order reduction based on proper generalized decomposition[END_REF] [10].

h(t) L u(x) x y

Parametrization of the space of initial conditions

So far discrete problems have been considered, with two and one degrees of freedom, respectively. Thus, initial conditions were composed by a vector of four components (two displacements plus two velocities) at most. When moving to large-scale continuous problems discretized by finite elements, initial conditions constitute a vector field in the domain Ω. After finite element discretization in space, both the unknown field v and initial conditions, u 0 and u0 , will be represented by the usual vector of nodal values for these variables. Thus, complex meshes could imply maybe several hundreds of thousands of degrees of freedom (if not millions). The number of parameters in the formulation presented before for discrete problems becomes thus too high. PGD has shown success for problems up to dimension 100 [START_REF] Chinesta | A short review on model order reduction based on proper generalized decomposition[END_REF] and maybe more, but it is not realistic to dream with efficient solutions for any number of dimensions. Therefore it seems reasonable to somehow optimize the parametrization of the space of initial conditions. To this end an optimal (in the sense of POD) basis is here constructed, based upon the results for complete problems similar to the one at hand.

A POD-PGD approach

It is well-known that assuming the displacement or velocity fields to be interpolated by piece-wise polynomials constitutes a reasonable assumption, that is on the basis of the success of finite element methods. However, it is not the best choice if we known in advance some information about the behavior of the response of the system. This is on the root of the development of model order reduction techniques [START_REF] Karhunen | Uber lineare methoden in der wahrscheinlichkeitsrechnung[END_REF] [27] [START_REF] Lorenz | Empirical Orthogonal Functions and Statistical Weather Prediction[END_REF]. Instead, it is preferred to run some tests on similar geometries and loading conditions and to obtain the best (in some norm) global, Ritz-like, shape functions that best approximate the solution for these particular problems, in the hope to be also a good choice for the original problem. This rationale is on the basis of Proper Orthogonal Decomposition (POD) methods.

Here a mixed POD-PGD approach is envisaged: based upon the previous solution of similar problems and the use of POD decompositions, a suitable parametrization of the space of initial conditions is obtained. Results for a particular problem will be then projected onto this basis and re-injected into the (off-line, multi-dimensional) obtained PGD solution to advance in time. We thus obtain a sort of response surface, meta-model or computational vademecum, to be evaluated on-line, that takes as input parameters the projection of the previous time step displacement and velocity field onto the reduced basis parametrization of initial conditions for the subsequent time step. A sketch of the proposed method is shown in Fig. 7.

!"#$%&'() *+%,-.+/-0'$./'1,2$ u=u(x,u n ,u n-1 )

u n 3$u n-1 ! u n+1 ,u n $ 4$!2/

Numerical Examples

To show how the proposed technique behaves for some engineering problems we consider two problems. The first represents the propagation of a wave on a solid block, while the second represents the palpation of a human liver during an endoscopic surgery procedure [START_REF] Niroomandi | Real-time deformable models of non-linear tissues by model reduction techniques[END_REF] [36] [START_REF] Niroomandi | Accounting for large deformations in real-time simulations of soft tissues based on reduced-order models[END_REF].

Wave propagation on a block

To show the performance of the proposed technique, we have analyzed an example of the propagation of a wave within a solid block of a continuous elastic material. The block is 5 × 5 × 10 m, with Young's modulus 2.66 × 10 5 kN/m 2 , Poisson's ratio 0.33, density ρ = 2.0 t/m 3 ( P-wave velocity c 1 = 443.9 m/s, data is obtained from [START_REF] Otto Von | Iterative coupling of fem and bem in 3d transient elastodynamics[END_REF]). The block is discretized with a very poor tetrahedral mesh involving 45 nodes, as shown in Fig. 8. In this example we do not care about the accuracy of the model itself, but rather to compare with state of the art numerical techniques. The block is subjected to an impulsive loading h(t) = 4.0 • H(t) kN/m 2 on the top face, where H(t) stands for the Heaviside function. At all other surfaces, the normal displacements are restraint to zero.

The problem is compared to a FE solution obtained by employing HHT integration [START_REF] Hilber | Improved numerical dissipation for time integration algorithms in structural dynamics[END_REF] with the following parameters: α = -0.05, β = 1 4 (1 -α) 2 and γ = 1 2 -α. This same solution is employed to extract the POD modes that will parameterize the space of initial conditions, see Fig. 9. The simulation employed a fixed time step of t = 0.00125 s, which corresponds more or less to the desired feedback rates for haptic devices and also provides good results for this particular problem. The displacement field at a typical time increment (in this case, the 201-th one) is shown in Fig. 10. Remark 2. Note that the discretization of the space of initial conditions must be done with an appropriate mesh. Too-coarse discretizations have led in our experience to meaningless results due to interpolation errors. For instance, in the wave propagation problem, results showing important, non-physical tractions where only compressive forces were applied have been found. It is important to estimate well the maximum amplitude that initial conditions could reach during the simulation (so as to avoid to perform extrapolations), but also a judicious choice of the meshes. Here, 100 nodes have been considered for the mesh (see Figs. 12 and13).

In Figs. [START_REF] Gonzalez | Proper generalized decomposition based dynamic data driven inverse identification[END_REF] and 16 a comparison is made between the reference solution and the solution obtained with a 3-degree-of-freedom solution and one employing 8 degrees of freedom, respectively, for the parametrization of the space of initial conditions. Both employed an explicit, forward Euler, time integration.

It is noteworthy to mention that these results, for a chosen time step of t = 0.00125 s. are obtained between 0.000992 and 0.001325 s. This enables the use of this technique for real-time haptic response purposes, where feedback rates between 0.5 -1.0 kHz are needed for a realistic touch sensation. These results are obtained by using Matlab running on a Windows-operated HP laptop. While the complete simulation took approximately 8 seconds, the same simulation in the Abaqus software takes 16 minutes and 35 seconds (thus making it 125 times faster than the mentioned commercial software).

If, on the contrary, we employ an energy-momentum conserving scheme such as that by Bathe [START_REF] Bathe | Conserving energy and momentum in nonlinear dynamics: A simple implicit time integration scheme[END_REF] results improve very much, and a great stability is conserved even for large integration times (10 seconds in the example depicted in Fig. 17).

Including non-linear strain measures

Extending the before presented technique to non-linear (hyper-)elastodynamics does not incorporate any conceptual difficulty, albeit some details deserve a comment. In [START_REF] Niroomandi | Real-time simulation of biological soft tissues: a pgd approach[END_REF] and [START_REF] Niroomandi | Model order reduction in hyperelasticity: a proper generalized decomposition approach[END_REF] the authors have presented and analyzed two different approaches for the efficient linealisation of hyperelastic constitutive models in the context of quasi-static PGD models. Notably, these works present an explicit and an implicit (respectively) approach to this problem that do not imply the need to re-construct the full tangent stiffness matrix of the complete model. This constitutes an active field of research in the framework of model order reduction methods, since non-linearities usually destroy all the gains obtained in reducing the model, if tangent matrices need to be assembled.

Here, an explicit approach, similar to that in [START_REF] Niroomandi | Real-time simulation of biological soft tissues: a pgd approach[END_REF] has been implemented. We consider again the problem of the wave propagation in Section 6.1 and, without modifying the POD basis employed to for the Green-Lagrange strain tensor, E, in the form

E = 1 2 (F T F -1) = γ l (u) + γ nl (u, u),
where F = ∇u + I is the gradient of deformation tensor and

γ l (u) = 1 2 (∇(u T ) + ∇(u)), (7a)
γ nl (u, u) = 1 2 ∇(u T )∇(u). (7b)
Hyperelastic materials are base on the assumption of a particular strain-energy function, Ψ. Then the second Piola-Kirchhoff stress tensor S can thus be obtained by

S = ∂Ψ ∂E ,
that is a symmetric tensor and is related to the first Piola-Kirchhoff stress tensor, P , by P = F S. The Kirchhoff-Saint Venant model is characterized by the energy function given by

Ψ = λ 2 (tr(E)) 2 + µE : E,
where λ and µ are Lame's constants. The second Piola-Kirchhoff stress tensor can be obtained by

S = ∂Ψ(E) ∂E = C : E,
in which C is the fourth-order constitutive (elastic) tensor.

Although it is very well known that Kirchhoff-Saint Venant model shows instabilities in compression, it provides with a very neat formulation to compare with the linear, small strains case. In Fig. 18 a comparison is made between the already obtained (via HHT time integration scheme) reference, linear solution and the one considering non-linear strain measures (the so-called Kirchhoff-Saint Venant already described). As expected, non-linear strains, when considering all the remaining parameters of the problem unchanged, only provides with some bigger displacements in the solution, but an overall similar behavior of the block. Considering hyperelastic models does not imply any further delay in the response feedback given by the method here presented, since the number of on-line operations is exactly the same as for the case of linear elastic models. Of course, the off-line part of the method is somewhat more involved, but it is solved only once and stored, as repeatedly mentioned. Note that the inclusion of non-linear strain measures induces higher displacements, as expected, but an overall similar response in the problem. Solution obtained with one only mode in the POD basis of the space of initial conditions.

Palpation of a liver

As mentioned before, computational surgery is one of the fields where real-time requirements are more astringent [13] [1]. Surgery planning provides results corresponding to large periods after the operation in a reasonable amount of time so as to allow the surgeon to take decisions in the operation room. But, notably, minimally invasive surgery simulation and training systems, equipped with haptic (force feedback) devices require us to provide responses at kHz rates. The liver is the biggest gland in the human body, after the skin. As in [START_REF] Niroomandi | Model order reduction in hyperelasticity: a proper generalized decomposition approach[END_REF], liver geometry has been obtained from the SOFA project [START_REF] Allard | SOFA an Open Source Framework for Medical Simulation[END_REF] and post-processed in order to obtain a mesh composed by 2853 nodes and 10519 tetrahedra, whose geometry is shown in Fig. 19. The anterior surface of the liver is considered free, while the posterior one was assumed to be supported over different organs (it is connected to the diaphragm by the coronary ligament, for instance). The inferior vena cava travels along the posterior surface, and the liver is frequently assumed clamped a that location. Although the assumed boundary conditions are not strictly correct from a physiological point of view, our main interest is to show that the model can be solved under real-time constraints with reasonable accuracy. Here, the human liver is considered, for simplicity, as a linear elastic material with E = 0.17 MPa and ν = 0.48 [START_REF] Delingette | Soft tissue modeling for surgery simulation[END_REF].

Simulating a palpation, a ramp load of 5 N is applied at a particular point of the liver surface during a period of 0.25 seconds, and then released during other 0.25 seconds. The liver is then left vibrating free. Even if the liver tissue is well-known to posses some kind of viscoelastic properties, these have been neglected. The purpose of this example is not to obtain an extremely realistic simulation from a physiological point of view, but to show the performance of the technique. In particular, the influence of the number of modes chosen to parametrize the space of initial conditions on the long-term behaviour of the solution. To this end, a reference solution has been computed by employing an HHT time integrator and standard finite elements. POD modes have been extracted from this reference solution to construct the basis for the combined PGD-POD integrator. Other approaches, such as non-linear manifold learning [START_REF] Millán | Nonlinear manifold learning for model reduction in finite elastodynamics[END_REF] could also be employed and most likely would improve these results, but for the time being POD is considered as a suitable method to construct the basis.

It can be noticed that, in the idealized situation of absence of any type of damping, after the load release, the liver continues vibrating indefinitely. A PGD-POD solution has been computed by employing 1, 3 and 7 modes in the basis of the space of initial conditions. It can be noticed, from Fig. 20 that increasing the number of modes, as expected, provides converging results towards the reference solution. In fact, it can be noticed how the reduction in the number of degrees of freedom, very much like in classical model order reduction, eliminates high frequencies from the solution (those with the lowest content in energy) and therefore oscillations around the reference solution can be observed. The richer the basis is, the smaller the amplitude of these oscillations and the smaller their period.

Conclusions

In this paper a combined PGD-POD technique has been presented for the solution of solid dynamics problems. By formulating the problem as a parametric one, where initial conditions are the parameters, the proposed method is able to generate a sort of "black-box" integrator that takes displacement values at previous time steps as inputs and gives the solution for the current one. Within this black box any time integrator can be used. The technique here presented is aimed at solving the parametric problem, not at proposing new integrators. This appealing approach presents, however, some difficulties. One is the enormous amount of degrees of freedom generated if every nodal displacement representing initial conditions is taken as a parameter of the formulation. To this end, an optimal bests, from the POD point of view, is used by utilizing results of similar problems computed off-line. But even then a parametric problem should be solved that may involve tens of parameters. To overcome the curse of dimensionality provoked by this fact, a PGD approach has been employed. Thus, by computing off-line a general solution (very much like the spirit of a response surface or meta-model, but computed without the need for prior computer experiments) a kind of computational vademecum is obtained that can be emptied on-line at extremely high feedback rates, on the order of one kHz.

Previous works of the authors have shown how to include within this general formulation parameters such as load position, value and orientation, and how to efficiently solve non-linearities arising from material constitutive equations [35] [37]. This work on the dynamics of the problem somehow closes one fundamental aspect: that related with the interactivity of the simulation with the user, thus configuring a true dynamic data driven application system (DDDAS). Some fundamental aspects remain open, however. For instance, within the computational surgery applications, surgery cutting is a fundamental aspect to achieve realism in the simulation. Within this framework, an algorithm should be developed able to predict the response of the organ to a surgical cut placed at any point of the model and with arbitrary length. This constitutes one of our current efforts of research.

Figure 1 .

 1 Figure 1. Frame structure considered in Section 4.2.

Figure 2 .

 2 Figure 2. Response of the frame structure. PGD results vs. analytical solution.

Figure 3 .

 3 Figure 3. Monodimensional functions used to approximate the response of the frame in Fig. 1 subjected to a step load. (a) F i(x). Note that the frame was modeled with the help of one single bar, discretized into two finite elements, and hence the nine degrees of freedom in the abscisa. Only the two non-clamped horizontal displacements are free (degrees of freedom 8 and 9 in the figure). (b) and (c) represent, respectively, the dependence of the solution on u n 1 and u n-1 1

Figure 4 .Figure 5 .

 45 Figure 4. Geometry of the doubly-clamped beam problem.

Figure 6 .

 6 Figure 6. Comparison of the results of the beam in Section. 4.3. Explicit vs. Energy-Momentum conserving integration proposed in [8].

Figure 7 .

 7 Figure 7. Sketch of the proposed algorithm. A POD basis is employed for the parametrization of the space of initial conditions of the problem (for forward Euler schemes, for instance, this means in practice that the solution deepens on the displacement at time steps n and n -1, respectively).

Figure 8 .Figure 9 .Figure 10 .

 8910 Figure 8. FE model for the problem of D-wave propagation. Points A and B are used to represent the solution in Figs. 15 and 16.

Figure 11 .

 11 Figure 11. Eight most important spatial modes F i(x) for the PGD solution of the wave propagation problem. Contour levels represent the value of the y component of the modes.

Figure 12 .

 12 Figure 12. Eight most important PGD modes Gi(u n ) for the wave propagation problem. Note that the function is three-dimensional, so that its components are depicted in (a), (b) and (c), respectively.

Figure 13 .Figure 14 .

 1314 Figure 13. Eight most important PGD modes Hi(u n-1 ) for the wave propagation problem. Note that the function is three-dimensional, so that its components are depicted in (a), (b) and (c), respectively.

Figure 15 .

 15 Figure 15. Wave propagation problem. Solution obtained with a three-degree-of-freedom parametrization of the space of initial conditions.

Figure 16 .

 16 Figure 16. Wave propagation problem. Solution obtained with an eight-degree-of-freedom parametrization of the space of initial conditions and an explicit integration in time.

Figure 17 .

 17 Figure 17. Wave propagation problem. Solution obtained with an eight-degree-of-freedom parametrization of the space of initial conditions and employing the energy momentum conserving scheme by Bathe[START_REF] Bathe | Conserving energy and momentum in nonlinear dynamics: A simple implicit time integration scheme[END_REF].

Figure 18 .

 18 Figure 18. Wave propagation problem. Comparison of the (linear elastic) HHT reference solution with the Kirchhoff-Saint Venant one. Note that the inclusion of non-linear strain measures induces higher displacements, as expected, but an overall similar response in the problem. Solution obtained with one only mode in the POD basis of the space of initial conditions.

Figure 19 .

 19 Figure 19. Geometry of the finite element model for the liver.

Figure 20 .

 20 Figure 20. Response of the liver for a peak load. Reference (FE) results, and PGD results with basis composed by 1, 3 and 7 modes.

† The use of direct time integration is mandatory if the objective of real-time simulation is to provide the analyst with an interactive tool. This is the case for many haptic applications, where space-time formulations as the one in[START_REF] Boucinha | Space-time proper generalized decompositions for the resolution of transient elastodynamic models[END_REF] provide with the solution of the entire time interval considered in the simulation, and therefore need to know a priori the entire loading history.

Contract/grant sponsor: Spanish Ministry of Economy and Innovation; contract/grant number: CICYT-DPI2011-27778-C02-


ANAEROBIC DIGESTION OF VINASSE AND CFD MODELLING APPROACH

CONTEXT AND OBJECTIVES OF THE WORK

• CHARACTERISATION of the substrates with the PHYSICO-CHEMICAL and RHEOLOGICAL properties ; • MODELLING and OPTIMISATION of WET and DRY anaerobic digestion ;

• CFD MODELLING of the fluids flow with the LARGE-EDDY SIMULATIONS (LES) approach with Boussinesq approximation ;

MATERIALS AND METHODS

RESULTS, DISCUSSIONS AND EXPECTED RESULTS

CONCLUSION AND PERSPECTIVES REFERENCES

Hélène Caillet 1 (helene.Caillet@univ-reunion.fr), Alain Bastide 1 (alain.bastide@univ-reunion.fr), Graziella Chuppa 2 (graziella.tostain@qualitropic.fr), , Thomas Petit 1,2 (Thomas.petit@univ-reunion.fr), Laetitia Adelard 1 (laetitia.adelard@univ-reunion.fr).

1 PIMENT Laboratory, Saint-Denis, Reunion Island, France 2 LCSNSA Laboratory, Sainte-Denis, Reunion Island, France  EXPERIMENTATIONS:

Main parameters Model validation Main outcomes

• Coupling of the two models ;

• Considerations of the rheological properties of the substrates ;

• Mechanical agitation optimisation.

 CFD MODELLING:

 EXPERIMENTATION RESULTS:

0 i i u x   
• Mass and momentum conservation equations:

      i j ij ij i z j i j j u u u P g e t x x x x                      2 . T i i j j D u t x x x            
• Heat transfer equation:

2 ij t ij S    
• SGS stress tensor:

We use the same geometry of Wu (2012). The internal cylinder is the sliding zone which contains the impeller.

• Dynamic Lagrangian turbulence model.

Figure

  Figure 1: Pilot (10L)

•

  Development of 3D SINGLE-PHASE model considering MECHANICAL AGITATION ;

m -3 Dynamic viscosity Pa.s

	Vinasse 1044.69	0.001 009 7
	Sludge 999.66	0.065

• Carry out CO-DIGESTION experimentations ; • Carry out PILOT tests at different TS rates ; • Develop MULTIPHASE model to predict the SPREAD OF MICROORGANISMS ; • COUPLING of the CFD model to the ADM1 model ; • OPTIMISATION of the MECHANICAL STIRRING of the process depending on the PHYSICAL PROPERTIES of the digestion medium. Fluid Density kg.

Table 1 :

 1 Physical parameters of Newtonian fluids

Table 2 :

 2 Physical parameters of sewage sludges in pseudo-plastic model

	TS %	T °C	k Pa.s n	n	Density kg.m -3	η min Pa.s	η max Pa.s	𝛄 s -1	J	C p

.kg -1 .°C -1

										λ
										W.m -1 .°C -1
	2.5	35	0.042	0.710	1000.36	0.006 0.008	4186.78	4186.78	0.6171
	5.4	35	0.103	0.639	1000.72	0.009 0.025	50-702	4185.94	0.6169
	7.5	35	0.525 0.533	1001.00	0.03	0.17	11-399	4185.33	0.6168
	9.1	35	1.052	0.467	1001.31	0.07	0.29	11-156	4184.87	0.6167
	12.1	35	5.885 0.367	1001.73	0.25	2.93	3-149	4184.00	0.6165
	15	17-24	31.3	0.30	973.0	3.35	44.24 0.61-24.37	4183.17	0.6163
	20	17-24	35.4	0.29	1091.8	3.69	48.6	0.64-24.14	4181.72	0.6159

Figure 3: Cumulative methane production average as a function of time of vinasse in BMP (left) and pilot (right) tests

Figure 2: Digester geometry (mm) B. Wu, "Large eddy simulation of mechanical mixing in anaerobic digesters," Biotechnol. Bioeng., vol. 109, no. 3, pp. 804-812, Mar. 2012. C. Holliger et al., "Towards a standardization of biomethane potential tests," Water Sci. Technol., vol. 74, no. 11, pp. 2515-2522, Dec. 2016. H. Landry, C. Laguë, and M. Roberge, "Physical and rheological properties of manure products," Appl. Eng. Agric., vol. 20, no. 3, p. 277, 2004. V. A. Vavilin and I. Angelidaki, "Anaerobic degradation of solid material: Importance of initiation centers for methanogenesis, mixing intensity, and 2D distributed model," Biotechnol. Bioeng., vol. 89, no. 1, pp. 113-122, Jan. 2005.

Table 3 :

 3 Physicochemical characteristics

Pilot test:

Temperature (T), consistency coefficient (k), power-law index (n), viscosity (η), shear rate ( 𝛾), specific heat (C P ), thermal conductivity (λ).

Velocity (u), fluid density (ρ), pressure (P), gravity (g), thermal expansion coefficient (β), temperature (θ), thermal diffusion coefficient (D T ), turbulent viscosity (μ t ), strain rate ( 𝑆 𝑖𝑗 ).

•

The vinasse is Newtonian ;

The sludge is non-Newtonian and modelled with the Power-Law:

The physical properties of liquid manure are used for the physical characterisation of sludge ; • Studies have shown that the viscosity and shear stress increase exponentially as total solid (TS) increases.

•

Hack Lange tests are used for the chemical oxygen demand (COD), the total carbon organic (TOC), the ammonium, the alkalinity and the Kjeldahl azote ;