J Girardot

• E Prulière

Elastic calibration of a discrete domain using a proper generalized decomposition

Keywords: Calibration, Lattice, Beam particle model, Discrete element method, Proper generalized decomposition

Current discrete/lattice methods can simulate a continuous mechanical behavior thanks to a network of bonds. The main drawback of these approaches is the need of a calibration process to link the emerging behavior of the structure and the parameters of the local mechanical bond. It is proposed in this work to use a fast and recent reduction model technique to build once and for all an exhaustive data chart and thus to avoid the calibration process. The proper generalized decomposition technique was used to build a parametric analysis in the case of a lattice beam structure. The results were in the range of the current calibration values found in the literature and extend it by giving a global calibration curve. They also allowed to discuss about the influence, in this specific case of lattice-beam structure, of the density of beams in terms of number of discrete elements and connectivity.

Introduction

The discrete approach was originally developed for granular media [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF]. It consists in evaluating contact forces on a multitude of discrete and generally spherical elements. The monitoring of the position of these elements over time is then carried out by integrating the accelerations, in the case of the smoothed contact dynamic [START_REF] Jean | Non-smooth contact dynamics approach of cohesive materials[END_REF]. This approach has been extended to continuous and homogeneous media by several authors, [START_REF] Delaplace | Discrete modeling of cracking of brittle materials in large relative motion and localization problem[END_REF][START_REF] Jebahi | Discrete element method to model 3D continuous materials[END_REF][START_REF] Pompe | Statistical models for the fracture of disordered media[END_REF] among others, by means of cohesive links. In these different works, these links are positioned thanks to their nodal extremities. The position of these nodes is obtained through different algorithms that can be grouped in two categories :

-algorithms based on a problem of granular stacking of discrete elements [START_REF] André | Discrete element method to simulate continuous material by using the cohesive beam model[END_REF];

B J. Girardot jeremie.girardot@ensam.eu E. Prulière etienne.pruliere@ensam.eu -algorithms based on a random spatial dispersion of a regular grid of nodes [START_REF] Vassaux | Beam-particle approach to model cracking and energy dissipation in concrete: identification strategy and validation[END_REF].

The final domain is therefore a so-called discrete medium, but with a continuous and homogeneous emerging behavior. The advantage of this type of representation is its ability to easily represent (in a numerical sense) fracture and/or contact mechanisms.

A calibration is nevertheless necessary to link the parameters of the cohesive bond to those of the desired homogeneous material [START_REF] André | The GranOO workbench, a new tool for developing discrete element simulations, and its application to tribological problems[END_REF][START_REF] Delaplace | Discrete 3D model as complimentary numerical testing for anisotropic damage[END_REF]. With an elastic bond behavior based on a Bernoulli beam kinematic, the calibration process links the Young's modulus, the Poisson's ratio and the radius of the beams (assumed with cylindrical shape) to the Young's modulus E and the Poisson's ratio ν of the material to simulate. It was indeed observed that the micro-elastic parameters needed for the beams are very different from the real ones at the macro-scale. For example, to represent a homogeneous elastic material with the following parameters E M = 30 GPa and ν M = 0.3, the micro-parameters are E µ = 200 GPa, ν µ = 0.4 and a radius ratio of 0.4 (for a coordination number of 6.2 using the same packing technique used in [START_REF] André | The GranOO workbench, a new tool for developing discrete element simulations, and its application to tribological problems[END_REF]).

At present, this calibration step is carried out at each new simulation work. Thanks to numerous parametric simulations, the recent work in [START_REF] Nguyen | Analytic laws for direct calibration of discrete element modeling of brittle elastic media using cohesive beam model[END_REF] gives some data charts to facilitate the calibration process but is limited by the range of the chosen parameters, and the same parametric study needs to be done for out-of-range input parameters or behavior (currently, only the elastic brittle behavior is taken into account in the literature dealing with calibration processes with also some attempts for orthotropic materials [START_REF] Curti | How to model orthotropic materials by the discrete element method (DEM): random sphere packing or regular cubic arrangement?[END_REF]).

So for a fast and global calibration, the point is to get the macroscopic response for all possible microscopic parameters. But this requires the resolution of a multi-parametric problem in high dimension.

