Numerical modelling of fluids mixing,heat transfer and non‐equilibrium redox chemical reactions in fluid‐saturated porous rocks

Non‐equilibrium redox chemical reactions of high orders are ubiquitous in fluid‐saturated porous rocks within the crust of the Earth. The numerical modelling of such high‐order chemical reactions becomes a challenging problem because these chemical reactions are not only produced strong non‐linear source/sink terms for reactive transport equations, but also often coupled with the fluids mixing, heat transfer and reactive mass transport processes. In order to solve this problem effectively and efficiently, it is desirable to reduce the total number of reactive transport equations with strong non‐linear source/sink terms to a minimum in a computational model. For this purpose, the concept of the chemical reaction rate invariant is used to develop a numerical procedure for dealing with fluids mixing, heat transfer and non‐equilibrium redox chemical reactions in fluid‐saturated porous rocks. Using the proposed concept and numerical procedure, only one reactive transport equation, which is used to describe the distribution of the chemical product and has a strong non‐linear source/sink term, needs to be solved for each of the non‐equilibrium redox chemical reactions. The original reactive transport equations of the chemical reactants with strong non‐linear source/sink terms are turned into the conventional mass transport equations of the chemical reaction rate invariants without any non‐linear source/sink terms. A testing example, for some aspects of which the analytical solutions are available, is used to validate the proposed numerical procedure. The related numerical solutions have demonstrated that (1) the proposed numerical procedure is useful and applicable for dealing with the coupled problem between fluids mixing, heat transfer and non‐equilibrium redox chemical reactions of high orders in fluid‐saturated porous rocks; (2) the interaction between the solute diffusion, solute advection and chemical kinetics is an important mechanism to control distribution patterns of chemical products in an ore‐forming process; and (3) if the pore‐fluid pressure gradient is lithostatic, it is difficult for the chemical equilibrium to be attained within permeable cracks and geological faults within the crust of the Earth. Copyright © 2005 John Wiley & Sons, Ltd.     

Weak forms for modelling of rotationally symmetric,multilayered structures,including anisotropic poro‐elastic media

A weak form of the anisotropic Biot's equation represented in a cylindrical coordinate system using a spatial Fourier expansion in the circumferential direction is presented. The original three dimensional Cartesian anisotropic weak formulation is rewritten in an arbitrary orthogonal curvilinear basis. Introducing a cylindrical coordinate system and expanding the circumferential wave propagation in terms of orthogonal harmonic functions, the original, geometrically rotationally symmetric three dimensional boundary value problem, is decomposed into independent two‐dimensional problems, one for each harmonic function. Using a minimum number of dependent variables, pore pressure and frame displacement, a computationally efficient procedure for vibro‐acoustic finite element modelling of rotationally symmetric three‐dimensional multilayered structures including anisotropic porous elastic materials is thus obtained. By numerical simulations, this method is compared with, and the correctness is verified against, a full three‐dimensional Cartesian coordinate system finite element model. Copyright © 2012 John Wiley & Sons, Ltd.     

Stable element‐free Galerkin solution procedures for the coupled soil–pore fluid problem

It is always difficult to solve the coupled soil–pore fluid problem when the soil is saturated and impermeable, because this situation often results in intensive oscillations of the solutions. This topic has been discussed widely in the field of the finite element method but rarely by meshless methods. The Element‐Free Galerkin (EFG) method has outstanding advantages of solving this problem, based on the fact that its interpolation function can be constructed flexibly by the nodes in the compact domain. In this study, an EFG numerical model together with a stabilization technique is proposed to obtain stable solutions, especially for the saturated and impermeable soil. Close agreement between computational results and analytical solutions for one‐and two‐dimensional examples shows that the proposed numerical model not only provides highly accurate solutions for the saturated soil with relatively high permeability, but also eliminates the oscillations of the solutions very effectively for the saturated and impermeable soil. Furthermore, the influences of the hydraulic anisotropy, a typical property of two‐dimensional problems, on the proposed stabilization technique are discussed. It is suggested that the optimal distribution of the nodes of essential variables can be designed, according to the relative importance of the permeability in different directions. Copyright © 2010 John Wiley & Sons, Ltd.     

