A pseudo-compressibility method for solving inverse problems based on the 3D incompressible Euler equations

Abstract

A numerical technique to solve the three-dimensional inverse problems that arise in aerodynamic design is presented. The approach, which is well established for compressible flows, is extended to the incompressible case via artificial compressibility preconditioning. The modified system of equations is integrated with a characteristic-based Godunov method. The solution of the inverse problem is given as the steady state of an ideal transient during which the flow field assesses itself to the boundary conditions, which are prescribed as design data, by changing the boundary contour. The main aspects of the Eulerian-Lagrangian numerical procedure are illustrated and the results are validated by comparing with theoretical solutions and experimental results.

Keywords:

• inverse problems
• aerodynamic design
• pseudo-compressibility methods

1 Introduction

Aerodynamic design consists of finding the most convenient shape for a fluid-immersed body that realizes some expected aerodynamic performances. There is a certain freedom in approaching aerodynamic design. Several routes were exploited in the past by researchers and engineers, with the same goals in mind. The most-followed approach iterates on a sequence of direct analyses. Starting from an initial configuration, the geometry is modified and the aerodynamic effects are either experimentally or numerically evaluated. Thanks to designer experience and intuition this process is driven toward the expected performances. In Computational Fluid Dynamics (CFD), the automation of this approach is known as the classical shape optimization.[1–3] This strategy defines a functional to measure the merit or cost of a certain aerodynamic solution, and a number of geometric control parameters that are used to modify the solution. Then, the functional gradient relative to the controls is evaluated numerically and its information is used to update the control parameters and to march towards the functional extremum. Several methodologies have been proposed in literature to evaluate the functional gradient and also to impose fluid-dynamic and geometrical constraints. Another way of solving the aerodynamic problem is based on the formulation of an inverse problem. In this case, the aerodynamic surface is unknown and some flow features are given on it. The solution of the inverse problem enables designers to obtain the geometry of an aerodynamic device satisfying the governing equations with all the given boundary conditions and other desirable flow constraints. One typical inverse problem is that of finding the airfoil geometry, given the flight speed and the pressure distribution on the profile.[4] The solution of the inverse problem is often a complex task, if compared to the direct analysis approach. Nevertheless, it offers the main advantage of requiring a fewer number of flow-field evaluations to determine the enquired geometry, and also it allows for an easier imposition of the flow constrains (e.g. no flow separation and adverse pressure gradient control [5, 6]). Approaches to the inverse problem solution are based on the potential flow theory and conformal mapping techniques,[4, 7, 8] on stream-function base formulations,[9–11] on boundary elements replacing the body surface [12] and on direct design methods.[13–16] Several examples of inverse problem solution methodologies are presented in [5, 17] and references therein.

A generalizable way of solving aerodynamic inverse problems identifies the problem solution as the steady state of an unsteady evolution during which the flow field tries to accommodate itself to the design data, which are prescribed as boundary conditions on the unknown surface.[18] To do this, the surface is allowed to move. To give an idea of the physical scenario involved let us consider, for instance, the design of an airfoil. First, an arbitrary initial geometry is guessed and the related flow field at the selected flight Mach number is evaluated. Then, the target pressure distribution along the profile is imposed. In general, the boundary conditions will not be in agreement with the inner flow field. If we let the profile surface move while it remains impermeable to the flow, walls start moving to decrease the pressure gap. If the system reaches a steady state, and this feature relies on the existence of the steady solution for the Euler equations, then the enquired geometry and the corresponding flow field are found. A similar approach for the determination of the unknown shape is the elastic membrane motion concept [19] or the transpiration velocity.[20–22]

One drawback of inverse problems is that they may be ill posed (according to the Hadamard’s criteria), or they can lead to unacceptable solutions from a designer point of view. If certain wall pressure distributions are required on airfoils, the result is an open or self-intersecting profile. Lighthill discovered the solvability conditions that have to be respected by pressure distributions within an incompressible potential flow model,[4] whereas the conditions for compressible flows and other issues were investigated in [8, 23, 24] and references therein. The main advantage of numerical optimization over inverse problem is that the first allows the maximization or minimization of global quantities, such as lift or drag, in the presence of constraints, whereas for inverse problems the design is limited to the selection of the pressure or velocity distribution on the boundary,[25] which is given on the basis of designer experience and therefore, somewhat arbitrarily. In addition, a limited control is possible on the final geometry. Optimization techniques based on the adjoint method can be adopted to drive inverse problems towards the maximization or minimization of target functionals.[6] From this point of view, the numerical solution of an inverse problem becomes attractive, as an alternative route to shape optimization and to automated design.

In the present work, a time-dependent technique of solving inverse problems [18, 24] is extended to the aerodynamic design in three-dimensional incompressible flows. The numerical method is based on the artificial compressibility approach to make the system of equations virtually hyperbolic.[26, 27] The choice of this preconditioning method allows us to easily set up an explicit, Godunov-type numerical scheme [28] and to reuse well-established techniques of solving the compressible Euler and Navier Stokes equations, as well as the related strategies of code optimization and parallelization. Although our investigation has been here limited to inviscid flows, the general methodology has been already extended to viscous flows and validated for the case of two-dimensional and axis-symmetric compressible flows.[29]

The plan of the paper is as follows: the mathematical model and the numerical technique are explained in the context of a characteristic-based finite volume method; then the equation of motion of impermeable surfaces is derived in either cartesian or cylindrical coordinates; finally, some numerical examples are presented and the accuracy of the method is studied through comparisons against theoretical solutions and experimental results.

2 The mathematical model

The pseudo-compressibility method is a preconditioning technique of overcoming the numerical issues related to the uncoupling of pressure and velocity field in the incompressible Euler equations. The governing equations are made artificially hyperbolic by adding a density time-derivative to the continuity equation.[26–28] Based on physical analogies, the time derivative of density is then linked to the pressure by introducing an artificial compressibility factor $c$.

In a Cartesian frame of reference, the system of modified Euler equations is written as1 $$\begin{array}{c}\hfill \frac{\mathit{\partial }\mathit{U}}{\mathit{\partial }t}+\mathbf{\Theta }\left(\frac{\mathit{\partial }\mathit{F}}{\mathit{\partial }x}+\frac{\mathit{\partial }\mathit{G}}{\mathit{\partial }y}+\frac{\mathit{\partial }\mathit{H}}{\mathit{\partial }z}\right)=0\end{array}$$1 where2 $$\begin{array}{c}\hfill \mathit{U}=\left\{\begin{array}{c}p\\ u\\ v\\ w\end{array}\right\}\mathbf{\Theta }=\left[\begin{array}{cccc}{c}^2& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]\end{array}$$2 and3 $$\begin{array}{c}\hfill \mathit{F}=\left\{\begin{array}{c}u\\ p+u^2\\ uv\\ uw\end{array}\right\}\mathit{G}=\left\{\begin{array}{c}v\\ uv\\ p+v^2\\ vw\end{array}\right\}\mathit{H}=\left\{\begin{array}{c}w\\ uw\\ vw\\ p+w^2\end{array}\right\}\end{array}$$3 as usual, $p$ is the pressure, $\mathit{q}=\{u,v,w\}$ is the velocity vector. All the variables are normalized to the reference length $l_{ref}$, density ${\mathit{\rho }}_{ref}={\mathit{\rho }}_o$ and pressure $p_{ref}$.

