Alexis Burr

Wendy Noël

Pierrick Trecourt

Mathieu Bourcier

Fabien Gillet-Chaulet

Armelle Philip

Christophe L Martin

The anisotropic contact response of viscoplastic monocrystalline ice particles

Keywords: Ice, Contact mechanics, Anisotropy, Digital Image Correlation, Finite Element Method

Introduction

It has been suggested that dislocation creep is the main deformation mechanism responsible for the densication of rn at low strain-rates (<10 -5 sec -1 ) [START_REF] Salm | Mechanical properties of snow[END_REF][START_REF] Ebinuma | Experimental studies on densication and pressuresintering of ice[END_REF]. It is taken into account in densication models of rn for which den-Email address: christophe.martin@grenoble-inp.fr (Christophe L. Martin) Preprint submitted to Acta Materialia 24 avril 2017 sity is in the range 0.55 -0.84 gcm -3 (relative density in the range 0.60-0.92, ρ ice =0.917 gcm -3 ) [START_REF] Arzt | The inuence of an increasing particle coordination on the densication of spherical powders[END_REF][START_REF] Arnaud | Modelling of the densication of polar rn : characterization of the snow-rn transition[END_REF]. However, some other mechanisms have been also proposed for rn densication. These mechanisms are grain rearrangement, grainboundary sliding, sublimation, recrystallization, volume and surface diusion [START_REF] Herron | Firn densication : an empirical model[END_REF][START_REF] Alley | Firn densication by grain-boundary sliding : a rst model[END_REF][START_REF] Theile | Simulating creep of snow based on microstructure and the anisotropic deformation of ice[END_REF].

The relative importance of these mechanisms is still an open question and depends on the densication stage considered. A common feature of these proposed mechanisms is that they do not take into account the anisotropic viscoplastic behaviour of an ice crystal . However, Rolland du Roscoat et al. [START_REF] Du Roscoat | Analysis of Snow Microstructure by Means of X-Ray Diraction Contrast Tomography[END_REF], using Diraction Contrast Tomography on in-situ compression experiments, have demonstrated its importance in the intragranular deformation of large-grained samples during densication. The hexagonal structure of ice makes dislocations slip mainly on basal planes. This leads to a strong viscoplastic anisotropy : for a given applied stress, a crystal of ice deforms approximately 10 3 times faster in the basal plane than in any other planes [START_REF] Duval | Rate-Controlling Processes in the Creep of Polycrystalline Ice[END_REF][START_REF] Petrenko | Physics of Ice[END_REF].

Additionally, in the low range of density, the granular nature of snow was recognized very early on by analytical models that represented snow as a collection of grains bonded by solid necks [START_REF] Brown | A volumetric constitutive law for snow based on a neck growth model[END_REF]. More recently, Theile et al. [START_REF] Theile | Simulating creep of snow based on microstructure and the anisotropic deformation of ice[END_REF] coupled tomography with Finite Element Modeling, meshing each grain of ice by a beam element. Ice anisotropy was taken into account by the use of orientationdependent material parameters in the ow law. The use of beams was justied by the density range (0.20 to 0.35 gcm -3 ) studied by the authors, which qualies the material as low-density snow. In this range of densities, the Discrete Element Method (DEM) is also a tool of choice since it allows the discrete nature of snow to be accounted for. Hagenmuller et al. [START_REF] Hagenmuller | Microstructure-based modeling of snow mechanics : A discrete element approach[END_REF], also proposed to couple tomography reconstructed volumes of snow with simulations. Representing ice particles as rigid clumps of spheres, they showed that the compression behavior of snow is mainly controlled by the actual relative density. Snow densication may also be tackled using homogenization techniques as in [START_REF] Wautier | Numerical homogenization of the viscoplastic behavior of snow based on X-ray tomography images[END_REF], where again X-ray tomography images were used as initial microstructures.

At the other extreme of density, for deep dense ice core, the coupling bet-ween anisotropy (fabrics) and the state of stress, and their strong impact on ice ow has received much attention (see for example [START_REF] Durand | Change in ice rheology during climate variations -implications for ice ow modelling and dating of the EPICA Dome C core[END_REF][START_REF] Martín | On the eects of anisotropic rheology on ice ow, internal structure, and the agedepth relationship at ice divides[END_REF][START_REF] Montagnat | Measurements and full-eld predictions of deformation heterogeneities in ice[END_REF]). In particular, the eect of dynamic recrystallization that accommodates deformation is well documented [START_REF] Montagnat | Rate controlling processes in the creep of polar ice, inuence of grain boundary migration associated with recrystallization[END_REF].

In between these two extremes, in the intermediate range of relative density (0.60-0.90), the eect of ice anisotropy on the densication of rn has received much less attention. In that density range, rn may be considered as a material that consists of crystallites of ice in contact with each other. Firn bears some important granular characteristics such as grain rearrangement, grain rotation and macroscopic densication under compressive stresses. It may be advantageously modeled as a set of discrete grains that interact with each other through their contacts. The modeling of such a granular packing requires a good understanding of the contact forces that arise upon densication. However, a reliable contact law is still in need to account for the mutual indentation of ice particles.

The main ingredients of such a contact law should incorporate viscoplasticity, ice anisotropy, and size eects. Ideally, the contact law should be simple enough to allow for its implementation into an analytical model or a DEM code. Thus, the improvement of current models for rn densication requires the contact mechanics between two grains of ice to be considered in more details, with special attention to anisotropy.

The aim of this paper is to provide a rm experimental basis for the derivation of a simple contact model for ice particles. The article is organized as follows. In section 2, we present a contact model based on indentation theory that incorporates explicitly the viscoplastic anisotropy of an ice crystal. In section 3, the diametrical loading of one or two monocrystalline ice cylinders provides the necessary data to understand the phenomenology of ice particle indentation with particular consideration for anisotropy. This section also presents a partial validation of the model proposed in section 2. Section 4 uses a 2D FEM simulation in which a transverse isotropic law is used to nalize the validation of the analytical model.

Contact model

In this paper, we report on experiments focusing on the indentation of monocrystalline ice cylinders and in particular on their anisotropic response. The c axis, which denes the cristallographic orientation of an ice crystallite, is perpendicular to the basal plane of the ice hexagonal structure. The elastic behavior of ice is itself anisotropic. The eective Young's modulus varies when it is evaluated about axes that are rotated away from the c axis direction. Still, the maximum value (along the c axis) is only 30% larger than the minimum [START_REF] Fletcher | The chemical physics of ice[END_REF]. At the rather high temperatures considered here, elasticity may be neglected. Indeed, ice monocrystals undergo plastic deformation when a shear stress as low as 0.02 MPa is imposed on the basal plane [START_REF] Schulson | Creep and Fracture of Ice[END_REF][START_REF] Chevy | Creep behaviour of ice single crystals loaded in torsion explained by dislocation cross-slip[END_REF]. Easy basal slip can occur in the three potential glide directions on the basal plane. Thus the analysis of experimental results uses a constitutive equation that accounts only for the viscoplastic anisotropy of ice crystals [START_REF] Meyssonnier | A model for the tangent viscous behaviour of anisotropic polar ice[END_REF] by modeling ice as a transversely isotropic medium. It has been shown that this simple model captures the essential features of deforming ice single crystals [START_REF] Mansuy | Localization of deformation in polycrystalline ice : Experiments and numerical simulations with a simple grain model[END_REF]. Further simplication is achieved in this model by assuming a Newtonian behaviour for the ice crystallite. The literature indicates power-law exponents less than 2.5 for single crystals [START_REF] Duval | Rate-Controlling Processes in the Creep of Polycrystalline Ice[END_REF][START_REF] Petrenko | Physics of Ice[END_REF][START_REF] Chevy | Creep behaviour of ice single crystals loaded in torsion explained by dislocation cross-slip[END_REF][START_REF] Higashi | Plastic Yielding in Ice Single Crystals[END_REF].

