Patent Publication Number: US-7899654-B2

Title: Hybrid front tracking algorithm for solving single phase fluid equations with a moving boundary on a quadrilateral grid

Description:
CROSS REFERENCE TO RELATED APPLICATIONS 
     This Application is related to application Ser. No. 11/777,076 entitled “A hybrid method for enforcing curvature related boundary conditions in solving one-phase fluid flow over a deformable domain,” filed on Jul. 12, 2007. This application is herein incorporated by reference in its entirety for all purposes. 
     BACKGROUND 
     Inkjet technology is traditionally used in a general home/office printing uses. In order to increase the resolution and to explore further uses of the inkjet technology, it is important to know how the fluid droplet from the inkjet head will finally settle when they dry. 
     The present invention is directed towards the study of single-phase fluids. A single-phase fluid, as used in the present invention refers to a pure fluid or a mixture. A boundary encapsulates this single-phase fluid, and over time the boundary may evolve and move through space. The solution or fluid domain is restricted to an area delimited by the boundary. For the purposes of the present invention, the effect of an interface, if any that does form between elements of such a mixture may be ignored when considering the motion of the boundary. 
     Motion of the boundary may be governed by one or more equations. The boundary may be defined by a set of markers. For some systems the equations may include terms related to the curvature of the boundary. Prior art methods have calculated the curvature of the boundary from a cubic spline fitted to the markers that define the boundary. Small errors in the positions of the marker can cause errors in the second derivatives, which can create variations in the curvature. Under certain conditions, this variation may increase in magnitude as the boundary evolves in time. This can lead to numerically unstable simulations, which exhibit unrealistic oscillations. 
     Prior art methods have attempted to address the problems that arise from defining the boundary using markers by employing multi-phase fluid simulation methods (e.g., front tracking methods). 
     Problems may arise when handling the boundary between the two fluids. Each fluid is characterized by a set of system variables (i.e., density, viscosity, and pressure). A jump condition necessarily exists at the boundary. This jump condition is represented by a step wise change in the system variables and results in a discontinuity in system variables and their derivatives at the boundary. This discontinuity can prevent numerical solutions from converging. Further, to accommodate the complex geometry, a quadrilateral grid is needed to be set aligned with the geometry of the bank structure. 
     SUMMARY 
     An embodiment of the present invention may be a system or method for simulating the flow of single-phase fluid flow. Markers represent a moving fluid boundary of the single-phase fluid at a first point in time. The moving fluid boundary separates a simulation space into a fluid space and a vacuum space. The single-phase fluid inhabits the fluid space. A signed distance function is evaluated at points surrounding the moving fluid boundary based upon markers. The curvature of the boundary is evaluated near the markers in the vacuum space based on the signed distance function. The velocity of the fluid and pressure jump conditions are calculated based upon the curvature of the boundary. The position of the markers is updated at a second point in time based on the velocity of the fluid. 
     In one embodiment, a computer implemented method for simulating a final pattern of a droplet of a fluid having a plurality of fluid properties is disclosed. The method includes representing a position of a moving fluid boundary at a first point in time with a plurality of fluid markers, wherein the moving fluid boundary provides a boundary for a single-phase fluid. The moving fluid boundary separates a simulation space into a fluid space and a vacuum space, and the single-phase fluid inhabits the fluid space. The method further includes evaluating a plurality of indicator functions defined on a quadrilateral grid for a velocity field, a pressure, and the plurality of fluid markers. The quadrilateral grid is formed by a plurality of points inside the fluid space and the vacuum space. The method also includes calculating surface tension at the moving fluid boundary through a level set method, solving a plurality of flow equations through a projection method, and representing the position of the moving fluid boundary at a second point in time by updating the positions of the fluid markers based on the plurality of velocity values. The fluid properties and the shape of the final pattern is then stored on a computer readable media. 
     In another embodiment, a system for simulating a final pattern of a droplet of a fluid is disclosed. The system includes a computer readable media to provide storage for programming instructions. The programming instructions including programming instructions to implement an indicator function evaluation module, a fluid boundary monitor module, a flow equation solver module, and a fluid properties display and storage module. A processor coupled to the computer readable media to execute the programming instructions is also included. A display is provided to plot the final pattern of the droplet and a plurality of fluid properties of the fluid. The fluid boundary monitor module marks a moving fluid boundary of the droplet. The indicator function evaluation module solves a plurality of indicator functions defined on a quadrilateral grid for a velocity field, a pressure, and the plurality of fluid markers. The flow equation solver module solves a plurality of flow equations through a projection method for a surface location on a fluid boundary of the droplet. 
     In yet another embodiment, a computer readable media having program instructions for simulating a final pattern of a droplet of a fluid having a plurality of fluid properties is provided. The computer readable media includes program instructions for representing a position of a moving fluid boundary at a first point in time with a plurality of fluid markers. The moving fluid boundary provides a boundary for a single-phase fluid. The moving fluid boundary separates a simulation space into a fluid space and a vacuum space, and the single-phase fluid inhabits the fluid space. Program instructions for evaluating a plurality of indicator functions defined on a quadrilateral grid for a velocity field, a pressure, and the plurality of fluid markers, are also provided. The quadrilateral grid is formed by a plurality of points inside the fluid space and the vacuum space. Further, program instructions for calculating surface tension at the moving fluid boundary through a level set method are provided. The computer readable media further includes program instructions for solving a plurality of flow equations through a projection method, program instructions for representing the position of the moving fluid boundary at a second point in time by updating the positions of the fluid markers based on the plurality of velocity values, and program instructions for storing the plurality of fluid properties and a shape of the final pattern on a computer readable media. 
     The advantages of the embodiments of the present invention are numerous. Most notably, since the embodiment of the present invention uses a quadrilateral grid for computation of various properties such as pressure and velocity, the simulation is realistic. In real life applications, the droplet is located in a more complex bank structure. To accommodate the complex geometry, a more complicated quadrilateral grid is needed to be set aligned with the geometry of the bank structure. 
     Other aspects and advantages of the present invention will become apparent from the following detailed description, taken in conjunction with the accompanying drawings, illustrating by way of example the principles of the present invention. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
       The present invention can be best understood by reference to the following description taken in conjunction with the accompanying figures, in which like parts may be referred to by like numerals. 
         FIG. 1  illustrates an exemplary use of the inkjet printing technology in the fabrication of Liquid Crystal Displays, in accordance with one embodiment of the present invention. 
         FIG. 2  illustrates a droplet on a non-uniform surface, in accordance with one embodiment of the present invention. 
         FIG. 3A  illustrates an exemplary quadrilateral grid layout, in accordance with one embodiment of the present invention. 
         FIG. 3B  and  FIG. 3C  illustrate exemplary quadrilateral grids and portions of the droplet to evaluate velocities extrapolations, in accordance with one embodiment of the present invention. 
         FIG. 4A  illustrates an exemplary quadrilateral grid of computational cells in which the pressure is evaluated at the center of each computational cell, in accordance with one embodiment of the present invention. 
         FIG. 4B  illustrates an exemplary quadrilateral grid for calculating curvatures using a level set method, in accordance with one embodiment of the present invention. 
         FIG. 5  illustrates an exemplary quadrilateral grid, the fluid markers, the fixed boundary, the deformable boundary, and the virtual boundary, in accordance with one embodiment of the present invention. 
         FIG. 6  illustrates a control flow of the process of simulating a final pattern of a droplet, in accordance with one embodiment of the present invention. 
         FIG. 7  illustrates an evolution of a droplet by surface tension, in accordance with one embodiment of the present invention. 
         FIG. 8  illustrates a system for simulating a final pattern of a droplet, in accordance with one embodiment of the present invention. 
     
    
    
     The figures are provided in order to provide a thorough understanding of the present invention. The figures should not be construed as limiting the breath of the invention in any manner. 
     DETAILED DESCRIPTION 
     The present invention relates to systems and methods for simulating a single-phase fluid flow enclosed by a deformable boundary while enforcing curvature related boundary and pressure jump conditions. The solution or fluid domain is restricted to an area delimited by the boundary. Further, the present invention is directed towards creating numerically stable simulations. 
     Inkjet technology is conventionally used printing applications such as printing on papers, etc. However, tremendous amount of resources can be saved by applying the inkjet technology in industrial printing process. For example, the inkjet technology can be used in semiconductor fabrication, fabricating Liquid Crystal Display (LCD), Thin film transistors (TFT), Organic light emitting diode (OLED) devices, solar cells, etc. 
     In the inkjet printing process, small droplets containing the solute of desired properties are deposited to the target area. After the droplets reaching the targeted area, the solvent evaporates and only the solute is left forming a final print pattern. The final pattern of the solute deposit directly determines the desired product quality. 
     In order to achieve a better quality of the final product, it is important to understand how the final pattern is formed in realistic environment. Embodiments of the present invention provide simulating the final pattern for a given solute in a given solvent under certain conditions such as pressure, the structure on which the droplets would rest, etc. 
     It can be challenging work to solve single-phase fluid flow problems in which the boundary conditions include time dependent deformable domains (e.g., a moving surface). Generally, only a few specific cases with a time dependent domain can be solved analytically. Most of the solutions to such problems can only be obtained through numerical methods. 
     Since the boundary of the domain is evolving with time, a time dependent solution needs to be calculated. This adds a further level of complexity to the numerical solution. The boundary conditions of such a problem may involve curvature related terms. How to treat these boundary conditions correctly in the numerical algorithm directly affects the accuracy and numerical stability of the numerical solution. With this overview in mind, the following figures will illustrate examples, structures, and functionality of methods and systems for simulating a final pattern of a droplet. 
       FIG. 1  illustrates an exemplary application of the use of inkjet printing technology in industrial printing. In one embodiment, the inkjet technology can be used for fabrication of liquid crystal display (LCD) devices  50 . Traditionally, LCDs  50  are fabricated using CMOS fabrication technology which requires several processing steps such as masking, unmasking, etc., to fabricated Red (R), Green (G), and Blue (B) pixels  52  on a non-uniform surface. With the use of the inkjet technology, the RGB pixels can be printed on the surface in one step. As well known, the inkjet technology uses small droplets of a fluid for printing. The inkjet printing head deposits small droplets in a desired pattern on the surface. Until the drops dry, shapes of these drops change continuously due to fluid evaporation, surface tension, fluid properties of the fluid, and properties of the surface, etc. Therefore, the final pattern of the drops may change into a shape and size that may not be desirable. For example, in the LCD fabrication example, if droplets for R and G are printed too far away from each other, a lower LCD resolution will result. However, if these droplets are places too close, the drops may overlap when dried. Therefore, it is important to simulate the final pattern of a drop for a better design of the devices. In one embodiment, the simulation of the final shape of the droplet can be repeated several times by adjusting the fluid properties of the droplet until a desired final shape is obtained. It should be noted that  FIG. 1  simply illustrates one of many uses of the inkjet technology in industrial applications. Many other applications are possible if a time dependent shape of the droplet on a complex bank structure can be simulated. 
       FIG. 2  illustrates a droplet of a single phase fluid on a complex bank structure surface  108 . A liquid  101  forms a drop on a substrate  108  in a vacuum  106 . The simulation boundary consists of a fixed boundary  104  and a deformable boundary  102 . The fixed boundary  104  separates the liquid  101  from the substrate  108  while deformable boundary  102  separates the liquid  101  from the vacuum  106 . In the present invention, the deformable boundary  102  of the domain is represented by a set of markers. 
     In order to achieve a better quality of the final product, it is important to understand how the final pattern of the droplet is formed in a realistic environment and the major factors affecting the final pattern of the droplet. Further, it is also important to understand as to how to control the fluid properties and other parameters that affect the final pattern to achieve the desired pattern. Although simulations in a rectangular domain provide some insights into achieving a desired final pattern, in a realistic setup, the droplet is located in a more complex bank structure  104 . To accommodate the complex geometry, a quadrilateral grid is needed to be set aligned with the geometry of the bank structure  104 . 
     In one embodiment, a hybrid method is used for enforcing the curvature related boundary conditions for solving a single phase fluid flow problem on a quadrilateral grid. In this method, the computational grid is assumed to be quadrilateral for both finite difference and finite element methods. The number of nodes that needs to be solved vary due to the change in computational domain over time. The boundary of the domain is explicitly tracked by a number of marker points on the fluid surface. The marker points move with respect to the background quadrilateral grid. Further, the pressure jump condition is enforced along this moving boundary. 
     Governing Equations 
     The liquid drop of an incompressible fluid sits on a substrate surface with two contact points. In one embodiment the liquid drop has a low Reynolds number, i.e., the ratio of inertial forces to viscous forces. In this embodiment, the creeping flow equation (1A) is solved with incompressibility condition (1B) enforced. The creeping flow equation (1A) and the incompressibility condition (1B) may be considered the governing equations for this problem. 
                     ρ   ⁢       ∂   u       ∂   t         =       -     ∇   p       +     μ   ⁢       ∇   2     ⁢   u                 (     1   ⁢   A     )                 ∇     ·   u       =   0           (     1   ⁢   B     )               
