P.-A Boucard

D Odièvre

F Gatuingt

A parallel and multiscale strategy for the parametric study of transient dynamic problems with friction

Keywords: parametric study, multiscale computational method, transient dynamics, domain decomposition, contact, friction, parallel processing

The objective of this work is to develop an efficient strategy for the parametric study of dynamic problems involving contacts with friction. Our approach is based on the multiscale LATIN method with domain decomposition. This is a mixed method which deals with the forces and velocities at the interfaces between the different subdomains simultaneously. We propose to take advantage of the capability of the multiscale LATIN method, called the multiparametric strategy, to reuse the solution of a given problem in order to solve similar problems. This strategy has already been applied successfully to a variety of static problems; here, it is extended to dynamics. First, we present the multiscale strategy in dynamics. Then, we show how the multiparametric strategy can be extended to dynamics. We illustrate the capabilities of the method through an academic 3D example and the simulation of a bolted joint.

INTRODUCTION

Modeling and simulation play an important role in engineering and design and give rise to multiple problems, particularly in dynamics when dealing with large assemblies with connections. These connections have a significant impact on the dimensioning process because they are subject to highly nonlinear local phenomena (contacts with friction), particularly in rapid transient dynamics problems, whose proper representation requires very fine meshes [START_REF] Hallquist | Sliding interfaces with contact-impact in large-scale lagrangien computations[END_REF]. In the case of structural assemblies, one's knowledge of the friction coefficients is especially limited. Although the incorporation of a system's parametric uncertainties into a model represents quite a challenge for structural engineers, without such information the structural response could not be analyzed accurately. It is much more difficult to explore this design space taking these parametric uncertainties into account than to solve the general parametric field problem because one must take into consideration nonlinear structural behavior (in our case, unilateral contact with friction). In order to do that, it is necessary to calculate the response of the structure for each set of values of the design parameters [START_REF] Benaroya | Finite element methods in probalistic structural analysis: A selective review[END_REF]. In this case, the design parameters are the friction coefficients and the prestresses. Upon completion of these calculations, the response parameters (maximum stresses, reaction forces...) are obtained and associated with the values of the input parameters. Then, the design space can be defined and explored by sampling parameter values within defined limits. In this context, the choice of an efficient computational method is of vital importance. The perturbation method is one of the main approaches used to obtain such structural responses ( [START_REF] Benaroya | Finite element methods in probalistic structural analysis: A selective review[END_REF][START_REF] Macias | Stochastic Finite Elements and fiability. Application to fracture mechanics[END_REF]; [START_REF] Liu | Transient probabilist systems[END_REF][START_REF] Qiu | Parameter perturbation method for dynamic responses of structures with uncertainbut-bounded parameters based on interval analysis[END_REF] for the dynamic case). In recent years, many researchers have focused on the stochastic finite element method, which is capable of dealing with the parametric uncertainties mentioned above [START_REF] Kleiber | Basic perturbation technique and computer implementation[END_REF]. Some of these methods have been developed in dynamics, such as the Neumann method [START_REF] Zhao | Neumann dynamic stochastic finite element method of vibration for structures with stochastic parameters to random excitation[END_REF] and the spectral method [START_REF] Ghanem | Stochastic finite element : a spectral approach[END_REF]. The POD method [START_REF] Chatterjee | An introduction to the proper orthogonal decomposition[END_REF] can also be used as it allows to build a basis for several sensitivity analysis [START_REF] Ryckelynck | On the a priori model reduction : Overview and recent developments[END_REF]. The objective of the work presented here is to develop a strategy which is suitable for problems involving multiple resolutions, which we will refer to as multiparametric problems. Some multiparametric strategies have been proposed for optimization problems [START_REF] Kim | Multi-resolution multi-scale topology optimization: a new paradigm[END_REF][START_REF] Boucard | A multiscale strategy for structural optimization[END_REF]. These lead to an efficient implementation of the multiresolution optimization technique. In this paper, we consider the case of assemblies of elastic structures connected through frictional contacts and dynamic loading conditions. Contact problems are characterized by constraints such as nonpenetration conditions, and an a priori unknown active area of contact (the area where contact indeed occurs). These conditions lead to stiff nonlinear systems of equations. There are several possible approaches for the resolution of static contact problems [START_REF] Kikuchi | Contact problems in elasticity: a study of variational inequalities and finite element methods[END_REF][START_REF] Zhong | Static contact problems : a review[END_REF][START_REF] Wriggers | Finite element algorithms for contact problems[END_REF]. In most of these approaches, the numerical methods used to enforce the contact constraints are either Lagrange multiplier methods or penalty methods [START_REF] Wriggers | Penalty and augmented Lagrangian formulations for contact problems[END_REF]. Penalty methods [START_REF] Kikuchi | Penalty/finite element approximations of a class of unilateral contact problems[END_REF][START_REF] Armero | Formulation and analysis of conserving algorithms for dynamic contact / impact problems[END_REF] are closely related to the regularization of the contact constraints and are usually formulated in terms of the displacement variables; therefore, they are primal methods. They enable the treatment of contact as a material behavior. Penalty methods can be subject to various numerical difficulties, especially ill-conditioning, when a too large or too small penalty parameter is introduced. Lagrange multiplier methods are dual methods in which the multipliers, which represent the contact reaction forces, are introduced in order to enforce the nonpenetration conditions exactly. Augmented Lagrange multiplier methods [START_REF] Arora | Multiplier methods for engineering optimization[END_REF][START_REF] Klarbring | Mathematical programming and augmented lagrangian methods for frictional contact problems[END_REF][START_REF] Radi | Treatment of the frictional contact via a Lagrangian formulation[END_REF][START_REF] Chabrand | Various numerical methods for solving unilateral contact problems with friction[END_REF] lead to mixed formulations involving both displacement and force unknowns. The numerical resolution schemes underlying both Lagrange multiplier methods and augmented Lagrange multiplier methods are often related to the Uzawa algorithm [START_REF] Arrow | Studies in nonlinear programming[END_REF][START_REF] Simo | An augmented lagrangian treatment of contact problems involving friction[END_REF][START_REF] Champaney | Modular analysis of assemblages of three-dimensional structures with unilateral contact conditions[END_REF].

