Thermolab: A Thermodynamics Laboratory for Nonlinear Transport Processes in Open Systems

We developed a numerical thermodynamics laboratory called “Thermolab” to study the effects of the thermodynamic behavior of nonideal solution models on reactive transport processes in open systems. The equations of the state of internally consistent thermodynamic data sets are implemented in MATLAB functions and form the basis for calculating Gibbs energy. A linear algebraic approach is used in Thermolab to compute Gibbs energy of mixing for multicomponent phases to study the impact of the nonideality of solution models on transport processes. The Gibbs energies are benchmarked with experimental data, phase diagrams, and other thermodynamic software. Constrained Gibbs minimization is exemplified with MATLAB codes and iterative refinement of composition of mixtures may be used to increase precision and accuracy. All needed transport variables such as densities, phase compositions, and chemical potentials are obtained from Gibbs energy of the stable phases after the minimization in Thermolab. We demonstrate the use of precomputed local equilibrium data obtained with Thermolab in reactive transport models. In reactive fluid flow the shape and the velocity of the reaction front vary depending on the nonlinearity of the partitioning of a component in fluid and solid. We argue that nonideality of solution models has to be taken into account and further explored in reactive transport models. Thermolab Gibbs energies can be used in Cahn-Hilliard models for nonlinear diffusion and phase growth. This presents a transient process toward equilibrium and avoids computational problems arising during precomputing of equilibrium data

A spectral/finite difference method for simulating large deformations of heterogeneous, viscoelastic materials

A numerical algorithm is presented that simulates large deformations of heterogeneous, viscoelastic materials in two dimensions. The algorithm is based on a spectral/finite difference method and uses the Eulerian formulation including objective derivatives ofthe stress tensor in the rheological equations. The viscoelastic rheology is described bythe linear Maxwell model, which consists of an elastic and viscous element connected inseries. The algorithm is especially suitable to simulate periodic instabilities. The derivatives in the direction of periodicity are approximated by spectral expansions, whereas the derivatives in the direction orthogonal to the periodicity are approximated by finite differences. The 1‐D Eulerian finite difference grid consists of centre and nodal points and has variable grid spacing. Time derivatives are approximated with finite differences using an implicit strategy with a variable time step. The performance of the numerical code is demonstrated by calculation, for the first time, of the pressure field evolution during folding of viscoelastic multilayers. The algorithm is stable for viscosity contrasts up to 5 × 105, which demonstrates that spectral methods can be used to simulate dynamical systems involving large material heterogeneities. The successful simulations show that combined spectral/finite difference methods using the Eulerian formulation are a promising tool to simulate mechanical processes that involve large deformations, viscoelastic rheologies and strong material heterogeneitie

Spontaneous thermal runaway as an ultimate failure mechanism of materials

The first theoretical estimate of the shear strength of a perfect crystal was given by Frenkel [Z. Phys. 37, 572 (1926)]. He assumed that as slip occurred, two rigid atomic rows in the crystal would move over each other along a slip plane. Based on this simple model, Frenkel derived the ultimate shear strength to be about one tenth of the shear modulus. Here we present a theoretical study showing that catastrophic material failure may occur below Frenkel's ultimate limit as a result of thermal runaway. We demonstrate that the condition for thermal runaway to occur is controlled by only two dimensionless variables and, based on the thermal runaway failure mechanism, we calculate the maximum shear strength $\sigma_c$ of viscoelastic materials. Moreover, during the thermal runaway process, the magnitude of strain and temperature progressively localize in space producing a narrow region of highly deformed material, i.e. a shear band. We then demonstrate the relevance of this new concept for material failure known to occur at scales ranging from nanometers to kilometers.Comment: 4 pages, 3 figures. Eq. (6) and Fig. 2a corrected; added references; improved quality of figure

An explicit GPU-based material point method solver for elastoplastic problems (ep2-3De v1.0)

We propose an explicit GPU-based solver within the material point method (MPM) framework using graphics processing units (GPUs) to resolve elastoplastic problems under two- and three-dimensional configurations (i.e. granular collapses and slumping mechanics). Modern GPU architectures, including Ampere, Turing and Volta, provide a computational framework that is well suited to the locality of the material point method in view of high-performance computing. For intense and non-local computational aspects (i.e. the back-and-forth mapping between the nodes of the background mesh and the material points), we use straightforward atomic operations (the scattering paradigm). We select the generalized interpolation material point method (GIMPM) to resolve the cell-crossing error, which typically arises in the original MPM, because of the C0 continuity of the linear basis function. We validate our GPU-based in-house solver by comparing numerical results for granular collapses with the available experimental data sets. Good agreement is found between the numerical results and experimental results for the free surface and failure surface. We further evaluate the performance of our GPU-based implementation for the three-dimensional elastoplastic slumping mechanics problem. We report (i) a maximum 200-fold performance gain between a CPU- and a single-GPU-based implementation, provided that (ii) the hardware limit (i.e. the peak memory bandwidth) of the device is reached. Furthermore, our multi-GPU implementation can resolve models with nearly a billion material points. We finally showcase an application to slumping mechanics and demonstrate the importance of a three-dimensional configuration coupled with heterogeneous properties to resolve complex material behaviour.</p

Melt Migration and Chemical Differentiation by Reactive Porosity Waves

