Equal order method for fluid flow simulation

Abstract

The present invention provides a method for solving the Navier-Stokes equation of viscous, incompressible laminar flows with moving free surfaces in complex domains. The method uses fixed mesh control volume-finite element techniques to track the flow. A gauss point velocity vector is defined as the average of its nodal counterpart, such that the gauss point velocity vector is constant over the element. The gauss point velocity vector is then inserted into the continuity constraint to form the Poisson pressure equation for solving the pressure field. The solution to the Poisson pressure equation is unique, the common checker-board problem is therefore eliminated. The corrected pressure field is substituted into the momentum equations, so that the resulting velocity field satisfies the continuity equation. Since velocity and pressure are evaluated at the same order, the global mass conservation can be evaluated to machine round-off tolerances.

Description

FIELD OF THE INVENTION

This invention relates to a numerical method for solving the Navier-Strokes equation and has application in computer simulation of viscous, incompressible laminar fluid flow; in particular, it relates to a method using finite element-control volume techniques, wherein the velocity and pressure of the fluid are evaluated at the same order to preserve accuracy of solution without sacrificing computation efficiency.

BACKGROUND OF THE INVENTION

Modeling incompressible, viscous flow generally involves solving the Navier-Stokes equation to obtain the velocity field. The Navier-Stokes equation describes the fundamental law of mass conservation and momentum balance of incompressible viscous flow. For a given pressure field, there is no particular difficulty in solving the momentum equation. By substituting the correct pressure field into the momentum equation, the resulting velocity field satisfies the continuity constraint. The difficulties in calculating the velocity field lies when the pressure field is unknown. The pressure field is indirectly specified via the continuity equation. Yet, there is no obvious equation for obtaining pressure.

The difficulty associated with the determination of pressure has led to methods that eliminate pressure from the momentum equations. For example, in the well-known “stream-function/vorticity method” described by others and Gosman et. al. (In Heat and Mass Transfer in Recirculating Flows , Academic, New York, 1969.) The stream-function/vorticity method has some attractive features, but the pressure, which has been so cleverly eliminated, frequently happens to be a desired result or even an intermediate outcome required for the calculation of density and other fluid properties. In addition, the method cannot be easily extended to three-dimensional situations. Because most practical problems are three-dimensional, a method that is intrinsically restricted to two dimensions suffers from a serious limitation.

Numerical methods have been developed for solving the momentum equation and obtaining the pressure field. The SIMPLE (Semi-Implicit Method for Pressure-Linked Equations by Patenkar, in Numerical Heat Transfer and Fluid Flow , McGraw Hill, 1980), and SIMPLER (SIMPLE Revised) algorithms are well-known examples of the numerical methods. The SIMPLE and SIMPLER methods discretize the flow domain, and transform the momentum equation into a set of algebraic equations. The pressure-gradient term is represented in the discretized momentum equations by the pressure drop between two grid points. One consequence of this representation is the non-uniqueness of the solution of the pressure field. As long as the pressure values at alternate grid points are uniform, the discretized form of the momentum equations is incapable of distinguishing a checker-board pressure field from a uniform pressure field. However, the potential to result in a checker-board pressure field should make the solution unacceptable. Another consequence of this representation is that pressure is, in effect, taken from a coarser grid than the one actually employed. This practice diminishes the accuracy of the solution. It has been recognized that the discretizing equation does not demand that all the variables be calculated on the same grid points. The SIMPLE and SIMPLER methods resolve the checker boarding of pressure by using a “staggered” grid for the velocity component. In the staggered grid, the velocity components are calculated on a grid staggered in relation to the normal control volume around the main grid points. A computer program based on the staggered grid must carry all the indexing and geometric information about the locations of the velocity components and must perform interpolation for values between grid points. Furthermore, a “staggered” grid is difficult to apply to a flow domain that is meshed into irregularly shaped finite elements.