There are several strategies to build parametric response of a system, mainly relying on meta-modeling (response surfaces, manifold learning, …). This kind of methods has in common the need of several problem solutions for many sets of parameters. These solutions are used to interpolate or to approximate the multi-dimensional solution. By principle, meta models do not rely on physical models, all the physics being embedded in the entry data.

Another approach for solving multi-dimensional problems directly from the physical equations has been developed during the last years. This approach uses the proper generalized decomposition (PGD) which is an efficient tool to solve a problem for all parameters in only one step.

Originally, the PGD was devised as an a priori model reduction method. The drawback of standards a posteriori model reduction techniques, such as the proper orthogonal decomposition (POD), is that they rely on some reduced basis that are difficult to obtain. The reduced basis is generally computed by costly finite element simulations [START_REF] Fahl | Computation of pod basis functions for fluid flows with lanczos methods[END_REF]. Two approaches have been developed to get the reduced solution of a problem without assuming the reduced basis (a priori model reduction):

-The a priori hyper-reduction method (APHR), which consists in enriching the reduced basis on the fly when the norm of the residual becomes higher than a tolerance value [20] -In the LATIN method, where the spatial reduced basis and the time evolution are computed at the same time [START_REF] Ladeveze | On a multiscale computational strategy with time and space homogenization for structural mechanics[END_REF]. The PGD approach is similar to the strategy used in the LATIN method.

Both of these methods were originally designed to solve space/time problems. After that, the PGD has been developed to solve the Fokker-Plank equation encountered in the statistical description of complex fluids based on the Kinetic theory with a drastic reduction of the computational cost [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].

Then, it has been used in a wide range of applications like, for instance, the mechanical behavior of thin structures [START_REF] Vidal | Explicit solutions for the modeling of laminated composite plates with arbitrary stacking sequences[END_REF]or to simulate an electromagnetic coupling [START_REF] Henneron | Proper generalized decomposition method applied to solve 3-D magnetoquasi-static field problems coupling with external electric circuits[END_REF]. In particular, the PGD has proven very efficient to solve parametric equations [START_REF] Nouy | A generalized spectral decomposition technique to solve a class of linear stochastic partial differential equations[END_REF]. In [START_REF] Pruliere | On the deterministic solution of multidimensional parametric models using the proper generalized decomposition[END_REF], an original use of the PGD is proposed to perform a parametric analysis of a mechanical problem. The PGD opens the possibility to build in only one simulation, the solution of a problem for a wide number of parameter sets. This hyper-solution can be seen as a computational vademecum and can be used for real-time applications and fast design procedures [START_REF] Chinesta | PGD-based computational vademecum for efficient design, optimization and control[END_REF]. Based on these former works, we propose to develop a model reduction approach to build in very short computation times a complete parametric analysis of a discrete domain. This study is a proof of concept and will be limited to:

-cohesive links of the elastic Bernoulli beam type -a nodal positioning strategy of granular type based on the work of André et al. [START_REF] André | Discrete element method to simulate continuous material by using the cohesive beam model[END_REF].

After detailing the parametric PGD implementation applied to a network of beams, the calibration results will then be compared to a literature reference. Finally, the influence on the calibration curves of the cohesive bond density within the discrete domain will be analyzed.

Domain generation

Granular packing

In this work, the generation of the domain is based on the work of André et al. [START_REF] André | The GranOO workbench, a new tool for developing discrete element simulations, and its application to tribological problems[END_REF][START_REF] André | Discrete element method to simulate continuous material by using the cohesive beam model[END_REF]. A granular packing is achieved in a three-dimensional cylindrical geometry in order to obtain the position of the discrete elements respecting classic granular filling properties such as the coordination number which is close to 6.2 in this case. These granular arrangement properties ensure a homogeneous positioning of the elements in a 3D volume. Then, a cohesive link is inserted between each particles in contact (i.e., when an interpenetration is recorded). Figure 1 shows the final result of this operation to create the positioning of discrete elements (whose radius is randomly drawn according to a Gaussian distribution) with the resulting network of cohesive links.

Beam model

A simple linear homogeneous beam based on the Navier-Bernoulli model is considered. The beam local coordinates (x, y, z) are defined such as the beam axis is set along the x axis. The displacement and rotation of the beam section in local coordinates are noted, respectively, u, v, w, θ x , θ y , and θ z . The displacement vector u is defined as: The beam is assumed to be a perfect cylinder whose length and diameter are noted L and D. The quadratic momentum following y and z is:

u = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ u v w θ x θ y θ z ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ . (1)
I yy = I zz = I = π D 4 32 . (2)
The quadratic polar momentum is then J = I /2, and the cross section area is S = π 4 D 2 . The material constitutive equation is linear elastic with an elastic coefficient E and a Poisson ratio ν. The second Lame coefficient is then:

µ = E 2(1 + ν) . (3)
The weak form of the beam equilibrium equation without inertia is given by:

L 0 ES du * dx du dx + µI 2 dθ * x dx dθ x dx +EI dθ * y dx dθ y dx + EI dθ * z dx dθ z dx dx = W ext (u * ), (4)
where W ext (u * ) is the virtual work associated with the external loads as a function of the virtual displacements u * , θ * x , θ * y and θ * z . Neglecting the shear deformation, the rotational degrees of freedom are related to the displacements one by θ y = dw dx and θ z = dv dx . The parametric PGD formulation will be described in the next section.

Tensile test

A monotonic pure tensile test is applied on the cylindrical discrete sample. Displacements on x-direction are applied on both extremities as 0 on the left side and 0.1 on the right to reach 10% of global deformation. We can notice then, since the model is linear, no large strain effect is considered and the global enforced strain has no effect on the calibration. This simple uniaxial tension test is generally used to identify the apparent Young's modulus and Poisson's ratio of the whole sample. To do so, mean elongation and radius evolutions are post-processed at the end of each simulations.

The sample used for the first step of this numerical work is made of a network of 33103 beams. This number comes from an approximate number of discrete elements (meaning undeformable spheres) of 10,000 which was defined as a converged density mesh [START_REF] André | Discrete element method to simulate continuous material by using the cohesive beam model[END_REF].

Figure 2 shows the displacement magnitude on each beam of the network on the sample. A typical emerging displacement field is observed for a homogeneous cylindrical sample loaded in uniaxial tension.

Proper generalized decomposition for a parametric Bernoulli beam

Parametric proper generalized decomposition approximation

