C Dereure

R Hall

F Baber

D Loison

E Lescoute

J.-P Guin

J.-C Sangleboeuf

L Berthe

M Nivard

I Guven
email: iguven@vcu.edu

Laser induced dynamic fracture of fused silica: Experiments and simulations

Keywords:

Fused silica samples were subjected to laser induced shock loading. Laser flux was varied in order to obtain different amounts and characteristics of damage in the samples. Three dimensional damage and fracture maps of two identical samples impacted by high and low laser flux values were obtained using both optical microscopy and X-ray computed micro-tomography. Three prevalent fracture and damage patterns were identified. Peridynamic approach was used to simulate the laser impact conditions on the samples in order to explain the causes of the observed fracture and damage morphologies. A proprietary shock physics code, ESTHER, was used to calculate the transient kinetic energy imparted to the samples based on the experimental laser flux values. The kinetic energy values were then integrated over time and provided target values to match for the peridynamic impact conditions. The main fracture patterns were captured by peridynamic simulations with reasonable quantitative accuracy. Explanations for initiation and propagation of each of the fracture patterns were presented based on the peridynamic dynamic fracture simulations. Limitations of the computational approach and recommendations for future work is provided.

Introduction

Protection against hypervelocity impacts is a major issue in many fields of applications, among which are the spacecraft industry or the high-power laser facilities. Indeed, meteoroids and space debris can impact satellites at several thousands of kilometres per second, possibly damaging key equipment including optics, solar panels and electronics [START_REF] Cour-Palais | Career in applied physics: apollo through space station[END_REF]. Comparably, some instruments equipping the Laser Megajoule (France) or the National Ignition Facility (USA) can be hit by numerous shrapnel and high velocity debris originating from the target [START_REF] Tobin | Using silica aerogel to characterize hypervelocity shrapnel produced in high power laser experiments[END_REF].

Among the materials potentially subjected to hyper-velocity impact conditions, glasses occupy an important place, as they are widely used in the aforementioned industries. Fused silica (SiO 2 ) in particular is the basis of many of these glasses. However, studying glass behavior under direct projectile impact loading is a difficult task. Among the difficulties is the brittle nature of glass materials leading to fragmentation and pulverization of the sample. Laser-induced shocks can be used as a surrogate to generate the similar loading conditions [START_REF] Nebolsine | Laser simulation of hypervelocity impact[END_REF][START_REF] Pirri | Theory for laser simulation of hypervelocity impact[END_REF], but with the advantage of allowing to recover the sample for post-mortem analysis. In both cases, the loading generates a shockwave in the material, followed by an unloading wave (also called rarefaction wave).

