Cell-centered Finite Volume Method Research Articles

A very high-order flux reconstruction approach coupled to the MPFA-QL finite volume method for the numerical simulation of oil-water flows in 2D petroleum reservoirs

In recent years, high-order methods for the discretization of the advection dominated saturation equation, and that results from the modelling of the 2D Oil-Water displacement through porous formations, have received renewed attention from the CFD community, as a suitable alternative to the traditional first-order upwind (FOU) finite volume (FV) scheme, which is the standard in industry. A key aspect of high-order methods applied to numerical simulation of petroleum reservoirs is the lower computational cost for a comparable degree of accuracy against their lower-order counterparts. In this context, in the present paper, for the first time in literature, we propose the use of the very high order Flux Reconstruction (FR) formulation for the discretization of the advective terms in the saturation equation using unstructured quadrilateral meshes. To solve the elliptic pressure equation, we use the recently developed Multipoint Flux Approximation method with a Quasi-Local Stencil (MPFA-QL). This method is very robust and able to handle highly heterogeneous and anisotropic permeability tensors. The system of equations is solved in a segregated manner by using the classical Implicit Pressure Explicit Saturation (IMPES) procedure. To properly couple the cell-centered MPFA-QL method with the FR formulation, it is necessary to obtain an adequate velocity reconstruction throughout the control volumes of the mesh. As the cell-centered finite volume method naturally delivers fluxes across control surfaces of the primal grid, then, we use a reconstruction operator based on the first-order Brezzi–Douglas–Marini (BDM) space, and the Piola transformation, to get the complete knowledge of the velocity field throughout the domain. To eliminate spurious oscillations around the discontinuities, which are typical in high-order methods, we have used a hierarchical multidimensional limiting process (hMLP). To compute fluxes over the control surfaces we have used the Roe-E approximate Riemann solver. Finally, through some representative examples, we show the accuracy, efficiency and robustness of our new high order numerical methodology to model oil-water displacements in highly heterogeneous and anisotropic petroleum reservoirs.

Numerical investigation of full helicopter with and without the ground effect

In the present work, the aerodynamic performance of the full helicopter PSP in hover flight is investigated using a simplified concept of multiple reference frame (MRF) technique in the context of high-order Monotone Upstream Centred Scheme for Conservation Laws (MUSCL) cell-centred finite volume method. The predictions were obtained for two ground distances and several collective pitch angle at tip Mach number of 0.585. The calculations were made for both out-of-ground-effect (OGE) and in-ground-effect (IGE) cases and compared with experimental data in terms of pressure distribution and integrated thrust and torque and vortex system.

A Cell-Centered Semi-Lagrangian Finite Volume Method for Solving Two-Dimensional Coupled Burgers’ Equations

A cell-centered finite volume semi-Lagrangian method is presented for the numerical solution of two-dimensional coupled Burgers’ problems on unstructured triangular meshes. The method combines a modified method of characteristics for the time integration and a cell-centered finite volume for the space discretization. The new method belongs to fractional-step algorithms for which the convection and the viscous parts in the coupled Burgers’ problems are treated separately. The crucial step of interpolation in the convection step is performed using two local procedures accounting for the element where the departure point is located. The resulting semidiscretized system is then solved using a third-order explicit Runge-Kutta scheme. In contrast to the Eulerian-based methods, we apply the new method for each time step along the characteristic curves instead of the time direction. The performance of the current method is verified using different examples for coupled Burgers’ problems with known analytical solutions. We also apply the method for simulation of an example of coupled Burgers’ flows in a complex geometry. In these test problems, the new cell-centered finite volume semi-Lagrangian method demonstrates its ability to accurately resolve the two-dimensional coupled Burgers’ problems.

A vertex-based reconstruction for cell-centered finite-volume discretization on unstructured grids