The viscoplatic response of the crystal submitted to a uniaxial stress test with axial stress σ and axial strain-rate ε writes [START_REF] Meyssonnier | A model for the tangent viscous behaviour of anisotropic polar ice[END_REF] :

σ = 2µ basal (4α -1) f (θ) -1 εn (1)
where µ basal is a material viscosity parameter for shear in the basal plane, θ is the angle between the c axis and the loading axis (Fig. 1), and f (θ) is a π 2 periodic function :

f (θ) = 3 cos 4 θ + 2α sin 4 θ + 2 (4α -1) -β β sin 2 θ cos 2 θ -cos 2 θ ( 2
)
where α is the ratio of the axial viscosity along the c-axis to the axial viscosity normal to the c-axis ; β is the ratio of the shear viscosity parallel to the basal plane to the shear viscosity in the basal plane (µ basal ). Axial viscosities are similar along the c-axis and the basal plane [START_REF] Mansuy | Localization of deformation in polycrystalline ice : Experiments and numerical simulations with a simple grain model[END_REF], thus α is approximated in the remaining of this paper to unity leading to a simple expression for f (θ) :

f (θ) = 2 1 + 3 1 -β β sin 2 θ cos 2 θ (3)
A value β = 0.01 was proposed in [START_REF] Mansuy | Localization of deformation in polycrystalline ice : Experiments and numerical simulations with a simple grain model[END_REF] to t, within the approximation of a Newtonian behaviour (n = 1), the experimental data obtained on a monocrystalline inclusion embedded in a ne-grained isotropic ice matrix. Whatever its initial crystallographic orientation, the shape of the initially circular inclusion, became ellipsoidal as axial strain increased. The β = 0.01 value allowed for a reasonable agreement between observed and computed evolutions of the deformed ellipse aspect ratio and of the crystallographic orientation.

Figure 1b shows that the minima of f (θ) are located at θ = 0 and θ = π/2 while the maximum arises for θ = π/4. In other words, the creep rate is maximum when the c-axis is oriented at θ = π/4 with the axial loading axis. For a given applied stress, Eq. [START_REF] Arzt | The inuence of an increasing particle coordination on the densication of spherical powders[END_REF] shows that the maximum creep rate is 3 + β 4β times the creep rate of the least favorable c orientation (θ = 0, π/2). For β = 0.01, this ratio is approximately 75, indicating the large anisotropy experienced by ice crystallites.

Indentation laws are generally written for spherical particles, taking advantage of spherical symmetries for the derivation. For viscoplastic isotropic materials, Storåkers et al. [START_REF] Storåkers | The viscoplastic compaction of composite powders[END_REF] have for example proposed analytical contact laws that give the indentation load on a spherical body as a function of the viscoplastic material properties. More importantly for the problem treated here, these authors have resolved the contact problem between two dissimilar viscoplastic bodies i and j with constitutive equations of the form :

σ = σ i εn (4)
where the material parameter σ i pertains to material i but the creep index n is identical for the two materials. Based on indentation theory and on von Mises isotropic ow theory, the indentation load, N v , writes for two spherical particles of radius R i and R j (R * = R i R j / (R i + R j )) [START_REF] Storåkers | The viscoplastic compaction of composite powders[END_REF][START_REF] Storåkers | Similarity analysis of inelastic contact[END_REF] :

N v = χ σ -1 n i + σ -1 n j -n (2R * ) (1-n 2 ) h (1-n 2 ) ḣn ( 5
)
where h is the indentation and ḣ is the normal rate of approach of the two particles and χ is a geometric factor with no unit. N v is positive in compression ḣ > 0 . For two particles i and j, which c axis make angles θ i and θ j with the contact normal, the formulation of Storåkers et al. applied to Eq. ( 1) leads to :

N v = 2χµ basal (4α -1) (f (θ i ) n + f (θ j ) n ) -n (2R * ) (1-n 2 ) h (1-n 2 ) ḣn (6)
Eq. ( 6) is the general contact law of two particles which anisotropic behaviour is given by Eq. ( 1). Now specializing the problem to a Newtonian material and for two particles of identical radius R and a small indentation (h R), the contact radius, a (Fig. 1), may be approximated by :

a 2 = Rh (7)
Introducing the area of contact, A = πa 2 , Eq. ( 6) rewrites (α = 1), :

N v = 6χµ basal (f (θ i ) + f (θ j )) -1 A π 1/2 ḣ ( 8
)
The aim of this work is to analyze the applicability of Eqs. ( 3) and ( 8) to the indentation problem of two ice particles having arbitrary c orientations. The general formulation of Storåkers et al. [START_REF] Storåkers | Similarity analysis of inelastic contact[END_REF] applies to an indenter of arbitrary prole. However, Eq. ( 6) has only been derived for axisymmetric contact conditions (typically spheres or at punches). From an experimental point of view, it is too complex a task to machine or form spherical monocrystalline ice particles.

Instead, for practical reasons, we opt for cylindrical geometry. The full derivation of Storåkers et al. formulation for the problem of two indenting cylinders is beyond the scope of this paper (see [START_REF] Johnson | Contact Mechanics[END_REF] for the elastic problem of two identical indenting cylinders). Still, we aim at evaluating the suitability of Eq. ( 8)

to the 2-dimensional problem of two cylinders indenting each other and forming a contact of radius a and area of contact A = 2aL with L the cylinder length.

The rationale for cylinders also rests on the pluses of a 2D conguration, which allows the propagation of basal slips to be visualized by polarized light (section 3.2), and strains to be quantied by Digital Image Correlation (section 3.3).

These observations would be much more dicult with spheres. Finite Element simulations using the open source code Elmer/Ice [27] complement the experimental investigation and allow for a direct comparison with Eq. ( 8). These simulations are carried out on a cylindrical geometry (section 4).

Experimental indentation of ice cylinders

Methods

Diametrical loading tests were carried out either on a single cylinder or on two cylinders (Fig. 2). Cylinders (diameter : 20 mm ± 0.05 mm, length : 8 mm ± 0.05 mm) were made of single crystals of ice, which were slowly grown in a crystallizer on top of a cooling device [START_REF] Baruchel | Synchrotron Xray imaging for crystal growth studies[END_REF]. The c axis is parallel to the planar faces of the cylinder (the out-of-plane scatter is ±2 • ) to ensure plane stress conditions. The out-of-plane strain (ε xx ) is negligible. Various θ angles were selected to investigate the eect of the crystal orientation on the viscoplastic contact behavior. The creep test set up is detailed in Mansuy et al. [START_REF] Mansuy | Modelling the ice single-crystal viscoplastic behaviour[END_REF][START_REF] Mansuy | Localization of deformation in polycrystalline ice[END_REF] and is sketched in gure 2. A constant load was applied on the diameter of the samples and the maximum applied load (109 N ± 1 N) was small enough to avoid fracture. Tests lasting typically ve hours were carried out in a cold room at -10

