Pascal Aubertin

Julien Réthoré

René De Borst

A MOLECULAR DYNAMICS / EXTENDED FINITE ELEMENT METHOD FOR DYNAMIC CRACK PROPAGATION

Keywords: multiscale methods, molecular dynamics, extended finite element method, fracture, crack propagation

A multiscale method is presented which couples a molecular dynamics approach for describing fracture at the crack tip with an extended finite element method for discretizing the remainder of the domain. After recalling the basic equations of molecular dynamics and continuum mechanics the discretization is discussed for the continuum subdomain where the partition-of-unity property of finite element shape functions is used, since in this fashion the crack in the wake of its tip is naturally modelled as a traction-free discontinuity. Next, the zonal coupling method between the atomistic and continuum models is described, including an assessment of the energy transfer between both domains for a one-dimensional problem. Finally, a two-dimensional computation is presented of dynamic fracture using the coupled model.

INTRODUCTION

Modern research into fracture commences with the seminal work of Griffith [START_REF] Griffith | The phenomena of rupture and flow in solids[END_REF]. Later, Irwin [START_REF] Irwin | Analysis of stresses and strains near the end of a crack traversing a plate[END_REF] and Rice [START_REF] Rice | A path independent integral and the approximate analysis of strain concentration by notches and cracks[END_REF] established the relation between the stress intensity factors and the energy release rate, and gave linear elastic fracture mechanics a firm basis. However, linear elastic fracture mechanics only applies to crack-like flaws in an otherwise linear elastic solid and when the singularity associated with that flaw is characterized by a nonvanishing energy release rate. The fracture and any dissipative processes must also be confined to a small region in the vicinity of the crack tip.

When the region in which the separation and dissipative process take place is not small compared to a structural dimension, but any nonlinearity is confined to a surface emanating from a classical crack tip, i.e. one with a non-vanishing energy release rate, cohesive zone models as introduced by Barenblatt [4] and Dugdale [START_REF] Dugdale | Yielding of steel sheets containing slits[END_REF] apply. The cohesive zone approach was extended by Hillerborg et al. [START_REF] Hillerborg | Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements[END_REF] and Needleman [START_REF] Needleman | A continuum model for void nucleation by inclusion debonding[END_REF] to circumstances where: (i) an initial crack-like flaw need not be present or, if one is present, it need not be associated with a non-vanishing energy release rate; and (ii) non-linear deformation behaviour may occur over an extended volume. Initially, cohesive-zone models were incorporate in finite element methods via special-purpose interface elements [START_REF] Rots | Smeared and discrete representations of localized fracture[END_REF][START_REF] Xu | Numerical simulations of fast crack growth in brittle solids[END_REF], but more recently, partition-of-unity finite element methods have shown to very amenable to the incorporation of cohesive-zone models, e.g. [START_REF] Wells | Discontinuous analysis of softening solids under impact loading[END_REF][START_REF] Wells | A consistent geometrically non-linear approach for delamination[END_REF][START_REF] De Borst | Numerical aspects of cohesive zone models[END_REF]. In particular, they naturally enable crack propagation, also in dynamics [START_REF] Duarte | A generalized finite element method for the simulation of three-dimensional crack propagation[END_REF][START_REF] Réthoré | An energy conserving scheme for dynamic crack growth with the extended finite element method[END_REF][START_REF] Menouillard | Efficient explicit time stepping for the extended finite element method[END_REF][START_REF] Remmers | A cohesive segments method for the simulation of crack growth[END_REF][START_REF] Remmers | The simulation of dynamic crack propagation using the cohesive segments method[END_REF] and in multi-phase continua [START_REF] Réthoré | A two-scale approach for fluid flow in fractured porous media[END_REF][START_REF] Réthoré | A two-scale model for fluid flow in an unsaturated porous medium with cohesive cracks[END_REF].

In spite of the power of the cohesive-zone approach, and its wide applicability on a range of scales, it remains a phenomenological approach. Probably, quantum mechanics is physically the most appropriate theory to describe fracture, but the difficulties to relate quantum mechanics to continuum mechanics, e.g. via Density Functional Theory [START_REF] Parr | Local density functional theory of atoms and molecules[END_REF][START_REF] Car | Unified approach for molecular dynamics and density-functional theory[END_REF] presently seem insurmountable. One scale of observation higher is to use Molecular Dynamics to describe fracture processes from a more fundamental physics point of view. Indeed, researchers have recently used this approach to describe fracture, e.g. [START_REF] Abraham | Simulating materials failure by using up to one billion atoms and the world's fastest computer: work-hardening[END_REF][START_REF] Zhou | Three-dimensional fracture via largescale molecular dynamics[END_REF][START_REF] Miller | Quasicontinuum models of fracture and plasticity[END_REF]. A disadvantage of the approach is that it is computationally demanding. For this reason multi-scale approaches have been introduced, in fracture [START_REF] Kolhoff | Crack propagation in bcc crystals studied with a combined finite element and atomistic model[END_REF], as well as in plasticity [START_REF] Shilkrot | Coupled atomistic and discrete dislocation plasticity[END_REF], in which only a part of the body is analysed using molecular dynamics, while the remaining part of the body is modelled using continuum mechanics and discretized using a finite element method. This manuscript furthers along this line and combines molecular dynamics for modelling the fracture process at the crack tip with an extended finite element method (XFEM), where the partition-of-unity property of the polynomial shape functions is exploited to model the crack in the wake of the tip as a traction-free discontinuity. It is noted that recently another approach has been published that couples atomistics and extended finite elements [START_REF] Gracie | Concurrently coupled atomistic and XFEM models for dislocations and cracks[END_REF], but the current paper makes a further advancement in that it includes dynamic loadings.

A major issue in multi-scale approaches as discussed above is the accurate coupling of both domains, especially when different descriptions are assumed on either domain. While the coupling can, in principle, either be achieved at a discrete interface, or on a zone of a finite size (overlap or zonal coupling), it is believed that zonal approaches, which include the Arlequin method [START_REF] Dhia | The Arlequin method as a flexible engineering design tool[END_REF], the bridging domain method [START_REF] Xiao | A bridging domain method for coupling continua with molecular dynamics[END_REF][START_REF] Zhang | A bridging domain and strain computation method for coupled atomistic-continuum modelling of solids[END_REF][START_REF] Guidault | On the L2 and the H1 couplings for an overlapping domain decomposition method using Lagrange multipliers[END_REF][START_REF] Xu | Conservation properties of the bridging domain method fod coupled molecular/continuum dynamics[END_REF], discreteto-continuum bridging [START_REF] Fish | Discrete-to-continuum bridging based on multigrid principles[END_REF], the discontinuous enrichment method [START_REF] Farhat | The discontinuous enrichment method for multiscale analysis[END_REF], and bridging scale decomposition [START_REF] Wagner | Coupling of atomistic and continuum simulations using a bridging scale decomposition[END_REF][START_REF] Farrell | Implementation aspects of the bridging scale method and application to intersonic crack propagation[END_REF] enable a more gradual transition from one domain to the other. The ability of a gradual transition is especially important for highly dissimilar domains and when wave propagation phenomena are considered, where preservation of the energy and avoiding spurious reflections when a wave exits one domain and enters the other can become an issue. Inspired by earlier work by Ben Dhia and Rateau [START_REF] Dhia | The Arlequin method as a flexible engineering design tool[END_REF] and Xiao and Belytschko [START_REF] Xiao | A bridging domain method for coupling continua with molecular dynamics[END_REF] we have chosen a weak coupling between the models in the two adjacent domains.

This paper is organized as follows. First, we briefly list the equations of molecular dynamics. This is followed by a succinct recapitulating of the governing equations of continuum mechanics, both in the strong and the weak forms. The discretization of the continuum subdomain is carried out using the extended finite element method, where the partition-of-unity property of finite element shape functions is used to model the tractionfree discontinuity in the wake of the crack tip. Subsequently, it is discussed how both domains can be coupled, see also [START_REF] Aubertin | Energy conservation of atomistic/continuum coupling[END_REF] and an analysis is presented of the energy conservation properties of the coupling scheme. The paper concludes with a full two-dimensional coupled analysis of dynamic crack propagation which shows multiple branching and suggests the formation of dislocations, which shows similarities with recent simulations on dynamic fracture using cohesive-zone models [START_REF] Xu | Numerical simulations of fast crack growth in brittle solids[END_REF][START_REF] Remmers | The simulation of dynamic crack propagation using the cohesive segments method[END_REF].

MOLECULAR DYNAMICS

For the discrete domain, i.e. Ω m , we build a grid of N a atoms, and, accordingly, the initial value problem in this domain can be written as:

For 1 ≤ i ≤ N a (t) and t ∈ [0; T ], given the initial conditions d (0) , ḋ (0) ,f i n d (d, f ) ∈D ad ×F ad such that:

m i di = f i (1)
with m i the mass of atom i and:

D ad = d = d i (t) 1≤i≤N a , ∀t ∈ [0,T] F ad = f = f i (t)=-∇ ∇ ∇ i U d (t) 1≤i≤N a , ∀t ∈ [0,T] (2)
from where it transpires that the interatomic forces are derived from a potential energy U . d and f assemble the discrete displacements d i and forces f i of the individual atoms, respectively. The internal energy of the discrete domain can be viewed as the sum of each atomic contribution U j :

