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.

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.

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.

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.

Numerical modelling of chemical effects of magma solidification problems in porous rocks

The solidification of intruded magma in porous rocks can result in the following two consequences: (1) the heat release due to the solidification of the interface between the rock and intruded magma and (2) the mass release of the volatile fluids in the region where the intruded magma is solidified into the rock. Traditionally, the intruded magma solidification problem is treated as a moving interface (i.e. the solidification interface between the rock and intruded magma) problem to consider these consequences in conventional numerical methods. This paper presents an alternative new approach to simulate thermal and chemical consequences/effects of magma intrusion in geological systems, which are composed of porous rocks. In the proposed new approach and algorithm, 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 the proposed mass source and physically equivalent heat source. The major advantage in using the proposed equivalent algorithm is that a fixed mesh of finite elements with a variable integration time-step can be employed to simulate the consequences and effects of the intruded magma solidification using the conventional finite element method. The correctness and usefulness of the proposed equivalent algorithm have been demonstrated by a benchmark magma solidification problem. Copyright (c) 2005 John Wiley & Sons, Ltd.

A numerical study of effect of grain boundaries on elastic and plastic properties in nanocomposite materials

Nanocomposite materials have received considerable attention in recent years due to their novel properties. Grain boundaries are considered to play an important role in nanostructured materials. This work focuses on the finite element analysis of the effect of grain boundaries on the overall mechanical properties of aluminium/alumina composites. A grain boundary is incorporated into the commonly used unit cell model to investigate its effect on material properties. By combining the unit cell model with an indentation model, coupled with experimental indentation measurements, the ''effective'' plastic property of the grain boundary is estimated. In addition, the strengthening mechanism is also discussed based on the Estrin-Mecking model.

Modelling of cold-formed purlin-sheeting systems .1. Full model

Purlin-sheeting systems used for roofs and walls commonly take the form of cold-formed channel or zed section purlins, screw-connected to corrugated sheeting. These purlin-sheeting systems have been the subject of numerous theoretical and experimental investigations over the past three decades, but the complexity of the systems has led to great difficulty in developing a sound and general model. This paper presents a non-linear elasto-plastic finite element model, capable of predicting the behaviour of purlin-sheeting systems without the need for either experimental input or over simplifying assumptions. The model incorporates both the sheeting and the purlin, and is able to account for cross-sectional distortion of the purlin, the flexural and membrane restraining effects of the sheeting, and failure of the purlin by local buckling or yielding. The validity of the model is shown by its good correlation with experimental results. A simplified version of this model, which is more suitable for use in a design environment, is presented in a companion paper. (C) 1997 Elsevier Science Ltd.

Modelling of cold-formed purlin-sheeting systems .2. Simplified model

A number of theoretical and experimental investigations have been made into the nature of purlin-sheeting systems over the past 30 years. These systems commonly consist of cold-formed zed or channel section purlins, connected to corrugated sheeting. They have proven difficult to model due to the complexity of both the purlin deformation and the restraint provided to the purlin by the sheeting. Part 1 of this paper presented a non-linear elasto plastic finite element model which, by incorporating both the purlin and the sheeting in the analysis, allowed the interaction between the two components of the system to be modelled. This paper presents a simplified version of the first model which has considerably decreased requirements in terms of computer memory, running time and data preparation. The Simplified Model includes only the purlin but allows for the sheeting's shear and rotational restraints by modelling these effects as springs located at the purlin-sheeting connections. Two accompanying programs determine the stiffness of these springs numerically. As in the Full Model, the Simplified Model is able to account for the cross-sectional distortion of the purlin, the shear and rotational restraining effects of the sheeting, and failure of the purlin by local buckling or yielding. The model requires no experimental or empirical input and its validity is shown by its goon con elation with experimental results. (C) 1997 Elsevier Science Ltd.

Elastic moduli of model random three-dimensional closed-cell cellular solids

Most cellular solids are random materials, while practically all theoretical structure-property results are for periodic models. To be able to generate theoretical results for random models, the finite element method (FEM) was used to study the elastic properties of solids with a closed-cell cellular structure. We have computed the density (rho) and microstructure dependence of the Young's modulus (E) and Poisson's ratio (PR) for several different isotropic random models based on Voronoi tessellations and level-cut Gaussian random fields. The effect of partially open cells is also considered. The results, which are best described by a power law E infinity rho (n) (1<n<2), show the influence of randomness and isotropy on the properties of closed-cell cellular materials, and are found to be in good agreement with experimental data. (C) 2001 Acta Materialia Inc. Published by Elsevier Science Ltd. All rights reserved.

Earth pressure on cantilever walls at design retained heights

There are many methods for the analysis and design of embedded cantilever retaining walls. They involve various different simplifications of the pressure distribution to allow calculation of the limiting equilibrium retained height and the bending moment when the retained height is less than the limiting equilibrium value, i.e. the serviceability case. Recently, a new method for determining the serviceability earth pressure and bending moment has been proposed. This method makes an assumption defining the point of zero net pressure. This assumption implies that the passive pressure is not fully mobilised immediately below the excavation level. The finite element analyses presented in this paper examine the net pressure distribution on walls in which the retained height is less, than the limiting equilibrium value. The study shows that for all practical walls, the earth pressure distributions on the front and back of the wall are at their limit values, Kp and K-a respectively, when the lumped factor of safety F-r is less than or equal to2.0. A rectilinear net pressure distribution is proposed that is intuitively logical. It produces good predictions of the complete bending moment diagram for walls in the service configuration and the proposed method gives results that have excellent agreement with centrifuge model tests. The study shows that the method for determining the serviceability bending moment suggested by Padfield and Mair(1) in the CIRIA Report 104 gives excellent predictions of the maximum bending moment in practical cantilever walls. It provides the missing data that have been needed to verify and justify the CIRIA 104 method.

Design, modelling and analysis of a six component force balance for hypervelocity wind tunnel testing

A combination of modelling and analysis techniques was used to design a six component force balance. The balance was designed specifically for the measurement of impulsive aerodynamic forces and moments characteristic of hypervelocity shock tunnel testing using the stress wave force measurement technique. Aerodynamic modelling was used to estimate the magnitude and distribution of forces and finite element modelling to determine the mechanical response of proposed balance designs. Simulation of balance performance was based on aerodynamic loads and mechanical responses using convolution techniques. Deconvolution was then used to assess balance performance and to guide further design modifications leading to the final balance design. (C) 2001 Elsevier Science Ltd. All rights reserved.

Analysis of a circular patch antenna radiating in a parallel-plate radial guide

A field matching method is described to analyze a recessed circular cavity radiating into a radial waveguide. Using the wall impedance approach, the analysis is divided into two separate problems of the cavity and its external environment. Based on this analysis, a computer algorithm is developed for determining wall admittances as seen at the edge of the patch in the cavity, the radial admittance matrix for the two-probe feed arrangement, and the input impedance as observed from the coaxial line feeding the cavity. This algorithm is tested against the general-purpose Hewlett-Packard finite-element High Frequency Structure Simulator as well as against measured results. Good agreement in all considered cases is noted.