A CFD Tutorial in Julia: Introduction to Laminar Boundary-Layer Theory

Abstract

Numerical simulations of laminar boundary-layer equations are used to investigate the origins of skin-friction drag, flow separation, and aerodynamic heating concepts in advanced undergraduate- and graduate-level fluid dynamics/aerodynamics courses. A boundary-layer is a thin layer of fluid near a solid surface, and viscous effects dominate it. Students must understand the modeling of flow physics and implement numerical methods to conduct successful simulations. Writing computer codes to solve equations numerically is a critical part of the simulation process. Julia is a new programming language that is designed to combine performance and productivity. It is dynamic and fast. However, it is crucial to understand the capabilities of a new programming language before attempting to use it in a new project. In this paper, fundamental flow problems such as Blasius, Hiemenz, Homann, and Falkner-Skan flow equations are derived from scratch and numerically solved using the Julia language. We used the finite difference scheme to discretize the governing equations, employed the Thomas algorithm to solve the resulting linear system, and compared the results with the published data. In addition, we released the Julia codes in GitHub to shorten the learning curve for new users and discussed the advantages of Julia over other programming languages. We found that the Julia language has significant advantages in productivity over other coding languages. Interested readers may access the Julia codes on our GitHub page.

1. Introduction

Computational fluid dynamics (CFD) simulation is one of the vital steps of the design of a product that includes fluid motion. Since fluid dynamics are extremely complex and the motion equations have nonlinear terms, usage of numerical approaches is inevitable to simulate or predict fluid motion. Predictions can also be done with experiments. However, it is usually costlier than a regular CFD simulation. On the other hand, fundamental knowledge about the flow, that is planned to simulate, is necessary because of the required numerical approach decision. Learning the canonical flows well is extremely important in this point because most of the complex flow consists of a combination of a couple of canonical flows. An airfoil CFD simulation can consist of boundary-layer flow over the smooth part of the airfoil, mixing layer flow where the tail ends, and blunt body flow in the wake region. One airfoil simulation consists of three different canonical flows. A full understanding of the canonical flows is extremely important to model an accurate airfoil simulation.

Prandtl [1] stated that some of the terms in the Navier–Stokes equations can be neglected for the boundary-layer flows. As a result of this assumption, well-known boundary-layer equations arose. These approaches are still valid after a century and the resultant system of equations inspired lots of researchers in their studies. One of them was Blasius who is the Ph.D. student of Prandtl. Blasius [2] worked on the same problem as Prandtl did. However, he aimed to overcome the enigma of turbulence by considering the phenomenon of boundary-layer flow explained by Hager [3]. He further simplified the boundary-layer equations for a flat plate. He assumed that the flow is parallel. In other words, the velocity component in the parallel direction is not zero and the velocity in the transverse direction is zero. Moreover, he represented the resultant system of equations with a third-order ordinary differential equation, which is known as the Blasius similarity solution.

The Falkner–Skan similarity solution is another laminar similarity solution. Cebeci [4] explains the Falkner–Skan equations named after V. M. Falkner and Sylvia W. Skan. Falkner and Skan generalized the Blasius similarity solution for non-parallel flows, such as wedge flows and corner flows. The resultant equation can be used to predict the boundary-layer thickness for wedge flow. It has to be noted that these assumptions are available in the laminar region. Once the turbulence occurs, both of these similarity solutions will be inaccurate. In the incompressible region, the Falkner–Skan equation can be used without additional equation; however, after the compressibility limit, the temperature effect must be introduced to the system as an additional equation. If the Falkner–Skan equation is solved with the energy equation, the boundary-layer profile can be obtained in the compressible region.

One of the other Ph.D. students of Prandtl is Karl Hiemenz who worked on Hiemenz flow which is a type of stagnation point flow. Hiemenz [5] formulated and calculated the stagnation point problem as explained in the Schlichting [6]. The problem was a special case of the Falkner–Skan similarity solution. Howarth [7] also worked on the same problem and concluded similar results. Another Ph.D. student of Prandtl, Fritz Homann [8] worked on the same problem for axisymmetric bodies as Schlichting [6] explains. Homann’s similarity formulation was for a sphere while Hiemenz’s similarity formulation was for a cylinder. Both similarity solutions are being widely used for stagnation flows.

The aforementioned papers are the origins of the boundary-layer theory. As it is mentioned, the boundary-layer theory is the origin of many problems, such as the laminar to turbulent boundary-layer transition and flow separation. The main concern of these researches is to increase the performance of the vehicle because the transition increases the heat transfer and the vehicle requires a better thermal protection system at high speeds due to increased heat transfer rate. On the other hand, flow separation may lead to a lift force loss on the wing as a result, the performance of aircraft decreases. Since the focus of the present paper is the fundamentals of laminar boundary-layer theory, the details will not be provided about advanced researches. However, readers who are interested in details may check [9,10,11,12] for subsonic boundary-layer transition, Refs. [13,14,15,16,17,18,19,20,21] for supersonic/hypersonic boundary-layer transition, and Refs. [22,23,24] for flow separation. The other researches where boundary-layer flow is involved are [25,26,27,28,29,30].

Understanding the aforementioned fundamental flows are crucial for a senior undergraduate student or a graduate student in order to simulate more complex flows. However, modeling these equations in a computer environment requires more than knowledge about canonical flows instead it requires knowledge about programming languages as well. There are several learning modules and papers [31,32,33,34,35] for computational fluid dynamics simulation coding; however, there is not enough publication for Julia language [36]. It is a relatively new coding language among the other coding languages such as Fortran, C/C++, Python, and MATLAB but it is getting popular fast because it is trying to fill the gap between the language of the state-of-art CFD codes, Fortran and C/C++, and straight-forward/user-friendly languages, Python and MATLAB. Most of the state-of-art CFD codes are written in Fortran and C/C++ because it is so fast and random access memory usage (RAM) can be reduced drastically. On the other hand, MATLAB and Python have user-friendly and easy syntax which makes them favorites of students, with a price, which is the speed of computation. Julia is a language that tries to fill this gap between two different coding styles. It has a user-friendly environment, while also being fast. It provides a working space that helps users to write clear, high-level, generic code that resembles mathematical formulas. Julia’s ability to combine high-performance with productivity makes it a great choice for researchers working in different areas.

In this paper, fundamental flow problems such as Hiemenz flow, Homann flow, Blasius flow, and Falkner–Skan flow will be derived from scratch and modeled in the Julia environment. The finite-difference discretization in space with Thomas algorithm as linear system solver is used. Outputs of the code are provided and they are compared with the literature. Additionally, the advantages of the Julia language over other languages will be discussed. The computer codes and the implementation instructions will help students to understand the fundamental flows which will provide insight for student’s course work and researches. Using these examples, they can solve more complex flow types and develop their own codes. We make all these codes available on GitHub and they are accessible to everyone. We provide installation instruction for Julia and the required packages in Appendix A. The GitHub link of the codes also can be found in Appendix A.

2. Laminar Boundary-Layer Theory