• C (±1 • C).
Ice samples are placed between two plates of glass to avoid sublimation and to ensure plane strain conditions. Horizontal black marks allow the samples to be correctly positioned during the tests. Crossed polarizers were used to follow the evolution of basal slips. Images were recorded and analyzed by ImageJ [START_REF] Schindelin | Fiji : an open source platform for biological image analysis[END_REF] to compute logarithmic strains and rotations. The error on strain was evaluated to be of the order of 0.3% for the single cylinder conguration, and less than 0.1% for the two-cylinders conguration. For some tests with two cylinders, both cylinders were truncated to mimic indented particles and to limit rotation of the c axis during loading. The initial conguration is characterized by the ratio a0 R , where 2a 0 is the initial length of the truncated contact (Fig. 2b). Some of the two-cylinders conguration experiments were analyzed by Digital Image Correlation (DIC) to quantify the strain eld. The surface of ice samples was slightly scratched with glass paper and a speckle (shoe polish) was applied to obtain a proper grey level and a sucient contrast for the DIC analysis. A correct light-level was obtained by using two LED lamps with the advantage of limiting brilliance, enabling sucient contrast and avoiding shadows. This method was already successfully used on polycrystalline ice by Grennerat et al. [START_REF] Grennerat | Experimental characterization of the intragranular strain eld in columnar ice during transient creep[END_REF] and recently by Chauve et al. [START_REF] Chauve | Strain eld evolution during dynamic recrystallization nucleation; A case study on ice[END_REF]. Strain elds are extracted using the CMV software [START_REF] Allais | Experimental characterization of the local strain eld in a heterogeneous elastoplastic material[END_REF] using a plane strain hypothesis. Details on strain calculations can be found in [START_REF] Bourcier | Multiscale experimental investigation of crystal plasticity and grain boundary sliding in synthetic halite using digital image correlation[END_REF]. Integration schemes number 2 or 3 were chosen [START_REF] Allais | Experimental characterization of the local strain eld in a heterogeneous elastoplastic material[END_REF][START_REF] Bourcier | Multiscale experimental investigation of crystal plasticity and grain boundary sliding in synthetic halite using digital image correlation[END_REF] to ensure a correct identication of local strains with good accuracy.

Biquintic subpixel optimization was performed to improve the nal correlation output.

Diametrical loading on a single cylinder

Polarised light allows for the observation and the identication of the deformation by basal slips. A color evolution near the contact is detected when deformation initiates in the sample (Fig. 3). The deformation is strongly heterogeneous. However, because test conditions are symmetric, heterogeneities are also symmetric with respect to the cylinder center. The strain is localized in slip bands along basal planes. Two perpendicular lines appear for each contact point in red, contrasting with the undeformed green color. Slip lines are propagating through the crystal perpendicularly to the c-axis direction (Fig. 3). The bottom contact in Images (A2) and (A3) could indicate an additional slip line parallel to the c-axis direction. However, a closer examination reveals that these patterns are in fact short slip lines that develop perpendicular to the c-axis. These lines fully develop in images (A4) and (A5).

Since the basal plane provides only two independent systems, additional degrees of freedom are necessary to accommodate arbitrary loading. In ice, non-basal planes are seldom activated, thus requiring other mechanisms [START_REF] Schulson | Creep and Fracture of Ice[END_REF]. For example, recent studies [START_REF] Chauve | Strain eld evolution during dynamic recrystallization nucleation; A case study on ice[END_REF][START_REF] Chauve | Investigation of nucleation processes during dynamic recrystallization of ice using cryo-EBSD[END_REF] emphasized that recrystallisation could occur even at low strain (less than 1 %) or low stress (0.5 MPa) at grain boundaries and triple junctions of polycrystalline ice. Moreover, despite being commonly observed in hexagonal metals, twinning is excluded in ice. Kink and bending bands on the other hand can appear and were observed by Mansuy et al. on polycrystalline ice deformed with our experimental set up [START_REF] Mansuy | Localization of deformation in polycrystalline ice[END_REF][START_REF] Lebensohn | Modeling viscoplastic behavior and heterogeneous intracrystalline deformation of columnar ice polycrystals[END_REF]. In our loading conditions, the cylinder has few constraints and can deform easily by shear on basal planes except for specic values of θ (0 • and 90 • ). This peculiar condition is well illustrated by sample C in Fig. 3. For all other θ values the cylinder can deform freely. Melting is another candidate for strain accommodation. It may have occurred as large stress concentration may develop below the contact points. The melting point temperature decreases by about 0.1 • C for every 1 MPa [START_REF] Petrenko | Physics of Ice[END_REF]. A 100 MPa stress would thus be necessary to trigger melting under contacts at -10 • C.

In any case, all these mechanisms (recrystallisation, kink and bending bands, and melting) would have left a clear signature that was not observed with polarised light. Note that bright spots at contact points on sample B are linked to an image contrast enhancement processing. Hence, dislocation motion on {0001} planes should be the main strain accommodation mechanism in our experiment. This should mainly occur by screw dislocations with Burgers vector a 3 < 11 20 > as these are the least energetic dislocation lines [START_REF] Schulson | Creep and Fracture of Ice[END_REF] We observed a clear evolution of the orientation of the c-axis during the test (Fig. 3). Whatever the θ value, the c-axis rotates towards the loading direction in agreement with the compression of a single crystal [START_REF] Honeycombe | The plastic deformation of metals[END_REF]. For example, crystal (A) has rotated by approximately 10 • between its initial state (A1) and its (A5) state . In accordance with crystal anisotropy, strain is primarily limited when the c-axis is perpendicular or parallel to the loading direction as shown by sample (C). Signicant rotation is observed for samples with loads of the order of 100 N as demonstrated by Fig. 3. For these samples, basal glide propagates from the contact point (CP) through the whole cylinder. This is due to the heterogeneous stress eld around CP and to the strong ice crystal anisotropy. Note that the two samples corresponding to θ = 0 • and θ = 90 • lead to ∆θ = 0 • even for very long times (Fig. 3) .

Fig. 5 sketches the proposed main deformation mechanisms encountered for small and large initial θ angles. Cylinders are sheared by two localized bands to form three dierent fragments. The central fragment is contained between the contact points and the rst basal planes that shear (Fig. 3). The three fragments dier in size depending on initial θ angle. The central fragment contains all the planes that are gliding during the deformation. The other two fragments are simply pushed away by the rigid planes. For small θ values, two small caps quickly form. This is illustrative of sample (A) in Fig. 3. In that case, the rotation initiates early (Fig. 4a) but saturates further on since the small cap hinders rotation. For large θ values, two large caps gradually form. This is illustrative of sample (B) in Fig. 3. The rotation takes more time to initiate since dislocations need to propagate on a larger distance along the large cap. Once initiated, the rotation is facilitated by the geometry of the two large caps. Thus, the c-axis rotates towards the loading direction. In that nal conguration, the fragments get stuck and there is no more driving force for rotation. Therefore, the rotation reaches its maximum and stays constant with time.