Recently, a vertex-based spatial reconstruction for unstructured cell-centered finite volume method has been proposed, showing advantages in accuracy, convergence, and efficiency. However, this method was only applied for inviscid flows. Moreover, for shock-capturing computations, the conventional cell-based limiter was the only choice to suppress numerical oscillations. In this work, the vertex-based reconstruction is extended for the solution of viscous flows, while automatically avoiding the typical decoupling problem. Moreover, a WENO-type nonlinear weighting strategy is designed and combined with the vertex-based reconstruction method, resulting in a simple and cost-efficient alternative to conventional slope limiters. In addition, the present method is combined with the recently proposed iterative near-boundary treatment, ensuring linear exactness even in the vicinity of boundaries, without sacrificing the overall computational efficiency. A series of numerical test cases involving viscous flows and shock waves provide evidence of the superior performance of the present method. The advantages of the present method are especially prominent on high aspect-ratio irregular triangular grids, which is a beneficial property for solving realistic fluid dynamic problems at high Reynolds numbers.

A 3D cell-centered ADER MOOD Finite Volume method for solving updated Lagrangian hyperelasticity on unstructured grids

In this paper, we present a conservative cell-centered Lagrangian Finite Volume scheme for solving the hyperelasticity equations on unstructured multidimensional grids. The starting point of the present approach is the cell-centered FV discretization named EUCCLHYD and introduced in the context of Lagrangian hydrodynamics. Here, it is combined with the a posteriori Multidimensional Optimal Order Detection (MOOD) limiting strategy to ensure robustness and stability at shock waves with piecewise linear spatial reconstruction. The ADER (Arbitrary high order schemes using DERivatives) approach is adopted to obtain second-order of accuracy in time. This strategy has been successfully tested in a hydrodynamics context and the present work aims at extending it to the case of hyperelasticity. Here, the hyperelasticity equations are written in the updated Lagrangian framework and the dedicated Lagrangian numerical scheme is derived in terms of nodal solver, Geometrical Conservation Law (GCL) compliance, subcell forces and compatible discretization. The Lagrangian numerical method is implemented in 3D under MPI parallelization framework allowing to handle genuinely large meshes. A relatively large set of numerical test cases is presented to assess the ability of the method to achieve effective second order of accuracy on smooth flows, maintaining an essentially non-oscillatory behavior and general robustness across discontinuities and ensuring at least physical admissibility of the solution where appropriate. Pure elastic neo-Hookean and non-linear materials are considered for our benchmark test problems in 2D and 3D. These test cases feature material bending, impact, compression, non-linear deformation and further bouncing/detaching motions.

A New Nonlinear Two-Point Flux Approximation Method for Solving the Anisotropic Diffusion Equation with Reduced Violations of the Discrete Maximum/Minimum Principle

Summary A class of monotone cell-centered nonlinear finite-volume methods has been proposed in the past decade to solve the anisotropic diffusion equation. The nonlinear two-point flux approximation (TPFA) (NTPFA) method preserves the nonnegativity of the solution values but can violate the discrete maximum/minimum principle (DMP). To enforce DMP, the nonlinear multipoint flux approximation (NMPFA) method ought to be used. In this work, we propose a novel NTPFA method that can reduce the severity of DMP violations significantly compared with the standard NTPFA method. The new formulation uses conormal decomposition for the construction of the one-sided fluxes. To define the unique flux through a connection between two cells, we choose a convex combination of the two one-sided fluxes such that the absolute differences of the magnitudes of the two transmissibility terms associated with the two neighboring cells are minimized, thus bringing the discrete coefficient matrix closer to having the zero row-sum property. Numerical experiments are conducted to test the performance of the new NTPFA method. The results demonstrate that the new scheme has comparable convergence order for both the solution and the flux compared with the standard NTPFA method or the classical multi-point flux approximation (MPFA-O) method. Moreover, the new NTPFA formulation shows marked improvements over the standard NTPFA in terms of reducing DMP violations. However, depending on the specific problem configuration, our new NTPFA formulation can lead to a system of nonlinear equations that is more difficult to solve.

A finite volume penalty based segment-to-segment method for frictional contact problems