Shockwaves can generate permanent densification, which is a wellknown phenomenon in quasi-static loading ( [START_REF] Bridgman | Effects of very high pressures on glass[END_REF][START_REF] Rouxel | Poisson's ratio and the densification of glass under high pressure[END_REF]) but is still an ongoing work under dynamic conditions ( [START_REF] Arndt | Shock-wave densification of silica glass[END_REF][START_REF] Okuno | A Raman spectroscopic study of shock-wave densification of vitreous silica[END_REF]). Shockwaves can also generate damage, which has different shapes and positions in the shocked specimen depending on the shooting conditions. Some samples were damaged on the front face, others underwent spalling damage on the rear face, and some were also damaged in their bulk. Spalling damage ( [START_REF] De Rességuier | Spallation of glass materials under laser induced shocks[END_REF][START_REF] De Resseguier | A 1D model for glass dynamic behaviour under explosive loading[END_REF]) is caused by the dynamic tension induced inside a sample due to the interaction of the release wave reflected from the free surface with the incident unloading wave. It produces one or several fragments ejected from the free surface. Front-face damage on the other hand is produced by stress gradient on the laser-irradiated zone. Failure waves ( [START_REF] Kanel | Failure waves in shock-compressed glasses[END_REF][START_REF] Kanel | Phenomenological description of the failure waves in glasses[END_REF]) can also be observed in some cases of shock-loaded glass. Existing studies on the shock-induced damage in fused silica are mainly about crater formation [START_REF] Michel | Hypervelocity impacts on thin brittle targets: experimental data and SPH simulations[END_REF], failure waves and spalling [START_REF] Davis | Laser-Shock-Induced Spall and the Intrinsic Strength of Glass[END_REF], but give no insight on the initiation and propagation of damage inside the material.

In order to gain a better understanding of the shock behaviour of fused silica, a series of laser shock experiments were performed using the ELFIE facility of the Laboratoire d'Utilisation des Lasers Intenses (LULI) at the French École Polytechnique. Fused silica samples were subjected to laser induced shock and subsequent damage. The shockwaves from laser impacts induced various types of modifications into the silica samples [START_REF] Renou | Investigating ramp wave propagation inside silica glass with laser experiments and molecular simulations[END_REF]. Two otherwise identical samples were selected for further computational analyses; one was subjected to higher laser flux and the other to lower. The samples were then characterized using optical microscopy and X-ray computed micro-tomography (called "micro-CT" henceforth) for quantification of fractures and damages. A number of distinct damage modes were identified. Peridynamic (PD) approach [START_REF] Silling | Reformulation of elasticity theory for discontinuities and long-range forces[END_REF]- [START_REF] Madenci | Peridynamic theory and its applications[END_REF] was used to quantitatively analyze the damage initiation and propagation in the samples. Peridynamics is a non-local continuum mechanics approach particularly well-suited for exploration of dynamic fracture problems with complex patterns of multiple cracks [START_REF] Agwai | Predicting crack propagation with peridynamics: a comparative study[END_REF][START_REF] Kilic | Peridynamic Theory for Progressive Failure Prediction in Homogeneous and Heterogeneous Materials[END_REF]. Explanations for multiple distinct failure patterns were provided using the peridynamic simulations.

The manuscript is organized as follows: Section 2 describes the samples, experimental setup, and characterization methods. A brief introduction to the peridynamic approach, the model used and the description of how the laser impact load is approximated are given in Section 3. Experimental observations and simulation results are presented in Section 4. Detailed analyses of different failure modes present in the experiments through peridynamic simulations along with limitations of the current approach are given in Section 5. Finally, Section 6 covers the conclusions and outlook.

Material and methods

Samples

The samples consisted of fused silica (spectrosil®2000) parallelepipeds of approximately 8 × 8 mm 2 and 1.5 to 3 mm thick. Samples were cut using a diamond saw to these dimensions to be compatible with the loading chamber dimensions. All surfaces were polished with diamond disks followed by diamond suspensions down to one micron size particle. A thin layer (10-15 μm thick) of pure aluminium was vapor deposited on the front face of each sample as illustrated in Fig. 1.

The fused silica samples were coated with a thin layer of aluminium for generation of a plasma cloud. During the laser-aluminium interaction, a thin layer of material is ablated into a plasma cloud which expands towards the laser source, generating a shockwave in the silica sample. Compared with other shock-generating methods as plate impact, the benefit of laser impulses is that the shock pressure is only applied to the laser-impacted zone followed by its rapid damping during propagation, which allows to reach high pressures without damaging the glass in a way that it is no more recoverable for postmortem investigation. Two specific samples were chosen for further computational analyses. The identical samples were shocked by a picosecond laser impulse in different flux domains: One was shocked in the high flux domain (6.01 TW/cm 2 , later called "sample HF", and the other in the low flux domain (0.61 TW/cm 2 , later called "sample LF"). Their thicknesses was 1.84 mm.

Experimental setup

A shockwave is formed by the expansion of the plasma generated by a laser-matter interaction. In order to have both as much control and repeatability as possible a thin vapor deposited layer of pure aluminium is used for laser-matter interaction. Laser shots were done at the ELFIE national facility (LULI laboratory, Ecole polytechnique, Palaiseau, France). The ELFIE laser source has the following characteristics: a 1057 nm wavelength energy up to 50 J with an impulse whose half width at mid height is 600 ps, focused on a 1.27 mm diameter spot on the front face of the target. Special care was given to avoid any local over-intensity in the laser spot. The pressure against time loading profile on the front face was assessed from one-dimensional lasermatter simulations with Esther code ( [START_REF] Colombier | Hydrodynamic simulations of metal ablation by femtosecond laser irradiation[END_REF][START_REF] Bardy | Numerical study of laser ablation on Aluminium for shock waves applications development of a suitable model by comparison with recent experiments[END_REF]). The laser shot experiment takes place under UHV vacuum conditions. The samples were recovered from laser shot experiments for post-mortem analyses.

Characterization methods

Optical microscopy was performed using a Keyence VHX 5000 microscope. Micro-CT measurements were conducted with a SkyScan 1173 X-ray to have a complete description of damage in the silica samples impacted by laser. The micro-CT investigation was done with a voxel resolution of 6.5 μm and an angle step of 0.2°. During the scan, the X-ray tube voltage and current were set to 80 kV and 80 μA, respectively. The measurements were executed without any filter, and for each scan a flat field correction was applied. Also, in order to enhance the contrast, the HF sample was submerged in a dye penetrant solution for 24 h in primary vacuum conditions prior to the scan. This penetrant was composed of 15 g of zinc iodide, 2.5 ml of distilled water, 2.5 ml of isopropyl and 2.5 ml of Kodak photo solution. The micro-CT scan images were then reconstructed using the NRecon software and the post-processing was performed with CTvox and Data Viewer.

Computational simulation method: peridynamics

Overview

Traditional mechanical simulation formulations are limited in their approach of damage [START_REF] Agwai | Predicting crack propagation with peridynamics: a comparative study[END_REF]. Damage in a solid generates discontinuities such as crack faces that are not straight-forward to address under the classical continuum mechanics stress-strain formalism. Current problem under consideration requires a simulation method that can handle a large number of concurrently occurring discontinuities that may grow in complex patterns under dynamic conditions, without prior knowledge of their origins. In these regards, peridynamic formulation is wellsuited for this study, as it is designed to predict damage under dynamic conditions. Therefore peridynamic (PD) approach [START_REF] Silling | Reformulation of elasticity theory for discontinuities and long-range forces[END_REF][17][START_REF] Madenci | Peridynamic theory and its applications[END_REF] was used in Fig. 1. Illustration of the experimental setup and of the generation of a shockwave by means of a laser impulse. The sample is placed in a vacuum chamber in which it is impacted by the laser. A thin layer of aluminium is ablated by the laser impact, forming a plasma cloud that expands towards the laser source. The counter-reaction on the sample generates a very high pressure which in turn generates a shockwave.

this study to quantitatively analyze the damage initiation and propagation in laser-shocked silica glass. Peridynamic theory is a non-local reformulation of continuum mechanics using integral forms, which removes the mathematical artifacts of singularities arising from spatial derivatives when geometric discontinuities such as cracks or material interfaces are present in the solution domain. The bond-based implementation of peridynamics is used in this study, as introduced by Silling [START_REF] Silling | Reformulation of elasticity theory for discontinuities and long-range forces[END_REF]. A material point interacts with other material points in its vicinity of a finite radius (δ, horizon); not only the nearest neighbors, but also those that are farther are considered in the family that affects the deformation of a material point (Fig. 2). The mutual influence of material points beyond each other's horizon is considered to diminish. The horizon parameter is an important part of the non-local formulation, allowing a length parameter to be part of the material behavior thus facilitating multi-scale analyses. When the material system under consideration has an underlying heterogeneity such as inclusions or microstructure, the horizon must be sufficiently small such that these features are represented with reasonable accuracy. Additionally, in the absence of material heterogeneity considerations, crack branching during dynamic fracture requires special care for horizon value selection. An extensive discussion on this topic is included in [START_REF] Bobaru | The meaning, selection, and use of the peridynamic horizon and its relation to crack branching in brittle materials[END_REF]. The equation of motion for a material point x at time t is written in the following form:

∫ ∂ ∂ = -′ - ′ ′- + ρ t t t t u T x x x T x x x b x { ( , ; ) ( , ; )} ( , ) H 2 2 (1)
where ρ and b denote the mass density and body force vector, respectively. In Eq. ( 1), H is the neighborhood of point x in which the members of its family, x ′ , are present. The displacement vector field u at point x results from the interaction between the material point x and points x ′ due to the pairwise force state vectors T(x, t; xx ′ ) and T(x ′ , t; x ′x). Details of the formulation and their relationships to engineering material properties are discussed in great detail in [17,[START_REF] Madenci | Peridynamic theory and its applications[END_REF]. It is worth noting that the term (xx ′ ) represents the relative position vector between points x and x ′ . Similarly, a relative displacement vector is calculated in the form of (u ′u) in which u ′ is the displacement vector of the point x ′ . Using these definitions, one can define a scalar quantity, termed stretch, between two interacting material points x and x ′ as

= -′ -′ - - -′ -′ s x x u u x x x x |( ) ( )| |( )| |( )| (2)
which is akin to normal strain in classical continuum theory. The constitutive relationships between the bond force and stretch between material points have been derived for various material models (see [START_REF] Madenci | Peridynamic theory and its applications[END_REF][START_REF] Foster | Viscoplasticity using peridynamics[END_REF][START_REF] Silling | Peridynamic theory of solid mechanics[END_REF] among others). As the solution domain is subjected to external stimuli leading to its deformation, forces and stretches between material points arise. In the current work, a linear-flat relationship between the bond force and bond stretch is assumed as illustrated in Fig. 3. In this material model, bond force and stretch are linearly related to each other up to a certain force value after which the bond cannot carry additional load while able to further deform. The deformation of the bond may continue to a value that leads to the termination of the bond; the stretch value s cr at this state is called the critical stretch. Silling and Askari [START_REF] Silling | A meshfree method based on the peridynamic model of solid mechanics[END_REF] formulated a relationship between the fracture energy G 0 and bulk modulus K of the material, and its peridynamic critical stretch as

= s G Kδ 5 9 . cr 0 (3)
It should be noted that the critical stretch given in Eq. ( 3) is valid in the general three-dimensional formulation. When the two-dimensional bond-based peridynamic formulation is used, the critical stretch takes the following form (Section 6.1 in [START_REF] Madenci | Peridynamic theory and its applications[END_REF]):

= s πG Kδ 3 . cr 0 (4)
Whether three-or two-dimensional formulation is used, the formulae given in Eqs. (3) and (4) are obtained using idealized peridynamic bond breakage configurations, and meant for providing an informed value. Methods for evaluating accurate values of critical stretch in peridynamics is an active research area.

If a bond reaches the flat portion but not the critical stretch, the unloading takes place from its current stretch value following a path parallel to the original loading curve. The input material parameter that decides the force level where the flat portion starts is referred to as, for lack of a better term, yield stress (σ y ) while recognizing that the current material model is not a rigorous implementation of plastic behavior.

The solution domain is discretized into smaller volumes represented by a material point in the center. The current implementation of the bond-based peridynamics uses cubic volumes of the same size in the entire solution domain. As the domain deforms, bonds exceeding the critical stretch value break permanently, which leads to redistribution of forces to the nearby bonds. In turn, those bonds may also break under the added load. During the peridynamic simulation, the number of bonds connected to each material point is monitored at every time step. A scalar parameter, damage, is defined at every material point such that its value is found by the ratio of number of broken bonds to the initial number of bonds at the material point. At the beginning of a simulation, damage value is zero while a value of unity indicates a material point that is disconnected from all of its original family members. This scalar damage parameter is not sufficient for deciding a crack and crack faces are present in the domain. There must be a group of material points having similar damage values in the form of a continuous curve. Once such a continuity is present in the results, the question of what the value of damage should be that describes a crack face arises. In this study, a value of damage greater than 0.4 is assumed to indicate potential sites for a crack faces. The rationale for using a value smaller than 0.5 is as follows: Consider a material point on a crack face. It represents a small cubic volume with its sides equal to the grid spacing. Prior to crack formation, in its pristine state, this point has all of its bonds intact, emanating from it to all the other material points in its family (those points within the sphere with a radius of horizon). Once the deformation and failure of bonds take place, the material point becomes a crack face by most of its bonds on one side (crack side) being broken while the bonds on the other side remaining intact. Because the material point is assumed to be on the crack face, its bonds that are parallel to the crack plane are also likely to be intact. Therefore more than half of the bonds of a material point on the crack face are likely to be intact, which puts the damage value somewhere between 0.4 and 0.5. If a larger damage value (~0.5) is used to qualify for a material point to be on a crack face, the measured crack lengths would be smaller as compared to those of a damage value less than 0.5. It is important to note that this definition of crack face is a post-processing action and that it does not affect the peridynamic solution.

The material points that are disconnected through bond breakage may find themselves on opposing crack faces. This is likely to happen when internal fractures such as those that are the subject of this study are present. The newly formed crack faces must not interpenetrate. In the current implementation, at every time step, the immediate vicinity of each material point is checked for other points that are not connected to the material point under consideration. When this is detected, contact conditions take effect. A contact algorithm (Section 10.1.2 of [START_REF] Madenci | Peridynamic theory and its applications[END_REF]) that calculates the dynamic contact force between material points that are not connected by peridynamic bonds has been implemented. When bonds between material points break, those points interact with each other through this contact condition. The contact detection and force calculations are independent of whether the crack is classified as mode I, II or III. However, the contact algorithm does not include friction in the current implementation.

Peridynamic model

Peridynamic simulations of both HF and LF samples were performed by using a 2-D rectangular peridynamic grid of spacing 8 μm with the horizontal dimension of 6 mm and a vertical dimension of 1.824 mm, as sketched in Fig. 4, resulting in nearly 171,000 material points. A horizon value of 24.12 μm (3.015 times the grid spacing) was selected. For glass materials, the maximum horizon size was suggested by [START_REF] Bobaru | The meaning, selection, and use of the peridynamic horizon and its relation to crack branching in brittle materials[END_REF] to be within 0.5 mm to 1.5 mm. The physical dimensions of the problem under consideration as well as the fracture-related length parameters characterized in Section 4 are considered to be within length scales where the classical continuum mechanics is valid, as discussed by [START_REF] Bobaru | The meaning, selection, and use of the peridynamic horizon and its relation to crack branching in brittle materials[END_REF]. This means the grid spacing value in conjunction with the horizon selection, has little effect on the solution. However, the grid spacing value will dictate the smallest crack lengths that can be resolved in the solution. A convergence case study had been performed in the current work to confirm that the main fracture morphologies predicted by a range of grid spacing values were close in size and orientation; the grid spacing size of 8 μm was found to be sufficiently refined while remaining computationally efficient. The material points located on the cylinder's circumferential vertical boundary are constrained from movement in x and y directions while they are free to move up or down.

The rectangular peridynamic domain is entirely comprised of silica glass to match the experimental conditions. In the simulations, mass density of ρ = 2.2 g/cm 3 and bulk modulus of K = 36.3 GPa was used for the silica glass [START_REF] Pabst | Elastic properties of silica polymorphs -a review[END_REF]. In this study, the stiffness of the peridynamic bonds are identical under tension and compression. However, the bonds are not allowed to fail under compression; bond breakage occurs only under tension. Although the peridynamic material model that was described in the previous subsection is linear flat, the material behavior for fused silica is specified as linear-only by specifying an arbitrarily large value for yield stress, σ y = 24 GPa, that has no physical meaning. The critical stretch is calculated using Eq. ( 4) with the fracture energy ranging from 3.7 to 5 J/m 2 as tabulated in [START_REF] Lucas | Determining fracture toughness of vitreous silica glass[END_REF]. The fracture energy value of 3.76 J/m 2 was found to yield satisfactory results leading to a critical stretch value of s cr = 2.12 × 10 -3 .

Laser shock loading calculations

As shock pressure and particle velocity are linked through the Hugoniot curve [START_REF] Meyers | Dynamic Behavior of Materials[END_REF], the applied transient velocity was obtained from the laser-matter interaction simulation provided by the Esther code [START_REF] Colombier | Hydrodynamic simulations of metal ablation by femtosecond laser irradiation[END_REF][START_REF] Bardy | Numerical study of laser ablation on Aluminium for shock waves applications development of a suitable model by comparison with recent experiments[END_REF]. In order to obtain the transient shock pressure and velocity time curves the interaction between the simulated shockwave and the Al-SiO 2 interface was considered [START_REF] Davison | Fundamentals of Shock Wave Propagation in Solids[END_REF].

The one-dimensional materials A and B shown in Fig. 5, assumed to be perfectly bonded, are in thermodynamic state 0 (pressure P 0 = 0 and particle velocity u 0 = 0) before laser-matter interaction. After the lasermatter interaction on material A, due to the shockwave propagation, part of material A is in thermodynamic state 1 (shock pressure P 1 and u 1 ), while the rest of the system remains in state 0 as represented in Fig. 5a. When the shock reaches the A-B interface two equilibrium conditions must be satisfied:

= = P P and u u A B A B (5)
Because A and B are different materials having different densities (ρ A and ρ B ) and sound velocities (C A and C B ), they exhibit two different shock impedances (Z A and Z B , Z i = ρ i C i with i = A, B). This may lead to two different scenarios:

• If Z A > Z B then a shock is transmitted in material B and a release wave is reflected in material A.

• If Z A < Z B then a shock is transmitted in material B and another shock is reflected in material A.

In the current study, material A is vapor deposited aluminium, and material B is silica. Impedances of aluminium and silica are 14.84 kg/s/ cm 3 and 12.98 kg/s/cm 3 , respectively, therefore when the initial shockwave interacts with the interface a shock is transmitted in silica and a release wave is reflected in the aluminium.

The pressure and velocity values of the thermodynamic state 2 can be determined by utilizing the Hugoniot curve of each material as shown in Fig. 5b. Note that in Fig. 5b that HF stands for High-Flux and LF stands for Low-Flux. Material A goes from state 1 to state 2 by the reverse Hugoniot curve (H A -). This state 2 also propagates in material B, initially at state 0, by the Hugoniot curve (H B +). State 2 is defined by the intersection of (H A -) and (H B +).

As the Hugoniot curve of silica is known from literature data [START_REF] Mader | Los Alamos Series on Dynamic Material Properties[END_REF][START_REF] Sugiura | Dynamic response of fused quartz in the permanent densification region[END_REF][START_REF] Chhabildas | Shock loading behaviour of fused quartz[END_REF][START_REF] Alexander | Changes to the shock response of fused quartz due to glass modification[END_REF], from the shock pressures in the aluminium layer computed by the Esther code, transmitted shock pressure in the silica glass samples can be calculated. In the case of HF sample, the shock pressure reached is higher than the values existing in the literature (see Fig. 5b). Thus Hugoniot curve of the silica was extrapolated using a second degree polynomial in order to estimate the pressure level.

The transient velocity values for both samples HF and LF at the laser impact site were then estimated by Esther code as plotted in Fig. 6; each curve represents the material velocity at the impact side of the fused silica in the 1-D approach. The bond-based peridynamic approach utilized in the current work is not suitable to solve a problem involving changes in density, pressure and temperature as governed by equations of state. Therefore, mechanical impact was used as a surrogate mechanism to introduce the laser impact onto the fused silica. This approximation may be considered suitable as laser impulses were already used as surrogate to hypervelocity impacts [START_REF] Pirri | Theory for laser simulation of hypervelocity impact[END_REF]. The impactor was made of a near-rigid material with a density of ρ = 2.649 g/cm 3 , a yield stress of σ y = 800 MPa, and a bulk modulus of K = 548 GPa. Further, the impactor was not allowed to fail in an effort to minimize the different mechanisms of energy dissipation associated with the impactor; this simplifies for identifying and keeping track of how the energy is partitioned during the impact.

The required impactor velocity in peridynamic simulations was found by utilizing the concept of action [35] in conjunction with the Esther simulations. The action S is defined by

∫ ∫ = = - S L E E dt ( )dt c p ( 6
)
where t is time, L is the Lagrangian, E c is the kinetic energy and E p is the potential energy. The goal is to find mechanical impact conditions that transfers a certain amount of energy, dictated by Esther simulations, over time to the substrate. Therefore, the definition of action given in Eq. ( 6) is modified to include only the action of the kinetic energy. It is assumed that the particle velocity of the impact face predicted by Esther is uniformly observed across the laser spot. Based on the particle velocity and the duration of the loading, action of the kinetic energy for a single layer of material points underneath the laser spot is found; this quantity is denoted by S Esther HF and S Esther LF for high and low flux laser impact cases, respectively. Peridynamic impact simulations were then performed to match these target action of the kinetic energy values by changing the initial impactor velocities iteratively. In the peridynamic simulations, the kinetic energies of the single layer of material points underneath the laser spot were numerically integrated over time and compared to the target Esther values in each iteration for convergence. The Esther action of kinetic energy values were found to be: S Esther HF = 3.36 × 10 -12 Js and S Esther LF = 2.29 × 10 -13 Js.

A number of impactor shape profiles were tested to obtain those that captured most of the experimental damage pattern; Fig. 7 shows the impactor shapes tested for LF and HF laser impact cases. Uniformly distributed impactors (a cylinder) led to fracture patterns that are not satisfactory. Impactor profile that is comprised of a solid cylinder combined with a ring having a wider radius, shown in Fig. 8a, was found to capture experimental fracture morphologies observed in HF laser impact case. It is hypothesized that densification at the center of impact disperses mechanical energy simultaneously leading to shearing at the edge of the densified region, which in turn facilitates availability of more energy around the edge to propagate and damage the sample. Incorporation of densification in the model is the subject of another study. Similarly, for the LF impact, the impactor shape profile that captures most closely the experimentally observed damage is presented in Fig. 8b. The impactor profiles used in the peridynamic analyses for LF and HF are different. In actuality, the difference between the spatiotemporal distribution of plasma-induced surface pressure might be a significant contributing factor to the differences between the damage morphologies observed in LF and HF samples. In an earlier work [START_REF] Fabbro | Physics and applications of lasershock processing[END_REF] plasma images were captured through high speed photography involving different levels of laser intensity. At low laser intensity, a uniformly luminous plasma throughout the impact event was observed while the high intensity laser impact led to uniformly luminous plasma in the early stages of the impact event followed by a clear ring formation later. A few peridynamic impactor profiles were experimented with using these results as suggestions for a starting point leading to the formations HF given in Fig. 8 to be suitable for the problem under

Results

Experimental results

The post mortem analysis of the silica samples revealed that the shocked samples displayed a wide variety of damage types on both the front and the rear faces. Some had both the front and back faces shattered while some had damage only on their front face; a few had no damage at all. From this set of different experiments two geometrically nearly-identical samples were selected. The difference between the two tests was the energy level of the laser impact. One of the samples was impacted using high flux domain (6.01 TW/cm 2 , HF sample henceforth) and the other using low flux domain (0.61 TW/cm 2 , LF sample henceforth). Table 1 lists the numerical values of energy flux and thickness for each sample. The HF sample was damaged on both front and rear faces while LF sample was only damaged on the front face as presented in Fig. 10. The top row of Fig. 10 shows optical microscope images of the front face (Fig. 10a), rear face (Fig. 10b), magnified rear face (Fig. 10c), and side view (Fig. 10d) of HF sample. Similarly, the bottom row presents the front (Fig. 10e) and rear (Fig. 10f) face images. The diameter of the damaged areas on the front faces of samples HF (Fig. 10a) and LF (Fig. 10e) are 1.78 mm and 1.4 mm, respectively. The rear face of sample HF displays heavy spalling damage, with a 2.5 mm diameter circular area of ejected matter (Fig. 10b andc) with a whitened multi-fractured surface. It also has obvious signs of internal damage underneath the surface. There is no visible damage neither on the rear face of LF sample nor internally as the front face damage can be readily observed from the back surface (Fig. 10f). Table 1 recapitulates the characteristics of both samples HF and LF.

The X-Ray micro-CT scan results for sample HF visualized in 3D from different angles of view in Fig. 11, starting from a plan view (top left) progressively changing to a side view (bottom right). The morphology presented in Fig. 11 reveals several internal damage modes in a quantifiable clarity that is not possible using optical microscopy. A specific cross-sectional side view is shown in Fig. 12. Underneath both front and rear face damage zones there are two networks of cracks parallel to the faces, labelled as H-cracks to represent horizontal cracks. There is also a system of vertical and tilted cracks that may or may not connect with each other; which are labelled as V-cracks and T-cracks, respectively. It is likely that what are labelled as the (i) horizontal, (ii) vertical, and (iii) tilted cracks are referred to in literature as (i) lateral, (ii) mode II shear, and. (iii) Hertz cone-like cracks, respectively. Nonetheless, they are labelled as H, V and T cracks throughout this manuscript for clarity. In addition to the fractures through the thickness, the CT-scan showed the presence of radial cracks initiating along the boundary of the column-shaped zone and propagating towards the free side surfaces. In order to quantify the fracture morphology, a total of eight cross-sections of the HF sample have been examined. Fig. 13 shows the top/plan view of the X-ray micro-CT of the sample on the left with the orientations of the cross-sectional side views shown on the right side of the figure. In all of these cross-sectional views the horizontal, vertical and tilted cracks are present. The vertical central axis in all eight planes are identical. Inspection of the cross-sectional side views at eight planes shows the variability and non-uniformity of the

Peridynamic simulation results

Peridynamic simulations of HF and LF samples were performed using the grid and boundary conditions described in Section 3. The damage contours at the end of simulations for HF sample are presented in Fig. 14 with the corresponding samples of different planes from micro-CT shown earlier in Fig. 13. The peridynamic simulation captured the V-and T-cracks as well as the rear-face H-cracks. Although a considerable amount of front-face damage is predicted, the morphology of the damage is not a perfect match for the experimentally observed front-face damage. One plausible explanation for the mismatch is that the peridynamic simulation does not take into account the equations of state governing the shock conditions dominant in the immediate vicinity of the laser impact.

The damage contours for LF sample obtained through peridynamic simulation are shown in Fig. 15. The simulation predicts damage to the front (impacted) face of the sample with no other damage morphology elsewhere, consistent with the experimental observation.

It is important to note that the simulation results are sensitive to the impactor shape profiles used. Although the action of the kinetic energy values for different impactor shape profiles match each other (and also the Esther values as described in Section 3.3), the damage predictions may be drastically different. This is demonstrated in Fig. 16a for HF laser impact case where damage contour map resulting from a slightly different impactor profile is shown; the action of the kinetic energy value is the same as that of Fig. 14. Similarly for LF laser impact case (Fig. 16b) the solid cylinder impactor profile generates damage map that includes back face manifestation contrary to Fig. 15 and the experiments. The simulations using the impactor shape profiles in Fig. 7 produced damage predictions with varying differences against the prediction shown in Figs. 14 and15, which were deemed to be the most closely correlated with the experimental observations. These results underline the need for further experimental characterization of spatiotemporal variation of the surface pressure generated by laser impact.

Predictions against experiments

The experimentally observed and peridynamic predicted fracture morphologies of HF sample presented in Figs. 13 and14 are similar in terms of general shape. As previously mentioned, V-cracks, T-cracks and rear-face H-cracks are present in both with some differences. The main difference is that the front-face H-cracks are not reproduced by the simulation. Instead the simulated front-face damage is concentrated under the impactors, while some of the material points on the front-face between the impactors are almost intact. This is most likely due to the actual material behavior in the first few micrometers of depth along the front face governed by shock physics equation of state that the simulations do not explicitly model.

In order to quantify the damage and fracture characteristics, the horizontal, vertical and tilted cracks are highlighted in Fig. 17 on the cross-sectional views of micro-CT images presented in Fig. 13. The specific damage characteristics measured are sketched on Fig. 18. Definitions of the parameters shown in Fig. 18 are given in Table 2. Based on these damage characteristics, peridynamic predictions are compared against the experimental measurements as listed in Table 3. The average experimentally measured horizontal crack radius H1 is 1920 μm with the corresponding peridynamic prediction of 1875 μm. The distance between the two opposing starting points of the vertical cracks near the front, H4, is predicted as 2290 μm as opposed to the average experimental value of 2874 μm. Another important measure is the distance between the two opposing tilted cracks, H6, which is experimentally measured as 1558 μm on average compared to the peridynamic prediction of a range of 1247 -2376 μm. The percent difference between the experiments and predictions for H1, H4 and H6 are calculated as 2.4%, 20.3% and 34.4%, respectively. When calculating the percent difference that involves a range of a parameter, the arithmetic average of the minimum and maximum values of that parameter is used. Additionally, the rear face damage depth related distances H2 and H5 as well as the horizontal distance travelled by the rear face horizontal crack after it turns towards the free surface, H3, exhibit greater differences between the experiments and simulations. These distances are small compared to the other distances measured (H1, H4 and H6) therefore the differences between the experiments and simulations lead to higher percent difference values. Finally, there are two characteristic angles of interest measured between the vertical axis and (i) the vertical cracks (V1), and (ii) the tilted cracks (T1) as indicated in Fig. 18. The percent difference between the experiments and the simulations for V1 and T1 are calculated to be 33.9% and 18.3%, respectively.

The fracture paths identified in Fig. 17 in the eight vertical planes are combined and plotted on top of each other, shown in Fig. 19a, to establish the variability of the damage morphologies observed. In addition, the peridynamic predicted fracture paths are plotted in dark on top of the experimentally observed fracture paths plotted in lighter grey in Fig. 19b.

In the case of the LF sample, in both the experiment and simulation, the only damaged zone is located on the front-face, directly under the impact. The diameter is measured in experiment as 1400 μm as compared to the prediction by the peridynamic simulation as 1120 μm, leading to 25% difference.

The peridynamic approach was able to capture the majority of the experimentally observed damage morphologies qualitatively. The spalled region, the V-cracks, T-cracks and rear-face H-cracks of the HF shot were present, as well as the surface damage zone on the LF sample. In addition, the order of magnitudes of the characteristic dimensions are close to each other. Considering the complexity of the multi-mode damage manifestation in the sample, irregular sample geometry, and uncertainties involved in quantitative understanding of the energy imparted on the sample by the laser plasma explosion, the peridynamic predictions are satisfactory. Simulations of other samples tested under conditions different than the current HF and LF would improve the knowledge and understanding of the laser-matter interaction and subsequent damage initiation and propagation.

Discussion

Simulations performed using the peridynamic approach provide data for convenient visualisation of the wave propagation and of the damage extent in the grid at multiple time steps. In the following subsections, the shockwaves are visualized by plotting the contours of horizontal (x 1 ) and vertical (x 3 ) velocity components in combination with the damage contours. These velocity and damage overlapping plots allow exploration of how the propagation of the waves is correlated with the initiation and propagation of H-, V-and T-cracks.

The generation of the front-face damage on the simulation of LF sample happens during the impact. The shockwaves that are formed at this moment reach the rear-face, but their reflection do not generate rear-face damage observed in the HF simulation. Therefore, only the HF sample results are analyzed in the following subsections.

Initiation of vertical cracks (V-cracks)

Series of plots given in Fig. 20 show damage and x 3 vertical velocity contours under the ring impactor immediately after the beginning of the impact event. The cracks initiate near the edges of the impactor due to shearing in the vertical direction, and propagate downwards. As the elastic waves from multiple contact points between the impactor and the target intersect, a complex system of fractures with a vertical symmetry plane appears directly underneath the impactor. Once the elastic waves initiated by the cylindrical impactor at the center of the top surface arrive and intersect with the existing crack system, a single prominent V-crack appears and propagates downwards. The initial propagation of the V-crack stops approximately 280 ns after the impact.

Initiation of horizontal cracks (H-cracks)

When the initial shockwave hits the rear-face, damage starts to appear underneath the rear surface. Shock theory [START_REF] Meyers | Dynamic Behavior of Materials[END_REF] states that the interaction between the reflected wave and the incident release wave generates tension, which can itself create spalling. The rear-face fracture that appears on the simulation immediately after the reflection of the shockwave is therefore most likely due to this phenomenon. But, as

Table 2

Descriptions of the parameters defined in Fig. 18.

Parameter Description

H1

Radius of the rear face horizontal cracks.

H2

Depth of the rear face horizontal crack measured from the rear face surface.

H3

Horizontal distance between rear face horizontal crack turning point and its intersection with the rear face.

H4

Distance between the two opposing starting points of the vertical cracks near the front (impacted) face.

H5

Depth of the rear face spalled region measured from the rear face surface.

H6

Distance between the two opposing tilted cracks.

V1

The angle between the vertical cracks and the vertical axis measured near the front (impact) face. T1

The angle between the tilted cracks and the horizontal axis measured near the rear face.

Table 3

Quantified fracture characteristics of HF sample observed in the experiment and corresponding simulation predictions. Definitions of the parameters are given in Fig. 18 along with their descriptions in Table 2. The experimental rear-face damage depth is highly heterogeneous, therefore the range is given. 2 This number was calculated by considering the maximum values of both experiment and simulation. 3 The H6 value in the simulation varies between the V-cracks (lower) and the T-cracks (higher).

seen in Fig. 21, after the initial clear spall fracture, a network of thinner cracks appear above it, which are also due to the tension caused by the crossing of the waves. These additional cracks are the H-cracks that slowly grow into those earlier seen in Fig. 14.

Initiation of tilted cracks (T-cracks)

Peridynamic simulations indicate that the T-cracks are extensions of the horizontal cracks (H-cracks). Their initiation points are located at the horizontal extremities of the H-cracks. Fig. 22a shows contours of horizontal (x 1 ) velocity and damage, plotted on top of each other, at 761 ns, marking the beginning of the T-crack initiation. Observation of the wave propagation prior to initiation helps identifying the mechanism leading to T-cracks curving up; Fig. 22b shows the snapshots x 1 velocity and damage contours at 533, 663, 694, 728, 761 and 793 ns. These cracks appear when two specific waves cross each other: the shockwave initiated on the far right end of the ring-impactor and reflected on the side, and the release wave initiated on the edge of the central cylindrical impactor (as shown in Fig. 22a)). The wave reflection pattern at the V-crack initiation time makes it difficult to see the intersection of individual waves clearly. In order to verify this hypothesis, the locations of the aforementioned shock and release wave fronts originating at different locations were quantified.

The velocities of these waves were measured in the simulation. The shockwave front travels at C s = 5.16 km/s, and the release wave front at C r = 3.07 km/s. In order for the two waves to meet at the V-crack initiation site, they have to cover the distances indicated by arrows in Fig. 22a. These distances were found to be L s = 3.581 μm ± 0.027 (shockwave initiated on the far end of the ring-impactor and reflected on the side) and L r = 2.080 ± 0.018 μm (the release wave initiated on the edge of the cylindrical impactor). Based on the measured distances and wave velocities, the predicted times for the waves to arrive at the initiation site were found to be Δt s = 0.694 ± 0.005 μs for the shock wave to cover the distance L s , and Δt r = 0.678 ± 0.006 μs for the release wave to cover L r . These numbers are compared to the times between the impact and the moment when the T-cracks start curving up in the peridynamic simulations, which are observed to be Δt s, PD = 0.705 ± 0.010 μs for the shockwave, and Δt r, PD = 0.686 ± 0.010 μs for the release wave. The difference between the time to Tcrack initiation based on the velocity/distance relationship for the shock and release waves and the actual predicted value are within their uncertainty ranges.

Propagation of vertical and tilted cracks

The mechanisms described in Subsections 5.1 and 5.3 are only associated with the initiation of the vertical and tilted cracks. Indeed, the simulation shows that their growth rate is not constant. Instead their growth happens in bursts, with periods when their length is constant and periods when they propagate.

As illustrated in Fig. 23a, the V-cracks end their initial extension from t = 280 ns after impact until t = 1424 ns. A growth step is between t = 1424 ns and 1567 ns, after which their growth stops again. The behaviour of the T-cracks follows a similar pattern. Their initial growth stops at t = 793 ns after impact, until t = 923 ns. At this moment they start growing at a steady rate until t = 1242 ns. This is illustrated in Fig. 23b. After the latest times shown in Fig. 23, the wave reflection patterns become highly complicated preventing clear association of individual waves and crack propagation steps. But this progression in bursts indicate that the growths of vertical and tilted cracks are due to the meeting of lingering reflected waves. Because of the several cracks and the two-dimensional effects, there are numerous waves propagating in the material. It is hypothesized that meeting of lingering waves in certain directional configurations in the vicinity of an existing crack tip may create favourable conditions for propagation of the V and T-cracks. Fig. 20. Time series of a V-crack under the ring-impactor during its initiation. In these plots, the damage parameter is superimposed to the x 3 velocity of the material points. The indicated time is counted after the impact. Although velocities and damage are displayed using the same LUT the former appears as a half tone background while the damage and notably cracks do appear as a bright blue-green localized contrast, which makes them distinguishable. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Conclusion

The mechanisms suggested by the peridynamics simulations give a possible scenario for the generation of the V-cracks, T-cracks and rearface H-cracks. De Rességuier [START_REF] De Rességuier | Etude expérimentale et numérique du comportement de matériaux de type verre soumis à des ondes de choc[END_REF] also observed T-cracks and rear-face H-cracks in soda-lime glass shocked by laser impulse. The H-crack initiation mechanism suggested by [START_REF] De Rességuier | Etude expérimentale et numérique du comportement de matériaux de type verre soumis à des ondes de choc[END_REF] relates to the spalling phenomenon in glass. Spalling happens when the release waves 1 and 1′ in Fig. 24 cross each other at the x r position and generate tension stress superior to the spall stress σ r for a specific time Δt r , as represented on the diagrams of Fig. 24. When fracture occurs, the compression waves 2 and 2′ are formed, moving at the speed C L , and "relax" the tension stresses. In metals these waves are faster than the release waves, and therefore spalling damage is limited to the surrounding of x r (Fig. 24). But in glasses, the velocity of waves 1 and 1′ presented in [START_REF] De Rességuier | Etude expérimentale et numérique du comportement de matériaux de type verre soumis à des ondes de choc[END_REF] is close to the sound speed C L , so relaxation is more progressive. The entire grey area of Fig. 24 is then subjected to an almost identical traction history. This causes multiple cracking in this zone, and the generation of new compression waves that gradually relax the stresses.

This explanation is in accordance with the peridynamics simulation. Indeed, the H-cracks appear when the initial shockwave interacts with the rear-face. The theory presented in [START_REF] De Rességuier | Etude expérimentale et numérique du comportement de matériaux de type verre soumis à des ondes de choc[END_REF] and the simulation show the same phenomena at the same moment.

De Rességuier [START_REF] De Rességuier | Etude expérimentale et numérique du comportement de matériaux de type verre soumis à des ondes de choc[END_REF] also proposed an explanation for T-cracks, based on the "crack-branching" phenomenon. It is characteristic of brittle materials, and happens when the propagation speed reaches the maximum limit. If the traction stress continues to increase, then the crack separates in two secondary cracks, as illustrated in Fig. 25. They propagate at a lower speed and with a different angle as compared to the initial direction of the crack.

There is indeed a secondary crack in the CT-scan and on the final peridynamic damage that has approximately the same origin point than the T-cracks. But in the crack-branching phenomena the two secondary cracks appear simultaneously. Yet, in Fig. 23b, it can be clearly seen Fig. 21. Initiation of the rear-face H-cracks in the simulation. The damage parameter is superimposed to the x 3 velocity of the material points to see the shockwave. The shockwave front is clearly seen on the top plot, and the damage starts appearing on the middle plot (377 ns); the indicated time is counted after the impact. that the T-crack appears before the secondary crack. It should be noted that the explanations proposed by De Rességuier in [START_REF] De Rességuier | Etude expérimentale et numérique du comportement de matériaux de type verre soumis à des ondes de choc[END_REF] regarding crack branching refer to idealized situations where an already existing crack perpendicular to the direction of a monotonic and persistent loading branches. By the time T-cracks evolve from the H-cracks in the peridynamic simulations, the wave landscape is highly complex and transient, which is likely the case in-situ. While De Rességuier [START_REF] De Rességuier | Etude expérimentale et numérique du comportement de matériaux de type verre soumis à des ondes de choc[END_REF] provides valuable insight under highly idealized configurations, its direct quantitative application to the problem under consideration may be questionable. Further, in [START_REF] De Rességuier | Etude expérimentale et numérique du comportement de matériaux de type verre soumis à des ondes de choc[END_REF], De Rességuier did not observe extensive front-face damage in laser-impacted soda-lime glass. In his work, the front-face damage is limited to the irradiated surface and a very thin layer under it. This is different from the current observations in silica glass, where the front-face manifests damage even under low flux impulses.

Presence of failure waves was not encountered in experiments or in peridynamics simulations in the current work. According to Brar et al [START_REF] Bless | Failure Waves in Glass[END_REF] failure waves are specific waves that propagate in a stressed brittle material, where the material is intact ahead of the wave and comminuted behind it. They are slower than the shockwave (as illustrated in Fig. 26). Kanel et al ([12,[START_REF] Kanel | Transformation of shock compression pulses in glass due to the failure wave phenomena[END_REF]) studied this phenomenon and compared it to a detonation wave. But their exact origin and nature are unknown, and are subject to debate. As none of the shocked samples was recovered in a pulverised state it is certain that this phenomenon was not observed in the present study.

Overall, the 2-D peridynamics simulation is in good agreement with the experimental results found in this study and in [START_REF] De Rességuier | Etude expérimentale et numérique du comportement de matériaux de type verre soumis à des ondes de choc[END_REF]. It revealed that the origin of some of the internal damage (the H-cracks) are caused by the initial shockwave propagation and to the brittleness of the material. On the other hand, the T-cracks are caused by wave reflections on the sides, and are therefore linked to the geometry of the target. Therefore, according to the peridynamic simulations, targets larger than the ones used in this study may have T-cracks with different geometric characteristics or may not be present at all. This is the subject of future testing.

The peridynamic implementation used in this study has a number of limitations. First, the temperature effects present in shock events are not represented in the current peridynamic analyses. In addition, rate sensitivity of the material properties including critical stretch is not taken into account. Finally, in the current study, a material is defined by its density, bulk modulus, yield stress and critical stretch. Although some work has been performed in the area of implementation of equations of state in peridynamics [START_REF] Silling | Modeling Shockwaves and Impact Phenomena with Eulerian Peridynamics[END_REF], current work does not include such formulation. As a result, densification of the silica under shock loading is not represented in the model.

The representation of laser impact loading conditions in a mechanical simulation remains an open question. In the current work, laser conditions (power, spot size, duration) are used in a 1-D shock physics code to predict the impact surface particle velocity variation. Peridynamic approach was used to provide possible explanations on the specific mechanisms at play in initiation and propagation of different fracture modes. While maintaining relationship with the 1-D shock physics predictions, 2-D peridynamic simulations were able to capture most fracture morphologies observed in the experiments qualitatively, and to a certain degree quantitatively. In order to increase the fidelity of peridynamic approach in its application to laser impact damage problems, a number of additional items should be considered: (1) Insitu measurement of surface velocities, both the impact and free surfaces if possible, (2) exploration of surface pressure distribution, both in space and in time, during the laser impact, (3) implementation of equation of state behaviour for silica in peridynamic framework in order to be able to take into account the densification behavior, (4) incorporation of rate sensitivity of critical stretch and thermal effects in the peridynamic simulations, and ( 5) performing peridynamic simulations of silica samples impacted under different laser conditions.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Two-dimensional illustration of bond-based peridynamics. The horizon value in this case is δ = 3 × grid spacing. The central material point (in red) is linked to all the material points (in green) contained in a δ radius sphere centred on it. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 4 .

 4 Fig. 4. Two-dimensional peridynamic model domain. The small rectangular region (orange-shaded) on the top side represents the laser-impacted region.

Fig. 5 .

 5 Fig. 5. a) Principle of transmission and reflection of a shockwave when it interacts with an interface between two different materials. b) Transmitted shock and reflected release wave in samples High-Flux (HF) and Low-Flux (LF). State 2 in sample HF was determined through extrapolation of the silica glass Hugoniot curve based on literature data.

Fig. 6 .

 6 Fig.6. Transient velocity profiles that were estimated using ESTER code on the front-face for samples HF and LF.

Fig. 7 .

 7 Fig. 7. Impactor shape profiles explored in search of close correlation against the experimentally observed damage morphologies; left and right columns show LF and HF cases, respectively.

Fig. 8 .

 8 Fig. 8. Impactor shape profiles used in peridynamic simulations for (a) high flux (HF), and (b) low flux laser impact.

Fig. 9 .

 9 Fig. 9. Physical interpretation of the peridynamic impactor loading used in the PD simulations for HF laser impact case: (a) Peridynamic pressure distribution on the sample surface due to the impactors at an instant, (b) cross-sectional view of the sample and the pressure field, (c) equivalent peridynamic impactors to simulate the nonuniform transient pressure, and (d) schematic representation of time variation of surface pressure.

Fig. 10 .

 10 Fig. 10. Optical microscope pictures of the front (left) and rear (middle and right) faces of samples HF (top) and LF (bottom) laser impacts taken with a x100 zoom on pictures a), b), d), e), f) and a x50 zoom on picture c). Picture f) shows only the front-face damage by transparency through the glass, the rear-face being intact. The front and rear-face damages are highlighted by circles. On c) the radial cracks are outlined by dotted lines. The grey areas near the damaged zones on the front-faces on a) and e) are the remnants of the aluminium coating. The damage visible on e) is the internal damage seen by transparency through the side of the whole sample.