Three‐dimensional coupled heat,moisture, and air transfer in a deformable unsaturated soil

A three‐dimensional numerical model is presented for three‐phase flow (moisture, air, and heat) in a deformable partly saturated soil with deformation calculated via a non‐linear elastic theory. The present work is an extension of a two‐dimensional analysis presented by Thomas and He. The objective of this work is the solution of problems of greater geometric complexity. The mathematical formulation of this coupled problem consists of four governing equations, developed from the principles of mass and energy conservations as well as the stress equilibrium equation. Darcy's flow law is used to describe the motion of liquid and air in the porous medium, and a Philip and de Vries type vapour flow approach is employed in the formulation. A Galerkin finite element method coupled with a finite difference recurrence relationship is used to obtain simultaneous solutions to the governing equations where pore liquid, pore air pressures, temperature and displacements are the primary variables. The method allows the non‐linear nature of the soil parameters to be modelled. Three‐dimensional 20‐noded isoparametric elements are used to simulate different types of cases for the verification of the work. Results are presented of the application of the new model to four problems, two of which are isothermal and two heating simulations. The three‐dimensional nature of the results achieved is highlighted. Copyright © 1999 John Wiley & Sons, Ltd.     

A stabilized finite element formulation for monolithic thermo‐hydro‐mechanical simulations at finite strain

An adaptively stabilized monolithic finite element model is proposed to simulate the fully coupled thermo‐hydro‐mechanical behavior of porous media undergoing large deformation. We first formulate a finite‐deformation thermo‐hydro‐mechanics field theory for non‐isothermal porous media. Projection‐based stabilization procedure is derived to eliminate spurious pore pressure and temperature modes due to the lack of the two‐fold inf‐sup condition of the equal‐order finite element. To avoid volumetric locking due to the incompressibility of solid skeleton, we introduce a modified assumed deformation gradient in the formulation for non‐isothermal porous solids. Finally, numerical examples are given to demonstrate the versatility and efficiency of this thermo‐hydro‐mechanical model. Copyright © 2015 John Wiley & Sons, Ltd.     

A coupling extended multiscale finite element method for dynamic analysis of heterogeneous saturated porous media

A coupling extended multiscale finite element method (CEMsFEM) is developed for the dynamic analysis of heterogeneous saturated porous media. The coupling numerical base functions are constructed by a unified method with an equivalent stiffness matrix. To improve the computational accuracy, an additional coupling term that could reflect the interaction of the deformations among different directions is introduced into the numerical base functions. In addition, a kind of multi‐node coarse element is adopted to describe the complex high‐order deformation on the boundary of the coarse element for the two‐dimensional dynamic problem. The coarse element tests show that the coupling numerical base functions could not only take account of the interaction of the solid skeleton and the pore fluid but also consider the effect of the inertial force in the dynamic problems. On the other hand, based on the static balance condition of the coarse element, an improved downscaling technique is proposed to directly obtain the satisfying microscopic solutions in the CEMsFEM. Both one‐dimensional and two‐dimensional numerical examples of the heterogeneous saturated porous media are carried out, and the results verify the validity and the efficiency of the CEMsFEM by comparing with the conventional finite element method. Copyright © 2015 John Wiley & Sons, Ltd.     

Robin‐Neumann transmission conditions for fluid‐structure coupling: Embedded boundary implementation and parameter analysis

