Adaptive backward Euler time stepping with truncation error control for numerical modelling of unsaturated fluid flow

An automatic time stepping scheme with embedded error control is developed and applied to the moisture‐based Richards equation. The algorithm is based on the first‐order backward Euler scheme, and uses a numerical estimate of the local truncation error and an efficient time step selector to control the temporal accuracy of the integration. Local extrapolation, equivalent to the use of an unconditionally stable Thomas–Gladwell algorithm, achieves second‐order temporal accuracy at minimal additional costs. The time stepping algorithm also provides accurate initial estimates for the iterative non‐linear solver. Numerical tests confirm the ability of the scheme to automatically optimize the time step size to match a user prescribed temporal error tolerance. An important merit of the proposed method is its conceptual and computational simplicity. It can be directly incorporated into existing or new software based on the backward Euler scheme (currently prevalent in subsurface hydrologic modelling), and markedly improves their performance compared with simple fixed or heuristic time step selection. The generality of the approach also makes possible its use for solving PDEs in other engineering applications, where strong non‐linearity, stability or implementation considerations favour a simple and robust low‐order method, or where there is a legacy of backward Euler codes in current use. Copyright © 2001 John Wiley & Sons, Ltd.     

Time integration for spherical acoustic finite–infinite element models

In this paper, we analyse the numerical time integration of models of exterior acoustics. The major challenge lies in the instabilities that may arise from the infinite elements. In this paper we consider the special case of spherical infinite elements formulations, which have shown their relevance for industrial applications. We propose a method that combines Crank–Nicholson's method with a filtering step by the backward Euler method. The paper is illustrated with an example relevant to industry. Copyright © 2005 John Wiley & Sons, Ltd.     

Step size adjustment and extrapolation for time‐stepping schemes in non‐smooth dynamics

In this paper we use step size adjustment and extrapolation methods to improve Moreau's time‐stepping scheme for the numerical integration of non‐smooth mechanical systems, i.e. systems with impact and friction. The scheme yields a system of inclusions, which is transformed into a system of projective equations. These equations are solved iteratively. Switching points are time instants for which the structure of the mechanical system changes, for example, time instants for which a sticking friction element begins to slide. We show how switching points can be localized and how these points can be resolved by choosing a minimal step size. In order to improve the integration of non‐smooth systems in the smooth parts, we show how the time‐stepping method can be used as a base integration scheme for extrapolation methods, which allow for an increase in the integration order. Switching points are processed by a small time step, while time intervals during which the structure of the system does not change are computed with a larger step size and improved integration order. The overall algorithm, which consists of a time‐stepping module, an extrapolation module and a step size adjustment module, is discussed in detail and some examples are given. Copyright © 2008 John Wiley & Sons, Ltd.     

Stabilized finite element methods for vertically averaged multiphase flow for carbon sequestration

A computationally efficient numerical model that describes carbon sequestration in deep saline aquifers is presented. The model is based on the multiphase flow and vertically averaged mass balance equations, requiring the solution of two partial differential equations – a pressure equation and a saturation equation. The saturation equation is a nonlinear advective equation for which the application of Galerkin finite element method (FEM) can lead to non‐physical oscillations in the solution. In this article, we extend three stabilized FEM formulations, which were developed for uncoupled systems, to the governing nonlinear coupled PDEs. The methods developed are based on the streamline upwind, the streamline upwind/Petrov–Galerkin and the least squares FEM. Two sequential solution schemes are developed: a single step and a predictor–corrector. The range of Courant numbers yielding smooth and oscillation‐free solutions is investigated for each method. The useful range of Courant numbers found depends upon both the sequential scheme (single step vs predictor–corrector) and also the time integration method used (forward Euler, backward Euler or Crank–Nicolson). For complex problems such as when two plumes meet, only the SU stabilization with an amplified stabilization parameter gives satisfactory results when large time steps are used. Copyright © 2016 John Wiley & Sons, Ltd.     