Fig. 11 .

 11 Fig. 11. Rotation of a 3D view of the damage in the bulk of sample HF measured by micro-CT (images done by CTvox). Both front and rear face of the actual sample are not visible. Top left corner: view of the damage from the front surface rotating progressively to a side view in the lower right corner. (Scale bar shown for only the bottom right due to other images having perspective projectionthe left column is zoomed in while the right column has the same scale).

Fig. 12 .

 12 Fig. 12. Damage types and dimensions observed on sample HF micro-CT scan. A cross-sectional view in the x 2x 3 plane, with both front and rear-faces outlined. The front and rearface surface damage are visible in the centre of the picture and the V, T and H-cracks are labelled.

Fig. 13 .

 13 Fig. 13. Micro-CT cross-sections of HF sample. On the left, a top view (x 1x 2 plane) of the sample is shown with eight different planes passing through x 3 -axis. On the right, the side views of each of the eight vertical planes with the horizontal, vertical and tilted cracks are presented.

Fig. 14 .

 14 Fig. 14. Damage contours from the peridynamic simulations for the HF sample.

Fig. 15 .

 15 Fig. 15. In HF sample, (a) damage contours from peridynamic simulations, and (b) the fracture morphology observed using micro-CT.

Fig. 16 .

 16 Fig. 16. Damage contours from peridynamic simulations for (a) HF and (b) LF laser impact cases with different impactor shape profiles.

