A mass-conservative control volume-finite element method for solving Richards’ equation in heterogeneous porous media

We present a mass-conservative vertex-centred finite volume method for efficiently solving the mixed form of Richards’ equation in heterogeneous porous media. The spatial discretisation is particularly well-suited to heterogeneous media because it produces consistent flux approximations at quadrature points where material properties are continuous. Combined with the method of lines, the spatial discretisation gives a set of differential algebraic equations amenable to solution using higher-order implicit solvers. We investigate the solution of the mixed form using a Jacobian-free inexact Newton solver, which requires the solution of an extra variable for each node in the mesh compared to the pressure-head form. By exploiting the structure of the Jacobian for the mixed form, the size of the preconditioner is reduced to that for the pressure-head form, and there is minimal computational overhead for solving the mixed form.     

A new multipoint flux approximation method with a quasi-local stencil (MPFA-QL) for the simulation of diffusion problems in anisotropic and heterogeneous media

In this paper, we present a new linear cell-centered finite volume multipoint flux approximation (MPFA-QL) scheme for discretizing diffusion problems on general polygonal meshes. This scheme uses a quasi-local stencil, based upon the conormal decomposition, to approximate the control face flux when solving the steady state diffusion problem, being able to reproduce piecewise linear solutions exactly and it is very robust when dealing with heterogeneous and highly anisotropic media and severely distorted meshes. In our linear scheme, we first construct the one-sided fluxes on each control surface independently and then a unique flux expression is obtained by a convex combination of the one-sided fluxes. The unknown values at the vertices that define a control surface are interpolated by means of a linearity-preserving interpolation procedure, considering control volumes surrounding these vertices. To show the potential of the MPFA-QL scheme, we solve some benchmark using triangular and quadrilateral meshes and we compare our scheme with other numerical formulations found in literature.     

A new fractional finite volume method for solving the fractional diffusion equation

The inherent heterogeneities of many geophysical systems often gives rise to fast and slow pathways to water and chemical movement. One approach to model solute transport through such media is by fractional diffusion equations with a space–time dependent variable coefficient. In this paper, a two-sided space fractional diffusion model with a space–time dependent variable coefficient and a nonlinear source term subject to zero Dirichlet boundary conditions is considered.Some finite volume methods to solve a fractional differential equation with a constant dispersion coefficient have been proposed. The spatial discretisation employs fractionally-shifted Grünwald formulas to discretise the Riemann–Liouville fractional derivatives at control volume faces in terms of function values at the nodes. However, these finite volume methods have not been extended to two-dimensional and three-dimensional problems in a natural manner. In this paper, a new weighted fractional finite volume method with a nonlocal operator (using nodal basis functions) for solving this two-sided space fractional diffusion equation is proposed. Some numerical results for the Crank–Nicholson fractional finite volume method are given to show the stability, consistency and convergence of our computational approach. This novel simulation technique provides excellent tools for practical problems even when a complex transition zone is involved. This technique can be extend to two-dimensional and three-dimensional problems with complex regions.     

A superconvergent fitted finite volume method for Black–Scholes equations governing European and American option valuation

We develop a superconvergent fitted finite volume method for a degenerate nonlinear penalized Black–Scholes equation arising in the valuation of European and American options, based on the fitting idea in Wang [IMA J Numer Anal 24 (2004), 699–720]. Unlike conventional finite volume methods in which the dual mesh points are naively chosen to be the midpoints of the subintervals of the primal mesh, we construct the dual mesh judiciously using an error representation for the flux interpolation so that both the approximate flux and solution have the second‐order accuracy at the mesh points without any increase in computational costs. As the equation is degenerate, we also show that it is essential to refine the meshes locally near the degenerate point in order to maintain the second‐order accuracy. Numerical results for both European and American options with constant and nonconstant coefficients will be presented to demonstrate the superconvergence of the method. © 2014 Wiley Periodicals, Inc. Numer Methods Partial Differential Eq 31: 1190–1208, 2015     

Robust intersection of structured hexahedral meshes and degenerate triangle meshes with volume fraction applications

