Method for the numerical simulation of a physical phenomenon with a preferential direction

• US 7,239,990 B2
• Filed: 02/18/2004
• Issued: 07/03/2007
• Est. Priority Date: 02/20/2003
• Status: Active Grant

• Current Assignee: Robert Struijs
• Original Assignee: Robert Struijs
• Inventors: Struijs, Robert
• Primary Examiner(s): FREJD, RUSSELL WARREN
• Application Number:
  US10/779,685
• Publication Number:
  US 20040167757A1
• Time in Patent Office: 1,231 Days
• Field of Search: 703/2, 703/9, 345/473
• US Class Current: 703/2
• CPC Class Codes: G06F 2111/10 Numerical modelling
• CPC Class Codes: G06F 30/23 using finite element methods [FEM] or finite difference methods [FDM]

Abstract

The invention provides a method of simulating behavior of a flow interacting with an object. The invention improves the accuracy of the approximation of the spatial p^th derivative D_p, p≧1, and therefore reduces the cost of the simulation. It includes in the approximation values of the parameter at points which do not lay on grid lines passing through the point of computation P, and by using these values to optimize the approximation. The use of these additional points depends on a preferential direction, such as determined by the advection direction of the flow. The points used in the approximation extend beyond a unit cube on the computational grid. The simulation is in three or more dimensions. The numerical simulation produces output which can be used in the design or optimization of the object which interacts with the flow. In certain cases, it is the output by itself which is important.

Claims

1. A computation method for simulating behavior of a flow interacting with an object, comprising the steps of:

developing a simulated numerical representation in N dimensions, N≧

3, composed of a plurality of approximated values at a multitude of points in at least a part of space where a flow interacts with an object, the approximated values being of a physical parameter u of the flow to which is associated a velocity field {right arrow over (a)} which determines a preferential direction by executing a numerical scheme wherein at least one spatial p^th derivative D_p, p≧

1, of the parameter u is approximated at the points of the part of space, through the steps of

i) constructing a discrete N-dimensional grid constructed by N families of coordinate lines and using the constructed grid for the part of space,

ii) computing, in at least one point P of the grid, called the point of computation, an approximated value of D_p^A of D_p with an error ε

_n, by using values u_S of the parameter in a collection of grid points, called the stencil S, and computational functions, evaluated with the values u_S, the computational functions depending on the numerical framework in which D_p is expressed, and

iii) choosing the computational functions for the approximated value D_p^A in such a way that the approximated value D_p^A is optimized for the preferential direction,

wherein the stencil S contains at least one point situated outside all the coordinate lines passing through the point of computation P, and the stencil S contains at least a first point and a second point, the first point being defined by N first coordinate lines of the N families of lines, the second point being defined by N second coordinate lines of the N families of lines, and for at least one family N_ƒ

of the coordinate lines, the first coordinate line belonging to the family N_ƒ

is different from and not adjacent to the second coordinate line belonging to the same family N_ƒ

;

outputting the simulated numerical representation that simulates, for the part of space, the behavior of the flow interacting with the object; and

determining the behavior of the flow interacting with the object by using the output simulated numerical representation.

2. The method of claim 1, wherein for varying preferential directions, the approximation D_p^A depends continuously on the values u_S.

3. The method of claim 1, wherein the computing step of D_p^A comprises the steps of:

providing a local basis B({right arrow over (e)}₁,{right arrow over (e)}₂,{right arrow over (e)}₃, . . . ) of curvilinear coordinates which has the unit vector {right arrow over (e)}₁ along the preferential direction, and

choosing the computational functions so that a contribution to the error ε

_n of at least one pure or one mixed derivative as expressed in the local basis B is minimized, while using as a formulation of the values u_S of the parameter at each of the points of the stencil S, a truncated Taylor series expansion with respect to the point of computation P with an error, called the truncation error ε

_S.

4. The method of claim 3, wherein the computational functions are individual coefficients C_S and the approximated value D_p^A is a linear combination of values u_S of the parameter.

5. The method of claim 3, wherein the computing step of D_p^A comprises the steps of:

using an integral formulation for the computation of the derivative D_p with computational functions which are fluxes through a control volume,