Fig. 17 .

 17 Fig. 17. Horizontal, vertical and tilted cracks observed in micro-CT scan of HF sample highlighted with thicker lines shown in selected eight vertical planes passing through x 3 -axis. The orientation of the vertical planes are shown in the top view of the sample on the left.

Fig. 18 .

 18 Fig. 18. Specific parameters defined for quantification of damage characteristics.

Fig. 19 .

 19 Fig. 19. (a) Experimentally observed fracture paths in HF sample on eight vertical planes examined plotted on top of each other, and (b) the peridynamic fracture predictions (dark) plotted on top of the experimentally observed fractures (light grey).

Fig. 22 .

 22 Fig. 22. Mechanisms of the generation of the T-cracks. The damage parameter is superimposed to the x 1 velocity of the material points to enhance the traction on the horizontal axis. a) View of the whole target to show the general phenomenon. b) Detailed view of the moment when the T-crack starts to curve up. The indicated time is counted after the impact.

Fig. 23 .

 23 Fig. 23. Illustration of the non-constant growth rate of the V-cracks (a) and the T-cracks (b). The damage parameter is superimposed to the x 1 velocity of the material points to enhance the traction on the horizontal axis. All times are counted after the impact.

Fig. 24 .

 24 Fig. 24. Schematic representation of the spalling phenomenon (a) in metallic materials and (b) in glass-type materials according to [37].

