A HYBRID FINIT-ELEMENT_VOLUME-OF FLUID METHOD FOR SIMULATING FREE SURFACE FLOWS AND INTERFACES

 

 

Free surface flows and interfaces between two immiscible fluids or materials with different phases are observed in many natural and industrial processes.  Different numerical techniques are developed to simulate these flows.  However, due to the complexity of the problem, each technique is tailored to a particular category of flows.  For instance, boundary integral techniques 1-3  are mainly used for simulating  inviscid irrotational flows and the limiting case of zero Reynolds number.

 Finite Element Methods (FEM) and Finite Difference Methods (FDMS) are potentially

applicable to generalized Navier-Stokes equations, however, they have to be coupled with a technique to track the advecting fluid boundaries and interfaces.  The difficulty in the interface tracking is inherently related to the complexity of its topology.  Therefore,

techniques which can handle small surface deformations fail when applied to large interface distortions.  For the simulation of  the former category of flows (small surface deformations), FEM is more popular.

Here, the fluid boundary is described by a set of fixed 4 or deforming 5-12 mesh, the location of which is obtained by either an iterative procedure or the Lagrangian movement of the interface nodes.

This results in the simultaneous calculation of the position of the free surface and the field variables at the new nodal positions.

  Boundary fitted orthogonal coordinates 13-15 and Lagrangian techniques 16,17 are also used to follow the advecting liquid interfaces. These techniques are confronted with difficulties when applied to large surface deformations, surface breaking and merging.

 

More versatile and robust category of techniques for free surface flow modeling are the front tracking methods. Here, an extra set of parameters is used to trace the fluid boundaries.  Front tracking techniques are divided into two groups: surface-tracking

and volume tracking methods.  In general the former class gives a more accurate description of the free surface, but the latter class can handle complicated liquid regions more easily.  In surface tracking methods the position of the free surface is described

in a direct way; either by specifying a set of marker points located on the free surface 18 , or by introducing a height function which explicitly describes the free-surface position 19.

There are several problems associated with the surface tracking methods.

The main problem is that  the marker points will be nonuniformly distributed as the interface evolves. Also, relatively high computer storage is needed to maintain the interface continuous and smooth.

 

Volume tracking techniques define a tracer  to follow the whole fluid region.  The two commonly used techniques are MAC and VOF techniques.  In Marker And Cell (MAC) technique 20-30 hundreds of massless marker particles are added to the fluid.  These particles are then advected in Lagrangian sense using the average of Eulerian velocities in their vicinity. In Volume Of Fluid (VOF) technique 24-32 a volume fraction parameter, f is described for every of the Eulerian grid cells. A cell is assumed to be filled with liquid if f = 1,  it is considered empty if f = 0 and partially full if 0 < f < 1.  Different methods are developed to advect the volume fraction field and to reconstruct the fluid surface.

VOF based techniques can handle the most complex free surface flow problems.  Surface breaking and merging can be  treated with this technique since the flow field calculations are decoupled from free surface location identification.  Here, knowing the initial fluid

surface location and the velocity field, the surface is advected using the volume fraction field; a  surface is reconstructed based on the newly calculated volume fractions; and finally the conservation equations are solved for the new fluid domain to find the  velocity field.  This procedure is repeated throughout the advection of the fluid region.  There is generally no need for iteration between the flow field calculation and interface location identification. However, presently there is a major shortcoming to these techniques.  Although the interface itself can locate inside a cell, the governing equations for the field variables are only solved for the whole cell. This results in significant inaccuracies in the treatment of the interface viscous stresses and the surface tension forces.

The present paper describes a new technique for the simulation of free surface flows and interfaces by combining FEM with deforming interface mesh and a VOF based technique called FLAIR31 for surface advection. The combination of these two will solve the major shortcomings of each technique. The VOF based technique will allow the simulation of interfaces with large deformations, even with surface merging and breaking. In addition, the FEM with deforming interface mesh allows accurate implementation of the boundary conditions.  In section 2 the governing equations of free surface viscous flows are provided. Their finite element formulation is given in Section 3 and the mesh generation and advection techniques are described in Section 4 and Section5,  respectively. Section6  and 7 provide detailed description of the simulation of the surface breaking and surface merging by two example problems, and concluding remarks are made in Section8.

 

2.GOVERNING EQUATIONS

 

The problem considered here is the laminar flow of a viscous incompressible fluid. The assumptions of Newtonian fluid with constant properties are also applied. Therefore, flow field is governed by continuity and momentum (Navier-Stokes) equations:

                                       ,                                           (1)

 

 

               (2)

 

               (3)

 