Fig. 4b shows the strain evolution in the loading direction z. It follows a pattern which is very similar to the pattern of Fig. 4a. It indicates that for a caxis oriented parallel (θ = 0 • ) or normal (θ = 90 • ) to the loading axis, no strain is observed although a non-negligible load (50 N) is applied. We observed that those samples exhibited less than 1 % of axial strain even after more than 17 hours. Similarly, those sample exhibited negligible rotation. Fig. 4a and b also demonstrate the eect of the applied load intensity with increasing rotation and strain for increasing loads for a given θ value.

Diametrical loading on two cylinders

Tests carried out with the two-cylinders conguration are in qualitative agreement with the observations gathered on the single cylinder conguration.

For this conguration, we were able to evaluate the full strain eld by DIC on the face of the two cylinders (Fig. 6). In particular, it is possible to obtain the shear strain in the basal plane (ε bc ). The θ values of each sample were selected using polarized light outside the set-up. Once installed in the set-up, these values were conrmed by the observation of the basal plane angles, which are assumed to correspond to the maximum shearing. This hypothesis is validated by the very small values of the diagonal components of the strain tensor in the crystallographic axes reference ( b, c), such as ε cc as illustrated in Fig. 6a. These low eigenvalues of strain obtained by DIC are associated to the resolved shear stress needed to activate non-basal glide, which is 60 times greater than for the basal plane activation for the same strain rate [START_REF] Duval | Rate-Controlling Processes in the Creep of Polycrystalline Ice[END_REF][START_REF] Petrenko | Physics of Ice[END_REF].

As already observed in the single-cylinder conguration, Figs. 6b-d show strong heterogeneities in shear strain. These gures illustrate that shearing in localized bands becomes very large with increasing strain. For example, Fig. 6d shows that the maximum value of the shear strain in the basal plane ε bc at 200 min is more than 40%, while the macroscopic strain ε zz is only 10%.

We observed that for small strains the cylinder with the most disorientated c-axis (j cylinder) consistently exhibits localized shear strain around its contact with the other cylinder (Fig. 6b). For larger strains, heterogeneities that start from loaded contact points or from the extremities of the neck rapidly overthrow any other heterogeneities (Figs. 6c-d).

For the top cylinder (i), which exhibits a small θ angle, we observed that a strong shear band always develops from the contact with the rigid plane. The maximum shear strain is consistently larger in the small θ angle cylinder as compared to the maximum shear strain in the large θ angle cylinder. The width of the shear band in cylinder i is clearly larger than the one of the bottom cylinder j. This is in agreement with experiments involving only one cylinder where we observed large peripheral fragments for large θ angles. This conrms that the size of the shear bands depends primarily on the c-axis orientation, whatever the conguration. Fig. 7 groups the rotation ∆θ of the c-axis angle against time for samples tested either with point contact between the two cylinders ( a0 R = 0) or with nite contact size ( a0 R = 0.2, 0.32, 0.35). The imposed load for all samples is 54 N. Fig. 7 only plots ∆θ for the cylinder that experiences maximum rotation but the same qualitative conclusions can be drawn for others. Although error bars on ∆θ are rather large (Fig. 4a), machining an initial nite neck size tends to decrease the ability of cylinders to rotate. This may reect a more realistic state for which multiple contacts between ice crystals exist in the rn. Such contacts would have enlarged over time through sintering mechanisms and viscoplastic deformation to form necks with nite size.

Neck size eect

When rotation is hindered by an initially large neck size, Eq. ( 8) should allow for a simple relation between strain and load for a given set of angles θ i and θ j .

In this section, we show that experimental data on truncated cylinders support such a relation. For small axial strain ε zz , the initial area of contact A 0 = 2a 0 L is assumed to remain constant and the ratio of actual to initial indentation is :

h h 0 = 1 + 2 R h 0 ε zz (9)
with ε zz positive in compression. Integrating between times 0 and t for a constant load test (creep), Eq. ( 8) leads to :

ε zz a 0 R 1/2 = (f (θ i ) + f (θ j )) π 2 1/2 12χµ basal R 3/2 L 1/2 N v t (10)
Fig. 8 plots the product ε zz a0 R 1/2 against time for three tests with a0 R = 0.2, 0.32, 0.35 for a given load (53 N) and a given set of angles θ i (θ i = 30 • , θ j = 65 • ). Fig. 8 shows that the three tests fall on a master curve for small strains (or times). For times smaller than 80 min, the master curve exhibits a linear portion. Using the slope of this linear portion, Eq. ( 10), which is only valid for small strains, allows for an estimation of the product χµ basal . Experimental data points lead to χµ basal ≈ 4 10 12 Pa s. Although the exact value of the geometric parameter χ is dicult to ascertain for two cylinders, we note that for two indenting Newtonian spheres, χ ≈ 2 [START_REF] Storåkers | The viscoplastic compaction of composite powders[END_REF]. The viscosity µ basal at -10 • C may be estimated from the literature for the stress range of interest here at the contact (0.1 -10 MPa). It is of the order of 0.5 to 5.10 12 Pa s [START_REF] Duval | Rate-Controlling Processes in the Creep of Polycrystalline Ice[END_REF][START_REF] Mansuy | Localization of deformation in polycrystalline ice : Experiments and numerical simulations with a simple grain model[END_REF], and thus in good accordance with the value of the product χµ basal found here. In the following an approximate basal viscosity µ basal = 1 10 12 Pa s is chosen.

Comparison between 2D FEM simulations and experiments

The aim of this section is to complement the experimental validation of Eq. neglecting inertial terms, with the full transverse isotropic law of an ice crystal as proposed by Meyssonnier and Philip [START_REF] Meyssonnier | A model for the tangent viscous behaviour of anisotropic polar ice[END_REF]. It simplies to Eq. ( 1) for a uniaxial stress state of compression. In accordance with the analytical model, the values of parameters α and β are set to 1. and 0.01, respectively. Gillet-Chaulet et al. [START_REF] Gillet-Chaulet | Castelnau, A user-friendly anisotropic ow law for ice-sheet modelling[END_REF] have studied the eect of those parameters on the ow response of ice.

Details on the model and on the numerical methods can be found in [START_REF] Gillet-Chaulet | Flow-induced anisotropy in polar ice and related ice-sheet ow modelling[END_REF][START_REF] Gagliardini | Capabilities and performance of Elmer/Ice, a new-generation ice sheet model[END_REF].

Here, we use the initial angle orientation of the two c vectors dening the contacting cylinders as input, together with the geometry of the truncated contact. Cylinders were meshed using triangular elements with approximate size of 0.5 mm. No specic mesh renement was used close to the contact zone.