Reduced order and reduced intergration methods have been suggested for modeling fluid flow in complex fluid domains. The reduced order method suggests that pressure be evaluated in a coarser grid, e.g. at alternate nodal positions of the finite element mesh (Baliga and Patankar, A Contol Volume Finite-Element Method for Two Dimensional Fluid Flows and Heat Transfer, in Nutmerical Heat Transfer , 6, 245-261 (1980)). The reduced order method diminishes the accuracy of the pressure field, and therefore compromises the accuracy of the subsequent mass flux calculations. Additionally, the method has been shown to be inadequate for modeling a fluid flow with a moving free surface. The reduced integration method suggests that less than full integration should be used for the pressure term. This method has shown some success in specific applications, but there are many limitations to its general application. Furthermore, reduced integration methods also sacrifice accuracy of the solution.

There is therefore a need for a method to solve the momentum equations of incompressive, viscous fluid flows such that accuracy of solution is preserved without sacrificing computation efficiency.

SUMMARY OF THE INVENTION

Accordingly, one preferred object of the present invention is to provide a numerical method for solving the Navier-Stokes equation of viscous incompressible, laminar fluid flow and wherein the velocity and pressure of the fluid are evaluated at the same order to preserve accuracy of solution without sacrificing computation efficiency.

Another preferred object of the present invention is to provide a numerical method for solving the momentum equation and providing a solution of such accuracy that the global mass conservation can be maintained to machine round-off tolerances.

One embodiment of the present invention is a numerical method for solving the coupled momentum and mass conservartion equation. The method resides and operates within a computer processor. The method defines a gauss point velocity vector as the average of the nodal velocity vectors. In essense, the expression of gauss point velocity vector is in the same order as the pressure gradient over each element. Because the gauss point fluxing vector field is “centered” inside the pressure field defined at the nodes, no checker boarding will occur. Additionally, because the gauss point velocities and pressure fields are evaluated from the same number of nodal positions, the accuracy of the solution is preserved. The predicted gauss point fluxing velocity field conserves mass to machine round-off levels for element types in which the pressure gradient and velocities can be evaluated on a consistent basis. In two dimensions, 3 noded, linear triangular elements will conserve mass to machine round-off tolerances. In three dimensions, 4 noded, linear tetrahedral elements will conserve mass to machine round-off tolerances.

One aspect of the present invention is to provide a numerical method operating within a processor environment for solving the momentum equation in a computer simulation of a viscous, incompressible, laminar flow within a flow domain comprising: (a) discretizing the flow domain into a mesh of elements having nodes and fluxing surfaces; (b) forming momentum equations; (c) computing nodal hat velocities at each node; (d) calculating gauss point hat velocities and gauss point coefficient, wherein the gauss point hat velocities and gauss point coefficient are average value of nodal hat velocities and nodal coefficient; (e) updating velocity field; and (f) checking if velocity field has converged, and returning to the forming momentum equation if the velocity field has not converged.

Other objects, embodiments, forms, aspects, features, advantages, and benefits will become apparent from the description and drawings provided herein.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a schematic diagram of a machine to implement the present invention.

FIG. 2 is a flowchart schematically showing an embodiment of the method of the present invention.

FIG. 3 is an example of a triangular finite element similar to those used in the illustrated embodiment of the method of the present invention.

FIG. 4 is an example of a triangular finite element mesh including a control volume and associated sub-control volumes and fluxing surfaces, preferred by the illustrated embodiment of the method of the present invention.

FIG. 5 is an example of a tetrahedral finite element preferred by an embodiment of the method of the present invention.

DESCRIPTION OF THE PREFERRED EMBODIMENT

For the purposes of promoting an understanding of the principles of the invention, reference will now be made to the embodiment illustrated in the drawings and specific language will be used to describe the same. It will nevertheless be understood that no limitation of the scope of the invention is thereby intended. Any alterations and further modifications in the illustrated device, and any further applications of the principles of the invention as illustrated therein being contemplated as would normally occur to one skilled in the art to which the invention relates are also included.

