S Marguet

P Rozycki

L Gornet

A Rate Dependent Constitutive Model for Carbon-Fiber Reinforced Plastic Woven Fabrics

Keywords: composite woven fabrics, delayed damage mesomodel, viscoplasticity, identification, direct search method, dynamic loadings

This paper deals with the modelling until rupture of composite structures made of carbon-fiber/epoxy-resin woven fabrics submitted to dynamic loadings. The model is built at the mesoscale of the elementary ply. It takes into account the slightly nonlinear brittle behavior of the fibers under tensile solicitations, their nonlinear behavior in compression as well as the strongly nonlinear and irreversible behavior of the ply in shear. Strain rate effects are also introduced and special attention is paid to the objectivity of the model in the context of finite element calculation. Therefore the choice of a delayed damage mesomodel coupled with viscoplasticity is made. In order to identify the values of the parameters of the model, an optimization procedure based on a gradient-free direct search method has been developed. As a logical procedure to this study, the models ability to avoid strain localization and mesh dependence is then checked on simple uniaxial examples. The last part of this paper is devoted to structural calculation. The results of the simulations of both the impact on composite plate and the crushing of thin-walled tube demonstrate the capability of the model to reproduce observed physical phenomena.

INTRODUCTION

Nowadays composite materials take a growing importance in the design and manufacture of multihull oceanic sailing ships. Their ability to provide significant stiffness combined with limited weight perfectly matches the more and more constraining needs of competition. Designing the ships with reduced safety margins is also a good way to improve the performances measured in term of the "velocity/wind power" ratio. Still, during great oceanic races, structural failures caused by extreme climatic conditions lead lots of ships to give up the race. The knowledge of the material behavior under dynamic loadings such as slamming or wave impacts then appears to be determinant for the design of safe structures.

Carbon-fiber/epoxy-resin woven fabrics composite materials have been studied for many years. Their linear orthotropic elastic behavior is now well known as can be seen in [START_REF] Herakovich | Mechanics of Fibrous Composites[END_REF]. The use of homogenization procedures to find the equivalent mechanical properties of stacking sequences of these materials has allowed the development of huge numerical models dedicated to the design of the sailing ships [START_REF] Gornet | Interactions modèles expériences sur des âmes nids d'abeilles Nomex R : Application au design d'un voilier multicoque de course océanique[END_REF]. However, in a context of competition, the knowledge of the elastic properties of the constitutive material is not enough for an optimized design of the structures. Nonlinear effects which correspond to the local deteriorations of the composite materials and eventually lead to the failure of the whole ship have to be dealt with.

Since the end of the 80s, important work has been performed on the understanding and modelling of the damage effects. Four principal modes of degradation of the elementary ply were isolated:

• the tensile fragile rupture of fibers in tension;

• the local buckling of the fibers in compression (responsible of the observed nonlinearity); • the microcracking of the resin under shear or transverse loadings (which explains the nonlinearity and the irreversible strains); • and the decohesion between the fibers and the matrix.

To these modes can be added the delamination between the different staked layers of the composite which also plays an important role in the failure of macroscopic structures. Based upon these observations, several models written at the homogenized mesoscale of the elementary ply have been developed and especially by Ladevèze and his co-workers [START_REF] Ladevèze | Damage modelling of the elementary ply for laminated composites[END_REF][START_REF] Ladevèze | A mesomodel for localization and damage computation in laminates[END_REF][START_REF] Ladevèze | A computation damage mechanics approach for laminates: Identification and comparison with experimental results[END_REF]. In a first step, damage mesomodels have demonstrated their ability to take into account the loss of stiffness observed in experiments. The adjonction of plastic models with isotropic hardening has also enabled the representation of irreversible strains. Yet, in the context of finite element calculation, local softening damage models are subjected to the strain localization phenomenon which results in a state of damage dependent from the spatial discretization. To circumvent this issue, several strategies have been tested. Nonlocal models [START_REF] Bazant | Nonlocal continuum damage, localization instability and convergence[END_REF], gradient dependent models [START_REF] Zhang | A new algorithm for numerical solution of dynamic elasticplastic hardening and softening problems[END_REF] or delayed damage mesomodels [START_REF] Ladevèze | A mesomodel for localization and damage computation in laminates[END_REF] are some of them. The aim is to introduce an internal length linked, for the delayed damage mesomodel, with a maximal rate of damage. As a consequence, the damage cannot immediately reach its critical value, and so, information can be passed through the mesh. Numerous studies have then been performed in order to check the ability of delayed damage mesomodels to avoid strain localization [START_REF] Suffis | Damage model with delay effect: Analytical and numerical studies of the evolution of the characteristic damage length[END_REF]. By the way, the identification of the values of the parameters and the validity of such models on large strain rate ranges still remain a matter. More recently, a model resting on classical damage and plastic theories was developed by [START_REF] Rozycki | Contribution au développement de lois de comportement pour matÈriaux composites soumis à l'impact[END_REF] to deal with the strain rate dependency of the glass-fiber/epoxy-resin composite materials. Satisfying results were obtained but the mesh independence was not ensured. The object of this study is to propose a local mesomodel dedicated to carbon-fiber/epoxy-resin woven fabrics submitted to high strain rate loadings. An optimization procedure based on a gradient-free direct search method (the pattern search algorithm) is also implemented in order to perform the identification of the parameters of the model. To check the behavior of the model in finite element calculations, the example of uniaxial tensile tests on bars was taken. The ability of the model to reproduce the failure of composite structures was eventually shown thanks to the simulation of the impact on a woven fabric composite plate, and of the crushing of a thin-walled composite tube.

MODELLING OF THE FABRIC

The present model is written at the mesoscale of the layer which brings a good compromise between the scales of both the constitutive materials (fibers and resin) and the structure. Moreover, and thanks to the thin aspect of most of the composite structures, the assumption of plane stress is made. At last, the effects of temperature are neglected and only small strains are considered. From now on, the subscripts 1 and 2 stand for the orthogonal fiber directions while 3 stands for the direction normal to the woven fabric ply.

Thermodynamic Potential and State Laws

To deal with both the reversible and the irreversible effects, the classical split of the total train t is assumed:

