<< . .

( 23)

. . >>

velocity ¬eld is continuous across the interface. Using this, the jump
(2µux )x + (µ(uy + vx ))y + (µ(uz + wx ))z
= , (21.15)
ρ ρ
[p] ’ 2 [µ] u · N , v · N , w · N · N = σκ (21.10)
(µ(uy + vx ))x + (2µvy )y + (µ(vz + wy ))z
= , (21.16)
can be written for the pressure. Notice that this reduces to [p] = σκ when ρ ρ
the viscosity is continuous across the interface. Further derivations lead to
(µ(uz + wx ))x + (µ(vz + wy ))y + (2µwz )z
«  = , (21.17)
[µux ] [µuy ] [µuz ] ρ ρ
 [µvx ] [µvy ] [µvz ] 
in expanded form.
[µwx ] [µwy ] [µwz ]
The two-phase incompressible ¬‚ow equations are discretized in the
«« « T « 
0 0
u same manner as the equations for one-phase incompressible ¬‚ow. First,
= [µ]  v   T1   T1  an intermediate velocity ¬eld V is computed using equation (18.19).
w Then equation (18.21) is solved to ¬nd the pressure, which is used in
T2 T2
«  equation (18.20) to make the velocity ¬eld divergence free. Due to the dis-
continuous nature of several quantities across the interface, special care
+ N N  v  NT N
is needed in discretizing the viscous terms in equation (18.19) and in
discretizing the pressure and density in equations (18.21) and (18.20).

« T « « T
0 0 u
 T1   T1   v  N T N · ,
’ 
21.3 Viscous Terms
T2 T2

which is useful for discretizing the viscous terms, especially since the right- In the δ-function approach the density and the viscosity are numerically
hand side of this equation involves only derivatives that are continuous smeared out, so that they are continuous across the interface. The contin-
across the interface. Notice that all the quantities on the left-hand side of uous viscosity simpli¬es the jump conditions, allowing the viscous terms
this equation become continuous across the interface when [µ] = 0 forces to be discretized just as they were for one-phase incompressible ¬‚ow. The
the right-hand side to be identically zero. See [106] for details. only di¬erence is that ρ and µ are de¬ned by equations (21.1) and (21.2)
Since the velocity is continuous across the interface, the material respectively. Averaging of φ is preferred to averaging of other quantities.
derivative, or Lagrangian, acceleration is continuous as well, For example, the viscosity at xi+ 1 ,j,k is de¬ned as

φi,j,k + φi+1,j,k
Du Dv Dw
µi+ 1 ,j,k = µ , (21.18)
= = = 0. (21.12)
Dt Dt Dt 2

as opposed to
Since the Navier-Stokes equations (18.2), (18.3), and (18.4) are valid on
µ (φi,j,k ) + µ (φi+1,j,k )
both sides of the interface, these equations do not jump across the interface; µi+ 1 ,j,k = , (21.19)
i.e., 2

