LINEBURG

Starting with [D] = 0, we derive G

wr ’ M ’

wu = n3 , (22.14)

ρr ρu

ρVN ’ ρ(VN ’ D)

= 0, (22.6) where N = (n1 , n2 , n3 ) is the local unit normal. Similarly, reacted ghost

ρ

velocities are de¬ned on a band of ghost cells on the unreacted side of the

ρVN ’ M

interface and used in the discretization of the reacted ¬‚uid velocity.

= 0, (22.7)

ρ When solving equation (18.23),

and p

· ·V ,

= (22.15)

1 ρ

[VN ] = M , (22.8)

ρ for the scaled pressure p , the jump in pressure given by equation 22.11 as

where the last equation follows from [M ] = 0. It is more convenient to write 1 1

[p ] = ’ tM 2 ’ (22.16)

ρr ρu

1

V =M N (22.9)

ρ is accounted for using the techniques developed in Liu et al. [106] and Kang

et al. [91].

as a summary of equation (22.8) and [VT1 ] = [VT2 ] = 0. The dot product of

Figure 22.1 shows the time evolution of two initially circular ¬‚ame fronts

equation (22.9) and N results in equation (22.8), while the dot product of as they grow to merge together. Figure 22.2 shows a snapshot of the velocity

equation (22.9) and T1 or T2 results in [VT1 ] = 0 or [VT2 ] = 0, respectively. ¬eld, illustrating its discontinuous nature across the interface. Figure 22.3

Equation (22.3) can be rewritten as shows the time evolution of two initially spherical ¬‚ame fronts as they grow

to merge together.

M2

+p =0 (22.10) Recently, Nguyen et al. [118] extended this approach to model ¬re for

ρ

computer graphics. In Figure 22.4, the φ = 0 isocontour is used to render

or a typical blue ¬‚ame core. S is smaller for the larger blue core on the right.

Figure 22.5 illustrates the e¬ect of increased expansion as the density jump

1

[p] = ’M 2 , (22.11) is increased from left to right. The yellow ¬‚ame color is calculated using a

ρ

blackbody radiation model based on the temperature pro¬le of the hot gas

again using [M ] = 0. emitted at the ¬‚ame front. Figures 22.6 shows a ball catching on ¬re, and

Figure 22.7 shows a camp¬re.

22.3 Treating the Jump Conditions

Since the normal velocity is discontinuous across the interface, caution is

needed in applying numerical discretizations near the interface. For exam-

ple, when discretizing the unreacted ¬‚uid velocity near the interface, one

should avoid using values of the reacted ¬‚uid velocity. Following the ghost

¬‚uid methodology, a band of ghost cells on the reacted side of the interface

is populated with unreacted ghost velocities that can be used in the dis-

cretization of the unreacted ¬‚uid velocity. This is done using equation (22.9)

22.3. Treating the Jump Conditions 243 244 22. Low-Speed Flames

level set contour

1

0.9

0.8

0.7

0.6

0.5

t=0

0.4

t=.01

t=.02

0.3

t=.03

t=.04

0.2

t=.05

0.1

0

0 0.5 1 1.5

Figure 22.1. Time evolution of two initially circular ¬‚ame fronts as they grow to

merge together.

velocity field (t=.035)

1

0.9

0.8

0.7

0.6

0.5

0.4

0.3

0.2

Figure 22.3. Time evolution of two initially spherical ¬‚ame fronts as they grow

to merge together.

0.1

0

0 0.5 1 1.5

Figure 22.2. Discontinuous velocity ¬eld depicted shortly after the two ¬‚ame

fronts merge.

22.3. Treating the Jump Conditions 245 246 22. Low-Speed Flames

Figure 22.4. Typical blue cores rendered using the zero isocontour of the level set

function. (See also color ¬gure, Plate 21.)

Figure 22.6. A ¬‚ammable ball catches on ¬re as it passes through a ¬‚ame. (See

also color ¬gure, Plate 23.)

Figure 22.5. The density ratio of the unburnt to burnt gas is increased from

left to right, illustrating the e¬ect of increased expansion. (See also color ¬gure,

Plate 22.)

22.3. Treating the Jump Conditions 247