where  for 2D Cartesian andfor axisymmetric flows. Equations (1)-(3) are written in non-dimensional form using the following non-dimensional parameters:

 

                

 

                

 

,        

 

 

 

 

                                                                                                   

                                                                                                                                                                                                                                                                                                                                                                             where ,  and   represent the density, viscosity, and surface tension coefficient of the fluid, respectively. Lengths and velocities are normalized with the characteristic scales L, and U.

The boundary conditions are given by:

                                          On S1,                       (4)

 

            (5)

 

 

where S 1 and S 2 are parts of the boundary with Dirichlet and

Neumann boundary conditions, respectively.

 

Fig1

and  bar denote z- and r- components of the total surface traction and nz and nr

denote direction cosines of the unit outer normal to the surface S2 . On the free surface,  and  are the components of the surface tension which is inversely proportional to the radius of curvature of the surface, R'c. Therefore,

 

                          (6)

 

The non-dimensional radius of curvature of an axisymmetric surface is defined as

 

,

 

 

 

with R 1 and R 2 , * being the principal radii of curvature of the surface, defined as:

 

                               

 

No contact between the fluid and solid surface is involved in the problems considered in this paper, therefore the contact angle  is not modeled here.

 

 

3. FINITE ELEMENT FORMULATION

 

The selection of the finite element solution to the conservation  equations (1) - (3)

is subject to two major restrictions. First, because of the requirements of the advection technique which is basically designed for finite difference schemes the finite element solution should be in terms of primitive variables based on the linear quadrilateral elements.

Furthermore, due to the significance of the surface tension effects in free surface flows the model must be capable of handling the pressure, velocity, velocity gradient, and stress boundary conditions, directly.

 

A finite element model which complies with these requirements is the penalty function formulation, in which the continuity equation is absorbed in the momentum equation by representing the pressure as:

 

                   (7)

 

 

Where (0(10 9)) is a large number depending on  and Re33. Substituting the pressure equation (7) into the momentum equations (2) and (3), and applying the Galerkin's method the following matrix equations are obtained:

 

 

 

 

 

 

 (8)

 

Here

 

 

                                     (9)

 

          (10)

 


           (11)

 

 

          (12)

 

               (13)   

 

 

 

          (14)

 

 

          (15)

 

                                    (16)

 

where the Ni  are the shape functions for the isoparametric quadrilateral elements in the domain with the boundary .

Integrals (9)-(12) and (16) are evaluated by 2 * 2 Gauss-Legendre integration rule.

The [L]-matrix should be singular for the  penalty function approach to be successful.  A

 reduced Gauss integration is used to evaluate the [L]- matrix given by integrals (13)-(15).  The surface traction resulting from the surface tension effects is directly integrable over the free surface using equation (16).

Equation (8) is integrated in time using an implicit forward differencing

finite difference technique.

 

 

4. MESH GENERATION

 

Mesh generation for a finite element problem has proven to be an important part of the overall solution procedure.  A large effort is devoted to this issue in order to generate the most efficient mesh for every particular problem, using various techniques such as Laplace's equation34, and transfinite mapping.35-37  These techniques generate the interior nodes based on the predefined boundary nodes and in general there is no control on the location of the generated interior nodes. In case of highly distorted domains, it is necessary to divide the whole domain into simple subdomains. However, where

moving boundaries are involved the solution domain experiences successive evolutions in the course of time, resulting in the redefinition of subdomains. This can be tedious  for very complicated boundary evolutions, such as the cases encountered in breaking or merging of two liquid regions.

These limitations along with the requirements of the advection technique for moving  interfaces makes the use of automatic mesh generators inconvenient.

 

In the present work rectangular finite element mesh is used in order to utilize the existing interface advection techniques.  Implementation of the finite element method on the Eulerian grid for solving the moving boundary problems is well developed4,38 and it can be applied to the present model.

In order to describe the mesh generation method used in this work consider the closed curve located on the Eulerian grid shown in Fig 1. The master element chosen for the finite element mesh is the isoparametric linear quadrilateral element which is also suitable for the penalty function method used in the momentum solver. Hence, each full cell will be considered as one element and we begin numbering the elements with the first full cell located on the left of the lowest row, and then continuing toward the full cells on the right of this cell, and repeating the same procedure for the other rows.  The

intersecting points of grid lines are the nodes and they are numbered counterclockwise for each element starting from the lower left node,  as shown in Fig 1.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure2.

 