In this paper, we develop a specific method in order to deal with the problems arising from the fact that the nonlinearities, as described previously, are localized in the connections. The friction law we use is Coulomb's classical law. Our approach is based on a decomposition of the assembly into substructures and interfaces. The interfaces play the vital role of providing an easy and accurate way to model local nonlinearities, such as contact and friction. The problem is solved in each substructure using the finite element method associated with a time integration scheme. An iterative scheme based on the multiscale LArge Time INcrement (LATIN) method developed at LMT Cachan [START_REF] Ladevèze | Nonlinear Computational Structural Mechanics -New Approaches and Non-Incremental Methods of Calculation[END_REF][START_REF] Ladevèze | A micro/macro approach for parallel computing of heterogeneous structures[END_REF][START_REF] Ladevèze | A multiscale computational approach for contact problems[END_REF] is used for the global resolution. This approach leads to a very significant reduction in computation cost for dynamic problems [START_REF] Odievre | A parallel, multiscale domain decomposition method for the transient dynamic analysis of assemblies with friction[END_REF]. The single-scale LATIN method has already been used in dynamics in [START_REF] Boucard | A modular approach to 3D impact computation with frictional Int[END_REF][START_REF] Lemoussu | A 3D shock computational strategy for real assembly and shock attenuator[END_REF][START_REF] Gupta | Mesodynamics of 3D C/C : a dedicated numerical strategy[END_REF]. The multiscale LATIN method is a mixed method which deals with both velocities and forces at the interfaces simultaneously and solves a condensed macroscopic problem in order to accelerate the convergence of the numerical scheme. When only static cases are considered, the LATIN method (without the multiscale approach) can be derived by other means, such as in Lions [START_REF] Lions | On the schwarz alternating method iii. a variant for non-overlapping subdomains[END_REF] or Glowinski and Le Tallec [START_REF] Glowinski | Augmented lagrangian interpretation of the nonoverlapping schwartz alternating method[END_REF].

More specifically, our objective is to calculate a large number of design configurations, each corresponding to a set of values of all the parameters involved in the mechanical analysis. Normally, a full calculation is needed for each set of parameters. Here, as an alternative, we propose to take advantage of the capabilities of the LATIN method and reuse the solution of one problem (for one set of parameters) to solve similar problems (for the other sets of parameters). This approach has already been used for viscoplastic and large-strain problems in [START_REF] Boucard | A multiple solution method for non-linear structural mechanics[END_REF][START_REF] Boucard | A suitable computational strategy for the parametric analysis of problems with multiple contact[END_REF][START_REF] Champaney | Adaptive multi-analysis strategy for contact problems with friction. Application to aerospace bolted joints[END_REF] and is now also being used for structural identification problems [START_REF] Allix | A new multi-solution approach suitable for structural identification problems[END_REF]. It leads to a significant reduction in the total computation cost of determining the design space. The objective of the work presented here is the extension of this method, which has already been largely developed in statics and quasi-statics, to dynamics.

THE MULTISCALE LATIN METHOD

This multiscale domain decomposition method consists of three components: the decomposition of the spatial domain, the separation of the scales and a resolution algorithm. The main features of these three components are described below. The details of the method itself can be found in [START_REF] Ladevèze | A micro-macro and parallel computational strategy for highly heterogeneous structures[END_REF].