Two methods for calculating the volume and surface area of the intersection between a triangle mesh and a rectangular hexahedron are presented. The main result is an exact method that calculates the polyhedron of intersection and thereafter the volume and surface area of the fraction of the hexahedral cell inside the mesh. The second method is approximate, and estimates the intersection by a least squares plane. While most previous publications focus on non-degenerate triangle meshes, we here extend the methods to handle geometric degeneracies. In particular, we focus on large-scale triangle overlaps, or double surfaces. It is a geometric degeneracy that can be hard to solve with existing mesh repair algorithms. There could also be situations in which it is desirable to keep the original triangle mesh unmodified. Alternative methods that solve the problem without altering the mesh are therefore presented. This is a step towards a method that calculates the solid area and volume fractions of a degenerate triangle mesh including overlapping triangles, overlapping meshes, hanging nodes, and gaps. Such triangle meshes are common in industrial applications. The methods are validated against three industrial test cases. The validation shows that the exact method handles all addressed geometric degeneracies, including double surfaces, small self-intersections, and split hexahedra.     

A robust iterative scheme for finite volume discretization of diffusive flux on highly skewed meshes

A new approach for diffusive flux discretization on a nonorthogonal mesh for finite volume method is proposed. This approach is based on an iterative method, Deferred correction introduced by M. Peric [J.H. Fergizer, M. Peric, Computational Methods for Fluid Dynamics, Springer, 2002]. It converges on highly skewed meshes where the former approach diverges. A convergence proof of our method is given on arbitrary quadrilateral control volumes. This proof is founded on the analysis of the spectral radius of the iteration matrix. This new approach is applied successfully to the solution of a Poisson equation in quadrangular domains, meshed with highly skewed control volumes. The precision order of used schemes is not affected by increasing skewness of the grid. Some numerical tests are performed to show the accuracy of the new approach.     

Finite Volume Methods for Multi-Symplectic PDES

We investigate the application of a cell-vertex finite volume discretization to multi-symplectic PDEs. The investigated discretization reduces to the Preissman box scheme when used on a rectangular grid. Concerning arbitrary quadrilateral grids, we show that only methods with parallelogram-like finite volume cells lead to a multi-symplectic discretization; i.e., to a method that preserves a discrete conservation law of symplecticity. One of the advantages of finite volume methods is that they can be easily adjusted to variable meshes. But, although the implementation of moving mesh finite volume methods for multi-symplectic PDEs is rather straightforward, the restriction to parallelogram-like cells implies that only meshes moving with a constant speed are multi-symplectic. To overcome this restriction, we suggest the implementation of reversible moving mesh methods based on a semi-Lagrangian approach. Numerical experiments are presented for a one dimensional dispersive shallow-water system.     

Graded Meshes for Higher Order FEM

A singularly perturbed one-dimensional convection-diffusion problem is solved numerically by the finite element method based on higher order polynomials. Numerical solutions are obtained using S-type meshes with special emphasis on meshes which are graded (based on a mesh generating function) in the fine mesh region. Error estimates in the ε-weighted energy norm are proved. We derive an 'optimal' mesh generating function in order to minimize the constant in the error estimate. Two layer-adapted meshes defined by a recursive formulae in the fine mesh region are also considered and a new technique for proving error estimates for these meshes is presented. The aim of the paper is to emphasize the importance of using optimal meshes for higher order finite element methods. Numerical experiments support all theoretical results.     

Planar mesh refinement cannot be both local and regular

We show that two desirable properties for planar mesh refinement techniques are incompatible. Mesh refinement is a common technique for adaptive error control in generating unstructured planar triangular meshes for piecewise polynomial representations of data. Local refinements are modifications of the mesh that involve a fixed maximum amount of computation, independent of the number of triangles in the mesh. Regular meshes are meshes for which every interior vertex has degree 6. At least for some simple model meshing problems, optimal meshes are known to be regular, hence it would be desirable to have a refinement technique that, if applied to a regular mesh, produced a larger regular mesh. We call such a technique a regular refinement. In this paper, we prove that no refinement technique can be both local and regular. Our results also have implications for non-local refinement techniques such as Delaunay insertion or Rivara's refinement. Received August 1, 1996 / Revised version received February 28, 1997     

An incomplete LU-based family of preconditioners for numerical resolution of a shallow water system using a duality method—applications

In this paper, we present a family of preconditioners well adapted to the solution of linear problems that arise from a particular discretisation of shallow water equations in the flux form. The formulation of the shallow water equations used here is discretised in time using the method of characteristics and the Euler implicit method, and solved by a duality technique with automatic choice of parameters. The space discretisation is performed using the first-order Raviart-Thomas finite element. The family of preconditioners designed for solving the linear problems that appear at each time iteration greatly improves convergence and significantly reduces the CPU time needed to solve them.     

Mass Conservative Approximation of Single and Two Phase Flow in Porous Media

In this paper a new error estimate for a multi point flux approximation ( MPFA ) control volume method on triangular meshes for approximating flow in porous media is presented. The scheme is mass conservative and has shown to simulate reliably flows with discontinuous permeability tensors and on irregular grids which is the major challenge in hydrological engineering. (© 2008 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim)     

