A High-Accuracy Finite-Difference Technique for Treating the Convection-Diffusion Equation

1975 ◽  
Vol 15 (06) ◽  
pp. 517-531 ◽  
Author(s):  
D.D. Laumbach

Abstract The convection-diffusion (C-D) equation arises from the conservation equations used in thermal recovery and miscible flooding. When convection predominates, this equation is very difficult to predominates, this equation is very difficult to represent numerically. The difficulty arises due to the hyperbolic character assumed as the Peclet number becomes large; consequently, the method presented here aims at providing the highest order presented here aims at providing the highest order of accuracy within this limit. The rationale underlying the treatment is to cancel a portion of the error in the convection term with that in the accumulation term. Thus, the technique presented is referred to as the truncation cancellation procedure (TCP). procedure (TCP). The application of this technique results in a new finite-difference representation of the C-D equation that is correct to the fourth order when the Peclet number approaches infinity. In this limit, when the dimensionless time step equals the dimensionless spatial increment, the discretization is exact. For very small time steps the method reduces to one previously considered to be one of the best semi-implicit discretizations of the C-D equation. For larger time steps, it yields significantly better results. The technique is applied to a linear form of the C-D equation in an equal-size grid block system; however, it is expected that the method should give favorable results in nonlinear and variable grid systems.The method is a semi-implicit one based on a three-point spatial and two-level time approximation. Thus, in one dimension, a set of difference equations is obtained that can be treated by solving a simple tridiagonal matrix. The extension of the TCP to two dimensions using a five-spatial-point star is demonstrated through an alternating-direction implicit solution method. The numerical results presented show that the TCP discretization yields presented show that the TCP discretization yields an excellent approximation in both Cartesian and radial coordinates. Comparisons with exact analytic solutions, conventional numerical techniques, and other high-accuracy numerical methods attest to the method's superiority over other formulations based on two time levels and three spatial locations. Introduction The purpose of this paper is to present a high-accuracy, finite-difference formulation for the convection-diffusion (C-D) equation. The discretization method presented is applied to a linear form of the C-D equation, with the expectation that the demonstrated improvements over existing methods will apply also to nonlinear forms, particularly those that are weakly nonlinear. The method is developed for equal grid spacing; however, it should give favorable results for variable grid spacing provided the grid size does not change too abruptly. provided the grid size does not change too abruptly. The difficulties associated with obtaining sufficient accuracy from conventional numerical representations of this equation were outlined quite clearly by Peaceman and Rachford. Improved methods for treating the C-D equation have been presented by Stone and Brian, Garder et al., presented by Stone and Brian, Garder et al., Price et al., Lantz, and Chaudhari. Price et al., Lantz, and Chaudhari. A significant improvement in the numerical treatment of the C-D equation was achieved by Stone and Brian. Unfortunately, the extension of this method to two dimensions has not been clear. Also, their method still possesses some oscillatory behavior in the vicinity of large gradients in the dependent variable when convection is strongly predominant This is particularly true for large time predominant This is particularly true for large time steps. The method proposed by Garder et al. uses the method of characteristics. Here the diffusion calculation is, in effect, an explicit one that consequently imposes a stability time-step limitation. Also, calculations are somewhat complicated by the moving points that must be tracked. Price et al. motivated a number of high-accuracy discretizations through the use of Galerkin's method. The discretization that they obtained through the use of chapeau basis functions is shown below to coincide with the recommended discretization of Stone and Brian. Consequently, this difference form suffers from the same oscillatory behavior in the vicinity of sharp fronts for large time steps. SPEJ P. 517

Geophysics ◽  
2021 ◽  
pp. 1-76
Author(s):  
Chunli Zhang ◽  
Wei Zhang

The finite-difference method (FDM) is one of the most popular numerical methods to simulate seismic wave propagation in complex velocity models. If a uniform grid is applied in the FDM for heterogeneous models, the grid spacing is determined by the global minimum velocity to suppress dispersion and dissipation errors in the numerical scheme, resulting in spatial oversampling in higher-velocity zones. Then, the small grid spacing dictates a small time step due to the stability condition of explicit numerical schemes. The spatial oversampling and reduced time step will cause unnecessarily inefficient use of memory and computational resources in simulations for strongly heterogeneous media. To overcome this problem, we propose to use the adaptive mesh refinement (AMR) technique in the FDM to flexibly adjust the grid spacing following velocity variations. AMR is rarely utilized in acoustic wave simulations with the FDM due to the increased complexity of implementation, including its data management, grid generation and computational load balancing on high-performance computing platforms. We implement AMR for 2D acoustic wave simulation in strongly heterogeneous media based on the patch approach with the FDM. The AMR grid can be automatically generated for given velocity models. To simplify the implementation, we employ a well-developed AMR framework, AMReX, to carry out the complex grid management. Numerical tests demonstrate the stability, accuracy level and efficiency of the AMR scheme. The computation time is approximately proportional to the number of grid points, and the overhead due to the wavefield exchange and data structure is small.


1970 ◽  
Vol 10 (04) ◽  
pp. 418-424 ◽  
Author(s):  
J.P. Letkeman ◽  
R.L. Ridings