The unsteady solutions of system (11 $$\begin{array}{c}\hfill \frac{\mathit{\partial }\mathit{U}}{\mathit{\partial }t}+\mathbf{\Theta }\left(\frac{\mathit{\partial }\mathit{F}}{\mathit{\partial }x}+\frac{\mathit{\partial }\mathit{G}}{\mathit{\partial }y}+\frac{\mathit{\partial }\mathit{H}}{\mathit{\partial }z}\right)=0\end{array}$$1 ) does not have a physical meaning, but from the numerical point of view, the methods of solution for compressible flow become available. At the steady state, the solution of (11 $$\begin{array}{c}\hfill \frac{\mathit{\partial }\mathit{U}}{\mathit{\partial }t}+\mathbf{\Theta }\left(\frac{\mathit{\partial }\mathit{F}}{\mathit{\partial }x}+\frac{\mathit{\partial }\mathit{G}}{\mathit{\partial }y}+\frac{\mathit{\partial }\mathit{H}}{\mathit{\partial }z}\right)=0\end{array}$$1 ) also satisfies steady Euler equations. The proposed strategy of solving the inverse problem requires the numerical integration of a three-dimensional flow field on a time-dependent domain. Unsteady flow field are solved, in the context of artificial compressibility methods, by introducing a dual time stepping technique.[28] Dual time stepping is here unnecessary, since we are interested in over-relaxing the system to the steady state. Moreover, the time-dependent derivative of the metrics can be neglected without a significant loss of numerical stability.

According to the Gauss formula, an integral form of Equation (11 $$\begin{array}{c}\hfill \frac{\mathit{\partial }\mathit{U}}{\mathit{\partial }t}+\mathbf{\Theta }\left(\frac{\mathit{\partial }\mathit{F}}{\mathit{\partial }x}+\frac{\mathit{\partial }\mathit{G}}{\mathit{\partial }y}+\frac{\mathit{\partial }\mathit{H}}{\mathit{\partial }z}\right)=0\end{array}$$1 ) is4 $$\begin{array}{c}\hfill \frac{\mathit{\partial }}{\mathit{\partial }t}{\int }_{\mathrm{\Omega }}\mathit{U}d\mathrm{\Omega }+\mathbf{\Theta }{\int }_{\mathrm{\Sigma }}(\mathit{F}\mathit{i}+\mathit{G}\mathit{j}+\mathit{H}\mathit{k})·\mathit{n}d\mathit{\sigma }=0\end{array}$$4 where $\mathrm{\Sigma }$ is the boundary of the volume $\mathrm{\Omega }$ and $\mathit{n}$ the outward normal. Equation (44 $$\begin{array}{c}\hfill \frac{\mathit{\partial }}{\mathit{\partial }t}{\int }_{\mathrm{\Omega }}\mathit{U}d\mathrm{\Omega }+\mathbf{\Theta }{\int }_{\mathrm{\Sigma }}(\mathit{F}\mathit{i}+\mathit{G}\mathit{j}+\mathit{H}\mathit{k})·\mathit{n}d\mathit{\sigma }=0\end{array}$$4 ) is approximated using a finite volume technique by discretizing the $(x,y,z)$ space with hexahedral cells whose shape depends on time. The integration in time is carried out according to a Godunov type two-step scheme. A first-order Flux Difference Splitting (FDS) is used at the predictor step: the variables ($p$, $u$, $v$ and $w$) are assumed as an averaged, constant value inside each cell. The fluxes $\mathit{F}$, $\mathit{G}$ and $\mathit{H}$ are evaluated by solving the Riemann problems pertinent to the discontinuities that take place at the cell interfaces. At the corrector level, the second order of accuracy is achieved by assuming a linear behaviour of the primitive variables inside the cells, according to an essentially non-oscillatory (ENO) procedure.[30] Implicit residual smoothing is applied to increase the numerical stability. The resulting scheme is second-order accurate in space.[28] In next section, we treat extensively the characteristic-based method of imposing the unconventional boundary conditions used to solve the inverse problem. The solution of the Riemann Problem for internal cells and standard boundary conditions is resumed in Appendix 1.

2.1 Dynamics of materially moving boundaries.

Walls on which a flow variable, e.g. the pressure, is imposed, are treated as flexible and impermeable surfaces whose motion is tracked in time until a steady state is reached. In inviscid flows, the velocity vector must be tangent to the wall surface, which is considered a flow surface. It therefore moves materially within the fluid. In this context, the no-through-flow condition represents the kinematic constraint from which the motion of the wall surface is derived.

Let us consider a new Cartesian coordinate system $\mathbf{\Gamma }=(x_1,x_2,x_3)$. A generic surface moving in time can be represented, for instance, by equation $x_3=b(x_1,x_2,t)$. Without any loss of generality, let us assume that $b(x_1,x_2,0)=0$. Thank to the change of variable $B=x_3-b$ , the surface is now defined by the manifold5 $$\begin{array}{c}\hfill B(x_1,x_2,x_3,t)=0\end{array}$$5 This manifold moves materially within the fluid, that is6 $$\begin{array}{c}\hfill \frac{dB}{dt}=\frac{\mathit{\partial }B}{\mathit{\partial }t}+\mathit{q}·\mathrm{\nabla }B=0\end{array}$$6 and, therefore7 $$\begin{array}{c}\hfill \frac{\mathit{\partial }B}{\mathit{\partial }t}=-\mathit{q}·\mathrm{\nabla }B=u_3-u_1\frac{\mathit{\partial }b}{\mathit{\partial }{x}_1}-u_2\frac{\mathit{\partial }b}{\mathit{\partial }{x}_2}\end{array}$$7 where $\mathit{q}=(u_1,u_2,u_3)$ is the flow velocity expressed in the $\mathbf{\Gamma }$ reference system.

In the following, Equation (77 $$\begin{array}{c}\hfill \frac{\mathit{\partial }B}{\mathit{\partial }t}=-\mathit{q}·\mathrm{\nabla }B=u_3-u_1\frac{\mathit{\partial }b}{\mathit{\partial }{x}_1}-u_2\frac{\mathit{\partial }b}{\mathit{\partial }{x}_2}\end{array}$$7 ) is coupled to a properly defined Riemann problem. The resulting system is used to compute the fluxes at the boundary and to evaluate the wall velocity. The latter is then integrated in time to find the new wall contour.

For the inverse problem on a nozzle or diffuser configuration, as shown in Figure 1(a), the pressure distribution is imposed on the unknown boundary. Aerodynamic devices as wings or blades, instead, are composed by two surfaces that represent a closed contour. In this case, the inverse problem formulation must ensure a correct profile closure and non-overlapping conditions. Among possible approaches, in addition to suitable boundary conditions at infinity, one may:

(i)	prescribe the distribution of thickness and pressure jump along the chord and enquire for the geometry of the camber line;
(ii)	prescribe the distribution of thickness and pressure on one side of the profile and enquire for the geometry of the camber line;
(iii)	prescribe the pressure distribution around the profile and enquire for its geometry.

Case (iii) is the most challenging since, when selecting the pressure distribution, additional closure constraints must be imposed.[23] Moreover, the convergence rate is low because smaller time-steps are required and also boundary crossing must be checked at each step.

All of the mentioned cases can be grouped in two approaches in imposing the boundary conditions: the case of a surface on which a pressure distribution is imposed and the case of two aerodynamic surfaces subjected to a prescribed pressure jump. In the following sections, the kinematic constraint (77 $$\begin{array}{c}\hfill \frac{\mathit{\partial }B}{\mathit{\partial }t}=-\mathit{q}·\mathrm{\nabla }B=u_3-u_1\frac{\mathit{\partial }b}{\mathit{\partial }{x}_1}-u_2\frac{\mathit{\partial }b}{\mathit{\partial }{x}_2}\end{array}$$7 ) is accomplished to the flow-field boundary conditions in order to derive the fluxes and the velocity on the unknown surfaces.

Figure 1. Inverse problem on a three-dimensional nozzle. (a) Schematic view of the imposed boundary conditions; (b) wave pattern of the Riemann problem at the unknown aerodynamic surface on which the pressure is prescribed.

2.1.1 Boundary condition for moving surfaces with imposed pressure.