In the equation (1A), ρ is the density, u is the velocity, t is the time, p is the pressure, and μ is the dynamic viscosity. In one embodiment, it is assumed that the drop meets with outside vacuum and that the pressure in vacuum is zero. The boundary condition along the liquid-vacuum surface is:
 
p=σκ  (2)
 
In equation (2), p is the pressure inside the fluid, σ is the surface tension on the surface of the drop, and κ is the curvature of the surface. In this embodiment, a projection method is utilized to solve equations (1A), (1B), and equation (2).
 
     Since a quadrilateral grid is being used, equation (1A) is transformed from physical space X=(x, y) to transformed space  =(ξ,η). The governing equation (1A and 1B) changes to: 
                       ρ   ⁢       ∂   u       ∂   t         =         -     1   J       ⁢     T   T     ⁢       ∇   Ξ     ⁢   p       +       μ   J     ⁢       ∇   Ξ     ⁢     ·     [       1   J     ⁢     TT   T     ⁢       ∇   Ξ     ⁢   u       ]               ,           (     3   ⁢   A     )                   ∇   Ξ     ⁢     ·     (   Tu   )         =   0.           (     3   ⁢   B     )               
In equation (3A), T is the transformation matrix and J is the Jacobian. The Jacobian depends of the cell being referred to, in the quadrilateral grid. J and T are solved by the following equation (4):
 
                     J   =       det     ⁢     (           x   ξ           x   η               y   ξ           y   η           )         ,     T   =       (           y   η           -     x   η                 -     y   ξ             x   ξ           )     .               (   4   )               
