1. BACKGROUND OF THE INVENTION

1.1 Field of Invention

This invention relates to the numerical methods for simulating invicid subsonic flows and solving a class of inverse problems, which is the shape design of aerodynamic body in invicid subsonic flows.

1.2 Prior Art

Most of the existing CFD (Computational Fluid Dynamics) works utilize the Eulerian formulation to solve the fluid flow governing equations, which means in the Eulerian plane constituted by the Cartesian coordinates the computational grid is generated in advance and mainly based on the geometry constraints. The grid forms the computational cells. Across the cell interfaces, the convective flux exists since the fluid particles pass over. It is this convective flux that causes severe numerical diffusion in the numerical solutions because the numerical diffusion is directly associated with the error resulting from numerically approximating the convective term. Since the last century, the primary CFD efforts of algorithm researchers had been concentrated on developing more robust, accurate, and efficient ways to reduce the diffusion. In particular, the upwind methods had a great deal of success in solving fluid flows, because the upwind methods reasonably represent the characteristics of the convective flux. Typically, the Godunov method [1], solving the Riemann problem on the cell interfaces, gives the most accurate results. The FVS (Flux Vector Splitting) method [2], applying the eigenstructure of the equations to treat the flux term on the cell interfaces, is more efficient. However, the convective diffusion, as a result of the fundamental nature of the Eulerian formulation, still exists in those characteristic-based numerical schemes.

On the other hand, Lagrangian description to fluid flow states the motions and properties of the given fluid particles as they travel to different locations. For the governing equations of the fluid flow, such as the Euler equations in Lagrangian formulation, there must be one coordinate representing the streamlines, which can be the stream function numerically. Another coordinate could be the distance of the particles travel. This coordinate system constitutes the Lagrangian plane, where the computational grid points are literally fluid particles and the computational grid is always simply rectangular. In particular, since the particle paths coincide with the streamlines in steady flow and no fluid particles will cross the streamlines, there is no convective flux across the cell interfaces and the numerical diffusion can be thus minimized when the equations are solved in the Lagrangian plane.

One advantage from applying the Lagrangian formulation is when solving a class of inverse problems. The typical inverse problem of aerodynamics is posed by specifying the pressure distribution on the solid-wall of aerodynamic body and determining the geometry of this solid-wall that realizes this pressure distribution. If we solve the problems in the Eulerian plane, using such as the adjoint method [3], we have to firstly estimate the unspecified body shape, then to generate a grid around this shape, and to apply the numerical methods to the flows and find the pressure distribution on the body. Next, a very important and time-consuming step is to solve the adjoint equation to modify the geometry shape. As this process has to be repeated for several successive steps for shape modification until the final target geometry is reached, the computation time of this process is usually long. A Lagrangian formulation, as mentioned before, using the stream function as one coordinate, is more suitable to work on the unspecified geometry problem, since the unspecified shape (solid-wall boundaries) are always represented by streamlines, there no need to modify the computational grid like that in the adjoint method. In the Lagrangian plane, the computational grid will keep on the same in anytime, no matter the geometry of the body shape how to change, because any solid-wall surface in Lagrangian plane is a straight line with the constant values of stream function. The method solving the geometry shape design problem in Lagrangian plane comes to the optimal process (the most efficient process).

Although the Lagrangian formulation presents such excellent properties, it was limited in supersonic flows simulation before [4]. When the Euler equations are numerically solved in Lagrangian plane, they spatially evolve downstream, there is no any upstream information need to consider, which therefore perfectly follows the physical characteristics in supersonic flows. Recently the Lagrangian formulation had been applied on the aerodynamic shape design problems in two-dimensional supersonic flow field [5]. Until recently, for the subsonic or low speed flow domain, using the advantage of the stream function as one coordinate working on the inverse shape design problems is only for potential flows (invicid, irrotational, incompressible) and linearized compressible flows [6].

1.3 Objects and Advantages