Compact schemes for anisotropic reaction–diffusion equations with adaptive time step

Many problems in biology and engineering are governed by anisotropic reaction–diffusion equations with a very rapidly varying reaction term. These characteristics of the system imply the use of very fine meshes and small time steps in order to accurately capture the propagating wave avoiding the appearance of spurious oscillations in the wave front. This work develops a fourth‐order compact scheme for anisotropic reaction–diffusion equations with stiff reactive terms. As mentioned, the scheme accounts for the anisotropy of the media and incorporates an adaptive time step for handling the stiff reactive term. The high‐order scheme allows working with coarser meshes without compromising numerical accuracy rendering a more efficient numerical algorithm by reducing the total computation time and memory requirements. The order of convergence of the method has been demonstrated on an analytical solution with Neumann boundary conditions. The scheme has also been implemented for the solution of anisotropic electrophysiology problems. Anisotropic square samples of normal and ischemic cardiac tissue have been simulated by means of the monodomain model with the reactive term defined by Luo–Rudy II dynamics. The simulations proved the effectiveness of the method in handling anisotropic heterogeneous non‐linear reaction–diffusion problems. Bidimensional tests also indicate that the fourth‐order scheme requires meshes about 45% coarser than the standard second‐order method in order to achieve the same accuracy of the results. Copyright © 2009 John Wiley & Sons, Ltd.     

A numerical scheme based on compact integrated-RBFs and Adams–Bashforth/Crank–Nicolson algorithms for diffusion and unsteady fluid flow problems

This paper presents a high-order approximation scheme based on compact integrated radial basis function (CIRBF) stencils and second-order Adams–Bashforth/Crank–Nicolson algorithms for solving time-dependent problems in one and two space dimensions. We employ CIRBF stencils, where the RBF approximations are locally constructed through integration and expressed in terms of nodal values of the function and its derivatives, to discretise the spatial derivatives in the governing equation. We adopt the Adams–Bashforth and Crank–Nicolson algorithms, which are second-order accurate, to discretise the temporal derivatives. The performance of the proposed scheme is investigated numerically through the solution of several test problems, including heat transfer governed by the diffusion equation, shock wave propagation and shock-like solution governed by the Burgers' equation, and torsionally oscillating lid-driven cavity flow governed by the Navier–Stokes equation in the primitive variables. Numerical experiments show that the proposed scheme is stable and high-order accurate in reference to the exact solution of analytic test problems and achieves a good agreement with published results for other test problems.     

An optimal weighted upwinding covolume method on non‐standard grids for convection–diffusion problems in 2D

In this paper we develop an optimal weighted upwinding covolume method on non‐standard covolume grids for convection–diffusion problems in two dimensions. The novel feature of our method is that we construct the non‐standard covolume grid in which the nodes of covolumes vary in the interior of different volumes of primary grid depending on the local weighted factors and further on the local Peclet's numbers. A simple method of finding the local optimal weighted factors is also derived from a non‐linear function of local Peclet's numbers. The developed method leads to a totally new scheme for convection–diffusion problems, which overcomes numerical oscillation, avoids numerical dispersion, and has high‐order accuracy. Some theoretical analyses are given and numerical experiments are presented to illustrate the performance of the method. Copyright © 2006 John Wiley & Sons, Ltd.     

A time‐integration method for the viscoelastic–viscoplastic analyses of polymers and finite element implementation