U = j U j (d) (3)
and the force f i acting on atom i can be written as the sum of all elementary forces:

f i = - ∂U ∂d i = j =i f ij (4)
In order to limit the cost of computing such a force, we reduce the summation by only including so-called "nearest" neighbours, within a cut-off radius r c :

f i ≃ r ij <r c f ij (5)
where r ij is the interatomic distance.

In the subdomain Ω m the weak formulation becomes:

∀w * ∈ Ḋad,0 , given the initial conditions d (0) , ḋ (0) , find d ∈D ad such that: a m (d, w * )=0 (6)
with w * the test function, and

a m (d, w * )= N a i=1 m i di • w * i + N a i=1 ∇ ∇ ∇ i U (d) • w * i 3. CONTINUUM MODEL
Assuming small deformation gradients for simplicity, the governing equations in the continuum subdomain Ω M can be written in a standard manner as: For x ∈ Ω M (t) and t ∈ [0; T ] , given the initial conditions u (x, 0) , u (x, 0) , find (u,σ σ σ) ∈U ad ×S ad such that:

ρü =divσ σ σ + g d (7)
with ρ the mass density and u the continuum displacement vector, σ σ σ the stress tensor, and g d the body force vector applied in Ω M , subject to the boundary conditions

U ad = u = u (x,t) ∈ H 1 (Ω M ) 3 ; u = u d on ∂ 1 Ω , ∀t ∈ [0,T] S ad = σ σ σ = K : ∇u (x,t) ∈ L 2 (Ω M ) 6 ; σ σ σ • n = F d on ∂ 2 Ω , ∀t ∈ [0,T] (8)
where K is the fourth-order stiffness tensor, n the outward normal vector to ∂Ω 2 ,a n du d and F d the prescribed displacements and tractions at ∂Ω 1 and ∂Ω 2 , respectively.

To allow for a discretization of the continuum subdomain we next specify the weak formulation:

∀v * ∈ U ad,0 , given the initial conditions u (x, 0) , u (x, 0) , find u ∈U ad such that:

a M (u, v * )=l M (v * ) (9)
with v * the test function, and

a M (u, v * )= Ω M ρü • v * dΩ + Ω M ǫ ǫ ǫ (u):K : ǫ ǫ ǫ (v * )dΩ (10) l M (v * )= ∂ 2 Ω F d • v * dS + Ω M g d • v * dΩ (11)
4. COUPLING SCHEME

Coupling functions

In order to enable an efficient coupling between the two domains, a coupling zone is defined and a coupling function is computed which is used to obtain the global energy. First a coupling length L c is chosen such that will be the characteristic length of the coupling region. Subsequently, the patch of atoms, Ω m , is included in the continuum at a given position. Atoms at a distance ℓ ≤ L c from the Molecular Dynamics Box (MD-Box) boundary are considered to be in the coupling zone Ω c . The finite elements in this zone are named "coupling elements". Inside the MD-Box, where only the atomistic model applies, elements are removed. The resulting, discretized domain is shown in Figure 1. In principle, each atom must be coupled to the finite element discretization of the underlying continuum. It suffices to couple only a limited number of atoms to the underlying continuum without loss of accuracy. Since this directly affects the size of the coupling matrices that Fig. 1. Discretized domain with the coupling region. The dark, "pear-shaped" area is the domain where an MD calculation is carried out (Ω m ). The coupling region Ω c consists of elements that are surrounded by a bold line will be derived next, such a limited coupling markedly decreases the computational effort, and therefore, a bigger scale transition becomes possible at reasonable computational costs.

In the coupling zone Ω c , we enforce a velocity coupling in a weak sense, and we distribute the energy between both models via a partition of unity [START_REF] Dhia | The Arlequin method as a flexible engineering design tool[END_REF]. For this purpose we define the following functions, see Figure 2:

α :Ω M → [0, 1] (12)
β :Ω m → [0, 1] (13
) such that:

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ α(x)=1 for x ∈ Ω M \Ω c β(x)=1 for x ∈ Ω m \Ω c α(x)+β(x)=1 for x ∈ Ω c (14)
The displacement and velocity fields in the domains Ω M and Ω m have a different nature. In Ω M we have a continuum field, while in Ω m we have a discrete field, which is only defined at the geometrical points corresponding to the atoms. Therefore we construct a so-called "mediator space", denoted by M, on which we project the fields u and ḋ in order to be able to compare them. The nature of M is constrained by the discrete character of the atomistic field. Indeed, its displacements cannot be extrapolated outside the atomic positions if we want to maintain a physical interpretation at the fine scale. Accordingly, M has to be a subspace of the physical atomistic space. More precisely, we project the velocities using an operator Π on a discrete subset Ω c of the atomic positions included in

∀μ * ∈M ,c μ * , Π u -Π ḋ =<μ * , Π u -Π ḋ > M ( 15
)
with c the classical scalar product on M. The global equations are coupled via Lagrange multipliers and can subsequently be written as:

∀ (v * , w * ,μ μ μ * ) ∈ U ad,0 × Ḋad,0 ×M, given the initial conditions u (x, 0) , u (x, 0) , d (0) , ḋ (0) , find (u, d,λ λ λ) ∈U ad ×D ad ×M such that: a α,M (u, v * )+a β,m (d, w * )+c (λ λ λ, Πv * -Πw * )+c μ μ μ * , Π u -Π ḋ = l α,M (v * ) (16)
The modified forms a α,M , a β,m and l α,M take into account the weighting functions α(x) and β(x), see [START_REF] Aubertin | Energy conservation of atomistic/continuum coupling[END_REF] for details. In the atomistic subdomain the modified form a β,m that takes into account the distribution of the energy reads:

a β,m (d, w * )=w * • m β d -w * • f β (17)
with:

β i = β (d i ) , m β =[β i δ i,j m i ] , f β =[f β,i ] (18)

Discretized problem

In a manner which is by now standard the interpolation of each component of the displacement field is enriched with discontinuous functions in order to properly capture the traction-free discontinuity in the wake of the crack tip:

∀x ∈ Ω M ,u h (x)= i∈N M N i (x)ū i + i∈N cut N i (x)H Γ d ûi (19)
where N i are standard finite element shape functions supported by the set of nodes N M included in the discretized domain Ω M .N o d e si nN cut have their support cut by the discontinuity. They hold additional degrees of freedom ûi corresponding to the discontinuous function H Γ d defined by:

H Γ d (x)= x • n Γ d x • n Γ d (20)
with n Γ d the normal to the discontinuity Γ. Symbolically, eq. ( 19) can be written as

∀x ∈ Ω M , u h = N T U (21)
where the matrix N contains the standard interpolation polynomials N i (x) as well as the discontinuous function H Γ d , and the array U contains the displacement degrees-of-freedom ūi and ûi . The transition within the domain Ω M between the subdomain where the nodes are "enriched" and the part which has just the standard formulation does not affect the Molecular Dynamics computation other than through the coupling matrices.

With the latter symbolic notation the bilinear form a α,M and the linear form l α,M become:

a α,M (u h , v * h )=V * T M α Ü + V * T K α U (22) l α,M (v * h )=V * T F α ( 23
)
where the term that represents the body forces has been omitted for simplicity, and

M = Ω M ρN T NdΩ (24) K = Ω M ∇N T K∇NdΩ (25)
the mass and stiffness matrices, respectively. With the standard definition of the scalar product, the coupling term in the continuum can be discretized as follows:

c (λ λ λ, Πv * h )=V * T C M Λ Λ Λ=V * T F L M ( 26
)
with C M the continuum coupling matrix. The vector Λ Λ Λ contains the Lagrange multipliers and its size equals the Ω c subset cardinal times the dimension of the space considered. F L M can be regarded as a fictitious force due to the coupling via the Lagrange multipliers. This force has a non-zero value only in the coupling zone Ω c .

the Lagrange multipliers in the atomistic domain similar to that in the continuum domain:

f L m = C m Λ Λ Λ (27)
the weighted and coupled system ( 16) can be cast in a matrix-vector format:

V * T M α Ü + K α U + C M Λ Λ Λ + W * T m β d -f β -C m Λ Λ Λ + μ μ μ * T C T M U -C T m ḋ = V * T F α ( 28
)
Since this set must hold for any admissible (V * , W * ,μ μ μ * ) we finally obtain:

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ M α Ü + K α U = F α -C M Λ Λ Λ m β d = f β + C m Λ Λ Λ C T M U = C T m ḋ (29)
with (U, d, Λ Λ Λ) the set of unknowns. Details on the time integration scheme associated with this set of coupled ordinary differential equations are given in Ref. [START_REF] Aubertin | Energy conservation of atomistic/continuum coupling[END_REF]. Fig. 3. Mechanical energy (drawn line) and the work stored in the Lagrange multipliers (dashed line) when the prediction step is with the full stiffness and mass matrices.

Energy transfer between the atomistics and continuum domains