t = e + vp ( 1
)
where e stands for the elastic strain which is an internal variable and vp stands for the viscoplastic strain. As viscoplasticity only occurs in shear, where fibers cannot prevent viscoplastic flow, the different strain vectors take the form t = [ t To take into account the irreversible strains, the loss of rigidity, the different tensile and compressive behavior in fiber direction and the strain rate dependency, the thermodynamic potential is derived from the Helmholtz's volumic free energy ψ, and written as follows:

                              
where:

• ˙ e must be considered as a parameter, it allows the elastic modulus to increase with the strain rate as observed in experiments;

• 0 ii are the coefficients of the stiffness matrix associated with the undamaged material which is supposed to be the initial one: 0 11 ( ˙ 11 ) To complete the description of the state potential, the following functions have been introduced:

= E 11 •V 11 ( ˙ 11 ) 1-ν 12 ν 21 ; 0 22 ( ˙ 22 ) = E 22 •V 22 ( ˙ 22 )
• P + (σ t-t ) = 1 if σ t-t ≥ 0,
0 otherwise (this is an explicit pragmatic way to check for the tensile or compressive state of the material from the state of stress at one step time before the actual time t); 2 do the same for the nondiagonal terms of the stiffness matrix and allows to conserve its symmetry.

• P -(σ t-t ) = 1 if σ t-t < 0, 0 otherwise; • V ii ( ˙ e ii ) = (1 + v ii • ln( | ˙ e ii | ˙ ref ))
Here it has to be noticed that the introduction of all the preceeding functions preserves the symmetry and the positive definite property of the stiffness matrix.

The state laws can now be derived from Eq. ( 2) leading to the expressions of σ, the Cauchy's stress vector dual of the elastic strain, Y the thermodynamic forces associated with the damage state and R the stress linked to the isotropic hardening:

     σ = ∂ψ ∂ e Y = -∂ψ ∂d R = ∂ψ ∂ p (3)
Taking the volumic energy due to isotropic hardening under the form h( p) = Q β+1 • p β+1 (with Q and β material parameters), the preceeding expressions once developed, the state laws become:

                                               σ 11 = (1 -d 11
R = Q • p β (4)

Dissipation Potential and Evolution Laws

To obtain the evolution laws, the existence of a dissipation potential that is a function of the stress variables is assumed [START_REF] Lema Ître | Mechanics of solid materials[END_REF]. The following expression is adopted in order to deal with both the rate dependent irreversible effects and the degradations:

φ(σ; Y ; R; (d)) = φ vp (σ; Y ; R; (d)) + φ d (σ; Y ; R; (d)) (5)
where φ vp stands for the dissipated energy per unit volume due to viscoplasticity and φ d is the volumic energy dissipated during the microcracking of fibers and resin.

Introducing the Perzyna's viscoplastic theory [START_REF] Perzyna | Fundamental problems in viscoplasticity[END_REF], and remembering that irreversible strains only occur in shear, the viscoplastic dissipation potential is:

φ vp (σ 12 ; R; (d 12 )) = K m + 1 f (σ 12 ; R; (d 12 )) K m+1 + (6)
with f , the yield surface, K and m parameters that drive the rate sensitivity of the material, and d 12 the damage variable associ-ated to the microcracking of the resin in shear. The choice of such a potential is motivated by two different aspects:

• the dependence of the woven fabric to the loading rate (material aspect that should come from the viscoelastic nature of the epoxy resin); • and the viscoplastic regularization effects on strain localization phenomenon (numerical aspect).

The coupling between viscoplasticity and damage is realized through the introduction of the effective stress in the expression of the yield surface f :

f (σ 12 ; R; (d 12 )) = f ( σ12 ; R) = | σ12 | -R -σ y (7)
with σ12 = σ12 1-d 12 the effective shear stress which stands for the stress that should be applyied to the undamaged material to obtain the same state of deformation as the one obtained with σ 12 on the damaged material.

From Eq. ( 6) it is now possible to determine the evolution laws of both the inelastic strain vp 12 and the cumulated plastic strain p:

˙ vp 12 = ∂φ ∂σ 12 = λvp ∂ f ∂σ 12 ṗ = -∂φ ∂R = -λvp ∂ f ∂R (8)
where:

λvp = f K m + (9)
corresponds to a Lagrange's multiplier that drives the viscoplastic flow and is named viscoplastic multiplier. The function x + stands for the Heaviside's function: x + = x if x > 0 and 0 elsewere. If f ≤ 0 no viscoplastic flow can occur and then ˙ vp 12 and ṗ are both equal to zero. If f > 0 then viscoplasticity appears and the material flows in the direction normal to the elastic surface. To completely define the dissipation potential φ, the choice of the damage dissipation potential φ d has to be carried out. Assumption is made that a dissipation potential, depending only on the damage state, exists. The following evolution laws for the damage variables are then stated:

ḋii = 1 τ ii • 1 -e -α ii • z ii + if d < d c
ii ; 0 elsewhere, and ii = {11; 22; 12} [START_REF] Lema Ître | Mechanics of solid materials[END_REF] where α ii and τ ii are material parameters and z ii stand for the functions that drive the appearance of the damage and its rate.

In fiber direction,

z ii = √ Y ii - √ Y o ii √ Y c ii - √ Y o ii -d ii (ii = 11, 22) whereas in shear direction z 12 = ln(Y eq )-ln(Y o 12 ) ln(Y c 12 )-ln(Y o
The evolution equations of the damage variables d define a damage mesomodel as was introduced by Ladevèze and his coworkers [START_REF] Ladevèze | A mesomodel for localization and damage computation in laminates[END_REF]. It is particularly interesting for its ability to avoid the strain localization phenomenon which results in mesh dependence of the results (see [START_REF] Suffis | Damage model with delay effect: Analytical and numerical studies of the evolution of the characteristic damage length[END_REF] for more details). The material parameter τ ii represents the inverse of the maximum damage rate and is linked to an internal length supposed to be the thickness of the layer for the mesomodel. As a consequence, the damage has to be constant through the out-of-plan direction of a layer of the composite ply.

A last key point in the construction of the model is to check the positivity of the intrinsic dissipation, coming from the second law of thermodynamics, and definied by the equation:

D = σ T ˙ vp -Y T ḋ -R ṗ ( 11
)
As no simple expression can be obtained at this point, particular attention should be paid to the verification of the second thermodynamic law.

Implementation and Numerical Aspects

This section deals with the following problem: "given the state of the material on the whole model at time t n , known the step time δt n+1 = t n+1t n and the prescribed conditions t n+1 , the state of the material at time t n+1 is to be found."

The integration is performed in two major steps which are:

1. integration of the delayed damage model in fiber directions; 2. integration of the coupled viscoplastic/delayed damage mesomodel in shear thanks to a modified Radial Return Mapping algorithm [START_REF] Simo | Computational inelasticity[END_REF].

For the behavior in fiber direction, no special difficulty occurs so that the integration method is not detailed. The damage state at time t n+1 is simply computed implicitly by the mean of a Newton Raphson's procedure considering whether the type of sollicitation at time t n is tension or compression. At the end of this first part, the damage variables d t 11 and d t 22 are updated and are used to check for the failure of the fibers in tension. Under compressive state, a simple criterion based on a critical stress drives the rupture.

The next step deals with the integration of the coupled viscoplastic/delayed damage mesomodel. The system of equations to solve comes from the implicit backward Euler's scheme applied to the variables from strain space (d, ḋ, vp and ˙ vp stand for d 12 , ḋ12 ,

d n+1 = d n + δt n+1 ḋn+1 vp n+1 = vp n+1 + δt n+1 ˙ vp n+1 p n+1 = p n + δt n+1 ṗn+1 (12)
where the evolution laws give the last right hand term of Eqs. [START_REF] Simo | Computational inelasticity[END_REF]. To condense these expressions, it is introduced the vector of unknows = (d; vp ; p) T which leads to:

n+1 = n + t n+1 ˙ n+1 (13)
As the expression of ˙ is a non linear expression of , the solution of the system is not trivial. In order to solve it, an incremental procedure of Newton Raphson is once again employed. At the increment k of the iterative procedure a solution to the problem at time t n+1 , noted k n+1 , is supposed to be known (for the first iteration k = 0, the approximation 0 n+1 = n is used). This hypothesis leads to the formation of a residual vector:

Res k n+1 = k n+1 -n -t n+1 ˙ k n+1 ( 14
)
If the norm of the residual is greater than a given tolerance, the solution is not acceptable, and an iteration has to be performed. Linearizing the Eq. ( 14) conducts to the following expression of the increment in strain variables:

δ k n+1 = - ∂Res ∂ k n+1 -1
Res k n+1 [START_REF] Lewis | Direct search methods: then and now[END_REF] with [START_REF] Johnson | Computational methods for predicting damage in composite structures[END_REF] which clearly shows the coupling between viscoplasticity and dalayed damage. From the preceding two equations, it is then possible to infer a new estimation of the state of the material at time t n+1 with:

∂Res ∂ k n+1 =    1 -t n+1 ∂ ḋ ∂d -t n+1 ∂ ḋ ∂ vp -t n+1 ∂ ḋ ∂ p -t n+1 ∂ ˙ vp ∂d 1 -t n+1 ∂ ˙ vp ∂ vp -t n+1 ∂ ˙ vp ∂ p -t n+1 ∂ ṗ ∂d -t n+1 ∂ ṗ ∂ vp 1 -t n+1 ∂ ṗ ∂ p    k n+1
k n+1 = k n+1 + δ k n+1 and k = k + 1 ( 17
)
which leads to the formation of a new residual vector until the tolerance is reached. Once n+1 has been obtained, the state laws allow the computation of the stress variables and the next step can be performed. Finally, the model is implemented in the commercial finite element code Abaqus Explicit thanks to a user subroutine VUMAT dedicated to the material behavior.

PARAMETERS IDENTIFICATION OF THE MODEL

The identification procedure consists in determining, thanks to experimental data, the values of the mathematical constants involved in the constitutive relations of a model. To perform the identification, two major different methods depending on the complexity of the material and its associated model can be used. First, when the following constraining conditions can be met:

• possibility to associate a physical or phenomenological interpretation to a given parameter and a given test (Young modulus to account for linear elasticity during a uniaxial test . . . ); • possibility to make the distinction beween the different physical phenomena (inelastic strain, loss of rigidity, rate dependence . . . ); • ability to make simplifying assumptions;

• capacity to use different experimental curves in a sequential order to find the parameters; a sequential identification procedure can be developed. This method is quite empirical and gives very satisfying results for numerous models. However, in this context, the identification procedure is specific to a given material and its associated constitutive model and it is not possible to consider strongly coupled phenomena.

In our case, all the aspects of the behavior are imbricated and it is not a simple matter to isolate the influence of each parameter in the global rate dependent response of the specimen. One way to circumvent this difficulty is to consider the identification procedure as an optimization problem. It constitutes the second major class of methods.

The basic idea consists in minimizing the gap between experimental data and results coming from numerical simulations. In order to achieve this goal, the first step is to define the variables of comparison. For example, one can compare the experimental and numerical stress versus strain curves or, with recent optic mesurement methods, fields of displacements. The next step is then to create a "cost function" that estimates the error made by the model while simulating the observed phenomena. Such an indicator can be an error within the meaning of least squares. Finally, from the knowledge of some sets of parameters and of the quality of their fit, a new set has to be determined in order to improve the numerical solution.

In this study, the identification procedure is based on a Pattern Search Algorithm which comes from the direct search (gradientfree) methods and enables to find minima of functions. This algorithm has been chosen for three principal qualities:

• its ability to deal with numerous variables;

• its robustness (if an incompatible sets of parameter leading to non real cost result is encountered, the algorithm goes through the trap); • and for its capacity to check a wide area in the parameters space which limits the starting point dependence in comparison with other methods such as descent ones.

Another point of interest in this algorithm is that the gradient of the cost function is not needed. However, it has to be noticed that the miminum found is a local one and so, the solution obtained is a particular solution that is expected to give a satisfying behavior to the model. In what follows, our attention is focused on the identification of the model in shear direction which is the part involving the most coefficients. Similar developments are used to identify the model in fiber direction for both tensile and compressive behaviors.

Experimental Data Available

In this first study and while waiting for the begining of an imminent experimental campaign, "pseudo experimental data" are generated from the model described in [START_REF] Rozycki | Contribution au développement de lois de comportement pour matÈriaux composites soumis à l'impact[END_REF]. This model includes both plasticity and static damage under the same framework as the present one. Here it is really important to notice that the "pseudo experimental data" follows the mechanical behavior of glass-fiber/epoxy-resin woven fabrics. So, in all that follows, the identification will also be made for this kind of materials. However, the method developed below is quite general and will be applied soon on new experimental data to deal with carbon-fiber/epoxy-resin woven fabrics. A fundamental point in an identification procedure is to solicit all the aspects of the model as was shown by [START_REF] Fossum | Parameter estimation for an internal variable model using nonlinear optimization and analytical response sensitiv-ities[END_REF][START_REF] Mahnken | Parameter identification for viscoplastic models based on analytical derivatives of a least-squares functional and stability investigations[END_REF]. For composite woven fabrics, this involves irreversible strains, loss of stiffness, and rate dependence. To match these requirements, the generated "pseudo experimental data" represent cyclic shear tests for strain rates evolving from 1 • 10 -3 s -1 to 1 • 10 3 s -1 . When the strain exceeds 0.1 m.m -1 , the strain rate is multiplied by a factor 10 in order to test more efficiently the strain rate effects. The stress versus strain curves obtained are plotted on Figure 1.

The "pseudo experimental data" are then computed in tables with the time and the associated total strain and stress. For each given strain rate, a table named data X P ˙ = t; t ; σ XP ˙ is generated for comparison with numerical results.

Pattern Search Method

This section deals with the Pattern Search Algorithm. Let χ be a set of parameters that is part of ϒ the space of the admissible parameters delimited by lower and upper bounds. As presented below, ϒ is composed of all the parameters involved in the shear For each set of parameters the time t and the strain t are used to drive the numerical simulations. Tables of numerical results are then computed under the form data Num ˙ = t; t ; σ Num ˙ . It has to be noticed that the identification procedure is performed at the level of an integration point. By substraction, stress error vectors are introduced and calculated by:

= σ XP -σ Num (18)
where each composant of stands for the error on one point between the experimental and numerical values of stress. A global indicator of the quality of the identification is then computed from the norm of . It takes into account all the experimental data available. In the least squares means, the indicator || || can be defined as:

γ = || || = 1 2 . ( 19
)
To summarize, the problem can be expressed as: "Find χ the set of parameters whose cost function γ is minimum with χ ∈ ϒ the space of the admissible parameters."

To solve this problem, the Pattern Search Algorithm presented in Figure 2 is employed.

Starting from a given set of parameters χ k , the algorithm computes the cost and then prospects for a better candidate χ k+1 by forming a "mesh" of points generated by the translation of χ k in each direction e i of the space of the parameters at a distance of a k i . If a better candidate appears, then a dilatation phase starts around χ k+1 in the meaning that the next increment will be larger than the previous one (||a k+1 || > ||a k ||). If there is no better candidate, the algorithm contracts its search by decreasing the distance of prospect (||a k+1 || < ||a k ||). The procedure goes on until a tolerance on ||a|| or a maximum number of iterations is reached. For more information, consult [START_REF] Lewis | Direct search methods: then and now[END_REF].

Implementation of the Method

The identification procedure is performed with the commercial software Matlab and the Genetic Algorithm Toolbox that includes a Pattern Search Algorithm. To conduct the procedure to success, some simplifications and assumptions have to be made on the set of parameters to identify. First of all, the huge vector of parameters: can be reduced with special considerations:

• G 12 the elastic mudulus and v G 12 its associated parameter that drives the strain rate dependence can be identified separately from the other parameters by considering only the linear reversible part of the experimental data;

• τ 12 that upper bounds the damage rate and α 12 can be estimated "à priori" from the litterature [START_REF] Ladevèze | A mesomodel for localization and damage computation in laminates[END_REF][START_REF] Suffis | Damage model with delay effect: Analytical and numerical studies of the evolution of the characteristic damage length[END_REF] since there is, to our knowledge, no simple experimental test that enables to give a value to this parameters.

Finally, and with the preceeding assumptions, the set of parameters to identify becomes:

χ = σ y ; K ; m; Q; β; Y o 12 ; v Y o 12 ; Y c 12 ; v Y c 12 ; w Y c 12 T
A huge campaign of identification has been performed on several sets of data. The upper and lower bounds of the parameters were determined upon litterature and physical considerations: First of all, it is visible that two curves differ from the others: ξ = 0 and ξ = 0.75. Next, for the initial sets of parameters located in the middle of the interval (ξ = 0.25 and ξ = 0.5), the final sets of parameters obtained after the optimization procedure are close to the ones from the others. Concerning the parameters of the delayed damage mesomodel, it appears that

Y o 12 , v Y o 12 , Y c 12 , v Y c
12 and w Y c 12 all lie between their lower bound and 12 percent of their upper bound. The relative independence of these results from the starting point, and the satisfying global costs, lead us to think that the adequate parameters for the delayed damage mesomodel are close to the one determined in this first step of identification. The upper bounds of these parameters are then reduced to 15 percent of their initial value in the next steps of optimization. Now, regarding the part of the model that deals with irreversible strains and strain rate dependence (viscoplasticity), the conclusions are more delicate to draw. For β and Q, and if we do not consider the curves ξ = 0.75 and ξ = 0, the results are quite homogeneous and close to the identified values obtained in [START_REF] Rozycki | Contribution au développement de lois de comportement pour matÈriaux composites soumis à l'impact[END_REF] that were used to generate the "pseudo experimental data." As β and Q drive the isotropic hardening phenomenon, the capability of the Pattern Search Algorithm to find the same values starting from very different points is quite interesting. For K , the Pattern Search Algorithm always leads nearby the lower bound. Except for ξ = 0.75, the same results are obtained for σ y which actually means that there is no elasticity, and that viscoplasticity occurs immediately. Last, m seems to be limited by the imposed upper bound.

From the preceeding considerations, three conclusions can be drawn:

• the parameters of the delayed damage mesomodel are almost identified (even if the model is coupled, various sets for the "viscoplastic part" do not seem to affect them); • the isotropic hardening parameters β and Q have been determined, and their value is similar to the one used to generate the "pseudo experimental curves"; FIG. 4. Comparison between experimental and numerical curves for both shear and fiber directions.

• efforts have to be made to finalize the optimization process on the viscoplastic part of the model.

After this first step, the bounds of the almost identified parameters are severly reduced, and the Pattern Search Algorithm is launched again with the objective to specify the values of the parameters of the viscoplastic part of the model. Several intervals of admissible parameters are generated with higher upper bound for m and lower upper bound for K .

Results

Finally, after a last step of optimization with a descent method to ensure that the minimum found is really a local minimum, the obtained results can be observed in the Figure 4 for shear direction (left side) and fiber direction (right side). The associated cost is equal to 34.

They clearly highlight the capacity of the optimization procedure to determine the parameters of the model. Moreover, the order of magnitude of the identified parameters seems to be physically reasonable (although the initial bounds were rather broad), which can give more confidence in the method. The Pattern Search Algorithm has demonstrated its ability to help with the identification of models dealing with strongly coupled physical phenomena. It has led to usable sets of parameters for the model in shear and fiber directions as summarized in Tables 1 and2 below.

The model presented in this paper allows to perform simulations on composite woven fabric structures upon six decades of strain rate.

APPLICATIONS

In the two first parts of this section, the case of the dynamic extension of a bar is studied to check the ability of the model to avoid strain localization and mesh dependence. This unidirectional example allows to differentiate the behavior of the elementary ply in fiber direction (first part) and in shear direction (second part). In the two following parts, structural calculations are presented. The third part deals with the impact of carbonfiber/epoxy-resin woven fabric plates, and the fourth is related to the dynamic crushing of composite tubes.

Tensile Simulations on a Bar-Behavior in Fiber Direction

The bar considered in this study measures 100 mm length for a cross section area of 1 mm 2 . One extremity is clamped whereas a constant velocity is imposed on to the other one. Mesh dependence is investigated through the use of three different meshes involving 100, 200 and 400 linear truss elements. The equation of motion is integrated through the time with a standard explicit central difference scheme. The material behavior is implemented in the finite element code Abaqus via a VUMAT subroutine. The material parameters employed for the model in fiber direction are given in Table 1.

Many numerical simulations have been performed, and convergence of the material behavior integration scheme was achieved in all cases with satisfying accuracy in fiew iterations. The left side of the Figure 5 presents the profile of the damage obtained after 2.5 • 10 -5 s along the bar for different meshes and for a prescribed velocity of 5 • 10 5 mm.s -1 (which corresponds to an average macroscopic strain rate of 5000 m • (m.s) -1 ). On the right side, the damage distribution along the bar for a mesh of 100 elements is ploted for various strain rates. The first observation that can be made about the preceeding curves is that the mesh objectivity of the results is ensured. Indeed, a modification of the mesh size of a factor 4 does not influence the response of the model in terms of damage distribution. Moreover, the curves obtained for various strain rates show that the fully damaged zone is not confined into a single element when the mesh is fine enough to capture the internal length associated with the delayed damage mesomodel and when the strain rate is high enough to make the damage propagate. It also has to be added that the positivity of the dissipation was verified in all cases.

In conclusion, one can say that:

• the model is mesh objective;

• spurious localization phenomenon is avoided.

Hence, the introduction of strain rate dependent moduli in the formulation of the delayed damage mesomodel has not removed the interesting properties of the initial model but has, on the contrary, enlarged its domain of validity in terms of strain rates. It can be noted that the integration step time has to be chosen with care, which means that:

• the elastic wave should not propagate to more than one element at a step time of integration t < min(l e )

Emax ρ

(with l e the characteristic lengths of the elements of the model); • the damage should not pass from 0 to 1 in one step time (which is traduced by the condition ḋmax t < 1 or written differently: δt < τ).

Finally, these considerations lead to:

t = min   min(l e ) E max ρ ; τ   (20)
and within these conditions, the Newton-Raphson procedure always converges in less than 5 or 6 iterations depending on the mesh size.

Tensile Simulations on a Bar-Behavior in Shear Direction

This part goes on with the preceeding example to study the behavior of the elementary ply submitted to shear loadings. The parameters used for the calculations are given in Table 2. The distribution of the damage along the bar is ploted in Figure 6. The influence of mesh size is presented for an imposed velocity taken equal to 5 • 10 5 mm.s -1 on the left and to 5 • 10 7 mm.s -1 on the right.

The coupling between viscoplasticity and delayed damage mesomodel disturbs neither the independence of the results from the mesh size nor the absence of localization phenomenon. Moreover, to check the influence of the step time on the response of the model, several calculations have been performed. The left curves of Figure 7 below focus on the distribution of the damage along the bar for a given loading but for different step times. As a result, it can be seen that a modification of the step time of a factor 100 does not modify the response of the model. The integration of the constitutive laws of behavior then appears to be sufficiently reliable to face more sophisticated simulations.

Eventually, and in order to show the mesh objectivity of the model on a "viscoplatic" variable, the right side of Figure 7 plots the repartition of the cumulated plastic strain among the bar. In conclusion, the uniaxial example of the extension of a bar has demonstrated the capacity of the model to give similar results whatever the meshes and the step times are. This fundamental intrinsic characteristic of the model allows to predict the behavior of structures independently from the spatial and temporal discretizations.

To go further on the analysis of the model, two structural applications are now presented.

Impact of a Composite Plate

The simulation of the impact of a composite plate is inspired from the work of Johnson et al. [START_REF] Johnson | Computational methods for predicting damage in composite structures[END_REF] performed within the German Aerospace Center (DLR). A square plate, of thickness 4.6 mm and length 300 mm, lying on a rigid square frame of length 250 mm, is impacted normaly to its surface by a sphere of mass 21 kg launched at an initial velocity of 6280 mm.s -1 . The layup sequence is [ 0 45 0 45 0 45 0 45 ] where each ply is a woven fabric with the properties described in Tables 1 and2. The Poisson's ratio ν 12 is taken equal to 0.13.

It is important to notice that the experiment was performed on a carbon-fiber/epoxy-resin woven fabric whereas the simulations are done with a glass-fiber/epxoy-resin material. As a consequence, the model presented in this paper will be qualitatively checked to show its ability to reproduce observed physical phenomena.

In the finite element model, the plate is meshed with the multilayer and reduced integration shell elements S4R available in Abaqus Explicit. These elements are 4 nodes linear elements with 6 degrees of freedom per node. Several spatial discretizations, leading to similar results (for fine enough meshes), are tested. Both the impactor and the square frame are modelled as rigid solids. The contact beween these two solids and the plate is considered to be a hard normal contact with no frictional sliding. All the degrees of freedom of the square frame are constrained FIG. 7. Damage and cumulated plastic strain distribution along the bar-influence of the step time and of the mesh. FIG. 8. Damage state in the direction 1 of the fibers after impact and evolution of the reaction force during the impact.

whereas the impactor is only free to translate normally to the glass-fiber/epoxy-resin woven fabric plate. The calculations are performed up to a time of 20 ms so that the penetration of the impactor can be observed. Figures 8 below plot the state of damage on the upper and lower sides of the plate in fiber direction 1, and the evolution of reaction force through the time for two initial velocities of 6280 and 2330 mm.s -1 (in which case, the impactor rebounds on the plate).The state of damage presented here, for a time t equal to 3 ms, corresponds to the begining of the failure of the plate as can be seen on the right curve with the abrupt decrease of the reaction force. The left figure clearly shows that the lower side of the plate is more damaged than the upper one. This numerical result is correlated with experimental observations. Moreover, as in Figures 9, the global kinematics of damage is well predicted by the model. A macroscopic crossshaped crack propagates under the impactor and finally causes the failure of the plate.

Axial Dynamic Crushing of a Composite Tube

This last structural example studies the dynamic crushing of a square composite thin-walled tube. Its interest lies in the preponderant sollicitation of the fibers in compression. Qualitative comparisons are made with the experimental results of Mamalis et al. [START_REF] Mamalis | On the response of thin-walled CFRP composite tubular components subjected to static and dynamic axial compressive loadings: Experimental[END_REF].

The section of the tube measures 100 by 100 mm and its height is also equal to 100 mm. The staking sequence is [ 0 ] 14 with layers made of carbon-fiber/epoxy-resin woven fabric. The total thickness is 3.5 mm. The thin-walled tube lies on a rigid ground and is impacted by a mass of 1 ton launched at an initial velocity of 5400 mm.s -1 . The modelling of this experiment is done in the same manner as before. Several discretizations of the tube have been checked, all leading to similar results. Ground and impactor were modelled as rigid solids and contact conditions were imposed on the tube. Figure 10 presents the macrocrack due to the failure of fibers in compression on FIG. 9. Fully damaged zone after impact-experimental carbon-fiber/epoxy-resin on the left-numerical-glass-fiber/epoxy-resin on the right. the left, and the evolution of the force applied by the tube on the impactor on the right. The calculation was stopped just after failure. The global collapse mode of the structure, a propagation of a mid-length macroscopic crack, is in agreement with the experimental observations. Moreover, the brittle failure of the thin-walled composite tube after a crushing displacement of 2 mm is in a correct order of magnitude. These results, despite the fact they were obtained from a model identified on glassfiber/epoxy-resin experimental data, are quite encouraging.

CONCLUSION

A model dedicated to carbon-fiber/epoxy-resin woven fabrics was presented in this paper. The brittle failure of the fibers under tensile loading, their nonlinear response due to the microbuckling in compression, the loss of stiffness and the irreversible strains observed in shear were all taken into account, and that, for a large range of strain rates. Moreover, special attention was paid to the construction of the model in order to avoid both the mesh dependency and the strain localization shortcomings that are often met in finite element simulations associated with softening material behaviors. After a successful step of identification based on a particular direct search method of optimization, calculations of uniaxial tensile tests on bars have shown the robustness of the model. Two structural applications have then been presented:

• the impact on a composite plate;

• and the crushing of a thin-walled tube.

In the first example, the model has predicted the kinematics of the failure of the plate with great accuracy. In the second example, the brittle rupture of the tube at mid-length was founded, in perfect agreement with the experimental observations.

To recall with the design of huge sailing ships, having a model able to reproduce various degradation phenomena is essential to ensure à priori the integrity of the structure submitted to dynamic loadings. Moreover, as numerical models of these boats involve numerous degrees of freedom, it appears really inter-esting to have reliable predictions whatever the meshes used for spatial discretization and the time steps of integration are. As encouraging results were obtained on the modelling of the glass-fiber/epoxy-resin woven fabrics, it can be said that the development of the model looks promising.

The next step of this study will be an experimental campaign on carbon-fiber/epoxy-resin woven fabrics. Other outlooks for this study deal with the modelling of the Nomex honeycomb core which connects the two composite skins in the panels used in nautical construction. First developments have been made in [START_REF] Gornet | Interactions modèles expériences sur des âmes nids d'abeilles Nomex R : Application au design d'un voilier multicoque de course océanique[END_REF] concerning the linear elastic behavior until rupture of these structural materials. Further ones, dealing with the rate effects, will come next. A delaminating interface between the skins and the core will also be investigated in order to be more accurate on the design of oceanic ships subjected to dynamic loadings.

12 ]

 12 T and, for irreversible strains, vp = [0 0 vp 12 ] T .

FIG. 1 . 12 Tϒ

 112 FIG.1. "Pseudo experimental data" for identification in shear direction.

FIG. 3 .

 3 FIG. 3. Cost and associated number of iterations obtained for various initial sets of parameters.

χFigure 3

 3 Figure3shows the final set of parameters and its associated cost versus the position ξ of the initial set of parameters χ 0 from the bounds. ξ = 0 means that χ 0 = χ lower and ξ = 1 means that χ 0 = χ upper . The calculations are performed for ξ = {0; 0.25; 0.5; 0.75; 1} with G 12 = 6396 MPa, v G 12 = 0.015 (obtained from a first step of minimization upon the linear reversible parts of the experimental curves), τ 12 = 1 • 10 -6 m 2 .(m -2 .s -1 ) and a reference strain rate˙ ref = 1 • 10 -3 m.(m.s) -1 .First of all, it is visible that two curves differ from the others: ξ = 0 and ξ = 0.75. Next, for the initial sets of parameters located in the middle of the interval (ξ = 0.25 and ξ = 0.5), the

FIG. 5 .

 5 FIG. 5. Damage along the bar for the model in fiber direction-influence of the mesh and of the strain rate.

FIG. 6 .

 6 FIG.6. Damage along the bar for the model in shear direction-influence of the mesh and of the strain rate.

FIG. 10 .

 10 FIG. 10. Simulation of the crushing of a composite tube-mid-length collapse mode.

  12

TABLE 1

 1

		Parameters for the model in fiber direction
		Elastic properties		Damage properties
	E 11 v E 11 ρ 11	30780 MPa 0.0211 1400 kg.m -3 3.64 • 10 -4	Y o 11 v Y o 11 Y c 11 v Y c 11 vv Y c 11	2.523 MPa 0.041 275.38 MPa -0.018 2.07 • 10 -4
			α 11	1
			τ 11	1 • 10 -5 m 2 • (m -2 .s -1 )

)d 12 with Y o ii the yield thermodynamic force that initates the damage and Y c ii a critical thermodynamic force.


Laurent Gornet

Steven Marguet

Gilles Marckmann

Modeling of Nomex honeycomb cores, linear and nonlinear behaviours

Keywords: Nomex R honeycomb core, homogenization, failure, nonlinear behavior

de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Composite sandwich panels are extensively used in lightweight constructions. There excellent stiffness to weight ratio make them particularly interesting for many applications such as aircraft and aerospace structures or oceanic race sailing boats. In the latter case, typical sandwich panels are made of two carbon fibers and epoxy resin skins, separated with a Nomex R honeycomb core as illustrated in Figure 1. The understanding of the mechanical behavior of the constituents up to failure is extremely important to ensure safe design of the structures.

Several damage mechanisms contribute to sandwich structure failure. Matrix cracking often appears at first in laminates. Then, the micro-cracks coalesce and initiate delaminations and large scale crack propagation. Interfaces between core and laminates debond or core crushing lead the sandwich to failure [START_REF] Gornet | Interactions modèles expériences sur des âmes nids d'abeilles Nomex R : Application au design d'un voilier multicoque de course océanique[END_REF]. In the context of finite element model, taking into account all these phenomena while keeping an admissible computational cost is not a simple matter. A classical simplifying assumption is to consider brittle linear elastic mechanical behaviors for both the skins and the core. Associated with anisotropic stress failure criteria such as Tsaï-Hill, Tsaï-Wu or Hoffman [START_REF] Tsa Ï | A general theory of strentgh for anisotropic materials[END_REF][START_REF] Hoffman | The brittle strength of orthotropic materials[END_REF], it then becomes possible to make calculations on the whole structure. As an illustration, Figure 2, presents a mesh of the multihull race sailing boat code B1 which is now well known as Orange II. The related finite element model involves multilayered thick shell elements which allow us to deal with both the skins and the core. Each layer is defined in terms of its orientation, its elastic moduli and a convenient stress failure criteria. With the assumption that no delamination occurs, a first design is then carried out. But in front of hard loadings, which can be wave impacts, slammings or violent swelling of veils, for example, the description of the damage phenomena has to be more precise. It results in more complicated models which involve more parameters can lead to increasing computational costs.

The aim of this study is to propose tools dedicated to the calculation of multihull race sailing boats. The focus is on the mechanical behavior of Nomex R honeycomb core. In the first section, strain based periodic homogenization technic is employed to determine the elastic constants that traduce the linear part of the mechanical behavior of the honeycombs. Next, starting from experimental observations, strengths of Nomex R honeycomb core are evaluated thanks to Euler's buckling theory. Comparisons between finite element analysis (performed within Cast3M, the finite element code from CEA) and data coming from manufacturers are presented. In the third section, a unidirectional model dedicated to the nonlinear compressive behavior of Nomex R honeycomb core is proposed. The identification technique is detailed and discussion is made on the results and on the extension of the model in full three dimensions.

ELASTIC PROPERTIES OF NOMEX R HONEYCOMB CORES

This section focuses on the method developped to compute the effecive elastic properties of Nomex R honeycomb cores. The motivation for this study comes from the fact that the only available data from manufacturers of Nomex R honeycomb cores [START_REF]HexWebTB Honeycomb Attributes and properties[END_REF][START_REF]Euro-composite, Nomex R cores[END_REF] are relative to the out-of-plane behavior which involves: the compressive elasticity modulus, the out-of-plane shear moduli and the corresponding strengths. These data are obviously insufficient to perform three dimensional calculations. This paper is in the continuity of previous studies about the plane behavior FIG. 2. Finite element mesh of code B1 multihull race sailing boat (Multiplast).

of Nomex R honeycomb cores [START_REF] Gellatley | The shear modulus of foil honeycomb cores[END_REF][START_REF] Gibson | Cellular solids structures and properties[END_REF] and on other ones, dealing with cellular solids [START_REF] Hohe | Effective elastic properties of triangular grid structures[END_REF][START_REF] Aminanda | Experimental analysis and modeling of the crushing of honeycomb cores[END_REF][START_REF] Grediac | A finite element study of the transverse shear in honeycomb cores[END_REF].

Strain Based Periodic Homogenization

Since modeling the true geometry of Nomex R honeycomb cores on the whole finite element model is unrealistic for reason of computational cost, another approch has to be implemented. Strain based periodic homogenization technic is a powerful tool to model the periodic structures that are honeycombs. The basic idea is to replace a periodic medium of complex geometry by a continuous medium with an equivalent mechanical behavior. The first step of the approach is to define a Representative Volume Element (RVE) which can, by series of translations, reproduce the whole discontinuous geometry [START_REF] Suquet | Plasticité et homogénéisation[END_REF]. In the case of this study, two geometries presented in Figure 3 are taken into account: hexagonal and over-expanded rectangular unit cells of honeycombs.

By convention, and in all that follows, the e 1 (strong) and e 2 (weak) axis are in plane axis whereas e 3 is the out-of-plane axis. O is the origin of basis = (O, e 1 , e 2 , e 3 ). Starting from the center of a unit cell, two orthogonal plans of symetry can be definied: (O, e 1 , e 3 ) and (O, e 2 , e 3 ). As a result, Nomex R honeycomb cores exibit an orthotropic mechanical behavior.

The link between the micro scale of the RVE and the macro scale of the equivalent continuous medium is done thanks to Hill Mandel's theorem. It postulates the equality between the mean elastic energy on the RVE and the energy in the equivalent medium, which reads:

: E = 1 V V σ(x) : (x)dV (1)
where:

• σ(x): is the stress at point (x) in the RVE;

• (x): is the strain measure at the same point;

• V : the volume of the RVE (dV = dx 1 dx 2 dx 3 is the elementary volume);

• : is the homogeneous macro stress of the equivalent continuous medium;

• E: is the homogeneous macro strain loading.

By introducing the mechanical behaviors at the micro scale (known behavior of Nomex R paper) and at the macro scale (unknown behavior of the equivalent medium), Eq. (1) deals to:

           σ(x) = K micro (x) : (x) = K micro : (x) = K MACRO : E E : K MACRO : E = 1 V V (x) : K micro : (x)dV (2)
with:

• K micro : elasticity tensor representative of the mechanical behavior of Nomex R paper; • K MACRO : elasticity tensor of the mechanical behavior of equivalent medium.

A macroscopic strain field E imposed is then imposed at the micro level of the RVE by prescribing a displacements field: u(x) = E imposed x, ∀x ∈ V . On the boundary of the RVE, periodicity conditions also have to be taken into account. If ∂V is the boundary of the RVE, it can be decomposed into ∂V σ which is the boundary with imposed stress and ∂V the boundary with imposed strain. As the study is strain based, no stress is imposed on ∂V σ , but on ∂V , periodicity conditions are introduced by writing that the field of displacements is the sum of two terms:

• the field of displacements due to the imposed macro strain field E imposed ; • a periodic field of displacement that links the opposite points on the boundary of the RVE.

If A and B are two opposite points on the boundary of the RVE, periodicity relations are:

     u(A) = E imposed x A + u periodicit y (x A ) u(B) = E imposed x B + u periodicit y (x B ) u(x A ) -u(x B ) = E imposed (x A -x B ) (3)
since u periodicit y (x A ) = u periodicit y (x B ). To summarize the previous considerations, the elementary mechanical problem to solve is now introduced within Eq. ( 4):

                   div(σ(x)) = 0 inV (equilibrium) σ(x) = K micro : (x) inV (behavior) (x) = grad S (u(x)) inV (kinematics) σ(x) = 0 on∂V σ (boundary conditions) u(x) = E imposed x on∂V (boundary conditions) +u periodicity (x) (4)
It has to be noted that, in finite element analysis, the last equation of Eq. ( 4) is introduced with linear relations between displacement degrees of freedom of the opposite nodes located on ∂V . Last step is the choice of the macroscopic strain field E imposed . By reminding us that due to symetries, Nomex R honeycomb cores present an orthotropic mechanical behavior, the elasticity tensor is:

                                                                         E = K -1

M AC RO

: or, with engineering notations:

          E 11 E 22 E 33 √ 2E 12 √ 2E 23 √ 2E 31           =                  1 E 1 -ν 12 E 1 -ν 13 E 1 0 0 0 -ν 21 E 2 1 E 2 -ν 23 E 2 0 0 0 -ν 31 E 3 -ν 32 E 3 1 E 3 0 0 0 0 0 0 1 2G 12 0 0 0 0 0 0 1 2G 23 0 0 0 0 0 0 1 2G 31                  ×            11 22 33 √ 2 12 √ 2 23 √ 2 31            (5)
Nine macroscopic loadings are then applied to the RVE in order to compute the macroscopic behavior of the equivalent continuous medium. They are presented in Figures 4 , 5 and6 for over-expanded rectangular cells.

The stored elastic energies of these finite element models are then evaluated and, by using Hill Mandell's theorem (Eq. 2), the characteristics of the equivalent continuous medium are computed.

Results, Effective Elastic Properties

The mechanical properties of numerous geometries of Nomex R cores have been derived from the method presented above. The microscopic mechanical behavior of Nomex R paper has been evaluated by rescaling a given out-of-plane shear modulus determined by periodic homogenization technique with the value given by manufacturers [START_REF]HexWebTB Honeycomb Attributes and properties[END_REF][START_REF]Euro-composite, Nomex R cores[END_REF]. By this way, with only one set of data from the manufacturer, it is possible to estimate the whole three dimensional mechanical behavior of a Nomex R honeycomb core. Table 1 presents the elastic properties for three kinds of over-expanded rectangular honeycombs of cell size 4.8 mm, sheet thickness 51 µm and densities: 29 kg.m -3 , 48 kg.m -3 and 64 kg.m -3 . Very good agreement is found with available data concerning the out-of-plane properties.

Next, the evolutions of the elastic moduli in function of the core density are plotted on Figure 7 for two types of honeycombs which are: 4.8 mm over-expanded rectangular and hexagonal cells. An increasing density of the core appears to have significant effect on the stiffness of the structure.

With this first study, relying on the strain based periodic homogenization technique, it becomes possible to determine the full three dimensional mechanical behavior of any Nomex R honeycomb core, any in the sense of any goemetry, any cell size, any density of Nomex R sheet and any thickness. This aspect is particularly interesting since available catalogues only provide the out-of-plane moduli of these materials for a given thickness of 12.7 mm.

A pragmatic approach to evaluating the strengths of honeycomb cores is presented in next section.

STRENGTHS OF NOMEX

R HONEYCOMB CORES

Experimental Observations

In the approach presented here, the determination of the strengths of Nomex R honeycomb cores relies on the assumption that the buckling phenomenon drives the failure of the structural materials. To justify this strong hypothesis, several tests are performed on out-of-plane compressive and shear specimens which are presented on Figure 8 The displacement is measured by the magnetic captor that drives the test, and the force is obtained with a 25 kN cell force. In the case of compressive tests the imposed velocity is 0.04 mm.s -1 equivalent to a strain rate of 0.001 s -1 . For simple shear tests, the imposed velocity is 0.16 mm.s -1 which also leads to an average strain rate of 0.001 s -1 .

Figure [START_REF] Aminanda | Experimental analysis and modeling of the crushing of honeycomb cores[END_REF] plots the crushing stress in function of crushing strain response of two over-expanded rectangular Nomex R honeycomb cores.

The photograph in Figure 9 has been taken just after the appearance of a global mode of instability. It shows the formation of a macroscopic buckling ply on the whole block of

Nomex

R core and is directly related to the brutal loss of stiffness visible on the response stress curve. The link between the failure of Nomex R honeycomb core and the buckling of its walls is clearly demonstrated. Figure 10 does the same as Figure 9 but for simple shear specimens. W (or e 1 ) and L (or e 2 ) stress versus strain curves are plotted on the graph. After a slightly nonlinear reversible response, a subit loss of stiffness is observed. For W direction, this decrease in stress appears in a range of strain evolving from approximatly 0.01 and 0.012 m.m -1 . The photograph on the upper left draws the lower right block of Nomex R just before the brittle failure (presented on the other photograph). At this stage, no macro crack is visible and moreover, a global bending mode appears on the whole block. This leads us to think that the loss of stiffness is not related to a brittle failure phenomenon but to the apparition of an instability which is the first buckling mode under simple shear loadings. The conclusion is that, for compressive loadings, in simple shear the failure of the structural material is due to the buckling phenomenon that is prior to the brittle failure of the block.

In the next part, and starting from this remark, a method based on Euler's buckling analysis is presented. Critical loads are determined for each kind of loading, and thanks to the homogenized behavior, the failure stress surface is constructed.

Euler's Buckling Analysis

The objective of this section is to evaluate the ultimate stress that a Nomex R honeycomb core can carry without failure. Euler's linear buckling theory is used. The basic idea is to impose an elementary unit loading on the RVE, and to evaluate the multipliers of this applied loading which result in incompatible stiffness matrix of the finite element model [START_REF] Bathe | Finite element procedures[END_REF]. For the unloaded structure, the stiffness matrix can be computed:

S 0 = V 0 B T (x) : K micro : B(x)dV (6)
where:

• S 0 : is the stiffness matrix of the unloaded structure obtained by integration on the undeformed structure; • B(x): is the gradient of the interpolation matrix.

After application of an elementary unit loading (by imposing a marcoscopic strain field E unit corresponding to out-of-plane compressive or simple shear loading for example), the stiffness matrix can be evaluated on the deformed shape:

S d = V d B T (x) : K micro : B(x)dV ( 7
)
The idea is to find a multiplier of the prescribed loading that results in a stiffness matrix S τ that has nul eigenvalue which traduces an instability:

S τ u τ = αu τ
where u τ = 0 [START_REF] Hohe | Effective elastic properties of triangular grid structures[END_REF] with:

• u τ : the matrix of the eigenvectors of S τ ;

• α: the diagonal matrix of the eigenvalues of S τ .

The main assumption in Euler's linear buckling analysis is to look for an incompatible stiffness matrix S τ under the form:

S τ = (I -λ)S 0 + (λ)S d (9)
where:

• I : is the identity matrix;

• λ: is the diagonal matrix of the load multipliers.

By introducing Eq. ( 9) into Eq. ( 8) , and by solving for nul eigenvalues: S τ u τ = αu τ = 0 u τ with u τ = 0, it can be obtained:

         [S 0 + λ -1 (I -λ)S d ]u τ = 0 whereu τ = 0 ⇒ [S 0 + γ 2 S d ]u τ = 0 withγ 2 = λ -1 (I -λ) ⇔ det S 0 + γ 2 i S d = 0 withγ i corresponding to the ithmode (10)
Equation ( 10) can then be solved for γ and next for λ the diagonal matrix of the eigenvalues which immediately gives the load multipliers that lead to instability, i.e., buckling. Once λ determined, the shape of the eigenvectors u τ can be evaluated with Eq. ( 8) and, most importanting, the macroscopic load applied on the RVE that leads to buckling. In other words, if the RVE is loaded with a macroscopic unit strain field E unit (corresponding to a pure stress loading unit on the RVE), the ith loading multiplier λ i associated with the ith buckling mode enables the calculation of the critical strain field: E c,i = λ i E unit . Thanks to the mechanical behavior of the equivalent medium obtained by homogenization technique, it is possible to determine the ultimate stress c,i corresponding to the ith mode of buckling:

c,i = K M AC RO : E c,i (11)
where:

• c,i has only one non nul value c,i in adequation with the imposed strain field E c,i .

The failure stress limits are evaluated on the first buckling mode (the one of lowest energy).

Results, Stress Failure Criterion

To determine the full three dimensional stress failure limits of Nomex R honeycomb cores, the method presented below has been applied with the following loadings:

• E unit such as E i j = 1 if i = j, 0 elsewhere, with (i; j) = {1; 2; 3} 2 (pure shear loadings);

• E unit such as E unit = K -1

M AC RO

: unit with i j = 1 if i = j, 0 elsewhere and (i; j) = {1; 2; 3} 2 (equivalent to uniaxial pure compressive tests in direction e i= j ).

The deformed shapes of the RVE submitted to out-of-plane shearing are plotted in Figure 11. A comparison between computed ultimate stress values and data available from catalogue of manufacturer Euro-composite [START_REF]Euro-composite, Nomex R cores[END_REF] is given in Table 2 for 4.8 mm cell size, 29 kg.m -3 density and 51 µm of thickness paper honeycombs.

Satisfactory correlation is founded on the out-of-plane ultimate shear stress values. A more important difference appears in compressive ultimated stress, difference which might be due to the boundary conditions on the skins of the RVE that can't prevent rotations of nodes and results in a lowest energy buckling mode (under evaluation of ultimate stress).

NONLINEAR COMPRESSIVE BEHAVIOR

As mentioned in introduction, hard loadings often occur during trans-oceanic races, inducing numerous abandons due to structural failure. To safely design the boats in front of these loadings, it appears necessary to understand the physics of the failure process and to be able to reproduce them. It means understanding how skins and core react not only in their linear domain but also in their nonlinear part. For example, when a wave impacts the hull of a boat, the energy is transmitted to the inside skin honeycomb core (or throught the foam, depending on the localization of the impact). Predicting the response of the core under out-of-plane loadings is then essential to make the calculation of the sandwich panel. The design of shock absorbers is another application that make intensive use of nonlinear models to enable the evaluation of the absorbed energy during the impact.

Compressive Tests on Sandwich Specimens

In this section, a campaign of compressive tests on sandwiches made of Nomex R honeycomb core and carbon fibers and epoxy resin composite skins is presented. The mechanical behavior of two kinds of Nomex Two TR-50 carbon fibers and epoxy resin composite skins of thickness 0.68 mm and stacked at [45-45] are stacked on each side of the core to complete the specimen. The tests are carried out on the same tensile testing machine as presented before. The specimens are placed between two plates of compression and crushed to a maximum strain of: 0.085 m.m -1 with an imposed velocity of 0.04 mm.s -1 for the first type of honeycombs and 0.08 m.m -1 with an imposed velocity of 0.025 mm.s -1 for the second type of honeycombs. This corresponds to a compressive quasi-static strain rate of 0.001 s -1 in both cases. Moreover, cyclic loadings and unloadings are imposed to the first type of specimens.

The results are summarized below for compressive tests on hexagonal honeycombs. Figure 12 presents the response of the hexagonal Nomex R honeycomb cores. The fully crushed specimen is also visible on the graph. For compressive tests on over-expanded rectangular honeycombs, results are available in Figure 9.

Four major crushing steps appear from the stress versus strain curves:

• First step: the linear growth of the stress versus strainit corresponds to the reversible elastic response of the Nomex R honeycomb cores and occurs at low level of strain (less than 0.05 m.m -1 ); • Second step: the brutal decrease of stress-this phenomenon is linked to the apparition of a huge macroscopic buckling ply orthogonal to the direction of the load as clearly visible on Figure 9, here, failure stress is arround 3.6 MPa; • Third step: the slightly rising plateau-it is due to the successive formation of folds on the walls;

• fourth step: the increasing stiffness of the specimen response-it traduces the fact that there is no longer formation of buckling ply, only densification occurs at this point.

The last observation coming from the stress versus strain curves is that the tests are very repetitive which give more confidence in the data. In what follows a model dedicated to the crushing behavior of Nomex R honeycomb cores is presented. This model tries to represent the response of the strucutral material that is Nomex R honeycomb core in the sense of an homogenized macroscopic behavior.

Modeling of the Compressive Behavior

The unidirectional model presented here relies on the thermodynamics of irreversible process [START_REF] Lemaître | Mechanics of solid materials[END_REF]. The Helmholtz's volumic free energy potential is chosen of the following form:

2ψ ( e ; d; p) = e E (1 -d) e if d < d c 2ψ ( e ; d; p) = e E e + 2h ( p) elsewhere (12)
where the following variables are respectively:

• ψ: the Helmholtz's volumic free energy [J.m -3 ];

• e : the recoverable part of the total strain, i.e., elastic strain [m.m -1 ]; • E: the elasticity modulus of the undamaged Nomex The state laws can now be derived from the Helmholtz's volumic free energy leading to the expressions of σ, the Cauchy's stress, Y the thermodynamic force associated to d and R the stress linked to nonlinear isotropic hardening. If d < d c then :

     σ = ∂ψ ∂ e = E(1 -d) e [Pa] Y = - ∂ψ ∂d = 1 2 e E e [J.m -3 ] (13) Else, if d ≥ d c then :        σ = ∂ψ ∂ e = E( t -p ) [ P a ] R = ∂ψ ∂ p = ∂h( p) ∂ p = (Q 1 + Q 2 p -p y β + ) [Pa] (14)
with:

• t : the total strain [m.m -1 ];

• p : the plastic strain which stands for the inelastic strains observed in experiments [m.m -1 ]; • p y : the cumulated plastic strains at the end of the plateau regime [m.m -1 ];

• Q 1 , Q 2 and β material constants ([Pa] for Q 1 and Q 2 ,
no unit for β).

Once the state laws are obtained, a dissipation potential has to be defined to determine the evolution laws. This potential is related to all the dissipative phenomena, which means, in the case of this study, to damage and to plasticity with nonlinear isotropic hardening:

φ (σ; Y ; R; (d)) = φ p (σ; Y ; R; (d)) + φ d (σ; Y ; R; (d)) [J.m -3 ] (15)
where:

• φ p (σ; Y ; R; (d)) = f (σ; Y ; R; (d)) = | σ 1-d |-σ y -R: is the plastic dissipation potential; • φ d (σ; Y ; R; (d))
: is the damage dissipation potential; with:

• σ y : the plastic yield stress [Pa].

Here it can be noted that the model is constructed under the assumption of associated plasticity since the dissipation potential is taken to equal the plastic flow surface f which delimits the recorvable states from the irreversible ones. The coupling between plasticity and damage is made as proposed by P. Ladevèze and E. Le Dantec [START_REF] Ladevèze | Damage modelling of the elementary ply for laminated composites[END_REF] with the introduction of effective stress | σ 1-d | which is the stress that might be applied on an undamaged specimen to obtain the same state of strain as that obtained by imposing σ on the damaged specimen. In what concerns the damage dissipation potential, it is supposed to exist, to be convex in the space of stress variables, to be positive and to take a zero value at its origin (φ d (0; 0; 0; 0) = 0).

By derivation, the evolution laws are then obtained:

           ḋ = 1 τ 1 -e -α √ Y - √ Yo √ Yc -d + if d < d c , 0 elsewhere ˙ p = λp ∂ f ∂σ if d ≥ d c , 0 elsewhere ṗ = -λp ∂ f ∂R if d ≥ d c , 0 elsewhere (16)
with:

• λp : the plastic multiplier that must check the consistency condition λp ḟ = 0:

• f < 0: elastic domain;

• f = 0 and λp > 0: plastic loading;

• f = 0 and λp = 0: plastic unloading; • Y o : the threshold that drives the initiation of damage [J.m -3 ]; • Y c : a critical energy release rate [J.m -3 ];

• α and τ : material parameters (no unit for α, [s.m 2 .m -2 ] for τ ).

The interest of using a delayed damage model has been highlighted in many works such as those of P. Ladevèze [START_REF] Ladevèze | A mesomodel for localization and damage computation in laminates[END_REF] and A. Suffis [START_REF] Suffis | Damage model with delay effect: Analytical and numerical studies of the evolution of the characteristic damage length[END_REF]. It introduces an internal length that results in the apparition of an upper bound for the damage rate ḋ which is equal to 1 τ . As a result, the localization phenonenon of damage into a single element (or spatial discretization dependence of the results in finite element analysis) is avoided.

Equation [START_REF] Simo | Computational Inelasticity," a-b[END_REF] gives the expression of the plastic multiplier:

λp = f ∂ σ f E∂ σ f + ∂ R f ∂ p R∂ R f (17)
The last point is the distinction between tensile and compressive behavior. It is done by a test on the sign of the stress at time t n . If σ (t n ) is positive, the tensile behavior (brittle elastic on stress failure criterion [START_REF] Tsa Ï | A general theory of strentgh for anisotropic materials[END_REF]) is integrated, if not, the compressive behavior presented above is used. Figure 13 illustrates the different steps.

The delayed damage model associated with linear elasticity is used until the damage variable d has reached d c a critical damage level. Then, at this stage (d ≥ d c ), the following variables are updated:

• the damage variable is locked: d = d c ;

• the yield stress is evaluated: σ y = E (1dc) e (with e = t at this point);

• inelastic strains are computed after elastic unloading: p = d c t ;

• and elastic strains are updated: e = (1d c ) t .

The stress is then computed: σ = E e = E (1d c ) t with the mechanical behavior relation (it is equal to σ y at the mo-ment when instability appears). Next, a usual return mapping plasticity algorithm [START_REF] Simo | Computational Inelasticity," a-b[END_REF] is used to update variables. Figure 13 illustrates the two major steps of the model.

Identification and Results

To determine the values of the material parameters involved in the model, an optimization procedure coupled with classical empirical rules is developed. Hence, the elasticity modulus E is measured as usual on the compressive stress versus compressive strain curve as the director coefficient in the elastic step of the response. Then the critical damage d c that initiates the plasticity regime, comes from the apparent elasticity modulus at the starting point of the regime plateau (see Figure 13). The yield cumulated plastic strain p y is defined as the inelastic strain that appears after elastic unloading just before the densification step.

For the other coefficients of the model:

• α, τ , Y o and Y c : for the damage model;

• Q 1 , Q 2 and β: for the nonlinear isotropic behavior;

an optimization procedure based on a direct search method: the pattern search algorithm, is employed [START_REF] Lewis | Direct search methods: then and now[END_REF]. The only restriction is the choice of the value of τ (which is still an open question) that is the upper bound of the damage rate. Its value is fixed at τ = 10 -5 s -1 as usualy found in literature [START_REF] Suffis | Damage model with delay effect: Analytical and numerical studies of the evolution of the characteristic damage length[END_REF]. The pattern search algorithm is presented below:

1. initialize: k = 0; a direction of search: η k in the space of parameters; a range of search: r k ; the initial guess of parameters: k ; 2. starting from the initial set of parameters 0 = [α; Y o ; Y c ; Q 1 ; Q 2 ; β] 0 , compute its cost function ϒ 0 (defined later) to estimate the quality of this set; 3. next, make a translation in direction η k of range r k in the space of the parameters to find a new candidate k+1 and evaluate its associated cost function ϒ k+1 ;

4. make a test: a) if ϒ k+1 ≤ ϒ k then, the direction of search is a promizing direction, and so increase the range of search r k+1 > r k without changing the direction of search η k+1 = η k ; go to step 3; b) if ϒ k+1 > ϒ k then, the direction of search is not a suitable direction, and so change η k+1 = η k (with an arbitrary empirical rule to test in all the directions of the parameters space); go to step 3 without changing the range of change r k+1 = r k until a better candidate is has been found; if such a candidate can't be found, decrease the range of search r k+1 < r k and go to step 3.