23

Heat Flow

23.1 Heat Equation

Starting from conservation of mass, momentum, and energy one can derive

ρet + ρV · ·V = · (k T ) ,

e+p (23.1)

where k is the thermal conductivity and T is the temperature. Assuming

that e depends on at most temperature, and that the speci¬c heat at con-

stant volume,cv is constant leads to e = eo + cv (T ’ To ), where eo is the

internal energy per unit mass at some reference temperature To (see, for

example, Atkins [10]). This and the incompressibility assumption · V = 0

simplify equation (23.1) to

ρcv Tt + ρcv V · · (k T ) ,

T= (23.2)

which can be further simpli¬ed to the standard heat equation

· (k T )

ρcv Tt = (23.3)

Figure 22.7. Camp¬re with realistic lighting of the surrounding rocks. (See also

color ¬gure, Plate 24.)

by ignoring the e¬ects of convection, i.e., setting V = 0.

Applying explicit Euler time discretization to equation (23.3) results in

T n+1 ’ T n 1

· (k T n ) ,

= (23.4)

t ρcv

where either Dirichlet or Neumann boundary conditions can be applied on

the boundaries of the computational domain. Assuming that ρ and cv are

250 23. Heat Flow 23.3. Poisson Equation 251

constants allows us to rewrite this equation as of the two-dimensional outline depicted in Figure (23.1). If one takes the

rather simple approach of embedding this complicated domain in a uniform

T n+1 ’ T n ˆ

· k Tn

= (23.5) Cartesian grid, then a level set function can be used to de¬ne the boundary

t

of the irregular region. The heat equation 23.3 can then be solved with, for

ˆ example, Dirichlet T = g(x, t) boundary conditions applied to the bound-

with k = k/(ρcv ). Standard central di¬erencing can be used for the spatial

ary where φ = 0. More complicated boundary conditions can be used as

derivatives and a time step restriction of

well.

2 2 2

ˆ ¤1

tk + + (23.6) The spatial derivatives are computed with the aid of the given values of

( x)2 ( y)2 ( z)2

T = g(x, t) on the interface. When using explicit Euler time discretization,

is needed for stability. the time-step restriction needed for stability becomes

Implicit Euler time discretization

2 2 2

ˆ ¤ 1,

tk + + (23.11)

T n+1 ’ T n x)2 y)2 (θ3 z)2

