Inverse method for viscous flow design using stream-function coordinates

In the paper, a new inverse method for viscous 2D laminar flows is developed. The method is based on incompressible Navier–Stokes equations transformed to the stream-function coordinate system (von Mises coordinates). The flow design problem with appropriate boundary conditions is formulated and solved numerically. The geometrical shape of the boundary is obtained through the integration along streamlines. The method may be coupled with a flow analysis solver to model the influence of known parts of geometry. Results for two analytically solvable cases (the Poiseuille and the Jeffery–Hamel flows) are presented. Then, the foil design problem is considered as an example. Potential applications and developments towards axisymmetric and 3D flows are discussed.


Introduction
In many engineering problems, it is desired to design the flow system (e.g., the nozzle wall geometry or a blading system) which satisfies prescribed conditions (Fig. 1). Some of these conditions do not refer directly to geometric parameters, but to flow quantities (for example: the pressure distribution). Yet, since the geometry is not known at this stage, the design process is brought to the field of inverse methods. Most of such methods in fluid dynamics rely on some model simplifications (which affects accuracy) or iterative optimisation of direct problems (which is computationally expensive) [1]. Consequently, there is a need of inverse methods to circumvent those limitations.
With the growth of computer power, the optimisation techniques [1,2] have gained attention as less restrictive, general purpose methods. In parallel, inverse methods are still developed using such concepts as bladeto-blade or through-flow models [3,4]. Mixed direct-inverse techniques, like moving boundary/transpiration concepts, have been proposed even for viscous flows [5,6] in recent years. However, the moving boundary method needs to handle time-consuming remeshing.
The history of stream-function coordinates (SFC, also known as the von Mises coordinates) goes back to the boundary layer analysis in von Mises' work [7] in 1927. For direct (analysis) problems, the SFC method was used for potential and viscous flows in curved geometries [8][9][10][11] and even for multiphase flows [12]. The SFC technique has also been adapted to inverse methods of various 2D and 3D design problems for inviscid flows [8,10,[13][14][15]. Keller [16] presented a 3D inverse method with possible extensions to viscous flows. However, Scascighini et al. [17] proved that such an approach is applicable only to the so-called lamellar flows. This is due to the use of natural coordinates [18] which exist if and only if the velocity field satisfies the condition u · (∇ × u) = 0 (where u is the velocity vector). Yet, the assumption about the flow being lamellar is very restrictive. It is fully satisfied for 2D flows or irrotational 3D flows. Also, notable work was done for the approximate treatment of viscous effects in 2D and quasi-3D flows [17,19] using the integral boundary layer approach, suitable for inertia-dominated (high Reynolds number) flows.
A new, fully inverse method (without any solution of the direct problem) for 2D incompressible viscous flows, without explicit remeshing, is presented in the paper. The method is based on the Navier-Stokes equations and is suitable for laminar steady flows. In the direct problem (flow analysis), the system of governing equations is solved to obtain fields of flow variables (velocity, pressure, etc.) for a given shape of the boundary, initial and boundary conditions (e.g., no-slip wall, inlet, outlet, etc.). Yet, in inverse problems, at least a part of the geometry is unknown. One can replace some of geometrical variables by fluid flow variables. That will affect the boundary condition which must be changed, too. One of the possibilities is the replacement of one geometrical variable through the transformation of equations to the stream-function coordinates. After the transformation, one geometrical variable (here, y) is replaced by a flow variable (here, the stream-function ψ).
To the best knowledge of the authors, the present paper is a first attempt to use the idea of stream-function coordinates to the inverse problem for viscous flows with the exact treatment of diffusion terms.
The essence of the approach and details of the mathematical model are described in Sect. 2. Since the equation system after the transformation remains non-linear, a suitable numerical method is needed to solve it. The computational domain in x − ψ coordinates is rectangular and can be easily meshed with a regular grid. The space discretisation of the equations is done here by the finite difference method. A pseudo-time stepping is performed to obtain the steady-state solution. Details of the numerical scheme are described in Sect. 3. In Sect. 4, results and comparisons for a selection of design problems are presented. The numerical stability and accuracy are tested on several grids. Geometrical shapes for different boundary pressure distributions are determined. Next, some practical considerations are discussed in Sect. 5 on the basis of foil design. Concluding remarks are promising as far as the use of the method for various kinds of design problems is concerned, including the implementation of the present method to preliminary blade design process.

