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.     

High‐order predictor–corrector algorithms

New predictor–corrector algorithms are presented for the computation of solution paths of non‐linear partial differential equations. The predictors and the correctors are based on perturbation techniques and Padé approximants. This extends the Asymptotic Numerical Method (ANM), which is an efficient high‐order continuation technique without corrector. The efficiency and the reliability of the new technique are assessed by several examples within thin shell theory and Navier–Stokes equations. Many variants have been tested to establish an optimal algorithm. Copyright © 2002 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.     

Remarks on the interpretation of current non‐linear finite element analyses as differential–algebraic equations

For the numerical solution of materially non‐linear problems like in computational plasticity or viscoplasticity the finite element discretization in space is usually coupled with point‐wise defined evolution equations characterizing the material behaviour. The interpretation of such systems as differential–algebraic equations (DAE) allows modern‐day integration algorithms from Numerical Mathematics to be efficiently applied. Especially, the application of diagonally implicit Runge–Kutta methods (DIRK) together with a Multilevel‐Newton method preserves the algorithmic structure of current finite element implementations which are based on the principle of virtual displacements and on backward Euler schemes for the local time integration. Moreover, the notion of the consistent tangent operator becomes more obvious in this context. The quadratical order of convergence of the Multilevel‐Newton algorithm is usually validated by numerical studies. However, an analytical proof of this second order convergence has already been given by authors in the field of non‐linear electrical networks. We show that this proof can be applied in the current context based on the DAE interpretation mentioned above. We finally compare the proposed procedure to several well‐known stress algorithms and show that the inclusion of a step‐size control based on local error estimations merely requires a small extra time‐investment. Copyright © 2001 John Wiley & Sons, Ltd.     

A new element‐by‐element method for trajectory calculations with tetrahedral finite element meshes

A new method to track massless particles in a three‐dimensional flow field is presented. The method is based on an element‐by‐element approach coupled with a predictor–corrector shooting scheme and does not use any time step. By analogy with time‐dependent schemes, the number of shootings is related to an equivalent number of time steps. The method has been implemented in a finite element framework using unstructured tetrahedral finite element meshes. However, it is general enough so that it can be implemented in finite difference and finite volume frameworks as well. It has been tested on a variety of flow systems namely: the Poiseuille flow in an empty circular pipe, the rotating flow in a stirred tank, the shear flow in a square tank and the flow through a static mixer. Accuracy has been found to depend on the accuracy of the velocity computation, the number of points per element and the level of mesh refinement. Copyright © 2006 John Wiley & Sons, Ltd.     

Boltzmann schemes for continuum gas dynamics

Many problems arising in the aerodynamic design of aerospace vehicles require the numerical solution of the Euler equations of gas dynamics. These are nonlinear partial differential equations admitting weak solutions such as shock waves and constructing robust numerical schemes for these equations is a challenging task. A new line of research called Boltzmann or kinetic schemes discussed in the present paper exploits the connection between the Boltzmann equation of the kinetic theory of gases and the Euler equations for inviscid compressible flows. Because of this connection, a suitable moment of a numerical scheme for the Boltzmann equation yields a numerical scheme for the Euler equations. This idea called the “moment method strategy” turns out to be an extremely rich methodology for developing robust numerical schemes for the Euler equations. The richness is demonstrated by developing a variety of kinetic schemes such as kinetic numerical method, kinetic flux vector splitting method, thermal velocity based splitting, multidirectional upwind method and least squares weak upwind scheme. A 3-D time-marching Euler code calledbheema based on the kinetic flux vector splitting method and its variants involving equilibrium chemistry have been developed for computing hypersonic reentry flows. The results obtained from the codebheema demonstrate the robustness and the utility of the kinetic flux vector splitting method as a design tool in aerodynamics. The work presented in this paper is based on the research work done by several graduate students at our laboratory and collaborators from research and development organizations within the country.     

A novel ‘optimal’ exponential‐based integration algorithm for von‐Mises plasticity with linear hardening: Theoretical analysis on yield consistency,accuracy, convergence and numerical investigations

In this communication we propose a new exponential‐based integration algorithm for associative von‐Mises plasticity with linear isotropic and kinematic hardening, which follows the ones presented by the authors in previous papers. In the first part of the work we develop a theoretical analysis on the numerical properties of the developed exponential‐based schemes and, in particular, we address the yield consistency, exactness under proportional loading, accuracy and stability of the methods. In the second part of the contribution, we show a detailed numerical comparison between the new exponential‐based method and two classical radial return map methods, based on backward Euler and midpoint integration rules, respectively. The developed tests include pointwise stress–strain loading histories, iso‐error maps and global boundary value problems. The theoretical and numerical results reveal the optimal properties of the proposed scheme. Copyright © 2006 John Wiley & Sons, Ltd.     