as was done in equation (18.16).
· „)
p( If one prefers to keep the density and viscosity sharp across the interface,
Vt + V · ’ ’ g = 0.
V+ (21.13)
the sign of the level set function can be used to determine µ as µ’ or µ+
ρ ρ
21.3. Viscous Terms 233 234 21. Two-Phase Incompressible Flow

and to determine ρ as ρ’ or ρ+ . Consider the case where µ and ρ are both Suppose that φL ¤ 0 and φM > 0, so that the interface lies between the
independently spatially constant on either side of the interface, allowing associated grid points. Then
the simpli¬cation of the viscous terms to |φL |
θ= (21.26)
|φL | + |φM |
µ (uxx + uyy + uzz )
, (21.20)
ρ can be used to estimate the interface location by splitting this cell into
two pieces of size θ x on the left and (1 ’ θ) x on the right. At the
interface, we denote the continuous velocity by uI and calculate the jump
µ (vxx + vyy + vzz )
, (21.21) as JI = θJM + (1 ’ θ)JL . Then we discretize the jump condition [µux ] = JI
uM ’ uI uI ’ uL
µ (wxx + wyy + wzz ) ’ µ’
µ+ = JI , (21.27)
, (21.22) (1 ’ θ) x θx
solving for uI to obtain
· V = 0. Since the velocities are continuous, their ¬rst
with the aid of
µ+ uM θ + µ’ uL (1 ’ θ) ’ JI θ(1 ’ θ) x
derivatives can be computed directly. However, the jump conditions in uI = , (21.28)
µ+ θ + µ’ (1 ’ θ)
equation (21.11) are needed to compute the second derivatives.
The right-hand side of equation (21.11) needs to be computed in order to so that we can write
evaluate the jumps across the interface. First, the continuous velocity ¬eld is uM ’ uI uM ’ u L µ JI θ
(µux )L = µ+ =µ
ˆ + , (21.29)
averaged from the MAC grid to the cell centers. Then central di¬erencing
(1 ’ θ) x x
is used to compute the ¬rst derivatives at each cell center. These ¬rst
derivatives are multiplied by the appropriate components of the normal
µ+ µ’
and tangent vectors to obtain a numerical estimate for the right-hand side
ˆ (21.30)
of equation (21.11), which we denote by the matrix J. Since J is a spatially µ+ θ + µ’ (1 ’ θ)
continuous function, spatial averages can be used to de¬ne J elsewhere. For
de¬nes an e¬ective µ. Similarly, if φL > 0 and φM ¤ 0, then
example, Ji+ 1 ,j,k = (Ji,j,k + Ji+1,j,k ) /2.
uM ’ uI uM ’ uL µ JI θ
Once J has been computed, the second derivatives are computed using (µux )L = µ’ ’

ˆ , (21.31)
(1 ’ θ) x x
techniques similar to those developed in [106]. For example, consider the
discretization of µuxx at xi+ 1 ,j,k using uM = ui+ 1 ,j,k and its neighbors where
2 2
uL = ui’ 1 ,j,k and uR = ui+ 3 ,j,k . We need averaged values of φ and J 11 µ’ µ +
2 2
µ= ’
ˆ (21.32)
(the appropriate scalar entry of J) at the same three spatial locations as µ θ + µ+ (1 ’ θ)
the u terms. If φL , φM , and φR are all greater than zero, we de¬ne
de¬nes an e¬ective µ.
uM ’ uL In similar fashion, if φR > 0 and φM ¤ 0, then
(µux )L = µ+ (21.23)
x |φR |
θ= (21.33)
|φR | + |φM |
is used to estimate the interface location with (1 ’ θ) x on the left and
uR ’ uM
(µux )R = µ+ , (21.24) θ x on the right. Then JI = θJM + (1 ’ θ)JR is used to discretize the
jump condition as
arriving at
uR ’ uI uI ’ uM
’ µ’
µ+ = JI , (21.34)
(µux )R ’ (µux )L (1 ’ θ) x
(µuxx )i+ 1 ,j,k = (21.25)
resulting in
µ’ uM θ + µ+ uR (1 ’ θ) ’ JI θ(1 ’ θ) x
in the standard fashion. A similar discretization holds when all three φ
uI = (21.35)
values are less than or equal to zero. µ’ θ + µ+ (1 ’ θ)
21.4. Poisson Equation 235 236 21. Two-Phase Incompressible Flow

and and solve for pI as
uI ’ uM uR ’ u M µ JI θ
ˆ β + pk+1 θ + β ’ pk (1 ’ θ) ’ bθ(1 ’ θ) x
(µux )R = µ’ ’

ˆ , (21.36) pI = , (21.43)
(1 ’ θ) x x β + θ + β ’ (1 ’ θ)
where so that approximations to the derivatives on the left and right sides of the
’+ interface can be written as
ˆ (21.37)
µ’ θ + µ+ (1 ’ θ) pI ’ p k ˆ pk+1 ’ pk ’ βb(1 ’ θ)
β’ =β (21.44)
θx x
de¬nes an e¬ective µ. If φR ¤ 0 and φM > 0, then
uI ’ uM uR ’ u M µ JI θ
(µux )R = µ+ =µ
ˆ + , (21.38) ˆ
(1 ’ θ) x x pk+1 ’ pI ˆ pk+1 ’ pk βbθ
β+ =β + , (21.45)
(1 ’ θ) x x
µ+ µ’ where
µ= +
ˆ (21.39)
µ θ + µ’ (1 ’ θ) ˆ
β= (21.46)
β + θ + β ’ (1 ’ θ)
de¬nes an e¬ective µ.
de¬nes an e¬ective β.
The new equations for the unknowns pk and pk+1 are then
21.4 Poisson Equation (pk+1 ’a)’pk pk ’pk’1
ˆ b(1’θ)
’ ’ βk’ 1
β β+
x x
= fk (21.47)
Consider solving
· (β p) = f (21.40)
pk+2 ’pk+1 uk+1 ’(pk +a)
ˆ bθ
βk+ 3 +
for the pressure p with speci¬ed jump conditions of [p] = a and [βpn ] = b β’
x x
= fk+1 , (21.48)
across the interface. In the context of equation (18.21), β = 1/ρ and x
f = ( · V )/ t. When the δ-function method is used, this equation is where we add the a term to correct for the fact the pressure is discontinuous
straightforward to solve, since a and b are set to 0, and equation (21.1) across the interface as well. Note that this correction was not necessary
is used to de¬ne a smeared-out continuous value for ρ. In this case, one in treating the viscosity, since the velocity ¬eld is continuous across the
obtains a numerically smeared-out pressure pro¬le that does not include interface. These new equations for pk and pk+1 can be rewritten in standard
surface-tension forces. As discussed above, the forcing function de¬ned form as
in equation (21.3) can be added to the momentum equations to recover
pk+1 ’pk pk ’pk’1
ˆ ’ βk’ 1
β ˆ ˆ
surface-tension e¬ects. On the other hand, if one wants to model the βb(1 ’ θ)
x x
= fk + + (21.49)
surface-tension e¬ects directly, the jump conditions for the pressure cannot ( x)2 β+ x
be ignored.
First, consider the one-dimensional problem where a standard second-
pk+2 ’pk+1 pk+1 ’pk
order accurate discretization of ’β
βk+ 3 ˆ ˆ
βa βbθ
x x
= fk+1 ’ +’ (21.50)
pi+1 ’pi pi ’pi’1
’ βi’ 1
βi+ 1 ( x)2
x β x
x x
2 2
= fi (21.41)
emphasizing that this discretization yields the standard symmetric linear
system with βk+1/2 = β.
can be written for each unknown pi . Suppose that the interface is lo-
More generally, at each grid point i, we write a linear equation of the
cated between xk and xk+1 . As in the treatment of the viscosity term,
we discretize the jump condition [βpx ] = b, obtaining, for example,
pi+1 ’pi pi ’pi’1
’ βi’ 1
βi+ 1
pk+1 ’ pI pI ’ p k x x
’ β’
β+ = fi + F L + F R
2 2
= b, (21.42) (21.51)
(1 ’ θ) x θx x
21.4. Poisson Equation 237 238 21. Two-Phase Incompressible Flow

and assemble the system of linear equations into matrix form. Each βk+1/2 are already balanced out on the right-hand side by the appropriate jumps
is evaluated based on the side of the interface that xk and xk+1 lie on, and included in the V term.
a special β is used when xk and xk+1 lie on opposite sides of the interface. The resulting linear system of equations can still be solved using a PCG
Then if the left arm of the stencil (the line segment connecting xi and gradient with an incomplete Choleski preconditioner, just as in the case
xi’1 ) crosses the interface, a nonzero F L is de¬ned with correction terms of one-phase incompressible ¬‚ow. However, one needs to use caution when
for [p] = a and [βpn ] = b. Likewise, if the right arm of the stencil (the plugging the resulting pressure into equation (18.20), since the pressure
line segment connecting xi and xi+1 ) crosses the interface, a nonzero F R is discontinuous across the interface. The pressure derivatives in equa-
is de¬ned with correction terms for [p] = a and [βpn ] = b. tion (18.20) should be computed in exactly the same fashion as they were
The multidimensional approach is treated in a dimension-by-dimension computed in solving equation (18.21); i.e., the correction terms are still
fashion. While the [p] = a jump condition is trivial to apply, some assump- needed.
tions are made in order to obtain a dimension-by-dimension approach for
[βpn ] = b. For example, in two spatial dimensions we assume that
[βpx ] = [βpn ]n1 (21.52)
[βpy ] = [βpn ]n2 , (21.53)
where n1 and n2 are the components of the local unit normal. Although
these equations are not generally true, adding n1 times the ¬rst equation to
n2 times the second equation leads to the correct [βpn ] = b jump condition
and numerically demonstrated convergence to the correct solution. The
errors in this approach can be characterized by adding t1 times the ¬rst
equation to t2 times the second equation to obtain [βpt ] = 0, implying that
the tangential derivative is incorrectly smeared out. The two-dimensional
application of the method consists in writing a linear equation of the form
pi+1,j ’pi,j pi,j ’pi’1,j
’ βi’1/2,j
βi+1/2,j x x
pi,j+1 ’pi,j pi,j ’pi,j’1
’ βi,j’1/2
βi,j+1/2 y y
x y
= fi,j + F + F (21.54)
at each grid point, where F x = F L + F R and F y = F B + F T are obtained
by considering each spatial dimension independently using either [βpx ] =
[βpn ]n1 = bn1 or [βpy ] = [βpn ]n2 = bn2 , respectively.
Before using the above-described numerical method to solve equa-
tion (18.21), the jump condition given in equation (21.10) needs to be
computed. This can be done with standard central di¬erencing of the av-
eraged cell-centered velocities, analogous to the way that J was computed
in discretizing the viscous terms. Note that we can set [px /ρ] = [py /ρ] =
[pz /ρ] = 0 in spite of the nonzero jumps in these quantities. Since the full
equations are continuous across the interface, one can take the divergence
of the full equations to derive equation (18.21) without the need for correc-
tion terms. The jumps in the derivatives of the pressure in equation (18.21)
240 22. Low-Speed Flames

to the discontinuity of the density, viscosity, and pressure). Delta-function
formulations smear out this velocity jump, forcing a continuous velocity
¬eld across the interface. This can be problematic, since this numerical
smearing adds a compressible character to the ¬‚ow ¬eld near the inter-
22 face. The divergence-free condition is not exactly satis¬ed in the separate
subdomains. In addition, di¬culties arise in computing the interface veloc-
ity, which depends on the local velocity of the unreacted material. Near the
Low-Speed Flames interface, the velocity of the unreacted material contains large O(1) numer-
ical errors where it has been nonphysically forced to be continuous with the
velocity of the reacted material. Partial solutions to these problems were
proposed by Helenbrook et al. [84], where the authors were able to remove
the numerical smearing of the normal velocity, obtaining a sharp interface
pro¬le. This method works well as long as ¬‚ame fronts remain well sepa-
rated with moderate curvature; see Helenbrook and Law [83]. This method
cannot treat merging ¬‚ame fronts or individual fronts with relatively high
curvature. These drawbacks were recently overcome by Nguyen et al. [119]
who extended the work of Kang et al. [91] to treat this problem.