The energy transfer for the coupling between the continuum and atomistic domains is studied by means of a longitudinal bar discretized with finite elements, and containing an atomistic region, with a coupling zone on both sides. The bar is submitted to a traction wave, which is enforced by displacing the left-most 20 elements in the initial configuration. The right-hand end is free. The whole domain is 59.142528 × 10 -9 m long and 100 atoms have been put in the atomistic domain. The interatomic distance is r e =0 .1234708 × 10 -9 m, and the finite element size is h = r e . We use a Lennard-Jones potential as constitutive model for the atoms, with a =3 2 .043529 × 10 -21 J and a mass m =0.0016599 × 10 -24 g. The elastic material properties for the finite element model have been derived from the atomic properties [START_REF] Aubertin | Coupling of Atomistic and Continuum Models: Dynamic Crack Propagation[END_REF]. The computation continues for 2000 time steps with ∆t =1× 10 -15 s, which amounts to 95% of the critical time step.

Figure 3 shows the total mechanical energy during the computation. We observe some fluctuations each time the wave crosses a coupling zone. In fact, work is stored by the Lagrange multipliers, and subsequently put back in the mechanical system when the wave exits the coupling zone. We observe that the work of the Lagrange multipliers is complementary to the mechanical energy, and the energy balance is therefore satisfied. From energy plots of Figure 4 we observe that for different coupling lengths (i) the total energy is preserved, and that (ii) the energy correctly passes from one domain to the other when the wave traverses the coupling zones.

DYNAMIC FRACTURE

We now proceed with a two-dimensional simulation of dynamic fracture. A copper single crystal is considered in its (111) plane, so that the two-dimensional lattice is hexagonal. The Lennard-Jones potential is used in the molecular dynamics simulation with parameters from [START_REF] Agrawal | Predicting trends in rate parameters for self-diffusion on FCC metal surfaces[END_REF] ), which corresponds to a mass density ρ = 1865.250812586 × 10 3 kgm -3 . In the present study, the temperature has not been taken into account, since the focus is on the coupling of Molecular Dynamics to an (extended) finite element method for crack propagation. When extending the methodology to explicitly include the temperature a "thermal equilibrium" has to be achieved in addition to the mechanical equilibrium. This can for instance be done using the Nose-Hoover thermostat method [START_REF] Nose | Constant-temperature molecular dynamics[END_REF].

The computational domain is 100 nm long and 77.5 nm wide with an initial crack of 10 nm. The finite element mesh consists of 1221 quadrilateral elements and 4868 nodes. The element size is about 10 times the interatomic distance. 36635 atoms are put in the vicinity of the crack tip. The width of the coupling domain is approximatively 3 nm and 33% of the atoms in this region hold Lagrange multipliers. For computational reasons, the results that are presented from now on have been obtained by only including the first neighbours in the atomistic interactions. Indeed, the equilibration techniques and the updates are expensive when we take into account many neighbours. Simulations on a smaller scale have indicated that the results are close to those obtained with more neighbours.

The test consists of applying a velocity on the top and bottom edges of the specimen. The prescribed velocity in this example is V p =47.4 ms -1 .T h et i m es t e pi s15.811388 fs (10 -15 s). Figure 5 gives the displacement field after the crack has propagated, at t = 33.209 ps. Close inspection reveals that, going from the atomistic zone to the continuum domain, no spurious wave reflections occur, which indicates that the coupling algorithm works properly, and is in agreement with results for linear elasticity presented in [START_REF] Aubertin | Energy conservation of atomistic/continuum coupling[END_REF]. We also see that, even though the prescribed loading corresponds to Mode-I and the atomistic lattice is a perfect single crystal, the crack path is not that which we would expect to obtain with classical continuum methods. We observe crack branching, the occurrence of dislocations, and locally, mixed-mode behaviour.

CONCLUDING REMARKS

A numerical approach has been proposed for combining a molecular dynamics method and a finite element method that exploits the partition-of-unity property of finite element shape functions (extended finite element method). The aim is to simulate dynamic fracture in an efficient manner on basis of elementary physical principles. To this end the zone around the crack tip is modelled using molecular dynamics. Around this so-called Molecular Dynamics Box a continuum mechanics approach is adopted, with the finite element method used for discretization. The partition-of-unity property of the finite element shape functions is exploited to model the crack in the wake of its tip as a traction-free discontinuity. The coupling between the continuum and molecular dynamics zones has a zonal character where the energy is partitioned over both models and a weak velocity coupling is enforced. In this manner, spurious reflections are avoided and energy is conserved when a wave travels from one zone into another.

Fig. 2 .

 2 Fig. 2. Partition of unity for the energy distribution

Fig. 4 .

 4 Fig.[START_REF] Barenblatt | The mathematical theory of equilibrium cracks in brittle fracture[END_REF]. Energy plots for the finite element -Molecular Dynamics coupling. In subfigures (a), (c) and (e) the drawn line is the energy in the first (left) continuum, the dashed line is the energy in the second (right) continuum, and the dash-dotted line represents the energy stored in the atoms. The bold drawn line is the total mechanical energy.

Fig. 5 .

 5 Fig. 5. Deformed configuration at t =33.209 ps

  : a =0 .415 eV and b =0 .2277 nm. The Young's modulus E

	Mechanical energy (zJ)				Total energy drift
	0.012						3.0e-11
							2.5e-11
	0.010
							2.0e-11
	0.008						1.5e-11
							1.0e-11
	0.006
							5.0e-12
	0.004						0.0e+00
							-5.0e-12
	0.002
							-1.0e-11
	0.000						-1.5e-11
	0.00	0.02	0.04	0.06	0.08	0.10	0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20
			Time (ps)			Time (ps)
	(a) Energy transfer for L c consisting	(b) Energy drift for L c consisting of
	of 2 elements					2e l e m e n t s
	Mechanical energy (zJ)				Total energy drift
	0.012						1e-11
	0.010						0e+00
	0.008						-1e-11
	0.006						-2e-11
	0.004						-3e-11
	0.002						-4e-11
	0.000						-5e-11
	0.00	0.02	0.04	0.06	0.08	0.10	0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20
			Time (ps)			Time (ps)
	(c) Energy transfer for L c consisting	(d) Energy drift for L c consisting of
	of 5 elements					5e l e m e n t s
	Mechanical energy (zJ)
	0.012
	0.010
	0.008
	0.006
	0.004
	0.002
	0.000
	0.00	0.02	0.04	0.06	0.08	0.10
			Time (ps)


Kai Wei
email: weikai@hnu.edu.cn

Fuhao Mo
email: fuhaomo@hnu.edu.cn

Ming Mei

Yujia He

Xujing Yang

2-2 twill woven fabric under shear pre-tension coupling

Keywords: Fabrics, Composite, Picture frame, Shear-tension coupling, Forming defect

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

Introduction

Due to the high strength/stiffness, low density and high manufacture flexibility, continuous fiber fabric reinforced resin matrix composites are widely used in vehicle engineering, especially in new energy vehicles for weight saving consideration [START_REF] Liu | Lightweight design of carbon twill weave fabric composite body structure for electric vehicle[END_REF][START_REF] Mo | Coupling investigation on tensile and acoustic absorption properties of lightweight porous LGF/PP composite[END_REF][START_REF] Wei | A viscoelastic model of compression and relaxation behaviors in preforming process for carbon fiber fabrics with binder[END_REF][START_REF] Wei | Mechanical analysis and modeling of metallic lattice sandwich additively fabricated by selective laser melting[END_REF].

It is well known that the preforming of the continuous fiber fabric is an important step in manufacture process of the composites, and the defects such as wrinkle in fabrics generated in the manufacture process have significant negative influence on the overall mechanical properties [START_REF] Do | Effects of the wrinkled fabric preform on the bending properties of resin transfer molded composites[END_REF][START_REF] Wei | Mechanical properties of Invar 36 alloy additively manufactured by selective laser melting[END_REF].

As indicated in reported literature, tensile, in-plane shear and bending are the primary deformation modes of the continuous fiber fabric in the forming process. Specially, in-plane excessive shearing, which is mainly related to the initiation of wrinkles defect [START_REF] Boisse | Simulation of wrinkling during textile composite reinforcement forming. Influence of tensile, in-plane shear and bending stiffnesses[END_REF][START_REF] Tephany | Development of an experimental bench to reproduce the tow buckling defect appearing during the complex shape forming of structural flax based woven composite reinforcements[END_REF][START_REF] Wang | Experimental and numerical analyses of manufacturing process of a composite square box part: Comparison between textile reinforcement forming and surface 3D weaving[END_REF][START_REF] Boisse | The bias-extension test for the analysis of in-plane shear properties of textile composite reinforcements and prepregs: a review[END_REF][START_REF] Zhu | An experimental study of in-plane large shear deformation of woven fabric composite[END_REF][START_REF] Hosseini | Identifying the distinct shear wrinkling behavior of woven composite preforms under bias extension and picture frame tests[END_REF][START_REF] Guzman-Maldonado | Experimental and numerical analysis of wrinkling during forming of multi-layered textile composites[END_REF], and significantly affects the fabric permeability [START_REF] Kim | Prediction of the vacuum assisted resin transfer molding (VARTM) process considering the directional permeability of sheared woven fabric[END_REF][START_REF] Pierce | Permeability characterization of sheared carbon fiber textile preform[END_REF] in liquid composite molding (LCM) process, is considered as the most dominant deformation mode [START_REF] Gereke | Experimental and computational composite textile reinforcement forming: A review[END_REF][START_REF] Bussetta | Numerical forming of continuous fibre reinforced composite material: A review[END_REF][START_REF] Krogh | Picture-frame testing of woven prepreg fabric: An investigation of sample geometry and shear angle acquisition[END_REF][START_REF] Cao | Characterization of mechanical behavior of woven fabrics: Experimental methods and benchmark results[END_REF]. Thus, in order to explore the shear behavior of the fabrics, the picture frame and the bias extension tests are widely reported in literature [START_REF] Boisse | The bias-extension test for the analysis of in-plane shear properties of textile composite reinforcements and prepregs: a review[END_REF][START_REF] Zhu | An experimental study of in-plane large shear deformation of woven fabric composite[END_REF][START_REF] Hosseini | Identifying the distinct shear wrinkling behavior of woven composite preforms under bias extension and picture frame tests[END_REF]. Fig. 1 illustrates the basic principles of these two testing protocols. Compared with the bias extension test, the superiority of the picture frame test is that it can provide relative more uniform shear deformation field without slippage of the yarns [START_REF] Kashani | The Role of Intra-Yarn Shear in Integrated Multi-Scale Deformation Analyses of Woven Fabrics: A Critical Review[END_REF][START_REF] Montazerian | A frameless picture frame test with embedded sensor: Mitigation of imperfections in shear characterization of woven fabrics[END_REF].