computing in at least one volume an approximated value D_p^A of D_p with an error ε

_n, where the approximated value is a function of the flux formulation employed.

6. The method of claim 3, wherein the computing step of D_p^A comprises the steps of:

subdividing the part of space into elements, containing nodes at which the approximation of the physical parameter u is stored,

defining basis computational functions φ

, called interpolation functions, on the elements, where the interpolation functions are used to approximate the physical value u on the element,

computing the integral of the derivative D_p on the element, using a test computational function ψ

, called weighting function,

computing in at least one element an approximated value D_p<sup>φ

,ψ</sup>

of D_p with an error ε

_n where the approximated value is a function of the interpolation functions φ and

the weighting function ψ

.

7. The method of claim 3, wherein the computing step of D_p^A comprises the steps of:

subdividing the part of space into elements, containing nodes at which the approximation of the physical parameter u is stored,

defining basis computational functions φ

, called interpolation functions, on the elements, where the interpolation functions are used to approximate the physical value u on the element,

computing the integral I_el of the derivative D_p on volumes, and distributing parts α

_iI_el to nodes i, where α

_i represent computational functions called distribution coefficients.

8. The method of claim 3, wherein the computing step of D_p^A comprises the steps of:

subdividing the part of space into elements, containing nodes at which the approximation of the physical parameter u is stored,

defining basis computational functions φ

, called interpolation functions, on the elements, where the interpolation functions are used to approximate the physical value u on the element,

computing fluxes ƒ

at the surfaces of volumes, and distributing parts α

_iƒ

to nodes i, where α

_i represent computational functions called distribution coefficients.

9. The method of claim 1, wherein the computing step of D_p^A comprises the steps of:

using a representation of the numerical solutions in Fourier components, and

choosing in the approximated value D_p^A the computational functions in such a way that the Fourier components are optimized for certain directions which are related to the velocity {right arrow over (a)}, while using the values u_S of the parameter at each of the points of the stencil S in the Fourier representation.

10. The method of claim 9, wherein the computational functions are individual coefficients C_S and the approximated value D_p^A is a linear combination of values u_S of the parameter.

11. The method of claim 9, wherein the computing step of D_p^A comprises the steps of:

using an integral formulation for the computation of the derivative D_p with computational functions which are fluxes through a control volume,

computing in at least one volume an approximated value D_p^A of D_p with an error ε

_n, where the approximated value is a function of the flux formulation employed.

12. The method of claim 9, wherein the computing step of D_p^A comprises the steps of:

subdividing the part of space into elements, containing nodes at which the approximation of the physical parameter u is stored,

defining basis computational functions φ

, called interpolation functions, on the elements, where the interpolation functions are used to approximate the physical value u on the element,

computing the integral of the derivative D_p on the element, using a test computational function ψ

, called weighting function,

computing in at least one element an approximated value D_p<sup>φ

,ψ</sup>

of D_p with an error ε

_n, where the approximated value is a function of the interpolation functions φ and

the weighting function ψ

.

13. The method of claim 9, wherein the computing step of D_p^A comprises the steps of:

subdividing the part of space into elements, containing nodes at which the approximation of the physical parameter u is stored,

defining basis computational functions φ

, called interpolation functions, on the elements, where the interpolation functions are used to approximate the physical value u on the element,

computing the integral I_el of the derivative D_p on volumes, and distributing parts α

_iI_el to nodes i, where α

_i represent computational functions called distribution coefficients.

14. The method of claim 9, wherein the computing step of D_p^A comprises the steps of:

subdividing the part of space into elements, containing nodes at which the approximation of the physical parameter u is stored,

defining basis computational functions φ

, called interpolation functions, on the elements, where the interpolation functions are used to approximate the physical value u on the element,

computing fluxes ƒ

at the surfaces of volumes, and distributing parts α

_iƒ

to nodes i, where α

_i represent computational functions called distribution coefficients.

15. The method of claim 1, wherein the N-dimensional grid is expressed in a coordinate system which is chosen from the group consisting of:

