Finite Volume Discretization Schemes Research Articles

A Kinetic Finite Volume Discretization of the Multidimensional PIDE Model for Gene Regulatory Networks

In this paper, a finite volume discretization scheme for partial integro-differential equations (PIDEs) describing the temporal evolution of protein distribution in gene regulatory networks is proposed. It is shown that the obtained set of ODEs can be formally represented as a compartmental kinetic system with a strongly connected reaction graph. This allows the application of the theory of nonnegative and compartmental systems for the qualitative analysis of the approximating dynamics. In this framework, it is straightforward to show the existence, uniqueness and stability of equilibria. Moreover, the computation of the stationary probability distribution can be traced back to the solution of linear equations. The discretization scheme is presented for one and multiple dimensional models separately. Illustrative computational examples show the precision of the approach, and good agreement with previous results in the literature.

The Traffic Reaction Model: A kinetic compartmental approach to road traffic modeling

In this work, a family of finite volume discretization schemes for LWR-type first order traffic flow models (with possible on- and off-ramps) is proposed: the Traffic Reaction Model (TRM). These schemes yield systems of ODEs that are formally equivalent to the kinetic systems used to model chemical reaction networks. An in-depth numerical analysis of the TRM is performed. On the one hand, the analytical properties of the scheme (nonnegative, conservative, capacity-preserving, monotone) and its relation to more traditional schemes for traffic flow models (Godunov, CTM) are presented. Finally, the link between the TRM and kinetic systems is exploited to offer a novel compartmental interpretation of traffic models. In particular, kinetic theory is used to derive dynamical properties (namely persistence and Lyapunov stability) of the TRM for a specific road configuration. Two extensions of the proposed model, to networks and changing driving conditions, are also described.

Assessment of a Hybrid Eulerian–Lagrangian CFD Solver for Wind Turbine Applications and Comparison with the New MEXICO Experiment

In this paper, the hybrid Lagrangian–Eulerian solver HoPFlow is presented and evaluated against wind tunnel measurements from the New MEXICO experiment. In the paper, the distinct solvers that assemble the HoPFlow solver are presented, alongside with details on their mutual coupling and interaction. The Eulerian solver, MaPFlow, solves the compressible Navier–Stokes equations under a cell-centered finite-volume discretization scheme, while the Lagrangian solver uses numerical particles that carry mass, pressure, dilatation and vorticity as flow markers in order to represent the flow-field by following their trajectories. The velocity field is calculated with the use of the decomposition theorem introduced by Helmholtz. Computational performance is enhanced by utilizing the particle mesh (PM) methodology in order to solve the Poisson equations for the scalar potential ϕ and the stream function ψ→. The hybrid solver is tested in 3-D unsteady simulations concerning the axial flow around the wind turbine (WT) model rotor tested in the New MEXICO experimental campaign. Simulation results are presented as integrated rotor loads, radial distribution of aerodynamic forces and moments and pressure distributions at various span-wise positions along the rotor blades. Comparison is made against experimental data and computational results produced by the pure Eulerian solver. A total of 5 PM nodes per chord length of the blade section at 75% have been found to be sufficient to predict the loading at the tip region of the blade with great accuracy. Discrepancies with respect to measurements, observed at the root and middle sections of the blade, are attributed to the omission of the spinner geometry in the simulations.

A physics-informed deep convolutional neural network for simulating and predicting transient Darcy flows in heterogeneous reservoirs without labeled data

The physics-informed neural network (PINN) is a general deep learning framework for simulating physical processes and surrogate modeling without labeled data. The basic idea is to formulate the loss function according to the governing PDEs such that the neural network (NN) can be trained to minimize the PDE residual along with other misfits such as initial and boundary conditions. Following PINN, various networks have been developed for simulating steady and transient flows with or without labeled data. However, according to literature review, it is still not clear how to use NNs to simulate transient Darcy flows in highly heterogeneous reservoir models with source/sink terms in the absence of labeled data. In the current study, a physics-informed deep convolutional neural network (PIDCNN) architecture for simulating and predicting such flows is presented. Convolutional neural network is found to be more efficient than fully-connected neural network since 2D variables can be regarded as images. The finite volume discretization scheme is adopted to build the loss function to approximate the PDE residual such that flux continuity between neighboring cells of different properties can be implemented conveniently using the two-point flux approximation. Test cases are used to show that PIDCNN can accurately simulate transient Darcy flows in homogeneous and heterogeneous reservoirs. Further, it is demonstrated that PIDCNN can be trained as a surrogate to predict the transient flow fields of reservoir models not included in training. In addition, the CNN structure in the current study can be trained as a surrogate with labels for a particular output for better accuracy. A workflow is presented to demonstrate that CNN can be trained as an accurate surrogate for production rates using labels generated by the PIDCNN-based solver such that the entire workflow is external-label-free.

