Numerical modelling of double diffusion driven reactive flow transport in deformable fluid-saturated porous media with particular consideration of temperature-dependent chemical reaction rates

Numerical methods ave used to solve double diffusion driven reactive flow transport problems in deformable fluid-saturated porous media. in particular, thp temperature dependent reaction rate in the non-equilibrium chemical reactions is considered. A general numerical solution method, which is a combination of the finite difference method in FLAG and the finite element method in FIDAP, to solve the fully coupled problem involving material deformation, pore-fluid flow, heat transfer and species transport/chemical reactions in deformable fluid-saturated porous media has been developed The coupled problem is divided into two subproblems which are solved interactively until the convergence requirement is met. Owing to the approximate nature of the numerical method, if is essential to justify the numerical solutions through some kind of theoretical analysis. This has been highlighted in this paper The related numerical results, which are justified by the theoretical analysis, have demonstrated that the proposed solution method is useful for and applicable to a wide range of fully coupled problems in the field of science and engineering.

Numerical modelling of single-layer folding: clarification of an issue regarding the possible effect of computer codes and the influence of initial irregularities

The influence of initial perturbation geometry and material propel-ties on final fold geometry has been investigated using finite-difference (FLAC) and finite-element (MARC) numerical models. Previous studies using these two different codes reported very different folding behaviour although the material properties, boundary conditions and initial perturbation geometries were similar. The current results establish that the discrepancy was not due to the different computer codes but due to the different strain rates employed in the two previous studies (i.e. 10(-6) s(-1) in the FLAC models and 10(-14) s(-1) in the MARC models). As a result, different parts of the elasto-viscous rheological field were bring investigated. For the same material properties, strain rate and boundary conditions, the present results using the two different codes are consistent. A transition in Folding behaviour, from a situation where the geometry of initial perturbation determines final fold shape to a situation where material properties control the final geometry, is produced using both models. This transition takes place with increasing strain rate, decreasing elastic moduli or increasing viscosity (reflecting in each case the increasing influence of the elastic component in the Maxwell elastoviscous rheology). The transition described here is mechanically feasible but is associated with very high stresses in the competent layer (on the order of GPa), which is improbable under natural conditions. (C) 2000 Elsevier Science Ltd. All rights reserved.

Prediction of rock alteration patterns: a potential tool in mineral exploration

Hydrothermal alteration of a quartz-K-feldspar rock is simulated numerically by coupling fluid flow and chemical reactions. Introduction of CO2 gas generates an acidic fluid and produces secondary quartz, muscovite and/or pyrophyllite at constant temperature and pressure of 300 degrees C and 200 MPa. The precipitation and/or dissolution of the secondary minerals is controlled by either mass-action relations or rate laws. In our simulations the mass of the primary elements are conserved and the mass-balance equations are solved sequentially using an implicit scheme in a finite-element code. The pore-fluid velocity is assumed to be constant. The change of rock volume due to the dissolution or precipitation of the minerals, which is directly related to their molar volume, is taken into account. Feedback into the rock porosity and the reaction rates is included in the model. The model produces zones of pyrophyllite quartz and muscovite due to the dissolution of K-feldspar. Our model simulates, in a simplified way, the acid-induced alteration assemblages observed in various guises in many significant mineral deposits. The particular aluminosilicate minerals produced in these experiments are associated with the gold deposits of the Witwatersrand Basin.

A numerical study of flexural buckling of foliated rock slopes

The occurrence of foliated rock masses is common in mining environment. Methods employing continuum approximation in describing the deformation of such rock masses possess a clear advantage over methods where each rock layer and each inter-layer interface (joint) is explicitly modelled. In devising such a continuum model it is imperative that moment (couple) stresses and internal rotations associated with the bending of the rock layers be properly incorporated in the model formulation. Such an approach will lead to a Cosserat-type theory. In the present model, the behaviour of the intact rock layer is assumed to be linearly elastic and the joints are assumed to be elastic-perfectly plastic. Condition of slip at the interfaces are determined by a Mohr-Coulomb criterion with tension cut off at zero normal stress. The theory is valid for large deformations. The model is incorporated into the finite element program AFENA and validated against an analytical solution of elementary buckling problems of a layered medium under gravity loading. A design chart suitable for assessing the stability of slopes in foliated rock masses against flexural buckling failure has been developed. The design chart is easy to use and provides a quick estimate of critical loading factors for slopes in foliated rock masses. It is shown that the model based on Euler's buckling theory as proposed by Cavers (Rock Mechanics and Rock Engineering 1981; 14:87-104) substantially overestimates the critical heights for a vertical slope and underestimates the same for sub-vertical slopes. Copyright (C) 2001 John Wiley & Sons, Ltd.