Equation (4) is used to solve the governing equations in a transformed space.
 
     Numerical Algorithm 
     In one embodiment, a numerical procedure for solving single-phase incompressible fluid  101  problems with a deformable boundary  102  (i.e. time dependent deformable solution domain) may use markers to track the boundary. The numerical procedure may include deriving indicator functions defined on a quadrilateral grid to solve for a velocity field, a pressure field, and the marker positions. The numerical procedure may also include deforming the boundary  102  if there are velocity boundary conditions on the boundary or for other reasons. The numerical procedure further includes calculating the surface tension term at the boundary  102  using a level set approach. It may be noted that the method of calculating surface tension term using the level set method is disclosed in the U.S. patent application Ser. No. 11/777,076, which is incorporated herein by reference, hence, the details of such calculations are being omitted. 
     Further, the numerical procedure also includes using a projection method to solve a set of flow equations given a particular location of the deformable boundary  102 . Because the fluid properties are unknown outside the fluid domain, to solve the flow equations, the velocity field at a time step is extrapolated to grid locations near but outside the surface or boundary  102 .  FIG. 3A  illustrates a cell  131  of the quadrilateral grid. The velocities and level set values are defined at cell centers  132 . Pressure and curvature values are defined at the grid points  133 . Referring now to  FIG. 3B  and  FIG. 3C , if a point (x,y) i,j  at which the velocity is unknown and needs to be determined is outside the computational domain (i.e., outside the boundary  102 ), the velocity at point (x,y) i,j  can be determined by applying one of the following rules. If one of the points that is adjacent to point (x,y) i,j  in the quadrilateral grid  130  is inside the fluid domain, then the velocity at point (x,y) i,j  is set to equal the velocity at the adjacent point in the fluid domain. Hence, as illustrated in  FIG. 3B , the velocity at point (x,y) i,j  will be set equal to the velocity at point (x,y) i,j−1  because point (x,y) i,j−1  is adjacent to point (x,y) i,j  and also point (x,y) i,j−1  is inside the fluid domain. 
     Referring now to  FIG. 3B . If three points that are adjacent to point (x,y) i,j  are in the fluid domain (i.e., inside the fluid, hence, the velocities at these points are known), the velocity at point (x,y) i,j  (i.e., point four in the calculations below) can be calculated by extrapolation. Assuming the coordinates of the points in three dimension are (x i , y i , z i ), where i=1, 2, 3, and 4 (representing four points in question). Since passing through any three points, a plane can be determined using equation (5): 
                                x   -     x   1             y   -     y   1             z   -     z   1                 x   -     x   2             y   -     y   2             z   -     z   2                 x   -     x   3             y   -     y   3             z   -     z   3                  =   0           (   5   )               
