ATOMIZAÇÃO_phd3

7/27/2019 ATOMIZAO_phd3

1/23

Chapter 3

NUMERICAL IMPLEMENTATION OF THE

VOLUME OF FLUIDCONTINUOUS

SURFACE FORCE METHOD

The purpose of computing is insight, not numbers. This motto isoften thought to mean that the numbers from a computing machine

should be read and used, but there is much more to the motto. Thechoice of the particular formula, or algorithm, influences not only thecomputing but also how we are to understand the results when theyare obtained . . . Thus computing is, or at least should be, intimatelybound up with both the source of the problem and the use that is going

to be made of the answers it is not a step to be taken in isolation

from reality.

R. W. Hamming (1986)

The solutions to complex free surface flow problems involve the solution

of the equations of motion and continuity for the two fluids subject to specified

boundary conditions and their numerical solution is briefly reviewed in section 1.3.

In this dissertation we use the Eulerian Volume of Fluid (VOF) method (Hirt

and Nichols, 1981), in which a marker function convected by the flow is used to

track the interface. The major incentive for using the VOF method is that it

allows for the description of highly complex interfaces such as those encountered

in multiphase flows. Reasonable accuracy is attainable with elemental control

volume balances, yet the method is relatively simple to implement. This algorithm

is publicly available in a two-fluid, 2-D program called SOLA-VOF (Nichols et al.,

1980). A 3-D version called FLOW-3D (Hirt, 1988) is also available commercially.

More recently a one-fluid program, RIPPLE, which implements the Continuum

33


2/23

34

Surface Force (CSF) algorithm and incorporates various improvements in the one-

fluid VOF algorithm, has been introduced by Kothe et al. (Kothe et al., 1991;

Brackbill et al., 1992). In this chapter we describe how we have combined the two-

fluid capability of the SOLA-VOF algorithm with the free surface implementation

of the CSF algorithm to investigate problems that involve transient 2-D free surface

flows with two immiscible fluids.

The source code for SOLA-VOF has been extensively modified and the

resulting code tested for consistency and accuracy with numerous test problems. It

has been rewritten to accommodate the addition of the second fluid with arbitrary

viscosity, to correct the incorporation of the viscous stress terms, and to add a

superior surface tension algorithm. The original SOLA-VOF 2-D program (Nichols

et al., 1980) is well suited for high Reynolds number flows, including those involving

free surfaces. Among the latter, however, it is better suited for gas-liquid than

for liquid-liquid systems, and relaxing this limitation has been an important part

of our efforts. Thus, we have implemented extensive modifications in the original

program. Both planar and axisymmetric 2-D flows can be simulated.

This chapter outlines the numerical methods used in combining the VOF(Hirt and Nichols, 1981) and CSF (Brackbill et al., 1992) algorithms. Although

the VOF references describe the algorithm for the iteration procedure, very little

justification is provided there. This chapter addresses this issue.

3.1 Governing Equations in a Cylindrical Coordinate System

It is assumed that the flow in each phase is axisymmetric, viscous, and

incompressible. The continuity equation (2.30) given in cylindrical, axisymmetric

coordinates (r, z) as:

u

r+

v

z+

u

r= 0 (3.1)


3/23

35

where (u, v) are the radial, axial components of the velocity field respectively. The

dynamic momentum equation (2.63) becomes:

u

t + u

u

r + v

u

z

=

p

r +

rr

r +

zr

z +

rr

r + gr +

F

r (3.2)

v

t+ u

v

r+ v

v

z

=

p

z+

rzr

+zzz

+rzr

+ gz + F

z(3.3)

where p is the pressure, gr, gz are the radial and axial components of the gravita-

tional acceleration, rr, zr , rz, zz are the components of the Newtonian stress

tensor from equations (2.32) and (2.35):

rr

= 2u

r,

zr=

rz= v

r+

u

z ,

zz= 2

v

z(3.4)

and is the curvature of the liquid-liquid interface (Kothe et al., 1991; Brackbill

et al., 1992) from equation (2.56),

= (n) (3.5)

where the unit normal (directed into fluid 2),

n = n|n|

(3.6)

is derived from a normal vector,

n = F (3.7)

The evolution equation (2.73) for the fluid function marker field, F, is:

F

t + u

F

r + v

F

z = 0 (3.8)

The density and viscosity fields are obtained from:

= 1 (1 F) + 2F (3.9)


4/23

36

= 1 (1 F) + 2F (3.10)

