A Wautier

C Geindreau

F Flin

Linking snow microstructure to its macroscopic elastic properties: a numerical homogenization method and its application to 3D images from X-ray tomography

Keywords: homogenization, KUBC, elasticity, microtomography, finite elements, snow

may come

Linking snow microstructure to its macroscopic elastic stiffness tensor: A numerical homogenization method and its application to 3-D images from X-ray tomography

Antoine Wautier, Christian Geindreau,

Introduction

Snow is a highly porous material that exists close to its melting point and exhibits a complex microstructure constantly evolving in time. Because overall mechanical properties of snow are strongly influenced by its density [START_REF] Mellor | A review of basic snow mechanics[END_REF] and the topology of its microstructure [START_REF] Shapiro | Snow mechanics: review of the state of knowledge and applications[END_REF], a correct multi-scale modeling of the mechanical properties of snow is of great interest when it comes to avalanche risk forecasting.

Recently, the generalization of X-ray tomography and its application to snow [START_REF] Brzoska | 3D visualization of snow samples by microtomography at low temperature[END_REF][START_REF] Schneebeli | Numerical simulation of elastic stress in the microstructure of snow[END_REF][START_REF] Kaempfer | A microstructural approach to model heat transfer in snow[END_REF][START_REF] Flin | The temperature-gradient metamorphism of snow: vapour diffusion model and application to tomographic images[END_REF][START_REF] Srivastava | Observation of temperature gradient metamorphism in snow by x-ray computed microtomography: measurement of microstructure parameters and simulation of linear elastic properties[END_REF][START_REF] Pinzer | Vapor flux and recrystallization during dry snow metamorphism under a steady temperature gradient as observed by time-lapse micro-tomography[END_REF] enabled the description of its complex geometry down to the micro-scale. Nowadays, good databases of 3D images for the different snow types described in the international classification [START_REF] Fierz | TheInternational Classificationi for Seasonal Snow on the Groun[END_REF] are available [START_REF] Calonne | 3-d image-based numerical computations of snow permeability: links to specific surface area, density, and microstructural anisotropy[END_REF][START_REF] Löwe | A general treatment of snow microstructure exemplified by an improved relation for thermal conductivity[END_REF].

With precise description of the ice skeleton of snow, its complex behavior can be investigated thanks to the fairly well established physical and mechanical properties of ice [START_REF] Schulson | Creep and fracture of ice[END_REF]. From this point, the macroscopic behavior of snow considered as a homogeneous material can be up-scaled thanks to the use of techniques derived from the homogenization theory based on the assumption of the separation of scales [START_REF] Dormieux | Introduction à la micromécanique des milieux poreux[END_REF][START_REF] Auriault | Homogenization of coupled phenomena in heterogenous media[END_REF]. In recent works, the combination of X-ray tomography imaging, finite elements techniques and ever increasing computing power was used to bridge the gap between the topology of the ice skeleton of snow and its mechanical behavior [START_REF] Schneebeli | Numerical simulation of elastic stress in the microstructure of snow[END_REF][START_REF] Pieritz | From snow xray microtomograph raw volume data to micromechanics modeling: first results[END_REF][START_REF] Srivastava | Observation of temperature gradient metamorphism in snow by x-ray computed microtomography: measurement of microstructure parameters and simulation of linear elastic properties[END_REF][START_REF] Hagenmuller | Numerical simulation of microstructural damage and tensile strength of snow[END_REF][START_REF] Chandel | Micromechanical analysis of deformation of snow using x-ray tomography[END_REF][START_REF] Schleef | Hot-pressure sintering of low-density snow analyzed by x-ray microtomography and in situ microcompression[END_REF]. However while Srivastava et al.

D R A F T

May 26, 2015, 7:39am D R A F T

[2010] used an approach rigorously based on the homogenization theory [START_REF] Van Rietbergen | Direct mechanics assessment of elastic symmetries and properties of trabecular bone architecture[END_REF], most of the up-scaled snow properties found in the literature were obtained thanks to numerical simulations of "real" experiments assuming an isotropic macroscopic behavior [START_REF] Hagenmuller | Numerical simulation of microstructural damage and tensile strength of snow[END_REF]?]. Moreover, while snow microstructure is generally anisotropic [START_REF] Calonne | 3-d image-based numerical computations of snow permeability: links to specific surface area, density, and microstructural anisotropy[END_REF][Calonne et al., , 2014]], most of the constitutive equations for snow mechanical behavior given in the literature are only unidimensional [START_REF] Mellor | A review of basic snow mechanics[END_REF].

In the wake of the study performed by [START_REF] Calonne | Numerical and experimental investigations of the effective thermal conductivity of snow[END_REF][START_REF] Calonne | 3-d image-based numerical computations of snow permeability: links to specific surface area, density, and microstructural anisotropy[END_REF]Calonne et al. [ , 2014] ] on the homogenization of the transfer properties of snow, the aim of this paper is to propose a rigorous multi-scale modeling of the elastic mechanical behavior of snow and to reckon the relative influence of the density and the topology of the microstructure on its macroscopic properties. Quantitative results are given for a very wide range of densities (from 100 to 600 kg/m 3 ) and the influence of different types of metamorphism is estimated.

In order to do so, boundary value problems derived from the homogenization theories were solved numerically on 3D snow images obtained by micro-tomography. Depending on the choice of the boundary value problem, the sample's representative elementary volume (REV) required to deduce the macroscopic mechanical properties varies. Three boundary value problems described in the homogenization theories are considered to give relatively small REV. From the biggest to the smallest ones [START_REF] Kanit | Determination of the size of the representative volume element for random composites: statistical and numerical approach[END_REF], these three approaches are named SUBC (statically uniform boundary conditions), KUBC (kinematically uniform boundary conditions) and PBC (periodic boundary conditions). In the KUBC approach, a macroscopic homogeneous strain is imposed on the boundary. Conversely, in the SUBC approach a macroscopic homogeneous stress is imposed on the

D R A F T May 26, 2015, 7:39am D R A F T X -5
boundary. In the PBC approach, mix stress/strain boundary conditions are applied by enforcing the periodicity of the displacement field over the three directions of space and the periodicity of the normal stress across the sample boundaries. If the PBC approach is supposed to be the best approach in terms of convergence with respect to the size of the REV [START_REF] Kanit | Determination of the size of the representative volume element for random composites: statistical and numerical approach[END_REF], its application on a non-periodic porous microstructure is not straightforward. This requires for instance to wrap the sample into a virtual bounding box or to fill up the porosity with a soft material. In order to avoid the introduction of artifacts, the KUBC approach was retained.

KUBC numerical homogenization procedure

The whole numerical homogenization procedure can be decomposed into four steps that are summarized in Figure 1 and detailed in the following subsections.

Converting images into 3D meshes (step 1)

Thanks to the use of the Matlab open-source toolbox iso2mesh [START_REF] Fang | Tetrahedral mesh generation from volumetric binary and grayscale images[END_REF],

a tetrahedral mesh was fitted onto the original 3D binary images of typical size 600 3 voxels. The toolbox was used with the cgalmesh utility option and adaptative meshes of controlled density and quality were created from binary snow images. The resulting meshes were composed of a typical number of 2 millions of elements and can be visualized in Figure 1.

Boundary value problem definition in Abaqus (step 2)

Under the assumption of the separation of scales, it is possible to deduce macroproperties of a given material by performing numerical tests on a representative elemen-D R A F T May 26, 2015, 7:39am D R A F T tary volume (REV), i.e. a volume that is "big enough" with respect to the size of the micro-structure to be considered as representative of the whole material and the physical phenomena. Then, thanks to Hill's lemma, the macroscopic stress tensor Σ and strain tensor E are deduced from their microscopic counterparts (σ and ε) by computing their mean values over the REV (V ).

E = ε = 1 |V | V ε dV, Σ = σ = 1 |V | V σ dV
Because of Green's formula, the mean stress and strain can either be averaged on the sample's volume or on its boundary. For instance, if ε is kinematically admissible to a continuous displacement field u(x), the equivalent strain can be computed with the equivalent formula:

E = 1 |V | ∂V u ⊗ s n dS
where ⊗ s is the symetrized tensorial product and n the outward normal vector on the boundary ∂V . With this localization property, a savvy choice in the boundary conditions enables the control of the macro stress or strain that is imposed on a given sample.

And indirectly, this choice in the boundary conditions governs the size of the REV. In the present work, the KUBC approach was retained and practically implemented thanks to the use of the Abaqus plugin HomTools developed by [START_REF] Lejeunes | Une toolbox abaqus pour le calcul de propriétés effectives de milieux hétérogènes[END_REF]. The plugin enables the automatic definition of the linear relationship on ∂V between the displacements u and the macroscopic homogeneous strain E defined as u = E • x for

x ∈ ∂V (Figure 1). From a microscopic point of view, the ice skeleton is constituted of an assemblage of mono-crystalline particles. While some recent studies [START_REF] Riche | Evolution of crystal orientation in snow during temperature gradient metamorphism[END_REF] show that ice grain orientation in snow may be anisotropic under specific circumstances, in the following the poly-crystalline ice skeleton is considered as an homogeneous isotropic elastic material.

Thus, the ice stiffness tensor C ice is defined as :

σ = C ice : ε = E 1 + ν ε + ν 1 -2ν Tr(ε)1 ,
where E and ν are the Young modulus and the Poisson ratio respectively, and 1 is the second order identity tensor. The ice elastic properties presented in [START_REF] Schulson | Creep and fracture of ice[END_REF]; [START_REF] Chandel | Micromechanical analysis of deformation of snow using x-ray tomography[END_REF] show that E ranges from 0.2 GPa to 9.5 GPa. In order to leave aside the resulting huge uncertainties, all the results in this paper are presented in a dimensionless form by setting E equal to one. Based on the literature, the Poison's ratio ν is chosen equal to 0.3.

Computing the homogenized stiffness tensor (step 4)

The homogenization problem to be solved on the REV (V ) is written:

               div σ = 0 for x ∈ V u = E • x for x ∈ ∂V ε = 1 2 ∇ u + t ∇ u for x ∈ V σ = C(x) : ε for x ∈ V , (1)
where C(x) is the local stiffness tensor equal to 0 in the air phase and C ice in the ice phase.