We are looking for the solution of Eq. ( 4 To get the solution of the parametric problem, several strategies can be carried out. As mentioned in the introduction, we choose to use the parametric PGD defined in [START_REF] Pruliere | On the deterministic solution of multidimensional parametric models using the proper generalized decomposition[END_REF]. The idea is to take the parameters as additional coordinates. We are then looking for an approximation of the global solution using this simple separated expression :

u(x, E,ν,D) ≈ n i=0 F i (x) × G i (E) × H i (ν) × K i (D), (5)
where G i , H i and K i are scalar functions and F i are vector functions for all integer i ∈[0, n]. F i is defined by: where each component of F i is a scalar function defined on x.

F i = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ F u i F v i F w i F θ x i F θ y i F θ z i ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ . ( 6
)

Weak form of the global problem

We want to find the solution for every values of parameters in only one simulation. Since the parameters are treated as additional coordinates, the problem is formulated in 6 dimensions: 3 dimensions for space (x, y and z) and 3 dimensions for parameters (E, ν and D). A beam in itself is modeled by a 1D domain, but the DEM mesh is represented by a dense 3D lattice structure. The weak formulation ( 4) is then modified to account for the new parametric coordinates:

E max E min ν max ν min D max D min L 0 → E π D
dx → dxdDdνdE = W ext (u * ). (7)
The well-known curse of dimensionality makes the solving of the 6-dimensional problem very costly because of the explosion of the number of degrees of freedom (10 12 nodes if we use 100 nodes per dimension). The use of the PGD approximation gives a solution to this very high number of degrees of freedom. As the displacement is assumed on a separated form, a kinematically admissible virtual displacement will also be approximated using a separated form:

u * (x, E,ν,D) = F * (x) × G * (E) × H * (ν) × K * (D). (8)
Introducing equations ( 5) and [START_REF] Delaplace | Discrete 3D model as complimentary numerical testing for anisotropic damage[END_REF] in the weak form equation [START_REF] Curti | How to model orthotropic materials by the discrete element method (DEM): random sphere packing or regular cubic arrangement?[END_REF] and using the definition of µ in (3)gives: [START_REF] Delaplace | Discrete modeling of cracking of brittle materials in large relative motion and localization problem[END_REF] with the function Ψ defined by:

n i=0 Ψ(F * , G * , H * , K * , F i , G i , H i , K i ) = W ext (u * ),
Ψ(F * , G * , H * , K * , F i , G i , H i , K i ) = L 0 dF u * dx dF u i dx dx E max E min EG * G i dE ν max ν min H * H i dν E max E min π D 2 4 K * K i dD + L 0 dF θ x * dx dF θ x i dx dx E max E min EG * G i dE ν max ν min 1 2(1 + ν) H * H i dν E max E min π D 4 64 K * K i dD + L 0 dF θ y * dx dF θ y i dx dx E max E min EG * G i dE ν max ν min H * H i dν E max E min π D 4 32 K * K i dD + L 0 dF θ z * dx dF θ z i dx dx E max E min EG * G i dE ν max ν min H * H i dν E max E min π D 4 32 K * K i dD . ( 10
)

Iterative solving algorithm

Several algorithms exist to find the PGD solution of the problem defined above. The simplest algorithm based on a greedy approach combined with an alternate direction fixed point method has brought very satisfactory results in terms of precision and computational cost [START_REF] Nouy | A generalized spectral decomposition technique to solve a class of linear stochastic partial differential equations[END_REF].

We assume that the solution Eq. ( 5) is known until the iteration n and we want to enrich the solution with a new term n + 1. Of course, nothing is known at the first iteration, and the solution is built iteratively.

The displacement is then:

u(x, E,ν,D) = n i=0 F i (x) × G i (E) × H i (ν) × K i (D) +F n+1 (x) × G n+1 (E) × H n+1 (ν) × K n+1 (D) (11)
And the weak formulation can be rewritten as:

Ψ(F * , G * , H * , K * , F n+1 , G n+1 , H n+1 , K n+1 ) = W ext (u * ) - n i=0 Ψ(F * , G * , H * , K * , F i , G i , H i , K i ) (12)
We start from random values of the unknown F n+1 , G n+1 , H n+1 and K n+1 that verify the boundary conditions, and we update alternatively F n+1 , G n+1 , H n+1 and K n+1 knowing the others.

When updating F n+1 , the other functions G n+1 , H n+1 and K n+1 are assumed known and the virtual field can be then expressed by:

u * (x, E,ν,D) = F * (x) × G n+1 (E) × H n+1 (ν) × K n+1 (D) (13)
The weak form becomes:

Ψ(F * , G n+1 , H n+1 , K n+1 , F n+1 , G n+1 , H n+1 , K n+1 ) = W ext (u * ) - n i=0 Ψ(F * , G i , H i , K i , F i , G i , H i , K i ) (14)
All the integrals in Ψ [see Eq. [START_REF] Fahl | Computation of pod basis functions for fluid flows with lanczos methods[END_REF]] can be computed numerically excepted the ones related to x. The problem is therefore reduced to a simple weak problem on x.Asx is a local coordinate, the weak formulation for many beams implies a change of base and a finite element assembly in global coordinates.

The system to solve is then built and solved following the standard finite element method. Now, we want to update, for instance, G n+1 , the other functions F n+1 , H n+1 and K n+1 being known. The weak form becomes then:

Ψ(F n+1 , G * , H n+1 , K n+1 , F n+1 , G n+1 , H n+1 , K n+1 ) = W ext (u * ) - n i=0 Ψ(F i , G * , H i , K i , F i , G i , H i , K i ) (15)
All the integrals in Ψ [see Eq. ( 10)] can be computed numerically except the ones related to E. It remains only a simple 1D problem on G n+1 that is solved very efficiently since there is no derivative depending on E.

The convergence of the PGD algorithm is reached when the last computed term (index n) becomes negligible compared to the most significant one, generally the first (index 0). The error is then computed using:

Err = ||F n || × ||G n || × ||H n || × ||K n || ||F 0 || × ||G 0 || × ||H 0 || × ||K 0 || (16)
where ||.|| is a suitable norm. In practice, every field is discretized, and then, the Euclidean norm is the natural choice. More details on the PGD, and in particular, how to apply the boundary conditions are given in [START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modelling of complex fluids: Part II: transient simulation using space-time separated representations[END_REF].

Results

All the convergence errors for the iterative solving algorithm in the simulations of this work are shown in Fig. 3. These results show good convergence for an error lower than 10 -5 (the convergence error threshold) reaching from 20 to 40 iterations.

For the next sections, all the result of the PGD calibration analysis is represented on two curves. (This is the dual representation of a calibration chart.) It was chosen to represent the ratio of the Young's moduli of the beams and the homo-Fig. 4 Comparison of the calibration curves with results from Andre [START_REF] André | Discrete element method to simulate continuous material by using the cohesive beam model[END_REF] for 33103 beams Fig. 5 Influence of the micro-Poisson's ratio on the calibration Fig. 6 From a coarse to a finer mesh used for the convergence analysis geneous material as well as the homogeneous Poisson's ratio as a function of the ratio (noted radius ratio r ) of the mean radius of the discrete elements in the granular arrangement on the beam mean radius.

Validation

In Fig. 4, calibration results of four different samples obtained by granular packings with the same input parameters (mean number of discrete elements, final coordination numbers) are plotted and compared to the already known calibration analysis in [START_REF] André | Discrete element method to simulate continuous material by using the cohesive beam model[END_REF]. In this case, a special attention was given to use similar post-processing measurement of the global elongation and radius variation after the tensile test than in this previous work. The corresponding four calibration curves are showing the same tendencies that actual discrete simulations show. This also validates the non-influence of the random position of the beams regarding the emerging elastic behavior: The four samples give the same calibration curves. Furthermore, Fig. 5 also shows the non-influence of the microscopic Poisson's ratio of each beam on the calibration chart. Here, only three values of micro-Poisson's ratio are plotted, but the global PGD analysis contains a space discretization of the coordinates ν of 100 points (from 0.01 to 0.9) and has again the same trends. The value 0 was not taken into account in the study as the standard Bernoulli beam kinematic is an hypothesis of this work.

Mesh density influence

In this part, the density of the mesh, i.e., the number of beams in the same volume sample, on the calibration chart is analyzed. Different meshes were used starting with a coarse mesh of 1830 beams and finishing a fine mesh of 33,103 beams (see Fig. 6). An influence of these different meshes on the two calibration parameters is observed in Fig. 7. Concerning the ratio of the modulus, no influence is observed, whereas to obtain the Poisson effect of the emergent structure, an important influence is observed but tends to converge to the mesh with 33,103 beams, which is also a well-known results as the equivalent number of discrete elements is around 10,000.

Coordination number influence

An analysis of the influence of the coordination number (see Fig. 8) has been performed. This parameter corresponds to the connectivity of a bond network and can have a strong influence on the structure behavior. Thus, Fig. 9 shows that the calibration of a lattice structure is also influenced by the nature of the mesh construction, and significantly, both for the ratio of the modulus, where a vertical shift of the curves is observed but also for the Poisson effect where even the trends are changing for small values of radius ratio.

Conclusion

The use of a reduction technique applied to the discrete medium calibration problem is relevant. The very short computation times made it possible to confirm some of the results already known in the literature, but above all to provide the new calibrations charts for homogeneous and elastic behaviors taking into account the discrete nature of the domain. The introduction of nonlinearities such as viscosity and plasticity is the next step of this work, consisting on building new weak formulations in the PGD paradigm.

As a final conclusion of this work, we can say that the reader is highly invited to use the proposed method to perform his own calibration study regarding his own simulation, meaning the use of different discrete bonds or a different protocol to position the nodes. It must be noticed that this computational approach could also be very helpful in lattice architecture materials such as the ones made with 3D printers. Perspectives of this work can be then numerous.

Compliance with ethical standards

Fig. 1

 1 Fig. 1 One of a cylinder sample (3D orthographic view) used in the present work meshed with a granular packing of spheres (left) and its corresponding bond network (right)

  ) for values of E on the interval [E min , E max ], for values of ν on the interval [ν min ,ν max ] and for values of D in the interval [D min , D max ]. E, ν and D are assumed to be constant along the beam.

Fig. 2

 2 Fig. 2 Magnitude displacement field of a 3D cylindrical 33,103 beams network under a global tensile loading

Fig. 3

 3 Fig. 3 PGD error during each simulations presented in this work

Fig. 7 Fig. 8 Fig. 9

 789 Fig.[START_REF] Curti | How to model orthotropic materials by the discrete element method (DEM): random sphere packing or regular cubic arrangement?[END_REF] Convergence results on the number of beams on the calibration curves

Conflict of interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.