Blasius boundary layer

A stationary, uniform flow passing over a vertical semi-infinite flat plate is analysed in this section. As it is shown in Fig. d.1.1, the fluid is flowing in the $ x $-direction at the constant velocity $ U_\inf $.

Fig. d.1.1. Blasius boundary layer along a semi-infinite flat plate.

The flow is considered to be an incompressible, homogeneous, Newtonian, 2-D steady-state flow, and the governing equations can be written as:
\begin{equation} \left\lbrace{ \begin{array}{l} \left(\textbf{u} \cdot \nabla\right) \textbf{u} - \nu\nabla^2\textbf{u} + \frac{1}{\rho} \nabla p = 0\\ \nabla \cdot \textbf{u} = 0 \end{array} }\right. \end{equation} The viscosity effects are important just inside a "thin" layer near the body wetted surface called boundary layer. The thickness of this layer is shown as $ \delta_x $ and:
\begin{equation} \frac{\delta _x}{x} \ll 1 \end{equation} should hold for every amount of $ x $, where $ x $ is the length from the beginning of the body.

To ensure about the presence of a boundary layer, it is necessary not to have so low Reynolds numbers and also to avoid flow separation.

Prandtl's 2-D boundary layer model

According to Prandtl's hypothesis and using continuity and momentum equation, it yields:
\begin{equation} \left\lbrace{ \begin{array}{l} U_y \sim U_x\frac{\delta_x}{x} \ll U_x \\ \ \\ \left(\textbf{u}\cdot \nabla \right)\textbf{u} \sim \nu \nabla^2 \textbf{u} \sim {U_x^2}/{x} \ \ \Longrightarrow \ \ \delta_x\sim \sqrt{\frac{\nu x}{U_x}}\\ \ \\ \frac{\partial p}{\partial y} \ negligible \end{array}}\right. \end{equation} Here the Prandtl's 2-D boundary layer model is presented:
\begin{equation} \left\lbrace{ \begin{array}{l} \partial_x u_x + \partial_y u_y = 0 \\ \ \\ u_x \partial_x u_x + u_y \partial_y u_x = \frac{1}{\rho}\partial_x p^{ext} + \nu\partial_{yy} u_x\\ \end{array}}\right. \end{equation} The $ y $-momentum equation is introduced inside the previous system through the pressure:
\begin{equation} \frac{\partial p}{\partial y} \ \text{negligible} \Longrightarrow p = p^{ext} \end{equation} and the $ x $-momentum equation is simplified in $ \nabla^2 u_x $:
\begin{equation} \delta_x \ll x \Longrightarrow \frac{\partial^2 u_x}{\partial x^2} \ll \frac{\partial^2 u_x}{\partial y^2} \end{equation}

Blasius boundary layer model

Now the Blasius boundary layer model can be introduced as follows:
\begin{equation} \left\lbrace{ \begin{array}{l} \partial_x u_x + \partial_y u_y = 0 \\ \ \\ u_x \partial_x u_x + u_y \partial_y u_x = \nu\partial_{yy} u_x\\ \end{array}}\right. \end{equation} Since it is a parabolic PDE, we can rewrite the system as:
\begin{equation} \frac{\partial [\textbf{u}(\textbf{x},t)]}{\partial t} + \textbf{L}[\textbf{u}(\textbf{x},t)] = \textbf{f}(\textbf{x},t) \end{equation} in which $ \textbf{L} $ is the elliptic operator, $ \textbf{u} $ is the solution of the problem, and $ \textbf{f} $ is the known function.

The solution admits just a characteristic line which defines the dependence domain. Along the parabolic direction $ t $, just one of the two edges needs a boundary condition and it represents the initial condition. However, each point of \textbf{x} edges needs a boundary condition. So the Blasius boundary layer can be illustrated as follows:
\begin{equation} \underbrace{u_x \partial_x u_x}_{\partial_t \textbf{u}(\textbf{x},t)} + \underbrace{u_y \partial_y u_x - \nu\partial_{yy} u_x}_{\textbf{L}[\textbf{u}(\textbf{x},t)]} = 0 \end{equation} with the initial condition:
\begin{equation} x = 0: u_x = U_\infty, \ u_y = 0 \end{equation} and boundary conditions:
\begin{equation} \left\lbrace{ \begin{array}{l} y = 0: u_x = 0, \ u_y = 0\\ \ \\ y \rightarrow \infty: u_x = U_\infty\\ \end{array}}\right. \end{equation}