Understanding canonical flows is crucial to understanding more complex flows such as flow over a wing or an airfoil. Hosseini et al. [37] studied flow over a wing section with a direct numerical simulation (DNS) study. They created a great video about how flow is formed over the aircraft. The video shows the development of the boundary-layer, how the trailing edge looks like a mixing-layer flow, and how Karman vortex street type of wake structures occurs. If the pre-study work is examined, it can be seen that the mesh structure is built on the flow prediction. Using a finer mesh on the critical regions is a key point of the CFD and it requires predictions about the possible flow behavior. The boundary-layer is the origin of many engineering problems in aerodynamics, including wing stall, the skin friction drag on an object, and the heat transfer that occurs in high-speed flight. In this present paper, boundary-layer theory will be examined under two subsections, which are laminar boundary-layer problems and stagnation point problems. In this present paper, fundamental fluid dynamics problems related to boundary-layer theory will be derived and implemented in a relatively new programming language, Julia. The contribution is employing the Julia language. The governing equations and solution methodologies are already published in the literature. The interested reader should refer to the additional references [6,38] for detailed derivations. The manuscript may enable students to adopt the programming language easily.

2.1. Laminar Boundary-Layer Flow Problems

Velocity distribution over a flat plate can be represented with a similarity solution. In this subsection, Blasius and Falkner–Skan similarity solutions will be derived from scratch and they will be solved numerically in the Julia environment. The Julia codes will be available to shorten the learning curve.

2.1.1. Blasius Flow Problem

Schematic description of flow over a flat plate, in other words, Blasius boundary-layer flow can be illustrated as in Figure 1. The governing equations for Blasius flow are boundary-layer equations. These equations can be obtained by non-dimensionalizing the two-dimensional Navier–Stokes equations that are:

(1)$\begin{array}{cc}\hfill \frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}& =0\hfill \end{array}$

(2)$\begin{array}{cc}\hfill u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}& =-\frac{1}{\rho }\frac{\partial p}{\partial x}+\nu \left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right)\hfill \end{array}$

(3)$\begin{array}{cc}\hfill u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}& =-\frac{1}{\rho }\frac{\partial p}{\partial y}+\nu \left(\frac{{\partial }^{2}v}{\partial x^2}+\frac{{\partial }^{2}v}{\partial y^2}\right),\hfill \end{array}$

where u and v are the velocities in the x-direction and y-direction, respectively. p is the pressure, $\nu$ is the kinematic viscosity, $\rho$ is the density. If Equations (1)–(3) are non-dimensionalized with:

(4)$\begin{array}{c}\hfill u^{*}=\frac{u}{U_{\infty }},v^{*}=\frac{v}{U_{\infty }},p^{*}=\frac{p}{\rho U_{\infty }^2},\end{array}$

(5)$\begin{array}{c}\hfill x^{*}=\frac{x}{L},y^{*}=\frac{y}{L},\end{array}$

where L is the plate length and $U_{\infty }$ is the free-stream velocity. The two-dimensional, incompressible non-dimensional Navier–Stokes equations can be shown as:

(6)$\begin{array}{cc}\hfill \frac{\partial u^{*}}{\partial x^{*}}+\frac{\partial v^{*}}{\partial y^{*}}& =0\hfill \end{array}$

(7)$\begin{array}{cc}\hfill u^{*}\frac{\partial u^{*}}{\partial x^{*}}+v^{*}\frac{\partial u^{*}}{\partial y^{*}}& =-\frac{\partial p^{*}}{\partial x^{*}}+\frac{1}{R{e}_{\infty }}\left(\frac{{\partial }^2{u}^{*}}{\partial x^{{*}^2}}+\frac{{\partial }^2{u}^{*}}{\partial y^{{*}^2}}\right)\hfill \end{array}$

(8)$\begin{array}{cc}\hfill u^{*}\frac{\partial v^{*}}{\partial x^{*}}+v^{*}\frac{\partial v^{*}}{\partial y^{*}}& =-\frac{\partial p^{*}}{\partial y^{*}}+\frac{1}{R{e}_{\infty }}\left(\frac{{\partial }^2{v}^{*}}{\partial x^{{*}^2}}+\frac{{\partial }^2{v}^{*}}{\partial y^{{*}^2}}\right).\hfill \end{array}$

Star superscript (${\left[\right]}^{*}$) corresponds to non-dimensional variable. In a boundary-layer flow, some variables are smaller than others. For example, $u^{*}=O\left(1\right)$ and $x^{*}=O\left(1\right)$ so $\frac{\partial u^{*}}{\partial x^{*}}=O\left(1\right)$. In the continuity equation, two terms must be in the same order to obtain 0 so $\frac{\partial v^{*}}{\partial y^{*}}=O\left(1\right)$. It can be seen that from Figure 1, $y=O\left({\delta }^{*}\right)$ where ${\delta }^{*}$ is the non-dimensionalized boundary-layer thickness. From here, it can be concluded that $v=O\left({\delta }^{*}\right)$. It has to be noted that ${\delta }^{*}<<1$. If the same approach is applied to momentum equations, the final system of equations will be:

(9)$\begin{array}{cc}\hfill \frac{\partial u^{*}}{\partial x^{*}}+\frac{\partial v^{*}}{\partial y^{*}}& =0\hfill \end{array}$

(10)$\begin{array}{cc}\hfill u^{*}\frac{\partial u^{*}}{\partial x^{*}}+v^{*}\frac{\partial u^{*}}{\partial y^{*}}& =-\frac{\partial p^{*}}{\partial x^{*}}+\frac{1}{R{e}_{\infty }}\frac{{\partial }^2{u}^{*}}{\partial y^{{*}^2}}\hfill \end{array}$

(11)$\begin{array}{cc}\hfill \frac{\partial p^{*}}{\partial y^{*}}& =0.\hfill \end{array}$

For an inviscid flow, where the Reynolds number is so high, the viscous term can be neglected. The pressure is constant as obtained in the Equation (11) and $v=0$. Equation (10) will be:

(12)$-\frac{d{P}_e^{*}}{d{x}^{*}}=U_e^{*}\frac{d{U}_e^{*}}{d{x}^{*}},$

where $U_e$ and $P_e$ are the inviscid velocity and pressure. Finally, the two-dimensional boundary-layer equations in the dimensionless form can be written as:

(13)$\begin{array}{cc}\hfill \frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}& =0\hfill \end{array}$

(14)$\begin{array}{cc}\hfill u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}& =U_e\frac{d{U}_e}{\partial x}+\nu \frac{{\partial }^{2}u}{\partial y^2}.\hfill \end{array}$

It has to be noted that star superscript is suppressed to avoid confusion. In order to solve this problem, similarity transformation can be used. The similarity parameter ($\eta$) can be chosen as:

(15)$\eta =\frac{y}{\sqrt{\frac{\nu x}{U_{\infty }}}}.$

After that, velocity components can be shown as:

(16)$\begin{array}{cc}\hfill u(x,y)& =U_{\infty }F\left(\eta \right)\hfill \end{array}$

(17)$\begin{array}{cc}\hfill v(x,y)& =\frac{G\left(\eta \right)}{\sqrt{\frac{x}{U_{\infty }\nu }}}.\hfill \end{array}$

If Equations (13) and (14) are rearranged with the given velocities, the boundary-layer equations will be:

(18)$\begin{array}{cc}\hfill G^{\prime }& =\frac{\eta }{2}{F}^{\prime }\hfill \end{array}$

(19)$\begin{array}{cc}\hfill F^{\prime \prime }& =-\frac{1}{2}\eta F{F}^{\prime }+G{F}^{\prime }.\hfill \end{array}$

In order to make the next procedure easier, it can be assumed that $F=f^{\prime }$. After that, Equation (18) can be integrated and put into Equation (19). The final ordinary differential equation can be obtained as:

(20)$2f^{\prime \prime \prime }+f{f}^{\prime \prime }=0.$

The boundary conditions of the system will be:

(21)$\begin{array}{cc}\hfill u(x,0)& =U_{\infty }F(\eta =0)=0⟶F\left(0\right)=f^{\prime }\left(0\right)=0\hfill \end{array}$

(22)$\begin{array}{cc}\hfill v(x,0)& =\frac{G(\eta =0)}{\sqrt{\frac{x}{U_{\infty }\nu }}}=0⟶G\left(0\right)=f\left(0\right)=0\hfill \end{array}$

(23)$\begin{array}{cc}\hfill u(x,y\stackrel{}{\to }\infty )& =U_{\infty }F(\eta \stackrel{}{\to }\infty )=1⟶F(\infty )=f^{\prime }(\infty )=f(\infty )=1.\hfill \end{array}$

The given assumptions reduced the system of second-order partial differential equations into a third-order ordinary differential equation. The computational approach will be examined in the Julia framework.

2.1.2. Numerical Solution of Blasius Flow Problem

In this computational section, the solution of the third-order ordinary differential equation will be transformed into two equations and they will be solved with Thomas algorithm [39] which is one of the best methods for the tridiagonal matrices. It has to be noted that the same equations can be solved with any other methods such as Runge-Kutta, Runge-Kutta-Fehlberg, compact finite difference, high-order finite-difference; however, in this present paper, authors used the finite difference method for spatial discretization with the Thomas algorithm for the linear system solution. If it is assumed that $f^{\prime }=h$ and $f=p$ in Equation (20), the system of equations will be:

(24)$\begin{array}{cc}\hfill 2h^{\prime \prime }+p{h}^{\prime }=& 0\hfill \end{array}$

(25)$\begin{array}{cc}\hfill p^{\prime }-h=& 0,\hfill \end{array}$

with the following boundary conditions:

(26)$\begin{array}{c}\hfill p\left(0\right)=0\end{array}$

(27)$\begin{array}{c}\hfill h\left(0\right)=0\end{array}$

(28)$\begin{array}{c}\hfill h(\infty )=1.\end{array}$

If the second-order central finite difference method is applied to h variable:

(29)$\begin{array}{cc}\hfill h_n^{\prime \prime }=& \frac{h_{n+1}-2h_n+h_{n-1}}{{(\Delta \eta )}^2}+O\left({(\Delta \eta )}^2\right)\hfill \end{array}$

(30)$\begin{array}{cc}\hfill h_n^{\prime }=& \frac{h_{n+1}-h_{n-1}}{2\Delta \eta }+O\left({(\Delta \eta )}^2\right).\hfill \end{array}$

In this paper, details of the finite difference and derivation of it from Taylor’s series will not be covered. If the reader is curious about the derivation of them, the book of Moin [40] can be checked. If the finite difference approach for h is substituted into the Equation (24), the final equation will be:

(31)$A_n{h}_{n+1}+B_n{h}_n+C_n{h}_{n-1}=0,$

where the $A_n$, $B_n$, and $C_n$ are:

(32)$\begin{array}{cc}\hfill A_n=& \frac{2}{\Delta {\eta }^2}+\frac{p_n}{2\Delta \eta }\hfill \end{array}$

(33)$\begin{array}{cc}\hfill B_n=& \frac{-4}{\Delta {\eta }^2}\hfill \end{array}$

(34)$\begin{array}{cc}\hfill C_n=& \frac{2}{\Delta {\eta }^2}-\frac{p_n}{2\Delta \eta },\hfill \end{array}$

with the boundary conditions for the system the tridiagonal system ($Ax=b$) can be shown as:

(35)$\underset{A}{\underbrace{\left[\begin{array}{cccccc}1& 0& 0& 0& \cdots & 0\\ A_2& B_2& C_2& 0& \cdots & 0\\ ⋮& \ddots & \ddots & \ddots & & ⋮\\ 0& 0& A_{N-2}& B_{N-2}& C_{N-2}& 0\\ 0& 0& 0& A_{N-1}& B_{N-1}& C_{N-1}\\ 0& 0& 0& 0& 0& 1\end{array}\right]}}\underset{x}{\underbrace{\left[\begin{array}{c}{h}_1\\ h_2\\ ⋮\\ h_{N-2}\\ h_{N-1}\\ h_N\end{array}\right]}}=\underset{b}{\underbrace{\left[\begin{array}{c}0\\ 0\\ ⋮\\ 0\\ 0\\ 1\end{array}\right]}}$

It has to be noted that the first and last rows of the matrix are known and they are coming from the boundary conditions. In order to solve this problem, the Thomas algorithm can be applied. If relation between $h_n$ and $h_{n-1}$ is:

(36)$\begin{array}{cc}\hfill h_n=& G_n+H_n{h}_{n+1}\hfill \end{array}$

(37)$\begin{array}{cc}\hfill h_{n-1}=& G_{n-1}+H_{n-1}{h}_n.\hfill \end{array}$

If Equation (37) is substituted into Equation (31), the final equation will be:

(38)$h_n=\frac{-C_n{G}_{n-1}}{B_n+C_n{H}_{n-1}}+\frac{-A_n}{B_n+C_n{H}_{n-1}}{h}_{n+1}.$

If Equations (36) and (38) are compared, $G_n$ and $H_n$ coefficients can be found as:

(39)$\begin{array}{cc}\hfill G_n=& \frac{-C_n{G}_{n-1}}{B_n+C_n{H}_{n-1}}\hfill \end{array}$

(40)$\begin{array}{cc}\hfill H_n=& \frac{-A_n}{B_n+C_n{H}_{n-1}}.\hfill \end{array}$

The problem here is that one has to know $G_1$ and $H_1$ in order to start the calculation. $H_1$ can be assumed as 0 and the $G_1$ can be calculated as 0 from the boundary conditions. Once $G_1$ and $H_1$ are calculated, $G_n$ and $H_n$ can be calculated from Equations (39) and (40). The given formulations can be implemented as it is shown in Listing 1.

The problem in Listing 1 is the usage of the p and h values which are not defined yet. In order to start to iteration, the initial assumption for h must be given. Most of the time, the linear assumption is the best assumption. Figure 2 shows the schematic of the profile transformation after some iterations. Once the initialization of h is done, p can be calculated from the integration of $p^{\prime }=h$. The discrete integration of p can be written as:

(41)$p_n=p_{n-1}+{\int }_{{\eta }_{n-1}}^{{\eta }_n}hd\eta .$

Listing 1. Implementation of Thomas algorithm for Blasius profile.

If the integration is done numerically with the trapezoidal rule [41], the initial guess for h and p can be calculated in Julia as shown in Listing 2.

After the initialization is done, Listing 1 can be run with p calculation of Listing 2, until the change in h is smaller than an arbitrary parameter $ϵ$ which can be taken as $1\times {10}^{-8}$. The final profile will be as shown in Figure 3. The profile is also compared with Schlichting [6] in order to validate the results. It has to be noted that, Schlichting used $\eta =\frac{y}{\sqrt{\frac{2\nu x}{U_{\infty }}}}$ as similarity coordinate so the figure is plotted according to this nondimensionalization.

Listing 2. Implementation of initial velocity guess and initial p for Blasius profile.

2.1.3. Falkner–Skan Flow Problem