Without loss of generality, we consider the aerodynamic design of a subsonic nozzle configuration, as shown in Figure 1. Being the upper wall, the unknown surface, on it we assign the pressure distribution $p_w(x,y)$. On that boundary, an half-Riemann problem has to be solved to derive the fluxes at the cell interface. The resulting wave pattern is shown in Figure 1(b). The Riemann invariants across fields (a) and (c) are8 $$\begin{array}{c}\hfill R_{3c}=R_{3a},R_{2v_c}=R_{2v_a}\to {\tilde{v}}_c={\tilde{v}}_a,R_{2w_c}=R_{2w_a}\to {\tilde{w}}_c={\tilde{w}}_a\end{array}$$8 where $\tilde{u},\tilde{v},\tilde{w}$ are the normal, tangent and binormal components of the velocity vector in the local frame of reference on the cell interface. The pressure level of field (c) is assigned as design variable $p_c=p_w(x,z)$ . The velocity ${\tilde{u}}_c$ is deduced from the first of conditions (88 $$\begin{array}{c}\hfill R_{3c}=R_{3a},R_{2v_c}=R_{2v_a}\to {\tilde{v}}_c={\tilde{v}}_a,R_{2w_c}=R_{2w_a}\to {\tilde{w}}_c={\tilde{w}}_a\end{array}$$8 ) and is given by9 $$\begin{array}{c}\hfill {\tilde{u}}_c=\frac{{(\tilde{u}-\tilde{a})}_a{p}_w-R_{3a}}{c^2}\end{array}$$9 being $\tilde{a}=\sqrt{{\tilde{u}}^2+c^2}$ the pseudo-velocity of sound.

Once the field (c) is known, the surface velocity is computed from Equation (77 $$\begin{array}{c}\hfill \frac{\mathit{\partial }B}{\mathit{\partial }t}=-\mathit{q}·\mathrm{\nabla }B=u_3-u_1\frac{\mathit{\partial }b}{\mathit{\partial }{x}_1}-u_2\frac{\mathit{\partial }b}{\mathit{\partial }{x}_2}\end{array}$$7 ) as10 $$\begin{array}{c}\hfill \frac{\mathit{\partial }b}{\mathit{\partial }t}={\tilde{u}}_c\left(n_z-\frac{\mathit{\partial }b}{\mathit{\partial }x}{n}_x-\frac{\mathit{\partial }b}{\mathit{\partial }y}{n}_y\right)+{\tilde{v}}_c\left(m_z-\frac{\mathit{\partial }b}{\mathit{\partial }x}{m}_x-\frac{\mathit{\partial }b}{\mathit{\partial }y}{m}_y\right)+{\tilde{w}}_c\left(l_z-\frac{\mathit{\partial }b}{\mathit{\partial }x}{l}_x-\frac{\mathit{\partial }b}{\mathit{\partial }y}{l}_y\right)\end{array}$$10 The same procedure holds for case (ii) with minor modifications. Once either the upper or lower surface of the profile has been selected, fluxes and boundary contour velocity are evaluated by Equations (88 $$\begin{array}{c}\hfill R_{3c}=R_{3a},R_{2v_c}=R_{2v_a}\to {\tilde{v}}_c={\tilde{v}}_a,R_{2w_c}=R_{2w_a}\to {\tilde{w}}_c={\tilde{w}}_a\end{array}$$8 )–(1010 $$\begin{array}{c}\hfill \frac{\mathit{\partial }b}{\mathit{\partial }t}={\tilde{u}}_c\left(n_z-\frac{\mathit{\partial }b}{\mathit{\partial }x}{n}_x-\frac{\mathit{\partial }b}{\mathit{\partial }y}{n}_y\right)+{\tilde{v}}_c\left(m_z-\frac{\mathit{\partial }b}{\mathit{\partial }x}{m}_x-\frac{\mathit{\partial }b}{\mathit{\partial }y}{m}_y\right)+{\tilde{w}}_c\left(l_z-\frac{\mathit{\partial }b}{\mathit{\partial }x}{l}_x-\frac{\mathit{\partial }b}{\mathit{\partial }y}{l}_y\right)\end{array}$$10 ) and its shape integrated in time. Concerning the other surface, since Equation (1414 $$\begin{array}{c}\hfill {\left(\frac{\mathit{\partial }b}{\mathit{\partial }t}\right)}_{up}={\left(\frac{\mathit{\partial }b}{\mathit{\partial }t}\right)}_{lw}\end{array}$$14 ) holds, at each step the geometry is deduced from the thickness distribution, while the fluxes can be computed using the moving wall boundary condition.

2.1.2 Boundary conditions for inverse problem on surfaces with imposed pressure jump.

In this section, we derive the boundary conditions to solve the inverse problem over three-dimensional profiles according to formulation (i). We consider the design of a stator blade, but the same procedure still holds, with simplifications, for 3D wings. The blade boundary is split into the upper and lower surface and represented by the equations ${\mathit{\vartheta }}_u=b_u(x,r,t)$ and ${\mathit{\vartheta }}_d=b_d(x,r,t)$ in a cylindrical frame of reference $(x\mathit{i},r\mathit{k},\mathit{\vartheta }\mathit{\eta })$.

The upper $(up)$ and lower $(lw)$ blade surfaces represent two boundaries on which two half-Riemann problems have to be formulated to impose boundary conditions. These Riemann problems are also coupled together by the prescribed pressure jump distribution and by a kinematic condition.

If we ideally cross the blade along the $\mathit{\vartheta }$-direction, we can identify states ($a$) and ($b$) as the flow conditions close to the lower and upper surfaces (see Figure 2). On the lower surface, the Riemann invariants across fields ($a$) and ($c$) $\equiv$ ($lw$) are11 $$\begin{array}{c}\hfill R_{3c}=R_{3a},R_{2v_c}=R_{2v_a}\to {\tilde{v}}_c={\tilde{v}}_a,R_{2w_c}=R_{2w_a}\to {\tilde{w}}_c={\tilde{w}}_a\end{array}$$11 On the upper surface, across fields ($b$) and ($d$) $\equiv$ ($up$), we also have12 $$\begin{array}{c}\hfill R_{1d}=R_{1b},R_{2v_d}=R_{2v_b}\to {\tilde{v}}_d={\tilde{v}}_b,R_{2w_d}=R_{2w_b}\to {\tilde{w}}_d={\tilde{w}}_b\end{array}$$12 The wall pressure is constrained, since the pressure jump is a given function $\mathrm{\Delta }p(x,r)$, that is,13 $$\begin{array}{c}\hfill p_{up}-p_{lw}=\mathrm{\Delta }p(x,r)\end{array}$$13 Moreover, since blade thickness is assigned, the velocities of the upper and lower surfaces have to be the same14 $$\begin{array}{c}\hfill {\left(\frac{\mathit{\partial }b}{\mathit{\partial }t}\right)}_{up}={\left(\frac{\mathit{\partial }b}{\mathit{\partial }t}\right)}_{lw}\end{array}$$14 To impose this kinematic constraint we have to derive the analogous of Equation (77 $$\begin{array}{c}\hfill \frac{\mathit{\partial }B}{\mathit{\partial }t}=-\mathit{q}·\mathrm{\nabla }B=u_3-u_1\frac{\mathit{\partial }b}{\mathit{\partial }{x}_1}-u_2\frac{\mathit{\partial }b}{\mathit{\partial }{x}_2}\end{array}$$7 ) in the cylindrical frame, which is15 $$\begin{array}{c}\hfill \frac{\mathit{\partial }b}{\mathit{\partial }t}=-c_x\frac{\mathit{\partial }b}{\mathit{\partial }x}-c_r\frac{\mathit{\partial }b}{\mathit{\partial }r}+\frac{c_{\mathit{\theta }}}{r}\end{array}$$15 The Euler equations are integrated using formulation (11 $$\begin{array}{c}\hfill \frac{\mathit{\partial }\mathit{U}}{\mathit{\partial }t}+\mathbf{\Theta }\left(\frac{\mathit{\partial }\mathit{F}}{\mathit{\partial }x}+\frac{\mathit{\partial }\mathit{G}}{\mathit{\partial }y}+\frac{\mathit{\partial }\mathit{H}}{\mathit{\partial }z}\right)=0\end{array}$$1 ), so that the evaluation of the fluxes is done on a Cartesian local frame. Therefore, three reference frames are used:

(a)	an absolute, Cartesian frame 1 $(x\mathit{i},y\mathit{j},z\mathit{k})$ on which the computational domain is described;
(b)	a local relative frame 2 $({\mathit{\xi }}_1\mathit{n},{\mathit{\xi }}_2\mathit{m},{\mathit{\xi }}_3\mathit{l})$ on which the fluxes at the cell interfaces are evaluated;
(c)	a cylindrical frame 3 $(x\mathit{i},r\mathit{k},\mathit{\vartheta }\mathit{\eta })$ on which the blade motion is represented.

We assume that the axial and radial directions coincide with the Cartesian $x$ and $z$-directions, respectively. First, the velocity components $c_x$, $c_{\mathit{\vartheta }}$, $c_r$ in the cylindrical frame 3 are expressed as a function of the velocity components $\tilde{u}$, $\tilde{v}$, $\tilde{w}$ in the local frame 2 as follows,16 $$\begin{array}{c}\hfill c_x=M_1\tilde{u}+K_1,c_{\mathit{\theta }}=M_2\tilde{u}+K_2,c_r=M_3\tilde{u}+K_3\end{array}$$16 with17 $$\begin{array}{c}\hfill \begin{array}{cc}{M}_1=n_x\hfill & K_1=m_x\tilde{v}+l_x\tilde{w}\hfill \\ M_2=n_{y}cos\mathit{\theta }-n_{z}sin\mathit{\theta }\hfill & K_2=(m_{y}cos\mathit{\theta }-m_{z}sin\mathit{\theta })\tilde{v}+(l_{y}cos\mathit{\theta }-l_{z}sin\mathit{\theta })\tilde{w}\hfill \\ M_3=n_{y}sin\mathit{\theta }+n_{z}cos\mathit{\theta }\hfill & K_3=(m_{y}sin\mathit{\theta }+m_{z}cos\mathit{\theta })\tilde{v}+(l_{y}sin\mathit{\theta }+l_{z}cos\mathit{\theta })\tilde{w}\hfill \end{array}\end{array}$$17 and $\mathit{n}=(n_x,n_y,n_z)$,   $\mathit{m}=(m_x,m_y,m_z)$,   $\mathit{l}=(l_x,l_y,l_z)$ are orthogonal unit vectors. Then, substituting Equations (1515 $$\begin{array}{c}\hfill \frac{\mathit{\partial }b}{\mathit{\partial }t}=-c_x\frac{\mathit{\partial }b}{\mathit{\partial }x}-c_r\frac{\mathit{\partial }b}{\mathit{\partial }r}+\frac{c_{\mathit{\theta }}}{r}\end{array}$$15 ) and (1616 $$\begin{array}{c}\hfill c_x=M_1\tilde{u}+K_1,c_{\mathit{\theta }}=M_2\tilde{u}+K_2,c_r=M_3\tilde{u}+K_3\end{array}$$16 ) into Equation (1414 $$\begin{array}{c}\hfill {\left(\frac{\mathit{\partial }b}{\mathit{\partial }t}\right)}_{up}={\left(\frac{\mathit{\partial }b}{\mathit{\partial }t}\right)}_{lw}\end{array}$$14 ), we have18 $$\begin{array}{c}\hfill {\tilde{u}}_d{E}_1-{\tilde{u}}_u{E}_2=E_3\end{array}$$18 where19 $$\begin{array}{cc}\hfill E_1& {\displaystyle ={\left[\frac{M_2}{r}-M_1\frac{\mathit{\partial }b}{\mathit{\partial }x}-M_3\frac{\mathit{\partial }b}{\mathit{\partial }r}\right]}_{lw}}\hfill \\ \hfill E_2& {\displaystyle ={\left[\frac{M_2}{r}-M_1\frac{\mathit{\partial }b}{\mathit{\partial }x}-M_3\frac{\mathit{\partial }b}{\mathit{\partial }r}\right]}_{up}}\hfill \\ \hfill E_3& {\displaystyle ={\left[\frac{K_2}{r}-K_1\frac{\mathit{\partial }b}{\mathit{\partial }x}-K_3\frac{\mathit{\partial }b}{\mathit{\partial }r}\right]}_{up}-{\left[\frac{K_2}{r}-K_1\frac{\mathit{\partial }b}{\mathit{\partial }x}-K_3\frac{\mathit{\partial }b}{\mathit{\partial }r}\right]}_{lw}}\hfill \end{array}$$19 Finally, by collecting the conditions (1111 $$\begin{array}{c}\hfill R_{3c}=R_{3a},R_{2v_c}=R_{2v_a}\to {\tilde{v}}_c={\tilde{v}}_a,R_{2w_c}=R_{2w_a}\to {\tilde{w}}_c={\tilde{w}}_a\end{array}$$11 )–(1414 $$\begin{array}{c}\hfill {\left(\frac{\mathit{\partial }b}{\mathit{\partial }t}\right)}_{up}={\left(\frac{\mathit{\partial }b}{\mathit{\partial }t}\right)}_{lw}\end{array}$$14 ), we are led to the following equation system for the two unknown states ($up$) and ($lw$)20 $$\begin{array}{c}\hfill \left[\begin{array}{cccc}1& -1& 0& 0\\ 0& {(\tilde{u}-\tilde{a})}_b& 0& -c^2\\ {(\tilde{u}+\tilde{a})}_a& 0& -c^2& 0\\ 0& 0& E_1& -E_2\end{array}\right]\left\{\begin{array}{c}{p}_{up}\\ p_{lw}\\ {\tilde{u}}_{up}\\ {\tilde{u}}_{lw}\end{array}\right\}=\left\{\begin{array}{c}\mathrm{\Delta }p\\ R_{3a}\\ R_{1b}\\ E_3\end{array}\right\}\end{array}$$20 whose solution is21 $$\begin{array}{cc}\hfill p_{up}& {\displaystyle =\frac{(e_2{R}_{1b}-c_b^2{E}_3)-e_1(R_{3a}+\mathrm{\Delta }p{(\tilde{u}-\tilde{a})}_a)}{e_2{(\tilde{u}+\tilde{a})}_b-e_1{(\tilde{u}-\tilde{a})}_a}}\hfill \\ \hfill p_{lw}& =p_{up}-\mathrm{\Delta }p\hfill \\ \hfill {\tilde{u}}_{lw}& {\displaystyle =\frac{{(\tilde{u}-\tilde{a})}_a{p}_d-R_{3a}}{c^2}}\hfill \\ \hfill {\tilde{u}}_{up}& {\displaystyle =\frac{{(\tilde{u}+\tilde{a})}_b{p}_u-R_{1b}}{c^2}}\hfill \end{array}$$21 The new shape is obtained by integrating in time the surface evolution equation (1515 $$\begin{array}{c}\hfill \frac{\mathit{\partial }b}{\mathit{\partial }t}=-c_x\frac{\mathit{\partial }b}{\mathit{\partial }x}-c_r\frac{\mathit{\partial }b}{\mathit{\partial }r}+\frac{c_{\mathit{\theta }}}{r}\end{array}$$15 ). Let us track the motion of the upper blade surface first. From the knowledge of the ($up$) state the corresponding flow velocities ${(c_x,c_r,c_{\mathit{\vartheta }})}_{up}$ are computed by equation (1616 $$\begin{array}{c}\hfill c_x=M_1\tilde{u}+K_1,c_{\mathit{\theta }}=M_2\tilde{u}+K_2,c_r=M_3\tilde{u}+K_3\end{array}$$16 ). The derivative $\mathit{\partial }b/\mathit{\partial }x$ and $\mathit{\partial }b/\mathit{\partial }r$ are approximated by an up-wind reconstruction.[18] Equation (1515 $$\begin{array}{c}\hfill \frac{\mathit{\partial }b}{\mathit{\partial }t}=-c_x\frac{\mathit{\partial }b}{\mathit{\partial }x}-c_r\frac{\mathit{\partial }b}{\mathit{\partial }r}+\frac{c_{\mathit{\theta }}}{r}\end{array}$$15 ) is then integrated numerically by a classical fourth-order Runge-Kutta scheme and the new shape of the upper blade surface is obtained. Finally, the lower surface can be either evaluated by a similar integration or deduced from the knowledge of the thickness distribution.