Abstract The numerical simulation of coning behavior bas been one of the most difficult applications of numerical analysis techniques. Coning simulations have generally exhibited severe saturation instabilities in the vicinity of the well unless time-step sizes were severely restricted. The instabilities were a result of using mobilities based on saturations existing at the beginning of the time step. The time-step size limitation, usually the order of a few minutes, resulted in an excessive amount of computer time required to simulate coning behavior. This paper presents a numerical coning model that exhibits stable saturation and production behavior during cone formation and after breakthrough. Time-step sizes a factor of 100 to 1,000 times as large as those previously possible may be used in the simulation. To ensure stability, both production rates and mobilities are extrapolated production rates and mobilities are extrapolated implicitly to the new time level. The finite-difference equations used in the model are presented together with the technique for incorporating the updated mobilities and rates. Example calculations which indicate the magnitude of the time-truncation errors are included. Various factors which affect coning behavior are discussed. Introduction The usual formulation of numerical simulation models for multiphase flow involves the evaluation of flow coefficient terms at the beginning of a time step and assumes that these terms do not change over the time step. These assumptions are valid only if the values of pressure and saturation in the system do not change significantly over the time step. The design of a finite-difference model to evaluate coning behavior of gas or water in a single well usually results in a model which uses radial coordinates. A two-dimensional single-well model is illustrated in Fig. 1. This type of model will often produce finite-difference blocks with pore volumes less than 1 bbl near the wellbore while producing large blocks with pore volumes greater producing large blocks with pore volumes greater than 1 million bbl near the external radius. If one chooses to use a reasonable time-step size of, say, 1 to 10 days, then normal well rates would result in a flow of several hundred pore volumes per time step through blocks near the wellbore. Therefore the assumption that saturations remain constant, for the purpose of coefficient evaluation, is not valid. Welge and Weber presented a paper on water coning which recognized the limitation of using explicit coefficients and applied an arbitrary limitation on the maximum saturation change over a time step. While this method is workable for a certain class of problems, it is not rigorous and is not generally applicable. In 1968, Coats proposed a method to solve the gas percolation problem which is similar in that it also results from explicit mobilities. This proposal involved adjusting the relative permeability to gas at the beginning of the time step so that an individual block would not be over-depleted of gas during a time step. This method is not conveniently extended to two dimensions nor to coning problems where a block is voided many times during a time step. Blair and Weinaug explored the problems resulting from explicitly determined coefficients and formulated a coning model with implicit mobilities and a solution technique utilizing Newtonian iteration. While this method is rigorous, achieving convergence on certain problems is difficult and, in many cases, time-step size is still severely restricted. In addition to the problems resulting from explicit flow-equation coefficients in coning models, the specification of rates requires attention to ensure that the saturations remain stable in the vicinity of the producing block. SPEJ P. 418


2017 ◽  
Vol 2017 ◽  
pp. 1-14 ◽  
Author(s):  
Minjie Xu ◽  
Kai Fu ◽  
Xianqing Lv

We propose combining the adjoint assimilation method with characteristic finite difference scheme (CFD) to solve the aerosol transport problems, which can predict the distribution of atmospheric aerosols efficiently by using large time steps. Firstly, the characteristic finite difference scheme (CFD) is tested to compute the Gaussian hump using large time step sizes and is compared with the first-order upwind scheme (US1) using small time steps; the US1 method gets E2 error of 0.2887 using Δt=1/450, while CFD method gets a much smaller E2 of 0.2280 using a much larger time step Δt=1/45. Then, the initial distribution of PM2.5 concentration is inverted by the adjoint assimilation method with CFD and US1. The adjoint assimilation method with CFD gets better accuracy than adjoint assimilation method with US1 while adjoint assimilation method with CFD costs much less computational time. Further, a real case of PM2.5 concentration distribution in China during the APEC 2014 is simulated by using adjoint assimilation method with CFD. The simulation results are in good agreement with the observed values. The adjoint assimilation method with CFD can solve large scale aerosol transport problem efficiently.


Geophysics ◽  
1976 ◽  
Vol 41 (1) ◽  
pp. 62-78 ◽  
Author(s):  
Irshad R. Mufti

Resistivity surveying is commonly done by using a point‐source dipole. Consequently, a finite‐difference evaluation of apparent resistivity curves implies the use of three‐dimensional simulation models which necessitate prohibitive computer costs. However, if we assume variation of resistivity only in two dimensions and use a line‐source dipole for setting up the finite‐difference model of a given structure, the potential field can be evaluated easily. A discrete version of the resistivity problem in two dimensions, which takes into account nonuniform grid spacing, is presented as a system of self‐adjoint difference equations. Since the iterative solution of such a system does not require grid spacing to be less than a certain critical value, it was successfully used for the development of fast‐convergence finite‐difference models. By examining in detail the characteristics of the matrix associated with the evaluation of the potential field, it is demonstrated that the proposed modeling procedure will remain stable for all conceivable geometries and resistivity distributions. It was used for the investigation of certain models for which the corresponding results could also be computed analytically. A direct superposition of results obtained in the two cases shows that they are virtually identical. By making use of the reciprocity theorem, a computational short‐cut, which provides the evaluation of vertical sounding curves for a line‐source dipole in a single step, is put forward. Special problems related to the optimization of acceleration parameters as well as the estimation of the potential function along the subsurface boundaries of the model are discussed. It is concluded that by surrounding the model by a termination strip of very large effective width, either Neumann‐ or Dirichlet‐type boundary conditions can be used for simulating a semiinfinite medium without introducing signficant errors in the results.


Author(s):  
Mani Mehra ◽  
Kuldip Singh Patel ◽  
Ankita Shukla

AbstractIn this article, compact finite difference approximations for first and second derivatives on the non-uniform grid are discussed. The construction of diffusion wavelets using compact finite difference approximation is presented. Adaptive grids are obtained for non-smooth functions in one and two dimensions using diffusion wavelets. High-order accurate wavelet-optimized compact finite difference (WOCFD) method is developed to solve convection–diffusion equations in one and two dimensions on an adaptive grid. As an application in option pricing, the solution of Black–Scholes partial differential equation (PDE) for pricing barrier options is obtained using the proposed WOCFD method. Numerical illustrations are presented to explain the nature of adaptive grids for each case.


Sign in / Sign up

Export Citation Format

Share Document