Abstract:
A system and method for determining fluid flow of compressible and non-compressible liquids is described. The system may include input means for receiving a model of an object defined as a plurality of cells having a plurality of nodes P, and a processor coupled to a memory. The processor may be configured for: discretizing a partial differential equation (PDE) corresponding to the received model; for each node P: (i) locating all neighbouring cells that share the node P, each of the neighbouring cells having a cell center; (ii) approximating the PDE at the cell center of the neighbouring cells using the discretized PDE; and (iii) updating a solution of the PDE at the node P by using the approximated discretized PDE at all the neighbouring cell centers; and iteratively updating the solution for all the nodes P from an initial guess until a convergence criterion is satisfied.

Description:
RELATED APPLICATIONS 
       [0001]    This application claims the benefit under 35 USC §119(e) to U.S. Provisional Application No. 61/457,589, filed on Apr. 26, 2011, the contents of which are hereby incorporated by reference in their entirety. 
     
    
     FIELD OF THE INVENTION 
       [0002]    This invention relates to a system and method for modelling internal and/or boundary conditions, as for example to model or determine fluid flows in, around and/or across objects or structures and, in particular, for both compressible and non-compressible liquids. 
       BACKGROUND OF THE INVENTION 
       [0003]    Computer methods and algorithms can be used to analyze and solve complex systems involving various forms of fluid dynamics having inputted boundary conditions. For example, computer modelling may allow a user to simulate the flow of air and other gases over an object or model the flow of fluid through a pipe. Computational fluid dynamics (CFD) is often used with high-speed computers to simulate the interaction of one or more fluids over a surface of an object defined by certain boundary conditions. Typical methods involve large systems of equations and complex computer modelling and include traditional finite difference methodology, cell-centered finite volume methodology and vertex-centered finite volume methodology. 
       Traditional Finite Difference Methodology 
       [0004]    Traditional Finite Difference Methodology (TFDM) requires a structured grid system, a rectangular domain and uniformed grid spacing. TFDM cannot be applied on a mesh system with triangular cells (elements). Rather, cells must be quadrilateral (2D) and cannot be polygonal (i.e., number of sides=4). In 3D, cells must be rectangular cubes. 
         [0005]    TFDM typically requires the use of coordinate transformations (i.e., grid generation) for curvilinear domains, to map the physical domain to a suitable computational domain. In addition, there may be a need to use a multiblock scheme if the physical domain is too complicated. Partial differential equations (PDEs) must be transformed to the computational domain. 
         [0006]    Traditional Finite Difference Methodology is typically difficult to deal with in complicated grid arrangements. Special treatment may be required near boundaries of the domain (e.g., in staggered grid systems or for higher-order schemes). Even with coordinate transformations, highly irregular domains may create serious difficulties for accuracy and convergence due to numerical discontinuities in the transformation metrics. Cell-Centered Finite Volume Methodology/Vertex-Centered Finite Volume Methodology 
       Cell-Centered Finite Volume Methodology (CCFVM) and Vertex-Centered Finite 
       [0007]    Volume Methodology (VCFVM) achieve greater flexibility in grid arrangement. Cells can be polygonal (e.g., triangular) in 2 Dimensional space or polyhedral (e.g., tetrahedral, prismatic) in 3 Dimensional space. With CCFVM/VCFVM there is no need for coordinate transformations to a computational domain. Rather, all calculations can be done in physical space. As well, grid smoothness is not an issue. Cell-centered schemes evaluate the dependent variable at the centroid of each cell. Vertex-centered (or vertex-based) schemes evaluate the dependent variable at the vertices of each cell. 
         [0008]    With CCFVM/VCFVM, inaccuracies due to calculation of fluxes across cell faces may be difficult to deal with. In addition, there are difficulties associated with treatment near boundaries for higher-order schemes, and accuracy and convergence issues associated with cells that are severely skewed or have a high aspect ratio. 
         [0009]    Accordingly, current computer modelling schemes are limited in the form of objects they can model and require different models and algorithms for different fluid applications, such as between compressible and non-compressible fluids. 
       SUMMARY OF THE INVENTION 
       [0010]    It is an object of this invention to provide a better method and system for determining and/or modelling boundary conditions, as for example, to determine or compute fluid dynamics of compressible and non-compressible liquids in, around or across objects. In one particular embodiment, it is an object of this invention to provide a better method and system to compute the fluid dynamics of compressible liquids in aeronautical application, the aeronautical applications having certain boundary conditions. 
         [0011]    Furthermore, another object of this invention to provide a better method and system for computing the fluid dynamics of non-compressible liquids within a pipe or transport mechanism, the pipe or transport mechanism having certain boundary conditions. 
         [0012]    The inventors have appreciated that if the solution domain can be discretized into a smooth structured grid, the FDM is better than the FVM (or the FEM) due to its efficiency. In particular, an FDM method requires less memory and has better stability. Furthermore, a system and method relying on an FDM has better convergence properties. 
         [0013]    In one aspect, the present invention resides in a system for modelling internal and/or boundary conditions and more preferably, by determining fluid flow of compressible and non-compressible liquids, as for example, in, around or across an object or structure. 
         [0014]    In one cell-centered finite different approach described hereafter lies in the discovery and development of a unified scheme for the numerical solution of Partial Differential Equations (PDEs), irrespective of their physical origin, which is solely based on the finite difference method, but is implemented in an innovative fashion that allows the use of an arbitrary mesh topology. Thus, the CCFDM enjoys the simplicity and strength of the traditional FDM, and the power and flexibility of the FVM and FEM. 
         [0015]    The resulting program may evolve into entirely new multiphysics computational continuum mechanics software, or replace the core numerical processing component of some existing software packages with a much simpler and more efficient algorithm which permits a natural and seamless coupling of fluid and solid interaction, allows for a more precise analysis of accuracy and produces faster, more accurate and more reliable results. 
         [0016]    The system may include input means for receiving a model of an object defined as a plurality of cells having a plurality of nodes P and a processor coupled to a memory. The processor may be configured for implementing the steps of discretizing a partial differential equation corresponding to the received model of the object; for each node P in the plurality of nodes P: (i) locating all neighbouring cells that share the node P, each of the neighbouring cells having a cell center; (ii) approximating the partial differential equation at the cell center of each of the neighbouring cells using the discretized partial differential equation; and (iii) updating a solution of the partial differential equation at the node P by using the approximated discretized partial differential equation at all the neighbouring cell centers; and iteratively updating the solution for all the nodes P from an initial guess until a convergence criterion is satisfied. 
         [0017]    In another aspect, the present invention resides in a computer-implemented method for approximating a partial differential equation for determining fluid flow of compressible and non-compressible liquids. The method comprising: discretizing the partial differential equation; receiving a model of the object defined as a plurality of cells having a plurality of nodes P; for each node P in the plurality of nodes P: (i) locating all neighbouring cells that share the node P, each of the neighbouring cells having a cell center; (ii) approximating the partial differential equation at the cell center of each of the neighbouring cells using the discretized partial differential equation; and (iii) updating a solution of the partial differential equation at the node P by using the approximated discretized partial differential equation at all the neighbouring cell centers; and iteratively updating the solution for all the nodes P from an initial guess until a convergence criterion is satisfied. 
         [0018]    In yet another aspect, the present invention resides in a computer readable medium having instructions stored thereon that when executed by a computer implement a method for approximating a partial differential equation for determining fluid flow of compressible and non-compressible liquids. The method may include discretizing the partial differential equation; receiving a model of an object defined as a plurality of cells having a plurality of nodes P; for each node P in the plurality of nodes P: (i) locating all neighbouring cells that share the node P, each of the neighbouring cells having a cell center; (ii) approximating the partial differential equation at the cell center of each of the neighbouring cells using the discretized partial differential equation; and (iii) updating a solution of the partial differential equation at the node P by using the approximated discretized partial differential equation at all the neighbouring cell centers; and iteratively updating the solution for all the nodes P from an initial guess until a convergence criterion is satisfied. 
         [0019]    Further and other features of the invention will be apparent to those skilled in the art from the following detailed description of the embodiments thereof. 
     
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
         [0020]    Reference may now be had to the following detailed description taken together with the accompanying drawings, in which: 
           [0021]      FIG. 1  shows schematically a computer system used in the modelling of flow dynamics in accordance with an embodiment of the present invention; 
           [0022]      FIG. 2A  shows a comparison of different modeling approaches of an embodiment of the present invention in comparison to prior art methods; 
           [0023]      FIG. 2B  shows a comparison of different modeling approaches of an embodiment of the present invention in comparison to prior art methods; 
           [0024]      FIG. 2C  shows a comparison of different modeling approaches of an embodiment of the present invention in comparison to prior art methods; 
           [0025]      FIG. 2D  shows a comparison of different modeling approaches of an embodiment of the present invention in comparison to prior art methods; 
           [0026]      FIG. 2E  shows a comparison of different modeling approaches of an embodiment of the present invention in comparison to prior art methods; 
           [0027]      FIG. 3A to 3C  show examples of grid arrangements on which CCFDM, CCFVM and VCFVM can be applied, but TFDM cannot be applied; 
           [0028]      FIGS. 4 and 4A  to  4 D show schematically a generic node P and exemplary two dimensional models, in accordance with an embodiment of the present invention; 
           [0029]      FIG. 5  shows a transformation of a cell center from a physical space (x, y) to a computational space (ξ, η) in accordance with an embodiment of the present invention; 
           [0030]      FIG. 6A  shows a flowchart for calculating flow dynamics in a modeled system in accordance with an embodiment of the present invention; 
           [0031]      FIG. 6B  shows a flowchart for calculating flow dynamics in a modeled system in accordance with another embodiment of the present invention; 
           [0032]      FIGS. 7A to 7D  show different representational models for a unit square in accordance with an embodiment of the present invention; 
           [0033]      FIG. 8A  and  FIG. 8B  show a comparison between a solution for the Laplace equation on a unit square with Dirichlet boundary conditions using an embodiment of the present invention with the exact solution; 
           [0034]      FIGS. 9A to 9C  show a comparison of CCFDM, FEM and the exact solution along horizontal lines y=0.05, 0.5 and 0.95, in accordance with an embodiment of the present invention; 
           [0035]      FIG. 10A  and  FIG. 10B  show different representational models for a unit disk in accordance with an embodiment of the present invention 
           [0036]      FIGS. 11A to 11C  show a comparison of a Laplace equation on a square with Dirichlet and Neumann boundary conditions using an embodiment of the present invention with a TFDM solution; 
           [0037]      FIGS. 12A to 12E  show a comparison of a Laplace equation on a hexagonal ring with Dirichlet and Neumann boundary conditions, in accordance with an embodiment of the present invention; and 
           [0038]      FIG. 13  shows a generic three dimensional cell for a three dimensional model, in accordance with an embodiment of the present invention. 
       
    
    
     DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS 
       [0039]      FIG. 1  illustrates schematically a computer system  10  on which a preferred embodiment of the present invention may be implemented. The computer system  10  includes a system bus  12  for communicating information, and a processor  16  coupled to the bus  12  for processing information. 
         [0040]    The computer system  10  further comprises a random access memory (RAM) or other dynamic storage device  25  (referred to herein as main memory), coupled to the bus  12  for storing information and instructions to be executed by processor  16 . Main memory  25  may also be used for storing temporary variables or other intermediate information during execution of instructions by the processor  16 . The computer system  10  may also include a read only memory (ROM) and/or other static storage device  26  coupled to the bus  12  for storing static information and instructions used by the processor  16 . 
         [0041]    A data storage device  27  such as a magnetic disk or optical disc and its corresponding drive may also be coupled to the computer system  10  for storing information and instructions. The computer system  10  can also be coupled to a second I/O bus  18  via an I/O interface  14 . A plurality of I/O devices may be coupled to the I/O bus  18 , including a display device  24 , an input device (e.g., an alphanumeric input device  23  and/or a cursor control device  22 ), and the like. A communication device  21  is used for accessing other computers (servers or clients) via an external data network (not shown). The communication device  21  may comprise a modem, a network interface card, or other well-known interface devices, such as those used for coupling to Ethernet, token ring, or other types of networks. 
       Embodiment #1 
     Compressible Fluids 
       [0042]    In a preferred embodiment, computer system  10  is used in conjunction with to model the fluid dynamics of an aircraft landing gear through the air, such as during aircraft landing and in flight. In use of the system  10 , a computer-aided designed (CAD) representation of the aircraft landing gear and supporting structure is inserted into the model. 
         [0043]    Airflow, as a compressible fluid, may be constrained by initial conditions entered as part of the model or taken from sensors from real-world applications. The airflow may be modeled as a partial differential equation, as known in the art of fluid dynamics. For example, data from temperature and speed sensors may be included in the model, taken from real-world applications. 
         [0044]    Once the boundary conditions and initial conditions have been inputted, the profile is input into the system  10  of the present invention and when the solution converges to a steady state, the solution is outputted. The solution may describe the flow of compressible fluid for the specific boundary conditions and initial conditions inputted into the model. 
         [0045]    The system advantageously allows a user to determine and analyze the turbulence in the compressible fluid caused by the different aircraft components passing through the airflow. The steady state output can be used to identify and analyze different flow regimes, such as laminar flow and turbulent flow including eddies, vortices and other flow instabilities. In addition, the behaviour of the fluid about the boundary layer is also outputted. In particular, the noise of the flow over the aircraft component can be modeled including the frequency of any noise created. 
         [0046]    It should be understood that the system  10  is capable of modeling any type of compressible fluid through a wide variety of applications, as further discussed below. Besides modelling the air passing over an aircraft component, other applications may include engine design, wind-tunnel effects and other airflow applications. In addition, the compressible fluid may be in a confined space, such as within a tunnel, or in a non-confined space, such as in flight. 
       Embodiment #2 
     Non-Compressible Fluids 
       [0047]    The above-described computer system  10  can also be used to model the fluid dynamics of a non-compressible fluid through a defined space. For example, in a preferred embodiment, the computer system  10  can model a fluid such as water through a pipe or other transport mechanism. 
         [0048]    As with the compressible embodiment, described above, a computer-aided designed (CAD) representation of the pipe is inserted into the simulation. Typical boundary conditions may be represented in the model. 
         [0049]    The system  10  then models the flow of the non-compressible fluid, i.e. water or gas, through the pipe in successive stages. The non-compressible fluid may be further defined by its initial conditions or parameters. For example, the non-compressible fluid may include particulate matter and have a specific viscosity. The non-compressible fluid may be constrained by initial conditions entered as part of the computer simulation or taken from sensors from real-world applications. These parameters may be inserted into the partial differential equation (PDE) used to model the compressible fluid flow. For example, flow and temperature data from real-world flow-analysis may be inputted automatically into the simulation. 
         [0050]    Once the solution of the system has converged to a steady state, the solution is outputted. The output data transformed into a usable format for describing the flow of the non-compressible fluid for the specific boundary conditions and initial conditions inputted into the simulation. 
         [0051]    The solution to the simulation advantageously allows a user to determine and analyze the turbulence in the non-compressible fluid caused by the boundary conditions (i.e. the pipe). The steady state output can be used to identify and analyze different flow regimes, such as laminar flow and turbulent flow including eddies, vortices and other flow instabilities. In addition, the behaviour of the fluid about the boundary layer is also outputted. Furthermore, the simulation may model the aggregate (i.e. the particulate matter) in the fluid and the Reynolds Number (Re), as would be known to person skilled in the art. 
         [0052]    It should be understood that the system  10  is capable of modeling any type of non-compressible fluid through a wide variety of applications. Besides simulating the flow of fluid passing through a pipe, other applications may include oil and gas applications and hydraulics. 
       Comparison of Cell-Centered Finite Difference Methodology and Prior Art Methods 
       [0053]    A improved method of solving partial differential equations (PDEs) in accordance with the present invention is now described. 
         [0054]      FIGS. 2A to 2E  show grid comparisons of Traditional Finite Difference Methodology (TFDM), Cell-Centered Finite Difference Methodology (CCFDM) in accordance with the present invention, Cell-Centered Finite Volume Methodology (CCFVM) and Vertex-Centered Finite Volume Methodology (VCFVM) across selected illustrated grid arrangements are illustrated. As shown, the four different methodologies: TFDM, CCFDM, CCFVM and VCFVM can be used to solve the illustrated grid arrangement. 
         [0055]    The TFDM requires that the PDE be applied at the node P. The TFDM is configured to use finite difference approximations for the derivatives in the PDE to “discretize” the equation. 
         [0056]    For the CFDM also shown in  FIG. 2A , the PDE is applied at cell centers  1 ,  2 ,  3  and  4 . The stencil is confined to each cell described. Furthermore, finite different approximations are used for the derivatives in the PDE to “discretize” the equation. 
         [0057]    For the CCFVM also shown in  FIG. 2A , the PDE is written in integral form and applied to each cell. The divergence theorem is used to write volume integrals as line integrals around the control volumes. In this manner, the line integrals can be approximated. 
         [0058]    For the VCFVM also shown in  FIG. 2A , the PDE is also written in integral form and applied to the control volume (CV) around P. The divergence theorem is used to write volume integrals as line integrals around the control volume. In this manner, the line integrals can be approximated. 
         [0059]      FIG. 2B  shows additional grid arrangements where all four methodologies may be used. For the TFDM, the grid arrangement must first be mapped to a unit square, as shown. CCFDM, CCFVM and VCFVM may be applied to provide solutions, as described in the previous case shown in  FIG. 2A . 
         [0060]      FIG. 2C  illustrates a grid arrangement on which CCFDM, CCFVM and VCFVM can be applied, but TFDM cannot. 
         [0061]      FIG. 2D  illustrates a grid arrangement where the four different methodologies can be used to solve the illustrated grid arrangement. However, in order to use TFDM, the grid arrangement must first be transformed using multiblock methodology and transformations, as shown in the figure. CCFDM, CCFVM and VCFVM may be applied to provide solutions, as described in the previous case shown in  FIG. 2A . 
         [0062]      FIG. 2E  illustrates curvilinear grid arrangement in which the four different methodologies can be used to solve a PDE. As shown, both TFDM and CCFDM first require coordinate transformations to map the solution domain to a rectangular plane. Once transformed, TFDM and CCFDM may be applied, as before. For CCFVM and VCFVM, no transformations are required. 
         [0063]      FIGS. 3A to 3C  illustrates additional examples of grid arrangements on which CCFDM, CCFVM and VCFVM can be applied, but TFDM cannot. 
         [0064]    When referring to a perfect solution, a problem or grid configuration is selected which has an absolute mathematical solution. The applicant&#39;s method is applied to the same problem or grid and then the results of the two calculations are compared. 
       Cell-Centered Finite Difference Method 
       [0065]    To implement the system of the present invention, an improved Cell-Centered Finite Difference Method (CCFDM) has been proposed. 
         [0066]    To illustrate the applicant&#39;s new numerical approximation process, in one example, a given Partial Differential Equation (PDE), or system of PDEs to be solved on a mesh arrangement containing elements (or cells), is shown in  FIG. 4 . Given the geometry of each cell, i.e., knowing the Cartesian coordinates of the cell vertices, the location of the cell centroids cc 1 , cc 2 , cc 3 , cc 4 , cc 5 , cc 6 , cc 7  is determined. Then, a finite difference stencil is constructed for each cell. This stencil has the unique feature that it is confined to the cell, intersecting the boundary edges of each cell at points w, e, s and n. 
         [0067]    For example, by examining the differencing stencil in cell # 1 , the distances from cc 1  to e and w are shown as not equal. Similarly, the distances from cc 1  to s and n are not equal. This inequality will degrade the accuracy of any central difference formula about the point cc 1 . To overcome this problem, 1D mappings are used from x to ξ and from y to η such that the line segment ‘w-cc 1 -e’ is mapped to a line segment −1≦ξ≦1 where cc 1  is mapped to ξ=0. A similar mapping is used to map the line segment ‘s-cc 1 -n’ to −1≦η≦1, as is shown in  FIG. 5 . 
         [0068]    The PDE, which will be applied at the cell centroid cc 1 , must also be transformed to the computational space. Consider, for example, the model elliptic equation (Poisson eqn.): 
         [0000]    
       
         
           
             
               
                 
                   
                     ∂ 
                     2 
                   
                    
                   T 
                 
                 
                   ∂ 
                   
                     x 
                     2 
                   
                 
               
               + 
               
                 
                   
                     ∂ 
                     2 
                   
                    
                   T 
                 
                 
                   ∂ 
                   
                     y 
                     2 
                   
                 
               
             
             = 
             
               f 
                
               
                 ( 
                 
                   x 
                   , 
                   y 
                 
                 ) 
               
             
           
         
       
     
         [0069]    Under the 1D mappings x=x(ξ), y=y(η), this equation transforms to 
         [0000]    
       
         
           
             
               
                 
                   1 
                   
                     x 
                     ′2 
                   
                 
                  
                 
                   
                     
                       ∂ 
                       2 
                     
                      
                     T 
                   
                   
                     ∂ 
                     
                       ξ 
                       2 
                     
                   
                 
               
               - 
               
                 
                   
                     x 
                     ′′ 
                   
                   
                     x 
                     ′3 
                   
                 
                  
                 
                   
                     ∂ 
                     T 
                   
                   
                     ∂ 
                     ξ 
                   
                 
               
               + 
               
                 
                   1 
                   
                     y 
                     ′2 
                   
                 
                  
                 
                   
                     
                       ∂ 
                       2 
                     
                      
                     T 
                   
                   
                     ∂ 
                     
                       η 
                       2 
                     
                   
                 
               
               - 
               
                 
                   
                     y 
                     ′′ 
                   
                   
                     y 
                     ′3 
                   
                 
                  
                 
                   
                     ∂ 
                     T 
                   
                   
                     ∂ 
                     η 
                   
                 
               
             
             = 
             
               
                 f 
                 _ 
               
                
               
                 ( 
                 
                   ξ 
                   , 
                   η 
                 
                 ) 
               
             
           
         
       
     
         [0070]    If one uses 3-point central differencing to approximate the partial derivatives in this equation, then the resulting difference equation can be written as: 
         [0000]        a   cc   T   cc   =a   w   T   w   a   e   T   e   +a   s   T   s   +a   n   T   n   −f   cc ; 
         [0000]    where the coefficients are expressed in terms of the physical Cartesian coordinates of the w, e, s and n points. This equation can be solved iteratively for the value of T at the cell centroid, assuming we have previous iteration values for T w , T e T s  and T n . 
         [0071]    Step  1 : Create a mesh for the region of interest. Label all nodes N 0 , N 1 , N 2 , etc. ( FIG. 4A ). Establish a fixed reference frame Oxy. Line segments N 1 -N 2 , N 2 -N 3 , - - - , N 6 -N 7  form the interface boundary curve between the solid material (solid region) and the fluid material (fluid region) depicted in this model. The mesh in the model can be arbitrary or user influenced, e.g., the user can apply a finer mesh (smaller size cells) in the areas of the model where variables have high gradients. The finer mesh will result in higher resolution in those areas. 
         [0072]    Step  2 : Select any node in the mesh, and determine the cells sharing that node. For example, in the diagram below, P is a node in the solid region and Q is a node in the fluid region. The cells surrounding P are P-N 1 -N 4 , P-N 4 -N 6 , P-N 6 -N 7 -N 8 -N 9 , etc. (as shown in  FIG. 4B ). The cells surrounding Q are Q-N 17 -N 24 -N 25 , Q-N 25 -N 26 -N 19 , Q-N 19 -N 6 -N 5  and Q-N 5 -N 4 -N 17 . 
         [0073]    Step  3 : For each cell surrounding P (or Q), determine the coordinates of the cell centroids cc 1 , cc 2 , etc. ( FIG. 4C ) 
         [0074]    Step  4 : Within each cell surrounding P (or Q), create a stencil centred at the cell centroid with arms parallel to the x, y, coordinate directions defined by the fixed reference frame, intersecting the cell faces at points w, e, s and n. For example, for node P refer to the cell formed by nodes P-N 13 -N 1  with cell centre cc 1 . For node Q refer to the cell formed by nodes Q-N 17 -N 24 -N 25 . As an alternative to using cell centroids in Steps  3  and  4 , it is possible to determine the coordinates of the point cc′ in the cell which has the property that the length of the line segments w-cc′ and cc′-e are equal and the length of the line segments s-cc′ and cc′-n are equal ( FIG. 4D ). 
         [0075]    Step  5 : For each cell surrounding P (or Q), determine the coordinates of the face intersection points w, c, s and n. 
         [0076]    Step  6 : Repeat Steps  2 - 5  for all nodes in the mesh. 
         [0077]    Step  7 : Select a node P in the mesh at which the dependent variable (T) is to be evaluated, and collect all the cells surrounding P. This node P may be in the solid region, in the fluid region, or on the interface boundary curve. 
         [0078]    Step  8 : For each cell surrounding node P, apply the appropriate mathematical equation (e.g., PDE for solids, or PDE for fluids), defined by the medium in which the cell lies, at the cell centre. Approximate the continuous derivatives in the mathematical equations by standard finite difference formulae, applied on the stencils created in Step  4 , to formulate a discrete approximation to the continuous equations. For each cell, this will result in a finite difference equation of the form 
         [0000]        a   cc   T   cc   +a   w   T   w   +a   e   T   e   +a   s   T   s   +a   n   T   n   =S   cc   (1)
 
         [0000]    if the cell is a solid cell, and of the same mathematical form 
         [0000]        a   cc   T   cc   +a   w   T   w   +a   e   T   e   +a   s   T   s   +a   n   T   n   =S   cc   (1)
 
         [0000]    if the cell is a fluid cell. In these equations the subscripts cc, w, etc., refer to the cell centre, face intersection point w, etc. The coefficients a cc , a w , a c , a s , a n , and the source term S cc  in equations (1) and (2) are not the same. These quantities depend on the nature of the continuous model equation (i.e., whether describing the solid motion or the fluid motion), the differencing scheme used, the cell topology and the coordinates of the face intersection points. Thus, in particular, the physical attributes of the medium, such as thermal conductivity, density, Young&#39;s modulus, Poisson&#39;s Ratio or modulus of elasticity for a solid cell, or such as kinematic viscosity, density, thermal conductivity or specific heat for a fluid cell, are embedded in these coefficients. From the computer&#39;s perspective, for each cell these coefficients are fixed constants and the solution process is identical, regardless of whether the cell is solid or fluid. 
         [0079]    Step  9 : The quantities T w , T e , T s  and T n  in equation (1) or (2) are approximated using an appropriate interpolation scheme based on neighbouring nodal and/or centroid values. These terms are taken to the right-hand side of the equation, and equation (1) or (2) is now approximated by 
         [0000]        a   cc   T   cc   =S   cc   −a   w   T   w   *−a   e   T   e   *−a   s   T   s   *−a   n   T   n *  (3)
 
         [0000]    where the superscript * refers to the approximate value obtained from the interpolation above. 
         [0080]    Step  10 : Equation (3) is solved for the quantity T cc : 
         [0000]    
       
         
           
             
               
                 
                   
                     T 
                     cc 
                   
                   = 
                   
                     
                       
                         S 
                         cc 
                       
                       - 
                       
                         
                           a 
                           w 
                         
                          
                         
                           T 
                           w 
                           * 
                         
                       
                       - 
                       
                         
                           a 
                           e 
                         
                          
                         
                           T 
                           e 
                           * 
                         
                       
                       - 
                       
                         
                           a 
                           s 
                         
                          
                         
                           T 
                           s 
                           * 
                         
                       
                       - 
                       
                         
                           a 
                           n 
                         
                          
                         
                           T 
                           n 
                           * 
                         
                       
                     
                     
                       a 
                       cc 
                     
                   
                 
               
               
                 
                   ( 
                   4 
                   ) 
                 
               
             
           
         
       
     
         [0081]    Step  11 : Repeat Steps  8 - 10  for each cell surrounding P, obtaining the value of T at all surrounding cell centres. 
         [0082]    Step  12 : Determine the value of T at node P by interpolation of the surrounding cell centre values. 
         [0083]    Step  13 : Select a new node P in the mesh and repeat Steps  8 - 12 . Continue until all nodes in the mesh have been updated. This completes one sweep of the mesh. 
         [0084]    The solution process described above] is iterative. Nodal values are repeatedly updated until some prescribed convergence criterion is satisfied. 
       Partial Differential Equations Solution Procedure 
       [0085]    The present system thus provides a preferred partial differential equations procedure shown in the process algorithm of  FIGS. 6A and 6B . P is a typical node in the domain at which the dependent variable is to be evaluated. The PDE solution procedure is as follows: 
         [0000]    a. find all the cells that share the current node (i.e. node P).
 