The illustrated embodiment of the present invention is a numerical method applicable for solving the momentum and mass conservation equations to obtain the velocity field and pressure distribution of incompressible fluid flow in a complex domain. The method can provide accurate solution that the global mass conservation can be maintained to the truncation limit of the computation device (machine round-off levels). An embodiment of the method is designed for incorporation in computer simulations of molten material flows.

The preferred embodiment of the present invention is designed to operate within a processor based computation device. Referring now to FIG. 1 , there is shown an apparatus for carrying out the preferred embodiment of the invention. A conventional computer system 50 , including a processor 60 including ROMs, RAMs, operatively connected by cables to a display 52 , a keyboard 54 , a mouse 56 , a printer 58 , and a storage device 62 , though a variety of connection means and input and output devices may b e substituted without departing from the invention. Processor 60 , in cooperation with an operation system, executes the codes of the present method, receives and stores data from the input devices, performs computations, and transmits data to the output devices. The illustrated embodiment of the method operates within a Silicon Graphic system using UNIX operating system. Other systems having similar computation capability, memory capacity and graphic displays can also be used.

The present invention is developed based on a discretized mesh of the flow domain. Any mesh may be employed, for example, a regular grid created according to the finite difference method or an irregular mesh created according to a finite element method. It is preferred that the grid can be split into uniquely defined control volumes and each of these control volumes can be associated with a unique node or grid point. It is also preferred that the flow domain is discretized into a finite element mesh. The irregular finite elements provide more flexibility in fitting irregular domains and in providing local grid refinement. In addition, the method is applicable to either moving or fixed finite element meshes. A fixed mesh is preferred.

Referring now to FIG. 2 which shows a flow-chart of an embodiment of the method of the present invention.

(a) Creation of the Finite Element Mesh, 101 .

In 101 , the flow domain is discretized or meshed into discrete control volumes to solve the governing equations. A finite element mesh of the flow domain may be created by any conventional procedures. For example, U.S. Pat. No. 5,010,501 issued on Apr. 24, 1991 to Arakawa and titled “Three-Dimensional Geometry Processing Method and Apparatus Therefor”, and U.S. Pat. No. 5,553,206 issued on Sep. 3, 1996 to Meshkat and titled “Method and System for Producing Mesh Representations of Objects”, both disclose suitable procedures.

The illustrated embodiment of the method is based on a mesh composed of a three-node, two dimensional, triangular finite element similar to that shown in FIG. 3 and 4 . Referring now to FIG. 3 , the triangle finite element 10 has three vertex nodes 12 . Thermal physical parameters of the flow, e.g. temperature, velocity and pressures are evaluated at nodal (grid) points 12 of element 10 . Particularly for this invention, gauss point fluxing vectors P g and V g are also computed.

Finite element 10 can be divided into three sub-control volumes (SCVs) 18 by joining the mid-side points 14 to the centroid 16 . Each of sub-control volumes 18 is bounded by four fluxing surfaces 30 where material can flow from one sub-control volume to another. Typically, a number of sub-control volumes 18 surround a vertex node 12 and conveniently defines a control volume (CV) for that node. The sub-control volumes 18 that are attached to a particular node i, are collectively known as the control volume of node i. For example, as shown in FIG. 4 , vertex node 22 is surrounded by six sub-control volumes 23 - 28 and the shaded area represents control volume 20 of vertex node 22 . While it is illustrated that control volume 20 is defined by the six sub-control volumes 23 - 28 attached to vertex node 22 , control volumes can be composed of one or more sub-control volumes depending on the geometry of the flow domain. For a given mesh, a flexible combination of the SCVs is much better than solely relying on CVs in defining the irregular shape of the free surface of the flowing fluid.