Besides the exploration of the shear experiments, actually, a number of works reported that applying pre-tension to fabrics through such as blank holder [START_REF] Rashidi | Passive control of wrinkles in woven fabric preforms using a geometrical modification of blank holders[END_REF][START_REF] Matveev | A novel criterion for the prediction of meso-scale defects in textile preforming[END_REF][START_REF] Tavana | Meso/macro-scale finite element model for forming process of woven fabric reinforcements[END_REF][START_REF] Nezami | Analyses of interaction mechanisms during forming of multilayer carbon woven fabrics for composite applications[END_REF]and spring-load clamps [START_REF] Chen | Formability optimisation of fabric preforms by controlling material draw-in through in-plane constraints[END_REF][START_REF] Harrison | Press forming a 0/90 cross-ply advanced thermoplastic composite using the double-dome benchmark geometry[END_REF] was an effective way to suppress wrinkles caused by excessive shearing. Thus, recently, several studies discussed the shear pre-tension coupling properties of the fabrics through the picture frame testing devices [START_REF] Willems | Picture frame shear tests on woven textile composite reinforcements with controlled pretension[END_REF][START_REF] Launay | Experimental analysis of the influence of tensions on in plane shear behaviour of woven composite reinforcements[END_REF][START_REF] Nosrat-Nezami | Characterisation of the shear-tension coupling of carbon-fibre fabric under controlled membrane tensions for precise simulative predictions of industrial preforming processes[END_REF][START_REF] Kashani | Analysis of a two-way tension-shear coupling in woven fabrics under combined loading tests: Global to local transformation of non-orthogonal normalized forces and displacements[END_REF].

For instance, Willems et al. [START_REF] Willems | Picture frame shear tests on woven textile composite reinforcements with controlled pretension[END_REF] revealed that the shear resistance increased with the enlarged pre-tension. Launay et al. [START_REF] Launay | Experimental analysis of the influence of tensions on in plane shear behaviour of woven composite reinforcements[END_REF] confirmed that the shear force was affected by the membrane pre-tension induced by the different clamping of fabric on the picture frame devices. Nosrat-Nezami et al. [START_REF] Nosrat-Nezami | Characterisation of the shear-tension coupling of carbon-fibre fabric under controlled membrane tensions for precise simulative predictions of industrial preforming processes[END_REF] utilized air cylinders to apply constant pre-tension force, and indicated that the shear force was affected by the membrane pre-tension in the whole shearing process. Haghi Kashani et al. [START_REF] Kashani | Analysis of a two-way tension-shear coupling in woven fabrics under combined loading tests: Global to local transformation of non-orthogonal normalized forces and displacements[END_REF] reported that the shear force of the fabric inclined with the increasing pre-tension provided by the displacement control at the initial shear angle. While, the distinctions of the shear force under different pre-tensions were not obvious at high shear angle due to the relaxation of the pre-tension on the fabric.

It should be well noted that the above studies on shear pre-tension coupling behavior of the fabric mainly focus on the characterization of macro performance including such as shear force versus angle curves and shear deformations. Actually, the evolution of mesoscopic feature of fabric such as yarn width is a key factor in the formation of wrinkle defects in the shear deformation. While, up to now, the studies on the yarn width evolution are carried out in the shear test without applying pre-tension [START_REF] Zhu | An experimental study of in-plane large shear deformation of woven fabric composite[END_REF][START_REF] Hosseini | Identifying the distinct shear wrinkling behavior of woven composite preforms under bias extension and picture frame tests[END_REF]. That is the effects of applying pre-tension on the yarn width evolution in the shear deformation have not been revealed yet.

Thus, the objective of this paper is to explore the effects of the pre-tension on the shear behavior of a typical kind of fabric, namely, the comingled glass-polypropylene 2-2 twill woven fabric in view of both macro and mesoscopic scales. As the necessary experimental basis, a self-developed picture frame fixture, which included a relatively concise approach to apply pre-tension loading, was originally designed. Then, the shear experiments under different levels of pre-tension were conducted, and the macroscopic shear force and shear deformation of the fabric were discussed. Finally, the mesoscopic features in shear pre-tension coupling deformation, including undulation of fabric surface and variation of yarn width, were analyzed and discussed.

Then, the devised picture frame fixture and the exploration of the shear pre-tension coupling behavior provide an important experimental basis to the purpose of suppressing forming defect.

Experiments

Fabric and sample preparation

As shown in Fig. 2a, the comingled glass-polypropylene 2-2 twill woven fabric, which is widely used in vehicle engineering, is selected as a representative for shear pre-tension coupling tests in this work. The fabric is composed of 60% E-glass fiber and 40% polypropylene fiber. It can be seen that the woven tightness of the selected fabric is large, and the basic specification and geometrical properties are listed in Table 1. As illustrated in Fig. 2b, the cross shape sample was designed and prepared before tests. The shear region in the middle of the sample was 100 mm × 100 mm.

Considering the effects of the local bending of yarns closed to the clamp induced by the clamp design on the pure shear deformation of the shear region, the rovings, which were adjacent to the shear region, were removed to reduce boundary effect according to the reported literature [START_REF] Zhu | An experimental study of in-plane large shear deformation of woven fabric composite[END_REF]. Additionally, the distance between the shear region and the clamp was enhanced to be 55mm to reduce the boundary effect on the center shear region. In order to measure the full-field deformation by the digital image correlation (DIC) method, in the pre-treatment, black dots were uniformly marked on the surface of shear region. Currently, the shear pre-tension coupling behavior was investigated by the developed picture frame fixture, which is composed of the traditional picture frame fixture and the applying pre-tension device [START_REF] Launay | Experimental analysis of the influence of tensions on in plane shear behaviour of woven composite reinforcements[END_REF][START_REF] Nosrat-Nezami | Characterisation of the shear-tension coupling of carbon-fibre fabric under controlled membrane tensions for precise simulative predictions of industrial preforming processes[END_REF][START_REF] Kashani | Analysis of a two-way tension-shear coupling in woven fabrics under combined loading tests: Global to local transformation of non-orthogonal normalized forces and displacements[END_REF]. The current applying pre-tension devices mainly utilized the cylinder or motor to provide pre-tension, and measured the pre-tension value by the sensor, which usually require various equipments, and are costly. Thus, in this work, similarly, a picture frame fixture was originally developed by mechanical structural design to test the shear pre-tension coupling behavior of the fabrics. Especially, different from the reported picture frame fixtures [START_REF] Launay | Experimental analysis of the influence of tensions on in plane shear behaviour of woven composite reinforcements[END_REF][START_REF] Nosrat-Nezami | Characterisation of the shear-tension coupling of carbon-fibre fabric under controlled membrane tensions for precise simulative predictions of industrial preforming processes[END_REF][START_REF] Kashani | Analysis of a two-way tension-shear coupling in woven fabrics under combined loading tests: Global to local transformation of non-orthogonal normalized forces and displacements[END_REF], here, the designed picture frame fixture with the springs is relatively concise and convenient to apply pre-tension. As shown in Fig. 3, the fixture includes two parts: a tension exerting and a picture frame for shear exertion. The picture frame part consists of four high strength aluminum (A7075) bars, which were connected via four hinge joints.

The hinge joints were machined with high precision to reduce friction. The high strength aluminum was used to avoid the intrinsic deformation of the fixture and to minimize the weight for convenient operation.

The tension exerting part contains adjusting bolts, clamp plates, springs for exerting force, guide pillars and traction bars. The clamp plates and the traction bars were connected via four springs with a spring coefficient of 3.3 N/mm. The spring tension can be adjusted through the regulation bolt. A measure ruler is engraved on the side of tensioning bar to conveniently obtain a specific elongation of the springs (see Fig. 3). Two guide pillars were set to ensure smooth movement of the traction bars and clamp plates, in order to ensure uniform and parallel force. Fabric slip in test should be definitely avoided, thus, a rubber blanket was attached inside the clamp plates, and a small circular groove was machined on the inside surface of the top clamp plate. An alternative connector manufactured by monomer casting nylon was designed to fit the adapters of different testing machines. According to above design, this picture frame fixture can be used to explore the effects of different pre-tension on the shear behavior of the fabric under a large range of testing shear angles.