This algorithm, quite simple, has several interesting advantages:

• its ability to deal with numerous parameters;

• its robustness (if an incompatible set of parameters is met leading to an unreal cost function for example, the algorithm changes its direction of search); • its capacity to check a wide area in the space of parameters.

However, it has some non negligible drawbacks such as a starting point dependency and a tendency to converge on a local minimum. In this study, the drawbacks are not a limiting point since even if the solution does not correspond to the global minimum, it has no real effect on the ability to perform calculations with the model and on the capacity to match experimental data. Many starting points were checked (initial guess 0 ), leading for an important part of them, to similar results in terms of set of parameters. To finish with the presentation of the method, the cost function is defined as an error in the sense of least squares between experimental data and computed values:

ϒ ( k ) = i data i -computedvalues i ( k ) 2 ( 18
)
at each point i of the curves. After the step of optimization all the parameters values are available. They are summarized in Table 3 for the over expanded rectangular honeycombs: With these parameters, the crushing stress versus crushing strain curve computed is compared to experimental data in Figure 14: From a general point of view, the computed values are very similar to the experimental data that were used to make the identification. The four steps of the failure of Nomex R honeycomb cores: elasticity, apparition of the instability, plateau and densification are reproduced with a satisfactory accuracy. However, at low level of strain, just after the apparition of the instability the unloading is not well represented, it is stiffer with the model than in experiment. This can be explained by the instantaneous change of mechanical behavior in the model, change of mechanical behavior which might be softer for real. Another small drawback of the model appears at very high level of strains. Experiments exibit an unloading stiffer than the initial elasticity modulus. This is due to the honeycomb sheets which are under compressive loading between the two skins of the sandwich. The model is not able to deal with this phenomenon.