This paper presents a new contact boundary condition for finite volume simulations of frictional contact problems involving geometrical and material non-linearities. Deformation of bodies in contact is described by the updated Lagrangian form of the momentum equation which is discretised in space using the cell-centred finite volume method. The proposed contact boundary condition is based on the finite volume implementation of the penalty based segment-to-segment contact force calculation method in which normal contact pressure, governed by a penalty law, is integrated across the discretised contact surface, enforcing contact constraints in an integral manner. Such an approach, as opposed to the pointwise contact force calculation algorithm, allows for more accurate and more robust treatment of the contact area edge, simultaneous calculation of contact forces on both contact surfaces as well as smoother contact force during large sliding. The proposed numerical method is tested on challenging mechanical contact problems showing very good agreement with the available benchmark results.

The Fictitious Domain Method with Sharp Interface for Elasticity Systems with General Jump Embedded Boundary Conditions

In the framework of the fictitious domain methods with immersed interfaces for the elasticity problem, the present contribution is to study and numerically validate the jump-integrated boundary conditions method with sharp interface for the vector elasticity system discretized by a proposed finite volume method. The main idea of the fictitious domain approach consists in embedding the original domain of study into a geometrically larger and simpler one called the fictitious domain. Here, we present a cell-centered finite volume method to discretize the fictitious domain problem. The proposed method is numerically validated for different test cases. This work can be considered as a first step before more challenging problems such as fluid-structure interactions or moving interface problems.

A Low Dissipative and Stable Cell-Centered Finite Volume Method with the Simultaneous Approximation Term for Compressible Turbulent Flows

A simultaneous-approximation term is a non-reflecting boundary condition that is usually accompanied by summation-by-parts schemes for provable time stability. While a high-order convective flux based on reconstruction is often employed in a finite-volume method for compressible turbulent flow, finite-volume methods with the summation-by-parts property involve either equally weighted averaging or the second-order central flux for convective fluxes. In the present study, a cell-centered finite-volume method for compressible Naiver–Stokes equations was developed by combining a simultaneous-approximation term based on extrapolation and a low-dissipative discretization method without the summation-by-parts property. Direct numerical simulations and a large eddy simulation show that the resultant combination leads to comparable non-reflecting performance to that of the summation-by-parts scheme combined with the simultaneous-approximation term reported in the literature. Furthermore, a characteristic boundary condition was implemented for the present method, and its performance was compared with that of the simultaneous-approximation term for a direct numerical simulation and a large eddy simulation to show that the simultaneous-approximation term better maintained the average target pressure at the compressible flow outlet, which is useful for turbomachinery and aerodynamic applications, while the characteristic boundary condition better preserved the flow field near the outlet.

Finite-volume flux reconstruction and semi-analytical particle tracking on triangular prisms for finite-element-type models of variably-saturated flow

Consistent particle tracking relies on conforming velocity fields that ensure local mass conservation on elements. Cell-centered finite-volume and mixed finite-element methods result in conforming velocity fields but this is not the case for continuous Galerkin methods, such as the standard finite element method (FEM). Nonetheless standard FEM is often used for subsurface flow modeling because it yields a continuous approximation of hydraulic heads, and it naturally handles unstructured grids and full material tensors. Acknowledging these advantages and the wide-spread use of finite-element-type simulations, we present a postprocessing method that reconstructs a cell-centered finite-volume solution from a finite-element-type solution of the variably-saturated subsurface flow equation to obtain conforming, mass-conservative fluxes. Using the linear average velocity field derived from these fluid fluxes, we employ element-wise analytical solutions for triangular prisms to compute particle trajectories and associated travel times. As a result, we can compute consistent particle trajectories for variably-saturated flow solutions generated by node-centered methods, such as finite element or finite difference methods, that do not yield conforming velocity fields. Our flux reconstruction solves a linear elliptic problem whose size is on the order of the number of elements, which is computationally much faster than solving the initial, non-linear transient variably-saturated flow equation. Compared to other postprocessing schemes, our flux reconstruction is numerically stable, fast to compute, and does not induce severe numerical artifacts when applied to heterogeneous domains with strongly varying velocities. However, these advantages come with a comparably high coding effort and the necessity of solving a global system of equations. We show that the results of our flux reconstruction are close to the node-centered primal solution for variably saturated three-dimensional flow with heterogeneous material properties.

