EWEA 2011 Scientific Proceedings
22
9
that is part of the above mentioned project. TheDG method [1, 2, 3, 15, 14] is a higly compact finite-element projection method that provides a practical framework for development of a high-order method using unstruc-tured grids. Furthermore, compared to Finite-Difference methods it has the advantage of not requiring smooth meshes.
2 Governing equations
The integral boundary layer equations (IBL) can be de-rived starting from the unsteady, two-dimensional bound-ary layer equations with a fundamental assumption which is that the only effect of the boundary layer and the wake is to displace the inviscid flow away from the phys-ical body to create an effective displacement thickness. This assumption is valid for most aerodynamic flows of interest. Integrating the following equation for n = 0, 1 with respect to the transverse direction (normal to the no-slip surface) over the boundary layer thickness, we can derive the n-th moment of boundary layer equation:
momentum equation × (n + 1)un (1)
−continuity equation × (un+1 e − un+1),
in which the free-stream velocity ue= ue(x, t)is
pre-sumably known from a potential-flow analysis. In the current study the 0-th (momentum integral) and 1-st (kinetic-energy integral) moments of momentum equations are used to obtain the global quantities such as displacement, momentum and energy thicknesses of the boundary layer. Furthermore, unsteady laminar and turbulent closure models are used as the third equation to close the system of equations. The equation set is written in the conservation form and the system is shown to be hyperbolic.
The integral boundary layer equations can be written
G(u) = ueθ ueδk ueCτ , (4) S(u) = Cf 2ue− (δ∗+ θ)∂ue∂x − δ ∗1 ue ∂ue ∂t CDue− 2δk ∂ue∂x − 2θue1∂ue∂t Cτue− Cτ∂ue∂x − 2 US 1 ue ∂ue ∂t , (5)
with, δ∗the displacement thickness, θ the momentum
thickness, δkthe kinetic energy thickness, C τthe shear
stress coefficient, Cfthe friction coefficient, CDviscous
diffusion coefficient and USthe non-dimensional slip
ve-locity. It is assumed that the velocity distribution at the edge of the boundary layer, ue, is prescribed.
2.1 Closure relations
The equation set given in equations (2) through (5) form an unclosed set of equations since there are more un-knowns then the number of equations. The system can be closed by modeling some of the unknowns in terms of other unknowns. These models, so-called closure re-lations, can be derived by using experimental data or analytical solutions of representative test cases under certain assumptions [25, 20, 6]. i.e. by defining the shape factor, H for the displacement thickness δ∗and
the shape factor, H∗for the kinetic energy thickness δk
as follows: H =δ ∗ θ, H ∗ =δ k θ. (6)
The closure models can be defined for kinetic energy shape factor, H∗(H), friction coefficient C
f(H, θ),
vis-cous diffusion coefficient CD(H, θ, H∗(H)), the slip
ve-locity US(H, H∗(H))and the shear stress coefficient
Cτ(H, H∗(H), US). For detailed closure set definitions
for laminar and turbulent flow conditions, please see ref-erences [13, 6, 21].
One should note that all the closure relations given in the literature mentioned above are derived for the In this study the two-dimensional, unsteady integral
boundary layer equations are solved numerically to-gether with a closure set for laminar and turbulent flows. A high-order discontinuous Galerkin method is used for the spatial discretization and a multi-stage Runge-Kutta scheme is employed for the time integration. Numerical results show good agreement with the literature up to a separation point for steady problems.
Keywords: integral boundary layer equations, closure models, discontinuous Galerkin method.
1 Introduction
U∞
inviscid region
viscous region
Figure 1:Viscous and inviscid regions of the flow field.
For the design of wind turbine blades an accurate pre-diction of the loading on the blade surface is necessary. With increasing blade sizes the flexibility of the blades increases and unsteady effects are increasingly impor-tant to model in the prediction methods. Currently avail-able aerodynamic prediction tools that use engineering models are fast but have the disadvantage of being inac-curate. On the other hand (commercially available) CFD tools can achieve desired accuracy though they are far more costly in computational time when compared to the engineering tools. In order to predict the loading on the wind turbine blades we do not need to resolve the flow field in detail but rather we are interested in certain in-tegral quantities of the boundary layer and the effects on the pressure distribution. These integral values
de-dimensions. For the global description of the boundary layer these integral values are obtained by integrating the boundary layer equations with respect to the direc-tion normal to the no-slip surface, z, over the boundary layer thickness that leads to the integral boundary layer equations.
The rotor aerodynamics simulation code under devel-opment at ECN in the RotorFlow project is a combination of a panel method flow solver for the unsteady, incom-pressible, inviscid external flow (outer region) and an integral boundary layer solver for the unsteady viscous flow near the blade surface (inner region) (see figure (1)). The strong interaction between these two flow regions in separated flows will be accounted for by a so-called viscous-inviscid interaction scheme. Within this study the numerical solution of the unsteady integral boundary layer equations is presented.
Integral boundary layer equations have been used widely for the global description of the flow [23, 18, 24] especially in engineering applications for aircraft aerody-namics. The methods based on these equations have the advantage of lowering the space dimension by one, and as a consequence there is no need for a volume grid. This leads to a reduction in computational costs and input efforts. At the start of the CFD era these methods therefore were used as airfoil analysis tools and even today some successful applications are based on this approach [7]. Integral boundary layer equations have been commonly discretized with Finite-Difference schemes [19, 6]. More recently also Finite-Element and Finite-Volume discretizations are employed [12]. The in-tegral boundary layer equations have been analyzed in detail and a lot of research has been carried out on lami-nar and turbulent closure relations and lamilami-nar-turbulent transition models.
Currently, in wind turbine aerodynamic simulations we face a similar problem where we would like to perform
EWEA 2011 Scientific Proceedings
23
0
steady flow cases and unsteady effects are not consid-ered. Within this study, for unsteady equation set given in equation (2), the steady closure relation set is consid-ered.
An analysis of the eigenvalues of the equation set (2) together with the chosen closure laminar set shows that the eigenvalues of the characteristic equation are both real and positive for the smooth part of the flow. At the point of separation one of the characteristics become zero and then becomes negative downstream of this point where inverse flow occurs. At the point of separa-tion the system matrix (Jacobian matrix) become singu-lar which has already been noted by Goldstein [8]. The analysis has shown that at the point of separation the shape factor, H = 4.1908
3 Numerical Method
TheIBLequations given in equation (2) together with the closure set are solved with a high-order discontinuous Galerkin (DG) method. Equation (2) can be represented as follows:
L(u) = S(u) (7)
with, L the partial differential operator and u is the vector of unknowns. We consider a solution u(·, t) such that for each time t ∈ It, u(·, t) belongs to the function space U
of the form
u(·, t) ∈ U3, U ≡ L2(Ω), (8)
where L2(Ω)denotes a Hilbert space of all square
inte-grable functions on Ω with an associated inner product defined by:
f, g Ω≡
Ω
f (x)g(x)dΩ, f, g ∈ L2(Ω). (9)
The weak formulation of theIBLequation can now be written as
L(u(·, t)), v = s, v , ∀v ∈ U3. (10)
In order to discretise theIBLwe divide the solution do-main Ω into non-overlapping elements Ωjsuch that
¯ Ω = Ne j=1 ¯ Ωj, (11)
where Ωj= Ωj∪∂Ωjis the closure of Ωjand the
bound-ary ∂Ωjbelongs to at most two elements and Ne
de-notes the number of elements. We consider an approxi-mate solution uh(·, t)to the solution u(·, t) in the
follow-ing form
uh(·, t) ∈ U3
h, Uh= span{bjk} ⊂ U, (12)
where Uhis a finite-dimensional subspace of U. The
functions {bjk}are linearly independent basis functions
defined such that
bjk(x) = 0, x∈ Ω/ j, (13)
bjk≡
bjk(x), x∈ Ωj,
0, x∈ Ω/ j. (14)
The basis functions are continuous in Ωj and k =
0, 1, .., Mis the index of the polynomials where the upper limit is defined as;
M (p, d) =1 d! d l=1 (p + l), (15)
with d the number of space dimensions and p the highest degree of the polynomials used. We consider the local approximate solution uh(·, t), in Ωj, of the solution u(·, t)
as an expansion on to the local basis set {bjk}
uh|Ωj(x, t) = ujk(t)bjk(x), ujk(t) ∈ L2(It), (16)
where, ujkare the solution expansion coefficients or the
degrees of freedom for the solution on Ωjand functions
of time only in this semi-discrete approach.
We approximate the weak formulation (Eq. (10)) by re-placing the solution u(·, t) with the approximate solution uh(·, t)
L(uh(·, t)), vh = s, vh , ∀vh∈ Uh3. (17)
Since Eq. (17) holds for any function vh∈ Uh3we can
replace vhby bjm∈ Uh3to get
L(uh(·, t)), bjm = s, bjm , ∀j ∈ (1, 2, .., N e),
∀m ∈ (0, 1, .., M ). (18) Inserting Eq. (18) into Eq. (7) and integrating over Ω leads to Ωj L(uh(x, t))bjmdΩ = Ωj sbjmdΩ. (19)
Integration by parts and applying Gauss’ theorem to the third term gives,
Ωj ∂ujk ∂t bjkbjmdΩ − Ωj fji∂bjm ∂xidΩ + ∂Ωj bjmfjinjidS = Ωj sbjmdΩ. (20)
where njidenotes the unit outward normal vector on
∂Ωj. For each j, Eq. (20) contains only the unknowns
ujk(t), k = 0, .., M, giving rise to Nesets of (M + 1)
or-dinary differential equations for the functions ujk(t). The
coupling of the functions ujk(t)in neighboring elements
is achieved by replacing the (normal component of the) flux fjiin the surface integral term by a numerical flux
fjinji|∂Ωj= h(¯uj, ¯ul, nj) (21) where ”l” is the index of the element adjacent to ele-ment ”j” and ¯uj(x, t) ≡ ujk(t)¯bjk(x). With Eq. (16) and
Eq.(21), Eq. (20) is recast as: Ωj ∂ujk ∂t bjkbjmdΩ − Ωj fji∂bjm ∂xidΩ + ∂Ωj bjmhdS = Ωj sbjmdΩ, ∀m, j. (22)
Analysis of the discontinuous Galerkin method proved a rate of convergence of at least hpfor general
triangu-lations and of hp+1for Cartesian grids employing basis
polynomials up to order p, where h is a length scale that represents the size of elements [11, 10, 16, 17].
At any interface between two elements, since the so-lution is allowed to be discontinuous, there is a left state and a right state leading to a Riemann problem which is represented by the flux vector as shown in equation (21). Solving the Riemann problem will provide the cou-pling and handle the discontinuity at element interfaces. Various kinds of flux formulas have been proposed and used in the literature to approximate the Riemann prob-lem. In this study we will consider the Lax-Friedrich flux formula of the form:
h(¯uj, ¯ul, nj) = 1
2{f (¯uj) + f (¯ul) −θ|a|max(¯ul− ¯uj)}, θ ≥ 0, (23)
where, |a|maxis the maximum (absolute value) of the
eigenvalues of the system matrix.
For the linear problems the terms of the equation (22) can be integrated exactly once for all. For the nonlinear problem considered here the boundary integral term and the volume integral term in Eq. (22) can be evaluated numerically by applying numerical quadrature formulas of the required order [4].
The basis functions are defined on each element ˆΩ, in the computational space. The set {bk}is complete in the
sense that it spans ˆPp( ˆΩ), the space of all polynomials
on ˆΩ with real coefficients and with a degree ≤ p: ˆ
Pp( ˆΩ) =span{b0, b1, . . . , bM}. (24)
In this study, the second kind of Chebyshev polynomials are used as basis functions.
The compactness of the method results in an easy im-plementation of the boundary conditions. The boundary conditions can be implemented by prescribing the exact external solution or by reformulating the Lax-Friedrichs flux in terms of the interior solution and the physical boundary conditions.
The time integration is performed by an explicit multi-stage Kutta algorithm [5]. Multi-multi-stage Runge-Kutta methods can be used to have a higher order time integration for every time step. High order methods are used to minimize the temporal errors and are required to have a stable explicit method for reasonable time steps. Stability requirement of the time integration is given by the CFL condition and additionally is restricted by the de-gree of the approximation polynomial used as the basis function.
4 Numerical results
4.1 Linear model problem
First of all theDGmethod is applied to a problem in which the analytical solution is known, in order to verify and in-vestigate the accuracy of the method. As a model prob-lem the following form of linear transport equation is con-sidered:
∂u ∂t+ a
∂u
∂x= 0, (25)
where a is a constant and u is the unknown variable. Periodic boundary conditions are assumed and an initial
EWEA 2011 Scientific Proceedings
23
1
Figure 4: Comparison of the L2error norms, for the
explicit RK-DG method for given polynomial degrees of the basis functions. Linear transport equation with initial condition u0(x) =sin(2πxL)and Courant number, C =
0.1.
shown in (Figure 5). Looking at (Figure 4) and (Figure 5) together, it can be seen that the method for p = 3 for theDOF=48 runs about 44= 256times faster then the
method for p = 0 forDOF=768 and gives about 103times
less error.
4.2 Flow over a flat plate
TheIBLequations are solved for the flow over a flat plate in which the velocity at the edge of the boundary layer is prescribed as ue= U∞. The exact solution for the
lam-inar flow over a flat plate is given by the solution of the Blasius equation [24]. Comparing the results of the nu-merical simulation to the mentioned exact solution would not be fair since theIBLequations include the closure re-lations which already introduces some error compared to the exact solution. An alternative is to find an exact solution to theIBLequations including the closure rela-tions for the given edge velocity profile. For the details of
Figure 5: Comparison of the CP U times per period,
for the explicit RK-DG method for given polynomial de-grees of the basis functions. Linear transport equation with initial condition u0(x) =sin(2πxL)and Courant
num-ber, C = 0.1.
such a solution please see Van den Boogaard [21]. In (Figure 6) a comparison of the momentum thick-ness of the numerical simulation with the exact solution is shown for a laminar flow over a flat plate for Re = 105.
A comparison of the numerical simulation with the
0 0.0005 0.001 0.0015 0.002 0.0025 0 0.2 0.4 0.6 0.8 1 θ x analytic DG
Figure 6:Comparison of the numerical result with the
analytical solution for the momentum thickness for a lam-inar flow over a flat plate for Re = 105.
ature for a turbulent flow over a flat plate is shown in (Figure 7) for the displacement thickness. The Reynolds number for the turbulent flow is Re = 107.
In (Figure 2) a comparison of the numerical solution of the model problem given in equation (25) with theDG method to the analytical solution is presented for a con-tinuous initial condition u0(x) = sin(2πx/L)for various
order of accuracies. The number of elements, N, used in the spatial domain is taken to be 10 and the discontinuity between the element interfaces are shown. When the order of the approximation polynomials is taken as p = 0 (first-order of accuracy in space) the jump in the discon-tinuous solution between each neighboring element is clearly larger leading to a high dissipation as expected. Please note that the dissipation could of course be min-imized for p = 0 by using more number of elements in spatial direction but in (Figure 2) the advantage of using high-order approximation is stressed. A more fair
com--1 -0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 u x/L analytical p=0 -1 -0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 u x/L analytical p=1 -1 -0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 u x/L analytical p=2 -1 -0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 u x/L analytical p=3
Figure 2:Comparison of the numerical simulation with
the analytical solution for various polynomial degrees for a constant number of elements, N=10. Simulation time of 1 period, u0(x) =sin(2πxL)and C = 0.1
.
parison between the order of accuracies can be made by
related with the . The spatial order of accuracy is
-1 -0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 u x/L analytical p=0 -1 -0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 u x/L analytical p=1 -1 -0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 u x/L analytical p=2 -1 -0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 u x/L analytical p=3
Figure 3:Comparison of the numerical simulation with
the analytical solution for various polynomial degrees for a constant number of degrees-of-freedom, DOF=24. Simulation time of 1 period, u0(x) =sin(2πxL)and C =
0.1
.
calculated by employing the L2-norm of the form:
L2= 1 L L 0 [u(x, t) − uexact(x, t)]2dx 1/2 . (29)
In (Figure 4) the L2-norm of the error as a function of DOFis shown. It is known that for aDGmethod employ-ing basis functions up to order p, the rate of convergence is hp+1
2in general (e.g.[9]).Here it is shown that the
cur-rent method is converging at a rate of hp+1for p = 1, 2
and 3 and with a rate slightly higher than hp+1 2for p = 0.
It is remarkable that the line of the order p method is sit-uated above the one for the order-(p − 1) method for any pconsidered. It is also observed that in case less points are used to evaluate the L2-norm, e.g. common grid
EWEA 2011 Scientific Proceedings
23
2
0 0.0005 0.001 0.0015 0.002 0.0025 0 0.2 0.4 0.6 0.8 1 δ∗ x DG 1/7 power law Green SekarFigure 7:Comparison of the numerical result with the
literature for the displacement thickness for a turbulent flow over a flat plate for Re = 107.
4.3 Flow over an airfoil
Another test case forIBLequations considered is the laminar and turbulent flows over NACA airfoils. The air-foils chosen are 0009 for laminar flow and NACA-0012 for turbulent flow. In both laminar and turbulent flow cases the boundary layer edge velocity distributions are extracted from a simulation performed by XFOIL [7] for an angle of attack of 0 degrees. and results are also compared to the XFOIL results. The unsteady compo-nent of the solution is used to converge to the steady problem.
Laminar flow over a NACA-0009 profile is performed for a reynolds number, Re = 105. In theDGsimulations
only 5 elements (N=5) are used in the solution domain and the maximum degree of the polynomial basis func-tions is set to 3 leading a 4th-order accurate method in space. A 4-stage Runge-Kutta method is used for the time integration.
In (Figure 8) and (Figure 9) comparisons of the numer-ical simulation with the XFOIL results are shown for dis-placement thickness and momentum thickness respec-tively. The numerical simulation is in very good agree-ment with the XFOIL results up to the point of separa-tion. One should note that both the presented numeri-cal method and XFOIL use the same set of closure re-lations. There is a discrepancy between two solutions after the point of separation where the shape factor, H, is about 4. As discussed earlier, for the steady
prob-lems there is a singularity in the solution at the sepa-ration point where one of the eigenvalues of the equa-tion system becomes negative. As menequa-tioned before, within this study the boundary layer edge velocity pro-file is prescribed and the integral boundary layer solver is not coupled with the potential solver with an viscous-inviscid interaction scheme. Because of the lack of this coupling theIBLequations are not able to go beyond the separation point. Results for a turbulent flow over
Figure 8: Comparison of displacement thickness for
laminar flow over NACA-0009 airfoil. Re = 105.
Figure 9:Comparison of momentum thickness for
lam-inar flow over NACA-0009 airfoil. Re = 105.
a NACA-0012 airfoil are shown in figures (10 and 11) for displacement thickness and momentum thickness re-spectively for a Reynolds number of, Re = 5 · 105. Also
for the turbulent flow current implementation of theIBL equations without viscous-inviscid coupling is not able to predict the solution correctly after point of separation.
0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02 0 0.2 0.4 0.6 0.8 1 δ∗ x Xfoil DG
Figure 10:Comparison of displacement thickness for
turbulent flow over NACA-0012 airfoil. Re = 5 · 105.
0 0.0005 0.001 0.0015 0.002 0.0025 0 0.2 0.4 0.6 0.8 1 θ x Xfoil DG
Figure 11:Comparison of momentum thickness for
tur-bulent flow over NACA-0012 airfoil. Re = 5 · 105.
4.4 Laminar flow over a cylinder
The laminar flow in the boundary layer of a cylinder started impulsively from rest is simulated with the follow-ing edge velocity profile:
ue(x) = V∞2sin(x). (30)
A simulation is performed for the flow over a cylinder with unit radius, for which the upper half has a length L = π.
Kinematic viscosity is taken unity as well, ν = 1. This unphysical exercise is suitable to demonstrate the un-steady behavior of the boundary layer since the simula-tion will not be able to reach a steady state condisimula-tion due to the singularity at the separation point. In (Figure 12) a comparison is shown between the numerical simula-tion and the results obtained by by van Dommelen and Shen [22] for the displacement thickness. In the men-tioned study, van Dommelen and Shen have used a field method for the boundary layer equations and Lagrangian coordinates to find an accurate solution. Please note that the current set of closure relations are not designed to capture large separation as present in the current ex-ercise which leads to a larger difference when compared to the literature especially after the point of separation.
Figure 12:Comparison of displacement thickness
be-tween numerical simulation and the literature for polyno-mial order p = 3.
5 Conclusions
In this study, the numerical solution of two-dimensional, unsteady integral boundary layer equations (IBL) to-gether with a closure set is performed using discontin-uous Galerkin method for laminar and turbulent flows. The boundary layer edge velocity distribution is assumed to be prescribed. TheIBLequation set together with the closure model is shown to be hyperbolic and for the steady, laminar set of equations the point of separation is determined via characteristics analysis. The high-order
EWEA 2011 Scientific Proceedings
23
3
neous generation of the singularity in a separating laminar boundary layer. Journal of Computational Physics, 38:125–140, 1980.
[23] T. von K´arm´an. On laminar and turbulent friction. Technical Report TM-1092, NACA, 1946. [24] F. White. Viscous Fluid Flow. McGraw-Hill, 2nd
edition, 1991.
[25] D. Whitfield. Analytical description of the complete
two-dimensional turbulent boundary layer velocity profile. AIAA journal, 17:1145–1147, 1979.
unless a potential flow solver is coupled via a viscous-inviscid interaction scheme. Finally the unsteady simula-tion is carried out for a flow over a cylinder started impul-sively from rest to demonstrate the unsteady capability of the current implementation.
References
[1] H. Atkins and C. Shu. Quadrature-free
implemen-tation of discontinuous Galerkin method for hyper-bolic equations. AIAA journal, 36:775–782, 1998.
[2] F. Bassi and S. Rebay. High-Order Accurate
Dis-continous Finite Element Solution of the 2D Eu-ler Equations. Journal of Computational Physics,
138:251–285, 1997.
[3] B. Cockburn, S. Hou, and C. Shu. TVB
Runge-Kutta Local Projection Discontinuous Galerkin Fi-nite Element Method for Conservation Laws II: Gen-eral Framework. Mathematics of Computation,
52(186):411–435, 1989.
[4] B. Cockburn, G. Karniadakis, and C. Shu.
Discon-tinuous Galerkin Methods; theory, computation and applications, volume 11 of Lecture Notes in Compu-tational Science and Engineering. Springer-Verlag,
2000.
[5] B. Cockburn and C. Shu. Runge-Kutta
discontin-uous Galerkin methods for convection-dominated problems. Journal od Scientific Computing,
16(3):173–261, 2001.
[6] M. Drela. Two-dimensional transonic aerodynamics
design and analysis using the Euler equations. PhD
thesis, MIT, 1985.
[7] M. Drela. XFOIL: An analysis and design system for
low Reynolds number airfoils. Low Reynolds
Num-ber Aerodynamics. Springer-Verlag, 1989.
Discontinuous Galerkin Method for a scalar hy-perbolic equation. Mathematics of Computation,
46(176):1–26, 1986.
[11] P. Lesaint and P. Raviart. chapter On a finite el-ement method for solving the Neutron Transport equation. Mathematical Aspects of Finite Elements in Partial Differential Equations. Academic Press, 1974.
[12] B. Mugal. Integral methods for three-dimensional
boundary-layers. PhD thesis, MIT, 1998.
[13] B. Nishida. Fully simultaneous coupling of the full
potential equation and the integral boundary layer equations in three dimensions. PhD thesis, MIT,
1997.
[14] H. ¨Ozdemir. High-order discontiuous Galerkin method on hexahedral elements for aeroacoustics.
PhD thesis, University of Twente, 2006. [15] H. ¨Ozdemir, R. Hagmeijer, and H.
Hoeijmak-ers. Verification of the higher order
Discontin-uous Galerkin Method on hexahedral elements. Comptes Rendus Mecanique, 333:719–725, 2005.
[16] T. Peterson. A note on the convergence of
Dis-continuous Galerkin Method for a scalar hyperbolic equation. SIAM J. Numer. Analysis, 28:133–140,
1991.
[17] G. Richter. An optimal-order error estimate for the
Discontinuous Galerkin Method. Mathematics of Computations, 50:75–88, 1988.
[18] L. Rosenhead. Laminar boundary layers. Dover, New York, 1966.
[19] T. Swafford. Analytical approximation of two-dimensional separated turbulent boundary layer ve-locity profiles. AIAA Journal, 26:923–926, 1983.
[20] T. Swafford. Three-dimensional, time-dependent,