Since point four is also on this plane, the values at point four can be calculated through extrapolation using equation (6):
 
                                  x   4     -     x   1               y   4     -     y   1               z   4     -     z   1                   x   4     -     x   2               y   4     -     y   2               z   4     -     z   2                   x   4     -     x   3               y   4     -     y   3               z   4     -     z   3                  =   0.           (   6   )               
In equation (6), z 4  can be calculated using equation (7). z 4  could be the horizontal or vertical velocity.
 
     
       
         
           
             
               
                 
                   
                     
                       z 
                       4 
                     
                     = 
                     
                       
                         
                           
                             z 
                             1 
                           
                           ⁢ 
                           A 
                         
                         - 
                         
                           
                             z 
                             2 
                           
                           ⁢ 
                           B 
                         
                         + 
                         
                           
                             z 
                             3 
                           
                           ⁢ 
                           C 
                         
                       
                       
                         A 
                         - 
                         B 
                         + 
                         C 
                       
                     
                   
                   , 
                   
                     
 
                   
                   ⁢ 
                   where 
                   , 
                   
                     
 
                   
                   ⁢ 
                   
                     A 
                     = 
                     
                        
                       
                         
                           
                             
                               
                                 x 
                                 4 
                               
                               - 
                               
                                 x 
                                 2 
                               
                             
                           
                           
                             
                               
                                 y 
                                 4 
                               
                               - 
                               
                                 y 
                                 2 
                               
                             
                           
                         
                         
                           
                             
                               
                                 x 
                                 4 
                               
                               - 
                               
                                 x 
                                 3 
                               
                             
                           
                           
                             
                               
                                 y 
                                 4 
                               
                               - 
                               
                                 y 
                                 3 
                               
                             
                           
                         
                       
                        
                     
                   
                   , 
                   
                     
 
                   
                   ⁢ 
                   
                     B 
                     = 
                     
                        
                       
                         
                           
                             
                               
                                 x 
                                 4 
                               
                               - 
                               
                                 x 
                                 1 
                               
                             
                           
                           
                             
                               
                                 y 
                                 4 
                               
                               - 
                               
                                 y 
                                 1 
                               
                             
                           
                         
                         
                           
                             
                               
                                 x 
                                 4 
                               
                               - 
                               
                                 x 
                                 3 
                               
                             
                           
                           
                             
                               
                                 y 
                                 4 
                               
                               - 
                               
                                 y 
                                 3 
                               
                             
                           
                         
                       
                        
                     
                   
                   , 
                   
                     
 
                   
                   ⁢ 
                   
                     C 
                     = 
                     
                       
                          
                         
                           
                             
                               
                                 
                                   x 
                                   4 
                                 
                                 - 
                                 
                                   x 
                                   1 
                                 
                               
                             
                             
                               
                                 
                                   y 
                                   4 
                                 
                                 - 
                                 
                                   y 
                                   1 
                                 
                               
                             
                           
                           
                             
                               
                                 
                                   x 
                                   4 
                                 
                                 - 
                                 
                                   x 
                                   2 
                                 
                               
                             
                             
                               
                                 
                                   y 
                                   4 
                                 
                                 - 
                                 
                                   y 
                                   2 
                                 
                               
                             
                           
                         
                          
                       
                       . 
                     
                   
                 
               
               
                 
                   ( 
                   7 
                   ) 
                 
               
             
           
         
       
     
     After extrapolating the velocity field, momentum equations are solved to obtain an intermediate velocity field. The intermediate velocity field is extrapolated, using the method of extrapolation described above, to grid locations near but outside the deformable boundary. 
     Solving of the flow equations also requires solving the pressure Poisson equation derived from the incompressibility condition of the fluid, and simultaneously enforcing the pressure jump boundary condition along the free surface using a finite element approach. The Poisson equation (8): 
                         ∇   2     ⁢   p     =       ρ     Δ   ⁢           ⁢   t       ⁢     ∇     ·     u   *             ⁢     
     ⁢       u   *     =       (     1   +         Δ   ⁢           ⁢   t     ρ     ⁢   μ   ⁢           ⁢     ∇   2         )     ⁢   u               (   8   )               