Fig. 25 .

 25 Fig.[START_REF] Silling | Peridynamic theory of solid mechanics[END_REF]. Illustration of the crack-branching phenomenon, which occurs in brittle materials submitted to very high tension stresses. The crack starts on the left edge, and propagates towards the right. It then separates when its propagation velocity reaches a certain limit.

Fig. 26 .

 26 Fig. 26. Distance-time diagram (taken from[START_REF] Bless | Failure Waves in Glass[END_REF]) illustrating the way a failure wave propagates and impacts the other waves in the material.

Table 1

 1 Sample and test characteristics of the two experimental shots that were chosen for computational analysis.

	Sample	Thickness	Energy flux	Front Face damage diameter	Rear Face damage diameter
	HF	1.839 mm ± 0.01	6.012 TW/cm 2	1.78 mm	2.50 mm
	LF	1.835 mm ± 0.01	0.607 TW/cm 2	1.40 mm	None

Acknowledgements

The authors gratefully acknowledge on the one hand the French Agence Nationale de la Recherche for its financial support through the program GLASS ANR-14-CE07-0020 and on the other hand the French ministry for research for a ½ PhD grant support. We would like to acknowledge Luli laboratory and specifically the team support of the ELFIE facility for their professionalism and their kindness. We are also in debt for H. Orain for his help regarding the preparation of the glass samples.