The flow model
Consider two-dimensional, incompressible viscous flow governed by the Navier-Stokes equations in a Cartesian coordinate system. The system has been arranged in a non-dimensional manner (with suitable scales): where p is the pressure and Re is the Reynolds number.

Coordinates transformation
Let us introduce the following coordinate transformation [10]: where x, y are Cartesian coordinates and ψ is the stream-function defined as so that where u x , u y are components of velocity vector u in Cartesian coordinates. Above, x corresponds to x direction in the new coordinate system. The spatial derivatives are transformed by using the chain rule Substituting (4) into (5) and taking into account that ∂ x /∂ x = 1 and ∂ x /∂ y = 0, the transformation takes the following form: Equation (6) represents the resulting transformation to a new coordinate system. The back transformation may be derived from the stream-function definition: Then, we can obtain the streamline shape through integration from inlet to outlet. However, the SFC method can only handle cases with no recirculation regions (no back-flow); there, u x = 0 at a certain point of the streamline, and the transformation in Eq. (5) will become singular as the result of Jacobian properties: So, the transformation is not directly invertible when u x = 0. Such a situation also appears at the walls due to the no-slip and impermeability conditions. The right hand side of Eq. (7) has the indeterminate form 0/0. To overcome this issue, we propose to use in the near-wall region the formula derived from Eq. (4) that may be integrated down to the wall (here ψ = ψ w ) from the nearest available streamline ψ = ψ w + ψ: where y w is the wall shape and y ψ w + ψ is the in-field streamline. Either, the above integral (9) does not behave well at the wall, since the velocity goes to zero due to no-slip boundary condition. So the function under the integral has a singularity on the boundary. This issue is resolved in Sect. 3. Additionally, the problem with recirculating flow may arise, where the multivaluedness of the streamfunction transformation appears. The solution process, based on pseudo-time marching (described in the next section), will diverge in that case, and no solution will be found. Also, the solver may detect existing separation during the calculation (by checking that velocity in the field is restricted to u x > 0). The possible appearance of a separation region is strictly related to the prescribed wall pressure distribution, mainly when an adverse pressure gradient is set. It is the user who should set such a boundary condition that will not produce the flow separation. Methods of choosing the pressure distribution were studied before [20,21]. This will be addressed in Sect. 5. Also, the problem may be solved using multiple regions of computation. However, this approach is far from trivial and as such is outside the scope of the present paper.

Transformation of the equation system
One may transform Eq. (1) to the stream-function coordinates using Eq. (6). This gives the following system: After the transformation, the convective term is simplified, but the diffusive part becomes much more complicated. However, in typical channel flow cases, one of the flow directions is dominant. The length scales satisfy L y L x , and velocity scales behave as U y U x , and then only diffusive terms u∂/∂ψ(u∂u/∂ψ) and u∂/∂ψ(u∂v/∂ψ) are dominant. This may be used for further simplifications. However, in the present paper, we will focus on the complete equation system (10).

Inverse problem
In the inverse problem, one wants to obtain the unknown shape y w (x) of the wall(s). First, we use Eq. (10) to solve the steady velocity and pressure fields in the (x, ψ) stream-function coordinates domain. Next, using Eqs. (7) and (9), we obtain the y shape of each streamline (including the wall).
As a closure of the equations system, standard boundary conditions may be applied except for the pressure at the wall streamline. The Dirichlet boundary conditions at the inlet and outlet are used for velocity and pressure. At the wall, the no-slip and impermeability conditions for velocity are set up.
Many solution approaches to the flow analysis problem use some boundary conditions for pressure at the wall, for example, of the Neumann type. Nonetheless, the analysis problem itself is uniquely defined without this boundary condition. The pressure distribution at the wall is then an explicit result of the method, directly connected to the wall shape. In the synthesis (design) problem, however, the shape is not known a priori and makes part of the solution. But, the shape handles the information how the velocity behaves near the wall boundary. In order to solve flow equations without a known boundary shape, something should be known about the pressure (similar to problem in [22]). This idea is readily understood if we consider the design of a nozzle using a very simple 1D flow model. The mean velocity u(x) and the sought cross-section area S(x) are linked by u(x)S(x) = const, whereas Bernoulli's principle connects them to pressure p(x) as follows: u 2 (x)/2+ p(x) = const. Consequently, the nozzle shape S(x) and velocity u(x) directly depend on the prescribed pressure distribution. For the presented inverse method, we assume that the Dirichlet pressure boundary conditions are known.