The Falkner–Skan similarity solution can be considered as a family of the similarity solutions since it is a general solution that includes Blasius flow, Hiemenz flow (see Section 2.2.1), and more. Falkner–Skan equation cannot represent Homann flow (see Section 2.2.3) because Homann flow is an axisymmetric flow. On the other hand, Falkner–Skan flow is a two-dimensional flow. The name ‘similarity solution’ arises from the solutions at two arbitrary stations which are related to one another by means of a scale factor. A flow that can be represented with boundary-layer Equations (13) and (14), and satisfies the following equation:

(42)$\frac{u\left(x_1\right),\frac{y}{g\left(x_1\right)}}{U_e\left(x_1\right)}=\frac{u\left(x_2\right),\frac{y}{g\left(x_2\right)}}{U_e\left(x_2\right)},$

can be considered as self-similar. If boundary-layer Equations (13) and (14) are considered and similarity a transformation is assumed as:

(43)$u(x,y)=U_e\left(x\right)f^{\prime }\left(\xi \right),\xi =\frac{y}{g\left(x\right)}.$

If the continuity equation of boundary-layer Equation (13) is modified with the given transformation, the final equation will be:

(44)$\frac{\partial v}{\partial \xi }=-U_e^{\prime }{f}^{\prime }g+U_e{f}^{\prime \prime }{g}^{\prime }\xi .$

If Equation (44) is integrated over $\xi$ to find the velocity v, the v velocity will be:

(45)$v(x,\xi )=-{\left(U_{e}g\right)}^{\prime }f+U_e\xi g^{\prime }{f}^{\prime }+H\left(x\right).$

$H\left(x\right)$ can be calculated from the boundary conditions on the wall where $\xi$ is zero and correspondingly, velocities and $f^{\prime }\left(0\right)$ are zero. $H\left(x\right)$ will be 0 if $f\left(0\right)$ is chosen as 0. If calculated derivatives and velocities are substituted into Equation (14), after some simplification, the final equation will be:

(46)$f^{\prime \prime \prime }+\frac{g}{\nu }{\left(U_{e}g\right)}^{\prime }f{f}^{\prime \prime }+\frac{g^2{U}_e^{\prime }}{\nu }(1-f^{\prime })=0.$

Since f must be the function of $\xi$ only and must not be a function of x, the coefficients of second and third terms must be constant. The final Falkner–Skan equation for the family of self-similar solutions:

(47)$f^{\prime \prime \prime }+\alpha f{f}^{\prime \prime }+\beta (1-f^{\prime 2})=0,$

where $\alpha =\frac{g}{\nu }{\left(U_{e}g\right)}^{\prime }$ and $\beta =\frac{g^2{U}_e^{\prime }}{\nu }$. The $\alpha$ and $\beta$ can be further solved for the velocity scale and length scale. If the following consideration is applied:

(48)$\begin{array}{cc}\hfill {\alpha }^{*}=& g{\left(U_{e}g\right)}^{\prime }\hfill \end{array}$

(49)$\begin{array}{cc}\hfill {\beta }^{*}=& g^2{U}_e^{\prime }.\hfill \end{array}$

If ${\left(U_e{g}^2\right)}^{\prime }$ is written in terms of ${\alpha }^{*}$ and ${\beta }^{*}$, the obtained equation integrated with respect to x, the resultant equation will be:

(50)$U_e{g}^2=(2{\alpha }^{*}-{\beta }^{*})x+c.$

The constant of the integration represents a shift in the origin on x. Hence it doesn’t affect the result and it also can be calculated from the stagnation point where $x=0$ and $U_e=0$ as a result, $c=0$. If relation of ${\beta }^{*}=g^2{U}_e^{\prime }$ is divided by Equation (50) and integrated with respect to x, $U_e$ can be calculated as:

(51)$U_e=c_1{x}^m,$

where $m=\frac{{\beta }^{*}}{2{\alpha }^{*}-{\beta }^{*}}$ and $c_1$ is a positive or negative constant which depends on the sign of $U_e$. It can be concluded from these calculations that similar solutions exist when the inviscid velocity is proportional to x raised to some power. Next, Equation (50) can be used by taking $c=0$ to calculate the g which is:

(52)$g=\sqrt{\frac{2{\alpha }^{*}}{c_1(1+m)}{x}^{1-m}}.$

Self-similar boundary-layers occur when the external velocity is the simple power law ($U_e=U_0{(x/L)}^m$), where the arbitrary constants $U_0$ and L have the same sign as U and x. The similarity variable for these kinds of flows can be written as:

(53)$\eta =\frac{y}{\delta }=\frac{y}{\sqrt{\pm \frac{\nu x}{U_e}}}=\frac{y}{\sqrt{\pm \frac{\nu L}{U_0}{\left(\frac{x}{L}\right)}^{1-m}}}.$

When $U_e$ and x have the same signs, the Falkner–Skan equation can be written as:

(54)$f^{\prime \prime \prime }+\frac{1}{2}(m+1)f{f}^{\prime \prime }+m(1-f^{\prime 2})=0.$

The two arbitrary constants $\alpha$ and $\beta$ have been reduced to one constant m by fixing the scale for the function $\delta \left(x\right)$. The boundary conditions of the equation are:

(55)$\begin{array}{c}\hfill f\left(0\right)=0\end{array}$

(56)$\begin{array}{c}\hfill f^{\prime }\left(0\right)=0\end{array}$

(57)$\begin{array}{c}\hfill f^{\prime }(\eta \stackrel{}{\to }\infty )=1.\end{array}$

If the Falkner–Skan Equation (54) is carefully examined, it can be seen that when the constant m is 0, the equation will be Blasius flow. Moreover, if the constant m is 1, the equation will be Hiemenz flow (see Section 2.2.1). This important point is stated before and it is emphasized one more time after the derivation. Another great usage of the Falkner–Skan equation is to simulate the boundary-layer over a wedge with half-angle $\theta =\frac{m\pi }{m+1}$ when the m is between 0 and 1. If the m is in between 1 and 2, the Falkner–Skan equation will solve a corner flow with $\theta >\frac{\pi }{2}$. The visual schematic of four different physical flows that can be calculated from the Falkner–Skan equation can be seen in Figure 4. One interesting point of the equation is that when $m=-0.0904$, the obtained profile will have zero-shear at the wall which corresponds to the verge of the separation point for all x stations.

2.1.4. Numerical Solution of Falkner–Skan Flow Problem

In this section, the Falkner–Skan equation will be solved with the Thomas algorithm and central finite difference scheme. In order to do that, the third-order ordinary differential equation should be reduced to second-order and first-order differential equations as it is done in Blasius flow. If it is assumed that $f^{\prime }=h$ and $f=p$ in Equation (54), the system of equations will be:

(58)$\begin{array}{cc}\hfill h^{\prime \prime }+\frac{1}{2}(m+1)p{h}^{\prime }+m(1-h^2)=& 0\hfill \end{array}$

(59)$\begin{array}{cc}\hfill p^{\prime }-h=& 0,\hfill \end{array}$

where the boundary conditions of the system are:

(60)$\begin{array}{c}\hfill p\left(0\right)=0\end{array}$

(61)$\begin{array}{c}\hfill h\left(0\right)=0\end{array}$

(62)$\begin{array}{c}\hfill h(\infty )=1\end{array}$

If finite difference scheme is applied to h variable, the final system of equations will be:

(63)$A_n{h}_{n+1}+B_n{h}_n+C_n{h}_{n-1}+D_n=0,$

where the $A_n$, $B_n$, $C_n$, and $D_n$ are:

(64)$\begin{array}{cc}\hfill A_n=& \frac{1}{\Delta {\eta }^2}+(m+1)\frac{p_n}{4\Delta \eta }\hfill \end{array}$

(65)$\begin{array}{cc}\hfill B_n=& \frac{-2}{\Delta {\eta }^2}-m{h}_n\hfill \end{array}$

(66)$\begin{array}{cc}\hfill C_n=& \frac{1}{\Delta {\eta }^2}-(m+1)\frac{p_n}{4\Delta \eta }\hfill \end{array}$

(67)$\begin{array}{cc}\hfill D_n=& m.\hfill \end{array}$

The relation between $h_n$ and $h_{n-1}$ is taken as it is done in Blasius flow. When Equation (37) is substituted into Equation (63), the final $G_n$ and $H_n$ will be:

(68)$\begin{array}{cc}\hfill G_n=& \frac{-C_n{G}_{n-1}+D_n}{B_n+C_n{H}_{n-1}}\hfill \end{array}$

(69)$\begin{array}{cc}\hfill H_n=& \frac{-A_n}{B_n+C_n{H}_{n-1}}.\hfill \end{array}$

The implementation of these variables in the Julia environment can be seen in Listing 3.

Listing 3. Implementation of the Thomas algorithm for the Falkner–Skan profile.

The final profiles for the varying m values can be seen in the Figure 5. The profiles are also compared with Schlichting [6] in order to validate the results. It has to be noted that, Schlichting used $\eta =\sqrt{\frac{m+1}{2}}\frac{y}{\sqrt{\frac{\nu x}{U_{\infty }}}}$ as similarity coordinate so the figure is plotted according to this nondimensionalization.

2.2. Stagnation Point Flow Problems

Boundary-layer velocity distribution over a wall that is perpendicular to the flow velocity vector can be represented with a similarity solution. In this subsection, Hiemenz and Homann similarity solutions will be derived from scratch and they will be solved numerically in the Julia environment. The Julia codes will be available to shorten the learning curve.

2.2.1. Hiemenz Flow Problem

This problem is that of a fluid flow that is parallel to the y-axis in the far-field impinging on a wall that coincides with the x-axis. The flow which is perpendicular to a cylinder can be assumed as the Hiemenz flow around the stagnation point. The schematic description of the Hiemenz flow can be seen in Figure 6. In the Hiemenz flow, viscous forces away from the wall become so small in comparison with the inertia forces, particularly when the Reynolds number is large. In this case, inviscid irrotational flow assumption ($\xi =\nabla \times U=0$) can be done. The velocity can be represented with a scalar function, $\varphi$, which is the velocity potential. The velocities in the x-direction and y-direction can be written as:

(70)$u=\frac{\partial \varphi }{\partial x},v=\frac{\partial \varphi }{\partial y}.$

If the dimensional Navier–Stokes equations (Equations (1)–(3)) are considered for this flow as well, and the continuity equation becomes:

(71)$\frac{{\partial }^2\varphi }{\partial x^2}+\frac{{\partial }^2\varphi }{\partial y^2}=0.$

On the wall, in the absence of viscosity, the flow can slip in the x-direction but in y-direction the velocity must be zero ($\frac{\partial \varphi }{\partial y}{|}_{y=0}=0$) because of the no-penetration boundary condition. The inviscid flow solution for the potential function was found to be:

(72)${\varphi }_{(x,y)}=\frac{a}{2}(x^2-y^2),$

where a is a constant that depends on the freestream flow and the body shape. This solution satisfies the governing equations for inviscid, irrotational flow, and the boundary conditions which can be shown as:

(73)$u=\frac{\partial \varphi }{\partial x}=ax,v=\frac{\partial \varphi }{\partial y}=-ay.$

Since potential the function is obtained, the stream-function can be calculated from the potential function. The velocity components, in terms of stream-function, ${\psi }_{(x,y)}$:

(74)$u=ax=\frac{\partial \psi }{\partial y},v=-ay=-\frac{\partial \psi }{\partial x}.$

If Equation (74) is integrated for x and y, the stream function can be found as:

(75)$\psi =axy+c.$

In order to find the pressure in the inviscid flow, Bernoulli equation [42] can be used. The pressure from the Bernoulli equation is:

(76)$p_0-p=\frac{a^2\rho }{2}(x^2+y^2),$

where $p_0$ is the stagnation pressure and $\rho$ is the density. So far, flow is assumed as inviscid flow; however viscous forces will modify the inviscid solution, in particular in the y-direction. Hence, the viscous velocity components can be assumed as:

(77)$\begin{array}{cc}\hfill u(x,y)=& xg\left(y\right)\hfill \end{array}$

(78)$\begin{array}{cc}\hfill v(x,y)=& -f\left(y\right)\hfill \end{array}$

If Equations (77) and (78) are substituted in the Navier–Stokes equations (Equations (1)–(3)), the final system of equations will be:

(79)$\begin{array}{cc}\hfill g=& f^{\prime }\hfill \end{array}$

(80)$\begin{array}{cc}\hfill x{g}^2-xf{g}^{\prime }=& -\frac{1}{\rho }\frac{\partial p}{\partial x}+\nu x{g}^{\prime \prime }\hfill \end{array}$

(81)$\begin{array}{cc}\hfill f{f}^{\prime }=& -\frac{1}{\rho }\frac{\partial p}{\partial y}+\nu f^{\prime \prime }.\hfill \end{array}$

Once Equation (79) is substituted into Equations (80) and (81), two non-linear ordinary differential equations can be obtained as:

(82)$\begin{array}{cc}\hfill x{f}^{\prime 2}-xf{f}^{\prime \prime }=& -\frac{1}{\rho }\frac{\partial p}{\partial x}+\nu x{f}^{\prime \prime \prime }\hfill \end{array}$

(83)$\begin{array}{cc}\hfill f{f}^{\prime }=& -\frac{1}{\rho }\frac{\partial p}{\partial y}+\nu f^{\prime \prime }\hfill \end{array}$

It has to be noted that these two equations are coupled, which means they have to be solved together. If the x-momentum and y-momentum are solved for pressure, it can be seen that the pressure will be in the form of:

(84)$p(x,y)=\frac{1}{2}{x}^2\left(constant\right)+H\left(y\right),$

where $H\left(y\right)$ is a function depends on only y. On the other hand, if the inviscid pressure compared with the viscous pressure as $y\stackrel{}{\to }\infty$, the pressure gradient in x-direction and y-direction can be found as:

(85)$\begin{array}{cc}\hfill \frac{\partial p}{\partial x}=& -a^2\rho x\hfill \end{array}$

(86)$\begin{array}{cc}\hfill \frac{\partial p}{\partial y}=& -\frac{a^2\rho }{2}{F}^{\prime }\left(y\right).\hfill \end{array}$

If the pressure gradients are substituted into Equations (82) and (83), the resultant system of ordinary differential equations can be found as:

(87)$\begin{array}{cc}\hfill f^{\prime 2}-f{f}^{\prime \prime }=& a^2+\nu f^{\prime \prime \prime }\hfill \end{array}$

(88)$\begin{array}{cc}\hfill f{f}^{\prime }=& \frac{a^2}{2}{F}^{\prime }-\nu f^{\prime \prime }.\hfill \end{array}$