With the knowledge of the microscopic stress, its homogenized macroscopic counterpart can be computed, and the linear relationship between Σ and E gives the homogenized

D R A F T May 26, 2015, 7:39am D R A F T stiffness tensor C hom : Σ = 1 |V | V σ dV = C hom : E.
The homogenized stiffness tensor C hom exhibits only 21 independent components in the most general case. By setting 5 out of 6 of the independent components of E equal to 0 (E = δ ij e i ⊗ s e j where δ ij is the Kronecker delta), only 6 numerical simulations are needed in order to fully compute C hom . However, based on the features of the compliance matrices computed on all the 3D images, it appears reasonable to model snow as an orthotropic material (i.e. a material with two perpendicular planes of symmetry) of principal axes (e 1 ,e 2 ,e 3 ), e 3 being along the direction of gravity. By doing so, the 21 independent components of C hom are reduced to 9: three Young moduli, three shear moduli and three Poisson ratios (because of the symmetry of the compliance matrix) and the inverse of C hom , i.e. the compliance matrix, is written :

       E 11 E 22 E 33 2E 23 2E 13 2E 12        =          1 / E 1 -ν 12/ E 1 -ν 13/ E 1 0 0 0 -ν 21/ E 2 1 / E 2 -ν 23/ E 2 0 0 0 -ν 31/ E 3 -ν 32/ E 3 1 / E 3 0 0 0 0 0 0 1 / G 23 0 0 0 0 0 0 1 / G 13 0 0 0 0 0 0 1 / G 12          •        Σ 11 Σ 22 Σ 33 Σ 23 Σ 13 Σ 12        . (2)
With the use of the Abaqus plugin HomTools, the macroscopic homogeneous stress Σ is automatically computed giving thus a straightforward access to the homogenized stiffness tensor (Figure 1).

Convergence analysis

Because homogenization theory is based on the assumption of the separation of scales, the fundamental notion of representative elementary volume (REV) is discussed in this section, and the integral range theory [START_REF] D R A F T Lantuejoul | Ergodicity and integral range[END_REF] is used in order to link REV and

D R A F T May 26, 2015, 7:39am D R A F T X -9
the correlation lengths presented by [START_REF] Löwe | A general treatment of snow microstructure exemplified by an improved relation for thermal conductivity[END_REF]. A sensitivity analysis is performed on a given snow sample in order to quantitatively define a dimensionless REV size.

Dimensionless volume parametrization

According to the integral range theory developed by Lantuejoul [1991], Kanit et al.

[2003] proposed a method for determining the typical REV size of a random heterogeneous material at a given precision. This prediction is based on the fact that the variance

D 2 Z (V )
of a given physical or mechanical property Z decreases as a power law of the ratio between the sample volume V considered and the integral range A 3 of the material [START_REF] D R A F T Lantuejoul | Ergodicity and integral range[END_REF]. For a porous material of porosity φ, this relation reads:

D 2 Z (V ) = φ(1 -φ)(Z m -Z air ) 2 A 3 V α Z (3)
where Z m -Z air is the property contrast between air and material phases, and α Z is an exponent smaller than one that depends both on the property and the boundary value problem considered. The bigger the edge effect, the smaller this exponent. It is noteworthy to notice that the integral range A 3 is a geometrical parameter that depends on the microstructure [START_REF] D R A F T Lantuejoul | Ergodicity and integral range[END_REF][START_REF] Kanit | Determination of the size of the representative volume element for random composites: statistical and numerical approach[END_REF]. If the two-point correlation function is denoted by

C(X, h) = p{x ∈ X, x + h ∈ X},
the integral range is defined as

A 3 = 1 C(X, 0) -C(X, 0) 2 R 3 (C(X, h) -C(X, 0) 2 ) dh. ( 4
)
By dividing the sample volume by the integral range of the snow considered, it is possible to restrict the volume convergence analysis of the homogenized mechanical properties

D R A F T May 26, 2015, 7:39am D R A F T
of snow to a single arbitrary snow sample and then define a general size for the REV as a multiple (n) of the integral range for any type of snow thanks to equation (3), i.e V = n A 3 . In practice, the numerical estimate of the integral range appears to be fairly well approximated by the product of the three correlation lengths ( 1 , 2 , 3 ) along the axes e 1 , e 2 and e 3 , which can be computed following the approach proposed in Löwe et al.

[2013].

Volume convergence analysis

For practical reasons, the volume convergence analysis was not performed in terms of variance but the dimensionless parametrization presented in the previous subsection was used. In order to evaluate the optimal size of the REV of snow, i.e. the optimal value of n, a volume convergence analysis has been performed on a rounded grain (RG) sample of a density of 285 kg.m -3 and of a typical correlation length = 3 √

1 2 3 3 √ A 3 of 70 µm. Cubic sub-volumes of size V = n A 3 = (m ) 3
were extracted from the raw binary tomographic image and were all meshed with iso2mesh with the same meshing options.

The homogenized stiffness and compliance matrices were then computed following the approach presented in 2.4. The Figure 1 shows the evolution of the Young moduli with respect to m. A similar plot was also obtained concerning the shear moduli.

This figure highlights the convergence of the young moduli towards a constant value with increasing m. A compromise between a good precision and a reasonable simulation time (approximately 10 hours) has been made and a REV of typical length of 30 was retained to perform all the simulations. The Figure 1 also shows that the volume convergence is slower for the Young moduli of the sample than for its density. This is to be put in parallel

D R A F T May 26, 2015, 7:39am D R A F T
with the theoretical expression of the variance given in the equation ( 3) and justifies the approach followed in this paper to not perform the convergence analysis on geometrical properties only.

Mesh convergence analysis

In parallel, a mesh convergence analysis has been performed on the sub-volume of a typical length scale of 30 . The evolution of the Young moduli with respect to the mesh size is given in Figure 1. A very similar plot was also obtained concerning the shear moduli. Again, a compromise between a good precision and a reasonable simulation time has been made and the meshing parameters used for the previous volume convergence analysis were kept.

Correlation between density and homogeneous elastic properties

The homogeneous elastic properties of a large variety of snow types were computed following the above methodology on 29 3D images from the database of the snow research center (CEN) (see e.g. [START_REF] Calonne | 3-d image-based numerical computations of snow permeability: links to specific surface area, density, and microstructural anisotropy[END_REF]). The numerical data used in this study are available in the supplementary material linked to this paper. average value is plotted together with the interval that contains the six Poison's ratio presented in (2). In order to visualize the difference in microstructure of the different samples, color is used in order to highlight the different types of snow according to the international classification [START_REF] Fierz | TheInternational Classificationi for Seasonal Snow on the Groun[END_REF].

In a log-log scale, the Young and shear moduli dependence with respect to density appears to be relatively linear throughout the wide range tested (from 102 to 544 kg.m -3 ).

A power law fit gives the first order dependence relationships with good regression coefficients (R 2 ):

       E snow E ice = 0.78 ρ snow ρ ice 2.33
, R 2 = 0.96

G snow G ice = 0.96 ρ snow ρ ice 2.55 , R 2 = 0.98 , for ρsnow ρ ice ∈ [0.1; 0.6]. (5)
The numerical points reported in Figure 2 were compared with the prediction of two analytical models and two energy bounds. The Young and shear moduli corresponding to the prediction of the classical Reuss's bound and the Hashin and Shtrikman upper bound for isotropic material [START_REF] Hashin | A variational approach to the theory of the elastic behaviour of multiphase materials[END_REF] were computed on the range of densities considered. The isotropic differential and Mori-Tanaka schemes, which can still be used for high porosity values, were also simulated. The differential scheme [START_REF] Zimmerman | Elastic moduli of a solid containing spherical inclusions[END_REF] consists of an iterative application of the dilute scheme that implies the resolution of an Eshelby problem where the true macroscopic strain is applied on the boundary. The Mori-Tanaka scheme [START_REF] Mori | Average stress in matrix and average elastic energy of materials with misfitting inclusions[END_REF] consists of an Eshelby problem where an implicit macroscopic strain is applied on the boundary in order to respect the mean strain distribution between matrix and inclusions. Eventually, in order to evaluate the consistency of our implementation of the KUBC homogenization process to snow, the numerical

D R A F T May 26, 2015, 7:39am D R A F T
estimates of the Young modulus in ? [START_REF] Srivastava | Observation of temperature gradient metamorphism in snow by x-ray computed microtomography: measurement of microstructure parameters and simulation of linear elastic properties[END_REF] and [START_REF] Schneebeli | Numerical simulation of elastic stress in the microstructure of snow[END_REF] were reported in Figure 2 for comparison.

Discussion

As visible in Figure 2, the parametrization (5) of the elastic properties of snow by a power law gives reasonable results for the prediction of the Young and shear moduli of snow on the whole range of densities explored. It is interesting here to draw a parallel with the analytical models classically used for the elastic properties of cellular materials [START_REF] Knackstedt | Elastic and transport properties of cellular solids derived from three-dimensional tomographic images[END_REF][START_REF] Gibson | Cellular solids: structure and properties[END_REF]. Indeed, for skeleton cellular materials locally loaded in bending, the relative Young modulus of the porous material is given by a power law of exponent n = 2 of its relative density. However for membrane cellular materials locally loaded in bending, the relative Young modulus is given by a power law of exponent n = 3 of its relative density. The exponent of the fitted law ( 5) is consistent with the predicted behavior of these materials also characterized by a very high porosity. Under the hypothesis that the ice grains are loaded in bending too, the macroscopic behavior of snow is closer to a skeleton behavior than a membrane one.

Concerning the comparison with the differential and Mori-Tanaka schemes in Figure 2, the differential scheme is found to overestimate the values of the Young modulus but gives a relatively good prediction of the exponent of the power law linking relative Young modulus and relative density. For the shear modulus, the differential scheme is not satisfactory anymore. As for Poisson's ratio in Figure 3 schemes considered were initially developed in order to model inclusion/matrix materials and it is not surprising that they cannot account for the connected porosity of snow.