Figure 2. Inverse problem with pressure jump imposed across a blade. Representation of the wave pattern of the upper ($up$) and lower ($lw$) half-Riemann problems.

2.2 Remarks on problem well-posedness

One drawback of inverse problems is that they may be ill-posed. While appropriate boundary conditions can be defined easily in the case of the inverse design of channels, a special treatment is required when dealing with more complex aerodynamic surfaces. For instance, if certain wall pressure distributions are required on airfoils, the result is an open or self-intersecting profile. Lighthill discovered the solvability conditions with respect to pressure distributions within an incompressible potential flow model,[4] whereas the similar integral conditions for compressible flows and other issues were investigated in Refs. [8, 23, 24]. The way of overcoming these issues depends on the strategy used to impose the closure of the body contour. This problem has already investigated in the proposed literature and can be addressed in two ways: one can satisfy the closure conditions or reformulate the inverse problem in order to automatically satisfy such constraints.[24]

As far as the proposed numerical approach is concerned, the problem we formulated to track the unsteady motion of the aerodynamic surface is always well-posed. Nevertheless, if the underlying inverse problem is ill-posed, the procedure simply does not converge to a steady solution. In other cases, the procedure may fail because of lack on problem uniqueness. If more than one solution could be admissible for the same inverse problem, the numerical procedure may be triggered towards one of the possible configurations. Often the problem uniqueness is recovered by changing a boundary condition.[24, 29] Finally, we remark that to the hyperbolic Equation (77 $$\begin{array}{c}\hfill \frac{\mathit{\partial }B}{\mathit{\partial }t}=-\mathit{q}·\mathrm{\nabla }B=u_3-u_1\frac{\mathit{\partial }b}{\mathit{\partial }{x}_1}-u_2\frac{\mathit{\partial }b}{\mathit{\partial }{x}_2}\end{array}$$7 ) is associated an initial condition $b(x_1,x_2,t)=b(x_1,x_2,0)$ that defines a locus on the unknown aerodynamic surface that does not change during integration. Among the possible choices we selected one of the edges, on that the surface appears to be fastened ‘as a flag’ during its transient motion. For blades, we selected the leading edge locus. Other choices (e.g. the trailing edge locus or a generic 3D stacking line) should be regarded as different initial conditions for the problem that, in general, lead to a different solution. This topic, anyway, belongs to the field of optimal design rather than to the inverse problem solution itself.

3 Numerical results

In next sections, the accuracy of the method is addressed. First, the numerical results of the inverse procedure is validated against two analytical solutions of the Euler equations. Then, accuracy and reliability of the procedure are checked for the more complex case of blade design. The inverse procedure is validated by using numerical results and experimental data (for the E/TU-1 Test Case).

Figure 3. The inverse problem on a three-dimensional nozzle. (a) Initial nozzle geometry with a wavy upper boundary.(b) Final nozzle geometry as computed by the inverse method. A slice with pressure isolines is also depicted. (c) Comparisons of the wall geometry and of the static pressure distributions with the corresponding target analytic values. (d) Time history of the $L^2$-norm residual and of maximum wall velocity.

3.1 Three-dimensional nozzle

The inverse design of a nozzle has been chosen as test case for the inverse problem formulation described in Section 2.1.1. The pressure distribution imposed on the unknown aerodynamic surface is derived from the potential flow theory. The planar solution is extended in three dimensions by applying periodic boundary conditions in the transversal direction $z$.

A converging–diverging nozzle is obtained from the flow field of a standing vortex on the half plane. The flow can be represented by the complex potential $W$22 $$\begin{array}{c}\hfill W(\mathit{\zeta })=u_{\infty }\mathit{\zeta }-i\mathit{\kappa }log\left(\frac{\mathit{\zeta }-ia}{\mathit{\zeta }+ia}\right),\mathrm{\Psi }(\mathit{\zeta })=\mathfrak{I}(W(\mathit{\zeta }))\end{array}$$22 where $\mathit{\zeta }=z+iy$ is the position, $u_{\infty }$ is the asymptotic velocity, $\mathit{\kappa }$ is the vortex strength and $a$ is the distance of the vortex from the $x$-axis. The two-dimensional nozzle is confined by the streamlines $\mathrm{\Psi }(0)$ and $\mathrm{\Psi }(ib)$, where $0<b<a$. The velocity field can be deduced from Equation (2222 $$\begin{array}{c}\hfill W(\mathit{\zeta })=u_{\infty }\mathit{\zeta }-i\mathit{\kappa }log\left(\frac{\mathit{\zeta }-ia}{\mathit{\zeta }+ia}\right),\mathrm{\Psi }(\mathit{\zeta })=\mathfrak{I}(W(\mathit{\zeta }))\end{array}$$22 ) as23 $$\begin{array}{c}\hfill u-iw=u_{\infty }-i\mathit{\kappa }\left(\frac{1}{\mathit{\zeta }-ia}-\frac{1}{\mathit{\zeta }+ia}\right)\end{array}$$23 while the pressure field is obtained from Bernoulli equation, once the total pressure level $P_{in}^o$ at inlet is fixed. The numerical test has been performed on a $90\times 30\times 20$ grid for the case with $b=0.7164$ , $u_{\infty }=0$, $\mathit{\kappa }=-1$. We also set $c=3$. With the scope of generating a fully 3D transient flow, we started from an initial geometry characterized by a wavy upper wall, as shown in figure 3(a). The theoretical values of total pressure $P_{in}^o=1$, the flow angles at inlet $\mathit{\alpha }=0$ and $\mathit{\beta }=0$, the static pressure distribution at the outlet $p_{ex}(y,z)$ and the target static pressure $p_{up}=p(\mathrm{\Psi }(ib))$ along the upper wall are imposed as boundary conditions. As final solution, the projection of the upper wall geometry on the $(x,y)$-plane is expected to match the theoretical curve $\mathrm{\Psi }(ib)$. Initial and final geometry are shown in Figure 3(a) and (b). The theoretical solution is obtained with high accuracy. As visible in Figure 3(c), the computed wall contour and the pressure distribution are almost superposed to the theoretical solution. The convergence history of the $L^2$-norm of the shape error and the maximum wall velocity are shown in Figure 3(d). Since Euler equations do not introduce dissipation or diffusion, the transient depends on how fast the propagating perturbations cancel out due to inertial effects and on how quickly they are absorbed by the flow-field boundaries. The full computation has required an half hour on a 8-core Intel i7 desktop computer.

Figure 4. Schematic representation of the free vortex flow and of the selected stream tube (a). Initial (b) and final (c) geometry of the boundary surface.

Figure 5. The free vortex flow. (a) A longitudinal and a cross-section of the pressure field obtained as solution of the inverse problem. (b) Convergence history of the $L^2$-norm of the shape error.

3.2 Free vortex flow

