1. Introduction
Attributes of snow (3-D) microstructural geometry have a dominant influence on its mechanical behavior. Geometric attributes like amount, size and shape distribution and spatial arrangement of 3-D microstructural features such as pores (or grains, inclusions, precipitates, etc.) (Reference Adams, Miller and BrownAdams and others, 2001) have a definitive effect on failure properties of snow like shear strength, fracture toughness, fatigue and creep resistance. Therefore, it is of interest to incorporate quantitative and qualitative description of actual 3-D microstructure in the micromechanical analysis of snow and ice.
In the past, snow 3-D micromechanical calculations have been performed on idealized microstructures having very simple geometry (e.g. one spherical inclusion in an infinite matrix) using analytical techniques or a voxel-based finite-element (FE) model approach, where the digitized voxels are converted to finite elements, resulting in a very large number of equations (Reference GarbocziGarboczi, 1998). Obviously, such models ignore the complex geometry and the physics of real 3-D microstructures. For example, the effect of local spatial arrangements of microstructural features like grain surface shape on the micro stress or strain distribution is completely ignored or overestimated in such calculations. Another prime source of difficulty is that data on mechanical properties like Young’s modulus and the bulk deformation behavior of the material have usually been organized and presented as functions of the snow density (Reference MellorMellor, 1975). However, the literature shows that snow density alone is not a reliable indicator of these properties. Instead, for a given temperature and loading condition, the response to load depends primarily on the bonding, microstructure and the geometric characteristics of the grains (Reference Schweizer, Schneebeli, Fierz and FöhnSchweizer and others, 1995).
The real micro geometry of snow can be explored by new X-ray microtomography techniques (Reference Coléou, Lesaffre, Brzoska, Ludwig and BollerColéou and others, 2001). This makes it possible to use these data to numerically model the micromechanical behavior. We describe the first developments in the analysis of the microstructure volume data to model and simulate the micromechanics of snow using image analysis and volumetric surface extraction techniques.
Section 2 of this paper explains the material and methods used to obtain the snow 3-D data and the development of the code for a 3-D geometry model and FE analysis simulation adapted to these scales. Section 3 introduces a hypothetical compressive loading simulation on a cylindrical 3-D snow sample geometry, presenting the constitutive equations, assumptions and boundary conditions.
2. Material and Methods
2.1. Snow sample modeling: metamorphism experiment
A 3 month long experiment of isothermalmetamorphism at –2°C was run at Col de Porte, Chartreuse mountain, French Alps, to provide tomographic 3-D data for this work and for the development and validation of metamorphism models (Reference Flin, Brzoska, Lesaffre, Coléou and PieritzFlin and others, 2003).
A 0.5 m×1m × 0.12 mthick slab of recently fallen snow was collected in the field on 16 January 2002, 12 hours after the snowfall (air temperature –1°C). It was stored inside a cold room (-2 ± 0.2°C) for 3 months in a closed Styrodur (extruded rigid polystyrene foam) box to prevent sublimation. All manipulations were done in the cold room at this temperature. A 3 cm diameter core was sampled at increasing intervals, ranging from 24 hours at the beginning to 1 week at the end of the experiment. All samples were taken at mid-height of the slab, always 410 cm away from already sampled regions of the slab.
Once sampled, each core was impregnated with 1-chloronaphthalene (90% purity, practical melting range after raw fractioned crystallization –15 to –20°C), then allowed to freeze and stored in a refrigerator at 20°C. Owing to the stiffness of its cyclic molecule, this low-toxicity filler is readily machined close to its melting point; moreover, its X-ray absorption properties provide good contrast to distinguish ice, filler and remaining air bubbles on 9 mm wide samples at 18–20 keV. At the end of the metamorphism experiments, the cold room was set to –25°C; each strengthened core was machined into a cylinder 9 mm in diameter and 9mm in height, and sealed into a gas-tight sample holder made of 0.2 mm thick clear plastic. Secured samples were stored until the beginning of the tomography in a refrigerator at -50°C.
2.2. The microtomography X-ray data acquisition
The principle of X-ray absorption tomography is to reconstruct the spatial distribution of the linear attenuation coefficient within the object from X-ray projections recorded at different angular settings of the sample. Each element in the recorded projection corresponds to a line integral of the attenuation coefficient along the beam path. Via a standard reconstruction algorithm (filtered back-projection) (Reference HermanHerman, 1980) the 3-D distribution of the attenuation coefficient can be calculated by combining the information for the different angular positions. In our case, 1500 projection images spanning 180° rotation around the vertical axis were recorded for each sample. The reconstructions were performed offline using standard filtered back-projection software available at the European Synchrotron Radiation Facility (ESRF), Grenoble, France.
The experiments were carried out at the ID19 beamline of the ESRF. ID19 is devoted to high-resolution imaging, and features an energy-tunable, parallel X-ray beam of up to 14 mm × 40 mm cross-section (energy range 7–100 keV, divergence1 μrad). The cryostat (Reference Coléou, Lesaffre, Brzoska, Ludwig and BollerColéou and others, 2001) was mounted on a precision mechanics sample stage, and the projections (Fig. 1) were recorded with a fast, high-resolution detector system. The latter consists of a powder phosphor screen which is coupled via light optics to a cooled 20482 charge-coupled device (CCD) camera. For this experiment an effective pixel size of 5 μm best met our needs to: (i) resolve the shape of individual snow grains; and (ii) maintain a sufficiently large field of view to observe a statistically significant volume.
Fig. 1. Vertical X-ray projection image for a snow sample (a), and its reconstructed central two-dimensional section (b) from the standard filtered back-projection software.
2.3. The 3-D volume data generation
Image analysis is performed using classical gray-level imaging and morphological algorithms in two and three dimensions. The aim of the procedure is to convert the gray-level images from the ESRF, containing three phases (ice, pore filler and remaining air bubbles), into a 3-D binary image where ice is distinguished from pore phase. The initial step is a volume image normalization by the first-neighbors 3-D averaging. In each plane, air bubbles are detected and their gray level is set to the average value of the pore space. The location of air bubbles is the union of the three following areas: the darkest one, the palest one and high-gradient zones owing to phase contrast during shots. Small holes inside bubbles are removed with a closing and hole-fill procedure (Reference SerraSerra, 1988). Noise speckles are removed using a two-dimensional median filter. Then an automatic threshold can be performed before a last 3-D smoothing of the image. Finally, the similarity of the binary image overlaid with the gray-level one is visually checked for each plane, and hand-made corrections may be made. The final result is a smoothed black-and-white 3-D image (Fig. 2).
Fig. 2. Three binary volumetric images from the isothermal metamorphism experiment, showing progressive change in morphology (resolution 5 μm, 6003 voxel volume) 15 hours (a) 470 hours (b) and 2026 hours (c) after snowfall.
Dry-snow evolution under low temperature gradient and wet-snow evolution are governed by local curvature (Reference ColbeckColbeck, 1997). A model has been developed for computing curvatures on discrete 3-D images (Reference Brzoska, Lesaffre, Coléou, Xu and PieritzBrzoska and others, 1999), the main results of which for the sample 15 hours (porosity 0.892; density 98 kgm–3) and the sample 2026 hours (porosity 0.715; density 260 kgm–3) are shown in Figure 3.
Fig. 3. Histograms of two curvature distributions: 15 and 2026 hours after snowfall (Fig. 2).
The pre-processor surface extraction capabilities allow the study of small grain geometries directly from the volumetric structure, such as the complex internal crystal channel (Fig. 4). This crystal (Fig. 4b; volume resolution 5 μm) was extracted by the “Tri-Linear Level Set Algorithm” (Reference SethianSethian, 1996) from the 15 hour binary microcomputed tomography (M-CT) image (Fig. 3), and compared to existing optical crystal images (Fig. 4a). The complex internal channels are observed by 3-D software rendering transparency from the external surface (Fig. 4b).
Fig. 4. A small ice crystal (a) observed by an optical microscopy, compared to a (b) small crystal surface (triangles surface) extracted by the pre-processor from volume 15 hours (binary M-CT image) and its (c) complex internal geometry (small channels, observed by surface transparency).
2.4. The RCFEA package
The new 3-D microstructure geometry modeling and FE analysis simulation package (RCFEA) was proposed as a complete environment to allow the execution of all steps from the 3-D raw binary image (or gray-value image) from a digital reconstruction (e.g. optical-computer tomography (O-CT) or X-ray microtomography (M-CT)) to the micromechanics simulation (thermal or mechanical). Like other developments (Reference RaabeRaabe, 1998), the package has different applications grouped in three main modules: pre-processor, processor and post-processor. All applications are developed with C++ technique standard code language using an object-oriented programming technique (OOP). These codes are portable and independent of the operational system.
2.4.1. The pre-processor
The pre-processor module is responsible for the control of all steps from the raw binary or gray volumetric image (input file from O-CTor M-CT) to the output description file for the simulation (processor module input). These main operations are:
-
(i) volumetric filters for identification of isolated objects, simplification, smooth, data shrinking;
-
(ii) surface triangles extraction (Reference SethianSethian, 1996) by linear and non-linear level set algorithms and raw voxel surface triangles;
-
(iii) triangles surface filter, decimation (simplification− triangles reduction) and quality control to filter overhead elements, holes and inclusions;
-
(iv) standard files output conversion to mesh generation;
-
(v) mesh interpretation for identification and labeling of different kind of surface nodes;
-
(vi) simulation and control output data file to allow operating the processor module (FE code and solver).
2.4.2. The processor
The processor module is an adapted implicit–explicit FE-analysis code where the computational domain is divided into FEs. The elements connect at nodes to form a 3-D domain. All fields (velocities, stresses, etc.) are assumed to be continuous over the domain. This differs from normal FE programs where quantities (e.g. stresses) can show jumps over element borders. The continuity requirement provides a numerical scheme that is more diffuse but also more stable. It has the further advantage that all fields can be stored in the nodes and not in integration points as in normal FE programs. This leads to a relatively easy implementation of remeshing and refining techniques. Only first order in time equations are solved. Time derivatives are approximated with Euler backward time discretization.
The main program capabilities are:
-
(i) differential equations for convection–diffusion, heat transfer, elasticity (isotropy and transverse isotropy), elasto-plasticity (Von Mises, Mohr–Coulomb, Gurson, etc.), hypo-plasticity (VonWolfersdorf, cohesion, intergranular strains, pressure-dependent initial void ratio) and damage;
-
(ii) parallelization (message-passing interface base) and iterative linear equations solver (diagonal preconditioned biconjugate gradient solver).
2.4.3. The post-processor
The post-processor module implements the interpretation and data meaning over the processor output results. The input files are the mesh and the simulated result file. The user can choose and visualize (by using a standard interface) the surface triangles, mesh tetrahedrals, mesh deformation, stress, strain, solver residues, node velocities, etc.
3. Snow Micromechanical Simulation
The hypothetical experiment is the uniaxial compression (Reference PopovPopov, 1968). For simplification, temperature and loading rate are assumed to be constant and the grains in the sample are initially bonded into a chain structure. The overall effect of the deformation is to strengthen the structure and increase the density of the snow. Concurrently, the bonding changes so the mechanical properties of the snow can vary through a wide range of values, depending on the deformation path. The macroscopic deformation of snow reflects the accumulated deformation on the scale of the grain. Usually, if tests on natural snow are of short duration, then strains are small and changes in bond structure limited. In such cases, the constitutive equations for linear-elastic, viscous or viscoelastic materials can be used to interpret the test results. In these cases, linear relationships are probably applicable over a relatively large range of stresses.
For our purposes, we only need to consider ideal elastic solids under isothermal conditions. An ideal elastic solid is defined by the property that the strain at any time depends only on the instantaneous magnitude of the stress and is independent of the stress history. Further, if the stress is removed from a deformed sample on an ideal elastic material, then the strain disappears, the sample returns to its original state and all the strain energy stored during deformation is recovered. However, determining a value of Young’s modulus (E) for snow from a static test in uniaxial loading includes the assumption that the stress–strain relation is linear, because the experimental data are fit to the one-dimensional form of Hooke’s law, σ = Eε, where σ is stress and ε is strain in the same direction. Similar arguments apply to the other elastic constants (the shear and bulk moduli, Poisson’s ratio, etc.).
3.1. Pre-processing: surface extraction and mesh generation
The snow data used to simulate the compression experiment is a subsample (3003 voxels) (Fig. 5a) from the sample 2046 hours after the snowfall. The snow consists of rounded grains, allowing the assumption that bonds and grains have the same crystal micromechanical properties as pure monocrystal ice (see section 3.2). The main steps in the pre-processor module (from the binary geometric volume data to the processor control file) are:
Fig. 5. The 3003 voxels: (a) subsample geometric volume data and (b) the numerical pre-processed surface cylinder for simulation.
-
(i) 3-D Gaussian filter to smooth and eliminate the small artifacts like 1–2 voxel bonds;
-
(ii) binary resolution shrinking (26 to 1 voxel filter threshold to minimize triangles surface elements);
-
(iii) non-connected small voxel objects filter to eliminate the top–bottom non-connected elements;
-
(iv) addition of top and bottom artificial plates (thickness: 17.5% in total volume voxel height for each plate) to uniform the uniaxial compression strength over the structure;
-
(v) cylinder data extraction to uniform radial stress distribution over the sample (Fig. 5b);
-
(vi) trilinear level set surface extraction (geometric surface description by triangles) (Reference SethianSethian, 1996);
-
(vii) surface smooth and triangles decimation;
-
(viii) filter small surface objects (grains internal porous) and export the surface file.
To perform the FE analysis, the reconstruction of an appropriatemesh in the microstructural region of interest is essential (Reference Frey and GeorgeFrey and George, 1998). Grains have complex shapes; they vary in size and orientation, and the surface curvature varies from one surface element to another, even on the same grain. Further, the distances between grains are not uniform; some grains are very closely packed, whereas others are not. To capture the effect of these geometric details on the local stress (or strain) distribution, the mesh size has to be very fine near the micro surfaces, and also in the regions of closely spaced pores. However, for computational efficiency it is not desirable to use a very fine mesh over the whole volume of the region of interest, as in a voxel-based FE model approach (Reference GarbocziGarboczi, 1998). A variable-size mesh has to be created in an interactive manner by tetrahedral inclusion from the triangulated surface geometry description.
From the exported surface file we generate the four-node tetrahedral mesh using a commercial software. The final mesh had (Fig. 6b) 30 540 nodes, 47530 surface triangles and 104731 tetrahedrals. It was interpreted by the preprocessor to recognize nodes (top and bottom surface nodes, surface inside nodes and volumetric inside nodes). This step allows the setting of physical domain properties and boundary conditions.
Fig. 6. The triangulated surface description: (a) extracted from the voxel volumetric data to the input simulation tetrahedral mesh (b).
3.2. Boundary conditions, simulation and post-processing
We assume the snow grains are ice grains, where the elastic behavior of ice is characterized by moderate anisotropy (Reference SchulsonSchulson, 1999). At temperatures near the melting point, Young’s modulus of single crystals varies by <30%, from 12 GPa along the least compliant direction (parallel to the c axis) to 8.6 GPa along the most compliant direction (inclined to both the c and a axes). Along directions within the basal plane, Young’s modulus is 10 GPa. For randomly oriented polycrystals, typical values of Young’s modulus and Poisson’s ratio are 9.0 GPa and 0.33 at –5°C.
The FE simulation is performed over a central subcylindrical sample inscribed in a cube of 3003 voxels of the 3-D geometry of the natural snow sample taken 2046 hours after snowfall (Fig. 5). The subsample size is lower than the representative elementary volume relevant to micromechanical properties.
The basic constitutive relationships and boundary conditions are:
-
(i) the constitutive equation used in the simulation is purely linear-elastic;
-
(ii) unidirectional top-down compression of 40 Pa with the immobile bottom surface;
-
(iii) four-node linear tetrahedron elements are used for meshing;
-
(iv) the FE model consists of 30 540 nodes, 47530 surface triangles (Fig. 6a) and 104731 tetrahedrals in a 3375 mm3 frame; the density of the FE mesh in the grains (Fig. 6b) is sufficiently high to represent steep local stress and strain gradients in the material.
In Figure 7, maximum normal stress is plotted. Peak stresses are found in the ice bonds (colored red). The maximum stress for this sample was 6–9 MPa in some small bonds. These preliminary results show that high stresses are reached in small bonds, which make a fundamental contribution to snow failure.
Fig. 7. Simulated compression stress over the cylindrical geometry (a, b) and the clipping mesh (c). The color scale shows the maximal normal stress in the range 0–1MPa (for visualization, greater values are plotted with the same maximum color). The maximum stress is 6–9MPa in some bonds (plotted by red color).
4. Discussion
The above results are only qualitative but they confirm the potential of this numerical method to perform microstructural analysis and simulate the micromechanical behavior of these complex materials.
The next step will be to increase the computational data domain to perform the simulations over a representative elementary volume for these scales. For that reason, we will improve the quality of the volumetric image techniques and the numerical simulation results, such as the 3-D image analysis and gray-level image threshold, surface extraction algorithm (Reference Flin, Brzoska, Lesaffre, Coléou and LamboleyFlin and others, 2001) and a special discrete Voronoi mesh generation (Reference Chassery and MontanvertChassery and Montanvert, 1991). Also, we plan in the future to validate the numerical procedure by experiments.
5. Conclusions
This paper introduced the first developments in the analysis of the microstructure volume data from the 3-D raw binary image from digital reconstruction (like O-CTor M-CT) to the micromechanical simulation using image analysis and volumetric surface extraction techniques. The algorithms and the numerical methods were introduced to simulate a hypothetical compression experiment on a real 3-D snow sample.
The micromechanical simulation results are only qualitative results, but they show the potential of using the microscopic geometric description of the medium, to study, for example, the grain-bond chain influence on the macroscopic properties of snow.
Acknowledgement
This research has been supported by the visiting scientist program fund of the Centre National de Recherches Météorologiques/Météo France. The financial support is gratefully acknowledged.