Numerical implementation
The SFC-transformed flow equations are solved by the time marching technique. As the divergence of the flow field is preserved only at the end of the solution process, this method is commonly known as the pseudo-time integration. Also, the artificial compressibility method was chosen to satisfy the continuity equation (with artificial compressibility coefficient β = 1): This brings us to a problem of three time-dependent equations. The equations are integrated numerically until a steady state is reached. The explicit Euler scheme was used for the test cases. The method is only first order accurate in time, but as the steady state is of interest, the pseudo-time integration scheme has no influence on the converged solution. For the solution of the flow fields in the pseudo-time, a stop criterion is chosen as a maximum acceptable change of a flow variable φ in pseudo-time, This criterion must be fulfilled for each of the fields u, v and p. We choose that each case should converge to at least t = 10 −5 .
The regular grid is used in the x direction of the (x, ψ) domain. Also, the uniform grid spacing is applied in ψ direction. The second order upwind scheme was used for the convective term in Eq. (10). The diffusive and pressure gradient parts in the momentum equation, as well as the continuity equation, were discretised with the second order central scheme.
Next, Eq. (7) is integrated by the Simpson rule. To integrate Eq. (9), one may expand u x near the wall into some chosen series u x (ψ) = i c i f i (ψ), with the desired accuracy, and integrate such a function analytically: The integral (9) has a singularity at the wall boundary and cannot be computed in a usual way. That is, we cannot use the expansion in the form that, applied to 1/u x (ψ), will lead to divergent numerical error. It will be increasing in the direction of ψ → ψ w (where u x (ψ w ) = 0). But one may note that in Cartesian coordinates the velocity profile along the wall may be expanded into a polynomial series: If we cut this expansion to and transform it to the stream-function coordinates, we obtain the following relation near the wall: On this basis, as firstly proposed by von Mises [7], the velocity field near the wall may be expanded in the Puiseux series: Concluding, the velocity field in the wall region is interpolated by series (18), using the velocity values on a streamline aligned to the wall, and then analytically integrated, cf. Eq. (9). So, assuming that ψ w = 0 (which implies c 0 = 0) and expanding velocity up to the second non-constant term u x (ψ) = c 1 √ ψ + c 2 ψ the formula (9) may be discretised: This discretisation was used in the presented inverse method. Also, please note that differential operators in stream-function coordinates are velocity dependent. Consequently, the discretisation affects the velocity, and as a result the discretisation itself. One may expect that the discretisation error will decrease better than the second order of accuracy.

Validation tests and results
As mentioned before, two test cases have been considered. For these cases, the boundary shape corresponding to the input velocity and pressure data is known exactly, so the design error can readily be measured. First, the method was validated by comparison with the analytical solution for the Poiseuille flow case. The convergence of the solution process is analysed. Next, the Jeffery-Hamel flow design problem-the linearly convergent nozzle-is considered. The difference between the numerical and exact design shape is analysed.