As a conclusion for this first unidirectional model dedicated to the nonlinear behavior of Nomex R honeycomb cores, it can be said that the obtained results are quite promising. The next step is to understand the coupling between out-of-plane compression and shear in order to propose a full three dimensional nonlinear behavior for Nomex R honeycomb cores.

CONCLUSION

In this work, the first point was to develop software dedicated to the determination of the effective elastic properties of any Nomex R honeycomb core (any geometry of cell, any cell cell size, any density, any thickness). This tool, relying on the strain based periodic homogenization technique, gave satisfactory results compared with manufacturer's data. The next point was the determination of the strengths of Nomex R honeycomb core. In agreement with experimental observations, the linear Euler's buckling analysis was introduced to evaluate the critical macroscopic strain field that leads to failure. Thanks to the homogenized mechanical behavior, the ultimate stress was computed, and, once again, good results have been obtained from data coming from manufacturers. The last point was a prospects on the nonlinear behavior of Nomex R honeycomb cores submitted to out-of-plane compressive loadings. A model based on the thermodynamics of irreversible effects, inclucing damage and plasticity was developped, and an identification technique proposed. Very good adequation was obtained between experimental data and computed values.

FIG. 1 .

 1 FIG. 1. Composite sandwich panel used for nautic construction: carbon fibres and epoxy matrix laminated skins separated by a Nomex R honeycomb core.