Decomposition into substructures and interfaces

An assembly is a set of substructures which communicate with one another through interfaces (see Figures 1(a) and 1(b)). Each interface represents a connection. The substructures and interfaces each have their own variables and equations (admissibility, equilibrium and behavior). The interface associated with two connected substructures Ω E and Ω E ′ is designated by Γ EE ′ . Each interface is a mechanical entity with its own variables and behavior, which depend on the type of connection. Many different types of connections, such as frictional contact, can be modeled with this approach. The interface variables consist of two force fields F E , F E ′ and two dual velocity fields W E , W E ′ (see Figure 1(b)). By convention, F E and F E ′ represent the action of the interface on the substructures, and W E and W E ′ are the velocities of the substructures viewed from the interface. Thus, the interface concept can be easily extended to the boundary, where the displacements, the velocities or the forces are prescribed. In our multiscale vision, the interfaces play the important role of separating the scales: the definitions of the microscopic and macroscopic fields are related to the interface quantities of the substructured problem and are specified prior to any discretization. Then the interface unknowns are divided into :

W E = W m E + W M E and F E = F m E + F M E
where W M and W m denote respectively the macro parts and the micro complements of the velocity field.

A comparison with the coarse problem introduced into the FETI method for the dynamic case by Farhat-Shen-Mandel [START_REF] Farhat | A scalable Lagrange multiplier based domain decomposition method for time-dependent problems[END_REF] can be found in [START_REF] Odievre | A parallel, multiscale domain decomposition method for the transient dynamic analysis of assemblies with friction[END_REF].

(a) The reference problem (b) Exchange between an interface and the adjacent substructures is the virtual space defined by

Ω E Γ EE f d ∂ 1 Ω F E F E W E W E U d F d Ω E ∂ 2 Ω Ω E Γ EE Ω E
U [0,T ] 0 = {u | u ∈ U [0,T ] , u = 0 on ∂ 1 Ω}. ε E (M, t)
is the strain field. The current state of the substructure is characterized by the stress field σ E (M, t), whose associated space is S [0,T ] . The mechanical problem to be solved within each substructure Ω E is: Find the evolutions of the displacement field u E (M, t) and stress field σ E (M, t) such that:

• Kinematic admissibility: ∀t ∈ [O, T ], u E ∈ U [0,T ] -Initial conditions: ∀M ∈ Ω E u E (t = 0) = U 0 E du E dt (t = 0) = V 0 E (1) -Boundary conditions: ∀t ∈ [O, T ], ∀M ∈ Γ EE ′ du E dt Γ EE ′ = W E u E ∂Ω 1 = U d F E ∂Ω 2 = F d (2)
• Equilibrium: ∀t ∈ [0, T ], ∀ u * ∈ U [0,T ] 0 , σ E ∈ S [0,T ] ΩE ρ d 2 u E dt 2 + f d u * dΩ + ΩE T r (σ E ε( u * )) dΩ = E ′ Γ EE ′ F EE ′ u * dΓ (3)
• Elastic behavior: ∀t ∈ [0, T ], ∀M ∈ Ω E σ E = K E ε(u E ) (4)
where K E is the Hooke's operator.

The problem at the interfaces:

The mechanical problem to be solved at each interface Γ EE ′ is:

Find the evolutions of the force fields F E (M, t), F E ′ (M, t) and velocity fields W E (M, t), W E ′ (M, t) such that:

• ∀t ∈ [O, T ], ∀M ∈ Γ EE ′ (F E , F E ′ ) = A Γ EE ′ (W E , W E ′ ) (5)
where the behavior is expressed as an evolution law A Γ EE ′ . This evolution law describes the behavior (perfect, contact, friction...) of the interface Γ EE ′ . More details can be found in [START_REF] Odievre | A parallel, multiscale domain decomposition method for the transient dynamic analysis of assemblies with friction[END_REF].

Resolution strategy: the LATIN method

The LATIN (LArge Time INcrement) method [START_REF] Ladevèze | Nonlinear Computational Structural Mechanics -New Approaches and Non-Incremental Methods of Calculation[END_REF] is a general, mechanics-based computational strategy for the resolution of time-dependent nonlinear problems which operates over the entire time-space domain. It has been applied successfully to a variety of problems [START_REF] Ladevèze | A multiscale computational approach for contact problems[END_REF][START_REF] Ladevèze | On a multiscale computational strategy with time and space homogenization for structural mechanics[END_REF][START_REF] Boucard | A modular approach to 3D impact computation with frictional Int[END_REF][START_REF] Lemoussu | A 3D shock computational strategy for real assembly and shock attenuator[END_REF].