can be solved for the pressure of the fluid  101  while enforcing an incompressibility condition of the fluid  101  and a boundary condition that is at least partially in motion. In equation (8) p is the pressure, ρ is the density, Δt is the time step, and {right arrow over (u)}* is the intermediate velocity vector. Also in equation (8) {right arrow over (u)} is the velocity vector and μ is the viscosity of the fluid  101 .
 
     As the boundary  102  deforms so do the marker positions. The markers may bunch up or spread out depending on the non-uniformity of the motion associated with the boundary  102 . New markers may be generated to keep the distance between markers within a given range (e.g. 0.2-0.8 times the cell size) or depending on the curvature of the boundary  102 . The new marker positions may be used in the next time step. 
     The boundary  102 , which is in motion, is represented explicitly by markers, which sit on the surface of the boundary  102 . These markers are not tied to a specific point in space, instead they move with the boundary  102  as the boundary  102  is moved by the local velocity of the fluid  101 , so do the markers. 
     The position of the boundary  102  is tracked explicitly in the time domain thus placing the present invention in the category of “front tracking” fluid flow simulation methods. Which should be distinguished from “front capturing” methods, such as the level set approach, in which the boundary is implied by the zero value contour of a level set function. Prior art methods have used front tracking methods to solve multi-phase fluid flow problems. In the past, front tracking methods have not been used to solve single-phase fluid flow problems. Front tracking methods were not used because a methodology did not exist to effectively enforce the boundary conditions, whereas, front capturing methods have problems conserving the mass of the fluid. The present invention conserves the mass of the fluid better than front capturing methods. The present embodiments address this drawback of the front tracking method when applied to single-phase problems. 
     Along the boundary  102 , curvature related boundary conditions often play an important role in determining local velocities of the fluid  101 . For example, when attempting to calculate the flow of a fluid  101 , the stress the fluid  101  experiences at the boundary  102  must meet certain conditions. An example of such a condition may include balancing the normal stress using the following equation (9).
 
[ p−n   i τ ij   n   j ]=σκ  (9)
 
In equation (9) the quantities in the square brackets [ ] represent a pressure jump condition across the moving boundary  102 . The pressure jump is defined as a difference between a pressure on a point immediately inside the boundary and immediately outside the boundary. The difference in pressure is caused by the surface tension on the boundary. Further, due to difference in curvature, the pressure immediately inside the boundary  102  is not the same at different points inside the boundary  102 . The symbol p is representative of the pressure, the symbol τ ij  is representative of the stress tensor, the symbols n i  and n j  are representative of the unit normal vector at the surface, σ is the surface tension coefficient of the fluid  101 , and κ is the local surface curvature. Enforcing curvature related boundary conditions can cause instabilities to occur in a numerical method because the curvature is related to the second derivative of the positions of the markers. In the prior art the curvature had been calculated based on the second derivative of a line fitted to the markers.
 
     In order to simplify the following discussion of the hybrid method, the stress term in equation (9) is ignored in the following discussion. Hence, the boundary condition becomes equation (10).
 
[p]=σκ  (10)
 
An individual skilled in the art would appreciate that the stress term from equation (9) may be included in the present invention without going beyond the spirit and scope of the present invention as recited in the claims.
 
     When the Poisson equation is solved for the pressure p, the pressure p may be evaluated at the grid points (i.e., corners of each cell) of the computational cell as shown in  FIG. 4A . The discrete pressure value at a point when the point is located outside the fluid domain enclosed by the surface or boundary  102 , is determined by the boundary condition, equation (10). 
     As illustrated in  FIG. 4A  in two space dimensions, in the discrete Poisson equation, p i+1,j+1  is set to p 0 +σκ (this value comes from the boundary condition equation (10) with the stress neglected), where p 0  is the pressure outside the fluid domain (or boundary). Next, the pressure is extrapolated to get a pressure at a point outside the surface  102  by solving the Poisson equation (3) derived by incompressibility using a finite element approach, equation (11). 
                       ∫   Ω     ⁢         Δ   ⁢           ⁢   t     ρ     ⁢       ∇   p     ·     ∇   Ψ       ⁢     ⅆ   x         =       ∫   Ω     ⁢         u   *     ·     ∇   Ψ       ⁢     ⅆ   x                 (   11   )               
Then, to get the pressure at a marker through extrapolation, a discretized Poisson equation (12) is solved for cells of the computational domain.
 
                     ∫     -   1     1     ⁢       ∫     -   1     1     ⁢           Δ   ⁢           ⁢   t       J   ⁢           ⁢   ρ       [       ∑     l   =   1     4     ⁢           p   ^     l     ⁡     (           N     l   ,   ξ             N     l   ,   η             )       ⁢     (           y   η           -     x   η                 -     y   ξ             x   ξ           )         ]     ·               [       ∑     k   =   1     4     ⁢           Ψ   ^     k     ⁡     (           N     k   ,   ξ             N     k   ,   η             )       ⁢     (           y   η           -     x   η                 -     y   ξ             x   ξ           )         ]     ⁢     ⅆ   ξ     ⁢     ⅆ   η       =       ∫     -   1     1     ⁢       ∫     -   1     1     ⁢       (           u   *           v   *           )     ·       ∑     k   =   1     4     ⁢           Ψ   ^     k     ⁡     (           N     k   ,   ξ             N     k   ,   η             )       ⁢     (           y   η           -     x   η                 -     y   ξ             x   ξ           )     ⁢     ⅆ   ξ     ⁢     ⅆ   η                               (   12   )               