FIG. 3 .

 3 FIG. 3. Hexagonal and over-expanded unit cells.

FIG. 4 .

 4 FIG. 4. Over-expanded rectangular unit cell submitted to simple stretches: E 11 , E 22 , and E 33 .

  FIG. 5. Over-expanded rectangular unit cell submitted to simple shear: E 12 , E 23 , and E 31 .

FIG. 6 .

 6 FIG. 6. Over-expanded rectangular unit cell submitted to biaxial stretches: E 11 & E 22 , E 22 & E 33 and E 33 & E 11 .

FIG. 10 .

 10 FIG. 10. Over-expanded rectangularhoneycomb cores submitted to simple shear tests.

FIG. 11 .

 11 FIG. 11. Deformed shapes of over-expanded rectangular RVE submitted to simple out-of-plane shear

  FIG.12. Hexagonal honeycomb core submitted to compressive test.

R

  honeycomb core [Pa]; • d: a damage variable [m 2 .m -2 ] (damage area for initial undamage unit area); • d c : a critical damage that traduces the initiation of the plateau regime; • h ( p): a function of the an accumulated plastic strain p associated with a nonlinear isotropic hardening [J.m -3 ].

FIG. 13 .

 13 FIG. 13. Compressive response of NomexR honeycomb cores.