Note that the interfacial surface forces are incorporated as body forces per

unit volume in the momentum equations (3.2) and (3.3) rather than as boundary

conditions. Instead of a surface tensile force boundary condition applied at a

discontinuous interface of the two fluids, a volume force is used which acts on fluid

elements lying within a transition region of finite thickness. The CSF formulation

makes use of the approach taken in numerical simulations that discontinuities

can be approximated, without increasing the overall error of approximation, as

continuous transitions within which the fluid properties vary smoothly from one

fluid to the other over a distance of O(h), where h is a length comparable to

the resolution of the computational mesh. Surface tension, therefore, is felt

everywhere within the transition region through the volume force included in the

momentum equations. A derivation is given in appendix A which shows that the

CSF formulation is equivalent in the limit of infinitesimal interfacial thickness h to

the classical description of these forces as boundary conditions at a discontinuous

interface, and supplements the previous derivation in the CSF references (Kothe

et al., 1991; Brackbill et al., 1992). It can be seen that this formulation is similar

to the approach taken in the finite element method in that the free surface forces

are included in the momentum integral residual equation (Georgiou et al., 1988;

Cuvelier and Schulkes, 1990; Malamataris and Papanastasiou, 1991). Also note

that interfacial surface tension is presumed constant, and surface tension gradient

induced stresses are not included in equations (3.2) and (3.3). Effects such as

surface tension gradient and surface dilatational and shear viscosities can, in

principle, be included, (Scriven, 1960) and may be important if surface active

agents are present.

Massless marker particles can be placed in the flow to follow fluid elements

if necessary for comparison with experimental results. Local velocities (up, vp)

of the marker particles are used to update the positions by use of the kinematic


5/23

37

relations:dr

dt= up,

dz

dt= vp (3.11)

For 2-D axisymmetric flows the streamfunction can be defined, as usual,

in order for the velocity to be divergence free (i.e., by satisfying the incompress-

ibility condition):

u =1

r

z, v =

1

r

r(3.12)

Then the streamfunction can be calculated from the following Poisson equation

(obtained by differentiating equations (3.12)), given the velocity field:

2

r2

+2

z2

= ru

z

v

r v (3.13)

These equations are to be solved subject to appropriate boundary conditions

discussed in section 3.6 (see, for example, chapter 5).

3.2 Momentum and Continuity Finite Differencing

The momentum equations are finite-differenced on a locally variable, stag-

gered mesh using the control volume approach, as illustrated in Figure 3.1. As

Figure 3.1 shows, the radial velocity ui+12

,j and axial velocity vi,j+12

are centered

at the right face and top face of each cell respectively, whereas the pressure, pi,j,

and marker function, Fi,j, are located at the center. The details can be found in

Nichols et al. (1980) and Hirt and Nichols (1981), while an introduction to control

volume approaches on staggered meshes can be found in Patankar (1980).

The equations (3.1) to (3.4) are solved for the variables u, v, and p at each

time step by an iterative method, a form of successive overrelaxation (SOR) that

will now be described. We may rewrite the equations explicitly in velocity and

implicitly in pressure for the (i, j)th cell with an Euler scheme from the currenttime n to the new time n + 1, accurate to first-order in time:

u

r+

v

z+

u

r

n+1= 0 (3.14)


6/23

38

zj

ri

Fictitious Cell Boundary

PhysicalBoundary

vi,j+1

pi,j

Fi,j

ui+1,j2

2

Figure 3.1: Staggered mesh.


7/23

39

un+1i+1

2,j

= uni+1

2,j

t

p

r

n+1

+t uur

vu

z

+1

rr

r

+zr

z

+rr

r+ gr +

F

r

n(3.15)

vn+1i,j+1

2

= vni,j+1

2

t

p

z

n+1

+t

u

v

rv

v

z+

1

zrr

+zzz

+zrr

+ gz +

F

z

n (3.16)

nrr = 2

urn

, nzr =

vr +

u

zn

, nzz = 2

vzn

(3.17)

The viscous terms in the momentum equations are calculated using central

finite differences accurate to second-order, taking into account the variable mesh.

In the original SOLA-VOF, the viscosity was assumed constant throughout the

entire domain and the viscous terms were calculated as in the Navier-Stokes

equations. The viscous term of the form 2v were then finite differenced.

However, a test problem showed that the tangential stresses were not equal atthe interface. In the present work the viscous terms were finite differenced using