Once the full cells are completely transformed into elements, the wet cells must be employed to define the rest of the mesh.  Depending on the intersection points of the interface line with the faces of the cell, nine possible cases for a wet cell can be recognized.  These cases are shown in Fig. 2 with the wet cell located at the center of a

3 * 3 cell unit. Note, that the unit is first rotated so that the heavy-side (the side which contains more mass than the other three sides) is on the lower right corner  as shown in Fig. 2. These cases can be distinguished by examining the four cells located on the middle of each side of the unit; namely, cells I-IV as shown in Fig. 2.  All the possible cases are summarized in Table 1. In cases (5) to (8) of Fig. 2 the four nodes are well defined, while in cases (1) to (4) a fourth node has to be defined between the two interface  nodes. Case (9) includes five nodes and cannot be considered as one element. Therefore, an extra node is defined between the two nodes which are located on the free surface.  This transforms the wet cell into two elements.

 

The above procedure may result in very small elements on the free surface adjacent to larger rectangular elements. These small elements are encountered when the free surface crosses an interface cell very close to its interior node(s) as seen in Fig. 3(a). Therefore, the mesh should be altered in this region to redefine the spacing between the nodes. This can be accomplished by moving the interior nodes of the surface cells away from the interface. Figure 3b shows the altered mesh which now has a better aspect ratio for the surface elements. The inward displacement of the interior node is proportional

to its distance from the interface, the closer the node is to interface the more it is displaced.

 

 

 

 

 

 

 

 

(a)

 

 

 

 

 

 

 

 

Figure 3                                     (b)

 

 

 

Table I. Case-distinguishing criteria for wet cells shown in Figure 2

 

Case

Cell I

Cell II

Cell III

Cell IV

(1)

Wet

Wet

Empty

Empty

(2)

Wet

Full

Empty

Empty

(3)

Full

Wet

Empty

Empty

(4)

Full

Full

Empty

Empty

(5)

Wet

Full

Empty

Wet

(6)

Full

Full

Empty

Wet

(7)

Full

Wet

Wet

Empty

(8)

Full

Full

Wet

Empty

(9)

Full

Full

Wet

Wet

 

 

 

5. SURFACE ADVECTION

 

 

Once the velocities are obtained the interface has to be advected. We have implemented the Flux Line-segment model for Advection and Interface Reconstruction (FLAIR) originally developed by Ashgriz and Poo31 for 2-$D$ Cartesian advection, and later extended by Mashayek and Ashgriz32 to Axisymmetric flows.

In this technique a volume fraction, , is defined for each cell such that: for the cells which are filled with liquid; for the empty cells; and  for the interface or wet cells.  Therefore, the liquid domain is discretized and represented by a set of volume fractions which we will refer to as the f- field.  FLAIR technique is designed  to advect the f- field by knowing the surface velocities and to reconstruct the surface topology based on the new $f$-field.

The principle of the surface advection by FLAIR is illustrated in Fig. 4.  It is assumed that the surface can be represented by a set of line-segments drawn at the boundary of the two neighboring cells.

Therefore, each line-segment can be described by: y=ax+b, where constants a and b can be determined based on the known volume fractions31  and .  Once the surface is

reconstructed by the combination of all the surface line-segments, it is advected by calculating the mass flux in and out of each cell.  The mass flux in each time step  moving from one cell to the neighboring cell can be calculated by using the velocity at the boundary of the two cells, ui,j , as follows:

 

 

 

 

Here x2 – x1 represent the distance the mass is advected in time  i.e. x2 - x1= ui,jTherefore, the mass flux in time step  for a surface with the orientation shown in Fig. 4 and x2 = 1 is

 

 

Other surface orientations and advection in the vertical direction can be considered in the similar fashion31,32 . The new cell volume fraction is then obtained by adding

the old value to the net in-flux of the volume fraction to each cell.

 

This advection technique has been basically designed for the finite difference staggered grids where the velocities are defined at the middle of each cell face.

However, finite element solution of the flow provides the velocities at the nodes of the element which are not suitable for this advection technique. Mashayek and Ashgriz32 have shown that the advection criterion can be deduced from satisfying the continuity of the flow through the faces of each cell during every time step. We use the same idea and calculate a mean velocity from the known velocities at the two nodes attached to the ends of one cell face, such that the same mass of fluid passes through that face.

 

Consider the axial direction in an axisymmetric flow.  Each cell in r-z plane is mapped on an isoparametric master element in  plane, as shown in Fig. 5. The mean velocity,

 across the face 2-3 is defined by:

                    (17)

 

 

 

