A family of multi‐point flux approximation schemes for general element types in two and three dimensions with convergence performance

A family of flux‐continuous, locally conservative, control‐volume‐distributed multi‐point flux approximation (CVD‐MPFA) schemes has been developed for solving the general geometry‐permeability tensor pressure equation on structured and unstructured grids. These schemes are applicable to the full‐tensor pressure equation with generally discontinuous coefficients and remove the O(1) errors introduced by standard reservoir simulation schemes when applied to full‐tensor flow approximation. The family of flux‐continuous schemes is characterized by a quadrature parameterization. Improved numerical convergence for the family of CVD‐MPFA schemes using the quadrature parameterization has been observed for structured and unstructured grids in two dimensions. The CVD‐MPFA family cell‐vertex formulation is extended to classical general element types in 3‐D including prisms, pyramids, hexahedra and tetrahedra. A numerical convergence study of the CVD‐MPFA schemes on general unstructured grids comprising of triangular elements in 2‐D and prismatic, pyramidal, hexahedral and tetrahedral shape elements in 3‐D is presented. Copyright © 2011 John Wiley & Sons, Ltd.     

Operator-splitting methods for the incompressible Navier-Stokes equations on non-staggered grids. Part 1: First-order schemes

New implicit finite difference schemes for solving the time-dependent incompressible Navier-Stokes equations using primitive variables and non-staggered grids are presented in this paper. A priori estimates for the discrete solution of the methods are obtained. Employing the operator approach, some requirements on the difference operators of the scheme are formulated in order to derive a scheme which is essentially consistent with the initial differential equations. The operators of the scheme inherit the fundamental properties of the corresponding differential operators and this allows a priori estimates for the discrete solution to be obtained. The estimate is similar to the corresponding one for the solution of the differential problem and guarantees boundedness of the solution. To derive the consistent scheme, special approximations for convective terms and div and grad operators are employed. Two variants of time discretization by the operator-splitting technique are considered and compared. It is shown that the derived scheme has a very weak restriction on the time step size. A lid-driven cavity flow has been predicted to examine the stability and accuracy of the schemes for Reynolds number up to 3200 on the sequence of grids with 21 × 21, 41 × 41, 81 × 81 and 161 × 161 grid points.     

On the characterization of grid density in grid refinement studies for discretization error estimation

This paper discusses the estimation of discretization errors on the basis of power series expansions for grid sets that are not geometrically similar, that is, grids not exhibiting a constant grid refinement ratio for the entire computational domain. Simple test cases with structured and unstructured grids are used to demonstrate that reliable error estimates on the basis of power series expansions can be made if the grids are refined systematically. However, if the grid refinement ratio is not constant in the complete domain, the definition of the typical cell size is not obvious, and the observed order of accuracy may not be equal to the expected theoretical order of the discretization. Some alternatives for the definition of the typical cell size are tested. In these tests, the error estimation does not show a significant effect of the definition of the typical cell size even for some cases with data sets clearly outside the ‘asymptotic range’. For non‐geometrically similar grids, the best estimates of the observed order of accuracy are obtained with the typical cell size on the basis of the mode of the cell size (the cell size that occurs more often in a given grid). Copyright © 2012 John Wiley & Sons, Ltd.     

Invariant discretization of the incompressible Navier-Stokes equations in boundary fitted co-ordinates

The discretization of the incompressible Navier-Stokes equation on boundary-fitted curvilinear grids is considered. The discretization is based on a staggered grid arrangement and the Navier-;Stokes equations in tensor formulation including Christoffel symbols. It is shown that discretization accuracy is much enhanced by choosing the velocity variables in a special way. The time-dependent equations are solved by a pressure-correction method in combination with a GMRES method. Special attention is paid to the discretization of several types of boundary conditions. It is shown that fairly non-smooth grids may be used using our approach.     

Two-grid methods for semi-linear elliptic interface problems by immersed finite element methods

In this paper, two-grid immersed finite element(IFE) algorithms are proposed and analyzed for semi-linear interface problems with discontinuous diffusion coefficients in two dimension. Because of the advantages of finite element(FE) formulation and the simple structure of Cartesian grids, the IFE discretization is used in this paper. Two-grid schemes are formulated to linearize the FE equations. It is theoretically and numerically illustrated that the coarse space can be selected as coarse asH= O(h~(1/4))(orH=O(h~(1/8))), and the asymptotically optimal approximation can be achieved as the nonlinear schemes. As a result, we can settle a great majority of nonlinear equations as easy as linearized problems. In order to estimate the present two-grid algorithms, we derive the optimal error estimates of the IFE solution in theL pnorm. Numerical experiments are given to verify the theorems and indicate that the present two-grid algorithms can greatly improve the computing efficiency.     