where u* and v* are the intermediate fluid velocities that were calculated previously, ρ is the fluid density, and Δt the time step size. If more than one neighboring locations are “out of phase”, for each out-of-phase location, {circumflex over (p)} l , we enforce {circumflex over (p)} l =σκ and move the corresponding terms to the right hand side of the linear equations.
 
     As mentioned earlier, the pressure at various points inside the surface  102  is different due to the surface curvature. The local curvature is a measure of the deviation of the fluid boundary  102  from a straight line. Given a boundary  102 , there is a normal vector {right arrow over (n)} on the boundary  102  that is perpendicular to the boundary. To measure the deviation of the boundary one can take the differential of the normal vector {right arrow over (∇)}·{right arrow over (n)} along the boundary  102 . To calculate the curvature κ of the surface  102 , first, the marker positions are utilized to generate a signed distance function, φ, over the whole domain or over a narrow band around the surface  102 . Since the marker positions at a particular time step are known, the value of signed distance function, φ, can be evaluated at any point inside the computational domain and in a narrow band around the surface  102 . 
     As indicated earlier, the marker positions are used to represent the boundary  102 . These marker positions are also used to generate a signed distance function φ. The signed distance function φ may be evaluated over the whole domain or over a narrow band around the boundary  102 . The signed distance function may be virtually extended into the region of space on both sides of the fixed boundary  104 . Two markers which represent the boundary  102  closest to the fixed boundary  104  may be used to linearly extend the boundary  102  and create a virtual boundary  104 . Instead of using the marker positions to generate the normal vector, the signed distance function may be used as in equation (13): 
     
       
         
           
             
               
                 
                   
                     n 
                     → 
                   
                   = 
                   
                     
                       
                         ∇ 
                         → 
                       
                       ⁢ 
                       ϕ 
                     
                     
                        
                       
                         
                           ∇ 
                           → 
                         
                         ⁢ 
                         ϕ 
                       
                        
                     
                   
                 
               
               
                 
                   ( 
                   13 
                   ) 
                 
               
             
           
         
       
     
     Since the value of φ can be evaluated at the corner points of each cell, the curvature value for each cell is used to approximate the average curvature along the surface  102  in the cell. This is a zero th  order approximation. In this way, even the markers on the surface moving close to each other as the surface evolves, the curvature term remains relatively smooth. 
     Next, φ is used to calculate the curvature value at the cell corners. This calculation is performed by employing the centered difference scheme. Since the curvature, 
                     κ   =       ∇     ·       ∇   ϕ            ∇   ϕ                =       1   J     ⁢       ∇   Ξ     ⁢     ·     (     T   ⁢         T   T     ⁢       ∇   Ξ     ⁢   ϕ                T   T     ⁢       ∇   Ξ     ⁢   ϕ                )               ,           (   14   )               
is defined at the same location as the grid point, the curvature at point (i,j) can be calculated using equation (15).
 
                         κ     i   ,   j       =         [       ∇   Ξ     ⁢     ·     (     T   ⁢         T   T     ⁢       ∇   Ξ     ⁢   ϕ                T   T     ⁢       ∇   Ξ     ⁢   ϕ                )         ]       i   ,   j       =       K       i   +     1   2       ,   j       -     K       i   -     1   2       ,   j       +     K     i   ,     j   +     1   2           -     K     i   ,     j   -     1   2                 ,     
     ⁢   where   ,     
     ⁢       K       i   +     1   2       ,   j       =         (     T   ⁢         T   T     ⁢       ∇   Ξ     ⁢   ϕ                T   T     ⁢       ∇   Ξ     ⁢   ϕ                )         i   +     1   2       ,   j       =       T       i   +     1   2       ,   j       ⁢           T       i   +     1   2       ,   j     T     ⁡     (       ∇   Ξ     ⁢   ϕ     )           i   +     1   2       ,   j                T       i   +     1   2       ,   j     T     ⁢       ∇   Ξ     ⁢   ϕ                      ⁢     
     ⁢       T       i   +     1   2       ,   j       =       1   2     ⁢     (       T     i   ,   j       +     T     i   ,     j   -   1           )         ⁢     
     ⁢         (       ∇   Ξ     ⁢   ϕ     )         i   +     1   2       ,   j       =       (             1   2     ⁢     (       ϕ       i   +   1     ,   j       -     ϕ       i   -   1     ,   j       +     ϕ       i   +   1     ,     j   -   1         -     ϕ       i   -   1     ,     j   -   1           )                   ϕ     i   ,   j       -     ϕ     i   ,     j   -   1                 )     .               (   15   )               
Similarly the forms for
 
               K       i   -     1   2       ,   j       ,     K     i   ,     j   +     1   2           ,     and   ⁢           ⁢     K     i   ,     j   -     1   2                   
can be derived.
 