where the velocity and radius are defined as:

 

            (18)

 

 

with  Ni being the shape functions for the isoparametric quadrilateral elements, defined as follows:

 

      (19)

 

 

 

 

Figure 5.

Then:

 

 

 

The side 2-3 of the cell is mapped onto the side 2-3 of the master element where  and  Therefore:

 

 

After substituting for N2, N3, and their derivatives we find that:

 

      (20)

 

The integral on the left hand side of equation (17) can be written as:

 

 

                         (21)

 

Combining equations (17), (20), and (21):

 

         (22)

Applying the same analysis on radial direction shows that is just the average of the v- velocities at the ends of the cell face:

 

                 (23)

In order to clarify the connection between the various steps employed in the FEM-VOF technique, consider the  oscillation of a liquid drop as shown in Fig. 6.  Owing to the symmetry a quarter of the drop, as depicted by the shaded area in Fig. 6(a), is analyzed. By assigning the value of the volume fraction f to every cell this area is converted to a set of numbers as shown in Fig. 6(b).

 The interface shown in Fig. 6(c) is reconstructed based on this f-field by using the FLAIR algorithm31,32 .Then, by implementing the method explained in Section 4, the FEM mesh is generated as shown in Fig. 6(d) with thick lines. After the momentum equations are solved the velocities are obtained at the nodal points of the mesh. These velocities which are shown in Fig. 6(e), are then transformed to the velocity field depicted in Fig. 6(f) by applying equations (22) and (23). In the next step the new f-field is obtained through advecting the old f-field by implementing FLAIR.

The new f-field is shown in Fig. 6g and it is used to reconstruct the new interface  depicted in Fig. 6(h). This will complete the sequences employed in the first time step and in continuation the steps of Figs. 6(d) to 6(h) will be repeated for each time step.

 

 

                 

 

 

 

 

 

 

 

 

 

 

          

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 



 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 6.

 

The FEM along with the VOF method for interface advection has several important advantages. The first advantage of FEM-VOF is in modeling of the surface breaking and merging.  On the fluid regions where surface breaking or merging is about to occur the accuracy of the solution can be greatly enhanced if very small fluid thicknesses at the breakup point or very small distances between the two merging surfaces are considered.  This can be easily achieved by FEM-VOF method since the interface boundary condition can be applied accurately.

The previous VOF techniques combined with the finite difference method (FDM-VOF) have to resort to interpolation technique for the implementation of the interface boundary conditions. Since the liquid surface is continuously moving, new cells are being formed.

Therefore, in order to complete the boundary conditions needed for the momentum equations new velocity components have to be defined for the newly  formed interface cells.  This is achieved by applying the continuity equation to these interface cells 31. In many of the VOF  techniques  when the thickness of the liquid comes within the order of one cell size this method fails since it results in double definition of the boundary velocities. Therefore, the criterion for the liquid breakup is set to be when the local liquid

thickness reaches one cell size.

Another advantage of the FEM-VOF over FDM-VOF is that the former defines a computational element for the part of the wet cell which actually contains fluid, while in

the FDM-VOF techniques the wet cells are treated in the same way as the full cells using the whole rectangular cell for the discretization of the momentum equation.  This allows accurate calculation of the boundary conditions without resorting to time consuming interpolations as in FDM-VOF.

 By implementing the FEM in momentum solver this technique is capable of handling the pressure and velocity boundary conditions accurately. The techniques which are based on the FDM are unable of applying these boundary conditions directly. In the most recent improvements which are presented by Unverdi and Tryggvason39 and Brackbill et al.40

the surface tension effects are considered as a body force distributed in the cells close to the interface. This may result in large errors when the thickness of the fluid is small, which is unavoidable in breaking and merging of liquids.  In the present technique, the interface is treated as a discontinuity which has better resolution in comparison to the techniques that smear the interface over some finite region.

Finally, FEM-VOF significantly increases the computational efficiency, because a smaller number of cells will be sufficient to model large surface deformations.  The number of cells in each direction used in this technique is two and a half times smaller than what is used in the FDM-VOF technique of Ashgriz and Poo31 in order to have the

same accuracy.  The problem of area change reported by Ashgriz and Poo31 is almost completely resolved by keeping the divergence of the velocity less than 10-9 for each element in each time step.

The momentum solver is fully implicit and therefore imposes no significant limit on the selection of time increment. Consequently, the time increment will be determined by the advection technique which is less restricted than the common explicit or semi-implicit momentum solvers. This will result in substantial decrease in the total number of time steps and will increase the overall efficiency and accuracy of the method.

 

 