Partitioned procedures are appealing for solving complex fluid‐structure interaction (FSI) problems, as they allow existing computational fluid dynamics (CFD) and computational structural dynamics algorithms and solvers to be combined and reused. However, for problems involving incompressible flow and strong added‐mass effect (eg, heavy fluid and slender structure), partitioned procedures suffer from numerical instability, which typically requires additional subiterations between the fluid and structural solvers, hence significantly increasing the computational cost. This paper investigates the use of Robin‐Neumann transmission conditions to mitigate the above instability issue. Firstly, an embedded Robin boundary method is presented in the context of projection‐based incompressible CFD and finite element–based computational structural dynamics. The method utilizes operator splitting and a modified ghost fluid method to enforce the Robin transmission condition on fluid‐structure interfaces embedded in structured non–body‐conforming CFD grids. The method is demonstrated and verified using the Turek and Hron benchmark problem, which involves a slender beam undergoing large transient deformation in an unsteady vortex‐dominated channel flow. Secondly, this paper investigates the effect of the combination parameter in the Robin transmission condition, ie, α_f, on numerical stability and solution accuracy. This paper presents a numerical study using the Turek and Hron benchmark problem and an analytical study using a simplified FSI model featuring an Euler‐Bernoulli beam interacting with a two‐dimensional incompressible inviscid flow. Both studies reveal a trade‐off between stability and accuracy: smaller values of α_f tend to improve numerical stability, yet deteriorate the accuracy of the partitioned solution. Using the simplified FSI model, the critical value of α_f that optimizes this trade‐off is derived and discussed.     

Stabilized low‐order finite elements for strongly coupled poromechanical problems

Finite element solutions of poromechanical problems often exhibit oscillating pore pressures in the limits of low permeability, fast loading rates, coarse meshes, and/or small time step sizes. To suppress completely the pore pressure oscillations, a stabilized finite element scheme with a better performance on monotonicity is proposed for modeling compressible fluid‐saturated porous media. This method, based on the polynomial pressure projection technique, allows the use of linear equal‐order interpolation for both displacement and pore pressure fields, which is more straightforward for both code development and maintenance compared to others. By employing the discrete maximum principle, a proper stabilization parameter is deduced, which is efficient to guarantee the monotonicity and optimal in theory in the 1‐dimensional case. An appealing feature of the method is that the stabilization parameter is evaluated in terms of the properties of porous material only, while no mesh or time step size is involved. Through comparing the numerical simulations with the analytical benchmarks, the efficiency of the proposed stabilization scheme is confirmed.     

A posteriori finite element bounds to linear functional outputs of the three‐dimensional Navier–Stokes equations

An implicit a posteriori finite element error estimation method is presented to inexpensively calculate lower and upper bounds for a linear functional output of the numerical solutions to the three‐dimensional Navier–Stokes (N–S) equations. The novelty of this research is to utilize an augmented Lagrangian based on a coarse mesh linearization of the N–S equations and the finite element tearing and interconnecting (FETI) procedure. The latter approach extends the a posteriori bound method to the three‐dimensional Crouzeix–Raviart space for N–S problems. The computational advantage of the bound procedure is that a single coupled non‐symmetric large problem can be decomposed into several uncoupled symmetric small problems. A simple model problem, which is selected here to illustrate the procedure, is to find the lower and upper bounds of the average velocity of a pressure driven, incompressible, steady Newtonian fluid flow moving at low Reynolds numbers through an endless square channel which has an array of rectangular obstacles. Numerical results show that the bounds for this output are rigorous, i.e. always in the asymptotic certainty regime, that they are sharp and that the required computational resources decrease significantly. Parallel implementation on a Beowulf cluster is also reported. Copyright © 2004 John Wiley & Sons, Ltd.     

A fully coupled finite element analysis of water‐table fluctuation and land deformation in partially saturated soils due to surface loading