Development of an immersed boundary-multiphase lattice Boltzmann flux solver with high density ratio for contact line dynamics

Interaction between a two-phase fluid and a structure involving contact line dynamics is a common phenomenon. In this paper, we aim to develop a fluid–solid coupling model that can study contact line dynamics in the case of a high density ratio between the two fluids. The fluids are treated using a multiphase lattice Boltzmann flux solver (MLBFS) that uses the cell-centered finite volume method to obtain macroscopic flow variables, and the interface fluxes are reconstructed locally by the standard lattice Boltzmann method (LBM) solutions. This approach retains the advantages of the original LBM while being more flexible in handling nonuniform grids and external force terms. The immersed boundary method (IBM) is an effective method for processing structural information, and here, the implicit boundary-condition-enforced IBM is used to accurately satisfy the Dirichlet boundary condition (no-slip boundary). Moreover, the Neumann boundary condition is deemed to represent the contribution from the structure boundary flux and is incorporated into the IB-MLBFS. The developed IB-MLBFS is verified by several test cases, including contact line motion of a two-phase fluid along a circular cylinder and droplet spreading on a flat plate, where both equilibrium results and dynamic process are correctly reproduced for different density ratios and wettability conditions. Furthermore, based on the IB-MLBFS established here, the contact line dynamics of a two-phase fluid between two square cylinders or two circular cylinders is studied. The effects of distance, structure size, and wettability on the interface state and the contact angle are studied in detail. The robustness of the proposed model is verified.

Simple multiple reference frame for high-order solution of hovering rotors with and without ground effect

In the present work, the aerodynamic performance of the Caradonna and Tung and S-76 in hover were investigated using a simplified concept of multiple reference frame (MRF) technique in the context of high-order Monotone Upstream Centred Scheme for Conservation Laws (MUSCL) cell-centred finite volume method. In the present methodology, the frame of reference is defined at the solver level by a simple user input avoiding the use of mesh interface to handle the intersections between frames of reference. The calculations were made for both out-of-ground-effect (OGE) and in-ground-effect (IGE) cases and compared with experimental data in terms of pressure distribution, tip-vortex trajectory, vorticity contours and integrated thrust and torque. The predictions were obtained for several ground distances and collective pitch angle at tip Mach number of 0.6 and 0.892.

Wind Farm Fault Detection by Monitoring Wind Speed in the Wake Region

A novel concept of wind farm fault detection by monitoring the wind speed in the wake region is proposed in this study. A wind energy dissipation model was coupled with a computational fluid dynamics solver to simulate the fluid field of a wind turbine array, and the wind velocity and direction in the simulation were exported for identifying wind turbine faults. The 3D steady Navier–Stokes equations were solved by using the cell center finite volume method with a second order upwind scheme and a k−ε turbulence model. In addition, the wind energy dissipation model, derived from energy balance and Betz’s law, was added to the Navier–Stokes equations’ source term. The simulation results indicate that the wind speed distribution in the wake region contains significant information regarding multiple wind turbine faults. A feature selection algorithm specifically designed for the analysis of wind flow was proposed to reduce the number of features. This algorithm proved to have better performance than fuzzy entropy measures and recursive feature elimination methods under a limited number of features. As a result, faults in the wind turbine array could be detected and identified by machine learning algorithms.

DEIM reduced order model constructed by hybrid snapshot simulation