Poiseuille channel flow
A steady laminar flow through a straight channel was chosen as a first test case. Its analytical solution is well known, so one can easily estimate the error of the method. For channel flow, the steady Navier-Stokes equations reduce to the following viscous relation (no convective terms): resulting in a quadratic velocity distribution (in y direction) and a linear pressure drop related to viscous stresses. However, after the transformation to the stream-function coordinates, the velocity profile across the channel (in ψ direction) is given by As the transformation is velocity dependent, the velocity itself influences the stream-function coordinates distribution. Due to this, the profile (21) is steeper near the wall, where the velocity is lower (Fig. 2). That may introduce instabilities into the solution process. The respective analytical solution will be used as the inlet and the wall boundary conditions for the inverse design problem of the straight channel. The Poiseuille flow channel was designed following the methodology just presented. Velocity and pressure fields were obtained as the solution of (10). The velocity difference (between numerical and theoretical values) taken at x = 0.5 may be useful to analyse the error of the solution field. So, let us define the measure across the flow channel: Additionally, we may analyse the convergence dependent on grid density using the total error: Also, from the analytical solution, it is known that the channel walls should be parallel to the x direction. The difference between the known (theoretical) and numerically designed shape represents an error of the inverse method itself: The order of the discretisation error in stream-function coordinates is dependent on the velocity fields (as the differential operators are functionals of the velocity). As a result, the total error u,total converges at a non-constant rate (see Fig. 4). Moreover, Fig. 3 illustrates that the error u increases near the wall. The velocity u x near the wall boundary is very steep in ψ direction. Due to this, the error of discretisation of velocity is higher near the wall. This proves that the error is related to a singularity at the wall. As seen in Fig. 4, the total error y of the designed shape decreases linearly with increasing grid density. The maximal estimated error is satisfactorily small.
The pressure drop in Poiseuille flow (cf. Sect. 4.1) is connected to viscous stress. If the pressure change along the wall is different from the one in the Poiseuille flow (for a given inlet velocity profile), the convective term will participate in it. These possibilities bring us to a problem of curvilinear channel flow.

Jeffery-Hamel flow
The classical Jeffery-Hamel flow occurs between two divergent/convergent planes. The problem can be reduced to a non-linear ordinary differential equation, solvable by numerical methods (the analytical solution exists in the form of elliptic integrals). The pressure distribution of the Jeffery-Hamel analytical solution at the prescribed angle α was used as the Dirichlet boundary condition along the walls. The inlet and outlet conditions of the inverse problem were calculated from the analytical solution, too. As a result, the geometry of the convergent nozzle and flow fields were obtained (Fig. 5).
The difference between numerically designed and exact shape, Eq. (24), was calculated for various grid sizes, as shown in Fig. 6. The absolute value of the difference is higher than in Poiseuille flow, as more terms participate in the solution. However, the error is at an acceptable level, and the convergence rate is constant.
As another flow case, we also computed the curvilinear convergent and divergent nozzle design problems with, respectively, a linear pressure decrease or increase in the main flow direction. The resulting geometries were then input to a commercial flow solver and compared with our pressure distributions. The results (not shown) are very much similar to the case of Jeffery-Hamel flow.