Momentum equations are decoupled in this system of equations so once, $f\left(y\right)$ is calculated from Equation (87), $F\left(y\right)$ in the Equation (88) can be solved. The boundary conditions of the system are:

(89)$\begin{array}{c}\hfill u(x,0)=0\stackrel{}{\to }{f}_{(y=0)}^{\prime }=0\end{array}$

(90)$\begin{array}{c}\hfill v(x,0)=0\stackrel{}{\to }{f}_{(y=0)}=0\end{array}$

(91)$\begin{array}{c}\hfill u(x,y\stackrel{}{\to }\infty )=ax\stackrel{}{\to }{f}_{\infty }^{\prime }=a\end{array}$

(92)$\begin{array}{c}\hfill v(x,y\stackrel{}{\to }\infty )=-ay\stackrel{}{\to }{f}_{\infty }^{\prime \prime }=0\end{array}$

Additionally, $F(y\stackrel{}{\to }\infty )=y^2$ can be obtained from y-momentum equation as $y\stackrel{}{\to }\infty$. The final Equations (87) and (88) can be solved with the given boundary conditions; however, it is possible to get rid of the dependence on a. In order to do that, affine transformation [43] can be used. If it is assumed that $f\left(y\right)=A\phi \left(\eta \right)$ and $y=B\eta$ where A and B are constant to be determined, the f derivative functions can be calculated as:

(93)$\begin{array}{cc}\hfill f^{\prime }=& \frac{A}{B}{\phi }^{\prime }\hfill \end{array}$

(94)$\begin{array}{cc}\hfill f^{\prime \prime }=& \frac{A}{B^2}{\phi }^{\prime \prime }\hfill \end{array}$

(95)$\begin{array}{cc}\hfill f^{\prime \prime \prime }=& \frac{A}{B^3}{\phi }^{\prime \prime \prime }\hfill \end{array}$

If Equations (93)–(95) substituted into the Equations (87) and (88) and assume that $\frac{\nu }{AB}=1$ and $\frac{a^2{B}^2}{A^2}=1$, the A and B coefficients can be calculated as:

(96)$\begin{array}{cc}\hfill A=& \sqrt{a\nu }\hfill \end{array}$

(97)$\begin{array}{cc}\hfill B=& \sqrt{\frac{\nu }{a}}\hfill \end{array}$

The given transformation allows to reduce two ordinary differential equations into one third-order ordinary differential equation as:

(98)${\phi }^{\prime \prime \prime }+\phi {\phi }^{\prime \prime }-{\phi }^{\prime 2}+1=0$

where the boundary conditions are:

(99)$\begin{array}{cc}\hfill f^{\prime }\left(0\right)=& 0\stackrel{}{\to }{\phi }^{\prime }\left(0\right)=0\hfill \end{array}$

(100)$\begin{array}{cc}\hfill f\left(0\right)=& 0\stackrel{}{\to }\phi \left(0\right)=0\hfill \end{array}$

(101)$\begin{array}{cc}\hfill f^{\prime }(\infty )=& a\stackrel{}{\to }{\phi }^{\prime }(\infty )=1\hfill \end{array}$

and the velocities:

(102)$\begin{array}{c}\hfill u=ax{\phi }^{\prime }\left(\eta \right)\end{array}$

(103)$\begin{array}{c}\hfill v=-\sqrt{a\nu }\phi \left(\eta \right)\end{array}$

the resultant equation along with the boundary conditions can be solved with different numerical approaches such as Runge-Kutta, Runge-Kutta-Fehlberg, compact finite difference, high-order finite difference; however, in this present paper, the Thomas algorithm with finite difference discretization will be used to solve the Hiemenz profile in the Julia framework.

2.2.2. Numerical Solution of Hiemenz Flow Problem

The computational approach for the Hiemenz flow is similar to the Blasius solution since the final ordinary equation of the Hiemenz flow (Equation (98)) is similar to the Blasius similarity solution (Equation (20)). If it is assumed that ${\phi }^{\prime }=h$ and $\phi =p$ in Equation (98), the system of equations will be:

(104)$\begin{array}{cc}\hfill h^{\prime \prime }+p{h}^{\prime }-h^2+1=& 0\hfill \end{array}$

(105)$\begin{array}{cc}\hfill p^{\prime }-h=& 0,\hfill \end{array}$

where the boundary conditions of the system are:

(106)$\begin{array}{c}\hfill p\left(0\right)=0\end{array}$

(107)$\begin{array}{c}\hfill h\left(0\right)=0\end{array}$

(108)$\begin{array}{c}\hfill h(\infty )=1.\end{array}$

If finite difference scheme is applied to the h variable as in Blasius flow. The final system of equations will be:

(109)$A_n{h}_{n+1}+B_n{h}_n+C_n{h}_{n-1}+D_n=0,$

where the $A_n$, $B_n$, $C_n$, and $D_n$ are:

(110)$\begin{array}{cc}\hfill A_n=& \frac{1}{\Delta {\eta }^2}+\frac{p_n}{2\Delta \eta }\hfill \end{array}$

(111)$\begin{array}{cc}\hfill B_n=& \frac{-2}{\Delta {\eta }^2}-h_n\hfill \end{array}$

(112)$\begin{array}{cc}\hfill C_n=& \frac{1}{\Delta {\eta }^2}-\frac{p_n}{2\Delta \eta }\hfill \end{array}$

(113)$\begin{array}{cc}\hfill D_n=& 1.\hfill \end{array}$

The relation between $h_n$ and $h_{n-1}$ is taken as it is done in Blasius flow. When Equation (37) is substituted into the Equation (109), the final $G_n$ and $H_n$ coefficients can be found as:

(114)$\begin{array}{cc}\hfill G_n=& \frac{-C_n{G}_{n-1}+D_n}{B_n+C_n{H}_{n-1}}\hfill \end{array}$

(115)$\begin{array}{cc}\hfill H_n=& \frac{-A_n}{B_n+C_n{H}_{n-1}}.\hfill \end{array}$

In order to start the calculation, $H_1$ can be assumed as 0 and the $G_1$ can be calculated as 0 from the boundary conditions. Once $G_1$ and $H_1$ are calculated, $G_n$ and $H_n$ can be calculated from Equations (114) and (115). p can be calculated as it is done for Blasius flow:

(116)$p_n=p_{n-1}+{\int }_{{\eta }_{n-1}}^{{\eta }_n}hd\eta .$

The system of equations can be implemented in the Julia environment as it is shown in Listing 4. The linear profile assumption is taken as the initial condition of h (see Figure 2) as it is done for the Blasius solution. The reason $hp$ and $errorProfile$ variables are used in the Listing 4 is that the linear profile converges to the Hiemenz profile in each iteration so in order to check the difference between previous and present profiles, the solution vector from the previous iteration is copied and it is compared with the new solution vector. If the difference between these two solution vectors is less than $ϵ$, which is an arbitrary limit and can be taken as $1\times {10}^{-8}$, then it can be said that the solution has converged. It has to be noted that $ϵ$ can be taken as any number; however, if it is small, the solution will be more accurate. The final result of the Hiemenz flow can be seen in Figure 7. The results are also validated with White [44].

2.2.3. Homann Flow Problem