There is an apparent need to apply the advantages of the Lagrangian formulation, such as the minimum numerical diffusion in simulation and the optimal process in shape design. Some industrial applications, such as two-dimensional invicid subsonic flow simulations and design cases for nozzle and airfoil, are common, useful and necessary in the preliminary phase of product design. Numerically solving subsonic flows by using the strict Lagrangian concept becomes an excessive obstacle, physically because there exist the upstream-propagating waves in subsonic flows. Thus, the existence of a body located downstream is transmitted to the oncoming fluid particles via the waves so that the particles, sensing the influences of upstream-propagating, can change motion accordingly. A key to the success of a numerical Lagrangian method for subsonic flows lies in how to properly and instantaneously feed these upstream-propagating waves to the particles. Several objects of the present invention are: (a) to provide a Lagrangian formulation of the two-dimensional Euler equations with the zeroed-out convective flux along one coordinate to minimize numerical diffusion; (b) to provide the numerical methods that solve the Euler equations in Lagrangian formulation to find more accurate solutions; (c) to provide the numerical method that is optimal to solve the inverse shape design problems in subsonic flows.

2. DESCRIPTION OF THE INVENTION

2.1 The two-dimensional Euler equations in Lagrangian formulation

This invention firstly provides a coordinate transformation from the Eulerian variables (t, x, y) in the Eulerian plane to the Lagrangian variables (τ, λ, ξ) ϊη the Lagrangian plane, where the variable τ with the same dimension as time term t in the Euler plane, is introduced to function as the 'time-like' marching term for subsonic flows; the other two independent variables are the stream function ξ with dimension and the Lagrangian-distance λ with dimension [m] . The coordinate ξ shall coincide with the streamlines of the fluid particles and λ is specified as the location of different particles along streamlines. This invention stars from the Euler equations in differential conservation form for two-dimensional unsteady inviscid compressible flows in the Eulerian plane

df d¥ dG

-+— + -0 , (1) dt dx dy

where f is the conservation variables vector and F , G are the convective fluxes two Cartesian directions. In detail,

The variables t , u , v , p and p are respectively the time, Cartesian components of flow velocity V , fluid density and pressure. The total energy E and enthalpy H are

2 , 2

U + V

£ = l ( _M 2 + v ² )+-^ and H = E + - + r p

(3)

2 γ - 1 p p 2 γ - 1 p

In above equation, γ is the specific heat ratio. The Jacobian matrix of the transformation from coordinates (t, x, y) to (τ, λ, ξ) can be written as

0 0

d(t,x,y)

cos Θ U (4) sin θ V

and its determinant J becomes

where h is called Lagrangian behavior parameter, Θ is the flow direction angle and the two quantities called the Lagrangian geometry state variables U , V with the dimension

