LINEBURG

(2µux )x + (µ(uy + vx ))y + (µ(uz + wx ))z

px

condition

= , (21.15)

ρ ρ

[p] ’ 2 [µ] u · N , v · N , w · N · N = σκ (21.10)

(µ(uy + vx ))x + (2µvy )y + (µ(vz + wy ))z

py

= , (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

pz

« = , (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-

u

continuous nature of several quantities across the interface, special care

+ N N v NT N

T

(21.11)

is needed in discretizing the viscous terms in equation (18.19) and in

w

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

w

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

2

φi,j,k + φi+1,j,k

Du Dv Dw

µi+ 1 ,j,k = µ , (21.18)

= = = 0. (21.12)

2

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)

2

i.e., 2

as was done in equation (18.16).

T

· „)

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

ρ

as

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

where

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.

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 |

and

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

x

jump condition as

arriving at

uR ’ uI uI ’ uM

’ µ’

µ+ = JI , (21.34)

(µux )R ’ (µux )L (1 ’ θ) x

θx

(µuxx )i+ 1 ,j,k = (21.25)

x

2

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

and

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

µ+ µ’ 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

2

= fk (21.47)

x

Consider solving

and

· (β 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

2

= 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 ’ θ)

βa

x x

2

= fk + + (21.49)

surface-tension e¬ects directly, the jump conditions for the pressure cannot ( x)2 β+ x

x

be ignored.

and

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

2

= 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

x

ˆ

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,

form

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)

and

[β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

x

pi,j+1 ’pi,j pi,j ’pi,j’1

’ βi,j’1/2

βi,j+1/2 y 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,

p

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.

by

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.

Copyright Design by: Sunlight webdesign