the expressions (3.17), which allow for a variable viscosity in the interface region.

If the inertial terms are approximated by second-order accurate central

differencing, the algorithm is numerically unstable (Hirt, 1968). If they are

approximated by first-order accurate upwinding or donor-cell differencing,

the algorithm is stable but it might be overstable due to the presence of excessive

artificial viscosity (Hirt, 1968). In the SOLA-VOF method a linear combination

of central (HN = 0) and donor-cell (HN = 1) differencing is used to approximate

the inertial terms with controlling parameter HN (Hirt and Nichols, 1981). To

develop it consider using second-order accurate central differencing for the first


8/23

40

convective term in the radial momentum equation:

u

u

rn

= uni+12 ,j

xi+1

un

i+12

,jun

i12

,j

xi

+ xi

un

i+32

,jun

i+12

,j

xi+1

xi+1 + xi (3.18)

However, this will lead to a numerically unstable algorithm. We may use first-

order accurate upwinding where the derivative upwind of the point in question

is used. In the case of flow in the positive r direction (sgn(ui,j) = +1) this means:

uu

r

n

= uni+1

2,j

uni+1

2,jun

i12

,j

xi (3.19)

In the case of flow in the negative r direction (sgn(ui,j) = 1):

u

u

r

n= un

i+12

,j

un

i+32

,jun

i+12

,j

xi+1

(3.20)

These two equations are stable, but accurate to only first-order. We may combine

the above three equations into one for stability and accuracy:

uu

r

n

=

uni+1

2,j

xi+1

un

i+12

,jun

i12

,j

xi

+ xi

un

i+32

,jun

i+12

,j

xi+1

+HN sgn(ui,j)

xi+1

un

i+12

,jun

i12

,j

xi

xi

un

i+32

,jun

i+12

,j

xi+1

xi+1 + xi + HN sgn(ui,j) (xi+1 xi)(3.21)

where HN is the interpolation factor between the methods such that whenHN = 0 or HN = 1 central differencing or upwinding are recovered respectively.


9/23

41

3.3 Pressure Iteration

This section, which describes the pressure iteration at each time step,

follows the procedure described by Gerhart et al. (1992). Pressure does not

appear explicitly for the new time step in the continuity equation so it must be

recast so that pressure appears as a variable at the new time step. If we solve

the momentum equations, the velocity field does not necessarily satisfy continuity

because we would be using the old time step pressure field pn. Therefore we must

use continuity to correct the pressure field to the new values so that the velocitiesun+1, vn+1

will satisfy continuity at the new time level. The pressure is split

into an old time-level component and a correction such that:

pi,j pn+1i,j p

ni,j (3.22)

Then we can decompose the pressure derivatives:

p

r

n+1=

r

pn+1 pn + pn

=

p

r+

p

r

n(3.23)

pz

n+1=

z

pn+1pn + pn

= p

z+

pz

n(3.24)

Substituting these in the momentum equations gives:

un+1i+1

2,j

= uni+1

2,j

t

p

r

+t

1

p

r u

u

rv

u

z+

1

rrr

+zrz

+rrr

+ gr +

F

r

n(3.25)

vn+1i,j+1

2

= vni,j+1

2

t

pz

+t

1

p

z u

v

rv

v

z+

1

zrr

+zzz

+zrr

+ gz +

F

z

n(3.26)


10/23

42

Defining the following quantities:

Dn 1

r

r(ru) +

v

zn

(3.27)

Qnr

1

p

r u

u

rv

u

z+

1

rrr

+zrz

+rrr

+ gr +

F

r

n(3.28)

Qnz

1

p

z u

v

rv

v

z+

1

zrr

+zzz

+zrr

+ gz +

F

z

n(3.29)

we may rewrite the momentum equations as:

un+1i+1

2,j

= uni+1

2,j

t

p

r+Qnr (3.30)

vn+1i,j+ !

2

= vni,j+1

2

t

p

z+Qnz (3.31)

We may start the iterative method by calculating a first estimate of the velocities

with a fully explicit guess (p(1) = 0):

u(1)

i+12

,j= un

i+12

,j+ Qnr (3.32)

v(1)

i,j+12

= vni,j+1

2

+ Qnz (3.33)

For an improved guess we must include the pressure correction p(2):

u(2)

i+12

,j = un

i+12 ,j t1

p

r(2)

+Qn

r = u(1)

i+12

,j t1

p

r(2)

(3.34)

