A comparison of numerical methods for solving the continuity equation R. McDermott Building and Fire Research Laboratory National Institute of Standards and Technology Gaithersburg, MD 20899-8663, USA
[email protected]
Abstract
In this work we study two schemes for time advancement of the continuity equation. The first is a conservative, second-order Runge-Kutta method (RK2), and the second is the predictor-corrector scheme used by the Fire Dynamics Simulator (FDS). The results of both schemes are compared with an analytical solution. We may conclude that the FDS scheme is second-order accurate, but perhaps overly complex since it shows no improvement over the accuracy of the RK2 scheme for the test problem considered here.
1
Introduction
In this write-up we report the results of a verification test to confirm the order of accuracy for the numerical scheme used to time advance the continuity equation in the Fire Dynamics Simulator (FDS) [5]. We test the FDS scheme and a conservative, second-order Runge-Kutta scheme. The results of both schemes are compared to a nontrivial analytical solution. Let ρ(x, t) denote the fluid mass density and let Ui (x, t) denote the ith component of the fluid mass-average velocity. The continuity equation is ∂ρ ∂ρUi + = 0, ∂t ∂xi
(1)
∂ρ ∂Ui ∂ρ +ρ + Ui = 0. ∂t ∂xi ∂xi
(2)
or equivalently
NIST Technical Note 1489
12 June 2007
Below we first describe a conservative, second-order Runge-Kutta scheme (RK2) to solve the continuity equation. Then in Section 3 we describe the FDS predictor-corrector scheme. Next, we present an analytical solution for the continuity equation to be used in a convergence study to confirm the order of accuracy of the two schemes for a stationary velocity field (both should be second-order accurate in space and time). Results of the convergence study are presented in Section 5 and conclusions are drawn in Section 6.
2
Conservative RK2
Here we consider a 1D uniform grid with spacing h. Let xi = ih − h/2 denote the location of
n n the cell center. The staggered-grid velocity is located at the cell face, Ui+ 1 ≡ U(xi + h/2, t ). 2
The density is stored at the cell center, ρni ≡ ρ(xi , tn ). Given U(x, t) and the initial condition for the density ρ(x, 0), the first stage (denoted by the superscript (1)) of the conservative (divergence form) RK2 [2] time advancement of the continuity equation is (1)
ρi = ρni −
∆t n n n n ρi+ 1 Ui+ 1 − ρ 1 1U i− 2 i− 2 , 2 2 h
(3)
where the linear interpolation of the density is ρi+ 1 = 2
1 (ρi + ρi+1 ) , 2
on the “east” cell face, for example. The second (and final) stage is 1 n 1 (1) ∆t (1) n+1 (1) n+1 ρ − ρ ρ . ρn+1 = ρ + U − U 1 1 i i− 12 i− 12 2 i 2 i h i+ 2 i+ 2
3
(4)
(5)
FDS scheme
The scheme used by FDS is a modified version of MacCormack’s method (see, e.g., [3]) introduced by Denaro et al. [1], where the spatial differencing is upwind biased in the predictor step and downwind biased in the corrector step. Further, the FDS developers choose to expand the flux divergence before differencing. Hence, FDS discretizes (2) instead of (1). 2
The predictor step is given by (1)
ρi
= ρni −
∆tρni
n n Ui+ 1 − U i− 1 2
2
h
n n ρi+1 − ρni ρi − ρni−1 1 1 n n n n − ∆t {1 − ci+ 1 }Ui+ 1 + {1 + ci− 1 }Ui− 1 , 2 2 2 2 2 h 2 h
(6)
n where the local Courant number is cni+ 1 = Ui+ 1 ∆t/h. 2
2
The corrector step is 1 n (1) ρn+1 = ρ + ρ i i 2 i n+1 n+1 ∆t (1) Ui+ 12 − Ui− 12 − ρ 2 i h 2 3 2 3 (1) (1) (1) (1) ∆t 1 ρ − ρ 1 ρ − ρ i 5 i−1 5 n+1 4 i+1 n+1 4 i − {1 + cn+1 }Ui+ + {1 − cn+1 }Ui− . 1 1 i+ 12 i− 12 2 2 2 2 h 2 h (7)
4
Analytical solution
In this section we present a 1D analytical solution to (1) for a specified velocity field formed from the combination of a constant advection velocity a and a stationary compression wave sin(x) [4]. Let q ≡ ln ρ. With U(x) = a + sin(x), the 1D continuity equation can be written as ∂q ∂q + [a + sin(x)] + cos(x) = 0 . ∂t ∂x
(8)
A solution to (8) is q(x, t) = q(x0 [x, t], 0)
© + ln −a2 − cos bt + 2 arctan[γ(x, t)] + b sin bt + 2 arctan[γ(x, t)]
© − ln −a2 − cos 2 arctan[γ(x, t)] + b sin 2 arctan[γ(x, t)] , where b≡
√
−1 + a2 > 0 , 3
(9)
(10)
γ(x, t) = and
5
1 + a tan b
x0 [x,t] 2
,
(11)
¨ « b 1 + a tan[x/2] bt 1 x0 (x, t) = 2 arctan tan arctan − − . a b 2 a
(12)
Results
We consider a 1D periodic domain of length L = 2π. The domain is divided into N cells of width h = L/N. The center of the ith cell is xi . The analytical solution is periodic in time with period τ = 2π/b. The time step for the numerical solution is specified as ∆t = 0.5h/a, with a = 2. The initial density field is set to ρ(x, 0) = 1. The numerical solution is advanced in time to t ≈ τ /2. The number of time steps required is n = round(0.5τ /∆t). The total simulation time is then n∆t. The results of the numerical solution are compared with the analytical solution and the error, erms =
N 1 X [ρni − exp{q(xi , n∆t}]2 N i=1
!1/2 (13)
is reported in Figure 1 for a range of grid spacings. Both the conservative RK2 scheme and the FDS scheme are second-order accurate for this problem.
6
Conclusions
There are two conclusions that can be drawn from this study. The first is that the FDS scheme, which is loosely based on the scheme described in [1], is indeed second-order accurate for this problem which uses a constant velocity field. It is not clear whether the upwind/downwind biasing will affect the order of accuracy for a more realistic flow field. This is an issue worth studying, and should be considered in conjunction with the flux limiters used to ensure scalar realizability. The second conclusion (or least observation) that can be made is that the biasing in the FDS scheme neither degrades nor improves the accuracy for this test case. It is therefore 4
questionable as to whether there is any benefit in the added complexity of the scheme, especially when the scheme is not guaranteed to be conservative.
Acknowledgements This research was performed while the author held a National Research Council Research Associateship Award at the National Institute of Standards and Technology.
References [1] F. M. Denaro, F. S. Marra, and G. Continillo. Upwind treatment of source terms in the numerical simulation of reactive flows. In 16th International Colloquium on the Dynamics of Explosions and Reactive Systems, Cracow, Poland, August 3-8, 1997. [2] S. Gottlieb, C. W. Shu, and E. Tadmor. Strong stability-preserving high-order time discretization methods. SIAM Review, 43(1):89–112, 2001. [3] Harvard Lomax, Thomas H. Pulliam, and David W. Zingg. Fundamentals of Computational Fluid Dynamics. Springer, 2001. [4] R. McDermott. Analytical solutions to the continuity equation. NIST Technical Note 1488, 2007. [5] K. McGrattan, S. Hostikka, J. Floyd, H. Baum, and R. Rehm. ics Simulator (Version 5) Technical Reference Guide. http://fire.nist.gov/fds/, 2007.
5
Fire Dynam-
NIST Special Pub. 1018-5,
erms
−2
10
−3
10
−1
h
10
Figure 1: Continuity error erms given by (13) for a range of grid spacings h = L/N where N = {16, 32, 64, 128} for the test problem defined by (8). The analytical solution for the density field is given by (9). The initial density field is ρ(x, 0) = 1. The time step for the numerical solutions is specified as ∆t = 0.5h/a with a = 2. Thus, the grid spacing and time step are decreased proportionally. The dashed line represents O(h) and the solid line
represents O(h2 ). The results for the conservative RK2 scheme are denoted by the stars and the results for the FDS scheme are denoted by the squares. Both schemes are second-order accurate for this test problem which utilizes a constant velocity field.
6