A displacement is imposed on the upper boundary to obtain a compressive force (54 N) with no lateral constraint. Note that neither the cylinder geometry nor the crystallographic orientation evolve for a given simulation. Fig. 9 shows the comparison between the FEM simulated and the DIC computed strain-rate elds. Since c vectors are given as input to the Elmer/ice code, the orientation of the shear bands coincides with the experimental one, as it should. More interestingly, FEM simulations predict satisfactorily the position and the width of the shear bands as compared to the DIC results. In particular, the FEM simulations reproduce the signicant eect of the contact on the spreading of the shear band. For example, considering the bottom cylinder, the shear band starts from point A and develops leftward until it encounters contact point B. This shear band is relatively thick due to the combination of geometrical and orientation constraints. In contrast, the shear band on the top cylinder is much thinner due to the c axis being nearly perpendicular to the load axis. Note that for the two bands, the strain-rate value becomes null at the extremities of the contact. Since DIC strain rates are calculated using incremental calculation of a strain map over time, they are far less resolved than strain maps which are calculated from the initial conguration (Fig. 6). Thus, the comparison between FEM simulations and experimental strain rate is merely qualitative. However, the strain-rate magnitudes of FEM and DIC are in reasonable agreement for the chosen value µ basal = 1 10 12 Pa s .

A second set of FEM simulations was carried out to conrm the validity of the analytical formulation proposed in Eq. [START_REF] Du Roscoat | Analysis of Snow Microstructure by Means of X-Ray Diraction Contrast Tomography[END_REF]. FEM calculations were performed for two truncated cylinders with various pairs of (θ i , θ j ) angles for four values of a0 R (0.141, 0.2, 0.3, 0.4). Fig. 10 plots the adimensional rate of approach of the two cylinders obtained by FEM simulations :

h = A 1/2 0 µ basal N v ḣ (11)
Using Eq. ( 8), the analytical expression of h rewrites :

h = f (θ i ) + f (θ j ) 6χ π 1/2 (12)
Fig. 10 compares the analytical expression to the FEM simulations. In this gure, the angle θ j of the vector c that characterizes the bottom cylinder is set to 0 and π/4 while the angle θ i is varied between 0 and π/2 by π/20 steps, taking advantage of the π/2 periodicity of the f (θ) function. Note that χ was used as a tting parameter to minimize by a least square method the error between Eqs. The match between Eqs. ( 11) and ( 12) is far from perfect with maximum error for the smallest a 0 /R value and θ j = π/4 of the order of 20% and negligible errors for θ j = 0. Still the comparison is satisfactory in that the general shape of the f (θ) function is well retrieved.

Concluding remarks

We have shown that the viscoplastic contact response between monocrystalline ice particles exhibits strong anisotropy. This anisotropy directly originates from the hexagonal structure of ice. When rotation is not hindered, the phenomenology of the indentation is rather complex, involving the propagation of basal slips and the formation of ice 'fragments' within the particle during loading. Localization of shear bands is the major strain mechanism. Light polarization and DIC are well suited to capture these strong heterogeneities that arise mainly parallel to the basal planes. These strain heterogeneities are coupled with a rotation of the c axis towards the loading axis.

When rotation is limited by a large contact area, the indentation evolution of two particles with arbitrary c axis orientation is more tractable. In that case, the simple analytical model presented here provides a means of extrapolating the contact response of strongly anisotropic ice particles. The simple forms of Eq.

(3), and Eq. ( 8) have been partially validated by a combination of experiments on truncated cylinders and by 2D FEM calculations. In particular, the role of contact size and of the crystallographic disorientation between the particles are correctly rendered.

It should be clear that the Newtonian law used here for single ice crystals is only a simplifying assumption which may not reect the correct strain-rate dependency. Power-law exponent of 1.5 to 2.5 are typically mentioned in the literature [START_REF] Duval | Rate-Controlling Processes in the Creep of Polycrystalline Ice[END_REF][START_REF] Chevy | Creep behaviour of ice single crystals loaded in torsion explained by dislocation cross-slip[END_REF]. For polycrystalline ice, both Newtonian and non-Newtonian behaviors have been observed [START_REF] Doake | Flow law for ice in polar ice sheets[END_REF][START_REF] Gillet-Chaulet | In-situ quantication of ice rheology and direct measurement of the Raymond Eect at Summit, Greenland using a phase-sensitive radar[END_REF]. In any case, it would be interesting to generalize this study to non-Newtonian behavior and investigate the eect of n > 1 in Eq. ( 5) on the shape of the h versus θ i curves (Fig. 10).

As stated above, the proposed analytical model is only valid when rotation may be neglected. Fortunately, this assumption is appropriate in many cases that are relevant for rn at intermediate relative density (above 0.60). In that density range, the rn still exhibits a particulate structure with ice particles contacting each other. In the densifying rn, a representative ice particle is in contact with an increasing number of neighboring particles. Assuming a spherical shape, the coordination number should be larger than ve for a relative density larger than 0.6. This should safely ensure a limited amount of rotation at the contact scale. Thus, we believe that the simple equations proposed here should be useful for applications where the viscoplastic anisotropy of the mu- Comparison of h given by nite elements simulations (Eq. ( 11)) and by the anisotropic contact model (Eq. ( 12)) for θ j = 0 (red dotted curves) and θ j = π/4 (blue plain curves). All curves are π/2 periodic.

Fig. 4a quantitatively conrms the qualitative observations of Fig. 3 .

 3 Fig.4aquantitatively conrms the qualitative observations of Fig.3. Fig.4aplots the rotation ∆θ of the c-axis angle against time for all tested samples in

( 8 )

 8 by specically investigating the eect of the θ angle on the normal force at the contact. The Finite Element Method (FEM) with the Elmer/ice opensource code was used to simulate some of the experiments with two truncated cylinders. The model implemented in Elmer/Ice solves the Stokes equations, i.e.

( 11 )

 11 and (12) (µ basal = 1 10 12 Pa s). A value χ ≈ 16 was found to minimize the error. Even though Fig.10only shows the comparison for θ j = 0 and θ j = π/4, the t has been conducted on a much larger set of θ j values (11 θ i × 11 θ j =121 pairs) for the four a 0 /R values studied here.

Figure 1 :Figure 2 :

 12 Figure 1: (a) Schematic representation of the anisotropic contact law between two indented spheres. ( b i , c i ) is the local reference of the ice crystal i. a 0 is the initial contact radius (Eq. 7) (b) Evolution with the angle θ of orientation function f given by Eq. (2) for α = 1 and β = 0.01.

Figure 4 :Figure 5 :Figure 6 :

 456 Figure 4: Evolution of the (a) rotation angle ∆θ and (b) strain εzz against time for monocrystalline cylinders of ice (single cylinder conguration). Samples are submitted to various loads and dierent initial angles of the c-axis with respect to the load direction are tested. Strain is positive in compression. Error bars are shown only for sample B for clarity. Black markers A, B, C relate to observations in gure 3. Sample C2 is outside the time scale (56 h). The two samples corresponding to θ = 0 • and θ = 90 • lead to ∆θ = 0 • . 25

Figure 7 :aFigure 8 : 29 2

 7829 Figure 7: Evolution of the disorientation of θ with the indendation for the cylinder that turns the most in two-cylinders conguration experiments (load 54 N).One can notice that maximum rotations for cylinders with initial neck size is lower than without neck. In other terms, the blue lines come from experiments in which cylinders rotate independently of one another in point contact congurations, while the red lines show that rotation is hindered because of the contact plane.

Figure 9 :

 9 Figure 9: Strain rate eld εzz (s -1 ) at 80 minutes under 54 N for θ initial (i) = 75 • and θ initial (j) = 63 • . a) Finite Element simulation results with µ basal = 1 10 12 Pa s and b) DIC results εzz.The gage length for this test is 864 µm ± 2 µm and the resolution is 21.60 µm/pixel ± 0.05 µm/pixel (integration scheme 2[START_REF] Bourcier | Multiscale experimental investigation of crystal plasticity and grain boundary sliding in synthetic halite using digital image correlation[END_REF]).