In the case of linear elastic substructures considered here, the solution u E (M, t), σ E (M, t) can be calculated from the boundary values W E (M, t), F E (M, t). Thus, a solution s is characterized by the force and velocity fields on both sides of an interface. The solution of problem s ref is expressed as a set of time-dependent fields within each substructure and at the corresponding interfaces:

s ref = E s E s E = F E (M, t), W E (M, t)

Separation of the difficulties:

The LATIN approach is based on the idea of dealing with each difficulty separately in order not to have to solve a global problem and a nonlinear problem at the same time. The equations are divided into global linear equations and local nonlinear equations, so that s ref = A d ∩ Γ is the intersection of two subspaces:

• A d , the space of the solutions of the linear equations associated with the substructures Ω E : kinematic admissibility, equilibrium, elastic behavior and admissibility of the macroquantities; • Γ, the space of the solutions of the local equations related to the interfaces Γ EE ′ and expressing their behavior.

A two-step iterative strategy:

The LATIN method consists in seeking fields of Γ and A d alternatively along two search directions E + and E -, as shown in Figure 2. Each iteration involves two stages, called the local stage and the linear stage:

Local stage: Given s n = F E , W E ∈ A d , find s n+1/2 = F E , W E such that: s n+1/2 ∈ Γ (interfaces) s n+1/2 -s n ∈ E + (search directions) (6)
Linear stage: Given

s n+1/2 = F E , W E ∈ Γ, find s n+1 = F E , W E such that: s n+1 ∈ A d (substructures) s n+1 -s n+1/2 ∈ E - (search directions) (7) + - ŝ n+1/2 s n+1 s n s ref Γ A d E E Figure 2.

An iteration of the LATIN method

In the case of linear elastic substructures, which we are considering here, the search directions are defined as follows:

s n+1/2 -s n ∈ E + ⇐⇒ F E -F E = k 0 ( W E -W E ) (8)
s n+1 -s n+1/2 ∈ E -⇐⇒ F E -F E = -k 0 (W E -W E ) (9)
where k 0 is a scalar parameter of the method. For the dynamic cases being addressed here, the optimum value of k 0 for a one-dimensional problem was given in [START_REF] Lemoussu | A 3D shock computational strategy for real assembly and shock attenuator[END_REF] and is k 0 = √ ρE, where E is the Young's modulus and ρ the density. k 0 can be viewed as a local impedance of the material.

The convergence of the algorithm toward s ref is controlled using an error indicator η, which is a measure of the distance between the two solutions s n+1 and s n+1/2 (Figure 2).

The local stage: s n+1/2

The local stage consists in building s n+1/2 ∈ Γ knowing s n ∈ Ad. Then, ( s n+1/2 -s n ) must follow the search direction E + defined in Equation 8. We obtain a differential equations problem-type in each point of the interface that is function of its behavior introduced in Equation 5 [START_REF] Odievre | A parallel, multiscale domain decomposition method for the transient dynamic analysis of assemblies with friction[END_REF].

The linear stage: s n+1

The linear stage consists in building s n+1 ∈ Ad knowing s n+1/2 ∈ Γ:

• Macroadmissibility: in order to ensure the admissibility conditions of the macro variables F M E and W M E . • Search direction: the unknowns (W E , F E ) must follow the search direction (Equation 9).

• Equations associated with the substructures: the unknowns (W E , F E ) must verify the dynamic equilibrium and the elastic behavior (Equations 3 and 4).

These equations lead to the resolution of three independent problems. First one have to solved a "micro" problem, in each substructure, then a "macro" problem and at the end a second "micro" problem (see Table I and[30]).

The algorithm and its parallelization

The LATIN method associated with the mixed domain decomposition method is inherently parallelizable [START_REF] Champaney | Large scale applications on parallel computers of a mixed domain decomposition method[END_REF]. In our case, this strategy was programmed in C++ in the framework of the finite element platform developed by H. Leclerc [START_REF] Lecerc | Developments and experimentations on some iterative numerical strategies for solid mechanics on multicore and graphic processors[END_REF]. Libraries such as MPI (Message Passing Interface) for the transfer of information among machines were used in order to enable the use of PC-cluster types of architectures. In order to parallelize the strategy, the first step consists in allocating the substructures and interfaces among the different processors. This is done through the METIS libraries [START_REF] Karypsis | METIS", A Software Package for Partitioning Unstructured Graphs, Partitioning Meshes, and Computing Fill-Reducing Orderings of Sparse Matrices[END_REF], which enable the amount of data which must circulate among the processors to be minimized in order to avoid an excessive decline in speedup.

Initialize

Loop over the substructures (on each processor) Calculation of rigidity and mass matrices of substructure Factorization Calculation and factorization of the macroproblem (on a single processor) Loop over the interfaces (on each processor)