The results of [START_REF] Srivastava | Observation of temperature gradient metamorphism in snow by x-ray computed microtomography: measurement of microstructure parameters and simulation of linear elastic properties[END_REF] in Figure 2 computed on a limited range of densities appear to be of the same order of magnitude as the results presented in this paper. The results obtained by ? on a wide range of densities are not consistent with our results. The exponential fit proposed by ? is plotted on Figure 2 and a power law with a good regression coefficient of R 2 = 0.917 has been proposed to better fit the data :

E snow E ice = 5.73 ρ snow ρ ice 5.15
.

The major difference between our approach and those found in the literature lies in the boundary value problem used in the homogenization procedure. While we use KUBC boundary conditions, ? uses a uniaxial compression test that requires a much bigger size of VER. In order to deduced the Young modulus value from numerical simulations, ?

were forced to postulate a macroscopic isotropic behavior for all their samples. Further more, the resolution of the images used by ? is only a tenth of the resolution used to design our finite element meshes. These differences in the testing procedure might result in the observed differences in Figure 2.

The last comment to be made on the results shown in Figures 2 and3 is that many samples present an anisotropic behavior that cannot be captured by the density parametrizations of the macroscopic properties of snow given in (5). The degree of anisotropy of a given type of snow can be reckoned thanks to the definition of two indicators A(E) and A(G) such that :

A(E) = E 3 (E 1 + E 2 )/2 , A(G) = G 12 (G 23 + G 13 )/2 , (6)

Influence of the type of metamorphism on the mechanical properties

As already shown by [START_REF] Srivastava | Observation of temperature gradient metamorphism in snow by x-ray computed microtomography: measurement of microstructure parameters and simulation of linear elastic properties[END_REF] and [START_REF] Schneebeli | Numerical simulation of elastic stress in the microstructure of snow[END_REF] for mechanical properties or by Calonne et al. [2014] for thermal conductivity for instance, the metamorphism of snow, even when occurring at constant density strongly influences its macroscopic behavior. The full 3D homogenization approach implemented here enables the quantitative analysis of these anisotropic mechanical properties. In this section, three types of snow metamorphism are considered in order to reckon their influence on the mechanical properties of snow.

Temporal evolutions of snow microstructures

In Figure 4 the time evolution of the three Young moduli and the three shear moduli are plotted for each of the three temporal evolutions considered: a temperature gradient experiment, a wet snow metamorphism experiment and an isothermal experiment. In this Figure the subscript 3 corresponds to the vertical direction. In order to reckon the influence of variations in density with time, its time evolution is also plotted in Figure 4.

Temperature gradient experiment: The first temporal series considered consists of the seven micro-tomographic images of Calonne et al. [2014]. They were obtained during a temperature gradient experiment at -4˚C. A temperature gradient of 43 K.m -1 was applied to a sieved snow layer during 500 h and samples were regularly extracted

D R A F T May 26, 2015, 7:39am D R A F T
from the middle of the layer, impregnated with 1-chloronaphtalene and stored at -20˚C.

They were then imaged by microtomography.

On the corresponding plot of Figure 4, it can be seen that the general fluctuations of the relative moduli follow the small fluctuations in density observed between the different images. In addition to that, the development of an anisotropic mechanical behavior can be seen with the apparition of an increasing gap between the vertical Young modulus E 3 and the two other ones, and between the horizontal shear modulus G 12 and the two other ones.

Wet snow metamorphism experiment: The second temporal series considered consists of 5 micro-tomographic images taken from [START_REF] Calonne | Numerical and experimental investigations of the effective thermal conductivity of snow[END_REF]. They were obtained during a grain coarsening experiment. Snow porosity was saturated with water kept at 0˚C and the coarsening of this wet snow was monitored over 142 h by sampling several specimens and imaging them by X-ray microtomography.

After a fast evolving phase at the begining of the experiment, the Young and shear moduli evolve in an isotropic way.

Isothermal experiment: The third temporal series considered consists of the first 8 micro-tomographic images of [START_REF] Flin | Three-dimensional geometric measurements of snow microstructural evolution under isothermal conditions[END_REF]. 15 h after a snowfall, a natural snow layer was kept under isothermal conditions and the metamorphism of the snow microstructure was monitored over 3 months by sampling several specimens and imaging them by X-ray microtomography.

On the corresponding plot of Figure 4, the time evolution of both the Young and shear moduli follows the increasing density.

D R A F T May 26, 2015, 7:39am D R A F T

Discussion

Temperature gradient experiment: In Figure 4, the fast initial decrease in the values of the Young moduli might be related to the initial sublimation of the small ice bridges between ice grains [Calonne et al., 2014]. Because of the vertical orientation of the vapor fluxes, horizontal bridges are more likely to be broken than the vertical ones. This is to be linked with the lower initial decrease of the vertical Young modulus E 3 . With the lengthening of the ice crystals in the vertical direction combined with vertical air chimneys due to sublimation/deposition phenomena, more ice is involved in vertical compression or shear solicitations than in their horizontal counterparts. This microstructure evolution is thus consistent with the relative increase in the vertical Young modulus E 3 and the vertical shear moduli G 13 and G 23 and with the results of Calonne et al. [2014] concerning the physical and geometrical properties of snow (thermal conductivity, Gaussian curvature, ...). During the experiment, the anisotropy indicators A(E) and A(G) defined in (6) vary from respectively 0.89 and 1.08 at the begining of the experiment to 1.56 and 0.83 in the end. Overall, most of the moduli decrease from the initial to the final state, which is consistent with the well-known fact that temperature gradient metamorphism is a cause of weak snow layer formation.

These observations can be compared to the results obtained by [START_REF] Srivastava | Observation of temperature gradient metamorphism in snow by x-ray computed microtomography: measurement of microstructure parameters and simulation of linear elastic properties[END_REF] for a similar temperature gradient experiment. Their findings are consistent with ours except the fact that they observed an increase in density at the beginning of their experiment. This explains why they didn't observe any initial decrease of the Young and shear moduli.

D R A F T May 26, 2015, 7:39am D R A F T

Wet snow metamorphism experiment: In Figure 4, the first increase in density is to be related to residual liquid water that remained in the sample due to capillary forces and froze before the sample was completely dried. Then the Young and shear moduli remain constant and evolve in a very isotropic way (A(E) ∈ [1.00; 1.12] and A(G) ∈ [0.94; 1.00]). The observed stagnation can be quite puzzling because the microstructure is modified in a noticeable way during this experiment. However this can be explained by the fact that the ice bridges between ice grains are growing and, in the meantime, their number per volume unit is decreasing. Overall, such a metamorphism seems to have little impact on the snow mechanical properties as far as the liquid water has been drained before refreezing. Otherwise, the increase in density due to refreezing would significantly strengthen the structure which is consistent with the practical knowledge of melt freeze effect on the mechanical properties of snow. with time after a snowfall, providing that no temperature gradient metamorphism is at stake.

Conclusion and Outlooks

In this work, a numerical homogenization method that can capture the full 3D macroscopic mechanical behavior of snow was successfully developed. A rigorous use of the theoretical results of the homogenization theory was made in the definition of the boundary value problem to be solved. Contrary to most of the previous studies on the subject, the convergence analysis was performed in a dimensionless form and in terms of physical and not geometrical parameters.

Thanks to this powerful numerical procedure, a much wider range of density and a greater diversity of microstructures has been explored in this paper than in previous studies. Over the whole range of density tested, simple power laws have been proposed in order to link the Young and shear moduli of snow to its density.

Thanks to the study of three temporal series, the influence of the three main types of metamorphism on the macroscopic elastic properties was successfully reckoned. The numerical results obtained are in good agreement with the qualitative effects each metamorphism type has on the mechanical properties of snow.

In the future, the proposed methodology might be easily transposed to the study of other materials, provided that the 3D images of their microstructure is available and that the constitutive behavior of their constituents is known.

4. 1 .

 1 Density influence on the snow elastic properties Young moduli, shear moduli and Poisson's ratio are plotted with respect to the relative density of snow in Figures 2 and 3. Triple "T" shapes are used to show the relative locations of the Young and shear moduli. While the large and the small bars of the triple "T" correspond to horizontal moduli (respectively E 1 and E 2 or G 23 and G 13 ), the bottom of the triple "T" corresponds to vertical moduli (E 3 or G 12 ). For the Poison's ratio, the D R A F T May 26, 2015, 7:39am D R A F T

  the Mori-Tanaka scheme proves to give a relatively good estimate contrary to the differential scheme. The two analytical D R A F T May 26, 2015, 7:39am D R A F T X -14 WAUTIER ET AL.: 3D HOMOGENIZATION OF THE ELASTIC PROPERTIES OF SNOW

D

  R A F T May 26, 2015, 7:39am D R A F T all the samples considered, we have A(E) ∈ [0.56; 1.57] and A(G) ∈ [0.77; 1.47] while a perfectly isotropic material would exhibit anisotropy indicators equal to 1.

  Isothermal experiment: In Figure4, the densification of snow is quantified with the increase in density due to the creeping of snow under its own weight. Because of the relatively small thickness of the snow layer considered (about 10 cm), no huge anisotropy appears with respect to the privileged vertical direction of compaction. During all the experiment, the anisotropy indicators remain relatively constant with A(E) ∈ [0.85; 1.03] and A(G) ∈ [0.95; 1.28]. However, as it has been observed by[START_REF] Flin | Three-dimensional geometric measurements of snow microstructural evolution under isothermal conditions[END_REF], the snow grains tend to rearrange with their longest dimension along the horizontal direction. This might be related to the fact that the vertical Young modulus is systematically the lowest one after 600 h. Overall, such a metamorphism results in a noticeable densification of snow, and thus in an increase of its stiffness. Recent snow on top of the snowpack evolves D R A F T May 26, 2015, 7:39am D R A F T

Figure 1 .

 1 Figure 1. Four steps procedure used in order to transform 3D microtomograph images

Figure 2 .

 2 Figure 2. Young moduli (top) and shear moduli (bottom) dependence with respect

Figure 4 .

 4 Figure 4. Time evolution of the three Young moduli (middle) and the three shear

Acknowledgments. We thank S. Lejeunes and S. Bourgeois from the LMA for sharing with us the code of their plugin HOMTOOLS and for helping us to use it in our scripts.