Computer simulations of coupled problems in geological and geochemical systems

The finite element method is used to simulate coupled problems, which describe the related physical and chemical processes of ore body formation and mineralization, in geological and geochemical systems. The main purpose of this paper is to illustrate some simulation results for different types of modelling problems in pore-fluid saturated rock masses. The aims of the simulation results presented in this paper are: (1) getting a better understanding of the processes and mechanisms of ore body formation and mineralization in the upper crust of the Earth; (2) demonstrating the usefulness and applicability of the finite element method in dealing with a wide range of coupled problems in geological and geochemical systems; (3) qualitatively establishing a set of showcase problems, against which any numerical method and computer package can be reasonably validated. (C) 2002 Published by Elsevier Science B.V.

A director theory for visco-elastic folding instabilities in multilayered rock

A model for finely layered visco-elastic rock proposed by us in previous papers is revisited and generalized to include couple stresses. We begin with an outline of the governing equations for the standard continuum case and apply a computational simulation scheme suitable for problems involving very large deformations. We then consider buckling instabilities in a finite, rectangular domain. Embedded within this domain, parallel to the longer dimension we consider a stiff, layered beam under compression. We analyse folding up to 40% shortening. The standard continuum solution becomes unstable for extreme values of the shear/normal viscosity ratio. The instability is a consequence of the neglect of the bending stiffness/viscosity in the standard continuum model. We suggest considering these effects within the framework of a couple stress theory. Couple stress theories involve second order spatial derivatives of the velocities/displacements in the virtual work principle. To avoid C-1 continuity in the finite element formulation we introduce the spin of the cross sections of the individual layers as an independent variable and enforce equality to the spin of the unit normal vector to the layers (-the director of the layer system-) by means of a penalty method. We illustrate the convergence of the penalty method by means of numerical solutions of simple shears of an infinite layer for increasing values of the penalty parameter. For the shear problem we present solutions assuming that the internal layering is oriented orthogonal to the surfaces of the shear layer initially. For high values of the ratio of the normal-to the shear viscosity the deformation concentrates in thin bands around to the layer surfaces. The effect of couple stresses on the evolution of folds in layered structures is also investigated. (C) 2002 Elsevier Science Ltd. All rights reserved.

Large amplitude folding in finely layered viscoelastic rock structures

We analyze folding phenomena in finely layered viscoelastic rock. Fine is meant in the sense that the thickness of each layer is considerably smaller than characteristic structural dimensions. For this purpose we derive constitutive relations and apply a computational simulation scheme (a finite-element based particle advection scheme; see MORESI et al., 2001) suitable for problems involving very large deformations of layered viscous and viscoelastic rocks. An algorithm for the time integration of the governing equations as well as details of the finite-element implementation is also given. We then consider buckling instabilities in a finite, rectangular domain. Embedded within this domain, parallel to the longer dimension we consider a stiff, layered plate. The domain is compressed along the layer axis by prescribing velocities along the sides. First, for the viscous limit we consider the response to a series of harmonic perturbations of the director orientation. The Fourier spectra of the initial folding velocity are compared for different viscosity ratios. Turning to the nonlinear regime we analyze viscoelastic folding histories up to 40% shortening. The effect of layering manifests itself in that appreciable buckling instabilities are obtained at much lower viscosity ratios (1:10) as is required for the buckling of isotropic plates (1:500). The wavelength induced by the initial harmonic perturbation of the director orientation seems to be persistent. In the section of the parameter space considered here elasticity seems to delay or inhibit the occurrence of a second, larger wavelength. Finally, in a linear instability analysis we undertake a brief excursion into the potential role of couple stresses on the folding process. The linear instability analysis also provides insight into the expected modes of deformation at the onset of instability, and the different regimes of behavior one might expect to observe.

Mantle convection modeling with viscoelastic/brittle lithosphere: Numerical methodology and plate tectonic modeling