[m ² s ^■ kg ^'1 ] are defined as u = ^, v (6) δξ δξ

Finally, the two-dimensional time-dependent Euler equations (1) have been transformed into the ones in Lagrangian formulation in the Lagrangian plane, which include two set of partial differential equations in the conservation form, 3G,

+ ^■ (?) δτ δλ 3ξ

5G„

^■ + - (8) δτ ΘΛ

where f _L , ¥ _L and G _L have the same definition as f F , G in (2) but in Lagrangian plane,

(10)

with the notations of J in (5), and e = tg f- (1 1)

K =«V - vU , (12) 0 < /z≤ 1 . (13)

The equations (7) are called the Lagrangian physical conservation equations and (8) are the geometric conservation equations. In (9), we find that the first and last element in vector G _L are zero, which means that for the computing grid in Lagrangian plane there apparently are no flux across the cell interfaces along ξ -direction for the continuity and energy equations, so that the numerical diffusion is minimized. There are some other properties of the Euler equations in Lagrangian formulation. They are

/. The equations (7) satisfy the homogeneity property, which implies that the FVS method can be used to resolve the flux term.

2. The equations (7) are not with the rotational invariance property, which implies that it is better to use the different numerical schemes for the convective flux along the two coordinate directions to keep the characteristic non-smeared.

3. The partial differential equation systems (7) and (8) are both strongly-hyperbolic, which implies that it is the well-posed initial problem and they can be solved using the time-like marching scheme. 2.2 Eigenstructure and characteristics of the Euler equations in

Lagrangian formulation

For the physical conservation equations (7), the eigenstructure (eigen values and eigen vectors) can be found, according to the definition of the Jacobian matrix,

d¥ _r dG ,

(14) dQ " dQ

where Q = [p, u, v , p] ^T . Ones can find the four eigenvalues of the Jacobian matrix B , and the four corresponding right eigenvectors

(16) along λ -direction; and the four eigenvalues of the Jacobian matrix C

and the four corresponding right eigenvectors

in ξ -direction, where a is the speed of sound. The characteristic equations of the Lagrangian physical conservation equations (7) can be written as d ,p + , pa Vu— ² +v ² du = 0 , (19) along ¾ = ^ ₄ ^ (l - ^ ² + v ² + ¾)^U ² +V ² ;

Riemann invariants in Lagrangian plane in λ -direction become

the Riemann invariants in -direction are

The equations (15) ~ (22) will be used to construct the numerical scheme for solving the governing equations (7) and (8).

2.3 Construction of numerical scheme

To solve the two-dimensional Euler equations in Lagrangian formulation (7) and (8) in Lagrangian plane, the FVM (Finite Volume Method) is chosen as the spatial discretization. As stated before, the computational grid in Lagrangian plane is simple rectangular. A cell-center scheme will be applied, where the flow quantities are stored at the centers of the grid cells. Based on the properties of the homogeneity, rotation invariance and hyperbolicity of the equations (7) and (8), to keep the efficiency and accuracy of the calculations, this invention, combining the advantages of both the Godunov method and FVS method, supplies a new method of solution, the dimensional- splitting scheme with hybrid upwind flux operators, which means that the convective fluxes across the cell interfaces will be calculated by the FVS method in λ -direction and the Godunov method with local Riemann solver in -direction sequentially. The Strang dimensional-splitting, which is the second-order accuracy in time, is applied for the temporal discretization. Now, let's drop the subscript L with the Lagrangian sense for the convenience reason. 2.3.1 FVS method in λ -direction

Firstly rewrite equations (7) only along λ -direction,

di d¥ _n

— +— = 0. (23) δτ δλ

By applying the FVS (flux vector splitting) method, one obtains the solutions at all the cell (z, j) at the new time-step f - ¹ = f ^* , - - [L^ - ι ,Γ,, + ¾ - ¾ ] , (24) where n is the number of time-step, i,j are the cell indices, Arthe time-step, Αλ is the spatial step. The split matrix

V = K _Lcy A ⁺ K l andL- = K _icv A ^~ K ^;Lc ¹ v ' (25) with the split matrix of eigenvalues

and the matrix

1 1 1

Va Va

u u

J Vu ² +V ² Vu ² +v ²

Ua Ua

v- v + -

J Vu ² +V ² Vu ² +v ²

2

W + v Ka a w ² +v ² Ka a

(Uw+Vv)^ + -

J 2 Vu ² +v ² i- 2 Vu ² +v ² i-

(27) and the inverse of it reads

-1

K Lev 2.3.2 Godunov method with Riemann solver in ξ -direction

For the Godunov method, the key ingredient is the solution of the Riemann problem. This invention solves the Riemann problem along ξ -direction, which is called the Riemann problem across streamlines. Figure 1 shows the definition of the Riemann problem across streamlines. Rewrite 7) only in -direction,

where f _L and f _R are the conservation variables at the left and right state of the cell interface respectively. From now, the subscripts L are for the left state, R for the right, and * for the star for the Riemann problem. At the left or right state, shock and expansion waves could be. Between them, it is the star state (region), which can be divided into left and right star states again like shown in Figure 1. The solutions of the new time-step can be obtain from the discretized equations of (29),

The task is to find the values of the primitive variables p,p,u,v from f to get G at cell interfaces j ± 1/ 2 , where the primitive variables are the variables in the star region in the Riemann problem. A Riemann solver was built as the following procedure. Connecting the left and right states to the star state by integrating along the characteristics equations (19) and (20), then ones have

