Simulation of flow and transport in petroleum reservoirs involves solving coupled systems of advectiondiffusionreaction equations with nonlinear flux functions, diffusion coefficients, and reactions/wells. It is important to develop numerical schemes that can approximate all three processes at once, and to high order, so that the physics can be well resolved. In this paper, we propose an approach based on high order, finite volume, implicit, Weighted Essentially NonOscillatory (iWENO) schemes. The resulting schemes are locally mass conservative and, being implicit, suited to systems of advectiondiffusionreaction equations. Moreover, our approach gives unconditionally Lstable schemes for smooth solutions to the linear advectiondiffusionreaction equation in the sense of a von Neumann stability analysis. To illustrate our approach, we develop a third order iWENO scheme for the saturation equation of twophase flow in porous media in two space dimensions. The keys to high order accuracy are to use WENO reconstruction in space (which handles shocks and steep fronts) combined with a twostage RadauIIA RungeKutta time integrator. The saturation is approximated by its averages over the mesh elements at the current time level and at two future time levels; therefore, the scheme uses two unknowns per grid block per variable, independent of the spatial dimension. This makes the scheme fairly computationally efficient, both because reconstructions make use of local information that can fit in cache memory, and because the global system has about as small a number of degrees of freedom as possible. The scheme is relatively simple to implement, high order accurate, maintains local mass conservation, applies to general computational meshes, and appears to be robust. Preliminary computational tests show the potential of the scheme to handle advectiondiffusionreaction processes on meshes of quadrilateral gridblocks, and to do so to high order accuracy using relatively long time steps. The new scheme can be viewed as a generalization of standard cellcentered finite volume (or finite difference) methods. It achieves high order in both space and time, and it incorporates WENO slope limiting.
more »
« less
Positivitypreserving and energydissipative finite difference schemes for the Fokker–Planck and Keller–Segel equations
Abstract In this work we introduce semiimplicit or implicit finite difference schemes for the continuity equation with a gradient flow structure. Examples of such equations include the linear Fokker–Planck equation and the Keller–Segel equations. The two proposed schemes are firstorder accurate in time, explicitly solvable, and secondorder and fourthorder accurate in space, which are obtained via finite difference implementation of the classical continuous finite element method. The fully discrete schemes are proved to be positivity preserving and energy dissipative: the secondorder scheme can achieve so unconditionally while the fourthorder scheme only requires a mild time step and mesh size constraint. In particular, the fourthorder scheme is the first high order spatial discretization that can achieve both positivity and energy decay properties, which is suitable for long time simulation and to obtain accurate steady state solutions.
more »
« less
 Award ID(s):
 1913120
 NSFPAR ID:
 10329079
 Date Published:
 Journal Name:
 IMA Journal of Numerical Analysis
 ISSN:
 02724979
 Format(s):
 Medium: X
 Sponsoring Org:
 National Science Foundation
More Like this


The thermal radiative transfer (TRT) equations form an integrodifferential system that describes the propagation and collisional interactions of photons. Computing accurate and efficient numerical solutions TRT are challenging for several reasons, the first of which is that TRT is defined on a highdimensional phase space that includes the independent variables of time, space, and velocity. In order to reduce the dimensionality of the phase space, classical approaches such as the P$_N$ (spherical harmonics) or the S$_N$ (discrete ordinates) ansatz are often used in the literature. In this work, we introduce a novel approach: the hybrid discrete (H$^T_N$) approximation to the radiative thermal transfer equations. This approach acquires desirable properties of both P$_N$ and S$_N$, and indeed reduces to each of these approximations in various limits: H$^1_N$ $\equiv$ P$_N$ and H$^T_0$ $\equiv$ S$_T$. We prove that H$^T_N$ results in a system of hyperbolic partial differential equations for all $T\ge 1$ and $N\ge 0$. Another challenge in solving the TRT system is the inherent stiffness due to the large timescale separation between propagation and collisions, especially in the diffusive (i.e., highly collisional) regime. This stiffness challenge can be partially overcome via implicit time integration, although fully implicit methods may become computationally expensive due to the strong nonlinearity and system size. On the other hand, explicit timestepping schemes that are not also asymptoticpreserving in the highly collisional limit require resolving the meanfree path between collisions, making such schemes prohibitively expensive. In this work we develop a numerical method that is based on a nodal discontinuous Galerkin discretization in space, coupled with a semiimplicit discretization in time. In particular, we make use of a second order explicit RungeKutta scheme for the streaming term and an implicit Euler scheme for the material coupling term. Furthermore, in order to solve the material energy equation implicitly after each predictor and corrector step, we linearize the temperature term using a Taylor expansion; this avoids the need for an iterative procedure, and therefore improves efficiency. In order to reduce unphysical oscillation, we apply a slope limiter after each time step. Finally, we conduct several numerical experiments to verify the accuracy, efficiency, and robustness of the H$^T_N$ ansatz and the numerical discretizations.more » « less