The earth's tectonic plates are strong, viscoelastic shells which make up the outermost part of a thermally convecting, predominantly viscous layer. Brittle failure of the lithosphere occurs when stresses are high. In order to build a realistic simulation of the planet's evolution, the complete viscoelastic/brittle convection system needs to be considered. A particle-in-cell finite element method is demonstrated which can simulate very large deformation viscoelasticity with a strain-dependent yield stress. This is applied to a plate-deformation problem. Numerical accuracy is demonstrated relative to analytic benchmarks, and the characteristics of the method are discussed.

Analysis of steady-state heat transfer through mid-crustal vertical cracks with upward throughflow in hydrothermal systems

We conduct a theoretical analysis of steady-state heat transfer problems through mid-crustal vertical cracks with upward throughflow in hydrothermal systems. In particular, we derive analytical solutions for both the far field and near field of the system. In order to investigate the contribution of the forced advection to the total temperature of the system, two concepts, namely the critical Peclet number and the critical permeability of the system, have been presented and discussed in this paper. The analytical solution for the far field of the system indicates that if the pore-fluid pressure gradient in the crust is lithostatic, the critical permeability of the system can be used to determine whether or not the contribution of the forced advection to the total temperature of the system is negligible. Otherwise, the critical Peclet number should be used. For a crust of moderate thickness, the critical permeability is of the order of magnitude of 10(-20) m(2), under which heat conduction is the overwhelming mechanism to transfer heat energy, even though the pore-fluid pressure gradient in the crust is lithostatic. Furthermore, the lower bound analytical solution for the near field of the system demonstrates that the permeable vertical cracks in the middle crust can efficiently transfer heat energy from the lower crust to the upper crust of the Earth. Copyright (C) 2002 John Wiley Sons, Ltd.

Effects of hot intrusions on pore-fluid flow and heat transfer in fluid-saturated rocks

We use the finite element method to solve coupled problems between pore-fluid flow and heat transfer in fluid-saturated porous rocks. In particular, we investigate the effects of both the hot pluton intrusion and topographically driven horizontal flow on the distributions of the pore-flow velocity and temperature in large-scale hydrothermal systems. Since general mineralization patterns are strongly dependent on distributions of both the pore-fluid velocity and temperature fields, the modern mineralization theory has been used to predict the general mineralization patterns in several realistic hydrothermal systems. The related numerical results have demonstrated that: (1) The existence of a hot intrusion can cause an increase in the maximum value of the pore-fluid velocity in the hydrothermal system. (2) The permeability of an intruded pluton is one of the sensitive parameters to control the pore-fluid flow, heat transfer and ore body formation in hydrothermal systems. (3) The maximum value of the pore-fluid velocity increases when the bottom temperature of the hydrothermal system is increased. (4) The topographically driven flow has significant effects on the pore-fluid flow, temperature distribution and precipitation pattern of minerals in hydrothermal systems. (5) The size of the computational domain may have some effects on the pore-fluid flow and heat transfer, indicating that the size of a hydrothermal system may affect the pore-fluid flow and heat transfer within the system. (C) 2003 Elsevier Science B.V. All rights reserved.

An equivalent algorithm for simulating thermal effects of magma intrusion problems in porous rocks

An equivalent algorithm is proposed to simulate thermal effects of the magma intrusion in geological systems, which are composed of porous rocks. Based on the physical and mathematical equivalence, the original magma solidification problem with a moving boundary between the rock and intruded magma is transformed into a new problem without the moving boundary but with a physically equivalent heat source. From the analysis of an ideal solidification model, the physically equivalent heat source has been determined in this paper. The major advantage in using the proposed equivalent algorithm is that the fixed finite element mesh with a variable integration time step can be employed to simulate the thermal effect of the intruded magma solidification using the conventional finite element method. The related numerical results have demonstrated the correctness and usefulness of the proposed equivalent algorithm for simulating the thermal effect of the intruded magma solidification in geological systems. (C) 2003 Elsevier B.V. All rights reserved.

Anisotropic convection model for the earth's mantle