Application of meshfree collocation method to a class of nonlinear partial differential equations

In this work, an algorithm for the numerical solution of the generalized Hirota–Satsuma equations and Jaulent–Miodek equations based on meshless radial basis functions (RBFs) method using collocation points, called Kansa's method, is presented. Four model problems with six different initial conditions are considered for the computation. A fairly explicit scheme is used to approximate the solution. The comparison is made with the exact solutions of each problem of the generalized Hirota–Satsuma coupled Korteweg–de Vries equations. A system consisting highly nonlinear partial differential equations known as Jaulent–Miodek equations and generalized Hirota–Satsuma coupled modified-Korteweg–de Vries equations are considered for comparison with the work already published. The multiquadric RBF results are compared with homotopy perturbation method (HPM) and variational iteration method (VIM) to highlight the excellent performance of the method.     

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.     

Fourth‐order energy‐preserving locally implicit time discretization for linear wave equations

A family of fourth‐order coupled implicit–explicit time schemes is presented as a special case of fourth‐order coupled implicit schemes for linear wave equations. The domain of interest is decomposed into several regions where different fourth‐order time discretizations are used, chosen among a family of implicit or explicit fourth‐order schemes. The coupling is based on a Lagrangian formulation on the boundaries between the several non‐conforming meshes of the regions. A global discrete energy is shown to be preserved and leads to global fourth‐order consistency in time. Numerical results in 1D and 2D for the acoustic and elastodynamics equations illustrate the good behavior of the schemes and their potential for the simulation of realistic highly heterogeneous media or strongly refined geometries, for which using everywhere an explicit scheme can be extremely penalizing. Accuracy up to fourth order reduces the numerical dispersion inherent to implicit methods used with a large time step and makes this family of schemes attractive compared with second‐order accurate methods. Copyright © 2015 John Wiley & Sons, Ltd.     

Integrating a micromechanical model for multiscale analyses

The Chang‐Hicher micromechanical model based on a static hypothesis, not unlike other models developed separately at around the same epoch, has proved its efficiency in predicting soil behaviour. For solving boundary value problems, the model has now integrated stress‐strain relationships by considering both the micro and macro levels. The first step was to solve the linearized mixed control constraints by the introduction of a predictor–corrector scheme and then to implement the micro–macro relationships through an iterative procedure. Two return mapping schemes, consisting of the closest‐point projection method and the cutting plane algorithm, were subsequently integrated into the interparticle force‐displacement relations. Both algorithms have proved to be efficient in studies devoted to elementary tests and boundary value problems. Closest‐point projection method compared with cutting plane algorithm, however, has the advantage of being more intensive cost efficient and just as accurate in the computational task of integrating the local laws into the micromechanical model. The results obtained demonstrate that the proposed linearized method is capable of performing loadings under stress and strain control. Finally, by applying a finite element analysis with a biaxial test and a square footing, it can be recognized that the Chang‐Hicher micromechanical model performs efficiently in multiscale modelling.     

An optimized predictor–corrector scheme for fast 3D crack growth simulations

An optimized predictor–corrector scheme for the accelerated simulation of 3D fatigue crack growth is presented. Based on experimental evidence, it is assumed that the crack front shape ensures a constant energy release rate. Starting from a crack front satisfying this requirement a predictor step is performed. Usually, the new crack front does not fulfill the requirement of a constant energy release rate. Therefore, several corrector steps are needed. Within the new predictor–corrector scheme the history of crack growth is taken into account to reduce the number of corrector steps. The efficiency of the new scheme is shown on two numerical examples providing a speed up of a factor above three.     

A general conjugate‐gradient‐based predictor–corrector solver for explicit finite‐element contact

A Lagrange‐multiplier based approach is presented for the general solution of multi‐body contact within an explicit finite element framework. The technique employs an explicit predictor step to permit the detection of interpenetration and then utilizes a corrector step, whose solution is obtained with a pre‐conditioned matrix‐free conjugate gradient projection method, to determine the Lagrange multipliers necessary to eliminate the predicted penetration. The predictor–corrector algorithm is developed for deformable bodies based upon the central difference method, and for rigid bodies from momentum and energy conserving approaches. Both frictionless and Coulomb‐based frictional contact idealizations are addressed. The technique imposes no time‐step constraints and quickly mitigates velocity discontinuities across closed interfaces. Special attention is directed toward contact between rigid bodies. Algorithmic moment arms conserve the translational and angular momentums of the system in the absence of external loads. Elastic collisions are captured with a two‐phase predictor–corrector approach and a geometrically approximate velocity jump criterion. The first step solves the inelastic contact problem and identifies inactive constraints between rigid bodies, while the second step generates the necessary velocity jump condition on the active constraints. The velocity criterion is shown to algorithmically preserve the system kinetic energy for two unconstrained rigid bodies. Copyright © 1999 John Wiley & Sons, Ltd. This paper was produced under the auspices of the U.S. Government and it is therefore not subject to copyright in the U.S.     