Figure 10 :

 10 Figure10: Evolution of the adimensional indentation rate h against θ i for four values of a 0 /R. Comparison of h given by nite elements simulations (Eq. (11)) and by the anisotropic contact model (Eq. (12)) for θ j = 0 (red dotted curves) and θ j = π/4 (blue plain curves). All curves are π/2 periodic.


Aaron Lichtner

Denis Roussel

Daniel Röhrens

David Jauffres

Julie Villanova

Christophe L Martin

Rajendra K Bordia

Anisotropic sintering behavior of freeze-cast ceramics by optical dilatometry and discrete-element simulations

Keywords: Freeze-Cast, Anisotropic Pores, Optical Dilatometer, Sintering, Discrete Element Simulations

come

Introduction

Freeze-casting is a materials processing technique that can be used to create a variety of anisotropic porous microstructures for use in fields such as materials for energy conversion devices [START_REF] Lichtner | Processing of Hierarchical and Anisotropic Porosity LSM-YSZ Composites[END_REF], dentistry [START_REF] Miller | Freeze Casting of High Strength Composites for Dental Applications[END_REF], bioengineering [START_REF] Deville | Freeze casting of hydroxyapatite scaffolds for bone tissue engineering[END_REF], and even food science [START_REF] Nguyen | Fast Dispersible Cocoa Tablets: A Case Study of Freeze-Casting Applied to Foods[END_REF]. This technique is attractive because of the high level of control and wide range of properties one can achieve with relatively simple equipment and primarily "green" methods. A significant number of researchers have reported on the effects of processing and microstructure on properties [START_REF] Deville | A meta-analysis of the mechanical properties of ice-templated ceramics and metals[END_REF][START_REF] Deville | Freezing Colloids: Observations, Principles, Control, and Use: Applications in Materials Science[END_REF] but little has been published concerning the actual sintering process of these highly anisotropic, porous structures. It is nonetheless an important topic of fundamental and practical interest. Freeze casting provides an interesting case study of the impact of a strongly anisotropic initial green microstructure on the sintering behavior. The shrinkage during sintering is an important fundamental outcome of sintering theories and models. It is well known that the shrinkage is a strong function of the microstructure. Thus, it is expected that these highly anisotropic microstructures will have anisotropic shrinkage. This study provides a quantitative confirmation of this and provides particle scale insights on the causes of this anisotropic shrinkage. From a practical point of view, the manufacturing of sintered parts with designed geometry requires in-depth understanding and modeling of the shrinkage during sintering and its dependence on the green microstructure.

Farhangdoust et al. [START_REF] Farhangdoust | The effect of processing parameters and solid concentration on the mechanical and microstructural properties of freeze-casted macroporous hydroxyapatite scaffolds[END_REF] reported strains in both directions after sintering hydroxyapatite and found that their samples shrunk less in the freezing direction than in their transverse direction. They argue that there is more empty space, spatially speaking, in the transverse direction so sintering will occur to a higher extent. However, it is generally admitted and observed that shrinkage is more pronounced in the direction of pore elongation (i.e. the radial direction for a uniaxially pressed sample) [START_REF] Olevsky | Theory of sintering: from discrete to continuum[END_REF]. Based on this commonly held belief, a freeze-cast sample would shrink more in the freezing direction, in apparent contradiction with the observation reported by Farhangdoust et al..

In this paper, we report on the in-situ sintering of freeze-cast composite porous ceramics with aligned macropores for use as solid oxide fuel cell (SOFC) cathodes. The fabrication process, described in [START_REF] Lichtner | Processing of Hierarchical and Anisotropic Porosity LSM-YSZ Composites[END_REF], consists of directionally freezing a well-dispersed, aqueous ceramic slurry. Due to the directional temperature gradient and anisotropic crystallographic nature of ice, long continuous ice channels form within the suspension, consolidating the ceramic particles into the interlamellar spaces between ice crystals. Subsequent sublimation of the ice leads to a green body with macropores where the ice crystals had been and micropores between the ceramic particles of the green body.

Sintering is then typically performed to consolidate the green body enough to provide adequate strength while maintaining a high degree of reactive surface area [START_REF] Lichtner | Dispersion, connectivity and tortuosity of hierarchical porosity composite SOFC cathodes prepared by freeze-casting[END_REF][START_REF] Lichtner | Effect of Macropore Anisotropy on the Mechanical Response of Hierarchically Porous Ceramics[END_REF].

To investigate the effects of the anisotropic microstructure on sintering behavior, an opti-cal dilatometer (TOMMIplus, Fraunhofer ISC, Wertheim, Germany) was used. This is a hightemperature furnace equipped with two aligned quartz windows on either side of the sintering sample (Fig. 1). At one quartz window is a high-powered light source and at the other is a highspeed camera with a telecentric lens. The light source creates a silhouette of the sample that is captured by the camera. Software is employed to track the edges of the silhouetted sample as it expands or contracts over time [START_REF] Raether | Current state of in situ measuring methods for the control of firing processes[END_REF]. Since an entire face of the sintering sample can be observed, the strains in multiple directions can be tracked simultaneously making this technique ideal for anisotropic or irregularly shaped samples compared to conventional contact dilatometry.

3D Discrete Element Method (DEM) simulations were also performed to confirm the sintering behavior observed experimentally and support the discussion. Real 3D freeze-cast structures, obtained from high resolution synchrotron-based X-Ray tomography, were used [START_REF] Villanova | Multi-Scale 3D imaging of absorbing porous materials for solid oxide fuel cells[END_REF]. Freeze-cast samples were placed so that the freezing direction (longitudinal) was pointing up.

Methods

Optical Dilatometry

Freeze-cast specimens with 80 vol. % green-state total porosity were made as per [START_REF] Lichtner | Dispersion, connectivity and tortuosity of hierarchical porosity composite SOFC cathodes prepared by freeze-casting[END_REF][START_REF] Lichtner | Effect of Macropore Anisotropy on the Mechanical Response of Hierarchically Porous Ceramics[END_REF] [START_REF] Lichtner | Effect of Macropore Anisotropy on the Mechanical Response of Hierarchically Porous Ceramics[END_REF], ensured that only the steady-state anisotropic regions of the freeze-casts were used for dilatometry experiments. After sectioning, the specimens were measured precisely then heated to 500 • C at 2 • C min -1 and held there for 2 hours until the resin had completely burned out.

The fragile samples, now devoid of any organic residue, were placed in the optical dilatometer with their freezing direction pointing up. During testing, the samples were sintered to 1400 • C at a rate of 10 • C min -1 in air. The optical dilatometer camera captured an image each minute by averaging a total of 10 images to reduce noise.

Data Analysis

For each sample tested, the series of images produced by the optical dilatometer were thresholded into binary images. Then, using a plugin written for ImageJ, black pixels in both the vertical and horizontal directions were counted at predetermined locations giving both the longitudinal and transverse strains as a function of time and temperature. On average, both longitudinal and transverse strains were measured at three distinct locations for each sample. The curves were then fitted with a gaussian smoothing algorithm to reduce noise.

Discrete-Element Simulations