For modeling a three-dimensional problem, where a three dimensional finite would be required, a tetrahedral finite element is preferred. As shown in FIG. 5 , tetrahedral element 40 has four vertex nodes 42 and four surfaces. Finite element 40 can be divided into four sub-control volumes (SCVs) 48 by joining the surface centroid 44 (centroid at each surface of the tetrahedral element 40 ) to the volume centroid 46 (the centroid of the tetrahedron). Each of the sub-control volumes 48 is a hexahedron having six fluxing surfaces.

While triangular and tetrahedral finite element are illustrated and preferred, it should be understood that other geometric shaped finite elements or a combination thereof might be used. It is also understood that linear finite elements are preferred but are not required.

(b) Inputing Initial Parameters, 102 .

In 102 , the finite element matrix, the properties of the fluid and initial values of variable parameters, e.g., temperature, pressure, density and velocities, are input. These values will be updated as the iterative computation progresses.

(c) Form Momentum Equations, 103 .

The general equations governing the conservation of mass and momentum of the fluid flow can be expressed by equations (1)-(4)

Conservation of mass: $$ ∂ ρ ∂ τ + ∂ ( ρ ⁢   ⁢ u ) ∂ x + ∂ ( ρ ⁢   ⁢ v ) ∂ y + ∂ ( ρ ⁢   ⁢ w ) ∂ z ⁢   = 0 ( 1 )

Conservation of momentum: $$ ρ ⁢ ∂ u ∂ τ + ρ ⁡ ( V · ∇ u ) = - ∂ P ∂ x + ρ ⁢   ⁢ g x + ∇ ( μ ⁢   ⁢ ∇ u ) ( 2 ) ρ ⁢ ∂ v ∂ τ + ρ ⁡ ( V · ∇ v ) = - ∂ P ∂ y + ρ ⁢   ⁢ g y + ∇ ( μ ⁢ ∇ v ) ( 3 ) ρ ⁢ ∂ w ∂ τ + ρ ⁡ ( V · ∇ w ) = - ∂ P ∂ z + ρ ⁢   ⁢ g z + ∇ ( μ ⁢ ∇ w ) ( 4 )

where:

x, y, z=cartesian coordinate variables

u, v, w=magnitude of the velocity vector V in the x-, y-, and z-directions, respectively,

τ=time

ρ=density

P=pressure

V=flow velocity vector

μ=viscosity

g x , g y , g z =magnitude of gravity in the x-, y-, and z-directions, respectively.

The momentum equations (2)-(4) above describes the fundamental conservation law of incompressible viscous flow. Equations (2)-(4) are made up of several terms each representing a physical phenomenon pertaining to fluid flow. Examples of the terms include inertia, advection, diffusion, body force and pressure gradient. Generally, the terms pertinent to the problem being solved are selected and the remaining terms are grouped and collectively defined as the source term. Using control volume techniques, an elemental matrix can be computed to approximate the terms of interest. These matrices are then assembled into a global momentum matrix (discretization equations) to be solved iteratively.

The Galerkin method and also Schneider and Zedan (“Control-Volume-Based Finite Element Formulation of the Heat Conduction Equations”, Progress in Astronomical and Aeronautical Engineering , Vol. 86, pp. 305-327, 1983) have demonstrated the procedure for derivation of the discretized equation for a finite element mesh. In both derivations, a “shaped function” or profile assumption is used to describe how a dependent variable varies over an element.

In assembling the global matrix of the illustrated embodiment of the present invention, the methodology presented by Schneider and Zedan is followed. However, any other method for forming a global matrix may be used. In the illustrated embodiment, the terms for advection, diffusion, and a pressure gradient are selected. Each term is integrated over a sub-control volume (e.g. 18 of FIG. 3 ) and then assembled to form global equations. If the integrand is considered constant, then equation (5) results as follows: $$ ∫ v ⁢ c ⁢   ⁢ ⅆ V = c ⁢   ⁢ V ( 5 )