6.  SURFACE BREAKING

 

As mentioned earlier, in FEM-VOF the continuous movement of the interface line inside a wet cell results in the formation of different elements at different time steps.  There is no limit on the size of these elements, therefore the thickness of the liquid can be reduced

to a small fraction of the cell size.  The advantage of this capability is illustrated by the surface movement of an unstable capillary jet as shown in Fig. 7. The initially disturbed surface is indicated in Fig. 7(a).  We have used two cells for the initial radius and

16 cells for one wavelength of the jet. Implementing FEM-VOF technique allows us to accurately decrease the thickness of the jet as much as it is desired.  For instance, the curve depicted in Fig. 7(b) shows the jet surface which has a minimum radius of about 5% of the undisturbed initial jet radius. This technique is implemented to simulate  the breakup of a capillary jet.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 7.

 

   Rayleigh,41 for the first time, obtained a linear analytical solution for the instability of inviscid jets in vacuum. Later, the analytical analysis extended to the nonlinear inviscid jets by Yuen42 and linear viscous jets by Chandrasekhar43, however the nonlinear

instability of viscous jets has not been the subject of analytical investigations. The experimental and numerical studies cover a broader range and are reported by a large number of researchers such as the works done experimentally by Goedde and Yuen44 .and

Vassallo and Ashgriz45, and numerically by Mansour and Lundgren46.

 

As a test problem for the FEM-VOF technique we have considered the breakup of an infinitely long jet in vacuum, subject to a periodic surface disturbance in the form of:

 

 

where and k are the amplitude and the wave number, respectively. The undisturbed radius of the jet is selected as the length scale, and the time scale is obtained by equating the Weber number to one. In order to validate the results, we consider the breakup of a glycerin-water jet with k=0.262 and Re=2.045, which has been experimentally studied by Goedde and Yuen44.  The results of our solution with are shown in Fig. 8 in different snapshots. Initially, the jet begins to swell at $ and necks down at z=0. However, after a short time interval the point having minimum radius starts moving toward the points .

Finally, the jet breaks at  leaving a spherical drop and a cylindrical ligament.

 

Figure 9a shows the variation of the amplitude of the swell and the neck regions of the jet with time.  The nonlinear effects are clearly seen close to the breakup point, where the amplitude of the swell region becomes flat and the neck region reduces rapidly.  To evaluate the accuracy of the calculations, we have compared our numerically calculated  growth rates with Goedde and Yuen's 44 experimental data and Chandrasekhar's43 analytical results.

 

 

 

Figure 8.

 

The same procedure, which is used by Goedde and Yuen44 for the calculation of the growth rates from the experimental data, is also used here. The exponential growth rate, , is defined by: , where  is the wave amplitude with the value of  at . Then:

.

 

Therefore,  is the slope of the curve versus time when this curve is approximated by a line.  We have observed the same variation of the growth rate as reported by Goedde and Yuen.  The exponential growth rate is not constant at any particular point of the jet surface, however, the exponential growth rate of the difference between the swell and the neck is nearly constant (except close to the breakup point) as shown in Fig. 9(b).  The numerically calculated average growth rate for the swell point is , and for the neck point is .

The growth rate for the difference between the amplitude of the swell

and neck points is  .  The analytically calculated growth rate for a viscous jet43

 for the same Re number is .

 Goedde and Yuen's44 experiments for approximately the same conditions resulted in a growth rate of  Considering the experimental error, which was estimated44 to be about 11%, the coarse numerical mesh used in our simulation, and that the initial disturbance amplitudes are not the same (initial disturbance amplitude for the experimental results are not reported), the observed differences between the experimental and the numerical results are acceptable.

 

 

 

 

 

 

 

 

 

 

 

Figure 9.

 

 

7.   SURFACE MERGING

 

 

The ability to accurately implement the boundary condition becomes more physically attractive when the collision of the two liquid drops in non-vacuum surrounding is concerned. In this case the fluid trapped between the drops prevents merging and  flat regions are developed on the surface of the drops which are about  to collide. Thickness of the fluid between drops is decreased due to the pressure generated by drops and merging becomes possible if the initial momentum of drops is capable of producing a pressure high enough to rupture the trapped liquid film, otherwise drops are bounced apart. In a real collision process the thickness of the film  may be reduced to a level comparable to molecular size and still the merging be prevented. Therefore, the margin of bouncing/merging collisions is not accurately found by implementing techniques that

model the merging phenomenon in the level of cell size.

 