Finite element and boundary element contact stress analysis with remeshing technique

The fundamental part of the contact stress problem solution using a finite element method is to locate possible contact areas reliably and efficiently. In this research, a remeshing technique is introduced to determine the contact region in a given accuracy. In the proposed iterative method, the meshes near the contact surface are modified so that the edge of the contact region is also an element’s edge. This approach overcomes the problem of surface representation at the transition point from contact to non-contact region. The remeshing technique is efficiently employed to adapt the mesh for more precise representation of the contact region. The method is applied to both finite element and boundary element methods. Overlapping of the meshes in the contact region is prevented by the inclusion of displacement and force constraints using the Lagrange multipliers technique. Since the method modifies the mesh only on the contacting and neighbouring region, the solution to the matrix system is very close to the previous one in each iteration. Both direct and iterative solver performances on BEM and FEM analyses are also investigated for the proposed incremental technique. The biconjugate gradient method and LU with Cholesky decomposition are used for solving the equation systems. Two numerical examples whose analytical solutions exist are used to illustrate the advantages of the proposed method. They show a significant improvement in accuracy compared to the solutions with fixed meshes.     

Robust a posteriori error estimates for conforming discretizations of a singularly perturbed reaction-diffusion problem on anisotropic meshes

In this paper, we derive robust a posteriori error estimates for conforming approximations to a singularly perturbed reaction-diffusion problem on anisotropic meshes, since the solution in general exhibits anisotropic features, e.g., strong boundary or interior layers. Based on the anisotropy of the mesh elements, we improve the a posteriori error estimates developed by Cheddadi et al., which are reliable and efficient on isotropic meshes but fail on anisotropic ones. Without the assumption that the mesh is shape-regular, the resulting mesh-dependent error estimator is shown to be reliable, efficient and robust with respect to the reaction coefficient, as long as the anisotropic mesh sufficiently reflects the anisotropy of the solution. We present our results in the framework of the vertex-centered finite volume method but their nature is general for any conforming one, like the piecewise linear finite element one. Our estimates are based on the usual H(div)-conforming, locally conservative flux reconstruction in the lowest-order Raviart-Thomas space on a dual mesh associated with the original anisotropic simplex one. Numerical experiments in 2D confirm that our estimates are reliable, efficient and robust on anisotropic meshes.     

Monotone combined edge finite volume–finite element scheme for Anisotropic Keller–Segel model

In this article, a new numerical scheme for a degenerate Keller–Segel model with heterogeneous anisotropic tensors is treated. It is well‐known that standard finite volume scheme not permit to handle anisotropic diffusion without any restrictions on meshes. Therefore, a combined finite volume‐nonconforming finite element scheme is introduced, developed, and studied. The unknowns of this scheme are the values at the center of cell edges. Convergence of the approximate solution to the continuous solution is proved only supposing the shape regularity condition for the primal mesh. This scheme ensures the validity of the discrete maximum principle under the classical condition that all transmissibilities coefficients are positive. Therefore, a nonlinear technique is presented, as a correction of the diffusive flux, to provide a monotone scheme for general tensors. © 2014 Wiley Periodicals, Inc. Numer Methods Partial Differential Eq 30: 1030–1065, 2014     

An iterative method for the creation of structured hexahedral meshes over complex orography

In this paper we propose a technique for measuring the quality of hexahedral Cartesian meshes used to model meso-scale atmospheric circulation in 3D. It is used to verify the progress of a novel method for satisfying the Delaunay criterion for structured hexahedral meshes over complex orography with high gradients and wide gradient variability. Based on a simile with potential energy, the iterative method of mesh smoothing is shown to improve mesh quality with logarithmic convergence. The method is evaluated in a practical application in a specific geographic location.     

On using meshes of mixed element types in adaptive finite element analysis

Four different automatic mesh generators capable of generating either triangular meshes or hybrid meshes of mixed element types have been used in the mesh generation process. The performance of these mesh generators were tested by applying them to the adaptive finite element refinement procedure. It is found that by carefully controlling the quality and grading of the quadrilateral elements, an increase in efficiency over pure triangular meshes can be achieved. Furthermore, if linear elements are employed, an optimal hybrid mesh can be obtained most economically by a combined use of the mesh coring technique suggested by Lo and Lau and a selective removal of diagonals from the triangular element mesh. On the other hand, if quadratic elements are used, it is preferable to generate a pure triangular mesh first, and then obtain a hybrid mesh by merging of triangles.     