For the advection terms, a streamlined upwind operator is used, as represented by the following equation (6): $$ ∫ v ⁢ [ u ⁢ ∂ φ ∂ x + v ⁢ ∂ φ ∂ y + w ⁢ ∂ φ ∂ z ] ⁢   ⁢ ⅆ V = ∫ v ⁢ [ u s ⁢ ∂ φ ∂ s ] ⁢   ⁢ ⅆ V = u s ⁢ V ⁢   ⁢ Δ ⁢   ⁢ φ Δ ⁢   ⁢ s ( 6 )

where φ is any dependent variable, S is distance along a streamline, subscript ν denotes volume and subscript s denotes streamline function.

For the diffusion terms, Green's theorem is used to convert the volume integral to an integration of fluxes through the bounding surface, as modeled by equation (7): $$ ∫ v ⁢ ∇ α ⁢ ∇ φ ⁢   ⁢ ⅆ V = ∫ σ ⁢ α ⁢ ∇ φ ·   ⁢ ⅆ ω = ∫ σ ⁢ J ·   ⁢ ⅆ ω ( 7 )

where σ, ν denote surface and volume integral, respectively, and φ is any dependent variable and J is the flux contribution. The flux contribution for a face j from a node i, is expressed by equation (8) as follows: $$ J i ⁢   ⁢ j = α ⁡ [ ∂ N i ∂ r ⁢ ( a 11 ⁢ ω x ⁢   ⁢ j + a 21 ⁢ ω y ⁢   ⁢ j + a 31 ⁢ ω z ⁢   ⁢ j ) + ∂ N i ∂ t ⁢   ⁢ ( a 13 ⁢ ω x ⁢   ⁢ j + a 23 ⁢ ω y ⁢   ⁢ j + a 33 ⁢ ω z ⁢   ⁢ j ) + ∂ N i ∂ s ⁢ ( a 12 ⁢ ω x ⁢   ⁢ j + a 22 ⁢ ω y ⁢   ⁢ j + a 32 ⁢ ω z ⁢   ⁢ j ) ] ( 8 )

where:

N i is the shape function, a 11 -a 33 are the cofactors of the Jacobian matrix, and ω xj , ω yj , ω zj are the components of the area ω, and α is diffusivity.

After the assembly process, each term at the left-hand side of the momentum equations is represented by a coefficient matrix operating on a velocity component. The resulting assembled global momentum equations can be modeled by equations (9) and (10) as follows $$ ∑ j = 1 N ⁢   ⁢ a i ⁢   ⁢ j ⁢ u j = - ∫ A ⁢ ∂ P ∂ x ⁢   ⁢ ⅆ A ( 9 ) ∑ j = 1 N ⁢   ⁢ a i ⁢   ⁢ j ⁢ v j = - ∫ A ⁢ ∂ P ∂ y ⁢   ⁢ ⅆ A ( 10 )

where the a ij terms make up the influence coefficient matrix. Pulling the pressure gradients outside the area integrals, equations (9) and (10) is recast as equations (11)-(15) as follows: $$ u i = u ^ i - K i ⁢ ∂ P ∂ x ( 11 )

and $$ v i = v ^ i - K i ⁢ ∂ P ∂ y ( 12 )

where, $$ K i = 1 a i ⁢   ⁢ i ⁢ ∫ ⅆ A , ( 13 ) u ^ i = - 1 a i ⁢   ⁢ i ⁢ ∑ i ≠ j ⁢   ⁢ a i ⁢   ⁢ j ⁢ u j , ( 14 )

and $$ v ^ i = - 1 a i ⁢   ⁢ i ⁢ ∑ i ≠ j ⁢   ⁢ a i ⁢   ⁢ j ⁢ v j ( 15 )

and nodal hat velocities, û i and {circumflex over (v)} i are composed of neighbor nodal velocities u j and contain no pressure component. It should be understood that while the illustrated embodiment of the present invention is applied to a two-dimensional flow, the method is applicable to one- and three-dimensional flows.