where C _L = p _L a _L , C _R = p _R a _R , and f (w, v) =— ^{v ¾2 + y2} + 1 u ln(v + . (35) f (u, v) is called the combination function, which can be considered as the combination of the velocity component (u, v) in Lagrangian plane and have the same role as the velocity term with dimension [m I s] . The (35) has its polar form as i ,e) =— (tg0 + cos Θ ln((l + sin ø)) + VlnV (36) with V = u = VcosO and v = VsinO . The solutions of (31)~(34) present the primitive variables in the star state

To find the velocity components («» , v, ) either in left or right star state, the equation (38) need to be solved. Firstly, one can rewrite (38) as

with the known conditions in left and right states,

^C L ^{■ U} L' ^V L)+ ^C R ^■ 1( ^U R> ^V R)+[PL ^~ PR]

B (42) c _L + c

The combination function (35) can be further modified. One can define

tg0, =e, (43) where θ„ is the flow angle either in left or right star state, then u„ - V _* cos Θ _* , v» =V, sin Θ _* (44)

+ ε" VT+ ε

where V„ representing the velocity magnitude either in the left or right star state. Further, the equation (41) can be finalized as a function of ε

The work needs to do is numerically to solve the equation F(^)=0 with a given velocity V„ to find ε for θ„ . The equation (45) does have the differentiability and its first-order derivative of F is Numerical experiments show that for the non-dimensional velocity V _* < 5 (subsonic flows) and f e [- l, l] , F '(s) > 0 , which means that the function F(^) is monotone within this domain; in the mean time, f (- ε) · Ρ(ε) < 0 , so that the Newton-Raphson iteration as a root-finding algorithm is a good way to find the root ε for the equation (45) with a given V„ . To find the value of V„ , the non-linear waves (shear, shock and expansion waves) in Riemann problem have to be recovered with the approximate values of p„ , p„ _L , p„ _R . In this invention, the following relations are considered:

Rankine-Hugoniot relations across shocks

If a shock wave appears between one state (say left) and the corresponding star state, across this shock, the Rankine-Hugoniot relations can be used to the left state and star state. They are,

where μ _! , is the shock wave speed, J„ _L and E„ _L are J (transformation Jacobian) and E (total energy) in the star state respectively. From (46) and (49), one receives and

Velocity absolute value V. _tL or V. _tR can be substituted into (45) to receive 6, _L οτθ,

Enthalpy constants across expansion waves

If an expansion wave appears between one state (say left) and the corresponding star state, the connection between the two states can be found by considering the enthalpy constant across the waves. Ones have VL + (52) γ - \ ρ L 2 Y - p _* r from which the velocity magnitude can be found as

and

The selection of the non-linear waves is based on the wave-type conditions, in which the values of pressure in left, right and star states can be found to judge what non-linear wave could appearing in the Riemann problem. Assuming

ί» =™(ΐΊ. ¾ ) ; (55) /'max = vasx{p _L , p _R ), (56) which, as well as p, , form the following wave-type choosing conditions:

(1) If p _min < p, < P _max , which means in the Riemann problem one shock wave in the state with p _min and one expansion wave in the state with ^ .

(2) If p, > p _max , which means there are two shock waves in Riemann problem.

(3) If p, < p _min , which means there are two expansion waves the in Riemann problem.

After <¾ is obtained, we can find G _L according to (9) and u, _L , v, _L oru, _R , ν _*Λ from (44). Those values, as well as p, , can be used to compute the flux term in Godunov method (30). At the same time, the geometry conservation equations (8) can be updated

n+l/ 2 n+l/ 2

^U i,j+l/ 2 ^U i,j-\I 2

(57) n+l / 2 n+l/ 2

^V ! , /+l/ 2 ^{~ V} i -V 2

Although this step is an additional to the method in the Eulerian formulation, it is a trivial step and will no cost more computing time.

2.3.3 Boundary conditions

Two layers of cells should be put, which are called the ghost cells, surrounding the whole computational domain, to make it possible to extend any spatial discretization scheme beyond the boundaries. Solid-wall boundary condition

In the selected scheme, the Riemann problem at the solid- wall boundary in -direction has to be solved, where the left or right state in Riemann problem consists of a mirror image with respect to the solid-wall. In the solid-wall boundary conditions, where subscript L means the quantities p , p , u , , Θ in the ghost cells; while subscript R means any ones in the boundary cells. 0 _W is the angle of the solid-wall. Based on the solutions of the Riemann problem given in (37)-(40), on the solid-wall,

P _* = P , (58)

The pressure in the star state can be considered as the pressure on the solid-wall p _w

P _W = P _* - (60) If the pressure on the stationary wall (p _w ) had been specified, the specified pressure can be directly used as the pressure in the star state (region) in the Riemann solver,

P* = P _W - (61)

In/out flow boundary condition

For the subsonic inlet boundary, the first relation corresponds to the positive, incoming, characteristics. Hence the values at the boundary, indicated by a subscribe b, are obtained from the Riemann invariance in (21)~(22). They are

where u _∞ , a _∞ are the velocity and the speed of sound for the free stream; the subscript i refers to a value at an internal point. The second relation is associated with a numerical boundary condition and has to be estimated from inside the domain by an appropriate extrapolation. Hence, the velocity u _b is obtained by adding equation (62) and (63), then u _b = " ' ■ (64)

A quite similar way can be used to treat the subsonic outflow boundaries, with the difference that the one s ecified variable at the boundary is p _∞ (the pressure at infinite),

Boundary conditions for Lagrangian geometry state variables

Another boundary conditions need to implement are the Lagrangian geometry state variables U andV . For the solid-wall boundary,

V _z = cos(2^ - ¾) . (67)

For in/outlet boundaries, we specify the flow angle Θ =0, then,

U = 0 , (68) V = 1 . (69)

A flow chart presenting the numerical procedures for one step of updating from [Q]" _v to Q\" ^w ^ ^tn second-order accuracy in time and space, where Q = [p, u, v , /?] ^r are the physical parameters of the flow field, is given in Figure 2.

2.4 Method for solving a class of inverse problems

Since using the Lagrangian formulation, the solid-wall boundaries are always represented by streamlines, which will be the stream function in Lagrangian plane, where any streamlines represented by the computing grids are straight lines. Exploring the advantage of the Lagrangian formulation on a class of inverse problem is another purpose in this invention. From now, this class problem is defined as the inverse geometry shape design problem.

In Section2, we introduced the solutions of the local Riemann problem in ^-direction in Lagrangian plane and the solid-wall boundary condition. From (37) and (57), ones have

where p _R ,p _R ,a _R ,V _R ,6 _R are the known state variables in boundary and f is the combination function in triangle form given in (35) and (36). From (70), ones have

which can be numerically solved with the specified p _w to get the flow angle <¾ ίη ghost cells. Using the mirror boundary condition, the wall angle is

3, = !(<¾ ⁺ <¾) · (72)

The inverse geometry shape design means finding the solid-wall angle based on the specified pressure distribution on the solid-wall. The recovering of the solid-wall angle needs to solve the inversed Riemann problem here. An inverse Riemann solver is built as the following procedure. One can define

then we have the velocity components in ghost cells u _L = V _R cos0 _L (74)

+ ε" T+ ε

which can be substituted into the combination function in triangle form (71), which can be further modified and finalized as a function of ε

The work needs to do is numerically to solve the equation F (£ ^■ ) = () with a given physical parameters p _R ,p _R ,V _R , 0 _R and p _w to find ε for 6 _L . The equation (75) does have the differentiability and its first-order derivative of F in respect to ε is

Similar to the method for (45), the updated value ε in new step can be expressed as

The flow angle in ghost cells can be recovered by 6 _L = tg ^~x [p ^n+x The solid-wall angle 0 _W can be found by using equation (72).

The flow chart for solving the inverse geometry design problem is given in Figure 3. This method starts from the guessed flow filed physical values and solid-wall angle, the inversed Riemann solver with the specified pressure distribution applied on the solid- wall boundary will give the flow angles in ghost cells, then the wall angle can be obtained and substituted back to the flow solver for updating the physical variables. This procedure will continue until an iterated solid-wall angle within the convergence tolerance is reached. It is obvious that the convergent criteria is the solid-wall angle instead of any physical parameter, which means that the solid-wall angle become a flow field variable. The change of the boundary each iteration doesn't change the characteristics of the equations and the equations are still the strong-hyperbolic system; therefore, it is still a well-posed initial problem. Compared to the methods used in the Eulerian plane, such as the adjoint method, the present method needs neither the step of the grid generation for the updated geometry shape, nor the step of solving the adjoint equation. It obtains the convergent physical solutions in flow field and the geometry shape concurrently, which will be the one-shot and optimal method in the inverse geometry shape design problem. The total CPU time for this class problem will reduce the order of one.

3. DRAWINGS - Figures

Figure 1 the Riemann problem across streamline along -direction in Lagrangian plane

Figure 2 flow chart of one time-step for updating the flow field using the second-order Scheme: dimensional-splitting with hybrid flux operators. Wherein the reference numbers represent the equations or variables, ® [Q]" · , © [Q]" _{±l/ 2 >} © [F]"±i/ 2j >

1 / 3

+ 2 /3 _ jn + 1 / 3 _ ^

, © [Q]"±l/ 2,/ > ® [F]" _± ±;1/ 2,7 '

Αξ

Figure 3 flow chart of an inverse method for geometry shape design

Figure 4 the Lagrangian and Eulerian solutions for parabolic nozzle flow t _in = 0.5

(a) Grid in Lagrangian plane

(b) Grid in Eulerian plane

(c) Streamlines by the Lagrangian solution (d) Streamlines by the Eulerian solution

(e) Comparison of the pressure distributions by Lagrangian and exact solution on the bottom solid-wall boundary

(f) Grid built by the Lagrangian solution

Figure 5 calculated geometry shape of the solid- wall based on the prescribed pressure distribution

(a) Grid in Eulerian plane

(b) Pressure distribution on the solid wall

(c) Designed solid-wall

4. DETAILED DESCRIPTION - PREFEREED EMBODIMENT

4.1 Invicid subsonic internal flow passing a divergent nozzle

According to the numerical methods discussed above, a complete solution to the invicid subsonic internal flow by using the Euler equations in Lagrangian formulation will be performed. This example is about the subsonic flow passing through a two-dimensional parabolic divergent nozzle with length L and inlet height H _in , whose geometry is defined as the two parts of arabolic curves

where H _in -

nozzle inlet is M _in - 0.5 . The flow is the pure subsonic flow. The Euler equations in

Lagrangian formulation (7) and (8) will be solved by the numerical scheme presented in Section 2, where the Strang dimensional-splitting of second-order accuracy in time and hybrid flux operators with second-order upwind MUSCL extrapolation with minmod limiter are used. The computing parameters: CFL = 0.48; the Lagrangian behavior parameter h = 0.98 ; totally 60x20 uniform cells in Lagrangian plane. The computing results will be compared with those obtained from the JST method [8] based on FVM in Eulerian plane and the exact solutions. Now, denoting the solution obtained in Lagrangian plane by the invented method as the Lagrangian solutions and the results obtained from the Eulerian plane by the Eulerian formulation as the Eulerian solutions. In detail, the computations are taken as the following steps with referring to Figure 2, (1) Initialization. The flow field variables Q° = [/? ⁰ , w ⁰ ,v ⁰ , ?°] ^r are known as the initial conditions as well as U ⁰ ,V ⁰ , where Q° take the values at the infinite and

U ° =0; V ° =1. The superscript 0 means that the initial conditions of a flow problem are given at t = 0 ( r = 0 ) in x - j plane and λ - ξ plane. Subsequently, the conservation variables f° are available. Then an appropriate rectangular λ— ξ coordinate mesh with uniform spatial-step ( Δλ , Δξ ) is formed. The corresponding mesh in x - y plane is laid on according to

dx - cos θ άλ and dy - sin θ άλ (79)

(2) Calculating the time-step. Based on the variables at time-step /? , (n - 0,1,2,... ), an FL<0.5, calculate the time-step,

where S , T = cos θ ^η , where θ ^η is from previous time-step, then

J = UT -VS , (81)

(3) Using the FVS method in λ - direction to solve equations (23). The conservation variables are updated from f" . to f" ^{+1 2} by the known values u, v, p] " . and ll " V. " . For the second-order schemes, the MUSCL extrapolation with

TVD limiter has to be taken.

(4) Using the Godunov method in ξ - direction to solve equation (29). With all f" ^+1/2 and Q" ^+1/2 known at time step n {n - 0,1,2,... ), the Godunov method will be used to update the conservation variables. A complete 2-D local Riemann solver gives the primitive variables u, v, p] _iJ±V2 at the interfaces of the neighboring cells. For the second-order schemes, the MUSCL extrapolation with TVD limiter has to be taken to decide the initial values for the Riemann problem.