The complex problem of  collision of two viscous liquid drops is chosen to demonstrate the capability of the FEM-VOF method in handling the merging of two liquid surfaces.  Because of the complexity of this problem majority of the published literature in this area

have been limited to the experimental research - Ashgriz and Poo47 provide a detailed count of the previous works in this field up to 1990.  Limited  numerical simulations have been attempted by Poo and  Ashgriz48 using FLAIR algorithm.  However, their simulation was for the collision of two liquid cylinders or two dimensional drops.

 

Here we present the binary collision of 3D drops in axisymmetric coordinates.  This will limit us only to the head-on collisions.

Figures 10(a)-10(h) show the time evolution of collision of two drops with Re=50 and We=10 - Re and We are evaluated based on the initial drop diameter and the relative velocity between the two drops which is equal to 6 for the present calculations. The initial

condition is shown in Fig. 10a, where only two elements - one element on each side of the axis of symmetry - are set to touch. This condition is not necessary, and the drops could be allowed to travel towards each other from a finite distance.  The velocity vector on each node is plotted on all figures to show the flow field and indicate the number of the elements in the system.  Size of the vectors are scaled relative to the initial velocity and the ratio is given at the top-right corner of figures by x.

Note that the number of the elements can change as it has in Fig. 10(b).

 

 

 

 

 

 

 

Figure 10

 

 

A very interesting transient flow field is observed during the collision of the two drops. Shortly after the collision (t= 0.09, Fig. 10(b)) the flow at the back side of each drop slows down while the top and bottom sides of the drops continue their motions. This results in the eventual formation of a shape shown in Fig. 10d.  At the collision surface the stagnation type flow pushes the interface to the sides  generating a combined drop as shown in Fig. 10(f).

The inertia forces keep pushing the liquid to the sides and surface of the drop moves inward around the axis of symmetry as seen in Fig. 10(e). Therefore, the surface tension forces are no longer in favor of pushing the liquid away from the axis. The maximum

inward displacement of the surface is reached in situation shown in Fig. 10(f) which indicates the initiation of the backward flow along the axis. The region covered by the reversed flow expands towards the sides of the drop (Fig. 10(g)) until finally surface tension and viscous forces overcome the initial momentum of the drop and the fluid flow is completely reversed. The drop tends to become spherical with identical curvature everywhere along the boundary. Without the viscous effects the drop will oscillate

infinitely, however the presence of viscosity will finally bring it to rest. Competition of surface tension and viscosity during the intermediate stages of drop deformation may result in various modes of oscillations.

 

There is no reported experiments  for the low Re number drop collision to compare the numerical results with.

The Re numbers for the water drop collisions of Ashgriz and Poo47 are an order of magnitude larger than the ones reported here, and they could not be simulated. However, the qualitative comparisons show the same evolutionary process.  The two examples presented here use a very coarse mesh which results in some surface roughness.  The mesh size is kept large to show the actual interface lines and also the versatility of the

method.  The surfaces can be made smooth by reducing the mesh size.

 

 

 

8. CONCLUDING REMARKS

 

A new technique is developed for simulating the free surface flows with large deformations. The Finite Element Method (FEM) is used to solve for the flow equations and a Volume of Fluid (VOF) technique is implemented for surface determination and interface advection.

Linking between FEM and VOF is provided by developing a new mesh generation technique. Detailed comparison to Finite Difference based VOF techniques is carried out by considering two most commonly encountered problems in free surface flows, namely liquid breaking and merging.  As compared to the other FEM techniques currently in use to simulate free surface flows, FEM-VOF may be less accurate, however, the other FEM techniques have difficulty in modeling the breakup and merging of liquids during a single flow simulation. Up to the knowledge of the authors this is the first attempt made to combine finite element method with volume of fluid techniques. Therefore, modifications are still expected to be made in order to improve the versatility and applicability of the technique. Mesh generation, for instance, can be modified to obtain more efficient mesh by providing better aspect ratios for the surface elements. Also loss of accuracy due to the velocity interpolations might be decreased by developing a new advection technique that utilizes nodal velocities rather than facial velocities.

 

 

 

 

9. REFERENCES

 

1.    R.W. Yeung, ‘Numerical methods in free surface flows,’ Ann. Rev. Fluid Mech. 14, 395 (1982).

2.     A.S. Geller, S.H. Lee, and L.G. Leal, ‘The creeping motion of a spherical particle normal to a deformable interface,’ J.  Fluid Mech. 169, 27, (1986).

3.     J.M., Rallison,  and A. Acrivos, ‘A numerical study of the deformation