Smooth implicit hybrid upwinding for compositional multiphase flow in porous media

In subsurface multiphase flow simulations, poor nonlinear solver performance is a significant runtime sink. The system of fully implicit mass balance equations is highly nonlinear and often difficult to solve for the nonlinear solver, generally Newton(-Raphson). Although the physical problem is inherently nonlinear, discretization schemes introduce several strong nonlinearities that can cause Newton iterations to diverge or converge very slowly. This frequently results in time step cuts, leading to wasted iterations and computationally expensive simulations. Much literature has looked into how to improve the nonlinear solver through enhancements or safeguarding updates. In this work, we take a different approach; we aim to improve convergence with a smoother finite volume discretization scheme which is more suitable for the Newton solver.Building on recent work, we propose a novel total velocity hybrid upwinding scheme with weighted average flow mobilities (WA-HU TV) that is unconditionally monotone and extends to compositional multiphase simulations. Analyzing the solution space of a one-cell problem, we demonstrate the improved properties of the scheme and explain how it leverages the advantages of both phase potential upwinding and arithmetic averaging. This results in a flow subproblem that is smooth with respect to changes in the sign of phase fluxes, and is well-behaved when phase velocities are large or when co-current viscous forces dominate. Additionally, we propose a WA-HU scheme with a total mass (WA-HU TM) formulation that includes phase densities in the weighted averaging.The proposed WA-HU TV consistently outperforms existing schemes, yielding benefits from 5% to over 50% reduction in nonlinear iterations. The WA-HU TM scheme also shows promising results; in some cases leading to even more efficiency. However, WA-HU TM can occasionally also lead to convergence issues. Overall, based on the current results, we recommend the adoption of the WA-HU TV scheme as it is highly efficient and robust.

Consistency and convergence for a family of finite volume discretizations of the Fokker–Planck operator

We introduce a family of various finite volume discretization schemes for the Fokker–Planck operator, which are characterized by different Stolarsky weight functions on the edges. This family particularly includes the well-established Scharfetter–Gummel discretization as well as the recently developed square-root approximation (SQRA) scheme. We motivate this family of discretizations both from the numerical and the modeling point of view and provide a uniform consistency and error analysis. Our main results state that the convergence order primarily depends on the quality of the mesh and in second place on the choice of the Stolarsky weights. We show that the Scharfetter–Gummel scheme has the analytically best convergence properties but also that there exists a whole branch of Stolarsky means with the same convergence quality. We show by numerical experiments that for small convection the choice of the optimal representative of the discretization family is highly non-trivial, while for large gradients the Scharfetter–Gummel scheme stands out compared to the others.

Voxelized simulation of cerebral oxygen perfusion elucidates hypoxia in aged mouse cortex.

Departures of normal blood flow and metabolite distribution from the cerebral microvasculature into neuronal tissue have been implicated with age-related neurodegeneration. Mathematical models informed by spatially and temporally distributed neuroimage data are becoming instrumental for reconstructing a coherent picture of normal and pathological oxygen delivery throughout the brain. Unfortunately, current mathematical models of cerebral blood flow and oxygen exchange become excessively large in size. They further suffer from boundary effects due to incomplete or physiologically inaccurate computational domains, numerical instabilities due to enormous length scale differences, and convergence problems associated with condition number deterioration at fine mesh resolutions. Our proposed simple finite volume discretization scheme for blood and oxygen microperfusion simulations does not require expensive mesh generation leading to the critical benefit that it drastically reduces matrix size and bandwidth of the coupled oxygen transfer problem. The compact problem formulation yields rapid and stable convergence. Moreover, boundary effects can effectively be suppressed by generating very large replica of the cortical microcirculation in silico using an image-based cerebrovascular network synthesis algorithm, so that boundaries of the perfusion simulations are far removed from the regions of interest. Massive simulations over sizeable portions of the cortex with feature resolution down to the micron scale become tractable with even modest computer resources. The feasibility and accuracy of the novel method is demonstrated and validated with in vivo oxygen perfusion data in cohorts of young and aged mice. Our oxygen exchange simulations quantify steep gradients near penetrating blood vessels and point towards pathological changes that might cause neurodegeneration in aged brains. This research aims to explain mechanistic interactions between anatomical structures and how they might change in diseases or with age. Rigorous quantification of age-related changes is of significant interest because it might aide in the search for imaging biomarkers for dementia and Alzheimer’s disease.