Loop over time (∀t Then, the different operators which are specific to the substructures are constructed on each processor. During the iterative resolution phase, the first and second microproblems are solved. The macroproblem, which has not been parallelized yet, is solved on a single processor. Finally, the local stage is completely parallelized because the interfaces are distributed among the different processors.

∈ [0, T ]) c WE = c W E ′ = 0 b FE = b F E ′ = 0 Iteration n = 1,
The LATIN method consists in performing linear stages and local stages alternatively. The iterations concern the whole time interval, i.e. a solution over the whole time interval is calculated at each iteration of the method. In the linear stage, an incremental formulation is used to solve problems within the substructures over the whole time interval. Thus, only two time steps of a substructure's fields (displacement, velocity and acceleration) need to be stored: the current time step being calculated and the previous time step. Conversely, the interface quantities ( F , W , F, W ) must be stored over the whole time interval. The fields are saved after the last iteration of the method for post-processing purposes. Algorithm (Table I) shows the key steps of an iteration of the multiscale method.

THE STRATEGY FOR THE PARAMETRIC STUDY

The strategy proposed consists in calculating design configuration where each design configuration corresponds to a set of values of all the variable parameters (friction coefficients, prestresses) which are introduced into the mechanical analysis.

At each iteration, the LATIN method leads to an approximate solution to the problem over the whole time interval. Therefore, the trick is to reuse one approximation (associated to one set of values of all the design parameters) to find the solution to another design configuration (another set of the design parameters) similar to the one for which it was calculated in the first place. Our strategy for the parametric study uses the fact that the LATIN algorithm can be initialized with any solution (usually an elastic solution) provided that it verifies the admissibility conditions. Therefore, the key to our technique is to initialize the process associated with a new design configuration using the results of the calculation carried out on the preceding set of values of the design parameter (Figure 3). In this scheme, Γ i correspond to one set of parameters and Γ i+1 to another one. S i (equal to S ref for Γ i ) is the solution of the first set of parameters and correspond to the initialization of the LATIN algorithm for the new ones. In this manner, a first approximation of the solution to the new design with a strong mechanical content is immediately available from the start. Note that as the previous computation is expected to have converged, the solution of the local stage should be identical to the one provided by the linear stage. One can use any of them to initialize the algorithm.

In this particular case of elastic structures in contact, the interfaces play a vital role: they enable one to initiate the calculation on the new design configuration without having to save all data on the substructures as well as to search for the solution of the new design configuration with an initial solution well-suited to the target problem. In the best-case scenario, only a few iterations are necessary: the solution to the problem is obtained at low cost. If the solutions to the design configurations are close enough, the latter can still be derived at a significantly lower cost than by using a full calculation. For the parametric study presented herein, we just change the parameters between iterations. Thus, the new computation is initialized by the solution to the previous one. If the parameters change slowly, the two solutions are close and only a few iterations are needed to reach convergence in the new calculation.

Initialize

