A Baillis

E Limam

Buzaud

M Brun

A Combescure

C Baillis

A Limam

E Buzaud

Simulations of the penetration of limestone targets using two dimensional multimodal fourier analysis

Keywords: projectile impact, axisymmetric shell, non axisymmetric loading, multimodal Fourier analysis

, in order to extend it for simulating penetration of targets composed of geological materials by steel projectiles. A finite element code with COMU elements is coupled with analytical forcing functions based on the dynamic expansion of a spherical cavity to represent the target. This combined analytical and computational approach is employed to reproduce normal impact of ogive-nosed VAR 4340 steel projectiles striking limestone targets with small angles of pitch and yaw.

).

A c c e p t e d m a n u s c r i p t 1. INTRODUCTION

During the last two decades, analytical models for penetration mechanics have been intensively developed.

Those application oriented tools are based on an approximation of the main physical processes implied in the interaction. In the case of a projectile striking a target, it is generally assumed that the projectile remains rigid during the penetration. The forces acting on the outer surface of the projectile during the penetration process can be estimated through empirical relations (Young [START_REF] Young | Simplified Analytical Model of Penetration with Lateral Loading : user's guide[END_REF]), or alternatively be derived from the resolution of a simplified problem. For instance, the solution of a dynamic spherical cavity expansion problem can be used as an input to the penetration equations, by assuming that the stress on the outer surface of projectile through penetration is determined by the stress on the surface cavity. This approach has been developed and adapted to several materials such as soil, concrete or limestone in Refs [START_REF] Forrestal | Penetration into soil targets[END_REF][START_REF]Penetration into dry porous rock[END_REF][START_REF] Luk | Penetration into semi-infinite reinforced-concrete targets with spherical and ogival nose projectiles[END_REF][START_REF] Forrestal | Penetration of concrete targets with deceleration-time measurements[END_REF]. The dynamic cavity expansion problem is solved by considering appropriate constitutive laws for the target materials. Functional forms used for the target representation incorporate an unknown target resistance strength parameter that is obtained from penetration depth versus striking velocity data.

Fully three-dimensional penetration calculations with a purely computational approach (finite element method for both projectile and target in Lagrangian, Eulerian or Arbitrary Lagrangian-Eulerian formulations) enable full structural dynamics to be reproduced, but the time required prohibits any prospect for use in an overall design tool. Recently, a combined approach using finite element method for the projectile and analytical functions for the target has been proposed for the simulation of penetration of aluminium, limestone or concrete target by spherical and ogive-nosed steel projectiles as discussed in [START_REF] Warren | Simulations of the penetration of limestone targets by ogive-nose 4340 steel projectiles[END_REF] and in Refs. [START_REF] Warren | Penetration of limestone targets by ogive-nosed VAR 4340 steel projectiles at oblique angles: experiments and simulations[END_REF][START_REF] Warren | Simulations of the penetration of 6061-T6511 aluminium targets by spherical-nosed VAR 4340 steel projectiles[END_REF][START_REF] Warren | Penetration of 6061-T6511 aluminium targets by ogive-nosed VAR 4340 steel projectiles at oblique angles[END_REF][START_REF] Warren | Penetration into low-strength (23 MPa) concrete: target characterization and simulations[END_REF][START_REF] Danielson | A meshless treatment of three-dimensional penetrator targets for parallel computation[END_REF]. In addition to the costefficiency of this combined approach, this method eliminates the need for a contact algorithm between target and projectile, and avoids the problem due to excessive mesh distorsion for the target in Lagrangian formulation. The finite element used in these computation was a classical three-dimensional hexahedral continuum elements for the projectile modelling. In this paper, a two-dimensional quasi axisymmetric element is fitted to impact problems in order to enhance even more the cost efficiency of the combined analytical and computational approach.

The finite element proposed in this work is an advanced version of the quasi-axisymmetric element COMU developed by the French Research Institute for Nuclear Energy (Commissariat à l'Energie Atomique) to predict buckling of thin shell structures for design. Buckling process in thin shell structures is particularly sensitive to the imperfections of the shells. In the COMU element formulation, the imperfections and displacements of the

A c c e p t e d m a n u s c r i p t

shell structure are described by Fourier series on circumferences of the quasi-axisymmetric structure following an approach proposed by Wunderlich et al. ([15]). The behaviour of the structure submitted to static or dynamic loads is computed from a finite element analysis by selecting a series of Fourier modes, on which are decomposed all physical values such as displacements, velocities, accelerations and forces. Strains and stresses are expanded in Fourier series as well. Thanks to the modal description of physical values on circumferences of the structure, imperfections and non axisymmetric loads are taken into account. So, the formulation enables to deal with three-dimensional problems, such as projectile impacting complex targets, while keeping the advantage of the two-dimensional analysis with low computation costs. To summarize, the COMU element has the following main features: non axisymmetric imperfection, non axisymmetric load, large displacements, non linear behaviour, low computation costs.

The initial formulation of this element has been improved in order to reproduce large rotations, that projectile may experience during the penetration into the target. To this end, the initial formulation based on an updated Lagrangian formulation is changed into a total Lagrangian formulation. In the following, the general formulation of the advanced COMU element will be presented in detail by pointing out the differences between the previous formulation and the new one. The COMU element developed in a total Lagrangian formulation is employed to reproduce penetration experiments reported by Frew et al. [START_REF] Frew | Penetration Experiments With Limestone Targets and Ogive-Nose Steel Projectiles[END_REF] involving limestone targets and ogive-nosed VAR 4340 steel projectiles. The target model is based on Forrestal's semi-empirical closed-form expression for the final depth of penetration. Values of limestone target resistance obtained by fitting the depth of penetration versus striking velocity data curve, are given by the authors. For the reproduction of the loss of confinement due to the entrance cratering effects, the methodology of successive layers proposed by Warren [START_REF] Warren | Simulations of the penetration of limestone targets by ogive-nose 4340 steel projectiles[END_REF] is adopted.

Results from non linear dynamic analyses with advanced COMU elements in terms of depth of penetration and projectile deformation through penetration are compared with experimental data and numerical results obtained by Warren using a three-dimensional finite element for the projectile.

GENERAL FORMULATION OF STRAINS FOR A SHELL

General formulation of strains for a geometrically perfect shell

In this first section, a geometrically perfect shell is considered. In other words, the initial imperfections of the shell structure are neglected. , where is the inward normal of the shell, n s is the tangent vector along the meridian direction and t is the tangent vector along the circumference of the structure. We consider a conical shell, that is the curvature radius in the meridian plane ) , ( s n , is assumed as infinite. In the meridian plane, the orientation of the conical shell is given by the angle ( )

s e z , - = ϕ .
The displacement field in the local coordinates is :

{ β v w u q t = } ( 1 )
where is the axial tangential displacement, is the normal displacement to the shell, is the circumferential tangential displacement, and

u w v
β is the rotation of the shell around the tangential vector t .

In the cylindrical system ( ) , , θ e e e z r , the displacement field is denoted by :

{ β v w u q t = } ( 2 ) with        = = - = + = β β ϕ ϕ ϕ ϕ v v w u u w u w cos sin sin cos ( 3 )
The rotations of the shell around the tangent vector s and normal vector n are not taken into account.

The strain tensor is given by : ( ) (

q q q q q t t ∇ ∇ + ∇ + ∇ = . 2 1 2 1 ) ( ε ) ( 4 )
where the gradient operator is expressed on the perfect configuration.

The tensor can be decomposed into two parts, as follows :

) ( ) ( ) ( q q q Q L ε ε ε + = ( 5 )
where is the linear part of strains, and is the quadratic part of strains involved in buckling and large displacement analysis.

) (q L ε ) (q Q ε
Let us write the linear strains in the following vector form : For a flat shell with a infinite curvature radius in the meridian plane, the expressions of linear membrane deformations are :

A c c e p t e d m a n u s c r i p t

{ } { }       =                     = L Lm
      ∂ ∂ + - ∂ ∂ =       + + ∂ ∂ =       ∂ ∂ = s v r v u r r u r w v r s u Lm s Lm Lm ss ϕ θ ε ϕ ϕ θ ε ε θ θθ sin 1 2 1 sin cos 1 ( 7
)
and the linear changes of curvature variations are given by : [START_REF] Luk | Penetration into semi-infinite reinforced-concrete targets with spherical and ogival nose projectiles[END_REF] The shear strain in the thickness of the shell is expressed as :

        ∂ ∂ + - ∂ ∂ + ∂ ∂ ∂ - =         ∂ ∂ + ∂ ∂ - ∂ ∂ - =         ∂ ∂ - = ϕ θ ϕ ϕ ϕ θ χ θ ϕ ϕ θ χ
      + ∂ ∂ = β ε ξ s w L s 2 1 ( 9 )
The Kirchoff-Love relation, classically assumed for thin shells, states that linear shear strain s ξ ε can be neglected. Consequently, the component β of the displacement field is linked to the longitudinal derivative of the normal displacement. We have :

s w ∂ ∂ - = β ( 1 0 )
Thus, the first curvature term in (8) becomes :

s ss ∂ ∂ = β χ ( 1 1 )
For the quadratic part of the strain tensor, only membrane strains are taken into account, as classically assumed for thin shells. The non linear membrane strains are given by :

A c c e p t e d m a n u s c r i p t

            + + ∂ ∂ ∂ ∂ +       - ∂ ∂       ∂ ∂ +       ∂ ∂       - ∂ ∂ =               + + ∂ ∂ +       - ∂ ∂ +       - ∂ ∂ =               ∂ ∂ +       ∂ ∂ +       ∂ ∂ = ϕ ϕ θ ϕ θ ϕ θ ε ϕ ϕ θ ϕ θ ϕ θ ε ε θ θθ sin cos sin cos 2 1 sin cos sin cos 1 2 1 2 1 2 2 2 2 2 2 2 u w v s v v u s u s w v w r u w v v u v w r r s v s u s w Qm s Qm Qm ss ( 1 2 )
At a point M located in the thickness of the shell from a distance ξ to the neutral fibre, the total strain is given by :

{ } { } { } { } Qm L Lm T s T T ss T ε χ ξ ε ε ε ε ε θ θθ + + =           = 2 ( 1 3 ) with { }           = = Qm s Qm Qm ss Qm Q q θ θθ ε ε ε ε ε 2 ) ( ( 1 4
)

General formulation of strains for a geometrically imperfect shell

The previous section provided analytical expressions for linear and non linear strains as a function of the displacement field from a perfect axisymmetric structure. Let us now introduce an initial imperfection field.

The imperfect shell with geometrical imperfections, denoted by I, is related to a virtual perfect shell, denoted by P, through a displacement field, denoted by . When a load F is applied to the actual structure, the structure experiences a displacement field, denoted by .

q D q

The purpose is to express the strain tensor induced by the imperfection field and the incremental displacement on the perfect configuration P. The methodology, which has been first proposed by Marguerre for a general tridimensional problem, is applied here to the case of an axisymmetric structure with nonaxisymmetric imperfections of any size. On the imperfect configuration I, the strain tensor induced by the incremental displacement follows the previous expression (4).

D q

The strain tensor associated with the displacement field can be obtained on the perfect configuration P by introducing the initial displacement field :

q D ) ( ) ( ) , ( D q D q D ε ε ε - + = ( 1 5 )
After expanding the terms of the above expression, we obtain the following expression for the strain tensor expressed on the perfect configuration :

) , ( ) , ( ) , ( ) ( ) , ( q q q D D q q q D Q Q Q L NL ε ε ε ε ε + + + = ( 1 6 )

A c c e p t e d m a n u s c r i p t

where quadratic terms of the strain tensor are defined by the quadratic operator :

[ 2 1 2 1 . 2 1 ) , ( q q q q t Q ∇ ∇ = ε ] ] ( 1 7 )
the gradient operator of the displacement fields and being defined on the perfect configuration P.

1 q 2 q
It is worth noting that, at this point, no assumption on the size of initial imperfections and displacements is made.

The original version of the quasi-axymmetric COMU element assumes that the incremental displacement is small with respect to the initial imperfection , so that the quadratic term of displacement could be neglected.

q D ) , ( q q Q ε When considering an initial perfect axisymmetric structure, the linear and non linear strains coming from the displacements of the shell are given in Equations ( 7) to [START_REF] Warren | Penetration of 6061-T6511 aluminium targets by ogive-nosed VAR 4340 steel projectiles at oblique angles[END_REF] in local coordinates. In a classical finite element formulation, these strains can be linked to the displacements of the nodes of the shell element through divergence operators. After expanding the imperfection displacements in Fourier series, it will be shown in the following that the divergence operators for strains can be also expanded in Fourier series.

Expressions of the linear strains and non linear strains for the perfect shell will be further employed to derive the total strains for an imperfect shell. The general outline of the formulation is to express strain and stress Fourier coefficients for a given COMU element as a function of Fourier coefficients of the nodal displacements in cylindrical coordinates.

) (q L ε ( NL ε ) , ( q q Q ε ) , q D

FOURIER DECOMPOSITION OF STRAINS AND STRESSES FOR A GEOMETRICALLY

IMPERFECT SHELL

Decomposition of the displacement field into Fourier series

By referring to the perfect configuration in local coordinates, the displacement field can be expanded in a Fourier series:

[ ∑ = + + = N n A S n q n q q q 1 0 ) ( ) ( ( 1 8 )

A c c e p t e d m a n u s c r i p t

where is a vector composed of four Fourier coefficients on the axisymmetric mode, is associated with the symmetric harmonic (symmetric with respect to the circumferential angle and, is associated with the antisymmetric harmonic (symmetric with respect to the circumferential angle

0 q n ) (n q S 0 = θ ),
) (n q A n 2 π θ = ).
So, in local coordinates, the four modal displacement vectors are :

, ) 0 ( ) 0 ( ) 0 ( ) 0 ( ) 0 (             = = = = = = n n v n w n u n q β ( 1 9 ) , cos ) ( sin ) ( cos ) ( cos ) ( ) (             = θ β θ θ θ n n n n v n n w n n u n q S S S S S             = θ β θ θ θ n n n n v n n w n n u n q A A A A A sin ) ( cos ) ( sin ) ( sin ) ( ) (
In the following, for the sake of clarity, only symmetric harmonics will be considered. Formula for antisymmetric harmonics are deduced from those for symmetric harmonics by swapping

Geometry of the finite element and vector notations

The COMU element is a flat shell element, with two nodes i and j . The geometry of the shell element is defined by the longitudinal length , the radius at the middle of the element , the thickness , and the orientation angle

L m R e ) ,
( s e z -= ϕ . The shape functions are linear. In local coordinates, we have :

j i f x f x f ) 1 ( 2
1 ) 1 ( 2 1 + + - = ( 2 0 )
where and refer to one of degrees of freedom, that is the displacements ( , , v ) or the rotation (

i f j f u w β ),
at nodes i and j in local coordinates, respectively ; x is the finite element parametric variable along the longitudinal local direction s , equal to -1 at node i , 0 at the middle the element and +1 at node j .

The partial derivative along the longitudinal direction is given by :

x f L s f ∂ ∂ = ∂ ∂ 2 ( 2 1 )
where can be replaced with the displacements ( u , , v ) or the rotation ( f w β ).

The displacement field in local coordinates can be expressed in cylindrical coordinates according to (3).

In the cylindrical reference, let us write the displacement vector for the shell element on the n symmetric harmonic ; this one is composed of two parts associated with the two nodes i and j of the element under consideration :

A c c e p t e d m a n u s c r i p t

{ } { } { }       = ) ( ) ( ) ( n q n q n q S j S i S ( 2 2 )
where displacement vectors of nodes i and j on the n symmetric harmonic are written as :

{ } [ ] { } [ ] ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( n n v n w n u n q n n v n w n u n q S j S j S j S j S j t S i S i S i S i S i t β β = = ( 2 3 )
The axisymmetric { } ) 0 ( = n q follow an analogous expression.

Note that these vectors are composed of Fourier coefficients. In the cylindrical reference, the total displacement of the node i is obtained by recomposing the Fourier series :

             + = = + = = + = = + = = ∑ ∑ ∑ ∑ = = = = N n S i i i N n S i i i N n S i i i N n S i i i n n n n n v n v v n n w n w w n n u n u u 1 1 1 1 cos ) ( ) 0 ( sin ) ( ) 0 ( cos ) ( ) 0 ( cos ) ( ) 0 ( θ β β β θ θ θ ( 2 4 )
In the following, we will express the linear and non linear strains in the (i,j) element as a function of the modal

displacements { } ) 0 ( = n q and { } ) (n q S
in the cylindrical reference, by considering an element without geometric imperfection. The case of the geometrically imperfect element with an initial imperfection field will be later deduced from these expressions.

D

Fourier decomposition of linear strains for a geometrically perfect element

We seek to express the linear strains , in terms of modal displacement vectors

) (q L ε { } ) 0 ( = n q and { } ) (n q S in cylindrical coordinates.
Firstly, let us consider the local displacements on the axisymmetric mode :

) 0 ( = = n u u , ) 0 ( = = n w w , , ) 0 ( = = n v v ) 0 ( = = n β β
. The expressions of linear strains, given in Equations ( 7) to [START_REF] Warren | Simulations of the penetration of 6061-T6511 aluminium targets by spherical-nosed VAR 4340 steel projectiles[END_REF], are developed by using the shape functions of the element (20), the relations between cylindrical coordinates and local coordinates (3), and the relation on the longitudinal derivative (21), in order to express linear strains as a function of the axisymmetric displacement { } ) 0 ( = n q in cylindrical coordinates. Finally, the linear strains associated with the axisymmetric mode of displacement can be written as : The same methodology is applied to the local displacements on the n symmetric. The n symmetric local displacement is defined by :

{ } [ ][ ] { ) 0 ( 2 2 ) ( 0 = = = n q B q L s L L ss Lm s Lm Lm ss t L θ θθ θ θθ χ χ χ ε ε ε ε } ( 2 5 )
θ n n u cos ) ( = u , θ n n w w cos ) ( = , θ n n v v sin ) ( = , θ β β n n cos ) ( = , θ
{ } [ ]{ } [ ]{ } θ θ ε n n q n B n n q n B q S S s S S c L sin ) ( ) ( cos ) ( ) ( ) ( + = ( 2 6 )
All the matrices [ introduced in the above expressions can be found in [START_REF] Combescure | Etude de la stabilité non linéaire géométrique et non linéaire matériau des coques minces. Application aux coques de révolution avec imperfections soumises à des chargements complexes[END_REF].

] B By adding the contributions of all the harmonics of the displacement and grouping the cosine and sine terms, the complete expression of the linear strains in the (i,j) element is expressed as :

{ } { } { } { } [ ]{ } [ ]{ } ( ) [ ]{ } ( ) θ θ χ ε ε n n q n B n n q n B n q B q q S S s N n S S c L Lm L sin ) ( ) ( cos ) ( ) ( ) 0 ( ) ( ) ( 1 0 + + = =       = Β = ∑ = (27)

Fourier decomposition of non linear strains for a geometrically perfect element

As previously carried out for the linear strains, the purpose is to link the non linear strains to the modal displacement vectors of the two-nodes element.

) (q Q ε Let us rewrite the non linear strain vector, given in [START_REF] Warren | Penetration of 6061-T6511 aluminium targets by ogive-nosed VAR 4340 steel projectiles at oblique angles[END_REF], in the following product form :

{ } [ ]{ ) ( ) ( 2 1 2 ) , ( q q H q q Qm s Qm Qm ss Q Θ =           = θ θθ ε ε ε ε } ( 2 8 )
where the 6 components vector and 3-by-6 components matrix are given by :

{ } { } { } [ ] { } [ [ ] { } { } {           Θ Θ Θ Θ = Η       Θ Θ =                               + + ∂ ∂       - ∂ ∂       - ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ = Θ ) ( ) ( ) ( 0 0 0 0 0 0 ) ( ) ( ) ( ) ( sin cos 1 sin 1 cos 1 ) ( 1 2 2 1 2 1 q q q q q and q q u w v r v u r v w r s v s u s w q t t t t ϕ ϕ θ ϕ θ ϕ θ ] } (29)

A c c e p t e d m a n u s c r i p t

In order to derive the Fourier series of the non linear strains, we consider successively each modal displacements, following an analogous methodology as previously described.

Let us consider the axisymmetric local displacement. Using the relations ( 20), ( 3) and (21) (shape functions, transfer matrix and relation on the longitudinal derivative, respectively), the vector { } ) (q Θ can be written as a product of a new 6-by-8 matrix [ and the axisymmetric displacement, as follows :

]

0 G { } { } { } [ ] [ ]{ } ) 0 ( ) ( ) ( ) ( 0 2 1 = = Θ Θ = Θ n q G q q q t t ( 3 0 )
From the above relation, the expression of the matrix [ ]

) (q Η
on the axisymmetric mode is easily derived and is written as :

[ ] [ 0 ) ( A q = Η ] ( 3 1 )
If we consider the n symmetric mode, matrices [ ]

) (n G S c , [ ] ) (n G S s relate the vector { to the modal displacements : } ) (q Θ { } [ ]{ } [ ]{ } θ θ n n q n G n n q n G q S S s S S c sin ) ( ) ( cos ) ( ) ( ) ( + = Θ ( 3 2 )
From the above expressions, matrices for the n symmetric displacement is given by :

[ ) (q Η ] [ ] [ ] [ ] ( ) θ θ n n A n n A q S s S c sin ) ( cos ) ( ) ( + = Η ( 3 3 )

Complete expressions of the matrices [ and ]

G

[ ]

A can be found in the [START_REF] Combescure | Etude de la stabilité non linéaire géométrique et non linéaire matériau des coques minces. Application aux coques de révolution avec imperfections soumises à des chargements complexes[END_REF].

Finally, by taking into account all the harmonics of the local displacement in the (i,j) element, we derive the complete expression of the non linear strains :

{ } [ ] [ ] [ ] ( )[ ]{ } [ ] [ ] ( ) {       + +       + + = ∑ ∑ = = N p S S s S c N n S s S c Q p q p p G p p G q G n n A n n A A q q 1 0 0 1 0 ) ( sin ) ( cos ) ( sin ) ( cos ) ( 2 1 ) , ( θ θ θ θ ε } (34)
Products of two sine or cosine functions appear in the above expression. As a result, the Fourier series of non linear strains is truncated at an order corresponding to twice the order N of displacements.

Fourier decomposition of linear and non linear strains for a geometrically imperfect shell

The formulation of strains associated with a geometrically perfect shell enables the strains in the case of a geometrically imperfect shell to be expressed by substituting into the non linear part of the strain the imperfection field for the displacement .

D q

As with the displacement field , let us expand the imperfection field into a Fourier series: Note that, unlike the original formulation of the quasi axisymmetric COMU element, the imperfection field incorporates here the axisymmetric harmonic. As a consequence, when non linear analysis is performed, the geometry of the perfect structure, on which the formulation is based, is not updated any more at each time step according to the incremental displacement associated with the mode axisymmetric. In this new formulation, all the geometry changes are taken into account in the imperfection field. Thus, the integration strategy of dynamic equations proposed in this paper uses a total Lagrangian formulation instead of an original updated Lagrangian formulation.

After changing with in the Equation (28), the linear and non linear strains (given in ( 16)) become :

q D { } { } { } [ ]{ } [ ]{ } [ ]{ } ) ( ) ( 2
1 ) ( ) ( 2
1 ) ( ) ( 2
1 ) , ( ) ( ) ( D q q D q q q D q q NL L Θ Η + Θ Η + Θ Η = Β = ε ε ( 3 6 )
where corresponds to the Fourier series of the linear part of the strains, whose complete expression can be decomposed in Fourier series in terms of

{ ) (q Β } [ ] B matrices, as expressed in (27).
From (28), it is easily shown that we have :

[ ]{ } [ ]{ } ) ( ) ( ) ( ) ( D q q D Θ Η = Θ Η (37)
So, we obtain the reduced expression of the non linear strains :

{ } { ) ( ) 2 1 ( ) , ( q q D q D NL Θ       + Η = ε } } ( 3 8 ) { ) (q Θ is expressed in terms of products of matrices [ ] G by the displacement . q       + Η ) 2 1 ( q D
is expressed in terms of matrices [ ] A , which can be obtained from (29) holding for a perfect shell, by replacing the field displacement with the field

q q D 2 1
+
, that is, the sum of the imperfection field with the half of the displacement field.

Total Lagangian formulation

The computation follows the classical explicit scheme in a total Lagrangian formulation. The imperfection field is equal to the displacement from the initial configuration at the impact time to the actual configuration at the

∆ ] B { } ( ) ( q    Η ∆
The matrices [ , and are expressed here on the perfect initial configuration and are computed only at the beginning of the non linear analysis. The operators

[G]

[ ]

A depend on the field q D ∆ + 2
1
, that is the displacement at the beginning of the time step added to the half of the incremental displacement. So, operators must be computed at each time step of the explicit scheme.

[ ]

A
At a given time step, the incremental displacement q ∆ is deduced from the results on the previous time step defined by q , velocity and acceleration . The incremental linear and non linear strains are derived from the incremental displacement as follows :

q q { } { } { } ) ( ) 2 1 ) , ( ) ( q q D q D q NL L ∆ Θ    ∆ + = ∆ ∆ Β = ∆ ∆ ε ε ( 3 9 )
The vectors and matrices involved in the above formula have been developed on a modal Fourier basis in Equations ( 27) and (38), by introducing the matrices [ ]

B , [ ] A and [ ] G .
As discussed previously, it is important to note that the axisymmetric mode is incorporated into the imperfection field . This point is different from the original version of the quasi axisymmetric element, in which only non axisymmetric displacements are considered as imperfections. The original version based on an updated Lagrangian method is replaced by a total Lagrangian method.

D

Moreover, it is worth noting that the strain term [ ] ( )

{ q q ∆ Θ ∆ Η ) ( 2 1
} in the full expression of strains (39) is taken into account in this new formulation at contrast to the original version of COMU element, in which the incremental displacement q ∆ was assumed to be small with respect to the imperfection field . D

Linear elastic stress-strain law

In large displacement analysis, we assume in this study that the Hooke's law is the same if expressed on the imperfect configuration or on the perfect one.

Under a plane stress assumption, the stresses are easily obtained from the strains. At a point located in the thickness of the shell from a distance ξ to the neutral fibre, the linear stress vector, composed of membrane and bending components, is :

A c c e p t e d m a n u s c r i p t

{ } [ ] [ ] [ ] { ) ( 0 0 ) ( q H H q L Lb s Lb Lb ss Lm s Lm Lm ss t L ε ξ σ σ σ σ σ σ σ θ θθ θ θθ       = = } ] (40)
where the [ is the law matrix for an homogeneous, isotropic material under the plane stress assumption, defined by :

H [ ]             - - = 2 1 0 0 0 1 0 1 1 2 ν ν ν ν E H ( 4 1 )
E being the Young's Modulus and ν the Poisson's ratio.

The non linear stress contains only membrane components and is written as :

{ } [ ] [ ]{ } ) , ( 2 ) , ( q D H q D NL Qm s Qm Qm ss t NL ε σ σ σ σ θ θθ = = ( 4 2 )

ANALYTICAL MODEL AND NUMERICAL SIMULATION

The approach adopted in this work lied in coupling analytical expressions representing the target resistance with a finite element method using COMU elements in a classical explicit scheme. First, we summarize the target resistance model based on Forrestal's semi-empirical closed-form expression for the final depth of penetration.

Then, the target resistance model is incorporated into an explicit scheme involving COMU elements. The coupling technique between analytical forcing functions representing the target and the two-dimensional COMU modelling for the projectile is presented. Finally, the constitutive model used for the projectile is briefly described.

Target resistance model

The target resistance functions are implemented into the CalPen3D code (Sibeaud et al. [START_REF] Sibeaud | CalPen3D: A three-dimensional penetration code suitable for Hardened Target Defeat Modeling -Numerical and Experimental assessment[END_REF], Buzaud et al. [START_REF] Buzaud | An experimental investigation of corner effects resulting from vertical attack on hardened structures[END_REF]) designed to provide an accurate determination of the curvilinear trajectory of a weapon while travelling through a three-dimensional target composed of concrete, rock, soil and air blocks. It has been extensively validated against a large experimental database in terms of depth of penetration, trajectory and time-history deceleration.

CalPen3D has been developed in Java language, hence, it can be run on any operating system without compilation, provided the installation of the Java Virtual Machine.

A c c e p t e d m a n u s c r i p t

The algorithm of CalPen3D is essentially based on the computation and summation of forces generated at the centre of each elementary four nodes flat surfaces on the outer surface of the projectile. Normal stresses generated at the projectile-target interface are obtained by using the spherical cavity-expansion solution (Refs [START_REF] Forrestal | Penetration into soil targets[END_REF][START_REF]Penetration into dry porous rock[END_REF][START_REF] Luk | Penetration into semi-infinite reinforced-concrete targets with spherical and ogival nose projectiles[END_REF][START_REF] Forrestal | Penetration of concrete targets with deceleration-time measurements[END_REF]).

At each time step, the position of the centre of each element is checked to determine if a target element has been penetrated. In this case, local forces are applied at this position to resist the penetration of the projectile. Those forces are approximated assuming the expansion of a spherical cavity solution with the corresponding local normal velocity V . This results in a polynomial relation between the normal stress n n σ and the normal velocity at the centre of the element:

n V 2 0 n n V R ρ σ + = ( 4 3 )
where 0 ρ is the target density and R is the target strength resistance.

For limestone targets, Frew et al. [START_REF] Frew | Penetration Experiments With Limestone Targets and Ogive-Nose Steel Projectiles[END_REF] deduced from impact experiments involving ogive-nosed rod projectiles the following relation between the target resistance and projectile diameter: Post-test observations [START_REF] Frew | Penetration Experiments With Limestone Targets and Ogive-Nose Steel Projectiles[END_REF] have exhibited a conical cratering region that is approximately two diameters in depth. This cratering region is followed by a tunneling region that is approximately one projectile diameter. Therefore,

      + = a a k K R 2 2 0 ( 4 4

A c c e p t e d m a n u s c r i p t

the first region needs a special treatment to take into account the less resistance due to the crater formation. We use the same methodology as proposed by Warren [START_REF] Warren | Simulations of the penetration of limestone targets by ogive-nose 4340 steel projectiles[END_REF]. The two diameters in depth cratering region is subdivided into a suite of uniformly spaced layers with increasing strength. The normal stress acting on the projectile outer surface is given by : In this region, normal stresses are only allowed to act on the nose of the projectile and not the shank. In the tunneling region where the target material is confined and cannot be ejected out, unreduced normal stresses are applied to the whole outer surface of the projectile.

n i i σ ψ σ = ( 4

Non linear dynamic analyses

The analytical and numerical combined method as employed by Warren [START_REF] Warren | Simulations of the penetration of limestone targets by ogive-nose 4340 steel projectiles[END_REF] is adopted in this work in order to simulate ogive-nosed steel rod projectiles impacting limestone targets. An explicit dynamic scheme with COMU elements has been implemented in the CalPend3D code in order to couple finite element method for the projectile with the analytical forcing functions representing the target.

The equation of equilibrium governing the dynamic response of the shell structure modelled with COMU elements is :

int F F U M ext - = ( 4 6 )
where M is the mass matrix decomposed in Fourier series, U contains the Fourier coefficients of the nodal accelerations, ext F and contain the Fourier coefficients of the external and internal forces, which will be explained in the following. int F Non linear computations are here carried out using an explicit dynamic scheme. This method allows to use a diagonal mass matrix M (Belytschko et al. [START_REF] Belytschko | Nonlinear Finite Elements for Continua and Structures[END_REF]) , which components are associated with the Fourier harmonics taken into account in the computations.

All the vectors involved in the momentum conservation are decomposed into Fourier series. Thus, for a given i node, Fourier coefficients of the acceleration are organised as follows : axisymmetric mode (p=0), p symmetric mode followed by p antisymmetric mode, with

N p … 1 =
, N being the chosen truncature order for the displacement expansion in Fourier series. So, each node of the modelled axisymmetric structure has

N × + 8 4

A c c e p t e d m a n u s c r i p t

degrees of freedom. Note that a complete Fourier base truncated at the order N is used in this formulation, that is, all the harmonics of an order less than N are taken into account into the Fourier series.

P

In this work, the complete description of the projectile is provided by the Fourier coefficients displacement and velocity, contained in the global vectors U and U involved in the equilibrium equation. For each COMU node, P points regularly distributed along the circumference associated with the COMU node are taken into account.

The displacement and velocity at each point along the circumference are obtained by recomposing Fourier series.

So the position and velocity of each point are computed at each time step with respect to the initial referential, associated with the initial configuration of the structure at the impact time. The coordinates and velocities of these points in the target referential are then easily obtained by using a transfer matrix from the initial referential

COMU element coupling with then target resistance model

The two-dimensional finite element mesh with quasi-axisymmetric COMU elements to model the ogive-nosed projectiles in all of the simulations is illustrated in Fig. 1. The projectile is modelled with 10 elements for the nose and 20 elements in the shank part. The mesh is defined from a two-dimensional profile of 30 COMU nodes, located at the neutral fibre of the rod projectile. The forcing function due to the target resistance through the penetration process acts on the outer surface of the projectile. As discussed previously, displacements and velocities at a set of points regularly distributed along the circumference are known by recomposing Fourier series. As theses points are located at the neutral fibre, a technique is set up to reconstruct at each time step of the explicit scheme the actual configuration of the projectile outer surface. This is explained in the following.

For each circumference, we can deduce the point at the axis from the set of points along the circumference.

For a given circumference corresponding to a COMU node, we denote this point located at the axis Q . From each point Q located along the circumference at the neutral fibre, we can obtain an associated point ' located on the outer surface of the projectile from:

A Q 2 ' e u u u u + = ( 4 7 )

A c c e p t e d m a n u s c r i p t

where is the vector from the point at the axis Q to Q at the neutral fibre and ' u A u is the vector from to ' located at the outer surface, and e is equal to the thickness associated with the considered circumference. In the general case, this corresponds to the shell thickness. Here, the thickness corresponds to the shank radius of the projectile rod.

A Q Q With this method, the actual configuration of the outer surface of the projectile can be constructed at each time step. Fig. 2 illustrates the reconstruction technique of the outer surface from a deformed configuration at a given time during the penetration process: this configuration selected for illustration has been obtained in the initial impact referential at 260 µs for a 12.7 mm diameter projectile with a striking velocity of 457 m/s. This simulation will be presented in the section devoted to the numerical results. As it can be seen, the outer surface is finally composed of quadrilateral elementary flat surfaces.

On a given elementary quadrilateral side defined by 4 points composing the outer flat surface, the force acting on this side is given by :

( ) [ ]da n V R df I 2 0 . ρ ψ + = ( 4 8 )
where I V is equal to the mean velocity of the 4 points at the neutral fibre deduced by recomposing Fourier series of velocity, is the outward unit vector normal to the diagonals of the outer quadrilateral side, n ψ accounts for the cratering region, and da is the surface of the quadrilateral side. This force is then distributed to the 4 points at the medium fiber. The forces are then computed for the whole outer surface and are expanded along each circumference in Fourier series. Finally, we obtain the applied load vector involved in the equation of motion (46) whose components are Fourier coefficients of the external force decomposed on the projectile circumferences.

Constitutive model for VAR 4340 steel projectile

The constitutive model adopted is a simplified Johnson-Cook type model (assuming the material elastic-plastic without temperature influence), with a yield criterion for plasticity defined by : ( )

n e p B p f + - = σ σ σ ) , ( ( 4 9 )
where σ is the Von Mises stress, e σ is the initial value of the yield surface, p is the effective plastic strain, B and are the coefficients of the plastic law.

n

A c c e p t e d m a n u s c r i p t

The model parameters are chosen so as to fit in a best way the constitutive model used by Warren [START_REF] Warren | Simulations of the penetration of limestone targets by ogive-nose 4340 steel projectiles[END_REF]. As the model adopted in this work is a simple elastic-plastic model, we chose to fit the parameters against the elasticviscoplastic model used by Warren at a nominal strain rate of 0.001. Fig. 3 At a time step of the dynamic explicit scheme, the plasticity is obtained for the set of points along each circumference of the projectile. By using a very classical initial stress method, the increment of plastic stress is numerically obtained. The plastic stresses around the circumference are then expanded in Fourier series with an order equal to . It is important to underline that the order of the truncated stress Fourier series must be sufficiently high so as to correctly reproduce the new distribution of plastic stresses. Otherwise, important bias could be engendered in the modal description of the physical plastic stresses.

σ N σ N
These new Fourier coefficients of stresses are then used to calculate the internal force . Due to the modal formulation of the element, the integration strategy employed here is particularly effective. The nodal forces are computed by integrating analytically the stresses on the circumference. Integration along the meridional direction is numerical and uses one Gauss point at the middle of the COMU element. Integration procedure can be found in Refs. [START_REF] Combescure | Static and dynamic buckling of large thin shells : design procedure, computation tools, physical understanding of the mechanisms[END_REF][START_REF] Combescure | Etude de la stabilité non linéaire géométrique et non linéaire matériau des coques minces. Application aux coques de révolution avec imperfections soumises à des chargements complexes[END_REF]. As discussed previously, the external force is computed at each time step as well. A stabilisation force is added to avoid spurious vibrations. Indeed, Hourglass control is required because COMU element uses one integration point and under-integrates the element resulting in a rank deficiency which manifests itself into spurious zero energy modes that must be constrained. The stabilisation procedure used in this work is based on the transverse shear energy, which acts as a penalty that enforces the Kirchoff-Love constraint, given in [START_REF] Warren | Penetration of limestone targets by ogive-nosed VAR 4340 steel projectiles at oblique angles: experiments and simulations[END_REF]. Transverse shears are assumed to be small, and the material response to such deformation is assumed to be linear elastic. Details about this classic method can be found in Hughes et al. [START_REF] Hughes | A simple and efficient finite element for plate bending[END_REF] and Zienkiewicz et al. [20] works. Finally, the Fourier coefficients of forces are incorporated into the equation of motion (46).

A c c e p t e d m a n u s c r i p t

In the following, results from non linear dynamic analyses with two-dimensional COMU elements are compared with experimental data obtained by Frew et al. [START_REF] Frew | Penetration Experiments With Limestone Targets and Ogive-Nose Steel Projectiles[END_REF] and numerical results published by Warren [START_REF] Warren | Simulations of the penetration of limestone targets by ogive-nose 4340 steel projectiles[END_REF] using the same target model and three-dimensional constant strain hexahedral continuum elements. The purpose is to validate the low-cost numerical approach using the two-dimensional modal element COMU proposed in this paper.

The three dimensional outer surface reconstructed from the two-dimensional COMU profile is illustrated in Fig. 4. The same mesh is used in all of the simulations for the 25.4-mm-Diameter, 12.7-mm-Diameter and 7.1-mm-Diameter projectiles. [START_REF] Belytschko | Nonlinear Finite Elements for Continua and Structures[END_REF]. Along each circumference, 20 points regularly spaced are taken into account for the plasticity control. As explained in section 4.4, stress Fourier series must be sufficiently high so as to correctly reproduce the distribution of plastic stresses along circumferences. So we take an order of stress Fourier series equal to 10 that is 21 harmonics. All the simulations were run on a Pentium 4 (2.8 Ghz) computer with Windows system and the release 1.4 of the Sun Java 2 Platform, Standard Edition (J2SE) environment. Computation times were between 500 and 1200 CPUs depending on striking velocity. The computations times are approximately an order lower than that is reported with hexahedral elements (Warren [4], Warren et al. [START_REF] Warren | Penetration of limestone targets by ogive-nosed VAR 4340 steel projectiles at oblique angles: experiments and simulations[END_REF]). This is essentially due to the use of a twodimensional mesh, the analytical stress integration along circumferences, and the fact that only one element COMU is required in the thickness of the rod projectile. It is also important to note that computations times given here are related to the interpreted Java language which has not the performance of a compiled language like Fortran. Authors believe that a more rigorous benchmark would be more clearly in favour of COMU approach.

The 12.7-mm-Diameter projectiles

In Table 1, it can be seen that numerical results in terms of depth of penetration using COMU element are consistent with experimental data. The analytical functions representing the target resistance applied to the outer surface of the projectile through penetration are derived from depth of penetration versus striking velocity data.

A c c e p t e d m a n u s c r i p t

Therefore results validate the combined approach coupling finite element method and analytical functions to represent the target resistance.

Results obtained with COMU elements are also compared with numerical results using constant strain hexahedral continuum elements obtained by Warren. It is observed that depths of penetration are in good agreement with a maximum relative difference of 6 %. The same simulations as selected by Warren are chosen to compare the projectile deformations through penetration process. In Fig. 5, images at three specific times are

shown for an impact with a velocity of 459 m/s and 3.55° of angle of inclination. The projectile deformations predicted by COMU elements are in very good agreement with those from three dimensional simulations. It is shown that the tail bends noticeably to the right, then reverses to the left. When travelling through the tunneling region, the projectile seems to be straightened out until it comes to rest at approximately 550 µs.

With a higher velocity of 1134 m/s and a 2.15° angle of inclination, the travel of the projectile through the target is illustrated in Fig. 6. Projectile configurations at three specific times are compared with those from Warren's computations. The bending of the projectile trough penetration is significantly higher than the previous simulation. The final configuration, which is consistent with the final configuration obtained by Warren and the experimental configuration after extracting the projectile from the target (photograph published by Frew et al.

[3]), shows that the projectile is straightened out in the tunneling region. As it is has been underlined by Warren, COMU computations confirmed that projectile experiences strong bending at intermediate states although the recovered projectiles didn't show high levels of deformation. Note that depth of penetration predicted by COMU approach is higher than the 3D element approach (0.645 m and 0.608 m, respectively), but the projectile configurations at 90 µs and 150 µs with both 2D modal and 3D approaches are consistent.

The 25.4-mm-Diameter projectiles

Results in terms of depths of penetration are summarized in Table 1. The discrepancy between the simulations using COMU elements or 3D elements is very low with a relative difference less of 1 %. Inclination angles are low in all of these simulations so the trajectories through the limestone target are nearly straight. As the bending is also very reduced, projectile deformations are not shown for these simulations.

The 7.1-mm-Diameter projectiles

In Table 1 it is shown that discrepancy in terms of depth of penetration between the COMU simulation and 3D element simulation is less than 8 %, obtained with the lowest striking velocity. Otherwise the discrepancy

A c c e p t e d m a n u s c r i p t

remains of an order of 3 %. The simulation of a 7.1-mm-Diameter projectile striking the limestone target with a velocity of 787 m/s and 1.25° angle of inclination is illustrated in Fig. 7 at three specific times. The bending of the tail to the left and to the right is well reproduced by the proposed COMU approach. With a higher velocity equal to 1340 m/s and a 2.01° angle of inclination, the tails whips back and forth into the tunneling region as shown in Fig. 8. Nonetheless it is observed that the projectile final configuration is slightly more tilted with COMU elements than 3D elements.

A final representation of the projectile travel through the limestone target is presented in Fig. 9 for a striking velocity of 1365 m/s and a 2.56° angle of inclination. Projectile trajectory through the target is obtained by plotting in the same figure different configurations at 125 different times (step time of 5 µs). A zoom on the tunneling region permits to exhibit the strong whipping of the tail. As Warren's computations, COMU simulations underpredict to a less extent depths of penetration for impact velocities of 1340 m/s and 1365 m/s with high angles of inclination. The reason conjectured by the authors lies in the fact that current expansion algorithm applies a load on the shank whenever there is a component of velocity in the outward normal direction and does not account for moved materials in the tunneling region.

CONCLUSION

The purpose of this work is to set up an effective two-dimensional modal finite element approach for simulating the dynamic structural behaviour of the projectile during a three-dimensional interaction with an infrastructure target. To this end, a modified formulation of the quasi-axisymmetric element COMU, based on the expansion of all the physical quantities into Fourier series, is proposed. This new formulation of 2D modal element COMU developed according to a total Lagrangian formulation has been detailed.

As regards penetration problems, the explicit scheme using new COMU elements has been implemented into the CalPen3D code. This simulation software uses an analytical approach, based on a spherical cavity-expansion approximation. The time integration in accordance with an explicit scheme using modified COMU elements allows to couple the two approaches, that is the finite element method to represent the deformation of the projectile and analytical method to represent the target.

Two-dimensional modal analyses have been performed to simulate ogive-nosed rod projectiles impacting limestone targets. Comparisons between experimental data and numerical results obtained with hexahedral

A c c e p t e d m a n u s c r i p t

continuum elements to represent the projectile and the same target modelling enable to check the accuracy of the two-dimensional modal approach. Moreover, thanks to the two-dimensional modelling and the efficiency of the stress integration decomposed in Fourier series, the two-dimensional modal approach proposed in this paper is able to reduce significantly the computation time in comparison to a classical three-dimensional analysis. On account of its cost-efficiency, the quasi axisymmetric modal approach may be well suited for parametric studies in order to provide useful insight into the effects of several variables such as impact angles, velocity, material and geometric characteristics of the projectile and target.

Although the formulation of the COMU elements has been established for thin shell structures, we showed that COMU elements still provide accurate predictions for penetration problems involving steel rod projectiles.

However, a limitation of the approach lies in the modal description of localisation process like damage or cracking: if a local damage appears (generated, for example, by a contact with a reinforced steel) the Fourier description with a few harmonics fails to reproduce the damage pattern. Use of higher order harmonics is still possible, but at a computational cost that makes COMU less attractive.

  radius, the circumferential angle and the vertical coordinate, respectively. At a point M, the local reference

  membrane strains and { } L χ the linear changes of curvature.

  s c r i p t where the 6-by-8 matrix [ ] 0 B corresponds to the divergence operator associated with the axisymmetric mode.

  being the circumferential angle. By introducing two 6-by-8 matrices [ ] and a sine factor, respectively, we can express the linear strains in the element in terms of the n symmetric displacement :

D A c c e p t e d m a n u s c r i p t

  beginning of the time step under consideration. At a time step of the explicit scheme, the incremental displacement is .

  q

.

  ) in which K and are constants obtained from data fits, is the reference projectile diameter, and is the projectile diameter under consideration. Three sets of experiments had been carried out with 7.1, 12.7 and 25.4 mm projectile diameters. The constant values found by the authors are: The density of limestone target is equal to 2130 kg/m this work, the three sets of penetration experiments published by Frew et al. are reproduced using a combined analytical and finite element method with COMU elements. All the VAR 4340 Rc 45 steel projectiles have a total length-to-diameter ratio of 10 and 3.0 caliber-radius-head (CRH) node shapes. The shank diameters and masses for each of the three sets of experiments are 25.4 mm, 0.931 kg (reference projectile); 12.7 mm, 0.117 kg; and 7.1 mm, 0.020 kg. So the target resistance parameters given by the relation (44) are 693 MPa, 779 MPa and 914 MPa with 25.4 mm, 12.7 mm and 7.1 mm projectile diameters, respectively. With this method accounting for the target resistance, all the constitutive behaviour of the target along with any frictional resistance is lumped into the previous expression (43).

,

  i being the layer index.

(

  at impact) to the target referential. After determining the Fourier coefficients of forces and , the Fourier coefficients of the acceleration contained in the global vector U are easily obtained by resolving the equation of motion with the diagonal mass matrix ext F int F M .

  ext

F

.

  compares the constitutive model of Warren with the elastic plastic model. Model parameters obtained from curve fits are: The density, Young's modulus, and Poisson's ratio of the VAR 4340 Rc=45 steel are taken to be ,

Figure 1 .Figure 2 .Figure 3 .Figure 4 .

 1234 Figure 1. 2D modelling of the projectile with COMU elements

Figure 5 .Figure 6 .Figure 7 .Figure 8 .Figure 9 .

 56789 Figure 5. Simulation of a 12.7mm-Diameter projectile penetrating a limestone target at 459 m/s and 3.55° angle of inclination; left side Warren's results with 3D elements [4] and right side results with COMU elements

  Figure 1.

Figure 3 .

 3 Figure 3.

  Figure 4.

  Figure 5.

  Figure 6.

  Figure 7.

Figure 8 .

 8 Figure 8.

Table 1

 1 summarizes penetration experiments and numerical results. There was no angle of obliquity with these experiments, so pitch and yaw angles are incorporated into a single angle of inclination. All of the COMU simulations are carried out with a modal basis containing 5 Fourier harmonics to represent displacements, velocities and accelerations. Thus the modal basis contains the axisymmetric mode, mode 1 symmetric and antisymmetric, and mode 2 symmetric and antisymmetric. That corresponds to

	N	=	2	in

Table 1 . Penetration results : depth of penetration P experimental data, numerical results with 3D element (Warren [4]) and with COMU elements; Relative difference between the numerical results.

 1

	Velocity (m/s)	Inclination	Experimental P SCE/PRONTOD	COMU 2D P	Relative
		(degree)		3D P (m)		difference
	Diameter projectiles (12.7 mm)
	459	3.55	0.141	0.137	0.137	0 %
	608	1.41	0.232	0.230	0.229	-0.5 %
	853	0.00	0.362	0.418	0.417	-0.2 %
	956	0.30	0.523	0.508	0.506	-0.4 %
	1134	2.15	0.562	0.608	0.645	+6 %
	1269	0.00	0.812	0.812	0.808	-0.5 %
	1404	0.70	0.924	0.931	0.945	+1.5 %
	1502	0.70	1.017	1.035	1.048	+1.2 %
	Diameter projectiles (25.4 mm)
	407	0.64	0.260	0.261	0.261	0 %
	566	0.60	0.390	0.454	0.452	-0.4 %
	800	0.61	0.790	0.826	0.820	-0.7 %
	917	0.30	1.020	1.045	1.036	-0.9 %
	1177	0.30	1.500	1.579	1.563	-1 %
	Diameter projectiles (7.1 mm)
	497	0.75	0.067	0.086	0.079	-8 %
	597	0.50	0.105	0.110	0.107	-2.7 %
	787	1.25	0.165	0.174	0.170	-2.3 %
	795	0.50	0.178	0.180	0.176	-2.2 %
	1037	0.90	0.271	0.280	0.274	-2.1 %
	1060	0.56	0.294	0.295	0.287	-2.7 %
	1230	0.25	0.392	0.379	0.368	-2.9 %
	1340	2.01	0.437	0.389	0.400	+2.8 %
	1365	2.56	0.430	0.369	0.383	+3.8 %
	1516	1.42	0.516	0.516	0.500	-3.1 %

ACKNOWLEDGEMENTS

This work, carried out by the applied research office EC2-MS (Etudes, Conseils et Calculs en Mécanique des Structures), was funded by contract from Centre d'Etudes de Gramat, Délégation Générale pour l'Armement.


B Castanié
email: bcastani@cict.fr

C Bouvet

Y Aminanda

J.-J Barrau

P Thevenet

Modelling of low energy/low velocity impact on Nomex honeycomb sandwich structures with metallic skins

Keywords: Sandwich structures, low velocity/low energy impact, finite element nonlinear analysis, indentation

In the aircraft industry, manufacturers have to decide quickly whether an impacted sandwich needs repairing or not. Certain computation tools exist at present but they are very time consuming and they also fail to perfectly model the physical phenomena involved in an impact. In a previous publication, the authors demonstrated the possibility of representing the Nomex TM honeycomb core by a grid of nonlinear springs and have pointed out both the structural behaviour of the honeycomb and the influence of core-skin boundary conditions. This discrete approach accurately predicts the static indentation on honeycomb core alone and the indentation on sandwich structure with metal skins supported on rigid flat support. In this study, the domain of validity of this approach is investigated. It is found that the approach is not valid for sharp projectiles on thin skins. In any case, the spring elements used to model the honeycomb cannot take into account the transverse shear that occurs in the core during the bending of a sandwich. To overcome this strong limitation, a multi-level approach is proposed in the present article. In this approach, the sandwich structure is modelled by Mindlin plate elements and the computed static contact law is implemented in a non linear spring located between the impactor and the structure. Thus, it is possible to predict the dynamic structural A c c e p t e d m a n u s c r i p t response in the case of low-velocity/low-energy impact on metal-skinned sandwich structures.

A good correlation with dynamic experimental tests is achieved.

1-Introduction

Sandwich structures exhibit static properties such as high stiffness-to-weight ratio and high buckling loads which are of great importance in the aeronautics field. Nevertheless, the current applications on commercial airplanes remain mainly limited to secondary structures like control surface or floor panels. In the field of helicopters where stress levels are lower, full sandwich structures are already in flight. In fact, one of the main limitations is linked to a lack of knowledge on the effects induced by impact damages [START_REF] Abrate | Localized impact on sandwich structure with laminated facings[END_REF]. However, in service, such structures are exposed and are often impacted during taxiing manoeuvres or take-off. They can also be damaged by tool drops during maintenance operations [START_REF] Guedra-Degeorges | Damage tolerance of sandwich structures[END_REF] and [START_REF] Razi | Analysis of sandwich panels with multiple-site damage[END_REF]. In such circumstances, the aircraft manufacturer has to inform the users very quickly of whether the impacted structure needs a repair or not. The study presented in this article is a step toward the final objective, which is to provide the aerospace industry with an inverse method. In this method, a 3D picture of the damage is first made by the airline company. Then, based on the shape of the impact, the impactor's shape is found and both the impact damage and the residual strength would be computed. This article investigates low speed and low energy impact phenomena. Tests in industry and most publications agree that, in this type of impact and for aeronautical sandwiches with

A c c e p t e d m a n u s c r i p t

skins that not exceed 2.5 mm, there is an equivalence between dynamic tests and static indentations [START_REF] Ferri | Static indentation and low velocity impact tests on sandwich plates[END_REF], [START_REF] Williamson | Response mechanism in the impact of graphite/epoxy honeycomb sandwich panels[END_REF] and [START_REF] Herup | Low-velocity impact damage initiation in graphite/epoxy/nomex honeycomb-sandwich plates[END_REF]. However, static tests have the advantage of being easier to set up and show low dispersion. This hypothesis will also be assumed in this study. At this first stage, the study will focus on modelling impacts on sandwich structures with Nomex TM honeycomb core and metallic skins. Metallic skins are used in order to avoid the complex mechanisms of laminate skins damaging.

The indentation response of sandwich structures can be determined by using Meyer's law with experimentally correlated coefficients, spring-mass models or energy-balance models [START_REF] Abrate | Localized impact on sandwich structure with laminated facings[END_REF]. Analytical models of beams or plates on elastic foundations have also been developed by several authors [START_REF] Swanson | Design of sandwich structures under contact loading[END_REF]. The first approaches were based on elasticity theory, so their practical use remains limited [START_REF] Abrate | Localized impact on sandwich structure with laminated facings[END_REF]. Soden [START_REF] Soden | Indentation of composite sandwich beams[END_REF] and Olsson [START_REF] Olsson | Simplified theory for contact indentation of sandwich panels[END_REF] used a foundation model with a perfect plastic law which may be able to represent the core behaviour. The model of Soden is also geometrically nonlinear. However, because the core crushes [START_REF] Bernard | Impact resistance of composite sandwich plates[END_REF] and [START_REF] Tsotsis | Characterization of localized failure modes in honeycomb sandwich panels using indentation[END_REF] and the skin plasticizes, this type of approach also remains limited. The damage of the core also proves to be very important as it even occurs at the weakest energies and without any visible damage on laminated skins.

Certainly, understanding and modelling the crushing phenomenon are the key points of the general problem of impact on sandwich structures with honeycomb. Core crushing is a complex mechanical phenomenon characterized by the appearance of various folds and failures in the hexagonal structure. This phenomenon is known for its energy-absorbing capacities and has been analysed since 1963 by Mc Farland [START_REF] Mcfarland | Hexagonal cell structures under post-buckling load[END_REF]. The crushing of aluminium honeycomb core has been extensively studied by Wierzbicki et Al [START_REF] Wierzbicki | Crushing analysis of metal honeycomb[END_REF] [14] and, more recently, by Mohr et Al [START_REF] Mohr | Deformation-induced folding systems in thin-walled monolithic hexagonal metallic honeycomb[END_REF]- [START_REF] Mohr | Experimental Investigation on the Plasticity of Hexagonal Aluminum Honeycomb Under Multiaxial Loading[END_REF]. Constitutive equivalent models have been developed [START_REF] Xue | Constitutive model for quasi-static deformation of metallic sandwich cores[END_REF], [START_REF] Wang | Yield surfaces of various periodic metal honeycombs at intermediate relative density[END_REF] and have been applied successfully to experiments on large structures subjected to blast A c c e p t e d m a n u s c r i p t loads. However, complex multiaxial experiments on honeycomb are required to identify the crushing laws and, as pointed out by Vaziri [START_REF] Mohr | Quasi-static punch indentation of a honeycomb sandwich plate: experiments and modelling[END_REF], the models are not suitable for simulating the phenomenon locally. Besides, the developments apply to metallic honeycomb only. It appears that Nomex TM honeycomb has a more complex micromechanical behaviour. Matrix cracking at the surface and local detachment and tears are observed [START_REF] Tsotsis | Characterization of localized failure modes in honeycomb sandwich panels using indentation[END_REF] and [START_REF] Aminanda | Etude expérimentale et modélisation du comportement en compression des structures nid d'abeille[END_REF].

Conventionally, modelling of this phenomenon is either based on an equivalent material law, which could be of the perfectly plastic type [START_REF] Thomson | Skin wrinkling of Impact Damaged Sandwich Composite[END_REF] and [START_REF] Besant | Finite element modelling of low velocity impact of composite sandwich panels[END_REF], or on a discrete crushing law (with a peak) for honeycomb blocks [START_REF] Zheng | Damage analysis of honeycomb core composite sandwich plate subjected low velocity impact[END_REF] and [START_REF] Minguet | A model for predicting the behaviour of impact-damaged minimum gage sandwich panels under compression[END_REF]. In this area, Horrigan [START_REF] Horrigan | Modelling of crushing due to impact in honeycomb sandwiches[END_REF] used an isotropic continuum damage model, but this modelling was limited to small indentations, since the continuity of the damage model and the plastic stress flow did not represent the real damage in the core.

All the approaches presented have a common feature which is to consider the honeycomb as a homogeneous material. Recently, based on a phenomenological approach and a structural analysis of honeycomb core with "soft" materials like paper or Nomex TM , the authors have proposed a discrete modelling of the honeycomb [START_REF] Aminanda | Experimental Analysis and Modeling of the Crushing of Honeycomb Core[END_REF]. A crushing law, including the influence of the boundary conditions which take into account the interaction between the skin and the honeycomb, was proposed. It predicts the static indentation of metal-skinned sandwiches supported on rigid flat support [START_REF] Aminanda | Modélisation de l'indentation des structures sandwichs à peaux métalliques[END_REF]. Those two points will be recalled in the first part of this article and an investigation of the domain of validity of this approach will be developed. It is noticeable that this approach can also predict the residual dent and the extension of the damage during compression after impact testing [START_REF] Aminanda | Experimental and Numerical Study of Compression After Impact of Sandwich Structure with Metallic Skins[END_REF]. An original criterion of structural strength in compression after impact (CAI) based on the collapse of cells at the edge of the impacted area was also proposed. The main limitation of this model is the fact that the vertical grid of nonlinear springs modelling the honeycomb core obviously cannot take

A c c e p t e d m a n u s c r i p t

into account the transverse shear stress. This limitation will be solved by demonstrating that the local indentation and the flexural response can be superimposed in static tests on a simplysupported sandwich structure. Then the method will be extended to dynamic tests that finally provide a multi-level approach.

2-Justification and limitation of the discrete modelling of honeycomb core

In this part, the key points of the work presented in Aminanda [START_REF] Aminanda | Experimental Analysis and Modeling of the Crushing of Honeycomb Core[END_REF], [START_REF] Aminanda | Modélisation de l'indentation des structures sandwichs à peaux métalliques[END_REF] are first summarized. Several uniform compression tests were conducted on honeycomb block with three different materials (paper, Nomex TM and aluminium). Each block of honeycomb was also bonded to one or two skins to simulate two types of boundary conditions between the honeycomb and the compression plate. Large cells were made for the paper honeycomb which allowed the mechanism of deformation to be observed easily. For a specimen with one skin, local buckling on one side of the free surfaces of the honeycomb occurred for very low loads while the vertical edge of the cell remained straight (Figure n°1). In the case of Nomex TM honeycomb, because of the small dimensions of the cells, it was not so easy to observe the mechanism. Nevertheless, a finite elements model showed that this local buckling occur at 20 % to 50 % of the peak or maximum load. The phenomenon seems to be triggered by the simply-supported boundary conditions on this free surface which allow local rotation of the cell walls. This sensitivity is also found on the force-displacement curves (Figure n°2).

For two-skin Nomex and Paper specimens the peak force was pronounced and the difference with the maximum force of one-skin specimens was significant. However, the difference was smaller for aluminium honeycomb, a reduction of about 10 % being observed. The Young modulus of aluminium is thirty times that of Nomex and the greater stiffness of the wall can

A c c e p t e d m a n u s c r i p t

partially explain the reduction. Therefore, it seems that the sensitivity to boundary conditions must be taken into account only in the case of "soft materials".

Since local buckling occurs early, the honeycomb structure works in a postbuckling mode and an analogy can be made with stiffened structures under compression in aircraft (Figure n°3). The stiffener in honeycomb is the vertical edge that is formed by the intersection of 3 cell wall. When the skin buckles, the compressive stress in the skin cannot be more than the buckling stress and the compression excess is therefore taken up by the stiffener and the nearby skin, with an equivalent half width equal to fifteen times the skin thickness (15.e)

[30]- [START_REF] Barrau | Buckling and post-buckling of beams with flat webs[END_REF]. By analogy, for the honeycomb cell, the compression is mainly taken by the vertical edges since the buckling of the walls occurs earlier. Then the collapse of the stiffened structure corresponds to the buckling of the stiffeners. In the case of honeycomb there is no collapse but, rather folding. The previous analysis shows that only the cell edge plays an important role from a structural point of view. This reasoning leads to the hypothesis that Nomex TM honeycombs under a crushing force behave like a juxtaposition of vertical cell edges and it is possible to model them by a grid of vertical nonlinear springs located at the angles of the hexagons (see Figure n°4). Furthermore, the local rotations must be taken into account as they modify the crush law of cell edge from a "with-peak" law to a "no peak" law. So, the following form is proposed for a generalized crush law of a cell edge:

Where w i and i are the vertical depth of crush and the local rotation on the upper part of cell edge i. w j , j = 1 to 3 are the crushing depth of the adjacent edges.

The function F Crush (w i , i ) is plotted in Figure n°5. The following variation is proposed:

• if i = 0 the experimental "with peak" F wp law is used (see the "two-skins" curve ).

A c c e p t e d m a n u s c r i p t

• if i >= critical then a "no peak" law F np is used.

• if 0 < i < critical then an intermediate law F i is defined.

The crush law for one edge "with peak", F wp is obtained by dividing the experimental compression law by the number of edges. The law presented was obtained using experiments on a Nomex TM honeycomb block (HRH78, 1/8-3 of 48kg/m 3 and 15 mm thick) containing 100 cells (10 on each side). The specimens were prepared by carefully cutting the specimen to preserve the vertical edges on the outer side of the honeycomb block (inset in Figure n°2).

The total number of cell edges of this block was 240.

The threshold parameter critical was identified by performing indentation tests on sandwich structures [START_REF] Aminanda | Modélisation de l'indentation des structures sandwichs à peaux métalliques[END_REF]. For the range of sandwich structures used in the experiments, this value was about 2.3°. H(w i -w j -w limit ) is a Heaviside function with:

The influence of F shear was investigated and w Limit was determined by means of tests using cylindrical indenters which created shear stress in the walls located in the circumferential zone of the indenter as explained by Wierzbicki et al [START_REF] Wierzbicky | Impact energy absorption of sandwich plates with crushable core[END_REF]. Quasi-static tests using 3 different radii (8, 14.75 and 25 mm) of cylindrical indenters were conducted on the same honeycomb (see Figure n°6). Two tests were performed for each radius but only single-test curves are presented because the behaviour is hardly varied. The circumference of the cylindrical indenter can be observed to coincide either with a vertical edge (POINT 1, Figure n°6) or with a wall (POINT 2, Figure n°6). According to the position on the circumference, it can be assumed that either the edge is subjected to compression or the wall is under shear load. Since it is assumed that these walls are subjected to shear load, in consequence the shear force is also proportional to the indenter radius. This additional shear force, F Sup can be expressed as:

F Sup = 2 r.q shear
where r is the radius of the cylindrical indenter and q shear is the shear force per unit length in the circumferential zone. The unit shear force q shear can thus be obtained by dividing the difference between the computation and the test by 2 r. Calculations show that the unit force q shear is not very different for the three indenter radii (Figure n°8). However, the peak force is lower for the 8 mm radius possibly because of the greater dispersion due to the location of the indenter and the small number of walls sheared. This result globally confirms the hypothesis that the extra indentation force is taken up by peripheral shearing and is in agreement with Wierzbicki's results [START_REF] Wierzbicky | Impact energy absorption of sandwich plates with crushable core[END_REF].

The results also demonstrate that the additional force, due to shear only, occurs from a certain indentation depth limit: w Limit =0.32 mm. Below this threshold, the error induced by assuming that the honeycomb can be modelled simply by its vertical edges, which behave independently of each other and are subjected only to compression load, seems to be negligible. A simple geometrical analysis (Figure n°9) shows that the corresponding limit angle is 0.27 rad, about 15°. This calculation was performed assuming

A c c e p t e d m a n u s c r i p t

that the circumference of the indenter passed through the middle of the cell, which avoided having to take the reaction of the vertical edge into account. It was then verified a posteriori that, in the case of the indentation experiments using spherical indenters [START_REF] Aminanda | Experimental Analysis and Modeling of the Crushing of Honeycomb Core[END_REF], this limit was not reached. The limit angle varied between 0.007 rad (radius 57.25mm) and 0.025 rad (radius 16.25mm) which is well below the calculated threshold. To see if the angle obtained had a real physical meaning, a complementary test was carried out with a conical indenter of halfangle 18° (> 15°) at the summit. It was observed that the edge-independence hypothesis still remained valid as the correlation between calculation and test results was very good (Figure n°10), which suggests that the actual threshold is little higher.

In consequence, for the indentation of a sandwich structure, it will be possible to model the honeycomb by its vertical edges alone if the deformation of the skin is not too pronounced and the difference in crushing between two neighbouring edges does not exceed w Limit , which is generally the case. For all our studies, w i -w j was always less than w limit . Therefore, we did not try to determine the F shear law because it is only needed for sharp projectiles on thin skins.

In such cases, the skin is generally perforated which is out of the scope of our research.

3-Computation of the static contact law

Experiments and numerical studies were carried out on sandwich structures with metal skin and the same Nomex TM honeycomb core as described in the last paragraph. Metal skins were used in order to avoid the complex failure damage mode of composite laminated skin.

Experimental investigations

Quasi-static indentation test were carried out on 100 x 100 mm sandwich structures.

Brass skins of 0.1 and 1 mm thickness were bonded to the Nomex TM honeycomb with a layer

A c c e p t e d m a n u s c r i p t 3.2 Comparison between computation and experiment.

An implicit nonlinear finite element model was made (see Figure n°14). The software used was the SAMCEF TM code (by SAMTECH GROUP). Nonlinear springs were placed at the same location as the honeycomb cell angles. The generalized law was implemented using special features of the software. The metal skins were modelled by Mindlin plate elements for thin skin (0.1mm) specimens and by volume elements in the thickness for thick skins (1 mm).

In this way, the triaxial stress state of the metallic skin located directly below the indenter was modelled satisfactorily. A fine mesh was created in the contact area just below the indenter to generate a smooth contact law without any slope discontinuity. A sensitivity study of the mesh gave a convergence result if 5 volume elements through the thickness and 36 elements per cell in the contact area were used.

The local rotation i corresponded to the rotation of the upper node of the spring which also belonged to the skin. This local rotation i was obtained directly with Mindlin plate elements but it had to be computed in the case of volume elements for the thick skin (more details can be found in [START_REF] Aminanda | Modélisation de l'indentation des structures sandwichs à peaux métalliques[END_REF]). In both cases, the rotation was assumed to be:

where R x and R y are the local rotations of the interface nodes between the skin and the spring. The elastic-plastic behaviour laws for brass skins were obtained from conventional tensile tests. Taking advantage of the symmetry of the structure, only one quarter of the plate was modelled. The numerical simulation was limited to 2 mm of indentation, which largely exceeds the threshold value of detectability known as BVID (Barely Visible Impact Damage).

The results of computation were compared to the tests on sandwiches with thin skins (Figure n°12) and a good test/computation comparison was obtained for all three indenter radii. Globally, the undulations observed during the test were also found numerically and cor-

A c c e p t e d m a n u s c r i p t

responded to the drop in the load after the peak load of each vertical edges located at the circumference of the indented area. For the thin skin, the rotation always proved to be less than 2.3°. It should therefore be possible to use the simple law with peak F wp for all vertical edges. This is hardly surprising since the skin shapes itself perfectly around the indenter and the test was similar to the one performed with the indenter directly in contact with the honeycomb.

However, as lateral folding of the edges could occur with the skin in place, the correlation was good whatever the diameter of the indenter. For the thick skins, if the F wp law was not corrected to take into account the rotation, there was a difference of about 15% between the computational and the experimental results. This can be explained by the fact that the rigid skin did not follow the shape of the indenter when it bent. The bend caused rotations at the menisci before the edges involved had reached peak load [START_REF] Aminanda | Modélisation de l'indentation des structures sandwichs à peaux métalliques[END_REF]. When the rotation is taken into account, the comparison for sandwich structures using thick skin gives a globally acceptable result considering the test dispersion (Figure n°13). For the thick skin, the bending stiffness was 1000 times that of the thin skin and the deformed shape of the skins generated local rotation in the upper edges of the honeycomb. Finally, the finite element model proposed was able to compute the static contact law with acceptable accuracy.

Indentation of a structure.

In practice, impacted aeronautical sandwich structures are mostly simply supported or clamped but are never fully supported. From a practical point of view, it is important to know whether or not a coupling effect exists, between the indentation and the stress field generated by the loading. To study this phenomenon, three-point bending tests were undertaken on sandwich plate.

A c c e p t e d m a n u s c r i p t

Three Nomex TM honeycomb specimen measuring 220 x 100 mm with 1 mm thick brass skins were produced (Table n°1). Globally, the experimental procedure and manufacturing of the specimen were the same as described previously (see Figure n°15). The distance between the cylindrical steel supports was 200 mm and three displacement sensors (DC3 to DC5) were located on the lower skin to measure the deflection of the specimen. Two displacement sensors (DC1 and DC2) were also located on the upper face of the indenter. Three tests were completed with indenter radii of 21.75 mm, 30.125 mm and 57.25 mm. The depth of indentation was obtained by calculating difference in displacement between the points on the upper and lower skins just below the indenter ((DC1 + DC2)/2 -DC4) in Figure n°15). The force/indentation curves could then be directly compared with those obtained previously from the tests on rigid flat supports (Figure n°16 et Figure n°17) and it was seen that the force/indentation contact laws were superimposed. Experimentally, for the configuration tested, it appears that there is no nonlinear coupling and, consequently, there is a superimposition of global bending and local indentation effects.

When a honeycomb sandwich structure undergoes bending, the transverse shear of the core plays an important part in the deflection and has to be taken into consideration. It is obvious that the proposed modelling of the core by a series of nonlinear vertical springs cannot take this transverse shear stress into account. However, the bending and indentation effects seem to be uncoupled. So, to represent the bending effect of the sandwich, equivalent nodal compression loads were added laterally at the nodes located at the edge of honeycomb block and at the indented skin side (Figure n°18). The computation result of this model was compared to the numerical simulations discussed in paragraph 3.4 and, as shown in Figure n°18, the results were equivalent. This approach enables the practical problem of sandwich structure indentation to be represented, since the bending/indentation uncoupling is again

A c c e p t e d m a n u s c r i p t

present. Physically, this uncoupling can be attributed to the existence of a plasticized area at a very early stage, under the indenter area. This area becomes saturated in stress and insensitive to the loading increments on the sandwich skin. The only possible coupling must appear at the initiation of the indentation area, which must logically be earlier when the skin is loaded under lateral compression. However, this phenomenon was not observable for the configuration tested.

At this stage of the study, it is now possible to model the static indentation of metalskinned sandwich structures. In the finite element model, a grid of nonlinear springs represents the honeycomb. The effects of local rotations of the skin during indentation on the crushing behaviour of the honeycomb are included. The compression law is obtained by a basic compression test on a block of honeycomb.

5-Dynamic study

In the previous paragraph, the possibility of obtaining the static contact law numerically was demonstrated. Nevertheless, in industrial cases the geometry can be more complex with various shapes and different local stiffnesses. So, the objective of this part of the study is to determine whether the approach developed in the last paragraph is suitable to model the dynamic behaviour of sandwich structures with metallic skins under low-velocity and low energy impact.

For this purpose, dynamic impact tests were performed using mass drop test equipment (see Figure n°19). The sandwich plates used for the experiments were of the same type as the ones for the three point bending tests (100x220 mm 2 , core thickness 15 mm, and skin thickness 1 mm). The boundary conditions were also the same (see Figure n°15). The impactor having spherical tip with radii of 30.125 or 57.25 mm, hit the centre of the plate.

A c c e p t e d m a n u s c r i p t

The mass of the impactor were respectively 885 et 865 g. The impact speeds recorded were respectively 2,58 and 2,80 m/s. The impact energy was about 3 Joules. The impactor was equipped with a load cell and an accelerator to provide the deflection and the force during the impact. The redundancy of these two data was voluntary. Nevertheless, practically, the force signal gave less interference and it was the only information that was used to measure the force. The force/time and force/displacement curves are shown Figure n°21 and Figure n°22.

The complex law of the dynamic contact force should a priori be the superimposition of the dynamic response of the sandwich structure and the local indentation.

To model the dynamic test, two assumptions were made:

• The global behaviour of the structure did not depend on local response during impact. This hypothesis agrees with the local nature of the impact reported in the literature.

• The static/dynamic equivalence was assumed for the range of structures and impacts studied. Thus, it was possible to use the static contact law computed previously.

With these assumptions, a multi-scale approach could be proposed: the local indentation law was computed first. The only parameter necessary was the crushing law for the spring, which was obtained simply by a basic compression test on a block of honeycomb. As the local dynamic effects were neglected, a nonlinear spring was used to represent the contact law in the global model as shown . The compression law for this spring was the indentation law previously computed in paragraph 3. This proposition is similar to the approach of Choi [START_REF] Choi | Low-velocity impact analysis of composite laminates using linearized contact law[END_REF] for laminated plates. An implicit finite element model was made using the same software as in paragraph 3.2. The sandwich structure was modelled by Mindlin plate finite elements (Figure n°20). The materials of the structure were assumed to be linearly elastic. The mass density of (implicit predictor-corrector scheme) was selected [START_REF]Samcef Users Guide[END_REF].

In spite of the relative simplicity of the model used, the numerical simulations fit the experimental results correctly (Figure n°21 and Figure n°22). The maximum contact force has been found, which will provide the damage area by a return to the indentation model. The static-dynamic equivalence for this range of structures and impacts, which was mainly observed experimentally, has also been confirmed numerically. Besides, the possibility of making indentations on fully supported specimens has been demonstrated. This result also globally validates the approach and the hypotheses made. The simplicity of the model should be an advantage in dealing with complex structures and multi-impact phenomena.

5-Conclusions

A method has been developed to model low velocity/low energy impacts on metalskinned sandwich structures, and gives good correlation of contact laws. An analysis of the The proposed multi-level approach consists of three steps. Firstly, a basic compression test must be performed on a block of honeycomb to obtain the initial crushing law. Secondly, using this law in nonlinear springs, it is possible to obtain the contact law by a finite element model and a nonlinear static computation. Finally, with the hypothesis of neglecting the dynamic effect at the contact, a basic finite element spring mass model using a nonlinear contact law computed at the last step is able to model a dynamic test. As a good correlation is obtained, the hypothesis of equivalence between static indentation and dynamic test is validated. This basic approach could be useful to model complex structures under impact or multi-impact. It is important to note again that the impact simulation is complete and is, finally based only on a simple economical compression test on a block of honeycomb core. It avoids the use of indentation tests on the complete structure to identify the Meyer's law coefficient. It also shows that the phenomenon remains local and, for the range of structures studied, is independent of the boundary conditions and the dimension of the plate. It would

Vertical displacement (mm) Indentation force (N)

C o m p u ta ti o n C o m p u t a t i o n C o m p u t a t i o n

Tests

FigureA c c e p t e d m a n u s c r i p t

  Figure n°7 compares the contact laws obtained from experiment and computation on

A c c e p t e d m a n u s c r i p t

  of REDUX TM 312/5 glue. Tension tests were performed on brass skin specimens for the two different thicknesses. For the 0.1 mm specimen, the elastic modulus was 103100 MPa and the yield stress was about 433 MPa. For the 1 mm specimen the elastic modulus was 70400 MPa and the yield stress was about 104 MPa. Obviously, alloys differ. Then, several specimens were made, as listed in . Indentation of all specimens was done using the same INSTRON TM machine and the same spherical indenters having different radii (57.25, 30.125 and 21.75 mm). Tests were performed with speed of 0.5 mm/minute, which can be considered as quasistatic loading. The specimens were fully supported on a rigid metallic part (Figuren°11). The loads were measured by the machine's sensor but displacements were measured using a dial comparator positioned on the indenter. Three tests were completed for each type of sandwich structure (with thin or thick skin) and for each indenter. The experimental contact laws are partially presented in Figuren°12, (0.1 mm skin) and (1 mm skin). The contact laws obtained show a low dispersion on the results in the case of 0.1 mm skins. This dispersion is higher for the 1 mm skin specimen (about 10%). For the 0.1 mm thin skin, the aspect of the crush curves shows nothing out of the ordinary, with the exception of small undulations which are similar to those observed on indentation on honeycomb alone. For the 1 mm thick skin, the curve shows a very high stiffness at the beginning of the indentation (because the bending stiffness was 1000 times that of the 0.1 mm specimens) which becomes weaker afterwards. This change can be qualitatively attributed to the collapse of the first honeycomb cell under the indenter. Indeed, the change of curve slope is observed for a penetration value (0.29 mm) which corresponds to the results obtained previously on the uniform compression test on the honeycomb alone. At that time, a noise corresponding to the breaking of the honeycomb was heard.

A c c e p t e d m a n u s c r i p t

  the brass skin and the Nomex were respectively 8000 kg/m 3 and 48 kg/m 3 . The transverse moduli of the core were G yz =44 MPa, G xz =30 MPa and E zz =120 MPa.Following the assumptions and to avoid local dynamic effects, it was decided to increase the transverse shear modulus of the core artificially on the right of the indenter with the following values: G yz = G xz =5000 MPa. Numerical tests were performed to demonstrate the little influence of the area dimension and the modulus of the local reinforcement on the global response of the structure. In our case, the nonlinear spring law was obtained numerically. The linear stiffness of the spring was 2.67 kN/mm, the yield stress corresponded to a force of 0.4 kN and the plastic stiffnesses were 0.714 and 0.93 kN/mm for the 30.125 and 57.25 indenters respectively. The initial velocities were those measured (2,58 and 2,80 m/s respectively). For the dynamic computation, a Hilbert Hughes Taylor algorithm with automatic time stepping

A c c e p t e d m a n u s c r i p t

  crushing from a structural point of view has enabled us to propose an original way of modelling Nomex TM honeycomb core using a grid of nonlinear springs. In practice, the springs in implicit finite elements provide a faster and a more robust computation, especially when the stiffness varies and decreases suddenly as is found to occur in brittle materials such as Nomex TM . The local rotation of the upper surface of the honeycomb which interfaces with the skin plays a role in the initiation of vertical edge buckling. This interface effect between skin and honeycomb is then taken into account in the model. Direct application of this modelling enables the contact law to be computed when metal-skinned sandwiches are quasi-statically indented on a flat support. This approach gives highly accurate correlation of indentation tests on a flat support or under 3-point bending.

Figure N° 7 :

 7 Figure N° 7: Experimental contact law with cylindrical indenter and discrete model comparison.

Figure N° 8 :Figure N° 9 :

 89 Figure N° 8: Peripheral shear force per unit length.

Figure N° 10 :Figure N° 11 :

 1011 Figure N° 10: Comparison between calculation and indentation tests with conical indenter (Angle: 18°).

Figure N° 12 :

 12 Figure N° 12: Test/numerical model correlation in the case of sandwich with thin metal skins.

Figure N° 13 :

 13 Figure N° 13: Test/numerical model correlation in the case of sandwich with thick metal skins.

Figure N° 14 :n l i n e a r s p r i n g s N o n l i n e a r s p r i n g s A c c e p t e d m a n u s c r i p t 9 Figure N° 15 :

 14915 Figure N° 14: Finite element modelling of the indentation problem on sandwich with thin metal skins.

Figure N° 16 :Figure N° 17 :Figure N° 18 :Figure N° 19 :Figure N° 20 :

 1617181920 Figure N° 16 : Comparison of the force/crush experimental curves in the case of rigid support and 3-points bending, indenter R = 21.75 mm.

Figure N° 21 :Figure N° 22 :

 2122 Figure N° 21: Dynamic load-displacement law (a) impactor 30.125 mm, (b) impactor 57.25 mm.

A c c e p t e d m a n u s c r i p t

seem that the approach that has been developed can be used for other cellular cores made from other materials thanks to the similar and common crushing mechanism.

The use of metallic skin in this study enables a step by step approach to the modelling of the impact but this remains a limitation in practice since it is rare nowadays for sandwich structures to be made using metallic skin in industry. The real challenge will still be to couple this modelling approach with laminated skin. In this case, the plasticity behaviour will be replaced by determination of the damage state in both the honeycomb core and the skin and the possible coupling between the stress state of the global structure and the indentation phenomenon will have to be taken into account. Knowing the local state, it will be then possible to compute the residual strength by a second model.

6-References

A c c e p t e d m a n u s c r i p t A c c e p t e d m a n u s c r i p t A c c e p t e d m a n u s c r i p t

Tables.

.

SANDWICH MATERIALS

FIGURES.

Figure