An individual skilled in the art will appreciate that this may be expanded to three dimensions or other coordinate systems without going beyond the scope of the present invention as recited in the claims.
 
     Equation (15) can be used to calculate a curvature value κ i,j  at the center of a computational cell. This may be done by approximating the derivatives of the signed distance function using the centered difference method. Other difference methods may be used without going beyond the scope of the invention. This curvature value κ i,j  is an approximation of the averaged curvature of the fluid boundary in the cell or near the cell. 
       FIG. 4B  is an illustration of the discretized signed distance function φ i,j  along with the boundary  102 . As shown in  FIG. 4B  the discretized signed distance function φ i,j  is evaluated at the cell centers. In one embodiment, if (x,y) i,j , (x,y) i+1,j , (x,y) i,j+1 , (x,y) i+1,j+1  are the coordinates of the four corners, the cell center&#39;s coordinate is (x,y) c  can be calculated using the following equation: 
     
       
         
           
             
               
                 ( 
                 
                   x 
                   , 
                   y 
                 
                 ) 
               
               c 
             
             = 
             
               
                 1 
                 4 
               
               ⁢ 
               
                 ( 
                 
                   
                     
                       ( 
                       
                         x 
                         , 
                         y 
                       
                       ) 
                     
                     
                       i 
                       , 
                       j 
                     
                   
                   + 
                   
                     
                       ( 
                       
                         x 
                         , 
                         y 
                       
                       ) 
                     
                     
                       
                         i 
                         + 
                         1 
                       
                       , 
                       j 
                     
                   
                   + 
                   
                     
                       ( 
                       
                         x 
                         , 
                         y 
                       
                       ) 
                     
                     
                       i 
                       , 
                       
                         j 
                         + 
                         1 
                       
                     
                   
                   + 
                   
                     
                       ( 
                       
                         x 
                         , 
                         y 
                       
                       ) 
                     
                     
                       
                         i 
                         + 
                         1 
                       
                       , 
                       
                         j 
                         + 
                         1 
                       
                     
                   
                 
                 ) 
               
             
           
         
       
     
       FIG. 5  is an illustration of the boundary  102  represented by a series of markers.  FIG. 5  also show virtual markers below the fixed boundary  104 . The virtual markers are extended below the fixed boundary  104  by linearly extending the boundary  102  from the fluid side of boundary  104  to the non-fluid side of the boundary. Virtual level set values are calculated relative to the virtual markers. 
       FIG. 6  is an overview of a method for implementing one embodiment of the present invention. This method may be implemented in hardware, software or a combination of both. The method is used to simulate the single-phase fluid  101  as illustrated in  FIG. 2 . 
     The single-phase fluid forms at least one boundary with its environment. This boundary is movable in one embodiment. The single-phase fluid may also form a second stationary boundary. In operation  180 , the moving fluid boundary is marked with a set of markers at a first point in time. Each marker consists of a set of coordinates, which describe the position of a point on the boundary at a particular time. The coordinates may be in one, two, three, or more dimensions. The coordinate system may be Cartesian, polar, cylindrical or some other coordinate system. 
     Next, in operation  182  an indicator function defined on a quadrilateral grid for velocity field, pressure, and the marker points, is evaluated. The boundary is an N−1 dimensional surface that moves through an N dimensional space. The boundary divides the N dimensional space into two groups (a positive space and a negative space). The negative space may be the single-phase fluid or the vacuum or vice versa. The indicator function or the signed distance function may be a smooth continuous function, which is defined throughout the N dimensional space. The signed distance function has a positive value in the positive space and a negative value in the negative space. The absolute value of the signed distance function at a particular point in the N dimensional space is equal to the shortest distance from the particular point to the moving boundary  102 . This statement is also true for all points in the N dimensional space. 
     The discrete number of points that represent the signed distance function may be distributed on a first grid. The first grid may coincide with portions of a second grid upon which values representative of the pressure of the liquid are calculated. The first grid may cover only a portion of the N dimensional space. The first grid may be limited to grid cells, which are one, two, three, or four grid cells from the deformable boundary based upon the order of the approximations used throughout the method  178 . 
     Next, in operation  184 , the surface tension at the deformable or moving boundary is calculated. The surface tension is calculated using a level set method, which is a well known method (hence, details are being omitted). Next, in operation  186 , flow equations are solved for current surface location. Solving the flow equations includes (a) extrapolating the velocity field at time step n to grid locations near but outside the boundary (i.e., immediately outside the boundary). The solving of the flow equations further includes projecting the momentum equations to obtain an intermediate velocity field and extrapolating the intermediate velocity field to grid locations near but outside the boundary. The pressure Poisson equation derived from the incompressibility condition of the fluid is then solved. A pressure jump boundary condition along the free surface is simultaneously enforced using a finite element approach. The curvature of the surface is also calculated in the process of calculating the pressure jump. The solving of the flow equations further includes calculating new velocity, through process of interpolation, at the next time step n+1 using the pressure obtained by deriving the Poisson equation. This new velocity field is extrapolated at time step n+1 to grid locations near but outside the surface. 
     Next, in operation  188  the fluid boundary is marked at a second point in time. The free surface is moved with the velocity field at time step n+1. Then, in operation  190 , a determination whether a final pattern of the droplet has evolved is made. If the final pattern has been evolved, the fluid properties and final pattern of the droplet is stored on a computer readable media. Else, the process is repeated by moving the control to operation  182 . 
     Moving now to  FIG. 7 , the movement of the fluid boundary for an exemplary set of values of surface tension and viscosity of the liquid is illustrated.  FIG. 7  illustrates the deformable boundary at three points in time. An initially deformed drop is driven by surface tension. Initially the drop assumes an elliptic shape with a larger height than width. The surface tension, in this example, is set to 1 (i.e., σ=1). The viscosity of the liquid is also set to 1 (i.e., μ=1). The density, in this example, is set to 1 (i.e., ρ=1). Further, in this example, full slip boundary condition is enforced at the bottom boundary. As apparent from  FIG. 7 , the elliptic drop spreads with time, and gradually moves toward the steady state which is a half circle. 
     System 
     Having described the details of the embodiments of the present invention, an exemplary system  230 , which may be used to implement one or more aspects of the present invention, will now be described with reference to  FIG. 8 . As illustrated in  FIG. 8 , the system includes a processor  242  that provides computing resources and controls the computer. The processor (or CPU)  242  may be implemented with a microprocessor or the like, and may also include a graphics processor and/or a floating point coprocessor for mathematical computations. The system  230  also includes memory  244  which may be in the form of random-access memory (RAM) and read-only memory (ROM). The system  230  may also include an I/O interface  246  to enable data input/output operations as well as interfacing with a fluid preparation system  252 . The I/O interface  246  also provides interface to one or more sensors to input data such as pressure, temperature, viscosity, etc. The fluid preparation system  252  may receive simulation data from the system  230  to enable the fluid preparation system  252  to adjust fluid properties of the fluid that is to be used in the inkjet printing. 
     The system  230  further includes an indicator function evaluation module  232  to provide computational functionality for solving an indicator function defined on a quadrilateral grid for a velocity field, a pressure, and the marker points. A fluid boundary monitor module is also included to monitor and mark the deformable fluid boundary at different points of time. Further, a flow equation solver module  236  is provided to solve flow equations as described in the description of  FIG. 6  and in detailed discussion in the previous paragraphs. A fluid properties display and storage module  238  is provided to display the fluid properties to a user or to store various fluid properties and simulation data in a computer readable media. In one embodiment, all modules that are parts of the system  230  are coupled together through an internal bus  240 . 
     A number of controllers and peripheral devices may also be provided, as shown in  FIG. 8 . The system  230  may also include an external storage media  248  for storing and/or transferring simulation data. A display  250  is also provided to enable display of simulation data and patterns including the final pattern of the droplet, as well as fluid properties of the fluid. The computer readable media such as magnetic tape or disk, or an optical medium can be used to record programs of instructions for operating systems, utilities and applications which may include embodiments of programs that implement various aspects of the present invention. 
     In the illustrated system, all major system components can connect to a bus  240 , which may represent more than one physical bus. However, various system components may or may not be in physical proximity to one another. For example, input data and/or output data may be remotely transmitted from one physical location to another. In addition, programs that implement various aspects of this invention may be accessed from a remote location (e.g., a server) over a network. Such data and/or programs may be conveyed through any of a variety of machine-readable medium including magnetic tape or disk or optical disc, or a transmitter, receiver pair. 
     The present invention may be conveniently implemented with software. However, alternative implementations are certainly possible, including a hardware implementation or a software/hardware implementation. Any hardware-implemented functions may be realized using ASIC(s), digital signal processing circuitry, or the like. Accordingly, any “means” terms in the claims are intended to cover both software and hardware implementations. Similarly, the term “machine-readable medium” as used herein includes software, hardware having a program of instructions hardwired thereon, or a combination thereof. With these implementation alternatives in mind, it is to be understood that the figures and accompanying description provide the functional information one skilled in the art would require to write program code (i.e., software) or to fabricate circuits (i.e., hardware) to perform the processing required. 
     In accordance with further aspects of the invention, any of the above-described methods or steps thereof may be embodied in a program of instructions (e.g., software), which may be stored on, or conveyed to, a computer or other processor-controlled device for execution. Alternatively, any of the methods or steps thereof may be implemented using functionally equivalent hardware (e.g., application specific integrated circuit (ASIC), digital signal processing circuitry, etc.) or a combination of software and hardware. 
     It should be understood that the invention may employ various computer-implemented operations involving data stored in computer systems. These operations are those requiring physical manipulation of physical quantities. Usually, though not necessarily, these quantities take the form of electrical or magnetic signals capable of being stored, transferred, combined, compared, and otherwise manipulated. Further, the manipulations performed are often referred to in terms, such as producing, identifying, determining, or comparing. 
     Any of the operations described herein that form part of the invention are useful machine operations. The invention also relates to a device or an apparatus for performing these operations. The apparatus may be specially constructed for the required purposes, such as the carrier network discussed above, or it may be a general purpose computer selectively activated or configured by a computer program stored in the computer. In particular, various general purpose machines may be used with computer programs written in accordance with the teachings herein, or it may be more convenient to construct a more specialized apparatus to perform the required operations. 
     The invention may be practiced with other computer system configurations including hand-held devices, microprocessor systems, microprocessor-based or programmable consumer electronics, minicomputers, mainframe computers and the like. The invention may also be practiced in distributing computing environments where tasks are performed by remote processing devices that are linked through a network. 
     The invention can also be embodied as computer readable code on a computer readable medium. The computer readable medium is any data storage device that can store data, which can thereafter be read by a computer system. Examples of the computer readable medium include hard drives, network attached storage (NAS), read-only memory, random-access memory, CD-ROMs, CD-Rs, CD-RWs, DVDs, Flash, magnetic tapes, and other optical and non-optical data storage devices. The computer readable medium can also be distributed over a network coupled computer systems so that the computer readable code is stored and executed in a distributed fashion. 
     While the invention has been described in conjunction with several specific embodiments, it is evident to those skilled in the art that many further alternatives, modifications and variations will be apparent in light of the foregoing description. Thus, the invention described herein is intended to embrace all such alternatives, modifications, applications and variations as may fall within the spirit and scope of the appended claims.