Fig. 3 The self-developed picture frame fixture to test the shear pre-tension coupling behavior of the fabric.

Experimental procedure

In shear pre-tension coupling tests, this devised picture frame fixture was installed in a universal testing machine (MTS E45.105). Shear tests were conducted on five levels of pre-tension including 0 N/mm, 0.6 N/mm, 1.2 N/mm, 1.8 N/mm and 2.4 N/mm. The loading speed is kept constant of 20 mm/min. For every pre-tension level, at least three repetitions were tested, and the averaged value was obtained for subsequent analysis.

The scheme of shear deformation mechanism in the fabric subjected to pre-tension is illustrated in Fig. 4. The parameters to characterize shear behavior under pre-tension is introduced as follow: N F is the tensile force measured by the testing machine.  is the theoretical shear angle, and can be calculated by the geometric relationship:

1 2arccos( ) 2 2 2 2 frame d L   = -= - + . (1)
Generally, the normalized NS F is calculated by [START_REF] Cao | Characterization of mechanical behavior of woven fabrics: Experimental methods and benchmark results[END_REF]:

( ) 2 2cos 2 N frame NS fabric FL F L   =  . ( 2
)
Similar to the reported literature [START_REF] Launay | Experimental analysis of the influence of tensions on in plane shear behaviour of woven composite reinforcements[END_REF][START_REF] Kashani | Analysis of a two-way tension-shear coupling in woven fabrics under combined loading tests: Global to local transformation of non-orthogonal normalized forces and displacements[END_REF], due to the inherent characteristics of the devised picture frame fixture, it can be seen that the pre-tension T F is applied by the springs, and the amplitude of the pre-tension remains constant since the variation of the springs' length is kept to be constant in tests. The direction of the pre-tension is also constantly along the spring direction, and is not along the yarn direction in tests.

Thus, this work investigates the effect of different levels of pre-tension on the shear behavior of the fabric.

The shear deformation process of the fabric was captured by a CCD and a camera(Point Grey GS3-PGE-91S6M; 9 million pixels; Resolution: 3376 x 2704) with a frequency of 1 frame/s. Then, the full-field displacement in shear process and the variation of fabric surface under different pre-tensions were calculated by digital image correlation software Vic 2D and Vic 3D, respectively. In order to obtain the yarn width at different shear angle, the corresponding pictures of shear deformation under specify shear angle captured by CCD were extracted and imported into ImageJ software. Then, the yarn width can be measured by the pixels, and the pixel value is converted to actual length using the scale bar. Finally, the yarn width was calculated by averaged value of 9 points randomly selected on three yarns to avoid the measurement error. Then, the variation of yarn width in shear process was also obtained. The comparison of practical shear angles measured by DIC results with the theoretical prediction by Eq. ( 1) is given in Fig. 6. It is noted that, generally, the measured shear angle is in agreement with theoretically calculated values under different pre-tensions in the considered range. While, when the theoretical shear angle is larger than40°, the measured shear angle gradually becomes remarkably lower than theoretical ones, indicating that the shear deformation reaches a limit above this shear angle, and the deformation becomes non-pure in-plane shear behavior. This phenomenon can be also confirmed in Fig. 7. As shown in Fig. 7,in order to compare the shear behavior before and after the shear limit, the local mesoscopic structure features before shear limit are extracted at shear angle of 30°and 40°, and the corresponding features after shear limit are extracted at shear angle of 50° to clearly reflect the variations. The points a, o, b are the intersections of weft yarn edge line and warp yarn edge line. The line oa and ob coincides with the edge of yarn at shear angle of 30° and 40° for all pre-tension levels, which means the fabric presents pure in-plane shear deformation, and demonstrates the sample design effectively limits the boundary effects. While, the line oa or ob is misaligned with the edges of yarn when the shear angle reaches large value of 50°, indicating that the in-plane shear deformation is not pure. Actually, the non-pure in-plane shear deformation at high shear angle of fabric is originated from the lateral compression limitation of the yarn. Besides, the excessive lateral compression of the yarn induces the out of plane deformation, which can easily cause wrinkle defects in the forming process. Therefore, the non-pure in-plane shear deformation of fabric should be avoided in fabric forming process.

Coincidence

Coincidence Misalignment

Effects of pre-tension on shear responses

According to Eq. ( 2), the normalized shear forces versus shear angles are shown in Fig. 8. It can be seen that the normalized shear forces rapidly rise when shear angle is higher than 40°, suggesting that this angle is a critical angle to determine the pure shear deformation limitation of this woven fabric, and wrinkle will be developed due Here, it is important to compare the shear forces which are affected by different pre-tension levels. According to the literature [START_REF] Chen | Formability optimisation of fabric preforms by controlling material draw-in through in-plane constraints[END_REF], a coupling factor CF is defined as: , ,0 () () ()

NS P NS F CF F    = , (3)
where ,

NS P

F

is the normalized shear force with pre-tension, and ,0 NS F is the normalized shear force without pre-tension. The coupling factor CF under different pre-tension levels is characterized in Fig. 9. Before 40° shear angle, as increase the pre-tension levels, the amplitude and slope of the CF curves rise, indicating the enhancement of pre-tension can improve the ability of the fabric to resist shear deformation. On the other hand, with the evolution of the shear angle, the CF curves under different pre-tension levels gradually become consistent, and approach to 1. That is the pre-tension cannot affect the shear force when shear angle is larger than 40°.

Besides, as shown in Fig. 9, at the beginning of shear deformation, the CF value is dramatically enhanced due to the increase of static friction among yarns when applying a pre-tension, which was also reported by Zhu's work [START_REF] Zhu | Mechanical behaviors of Kevlar 49 fabric subjected to uniaxial, biaxial tension and in-plane large shear deformation[END_REF]. Then, considering the enhancement of interaction force among yarns can reduce the evolution of fabric surface due to the fabric weaving structure. Thus, in order to explore the effect of the pre-tension on the interaction force among yarns, the evolution of un-sheared fabric surface under different pre-tensions was calculated through the software Vic 3D as shown in Fig. 10. The Z-coordinate of the fabric surface decreases with the increasing pre-tension. The maximum amplitudes are 1.14 mm, 0.73 mm, 0.69 mm, 0.65 mm for pre-tension level of 0.6 N/mm, 1.2 N/mm, 1.8 N/mm, 2.4 N/mm, respectively, indicating the interaction force is enhanced.

Obviously, the variation of surface is not significant when pre-tension is higher than 1.2 N/mm. While, the interaction force should be still increased remarkably as described of the static friction enhancement among yarns in Fig. 9, suggesting the deformation of yarn becomes difficult.

Variation of yarn width

As shown in Fig. 11, generally, in the view of mesoscopic scale, there are three typical deformation stages in shear deformation of woven fabrics [START_REF] Zhu | An experimental study of in-plane large shear deformation of woven fabric composite[END_REF]. The yarns firstly rotate with each other until the yarn gaps disappear. Then, the yarns rotate without the gaps and suffer both lateral and longitudinal compression. Local buckling of the yarns will also arise. Finally, the yarns reach a compression limit (shear locking).Beyond the shear locking angle, the yarns primarily experience lateral compression, and the in-plane shear deformation is not pure. The yarn width continues to decrease, leading to the out of plane waviness of fibers, namely, the appearance of wrinkle defect. Thus, obviously, the variation of the yarn width in the pre-tension levels shows clear enhancements at the same shear angle compared with that with no pre-tension. For instance, the ( ) w  is increased by 0.07 at the shear angle of 42° for 0.6 N/mm pre-tension, indicating the previous conclusions with regard to the yarn evolution under no pre-tension are not appropriate for that of under specific pre-tension. Then, in order to explore the effect of pre-tension on the yarn width variation, the variation of yarn width ( )

w  
is presented in Fig. 12b, and the value is calculated by: ( ) ( ) ( )

0 w w w   =  - , ( 4
)
where the ( )

0 w  is the yarn width of un-sheared fabric. It is remarkable that ( ) w  
decreases with the increasing pre-tension, suggesting that the pre-tension suppresses the lateral compression deformation of yarn. Especially, the suppression effect becomes remarkable with the enlarged shear angle. Actually, the reduction of ( ) w   with the enlarged pre-tensions will decrease the deformation in thickness direction of yarn. Therefore, this result confirms that the pre-tension can suppress the wrinkle defect which is always induced by large deformation in thickness direction of yarn. It can be seen that the ( ) w   is hard to be distinguished when pre-tension is higher than 1.2N/mm, indicating that the effect of pre-tension to suppress wrinkle defect cannot be unlimitedly improved. These results also reflect that the serious wrinkle defects induced by excessive shear deformation can be controlled by applying pre-tension, but cannot be unlimitedly suppressed by solely improving the level of pre-tension in the practical forming process.