rectangular coordinates, spherical coordinates, cylindrical coordinates, parabolic cylindrical coordinates, paraboloidal coordinates, elliptic cylindrical coordinates, prolate spheroidal coordinates, oblate spheroidal coordinates, bipolar coordinates, toroidal coordinates, conical coordinates, confocal ellipsoidal coordinates and confocal paraboloidal coordinates.

16. The method of claim 1, wherein the N-dimensional grid is chosen from the group consisting of:

a grid with non-uniform mesh spacing, a grid which is moving, a grid which is deforming, a grid which is rotating, and a grid which is staggered, and any combination thereof.

17. The method of claim 1, wherein the spatial p^th derivative D_p is a pure derivative

$\frac{{\partial }^{}}{}$

p ⁢

u ∂

e i p .

18. The method of claim 1, wherein p=1.

19. The method of claim 1, wherein p=2.

20. The method of claim 1, wherein the spatial p^th derivative D_p is a mixed derivative

$\frac{{\partial }^{}}{}$

p ⁢

u ∂

e 1 p 1 ⁢

∂

e 2 p 2 ⁢

…

with p₁+p₂+ . . . =p.

21. The method of claim 3, computing an approximated value D₁^A of D₁, where the approximated value is denoted by

$\begin{array}{c}(\begin{array}{c}{\left(u_x\right)}_{i,j,k,}\end{array}\end{array}$

…

( u y ) i , j , k , ⁢

…

( u z ) i , j , k , ⁢

…

⋮

) , ( 27 )

wherein

$\begin{array}{cc}\begin{array}{c}{\left(u_x\right)}_{i,j,k,<br><br>\mathrm{\dots <br><br>}}=<br><br>\frac{1}{\mathrm{\Delta <br><br>}<br><br><br><br>x}<br><br>\{a_{-m}<br><br>u_{i-m,j,k,<br><br>\mathrm{\dots <br><br>}}+a_{-m+1}<br><br>u_{i-m+1,j,k,<br><br>\mathrm{\dots <br><br>}}+\mathrm{\dots <br><br>}+\\ <br><br>a_{n-1}<br><br>u_{i+n-1,j,k,<br><br>\mathrm{\dots <br><br>}}+a_n<br><br>u_{i+n,j,k,<br><br>\mathrm{\dots <br><br>}}+\mathrm{Tx}\},\end{array}<br><br><br><br>\begin{array}{c}{\left(u_y\right)}_{i,j,k,<br><br>\mathrm{\dots <br><br>}}=<br><br>\frac{1}{\mathrm{\Delta <br><br>}<br><br><br><br>y}<br><br>\{a_{-m}<br><br>\left(u_{i-m,j-m,k,<br><br>\mathrm{\dots <br><br>}}-u_{i-m,j,k,<br><br>\mathrm{\dots <br><br>}}\right)+\\ <br><br>a_{-m+1}<br><br>\left(u_{i-m+1,j-m+1,k,<br><br>\mathrm{\dots <br><br>}}-u_{i-m+1<br><br>j,k,<br><br>\mathrm{\dots <br><br>}}\right)+\mathrm{\dots <br><br>}+\\ <br><br>a_{n-1}<br><br>\left(u_{i+n-1,j+n-1,k,<br><br>\mathrm{\dots <br><br>}}-u_{i+n-1,j,k,<br><br>\mathrm{\dots <br><br>}}\right)+\\ <br><br>a_n<br><br>\left(u_{i+n,j+n,k,<br><br>\mathrm{\dots <br><br>}}-u_{i+n,j,k,<br><br>\mathrm{\dots <br><br>}}\right)+\mathrm{Ty}\},\end{array}<br><br><br><br>\begin{array}{c}{\left(u_z\right)}_{i,j,k,<br><br>\mathrm{\dots <br><br>}}=<br><br>\frac{1}{\mathrm{\Delta <br><br>}<br><br><br><br>z}<br><br>\{a_{-m}<br><br>\left(u_{i-m,j-m,k-m,<br><br>\mathrm{\dots <br><br>}}-u_{i-m,j-m,k,<br><br>\mathrm{\dots <br><br>}}\right)+\\ <br><br>a_{-m+1}(u_{i-m+1,j-m+1,k-m+1,<br><br>\mathrm{\dots <br><br>}}-\\ <br><br>u_{i-m+1,j-m+1,k,<br><br>\mathrm{\dots <br><br>}})+\mathrm{\dots <br><br>}+\\ <br><br>a_{n-1}<br><br>\left(u_{i+n-1,j+n-1,k+n-1,<br><br>\mathrm{\dots <br><br>}}-u_{i+n-1,j+n-1,k,<br><br>\mathrm{\dots <br><br>}}\right)+\\ <br><br>a_n<br><br>\left(u_{i+n,j+n,k+n,<br><br>\mathrm{\dots <br><br>}}-u_{i+n,j+n,k,<br><br>\mathrm{\dots <br><br>}}\right)+\mathrm{Tz}\},\\ <br><br>\mathrm{⋮<br><br>}\end{array}& \left(28\right)\end{array}$

