Notes on the Fire Dynamics Simulator (FDS) projection scheme R. McDermott Building and Fire Research Laboratory National Institute of Standards and Technology Gaithersburg, MD 20899-8663, USA
[email protected]
Abstract
The predictor-corrector method used for time advancement of the velocity field in the Fire Dynamics Simulator (FDS) is cast in the form of a projection method. We discuss two general classes of projection schemes – a divergence class and an advective class – and comment on their respective strengths and weaknesses. We then show that the FDS time advancement belongs to the advective class of projection methods.
1
Introduction
In this article we describe the Fire Dynamics Simulator (FDS) time integration scheme and show that this scheme can be viewed as a variable-density projection scheme. The material in this document is supplementary to the FDS Tech Guide [8] and the reader is encouraged to become familiar with that document before tackling the details presented here. The FDS5 scheme is second-order in time (for velocity) and enforces the velocity divergence constraint imposed by the thermal energy equation. This puts the scheme in the same class of projection methods as used by the Center for Computational Sciences and Engineering (CCSE) group at Lawrence Berkeley National Laboratory [1]. Similar to the CCSE scheme, but for different reasons, the FDS scheme does not enforce discrete conservation of momentum or kinetic energy. The CCSE scheme does not conserve momentum because NIST Technical Note 1485
22 August 2007
it is collocated. The FDS scheme does not conserve momentum because the Stokes form (as opposed to the divergence form) of the Navier-Stokes (NS) equations is discretized in space. Neither scheme discretely conserves kinetic energy, which requires a time-centered fully-implicit scheme [4, 7]. As will be shown, the FDS scheme is semi-implicit (explicit in velocity and implicit [and first-order] in pressure). In projection methods [2] the pressure is simply a scalar responsible for enforcing the constraints on the velocity field imposed by the continuity (mass conservation) equation. For variable-density flows, projection methods can be broken into two general classes: a divergence class and an advective class, based mainly on the form of the momentum and continuity equations which are solved. Let ρ be the fluid mass density. The ith component of velocity is ui ; the hydrodynamic pressure is p; τij is the viscous stress tensor∗ ; and gi is a body force per unit mass in the ith direction. Here and throughout this document summation is implied over repeated suffixes. However, bracketed suffixes and superscripts do not follow the summation convention. In the divergence class of projection schemes, the typical starting point is to write down the momentum and continuity equations in divergence form: ∂p ∂τij ∂ρui ∂ρui uj + =− − + ρgi , ∂t ∂xj ∂xi ∂xj
(1)
∂ρ ∂ρuj + = 0. ∂t ∂xj
(2)
In the advective class, the typical starting point is to write down the governing equations in advective form:
where ∗
Dui 1 ∂p 1 ∂τij =− − + gi , Dt ρ ∂xi ρ ∂xj
(3)
1 Dρ ∂uj =− , ρ Dt ∂xj
(4)
D( ) ∂( ) ∂( ) ≡ + uj is the material derivative. Dt ∂t ∂xj
For brevity we are only looking at the direct numerical simulation (DNS) version of the governing
equations, which is sufficient for the discussion of time integration. It should be noted that the large-eddy simulation (LES) form of the equations could easily be substituted here with no real change in the projection method.
2
While these sets of equations, {(1) and (2)} and {(3) and (4)}, are mathematically identical, the numerical schemes that result from their spatial and temporal discretizations can have drastically different properties. Generally, when using the divergence form, momentum is discretely conserved but the continuity equation is typically represented to first-order accuracy in time for a given stage. However, usually multiple stages are employed to increase the temporal order of accuracy. Also, with the divergence form it is difficult to directly incorporate the effects of thermal expansion on the velocity field without some sort of iterative correction. Hence, these schemes tend to be “loosely implicit” in that the iterations move the velocity field toward an implicit solution, but it is simply impractical to converge this solution to any reasonable tolerance. Examples of divergence-class projection methods can be found in the schemes used by Stanford [10], Darmstadt [5], and Utah [7]. The velocity field generated by advective-class projection schemes does not discretely conserve momentum. However, this is not an absolute requirement. In direct numerical simulation (DNS) and large-eddy simulation (LES), kinetic energy conservation is the key consideration and it is sufficient to conserve momentum to some order of accuracy (hopefully, second order). Unfortunately, no one has yet designed a “kinetic-energy preserving” projection scheme for variable-density flows and so most codes rely heavily on flux limiters and other sources of numerical viscosity for stability (which can seriously limit the order of accuracy). An advantage of the advective form is that the divergence form of the continuity equation is typically used to update the density with second-order time accuracy. The ideal gas law is used as the equation of state and is differentiated to obtain the left-hand side (LHS) of (4). The result requires (among other things) input from the thermal energy equation (or some other representation of the temperature; FDS obtains the temperature from the equation of state). Eq. (4) then is used as a divergence constraint for the velocity field. We will not go into the details of the divergence-class schemes. We will, however, look closely at the advective-class projection method used by FDS. It is worth noting that both classes can be written in the same general framework, which is why they are both considered
3
projection methods. To illustrate the framework, we consider a single-stage scheme, but the framework generalizes to multi-stage schemes, as we will see later. Here we give an outline of the framework. A given stage of any projection scheme can be derived as follows: 1. Write down the desired numerical update equation for the velocity (or momentum) assuming that the velocity field obeys the desired constraints. The hydrodynamic pressure on the right-hand side (RHS) of this equation (which we will denote P to distinguish it from the true pressure in the governing equations) has no real time location! Its only function is to enforce whatever constraints are to be applied to the velocity field. 2. Provide a guess for the pressure field. Note that P (x) = 0 for all x is a perfectly legitimate guess. (Note that, historically, setting P (x) = 0 resulted in “projection methods,” whereas using some other guess for the pressure resulted in “pressure correction methods.” There is really no good reason to distinguish between the two.) 3. Write down the velocity update equation using the guessed pressure. Call the resulting velocity U∗ (or somehow denote this field to indicate that it does not obey the desired mass-conservation constraints). 4. Subtract the uncorrected velocity (or momentum) equation from the desired equation. That is, subtract the result of item 3 from item 1. 5. Take the divergence of the result from item 4. This results in a linear equation for the pressure. It contains unknown terms that must be estimated: In the case of divergenceclass schemes, it contains the divergence of the mass flux at the new time location. Often, the time derivative of the density is substituted here via the continuity equation. However, this derivative is located at the new time location and any effort to estimate the derivative involves one-sided differences, unless an implicit scheme is employed. For advective-class schemes, the Poisson equation contains the divergence of the velocity at the new time location as an unknown. Here, (4) is used as an estimate. It is easier to couple the thermal constraints in time using this method. 4
6. Solve the Poisson equation for the pressure. 7. Using the pressure just obtained, project the uncorrected velocity onto the space of vector fields constrained by item 5. The type of scheme resulting from this framework is often referred to as a “semi-implicit” or “fractional step” method. Often (as is the case in FDS) the velocity is treated explicitly and the pressure is treated implicitly. However, it is also common to mix the time locations of the velocity and treat the advective terms explicitly and the diffusive terms implicitly [10]. The same general framework can be applied for deriving the projection method, but iteration within a given stage is required to reconcile the implicit treatment of both the diffusive terms and pressure. As mentioned previously, in practice it is common to only loosely converge the implicit solution (often by specifying the number of iterations [sometimes as few as two or three] instead of an error tolerance). The remainder of this document is organized as follows: In Section 2, we cast the FDS time integration scheme in terms of the above framework, without considering the Stokesform variation of the momentum equation. Then, in Section 3, we translate this into the notation of the FDS Technical Reference Guide [8]. The equivalence of the projection methods from the two sections is thus established. One should note that the overall “numerical method” also involves details related to the spatial discretization. These details are not considered here because they are irrelevant to the projection scheme. In Section 4 we summarize our findings.
2
A two-stage projection scheme
In this section we describe a two-stage projection method based on a second-order RungeKutta (RK2) time integration scheme [3]. We denote the cell-centered mass density by ρ, the ith component of the staggered-grid velocity by Ui , the cell-centered pressure by p, the viscous stress tensor by τij (off diagonal components are stored at vertices and diagonal components are cell-centered), and the cell-centered body force (taken to be uniform) in 5
the ith direction by gi . The semi-discrete form [6] of the advective form of the momentum equation can be written as
dUi 1 δp = − Fi + (xi ) , dt ρ δxi where
(5)
(x )
δU i j 1 δτij F i ≡ Uj + (xi ) + gi . δxj δxj ρ The overline, ( )
(6)
(xi )
, denotes an interpolation operation in the ith direction and the staggeredδ( ) . A detailed discussion of grid differential operator in the ith direction is denoted by δxi these operators can be found in [9]. Our goal is to integrate (5) from t(n) to t(n+1) ≡ t(n) + ∆t, under the constraint that the staggered-grid velocity obeys the discrete analog to (4) (The expression for the discrete divergence constraint is complex and presenting it here would distract from the focus of this work. The interested reader is referred to [1, 8].). The RK2 scheme (also known as “Modified Euler”) basically takes two full Forward Euler (FE) steps and then averages the current state with the final state. We denote the intermediate stages of the scheme using the superscripts (1) and (2). Thus, the FE stages take the solution through the sequence (n)
Ui
2.1
(1)
→ Ui
(2)
(n+1)
→ Ui , and the final solution is then given by Ui
=
1 2
(n)
Ui
(2)
+ Ui
.
Stage 1
We assume that the density has been updated prior to the velocity. A Forward Euler advancement is sufficient here,
(1)
ρ
(n)
=ρ
(n)
δ ρ(xj ),(n) Uj − ∆t δxj
.
(7)
Following the framework outlined in Section 1, the desired first stage of the scheme is a simple Forward Euler update, " (1) Ui
=
(n) Ui
− ∆t
(n) Fi
#
δp(1) + (xi ),(1) . δxi ρ 1
(8)
It is important to note that the pressure is not associated with a time location. The superscript is simply used to distinguish this value of the pressure from other stage values. 6
The next step in the framework is to write down the estimated velocity update with a guessed value for the pressure. We will use zero for the guess. Hence, we have (1),∗
Ui
(n)
= Ui
h
(n)
− ∆t Fi
i
.
(9)
The asterisk indicates that this velocity field does not satisfy the divergence constraint. We then subtract the estimated update from the desired update to obtain "
(1) Ui
=
(1),∗ Ui
#
δp(1) − ∆t (xi ),(1) . δxi ρ 1
(10) (1),∗
Now, given p(1) , Eq. (10) is the projection: it projects the velocity Ui
onto the space of
vector fields which satisfy the discrete divergence constraint. The linear system for the pressure is derived by taking the discrete divergence of (10), 1 − ∆t
(1)
(1),∗
δUi δU − i δxi δxi
δ = δxi
δp(1) ρ(xi ),(1) δxi 1
!
.
(11)
(1)
δU The term i is formed from (4) using values which have been updated to the t(1) = t(n) +∆t δxi time location. The linear system (11) has variable coefficients due to the inverse density multiplying the pressure gradient. Equation (11) is solved for p(1) and the result is used in the projection (10). This completes the first stage of the projection scheme.
2.2
Stage 2
Desired update for the second stage of the RK2 scheme is (n+1)
Ui
1 (n) (2) Ui + Ui , 2 " #! 1 δp(2) 1 (n) (1) (1) Ui + Ui − ∆t Fi + (xi ),(2) , = 2 δxi ρ
=
where
ρ(2)
(1)
δ ρ(xj ),(1) Uj (1) = ρ − ∆t δxj
7
(12)
.
(13)
The estimated update with a guessed pressure of zero is (n+1),∗
Ui
=
h i 1 (n) (1) (1) Ui + Ui − ∆t Fi . 2
(14)
Subtracting (14) from (12) yields the projection (n+1)
Ui
(n+1),∗
= Ui
−
∆t 1 δp(2) . 2 ρ(xi ),(2) δxi
(15)
Taking the divergence of (15) gives the linear system 2 − ∆t
(n+1)
δUi δxi
(n+1),∗
δU − i δxi
δ = δxi
δp(2) ρ(xi ),(2) δxi 1
!
.
(16)
(n+1)
δUi , is obtained from (4) using updated thermal δxi properties. For example, the updated density is given by the final RK2 stage, The divergence of the updated velocity,
ρ(n+1) =
1 (n) ρ + ρ(2) , 2
(17)
where ρ(2) is given by (13). Note that the factors of 2 in (15) and (16) cancel. Hence, one may omit the 2 in the linear system if one also omits the 2 in the projection. The projection and linear system are then given by (n+1)
Ui and 1 − ∆t
(n+1)
δUi δxi
(n+1),∗
= Ui
(n+1),∗
δU − i δxi
−
δp(2) . ρ(xi ),(2) δxi ∆t
δ = δxi
δp(2) ρ(xi ),(2) δxi 1
(18) !
,
(19)
respectively. A final comment on this particular form of the two-stage projection: One may notice that the desired update given by (12) is set up so that we first average the velocity value at t(2) and then project the update to obey the divergence constraint. Alternatively, we could have first projected the Stage 2 velocity and then averaged that result with the current velocity to achieve effectively the same RK2 update. However, as pointed out by Borodai in [7], the former strategy is to be preferred because it prevents the accumulation of mass 8
divergence errors (which can be more easily assessed in the incompressible case). We refer to this approach as the “average → project” approach, as opposed to the “project → average” approach.
3
The two-stage projection scheme in FDS notation
To facilitate comparison to the notation used in the FDS Technical Reference Guide [8] (hereafter referred to as “FDS5”) we will often refer to equation numbers from that document. Because the momentum equation is presented in Stokes’ form in FDS5, the form of the semi-discrete update is written
δH dUi = − Fi + dt δxi
.
(20)
Here the form of Fi is different than that given by (6), but this is inconsequential for the current discussion. Hence, we do not spell it out explicitly. Further, the scalar H is a “pressure-like” variable (hereafter referred to simply as pressure), which, in the FDS formulation, plays the role of enforcing the divergence constraint. The only difference between (20) and (5) is that in the FDS formulation the inverse density does not multiply grad H. As we have mentioned, the types of projection schemes presented here are usually referred to as “semi-implicit” schemes, because the pressure is viewed as being treated “implicitly.” To be precise, calling a scheme “implicit” means that certain parameters on the RHS are supposed to be located at time t(n+1) . Strictly speaking, this is not the case for the pressure, which has the sole purpose of enforcing mass conservation given the form of the stated desired numerical update for the momentum. Consider, for example, that the force term Fi could be evaluated at time level t(n) , or it could be evaluated using Adams-Bashforth extrapolation (which would involve values from time t(n−1) ); or, part of Fi could be treated implicitly (here in the strict sense of the word, such that terms at time t(n+1) are present). The resulting pressure for each case would be different. Hence, the main point is that pressure is implicit only in the loose sense that solving for the pressure requires inversion of a linear system. Pretending that the pressure is evaluated at a specific time location 9
results in a first-order representation of the pressure∗ . This is acceptable, however, because in the low-Mach equations the hydrodynamic pressure does not affect other thermodynamic quantities.
3.1
FDS5 Stage 1
Now, just as the pressure is not located at t(n+1) , it is also not located at t(n) , as is suggested by the notation used in FDS5. The notation in FDS5 is designed to present the RK2 scheme as a “predictor-corrector” scheme. Hence, what we referred to as the “desired Stage 1 (1)
velocity” in the previous section, denoted Ui , is referred to as the “estimated velocity” in (n+1)e
FDS5, denoted Ui
. The desired Stage 1 update in FDS5 is written as " (n+1)e Ui
=
(n) Ui
#
δH(n) + . δxi
(n) Fi
− ∆t
(21)
This is Eq. (4.8) in FDS5 and is analogous to (8) in Section 2. What is not explicitly stated in FDS5 is what we referred to in Section 2 as the “uncorrected velocity.” We will again denote the uncorrected velocity with an asterisk. In FDS5 notation we have, (n+1)e ,∗
Ui
h
(n)
= Ui
(n)
− ∆t Fi
i
.
(22)
Subtracting (22) from (21) yields the projection (n+1)e
Ui
(n+1)e ,∗
= Ui
− ∆t
δH(n) . δxi
(23)
Taking the divergence of (23) yields the Poisson equation for the pressure, 1 − ∆t
(n+1)e
δUi δxi
(n+1)e ,∗
δU − i δxi
=
δ 2 H(n) , δxi δxi
(24)
which is identical to Eq. (4.7) in FDS5. Even though the Poisson equation and the projection are here presented in a different order than that given in FDS5, the method is identical. This completes Stage 1 of the projection scheme. ∗
The pressure may be obtained to second order accuracy if the momentum equation is advanced using
Crank-Nicolson [4].
10
3.2
FDS5 Stage 2
Again, note that FDS5 denotes the pressure for the second stage by H(n+1)e , but this pressure is completely analogous to p(2) used in Section 2. The desired update for the second stage of the time integration is "
(n+1) Ui
1 δH(n+1)e (n) (n+1)e (n+1)e = Ui + Ui − ∆t Fi + 2 δxi
#!
.
(25)
The uncorrected update, assuming zero pressure, is (n+1),∗
Ui
=
h i 1 (n) (n+1)e (n+1)e Ui + Ui − ∆t Fi . 2
(26)
Subtracting (26) from (25) yields the projection (n+1)
Ui
(n+1),∗
= Ui
−
∆t δH(n) . 2 δxi
(27)
Taking the divergence of (27) yields the Poisson equation 2 − ∆t
(n+1)
δUi δxi
(n+1),∗
δU − i δxi
=
δ 2 H(n+1)e , δxi δxi
(28)
which is identical to Eq. (4.11) in FDS5. As before, one may omit the factors of 2 in (27) and (28). We note that the FDS5 projection scheme correctly takes the “average → project” approach, as discussed at the end of Section 2.2.
3.3
Summary of the FDS5 two-stage projection algorithm
In this section, we provide a step-by-step procedure for updating the velocity. We do this because the order in which the equations are presented in the derivation above is not the same as the order in which they are typically solved computationally. The algorithm presented below is equivalent to that given in FDS5 (including cost), though the storage of data is likely somewhat different. Stage 1: 1. Solve for the uncorrected velocity via (22). 11
2. Form and solve the Poisson equation (24). 3. Project the uncorrected velocity to the corrected velocity via (23). Stage 2: 1. Solve for the uncorrected velocity via (26). 2. Form and solve the Poisson equation (28), omitting the factor of two. 3. Project the uncorrected velocity to the corrected velocity via (27), omitting the factor of two.
4
Conclusions
The main purpose of this document is to show that the FDS5 time integration scheme can be viewed as a projection method. In Section 2, we derive a classical projection scheme, which is similar to that used by the CCSE group at Berkeley Labs [1]. In Section 3, we show that the FDS5 scheme is equivalent to the CCSE scheme with the exception of the density weighting in the Poisson equation† . Thus, indeed the FDS5 scheme is a projection method. We also note that the FDS5 scheme correctly utilizes the “average → project” approach, which eliminates the accumulation of divergence errors. The FDS5 scheme is second-order in velocity and first-order in pressure. Due to the use of the Stokes form of the momentum equation, the FDS5 scheme does not discretely conserve momentum, but does so to secondorder in space. Further, as with most all variable-density projection schemes, the FDS5 scheme does not discretely conserve kinetic energy and therefore relies on numerical viscosity for stability. †
This difference is due to the use of the Stokes form of the momentum equation in FDS5, which improves
the speed of the algorithm at the cost of a small time lag in the baroclinic torque.
12
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] J. Bell. AMR for low Mach number reacting flow. Lawrence Berkeley National Laboratory Paper LBNL-54351, 2004. [2] A. J. Chorin and J. E. Marsden. A Mathematical Introduction to Fluid Mechanics. Springer, third edition, 1990. [3] S. Gottlieb, C. W. Shu, and E. Tadmor. Strong stability-preserving high-order time discretization methods. SIAM Review, 43(1):89–112, 2001. [4] F. E. Ham, F. S. Lien, and A. B. Strong. A fully conservative second-order finite difference scheme for incompressible flow on non-uniform grids. J. Comp. Phys., 177:117–133, 2002. [5] A. Kempf. Large-Eddy Simulation of Non-Premixed Turbulent Flames. PhD thesis, Technischen Universit¨at Darmstadt, 2003. [6] Harvard Lomax, Thomas H. Pulliam, and David W. Zingg. Fundamentals of Computational Fluid Dynamics. Springer, 2001. [7] R. McDermott. Toward One-Dimensional Turbulence Subgrid Closure for Large-eddy Simulation. PhD thesis, The University of Utah, 2005. [8] K. McGrattan, S. Hostikka, J. Floyd, H. Baum, and R. Rehm. ics Simulator (Version 5) Technical Reference Guide. http://fire.nist.gov/fds/, 2007.
13
Fire Dynam-
NIST Special Pub. 1018-5,
[9] Y. Morinishi, T. S. Lund, O. V. Vasilyev, and P. Moin. Fully conservative high order finite difference schemes for incompressible flow. J. Comp. Phys., 143:90–124, 1998. [10] C. D. Pierce. Progress-Variable Approach for Large-Eddy Simulation of Turbulent Combustion. PhD thesis, Stanford University, 2001.
14
