LINEBURG


<< . .

 21
( 23)



. . >>

1 1
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

<< . .

 21
( 23)



. . >>

Copyright Design by: Sunlight webdesign