FIG. 14 .

 14 FIG. 14. Comparison between computed and experimental crushing curves.

TABLE 1

 1

			Computed elastic properties of Nomex	R honeycomb cores
		ECA-R 4.8 mm 29 kg.m -3 51 µm		ECA-R 4.8 mm 48 kg.m -3 51 µm
		FE analysis	Euro-composite		FE analysis	Euro-composite
		Minimal	Typical	Minimal	Typical	Minimal	Typical	Minimal	Typical
	ν 12	0.2606				0.2606
	ν 23	0.0219				0.0219
	ν 31	0.224					0.224
	E 1	44.5	72.5			101.1	128.7
	E 2	4.4	7.1				9.9	12.6
	E 3	79.4	129.5			180.5	229.8
	G 12	0.5	0.9				1.2	1.6
	G 23	8.6	14.1	9	14	19.6	25	18	24
	G 31	14.6	23.9	14	24	33.3	42.4	36	44

FIG. 7. Evolution of elastic properties with the density. FIG. 8. Compressive and simple shear specimens. FIG. 9. Over-expanded rectangularhoneycomb cores submitted to compressive tests.

TABLE 2

 2

	Ultimate stress of over-expanded rectangular
			honeycomb cores.
	ECA-R	FE (min)	Tests (min)	FE (typ)	Tests (typ)
	σ c 11 σ c 22 σ c 33 σ c 12 σ c 23 σ c 31	0.185 0.282 0.39 0.72 0.25 0.33	--0.6 -0.31 0.32	0.31 0.48 0.62 1.2 0.415 0.55	--0.85 -0.42 0.44
	throught the outside skin first, and then throught the Nomex	R