The free vortex flow is a special case of flow in radial equilibrium. [31] This axis-symmetric analytical solution of the Euler equations is used to validate the inverse procedure formulation (i) of Section 2.1.2. The test case is derived from the general solution of the flow field between two rotating cylinders, as depicted if Figure 4(a). Let us suppose that the inner cylinder rotates with angular velocity ${\mathrm{\Omega }}_1=k/R_1$, and the outer one with angular velocity ${\mathrm{\Omega }}_2=k/R_2$, where $R_1$, $R_2$ are the radii of the inner and outer cylinder, respectively, and $k$ is a given constant. The velocity field is defined as24 $$\begin{array}{c}\hfill c_x=u_o,c_{\mathit{\theta }}=k/r,c_z=0\end{array}$$24 where $u_o$ is also constant. Equation (2424 $$\begin{array}{c}\hfill c_x=u_o,c_{\mathit{\theta }}=k/r,c_z=0\end{array}$$24 ) represents an helicoidal motion with translation velocity $u_T=h\mathit{\omega }=u_o$, tangential velocity $v_T=\mathit{\omega }R$ and angular velocity $\mathit{\omega }=k/R^2$. By integration, the flow trajectories can be found as25 $$\begin{array}{c}\hfill X=X_o+u_{o}t,Y=Y_o+Rcos\mathit{\omega }t,Z=Z_o+Rsin\mathit{\omega }t\end{array}$$25 with $\mathit{\omega }t=\mathit{\vartheta }(t)-{\mathit{\vartheta }}_o$; in our case it is26 $$\begin{array}{c}\hfill {\displaystyle R=r,\mathit{\vartheta }-{\mathit{\vartheta }}_o=\frac{(X-X_o)k}{u_o{r}^2}}\end{array}$$26 In the numerical test, we tried to find the three-dimensional shape of the stream tube having a sector of circular crown as planar section at $x=x_0$ (see the grey zone in Figure 4(a). The analytical solution of the stream surface can be easily evaluated by using Equation 2626 $$\begin{array}{c}\hfill {\displaystyle R=r,\mathit{\vartheta }-{\mathit{\vartheta }}_o=\frac{(X-X_o)k}{u_o{r}^2}}\end{array}$$26 for each point of the stream tube contour at the axial station $x=x_0$. The test follows the footsteps of the previous case. An initial, arbitrary configuration is assumed and the boundary conditions are imposed according to the theoretical solution and to the inverse problem formulation. The total pressure $P^o=1$ and the theoretical flow angles are imposed at inlet. At the outlet, the static pressure is given by27 $$\begin{array}{c}\hfill {\displaystyle p_{ex}(r)=P^o-\frac{1}{2}\left(u_o^2+\frac{k}{r}\right)}\end{array}$$27 The domain boundaries at the outer and inner radii are treated as solid walls, while the lateral boundaries are the unknown surfaces we aim to resolve. As stream surfaces, they are impermeable, since the velocity vector must be tangent to them. Moreover, they cannot sustain any force, so that the pressure jump across these interfaces is null.

We have recovered the lateral surfaces shape by imposing a zero pressure jump condition and by allowing these boundaries to move according to the dynamic model of Section 2.1.2 until a steady state is reached. In the final configuration, the flow field is expected to coincide with the free vortex solution and the two moving boundaries with the corresponding flow surface, according to Equation (2626 $$\begin{array}{c}\hfill {\displaystyle R=r,\mathit{\vartheta }-{\mathit{\vartheta }}_o=\frac{(X-X_o)k}{u_o{r}^2}}\end{array}$$26 ). The initial guess and the final domain shape are shown in Figure 4. At the steady state, the free vortex flow is obtained and it is presented in Figure5(a), where the iso-pressure contours on a longitudinal and a transversal section are shown. The $L_2$-norm of the shape error is given in Figure 5(b). The absolute error is of the order of ${10}^{-5}$.

Figure 6. Inverse problem on a stator blade. (a–c) Evolution of the computational domain and of the blade geometry. Iso-pressure contours are also shown. (d) Comparison of the target blade load with the corresponding distributions computed by the inverse procedure and by the direct solver.

3.3 Inverse problem on a stator blade

The numerical experiment here proposed is a check of the full correspondence, at the steady state, between the inverse procedure and the direct solver solution on the same geometry. The blade geometry computed by the inverse procedure is introduced in a direct solver and the flow field is computed again, now by using standard boundary conditions. The two solutions are then compared. Arbitrarily, we assumed the following expression for the blade thickness28 $$\begin{array}{c}\hfill t\left(\frac{x}{c_a}\right)=\left[b_0+k{c}_a\left(1-\mathit{\xi }\right)\right]\sqrt{\mathit{\xi }\left(1-\mathit{\xi }\right)},\mathit{\xi }=\frac{x-x_l}{C_a}\end{array}$$28 where $x_l$, $x_t$ are the leading and trailing edge axial stations, respectively, and $C_a=x_t-x_l$ is the axial chord. The constants are set to $b_0=0.04$ and $k=0.12$. Similarly, the blade load is prescribed as29 $$\begin{array}{c}\hfill \mathrm{\Delta }p\left(\frac{x}{c}\right)=\mathit{\gamma }[1-cos(2\mathit{\pi }\mathit{\xi })]\end{array}$$29 with $\mathit{\gamma }=0.675$. The inlet and outlet boundary conditions are: $P_{in}^o=1$, $\mathit{\alpha }=0$, $\mathit{\beta }=0$ and $p_{ex}=0.99$, respectively. Starting from the initial guess of a blade with planar camber line, the transient dynamics of the system is integrated numerically until a steady state is reached. At each step, the blade geometry is recovered from the computed velocity of the blade surface. Only the dynamics of the blade suction side is computed. Pressure-side surface is deduced from the knowledge of thickness. A snapshot sequence of the blade moving surface is depicted in Figure 6(a–c). The transient resembles the motion of a flag on the wind, fastened on one side. More precisely, the grid describes a sliding motion on the tangential direction, since the axial chord of the blade does not change during the inverse process. At convergence, the blade geometry is found. By switching the code in direct mode, the flow field is computed with standard boundary conditions, and the blade load is evaluated. The solutions of the inverse and direct solver computations are in good agreement. The computed blade load $\mathrm{\Delta }p$, in direct and inverse modes and the target ${(\mathrm{\Delta }p)}_t$, defined in Equation (2929 $$\begin{array}{c}\hfill \mathrm{\Delta }p\left(\frac{x}{c}\right)=\mathit{\gamma }[1-cos(2\mathit{\pi }\mathit{\xi })]\end{array}$$29 ), are compared in Figure 6(d). A pressure mismatch, with respect to the target blade load distribution of Equation (2828 $$\begin{array}{c}\hfill t\left(\frac{x}{c_a}\right)=\left[b_0+k{c}_a\left(1-\mathit{\xi }\right)\right]\sqrt{\mathit{\xi }\left(1-\mathit{\xi }\right)},\mathit{\xi }=\frac{x-x_l}{C_a}\end{array}$$28 ), can be observed close to the blade leading edge. This is a well-known local effect caused by the poor resolution properties of the H-type grids close to the leading edge. A finer grid gives a more detailed flow but is meaningless for the test proposed, since the inverse and the direct solvers are acting on the same grid.

Table 1. Main parameters of E/TU-1 test-case.

Parameter	Value
Tip diameter	0.71 m
Hub diameter	0.566 m
Chord length	0.12 m
Axial chord length	0.087 m
Aspect ratio	0.6
Number of blades	21
Pitch to chord ratio (hub;midspan;tip)	0.705; 0.795; 0.885

Figure 7. The E/TU-1 Test Case. Initial (a) and final (b) blade configuration, as deduced at the steady state by the inverse procedure. The computational grid used and the pressure field are also visible.

3.4 E/TU-1 test case

The E/TU-1 test-case [32] focuses on a low speed, low aspect ratio, high turning annular nozzle guide vane (Table 1). The blades have a constant profile over the blade height and are untwisted. The severe inlet skew (e.g. ${51}^{\circ }$ at the hub) and the high turning angle (about ${65}^{\circ }$) generate a fully three-dimensional flow. We aim to compare a blade geometry computed by the inverse procedure with that given in Ref. [32]. Again, in solving the inverse problem, the formulation (i) is used. The blade thickness distribution was taken from [32], while the blade load has been deduced from the numerical results of Oktay et al. ([33], Figure A1).

In the initial geometry the effective thickness distribution has been superposed to a planar camber surface, as shown in Figure 7(a). A grid refinement study showed that a $90\times 32\times 20$ stretched grid was appropriate for describing the evolution of the computational domain.

Figure 8. The E/TU-1 Test Case. Distribution of the static pressure coefficient $C{p}_s$ along the blade at various radial positions.

Figure 9. The E/TU-1 Test Case. (a) Convergence history of the $L^2$-norm of the shape error and (b) evolution of the computed blade profile at the midspan (the real geometry is shown as ‘target’) and (c) computed blade profiles at different radii.

The numerical test follows the footsteps of the previous cases. The initial guess and the solution of the inverse problem are shown in Figure 7. The procedure has converged at slower rate than previous examples. The flow field around the thicker blade with rounded leading and trailing edges is resolved by the H-type grid with lower accuracy. It must be noted that also the target pressure distribution itself was extracted from a numerical simulation affected by the same accuracy problem.

The inverse procedure has required more grid points than Ref. [33] in the stream-wise direction ($90$ nodes instead of $60$) to ensure an adequate resolution of the blade contour motion. At each time step the mesh was recomputed from scratch, by maintaining the same stretching constants and domain discretization. As far as Euler equations are concerned, the computations remain quite affordable, even by using this rough remeshing approach.

In Figure 8, the final pressure distribution obtained along the blade is compared with that extracted from Ref. [33]. Although the blade load was enforced in the numerical test, the local pressure distribution is also in fair agreement with the reference data-set. The higher mismatch is observed at the tip radius (Figure 8). In that figure, it is also evident that both computations suffer the same discretization problem (a H-type grid is also used in Ref. [33]). Finally, the convergence history and the evolution of the blade geometry are shown in Figure 9.

4 Conclusions

A numerical method of solving three-dimensional inverse problems in incompressible flows has been presented. Details have been given for the derivation of the inverse problem over configurations of interest in the aerodynamic design as nozzles, diffusers, blades or wings. The numerical method is based on the pseudo-compressibility technique of preconditioning the unsteady Euler equations. Time accuracy is unnecessary, so that dual time stepping and the computation of the time derivatives of the metrics are avoided. Moreover, the characteristic-based formulation of the problem simplifies the parallelization of the code. The accuracy of the numerical procedure has been tested through comparisons with theoretical solutions, numerical computations performed with direct solvers and finally, test against experimental data. Discrepancies at the blade leading edge have been observed, caused by the H-type grid we used. Further extension of the inverse procedure to the Navier-Stokes equations, in a way similar to Ref. [29], is planned. An intermediate step required is the introduction of efficient over-relaxing techniques, to increase the convergence ratio. For the inverse design of aerodynamic profiles, the surface tracking algorithm should be modified to deal with C-type grid generation systems that allow for an accurate evaluation of the viscous stresses.

Appendix 1

Riemann invariants

In order to highlight the wave structure of the system, we switch to the quasi-linear form of equations (11 $$\begin{array}{c}\hfill \frac{\mathit{\partial }\mathit{U}}{\mathit{\partial }t}+\mathbf{\Theta }\left(\frac{\mathit{\partial }\mathit{F}}{\mathit{\partial }x}+\frac{\mathit{\partial }\mathit{G}}{\mathit{\partial }y}+\frac{\mathit{\partial }\mathit{H}}{\mathit{\partial }z}\right)=0\end{array}$$1 )A1 $$\begin{array}{c}\hfill \frac{\mathit{\partial }\mathit{U}}{\mathit{\partial }t}+{\mathit{F}}_{\mathit{U}}\frac{\mathit{\partial }\mathit{U}}{\mathit{\partial }x}+{\mathit{G}}_{\mathit{U}}\frac{\mathit{\partial }\mathit{U}}{\mathit{\partial }y}+{\mathit{H}}_{\mathit{U}}\frac{\mathit{\partial }\mathit{U}}{\mathit{\partial }z}=0\end{array}$$A1 whereA2 $$\begin{array}{c}\hfill {\mathit{F}}_{\mathit{U}}=\left[\begin{array}{cccc}0& c^2& 0& 0\\ 1& 2u& 0& 0\\ 0& v& u& 0\\ 0& w& 0& u\end{array}\right]{\mathit{G}}_{\mathit{U}}=\left[\begin{array}{cccc}0& 0& c^2& 0\\ 0& v& u& 0\\ 1& 0& 2v& 0\\ 0& 0& w& v\end{array}\right]{\mathit{H}}_{{}_{\mathit{U}}}=\left[\begin{array}{cccc}0& 0& 0& c^2\\ 0& w& 0& u\\ 0& 0& w& v\\ 1& 0& 0& 2w\end{array}\right]\end{array}$$A2 After integration of Equation (A1A1 $$\begin{array}{c}\hfill \frac{\mathit{\partial }\mathit{U}}{\mathit{\partial }t}+{\mathit{F}}_{\mathit{U}}\frac{\mathit{\partial }\mathit{U}}{\mathit{\partial }x}+{\mathit{G}}_{\mathit{U}}\frac{\mathit{\partial }\mathit{U}}{\mathit{\partial }y}+{\mathit{H}}_{\mathit{U}}\frac{\mathit{\partial }\mathit{U}}{\mathit{\partial }z}=0\end{array}$$A1 ), according to the Gauss formula, the fluxes are written in the formA3 $$\begin{array}{c}\hfill \mathit{C}={\mathit{F}}_{\mathit{U}}{n}_x+{\mathit{G}}_{\mathit{U}}{n}_y+{\mathit{H}}_{\mathit{U}}{n}_z\end{array}$$A3 where $\mathit{n}$ is the outward normal to the cell surface. As a property of hyperbolic systems, matrix $\mathit{C}$ is decomposed asA4 $$\begin{array}{c}\hfill \mathit{C}=\mathit{R}\mathit{D}\mathit{L}\end{array}$$A4 where $\mathit{R}$ and $\mathit{L}$ are the matrices of the right and left eigenvectors, respectively, and $\mathit{D}$ is a diagonal matrix that contains the eigenvalues of $\mathit{C}$A5 $$\begin{array}{c}\hfill \mathit{D}=\left[\begin{array}{cccc}\tilde{u}& 0& 0& 0\\ 0& \tilde{u}& 0& 0\\ 0& 0& \tilde{u}+\tilde{a}& 0\\ 0& 0& 0& \tilde{u}-\tilde{a}\end{array}\right]\end{array}$$A5 where $\tilde{u}$, $\tilde{v}$ and $\tilde{w}$ are the velocity components in the local frame of reference, and $\tilde{a}=\sqrt{{\tilde{u}}^2+c^2}$ is the pseudo-velocity of sound. The matrix $\mathit{R}$ isA6 $$\begin{array}{c}\hfill \mathit{R}=\left[\begin{array}{cccc}0& 0& c^2\tilde{a}& -c^2\tilde{a}\\ -n_z& -n_y& c^2{n}_x+u(\tilde{u}+\tilde{a})& c^2{n}_x+u(\tilde{u}-\tilde{a})\\ 0& n_x& c^2{n}_y+v(\tilde{u}+\tilde{a})& c^2{n}_y+v(\tilde{u}-\tilde{a})\\ n_x& 0& c^2{n}_z+w(\tilde{u}+\tilde{a})& c^2{n}_z+w(\tilde{u}-\tilde{a})\end{array}\right]\end{array}$$A6 and, by inversion ($\mathit{L}={\mathit{R}}^{-1}$ )A7 $$\begin{array}{c}\hfill {\displaystyle \mathit{L}=\left[\begin{array}{cccc}{\displaystyle \frac{\tilde{u}{n}_z-w}{{\tilde{a}}^2{n}_x}}& {\displaystyle -\frac{w\tilde{u}+n_z{c}^2}{{\tilde{a}}^2}}& {\displaystyle -\frac{(w\tilde{u}+n_z{c}^2)n_y}{{\tilde{a}}^2{n}_x}}& {\displaystyle \frac{(1-n_z^2)c^2+\tilde{u}(\tilde{u}-w{n}_z)}{{\tilde{a}}^2{n}_x}}\\ {\displaystyle \frac{\tilde{u}{n}_y-v}{{\tilde{a}}^2{n}_x}}& {\displaystyle -\frac{v\tilde{u}+n_y{c}^2}{{\tilde{a}}^2}}& {\displaystyle \frac{(1-n_y^2)c^2+\tilde{u}(\tilde{u}-v{n}_y)}{{\tilde{a}}^2{n}_x}}& {\displaystyle -\frac{(v\tilde{u}+n_y{c}^2)n_z}{{\tilde{a}}^2{n}_x}}\\ {\displaystyle -\frac{\tilde{u}-\tilde{a}}{2{\tilde{a}}^2{c}^2}}& {\displaystyle \frac{n_x}{2{\tilde{a}}^2}}& {\displaystyle \frac{n_y}{2{\tilde{a}}^2}}& {\displaystyle \frac{n_z}{2{\tilde{a}}^2}}\\ {\displaystyle -\frac{\tilde{u}+\tilde{a}}{2{\tilde{a}}^2{c}^2}}& {\displaystyle \frac{n_x}{2{\tilde{a}}^2}}& {\displaystyle \frac{n_y}{2{\tilde{a}}^2}}& {\displaystyle \frac{n_z}{2{\tilde{a}}^2}}\end{array}\right]}\end{array}$$A7 The Riemann invariants are obtained byA8 $$\begin{array}{c}\hfill \mathit{\partial }\mathit{W}=\mathit{L}\mathit{\partial }\mathit{U}={\left\{R_{2v},R_{2w},R_1,R_3\right\}}^T.\end{array}$$A8 After algebraic manipulations, the result isA9 $$\begin{array}{c}\hfill R_{2v}=\tilde{v},R_{2w}=\tilde{w},R_1=(\tilde{u}+\tilde{a})p-c^2\tilde{u},R_3=(\tilde{u}-\tilde{a})p-c^2\tilde{u}\end{array}$$A9

Figure A1. The Riemann Problem (RP). Flow variables (a) and wave pattern (c) at inner interfaces. Characteristics wave pattern of the half-Riemann problems at the inlet (b) and outlet (c).

Riemann problem for internal points

Analysing the wave system depicted in Figure A1, the following conditions hold across wave system I, between fields (a) and (c),A10 $$\begin{array}{c}\hfill R_{3c}=R_{3a},{\tilde{v}}_c={\tilde{v}}_a{\tilde{w}}_c={\tilde{w}}_a\end{array}$$A10 whereas, across wave system III, between fields (b) and (d), we haveA11 $$\begin{array}{c}\hfill R_{1d}=R_{1b},{\tilde{v}}_d={\tilde{v}}_b{\tilde{w}}_d={\tilde{w}}_b\end{array}$$A11 Wave system II represents a contact discontinuity, so thatA12 $$\begin{array}{c}\hfill {\tilde{u}}_c={\tilde{u}}_d,p_c=p_d\end{array}$$A12 the second relation holds since this surface cannot sustain any force. The solution of system (A10A10 $$\begin{array}{c}\hfill R_{3c}=R_{3a},{\tilde{v}}_c={\tilde{v}}_a{\tilde{w}}_c={\tilde{w}}_a\end{array}$$A10 )–(A12A12 $$\begin{array}{c}\hfill {\tilde{u}}_c={\tilde{u}}_d,p_c=p_d\end{array}$$A12 ) isA13 $$\begin{array}{cc}& p_c=p_d=\frac{R_{1b}-R_{3a}}{{(\tilde{u}+\tilde{a})}_b-{(\tilde{u}-\tilde{a})}_a}\hfill \end{array}$$A13 A14 $$\begin{array}{cc}& {\tilde{u}}_c={\tilde{u}}_d=\frac{{(\tilde{u}+\tilde{a})}_b{R}_{3a}-{(\tilde{u}-\tilde{a})}_a{R}_{1b}}{[{(\tilde{u}+\tilde{a})}_b-{(\tilde{u}-\tilde{a})}_a]c^2}\hfill \end{array}$$A14 A15 $$\begin{array}{cc}& {\tilde{v}}_c={\tilde{v}}_a\frac{1+sign({\tilde{u}}_c)}{2}+{\tilde{v}}_b\frac{1-sign({\tilde{u}}_c)}{2}\hfill \end{array}$$A15 A16 $$\begin{array}{cc}& {\tilde{w}}_c={\tilde{w}}_a\frac{1+sign({\tilde{u}}_c)}{2}+{\tilde{w}}_b\frac{1-sign({\tilde{u}}_c)}{2}\hfill \end{array}$$A16 where the Riemann invariants of fields (a),(b),(c) and(d) are approximated byA17 $$\begin{array}{cc}\hfill {\displaystyle R_{3a}}& ={(\tilde{u}-\tilde{a})}_a{p}_a-c^2{\tilde{u}}_a,\hfill \\ \hfill {\displaystyle R_{1b}}& ={(\tilde{u}+\tilde{a})}_b{p}_b-c^2{\tilde{u}}_b\hfill \\ \hfill {\displaystyle R_{3c}}& ={(\tilde{u}-\tilde{a})}_a{p}_c-c^2{\tilde{u}}_c,\hfill \\ \hfill {\displaystyle R_{1d}}& ={(\tilde{u}+\tilde{a})}_b{p}_d-c^2{\tilde{u}}_d\hfill \end{array}$$A17

Boundary conditions

Inlet conditions   At the inlet, a half-Riemann problem is formulated by replacing the missing signals with an equal number of boundary conditions. Three wave signals are unavailable at inlet, two travelling with velocity $\tilde{u}$ and one travelling with velocity $\tilde{u}+\tilde{a}$ (see Figure A1(b)). The former wave conveys signals $\tilde{v}$ and $\tilde{w}$, here replaced by imposing the flow angles $\mathit{\alpha }={tan}^{-1}(v/u)$ and $\mathit{\beta }={tan}^{-1}(w/u)$. The signal conveyed along $\tilde{u}+\tilde{a}$ is replaced by the condition of constant inlet total pressure $P^o$.

Outlet Conditions   In subsonic flows, ${\mathit{\lambda }}_1=\tilde{u}-\tilde{a}$ is always negative. The signal that travels along the characteristic with slope $\tilde{u}-\tilde{a}$ is missing at the outlet (see Figure A1(d)) and is replaced by the boundary condition $p=p_{exit}$.

Wall Boundary Condition   In inviscid flows, the no-through-flow condition is imposed at wall. When the wall is fixed, the wave configuration is represented by the scheme shown in Figure A1(d), where wave II is now travelling with velocity $\tilde{u}=0$ at the boundary. The characteristic with slope $\tilde{u}-\tilde{a}$ is missing, so that one boundary condition is required, this is written as ${\tilde{u}}_c=0$. If the wall is moving with velocity $V_p$, the above-mentioned conditions hold in the relative motion.