(5) Updating the geometry variables. Since the velocities [«, v] "^ ₊ ^ at the interface of cells are available, then the discretized geometrical conservation equations can be updated to new time step at once by (57), whose accuracy in time and space is decided by the accuracy of the velocities [ , v] " ^~ _+y2 , though the discretized equations are presented in the first-order accuracy. (6) Finding the physical primitive variable values at the new time step. Decoding from the conservation variables f _{t j} , = /Ί,/ ₂ 4 ]"^ to obtain the primitive variables [p,p, u, v] " ⁺¹ . The primitive variables read as

(6) Grid building. The last step in one time-step is to build the grid if needed. The grid points are given by

according to which each grid point represents the center of a fluid particle representing by a computational cell. The principle of the time-dependent

Lagrangian approach in this invention tells that along ξ - direction, the grid lines built here actually are the streamline themselves.

(7) Loop (2). After the numerical procedure for advancing one time-step is complete, to march forward further, one goes back to (2) and repeat (3)-(7), until the conservation variables or primitive variables go to convergence.

The Lagrangian solution starts from the rectangular grid with Lagrangian distance in λ - direction and stream function in ^-direction in Figure 4 (a). The streamlines shown in Figure 4(c) by Lagrangian solution apparently agreed well with Eulerian solution in Figure 4(d), and the pressure distributions on the solid-wall shown in Figure 4(e) totally agreed with the exact solution. The grid used in the Eulerian coordinates in Figure 4(b) has little difference with the grid built by the Lagrangian solution, where the lines along x-direction are streamlines and the lines along y-direction are not perpendicular to x- direction to preserve the length of the streamlines. It is worth to indicate that the final grids in Figure 4(f) evolve from the initial grid in Figure 4(a), and it is the result of the Lagrangian solution since the grid lines are the streamlines themselves. Although one does not need this grid in the computing, an evolving procedure can give a clear understanding of the evolution of the Lagrangian solution. 4.2 Inverse geometry shape design case