We thank PK. Srivastava for accepting to share with us the numerical data of Srivastava et al. [2010]. CNRM-GAME/CEN is part of the LabEx OSUG@2020 (ANR10 LABX56).

3SR lab is part of the LabEx Tec 21 (Investissements d'Avenir -grant agreement ANR11 LABX0030).


Albert Marouani
email: marouani@unice.fr

Keywords: University, Institutions, Prospective, Public goods, Markets, Regulation

The University as an object of reflection is itself a very fascinating institution where knowledge is generated, stored, transformed, taught, converted into truth or ideology. These many roles of university have been established in Oxford University, the Sorbonne and Bologna in the 11th century. What has changed? What are the invariants?

What are the roots of the university? How will or how should the university develop to fulfill its roles and to contribute to peace? To answer this question we will use the past as tools for prospecting the future and a balance between vision and reality. « As an institution, (the university) looks far into the past and far into the future, and is often at odds with the present » Clark Kerr (1994 p.14).

1.

Keep the memory of the core values of the university of origins: community, autonomy and freedom.

2.

Overcoming the challenges of the present.

3.

Reinventing the global university of the future.

The fundamentals values of the past.

The university has a long and persistent history that are still relevant today. In particular those related to the defense of intellectual autonomy, the dynamics of teaching combined with research, academic freedom and the nurturing of the humanist tradition. The Latin word "universitas" means an autonomous group which has been called "guild", the members of which has common interests. Among the many medieval guilds developed in a society of people, teachers and students who decided to look at the world in a sholarly fashion. They established the only guild that, to this day, has retained the label "universitas" to define its community. Prominent universities became "studium generale", not because of their teaching was wide, but because of their students came from many different places; their teachers also licensed to teach in different institutions.

According to Wilhelm von Humboldt, universities are institutions, centers of scholarship where teachers and students work together in the pursuit of knowledge. It is both the task of dealing with knowledge and its setting, and the dynamics of combining research and training that make the university a very special community, unique and specific. Both Humboldt and Whitehead, a century later, claim that knowledge, imagination and their dynamic combination are of fundamental value for the academic community that makes the university.

The message coming from different authors, German or American historians and philosophers, is quite clear: the university is a community of teachers and steadfastness is the acquisition, preservation and transmission of knowledge in the service of mankind. For instance Georges Gusdorf, a French philosopher, writes: "The essential aspects of the university as they come from their earlier history are community and interdisciplinarity.

The corporation unites masters and students in the enterprise of teaching; in the service of the intellectual values. I do believe that this idea of the university is not outdated and we have in the morning twilight of the 21st century to reinvent these fundamental values.

To challenge the present : The pressure of the global market

We can distinguish four different tradition in modern University:

• The Napoleonian model of the French Universities considers the university as a tool for society's development, closely linked to state needs. The universities are controlled centrally and their autonomy is a real challenge for all universities in France.

• The German model is under Humboldtian influence. The universities enjoyed considerable independence for their own administration. Academic freedom in research and teaching is a strong feature of this system and non-academic considerations do not interfere with scholarly work.

• The English model come under the sway of Newman's ideals. The university should offer the most suitable environment for learning by close interaction among students and tutors.

The university is governed by the academics themselves within a collegial community of learning.

• The American model is a mix of German and English models. The University has a basic mission to the community it belongs to. There is a clear distinction between undergraduate and graduate education. The university is responsible for the professional management of faculty members and the faculty of the faculty responsible for study programs and the organization of research. Despite the differences between the university and the world, we can notice a convergence towards the American model as an effect of the globalization process. All the Universities in the world emphasize the importance of academic freedom and the role of the academic community in running the institution. University tends to be concerned both for its fruitful advances and for the retention of antiquated traditions like collegial governance for instance. The Universities receive contradictory messages from civil society and industry: they are strongly incited to the frontiers of science by competing on scientific score listings, They are here to help you by working in the field of business, by providing first class future generations fundamental scientific knowledge, philosophy, history, political and moral and critical thought.

Facing these strong and conflicting demands made by a multitude of stakeholders, they seem to be uncertain about their roles, obligations and loyalties, searching for their path into the future, wondering which way to go. The professors and researchers know or feel that they are in their own right.

Higher Education is both a public and a private good. Higher education, as service industry is included in the WTO trade negotiations extending free trade to services. A growing number of universities are selling more education similar to any other commodity. But on the other hand, as a human capital higher education contributes to growth and competitiveness of a nation and can be considered as a public good. It is obvious that today science and knowledge, in the abstract are also a global good. However, there is a general tendency to trust the market. But I am convinced that public research universities can compete with the private ones. Even if many governments encourage public universities to adopt the market, they are private companies and thus depriving the academics power. It is not a commodity for a competitive market, but it is not a business activity. However, public universities must hear the demand of the civil and economic society.

Other tendencies must be considered to understand the evolution of the university: The emergence of the natural sciences and engineering within the university and their relevance for economic development to the gradual subordination of other university disciplines due to their overwhelming scientific and economic strength. Progressively they invaded and took over the performance and relevance of universities. The decline of humanities and social sciences is now a general concern for all universities in the world. There is a growing tension between a professional and an academic or theoretical emphasis in education and research as well as scientific research. This tension is between governments, industry and students who immediately succeed, and on the other side, the institutional academic and incentive schemes are driven both by the desire of High status and institutional ranking and by research and tenure plans.

In the context of globalization we can assume that their influence and domination over the universities and at the same time to reduce their commitments to universities by pushing them towards the market. This contradiction is now a dilemma for governments and politics at national level. Even if they believe in their policies, we doubt that they are going to the extent they think. In France for instance, the actual government asks the universities to care for the needs of their students in terms of employability, and indeed, universities have increased the relevance of their service dimension, in a visible and transparent way .

Today all the universities in the world find them in a paradoxical situation: they have a stencilled attitude, trust, consideration and respect. They have received a lot of rhetorical support from both governments and industry, which is of the utmost importance in the knowledge of the world, and is of fundamental importance to the world. new age of welfare, freedom and democracy. The universities are considered in this context as a major social institution and are expected to play a role in social roles. On the other hand is a little bit different. Universities cannot afford to be in a position to cope with the challenges of the modern world: environment, disease, economic crisis, inequality, poverty, social and urban violence, political and military conflict, etc. The Universals of Existing Economies, Existential Uncertainty, and Fiction. In other words the universities consider they are not moral nor the financial backing they feel the official rhetoric does imply.

Universities are able to demonstrate their dynamism, their sensitivity to stakeholders. They know that they are accountable to society and to the government that decides on the allocation of public funds. The culture of internal and external evaluation, accreditation and quality assurance is growing under pressure from public authorities and different stakeholders. Most universities around the world see the need to position themselves in a very competitive market. They are pushing for the production of more scientific articles, more registered patents. At the same time, they tend to brag about the number of students and graduates and become more independent of public funds by developing new channels of communication with civil society and in particular with the world of industry (academic foundations and partnership). Universities feel compelled to seek research and training funds from industry, mainly because of the results of vocational training and research.

There are many different types of universities. Clark Kerr's multiversity concept suggests that there are different types of universities and many types of universities around the world. However, despite their differences, they share something important. Each university has a unique local profile, a very specific environment and a particular mission in a sociocultural context. But we insist that academic freedom belongs to a social entity that produces, transforms, teaches, and studies science and engages in an ethical transfer of science and scholarship from one generation to the next. other. We are convinced that the increase of freedoms, once practiced within the institution, has repercussions on its members and all its stakeholders. We believe that universities are improving ethical standards, and that they are gaining in importance, and we believe that their authority is becoming more important. Academic freedom can be used and abused. We cannot reduce academic freedom to personal and individual issues. Universities have a moral obligation to ensure the protection and assistance of scholars to the victims of persecution and to fight against unethical methods leading to various forms of torture and manipulation of human rights. Universities are responsible for protecting their freedom and their voice. They need allies in civil society, political circles and the media. Security issues, national boundaries, marginalized groups of disciplines, and so on. Internal and external corruption undermines institutional autonomy and academic freedom. Academic freedom is not an easy ingredient of institutional politics, neither for politicians nor for all members of each university.

Academic freedom is a permanent challenge for all, more importantly, for the institutions themselves, who are the main providers of freedom in the political arena. Indeed, the very special freedoms needed to run a university need protection and care.

Reinventing the future: sharing, cooperation and solidarity

The future of the University (likely, or possible or preferable). What will be or what should be the future? Should we separate research from teaching? Should graduate studies be separated from undergraduate studies? Should universities evolve under the only pressure of globalization and market forces? Many claims that universities are slow and cannot keep pace with economic and social change. This may be the case, but the research needs a long-term perspective. Similarly, education and training must have time to mature. Could we receive in another way the services we receive today from universities? Could stakeholders find alternative partners elsewhere? In the Middle Ages, universities were weakly linked to a given territory and, when the situation became uncomfortable or dangerous, the university could easily migrate to another city. Today, in the age of the digital revolution, the university can be virtual. But "MOOCs" will never replace the vital need for physical proximity and human encounters on a campus as a particular place of study and located geographically in a space steeped in history.

Research, to a large extent, can actually be done outside the walls of universities and, in fact, is often done in public or private research units. But only universities can bring together large, creative and inventive scientific communities on a global scale. The university remains the best system for accumulating knowledge and disseminating new ideas through teaching. Doctoral research in universities ( or some High Schools) remains the most economical way to mobilize an intellectual workforce. Even if there are cultural differences between countries and universities, we must build the future of the University on the basis of common universal values based on a common ethic of sharing. The university's constant concern is to achieve universal knowledge, to transcend geographic and political boundaries and to affirm the need for different cultures to enrich one another.

According to Magna Charta Observatory:

1. The University should be an autonomous institution. To meet the needs of the world around him, his research and teaching must be morally and intellectually independent of all political authority and economic power. What are the common social missions that each university must provide?

We can identify four categories of university social mission within its main function: to cultivate learning through science education and research.

1. A cultural and academic mission including fundamental knowledge and a scientific culture.

2. An economic mission that involves professional training and research for innovation.

3. A political mission that ensures the education of free citizens through the dissemination of information and critical knowledge for democratic debates.

An egalitarian mission that ensures equal opportunities for all students, including through a liberal access policy or lifelong learning.

Academic values such as intellectual honesty, personal integrity, solidarity, cooperation, etc. should not be taken for granted in the construction of the university of tomorrow. One of the main concerns of the universe is to regain self-confidence and ensure the trust of society. Will be established in the field of research and education. All universities will learn to balance competition and cooperation ("coopetition").

In 2000, the student population of higher education reached more than 100 million worldwide. It increases by about 4% per year. If we consider that the world population increases by 1% per year, we should reach 200 million students worldwide by 2020. The number of universities should increase in the same proportion. The global academic landscape will change with a shift to the East (Asia), the South (Africa, Latin America) and the countries of the MENA region. The issues of interuniversity cooperation must be considered by European and American universities to support this movement and strengthen the research capacities of universities in the East and South. We believe that the University should be considered a global public good that must be sheltered and outside nationalist political conflicts.

A fortiori, racism, anti-Semitism, communitarianism, and culturalism under their most diverse and insidious forms have not, or should not, have their place in the university. These principles must be affirmed with force, especially in all the countries around the Mediterranean and the Middle East. A reminder of the thousand-year history and ideas that have guided the development of Mediterranean universities should be used to chart the way forward to strengthen solidarity and community among all universities in the Mediterranean basin without exclusion and discrimination.

Between a highly centralized and state-controlled public university system and a marketbased private higher education system, there is room for a decentralized higher education system. We have to return to the model of middle-age universities governed by communities of professors and students making decisions in a collegiate manner and benefiting from both public and private funding.

Conclusion

Universities must be proactive in designing their strategic directions. They must be actresses of their own destiny rather than passive subjects depending on the market or the state. Autonomy and academic freedom to be effective must be integrated into the democratic life of the university.

All the decision-makers of the university, all the actors of the university, inside and outside the institution, all those who care about the future of the university should defend individually and collectively the values academic and social issues such as equality, democracy, critical thinking, preservation of memory, the fight against ignorance, etc.

In my opinion, the "raison d'être" of the future of academic identity is also to contribute to peace between nations through academic freedom, institutional autonomy, international scientific cooperation between free institutions and free scholars. We firmly believe that freedom, the privilege of higher education institutions, is the way to achieve permanent peace among all Mediterranean countries.  Bok, D(2003). Universities in the marketplace:the commercialization of higher education. Princeton, N.J; Woodstock:Princeton University Press.

Bibliography

2.

  Education and research must be inseparable.3. Freedom of research and training is a fundamental requirement of the university.4. The rejection of intolerance and openness to dialogue must be shared by teachers, researchers, students and all actors.



  Brint, S. (2007). Can public research universities compete? in R.L Geiger, C.L.Colbeck,R.L Williams & C.K. Anderson (Eds), Future of the American public research university (pp.91-118). Rotterdam: Sense Publishers.  R.L Geiger, C.L.Colbeck,R.L Williams & C.K. Anderson (Eds), Future of the American public research university. Rotterdam: Sense Publishers.  Graham, G.J (2005). The insititution of intellectual values:realism and idealism in higher education. Exeter: Imprint Academic.  Gusdorf, G. (1964). L'université en question, Paris: Payot.  Humboldt, W.v. (1810/1970). On the spirit and the organisational framework of intellectual Insitutions in Berlin. Minerva, 8, 242-250.  Kerr, C. (1994) The Uses of the University (4th ed) Cambridge, MA: Harvard University Press.  Kohler, J., & J. Huber (2006). Higher education governance between democratic culture, academic aspirations and market forces. Strasbourg: Council of Europe Pub.  Morin, E. (2001). Seven complex lessons in education for the future. Paris: UNESCO Pub.  Neave, R. & T. Nybom (Eds), The european research university: an historical parenthesis: essays in honour of professor Stig Strömholm , Basingstoke: Palgrave Macmillan.  UNESCO/OECD World Education Indicators Programme (2005)  Watson, D. (2008) The university in the modern world:Ten lessons of civic and community engagement Education citizenship and social justice, 3(1), 43-55.  Whitehead, A.N. (1929/1949) The Aims of Education and other Essays. New York: New American Library.  Inventing Tomorrow's University, What is to Take the Lead? (2008.) An Essay of the Magna Charta Observatory. Bolognia University Press, Bologna.


André Galligo
email: galligo@unice.fr

Frédéric Lesage

Sebastian Minjeaud
email: sebastianminjeaud@unice.fr

Frédéric Lesage

Sebastian Minjeaud Modeling

S Barland

P Reynaud- Bouret

F Cauneau

; K Guillouzouic

U Kuhl

T Passot

F Planchon

Galligo André

Lesage Frédéric

Minjeaud Sebastian

Frederic Lesage
email: frederic.lesage@unice.fr

Modeling and Computation of a liquid-vapor bubble formation

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

Introduction

Boiling is a complex phase transitioning process in which a new liquid-vapor interface (bubbles) is created. It is an effective mechanism widely used in energy conversion industrial facilities. However, a precise descriptive/predictive model of bubble formations is needed to better understand its heat transfer characteristics. As an important first step, the question can be addressed from a quasi-static viewpoint. Indeed, at the interface, a classical conservative law governs the surface curvature of the generated bubble. This non-linear law links the difference of pressure and the surface tension effect. After normalization, it can be expressed at each point of the surface by the following equation:

κ 1 + κ 2 -Bz =constant;
where B is a Bond number, summarizing the context, z is the elevation and κ 1 + κ 2 is twice the Gaussian mean curvature. Therefore, the study of a growing bubble is amenable to solving geometric computations and non-linear differentiable equation. The main difficulty is to design criteria to detect when the growing bubble will detach from the substrate's nucleation site. Inspired by experimental observations, the second author was able to develop [START_REF] Frédéric | Experimental and numerical analysis of quasistatic bubble size and shape characteristics at detachment[END_REF]2] a successful computational strategy. The computational model implemented for the axis-symmetric case (which relies on 2D geometric computations), provides useful information, fully confirmed by experimental observations. Moving forward, we were able to speed up the computations and introduce the use of splines for representing the plane section of the surfaces.

In future work, we will address with similar tools the case for which the nucleation site is on an inclined plane; this geometry better models the one effectively used in industrial plants. Our next target is to better formalize our detachment criteria, in collaboration with an experimental team in Dublin, Ireland; then develop and implement a 3D model for a more general situation. Future studies will also include external stress terms such as cross flows and electrical fields in an effort to move towards realistic and accurate boiling models.

Young-Laplace Equation and Capillary Equation

We consider an infinitesimal section of a vapour formation interface defined by two principal radii of curvature in perpendicular planes, as illustrated in Fig. 1.

We suppose that the work done on the surface is considered to be equal to the work attributed to the pressure difference, as expressed by the Young-Laplace Equation:

∆ p = σ (κ 1 + κ 2 ) (1)
in which κ 1 and κ 2 are the curvatures corresponding to the principal radii of curvature in orthogonal planes, and σ is defined as the work per unit area required to produce an area variation of the interface. Applying the Young-Laplace pressure balance to the context of a single vapour formation event, we consider a surface that is axisymmetric. The surface can thus be described by a smooth curve C that defines the two-phase interface, as illustrated in Fig. 1. The horizontal coordinate is denoted by r and the vertical one, oriented downward from the apex, by z. We normalize the problem by imposing a foot radius of unity, in other words the point (η, 1) belongs to the curve C.

The curvature κ 1 is in the zr plane and the curvature κ 2 is in the plane orthogonal to the tangent unit vector T. Due to symmetry, the curvatures are equal at the apex origin, they are denoted by κ 0 = κ 1 = κ 2 . Let g be the gravitational constant and ∆ ρ the difference in density between the two phases, (1) implies the capillary equation:

2σ κ 0 -∆ ρgz = σ (κ 1 + κ 2 ), on C.
(2)

Introducing a dimensionless Bond number, Bo = ∆ ρr(η) 2 g/σ , we get the adimensionalized form of the capillary equation:

2κ 0 = κ 1 + κ 2 + Bo • z, on C.
(3)

Strategy for Solving the problem

By choosing a parametrization of the curve C and explicit expressions of the curvatures, the equation (3) can be expressed as a second order Ordinary Differential Equation (ODE) with a free parameter κ 0 (the Bond number being determined by the physical properties of the two fluids). In order to simplify the bubble detachment problem, we decompose the difficulty by fixing a positive value for the parameter η and we consider the corresponding subproblem with the following boundary conditions:

(0, 0) ∈ C ; (η, 1) ∈ C and dz dr = 0 at the apex.

Numerical procedures allow to compute a candidate curve solution for each value of the fixed parameter η. Roughly speaking, two boundary conditions are needed to determine the curve from the second order ODE (3) while the third one determines a value for the free parameter κ 0 . So, we obtain a familly of curves depending on the parameter η. Then we need a criterion to select the curve that best models the (physical) bubble detachment phenomena. We claim that the detachment profile is identified by noting that at a given height η, below a threshold, one smooth curve satisfies the boundary conditions (4) with at most one inflexion point. This is illustrated in Fig. 2.

The detachment height is then deduced through incremental increases in the parameter to maximize the height η till this property can be satisfied.

First numerical procedure

This first method is an established approach [START_REF] Frédéric | Experimental and numerical analysis of quasistatic bubble size and shape characteristics at detachment[END_REF]3]. The half (smooth) curve C, is divided in two parts: first a graph of a function r(z) near the apex, then a graph of a function z(r) till the foot. Expressing the curvatures κ 1 and κ 2 in each case, relying on differential geometry formulas, gives rise to a pair of ODEs:

2κ 0 = z (r) (z (r) 2 + 1) 3/2 + z (r) r 1 + z (r) 2 + Bo • z(r) (ODE1) and 2κ 0 = -r (z) (r (z) 2 + 1) 3/2 + 1 r(z) 1 + r (z) 2 + Bo • z (ODE2)
For each η, we first choose a value of κ 0 . We integrate (ODE1) with the initial boundary conditions z(0) = 0 and z (0) = 0. This is done till some value r 1 of r such that z (r 1 ) = 1. We denote by z 1 = z(r 1 ). Then we integrate (ODE2)

with the initial boundary conditions r(z 1 ) = r 1 and r (z 1 ) = 1. This second integration is done until z = η. The r(η) is computed, it will serve for a shooting method. Indeed, performing this process for different values of κ 0 , allows to find a minimum value of κ 0 for which r(η) = 1. This numerical approach was successfully implemented using the software Mathematica.

Second numerical procedure

The half (smooth) curve C is represented by a kind of arc-length parameterization (0 ≤ t ≤ 1): r(t), z(t) with |r (t)| 2 + |z (t)| 2 = λ 2 , depending on a new parameter λ . Again, we can express κ 1 and κ 2 , by similar differential geometry formulas, which provide a system of ODE with respect to the variable t:

2z (t)r(t)z (0) = -λ r(t)r (t) + λ z (t) 2 + Boλ 2 z(t)z (t)r(t), r (t) 2 + z (t) 2 = λ 2 .
(5)

The boundary conditions can be expressed as (r(0), z(0)) = (0, 0), (r(1), z(1)) = (1, η), z (0) = 0. Now, we can approximate the functions r(t) and z(t) by cubic spline functions on a subdivision of the interval (0, 1) and express the Boundary Value Problem as a polynomial system of equations, which can be solved by standard techniques. The unknowns of the system are the coefficients of the two splines r and z and the parameter λ .

This second numerical approach was implemented using the software Scilab and is quite efficient since, for a given Bond number, it converges in few iterations of the nonlinear solver. Obviously, the curves obtained by the two numerical methods agree very well. Numerical results obtained with the second approach using a regular subdivision of 50 cells are presented in Fig. 3: at left, for a Bond number equal to 0.0052, the detachment height is estimated to 15.47 whereas at right, for a Bond number equal to 0.047, it is estimated to 7.90.

This second approach is still amenable when we replace the maximization parameter η by the volume of the bubble, an a priori more relevant physical quantity. The details will be given in a forthcomming article.

Conclusion

Two numerical methods that solve the capillary equation are reported. The objective of the numerical treatments is to generate curves that model the interface of vapour formations growing from a nucleation site. The first method is an established approach and is used here to provide context and clarity to the solution procedure. The second method is a new refinement to the numerical solution of the Capillary equation. By implementing splines into the procedure, it reduces the computational expense and presents a stepping stone towards more accurate heat transfer predictions of single boiling events.

Fig. 1

 1 Fig. 1 Axisymmetric surface with smooth contour.

Fig. 2

 2 Fig. 2 Set of curves satisfying (3) with the boundary conditions (4) for fixed values of η

Fig. 3

 3 Fig. 3 Numerical solutions obtained for Bo = 0.0052 (left) and Bo = 0.047 (right)

and

  Mass Transfer, 64:53-69, 2013. 2. Frédéric J Lesage, James S Cotton, and Anthony J Robinson. Analysis of quasi-static vapour bubble shape during growth and departure. Physics of Fluids, 25(6):067103, 2013. 3. D. Gerlach, G. Biswas, F. Durst, and V. Kolobaric. Quasi-static bubble formation on submerged orifices. International Journal of Heat and Mass Transfer, 48(2):425 -438, 2005.

Acknowledgements

This work has been supported by the French government, through the UCA-JEDI Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-15-IDEX-01.


Chang-Sheng Wang

Thierry Roisnel

Pierre H Dixneuf

Jean-François Soulé
email: jean-francois.soule@univ-rennes1.fr

Access to 3-(2-Oxoalkyl)-azaspiro[4.5]trienones via Acid-Triggered Oxidative Cascade Reaction through Alkenyl Peroxide Radical Intermediate

Keywords: Radical reaction, Cascade, Heterocycles, Ketones, Peroxides

trienones bearing ketone side chains at the 3-position are prepared from N-alkylarylpropiolamides and ketones via oxidative 1,2difunctionalization of alkynes. The cascade sequence starts with the generation of alkenyl peroxide intermediates, which are obtained by addition of tertbutyl hydroperoxide to ketones in presence of a catalytic amount of a strong acid. Then, the ketone radical adds to alkynes, followed by spirocyclization and dearomatization process. This method represents a new example of difunctionalization of alkynes with simultaneous formation of two carbon-carbon single bonds and one carbon-oxygen double bond in one step.

Spirocyclic compounds are important motifs embedded in many natural products, bioactive compounds and functional molecules. [1] Spiro [4.5]trienones, which are also valuable intermediates in organic synthesis, are often prepared via multi-steps synthesis. [START_REF]For reviews on the synthesis of spiro[4,5]decatrienones and related spirocyclic compounds[END_REF] One of the most efficient strategy for the construction of multiple bonds in a single operation involves 1,2-difunctionalization of unsaturated molecules, [START_REF]For reviews on 1,2-difunctionalization of unsaturated molecules[END_REF] including photoredox radical process. [START_REF] Courant | For reviews on radical 1,2-difuntionalization of unsaturated compounds using photoredox systems[END_REF] Among them oxidative radicalmediated spirocyclization of N-arylpropiolamides has proven to be an attractive and straightforward method to access spiro [4.5]trienones in a single step (Scheme 1). For example, radical oxidation of thiols, [5] sulfonylhydrazides, [6] diaryl diselenide, [7] nitrite, [8] or silanes [9] gave formation of heteroatom-centered radicals, which added to alkynes followed by spirocyclization and dearomatization (Scheme 1). This procedure allowed the 1,2-difunctionalization of activated alkynes through the one-pot formation of C-heteroatom and C-C single bonds. However, there are only few reports focusing on the use of carboncentered radicals, albeit it offers a direct synthesis of polycyclic compounds with the construction of two C-C single bonds. As examples, Liu and Liang have successfully employed Langlois' reagent (CF3SO2Na) as trifluoromethyl radical precursor in the synthesis of 3-trifluoromethylated azaspiro [4.5]trienones. [10] Alkyl radicals -generated by hydrogen abstraction process of ether [11] or cycloalkane C(sp 3 )-H bonds [12] using Cu catalyst with tert-butyl hydroperoxide (tBuOOH)-have been also employed in those cascade reactions (Scheme 1).

Scheme 1. Previous Syntheses of Functionalized spiro [4.5]trienones via Oxidative Radical-Mediated Spirocyclization of N-Arylpropiolamides.

In 2014, Klussmann and co-workers reported acidmediated 1,2-difunctionalization of alkenes with ketones and tert-butyl hydroperoxide (tBuOOH) (Scheme 2a). [13] The reaction involves the generation of ketone radical via the decomposition of alkenyl peroxide A -generated by the addition of tBuOOH to

R 4 R 4 N R 3 O Ar R O N R 3 O Ar R Cat.
Additive CuCl H C(sp

Accepted Manuscript

ketone catalyzed by acid-through homolytic bond O-O cleavage. [START_REF] Klussmann | For a review on alkenyl peroxide chemistry[END_REF] Independently, Klussmann's and Xia's groups have employed ketone radicals in cascade reaction. Indeed, after addition of ketone radical to alkene unit of N-aryl acrylamides, a radical cyclization affords oxindoles in good yields (Scheme 2b). [15] Xia and coworkers also reported a similar procedure for the preparation of isoquinolinonediones (Scheme 2c). [16] To the best of our knowledge, alkenyl peroxide intermediates have never been employed in addition of ketone radicals into alkynes, yet. In our continuous efforts to employ carboncentered radicals in C-H bond functionalizations, [17] we decided to investigate the reactivity of ketones in presence of tBuOOH and activated alkynes such as N-arylpropiolamides in order to introduce at least two different functional groups into unsaturated molecules (Scheme 2d).

Scheme 2. Application of Alkenyl Peroxide Intermediates in Organic Synthesis.

We selected N-methyl-N,3-diphenylpropiolamide and pinacolone as model substrates (Table 1). In the presence of 10 mol% p-TsOH and tBuOOH (in decane solution) in acetonitrile at 60 ºC, 3-(3,3dimethyl-2-oxobutyl)-1-methyl-4-phenyl-1azaspiro [4.5]deca-3,6,9-triene-2,8-dione (1) is obtained in a poor yield (entry 1). When the reaction is carried out without acid catalyst or peroxide, no reactions occurred (entries 2 and 3). The presence of oxygen has a deleterious effect on this cascade reaction (entry 4). The use of a larger amount of ketone (10 equiv.) allowed to improve the yield in 1 to 23% (entry 5). A better yield of 52% in 1 is obtained using pinacolone as reactant and solvent (0.2 mol/L) (entry 6). Aqueous solution of tBuOOH and other oxidants such as (tBuO)2 or PhI(OAc)2 are not successful (entries 7-9). p-TsOH outperformed other catalyst such as H2SO4, AcOH and Cu(OTf) (entries 10-12). A lower catalyst loading resulted in a lower yield, whereas a higher loading of 20 mol% failed to improve the yield (entries 13 and 14). Finally, a prolonged reaction time allowed to improve the yield up to 76% after 48 h (entries 14-16). The structure of 1 has been confirmed by X-ray analysis. [18] Interestingly, previous method to access azaspiro [4.5]trienones bearing ketone side chains at the 3-position involves oxidative ipso-annulation of specific substrates (i.e. N-(pmethoxyaryl)propiolamides with α-carbonyl alkyl bromides). [19] Table 1 1 mL p-TsOH (10) -82 (76) a) Determined by GC-analysis using n-dodecane as internal standard and isolated yield was shown in parentheses. b) Reaction carried out without tBuOOH. c) Reaction carried out under air atmosphere. d) tBuOOH in water solution used instead of tBuOOH in decane. e) (tBuO)2 used instead

R 1 O R 2 R 3 R 1 O R 2 O R 3 OtBu R 1 R 2 O OtBu via A
a. Acid-Catalyzed Oxidative Radical Addition of Ketones to Alkenes (Klussmann) [13] p-TsOH (cat.)

tBuOOH CH 3 CN, 50 ºC R 1 O R 2 b. Acid-Catalyzed Oxidative Radical Cascade Addition of Ketones to
Alkenes for Synthesis of Oxindoles (Klussmann) [15a] and (Xia) [15b] p-TsOH or MsOH (cat.)

tBuOOH 50 or 60 ºC N R 3 N R 1 O O O R 3 R 2 R 4 R 4 R 1 O R 2
c. Acid-Catalyzed Oxidative Radical Cascade Addition of Ketones to Alkenes for Synthesis of Isoquinolinonediones (Xia) [16] MsOH (cat.)