A fully coupled numerical model is presented for the water‐table fluctuation and land deformation in partially saturated soils due to surface loading. This numerical model is developed based on the poroelastic governing equations for groundwater flow in deforming variably saturated porous media and the Galerkin finite element method. The numerical model is verified and validated against a one‐dimensional consolidation problem concerning surface loading on a soil column which has six different initial water‐table elevations. The numerical model is then applied to a two‐dimensional consolidation problem of surface loading on a partially saturated soil at a construction site. Results from the numerical simulations of both problems show that the water table fluctuates in the partially saturated soils, and the unsaturated zone above the water table has significant effects on the consolidation behaviour of the partially saturated soils under surface loading. Such effects are caused by the permanent absorption of a portion of the mechanical loading stress and the weak hydromechanical coupling between the solid skeleton deformation field and the groundwater flow field in the unsaturated zone due to its partial saturation. Copyright © 2000 John Wiley & Sons, Ltd.     

Improvements of explicit crack surface representation and update within the generalized finite element method with application to three‐dimensional crack coalescence

This paper presents improvements to three‐dimensional crack propagation simulation capabilities of the generalized finite element method. In particular, it presents new update algorithms suitable for explicit crack surface representations and simulations in which the initial crack surfaces grow significantly in size (one order of magnitude or more). These simulations pose problems in regard to robust crack surface/front representation throughout the propagation analysis. The proposed techniques are appropriate for propagation of highly non‐convex crack fronts and simulations involving significantly different crack front speeds. Furthermore, the algorithms are able to handle computational difficulties arising from the coalescence of non‐planar crack surfaces and their interactions with domain boundaries. An approach based on moving least squares approximations is developed to handle highly non‐convex crack fronts after crack surface coalescence. Several numerical examples are provided, which illustrate the robustness and capabilities of the proposed approaches and some of its potential engineering applications. Copyright © 2013 John Wiley & Sons, Ltd.     

Finite element modelling of saturated porous media at finite strains under dynamic conditions with compressible constituents

Two finite element formulations are proposed to analyse the dynamic conditions of saturated porous media at large strains with compressible solid and fluid constituents. Unlike similar works published in the literature, the proposed formulations are based on a recently proposed hyperelastic framework in which the compressibility of the solid and fluid constituents is fully taken into account when geometrical non‐linear effects are relevant on both micro‐ and macroscales. The first formulation leads to a three‐field finite element method (FEM), which is suitable for analysing high‐frequency dynamic problems, whereas the second is a simplification of the first, leading to a two‐field FEM, in which some inertial effects of the pore fluid are disregarded, hence the second formulation is suitable for studying low‐frequency problems. A fully Lagrangian approach is considered, hence all terms are expressed with reference to the material setting; the balance equations for the pore fluid are also expressed in terms of the chemical potential and the mass flux of the pore fluid in order to take the compressibility of the fluid into account. To improve the numerical response in the case of wave propagation, a discontinuous Galerkin FEM in the time domain is applied to the three‐field formulation. The results are compared with analytical and semi‐analytical solutions, highlighting the different effects of the discontinuous Galerkin method on the longitudinal waves of the first and second kind. Copyright © 2010 John Wiley & Sons, Ltd.     

A coupled discrete element‐finite difference approach for modeling mechanical response of fluid‐saturated porous materials

A model of fluid‐saturated poroelastic medium was developed based on a combination of the discrete element method and grid method. The developed model adequately accounts for the deformation, fracture, and multiscale internal structure of a porous solid skeleton. The multiscale porous structure is taken into account implicitly by assigning the porosity and permeability values for the enclosing skeleton, which determine the rate of filtration of a fluid. Macroscopic pores and voids are taken into account explicitly by specifying the computational domain geometry. The relationship between the stress–strain state of the solid skeleton and pore fluid pressure is described in the approximations of simply deformable discrete element and Biot's model of poroelasticity. The developed model was applied to study the mechanical response of fluid‐saturated samples of brittle material. Based on simulation results, we constructed a generalized logistic dependence of uniaxial compressive strength on loading rate, mechanical properties of fluid and enclosing skeleton, and on sample dimensions. The logistic form of the generalized dependence of strength of fluid‐saturated elastic–brittle porous materials is due to the competition of two interrelated processes, such as pore fluid pressure increase under solid skeleton compression and fluid outflow from the enclosing skeleton to the environment. Copyright © 2015 John Wiley & Sons, Ltd.     

