## Abstract

This paper presents a new numerical method to solve the equations of the asymptotic theory of separated flows. A number of measures was taken to ensure fast convergence of the iteration procedure, which is employed to treat the nonlinear terms in the governing equations. Firstly, we selected carefully the set of variables for which the nonlinear finite difference equations were formulated. Secondly, a Newton–Raphson strategy was applied to these equations. Thirdly, the calculations were facilitated by utilizing linear approximation of the boundary-layer equations when calculating the corresponding Jacobi matrix.

The performance of the method is illustrated, using as an example, the problem of laminar two-dimensional boundary-layer separation in the flow of an incompressible fluid near a corner point of a rigid body contour. The solution of this problem is non-unique in a certain parameter range where two solution branches are possible.

## 1. Introduction

One of the success stories in theoretical fluid dynamics in the second half of the twentieth century was the development of the theory of interaction between the boundary layer and external inviscid flow, usually referred to as ‘viscous–inviscid interaction’. It was first published by Stewartson & Williams (1969) and Neiland (1969) in relation to boundary-layer separation in supersonic flow, and by Stewartson (1969) and Messiter (1970), who studied incompressible flow near the trailing edge of a flat plate. Then, the theory of viscous–inviscid interaction, also known as the ‘triple-deck theory’, was applied to a variety of fluid flows, and by the mid-1970s, it had become clear that the viscous–inviscid interaction played a key role in many fluid dynamic phenomena. In addition to different forms of separation, it was shown to govern upstream influence in the supersonic and hypersonic boundary layers, development of different modes of instabilities, bifurcation of the solution and possible hysteresis in separated flows.

The incipience of separation in the boundary layer near a corner point of a rigid body contour was one of the first problems to which the theory was applied. The set of asymptotic equations governing the flow behaviour near the corner was formulated by Stewartson (1970), who also was able to construct analytically the solution of the linearized version of the problem. It describes a preseparated flow regime when the angle of rotation of the tangent to the body contour at the corner point *θ* is sufficiently small to cause only weak perturbations in the boundary layer. A numerical solution to the nonlinear problem was first produced by Ruban (1976), who observed the formation of a local separation region in the boundary layer on convex and concave corners. Similar calculations were also performed by Smith and Merkin (1982). The problem was later re-examined by Korolev (1992), who discovered that the solution is non-unique for the flow over convex corners. It appeared that with the same body geometry, two different solutions are possible; one with a ‘short’ separation region forming behind the corner point, and the other with a ‘long’ separation region.

In this paper, we shall utilize the corner problem as a test example for our numerical method.

## 2. Governing equations

Let us consider the two-dimensional flow of an incompressible fluid past a rigid body whose surface has a corner point, and near this point the body contour may be thought as being constructed of two straight lines forming an angle *θ* with one another (see figure 1). We shall denote the fluid velocity in the oncoming free stream by *U*_{∞} and the pressure by *p*_{∞}, and the characteristic length-scale of the body as *L*. The density of the fluid and dynamic viscosity coefficient are constant all over the flow-field and are denoted as *ρ* and *μ*, respectively. We assume that the Reynolds numberis large, and the angle *θ*=*Re*^{−1/4}*θ*_{0}, with *θ*_{0} being an order one quantity.

Under these conditions, the flow near the corner is governed by the interaction between the boundary layer and inviscid external flow. The region of interaction is known to occupy an O(*Re*^{−3/8}) vicinity of the corner point (see figure 1). It has a three-tiered structure being composed of the viscous near-wall sublayer (region 1), the main part of the boundary layer (region 2) and an inviscid potential flow (region 3) situated outside the boundary layer.

Let and be dimensional Cartesian coordinates, with measured from the corner point along the tangent to the body contour before the corner and in the normal direction. The dimensional velocity components in these coordinates are denoted by , respectively, and the pressure by . Using these variables, the solution in the viscous sublayer may be represented in the form (see Stewartson 1970)(2.1)Here, *λ* is the non-dimensional shear stress on the body surface in the boundary layer immediately upstream of the interaction region. Function *f*(*x*) defines the shape of the body. With *α*=*θ*_{0}/*λ*^{1/2} denoting the wall deflection angle, it may be written asIn order to avoid the singularity at *x* = 0, we shall assume (as in previous studies) that the corner is slightly rounded. Thus,(2.2)In our calculations, parameter *r* was typically taken to be *r* = 0.1.

Substituting (2.1) into the Navier–Stokes equations and setting *Re*→∞ results in the boundary-layer equations(2.3)They have to be solved with the no-slip condition on the body surface(2.4)and the matching conditions with the solution in the main part of the boundary layer (region 2)(2.5)and in the boundary-layer upstream of the interaction region (Ruban 1976)(2.6)Here, *η*=*y*/(−*x*)^{1/3} and function *g*(*η*) may be calculated by taking into account that(2.7)Constant *α* in (2.2) and (2.7) is the only non-dimensional parameter governing the flow behaviour.

The Hilbert integral of the linear aerofoil theory(2.8)closes the formulation of the interaction problem in (2.3)–(2.6). It may be deduced by analysing the response of the inviscid flow in region 3 to the displacement effect of the boundary layer.

## 3. Numerical technique

We introduce a non-uniform meshand denote the values of function *A*(*x*) at points *x*_{i} by *A*_{i}. Considered together, they constitute the grid function {*A*_{i}}. As pressure is also independent of *y*, it may be represented by the grid function {*p*_{i}}, whose elements are defined as *p*_{i}=*p*(*x*_{i}).