FURTHER DISCUSSION OF NUMERICAL ERRORS IN CFD

The methods of estimating numerical errors given in an earlier paper are extended in directions that make them useful in actual CFD applications. In particular, the method of estimating convergence error (the error due to insufficient iteration) is extended to allow the possibility of complex eigenvalues; an ad hoc method that can be applied to any case is also given. For the discretization error, which arises from the numerical approximation of the differential equation(s), methods that can be used on non-uniform drids are presented; they can be extended to unstructured grids as well. The utility of these methods is demonstrated for linear problems as well as solutions of the Navier-Stokes equations. The examples show that the estimation of errors is neither difficult nor expensive.     

Application of a perfectly matched layer to the nonlinear wave equation

We derive a perfectly matched layer-like damping layer for the nonlinear wave equation. In the layer, only two auxiliary variables are needed. In the linear case the layer is perfectly matched, but in the nonlinear case it is not. Well posedness is established for the linear case. We also prove various energy estimates which can be used as a starting point for establishing stability of more general cases. In particular, we are able to show estimates for a special type of nonlinearity.Numerical experiments that show the effectiveness of the layer are presented both for nonlinear and linear problems. In the computations, we use an eighth order summation-by-parts discretization in space and implement the boundary conditions using a penalty procedure. We present new stability results for this discretization applied to the second order wave equation in the case with Dirichlet boundary conditions.     

Anisotropy favoring triangulation CVD(MPFA) finite‐volume approximations

A family of flux‐continuous, control‐volume distributed multi‐point flux approximation schemes CVD (MPFA) have been developed for solving the general geometry‐permeability tensor pressure equation on structured and unstructured grids (Comput. Geo. 1998; 2 : 259–290, Comput. Geo. 2002; 6 : 433–452). The locally conservative schemes are applicable to the diagonal and full‐tensor pressure equation with generally discontinuous coefficients and remove the O(1) errors introduced by standard reservoir simulation schemes when applied to full‐tensor flow approximation. The family of flux‐continuous schemes is quantified by a quadrature parameterization. Improved numerical convergence for the family of CVD(MPFA) schemes for specified quadrature points has been observed for lower anisotropy ratios for both structured and unstructured grids in two dimensions. However, for strong full‐tensor anisotropy fields the quadrilateral schemes can induce strong spurious oscillations in the numerical solution. This paper motivates and demonstrates the benefit of using anisotropy favoring triangulation for treating such cases. Test examples involving strong full‐tensor anisotropy fields are presented in 2‐D and 3‐D, which show that the family of schemes on anisotropy favoring triangulation (prisms in 3‐D) yield well‐resolved pressure fields with little or no spurious oscillations. Copyright © 2010 John Wiley & Sons, Ltd.     

Fourier analysis of Schwarz domain decomposition methods for the biharmonic equation

Schwarz methods are an important type of domain decomposition methods. Using the Fourier transform, we derive error propagation matrices and their spectral radii of the classical Schwarz alternating method and the additive Schwarz method for the biharmonic equation in this paper. We prove the convergence of the Schwarz methods from a new point of view, and provide detailed information about the convergence speeds and their dependence on the overlapping size of subdomains. The obtained results are independent of any unknown constant and discretization method, showing that the Schwarz alternating method converges twice as quickly as the additive Schwarz method.     

Selection of a staggered grid for inertia‐gravity waves in shallow water

The problem of accuracy in propagating inertia‐gravity waves on Arakawa grids is investigated. It is shown that the sole analysis of spatial discretization and the recommendation of the B‐grid for coarse resolution models and C‐grid for high resolution models must be re‐analysed when time discretization is taken into account as well. For a chosen time discretization, a coarse C‐grid is shown for example, to increase precision when using larger time‐steps (up to the stability limit) whereas the precision of the B‐grid decreases. Here, an analysis of error for different grids in function of the space–time resolutions and computational costs is presented and some recommendations on the choice of the particular staggered grid for a given application are outlined. Copyright © 2002 John Wiley & Sons, Ltd.     

Assessment of adaptive and heuristic time stepping for variably saturated flow