Enhanced multiscale restriction-smoothed basis (MsRSB) preconditioning with applications to porous media flow and geomechanics

A novel method to enable application of the Multiscale Restricted Smoothed Basis (MsRSB) method to non M-matrices is presented. The original MsRSB method is enhanced with a filtering strategy enforcing M-matrix properties to enable the robust application of MsRSB as a preconditioner. Through applications to porous media flow and linear elastic geomechanics, the method is proven to be effective for scalar and vector problems with multipoint finite volume (FV) and finite element (FE) discretization schemes, respectively. Realistic complex (un)structured two- and three-dimensional test cases are considered to illustrate the method's performance.

Verification and Validation of openInjMoldSim, an Open-Source Solver to Model the Filling Stage of Thermoplastic Injection Molding

In the present study, the simulation of the three-dimensional (3D) non-isothermal, non-Newtonian fluid flow of polymer melts is investigated. In particular, the filling stage of thermoplastic injection molding is numerically studied with a solver implemented in the open-source computational library O p e n F O A M ® . The numerical method is based on a compressible two-phase flow model, developed following a cell-centered unstructured finite volume discretization scheme, combined with a volume-of-fluid (VOF) technique for the interface capturing. Additionally, the Cross-WLF (Williams–Landel–Ferry) model is used to characterize the rheological behavior of the polymer melts, and the modified Tait equation is used as the equation of state. To verify the numerical implementation, the code predictions are first compared with analytical solutions, for a Newtonian fluid flowing through a cylindrical channel. Subsequently, the melt filling process of a non-Newtonian fluid (Cross-WLF) in a rectangular cavity with a cylindrical insert and in a tensile test specimen are studied. The predicted melt flow front interface and fields (pressure, velocity, and temperature) contours are found to be in good agreement with the reference solutions, obtained with the proprietary software M o l d e x 3 D ® . Additionally, the computational effort, measured by the elapsed wall-time of the simulations, is analyzed for both the open-source and proprietary software, and both are found to be similar for the same level of accuracy, when the parallelization capabilities of O p e n F O A M ® are employed.

Hybrid finite volume and Monte Carlo method for solving multi-dimensional population balance equations in crystallization processes

A novel hybrid framework that combines our recently proposed event driven constant number Monte Carlo (MC) method and a finite volume discretization scheme is developed to solve population balance equation (PBE) involving nucleation, growth, and breakage events for a batch crystallization process. At each time step, the nucleation and growth terms of the PBE are solved using high resolution finite volume method (HRFV), and then the breakage term is solved using the MC approach. This integrated approach is especially advantageous for solving multi-dimensional PBE involving nucleation, growth, and breakage events. The efficiency of the proposed approach is demonstrated by solving both one-dimensional (1D) and two-dimensional (2D) PBE from published literature. Also, a set of seeded batch cooling crystallization of plate-like pyrazinamide crystals were performed in presence of ultrasound (US) and the proposed hybrid solution scheme was used successfully to model and simulate this system characterised by 2D PBE.

Nonlinear diffusion, boundary layers and nonsmoothness: Analysis of challenges in drift–diffusion semiconductor simulations

We study different discretizations of the van Roosbroeck system for charge transport in bulk semiconductor devices that can handle nonlinear diffusion. Three common challenges corrupting the precision of numerical solutions will be discussed: boundary layers, discontinuities in the doping profile, and corner singularities in L-shaped domains. We analyze and benchmark the error and the convergence order of finite difference, finite-element as well as advanced Scharfetter–Gummel type finite-volume discretization schemes. The most problematic of these challenges are boundary layers in the quasi-Fermi potentials near ohmic contacts, which can have a drastic impact on the convergence order. Using a novel formal asymptotic expansion, our theoretical analysis reveals that these boundary layers are logarithmic and significantly shorter than the Debye length.

Computational modeling of magnetoconvection: Effects of discretization method, grid refinement and grid stretching

The problem of natural convection in a laterally heated three-dimensional cubic cavity under the action of an externally imposed magnetic field is revisited. Flows at the Rayleigh number Ra=106 and the Hartmann number Ha=100, and three different orientations of the magnetic field are considered. The problem is solved using two independent numerical methods based on the second order finite-volume discretization schemes on structured Cartesian grids. Convergence toward grid-independent results is examined versus the grid refinement and near-wall grid stretching. Converged benchmark-quality results are obtained. It is shown that for convection flows with a strong magnetic field a steep, sometimes extremely steep, stretching near some of the boundaries is needed. Three-dimensional patterns and integral properties of the converged flow fields are reported and discussed. It is shown that the strongest magnetic suppression is yielded by the field directed along the imposed temperature gradient. The horizontal magnetic field perpendicular to the imposed temperature gradient stabilizes the main convection roll and leads to a flow with higher kinetic energy and heat transfer rate than in the non-magnetic case. Applicability of the quasi-two-dimensional model to natural convection flows in a box is discussed.