b. for each one of these cells;
 
i. calculate the cc coordinates, and the coordinates of w, s, e and n intersections.
 
ii. calculate T c  by weighted averaging between the two cc&#39;s that share e (i.e. cc 1  and cc 2 ). Similarly, evaluate T n , T w  and T s .
 
iii. evaluate T cc  from the discretized CCFDM form of the model equation.
 
c. update node P by weighted averaging from all adjacent cell centres.
 
         [0086]    The calculations start with an initial guess at P, which is then updated iteratively until the convergence criterion is satisfied. 
       EXAMPLES 
       [0087]      FIGS. 7A to 7D ,  8 A,  8 B and  9 A to  9 C illustrate graphically first example test cases with uniform structured mesh, clustered structured mesh, unstructured triangulated mesh, and unstructured refined mesh, as well as comparisons between the CCFDM solution, FEM and exact solutions. 
         [0088]    In  FIGS. 7A to 7D , a sample test case using the Laplace equation on a unit square with Dirichlet boundary conditions is illustrated. The Laplace equation for the sample test is described according to the following equations: 
         [0000]    
       
         
           
             
               
                 
                   
                     
                       ∂ 
                       2 
                     
                      
                     T 
                   
                   
                     ∂ 
                     
                       x 
                       2 
                     
                   
                 
                 + 
                 
                   
                     
                       ∂ 
                       2 
                     
                      
                     T 
                   
                   
                     ∂ 
                     
                       y 
                       2 
                     
                   
                 
               
               = 
               0 
             
             ; 
             
               ∀ 
               
                 
                   ( 
                   
                     x 
                     , 
                     y 
                   
                   ) 
                 
                 ∈ 
                 
                   
                     ( 
                     
                       0 
                       , 
                       0 
                     
                     ) 
                   
                   × 
                   
                     ( 
                     
                       1 
                       , 
                       1 
                     
                     ) 
                   
                 
               
             
           
         
       
       
         
           
             
               
                 T 
                  
                 
                   ( 
                   
                     x 
                     , 
                     0 
                   
                   ) 
                 
               
               = 
               
                 
                   T 
                    
                   
                     ( 
                     
                       x 
                       , 
                       1 
                     
                     ) 
                   
                 
                 = 
                 
                   
                     T 
                      
                     
                       ( 
                       
                         1 
                         , 
                         y 
                       
                       ) 
                     
                   
                   = 
                   0 
                 
               
             
             , 
             
               
                 T 
                  
                 
                   ( 
                   
                     0 
                     , 
                     y 
                   
                   ) 
                 
               
               = 
               1 
             
           
         
       
     
         [0089]    Different meshes are shown for the unit square. In  FIG. 7A , a uniform structured mesh is shown. In  FIG. 7B , a clustered structured mesh is shown. In  FIG. 7C , an unstructured (triangulated) mesh is shown. In  FIG. 7D , an unstructured refined mesh is shown. 
         [0090]      FIG. 8A  shows the CCFDM solution to the sample test case using the Laplace equation on a unit square, in accordance with an embodiment of the present invention. As compared to the exact solution shown in  FIG. 8B , the approximated version using CCFDM shows a solution which closely matches the exact solution. 
         [0091]      FIGS. 9A to 9C  show a comparison of CCFDM, FEM and the exact solution along horizontal lines y=0.005, 0.5 and 0.95; respectively. As shown, both the CCFDM and FEM provide solutions which closely match the exact solution. 
         [0092]    A second example was conducted with a coarse mesh and fine mesh for a Poisson equation on a unit disk with Dirichlet boundary conditions. The coarse mesh is shown in FIG.  10 A and the fine mesh is shown in  FIG. 10B . The second example is described according to the following equations: 
         [0000]    
       
         
           
             
               
                 
                   
                     
                       
                         ∂ 
                         2 
                       
                        
                       T 
                     
                     
                       ∂ 
                       
                         x 
                         2 
                       
                     
                   
                   + 
                   
                     
                       
                         ∂ 
                         2 
                       
                        
                       T 
                     
                     
                       ∂ 
                       
                         y 
                         2 
                       
                     
                   
                 
                 = 
                 
                   - 
                   1 
                 
               
               ; 
               
                 ∀ 
                 
                   
                     ( 
                     
                       x 
                       , 
                       y 
                     
                     ) 
                   
                   ∈ 
                   Ω 
                 
               
             
             , 
             
               
                 T 
                  
                 
                   ( 
                   
                     x 
                     , 
                     y 
                   
                   ) 
                 
               
               = 
               0 
             
             , 
             
               ∀ 
               
                 
                   ( 
                   
                     x 
                     , 
                     y 
                   
                   ) 
                 
                 ∈ 
                 
                   ∂ 
                   Ω 
                 
               
             
           
         
       
     
         [0093]    When run, the results were compared to the exact solution. In particular, the relative error (RE) and root mean square error (RMSE) for the above example is shown in Table 1: 
         [0000]    
       
         
               
             
               
               
               
               
             
               
               
               
               
             
           
               
                 TABLE 1 
               
             
             
               
                   
               
               
                 Error compared to Exact Solution 
               
             
          
           
               
                   
                 Mesh 
                 RE 
                 RMSE 
               
               
                   
                   
               
             
          
           
               
                   
                 Coarse 
                 3.43% 
                 6.115e−04 
               
               
                   
                 Fine 
                 2.03 
                 1.842e−04 
               
               
                   
                   
               
             
          
         
       
     
         [0094]      FIGS. 11A to 11C  illustrate third example results by the 21×41 grid shown in  FIG. 11A  with CCFDM solution shown in  FIG. 11B  and the TFDM solution shown in  FIG. 11C . 
         [0095]    The third example was conducted on a square with Dirichlet and Neumann boundary conditions according to the following equations: 
         [0000]    
       
         
           
             
               
                 
                   
                     
                       ∂ 
                       2 
                     
                      
                     T 
                   
                   
                     ∂ 
                     
                       x 
                       2 
                     
                   
                 
                 + 
                 
                   
                     
                       ∂ 
                       2 
                     
                      
                     T 
                   
                   
                     ∂ 
                     
                       y 
                       2 
                     
                   
                 
               
               = 
               0 
             
             ; 
             
               ∀ 
               
                 
                   ( 
                   
                     x 
                     , 
                     y 
                   
                   ) 
                 
                 ∈ 
                 
                   
                     ( 
                     
                       0 
                       , 
                       0 
                     
                     ) 
                   
                   × 
                   
                     ( 
                     
                       4 
                       , 
                       4 
                     
                     ) 
                   
                 
               
             
           
         
       
       
         
           
             
               
                 T 
                  
                 
                   ( 
                   
                     0 
                     , 
                     y 
                   
                   ) 
                 
               
               = 
               80 
             
             , 
             
               
                 T 
                  
                 
                   ( 
                   
                     4 
                     , 
                     y 
                   
                   ) 
                 
               
               = 
               0 
             
             , 
             
               
                 T 
                  
                 
                   ( 
                   
                     x 
                     , 
                     4 
                   
                   ) 
                 
               
               = 
               180 
             
             , 
             
               
                 
                   
                     ∂ 
                     T 
                   
                   
                     ∂ 
                     y 
                   
                 
                  
                 
                   ( 
                   
                     x 
                     , 
                     0 
                   
                   ) 
                 
               
               = 
               0 
             
           
         
       
     
         [0096]    As shown, CCFDM and TFDM provide similar numerical approximations. 
         [0097]      FIGS. 12A to 12E  illustrate graphically fourth example results and relative difference contours on hexagonal ring between CCFDM solution and the FEM solution. In  FIG. 12A , the hexagonal ring is illustrated. 
         [0098]    In  FIGS. 12B ,  12 C and  12 D, a coarse mesh, a clustered mesh near the discontinuities, and a fine mesh are illustrated; respectively. 
         [0099]    In  FIG. 12E , the relative differences between the CCFDM solution and the FEM solution are illustrated. As shown, the relative differences between these methods are low. However, the inventors have appreciated that the numerical approximation provided by CCFDM provides a simpler and faster approach using discretized partial differential equations. The approach requires less memory and may be adaptable to a wider range of applications. 
         [0100]    In particular, CCFDM may be used to solve both compressible fluids and non-compressible fluids in the same equation. 
       Extension of CCFDM Formulation to 3D 
       [0101]    The extension of the CCFDM to 3-dimensional problems is straightforward. Consider, for example, the tetrahedral cell shown in  FIG. 13 . Each face of this 4-faced volume element is triangular in shape. To simplify the discussion, the global Cartesian coordinate system is placed with its origin at one of the vertices of the tetrahedron OABC. Face OAB lies in the xy-plane, face OBC lies in the yz-plane and face OCA lies in the xz-plane. 
         [0102]    For 3 Dimensional problems, the preferred procedure is as follows: 
         [0000]    1. Given the coordinates of A, B and C, calculate the coordinates of the centroid cc of the cell.
 