The performance of improved initial estimates and ‘heuristic’ and ‘adaptive’ techniques for time step control in the iterative solution of Richards equation is evaluated. The so‐called heuristic technique uses the convergence behaviour of the iterative scheme to estimate the next time step whereas the adaptive technique regulates the time step on the basis of an approximation of the local time truncation error. The sample problems used to assess these various schemes are characterized by nonuniform (in time) boundary conditions, sharp gradients in the infiltration fronts, and discontinuous derivatives in the soil hydraulic properties. It is found that higher order initial solution estimates improve the convergence of the iterative scheme for both the heuristic and adaptive techniques, with greater overall performance gains for the heuristic scheme, as could be expected. It is also found that the heuristic technique outperforms the adaptive method under strongly nonlinear conditions. Previously reported observations suggesting that adaptive techniques perform best when accuracy requirements on the numerical solution are very stringent are confirmed. Overall both heuristic and adaptive techniques have their limitations, and a more general or mixed time stepping strategy combining truncation error and convergence criteria is recommended for complex problems. Copyright © 2006 John Wiley & Sons, Ltd.     

Geometrical interpretation of the multi‐point flux approximation L‐method

In this paper, we first investigate the influence of different Dirichlet boundary discretizations on the convergence rate of the multi‐point flux approximation (MPFA) L‐method by the numerical comparisons between the MPFA O‐ and L‐method, and show how important it is for this new method to handle Dirichlet boundary conditions in a suitable way. A new Dirichlet boundary strategy is proposed, which in some sense can well recover the superconvergence rate of the normal velocity. In the second part of the work, the MPFA L‐method with homogeneous media is studied. A systematic concept and geometrical interpretations of the L‐method are given and illustrated, which yield more insight into the L‐method. Finally, we apply the MPFA L‐method for two‐phase flow in porous media on different quadrilateral grids and compare its numerical results for the pressure and saturation with the results of the two‐point flux approximation method. Copyright © 2008 John Wiley & Sons, Ltd.     

Consistent discretization for simulations of flows with moving generalized curvilinear coordinates

We develop a consistent discretization of conservative momentum and scalar transport for the numerical simulation of flow using a generalized moving curvilinear coordinate system. The formulation guarantees consistency between the discrete transport equation and the discrete mass conservation equation due to grid motion. This enables simulation of conservative transport using generalized curvilinear grids that move arbitrarily in three dimensions while maintaining the desired properties of the discrete transport equation on a stationary grid, such as constancy, conservation, and monotonicity. In addition to guaranteeing consistency for momentum and scalar transport, the formulation ensures geometric conservation and maintains the desired high‐order time accuracy of the discretization on a moving grid. Through numerical examples we show that, when the computation is carried out on a moving grid, consistency between the discretized scalar advection equation and the discretized equation for flow mass conservation due to grid motion is required in order to obtain stable and accurate results. We also demonstrate that significant errors can result when non‐consistent discretizations are employed. Copyright © 2009 John Wiley & Sons, Ltd.     

Inverse finite element formulations for transient heat conduction problems

Transient heat conduction problems are normally simulated by the conventional consistent and lumped finite element methods. The discretization error and the physical reality violation in such problems are noticeable and unwanted responses are observed in the results when using the consistent formulations. Although in utilizing the lumped formulations, the violation of physical reality becomes reduced; however, emerging the discretization error would also become obvious to the degree of being quantifiable. In using the inverse finite element method without considering the element shape functions, the element matrices will be obtained by minimizing the governing equation and its generalized discretized corresponding formula. The results obtained by using this method indicates that the reduction in both the discretization error as well as the violation of physical reality would be realized.     

Spectral Element Discretization of the Circular Driven Cavity. Part IV: The Navier-Stokes Equations

This paper deals with the spectral element discretization of the Navier-Stokes equations in a disk with discontinuous boundary data, which is known as the driven cavity problem. The numerical treatment does not involve any regularization of these data. Relying on a variational formulation in the primitive variables of velocity and pressure, we describe a discretization of these equations and derive error estimates in appropriate weighted Sobolev spaces. We propose an algorithm to solve the nonlinear discrete system and present numerical experiments to verify its efficiency.     

Spectral approximation of the Navier-Stokes equations coupled with the heat equation

This paper is devoted to the approximation by spectral methods of a problem for viscous incompressible flows, where the Navier-Stokes equations are coupled with the heat equation. The numerical analysis of theses discretizations uses the discrete implicit function theorem of BREZZI, RAPPAZ and RAVIART that allows us to derive some error estimates.     