The present study introduces a time‐integration algorithm for solving a non‐linear viscoelastic–viscoplastic (VE–VP) constitutive equation of isotropic polymers. The material parameters in the constitutive models are stress dependent. The algorithm is derived based on an implicit time‐integration method (Computational Inelasticity. Springer: New York, 1998) within a general displacement‐based finite element (FE) analysis and suitable for small deformation gradient problems. Schapery's integral model is used for the VE responses, while the VP component follows the Perzyna model having an overstress function. A recursive‐iterative method (Int. J. Numer. Meth. Engng 2004; 59 :25–45) is employed and modified to solve the VE–VP constitutive equation. An iterative procedure with predictor–corrector steps is added to the recursive integration method. A residual vector is defined for the incremental total strain and the magnitude of the incremental VP strain. A consistent tangent stiffness matrix, as previously discussed in Ju (J. Eng. Mech. 1990; 116 :1764–1779) and Simo and Hughes (Computational Inelasticity. Springer: New York, 1998), is also formulated to improve convergence and avoid divergence. Available experimental data on time‐dependent and inelastic responses of high‐density polyethylene are used to verify the current numerical algorithm. The time‐integration scheme is examined in terms of its computational efficiency and accuracy. Numerical FE analyses of microstructural responses of polyethylene reinforced with elastic particle are also presented. Copyright © 2009 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.     

Unconditionally energy stable implicit time integration: application to multibody system analysis and design

This paper focuses on the development of an unconditionally stable time‐integration algorithm for multibody dynamics that does not artificially dissipate energy. Unconditional stability is sought to alleviate any stability restrictions on the integration step size, while energy conservation is important for the accuracy of long‐term simulations. In multibody system analysis, the time‐integration scheme is complemented by a choice of co‐ordinates that define the kinematics of the system. As such, the current approach uses a non‐dissipative implicit Newmark method to integrate the equations of motion defined in terms of the independent joint co‐ordinates of the system. In order to extend the unconditional stability of the implicit Newmark method to non‐linear dynamic systems, a discrete energy balance is enforced. This constraint, however, yields spurious oscillations in the computed accelerations and therefore, a new acceleration corrector is developed to eliminate these instabilities and hence retain unconditional stability in an energy sense. An additional benefit of employing the non‐linearly implicit time‐integration method is that it allows for an efficient design sensitivity analysis. In this paper, design sensitivities computed via the direct differentiation method are used for mechanism performance optimization. Copyright © 2000 John Wiley & Sons, Ltd.     

A finite volume method with a modified ENO scheme using a Hermite interpolation to solve advection diffusion equations

In this paper we have developed a finite volume ENO scheme, third‐order accurate, based on cell averages and a TVD Runge–Kutta time discretization to solve advection–diffusion equations in a two‐dimensional spatial domain. We have designed a special interpolating polynomial based on a modified ENO scheme and a Hermite procedure which avoids the excessive smearing in regions with sharpconcentration fronts and the overcompression effects produced by the modified ENO technique. Thesemodifications do not affect the non‐oscillatory philosophy since we compare divided differences inthe modified ENO scheme and in the evaluation of the Hermite polynomial derivatives. Numericalresults compare favourably with their respective analytical solutions. Copyright © 2001 John Wiley & Sons, Ltd.     

Dynamic inflation of non‐linear elastic and viscoelastic rubber‐like membranes

The present paper deals with the dynamic inflation of rubber‐like membranes.The material is assumed to obey the hyperelastic Mooney's model or the non‐linear viscoelastic Christensen's model. The governing equations of free inflation are solved by a total Lagrangian finite element method for the spatial discretization and an explicit finite‐difference algorithm for the time‐integration scheme. The numerical implementation of constitutive equations is highlighted and the special case of integral viscoelastic models is examined in detail. The external force consists in a gas flow rate, which is more realistic than a pressure time history. Then, an original method is used to calculate the pressure evolution inside the bubble depending on the deformation state. Our numerical procedure is illustrated through different examples and compared with both analytical and experimental results. These comparisons yield good agreement and show the ability of our approach to simulate both stable and unstable large strain inflations of rubber‐like membranes. Copyright © 2001 John Wiley & Sons, Ltd.     

A novel numerical method of high‐order accuracy for flow in unsaturated porous media