(d) Compute nodal hat velocities û i , {circumflex over (v)} i , 104 .

In 104 , the coefficient matrix a ij is first calculated from the estimated velocity field provided initially. Using these newly calculated values of the coefficient matrix and the nodal velocity values, u j and v j , from the neighboring sub-control volumes, the nodal hat velocities, û i and {circumflex over (v)} i , are computed.

(e) Compute the Gauss Point Hat Velocities and Gauss Point Coefficient, 105

In 105 , the gauss point hat velocity û g , {circumflex over (v)} g and the gauss point coefficient term K g are defined as the average value of their nodal counterparts, as provided by the following equations (16)-(18): $$ u ^ g = 1 N e ⁢ ∑ j = 1 N e ⁢   ⁢ u ^ j , ( 16 ) v ^ g = 1 N e ⁢ ∑ j = 1 N e ⁢   ⁢ v ^ j , ( 17 )

and $$ K g = 1 N e ⁢ ∑ j = 1 N e ⁢   ⁢ K j ( 18 )

and the gauss point velocities are expressed as: $$ u g = u ^ g - K g ⁢ ∂ P ∂ x ( 19 ) v g = v ^ g - K g ⁢ ∂ P ∂ y ( 20 )

In this expression for the gauss point fluxing velocities u g and v g , the pressure gradients ∂P/∂x, ∂P/∂y are constant over the fluxing surface. For consistency, the terms u g , v g and K g have to be constant over the fluxing surface as well. Traditionally, the interpolation of the gauss point terms would have been computed using the element shape functions, Ni. However, computations involving the shape functions would have given quantities like the gauss point velocities a linear variation across the element, while the pressure gradient is constant over the element. This inconsistency will result in computational schemes that can not conserve mass to machine round-off levels.

(f) Form and Solve the Poisson Pressure Equation, 106 .

In 106 , by inserting the expression for the gauss point fluxing velocities u g and v g into the continuity equation (21), $$ ∂ u ∂ x + ∂ v ∂ y = 0 ( 21 )

a Poisson equation for pressure is obtained. The Poisson pressure equation (22) is expressed as: $$ ▽ ⁢   ⁢ k ⁢   ⁢ ▽ ⁢   ⁢ p = ∂ u ^ ∂ x + ∂ v ^ ∂ y ( 22 )

Integrating this equation over an area A, the fluxing surface of the element, the pressure field can be expressed as: $$ ∫ A ⁢ ▽ ⁢   ⁢ k ⁢   ⁢ ▽ ⁢   ⁢ P ⁢ ⅆ A = ∫ A ⁢ [ ∂ u ^ ∂ x + ∂ v ^ ∂ y ] ⁢ ⅆ A ( 23 )

Using Green's theorem, these area integrals are converted to line integrals, as represented by the following equation (24): $$ ∫ σ ⁢ [ k g ⁢ ∂ P ∂ x ] ⁢ ⅆ y + [ k g ⁢ ∂ P ∂ x ] ⁢ ⅆ x = ∫ σ ⁢ u ^ g ⁢ ⅆ y - v ^ g ⁢ ⅆ x ( 24 )

where σ is the line path. Using the gauss point hat velocities and gauss point coefficient, û g , {circumflex over (v)} g , and k g previously calculated in 105 , the Poisson pressure equation can be assembled and solved to determine the pressure field. While it is preferred that the area integrals is converted into line integrals before the pressure field expression (23) is solved, other methods which solve the original area integrals, e.g. using a Galerkin formulation, may be used.

(g) Update Velocity Field, 107 .

The nodal velocities, u i and v i and the gauss point velocities, u g and v g , are computed by inserting the newly derived pressure field and û i , {circumflex over (v)} i , K i , û g , {circumflex over (v)} g , and K g in the momentum equations.

The discretization equations for other variables, e.g. temperature, are solved using the update velocity field if they influence the flow field. If a particular variable does not influence the flow field, it is preferable to calculate it after a converged solution for the flow field has been obtained.