22.2 Governing Equations
22.1 Reacting Interfaces We ignore viscous e¬ects and consider the equations for inviscid incom-
pressible ¬‚ow
In Chapter 21 the interface moved with the local ¬‚uid velocity only,
and individual ¬‚uid particles did not cross the interface. In this chap- Vt + V · V+ =0 (22.1)
ter we consider interfaces across which a chemical reaction is converting
one incompressible ¬‚uid into another. The interface moves with the lo- independently for each ¬‚uid. The interface velocity is W = DN , where D is
cal velocity of the unreacted ¬‚uid plus a reaction term that accounts for the normal component of the interface velocity de¬ned by D = (VN )u + S.
the conversion of one ¬‚uid into the other as material moves across the The “u” subscript indicates that the normal velocity is calculated using
interface. Consider an interface separating liquid and gas regions where
the velocity of the unreacted material only. This is important to note,
the liquid is actively vaporizing into the gaseous state. Juric and Tryg- since VN is discontinuous across the interface. The ¬‚ame speed is de¬ned
gvason [90] developed a front-tracking approach to this problem using a as S = So + σκ where, So and σ are constants and κ is the local curvature
δ-function formulation to treat the interface boundary conditions. Son and of the interface.
Dir [153] and Welch and Wilson [171] developed level-set-based and volume- Conservation of mass and momentum imply the standard Rankine-
of-¬‚uid-based (respectively) approaches to this same problem also using a Hugoniot jump conditions across the interface
δ-function formulation.
Another example of reacting interfaces occurs in premixed ¬‚ames. As- [ρ(VN ’ D)] = 0 (22.2)
suming an in¬nitely thin ¬‚ame front allows us to treat the ¬‚ame front as a ρ(VN ’ D)2 + p = 0 (22.3)
discontinuity separating two incompressible ¬‚ows. The unreacted material
as well as continuity of the tangential velocities, [VT1 ] = [VT2 ] = 0, as long
undergoes reaction as it crosses the interface, producing a lower-density
as S = 0. Note that S = 0 only in the case of a contact discontinuity (not
(higher-volume) reacted material. Qian et al. [135] devised a front-tracking
a ¬‚ame). Denoting the mass ¬‚ux in the moving reference frame (speed D)
approach to this problem using a δ-function formulation.
Typically, the density is discontinuous across the interface. Thus, mate-
rial must instantaneously expand as it crosses the interface, implying that
M = ρr ((VN )r ’ D) = ρu ((VN )u ’ D) (22.4)
the normal velocity is discontinuous across the interface as well (in addition
22.3. Treating the Jump Conditions 241 242 22. Low-Speed Flames

allows us to rewrite equation (22.2) as [M ] = 0. Here the “r” subscript to obtain
denotes a reacted material quantity. Substitution of D = (VN )u + S into 1 1
uG ur ’ M ’
= n1 , (22.12)
equation (22.4) yields u
ρr ρu
1 1
M = ’ρu S, (22.5) G
vr ’ M ’
vu = n2 , (22.13)
ρr ρu
which is a rather simple quantity for computations.

<< . .

( 23)

. . >>

Copyright Design by: Sunlight webdesign