and burst of a viscous drop in an extensional flow,’  J. Fluid Mech. 89,   191, (1978).

   4. K.J. Bathe and M.R. Khoshgoftaar, ‘Finite element free surface seepage analysis without mesh iteration,’ Int. J. Num. Analyt. Mech. Geomech., 3,   13, (1979).

   5. R.I. Tanner, R.E. Nickell and R.W. Bilger, ‘Finite element methods for the solution of some incompressible non-Newtonian fluid mechanics problems with free surfaces,’

 Comput. Meth. Appl. Mech. and Eng., 6, 155, (1975).

  6. H. Saito and L.E. Scriven, ‘Study of coating flow by the finite element  method,’

  J. Comput. Phys., 42, 53, (1981).

 

   7. H. Ettouney and R.A. Brown, ‘Finite element method for steady

solidification problems,’ J. Comput. Phys.,49, 118, (1983).

   8. R. Bonnerot,  and P. Jamet, ‘Numerical computation of the free boundary for the two-dimensional Stefan problem by space-time finite element,’ J. Comput. Phys., 25,   163, (1977).

   9. D.R. Lynch, ‘Unified approach to simulation on deforming elements

with application to phase change problems,’ J. Comput. Phys. 47,   387, (1982).

   10. P. Bach, and O. Hassager, ‘An algorithm for the use of the Lagrangian

specification in Newtonian fluid mechanics and applications to free  surface flow,’

 J. Fluid Mech. 152,   173, (1985).

   11. M.S. Sadeghipour, and F. Mashayek, ‘A finite element analysis of

the transient freezing of liquids flowing inside tubes,’ Proc. National Heat Transfer Conference, HTD-Vol. 109  , (1989).

   12. M.S. Engelman, and R.L. Sani, ‘Finite element simulation of incompressible

fluid flows with a free/moving surface,’ Recent Advances in Numerical Methods

in Fluids, 5  , 47, (1986).

   13. I.S. Kang,  and L.G. Leal, ``Numerical solution of axisymmetric, unsteady free-boundary problems at finite Reynolds number. 1.Finite-difference scheme and its application to the deformation of a bubble in a uniaxial straining flow,''  Phys. Fluids, 30,   1929, (1987).

   14. G. Ryskin, and L.G. Leal, ‘Numerical solution of the free boundary

problems in fluid mechanics. Part 1. The finite-difference technique,’ J. Fluid Mech., 148,   1,(1984).

   15. N.S. Asaithambi, ‘Computation of free surface flows,’ J. Comput. Phys., 73,   380, (1987).

   16. D.E. Fyfe, E.S. Oran,  and M.J. Fritts,  Naval Res. Lab. NRL Memorandum Report No. 6185, (1988).

   17. M.J. Fritts, and J.P. Boris, ‘The Lagrangian solution of transient