Numerical implementation of a shape memory alloy thermomechanical constitutive model using return mapping algorithms

Previous studies by the authors and their co‐workers show that the structure of equations representing shape Memory Alloy (SMA) constitutive behaviour can be very similar to those of rate‐independent plasticity models. For example, the Boyd–Lagoudas polynomial hardening model has a stress‐elastic strain constitutive relation that includes the transformation strain as an internal state variable, a transformation function determining the onset of phase transformation, and an evolution equation for the transformation strain. Such a structure allows techniques used in rate‐independent elastoplastic behaviour to be directly applicable to SMAs. In this paper, a comprehensive study on the numerical implementation of SMA thermomechanical constitutive response using return mapping (elastic predictor‐transformation corrector) algorithms is presented. The closest point projection return mapping algorithm which is an implicit scheme is given special attention together with the convex cutting plane return mapping algorithm, an explicit scheme already presented in an earlier work. The closest point algorithm involves relatively large number of tensorial operations than the cutting plane algorithm besides the evaluation of the gradient of the transformation tensor in the flow rule and the inversion of the algorithmic tangent tensor. A unified thermomechanical constitutive model, which does not take into account reorientation of martensitic variants but unifies several of the existing SMA constitutive models, is used for implementation. Remarks on numerical accuracy of both algorithms are given, and it is concluded that both algorithms are applicable for this class of SMA constitutive models and preference can only be given based on the computational cost. Copyright © 2000 John Wiley & Sons, Ltd.     

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.     

Numerical simulation of granular materials by an improved discrete element method

In this paper, we present an improved discrete element method based on the non‐smooth contact dynamics and the bi‐potential concept. The energy dissipated during the collisions is taken into account by means of restitution coefficients. The interaction between particles is modelled by Coulomb unilateral contact law with dry friction which is typically non‐associated: during the contact, the sliding vector is not normal to the friction cone. The main feature of our algorithm is to overcome this difficulty by means of the bi‐potential theory. It leads to an easy implement predictor–corrector scheme involving just an orthogonal projection onto the friction cone. Moreover the convergence test is based on an error estimator in constitutive law using the corner stone inequality of the bipotential. Then we present numerical simulations which show the robustness of our algorithm and the various possibilities of the software ‘MULTICOR’ developed with this approach. Copyright © 2004 John Wiley & Sons, Ltd.     

Modelling ion diffusion mechanisms in porous media

The main features of a numerical model aiming at predicting the drift of ions in electrolytic solutions are presented. The mechanisms of ionic diffusion are described by solving the Nernst–Planck system of equations. The electrical coupling between the various ionic fluxes is accounted for by the Poisson equation. Two algorithms using the finite element method for spatial discretization are compared for simple test cases. One is based on the Picard iteration method while the other is based on the Newton–Raphson scheme. Test results clearly indicate that the range of application is broader for the algorithm based on the Newton–Raphson method. Selected examples of the application of the algorithm to more complex 1‐D and 2‐D cases are given. Copyright © 1999 John Wiley & Sons, Ltd.     

A Jacobian‐free‐based IIM for incompressible flows involving moving interfaces with Dirichlet boundary conditions

In this paper, a finite difference marker‐and‐cell (MAC) scheme is presented for the steady Stokes equations with moving interfaces and Dirichlet boundary condition. The moving interfaces are represented by Lagrangian control points and their position is updated implicitly using a Jacobian‐free approach within each time step. The forces at the moving interfaces are calculated from the position of the interfaces and interpolated using cubic splines and then applied to the fluid through the related jump conditions. The proposed Jacobian‐free Newton–generalized minimum residual (GMRES) method avoids the need to form and store the matrix explicitly in the computation of the inverse of the Jacobian and betters numerical stability. The Stokes equations are discretized on a MAC grid via a second‐order finite difference scheme with the incorporation of jump contributions and the resulting saddle point system is solved by the conjugate gradient Uzawa‐type method. Numerical results demonstrate very well the accuracy and effectiveness of the proposed method. The present algorithm has been applied to solve incompressible Navier–Stokes flows with moving interfaces. Copyright © 2010 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.     

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.