A multigrid scheme for solving convection–diffusion-integral optimal control problems

The fast multigrid solution of an optimal control problem governed by a convection–diffusion partial-integro differential equation is investigated. This optimization problem considers a cost functional of tracking type and a constrained distributed control. The optimal control sought is characterized by the solution to the corresponding optimality system, which is approximated by a finite volume and quadrature discretization schemes and solved by multigrid techniques. The proposed multigrid approach combines a multigrid method for the governing model with a fast multigrid integration method. The convergence of this solution procedure is analyzed by local Fourier analysis and validated by results of numerical experiments.

Family of Runge–Kutta-Based Algorithms for Unsteady Incompressible Flows

A family of algorithms based on implicit Runge–Kutta time integration suitable for unsteady incompressible flow is developed. Spatially discretized momentum equations are integrated in time using implicit Runge–Kutta methods. An equation for pressure is formed by combining the discretized continuity equation with the discretized momentum equations. This family of algorithms, called IRK-SIMPLER, uses only exact discretized mass and momentum equations and requires no relaxation to converge. Two distinct IRK-SIMPLER variants are analyzed. The first variant solves the pressure equation once per time step and has the advantage of fewer computations per time step but is limited to first-order accuracy in time. The second variant centers on reformulating a pressure equation each Runge–Kutta stage by manipulating the momentum stage equations, and this pressure equation is solved multiple times within a time step. The second variant is shown to achieve a higher temporal order of accuracy. By analyzing the largest time step allowable and run time required to achieve time-accurate solutions, the most efficient algorithms in this family are delineated. For the chosen problems, IRK-SIMPLER is able to achieve as high as 101 times speedup compared to the traditional SIMPLER algorithm with Crank–Nicolson time integration using a finite volume discretization scheme.

Fast Iterative Solvers for Linear Systems Arising from Time-Dependent Space-Fractional Diffusion Equations

In this paper, we study the linear systems arising from the discretization of time-dependent space-fractional diffusion equations. By using a finite difference discretization scheme for the time derivative and a finite volume discretization scheme for the space-fractional derivative, Toeplitz-like linear systems are obtained. We propose using the approximate inverse-circulant preconditioner to deal with such Toeplitz-like matrices, and we show that the spectra of the corresponding preconditioned matrices are clustered around 1. Experimental results on time-dependent and space-fractional diffusion equations are presented to demonstrate that the preconditioned Krylov subspace methods converge very quickly.

Coupled fine-mesh neutronics and thermal-hydraulics – Modeling and implementation for PWR fuel assemblies

In this paper we present a fine-mesh solver aimed at resolving in a coupled manner and at the pin cell level the neutronic and thermal-hydraulic fields. Presently, the tool considers Pressurized Water Reactor (PWR) conditions. The methods and implementation strategy are such that the coupled neutronic and thermal-hydraulic problem is formulated in a fully three-dimensional (3D) and fine mesh manner, and for steady-state situations. The solver is built on finite volume discretization schemes, matrix solvers and capabilities for parallel computing that are available in the open source C++ library foam-extend-3.0. The angular neutron flux is determined with a multigroup discrete ordinates method (SN), solved by a sweeping algorithm. The thermal-hydraulics is based on Computational Fluid Dynamics (CFD) models for the moderator/coolant mass, momentum, and energy equations, together with the fuel pin energy equation. The multiphysics coupling is solved by making use of an iterative algorithm, and convergence is ensured for both the separate equations and the coupled scheme. Since all the equations are implemented in the same software, all fields can be directly accessed in such a manner that external transfer and external mapping are avoided. The parallelization relies on a domain decomposition which is shared between the neutronics and the thermal-hydraulics. The latter allows to exchange the coupled data locally on each CPU, thus minimizing the data transfer. The code is tested on a quarter of a 15×15 PWR fuel lattice. The results show that convergence is successfully reached, and correct physical behaviors of all fields can be achieved with a reasonable computational effort.

Simulation of Premixed Turbulent Combustion in a Gas Engine using the Immersed Boundary Method and a Level Set Flame Model