(h) Check for the Velocity Converging, 108 .

In 108 , an iterative process is said to have converged when further iterations will not improve the accuracy of the dependent variables. In practice, the iterative process is terminated when the desired accuracy is obtained. Although the computed velocity field at the gauss points conserves mass to machine (computation device) round-off limits at every iteration, the momentum equations are non-linear and therefore iterations are required to refine the momentum balance. In this technique, the flow equations are iterative solved until all velocity values converge to a user supplied tolerance. Numerical tests have revealed that the above algorithm, converges rapidly for highly viscous flows. When lower Reynolds numbers are encountered, solving the momentum equations as a complete equation set can aid in obtaining a more rapid convergence. While the above convergence criterion was chosen for the illustrated embodiment, other convergence criteria may be used without deviating from the spirit of the present invention.

If the velocity fields have not converged, then return to 103 of forming the momentum equations using the updated velocity field. If the velocity field has converged, the method is completed.

As a result of making the conservation of mass and momentum equations consistent at the element level, the predicted gauss point fluxing velocity field does indeed conserve mass to machine round-off levels.

Furthermore, pressure checker-boarding is eliminated, because the gauss point fluxing vector field is centered between the pressures defined at the nodes. Because these vectors are constant within an element, the right number of degrees of freedom is ensured, guaranteeing a unique solution to pressure and velocity.

It should be understood that the present invention of using gauss point velocity and pressure fields to improve accuracy and computation efficiency in solving the coupled mass and momentum equation is not limited to the methodology illustrated; it is applicable to other discretized forms of the momentum equations derived by other formulation, e.g. Galerkin's method or Newton-Raphson's method. In addition, the present invention is also applicable to methods which solve the equation implicitly, e.g. the penalty function method.

While the invention has been illustrated and described in detail in the drawings and foregoing description, the same is to be considered as illustrative and not restrictive in character, it being understood that only the preferred embodiment has been shown and described and that all changes and modifications that come within the spirit of the invention are desired to be protected.

Claims

1. A numerical method operating within a processor environment for solving the momentum equation in a computer simulation of a viscous, incompressible, laminar flow within a flow domain comprising:(a) discretizing said flow domain into a mesh of elements having nodes and fluxing surfaces; (b) forming momentum equations; (c) computing nodal hat velocities at each said node; (d) calculating gauss point hat velocities and gauss point coefficient, wherein said gauss point hat velocities and gauss point coefficient are average value of nodal hat velocities and nodal coefficient; (e) updating velocity field; (f) checking if velocity field has converged, and returning to said forming momentum equation if said velocity field has not converged.

2. The method of claim 1, wherein said forming momentum equation includes inputting matrix of said mesh, and initial values of variable parametes.

3. The method of claim 2, wherein said forming of momentum equations includes selecting terms for said momentum equation, integrating said terms over a sub-control volume, and assembling said terms to form global momentum equations represented as follows: ∑Nj=1⁢aij⁢uj=-∫A⁢∂P∂x⁢ⅆA∑Nj=1⁢aij⁢vj=-∫A⁢∂P∂y⁢ⅆAand wherein said global momentum equations each including pressure gradient within an area integral.

4. The method of claim 3, wherein said forming momentum equations further includes recasting said global momentum equations, after pulling said pressure gradients outside said area integral, as follows: ui=u^i-Ki⁢∂P∂x⁢ ⁢and⁢ ⁢vi=v^i-Ki⁢∂P∂ywhere, Ki=1aii⁢∫ⅆA,⁢u^i=-1aii⁢∑i≠j⁢aij⁢uj,and v^i=-1aii⁢∑i≠j⁢aij⁢vj.