wherein a<sub>−

m</sub>≠

0, a_n≠

0, m and n are given integers, m+n>

0, and m+n>

1 if m*n=0, the terms T_x, T_y, T_z, . . . represent the degrees of freedom which are used in the optimization of the approximated value D₁^A, and where indices (i, j, k, . . . ) define the point of computation P on the N-dimensional grid, and Δ

x, Δ

y, Δ

z, . . . denote the mesh spacings of the N-dimensional grid in each coordinate direction.

22. The method of claim 21, computing an approximated value of D₁^A of

$D_1=\frac{\partial }{}$

u ∂

e 1

by a discretization of order M, in which the terms

$\frac{{\partial <br><br>}^{M+1}<br><br>u}{\partial <br><br>e_2^{M_2}<br><br>\partial <br><br>e_3^{M_3}<br><br>\mathrm{\dots <br><br>}}$

with M₂+M₃+ . . . =M+1 are eliminated in the approximated value D₁^A.

23. The method of claim 21, computing an approximated value D₁^A of

$D_1=\frac{\partial }{}$

u ∂

e 1

by a discretization of order M, in which the terms

$\frac{{\partial <br><br>}^{M+1}<br><br>u}{\partial <br><br>e_1^{M_1}<br><br>\partial <br><br>e_2^{M_2}<br><br>\partial <br><br>e_3^{M_3}<br><br>\mathrm{\dots <br><br>}}$

with M₁+M₂+M₃+ . . . =M+1 and M₁<

M+1 are eliminated in the approximated value D₁^A in the case that {right arrow over (e)}₁ is along the x-axis or along diagonals.

24. The method of claim 1, computing an approximated value D₁^A of D₁, where the approximated value is denoted by

$\begin{array}{c}(\begin{array}{c}{\left(u_x\right)}_{i,j,k,}\end{array}\end{array}$

…

( u y ) i , j , k , ⁢

…

( u z ) i , j , k , ⁢

…

⋮

) , ( 29 )

wherein