v(2)

i,j+12

= vni,j+1

2

t

1

p

z

(2)+Qnz = v

(1)

i,j+12

t

1

p

z

(2)(3.35)


11/23

43

To solve for the pressure correction, substitute these into the continuity equation

which we require to be satisfied at the new iteration:

D(2) 1

r

r ru +v

z(2)

=1

r

r

ru

(1)

i+12

,j rt

1

p

r

(2)+

z

v(1)

i,j+12

t

1

p

z

(2)

=1

r

r

ru

(1)

i+12

,j

+

z

v(1)

i,j+12

t

r

r

r

p(2)

r

t

z

1

p(2)

z

= 0(3.36)

So that a Poisson-like equation results:

r

1

p(2)

r

+

z

1

p(2)

z

+

1

r

p(2)

r=

D(1)

t(3.37)

This equation must be solved for p(2). We may use finite differences taking into

consideration the variable mesh to obtain:

D(1)i,j

t =

2

ri p

(2)i+1,j p

(2)i,j

(r)i+1 + (r)i

p(2)i,j p

(2)i1,j

(r)i + (r)i1

+2

zj

p(2)i,j+1 p(2)i,j

(z)j+1 + (z)j

p(2)i,j p

(2)i,j1

(z)j + (z)j1

+2

ri

p(2)i+1,j p(2)i1,j

(r)i+1 + 2(r)i + (r)i1

(3.38)

where mesh-weighted averages are used to obtain the quantities:

(r)i+1 + (r)i = i+1,jri + i,jri+1 (3.39)

(z)j+1 + (z)j = i,j+1zj + i,jzj+1 (3.40)


12/23

44

i,j = 1 (1 Fi,j) + 2Fi,j (3.41)

i,j = 1 (1 Fi,j) + 2Fi,j (3.42)

We have to solve this (linear) system (3.38) for p(2)i,j . We may use an SOR

iterative procedure where the off-diagonal terms are neglected (e.g., p(2)i+1,j 0)

so that the system may solved directly for the new p(2)i,j :

p(2)i,j =

D(1)i,j

2tri

1

(r)i+1+(r)i+ 1(r)i+(r)i1

+ 2tzj

1

(z)j+1+(z)j+ 1(z)j+(z)j1

(3.43)

where is the relaxation factor and the old divergence is:

D(1)i,j =

u(1)

i+12

,j u

(1)

i12

,j

ri+

v(1)

i,j+12

v(1)

i,j12

zi+

u(1)

i+12

,j+u

(1)

i12

,j

2ri(3.44)

We may now update u(1)i,j , v

(1)i,j and p

(1)i,j using (3.34):

p(2)i,j = p

(1)i,j + p

(2)i,j (3.45)

u(2)i+1

2,j

= u(1)i+1

2,j

+ 2tp

(2)

i,j

(r)i+1 + (r)i(3.46)

u(2)

i12

,j= u

(1)

i12

,j

2tp(2)i,j

(r)i1 + (r)i(3.47)

v(2)

i,j+12

= v(1)

i,j+12

+2tp

(2)i,j

(z)j+1 + (z)j(3.48)

v(2)

i,j12

= v(1)

i,j12

2tp

(2)i,j

(z)j1 + (z)j(3.49)

Iteration is continued until maxDki,j < and the variables are ready for the next

time step. The SOR parameter must meet the criterion:

1 < 2 (3.50)


13/23

45

As is discussed in Ciarlet (1991), this parameter is adjusted to obtain the minimum

number of iterations by sensitivity studies depending upon the particular case of

interest.

3.4 Kinematic Condition

The kinematic condition for the advection of the VOF function can be

written in finite-difference form as:

Fn+1i,j = Fni,j t

u

F

r+ v

F

z

(3.51)

The convective terms may be rearranged to divergence form, minus a divergence

correction term:

uF

r+ v

F

z=

1

r

r(F ru) +

z(F v) F

1

r

r(ru) +

v

z

(3.52)

Substituting (3.52) into (3.51) we get the divergence part:

Fi,j = Fn

i,j t

1

r

r(rF u) +

z(F v)

n(3.53)

which is then updated with the correction:

Fn+1i,j = Fi,j + tFni,j

1

r

r(ru) +

v

z

n+1(3.54)

This sequence insures conservative advection of F (Kothe et al., 1991). Ordinarily

the correction in (3.54) would be zero due to continuity, but it has been found