In order to predict the yarn width evolution of the comingled glass-polypropylene 2-2 twill woven fabric under different pre-tension level, a model ( ) , T wF  considering pre-tension effect is proposed based on the reported work [START_REF] Zhu | An experimental study of in-plane large shear deformation of woven fabric composite[END_REF], and it is given as:

( ) ( ) ( ) , 0 cos T w F w A  =    , (5)
where the parameter A is the effect factor of the pre-tension. To simplify the equation, the concept of width variation ratio

( ) ( ) ,0 T r w F w
 = , which is the transformation of Eq. ( 5), is used here. As shown in Fig. 13a, the measured data of r versus shear angle is fitted by the model. The effect factor of parameters A under different pre-tension levels can be also obtained, and is listed in Table 2. The high value of correlation factor Adj (R2) demonstrates the fitting effectiveness of the model. That is the model can effectively capture the yarn width evolution. predict the variation of parameters A under the pre-tension range of 0.0-2.4N/mm. It can be seen that the parameter A is reduced with the increasing of pre-tension level.

While, the slope of parameter A curve decreases, indicating the effects of pre-tension on the yarn width variation gradually weaken. Then, substitute the exponential function into Eq. ( 5), the yarn width ( )

, T wF  when the fabric suffer shear deformation under pre-tension can be expressed as:

( ) ( ) ( ) ( ) ( 0
=   +    . (6)
Thus, in view of mesoscopic scale, Eq. ( 6) can capture the evolution of yarn width under different shear angle and especially under pre-tensions. That is the yarn width under shear pre-tension coupling can be predicted for the comingled glass-polypropylene 2-2 twill woven fabric. In addition, according to the Eq. ( 5), the parameter A has a non-linear dependence with the shear behavior when the shear angle is kept constant as shown in Fig 14a . While, as the effective interval of parameter A is limited by the applied pre-tension as shown in the Fig. 13b, the non-linear relationship is not significant in the effective interval of this work (Fig. 14b). With the decreasing of parameter A, the yarns width is enlarged, indicating the yarn width variation ( ) ( ) ( )

0 w w w   =  - is suppressed. While,

Discussion

Macroscopic shear performance

In this work, the shear force of the comingled glass-polypropylene 2-2 twill woven fabric shows an evident increase when subjected to pre-tension, as pre-tension enhances the friction between yarns. However, the improvement of shear force is not obvious when the shear angle is larger than 40°. This convergence phenomenon of curves in Fig. 8 is also reported in studies [START_REF] Matveev | A novel criterion for the prediction of meso-scale defects in textile preforming[END_REF][START_REF] Harrison | Press forming a 0/90 cross-ply advanced thermoplastic composite using the double-dome benchmark geometry[END_REF], and up to now, these works believe that this result is owing to the inherent design limitation of the test device, which causes a relaxation of the pre-tension force applied to the fabric force when the shear deformation is enlarged to specific degree. However, for the devised picture frame fixture in current work, the pre-tension provided by springs will not be relaxed in experiment. Actually, when the fabric reached a specific shear angle (about 40° in this work, i.e. stage 3), the gaps among yarns and the gaps among fibers disappeared, and the deformation approached to be no-pure in-plane shear as described in the Fig. 6 and 7. Besides, when the shear angle is larger than critical value, such as 40° in this work, the yarns rotation was largely limited, and the lateral compression was primary deformation mode of the fabrics. Thus, the shear deformation of the densified fabrics can be regarded as the deformation of an entity without gaps, making the shear force almost unaffected by pre-tension.

Mesoscopic deformation

As shown in Fig. 10 and 12, the pre-tension has significant effect on the mesoscopic geometric features of the fabric. However, when the pre-tension is above 1.2 N/mm, this effect is not remarkable. As shown in Fig. 15, the normal pressure between the surfaces of two yarns is enhanced when yarns are subjected to pre-tension, and will also result in a reduction in the cross-sectional height of the yarns. Thus, the undulation of the fabric surface will decrease. However, with the increasing of deformation degree of the yarn cross section, the degree of deformation difficulty is enlarged exponentially [START_REF] Wei | Preforming behaviors of carbon fiber fabrics with different contents of binder and under various process parameters[END_REF]. Thus, the effect of pre-tension on reduction of the undulation of fabric surface is not significant when the value of pre-tension is higher than 1.2N/mm (Fig. 10).As shown in Fig. 15, the fibers in the yarns will be arranged more closely as cross-sectional height of the yarn decreases, indicating that the ability of yarn to resist lateral compression is enhanced. Thus, the variation of the yarn width ( ) w   is reduced when subjected to pre-tension. While, the ( ) w   approaches to be constant when the pre-tension is larger than 1.2N/mm, since the cross-sectional height of yarn can be hardly changed more.

It is well know that the ( ) w   is a significant cause of wrinkle appearance, since the yarn width reduction will lead to the increase of the yarn thickness, and the nature of the wrinkle defect is the over increase of the yarn thickness [START_REF] Zhu | An experimental study of in-plane large shear deformation of woven fabric composite[END_REF]. That is the relative large

( )

w   will lead to the significant increase of the yarn thickness, and the wrinkle defect will arise more easily. In this work, when the pre-tension is applied on the yarns, we experimentally identify that at the same shear angle, the yarn width The effects of pre-tension on macro shear force and mesoscopic features, including undulation of fabric surface and variation of yarn width, were obtained and discussed.

Main conclusions are summarized as follows:

(1) A picture frame fixture was originally designed with a concise approach to apply pre-tension for testing shear pre-tension coupling behavior. The experimentally obtained shear deformation under applied pre-tension is uniform, demonstrating the effectiveness and reliability of the self-developed picture frame fixture.

(2) The critical shear angle of the comingled glass-polypropylene 2-2 twill woven fabric is experimentally confirmed to be 40°, and the in-plane shear deformation is not pure when shear angle is higher than 40°.The shear force is enhanced significantly by increasing pre-tension when the shear angle is lower than 40°.

While, the shear force can not be enhanced when shear angle is beyond 40°.

(3) The pre-tension suppresses the mesoscopic undulation of un-sheared fabric surface, suggesting the enhancement of coupling effect at the initiation of shear deformation due to the static friction. However, the suppression effect on the undulation is not obvious when the pre-tension is higher than 1.2N/mm. w   approaches to be constant when pre-tension is higher than 1.2N/mm.

In the future work, relative large shear angle of the picture frame fixture should be designed to induce obvious wrinkle defects, so that the correlations among the pre-tensions, yarn width evolution and wrinkle defects can be further revealed.

Fig. 1

 1 Fig. 1 Two protocols commonly used to characterize shear behavior of the fabrics: (a)

Fig. 2 (

 2 Fig. 2 (a): the comingled glass-polypropylene fabrics and (b): schematic diagram of

NSF

  is the normalized shear force in the fabric. T F indicates the pre-tension force applied to fabric. frame L is the distance between the adjacent rotatory joints. fabric L is the sample's clamping length. d represents movement distance of the picture frame device along N F .  is the angle between two bars.

Fig. 4 Fig. 5

 45 Fig. 4 Schematic diagram of the coupling between the shear and pre-tension

Fig. 5

 5 Fig. 5 X and Y displacement in shear deformation under a typical pre-tension of 2.4

Fig. 6

 6 Fig. 6 Comparison of the theoretically calculated shear angle with the measured shear

Fig. 7

 7 Fig. 7 Mesoscopic structures of the fabric under shear deformation with different

Fig. 8

 8 Fig. 8 Normalized shear force versus shear angle curves under different pre-tensions.

Fig. 9

 9 Fig.[START_REF] Wang | Experimental and numerical analyses of manufacturing process of a composite square box part: Comparison between textile reinforcement forming and surface 3D weaving[END_REF] The coupling factors with the evolution of shear angle under different levels of

Fig. 10

 10 Fig. 10 Contours plots of the measured distribution of the deformation of the fabric

Fig. 13b ,

 13b Fig. 13b, and the values are fitted by exponential function by using the least square

FrFig. 14 (

 14 Fig. 12 (a) Yarn width ( ) w  versus shear angle and (b) the variation of yarn width

Fig. 15 5 Conclusions

 155 Fig. 15 Diagram of the yarn deformation under pre-tension.

( 4 )

 4 The enlarged shear angle leads to the reduction of Yarn width. With the increasing pre-tension, the variation of yarn width ( ) w   is reduced, since the pre-tension enhances the ability to resist lateral deformation. While, the ( )

Table 1

 1 Specification and geometrical properties of the comingled glass-polypropylene fabrics.

	Weight	Area	Yarn linear	Thickness	Yarn	Yarn
	faction	density	density	(mm)	spacing	width
	(glass)	(g/m 2 )	(tex)		(mm)	(mm)
	60%	950	1870	1.4	4.9	4.17
	2

.2 Picture frame fixture

Table 2

 2 The values of correlation factor Adj (R2) and parameter A. obtained in this work versus pre-tension level are plotted as scatters in the

	FT(N/mm)

A

Acknowledgements

This work is supported by National Natural Science Foundation of China (Grant No. 51621004, 51875187 and 11972154), Hunan Province Science and Technology Plan (Grant No. 2019JJ40021), and Shenzhen Science Innovation Committee Foundation (Grant No. JCYJ20160530193728307). This work is also supported by the Fund of State Key Laboratory of Structural Analysis for Industrial Equipment under grant #GZ18106.

Data Availability Statement

All data obtained or analyzed in this work are included in this manuscript.

shear process, as a most important mesoscopic feature, is a significant index to reflect the different stages of shear deformation. Most importantly, the variation of the yarn with under different levels of pre-tension can clearly index the effect of shear pre-tension coupling effect. As the arrangement of this woven fabric yarns is tight, the variation of yarn width in the first stage of shear deformation, namely, small value of shear angle, is not pronounced. Consistent with the yarn width variation under condition of no pre-tension in the literature [START_REF] Zhu | An experimental study of in-plane large shear deformation of woven fabric composite[END_REF], with or without pre-tension, the yarn width shows a decreasing trend when the shear angle increases . While, the ( )


Costas Olga Sambataro

Costas Smaragdakis

Michael Taroudakis

A COMPARISON OF PROCESSING TECHNIQUES APPLIED TO TIME-FREQUENCY REPRESENTATION OF ACOUSTIC SIGNALS INTENDED FOR GEOACOUSTIC INVERSIONS

A comparison of processing techniques applied to time-frequency representation of acoustic signals intended for geoacoustic inversions

INTRODUCTION

Inverse problems concern the estimation of the parameters of a specific system using suitable mesurements of a measuraable quantityndirectly or indirectly associated with these parameters. In underwater acoustics, inversion techniques are used for estimating unknown parameters of an ocean environment. Typically these parameters conern the water column (sound speed profile) and/or the seabed (compressionaland shear velocities, densities, thickness of the sediment layers and attenuation). A propagation model is required to associate suitable physical or statistical "observables" required as input to the inverse prob-lem, which in turns are defined by some suitable transformation of the acoustic field which results in the appearance of the discrete characteristics of the acoustic signal required to define the observables. The starting point in the approach considered in this work is to use a single recording of the acoustic field due to a broad-band sound source, and perform a suitable time-frequency analysis of the recorded acoustic signal. The time-frequency representations show simultaneously the frequency and the time content of the signal and its energy distribution and may lead to suitable physical observables associated with the modal character of the acoustic field. The two most common time-frequency representations are the spectrogram and the scalogram [START_REF] Finn | Computational Ocean Acoustics. Modern Acoustics and Signal Processing[END_REF]. Spectrogram is obtained by taking the Short Time Fourier Transform (STFT) of the signal, while scalogram is obtained through the Continuous Walevet Transform (CWT). Due to its characteristics the second one can be proved more useful for analyzing realworld signals with complex time-frequency features.

In order to increase the resolution of the time-frequency distribution and to obtain a well separation between modes, several techniques have been proposed. Mode separation is required for many inversion procedures that exploit the modal content of the signal. In 2010 Bonnel [START_REF] Bonnel | Estimation of modal group velocities with a single receiver for geoacoustic inversion in shallow water[END_REF] proposed a signal processing method that involves time-frequency warping of the signal to improve the resolution of the time frequency distribution. This method has been applied in shallow water and in general provides a better separability of the modes than traditional spectrogram, facilitating the recognition of the modl dispersion in the environment under consideration. The second method reassigns the spectrogram and the scalogram and was first applied by Kodera, Gendrin, and de Villedary [START_REF] Kodera | A new method for the numerical analysis of time varying signals with small BT values[END_REF][START_REF] Kodera | Analysis of time-varying signals with small BT values[END_REF] and subsequently generalized by Auger and Flandrin [START_REF] Auger | Improving the readability of time-frequency and time-scale representations by the reassignment method[END_REF]. This method provides an improved identification of the modes due to the reallocation of the energy distribution in the time-frequency fomain [START_REF] Taroudakis | On the use of reassignment wavelet transform for mode identification[END_REF].

In the present work we perform both, the time-warping and the reassignment techniques in order to compare their efficiency in mode separation. We present here in brief the theoretical aspects of these two techniques and we apply them to both spectrogram and scalogram by means of a simulated and an actual acoustic signal. The simulated signal is generated by the MODE1 program developed at FO.R.T.H. provides the acoustic field in an oceanic waveguide given model parameters of the environment. In our simulated test case we have considered a Gaussian point source in a shallow water range-independent ocean environment. Finally, the considered actual signal is a denoised version of a raw measured signal due to a light-bulb implosion captured during the SW06 experiment [START_REF] Tang | Shallow Water' 06-A joint acoustic propagation/nonlinear internal wave physics experiment[END_REF].

DISPERSION NOTIONS

One of the main goals ofinversion techniques in acoustical oceanography is to exploit the information included in the ocean acoustic signals to determine the water column and/or the sea-bed geoacoustic parameters in which sound propagates. In such a scenario dispersion can play an important role as its character depends on the environment through which sound has propagated. According to the normal mode theory the complex acoustic field is represented at long ranges as a sum of terms each one corresponding to a single propagation mode.

p(r, z) = N m=1 Φ m (r)Ψ m (z), (1)
where N is the number of propagating modes (normal modes). The acoustic energy carried out by the signal propagates in wave packets and each of them have specific dispersion properties. Mode dispersion is characteristic of a specific acoustic environment and for this reason an appropriate processing of the acoustic measurements can provide information about the propagation medium, be proven as an efficient tool for inversion problems. Recall that each mode has its own travel time associated with its group velocity. Group velocity represents the energy transport velocity of a particular mode in the waveguide and it is defined by

v g,m = ∂ω ∂k m ω0 , (2)
where ω 0 is the central frequency of the signal bandwidth and k m is the wavenumber of the m-th mode.

The dispersion curves describe the relation of the angular frequency ω with the wavenumber of the mode m, k m and group velocity is the slope of the dispersion curves. This relation is non-linear and it is strictly connected to the fact that the sound channel is dispersive. The notion of group velocity has been associated with the arrival pattern of an underwater acoustic signal to identify the modal arrivals which have been used as the observables of an inverse problems exploiting the despersive characteristic of the sound channel for the estimation of the model parameters of the environment [START_REF] Taroudakis | Tomographic inversions in shallow water using modal travel time measurements[END_REF][START_REF] Bouchage | Fluctuations of the modal arrival times due to linear internal waves : Application to Inversion[END_REF]. Using modal arrivals the full dispersion characteristics of the environment are not exploited. For this reason spectrograms and scalograms representing a more detailed time-frequency analysis of the signal have been proposed for inversion purposed.

SPECTROGRAM AND SCALOGRAM

The spectrogram is generated by the Short Time Fourier Transform (STFT) of the signal s(t), with unit energy analysis window h(t) and given as follows:

STFT h (s; τ, ω) = ∞ -∞ s(t)h(t -τ )e -iωt dt (3)
and it is defined [START_REF] Finn | Computational Ocean Acoustics. Modern Acoustics and Signal Processing[END_REF] as the energy distribution associated with the STFT

SP h (s; τ, ω) = |STFT h (s; τ, ω)| 2 . ( 4
)
Scalograms are generated by Continuous Wavelet Transform (CWT). The latter provides a flexible time-frequency window which changes its amplitude automatically according to the frequency components and for this reason is particularly useful for analyzing signals at different scales.

The CWT of the signal s(t) with scale parameter a ∈ R and time-location parameter τ ∈ R is

CWT ψ (s; τ, a) = 1 |a| 1/2 ∞ -∞ s(t)ψ t -τ a dt, (5)
where ψ(t) is a continuous function in both the time domain and the frequency domain and it is the basic wavelet from which they are generated the translated and scaled windows. We remind that, in order that the definition of the CWT is valid, the basic wavelet ψ(t) must satisfy the admissibility condition

C ψ = ∞ -∞ |Ψ(ω)| 2 |ω| dω < ∞, (6)
where Ψ(ω) is the Fourier Transform of ψ(t). Finally, the scalogram is defined as

SC ψ (s; τ, α) = |CWT ψ (s; τ, α)| 2 . (7)

TIME-WARPING METHOD

Time warping [START_REF] Bonnel | Estimation of modal group velocities with a single receiver for geoacoustic inversion in shallow water[END_REF][START_REF] Bonnel | Geoacoustic inversion in a dispersive waveguide using warping operators[END_REF] consists a non-uniform resampling of the time coordinate of original signal. Warped time t and signal frequency in the warped domain f are defined in [START_REF] Godin | Normal mode dispersion and time warping in the coastal ocean[END_REF].

t = t 2 -t 2 r , f = f (t) 1 - t r t 2 , ( 8
)
where t r is a reference travel time and f (t) the instantaneous frequency of the signal arriving at the receiver at time t. For a range independent environment, like the one considered in this work, t r = r/c, where r is the horizontal distance between source and receiver and c is the sound speed in the water column. Here the signal is assumed to be emitted at t = 0.

The idea is to identify the modes in the warped spectrogram where they are well separated and then to perform the inverse time-warping (inverting eq. 8) to go back the original space. Examples of how spectrogram and scalogram appear after warping are shown in the figures associated with the applications that will be presented bellow.

REASSIGNMENT METHOD

The reassigned approach can be applied to both timefrequency representations(spectrogram and scalogram) ) [START_REF] Kodera | A new method for the numerical analysis of time varying signals with small BT values[END_REF][START_REF] Kodera | Analysis of time-varying signals with small BT values[END_REF]. It creates a modified version of these representations resulting in an increase of the signal energy components concentration so that the discrete energy pakets of the signal are better localized. The use of this method is justified by the fact that both scalogram and spectrogram are greatly influenced by the presence of cross-terms along the signal bandwidth and thist can lead to a disturbed visualization of the time frequency distribution. In particular, in this work we are interested in in the mode indentification so we search for a good separability between modes in the time frequency pattern.

We present here the basic elements of the reassigned method as they have been published in the relevant literture.

So far, we have defined spectrogram and scalogram as the squared modulus of a STFT and CWT, respectively, it has been shown in [START_REF] Auger | Improving the readability of time-frequency and time-scale representations by the reassignment method[END_REF] that these representations can be expressed as well as function of the Wigner-Ville Distribution (WVD) W (x; t, ω) of x(t)

W (x; t, ω) := ∞ -∞
x(t + t /2)x(t -t /2)e -iωt dt . [START_REF] Bouchage | Fluctuations of the modal arrival times due to linear internal waves : Application to Inversion[END_REF] In the next sub-sections we will give the expressions for the reassigned spectrogram and scalogram without going into details on their derivation as these can be found in the cited literature.

Reassigned spectrogram

Spectrogram as a function of WVD, and its reassigned version are defined by

SP h (s; t, ω) = ∞ -∞ ∞ 0 W (h; u, Ω)× × W (x; t -u, ω -Ω) dudΩ 2π . ( 10
) RSP h (s; t , ω ) = ∞ -∞ ∞ 0 SP h (s; t, ω)× × δ(t -t(s; t, ω))δ(ω -ω(s; t, ω)) dtdω 2π (11)
t(s; t, ω) = t - 1 SP h (s; t, ω) ∞ -∞ ∞ 0 u • W (s; u, Ω)× × W (x; t -u, ω -Ω) dudΩ 2π (12) ω(s; t, ω) = ω - 1 SP h (s; t, ω) ∞ -∞ ∞ 0 Ω • W (s; u, Ω)× × W (x; t -u, ω -Ω) dudΩ 2π , (13)
where t and ω are the coordinates of the center of gravity of the energy concentration.

Reassigned scalogram

Analogously, the following are valid for the scalogram

SC ψ (s; t, a) = ∞ -∞ ∞ -∞ W (ψ; u/a, ω 0 -aω)× × W (s; t -u, ω) dudω 2π (14) RSC ψ (s; t , a ) = ∞ -∞ ∞ -∞ (a ) 2 SC ψ (s; t, a)× × δ(t -t(s; t, a))δ(a -â(s; t, a)) dt da a 2 , ( 15
) t(s; t, a) = t - 1 SC ψ (s; t, a) ∞ -∞ ∞ -∞ u× × W (ψ; u/a, ω 0 -aω)W (s; t -u, ω) dudω 2π (16) ω(s; t, a) = ω - 1 SC ψ (s; t, a) ∞ -∞ ∞ -∞ ω× × W (ψ; u/a, ω 0 -aω)W (s; t -u, ω) dudω 2π , (17)
where t and ω are the coordinates of the center of gravity of the energy concentration as in the case of the spectrogram.

Here ω 0 is the central frequency of the basic wavelet's spetrum and â(t, α) = ω 0 /ω(s; t, a).

APPLICATION TO A SIMULATED SIGNAL

A simulated signal s(t) is calculated using the normalmode program MODE1, developed at FO.R.T.H. It refers to the environment shown in Figure 1, with source depth z s = 10 m, receiver depth z r = 30 m and distance between source and receiver (Range) 20 km. All the other parameters are resumed in table 1. The signal is modelled in the frequency domain using a Gaussian pulse of central frequency f 0 = 150 Hz and a bandwidth ∆f = 50 Hz. The simulated arrival of the signal in the time domain is shown in Figure 2 and the results using all the techniques discussed above, are collected in Figure 5. The warped signal is indicated with ws(t).

APPLICATION TO REAL SIGNAL

The techniques described in the previous sections are now applied to an actual signal, recorded in 2006 during the SW06 experiment carried out off the coast of New Jersey [START_REF] Tang | Shallow Water' 06-A joint acoustic propagation/nonlinear internal wave physics experiment[END_REF]. The acoustic signal was generated by the implosion of a glass light bulb, lowered into the water. The implosion occurred at a depth of approximately 22 m in a shallow water waveguide with almost range-independent characteristics. The water depth was 79.1 m and the receiver was placed 7 km away from the source. The environment of the SW06 experiment is shown in Figure 4. Table 2 presents both the known environmental parameters and the estimated parameters of the sea-bed obtained using [START_REF] Smaragdakis | Acoustic Signal Characterization using Hidden Markov Models With applications in Acoustical Oceanography[END_REF]. In order to obtain the theoretical dispersion curves for comparison with the dispersion characteristics obtained by the scalogram and spectrogram, these parameters have been used to create a simulated signal.

Recorded and simulated signal are compared in Figure 3. Note that the agreement of the recorded and simulated signals is not perfect. It is still adequate to justify the use of the sea-bed parameters as obtained by the inversion method presented in [START_REF] Smaragdakis | Acoustic Signal Characterization using Hidden Markov Models With applications in Acoustical Oceanography[END_REF]. Also, the agreement between the dispersion characteristics of the signal obtained by the scalogram and the spectrogram are in a good accordance with the dispersion curves obtained by the simulated signal. It is interesting to note that for this particular case the spectrogram using the warping operator gives a better representation with respect to scalogram obtained using the same operator. In contrary when the typical spectrogram or scalogram are computed the scalogram seems superior to the spectrogram. In all cases less uncertainty can be obtained after applying reassignment, as shown in Figure [START_REF] Taroudakis | On the use of reassignment wavelet transform for mode identification[END_REF].

CONCLUSIONS

In this work we presented a comparative study of methods concerning the identification of the modal structure of underwater acoustic signals in the time-frequency domain by means of the time-warping operator and the reassignment transform. Time-warping and reassignment are useful tools for improving the readability of time-frequency and time-scale representation of acoustic signals recorded in a shallow water environment. In order to compare the performance of both methods, we start from the comparison between the standard representations, spectrogram and scalogram, and then we apply both techniques separately and together to a simulated and a real signal. As we expected, since the scalogram is defined by means of the Wavelet transform, which represent the structures of the signal at different frequency scales, it generally performs better than the spectrogram. This is in particular true in the simulated signal. Although, in both cases scalogram does not reduce the uncertainty in mode identification, it offers a clearer separation between modes. An improved modal separation can be obtained by the time-warping method. As we can see in Figure 5 in the warped space, each mode is perfectly separated from the others. Furthermore, timewarping seems to work better for spectrograms than scalograms.

On the other hand, in order to improve the timefrequency resolution and to reduce the uncertainty, a reassignment technique is essential. In particular, it offers perfectly localized modes and gives very good results when it is applied to scalograms. As we can see in Figure 5 for the simulated signal and Figure 6 for the real signal, the agreement between scalograms and the analytical dispersion curves is almost perfect, while reassignment applied to spectrogram does not seem to offer significant advantages. These conclusions are valid for both standard and warped time-frequency representations. To summarize, time-warping techniques are useful tools when the modal structures of the signal are not well separated, while reassignment is an excellent tool to improve the timefrequency resolution of the signals. The use of the two techniques simultaneously can be proved useful when a better separation between modes by reducing the uncertainty in modal identification is required.

Figures 5 [

 5 Figures 5[(a),(c)] and 6[(a),(c)] show a comparison between spectrogram and scalogram for a simulated signal and the SW06 experiment recorded signal.

RangeFigure 1 :

 1 Figure 1: The environment for the simulated test case. The sound speed profile in water, sediment layer and substrate are c w , c sl , c s respectively.

Figure 2 :

 2 Figure 2: The signal for the simulated test case.

Figure 3 :Figure 4 :

 34 Figure 3: Comparison between recorded and simulated signals of the SW06 experiment.

Figure 6 :

 6 Figure 6: Time-frequency representations of the light-bulb implosion signal recorded during the SW06 experiment.

Table 1 :

 1 Environmental parameters for the simulated test case. The environment is shown in Figure2.

	Parameter	Value
	Range	20000 m
	cw(0)	1500 m/s
	cw(20)	1490 m/s
	cw(100)	1510 m/s
	Sediment thickness	20 m
	Sediment sound speed	1700 m/s
	Substrate sound speed	1800 m/s
	Sediment density	1500 kg/m 3
	Substrate density	1800 kg/m 3
	Sediment attenuation	0 nep/mkHz
	Substrate attenuation	0 nep/mkHz

Table 2 :

 2 Environment parameters used to create the simulated signal of the SW06 experiment.

10.48465/fa.2020.0351

e-Forum Acusticum, December 7-11, 2020

ACKNOWLEDGMENTS

We acknowledge support of this work by the project "Innovative Actions in Environmental Research and Development (PErAn)" (MIS 5002358) which is implemented under the "Action for the Strategic Development on the Research and Technological Sector" , funded by the Operational Programme "Competitiveness, Entrepreneurship and Innovation" (NSRF 2014-2020) and co-financed by Greece and the European Union (European Regional Development Fund).

s(t) spectrogram