Given {*A*_{i}}, equation (2.8) allows calculation of what we shall call the ‘inviscid pressure gradient’,On the other hand, for the same {*A*_{i}}, we can calculate the ‘viscous pressure gradient’,based on equations (2.3)–(2.6). The task is to find function {*A*_{i}} that satisfies the following implicit set of equations:In order to solve these equations Newton–Raphson linearization has been used, with an improved approximation written as *A*_{i}+δ*A*_{i}. The correction function {*δA*_{i}} is determined by solving the matrix equation

### (a) Calculation of the Hilbert integral

We express integral (2.8) in the form(3.1)and note that with (2.2) the second term on the right-hand side of (3.1) may be calculated analyticallyThen, the first integral in (3.1) was represented for each *x* = *x*_{i} in the formand the three integrals on the right-hand side of this equation were calculated analytically using the Taylor expansion for *A*(*s*), which was written asfor the first, second and third integrals, respectively. This reduces (3.1) toThis equation may be obviously expressed in the vector formHere, vector is composed of the values *β*(*x*_{i}) of the pressure gradient d*p*/d*x* at the mesh points {*x*_{i}}, and vector of the elements of the grid function {*A*_{i}}. The elementsof matrix may be calculated analytically by direct differentiation of the above expression for the pressure gradient. However, we found that this can be done more easily by making perturbations Δ*A*_{k} of *A*_{k} and calculating the corresponding perturbations of the pressure gradient. Owing to the linear relationship between these quantities, the values of Δ*A*_{k} can be chosen arbitrarily. In our calculations, we simply assumed Δ*A*_{k}=1.

### (b) Solution of the boundary-layer equations

The boundary-layer equations (2.3) were calculated marching from one position *x*_{i} to the next one *x*_{i+1}. At each cross-section of the boundary layer *x*_{i}, starting with *i*=3, we formulate a set of the operator equationsFor this purpose, we utilize the momentum equation from (2.3). Provided that *u*_{j,m}≥ 0, we have for *m*=2,…, *M*−2(3.2)whereNear the wall (i.e. for *m*=1), we write(3.3)To close this set of equations, we utilize the boundary condition (2.5) at the outer edge of the viscous sublayer, which yields(3.4)Here, and in (3.2) and (3.3), the *y*-component of the velocity vector is calculated using the continuity equation from (2.3). It is written in finite difference form as follows:This equation is solved for *v*_{i,m} and is used for progressively increasing *m*=2, …, *M* starting from the wall where, according to (2.4), *v*_{i,1}=0.

In order to ensure the stability and convergence of the numerical scheme, at each point (*x*_{i}, *y*_{m}) where *u*_{i,m}<0, we change (3.2) and (3.3) toHere, is the solution obtained in the course of the previous march through the flow-field, andThe main idea of the method is as follows. Given *A*_{i}, the above set of equations serves to determine the longitudinal velocity component *u*_{i,m}, *m*=2, …, *M*−1 and pressure gradient at position *x*_{i}. This then allows us to calculate the elements of the matrix(3.5)by repeating the calculations *N*−2 times. Each calculation should be performed with just one element *A*_{k} of the grid function {*A*_{i}}, perturbed by a small value Δ*A*_{k}. This method is significantly simpler in programming but almost as fast as the direct method described, for example, in Korolev *et al*. (2002).

We found that for calculating the matrix elements (3.5), it is convenient to perform Newtonian linearization analytically. For this purpose, we represent the fluid dynamic functions in the boundary layer as(3.6)Here, *u*′, *v*′, *p*′ and *A*′ are assumed small. Substituting (3.6) into (2.3) and (2.6), and neglecting squares of perturbations, we arrive at the linearized boundary-layer equationswhich were used to calculate (3.5).

## 4. Numerical results

The calculations were performed using a non-uniform mesh 251×151, and then repeated on a double‐sized mesh to check the accuracy of the results. A concentration of the mesh points was arranged near the corner where minimal step sizes were Δ*x*=0.005 and Δ*y*=0.02. In the separation region, we typically used Δ*x*=0.1, which was gradually reduced to Δ*x*=0.025 near the reattachment point.

The results of the calculations are summarized in figure 2, which shows the length of the separation region *l* as a function of the scaled angle *α*. We see that the solution is non-unique, that is, in a certain range of the variation of *α*, two solutions are possible for each value of *α*; one with a ‘short’ and another with a ‘long’ separation region.

This is illustrated by figure 3, where the streamline pattern is shown for the solution with the ‘long’ separation region at *α*=−5.4. For comparison, we also show in this figure (dashed line) the contour of the ‘short’ separation region given by the second solution.

In figure 4, a detailed comparison of the displacement function *A*(*x*), non-dimensional skin friction pressure gradient d*p*/d*x* and pressure *p*(*x*) is conducted for the two solutions. The solid lines show the solution with a ‘long’ separation region, and dashed lines the solution with a ‘short’ separation region.

Figures 5 and 6 show the flow behaviour for *α*=−4.67. At this value of the scaled angle *α*, there are also two solutions to the problem. One represents the flow, which remains attached to the corner surface. A detailed analysis of this solution has been performed elsewhere (e.g. Ruban 1976). The second solution, which is of interest here, contains a very large separation region. It is interesting to note that the behaviour of the pressure and skin friction (see figure 6) near the reattachment point is well in line with the theoretical predictions of Ruban (1980).

Finally, it should be noted that there exists a minimal value of *α* (see figure 2), which suggests that with increasing angle |*α*|, the flow on the corner should undergo an abrupt reconstruction when the local separation region suddenly bursts, with a large-scale separation (originating from the corner point) being formed instead.

## Acknowledgments

This research was partly supported by the Russian Foundation for Basic Research under grant 03-01-00918.

## Footnotes

One contribution of 19 to a Theme ‘New developments and applications in rapid fluid flows’.

- © 2005 The Royal Society