Nicolas Pirc

Florian Bugarin

Fabrice M Schmidt

Marcel Mongeau

Cooling channel optimization for injection molding

Keywords: BEM, heat transfer, injection molding, mold-cooling simulation, optimization, SQP PACS: 02.70.Pt

published or not. The documents may come

INTRODUCTION

Numerical simulations are more and more used for designing injection molds. The location of the cooling channels is a major element in the design of the mold because the cooling time can represents up to 70% of the injection cycle. Rapid prototyping processes such as layered design or selective laser sintering start to be used to manufacture injection molds. The advantage of these assembling processes is the possibility to obtain almost any desired shape of mold geometry. However; it is then more difficult to locate and to shape the cooling channels. Numerical simulation can be used to perform automatic optimization of the position and the shape of the cooling channels.

Simulation of Heat Transfer during Injection Mold

Numerical simulation is used to solve the heat balance equation and evaluate a cost (or objective) function related to part quality and/or productivity. Several numerical methods such as Finite Element Method (FEM) [1] or Boundary Element Method (BEM) [START_REF] Polynkin | international polymer processing[END_REF] can be used. Using BEM makes it possible to reduce by one the dimension of the mesh involved, meshing only the contour surface.

Mold-Cooling Optimization Approaches

The modeling of the heat transfer using BEM makes it possible to re-mesh only the surface of the channels during the optimization. An optimization method is used to modify the location and shape parameters of the channels and improve the cooling performance of the mold. L. Silva [START_REF] Silva | international polymer processing[END_REF] uses 3D FEM software to model injection molding cycle for complex geometries of molds. However, the computational burden of the 3D renders the integration of optimization difficult.

Mathey [START_REF] Mathey | AIP Conference Proceedings[END_REF] used 2D stationary and transient BEM simulations of injection molds coupled with an SQP algorithm (Sequential Quadratic Programming) [START_REF]Optimization Toolbox For Use With Matlab[END_REF] to improve mold injection cooling.

In this paper, we extend Mathey's approach in view of addressing 3D optimization. An other contribution is the improvement of the optimization formulation. We perform numerical simulations on a semi-industrial injection mold designed for the European project: Eurotooling 21. We optimize cooling channel locations on both 2D and 3D models. We attempt at minimizing the maximal cooling temperature of the plastic part subject to a temperature uniformity constraint.

STATIONARY HEAT TRANSFER PROBLEM

The variation of the temperature in the injection mold in production, after some cycle, is a quasistationary mode. After the average temperature of the mold during a cycle is stabilized (Fig. 1), we consider a stationary regime. We neglect the transitory oscillations of the temperature.

Solving Heat Balance Equation using Boundary Elements Method

The stationary heat conduction problem assuming constant thermal conductivity reduces to Laplace equation.

€ ΔT

= 0 (1)
Multiplying this equation by a weighting function T*, and using Green's theorem, we obtain the wellknown Somigliana's equation [START_REF] Brebbia | Boundary elements an introduction of course[END_REF]:

€ C.T P ( ) + T. ∇T * .n ( ) .dΓ Γ ∫ = ∇T.n ( ) .T * .dΓ Γ ∫ (2)
T is the temperature, and C is equal to 1 inside the domain and to 0.5 on its regular boundary. T* and q* denotes the fundamental solutions of the problem socalled Green's functions [START_REF] Brebbia | Boundary elements an introduction of course[END_REF]:

€ 2D T * = -1 2π lnr q * = -r.n 2πr        3D T * = 1 2πr q * = -r.n 2πr 2        (3)
After reorganization of equation ( 2), we obtain a non-symmetrical linear system.

Boundary Conditions

€ -λ ∇T .n = h T -T C ( ) (4)
On the cavity surface, the flux density is calculated from the cycle time and polymer properties (heat due to polymer crystallization is neglected) [START_REF] Jui- | Multi-objective optimization scheme for quality control in injection molding[END_REF].

Thermo Physical Parameters

All the parameters necessary for the numerical heat transfer simulations are referenced in Table 1.

Mold Cooling Optimization

Figure [START_REF] Silva | international polymer processing[END_REF] displays the coupling between the thermal solver and the optimization algorithm. Heat transfer computation is coupled with optimization method in order to modify automatically parameters. BEM simulation is performed and cost function is computed. The optimization method allows updating parameters with respect to constraints until a minimum of the cost function is reached. Matlab® SQP, based on the Newton method, is used for the optimization of continuous non-linear functions with continuous nonlinear constraints.

Cost Function and Constraints

SQP algorithm minimizes just one objective. One method can consider the cost function as the weighting of two objectives, but this method imposes to choose a weighting coefficient. We propose to use the first objective as an optimization criterion (5), and the second as a non-linear constraints (eq 6):