$\begin{array}{cc}\begin{array}{c}{\left(u_x\right)}_{i,j,k,<br><br>\mathrm{\dots <br><br>}}=<br><br>\frac{1}{\mathrm{\Delta <br><br>}<br><br><br><br>x}<br><br>\{a_{-m}<br><br>u_{i-m,j,k,<br><br>\mathrm{\dots <br><br>}}+a_{-m+1}<br><br>u_{i-m+1,j,k,<br><br>\mathrm{\dots <br><br>}}+\mathrm{\dots <br><br>}+\\ <br><br>a_{n-1}<br><br>u_{i+n-1,j,k,<br><br>\mathrm{\dots <br><br>}}+a_n<br><br>u_{i+n,j,k,<br><br>\mathrm{\dots <br><br>}}<br><br>\},\end{array}<br><br><br><br>\begin{array}{c}{\left(u_y\right)}_{i,j,k,<br><br>\mathrm{\dots <br><br>}}=<br><br>\frac{1}{\mathrm{\Delta <br><br>}<br><br><br><br>y}<br><br>\{a_{-m}<br><br>\left(u_{i-m,j-m,k,<br><br>\mathrm{\dots <br><br>}}-u_{i-m,j,k,<br><br>\mathrm{\dots <br><br>}}\right)+\\ <br><br>a_{-m+1}<br><br>\left(u_{i-m+1,j-m+1,k,<br><br>\mathrm{\dots <br><br>}}-u_{i-m+1<br><br>j,k,<br><br>\mathrm{\dots <br><br>}}\right)+\mathrm{\dots <br><br>}+\\ <br><br>a_{n-1}<br><br>\left(u_{i+n-1,j+n-1,k,<br><br>\mathrm{\dots <br><br>}}-u_{i+n-1,j,k,<br><br>\mathrm{\dots <br><br>}}\right)+\\ <br><br>a_n<br><br>\left(u_{i+n,j+n,k,<br><br>\mathrm{\dots <br><br>}}-u_{i+n,j,k,<br><br>\mathrm{\dots <br><br>}}\right)\},\end{array}<br><br><br><br>\begin{array}{c}{\left(u_z\right)}_{i,j,k,<br><br>\mathrm{\dots <br><br>}}=<br><br>\frac{1}{\mathrm{\Delta <br><br>}<br><br><br><br>z}<br><br>\{a_{-m}(u_{i-m,j-m,k-m,<br><br>\mathrm{\dots <br><br>}}-\left(u_{i-m,j-m,k,<br><br>\mathrm{\dots <br><br>}}\right)+\\ <br><br>a_{-m+1}<br><br>\left(u_{i-m+1,j-m+1,k-m+1,<br><br>\mathrm{\dots <br><br>}}-u_{i-m+1,j+1,k,<br><br>\mathrm{\dots <br><br>}}\right)+\mathrm{\dots <br><br>}+\\ <br><br>a_{n-1}<br><br>\left(u_{i+n-1,j+n-1,k+n-1,<br><br>\mathrm{\dots <br><br>}}-u_{i+n-1,j+n-1,k,<br><br>\mathrm{\dots <br><br>}}\right)\\ <br><br>a_n<br><br>\left(u_{i+n,j+n,k+n,<br><br>\mathrm{\dots <br><br>}}-u_{i+n,j+n,k,<br><br>\mathrm{\dots <br><br>}}\right)\},\\ <br><br>\mathrm{⋮<br><br>}\end{array}& \left(30\right)\end{array}$

wherein a<sub>−

m≠</sub>

0, a_n≠

0, m and n are given integers, m+n>

0, and m+n>

1 if m*n=0, and where indices (i, j, k, . . . ) define the point of computation P on the N-dimensional grid, and Δ

x, Δ

y, Δ

z, . . . denote the mesh spacings of the N-dimensional grid in each coordinate direction.

25. The method of claim 21, wherein

$\begin{array}{c}\begin{array}{c}{\left(u_x\right)}_{i,j,k,}\end{array}\end{array}$

…

= 1 Δ

⁢

⁢

x ⁢

{ 1 2 ⁢

( u i + 1 , j , k , ⁢

…

- u i - 1 , j , k , ⁢

…

) + Tx } , ( u y ) i , j , k , ⁢

…

= ⁢

1 Δ

⁢

⁢

y ⁢

{ 1 2 ⁢

u i + 1 , j + 1 , k , ⁢

…

- u i + 1 , j , k , ⁢

…

+ ⁢

⁢

u i - 1 , j , k , ⁢

…

- u i - 1 , j - 1 , k , ⁢

…

) + Ty } , ( u z ) i , j , k , ⁢

…

= ⁢

1 Δ

⁢

⁢

z ⁢

{ 1 2 ⁢

( u i + 1 , j + 1 , k + 1 , ⁢

…

- u i + 1 , j + 1 , k , ⁢

…

+ ⁢

⁢

u i - 1 , j - 1 , k , ⁢

…

- u i - 1 , j - 1 , k - 1 , ⁢

…

⁢

) + Tz } , ⁢

⋮

( 31 )

26. The method of claim 21, wherein

$\begin{array}{c}\begin{array}{c}\begin{array}{c}{\left(u_x\right)}_{i,j,k,}\end{array}\end{array}\end{array}$

…

= ⁢

1 Δ

⁢

⁢

x ⁢