2. Draw a line through cc parallel to the z-axis, extending it until it intersects two faces of the cell, at points n (on face ABC) and s (on face OAB) in the figure. Determine the coordinates of n and s.
 
3. Draw a line through cc parallel to the y-axis, extending it until it intersects two faces of the cell, at points w (on face OCA) and e (on face ABC) in the figure. Determine the coordinates of w and e.
 
4. Draw a line through cc parallel to the x-axis, extending it until it intersects two faces of the cell, at points f (on face ABC) and b (on face OBC) in the figure. Determine the coordinates of f and b.
 
5. Use three 1D mappings to map the non-uniform stencil in the physical domain to a computational stencil which has uniform spacing in each direction.
 
6. Apply the appropriate finite difference formulas at the cell centroid to discretize the governing PDEs.
 
7. Use interpolation formulae to evaluate the dependent variables at the points n, s, w, e, f and b.
 
8. Use the values obtained in #7 and the discretized equations in #6 to determine the values of the dependent variables at the cell centroid.
 
         [0103]    To determine the solution at a node in 3D space, all cells that share that node are first identified. The above procedure is applied to each of these cells to determine the values at the centroids of these cells. Then, a weighted average of the cell centroid values can be used to determine the nodal value. 
         [0104]    Embodiments of the invention may include various steps as set forth above. While described in a particular order, it should be understood that a different order may be taken, as would be understood by a person skilled in the art. Furthermore, the steps may be embodied in machine-executable instructions. The instructions can be used to cause a general-purpose or special-purpose processor to perform certain steps. Alternatively, these steps may be performed by specific hardware components that contain hardwired logic for performing the steps, or by any combination of programmed computer components and custom hardware components. 
         [0105]    Elements of the present invention may also be provided as a machine-readable medium for storing the machine-executable instructions. The machine-readable medium may include, but is not limited to, floppy diskettes, optical disks, CD-ROMs, and magneto-optical disks, ROMs, RAMs, EPROMs, EEPROMs, magnetic or optical cards, propagation media or other type of media/machine-readable medium suitable for storing electronic instructions. For example, the present invention may be downloaded as a computer program which may be transferred from a remote computer (e.g., a server) to a requesting computer (e.g., a client) by way of data signals embodied in a carrier wave or other propagation medium via a communication link (e.g., a modem or network connection). 
         [0106]    As well, the procedure described above can be implemented on any arbitrary cell topology, ie., any polyhedral shape, and any combination of cell shapes, referred to as hybrid meshes. 
         [0107]    The CCFDM method described above is designed to be applicable to a number of physical problems that can be mathematically modeled by partial differential equations with associated initial conditions (for time-dependent problems) and/or boundary conditions. These include, but are not limited to providing output data and/or the manual or automated computer modelling and/or control of at least the following potential applications:
       steady and unsteady fluid and gas flows   multi-component and multiphase fluid flows   solid mechanics, elasticity, stress analysis   heat conduction   fluid flow and heat transfer   scour simulations   sediment transport   electrostatics, electromagnetics   fluid-structure interaction   multiphysics simulations   cardiovascular flows   higher-order numerical schemes   direct numerical simulation of turbulence       
 
         [0121]    Although this disclosure has described and illustrated certain preferred embodiments of the invention, it is also to be understood that the invention is not restricted to these particular embodiments rather, the invention includes all embodiments which are functional, or mechanical equivalents of the specific embodiments and features that have been described and illustrated herein. Furthermore, the various features and embodiments of the invention may be combined or used in conjunction with other features and embodiments of the invention as described and illustrated herein. The scope of the claims should not be limited to the preferred embodiments set forth in the examples, but should be given the broadest interpretation consistent with the description as a whole. 
         [0122]    As used herein, the aforementioned acronyms shall have the following meanings:
   PDE—Partial Differential Equation   TFDM—Traditional Finite Difference Methodology   CCFDM—Cell-Centered Finite Difference Methodology   CCFVM—Cell-Centered Finite Volume Methodology   VCFVM—Vertex-Centered Finite Volume Methodology   CV—Control Volume   FEM—Finite Element Methodology