We construct finite volume schemes of very high order of accuracy in space and time for solving the nonlinear Richards equation (RE). The general scheme is based on a three‐stage predictor–corrector procedure. First, a high‐order weighted essentially non‐oscillatory (WENO) reconstruction procedure is applied to the cell averages at the current time level to guarantee monotonicity in the presence of steep gradients. Second, the temporal evolution of the WENO reconstruction polynomials is computed in a predictor stage by using a global weak form of the governing equations. A global space–time DG FEM is used to obtain a scheme without the parabolic time‐step restriction caused by the presence of the diffusion term in the RE. The resulting nonlinear algebraic system is solved by a Newton–Krylov method, where the generalized minimal residual method algorithm of Saad and Schulz is used to solve the linear subsystems. Finally, as a third step, the cell averages of the finite volume method are updated using a one‐step scheme, on the basis of the solution calculated previously in the space–time predictor stage. Our scheme is validated against analytical, experimental, and other numerical reference solutions in four test cases. A numerical convergence study performed allows us to show that the proposed novel scheme is high order accurate in space and time. Copyright © 2011 John Wiley & Sons, Ltd.     

An artificial compressibility method for viscous incompressible and low Mach number flows

Within the framework of the pressure‐based algorithm, an artificial compressibility method is developed on a non‐orthogonal grid for incompressible and low Mach number fluid flow problems, using cell‐centered finite‐volume approximation. Resorting to the traditional pseudo‐compressibility concept, the continuity constraint is perturbed by the time derivative of pressure, the physical relevance of which is to invoke matrix preconditionings. The approach provokes density perturbations, assisting the transformation between primitive and conservative variables. A dual‐dissipation scheme for the pressure–velocity coupling is contrived, which has the expediences of greater flexibility and increased accuracy in a way similar to the monotone upstream‐centered schemes for conservation laws approach. To account for the flow directionality in the upwinding, a rotational matrix is introduced to evaluate the convective flux. Numerical experiments in reference to a few well‐documented laminar flows demonstrate that the entire contrivance expedites enhanced robustness and improved overall damping properties of the factored pseudo‐time integration procedure. Copyright © 2008 John Wiley & Sons, Ltd.     

Finite element solution of free‐surface ship‐wave problems

An unstructured finite element solver to evaluate the ship‐wave problem is presented. The scheme uses a non‐structured finite element algorithm for the Euler or Navier–Stokes flow as for the free‐surface boundary problem. The incompressible flow equations are solved via a fractional step method whereas the non‐linear free‐surface equation is solved via a reference surface which allows fixed and moving meshes. A new non‐structured stabilized approximation is used to eliminate spurious numerical oscillations of the free surface. Copyright © 1999 John Wiley & Sons, Ltd.     

Hybrid BEM for the early stage of unsteady transport process

This paper presents a new approach for solving the early stage of 3D time‐dependent transport problems in non‐homogeneous fractured porous media in which the initial distribution of concentration presents discontinuous jump at the interface of two regions of transport properties differing in several orders of magnitude. The goal of this formulation is to overcome the problems of different scales and apparent large flux during early time by combining the 3D dual reciprocity boundary element method and a semi‐analytical solution of the time‐dependent advection–diffusion equation, employing a two‐level finite difference time integration scheme. Theoretical background, validation results and practical applications for the advection–diffusion equation in fractured and continuous porous media are reported. The results show advantages for relatively large‐scale models and complex geometries. Copyright © 2006 John Wiley & Sons, Ltd.     

A two‐dimensional stochastic algorithm for the solution of the non‐linear Poisson–Boltzmann equation: validation with finite‐difference benchmarks