The paper presents a theory for modeling flow in anisotropic, viscous rock. This theory has originally been developed for the simulation of large deformation processes including the folding and kinking of multi-layered visco-elastic rock (Muhlhaus et al. [1,2]). The orientation of slip planes in the context of crystallographic slip is determined by the normal vector - the director - of these surfaces. The model is applied to simulate anisotropic mantle convection. We compare the evolution of flow patterns, Nusselt number and director orientations for isotropic and anisotropic rheologies. In the simulations we utilize two different finite element methodologies: The Lagrangian Integration Point Method Moresi et al [8] and an Eulerian formulation, which we implemented into the finite element based pde solver Fastflo (www.cmis.csiro.au/Fastflo/). The reason for utilizing two different finite element codes was firstly to study the influence of an anisotropic power law rheology which currently is not implemented into the Lagrangian Integration point scheme [8] and secondly to study the numerical performance of Eulerian (Fastflo)- and Lagrangian integration schemes [8]. It turned out that whereas in the Lagrangian method the Nusselt number vs time plot reached only a quasi steady state where the Nusselt number oscillates around a steady state value the Eulerian scheme reaches exact steady states and produces a high degree of alignment (director orientation locally orthogonal to velocity vector almost everywhere in the computational domain). In the simulations emergent anisotropy was strongest in terms of modulus contrast in the up and down-welling plumes. Mechanisms for anisotropic material behavior in the mantle dynamics context are discussed by Christensen [3]. The dominant mineral phases in the mantle generally do not exhibit strong elastic anisotropy but they still may be oriented by the convective flow. Thus viscous anisotropy (the main focus of this paper) may or may not correlate with elastic or seismic anisotropy.

Convective instability of 3-D fluid-saturated geological fault zones heated from below

We conduct a theoretical analysis to investigate the convective instability of 3-D fluid-saturated geological fault zones when they are heated uniformly from below. In particular, we have derived exact analytical solutions for the critical Rayleigh numbers of different convective flow structures. Using these critical Rayleigh numbers, three interesting convective flow structures have been identified in a geological fault zone system. It has been recognized that the critical Rayleigh numbers of the system have a minimum value only for the fault zone of infinite length, in which the corresponding convective flow structure is a 2-D slender-circle flow. However, if the length of the fault zone is finite, the convective flow in the system must be 3-D. Even if the length of the fault zone is infinite, since the minimum critical Rayleigh number for the 2-D slender-circle flow structure is so close to that for the 3-D convective flow structure, the system may have almost the same chance to pick up the 3-D convective flow structures. Also, because the convection modes are so close for the 3-D convective flow structures, the convective flow may evolve into the 3-D finger-like structures, especially for the case of the fault thickness to height ratio approaching zero. This understanding demonstrates the beautiful aspects of the present analytical solution for the convective instability of 3-D geological fault zones, because the present analytical solution is valid for any value of the ratio of the fault height to thickness. Using the present analytical solution, the conditions, under which different convective flow structures may take place, can be easily determined.

Theoretical and numerical analyses of convective instability in porous media with temperature-dependent viscosity

Exact analytical solutions of the critical Rayleigh numbers have been obtained for a hydrothermal system consisting of a horizontal porous layer with temperature-dependent viscosity. The boundary conditions considered are constant temperature and zero vertical Darcy velocity at both the top and bottom of the layer. Not only can the derived analytical solutions be readily used to examine the effect of the temperature-dependent viscosity on the temperature-gradient driven convective flow, but also they can be used to validate the numerical methods such as the finite-element method and finite-difference method for dealing with the same kind of problem. The related analytical and numerical results demonstrated that the temperature-dependent viscosity destabilizes the temperature-gradient driven convective flow and therefore, may affect the ore body formation and mineralization in the upper crust of the Earth. Copyright (C) 2003 John Wiley Sons, Ltd.

Determination of mechanical properties of PECVD silicon nitride thin films for tunable MEMS Fabry-Pérot optical filters

This paper reports an investigation on techniques for determining elastic modulus and intrinsic stress gradient in plasma-enhanced chemical vapor deposition (PECVD) silicon nitride thin films. The elastic property of the silicon nitride thin films was determined using the nanoindentation method on silicon nitride/silicon bilayer systems. A simple empirical formula was developed to deconvolute the film elastic modulus. The intrinsic stress gradient in the films was determined by using micrometric cantilever beams, cross-membrane structures and mechanical simulation. The deflections of the silicon nitride thin film cantilever beams and cross-membranes caused by in-thickness stress gradients were measured using optical interference microscopy. Finite-element beam models were built to compute the deflection induced by the stress gradient. Matching the deflection computed under a given gradient with that measured experimentally on fabricated samples allows the stress gradient of the PECVD silicon nitride thin films introduced from the fabrication process to be evaluated.