desirable to include it numerically (Brackbill et al., 1992) because although it is

small, but is non-zero and of the order of t. In the VOF code, this equation is

utilized at the end of the pressure iteration using the new un+1, vn+1 values.

A brief description of the Hirt and Nichols (1981) advection algorithm,

which takes special care in preserving the discontinuous nature of F, is now

given; further details can be found in Hirt and Nichols (1981). The divergence


14/23

46

equation for F can be finite differenced for the radial advection term in terms of

an upstream donor(d) cell at (id,j) and a downstream acceptor(a) cell at (ia,j) (if

un+1i,j > 0,idm = i1, id = i,ia = i + 1, otherwise, ia = i,id = i + 1,idm = i +2):

Fid,j = Fn

id,j 1

rid

ri+12

t (F u)

xid(3.55)

Fia,j = Fn

ia,j +1

ria

ri+12

t (F u)

xia(3.56)

where the amount of F fluxed across the cell face in t is:

t (F u) = min

Fiad,j

un+1

i+12

,jt

+ CFr, Fid,jxid

(3.57)

with the correction factor:

CFr = max

Fidm,j Fiad,j

un+1i+12

,jt

Fidm,j Fid,jxid, 0

(3.58)

In these equations, the subscript iad = ia if the fluid is more vertical than

horizontal, otherwise iad = id. The min feature prevents the fluxing of more

of fluid 2 (F = 1) from the donor cell than it has to give, while the max feature

accounts for an additional flux, CFr, if the amount of fluid 1 (F = 0) to be fluxed

exceeds the amount available. A subsequent pass is made for the axial convective

term:

Fi,jd = Fn

i,jd t (F v)

zjd(3.59)

Fi,ja = Fn

i,ja +t (F v)

zja(3.60)

t (F v) = min

Fi,jad

vn+1

i,j+12

t

+ CFz, Fi,jdzjd

(3.61)

CFz = max

Fi,jdm Fi,jad vn+1i,j+1

2t Fi,jdm Fi,jdzjd, 0 (3.62)

Finally, the divergence correction is added:

Fn+1i,j = Fi,j + tFni,jD

n+1i,j (3.63)


15/23

47

3.5 Surface Force Terms

This section is a brief summary of the CSF calculation of the surface force

terms in the momentum equations (Kothe et al., 1991; Brackbill et al., 1992). The

surface force in the momentum equations is:

Fsv i,j = i,j ni,j (3.64)

and is the curvature of the liquid-liquid interface which is calculated as:

= (n) =1

|n|

n

|n|

|n| (n)

(3.65)

Brackbill et al. (1992) have rewritten the curvature in terms ofn and |n| (equation

(3.7)) instead ofn (equation (3.6)) to insure that principal contributions to finite-

difference approximation of come from the center of the transition region rather

than the edges.

The normal vector at cell corners is obtained by finite differencing, for

example at the top right corner of cell (i, j):

nr i+12 ,j+12 =(Fi+1,j+1 Fi,j+1) zj + (Fi+1,j Fi,j) zj+1

(zj + zj+1) ri+12

(3.66)

nz i+12

,j+12

=(Fi+1,j+1 Fi+1,j ) ri + (Fi,j+1 Fi,j) ri+1

(ri + ri+1) zj+12

(3.67)

The calculation of the first term in the curvature expression is (|n| =

n2r + n2z):

n

|n|

|n| =

nr|n|i,j

|n|

r i,j+

nz|n|i,j

|n|

z i,j=

nr|n|

2i,j

nrr

i,j

+

nrnz|n|2

i,j

nrz

+nzr

i,j

+

nz|n|

2i,j

nzz

i,j

(3.68)

The cell-centered partial derivatives follow from knowledge of n at the corners, for


16/23

48

example :

nrr

i,j

=

nr i+1

2,j+1

2+ nr i+1

2,j1

2

nr i12

,j+12

+ nr i12

,j12

2r

i

(3.69)

The cell-centered normal is the average at the corners:

nr i,j =1

4

nr i+1

2,j+1

2+ nr i1

2,j+1

2+ nr i+1

2,j1

2+ nr i1

2,j1

2

(3.70)

nz i,j =1

4

nz i+1

2,j+1

2+ nz i1

2,j+1

2+ nz i+1

2,j1

2+ nz i1

2,j1

2

(3.71)

The second term in the curvature expression is calculated:

(n)i,j =1

ri

r(rnr)i,j +

nzz

i,j