DEM simulations are based on the sintering model used in [START_REF] Martin | Discrete element modeling of metallic powder sintering[END_REF][START_REF] Martin | Evolution of Defects During Sintering: Discrete Element Simulations[END_REF]. This model originates from where grain boundary and surface diffusion are the major mechanisms of mass transport. In DEM, particles rearrange due to the action of sintering forces at contacts with neighbouring particles.

For this study it should be made clear that we are interested in the evolution of the freeze-cast microstructures with and without macropores, not in the kinetics of sintering.

Particles are modeled as spheres interacting dynamically with each other through their contacts.

The diffusion parameter ∆ b is modeled as:

∆ b = Ω kT δ b D b ( 1
)
where Ω is the atomic volume and

D b = D 0b exp(-Q b /RT ) is the diffusion coefficient for vacancy
transport in a grain boundary with thickness δ b and activation energy Q b . For two spherical particles of identical radius R, having a contact radius a s and particle overlap h, the normal force, N s between the particles is given by:

N s = πa 4 s 8∆ b dh dt - 9 8 πRγ s (2)
where the surface energy between the particles is denoted as γ s . Tangential contact forces are neglected in these simulations since we assume that contacts are large enough to oppose any rotation of particles.

The contact radius simply grows according to Coble's geometric model [START_REF] Coble | Initial sintering of alumina and hematite[END_REF]:

da s dt = R a s dh dt (3)
For a multimodal packing as modeled here, we define R * as the equivalent radius for two particles of radii R 1 and R 2 :

R * = R 1 R 2 /(R 1 + R 2 ), (4)
As suggested by Parhami et al. [START_REF] Parhami | A model for the sintering and coarsening of rows of spherical particles[END_REF], the model of Bouvard and McMeeking may be generalized for contacts between particles of different sizes by replacing R in equations ( 2) and ( 3) by 2R * .

For this study, two numerical structures were generated, one with realistic freeze-cast macropores and one without. To do so, two identical packings of particles were generated with 50 vol. % total porosity, typical for green bodies without macroporosity. A bimodal size distribution corresponding to the 40:60 vol% LSM:YSZ ratio was used with particle sizes of 0.5µm (YSZ) and 1µm (LSM).

A 3D freeze-cast structure was reconstructed from a synchrotron X-ray characterization of one of our samples [START_REF] Villanova | Multi-Scale 3D imaging of absorbing porous materials for solid oxide fuel cells[END_REF]. The voxel resolution for this image was approximately 75 nm, while the field of view was ∼ 90 3 µm 3 . The image was merged with one of these DEM particle packings to produce macropores in addition to the already present interparticle pores. The DEM structure from the synchrotron reconstructed image contained approximately 208,000 YSZ particles and 17,000 LSM particles. The average initial coordination number was 3.1 in the overall structure and 3.8 inside the particle-packed regions. The procedure to generate the DEM structure from the synchrotron reconstructed image has been described in more detail in [START_REF] Roussel | Strength of hierarchically porous ceramics : discrete simulations on X-ray nanotomography images[END_REF]. Both DEM structures, with and without macropores, were then sintered numerically until a density of 0.75 was reached within the particle-packed regions, hereafter called the walls. The intrawall density of the simulated structure was equivalent to that of a physical LSM-YSZ freeze-cast sintered to 1400 • C.

It is important to note that since the size of the simulated structures is so much smaller than the physical ones (typical experimental samples are tens of mm in size, numerical sample are only 90 µm per side), there are some important differences in their overall structures. Large scale freeze-casts are composed of colonies of lamellar pores which are randomly oriented in the x-y plane [START_REF] Lichtner | Effect of Macropore Anisotropy on the Mechanical Response of Hierarchically Porous Ceramics[END_REF]. The simulated freeze-casts contain only those pores within a single colony, consequently, they lose the isotropic nature of the transverse direction seen in the experimental samples. The result of this is that physical samples have two distinct microstructural and behavioral directions (longitudinal and transverse) while the DEM samples have three (longitudinal, transverse-align, and transverse-non-align). The behavior of the transverse direction in physical freeze-casts will be some linear combination of the transverse-align and transverse-non-align direction observed in the DEM samples.

Results

Experimental Sintering Behavior

Sintering studies of 40:60 vol. % LSM-YSZ show that densification starts to occur at approximately 1000 • C and continues until full densification is achieved at approximately 1500 • C. The vast majority of densification occurs in the 1100 -1300 • C range. The optical dilatometer used in this study was limited to a 1400 • C maximum, which despite not capturing the full sintering profile, does capture the majority of densification. Figure 3 shows the aggregated sintering profiles from multiple freeze-casts. One can see that the longitudinal and transverse strains begin to deviate at approximately 1200 • C, with the longitudinal direction shrinking to a greater extent, on average, than the transverse one.

Physical measurements of the porous freeze-casts before and after sintering to 1400 • C clearly show the anisotropic sintering behavior experienced by the freeze-casts in the two testing directions (Fig. 4). In the longitudinal direction, on average there was an additional shrinkage of 3.4 % compared to the tranverse direction. This is contradictory to the results of Farhangdoust et al. [START_REF] Farhangdoust | The effect of processing parameters and solid concentration on the mechanical and microstructural properties of freeze-casted macroporous hydroxyapatite scaffolds[END_REF] who claim that sintering strains are enhanced in the transverse direction, not the longitudinal one. In our experiments, cuboid samples were cut so as to exclude the initial, isotropic cellular structure found at the bottom that may have influenced the sintering behavior seen by Farhangdoust et al..

In addition, by doing so, the continuous skin that generally forms on the side of a freeze-casted cylinder is also removed. The presence of such a continuous skin may have promoted the transverse shrinkage seen by Farhangdoust [START_REF] Farhangdoust | The effect of processing parameters and solid concentration on the mechanical and microstructural properties of freeze-casted macroporous hydroxyapatite scaffolds[END_REF].

Simulated Sintering Behavior

DEM simulations were performed on particle packings with and without macroporous structures.

For simulations performed on structures with realistic freeze-cast macropores, plotting the strains in the x, y and z-directions one can see that there is one direction that shrinks to a lesser extent (∼ 3 %) than the other two (Fig. 5). As expected, this direction , y, is normal to the freezing-direction.

The microstructure along this direction (Transverse-non-align) is most similar to the transverse direction in physical samples which also saw lower strain values during sintering (Fig. 5).

Note, that the x-direction (transverse-align) sees equivalent strain values to the freezing direction (z), illustrating that within a colony of pores, the properties in the freezing-direction are similar to those in the transverse-align direction. The quantitative comparison between the macroscopically measured strain (Transverse: -0.155 ± 0.015, Longitudinal: -0.19 ± 0.01) and the DEM simulations (y : -0.1, x, z : -0.14) shows that DEM simulations underestimate strains. This may be due to the size of the DEM representative volume obtained from the X-ray tomography which may be too small to fully model the behavior seen in physical samples. Still, we note that the ∼ 90 3 µm 3 freeze-cast DEM structure captures approximately both the strain magnitude and the ratio between the transverse and the longitudinal strains.

As a comparison, the sintering strain experienced by the DEM packing without macropores is plotted alongside strain values from the sintered freeze-cast DEM structure (Fig. 6). We show here that the freeze-cast structure has an unequivocal effect on the sintering process.