Melt transport across the ductile mantle is essential for oceanic crust formation or intraplate volcanism. However, mechanisms of melt migration and associated chemical interaction between melt and solid mantle remain unclear. Here, we present a thermo-hydro-mechanical-chemical (THMC) model for melt migration coupled to chemical differentiation. We consider melt migration by porosity waves and a chemical system of forsterite-fayalite-silica. We solve the one-dimensional (1D) THMC model numerically using the finite difference method. Variables, such as solid and melt densities or MgO and SiO2 mass concentrations, are functions of pressure (P), temperature (T), and total silica mass fraction (urn:x-wiley:15252027:media:ggge22741:ggge22741-math-0001). These variables are pre-computed with Gibbs energy minimization and their variations with evolving P, T, and urn:x-wiley:15252027:media:ggge22741:ggge22741-math-0002 are implemented in the THMC model. We consider P and T conditions relevant around the lithosphere-asthenosphere boundary. Systematic 1D simulations quantify the impact of initial distributions of porosity and urn:x-wiley:15252027:media:ggge22741:ggge22741-math-0003 on the melt velocity. Larger perturbations of urn:x-wiley:15252027:media:ggge22741:ggge22741-math-0004 cause larger melt velocities. An adiabatic or conductive geotherm cause fundamentally different vertical variations of densities and concentrations, and an adiabatic geotherm generates higher melt velocities. We quantify differences between melt transport (considering incompatible tracers), major element transport and porosity evolution. Melt transport is significant in the models. We also quantify the relative importance of four porosity variation mechanisms: (a) mechanical compaction and decompaction, (b) density variation, (c) compositional variation, and (d) solid-melt mass exchange. In the models, (de)compaction dominates the porosity variation. We further discuss preliminary results of 2D THMC simulations showing blob-like and channel-like porosity waves

The problem of depth in geology: When pressure does not translate into depth

We review published evidence that rocks can develop, sustain and record significant pressure deviations from lithostatic values. Spectroscopic studies at room pressure and temperature (P-T) reveal that in situ pressure variations in minerals can reach GPa levels. Rise of confined pressure leads to higher amplitude of these variations documented by the preservation of α-quartz incipiently amorphized under pressure (IAUP quartz), which requires over 12 GPa pressure variations at the grain scale. Formation of coesite in rock-deformation experiments at lower than expected confined pressures confirmed the presence of GPa-level pressure variations at elevated temperatures and pressures within deforming and reacting multi-mineral and polycrystalline rock samples. Whiteschists containing garnet porphyroblasts formed during prograde metamorphism that host quartz inclusions in their cores and coesite inclusions in their rims imply preservation of large differences in pressure at elevated pressure and temperature. Formation and preservation of coherent cryptoperthite exsolution lamellae in natural alkali feldspar provides direct evidence for grain-scale, GPa-level stress variations at 680°C at geologic time scales from peak to ambient P-T conditions. Similarly, but in a more indirect way, the universally accepted' pressure-vessel' model to explain preservation of coesite, diamond and other ultra-high-pressure indicators requires GPa-level pressure differences between the inclusion and the host during decompression at temperatures sufficiently high for these minerals to transform into their lower pressure polymorphs even at laboratory time scales. A variety of mechanisms can explain the formation and preservation of pressure variations at various length scales. These mechanisms may double the pressure value compared to the lithostatic in compressional settings, and pressures up to two times the lithostatic value were estimated under special mechanical conditions. We conclude, based on these considerations, that geodynamic scenarios involving very deep subduction processes with subsequent very rapid exhumation from a great depth must be viewed with due caution when one seeks to explain the presence of microscopic ultrahigh-pressure mineralogical indicators in rocks. Non-lithostatic interpretation of high-pressure indicators may potentially resolve long-lasting geological conundrum

Spontaneous dissipation of elastic energy by self-localizing thermal runaway

Thermal runaway instability induced by material softening due to shear heating represents a potential mechanism for mechanical failure of viscoelastic solids. In this work we present a model based on a continuum formulation of a viscoelastic material with Arrhenius dependence of viscosity on temperature, and investigate the behavior of the thermal runaway phenomenon by analytical and numerical methods. Approximate analytical descriptions of the problem reveal that onset of thermal runaway instability is controlled by only two dimensionless combinations of physical parameters. Numerical simulations of the model independently verify these analytical results and allow a quantitative examination of the complete time evolutions of the shear stress and the spatial distributions of temperature and displacement during runaway instability. Thus we find that thermal runaway processes may well develop under nonadiabatic conditions. Moreover, nonadiabaticity of the unstable runaway mode leads to continuous and extreme localization of the strain and temperature profiles in space, demonstrating that the thermal runaway process can cause shear banding. Examples of time evolutions of the spatial distribution of the shear displacement between the interior of the shear band and the essentially nondeforming material outside are presented. Finally, a simple relation between evolution of shear stress, displacement, shear-band width and temperature rise during runaway instability is given.Comment: 16 pages, 7 figures. Extended conclusion; added reference

Relationship between tectonic overpressure, deviatoric stress, driving force, isostasy and gravitational potential energy

ISSN:0956-540XISSN:1365-246

Failure patterns caused by localized rise in pore-fluid overpressure and effective strength of rocks

In order to better understand the interaction between pore-fluid overpressure and failure patterns in rocks we consider a porous elasto-plastic medium in which a laterally localized overpressure line source is imposed at depth below the free surface. We solve numerically the fluid filtration equation coupled to the gravitational force balance and poro-elasto-plastic rheology equations. Systematic numerical simulations, varying initial stress, intrinsic material properties and geometry, show the existence of five distinct failure patterns caused by either shear banding or tensile fracturing. The value of the critical pore-fluid overpressure at the onset of failure is derived from an analytical solution that is in excellent agreement with numerical simulations. Finally, we construct a phase-diagram that predicts the domains of the different failure patterns and at the onset of failure