(3.72)

where:

1

ri

r(rnr)i,j =

ri+12

nr i+1

2,j+1

2+ nr i+1

2,j1

2

ri1

2

nr i1

2,j+1

2+ nr i1

2,j1

2

2riri

(3.73)

The face centered values of Fsv i,j at the right and top faces are required for the

momentum equations and are calculated from cell centered values:

Fs v r i+1,j =riFs v r i+1,j + ri+1Fsvr i,j

ri + ri+1(3.74)

Fsv z i,j+1 =zj Fsv z i,j+1 + zj+1Fsv z i,j

zj + zj+1(3.75)

The CSF method has the ability to use the smoothed or mollified VOF

function Fi,j introduced in equation (2.45) for the calculation of the curvature

i,j in the volume force equation (3.64), which is different from the unsmoothed

function Fi,j used to calculate the normal vector in equation (3.64). This enables

the algorithm to calculate a smoother curvature for accuracy and, we have found,

for decreased number of pressure iterations. The smoothed VOF function is


17/23

49

computed by convolving F with a B-spline of degree l (de Boor, 1967; Brackbill

et al., 1992), (l) (|x x| ; h), (with l 2 in this dissertation) where (l) = 0 only

for |x x| < (l + 1) h/2 = 3h/2. The smoothed VOF function is given by:

Fi,j =k

i,j=1

Fi,j (l)

xi,j xi,j; h

(l)

yi,j yi,j; h

(3.76)

where the sum gathers contributions from the nine values (for l = 2 in 2-D) ofFi,j

within the support of (2). In our case this formulation becomes simply:

Fi,j =9

16Fi,j +

3

32(Fi+1,j + Fi1,j + Fi,j+1 + Fi,j1)

+

1

64 (Fi+1,j+1 + Fi+1,j1 + Fi1,j+1 + Fi1,j1)

(3.77)

This formula may be applied iteratively by multiple passes through the mesh for

increased degrees of smoothing. Our experience has shown that one to three passes

are optimal, with most calculations carried out with one pass.

3.6 Boundary Conditions

The equations are solved with appropriate boundary conditions in axisym-

metric coordinates (r, z), with the boundary conditions as follows. The conditions

which we are concerned with are no-slip, free-slip, and continuative. For the solid

walls no-slip conditions are used:

u = v = 0 (3.78)

For example, at a bottom wall:

ui+12

,1 =ui+1

2,2z1

z2(3.79)

vi,32

= 0 (3.80)

where the z1 and z1 appear as interpolation factors in equation (3.79) so that

the radial velocity may be exactly zero on the boundary. For the axis of symmetry


18/23

50

at r = 0 (same as a free slip wall):

u = 0,v

r= 0 (3.81)

u32

,j = 0 (3.82)

v1,j+12

= v2,j+12

(3.83)

For a continuative boundary at the bottom of the mesh, for example, there is no

change in the axial direction:u

z=

v

z= 0 (3.84)

ui+12

,1 = ui+12

,2 (3.85)

vi,32

= vi,52

(3.86)

For problems involving the effects of wall adhesion characterized by a given

equilibrium or dynamic contact angle, , such as the evaluation of the interface

separating two stationary immiscible fluids, this condition can be imposed by

requiring that the unit normal to the interface at the wall contact line, n, satisfy

(Kothe et al., 1991):

n = nw cos + tw sin (3.87)

where nw and tw are the unit normal and tangent to the wall at the contact line

respectively.

In the jet simulations, special attention was given to force the F function to

intersect, or be pinned, to the nozzle lip as shown in Figure 5.1, while using the

Hirt-Nichols acceptor/donor method (Nichols et al., 1980). This was accomplished

by the simple expedient of not allowing fluid 2 to flux across the top of the nozzle

(if it so desired), thus imitating the experimentally observed condition. Note that

this procedure does not influence the steady-state solution reported in section 5.3

since the corresponding condition that activates it is not applicable there. Thus, in

practice, this condition was only necessary in order to stabilize the initial transient

during startup.


19/23

51

3.7 Numerical Stability

The numerical difference equations are subject to linear numerical stability

conditions that are detailed in Hirt and Nichols (1981). Material cannot move

more than one cell in one time step yielding the Courant condition:

t < min

riui+1

2,j

,zivi,j+1

2

(3.88)

Second, momentum must not diffuse more than about one cell in one time step:

1t

1

Documents

ATOMIZAÇÃO_phd3