Airfoil design-practical aspects
The convergent-divergent channel design problem is another important issue in turbomachinery, since this type of flows arises in the blade passages. The results presented here constitute a first step towards creating the 2D blade geometry in a flow channel (left for future work).
Two aspects need to be included into the consideration. First, the real inlet/outlet conditions may depend on parts of known geometry (like intake), but also may be influenced by the flow field in inverse design domain. So one may combine the flow analysis with geometry synthesis in a block-to-block manner, as shown in Fig. 7. On the other hand, for the purpose of inverse design of a convergent-divergent nozzle, one has to know a boundary pressure distribution at the walls. However, an additional constraint will appear for the blade channels that should have the same height at the inlet as well as at the outlet. Then, at the outlet, we have two boundary conditions (one for pressure p and one for geometry y) which makes the system overdetermined. In the blade design, this leads to the trailing edge closure problem, for which a pressure distribution is not known exactly at the designing stage. As shown by Mangler [20] and Lighthill [21], the pressure distribution is constrained by three integral relations to ensure a uniform free-stream velocity at infinity and to keep the foil closed at the trailing edge. Improper pressure distributions will lead to a blade profile with overly large or non-physical negative trailing edge thickness as shown in Fig. 8. Many authors addressed the problem in the past, and some  of them include an additional degree of freedom to the pressure boundary conditions (as we can see in [19]). A similar solution to that problem has been adopted here, as follows.
The requirement of an admissible (generally, positive) thickness at the trailing edge (in the meaning of blade design, cf. Fig. 7) is tantamount to saying that the flow domain as a whole has a convergent nozzle-like geometry. Otherwise, a negative thickness is equivalent to a divergent overall shape of the blading channel. In the inverse design process, one can note that a monotonically increasing pressure distribution will generally lead to divergent nozzle geometry, and a monotonically decreasing pressure distribution will result in a convergent nozzle or straight (as in Poiseuille flow) channel geometry. This brings the idea of prescribing the pressure distribution in the inverse method not only as a function of x direction but also as being dependent on the trailing edge thickness history τ te (t) during the iteration process. According to this idea, the trailing edge thickness τ te will be going to zero, as a result of adding a correction (a divergent/convergent update) to the originally assumed pressure distribution. So, at a given pseudo-time step, the pressure distribution will have the form: We propose now to prescribe the additional pressure update as a linear function of x multiplied by an unknown functional f [τ te ] related to the pseudo-time history of τ te and scaled by a suitably chosen constant c 0 : The functional f [τ te ] may be assumed as an integral over the pseudo-time: It is easy to see that if the iteration process with an updated pressure distribution leads to a zero trailing edge thickness, f will tend to be constant over the pseudo-time. That is because the thickness τ te = 0 will no longer affect the integral (27), and a steady solution will be obtained. So finally, the pressure distribution at the nozzle wall will be described as: The value of c 0 is chosen to assure and accelerate the convergence of the method. The constant was set to c 0 = 0.1 here. The final geometry of the convergent-divergent nozzle (the equivalent of a properly enclosed foil) designed in the above way is shown in Fig. 9. Generally, the foils need to have precise design in the leading edge region. However, there appears the limitation to avoid the singularity where u x = 0. This leads to a cusped leading edge. For practical purposes, the leading edge section of the foil will need to be replaced by a predesigned shape (e.g., a circular arc). The total CPU time overhead due to the correction procedure is negligible, which is important from the numerical point of view. Figure 10 shows the pressure distribution p w (x, τ te ) and its original set up p w (x) at the steady state. This shows that the final pressure distribution may be different from that initially prescribed, p w (x). This fact simply reflects the nature of inverse methods: for given flow parameters, the solution may not exist (as the system with prescribed p w (x) is overdetermined [sic!]).
The method has an additional advantage that the pressure distribution in the inverse problem may be assumed to be also a function of other variables like efficiency or dissipation coefficient. In that way, it will be straightforward to include the optimisation process within the inverse method pseudo-time stepping algorithm with only a minor increase of the design time.

Conclusions
The new inverse method for 2D viscous flow design was formulated and implemented. We explicitly excluded recirculating flow regions (normally, reflecting a bad design) and situations where the streamwise velocity component u x is zero (non-cusped leading edge). However, for the unavoidable treatment of walls (where, obviously, both velocity components are zero), we have proposed a technique to deal with the indeterminacy of backward transformation at the domain boundary as described in Sect. 2.1. To the best of the authors' knowledge, the presented method is the first inverse numerical solution of the Navier-Stokes equations in two dimensions in the SFC formulation with direct treatment of diffusive terms.
The simple pseudo-time stepping technique was adopted. The second order numerical scheme for space discretisation was implemented, however with a caveat of only linear convergence in certain cases (Sect. 3). Comparisons with analytically known cases of the Poiseuille channel and the Jeffery-Hamel flow served to Fig. 11 The stream-function ψ on λ surfaces in the three-dimensional model validate the approach. The total error of the method was kept at an acceptable level. Convergence of the solution process may be seen indirectly as a proof of existence of the inverse problem solution.
An extension of the method to a two-dimensional axisymmetric model is straightforward as well as an extension to compressible flow problems. In such a case, we need to change the stream-function formulation. For example, in compressible flow, the stream-function depends not only on the velocity field but also on the mass density ρ and its reference level ρ 0 :