This paper presents a methodology that combines the hybrid snapshot simulation (Bai and Wang in Int J Comput Methods 18:2050029, 2020) and the discrete empirical interpolation method (DEIM) to reduce the computational cost of constructing a DEIM-based nonlinear reduced order model (ROM) while preserving its accuracy. In distinct contrast to its traditional counterpart, the time span of the hybrid snapshot simulation is divided into multiple intervals, and a majority of the intervals are simulated by local ROMs, while a fraction by the full order model (FOM). To tackle the challenge associated with limited FOM data for ROM-DEIM construction, a new approach is proposed to reconstruct and enrich the snapshot data of nonlinear model terms using solution data in ROM intervals. A formulation and procedure based on the cell-centered finite volume method (FVM) scheme are also developed to take into account two types of nonlinearities in ROM-DEIM: the componentwise and the transport-related, non-componentwise. A global ROM-DEIM is produced immediately at the end of the snapshot simulation and can be used to accelerate online simulation for different scenarios. The proposed methodology demonstrates excellent computational performance. Specifically, the computational time of the snapshot simulation is reduced by \(\sim\) 50% without compromising ROM-DEIM accuracy (relative error less than 0.8% compared with FOM). The results prove feasibility of combining hybrid snapshot simulation and DEIM for robust and efficient ROM-DEIM development.

Nonlinear multigrid based on local spectral coarsening for heterogeneous diffusion problems

This work develops a nonlinear multigrid method for diffusion problems discretized by cell-centered finite volume methods on general unstructured grids. The multigrid hierarchy is constructed algebraically using aggregation of degrees of freedom and spectral decomposition of reference linear operators associated with the aggregates. For rapid convergence, it is important that the resulting coarse spaces have good approximation properties. In our approach, the approximation quality can be directly improved by including more spectral degrees of freedom in the coarsening process. Further, by exploiting local coarsening and a piecewise-constant approximation when evaluating the nonlinear component, the coarse level problems are assembled and solved without ever re-visiting the fine level, an essential element for multigrid algorithms to achieve optimal scalability. Numerical examples comparing relative performance of the proposed nonlinear multigrid solvers with standard single-level approaches—Picard’s and Newton’s methods—are presented. Results show that the proposed solver consistently outperforms the single-level methods, both in efficiency and robustness.

Development of a Massively Parallelized Fluid-Based Plasma Simulation Code With a Finite-Volume Method on an Unstructured Grid

A massively parallelized unstructured-grid plasma simulation code based on a multifluid plasma model has been developed and validated. The code uses a collocated cell-centered finite-volume method and is designed to simulate intermediate low-pressure and atmospheric-pressure (AP) plasma sources with complex machine geometries. One of the novel features of this code is that it is implemented in a highly flexible computational platform named ultrafast massively parallel processing (ultraMPP), which allows straightforward addition and integration of different partial differential equation (PDE) solvers in a self-consistent manner. As to the numerical methods, the Scharfetter–Gummel scheme is adopted to handle the drift-diffusion flux of electrons, whereas the Harten–Lax-van Leer (HLL)-type approximate Riemann solver is used to handle the convection terms of ion momentum equations. For discretization of the diffusion terms in an unstructured grid, the Taylor expansion is used to deal with the effects of nonorthogonality of the cells, and the cell-center gradient is calculated using a least-squares method. The simulation code with a local-field approximation (LFA) was validated for two cases of AP plasmas as well as a case of an argon capacitively coupled plasma generated in a Gaseous Electronic Conference (GEC) reference cell with local mean energy approximation (LMEA) at intermediate low pressure. The simulation results were found to be in good agreement with the previously published experimental and simulation data.

Higher order fully implicit cell-centered finite volume method for simulation of oil-water displacement in porous medium

In this paper a fully implicit cell-centered finite volume method with Higher-Order Upwind scheme (HOU) for convective flow approximation is presented for numerical solution of oil-water displacement problem. Main advantages of this scheme are the robustness, accuracy and rapid convergence of Newton technique, used for linearization obtained discrete equations. Using of HOU scheme allows receive second order approximation of flux in the regions of slow saturation changing and stable solution near to the shocks. Efficiency of proposed scheme is verified by numerical tests.

Scalable parallel finite volume lattice Boltzmann method for thermal incompressible flows on unstructured grids