AbstractAn efficient simulation approach for turbulent flame brush propagation is a level set formulation closed by the turbulent flame speed. A formulation of the level set equation with the corresponding treatment of the turbulent mass burning rate that is compatible with standard Finite Volume discretization schemes available in computational fluid dynamics codes is employed. In order to simplify and to speed up the meshing process in complicated geometries (here in gas engines) the immersed boundary method in a continuous formulation, where the forces replacing the boundaries are introduced in the momentum conservation equations before discretization, is employed. In our contribution, aspects of the numerical implementation of the level set flame model combined with the immersed boundary formulation in OpenFOAM are presented. First representative simulation results of a homogeneous methane/air mixture combustion in a simplified engine geometry are shown. (© 2014 Wiley‐VCH Verlag GmbH & Co. KGaA, Weinheim)

Isotropic finite volume discretization

Finite volume methods traditionally employ dimension by dimension extension of the one-dimensional reconstruction and averaging procedures to achieve spatial discretization of the governing partial differential equations on a structured Cartesian mesh in multiple dimensions. This simple approach based on tensor product stencils introduces an undesirable grid orientation dependence in the computed solution. The resulting anisotropic errors lead to a disparity in the calculations that is most prominent between directions parallel and diagonal to the grid lines. In this work we develop isotropic finite volume discretization schemes which minimize such grid orientation effects in multidimensional calculations by eliminating the directional bias in the lowest order term in the truncation error. Explicit isotropic expressions that relate the cell face averaged line and surface integrals of a function and its derivatives to the given cell area and volume averages are derived in two and three dimensions, respectively. It is found that a family of isotropic approximations with a free parameter can be derived by combining isotropic schemes based on next-nearest and next-next-nearest neighbors in three dimensions. Use of these isotropic expressions alone in a standard finite volume framework, however, is found to be insufficient in enforcing rotational invariance when the flux vector is nonlinear and/or spatially non-uniform. The rotationally invariant terms which lead to a loss of isotropy in such cases are explicitly identified and recast in a differential form. Various forms of flux correction terms which allow for a full recovery of rotational invariance in the lowest order truncation error terms, while preserving the formal order of accuracy and discrete conservation of the original finite volume method, are developed. Numerical tests in two and three dimensions attest the superior directional attributes of the proposed isotropic finite volume method. Prominent anisotropic errors, such as spurious asymmetric distortions on a circular reaction–diffusion wave that feature in the conventional finite volume implementation are effectively suppressed through isotropic finite volume discretization. Furthermore, for a given spatial resolution, a striking improvement in the prediction of kinetic energy decay rate corresponding to a general two-dimensional incompressible flow field is observed with the use of an isotropic finite volume method instead of the conventional discretization.

Accuracy analysis of unstructured finite volume discretization schemes for diffusive fluxes

Numerical errors form a major source of uncertainty in CFD simulations. This source of error, which is more significant on unstructured meshes, can be alleviated by a careful choice of discretization scheme. In the context of finite volume methods, there are many suggestions for the discretization of diffusive fluxes appearing in viscous flow simulations. In this paper, a wide range of discretization schemes commonly used for diffusive flux approximation in cell-centered unstructured finite volume solvers are compared in terms of discretization and truncation errors. In addition, a novel eigenanalysis tool is proposed to relate these two forms of numerical error to each other and to interpret the error behaviors obtained by each scheme. The error comparisons are performed on unstructured meshes consisting of both isotropic and anisotropic triangles. Our results suggest that adding a solution jump term to the baseline face gradient (determined as the average of two adjacent cell gradients) provides the most accurate approximation for diffusive fluxes. Also, this term produces sufficient damping to stabilize the discretization even on highly skewed meshes.

Two-dimensional depth-averaged modelling of dam-break flows over mobile beds

The paper is aimed at the description and validation of a novel two-dimensional depth-averaged simulation tool for highly unsteady discontinuous flows over complex time-evolving geometries. The conceptual model is developed within the shallow-flow framework and features non-equilibrium sediment transport. A fully conservative finite-volume discretization scheme is employed. The treatment of bed-slope source terms ensures that the scheme is well balanced and leads to correct energy losses in discontinuities. Existing laboratory data are used to validate the model and to discuss some of its embedded formulations. Blind-test comparison shows that the model is capable of describing the main dam-break flow and bed morphology features. The solution is shown to be sensitive to the formulations for sediment transport capacity and for adaptation length (Λ). A better functional relation of Λ with the Shields parameter is investigated. The quality of the solution is shown to be only marginally improved by the inclusion of turbulent stresses.