(θ1 (θ2

ˆ

· k T n+1

= (23.7)

t

where θ1 , θ2 , and θ3 are the cell fractions in each spatial dimension for cells

avoids this time step stability restriction. This equation can be rewritten cut by the interface with 0 < θi ¤ 1. Since the θi ™s can be arbitrarily small,

as leading to arbitrarily small time steps, implicit methods need to be used,

ˆ e.g., backward Euler or Crank-Nicolson. Then a linear system of equations

T n+1 ’ t · k T n+1 = T n (23.8)

can be solved for the unknowns Tin+1 . Since the coe¬cient matrix depends

ˆ on the details of the spatial discretization, a robust method for treating

· k T n+1

where the term is discretized using central di¬erencing.

the cut cells is crucial. Below, we outline how to do this for the simpler

For each unknown Tin+1 , equation (23.8) is used to ¬ll in one row of a

variable-coe¬cient Poisson equation, which has spatial derivatives identical

matrix, creating a linear system of equations. Since the resulting matrix is

to those of the heat equation.

symmetric, a number of fast linear solvers can be used (e.g., a PCG method

with an incomplete Choleski preconditioner; see Golub and Van Loan [75]).

Equation (23.7) is ¬rst-order accurate in time and second-order accurate in

23.3 Poisson Equation

space, and t needs to be chosen proportional to x2 , in order to obtain

an overall asymptotic accuracy of O( x2 ). However, the stability of the

Consider the variable-coe¬cient Poisson equation

implicit Euler method allows one to chose t proportional to x saving

dramatically on CPU time. The Crank-Nicolson scheme · (β(x) u(x)) = f (x), (23.12)

n+1 n

’T

T 1 1

ˆ ˆ

· k T n+1 + · k Tn where β(x) is positive and bounded below by some > 0. As above, consider

= (23.9)

t 2 2 an irregularly shaped domain (as in Figure 23.1) de¬ned by a level set

can be used to achieve second-order accuracy in both space and time with function on a Cartesian grid with Dirichlet u = g(x, t) boundary conditions

t proportional to x. For the Crank-Nicolson scheme, on the φ = 0 isocontour.

For simplicity consider the one-dimensional case (βux )x = f . Since

t t

ˆ ˆ

T n+1 ’ · k T n+1 = T n + · k Tn

(23.10) β and φ are known only at the grid nodes, their values between grid

2 2

nodes are de¬ned by the linear average of the nodal values, e.g., βi+ 1 =

is used to create a symmetric linear system of equations for the unknowns 2

(βi + βi+1 ) /2. In the absence of cut cells, the standard discretization

Tin+1 . Again, all spatial derivatives are computed using standard central

di¬erencing. ui+1 ’ui ui ’ui’1

’ βi’ 1

βi+ 1 x x

2 2

= fi (23.13)

x

23.2 Irregular Domains can be used to solve this problem. For each unknown ui , equation (23.13)

is used to ¬ll in one row of a matrix, creating a linear system of equations.

Instead of a uniform Cartesian domain, suppose we wish to solve equa- The resulting matrix is symmetric and can be solved with a number of fast

tion (23.3) on an irregularly shaped domain, for example in the interior linear solvers.

252 23. Heat Flow 23.3. Poisson Equation 253

and

2uI + 2θ2 ’ 2 ui + ’θ2 + 1 ui’1

uG = (23.17)

i+1

θ2 + θ

with constant, linear, and quadratic extrapolation respectively. Here θ ∈

12

[0, 1] is de¬ned by θ = (xI ’ xi ) / x, and it can be calculated as θ =

10

|φ|/ x, since φ is a signed distance function vanishing at xI . Since equa-

tions (23.16) and (23.17) are poorly behaved for small θ, they are not used

8

when θ ¤ x. Instead, ui is set equal to uI , which e¬ectively moves the

6

interface location from xI to xi . This second-order accurate perturbation

of the interface location does not degrade the overall second-order accuracy

4

of the solution obtained using equation (23.13) to solve for the remaining

2

unknowns. Furthermore, ui = uI is second-order accurate as long as the

solution has bounded ¬rst derivatives.

0

3 Plugging equation (23.17) into equation (23.13) gives an asymmetric

1 discretization of

2

0.5

ui ’ui’1

uI ’ui

0 ’

1

θx x

-0.5

= fi (23.18)

0 .5 (θ x + x)

-1

(when β = 1). Equation (23.18) is the asymmetric discretization used by

Chen et al. [43] to obtain second-order accurate numerical methods in the

· (β u) = f

Figure 23.1. Solution of the two-dimensional Poisson equation

context of solving Stefan problems. Alternatively, Gibou et al. [44] pointed

with Dirichlet boundary conditions. The circles are the computed solution, and

out that plugging equation (23.16) into equation (23.13) gives a symmetric

the solid line contour outlines the irregularly shaped computational domain.

discretization of

ui ’ui’1

uI ’ui

’ βi’ 1

βi+ 1 θx x

2 2

= fi (23.19)

x

Suppose that an interface point xI is located between two grid points xi

and xi+1 with a Dirichlet u = uI boundary condition applied at xI . Con- based on linear extrapolation in the cut cell. It turns out that this sym-

sider computing the numerical solution in the domain to the left of xI . metric discretization is second-order accurate as well. Moreover, since the

Equation (23.13) is valid for all the unknowns to the left and including discretization is symmetric, the linear system of equations can be solved

ui’1 , but can no longer be applied at xi to solve for ui , since the subdo- with a number of fast methods such as PCG with an incomplete Choleski

main to the left of xI does not contain a valid value of ui+1 . This can be preconditioner.

remedied by de¬ning a ghost value of uG at xi+1 and rewriting equation To see this, assume that the standard second-order accurate discretiza-

i+1

(23.13) as tion in equation (23.13) is used to obtain the standard linear system of

equations for u at every grid node except xi , and equation (23.14) is used

uG ’ui ui ’ui’1

to write a linear equation for ui , introducing a new unknown uG . The sys-

’ βi’ 1

i+1

βi+ 1 i+1

x x

2 2

tem is closed with equation (23.16) for uG . In practice, equations (23.16)

= fi (23.14) i+1

x

and (23.14) are combined to obtain equation (23.19) and a symmetric lin-

ear system. Solving this linear system of equations leads to well-determined

in order to solve for ui . Possible candidates for uG include

i+1

values of u at each grid node in the subdomain as well as a well-determined

value of uG (from equation (23.16)). Designate u as the solution vector

i+1

uG = u I (23.15)

i+1 containing all these values of u.

uI + (θ ’ 1) ui Next, consider a modi¬ed problem where a Dirichlet boundary condi-

uG = (23.16)

i+1

tion of ui+1 = uG is speci¬ed at xi+1 with uG chosen to be the value

θ i+1 i+1

254 23. Heat Flow 23.4. Stefan Problems 255

of uG from u (de¬ned above). This modi¬ed problem can be discretized velocity is W = DN , where D = (VN )u + S for some reaction speed S.

i+1

to second-order accuracy everywhere using the standard discretization in Here the “u” subscript denotes an unreacted material quantity. Including

equation (23.13) at every node except at xi , where equation (23.14) is used. the e¬ects of thermal conductivity, the Rankine-Hugoniot jump condition

Note that equation (23.14) is the standard second-order accurate discretiza- for conservation of energy is

tion when a Dirichlet boundary condition of ui+1 = uG is applied at xi+1 .

ρ(VN ’ D)2

i+1

+ p (VN ’ D) = k T · N ,

Thus, this new linear system of equations can be solved in standard fashion ρe + (23.21)

2

to obtain a second-order accurate solution at every grid node. The realiza-

tion that u (de¬ned above) is an exact solution to this new linear system where we have assumed that D = VN (i.e., S = 0), so that the tangen-

implies that u is a valid second-order accurate solution to this modi¬ed tial velocities are continuous across the interface. This equation can be

problem. rewritten as

Since u is a second-order accurate solution to the modi¬ed problem, u ρ2 S 2 1

p

+u

’ρu S = k T ·N

e+ (23.22)

can be used to obtain the interface location for the modi¬ed problem to

ρ2

ρ 2

second-order accuracy. The linear interpolant that uses ui at xi and uG ati+1

using the Rankine-Hugoniot jump condition for conservation of mass,

xi+1 predicts an interface location of exactly xI . Since higher-order accurate

[ρ(VN ’ D)] = 0. Assuming that the enthalpy per unit mass h = e + (p/ρ)

interpolants (higher than linear) can contribute at most an O(∆x2 ) pertur-

depends on at most temperature, and that the speci¬c heat at constant

bation of the predicted interface location the interface location dictated by

pressure cp is constant leads to h = ho + cp (T ’ To ), where ho is the en-

the modi¬ed problem is at most an O(∆x2 ) perturbation of the true inter-

thalpy per unit mass at some reference temperature To ; see [10]. This allows

face location, xI . Thus, u is a second-order accurate solution to a modi¬ed

us to rewrite equation (23.22) as

problem where the interface location has been perturbed by O(∆x2 ). This

makes u a second-order accurate solution to the original problem as well. ρ2 S 2 1

u

’ρu S [ho ] + [cp ] (TI ’ To ) + = k T ·N , (23.23)

Note that plugging equation (23.15) into equation (23.13) e¬ectively per-

ρ2

2

turbs the interface location by an O( x) amount, resulting in a ¬rst-order

where we have used the fact that the temperature is continuous across the

accurate algorithm.

interface, [T ] = 0, and labeled the interface temperature TI . It is convenient

When β is spatially varying, βi+1/2 in equation (23.19) can be determined

to choose the reference temperature To equal to the standard temperature

G G

from a ghost value βi+1 and the usual averaging βi+1/2 = βi + βi+1 /2,

at which the reaction takes place; e.g., in the case of freezing water To =

where the ghost value is de¬ned using linear extrapolation

Copyright Design by: Sunlight webdesign