tBuOOH 60 ºC N R 3 O R 4 O N R 1 O R 2 O O R 3 R 4 R 1 O R 2

Accepted Manuscript

of tBuOOH in decane. f) PhI(OAc)2 used instead of tBuOOH in decane. g) 24 h. h) 48 h.

With the best conditions in hand, we then studied the effect of the nitrogen substituent on this acidmediated oxidative radical cascade reaction (Scheme 3). From, non-substituted N,3-diphenylpropiolamide, no reaction occurred. The replacement of N-methyl by N-ethyl substituent has almost no effect as the azaspiro [4.5]trienones 3 is isolated in 79% yield. N,3-diphenylpropiolamide bearing a removal Nprotecting group, e. g., N-allyl or N-benzyl, also reacted with pinacolone to give the corresponding 3-(3,3-dimethyl-2-oxobutyl)-azaspiro [4.5]trienones 4 and 5 in 53% and 71% yield, respectively without the cleavage of the N-protecting group. Then, we investigated the reactivity of other ketones in this acid-mediated oxidative radical cascade reaction (Scheme 4). Acetone underwent cascade reaction with N-methyl-N,3-diphenylpropiolamide or N-ethyl-N,3-diphenylpropiolamide to afford the azaspiro [4.5]trienones 6 and 7 in good yields. Going from primary to secondary ketones, a decrease of reactivity is observed as the spirocyclization with diethyl ketone afforded the product 8 in only 49% yield. This lack of reactivity could be attributed to steric factors. Unsymmetrical 1-cyclopropylethan-1one showed a preference for radical formation at the primary carbon, giving the corresponding cyclized product 9 in 68% yield. It is important to note that even under these radical conditions the cyclopropyl group remains intact at the end of the reaction. This observation, suggested that the formation of carboncentered radical via the decomposition of alkenyl peroxide is faster at primary carbon than tertiary carbon. Other unsymmetrical ketones such as hexan-2-one or dodecan-2-one displayed the same trend of reactivity and selectively reacted at the primary carbon yielding the targeted compounds 10 and 11 in 76% and 63% yield, respectively. This spirocyclization cascade reaction is tolerant to ester group as from N-methyl-N,3-diphenylpropiolamide and methyl levulinate, a bio-sourced ketone, the cyclized product 12 is isolated in 52% yield. However, no reaction occurred when (substituted) acetophenones are employed as reactants owing such ketones have more pronounced radical cation characters (see supporting information for more details). [20] Scheme 4. Scope of the Ketones in Acid-Triggered Oxidative Radical Cascade Addition -Spirocyclization -Dearomatization Process.

To further explore the scope of this radical cascade spirocyclization reaction, various substituted Narylpropiolamides are investigated under the optimized reaction conditions with pinacolone as radical precursor (Scheme 5). Firstly, the impact of ortho-substituents on aniline unit of Narylpropiolamide is investigated. When the cascade reaction is performed with ortho-C-halo (halo = F, Cl, Br) substituted substrates, the azaspiro [4.5]trienones 13-15 are obtained in moderate to good yields without cleavage of the C-halo bonds allowing further transformations.

From ortho-phenyl substituted N-phenylpropiolamide the cascade reaction, 16 is obtained in 77% yield. The reaction is slightly sensitive to the steric effect. Indeed, from Narylpropiolamides, in which the aniline moiety bears two electron-donating groups (e. g., Et or iPr) at both ortho-positions, the corresponding cyclized products 17 and 18 are obtained in 42% and 54% yield, respectively. Bromo or even iodo atom at the metaposition of the aniline unit are also tolerated affording the desired products 19 and 20 in good yields. Notably, when the reaction is performed using a substrate with a methyl group attached to the triple bond, N-methyl-N-phenylbut-2-ynamide, the acid-

(0.2 mmol) (1 mL) R = H R = Et R = All R = Bn NR 79% 53% 71% N R O Ph R 1 O p-TsOH (10 mol%) tBuOOH (5 eq.) 60 ºC, 48 h N R O Ph R 1 O O (0.2 mmol) (1 mL) N R O Ph O O R = Me 6 81% R = Et 7 86% N O Ph O O 8 49% R 2 R 2 N O Ph O O 9 68% N O Ph nBu O O 10 76% N O Ph nDec O O 11 63% N O Ph O O 12 52% CO 2 Me
Accepted Manuscript mediated radical cascade spirocyclization with pinacolone failed to deliver the corresponding azaspiro [4.5] Interestingly, the reaction with N-arylpropiolamides containing an aniline unit substituted by a paramethoxy or a para-fluoro group afforded 1 in moderate yields (Scheme 6). Moreover, the reaction with substrates bearing a para-chloro or para-bromo substituent to the amide moiety also give 1, albeit in a poor yield (16-24%). Whereas, no reaction occurred with substrates bearing a para-methyl substituent.

We attributed these results to the addition of the tertbutylperoxy radical to the para-position of the phenyl group, forming a C-O bond and leading to cleavage of the C-OMe, or C-X (X = F, Cl, Br) bonds to yield MeOH, and HX, respectively, prior to final elimination of tBuOH to afford the desired products (Scheme 8). Several control experiments are carefully carried out to gain deep insight into the reaction pathway. Firstly, azaspiro- [4.5]trienone (1) is not observed when 2 equivalents of 2,2,6,6-tetramethyl-1-piperidinyloxy (TEMPO), a radical scavenger, is added into the reaction system. However, the ketone-TEMPO adduct 22 is detected by GC-MS analysis. This observation confirms that the cascade spirocyclization involves ketone radical species (Scheme 7a). From (d6)-acetone KIE values of 1.50-1.06 are determined from intermolecular competition reaction and from two parallel reactions. [21] These results indicate that the formation of carbon-centered radical via the decomposition of alkenyl peroxide is not the rate-determining step (Scheme 7b). Finally, in the presence of 1 equivalent of water-18 O, the azaspiro- [4.5]trienone (1) is obtained in only 45% yield an no 18 O-labeled product was detected by GC-MS analysis (Scheme 7c). This labeling experiment indicated that the oxidation of aniline mostly occurred via the action of tert-butylperoxy radical (tBuOO • ) rather than a water attack. [10] N O R 5

tBu O p-TsOH (10 mol%) tBuOOH (5 eq.) 60 ºC, 48 h On the basis of this study and previous reports, [START_REF] Courant | For reviews on radical 1,2-difuntionalization of unsaturated compounds using photoredox systems[END_REF][5][6][7][8][9][10][11] a plausible mechanism is proposed (Scheme 8). Firstly, in the presence of a strong acid (p-TsOH), alkenyl peroxide intermediate B, is formed after addition of tBuOOH to the ketone and dehydration. Then, the decomposition of alkenyl peroxide B through homolytic bond cleavage leads to ketone radical D and tBuO • , which latter reacts with tBuOOH to give t-BuOO • and tBuOH. [13a] A subsequent radical addition of C into the alkyne at the α-position of the C=O bond of N-arylpropiolamide generates vinyl radical E, which undergoes an intramolecular radical cyclization at the ipso-position to give F. F is then trapped by the tert-butylperoxy radical (tBuOO • ) to afford G. Finally, H undergoes deprotonation and the elimination of tBuOH to afford the desired azaspiro- [4.5]trienone derivative and tBuOH. In summary, we have developed metal-free conditions for the diflunctionalization of activated alkynes with ketones. The reaction enables radical cascade additionspirocyclizationdearomatization process initiated by the formation of ketone radicals formation from ketones and tert-butyl hydroperoxide under acidic conditions. This one-pot cascade reaction give a straightforward access to 3-(2oxoalkyl)azaspiro- [4.5]trienones, incorporating a spirocycle unit that is a common structural motif in many natural products and pharmaceuticals.

N O R 5 tBu O O (0.2 mmol) (1 mL) R 3 R 3 R 4 R 2 R 4 R 2 R 1

Experimental Section

General Procedure for Acid-Triggered Oxidative Radical Addition -Spirocyclization -Dearomatization Process: In an oven-dried 15 mL Schlenk tube, N, 3-diarylpropiolamide (0.2 mmol, 1.0 equiv.) and TBHP (1 mmol, 5.0 equiv.) (TBHP: 5.5M in the dodecane solution) were dissolved in degassed ketone (1 mL). Then, 4-methylbenzenesulfonic acid (0.02 mmol, 0.1 equiv.) was added under a stream of argon, and the reaction mixture was stirred at 60 ºC over 48 h. After evaporation of solvent, the product was purified using flash column chromatography on silica gel to afford desired azaspiro [4.5]

  d. Acid-Triggered Oxidative Radical Cascade Addition of Ketones toAlkynes for Synthesis of Azaspiro[4.5]trienones (This work) structure of 1x (equiv or mL)

Scheme 3 .

 3 Scheme 3. Effect of the N-Substituent of Narylpropiolamides in Acid-Triggered Oxidative Cascade Addition -Spirocyclization -Dearomatization Process

Scheme 6 .

 6 Scheme 6. Reactivity of N-Arylpropiolamides Bearing a para-Substituent on the Aniline Unit in Acid-Triggered Oxidative Radical Cascade Addition -Spirocyclization -Dearomatization Proces.

  Scheme 7. Control Experiments

Scheme 8 .

 8 Scheme 8. Proposed Mechanism for Acid-Triggered Oxidative Radical Addition -Spirocyclization -Dearomatization Process

  trienone 21. However, substituent on the aryl group attached to the triple bond are welltolerated. As example p-hexyloxy substituted aryl alkyne delivered 22 in 62% yield, and o-phenylsubstituted aryl alkynes gave the corresponding product 23 in 57% yield. Naphthalen-1-yl alkyne was also viable to furnish 24.

	Scheme 5. Scope of the N-arylpropiolamides in Acid-
	Triggered Oxidative Radical Addition -Spirocyclization -
	Dearomatization Process

Acknowledgements

We acknowledge the China Scholarship Council (CSC) for a grant to CSW.

COMMUNICATION


André Galligo
email: galligo@unice.fr

Jean Rajchenbach
email: jean.rajchenbach@unice.fr