An optimizing implicit difference scheme based on proper orthogonal decomposition for the two‐dimensional unsaturated soil water flow equation

An optimizing reduced implicit difference scheme (IDS) based on singular value decomposition (SVD) and proper orthogonal decomposition (POD) for the two‐dimensional unsaturated soil water flow equation is presented. An ensemble of snapshots is compiled from the transient solutions derived from the usual IDS for a two‐dimensional unsaturated flow equation. Then, optimal orthogonal bases are reconstructed by implementing SVD and POD techniques for the ensemble of snapshots. Combining POD with a Galerkin projection approach, a new lower dimensional and highly accurate IDS for the two‐dimensional unsaturated flow equation is obtained. Error estimates between the true solution, the usual IDS solution, and the reduced IDS solution based on POD basis are derived. Finally, it is shown by means of a numerical example using the technology of local refined grids that the computational load is greatly diminished by using the reduced IDS. Also, the error between the POD approximate solution and the usual IDS solution is proved to be consistent with the derived theoretical results. Thus, both feasibility and efficiency of the POD method are validated. Copyright © 2011 John Wiley & Sons, Ltd.     

Steady Flows in Unsaturated Soils Are Stable

The stability of steady, vertically upward and downward flow of water in a homogeneous layer of soil is analyzed. Three equivalent dimensionless forms of the Richards equation are introduced, namely the pressure head, saturation, and matric flux potential forms. To illustrate general results and derive special results, use is made of several representative classes of soils. For all classes of soils with a Lipschitz continuous relationship between the hydraulic conductivity and the matric flux potential, steady flows are shown to be unique. In addition, linear stability of these steady flows is proved. To this end, use is made of the energy method, in which one considers (weighted) L ²-norms of the perturbations of the steady flows. This gives a general restriction of the dependence of the hydraulic conductivity upon the matric flux potential, yielding linear stability and exponential decay with time of a specific weighted L ²-norm. It is shown that for other norms the ultimate decay towards the steady-solution is preceded by transient growth. An extension of the Richards equation to take into account dynamic memory effects is also considered. It is shown that the stability condition for the standard Richards equation implies linear stability of the steady solution of the extended model.     

Comparison of Iterative Methods for Improved Solutions of the Fluid Flow Equation in Partially Saturated Porous Media

Abstract. The Picard and modified Picard iteration schemes are often used to numerically solve the nonlinear Richards equation governing water flow in variably saturated porous media. While these methods are easy to implement, they are only linearly convergent. Another approach to solve the Richards equation is to use Newton's iterative method. This method, also known as Newton–Raphson iteration, is quadratically convergent and requires the computation of first derivatives. We implemented Newton's scheme into the mixed form of the Richards equation. As compared to the modified Picard scheme, Newton's scheme requires two additional matrices when the mixed form of the Richards equation is used and requires three additional matrices, when the pressure head-based form is used. The modified Picard scheme may actually be viewed as a simplified Newton scheme.Two examples are used to investigate the numerical performance of different forms of the 1D vertical Richards equation and the different iterative solution schemes. In the first example, we simulate infiltration in a homogeneous dry porous medium by solving both, the h based and mixed forms of Richards equation using the modified Picard and Newton schemes. Results shows that, very small time steps are required to obtain an accurate mass balance. These small times steps make the Newton method less attractive.In a second test problem, we simulate variable inflows and outflows in a heterogeneous dry porous medium by solving the mixed form of the Richards equation, using the modified Picard and Newton schemes. Analytical computation of the Jacobian required less CPU time than its computation by perturbation. A combination of the modified Picard and Newton scheme was found to be more efficient than the modified Picard or Newton scheme.     

An explicit finite volume scheme on staggered grids for the Euler equations: Unstructured meshes,stability analysis,and energy conservation

We set up a numerical strategy for the simulation of the Euler equations, in the framework of finite volume staggered discretizations where numerical densities, energies, and velocities are stored on different locations. The main difficulty relies on the treatment of the total energy, which mixes quantities stored on different grids. The proposed method is strongly inspired, on the one hand, from the kinetic framework for the definition of the numerical fluxes, and, on the other hand, from the discrete duality finite volume (DDFV) framework, which has been designed for the simulation of elliptic equations on complex meshes. The time discretization is explicit and we exhibit stability conditions that guaranty the positivity of the discrete densities and internal energies. Moreover, while the scheme works on the internal energy equation, we can define a discrete total energy which satisfies a local conservation equation. We provide a set of numerical simulations to illustrate the behavior of the scheme.