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]


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.



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


∂ρ ∂Ui ∂ρ +ρ + Ui = 0. ∂t ∂xi ∂xi


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.


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


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




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)


= ρni −



n n Ui+ 1 − U i− 1 2




 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


n where the local Courant number is cni+ 1 = Ui+ 1 ∆t/h. 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)


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


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



γ(x, t) = and


1 + a tan b

x0 [x,t] 2




– ¨ « ™ Œ b 1 + a tan[x/2] bt 1 x0 (x, t) = 2 arctan tan arctan − − . a b 2 a



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.



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.


Fire Dynam-

NIST Special Pub. 1018-5,









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.


A comparison of numerical methods for solving the ...

Jun 12, 2007 - solution. We may conclude that the FDS scheme is second-order accurate, but .... −a2 − cos bt + 2 arctan[γ(x, t)] + bsin bt + 2 arctan[γ(x, t)]. − ln.

76KB Sizes 0 Downloads 99 Views

Recommend Documents

a comparison of methods for determining the molecular ...
2009). The spatial and mass resolution of the cosmological simulation is matched exactly to the spatial and mass resolution of the simulation with the fixed ISM. Due to computational expense, the fully self- consistent cosmological simulation was not

A Comparison of Methods for Finding Steady- State ...
a converter is described as a linear time-invariant system. (LTI). ... resonant converters the procedure to find the solution is not as easy as it may .... New set of Ai.

Numerical Treatment for Solving Linear.pdf
Faculty of Science and Science Education at the University of Sulaimani, as partial. fulfillment of the requirements for the degree of Master of Science in.

A straightforward application of Monte Carlo (MC) method leads to the .... obviously not the case in the PC approach that requires the development of a modified.

Intrinsic Methods for Comparison of Corpora - raslan 2013
Dec 6, 2013 - syntactic analysis. ... large and good-enough-quality textual content (e.g. news ... program code, SEO keywords, link spam, generated texts,. . . ).

Comparison of Camera Motion Estimation Methods for ...
Items 1 - 8 - 2 Post-Doctoral Researcher, Construction Information Technology Group, Georgia Institute of. Technology ... field of civil engineering over the years.

Comparison of Training Methods for Deep Neural ... - Patrick GLAUNER
Attracted major IT companies including Google, Facebook, Microsoft and Baidu to make ..... Retrieved: April 22, 2015. The Analytics Store: Deep Learning.

Comparison of Four Methods for Determining ...
limit the access of lysostaphin to the pentaglycine bridge in the ... third glycines to the pentaglycine cross bridge. .... Lysostaphin: enzymatic mode of action.

A comparison of thinning methods in red pine
Abstract: Long-term replicated experiments that contrast thinning method (dominant thinning, .... ameter, using data from two long-term replicated studies in.

A numerical method for the computation of the ...
Considering the equation (1) in its integral form and using the condition (11) we obtain that sup ..... Stud, 17 Princeton University Press, Princeton, NJ, 1949. ... [6] T. Barker, R. Bowles and W. Williams, Development and application of a.

Evaluation of two numerical methods to measure the ...
singular measures, fractal dimension and fractional analysis play an impor- tant role in the ... d1, and the coefficient M that defines the size of the exploratory squares, ... possible approach, for example, is to use self-similar constructions). Mo