X Y Z X Y Z X Y Z
with Freeze-Cast Macropores No Macropores Strain

Discussion

Freeze-cast LSM-YSZ underwent anisotropic densification during sintering. Experimentally, samples shrunk to a greater extent (3 -4 %) longitudinally (freezing direction) compared with their transverse directions when sintered to 1400 • C. Discrete element sintering simulations of equivalent particle packings confirm that the presence of macropores contained within the freezecast responsible for the anisotropic behavior (Fig. 6).

To complement the anisotropic, macroscopic DEM strains seen in Fig. 5 we took strain measurements in between pairs of specified particles of the DEM numerical structure. These investigations provide additional understanding of how sintering evolves at the local level within freeze-cast structures. Fig. 7 shows the evolution of these local strains as a function of relative density. We chose to measure the strains in the walls in the x, y, z directions and across two walls (i.e. the strain in a macropore) in the y direction. Several points were chosen to evaluate the dispersion accompanying these strain measurements. Fig. 7 shows that the local strains in the y direction, be it in the walls, wall(y) , or in the macropores, macropore(y) , are on average, smaller than the strains in the wall in the x or in the z directions ( wall(x) , wall(z) ). This is consistent with the macroscopic strains that have been measured at the scale of the whole freeze-cast DEM structure (Fig. 5). Strains in Fig. 7 are of course distributed since they characterize strains in between pairs of discrete points in the structure. Interestingly, the strains in the y direction are similar in the solid wall and in the macropore, demonstrating that the freeze-cast structure shrinks homogeneously in the transverse-non-align direction, albeit to a lesser extent than in the transverse-align direction. ε macropore(y)

ε wall(x)
ε wall(z) As another level of investigation, we can examine the coordination number of particles in a DEM structure after sintering. Areas which are more densely packed will have densified more substantially than others. Looking at results of the discrete element sintering simulations in figure 8, one can see that the regions of highly coordinated particles run along the x and z-directions. These are both the directions where greater overall shrinkage values were observed. Comparatively, the y-direction is highly discontinuous, with only small particulate bridges connecting walls. In addition these bridges have experienced lower densification than the walls as revealed by their lower coordination numbers.

Tanaka et al. [START_REF] Tanaka | Packing structure of particles in a green compact and its influence on sintering deformation[END_REF] demonstrated that the anisotropy of the contact orientation network and its influence on sintering may be quantitatively rendered by studying 2D images of die-pressed compacts. Alternatively, using the full 3D information from our DEM structure, the stereographic projections of the contact orientation network can be obtained. Fig. 9 a shows that indeed in the center of the projection, there is a lack of contacts whose normals are preferentially oriented in the y direction (as indicated by the yellow color in the center of the pole figure). This lack of contacts with normals oriented in the y direction is not immediately obvious upon observing the stereographic projection. Still, it is confirmed by Fig. 9b, which shows the histogram of contact angles with respect to the y axis. As the contact angle increases from θ = 0, a gradual increase in the number of occurrences is clearly indicated.

The anisotropy shown in Fig. 9 is not pronounced because it is only the particles that are found along free surface of the walls that lead to a macroscopic contact anisotropy (the contacts within the walls themselves are still isotropically distributed). Still this anisotropy in contact orientation is sufficient to lead to the observed anisotropic shrinkage. It also confirms the qualitative information given by the coordination number map in Fig. 8.

Sintering is a mass-transport phenomena, it is driven by the system's desire to reduce surface area and controlled by mass diffusion. Particles that are well within the walls of a freeze-cast are surrounded in all directions by particles and so will sinter isotropically to reduce that energy.

Particles which find themselves running along the free surfaces of macropores are only partially surrounded by other particles (Fig. 8). Because of this they will experience a lower sintering driving force. In other words, the lack of contacts whose normals are oriented along the y direction (Fig. 9), will result in enhanced shrinkage along the the walls (x and z).

Conclusions

Sintering of freeze-cast structures resulted in significantly different relative levels of shrinkage in the longitudinal and transverse directions. As shown by optical dilatometry measurements and DEM simulations, the longitudinal direction shrinks more than the tranverse one ( ∼ 3 -4 % more). Localized strain measurements between pairs of particles within DEM simulations indicate that not only do freeze-cast samples shrink anisotropically on the macro-scale but do so at the scale of the walls and macropores themselves. Stereographic projection of the particle contact orientation network from the unsintered sample suggests that the freeze-casting process reduces the number of contacts with specific orientations. The relatively lower ratio of particle contacts with normals aligned along the transverse direction results in less shrinkage normal to the freezing direction Reduced densification of the particulate bridges that connect walls, as compared to the walls themselves, is confirmed by DEM results. The orientation of free surfaces and the fact that, spatially there is more void space in the transverse direction (macropores) favors anisotropic shrinkage. This study demonstrates that near netshape sintering of porous anisotropic parts might not be straightforward even for simple geometries. Thus, the successful forming of porous anisotropic ceramics requires a detailed investigation and further modeling of sintering behavior.

Figure 1 :

 1 Figure1: Schematic of the inner-workings of an optical dilatometer. Freeze-cast samples were placed so that the freezing direction (longitudinal) was pointing up.

T 1 Figure 2 :

 12 Figure 2: Schematic illustrating how strains are determined from images in multiple directions. The figure shoes the sample at the initial temperature, T 0 (black) and the sample at the sintering temperature, T 1 (grey).

Figure 3 :

 3 Figure 3: Aggregated sintering profiles from samples oriented with their freezing-direction up (longitudinal). Stillframes from a freeze-cast sintering run are included for reference. A minimum of 15 measurements were done for each curve.

Figure 4 :

 4 Figure 4: Sintering shrinkage from averaged physical measurements taken before and after sintering freeze-casts for sintering temperatures of 1400 • C in air. A minimum of 5 samples were utilized per treatment.

Figure 5 :

 5 Figure 5: (a) Simulated sintering of a freeze-cast DEM structure (∼ 90 3 µm 3 ) showing (b) anisotropic shrinkage. Also shown are the (c) x-y and (d) z-y cross-sections to illustrate the macrostucture and relevant shrinkage. Macropores are shown in white.

Figure 6 :

 6 Figure 6: Strain values in each direction at the final density of 0.75 for DEM simulated structures with and without macropores.

Figure 7 :

 7 Figure 7: Evolution of strains in the walls ( wall ) and in macropores ( macropore) at selected pairs of particles used as strain gauges.

Figure 8 :

 8 Figure 8: Coordination analysis of particles in a sintered DEM structure.

Figure 9 :

 9 Figure 9: Anisotropy of contact orientation. a) Stereographic projections of the contact orientation network in the (x, z) plane from the unsintered DEM structure. Colors indicate relative density in arbitrary units. b) Histogram of contact angles with respect to the y axis. The number of occurrences is normalized by 2π (cos (θ) -cos (θ + δθ)) for each bin.

Acknowledgments

This work was supported by the U.S. National Science Foundation (NSF) under grant no. 1008600 and l'Agence de la Recherche under grant no. 2010 BLAM 0931 01. RKB also acknowledges partial support from the NSF under grant no. 1502392. Fig. 9a was built with the help of Alexis Burr using the AITA python toolbox (https://github.com/ThomasChauve/aita).