5. The method of claim 4, wherein said calculating gauss point hat velocities and gauss point coefficient further includes expressing gauss point velocities, ug, vg, as: ug=u^g-Kg⁢∂P∂x⁢ ⁢and⁢ ⁢vg=v^g-Kg⁢∂P∂ywherein said gauss point hat velocities ûg and {circumflex over (v)}g, and gauss point coefficient Kg are element average quantities, such that: u^g=1Ne⁢∑Nej=1⁢u^j,v^g=1Ne⁢∑Nej=1⁢v^j,Kg=1Ne⁢∑Nej=1⁢Kj.

6. The method of claim 5, wherein said pressure gradients, ∂P/∂x and ∂P/∂y in said gauss point velocity expression are constant over an fluxing surface of an element, and said gauss point velocities, ug, vg, and gauss point coefficient Kg are constant over said fluxing surface of said element.

7. The method of claim 6, wherein said updating velocity field further includes forming a Poisson pressure equation, solving said Poisson pressure equation for a gauss point pressure field.

8. The method of claim 7, wherein said forming a Poisson pressure equation further includes inserting the expression for said gauss point fluxing velocities, ug, vg, into a continuity equation and expressing said Poisson pressure equation as: ▽⁢ ⁢k⁢ ⁢▽⁢ ⁢p=∂u^∂x+∂v^∂y.⁢

9. The method of claim 8, wherein said solving said Poisson pressure equation for a gauss point pressure further includes integrating said Poisson equation over an area of said fluxing surfaces of said elements, and forming said pressure field represented as: ∫A⁢▽⁢ ⁢k⁢ ⁢▽⁢ ⁢P⁢ⅆA=∫A⁢[∂u^∂x+∂v^∂y]⁢ⅆA.

10. The method of claim 9, wherein said solving said Poisson pressure equation further includes converting area integrals in said pressure field to line integrals via Green's theorem, and represented as: ∫σ⁢[kg⁢∂P∂x]⁢ⅆy+[kg⁢∂P∂x]⁢ⅆx=∫σ⁢u^g⁢ⅆy-v^g⁢ⅆx.

11. A method, comprising:providing a mesh of a flow domain of a fluid flow, the mesh including a number of finite elements with corresponding nodes; determining a number of momentum equilibrium equations, the equations each being determined as a function of pressure and velocity at a corresponding one of the nodes; defining velocity at a Gauss-point inside each of the elements to provide a Gauss-point velocity field, and pressure gradient over each of the elements, having the same order as the Gauss-point velocity field; calculating velocity and pressure at each of the nodes as a function of the Gauss-point velocity field; and repeating said determining, said defining, and said calculating until a convergence test is satisfied.

12. The method of claim 11 wherein said determining comprises inserting velocities at said Gauss-point inside each of the elements into a continuity equation forming a Poisson pressure equation, and solving said Poisson pressure equation for a pressure field.

13. The method of claim 12, further comprising determining a velocity field as a function of the pressure field.

14. The method of claim 11, wherein the finite elements are each of a three-node linear triangle form, and the velocity at the Gauss-point inside each of the elements is an average of velocities defined at three corresponding nodes.

15. A method for simulating a fluid flow within a flow domain, comprising:(a) defining the flow domain with a mesh of finite elements each having a corresponding number of nodes; (b) determining momentum relationships based on velocity and pressure defined at each of the nodes of the finite elements; (c) calculating a Gauss-point velocity inside each of the finite elements, wherein the Gauss-point velocity is computed from the velocity at each of the nodes of a corresponding one of the finite elements; (d) updating the velocity at each of the nodes as a function of said calculating; (e) repeating said determining, calculating, and updating until a convergence test is satisfied.

16. The method of claim 15, wherein said mesh is a fixed mesh.

17. The method of claim 15, wherein said finite elements are triangular elements.

18. The method of claim 15, wherein said finite elements are tetrahedral elements.

19. The method of claim 15, which includes determining a pressure field as a function of the Gauss-point velocity inside each of the finite elements.

20. The method of claim 15, which includes operating a computer to perform said defining, determining, calculating, updating, and repeating.