TABLE 3

 3 Parameters values for the model.

	E	128 MPa	Q 1 2.42 J.mm -3	α	10
	d c 0.87 m 2 •m -2 Q 2 450 J.mm -3	τ	10 -5 s -1
			β	0.54	Y o 0.06 J.mm -3
			p y	0.53 m.m -1	Y c 0.01 J.mm -3


Thomas Sauvaget

A FORMULA FOR ζ(2n + 1) AND SOME RELATED EXPRESSIONS

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Are all values of the Riemann zeta function irrational numbers when the argument is a positive integer ? This question goes back to the XVIIIth century when Euler published in 1755, n being a positive integer, that ζ(2n) = (-1) n+1 B2n(2π) 2n

2(2n)!

(where B 2n ∈ Q is an even Bernoulli number) and Lindemann proved in 1882 that π is transcendental (hence none of its powers is rational) [START_REF] Lindemann | Über die Zhal π[END_REF] 1 .

On the other hand, only in 1978 did Apéry [START_REF] Apéry | Irrationalité de ζ(2) et ζ(3)[END_REF] famously proved that ζ(3) is irrational. This was later reproved in a variety of ways by several authors, in particular Beukers [START_REF] Beukers | A note on the irrationality of ζ(2) and ζ(3)[END_REF] who devised a simple approach involving certain intergrals over [0, 1] 3 (which will be recalled in section 3). The reader should consult Fischler's very informative Bourbaki Seminar [START_REF] Fischler | Irrationalité de valeurs de zêta [d'après Apéry, Rivoal[END_REF] for more details and references. In the early 2000s, an important work of Rivoal [START_REF]La fonction zêta de Riemann prend une infinité de valeurs irrationnelles aux entiers impairs[END_REF] and Ball and Rivoal [START_REF] Ball | Irrationalité d'une infinité de valeurs de la fonction zêta aux entiers impairs[END_REF] determined that infinitely many values of ζ at odd integers are irrational, and the work of Zudilin [START_REF] Zudilin | One of the numbers ζ(5), ζ(7), ζ(9), ζ([END_REF] proved that at least one among ζ(5), ζ [START_REF] Brown | Irrationality proofs for zeta values, moduli spaces and dinner parties[END_REF], ζ [START_REF] Dupont | Odd zeta motive and linear forms in odd zeta values[END_REF] and ζ [START_REF] Gourdon | Irrationality proofs[END_REF] is irrational. Despite these advances, to this day no value of ζ(2n + 1) with 2n + 1 > 3 is known to be irrational.

One dimensional integral formulas for ζ(2n + 1) have been known for a long time, for instance the 1965 monograph of Abramowitz and Stegun [START_REF] Abramowitz | Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables[END_REF] gives:

ζ(2n + 1) = (-1) n+1 (2π) 2n+1 2(2n + 1)! 1 0 B 2n+1 (x) cot(πx)dx
While it bears a striking structural analogy with Euler's formula for ζ(2n), it is not obvious how one might try to prove or disprove that these numbers are irrational.

On the other hand, multidimensional integral formulas for ζ(2n + 1) are more recent: as mentionned by Baumard in his PhD Thesis [START_REF] Baumard | Aspects modulaires et elliptiques des relations entre multizêtas[END_REF], quoting Zagier [START_REF] Zagier | Values of Zeta Functions and their Applications[END_REF], it is Kontsevich in the early 1990s who found such a type of formula for Multiple Zeta Values, which in the case of simple zeta boils down, for any odd or even k, to:

ζ(k) = 1 0 dx 1 x 1 x1 0 dx 2 x 2 • • • x k-2 0 dx k-1 x k-1 x k-1 0 dx k 1 -x k
An erroneous claim of a proof of irrationality of all ζ(2n + 1) in the previous version of this manuscript (arXiv.org/03174v3) has been withdrawn. The author apologizes for that claim and wishes to thank the editorial board of PMB for pointing out the error and rejecting the paper. What remains of this work is not intented for publication anymore. The author is very grateful to the anonymous referee of an earlier version of this work (arXiv.org/03174v2) for (a) many comments that improved greatly the readability of the paper, (b) a technical suggestion which allowed the author to prove Theorem 2.5 (which had been stated earlier as a conjecture based on numerical observations), (c) inviting the author to resubmit a revised version. The author also would like to thank the numerous contributors to useful freely available online knowledge resources, in particular the arXiv, Wikipedia, the SagemathCloud, WolframAlpha, and Stack Exchange sites.

1 Lambert proved in 1761 that π is irrational [START_REF] Gourdon | Irrationality proofs[END_REF], but this does not imply that all powers of π are so.

This is easily proved by expanding the integrand in geometric series and integrating. This can be rewritten more simply as a multidimensional integral over a unit hypercube:

ζ(k) = [0;1] k dx 1 • • • dx k 1 -x 1 • • • x k
While this is much closer to the type of integrals that Beukers used, it is not clear how it might be adapted directly to prove that zeta is irrational at odd integers.

One should mention that Brown [START_REF] Brown | Irrationality proofs for zeta values, moduli spaces and dinner parties[END_REF] has in the past few years outlined a geometric approach to understand the structures involved in Beukers's proof of the irrationality ζ(3) and how this may generalize to other zeta values, see also the recent work of Dupont [START_REF] Dupont | Odd zeta motive and linear forms in odd zeta values[END_REF] on that topic.