Bernard Rousselet
email: bernard.rousselet@unice.fr

S Barland

P Reynaud-Bouret

F Cauneau

U Guillouzouic

T Kuhl

F Passot

Planchon

Compressions of a polycarbonate honeycomb

Introduction

Honeycombs are widely used material and it is important to understand the mechanism governing their responses to compressions. Here we consider a polycarbonate honeycomb with circular close-packed cells, (as shown in the following pictures) that we submitted to a in-plane quasi-static uniaxial compression, with different orientations. Since the common tangents to the circles of the plane initial configuration form an hexagonal mesh, there are two different natural directions of compression, namely when the vertical compression axis is either parallel or perpendicular to a tangent. In the first case, the circles are positioned in columns, in the second case they are in a staggered arrangement. We observed that, as expected, in the first case the deformation was homogeneous, compressing all the rows. However, in the second case we observed a localization along few horizontal rows. More surprisingly the second kind of deformations were also reversible, contradicting usual expectations.

To understand these phenomena, we first compared our observations with other published experiments, mainly the excellent work of S. Papka et S. Kyriakides [START_REF] Papka | In-plane crushing of a polycarbonate honeycomb[END_REF], which is accompanied with coherent simulations performed by a detailed computer model. While their aim was to describing the different steps of a complete crushing of a material, the first steps of their described experiments coincide with ours. So our first contribution is the observation of what happens when one releases the load after the localization is achieved, and the surprising fact is that the process is reversible.

Our second contribution is the description of a geometric model for the deformation of the circular close-packed cells and the corresponding hexagonal mesh, controlled by the evolution of the total bending energy of the mechanism. For that purpose, we propose an alternative decomposition of the material in an aggregate of curved triangles instead of the obvious aggregate of circles. We also develop spline approximate models of their geometry. All of this allows us to provide a rational explanation of the considered phenomena.

The last section will be dedicated to the presentation and the mathematical analysis of a simple mechanical model with two degrees of freedom. The target will be to somehow illustrate, with a much simpler model, the observed localization of the deformations which combines a kind of buckling with a rotation; in particular the process is reversible.

We thank the CNRS Fédération Doeblin, for its support during the preparation of this work.

Experiments

We submitted a circular polycarbonate honeycomb to a in-plane quasi-static uniaxial compression, along two compression directions either parallel or perpendicular to a cell tangent. We observed two completely different behaviors: from the one hand a homogeneous deformation and for the other hand a localization.

Figure 1 shows the first and last step of a compression (before a release of the load); we observe that the middle rows of the cells are strongly bent on the left or on the right and almost crushed. In other words the deformation localizes along the middle rows.

We also notice a quasi invariance by some horizontal translations. The very surprising fact is that when we stop the compression at this point and leave the material free to relax, it returns to its initial configuration.

In an important paper [START_REF] Papka | In-plane crushing of a polycarbonate honeycomb[END_REF], S. Papka et S. Kyriakides studied the crushing of the same kind of polycarbonate honeycombs that we are considering.

We refer to the figure 2 of page 243 of the article [START_REF] Papka | In-plane crushing of a polycarbonate honeycomb[END_REF] of S. Papka et S. Kyriakides which describes the sequence they observed through 8 successive images of deformed configurations, corresponding to response to a long compression. The picture labeled (0) to (3) are completely conform to our observations, with the same kind of localization phenomena (before we release the load).

S. Papka et S. Kyriakides also mentioned accurate measures of the geometric characteristics of the cells of the honeycomb, pointing out variations of wall thickness and "ovalisation", illustrated by a picture that we sketch in the left panel of Figure 2. These facts were also confirmed in an article [START_REF] Hu | The inhomogeneous deformation of polycarbonate circular honeycombs under in-plane compression[END_REF] by L.L. Hu, T.X. Yu, Z.Y. Gao and X.Q. Huang, who also provided a convincing photo taken with a microscope.

We emphasize that, on the picture, we can notice that near the stick points, the ovals coincide at an higher order than an usual tangent point. We based our modelization on this important observation. Finally, we also report that they provided a computer model of the compression process. Their problem was discretized using the software ABAQUS with quadratic beam elements; the simulation of crushing of a honeycomb was conform to the observations summarized in the figure cited above and allowed to describe precisely the process.

Geometric Modelization

We consider that the basic polycarbonate material of the curved walls is isotropic and hyper elastic which implies that it admits reversible elementary deformations.

We will assume that the length of each arc of curve between two stick points does not vary during a deformation and that the deformations keep the sections in the same plane, so all the study will be conducted in 2D. An invariance by horizontal translations is also assumed.

We consider that the glue between the "circular" packed cells play a key role that we formalize through the two following geometric hypothesis: During all the compression and release process we consider that the curvatures at the sticking points are zero. At the beginning of the compression, the 3 semitangents at the stick points of a curved triangle meet in a same point. This property remains true for small deformations or when the triangle admits a symmetry axis. Assuming this simplifying property, we can associate to the polycarbonate an abstract hexagonal mesh (made of small tangent segments) which resembles an hexagonal honeycomb and follows the deformation of the polycarbonate.

It makes sense to view the polycarbonate honeycomb not as an aggregate of packed circles, but as an aggregate of curved triangles (represented by Bezier splines) that can be deformed, provided that any two adjacent semi-tangents coincide. We now present an explanation of the unexpected phenomena observed in our experiments.

Let us emphasize, with the illustration of Figure 2, that the initial curved triangles are positioned differently relatively to the load, so that after the first step of the compression the remaining axial symmetries are different.

Let us consider in both cases, (stacks in columns or staggered) the effects of a compression.

Stack "in columns'

• The compression respects the horizontal symmetry of each curved triangle, of each oval, and of each hexagon formed by the tangents. • The stack of hexagons made by the tangents deform (reversibly) to a flattened stack of hexagons.

"Staggered" stack

• After a certain load the "horizontal" side of the triangle reaches a maximum allowed length, then either it stops deforming or it buckles and a vertex bends (randomly on the left or on the right), then the triangle rotates in that direction. • This behavior localizes along few rows where imperfections create a weakest resistance to buckling. • In this process these curved triangles eventually rotate by an angle of π 6 , and behaves in a way similar to the one described for the curved triangles in a stack "in columns".

The whole phenomena is illustrated on the sequence of four images shown in Figure 2.0.2, their ordering 1 to 4 correspond to those of Figure 2. To emphasize the described process, we chose three different colors to single out the behavior of each of the curved triangles surrounding an oval in a localizing row. Note that a central symmetry is clearly conserved.

When the load is released, each flattened curved triangle reacts elastically, in particular the ones in the upper and lower rows which did not buckle or rotate recover their original form. In the localized rows, the vertices which were forced to meet are anew separated, can relax, the process is then inverted and the curved triangles in the localization row rotate in the opposite direction.

A simple model for buckling deformations

In this section we present a model with 2 degrees of freedom made of 3 bars linked with 2 hinges; one bar is linked to a fixed hinge whereas the third one is linked to an hinge which can move horizontally and on which a horizontal load is applied. The system may be described with 3 angles of the bars with respect to the unloaded configuration φ 1 , φ 2 , φ 3 ; see the figure below. The work of the applied load is Pw; where the horizontal displacements of the hinges are:

u 0 = 0, u 1 = l 1 (1 -cos(φ 1 )), u 2 = l 1 (1 -cos(φ 1 )) + l 2 (1 -cos(φ 2 )), w = l 1 (1 -cos(φ 1 )) + l 2 (1 -cos(φ 2 )) + l 3 (1 -cos(φ 3 )).
The horizontal ones are: x 0 = 0, x 1 = l 1 sin(φ 1 ), x 2 = l 2 sin(φ 2 ), x 3 = 0 and we have:

x 2 -x 1 = l 3 sin(φ 3 ); φ 1 = arcsin( x 1 l 1 ), φ 2 = arcsin( x 2 l 2 ) and φ 3 = arcsin( x 2 -x 1
l 3 ). We assume that torsion springs are active on hinges A 1 and A 2 ; we denote with θ 1 , θ 2 the relative angle of rotation of each bar with respect to its neighbor;

θ 1 (x 1 , x 2 ) = φ 1 (x 1 , x 2 ) -φ 3 (x 1 , x 2 ), θ 2 (x 1 , x 2 ) = φ 2 (x 1 , x 2 ) + φ 3 (x 1 , x 2 ).
Assuming the torsion springs to be linear elastic, the strain energy is 1 2 (K 1 θ 2 1 + K 2 θ 2 2 ). Then the total energy involving the work of the applied load is

V (x 1 , x 2 ) = 1 2 [K 1 (θ 1 (x 1 , x 2 )) 2 + K 2 (θ 2 (x 1 , x 2 )) 2 ] -Pw(x 1 , x 2 ). (1)
With this energy, we get the equilibrium equation:

∂V ∂ x 1 = 0, ∂V ∂ x 2
= 0. To simplify the computations, we assume here that l 1 = l 2 = l 3 = 1 and that K1 = K2 = 1. Then by symmetry, the equilibrium can be reached only if x 2 = x 1 or x 2 = -x 1 .

For P < 1, then x = 0 is the only stable equilibrium. For 1 < P < 3, then x = 0 is an unstable equilibrium and two stable equilibriums appears for x 2 = x 1 and P = arcsin(x 1 )

y A 1 A 2 A 3 A 1 A 2 θ 1 x 1 θ 2 x 2 φ 1 φ 2 φ 3 A 0 A 3 P

Fig. 1

 1 Fig. 1 Compression of a circular honeycomb in a staggered stack

Fig. 2

 2 Fig. 2 Left panel: Sketches of the microscopic geometry near a stick point. Right panel: Non linear behavior of the response.

Fig. 3

 3 Fig. 3 Model of curved triangle, undeformed and deformed.

Fig. 4

 4 Fig. 4 Two distinct directions of conserved symmetry.

Fig. 5

 5 Fig. 5 Deformations and rotations around an oval.

x 1 i.e |x 1 |

 11 almost equal to (6(P -1). For P > 3, two stable equilibriums appears for x 2 = -x 1 and |x 1 | almost equal to (2(P -3)/3. It corresponds to the configuration shown in the next figure.

  x