{ 1 12 ⁢

( u i - 2 , j , k , ⁢

…

- 8 ⁢

u i - 1 , j , k , ⁢

…

+ ⁢

⁢

8 ⁢

u i + 1 , j , k , ⁢

…

- u i + 2 , j , k , ⁢

…

⁢

) + Tx } , ( u y ) i , j , k , ⁢

…

= ⁢

1 Δ

⁢

⁢

y ⁢

{ 1 12 ⁢

( u i - 2 , j - 2 , k , ⁢

…

- u i - 2 , j , k , ⁢

…

- ⁢

8 ⁢

u i - 1 , j - 1 , k , ⁢

…

+ 8 ⁢

u i - 1 , j , k , ⁢

…

+ ⁢

8 ⁢

u i + 1 , j + 1 , k , ⁢

…

- 8 ⁢

u i + 1 , j , k , ⁢

…

- ⁢

⁢

u i + 2 , j + 2 , k , ⁢

…

+ u i + 2 , j , k , ⁢

…

⁢

) + Ty } , ( u z ) i , j , k , ⁢

…

= ⁢

1 Δ

⁢

⁢

z ⁢

{ 1 12 ⁢

( u i - 2 , j - 2 , k - 2 , ⁢

…

- u i - 2 , j - 2 , k , ⁢

…

- ⁢

8 ⁢

u i - 1 , j - 1 , k - 1 , ⁢

…

+ 8 ⁢

u i - 1 , j - 1 , k , ⁢

…

+ ⁢

8 ⁢

u i + 1 , j + 1 , k + 1 , ⁢

…

- 8 ⁢

u i + 1 , j + 1 , k , ⁢

…

- ⁢

⁢

u i + 2 , j + 2 , k + 2 , ⁢

…

+ u i + 2 , j + 2 , k , …

⁢

) + Tz } , ⁢

⋮

( 32 )

27. The method of claim 21, wherein

$\begin{array}{c}\begin{array}{c}{\left(u_x\right)}_{i,j,k,}\end{array}\end{array}$

…

= 1 Δ

⁢

⁢

x ⁢

( 3 2 ⁢

u i , j , k , ⁢

…

- 2 ⁢

u i - 1 , j , k , ⁢

…

+ 1 2 ⁢

u i - 2 , j , k , ⁢

…

+ Tx ) , ( u y ) i , j , k , ⁢

…

= ⁢

1 Δ

⁢

⁢

y ⁢

( - 2 ⁢

( u i - 1 , j - 1 , k , ⁢

…

- u i - 1 , j , k , ⁢

…

⁢

) + ⁢

1 2 ⁢

( u i - 2 , j - 2 , k , ⁢

…

- u i - 2 , j , k , ⁢

…

⁢

) + Ty ) , ( u z ) i , j , k , ⁢

…

= ⁢

1 Δ

⁢

⁢

z ⁢

( - 2 ⁢

( u i - 1 , j - 1 , k - 1 , ⁢

…

- u i - 1 , j - 1 , k , ⁢

…

⁢

) + ⁢

1 2 ⁢

( u i - 2 , j - 2 , k - 2 , ⁢

…

- u i - 2 , j - 2 , k , ⁢

…

⁢

) + Tz ) , ⁢

⋮

( 33 )

28. A method according to claim 1, in three dimensions, obtaining the approximation D_p^A with order M of the derivative

$D_p=\frac{{\partial }^{}}{}$

p ⁢

u ∂

e 1 p 1 ⁢

∂

e 2 p 2 ⁢

∂

e 3 p 3

with p₁+p₂+p₃=p on a grid of given extent from the output of a computer program.

29. A method according to claim 28, wherein in the approximation D_p^A by a discretization of order M terms are eliminated, which are the terms

$\frac{{\partial }^{}}{}$

M + 1 ⁢

u ∂

e 2 M 2 ⁢

∂

e 3 M 3 ⁢

…

with M₂+M₃+ . . . =M+1, or the terms

$\frac{{\partial <br><br>}^{M+1}<br><br>u}{\partial <br><br>e_1^{M_1}<br><br>\partial <br><br>e_2^{M_2}<br><br>\partial <br><br>e_3^{M_3}<br><br>\mathrm{\dots <br><br>}}$