Similarity solution for Blasius boundary layer

Assuming:
\begin{equation} \frac{u_x(x,y)}{U_\infty}=\phi\left[{\frac{y}{\delta_x(x)}}\right] = \phi[\eta(x)] \end{equation} and expressing the derivatives with respect to $ \eta $:
\begin{equation} \left\lbrace{\begin{array}{l} \frac{\partial}{\partial x} = \frac{d}{d \eta} \frac{\partial \eta}{\partial x} = -\frac{y}{\delta_x^2(x)}\frac{d\delta_x}{d x}\frac{d}{d \eta} = - \frac{\eta}{\delta_x(x)}\frac{d\delta_x}{d x}\frac{d}{d \eta} \\ \ \\ \frac{\partial}{\partial y} = \frac{d}{d \eta} \frac{\partial \eta}{\partial y} = \frac{1}{\delta_x(x)}\frac{d}{d \eta} \end{array}}\right. \end{equation} we have:
\begin{equation} \Longrightarrow \left\lbrace{\begin{array}{l} \partial_x = - \frac{\eta}{\delta_x(x)}\frac{d\delta_x}{d x} d_\eta \\ \ \\ \partial_y = \frac{1}{\delta_x(x)} d_\eta \end{array}}\right. \end{equation} Applying the derivatives on continuity equation gives:
\begin{equation} - \frac{\eta}{\delta_x(x)}\frac{d\delta_x}{d x}U_\infty \frac{d\phi}{d \eta} + \frac{1}{\delta_x(x)}\frac{d u_y}{d \eta} = 0 \end{equation} Then it is possible to calculate $ u_y $ using direct integration (integration by part):
\begin{equation} u_y = \frac{d\delta_x}{d x}U_\infty \int_0^\eta{\eta}\frac{d\phi}{d \eta}d\eta = \frac{d\delta_x}{d x}U_\infty \left[{\eta\phi -\int_0^\eta{\phi d\eta}}\right] \end{equation} Now we can introduce a new variable $ f $:
\begin{equation} f(\eta) = \int_0^\eta{\phi d\eta} \end{equation} and the velocity values can be written in terms of $ f $ as follows:
\begin{equation} \left\lbrace{\begin{array}{l} \frac{u_x}{U_\infty} = \phi(\eta) = \frac{df}{d\eta} \\ \ \\ \frac{u_y}{U_\infty} = \frac{d\delta_x}{d x}\left({\eta \frac{df}{d\eta} -f}\right) \end{array}}\right. \end{equation} Now it is possible to rewrite the momentum equation in terms of $ f $:
\begin{equation} U_\infty\left[ \underbrace{-\frac{df}{d\eta}\frac{\eta}{\delta_x}\frac{d\delta_x}{d x}\frac{d^2 f}{d \eta^2}}_{a} + \frac{d\delta_x}{d x}\left(\underbrace{\eta \frac{d f}{d \eta}}_{b} - f\right)\frac{1}{\delta_x}\frac{d^2 f}{d \eta^2}\right] = \frac{\nu}{\delta_x^2}\frac{d^2 f}{d \eta^2} \end{equation} in which terms $ a $ and $ b $ can be eliminated from the equation, since $ a = -b(d\delta_x/dx) $. Then we have:
\begin{equation} Re_{\delta} \frac{d \delta_x}{d x}f f^{(2)} + f^{(3)} = 0 \end{equation} in which $ Re_\delta = \frac{\rho U_\infty \delta_x(x)}{\mu} $. Similarity solution implies constant coefficients. Therefore, we have:
\begin{equation} Re_{\delta} \frac{d \delta_x}{d x} = k \Longrightarrow \int_0^x{\delta_x \frac{d \delta_x}{d x} dx} = \frac{1}{2}\left[{\delta_x^2(x) - \delta_x^2(0)}\right] = \frac{1}{2}\delta_x^2(x) = \frac{k\mu x}{\rho U_\infty} \end{equation} in which the arbitrary value of $ k $ can be chosen as 0.5, then we have:
\begin{equation} f f^{(2)} + 2 f^{(3)} = 0 \ \ , \ \ \text{where} \ \delta_x = \sqrt{\frac{\nu x}{U_\infty}} \end{equation} Moreover, writing boundary conditions in terms of $ f $, we have:
\begin{equation} \left\lbrace{\begin{array}{l} y = 0 \ \ \eta = 0: \ \ \left\lbrace{\begin{array}{l} u_x = \frac{U_\infty}{\delta_x}\sqrt{\frac{\nu x}{U_\infty}}f^{1}(0) = U_\infty f^{1}(0) = 0\\ \ \\ u_y = -\frac{U_\infty}{2}\sqrt{\frac{\nu}{U_\infty x}}\left({f(0) - \eta f^{(1)}(0)}\right) = 0 \end{array}}\right. \\ \ \\ y \rightarrow \infty \ \ \eta \rightarrow \infty: \ \ U_\infty f^{(1)}(\eta) \rightarrow U_\infty \end{array}}\right. \end{equation} Therefore, the complete form of the similarity solution of Blasius boundary layer problem can be expressed as:
\begin{equation} \left\lbrace{\begin{array}{l} f f^{(2)} + 2 f^{(3)} = 0\\ \ \\ b.c. : \left\lbrace{ \begin{array}{ll} \eta = 0: & f^{(1)}(0) = 0, \ f(0) = 0\\ \ \\ \eta \rightarrow \infty: & f^{(1)}(\eta) \rightarrow 1 \end{array}}\right. \end{array}}\right. \end{equation} But still it is possible to make the equations simpler:
\begin{equation} \left\lbrace{\begin{array}{l} f \psi^{(1)} + 2 \psi^{(2)} = 0\\ \ \\ f^{(1)} = \psi\\ \ \\ b.c. : \left\lbrace{ \begin{array}{l} \eta = 0: \ \ \ \psi(0) = 0, \ f(0) = 0\\ \ \\ \eta \rightarrow \infty: \ \psi(\eta) \rightarrow 1 \end{array}}\right. \end{array}}\right. \end{equation} Now this problem can be solved easier, using a centered scheme discretisation:
\begin{equation} A(\textbf{g}) = 0: \ \left\lbrace{\begin{array}{l} f_i \frac{\psi_{i+1} - \psi_{i-1}}{2\Delta \eta} + 2 \frac{\psi_{i+1} - 2\psi_{i} + \psi_{i-1}}{\Delta \eta^2} = 0\\ \ \\ \frac{f_{i+1} - f_{i-1}}{2\Delta \eta} - \psi_i = 0\\ \ \\ \text{Dirichlet's b.c.} \end{array}}\right. \end{equation} The unknowns can be ordered as follows:
\begin{equation} \left\lbrace{\begin{array}{c} g_1 \\ g_2 \\ \vdots \\ g_{2N-1} \\ g_{2N} \end{array}}\right\rbrace = \left\lbrace{\begin{array}{c} \psi_1 \\ f_1 \\ \vdots \\ \psi_N \\ f_N \end{array}}\right\rbrace \ \ \ , \ \ \ A(\textbf{g}) = \left({\begin{array}{c} A_{\text{1st Eq.}} \\ A_{\text{2nd Eq.}} \\ A_{\text{1st Eq.}} \\ A_{\text{2nd Eq.}} \\ \vdots \end{array}}\right)(\textbf{g}) \end{equation} The algebraic system is nonlinear, so it is necessary to apply an iterative method to solve this system of ODE equations. Here the Newton's method is implemented:
\begin{equation} A(\textbf{g}) = \textbf{0} \ \Longrightarrow \ A(\textbf{g}) \approx A(\hat{\textbf{g}}) + [J]_{\hat{\textbf{g}}}(\textbf{g} - \hat{\textbf{g}}) = \textbf{0} \end{equation} \begin{equation} \Longrightarrow \ \textbf{g}^{k+1} \approx \textbf{g}^k - [J]^{-1}_{\textbf{g}^{k}}A(\textbf{g}^k) \end{equation} where $ \textbf{g}^{k} $ is the $ k $-th iteration step and $ [J] $ is the Jacobian matrix. Instead of solving the system in terms of $ \textbf{g} $, it is more efficient to do it considering the residuals:
\begin{equation} \textbf{res}^{k+1} = \textbf{g}^{k+1} - \textbf{g}^k \approx - [J]^{-1}_{\textbf{g}^{k}}A(\textbf{g}^k) \end{equation} Moreover, the set of equations can be named as follows:
\begin{equation} \left\lbrace{\begin{array}{ll} \text{Odd Eqs.} \ \ A_{2i-1} :& f_i \frac{\psi_{i+1} - \psi_{i-1}}{2\Delta \eta} + 2 \frac{\psi_{i+1} - 2\psi_{i} + \psi_{i-1}}{\Delta \eta^2} = 0\\ \ \\ \text{Even Eqs.} \ \ A_{2i} :& \frac{f_{i+1} - f_{i-1}}{2\Delta \eta} - \psi_i = 0 \end{array}}\right. \end{equation} and the Jacobian matrix can be built:
\begin{equation} [J] = \left\lbrace{\begin{array}{l} J_{2i-1,2i-3} = \frac{\partial A_{2i-1}}{\partial \psi_{i-1}} \ , \ J_{2i-1,2i-1} = \frac{\partial A_{2i-1}}{\partial \psi_{i}} \ , \ J_{2i-1,2i+1} = \frac{\partial A_{2i-1}}{\partial \psi_{i+1}} \\ \\ J_{2i-1,2i-2} = \frac{\partial A_{2i-1}}{\partial f_{i-1}} \ , \ J_{2i-1,2i} = \frac{\partial A_{2i-1}}{\partial f_{i}} \ , \ \ \ \ J_{2i-1,2i+2} = \frac{\partial A_{2i-1}}{\partial f_{i+1}} \\ \\ J_{2i,2i-3} = \frac{\partial A_{2i}}{\partial \psi_{i-1}} \ , \ \ \ \ \ J_{2i,2i-1} = \frac{\partial A_{2i}}{\partial \psi_{i}} \ , \ \ \ \ \ \ J_{2i,2i+1} = \frac{\partial A_{2i}}{\partial \psi_{i+1}} \\ \\ J_{2i,2i-2} = \frac{\partial A_{2i}}{\partial f_{i-1}} \ , \ \ \ \ \ J_{2i,2i} = \frac{\partial A_{2i}}{\partial f_{i}} \ , \ \ \ \ \ \ \ \ \ J_{2i,2i+2} = \frac{\partial A_{2i}}{\partial f_{i+1}} \end{array}}\right. \end{equation} Furthermore, in order to have more accurate results, it is recommended to use a non-uniform grid. Then the equations can be written as:
\begin{equation} \left\lbrace{\begin{array}{l} f_i \frac{\psi_{i+1} - \psi_{i-1}}{\eta_{i+1} - \eta_{i-1}} + 2 \frac{\frac{\psi_{i+1} - \psi_{i}}{\eta_{i+1} - \eta_{i}} - \frac{\psi_{i} - \psi_{i-1}}{\eta_{i} - \eta_{i-1}}}{\frac{\eta_{i+1} - \eta_{i-1}}{2}} = 0\\ \ \\ \frac{f_{i+1} - f_{i-1}}{\eta_{i+1} - \eta_{i-1}} - \psi_i = 0 \end{array}}\right. \end{equation} And the boundary conditions are illustrated in the matrix form as follows:
\begin{equation} \hspace{-0.5cm} \left\lbrace{\begin{array}{c} g_1\\ g_2\\ \vdots\\ A(g^k_i)\\ \vdots\\ g_{2N-1}-1 \\ -g_{2N-1}+g^{(1)}_{2N} \end{array}}\right\rbrace + \left( \begin{array}{cccccc} 1 & 0 & \ldots & 0 & 0 & 0 \\ 0 & 1 & \ldots & 0 & 0 & 0 \\ & & & \vdots & & \\ & & & [J]_{\textbf{g}^k} & & \\ & & & \vdots & & \\ 0 & 0 & \ldots &0& 1 & 0 \\ 0 & 0 & \ldots &-\frac{1}{\eta_N - \eta_{N-1}} & -1 & \frac{1}{\eta_N - \eta_{N-1}} \\ \end{array} \right) \left\lbrace\begin{array}{c} res_1\\ res_2\\ \vdots\\ res_i\\ \vdots\\ res_{2N-1} \\ res_{2N} \\ \end{array}\right\rbrace^{k+1} = \left\lbrace{\begin{array}{c} 0\\ 0\\ \vdots\\ 0\\ \vdots\\ 0 \\ 0 \end{array}}\right\rbrace \end{equation}

Numerical results

Fig. d.1.2 illustrates the stream-wise velocity as a function of $ \eta $, Fig. d.1.3 the wall-normal velocity (both scaled according to their magnitude) and in Fig. d.1.4 the Newton's method iterations residuals is shown.