Loop over the substructures (on each processor) Calculation of matrices KE, ME, ... Factorization Calculation and factorization of the macroproblem (on a single processor Loop over the interfaces (on each processor)

Loop over time (∀t

∈ [0, T ]) c WE = c W E ′ = 0 b FE = b F E ′ = 0 Definition of

the sets of values of all the design parameters

Loop over the sets of values of all the design parameters (k = 1, 2, ...)

Restoration of the quantities at the interfaces Iteration until convergence Linear stage (Table I) Local stage (Table I) Convergence test Save the interface solution for the k th parameter set Table II. Algorithm of the parametric study using the LATIN method

FIRST EXAMPLE: AN ACADEMIC PROBLEM

In order to test the method described previously, let us consider a simple 3D example consisting in the propagation of a compression wave in a bar composed of three parallelepipeds (Figure 4). Each parallelepiped is 30 mm-long and 15 mm-wide. The material properties are: Young's modulus 150 GP a, mass density 7, 800 kg/m 3 and Poisson's ratio 0.45 for the central parallelepiped, and Young's modulus 200 GP a, mass density 7, 800 kg/m 3 and Poisson's ratio 0.3 for the other two. Contact with friction is assumed to take place between theses parallelepipeds with friction coefficients f 1 for Interface 1 and f 2 for Interface 2 (Figure 4). First, the bar was subjected to a static loading in the form of a prescribed displacement applied to the upper parallelepiped (U = -0.1mm along the z-axis). This loading consisted of a normal pressure at each contact interface, which corresponds to a static prestress in the structure. Then, a pressure was applied laterally to the central parallelepiped going from zero to a maximum value of 50 M P a over 10 µs. This pressure resulted in slipping of the middle parallelepiped along the x axis by an amount which depends on the value of the two friction coefficients. For the purpose of this test, each parallelepiped was decomposed into several substructures. The duration of the test was 0.2 ms and we used 50 time steps of 4 µs each. The objective of this example was to study the influence of the friction coefficients on the slippage of the central parallelepiped. Each friction coefficient was given 11 different values from 0.1 to 0.2 in increments of 0.01. At the end of each simulation, the amounts of slippage at points A and B were saved (Figure 4). The parameters of the simulation (material properties, loading intensities, friction coefficients,...) were chosen such that the pressure on the central parallelepiped would induce slipping during the study and the parallelepipeds would stick at the end of the study. Figure 5 shows the displacement and velocity of point A (along x) obtained for f 1 = 0.1 and f 2 = 0.1, with the fields of the upper parallelepiped in solid black lines and the fields of the middle parallelepiped in dotted black lines. The slippage of point A is equal to the difference between the two displacement curves.

Each parallelepiped was divided into 16 substructures consisting of 10 × 10 × 10 8-node linear elements. The simulation was carried out using 16 processors. This example of a parametric study required 121 resolutions of a nonlinear problem meshed with 144,000 DOFS over 50 time steps. We studied the slippage of points A and B of the two contact interfaces (Figure 4). Figure 6 shows the response surfaces of the parametric study, i.e. the evolutions of the slippage at points A and B as functions of the two friction coefficients. These two surfaces illustrate the influence of the friction coefficients. For example, the slippage of point A depends on the friction coefficient f 1 alone, and this dependence is highly nonlinear. Moreover, above f = 0. For each value of f 1 starting with f 1 = 0.1, we applied the parametric strategy to perform the 11 simulations associated with each value of f 2 starting with f 2 = 0.1, i.e. we applied the multiparametric strategy 11 times for each value of f 1 . In order to optimize also the first calculation of each multiparametric strategy (f 2 = 0.1), we saved the solution of the previous simulation to initialize the next multiparametric strategy (for the next value of f 1 ). For this example, the convergence criterion of the iterative algorithm was set to 2.10 -4 . Figure 7 shows the evolution of the LATIN error indicator during the first multiparametric strategy (f 1 = 0.1). The number of iterations and the CPU time associated with each simulation (for each value of f 2 ) are shown in Table III. The error indicator increased every time the friction coefficient changed, but to a lesser extent than in the first initialization thanks to the multiparametric strategy. Thus, one can observe that the number of iterations required to achieve convergence decreased throughout the parametric study. With the multiparametric strategy, both the number of iterations and the computation cost decreased because the results of the successive computations were very similar (as shown by Figure 6). Table III. Numerical costs of the first multiparametric strategy (f1 = 0.1)

In Table III, the number of iterations decreases continuously. The gain is very important for the last values of f 2 (0.18...0.2). That mean that the improvement does not only depend on the method, but so on the nature of the solutions. Indeed the solutions for f 2 equal to 0.18.. 0.2, only involve stick statuses over the interface, which is easier to solve than mixed solutions with stick and slip statuses. This particular point of the strategy is highlight by testing a parametric study starting for f 2 from 0.2 to 0.1 by a decreasing way. Figure 8 shows the evolution of the LATIN error indicator during this new multiparametric strategy (f 1 = 0.1). The total number of iterations in this case is 753 (to compare with the 570 iteration of the first parametric study). In terms of CPU time, for this second study total time is 102.4 min (to compare with 71 min). We can also compare the gain: 2.9 for the second parametric study and 4.1 for the first one. This result shows that it is possible to improve the parametric strategy by seeking the best point for reinitializing the LATIN algorithm. Indeed, at a given time in the parametric study, one can have access to several already calculated solutions. Thus, in performing a new calculation, one can reasonably assume that among the previously calculated solutions one is likely to maximize the cost reduction for this new calculation. In this example the i th calculation is reinitialized from the results of calculation i -1, but one can easily imagine that it could be better if the i th calculation is reinitialized from the results of the closest calculation, in the sense of a distance in the parametric space [START_REF] Boucard | Multidisciplinary Design Optimization in Computational Mechanics[END_REF].

Error indicator

Iterations We compare the results with those obtained with the industrial explicit dynamic code LS-DYNA [START_REF] Hallquist | [END_REF]. To obtain the same quality in terms of contact and sliding at the interface, we had to use a finer mesh in LS-DYNA. For each parallelepiped we use 40x40x80 elements (the number of elements in each direction is doubled in comparison with the mesh used with our own software: 20x20x40 elements). Indeed, if we use the same mesh with LS-DYNA, the behavior of the interface is not verified exactly: the parallelepipeds are not in contact as they should be, they interpenetrate at the end of the calculation. We compare the results in terms of velocity of point A (figure 9). It is obvious that the solutions are similar and that a finer mesh is necessary for LS-DYNA to obtain the solution. The total computation time with LS-DYNA3D on one computation and one processor (f 1 = 0.1 and f 2 = 0.1) is about 7h and the memory used about 850 MBytes. In comparison, the computation time with our own software is about 4h and 50 minutes and the memory used is about 1800 MBytes. Therefore, the computational times associated with a solution quality are comparable. Our research software is not optimized in terms of memory occupation. Nevertheless, one limitation of the proposed approach is a limited number of time steps. We compare the memory used for several time discretization : 50, 100, and 200 time steps. For this 3 values, the memory used is respectively : 1800, 2400 and 3600 MBytes. Contrary to a classical incremental scheme for which the volume of computations roughly doubles when doubling the time steps, here, not only this number of operations doubles, but also the storage of the iterates on which the computations are performed. Moreover, when we double the number of time step, the total time of the analysis nearly double also.

f 2 = 0.12 f 2 = 0.13 f 2 = 0.1 .2 f 2 = 0.2 f 2 = 0.18 f 2 = 0.19 f 2 = 0.11 f 2 = 0.17 f 2 = 0.16 f 2 = 0.14 f 2 = 0.15

APPLICATION TO A BOLTED JOINT USED IN THE AEROSPACE INDUSTRY

Let us now consider the example of a bolted joint. Damping in such joints is sensitive to friction, pretension in the bolts, etc... These parameters are naturally scattered, and a full test campaign to evaluate their actual influence would involve a very large number of specimens. Such studies in quasi-statics have been reported in [START_REF] Caignot | Virtual testing for the prediction of damping in joints[END_REF]. This example concerns the propagation of a compression wave in the bolted joint. This 3D assembly is a joint between two sandwich composite structures consisting of two metal parts and three bolts (Figure 10). For the purpose of this example, the core and the composite material of the sandwich were assumed to be homogeneous, elastic and isotropic media. The interfaces between the sandwich composites and the metal parts were assumed to be perfect. Frictional contact was taken into account both at the interfaces between the bolts and the metal parts and at the interface between the two metal parts. First, we applied a static loading and prestresses in the bolts. A static calculation was performed in order to define the initial conditions of the dynamic analysis. For the dynamic calculation, the loading consisted of the prescribed velocity shown in Figure 10. The calculations were performed using a 450,000-DOF mesh for the whole assembly. In order to use our approach, we divided the different parts into several substructures. Figure 10 shows a portion of this decomposition. We chose to use substructures with similar numbers of DOFs in order to achieve a well-balanced load among the different processors and a good speedup of the parallel algorithm. We used 194 substructures with about 2,500 DOFs each. An implicit scheme with 39 time steps of 6 µs each was used for the time integration for a total duration of 240 µs. The objective of this example was to study the influence of the friction coefficients and of the prestress in the bolt on the slippage between the two metal parts. The friction coefficient was given 11 different values from 0.2 to 0.4 in increments of 0.02, and the prestress in the middle bolt was also given 11 different values from 0.02 mm and 0.04 mm in increments of 0.002 mm. The prestress in the bolt was represented as a relative displacement between the body and the head of the bolt. At the end of each simulation, we saved the amount of slippage of point A along x (Figure 10). Figure 11 shows the displacement U x and the velocity V x of point A along x obtained for a friction coefficient equal to 0.2 and a prestress equal to 0.02 mm, with the fields of the right-hand side in solid black lines and the fields of the left-hand side in dotted black lines. The slippage of point A is equal to the difference between the two displacement curves. This example of a parametric study involved 121 resolutions of a nonlinear problem with a 450, 000-DOF mesh and 39 time steps. The simulation was carried out using 16 processors.

Figure 12 shows the response surfaces of the parametric study: the evolution of the slippage of point A as a function of the friction coefficient and of the prestress at the beginning and at the end of the study. These two surfaces illustrate the influence of the parameters on the response of the structure. For example, the dependence of the slippage on these parameters at the beginning of the study is highly nonlinear, with the slippage being equal to zero for part of the variation range of the parameters.

For each value of the prestress starting with 0.02 mm, we applied the multiparametric strategy to calculate the 11 simulations associated with each value of the friction coefficient starting with 0.2, i.e. we applied the multiparametric strategy 11 times for each value of the prestress.

In order to optimize also the first calculation of each multiparametric strategy (f = 0.1), we saved the solution of that first simulation to initialize the next multiparametric strategy (for the next value of the prestress). Figure 13 shows the number of iterations carried out for each parameter set. While the first calculation required more than 700 iterations, each of the subsequent calculations took less than 100 iterations thanks to the multiparametric strategy.

Table V shows a comparison of the computation costs of the different strategies. The first CPU time (Direct LATIN, sequential) corresponds to 121 times the CPU time of the first calculation (700 iterations) carried out using a single processor, the second CPU time (Direct LATIN, parallel) corresponds to 121 times the CPU time of the first calculation carried out using 16 processors and the third CPU time corresponds to the multiparametric strategy carried out using 16 processors. For this example, the gain in CPU time obtained thanks to parallelism using 16 processors was 9.5 (1354/141.1). The gain in CPU time obtained thanks to the multiparametric strategy was 5.8 (comparison between direct LATIN parallel and ). Parallelism and the multiparametric strategy together enabled us to carry out the multiparametric study in a reasonable time compared to the direct approach, with a global gain of about 55.

CONCLUSION

We presented the extension of the multiparametric strategy based on the LATIN method to dynamic problems involving complex 3D structural assemblies. The LATIN method is based

Figure 1 .

 1 Figure 1. Decomposition of the reference problem into substructures and interfaces

Figure 3 .

 3 Figure 3. Schematic view of the reuse of a previous LATIN solution

FFigure 4 .

 4 Figure 4. The numerical example

2 ,Figure 5 .Figure 6 .

 256 Figure 5. Evolution of the displacement and velocity of point A (along x), upper parallelepiped in solid black lines and middle parallelepiped in dotted black lines

Figure 7 .

 7 Figure 7. Evolution of the error indicator for f1 = 0.1

Figure 8 .

 8 Figure 8. Evolution of the error indicator for f1 = 0.1

Figure 9 .

 9 Figure 9. Evolution of the displacement of point A (along x), upper parallelepiped in solid black lines and middle parallelepiped in dotted black lines

Figure 10 .

 10 Figure 10. The model of the joint

Figure 11 .

 11 Figure 11. Evolution of the displacement and velocity of point A along x, right-hand side in solid black lines and left-hand side in dotted black lines

Figure 12 .Figure 13 .

 1213 Figure 12. Response surface -slippage of point A at the beginning and at the end of the study

Table I .

 I Algorithm -The micro/macro LATIN method

	First microproblem:
	Solve the first micro problem at time step t
	Macroproblem (for a single processor):
	Determination of the macro part of the solution a time step t
	Loop over the substructures (on each processor)
	Second microproblem:
	Solve the second micro problem at time step t
	(operators are the same as those used in the first microproblem)
	Local stage -Loop over time (∀t ∈ [0, T ])
	Loop over the interfaces (on each processor)
	Determination of ( c WE, b FE) given (WE, FE)
	Convergence test

2, ... until convergence Linear stage -Loop over time (∀t ∈ [0, T ])

Loop over the substructures (on each processor)

Table V .

 V TableIVshows a comparison of the computation costs of the different strategies. The first CPU time (Direct LATIN, sequential) corresponds to 121 times the CPU time of the first calculation (194 iterations) carried out using a single processor, the second CPU time (Direct LATIN, parallel) corresponds to 121 times the CPU time of the first calculation carried out using 16 processors, and the third CPU time corresponds to the multiparametric strategy carried out using 16 processors. For this example, the gain in CPU time obtained thanks to parallelism using 16 processors was 11 (588/53.4). The gain in CPU time obtained thanks to the multiparametric strategy was 4.4 (comparison between the direct LATIN parallel and the Multiparametric LATIN parallel: 53.4/12.1). Parallelism and the multiparametric strategy together enabled us to carry out the multiparametric study in a reasonable time compared to the direct approach, with a global gain of about 48. Cost comparison of the different calculations

	Calculation	Cost (hrs) Cost (days)	Gain vs. Direct
	Direct LATIN (sequential, estimated)	588	24.5	LATIN sequential
	Direct LATIN (parallel, estimated)	53.4	2.2	11
	Multiparametric LATIN (parallel)	12.1	0.5	48

Int. J. Numer. Meth. Engng 20; :-Prepared using nmeauth.cls

on three components: a mixed decomposition of the structure (which provides significant modularity to the description of the problem), a multiscale description and an iterative resolution scheme. Using scale separation, a macroproblem is defined at the interface by condensing the microproblem, which accelerates the convergence of the iterative algorithm. This method is particularly well-suited to problems involving frictional contact. The strategy is also fully parallelized. Its implementation into a cluster architecture leads to good behavior of the parallel algorithm and enables one to treat complex real-life structural assemblies with very large numbers of degrees of freedom.

The multiparametric strategy proposed here is based on the LATIN method and, more specifically, on its capability to reuse the solution of a given problem in order to solve similar problems. Our first numerical example showed the very good behavior of the algorithm in the case of multiple resolutions in the analysis of dynamic problems with contact and friction.

The second example concerned a real 3D joint. For this assembly, parametric studies were carried out on the values of the connection parameters (friction coefficient, prestress). The two examples presented show that the algorithm can be very efficient numerically. The parametric studies were carried out about 50 times faster than with the direct approach. Thanks to the multiparametric strategy and parallelism of the LATIN method, it is possible to carry out such parametric studies in reasonable time.