€ Φ = max(T) (5)
€

and 0 ≤ (T i -T moy ) 2 i ∑ ≤ 4 (6)
The first criterion is the maximum of the temperature on the cavity surface [START_REF]Optimization Toolbox For Use With Matlab[END_REF].

The second criterion (6) improves the uniformity. For that, Park [START_REF] Wen-Hsien | integrated numerical simulation of injection molding using true 3D aproch[END_REF] proposes to minimize the variation of temperature distribution on the cavity surface compared to its average temperature. This function is defined as non-linear constraints. The final configuration provides an uniform and high average temperature, which give a very long cooling time.

Optimization Parameters

In this paper, we consider circular cooling channels in 2D and cylinder in 3D, 3 parameters are needed to describe each cooling channel. Our methodology can also account for optimizing non-linear geometries parameters. These processing parameters can define for instance the thermal regulation: temperature and flow rate of the coolant.

APPLICATION 1: RESULTS FOR A 2D CASE

In order to validate our approach, we consider the 2D injection geometry as described in Table 1 and2.

Geometrical Parameters

Our methodology allows for various industrials constraints such as the definition of forbidden zones where one cannot put the cooling channels (due for the presence of the ejectors, or to keep the cooling channel within the mold, or to avoid inter-channel collisions). δ is the possible displacement according to axis X and Y for the channels.

Optimization Parameters

Recall that the stationary BEM is used to evaluate the objective function at the mold cavity surface.

The optimization parameters for the locating of the i th cooling channel are (xi,yi) (see Figure 4). Thus, this 2D problem involves 16 optimization variables. As illustrated Figure 5, after only 3 optimizations iterations, the objective-function value is already within 5% from the final (optimized) value. The cost function is divided by the initial value to insure 0 ≤φ≤ 1. We observe that both temperature variance and temperature average decreased significantly. Figure 7 shows the temperature distribution inside the mold temperature. We now report computational results on a 3D plastic part whose features are displayed on Figure 8 (unit in mm). The forbidden zones and the input data are identical to those of the 2D case.

Meshing Parameters

The data for the mesh are reported in Table 7. The geometrical optimization parameters are now the coordinates of points P1 and P2 (see Figure 9). The shape of the channels is assumed to be cylindrical and to have a rigid-body motion. Indeed, P2 can then be expressed in terms of P1 coordinates and of the constant channel length L (Figure 9): Thus, the optimization parameters for positions of the i th cooling channel are xi, yi, zi. Our 3D problem therefore involves 24 optimization variables.

Sensibility of the Cost Function

As is illustrated Figure 10, the speed convergence of the cost function is fast. As it is illustrated in Figure 10, the form of the convergence of the cost function is fast. The value of the cost function is reduced by 90% in 4 iterations. SQP cannot respecting the constraints during optimization if the converge more quickly. This is an advantage in the industrials problems because an initial location of the channels, respecting all the P 1 (x i ,y i ,z i )

L P 2 (x i ,y i ,z i ) X Z Y P 2 = (X 1 ,Y 1 ,Z 1 +L) P 1 =(X 1 ,Y 1 ,Z 1 ) 112

CONCLUSION

The BEM has been used to solve heat transfer equation during the cooling step of the molding process, for a 2D and 3D problem. This model is coupled with SQP optimization algorithm to find the best geometrical parameters according to a cost function. SQP is a mono-objective algorithm but it is possible to define one of the objectives like non-linear constraints. It is crucial to study the parameterization of the channels to reduce optimization variable and thus CPU time. The use of the elliptical cooling channels should be further investigated for more complex shapes in future works.

FIGURE 1 .

 1 FIGURE 1. Schematic temperature evolution in injection mold during several cycles.

Figure 2

 2 Figure 2 displays the boundary conditions on the mold.

FIGURE 2 .

 2 FIGURE 2. Boundary conditions for the mold. The temperature of the coolant is T C and the heat transfer coefficient, h, is related to the coolant flow rate (via Colburn correlation coefficient).

FIGURE 3 .

 3 FIGURE 3. Heat transfer simulation coupled with optimization algorithm.

FIGURE 4 .

 4 FIGURE 4. Definition of the geometry parameters.

FIGURE 6 .

 6 FIGURE 6. Temperature at the surface of the 2D mold cavity before and after optimization.

Figure 6

 6 Figure6displays the temperature distribution along the mold cavity surface before and after optimization. We observe that both temperature variance and temperature average decreased significantly. Figure7shows the temperature distribution inside the mold temperature.

FIGURE 7 .

 7 FIGURE 7. Temperature distribution after optimization for the 2D approach.

FIGURE 8 .

 8 FIGURE 8. Schema of plastic part.

FIGURE 9 .

 9 FIGURE 9. Channel geometric parameters.

FIGURE 10 .

 10 FIGURE 10. Objective function versus iterations.

  Figure 11 displays the temperature gradient of the surface of mold cavity. The results are similar to 2D geometry, except for the boundary effect.

FIGURE 11 .

 11 FIGURE 11. Surface temperature distribution at the surface of the mold cavity.

TABLE 1 .

 1 Thermo physical parameters for the mold and polymer.

		Polymer (PP)	Mold (steel
			40cmd8s)
	λ [ W.m -1 .K -1 ]	0.63	34
	ρ [Kg.m -3 ]	891	7800
	Cp [J.Kg -1 .K -1 ]	2740	460

TABLE 2 .

 2 Process

		parameters
	T°C of injection	T°C of ejection	Time of injection =
	= 240°C	= 100°C	20 s

TABLE 3 .

 3 Numerical parameters for 2D geometry.

	Elements type	Linear
	Number of elements	1013
	Number of nodes	1014
	Number of internal nodes	489

TABLE 4 .

 4 Channels positions before and after optimization.

		Initial values	Final values
	Channel « i »	X	Y	X	Y
	1	0.135 0.018 0.129	0.014
	2	0.11	0.018 0.106	0.014
	3	0.067 0.018 0.063	0.014
	4	0.05	0.018 0.046	0.0198
	5	0.037 0.049 0.033	0.045
	6	0.067 0.049 0.063	0.05
	7	0.1	0.049 0.0984	0.0459
	8	0.14	0.049 0.136	0.0496

TABLE 5 .

 5 CPU timed for the direct 2D computing and the optimization.

	Meshing	Matrix A construction	Solving the system	Internals points computing
	0.75 s	1.2188 s	0.64063 s	11.0313 s

The Table

5

displays the CPU time repartition for a direct computing. The matrix construction represents 45% of the direct computing. FIGURE 5. Objective function versus iterations.

TABLE 6 .

 6 CPU time for the 2D optimization problem

	Optimization time	Iterations	Number of cost function evaluations
	1320 s	17	367

TABLE 7 .

 7 Meshing parameters for 3D geometry.

	Elements type	Linear
	Number of elements	3266
	Number of nodes	3204

TABLE 8 .

 8 channels positions before and after optimization.

		Initial values	Initial values
	Channel	X	Y	X	Y
	1	0.135	0.018	0.131	0.021
	2	0.11	0.018	0.101	0.02
	3	0.067	0.018	0.058	0.02
	4	0.05	0.018	0.046	0.019
	5	0.037	0.049	0.033	0.045
	6	0.067	0.049	0.063	0.05
	7	0.1	0.049	0.0984	0.0459
	8	0.14	0.049	0.136	0.0496
	Table 9 displays the CPU time for a direct
	computational. The matrix construction represent 80%
	of the CPU time for the direct computing.

TABLE 9 .

 9 CPU time for the direct 2D computing and the optimization.

	Meshing	Matrix A construction	Resolution system of the	CPU total
	0.75 s	72.21 s	13.295s	86.145 s

Table 10

 10

displays the CPU time optimization.

TABLE 10 .

 10 CPU time for the 2D optimization.

	Optimization time	Iterations	Evaluation's number of the cost function
	13.5 H	24	576

ACKNOWLEDGMENTS

This study was conducted within the frameworks of EUROTOOLING 21 (www.eurotoolin21.com, IP 505901-5). NOMENCLATURE M, P and C indexes refer respectively to the mold, the polymer part and the coolant. C: Heat capacity [J.Kg -1 .K -1 ] d: Cooling channel diameter [m] h: Heat transfer coefficient [W.m -2 .K -1 ] S: Surface [m 2 ] q: Temperature gradient [K.m -1 ] R: Thermal contact resistance [K.m 2 .W -1 ] Tc : Temperature of the coolant [s] T: Temperature [°C] V: Volume [m 3 ] λ: Thermal conductivity [W.m -1 .K -1 ] Φ: Heat flux [W.m 2 ] ρ: Mass density [Kg.m -3 ]


Maxime Bordival

Fabrice M Schmidt

Elian Yannick Le Maoult

E Coment

Measurement Of Thermal Contact Resistance Between The Mold And The Polymer For The Stretch-blow Molding Process

Keywords: stretch-blow molding process, thermal contact resistance, inverse problem PACS: 44.10.+i, 44.40.+a, 44.20.+b

may come

Measurement of thermal contact resistance between the mold and the polymer for the stretch-blow molding process

INTRODUCTION