Homann flow is similar to Hiemenz flow. The only difference is that Homann flow is an axisymmetric version of the Hiemenz flow. The same schematic (Figure 6) can represent this flow as well; however, in this flow, the circle is the projection of a sphere. On the other hand, it was the projection of a cylinder in Hiemenz flow. Derivation of the Homann similarity solution has the almost same procedure as well but it uses cylindrical coordinates instead of Cartesian coordinates. The velocity components of the flow can be shown as:

(117)$\begin{array}{cc}\hfill v_r=& v_r(r,z)\hfill \end{array}$

(118)$\begin{array}{cc}\hfill v_{\theta }=& 0\hfill \end{array}$

(119)$\begin{array}{cc}\hfill v_z=& v_z(r,z)\hfill \end{array}$

(120)$\begin{array}{cc}\hfill p=& p(r,z).\hfill \end{array}$

As axisymmetric assumptions, derivative with respect to $\theta$ ($\frac{\partial }{\partial \theta }$) and $v_{\theta }$ can be assumed as zero. The Navier–Stokes equations in cylindrical coordinates can be written as:

(121)$\begin{array}{cc}\hfill \frac{\partial v_r}{\partial r}+\frac{v_r}{r}+\frac{\partial v_z}{\partial z}=& 0\hfill \end{array}$

(122)$\begin{array}{cc}\hfill v_r\frac{\partial v_r}{\partial r}+v_z\frac{\partial v_z}{\partial z}=& -\frac{1}{\rho }\frac{\partial p}{\partial r}+\nu \left(\frac{{\partial }^2{v}_r}{\partial r^2}+\frac{1}{r}\frac{\partial v_r}{\partial r}-\frac{v_r}{r^2}+\frac{{\partial }^2{v}_r}{\partial z^2}\right)\hfill \end{array}$

(123)$\begin{array}{cc}\hfill v_r\frac{\partial v_z}{\partial r}+v_z\frac{\partial v_z}{\partial z}=& -\frac{1}{\rho }\frac{\partial p}{\partial z}+\nu \left(\frac{{\partial }^2{v}_z}{\partial r^2}+\frac{1}{r}\frac{\partial v_z}{\partial r}+\frac{{\partial }^2{v}_z}{\partial z^2}\right).\hfill \end{array}$

If the same procedures applied to Navier–Stokes in Cartesian coordinates for Hiemenz flow, are applied in cylindrical coordinates, the final potential function and stream function will be:

(124)$\begin{array}{cc}\hfill \varphi (r,z)=& k(\frac{r^2}{2}-z^2)\hfill \end{array}$

(125)$\begin{array}{cc}\hfill \psi (r,z)=& kz{r}^2,\hfill \end{array}$

where k is a constant which depends on the freestream flow and the body shape. The corresponding velocities are:

(126)$\begin{array}{cc}\hfill v_r=& \frac{\partial \varphi }{\partial r}=\frac{1}{r}\frac{\partial \psi }{\partial z}=kr\hfill \end{array}$

(127)$\begin{array}{cc}\hfill v_z=& \frac{\partial \varphi }{\partial z}=\frac{-1}{r}\frac{\partial \psi }{\partial r}=-2kz.\hfill \end{array}$

Once velocities are calculated, the pressure from the Bernoulli can be calculated as:

(128)$p_0-p=\frac{1}{2}\rho k^2(r^2+4z^2).$

The viscous velocity components and pressure can be calculated as:

(129)$\begin{array}{cc}\hfill v_r=& rg\left(z\right)\hfill \end{array}$

(130)$\begin{array}{cc}\hfill v_z=& -2f\left(z\right)\hfill \end{array}$

(131)$\begin{array}{cc}\hfill p_0-p=& \frac{1}{2}\rho k^2(r^2+F\left(z\right)).\hfill \end{array}$

Listing 4. Implementation of Thomas algorithm for Hiemenz profile.

After the affine transformation, the final ordinary differential equation will be:

(132)${\phi }^{\prime \prime \prime }+2\phi {\phi }^{\prime \prime }-{\phi }^{\prime 2}+1=0,$

where the boundary conditions of the equation is:

(133)$\begin{array}{c}\hfill {\phi }^{\prime }\left(0\right)=0\end{array}$

(134)$\begin{array}{c}\hfill \phi \left(0\right)=0\end{array}$

(135)$\begin{array}{c}\hfill {\phi }^{\prime }(\infty )=1\end{array}$

If one compares Hiemenz and Homann similarity solutions (see Equations (98) and (132)) the only difference is the 2 in the second term. In the same manner, the computational process will be the same except for two lines of code.

2.2.4. Numerical Solution of Homann Flow Problem

Since equations of Hiemenz and Homann similarity solutions are the same except one coefficient, procedures for the solution are also the same, except for two lines of code. If it is assumed that ${\phi }^{\prime }=h$ and $\phi =p$ in Equation (132), the system of equations will be:

(136)$\begin{array}{cc}\hfill h^{\prime \prime }+2p{h}^{\prime }-h^2+1=& 0\hfill \end{array}$

(137)$\begin{array}{cc}\hfill p^{\prime }-h=& 0,\hfill \end{array}$

where the boundary conditions of the system are:

(138)$\begin{array}{c}\hfill p\left(0\right)=0\end{array}$

(139)$\begin{array}{c}\hfill h\left(0\right)=0\end{array}$

(140)$\begin{array}{c}\hfill h(\infty )=1\end{array}$

If the finite difference scheme is applied to the h variable, the final system of equations will be:

(141)$A_n{h}_{n+1}+B_n{h}_n+C_n{h}_{n-1}+D_n=0,$

where the $A_n$, $B_n$, $C_n$, and $D_n$ are:

(142)$\begin{array}{cc}\hfill A_n=& \frac{1}{\Delta {\eta }^2}+\frac{p_n}{\Delta \eta }\hfill \end{array}$

(143)$\begin{array}{cc}\hfill B_n=& \frac{-2}{\Delta {\eta }^2}-h_n\hfill \end{array}$

(144)$\begin{array}{cc}\hfill C_n=& \frac{1}{\Delta {\eta }^2}-\frac{p_n}{\Delta \eta }\hfill \end{array}$

(145)$\begin{array}{cc}\hfill D_n=& 1.\hfill \end{array}$

The final $G_n$ and $H_n$ coefficients are the same as Hiemenz flow as well and they are:

(146)$\begin{array}{cc}\hfill G_n=& \frac{-C_n{G}_{n-1}+D_n}{B_n+C_n{H}_{n-1}}\hfill \end{array}$

(147)$\begin{array}{cc}\hfill H_n=& \frac{-A_n}{B_n+C_n{H}_{n-1}}.\hfill \end{array}$

The final code can be seen in Listing 5. If the Hiemenz code (4) and Homann code (5) are compared, the only difference is in the A and C and it is because of the finite difference approach for the $h^{\prime }$.

Listing 5. Implementation of Thomas algorithm for Homann profile.

The final solution profile can be seen in Figure 8. The results are also validated with White [44]. One can use these results to validate their own codes. Reference data will be shared on GitHub as well.

3. Discussion

Fast computational fluid dynamics solvers are crucial for engineers because the design process requires a lot of simulations to reach the final and optimized design. One of the most important factors that make a solver fast is the language itself. Writing almost the same script in different languages may give different solution central processing unit (CPU) times. For instance, a script that uses long and complex for loops in both Julia and Python environments will result in different solution times because Python is slower with for loops in general. The reasons that cause this slow behavior will not be covered in this paper since it is out of the scope of this paper. However, it is important to state these differences in order to decide which coding language is proper for the simulation that is planned.