In this work, we go along another path and prove polylogarithmic identities which then allow to write each ζ(2n + 1) as an alternating sign (-1) n+1 times a multiple integral over a n-dimensional unit hypercube of certain functions involving logarithms (rather than unsigned integrals over a (2n + 1)-dimensional unit hypercube as in the previously mentioned formulas). These functions are shown to have an interesting property: raising some of the variables to a power leads to a fractional multiple of ζ(2n + 1) that belongs to the interval ]0, ζ(2n+1)[. We also investigate two related families of integrals for n = 2, for which we establish closed-form expressions, but show that when used in Beukers's framework they fail to produce the irrationality of ζ [START_REF] Baumard | Aspects modulaires et elliptiques des relations entre multizêtas[END_REF]. Finally, we study another type of functions with better decay properties in a 1-dimensional integral related to the previous formulas.

It is rather curious that these precise identities and integrals seem not to have been considered before, despite their simplicity. A search through the literature did not return them (we have used the treatise of Lewin [START_REF] Lewin | Stuctural Properties of Polylogarithms[END_REF] as well as the relevant page on functions.wolfram.com [15]). Identities involving polylogarithms of different degrees are rather scarce, all the more so when all variables must be integers, and representations of ζ(s) as multiple integrals over bounded domains, including some that have been worked out very recently by Alzer and Sondow [START_REF] Alzer | A parameterised series representation for Apéry's constant ζ(3)[END_REF], only go as far as a double integral. The idea to consider the formulas presented below came to the author in a fortunate way after studying and trying to generalize an integral representation of

ζ(3) = 1 2 1 0 log(x) log(1-x) x(1-x)
dx established by Janous [START_REF] Janous | Around Apéry's constant[END_REF] (and mentioned by Alzer and Sondow, where the author first learned about it), while the idea of trying to prove the irrationality of ζ(5) was a reaction to a footnote in a section of the fine undergraduate book of Colmez [START_REF] Colmez | Éléments d'analyse et d'algèbre (2e édition)[END_REF] devoted to Nesterenko's proof of the irrationality of ζ(3).

Values of ζ at odd integers as multidimentional integrals on unit hypercubes

Recall that the polylogarithm function of order s ≥ 1 is defined for z

∈ {z ∈ C, |z| < 1} by Li s (z) := +∞ k=1 z k
k s (and is extented by analytic continuation to the whole complex plane).

The aim of this section is to establish the following results (which we could not locate in the literature).

Theorem 2.1. Let n be a positive integer. Then the value of the Riemann zeta function at odd integers is:

ζ(2n + 1) = (-1) n+1 [0;1] n n i=1 log(x i ) x i log 1 - n i=1 x i dx 1 • • • dx n Proof.
Let n be a positive integer, and for any integer 1 ≤ k ≤ n define D k,n to be the set of all ordered k-tuples j

1 < • • • < j k of distinct integers taken in {1, . . . , n}. So #D k,n = n k .
Define for any (x 1 , . . . ,

x n ) ∈]0, 1[ n (the open unit hypercube of dimension n) the function M n as M n (x 1 , . . . , x n ) := (-1) n+1 Li 2n+1 ( n u=1 x u ) + n i=1 (-1) i   J∈Dn-i,n j∈J log(x j )   Li n+i ( n u=1 x u )
For any positive integer n, the function M n is at least C n when x i = 0 for all i = 1, . . . , n. We are going to show that:

∂ n ∂x 1 ∂x 2 • • • ∂x n M n (x 1 , . . . , x n ) = n i=1 log(x i ) x i log 1 - n i=1
x i

Differentiating M n with respect to x 1 one finds :

∂M n ∂x 1 (x 1 , . . . , x n ) = (-1) n+1 Li 2n ( n k=1 x k ) x 1 + n i=1 (-1) i J∈D * (1) n-i,n 1 x 1 j∈J log(x j ) Li n+i ( n k=1 x k ) + J∈Dn-i,n j∈J log(x j ) Li n+i-1 ( n k=1 x k ) x 1
where D * (1) n-i,n denotes elements of the set D n-i,n where 1 is not in the (n -i)-uplet, so this is also exactly the set of ordered (n -i -1)-uplets of distinct elements taken in the set {2, . . . , n}, and so we have the inclusion D * ( 1)

n-i,n ⊂ D n-i-1,n .
Hence the telescopic cancellations witnessed in the case n = 4 occur between terms of two consecutive values of i. The remaining differentiations with respect the other variables ultimately lead to the desired expression. Thus we have an explicit antiderivative of

n i=1 log(xi) xi log (1 - n i=1
x i ) and it is clear from its expression that the generalized integral exists. Evaluating M n at (x 1 , . . . , x n ) = (1, . . . , 1) finishes the proof.

Remark 2.2. As remarked by the referee of a first version of this work:

(a) the formula in fact readily follows by using the entire series expansion -log (1

- n i=1 x i ) = +∞ k=1 1 k ( n i=1 x i ) k
(valid here as n i=1 x i ∈ (0, 1)) and a repeated use of Fubini's theorem with the integral of monomials

1 0 x k-1 i log(x i )dx i = -1
k 2 (though this would not provide the closed-form expression of the antiderivative) ;

(b) using the case n = 1 and the change of variable t = 1 -x one obtains

ζ(3) = 1 0 log(x) log(1 -x) x dx = 1 0 log(x) log(1 -x) 1 -x dx
and the formula of Janous stated in the introduction follows by adding those two integrals.

Corollary 2.3. Let n and r be positive integers. Then we have:

ζ(2n + 1) r 2n = (-1) n+1 [0;1] n n i=1 log(x i ) x i log 1 - n i=1 x i r dx 1 • • • dx n
Proof. This follows from a change of variable x → x r in the previous theorem 2 .

Remark 2.4. By summing over r ∈ N * this last formula we obtain:

ζ(2n)ζ(2n + 1) = [0;1] n n i=1 log(x i ) x i log φ n i=1 x i dx 1 • • • dx n
where φ is Euler's function defined for q ∈ [0, 1] by φ(q) := +∞ n=1 (1 -q n ).

2 This too was observed by the referee of the previous version of this work, the author had tediously proposed a proof along the lines of the theorem.

Theorem 2.5. Let n ≥ 2 and k ≥ 1 be integers. Then we have the following closed-form expressions: (i)

[0;1] 2 log(x) x log(y) y log(1 -(xy) 1 ) log(x) log(y)(xy) 2k+1 dxdy = 4ζ(6) (2k + 1) + 4ζ(5) (2k + 1) 2 + 4ζ(4) (2k + 1) 3 + 4ζ(3) (2k + 1) 4 + 4ζ(2) (2k + 1) 5 - 4 (2k + 1) 6 2k+1 n=1 1 n -4 6 j=1 1 (2k + 1) j 2k+1 i=1 1 i 7-j (ii) [0;1] 2 log(x) x log(y) y log(1 -(xy) 2 ) log(x) log(y)(xy) 2k+1 dxdy = 63ζ(6) 8(2k + 1) + 31ζ(5) 4(2k + 1) 2 + 15ζ(4) 2(2k + 1) 3 + 7ζ(3) (2k + 1) 4 + 6ζ(2) (2k + 1) 5 + 8 log(2) (2k + 1) 6 - 4 (2k + 1) 6 k i=1 1 2i + 1 -8 6 j=1 1 (2k + 1) j k i=0 1 (2i + 1) 7-j
Proof. (i) using the entire series expansion of log(1 -xy) and Fubini's theorem we have, twice integrating by parts: and the result follows by adding and substracting the first terms in each sum to make the zeta values appear.

(ii) the begining proceeds in a similar fashion to give:

[0;1] 2 log(x) x log(y) y log(1-(xy) 2 ) log(x) log(y)(xy) 2k+1 dxdy = +∞ u=1 -1 u 1 0 (log(x)) 2 x 2k+1+2u dx 1 0 (log(y)) 2 y 2k+1+2u dy = - +∞ u=1 4 u 1 (2k + 1 + 2u) 6 = - 4 (2k + 1) 6 +∞ u=1 1 2u - 2 2k + 1 + 2u + 4 (2k + 1) 5 +∞ u=1 2 (2k + 1 + 2u) 2 + 4 (2k + 1) 4 +∞ u=1 2 (2k + 1 + 2u) 3 + 4 (2k + 1) 3 +∞ u=1 2 (2k + 1 + 2u) 4 + 4 (2k + 1) 2 +∞ u=1 2 (2k + 1 + 2u) 5 + 4 (2k + 1) +∞ u=1 2 (2k + 1 + u) 6
To conclude this time, we need a few more steps. First we use the well-known relation between sums on even and odd indices in zeta values, i.e. for any positive integer m one has

ζ(m) = +∞ u=1 1 (2u) m + +∞ u=0 1 (2u+1) m so that +∞ u=1 1 (2u+1) m = 1 -1 2 m ζ(m) -1.
Second, the sum multiplied by -4 (2k+1) 6 is no longer finite, and we use the opposite of the alternating harmonic series: 6) in it, and had added a third statement which we realized later is in fact completely erroneous). We are very grateful to the referee of a previous version of this work for suggesting that these expressions might be established using the same strategy she/he had mentioned about her/his alternative proof of theorem 2.1 (entire series expansion and Fubini's theorem).

- +∞ u=1 (-1) u+1 u = -log(2
3. Proof of the irrationality of all ζ(2n + 1) for n ∈ N * First let us recall the standard irrationality criteria of Dirichlet :

Lemma 3.1. (Dirichlet, 1848) α ∈ R\Q ⇔ ∀ > 0 ∃p ∈ N ∃q ∈ Q such that |pα -q| < .
In 1979, shortly after Apéry presented his proof of the irrationality of ζ(2) and ζ(3), Beukers [START_REF] Beukers | A note on the irrationality of ζ(2) and ζ(3)[END_REF] found another proof using Dirichlet's criteria applied to some particular integral representations of those two numbers. We quickly summarize the strategy as follows (the author also benefited from the extremely clear slides of Brown [START_REF] Brown | Irrationality proofs for zeta values, moduli spaces and dinner parties[END_REF]).

• step 1: we have

1 0 1 0 -log(xy)
1-xy dxdy = 2ζ(3) and for any integer r ≥ 1 we have

1 0 1 0 -log(xy) 1-xy (xy) r dxdy = 2 ζ(3) -1 1 3 -• • • -1 r 3 ≤ 2ζ(3)
• step 2: by denoting a Legendre-type polynomial

P k (x) := 1 k! d dx k x k (1 -x k ) ∈ Z[X]
and using the previous step we have • step 6: we can bound uniformly for 0 ≤ x, y, z ≤ 1 one part of the integrand x(1-x)y(1-y)z(1-z)

I k := 1 0 1 0 -log(xy) 1-xy P k (x)P k (y)dxdy = A k +B k ζ(3)
1-(1-xy)z ≤ ( √ 2 -1) 4 < 1 2
(this is the reason for introducing P k rather than working with the integrand of step 1 where (xy) r can only be bounded by 1)

• step 7: by using most of the previous steps we find 0

< A k +B k ζ(3) d 3 k ≤ 2ζ(3)( √ 2 -1) 4k
• step 8: using now the information on the growth of d k , so of d 3 k too, we get 0

< |A k + B k ζ(3)| ≤ 4 5 k ,
which concludes the proof by Dirichlet's criteria.

Unfortunately in the ensuing years and decades no tweak to that strategy could be made to work for values of ζ at other odd integers. Vasilyev [START_REF] Vasilyev | On small linear forms for the values of the Riemann Zeta function at odd points[END_REF]

+ 7ζ(3) (2k + 1) 4 + 6ζ(2) (2k + 1) 5 + 8 log(2) (2k + 1) 6 = 31ζ(5) 4(2k + 1) 2 - 4 (2k + 1) 6 k i=1 1 2i + 1 -8 6 j=1 1 (2k + 1) j k i=1 1 (2i + 1) 7-j
Using the triangle inequality on the definition of I k , and determining by a routine study of critical points that the supremum of the absolute value | log(x) log(y)(xy) 2k+1 | occurs at x * = y * = e -1/(2k+1) , we get the inequality: 5 which, multiplying both sides by (2k + 1) 2 , can be rewritten as:

|I k | ≤ Sup{| log(x) log(y)(xy) 2k+1 | where 0 ≤ x, y ≤ 1} = 1 e 2 (2k+1) 2 < 1 (2k+1) 2 × [0;1] 2 log(x) x log(y) y log(1 -(xy))dxdy =ζ(5) + 4ζ (6)
(2k + 1) + 4ζ (4)
(2k + 1) 3 + 4ζ(3) (2k + 1) 4 + 4ζ(2) (2k + 1)
0 < 4ζ(5) - a k b k < ζ(5) + 4(2k + 1)ζ(6) + 4ζ(4) 2k + 1 + 4ζ(3) (2k + 1) 2 + 4ζ(2) (2k + 1) 3 where a k b k := (2k + 1) 2 4 (2k+1) 6 2k+1 n=1 1 n + 4 6 j=1 1 (2k+1) j 2k+1 i=1 1 i 7-j
> 0 is a certain irreducible rational number.

Similarly with J k we ultimately have:

0 < 31ζ(5) - c k d k < 4ζ(5) + 63(2k + 1)ζ(6) 2 + 30ζ(4) 2k + 1 + 28ζ(3) (2k + 1) 2 + 24ζ(2) (2k + 1) 3 + 32 log(2) (2k + 1) 4 where c k d k := (4(2k + 1) 2 ) 4 (2k+1) 6 k i=1 1 2i+1 + 8 6 j=1 1 (2k+1) j k i=0 1 (2i+1) 7-j
> 0 is a certain irreducible rational number.

In those two inequalities, to obtain irrationality of ζ(5) one would need that after multiplying all sides by b k (respectively d k ) the resulting right-hand side expression be a function that goes to 0 as k goes to infinity, which is not the case. This led the author to think about finding a function whose supremum on [0, 1] has a faster decay than O 1 (2k+1) 2 . For any positive integer m, we can show in a similar fashion as before the closed-form expression of the following one-dimensional integral (where s k,m and t k,m are positive integers resulting from the combined completions of the sums): While the right hand side does tend to 0 as k goes to infinity, contrarily to what was falsely claimed in a previous version of this work, this does not prove the irrationality of ζ(m + 2). For that, one would first need to multiply all sides by t k,m and that the resulting right-hand side be a function that goes to 0 as k goes to infinity. We have not computed t k,m explicitely, but it appears not to be the case.

  xy)1 ) log(x) log(y)(xy) 2k+1 dxdy =

d 3 n

 3 with A k , B k ∈ Z and d k := lcm(1, . . . , k) • step 3: by the Prime Number Theorem we have for any integer k ≥ 1 that d k < 3 k xy)z dz = -log(xy) 1-xy • step 5: by integration by parts one finds also that I k =

1 - 5 = 1 -

 151 could show that a direct generalization of Beukers's integral for ζ(5) can only show that one of ζ(3) and ζ(5) is irrational. In what follows we shall show that the expressions from the previous section do not seem to allow any progress on these matters. Define for positive integers k the two following sequences of numbers: (xy)) log(x) log(y)(xy) 2k+1 dxdy -4ζ(6) (2k + 1) + 4ζ(4) (2k + 1) 3 + 4ζ(3) (2k + 1) 4 + 4ζ(2) (2k + 1) (xy) 2 ) log(x) log(y)(xy) 2k+1 dxdy -

1 - 1 -

 11 x)(log(x)) m x 2k+1 dx = x)(log(x)) m x 2k+1 dxwe obtain, after using the fact that the supremum of |(log(x)) m x 2k+1 | occurs at x * = e -m/(2k+1) , and multiplying both sides by 2k + 1, that:0 < (m + 1)!ζ(m + 2) -(2k + 1)s k,m t k,m ≤ m m ζ(3) (2k + 1) m-1 e m + m j=1 (m + 1)!ζ(m + 2 -j) (2k + 1) j

  Remark 2.6. In the previous version of this work, this had been stated as a conjecture (in a much less legible manner since we had not yet recognized ζ(2), ζ(4) and ζ(

).