During a stretch-blow molding cycle, the PET preform is successively blown and maintained into contact with an aluminum mold, thanks to air pressure. The polymer is cooled down by the die, and it rigidifies inside the cavity. In such process, the solidification time of the part is controlled by the heat transfer between the plastic and the mold. The PET cooling rate affects not only the process efficiency, but also mechanical properties of the final bottle. For these reasons, the PET cooling rate prediction is of prime interest, and can not be carried out without a precise understanding of the heat transfer properties. Owing to the imperfect nature of contact between plastic part and mold inner surfaces, the heat transfer is affected by a Thermal Contact Resistance (TCR). The TCR is generally used in numerical software in order to impose the heat flux density boundary condition. Consequently its value is critical in the aim of achieving accurate simulations. However this parameter is often badly estimated. A few researchers have reported values of TCR at plastic-metal interfaces from steady state experiments [START_REF] Fletcher | Experimental Heat Transfer, Fluid Mechanics and Thermodynamics[END_REF]. However, the TCR in Stretch-Blow Molding Process (SBMP) is an unsteady phenomenon, and one cannot simply extend steady state results to such a transient process. Evolution of TCR with time was abundantly studied for different processes like casting die [START_REF] Dour | [END_REF][START_REF] Hamasaiid | MCWASP Conference Proceedings XI[END_REF], or injection-molding [START_REF] Bendada | [END_REF][5][START_REF] Massé | Les couplages thermomécaniques lors de la solidification de matériaux polymères[END_REF]. Some researchers have proposed analytical model allowing to link TCR to the pressure [START_REF] Massé | Les couplages thermomécaniques lors de la solidification de matériaux polymères[END_REF]. However in such processes, the nature of contact is quiet different since the material is molten. Moreover applied pressures and temperature ranges are strongly higher than for the SBMP. In the present work, we have investigated the efficiency of the method exposed in [START_REF] Dour | [END_REF][START_REF] Bendada | [END_REF] for the SBMP. For that, specific sensors have been developed in order to measure at the same time, the evolution of the TCR and the air pressure inside the preform. Main difficulties related to this measurement are the very short time of the phenomenon (around 0.5 s) as well as the small temperature increases which must be detected (typically some degrees for the mold).

MEASUREMENT STRATEGY AND PROCEDURE

An imperfect contact between two materials gives rise to the so-called Thermal Contact Resistance, characterized by a sharp drop in temperature at the interface. It may be defined per surface unit as:

ϕ Mold PET T T TCR - = (1)
where T PET and T Mold are respectively polymer and mold surface temperatures, and φ is the heat flux density exchanged. The TCR (in m².K.W -1 ) is also the inverse of the heat transfer coefficient h. In order to compute the TCR, the strategy adopted in this study is to determinate φ, T PET , and T Mold .

Measurement Procedure for the Polymer Surface Temperature

The strategy adopted for this measurement is to use an infrared pyrometer, inserted through the mold thickness, in middle-height. This type of measurement has the advantage of being non-intrusive. However, it can't be suitably performed without a precise knowledge of PET optical properties. In order to ensure that the pyrometer measures a surface temperature, the PET must behave like an opaque body. For that, the pyrometer spectral range must be chosen judiciously. Moreover the PET mean emissivity over the considered spectral range is required for the sensor calibration.

In this study, optical properties of PET Tergal F9 have been measured with a Perkin Elmer FT-IR spectrometer, over the range 2.5-25 µm. The method consists in determining both spectral specular reflection coefficient ρ λ , and transmission coefficient τ λ . Reflection measurements are performed on thick samples (2 mm for instance), while transmission ones are realized for different thicknesses. At the thermal equilibrium, the spectral emissivity ε λ is deduced form the Kirchhoff law:

λ λ λ τ ρ ε - - = 1 (2)
For an opaque body (τ λ = 0), the spectral emissivity can be easily determinated by measuring the reflection coefficient. Finally, the mean emissivity ε defined as following can be computed:

∫ ∫ = ∆ ∆ λ λ λ λ ε ε λ λ λ d T I d T I T PET PET PET ) ( ) ( ) ( o o (3)
where I λ °(T PET ) is Planck's Intensity [START_REF] Modest | Radiative Heat Transfer[END_REF] at PET temperature. Full details about this procedure are reported in [START_REF] Monteix | [END_REF]. Results are illustrated FIGURE 1. FIGURE 1A reveals that PET presents successively transparent and opaque regions between 1 µm and 20 µm. The range 8-14 µm is particularly interesting since the transmission coefficient is less than 5 %, even for relatively thin samples (0.24 mm). For this reason we have chosen to compute the mean emissivity over this range. Results are illustrated FIGURE 2A. We observe that for a thickness higher than 0.3 mm, the mean emissivity reaches an asymptotic value equal to 0.93. This characterizes an opaque body for which the emission is a surface phenomenon.

As a conclusion these results demonstrate that a pyrometer operating in the range 8-14 µm will measure the PET surface temperature, provided that the preform thickness is sufficiency large. For this study, a pyrometer Optris ® CTfast 8-14 µm has been chosen (FIGURE 2B). The sensor response time, given by the constructor, is equal to 17 ms, which is sufficiently small regarding to the typical time of contact (0.5s). The second sensor developed for this application (FIGURE 3A), is inserted at the same height position as the pyrometer, but on the opposite side of the mold (FIGURE 3B). Its external dimensions are the same as a Kistler pressure sensor to make its installation in dies easy. It is realized in the same aluminum as the mold in order to avoid thermal disturbances. Therefore, the heat flux is considered as unidirectional through the mold thickness, assumption necessary for the inverse method. The measurement is realized thanks to two K-type micro-thermocouples, brazed into the gauge at different depth, and perpendicularly to the heat flux direction. Critical design features that are essential to an accurate estimation of the heat flux are discussed in [START_REF] Dour | [END_REF]. The method is based on the linear inverse heat condution problem (Function Specification Method), fully described in [9].

Experimental Device

Trials were performed on the CROMeP blowing machine (FIGURE 4A), with APT_PACK 18.5 g preforms (FIGURE 4B), injected from PET Tergal F9. The preform heating is ensured by an infrared oven composed of six halogen lamps of nominal power 1kW. The air pressure inside the preform is measured using a Kulite LE 125 sensor (FIGURE 4B). Data are recorded thanks to a Nimtech acquisition system.

RESULTS AND DISCUSSION

During the blowing step (without any stretch rod), the air pressure inside the bottle follows typical variations. In the first time, the pressure increases sharply. As soon as the pressure is sufficient to blow the preform, the bottle volume increases and consequently the air pressure drops. While preform internal volume remains constant, the pressure reaches gradually its nominal value. This typical evolution of air pressure gives a good representation of the blowing kinematic, and remains a reference in order to understand variations of measured temperatures. FIGURE 5A illustrates the PET surface temperature. This one starts to increase before the beginning of the blowing. This is probably due to the own emission of the preform, which is already inside the mold. The pyrometer receives radiation from the preform, after different reflexions inside the cavity. At approximately 2 s we observe a second temperature rise of approximately 14 °C. This one is more difficult to explain, and may be due to different phenomena: 1-The temperature increases because of the polymer viscous dissipation, which appears during the stretching. 2-The temperature rise is due to an overestimating of the polymer emissivity, and in other words, to a measurement error. The thickness of the final bottle was measured in middle height. Its value is equal to 0.2 mm. Regarding to FIGURE 1A, the emissivity decrease is relatively small, and should not be sufficient to justify the temperature increase. As a conclusion, we can assume that the temperature increase is probably due to a combination of both phenomena. FIGURE 5B illustrates evolution of the mold surface temperature versus time, computed using the inverse method. We observe that for one blowing cycle, the temperature increase is around 1 °C. For one-shot bottle, this allows to verify the isotherm assumption generally adopted in numerical software to model heat transfer inside the mold. FIGURE 6A illustrates the heat transfer coefficient (inverse of the TCR) versus time. This one drops dramatically in the first second following the peak value, while the air pressure still increases. This phenomenon was observed by different investigators for other processes [START_REF] Dour | [END_REF][START_REF] Bendada | [END_REF]. It is generally admitted that this sudden rise of TCR (FIGURE 6B) is related to the air gap due to polymer shrinkage. However evolution of contact should be quiet different for the SBMP, since PET does not change of phase. A more probable assumption could be that the polymer cooling is slowed down by the air gap between the pyrometer and the bottle. The peak value of the heat transfer coefficient (averaged on five trials) is equal to 225 W.m -².K -1 , which corresponds to a TCR equal to 0.0044 m².K.W -1 . The standard deviation of this peak value is equal to 20 W.m -².K -1 (9% error), which indicates that the repetitiveness is good.

CONCLUSION AND PROSPECTS

Preliminary results demonstrate that response time of sensors is sufficiently fast to measure the variation of the heat transfer coefficient in a very short time (0.5 s), corresponding to real process conditions. Typical values of TCR that were determinated in this study should be helpful for more accurate simulation of the stretch-blow molding process. However measurements have been performed with constant process conditions. The air pressure as well as the preform temperature should have an important impact on the TCR. Future work will focus on the influence of such parameters.

FIGURE 1 .

 1 FIGURE 1. A: Spectral transmission coefficient of PET Tergal F9 for two sample thicknesses -B: Spectral reflection coefficient of PET Tergal F9 for a 2 mm thick sample

FIGURE 2 .

 2 FIGURE 2. A: Planck mean emissivity of PET Tergal F9 versus sample thickness -B: Pyrometer Optris ® CTfast 8-14 µm

FIGURE 3 .

 3 FIGURE 3. A: Heat flux density sensor -B: Location of sensors through the mold

FIGURE 4 .

 4 FIGURE 4. A: CROMeP blowing machine -B: APT_PACK 18.5 g preform and pressure sensor

FIGURE 5 .

 5 FIGURE 5. A: Measured PET surface temperature and air pressure versus time -B: Computed mold inner surface temperature and pressure versus time

FIGURE 6 .

 6 FIGURE 6. A: Heat transfer coefficient and air pressure versus time -B: Corresponding TCR versus time over the range 2.5 -5.3 s

ACKNOWLEDGMENTS

This work war carried out as a community STRP project named APT_PACK, contact n° STRP 505204-1. Authors would thank Logoplast Company for injecting preforms used in this work, as well as E. Coment and V. Lucin for their contribution.