A parallel cell-centered finite volume thermal lattice Boltzmann method (FV-TLBM) for two-dimensional unstructured grids is proposed to simulate thermal incompressible flows. The governing equations for the velocity and temperature fields are solved using the D2Q9 lattice model in the lattice Boltzmann method and the cell-centered finite volume method. In this paper, the advective fluxes in the governing equations are evaluated by a low-diffusion Roe scheme, and the gradients of the particle distribution functions are computed with the least-square approach. To simulate large-scale complex flows within a reasonable time, a parallel algorithm for the FV-TLBM on unstructured grids is devised. The present method is validated by numerical simulations of the natural convection in a square cavity and natural convection in a concentric annulus at different Rayleigh numbers. The results obtained by the proposed method agree well with previous studies. The parallel performance results show that the proposed algorithm has considerable scalability and the parallel efficiency is as high as 96.50% for a case with over 5 billion grid cells using 6000 processor cores.

Preconditioning methods for compressible flow CFD codes: Revisited

The present article aims to implement and investigate a different preconditioning method based in a three-dimensional in-house compressible CFD code that ensures the robustness and numerical stability to determine the flowfield considering low Mach number flow. The present preconditioning method involve two different methodologies developed by references [8, 32]. The CFD solver was developed to calculate the Euler and Navier-Stokes equations, numerically, for steady-state regime based on the cell-centered finite volume method (FVM) using Reynolds Averaged Navier-Stokes equations (RANS). The centered second-order scheme was used for the discretization of convective terms from momentum equations. The explicit second-order five-step Runge-Kutta scheme was employed for the time-marching procedure, using an implicit residual smoothing technique to enhance the numerical stability. A local preconditioning method was implemented due to its robustness in predicting low Mach number flows in a compressible CFD code environment. However, for low Mach number flows, near stagnation points, numerical perturbations were amplified generating a stiffness in the convergence rate, which provided an inaccurate solution and numerical stability degradation. Aiming to improve the preconditioning robustness, a flux function and a new limiter were applied to operate with the preconditioning technique based on a pressure sensor. Those corrections re-scale the eigenvectors and ensure the locality of the algorithm, which improves the numerical stability and guarantees the convergence for low-speed flows. The inviscid flow over a NACA 0012 airfoil for compressible and incompressible cases shown accurate and robust solutions. For a viscous flow over a flat plate case in the compressible and incompressible cases, the preconditioning technique purposed in this work supplied good numerical solution in agreement with the analytical solution.

A patient-specific numerical model of the ankle joint for the analysis of contact pressure distribution.

A patient-specific numerical model of the ankle joint has been developed using open-source software with realistic material properties that mimics the physiological movement of the foot during the stance phase of the gait cycle. The patient-specific ankle geometry has been segmented as a castellated surface using 3DSlicer from the computed tomography image scans of a subject with no congenital or acquired pathology; subsequently, the bones are smoothed, and cartilage is included as a uniform thickness extruded layer. A high-resolution Cartesian mesh has been generated using cfMesh. The material properties are assigned in the model based on the CT image Hounsfield intensities and compared to a sandwich-based material model. Gait data of the same subject was obtained and used to relatively position the tibia, talus, and calcaneus bones in the model. The stance phase of the gait cycle is simulated using a cell-centred finite-volume method implemented in open-source software OpenFOAM. The predicted peak contact pressures occur in the range of 4.85-5.53 MPa with average pressures in the range of 1.56-1.95 MPa, and the contact area ranges between 429 and 707.8 mm2 for the entire stance phase with the mid-stance phase predicting the maximum contact area. These predictions are in agreement with results from the literature. The effect of arthritis on the contact characteristics of the ankle joint has also been examined. A concentrated increase in pressure was predicted that could be manifested as pain, thereby leading to reduced motion in the ankle. The model, with continued development, has the capability to understand the effect of joint degradation and furthermore, could help provide a tool to predict the efficiency of therapeutic surgical procedures as well as guide the development of next generation ankle prostheses. The work would be made available in the University College Dublin depository (https://github.com/laxmimurali/anklejoint) as well as research gate once the article has been published.