An implicit–explicit solution method for electro‐osmotic flow through three‐dimensional micro‐channels

This paper presents a comprehensive finite‐element modelling approach to electro‐osmotic flows on unstructured meshes. The non‐linear equation governing the electric potential is solved using an iterative algorithm. The employed algorithm is based on a preconditioned GMRES scheme. The linear Laplace equation governing the external electric potential is solved using a standard pre‐conditioned conjugate gradient solver. The coupled fluid dynamics equations are solved using a fractional step‐based, fully explicit, artificial compressibility scheme. This combination of an implicit approach to the electric potential equations and an explicit discretization to the Navier–Stokes equations is one of the best ways of solving the coupled equations in a memory‐efficient manner. The local time‐stepping approach used in the solution of the fluid flow equations accelerates the solution to a steady state faster than by using a global time‐stepping approach. The fully explicit form and the fractional stages of the fluid dynamics equations make the system memory efficient and free of pressure instability. In addition to these advantages, the proposed method is suitable for use on both structured and unstructured meshes with a highly non‐uniform distribution of element sizes. The accuracy of the proposed procedure is demonstrated by solving a basic micro‐channel flow problem and comparing the results against an analytical solution. The comparisons show excellent agreement between the numerical and analytical data. In addition to the benchmark solution, we have also presented results for flow through a fully three‐dimensional rectangular channel to further demonstrate the application of the presented method. Copyright © 2007 John Wiley & Sons, Ltd.     

A Galerkin/least‐squares finite element formulation for consolidation

A Galerkin/least‐squares (GLS) finite element formulation for problem of consolidation of fully saturated two‐phase media is presented. The elimination of spurious pressure oscillations appearing at the early stage of consolidation for standard Galerkin finite elements with equal interpolation order for both displacements and pressures is the goal of the approach. It will be shown that the least‐squares term, based exclusively on the residuum of the fluid flow continuity equation, added to the standard Galerkin formulation enhances its stability and can fully eliminate pressure oscillations. A reasonably simple framework designed for derivation of one‐dimensional as well as multi‐dimensional estimates of the stabilization factor is proposed and then verified. The formulation is validated on one‐dimensional and then on two‐dimensional, linear and non‐linear test problems. The effect of the fluid incompressibility as well as compressibility will be taken into account and investigated. Copyright © 2001 John Wiley & Sons Ltd.     

Weak formulation of Biot's equations in cylindrical coordinates with harmonic expansion in the circumferential direction

A weak symmetric form of Biot's equation in cylindrical coordinates with a spatial Fourier expansion in the circumferential direction is presented. The solid phase displacement and the pore pressure are used as the dependent variables. The original three‐dimensional boundary value problem is here, due to the orthogonality of the harmonic functions and the rotationally symmetric geometry, decomposed into independent two‐dimensional problems, one for each harmonic function. This formulation provides a computationally efficient procedure for vibroacoustic finite element modelling of rotationally symmetric three‐dimensional multilayered structures including porous elastic materials. By numerical simulations, this method is compared with, and verified against, full three‐dimensional Cartesian coordinate system finite element models. Copyright © 2009 John Wiley & Sons, Ltd.     

A two‐scale approach for fluid flow in fractured porous media

A two‐scale numerical model is developed for fluid flow in fractured, deforming porous media. At the microscale the flow in the cavity of a fracture is modelled as a viscous fluid. From the micromechanics of the flow in the cavity, coupling equations are derived for the momentum and the mass couplings to the equations for a fluid‐saturated porous medium, which are assumed to hold on the macroscopic scale. The finite element equations are derived for this two‐scale approach and integrated over time. By exploiting the partition‐of‐unity property of the finite element shape functions, the position and direction of the fractures is independent from the underlying discretization. The resulting discrete equations are non‐linear due to the non‐linearity of the coupling terms. A consistent linearization is given for use within a Newton–Raphson iterative procedure. Finally, examples are given to show the versatility and the efficiency of the approach, and show that faults in a deforming porous medium can have a significant effect on the local as well as on the overall flow and deformation patterns. Copyright © 2006 John Wiley & Sons, Ltd.     