problems in hydrodynamics using a triangular grid,’

 J. Comput. Phys., 31,   173, (1979).

   18. I-L. Chern, J. Glimm, O. McBryan, B. Plohr, and S. Yaniv, ‘Front tracking for gas dynamics,’ J. Comput. Phys., 62,   83-110, (1986).

   19. A.E.P. Veldman, {National Aerospace Laboratory,  ‘Liquid sloshing under low-$g$ conditions: Mathematical model and basic numerical method,’ NLR Report TR 79057 U, The Netherlands,  (1979).

   20. J.E. Welch, F.H. Harlow, J.P. Shannon, and B.J. Daly,

‘The MAC Method: a computing technique for solving viscous, incompressible, transient fluid-flow problems involving free surfaces,’ Los Alamos Scientific Lab. Report No. LA-3425, (1966).

   21. F.H. Harlow, and J.F. Welch, ‘Numerical calculation of time-dependent

viscous incompressible flow of fluid with free surface,’ Phys. Fluids, 8,   2182, (1965).

   22. H. Miyata, and S. Nishimura, ‘Finite-difference simulation of nonlinear ship waves,’ J. Fluid Mech., 157,   327, (1985).

 

   23. H. Miyata, ‘Finite-difference simulation of breaking waves,’ J. Comput. Phys., 65,   179, (1986).

   24. C.W. Hirt, and B.D. Nichols, ‘Volume of Fluid (VOF) method for the

dynamics of free boundaries,’ J. Comput. Phys., 39,   201, (1981).

   25. B.D. Nichols, C.W. Hirt,  and R.S. Hotchkiss, ‘SOLA-VOF: A solution

algorithm for transient fluid flow with multiple free boundaries,’ Los Alamos Scientific Lab. Report No. LA8355, (1980).

   26. M.D. Torrey, L.D. Cloutman, R.C. Mjolsness,  and C.W. Hirt, ‘NASA-VOF2D: A computer program for incompressible flows with free surfaces,’ Los Alamos National Lab. Report No. LA-10612-MS, (1985).

   27. W.F. Noh, underbar {Methods in Computational Physics,   3, 

edited by   B. Adler, S. Fernbach and M. Rotenberg, p. 117, (1964).

   28. W.F. Noh, and P. Woodward, ‘SLIC (Simple Line Interface Calculation),’

in  Proceedings of the Fifth International Conference on Numerical Methods in Fluid Dynamics, edited by A.I. van de Vooren and P.J. Zandbergen  (lecture Notes in

Physics,  59, Springer-Verlag, NY, (1976)), pp. 330-340.

   29. A.J. Chorin, ‘Flame advection and propagation algorithm,’

 J. Comput. Phys., 35, 1, (1980).

   30. A.J. Chorin, ‘Curvature and solidification,’ J. Comput. Phys., 57,   472, (1985).

   31. N. Ashgriz, and J.Y. Poo, ‘FLAIR: Flux Line-Segment Model for Advection

and Interface Reconstruction,’ J. Comput. Phys., 93,   449,(1991).

   32. F. Mashayek and N. Ashgriz, ‘Advection of axisymmetric interfaces,’

Submitted to  Int. J. for Numer. Meth. in Fluids  , (1993).

   33. T.J.R. Hughes, W.K. Liu, and A. Brooks, `Finite element analysis of

incompressible viscous flows by the penalty function formulation' J. Comput. Phys., 30, 1-60 (1979).

   34. W.R. Buell, and B.A. Bush, ‘Mesh generation - A survey,’ ASME J. Industrial Engineering, Ser B, 95,   332, (1973).

   35. R. Haber, M.S. Shephard, J.F. Abel, R.H. Gallagher, and D.P. Greenberg,

‘A general two-dimensional, graphical finite element preprocessor utilizing discrete transfinite mapping,’ Int. J. Num. Meth.  Eng., 17,   1015, (1981).

   36. T.I-P. Shih, R.T. Bailey, H.L. Nguyen, and R.J. Roelke, ‘Algebraic grid generation for complex geometries,’ Int. J.  Num. Meth.  Eng., 13,   1, (1991)

   37. C. Chinnaswamy, B. Amadei, and T.h. Illangasekare,

‘A new method for finite element transitional mesh generation,’

 Int. J. Num. Meth.  Eng., 31, 1253, (1991).

    38. L. A. Crivelli, and S. R. Idelsohn,  ‘A temperature-based finite element solution for phase-change problems,’  Int. J.  Num. Meth. Eng., 23,   99, (1986).

    39. S. O. Unverdi, and G. Tryggvason, ‘A front-tracking method for viscous,

incompressible, multi-fluid flows,’

 J. Comput. Phys., 100,   25, (1992).

   40. J.U. Brackbill, D.B. Kothe, and C. Zemach, ‘A continuum method for

modeling surface tension,’ J. Comput. Phys., 100,   335,( 1992).

   41. W.S. Rayleigh, ‘On the instability of jets’,  Proc. London Math. Soc.   10, 4-13 (1879).

   42. M.C. Yuen, ‘Non-linear capillary instability of a liquid jet,’ J. Fluid Mech., 33  , 151-163, (1968).

   43. S. Chandrasekhar, underbar Hydrodynamic and Hydromagnetic Stability    Claredon, Oxford, (1961).

   44. E.F. Goedde, and M.C. Yuen, ‘Experiments on liquid jet instability,’ J. Fluid Mech., 40, 495-511, (1970).

   45. P. Vassallo, and N. Ashgriz, ‘Satellite formation and merging in

liquid jet breakup,’ Proc. Royal Soc. London A   433,

269-286, (1991).

    46. N.N. Mansour, and T.S. Lundgren, ‘Satellite formation in capillary

jet breakup,’ Phys. Fluids, A2(7),   1141-1144, (1990).

   47. N. Ashgriz and J.Y. Poo, ‘Coalescence and separation in binary

collisions of liquid drops,’ J. Fluid Mech., 221  , 183, (1990).

   48. J.Y. Poo, and N. Ashgriz, ‘Numerical simulation of collision of

two-dimensional drops,’ Proc. of  5th Annual Conference on

Liquid Atomization and Spray Systems, (1992).