This paper presents a two‐dimensional floating random walk (FRW) algorithm for the solution of the non‐linear Poisson–Boltzmann (NPB) equation. In the past, the FRW method has not been applied to the solution of the NPB equation which can be attributed to the absence of analytical expressions for volumetric Green's functions. Previous studies using the FRW method have examined only the linearized Poisson–Boltzmann equation. No such linearization is needed for the present approach. Approximate volumetric Green's functions have been derived with the help of perturbation theory, and these expressions have been incorporated within the FRW framework. A unique advantage of this algorithm is that it requires no discretization of either the volume or the surface of the problem domains. Furthermore, each random walk is independent, so that the computational procedure is highly parallelizable. In our previous work, we have presented preliminary calculations for one‐dimensional and quasi‐one‐dimensional benchmark problems. In this paper, we present the detailed formulation of a two‐dimensional algorithm, along with extensive finite‐difference validation on fully two‐dimensional benchmark problems. The solution of the NPB equation has many interesting applications, including the modelling of plasma discharges, semiconductor device modelling and the modelling of biomolecular structures and dynamics. Copyright © 2005 John Wiley & Sons, Ltd.     

A discontinuous Galerkin method for a hyperbolic model for convection–diffusion problems in CFD

This paper proposes a hyperbolic model for convection–diffusion transport problems in computational fluid dynamics (CFD). The hyperbolic model is based on the so‐called Cattaneo's law. This is a time‐dependent generalization of Fick's and Fourier's laws that was originally proposed to solve pure‐diffusive heat transfer problems. We show that the proposed model avoids the infinite speed paradox that is inherent in the standard parabolic model. A high‐order upwind discontinuous Galerkin (DG) method is developed and applied to classic convection‐dominated test problems. The quality of the numerical results is remarkable, since the discontinuities are very well captured without the appearance of spurious oscillations. These results are compared with those obtained by using the standard parabolic model and the local DG (LDG) method and with those given by the parabolic model and the Bassi–Rebay scheme. Finally, the applicability of the proposed methodology is demonstrated by solving a practical case in engineering. We simulate the evolution of pollutant being spilled in the harbour of A Coruña (northwest of Spain, EU). Copyright © 2007 John Wiley & Sons, Ltd.     

A contact‐stabilized Newmark method for dynamical contact problems

The numerical integration of dynamical contact problems often leads to instabilities at contact boundaries caused by the non‐penetration condition between bodies in contact. Even an energy dissipative modification (see, e.g. (Comp. Meth. Appl. Mech. Eng. 1999; 180 :1–26)), which discretizes the non‐penetration constraints implicitly, is not able to circumvent artificial oscillations. For this reason, the present paper suggests a contact stabilization in function space, which avoids artificial oscillations at contact interfaces and is also energy dissipative. The key idea of this contact stabilization is an additional L²‐projection at contact interfaces, which can be easily added to any existing time integration scheme. In case of a lumped mass matrix, this projection can be carried out completely locally, thus creating only negligible additional numerical cost. For the new scheme, an elementary analysis is given, which is confirmed by numerical findings in an illustrative test example (Hertzian two‐body contact). Copyright © 2007 John Wiley & Sons, Ltd.     

A robust and efficient substepping scheme for the explicit numerical integration of a rate‐dependent crystal plasticity model

This paper describes the development of efficient and robust numerical integration schemes for rate‐dependent crystal plasticity models. A forward Euler integration algorithm is first formulated. An integration algorithm based on the modified Euler method with an adaptive substepping scheme is then proposed, where the substepping is mainly controlled by the local error of the stress predictions within the time step. Both integration algorithms are implemented in a stand‐alone code with the Taylor aggregate assumption and in an explicit finite element code. The robustness, accuracy and efficiency of the substepping scheme are extensively evaluated for large time steps, extremely low strain‐rate sensitivity, high deformation rates and strain‐path changes using the stand‐alone code. The results show that the substepping scheme is robust and in some cases one order of magnitude faster than the forward Euler algorithm. The use of mass scaling to reduce computation time in crystal plasticity finite element simulations for quasi‐static problems is also discussed. Finally, simulation of Taylor bar impact test is carried out to show the applicability and robustness of the proposed integration algorithm for the modelling of dynamic problems with contact. Copyright © 2014 John Wiley & Sons, Ltd.