Sometimes, some other criteria, such as a user-friendly environment, might be the critical condition. The easy matrix–matrix multiplication and backslash linear system solver can be some of the user-friendly examples. In CFD, the usage of vectors and matrices is so common. Most of the time, it is required to multiply or add vectors or matrices with one another. Fortran, which is one of the fastest languages and also one of the fundamental languages in the CFD industry, does not have a built-in element-wise vector or matrix multiplication feature. This requires the use of for loops for each element-wise matrix and vector operation. As a result of this, every time, one needs to use a for loop to do these operations. This ends up with hard-to-read codes and also excessive usage of indices is another source of possible mistakes that will cause trouble during the debugging. However, in the Julia environment, it can be done in one line without any trouble. Listing 6 is showing the two for loop usage. Both of them use Julia syntax in order to prevent confusion; however, it is required to state that Fortran has a different syntax than this but, logically, the for loops are the same with Fortran logic.

Another user-friendly feature of Julia over Fortran is to plot a vector in the code with built-in functions. However, in Fortran, it is not possible to do it with built-in functions so one needs to extract the solution vector or matrix to an external file to visualize it. As it is stated before, coding language selection can be an important topic. In this section, the strong sides of the Julia environment will be stated. One of the user-friendly features of Julia is compact for loop syntax. Listing 7 shows the traditional way and compact way to use a for loop in the Julia environment.

Listing 6. Comparison of the element-wise and compact matrix-matrix multiplication.

Listing 7. Comparison of the traditional and compact way of for loop usage.

The readability of the code is also extremely important in order to explain the code to others or provide it as open-source code. Julia is so strong in this topic because it allows the usage of LATEX language within the code. It is one of the unique abilities of this language. If you are writing a formula that is full of Greek letters and regular letters, it is inevitable to use complex variable names. Comparison of the usage of Greek letters and regular formula writing can be seen in Listing 8.

Listing 8. Comparison of the Greek letter usage and conventional formula writing in Julia.

The last advantage of Julia stated in this present paper is the ability to do both dynamic and static RAM allocations. It may be important to manually allocate the variable sizes and types for optimized code. However, sometimes it is easier to use just variables without any initializing. Initializing differences in the Julia environment can be seen in Listing 9.

Listing 9. Different initializing methods in Julia environment.

As it is seen from the example scripts, Julia is a user-friendly, fast, open-source, and free language which can increase productivity drastically [36]. It is a great choice both for those who are new to coding and coding experts. Julia also can call the C, Fortran, and Python libraries so it is great for experienced engineers who think their previous code in other coding languages will be useless.

4. Conclusions

In the computational fluid dynamics industry, it is crucial to have some predictions about the flow that will be simulated. It helps to spot the location at which finer mesh is required. Fundamental knowledge about canonical flows is crucial in this point because most of the complex flow consists of a combination of a couple of canonical flows. For example, an airfoil CFD simulation can consist of boundary-layer flow over the smooth part of it, mixing layer flow where the tail ends, and blunt body flow in the wake region. In other words, one airfoil simulation consists of three different canonical flows. A full understanding on canonical flows is extremely important in order to simulate similar flows accurately and cheaply. In this paper, boundary-layer theory is introduced and boundary-layer flows are derived from scratch. The Blasius flow, Hiemenz flow, Homann flow, Falkner–Skan flow are the focus of this paper. Once the derivations of them are completed, derived forms are implemented in the Julia environment. In order to model the equations, a finite difference scheme for space discretization is used and the Thomas algorithm is used for the linear system solution. It has to be noted that other methods such as Runge-Kutta, Runge-Kutta-Fehlberg, compact finite difference, high-order finite difference can also be used to solve the given ODEs; however, in this present paper, authors preferred the finite difference method and Thomas algorithm. The authors further discussed the advantages of the Julia language over some other coding languages available in the literature. It is shown that Julia syntax is so straightforward, easy, and user-friendly. Strong sides of the Julia environment are stated with the given comparisons as well. The popularity of Julia may drastically increase in the near future because of its potential.

Author Contributions

Code generation, F.O. and K.K.; validation, F.O. and K.K.; writing—original draft preparation, F.O.; writing—review and editing, F.O. and K.K.; visualization, F.O. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

All the data used and generated in this study is available in the GitHub link provided in Appendix A.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

Julia setup files can be downloaded from their website (https://julialang.org/downloads/ (accessed on 1 June 2021)). The website also includes instructions on how to install Julia on Windows, Linux, and mac operating systems. Some of the useful resources for learning Julia are listed below:

• https://docs.julialang.org/en/v1/ (accessed on 1 June 2021)

• https://www.coursera.org/learn/julia-programming (accessed on 1 June 2021)

• https://www.youtube.com/user/JuliaLanguage/featured (accessed on 1 June 2021)

• https://www.youtube.com/user/ParallelComputingandScientificMachineLearning (accessed on 1 June 2021)

• https://discourse.julialang.org/ (accessed on 1 June 2021)

It is common to use external packages for Julia. In order to do that, Pkg, which is Julia’s built-in package manager, can be used. Once Julia is opened, Pkg can be activated with the “]” button in Windows. In Linux, calling “julia” in the terminal will open it. After that “Pkg.add(“Pluto”)” will trigger the setup process for that package. In here, we used Pluto as an example because, in GitHub, our codes are developed in the Pluto environment. After Pluto is installed. Pluto can be run with “Pluto.run()”. This command will open a new tab in the browser which you can run your Julia codes. After that, the “using Pluto” line must be placed to the top of the file. For “Plots” package, the commands will be “Pkg.add(“Plots”)” and “using Plots”. Since the Plots package does not have a GUI, there is not a command called “Plots.run()”.

Other than Pluto, JuliaPro which includes Julia and the Juno IDE (https://juliacomputing.com/products/juliapro/ (accessed on 1 June 2021)) can be used as an editor and compiler. This software contains a set of packages for plotting, optimization, machine learning, database, and much more. Pluto is appropriate for small scripts while JuliaPro is better for more complex codes. The GitHub link of the codes used in this paper is:

• https://github.com/frkanz/A-CFD-Tutorial-in-Julia (accessed on 1 June 2021)

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures

Figure 1. Schematic description of the flow over a flat plate.

Figure 2. Schematic description of the change in the velocity profile from initial guess to final solution.

Figure 3. The velocity distribution of the Blasius similarity solution obtained by given Julia code and data digitized from Schlichting [6].

Figure 4. The representation of flow types that can be calculated with the Falkner–Skan equation. (a) Blasius flow (m=0), (b) wedge flow (0<m<1), (c) Hiemenz flow (m=1), and (d) corner flow (1<m<2).

Figure 5. The velocity distribution of the Falkner–Skan similarity solution obtained by Julia code and data digitized from Schlichting [6] for varying m values.

Figure 6. Schematic description of the Hiemenz and Homann flow. Left side corresponds to general flow and the Right side is the extended vision of the dashed rectangle. The flow represents Hiemenz flow if the circle is the projection of a cylinder and Homann flow if the circle is the projection of a sphere.

Figure 7. The velocity distribution of Hiemenz similarity solution obtained by Julia code and data digitized from White [44].

Figure 8. The velocity distribution of Homann similarity solution obtained by Julia code and digitized from White [44].