with M₁+M₂+M₃+ . . . =M+1 and M₁<

M+1 in the case that {right arrow over (e)}₁ is along the x-axis or along diagonals.

30. A method according to claim 28, wherein the approximation of D_p^A is a first order accurate discretization.

31. The method of claim 3, computing in three dimensions an approximation of the second derivative

$D_2=\frac{{\partial }^{}}{}$

2 ⁢

u ∂

e 1 2 ,

where the approximation D₂^A is expressed in the terms u_xx, u_yy, u_zz, u_xy, u_yz and u_zx which are given by

$\begin{array}{cc}\begin{array}{c}{\left(u_{\mathrm{xx}}\right)}_{i,j,k}=\frac{1}{{\left(\mathrm{\Delta <br><br>}<br><br><br><br>x\right)}^2}<br><br>\left(u_{i+1,j,k}-2<br><br>u_{i,j,k}+u_{i-1,j,k}+\mathrm{Txx}\right),\\ {\left(u_{\mathrm{yy}}\right)}_{i,j,k}=\frac{1}{{\left(\mathrm{\Delta <br><br>}<br><br><br><br>y\right)}^2}<br><br>\left(u_{i,j+1,k}-2<br><br>u_{i,j,k}+u_{i,j-1,k}+\mathrm{Tyy}\right),\\ {\left(u_{\mathrm{zz}}\right)}_{i,j,k}=\frac{1}{{\left(\mathrm{\Delta <br><br>}<br><br><br><br>z\right)}^2}<br><br>\left(u_{i,j,k+1}-2<br><br>u_{i,j,k}+u_{i,j,k-1}+\mathrm{Tzz}\right),\\ \begin{array}{c}{\left(u_{\mathrm{xy}}\right)}_{i,j,k}=<br><br>\frac{1}{\mathrm{\Delta <br><br>}<br><br><br><br>x<br><br><br><br>\mathrm{\Delta <br><br>}<br><br><br><br>y}<br><br>\{\frac{1}{2}<br><br>(u_{i+1,j+1,k}-u_{i,j,+1,k}-u_{i+1,j,k}+2<br><br>u_{i,j,k}-\\ \begin{array}{c}\\ \end{array}<br><br><br><br>u_{i-1,j,k}-u_{i,j-1,k}+u_{i-1,j-1,k})+\mathrm{Txy}\},\end{array}\\ \begin{array}{c}{\left(u_{\mathrm{yz}}\right)}_{i,j,k}=<br><br>\frac{1}{\mathrm{\Delta <br><br>}<br><br><br><br>y<br><br><br><br>\mathrm{\Delta <br><br>}<br><br><br><br>z}<br><br>\{\frac{1}{4}<br><br>(u_{i+1,j+1,k+1}-u_{i+1,j,k+1}-u_{i+1,j+1,k}+\\ <br><br>u_{i+1,j,k}+u_{i,j+1,k+1}-u_{i,j,k+1}-u_{i,j+1,k}+2<br><br>u_{i,j,k}-\\ <br><br>u_{i,j-1,k}-u_{i,j,k-1}+u_{i,j-1,k-1}+u_{i-1,j-1,k-1}-\\ \begin{array}{c}\\ \end{array}<br><br><br><br>u_{i-1,j-1,k}-u_{i-1,j,k-1}+u_{i-1,j,k})+\mathrm{Tyz}\},\end{array}\\ \begin{array}{c}{\left(u_{\mathrm{zx}}\right)}_{i,j,k}=<br><br>\frac{1}{\mathrm{\Delta <br><br>}<br><br><br><br>z<br><br><br><br>\mathrm{\Delta <br><br>}<br><br><br><br>x}<br><br>\{\frac{1}{4}<br><br>(u_{i+1,j+1,k+1}-u_{i,j+1,k+1}-u_{i+1,j+1,k}+\\ <br><br>u_{i,j+1,k}+u_{i+1,j,k+1}-u_{i+1,j,k}-u_{i,j,k+1}+2<br><br>u_{i,j,k}-\\ <br><br>u_{i-1,j,k}-u_{i,j,k-1}+u_{i-1,j,k-1}+u_{i-1,j-1,k-1}-\\ <br><br>u_{i-1,j-1,k}-u_{i,j-1,k-1}+u_{i,j-1,k})+\mathrm{Tzx}\},\end{array}\end{array}& \left(34\right)\end{array}$