多边形网格上扩散方程新的单调格式

 本文在星形多边形网格上, 构造了扩散方程新的单调有限体积格式.该格式与现有的基于非线性两点流的单调格式的主要区别是, 在网格边的法向流离散模板中包含当前边上的点, 在推导离散法向流的表达式时采用了定义于当前边上的辅助未知量, 这样既可适应网格几何大变形, 同时又兼顾了当前网格边上物理量的变化. 在光滑解情形证明了离散法向流的相容性.对于具有强各向异性、非均匀张量扩散系数的扩散方程, 证明了新格式是单调的, 即格式可以保持解析解的正性. 数值结果表明在扭曲网格上, 所构造的格式是局部守恒和保正的, 对光滑解有高于一阶的精度, 并且, 针对非平衡辐射限流扩散问题, 数值结果验证了新格式在计算效率和守恒精度上优于九点格式.     

Une méthode de volumes finis pour les équations elliptiques du second ordre

We present a new finite volume method for approximating second order elliptic equations. This method has several advantages: — it allows large mesh distortions which occur in Lagrangian hydrodynamics calculations; — it is convenient for the approximation of crossed derivative terms as those which take place in magnetohydrodynamics will Hall effect; — it generalizes the finite difference method and the finite volume method using Delaunay-Voronoi meshes.     

High order upwind scheme for modelling turbulent shallow water flow in hydraulic structures

An unstructured finite volume model for quasi-2D free surface flow with wet-dry fronts and turbulence modelling is presented. The convective flux is discretised with either a an hybrid second-order/first-order scheme, or a fully second order scheme, both of them upwind Godunov's schemes based on Roe's average. The hybrid scheme uses a second order discretisation for the two unit discharge components, whilst keeping a first order discretisation for the water depth [2]. In such a way the numerical diffusion is much reduced, without a significant reduction on the numerical stability of the scheme, obtaining in such a way accurate and stable results. It is important to keep the numerical diffusion to a minimum level without loss of numerical stability, specially when modelling turbulent flows, because the numerical diffusion may interfere with the real turbulent diffusion. In order to avoid spurious oscillations of the free surface when the bathymetry is irregular, an upwind discretisation of the bed slope source term [4] with second order corrections is used [2]. In this way a fully second order scheme which gives an exact balance between convective flux and bed slope in the hydrostatic case is obtained. The k – ε equations are solved with either an hybrid or a second order scheme. In all the numerical simulations the importance of using a second order upwind spatial discretisation has been checked [1]. A first order scheme may give rather good predictions for the water depth, but it introduces too much numerical diffusion and therefore, it excessively smooths the velocity profiles. This is specially important when comparing different turbulence models, since the numerical diffusion introduced by a first order upwind scheme may be of the same order of magnitude as the turbulent diffusion. (© 2008 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim)     

Quadrilateral surface meshes without self-intersecting dual cycles for hexahedral mesh generation

Several promising approaches for hexahedral mesh generation work as follows: Given a prescribed quadrilateral surface mesh they first build the combinatorial dual of the hexahedral mesh. This dual mesh is converted into the primal hexahedral mesh, and finally embedded and smoothed into the given domain. Two such approaches, the modified whisker weaving algorithm by Folwell and Mitchell, as well as a method proposed by the author, rely on an iterative elimination of certain dual cycles in the surface mesh. An intuitive interpretation of the latter method is that cycle eliminations correspond to complete sheets of hexahedra in the volume mesh.

Although these methods can be shown to work in principle, the quality of the generated meshes heavily relies on the dual cycle structure of the given surface mesh. In particular, it seems that difficulties in the hexahedral meshing process and poor mesh qualities are often due to self-intersecting dual cycles. Unfortunately, all previous work on quadrilateral surface mesh generation has focused on quality issues of the surface mesh alone but has disregarded its suitability for a high-quality extension to a three-dimensional mesh.

In this paper, we develop a new method to generate quadrilateral surface meshes without self-intersecting dual cycles. This method reuses previous b-matching problem formulations of the quadrilateral mesh refinement problem. The key insight is that the b-matching solution can be decomposed into a collection of simple cycles and paths of multiplicity two, and that these cycles and paths can be consistently embedded into the dual surface mesh.

A second tool uses recursive splitting of components into simpler subcomponents by insertion of internal two-manifolds. We show that such a two-manifold can be meshed with quadrilaterals such that the induced dual cycle structure of each subcomponent is free of self-intersections if the original component satisfies this property. Experiments show that we can achieve hexahedral meshes with a good quality.