A discontinuous Galerkin finite element method for dynamic and wave propagation problems in non‐linear solids and saturated porous media

A time‐discontinuous Galerkin finite element method (DGFEM) for dynamics and wave propagation in non‐linear solids and saturated porous media is presented. The main distinct characteristic of the proposed DGFEM is that the specific P3–P1 interpolation approximation, which uses piecewise cubic (Hermite's polynomial) and linear interpolations for both displacements and velocities, in the time domain is particularly proposed. Consequently, continuity of the displacement vector at each discrete time instant is exactly ensured, whereas discontinuity of the velocity vector at the discrete time levels still remains. The computational cost is then obviously saved, particularly in the materially non‐linear problems, as compared with that required for the existing DGFEM. Both the implicit and explicit algorithms are developed to solve the derived formulations for linear and materially non‐linear problems. Numerical results illustrate good performance of the present method in eliminating spurious numerical oscillations and in providing much more accurate solutions over the traditional Galerkin finite element method using the Newmark algorithm in the time domain. Copyright © 2003 John Wiley & Sons, Ltd.     

Non‐planar mixed‐mode growth of initially straight‐fronted surface cracks,in cylindrical bars under tension,torsion and bending,using the symmetric Galerkin boundary element method‐finite element method alternating method

In this paper, the stress intensity factor (SIF) variations along an arbitrarily developing crack front, the non‐planar fatigue‐crack growth patterns, and the fatigue life of a round bar with an initially straight‐fronted surface crack, are studied by employing the 3D symmetric Galerkin boundary element method‐finite element method (SGBEM‐FEM) alternating method. Different loading cases, involving tension, bending and torsion of the bar, with different initial crack depths and different stress ratios in fatigue, are considered. By using the SGBEM‐FEM alternating method, the SIF variations along the evolving crack front are computed; the fatigue growth rates and directions of the non‐planar growths of the crack surface are predicted; the evolving fatigue‐crack growth patterns are simulated, and thus, the fatigue life estimations of the cracked round bar are made. The accuracy and reliability of the SGBEM‐FEM alternating method are verified by comparing the presently computed results to the empirical solutions of SIFs, as well as experimental data of fatigue crack growth, available in the open literature. It is shown that the current approach gives very accurate solutions of SIFs and simulations of fatigue crack growth during the entire crack propagation, with very little computational burden and human–labour cost. The characteristics of fatigue growth patterns of initially simple‐shaped cracks in the cylindrical bar under different Modes I, III and mixed‐mode types of loads are also discussed in detail.     

A semi‐Lagrangean time‐integration approach for extended finite element methods

Many computational problems incorporate discontinuities that evolve in time. The eXtendend Finite Element Method (XFEM) is able to represent discontinuities sharply on fixed arbitrary meshes, but numerical difficulties arise if these discontinuities move in time. We point out that this issue is crucial for interface problems with strongly discontinuous fields on fixed grids. A method using semi‐Lagrangean techniques is proposed to adequately handle time integration based on finite difference schemes in the context of the XFEM. The basic idea is to adapt previous numerical solutions to the current interface position by tracking back virtual Lagrangean particles to their previous positions, where an appropriate solution can be extrapolated from a smooth field. Convergence properties of the proposed method in time and space are thoroughly studied for two one‐dimensional model problems. Finally, the method is applied to the particularly challenging problem of premixed combustion, where the discontinuity appears at the flame front separating the burnt from the unburnt gases. A two‐dimensional and a three‐dimensional expanding flame demonstrates that the method is sufficiently accurate to retain the properties of the overall Nitsche‐type formulation for interface problems with embedded strong discontinuities. Copyright © 2014 John Wiley & Sons, Ltd.