wherein the terms T_xx, T_xy, T_xz, T_yy, T_yz and T_zz represent the degrees of freedom which are used in the optimization of the approximated value D₂^A and where indices (i, j, k) define the point of computation P on the three-dimensional grid, and Δ

x, Δ

y, Δ

z denote the mesh spacings of the three-dimensional grid in each coordinate direction.

32. The method according to claim 1, wherein

D_p^A=Σ

_nL_nD_p,n^A where L_n are constants and each D_p,n^A is a function of values u_S of the parameter in a collection of grid points, called stencil S_n, with individual computation functions, which depend on the numerical framework in which D_p,n^A is expressed, and wherein in the approximation D_p,n^A, the computation functions are chosen in such a way that the approximated value D_p,n^A is optimized for the preferential direction.

33. The method of claim 1, wherein

D_p^A=Σ

_nL_nD_p,n^A where L_n are limiting functions of the values u_S of the stencil S, and at least one D_p,n^A is a function of values u_S of the parameter in a collection of grid points, called stencil S_n, with individual computation functions, which depend on the numerical framework in which D_p,n^A is expressed, and wherein in the approximation D_p,n^A, the computation functions are chosen in such a way that said approximated value D_p,n^A is optimized for the preferential direction.

34. The method of claim 1, wherein the stencil S is chosen from the group consisting of:

upwind discretization stencils, centered discretization stencils, and discretization stencils which are a combination of at least one upwind discretization stencil and at least one centered discretization stencil.

35. The method of claim 1, wherein the numerical discretization is a non-linear discretization.

36. The method of claim 1, wherein the numerical scheme is chosen from the group consisting of:

the Lax-Wendroff scheme, the Lax-Friedrich scheme, the MacCormack scheme, the leap-frog scheme, the Crank-Nicholson scheme, the Stone-Brian scheme, the box scheme, Henn'"'"'s scheme, the QUICK scheme, the κ

scheme, the Flux Corrected Transport (FTC) scheme, the family of ENO schemes, schemes in the class of the Piecewise Parabolic Method (PPM), multi-level schemes, and schemes obtained with the fractional step method and variants thereof.

37. The method of claim 1, wherein the numerical scheme includes the discretization of a plurality of equations.

38. The method of claim 1 for the numerical simulation of physical phenomena which are modeled by the Navier-Stokes equations with equation(s) of state.

39. The method of claim 1 for the numerical simulation of physical phenomena which are modeled by the Euler equations with equation(s) of state.

40. The method of claim 1 for the numerical simulation of physical phenomena which are modeled by the magneto-hydrodynamic equations with equation(s) of state.

41. The method of claim 1 in combination with at least one model for the physical phenomena chosen from the group consisting of:

equations to model turbulence, equations to model chemical reactions, equations to model electromagnetic phenomena, equations to model multiphase flow and equations to model multiphysics phenomena.

42. The method of claim 1 for the numerical simulation of physical phenomena in combination with at least one acceleration technique chosen from the group consisting of:

local time-stepping, multi-grid, GMRES and preconditioning.

43. The method of claim 1 for the simulation of a physical phenomenon which includes a material flow.

44. The method of claim 1 for the simulation of a physical phenomenon which includes a material object interacting with material flow.

45. The method of claim 1, wherein the object is a vehicle.

46. The method of claim 1, wherein the object includes a rotating blade.

47. The method of claim 1, wherein the object is at least a part of the atmosphere of the earth.

48. The method of claim 1, wherein the flow includes oil.

49. The method of claim 1, wherein the flow includes combustion.

50. A data processing system programmed to implement a method according to claim 1.

51. A computer program that can be loaded in a data processing system so as to implement a method according to claim 1.

52. A digital storage computer-readable medium containing a stored computer program that is configured to be loadable in a data processing system so as to control the data processing system to implement a method according to claim 1.