Computational methods

7.11.1 Methods based on the momentum integral equation

In the general case with an external pressure gradient the momentum integral equation must be solved numerically. There are a number of ways in which this can be done. One method for laminar boundary layers is to use the approximate expressions (7.64a’,b’,c’) with Eqn (7.59) (or 7.98) rewritten as

d0/dx can be related to d<$/dx as follows:

dв d(6I) _ d6 dl _ d6 dl dA d<5

d. v dx dx dx dx dA d<$ dx

. dl dA = |, + SdAdi

It follows from Eqns (7.62d) and (7.64b’) respectively that

dA 262 dUc

=—- ^~ = 2A

d о v dx

d / П АЛ

dA“ “63 15 + 72У

So

(7.128)

where

. 4 ч r 2Л / l A 1 ^ _ “ 63 (l5 + 72

Thus Eqn (7.127) could be readily converted into an equation for ddi/dx by dividing both sides by E’i(A). The problem with doing this is that it follows from Eqn (7.64c’) that

Thus in cases where either 6 = 0 or t/e = 0 at x = 0, the initial value of Q, and therefore the right-hand side of Eqn (7.127) will be infinite. These two cases are in fact the two most common in practice. The former corresponds to sharp leading edges and the latter to blunt ones. To deal with the problem identified above both sides of Eqn (7.127) are multiplied by 2Ue6/v, whereby with use of Eqn (7.128) it becomes

where use of Eqns (7.64a’,b’,c’) gives

а л 3 , Л2 4 /37 Л A2 Е2(Л) -4 + т-—Л + — + —( — – — – —

To obtain the final form of Eqn (7.127) for computational purposes Fi(A)62dt/c/dx is added to both sides and then both sides are divided by Fi(A); thus

dZ _ Тг(Л) д______ 2__M

dx F (А) Fi(A) v

where the dependent variable has been changed from 6 to Z = 62Ue/v. In the usual case when Vs = 0 the right-hand side of Eqn (7.129) is purely a function of A. Note that

Since Ue is a prescribed function of x this allows both Л and 6 to be obtained from a value of Z. The other quantities of interest can be obtained from Eqn (7.64a’, b’, c’).

With the momentum-integral equation in the form (7.129) it is suitable for the direct application of standard methods for numerical integration of ordinary differ­ential equations.* It is recommended that the fourth-order Runge-Kutta method be used. It could be used with an adaptive stepsize control, the advantage of which is that small steps would be chosen in regions of rapid change, e. g. near the leading edge, while larger steps would be taken elsewhere.

In order to begin the calculation it is necessary to supply initial values for Z and A at x = 0, say. For a sharp leading edge 6 = 0 giving Z = A = 0atx = 0, whereas for a round leading edge Ue = 0 and Z = 0 but Л = 7.052, see Example 7.8 above, this value of A should be used to evaluate the right-hand side of (7.129) at x = 0.

Boundary-layer computations using the method described above have been carried out for the case of a circular cylinder of radius 1 m in air flowing at 20m/s (v = 1.5 x 10-5 m2/s). In Fig. 7.45, the computed values of A and momentum

Fig. 7.45

Complete FORTRAN subroutines for this are to be found in W. H. Press etal. (1992) loc. cit.

thickness are plotted against angle around the cylinder’s surface measured from the fore stagnation point. A fourth-order Runge-Kutta integration scheme was used with 200 fixed steps between x = 0 and x = 3. This gave acceptable accuracy. According to this approximate calculation the separation point, corresponding to Л = —12, occurs at 106.7°. This should be compared with the accurately computed value of 104.5° obtained using the differential form of the equations of motion. Thus it can be seen that the Pohlhausen method gives reasonably acceptable results in terms of a comparison with more accurate methods. In point of fact neither value given above for the separation point is close to the actual value found experimentally for a laminar boundary layer. Experimentally, separation is found to occur just ahead of the apex of the cylinder. The reason for the large discrepancy between theoretical and observed separation points is that the large wake substantially alters the flow outside the boundary layer. This main-stream flow, accordingly, departs markedly from the potential-flow solution assumed for the boundary-layer calcula­tions. Boundary-layer theory only predicts the separation point accurately in the case of streamlined bodies with relatively small wakes. Nevertheless, the circular cylinder is a good test case for checking the accuracy of boundary-layer computations.

Numerical solutions of the momentum integral equation can also be found using Thwaites’ method.[47] This method does not make use of the Pohlhausen approximate velocity profile Eqn (7.63) or Eqn (7.65). It is very simple to use and for some applications it is more accurate than the Pohlhausen method described above. A suitable FORTRAN program for Thwaites method is given by Cebeci and Bradshaw (1977).

Computational methods based on the momentum integral equation are also available for the turbulent boundary layer. In this case, one or more semi-empirical relationships are also required in addition to the momentum integral equation. For example, most methods make use of the formula for Q due to Ludwieg and Tillmann (1949), namely

/77 7) —10-268

Cf = 0.246 x Ю-°-678іЧ—J

A very good method of this type is due to Head (1958)7 This method is relatively simple to use but, nevertheless, performs better than many of the much more com­plex methods based on the differential equations of motion. A FORTRAN program based on Head’s method is also given by Cebeci and Bradshaw loc. cit.

In order to begin the computation of the turbulent boundary layer using Head’s method it is necessary to locate the transition point and to supply initial values of 6 and H = 8*16. In Section 7.7.7, it was shown that for the boundary layer on a flat plate Eqn (7.89) held at the transition point, i. e. there is not a discontinuous change in momentum thickness. This applies equally well to the more general case. So once the transition point is located the starting value for 6 in the turbulent boundary layer is given by the final value in the laminar part. However, since transition is assumed to occur instantaneously at a specific location along the surface, there will be a discontinuous change in velocity profile shape at the transition point. This implies a discontinuous change in the shape factor H. To a reasonable approximation

Нл = Hu — AH

where

ДЯ = 0.821 +0.114loglo(Re0t) Reet< 5 x 104 ДЯ= 1.357 Reei > 5 x 104