Next case is an example to show the detail how to find the geometry of a solid-wall by solving the two-dimensional Euler equation in Lagrangian formulation. It is an inverse shape design problem. The geometry of the example is a convergent-divergent nozzle. Since the prescribed pressure distribution on the geometrically-unspecified solid-wall need to be input in the inverse problem, one can firstly specify this nozzle's shape and find the pressure distribution, then use this pressure as the input to find the geometry shape by the invented method. The calculated shape should match with the original shape. The profile of the nozzle's solid-wall is given by a function

γ = Η _Λ β ^βχ - 1 < JC < 1 , (84) where H _th is taken to 0.1 meaning the 10% bump high of the inlet and β is taken 6.5 to control the smoothed convex part which counts about a third of the channel length. For an inlet Mach number of 0.5, it is a pure subsonic flow. The grid generation of this geometry in the Lagrangian plane is identical to that in the previous case shown in Figure 4(a). The computational grid in Euler plane is given in Figure 5(a). The Euler equations in Lagrangian formulation (7) and (8) will be solved by the numerical scheme presented in Section 2, where the Strang dimensional-splitting of second-order accuracy in time and hybrid flux operators with second-order upwind MUSCL extrapolation with minmod limiter are used. The computing parameters: CFL = 0.48; the Lagrangian behavior parameter h = 0.98 ; two kinds grid, totally 60x20 and 90X30 uniform cells in Lagrangian plane. To obtain the specified pressure distributions on the solid- wall, the flow field calculations in the Lagrangian and Eulerian plane firstly are implemented like in the previous case. Figure 5(b) shows the pressure distribution on the solid-wall by the Lagrangian and the Eulerian solution, which will be used as the prescribed pressure distribution in the inverse geometry shape design problem. Also, the adjoint method will be used to solve this inverse design problem to test the computing efficiency of the invented method. The computations are taken as the following steps, in detail:

(1) Initializing the flow field and the solid- wall angle. The initialed flow field variables are set as the same as in the previous case. The solid-wall angle is assumed as 6 _W = 0° . Since the flow angle on the whole flow field are assumed as zero, the flow angle in the ghost cells is zero 9 _L = 0° .

(2) Imputing the prescribed pressure distribution on the solid- wall p _w . (3) Solving the inversed Riemann problem across the streamline on the solid-wall boundary based on p _w , according to the equations (70)~(77) in Section 2.

(4) Calculating the solid-wall angle. The solid-wall angle can be calculate based on the mirror state condition, then we have 0 _W = !(<¾ + <¾ ) · (84)

(5) Checking the convergence of the solid-wall angle. The variations of the solid-wall angle is the difference of 0 _W between the time step n+1 and n

M _w = {e -ei). (85)

If S0 _W is smaller than the convergent tolerance, then the flow field and solid- wall angle come to the final solution ^ ⁺¹ . Otherwise, continue the next steps.

(6) Updating the flow field. The task of this step is to calculate the flow field physical parameters from [oj" _j to [(J]"* ¹ based the updated solid-wall angle 6^ ⁺¹ . This step includes the step (1) ~ (7) introduced in the previous case. After this step, the new flow field parameters p _R , p _R , V _R , 6 _R ,) are obtained.