Although the full form of the Rayleigh–Plesset (RP) equation more accurately depicts the bubble behavior in a cavitating flow than its reduced form, it finds much less application than the latter in the computational fluid dynamic (CFD) simulation due to its high stiffness. The traditional variable timestep scheme for the full form RP equation is difficult to be integrated with the CFD program since it requires a tiny time step at the singularity point for convergence and this step size may be incompatible with time marching of conservation equations. This paper presents two stable and efficient numerical solution schemes based on the finite difference method and Euler method so that the fullform RP equation can be better accepted by the CFD program. By employing a truncation bubble radius to approximate the minimum bubble size in the collapse stage, the proposed schemes solve for the bubble radius and wall velocity in an explicit way. The proposed solution schemes are more robust for a wide range of ambient pressure profiles than the traditional schemes and avoid excessive refinement on the time step at the singularity point. Since the proposed solution scheme can calculate the effects of the secondorder term, liquid viscosity, and surface tension on the bubble evolution, it provides a more accurate estimation of the wall velocity for the vaporization or condensation rate, which is widely used in the cavitation model in the CFD simulation. The legitimacy of the solution schemes is manifested by the agreement between the results from these schemes and established ones from the literature. The proposed solution schemes are more robust in face of a wide range of ambient pressure profiles.more » « less

Abstract We present a discontinuity aware quadrature (DAQ) rule and use it to develop implicit selfadaptive theta (SATh) schemes for the approximation of scalar hyperbolic conservation laws. Our SATh schemes require the solution of a system of two equations, one controlling the cell averages of the solution at the time levels, and the other controlling the spacetime averages of the solution. These quantities are used within the DAQ rule to approximate the time integral of the hyperbolic flux function accurately, even when the solution may be discontinuous somewhere over the time interval. The result is a finite volume scheme using the theta time stepping method, with theta defined implicitly (or selfadaptively). Two schemes are developed, selfadaptive theta upstream weighted (SAThup) for a monotone flux function using simple upstream stabilization, and selfadaptive theta Lax–Friedrichs (SAThLF) using the Lax–Friedrichs numerical flux. We prove that DAQ is accurate to second order when there is a discontinuity in the solution and third order when it is smooth. We prove that SAThup is unconditionally stable, provided that theta is set to be at least 1/2 (which means that SATh can be only first order accurate in general). We also prove that SAThup satisfies the maximum principle, and is total variation diminishing under appropriate monotonicity and boundary conditions. General flux functions require the SAThLF scheme, so we assess its accuracy through numerical examples in one and two space dimensions. These results suggest that SAThLF is also stable and satisfies the maximum principle (at least at reasonable CourantFriedrichsLewy numbers). Compared to the solutions of finite volume schemes using Crank–Nicolson and backward Euler time stepping, SAThLF solutions often approach the accuracy of the former, but without oscillation, and they are numerically less diffuse than the latter.more » « less

null (Ed.)In this paper we propose and analyze a finite difference numerical scheme for the PoissonNernstPlanck equation (PNP) system. To understand the energy structure of the PNP model, we make use of the Energetic Variational Approach (EnVarA), so that the PNP system could be reformulated as a nonconstant mobility H − 1 H^{1} gradient flow, with singular logarithmic energy potentials involved. To ensure the unique solvability and energy stability, the mobility function is explicitly treated, while both the logarithmic and the electric potential diffusion terms are treated implicitly, due to the convex nature of these two energy functional parts. The positivitypreserving property for both concentrations, n n and p p , is established at a theoretical level. This is based on the subtle fact that the singular nature of the logarithmic term around the value of 0 0 prevents the numerical solution reaching the singular value, so that the numerical scheme is always welldefined. In addition, an optimal rate convergence analysis is provided in this work, in which many highly nonstandard estimates have to be involved, due to the nonlinear parabolic coefficients. The higher order asymptotic expansion (up to third order temporal accuracy and fourth order spatial accuracy), the rough error estimate (to establish the ℓ ∞ \ell ^\infty bound for n n and p p ), and the refined error estimate have to be carried out to accomplish such a convergence result. In our knowledge, this work will be the first to combine the following three theoretical properties for a numerical scheme for the PNP system: (i) unique solvability and positivity, (ii) energy stability, and (iii) optimal rate convergence. A few numerical results are also presented in this article, which demonstrates the robustness of the proposed numerical scheme.more » « less