(7) Repeating (2).

Figure 5(c) shows the graphic comparison of the computed and original solid- wall shapes. From the results, we can numerically evaluate the precision of the resolution. The computed geometrical shape was found in good agreement with the original shape used as the input; the maximum error is limited within 0.5%. The totally CPU time is only tenth of that by the adjoint method.

5. Conclusions

This invention supplies a new technology to simulate the invicid subsonic flows. It comprises a transformation to transfer the two-dimensional Euler equations into the ones in Lagrangian formulation in Lagrangian plane, which will not only accept a simple rectangular grid, but also produce the minimum convective diffusion when the characteristic-based numerical methods are use to solve the equations. A two- dimensional Riemann solver is created to solve the Riemann problem across streamlines in Lagrangian plane. An advantage of using this new technology lies on solving a class of inverse problem of aerodynamic shape design. Solving the inverse Riemann problem on the solid-wall boundary greatly reduces the complicity of this kind problem, and the computing time reduces as much as the order of one. It gives the flow field solution and the geometry shape concurrently. So that, solving this class of inverse problem comes to the most efficient.

6. Reference

(1) Godunov, S.K.: A difference scheme for numerical computation discontinuous solution of hydrodynamic equations. Translated US Publ. Res. service, JPRS 7226, 1969.

(2) Van Leer, B.: Flux vector splitting for the 1990 's. Invited Lecture for the CFD Symposium on Aero propulsion, Cleveland, Ohio, 1990.

(3) Jameson, A.: Aerodynamic design via control theory. Journal of Scientific Computing, 3:233- 260, 1988.

(4) Loh, C.Y. and Liu, M.S.: A new Lagrangian method for three-dimensional steady supersonic flows. Journal of Computational Physics, 113, pp.224-248, 1994.

(5) Mateescu, D.: Analysis of aerodynamic problem with geometrically unspecified boundaries using an enhanced Lagrangian method. Journal of Fluids and Structures, 17, pp.603-626, 2003.

(6) Stanitz, D. John: A review of certain inverse methods for the design of duct with 2-or 3- dimensional potential flow. Appl. Mech. Rev. Vol. 41, No.6, June 1988.

(7) Jameson, A., Schmidt, W., Turkel, E.: Numerical solutions of the Euler equations by finite volume methods using Runge-Kutta time-stepping schemes. AIAA-paper, 81-1259., 1981.