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

: A boundary-layer is a thin ﬂuid layer near a solid surface, and viscous effects dominate it. The laminar boundary-layer calculations appear in many aerodynamics problems, including skin friction drag, ﬂow separation, and aerodynamic heating. A student must understand the ﬂow physics and the numerical implementation to conduct successful simulations in advanced undergraduate-and graduate-level ﬂuid dynamics/aerodynamics courses. Numerical simulations require writing computer codes. Therefore, choosing a fast and user-friendly programming language is essential to reduce code development and simulation times. Julia is a new programming language that combines performance and productivity. The present study derived the compressible Blasius equations from Navier–Stokes equations and numerically solved the resulting equations using the Julia programming language. The fourth-order Runge–Kutta method is used for the numerical discretization, and Newton’s iteration method is employed to calculate the missing boundary condition. In addition, Burgers’, heat, and compressible Blasius equations are solved both in Julia and MATLAB. The runtime comparison showed that Julia with for loops is 2.5 to 120 times faster than MATLAB. We also released the Julia codes on our GitHub page to shorten the learning curve for interested readers.


Introduction
Until the 19th century, scientists neglected the effects of viscosity in their hydrodynamic and aerodynamic calculations using potential flow theory. However, this assumption led to a contradiction between theoretical predictions and experimental measurements of drag force acting on a moving body, now known as the d'Alembert paradox [1]. Later, a revolutionary boundary-layer concept is introduced [2,3]. In this concept, the fluid flow over a surface is divided into two regions by the boundary-layer edge: an area between the surface and the boundary-layer edge dominated by the viscous effects and a region outside the boundary-layer edge where the viscous effects can be neglected. It enables a significant simplification of full Navier-Stokes equations.
The boundary-layer theory was first presented by Prandtl [4] in 1904, and it provides the solutions of velocity and temperature profiles within the boundary-layer by using approximations. One can obtain Blasius [5,6], Falkner-Skan [7], and compressible Falkner-Skan [2,8] solutions by using this approach. Researchers extensively use these solutions to validate the computational fluid dynamics (CFD) simulations. Moreover, in a CFD simulation, one can calculate the boundary-layer thickness in advance to estimate the required grid parameters to resolve the boundary-layer region. Understanding the fundamentals of boundary-layer theory is critical for engineers to solve today's aerodynamic design challenges.
It may be challenging to fully understand the fundamentals of the boundary-layer theory in undergraduate-and graduate-level boundary-layer courses. Most of the time, books skip or briefly mention some steps in the derivation of a system of equations. The visual abstract of the present paper which is mainly designed around 3 major points. Major points are further divided into smaller points which correspond to purposes/ideas of the particular major point. Each major point is connected to one another, which makes them complete.

Compressible Laminar Boundary-Layer
Compressibility effect in the boundary-layer requires additional calculations. Constant density assumption in incompressible speeds is no longer valid for the compressible boundary-layer. In compressible speeds, temperature and density change within the boundary-layer. It is crucial to capture the velocity, temperature and density variations in the boundary-layer to obtain accurate simulation results. One can estimate the number of element required to resolve the boundary-layer in the CFD simulation by using the boundary-layer theory. Compressible Blasius is also widely used for CFD validations in high-speed flows. In this section, compressible Blasius equations will be derived from scratch and implemented in the Julia environment. The contribution of this paper is employing the Julia language. The equations used in this paper are already in the literature [2,44]. The manuscript may enable students to adopt the programming language with easily and available GitHub codes, which may shorten the learning curve.

Compressible Blasius Equations
Incompressible Blasius solution is a similarity solution for a flat plate. The assumptions for the incompressible Blasius equations are given in our previous work [19]; interested readers can check the details from there. In the compressible region, the temperature effects must be taken into account for an accurate solution. In the incompressible region, the temperature and density changes are small enough to be neglected. In the compressible region, the temperature can increase drastically as a result; density decreases within the boundary-layer. For example, the temperature on the solid wall can reach 7 times the freestream temperature in Mach 6 flow over a wedge. If the freestream temperature is 300 K, the wall temperature will be around 2100 K. In order to compare the quantity, the melting point of titanium is around 1941 K [45]. This problem is still a challenge for aerospace applications in which high Mach numbers are involved.
The compressible Blasius equations can be derived from the compressible Navier-Stokes equations, which can be expressed in two spatial dimensions as: where ρ is the density, u and v are the velocities in xand ydirections, p is the pressure, µ is the dynamic viscosity, λ is the second viscosity coefficient, k is the thermal conductivity, T is the temperature, c p is the specific heat at constant pressure, and Φ is the dissipation function, which can be written as: In order to obtain the boundary-layer equations, dimensional analysis is required to neglect the variables that have smaller orders than others. The flat plate boundary-layer development is illustrated in Figure 2. In this flow, u velocity is related to freestream velocity and the order of magnitude is one. The x is related to plate length, so its order of magnitude is also one. The y distance is related to boundary-layer thickness, so it is in the order of δ which is the boundary-layer thickness. The density, ρ, is related to freestream density so its order of magnitude is also one. The magnitude of the v velocity can be calculated from the continuity equation, Equation (1). In order to get zero from this equation, all variables must be in the same order so v is in the order of δ as a result of this, ∂(ρv) ∂y = O(1). When the magnitude analysis is completed in the same manner, the boundary-layer equations can be obtained. It has to be noted that dynamic viscosity is in the order of δ 2 , pressure and temperature are in the order of one. The specific heat at constant pressure is in the order of one. The second viscosity coefficient, λ, can be taken as −2/3µ because of Stokes' hypothesis. Once the order of magnitude is obtained for each of the terms, some of the terms can be neglected because δ 1. The final system of equations in steady-state condition ( ∂ ∂t = 0) will be: The boundary-layer velocity profile is illustrated with a blue line. The black dot corresponds to the boundary-layer edge at that station. The density, temperature, and velocity at the boundary-layer edge are ρ e , T e , and u e , respectively. The boundary-layer thickness is defined with δ(x), which is the function of x.
Equation (7) can be expressed at the boundary-layer edge as: The variables are changing from the solid surface up to the boundary-layer edge. At the boundary-layer edge, they reach to freestream value for the corresponding variable and remain constant. The velocity change in the y-direction at the boundary-layer edge is zero ( ∂u ∂y | y=δ = 0), because it is constant at boundary-layer edge. Equation (8) indicates that the pressure gradient in the y-direction is zero, so pressure at the boundary-layer edge equals the pressure within the boundary-layer (p e = p). Equation (10) becomes: The velocity at the boundary-layer edge is equal to freestream velocity, which is constant in x-direction for a flat plate. In other words, edge velocity gradient in x-direction is zero ( ∂u e ∂x = 0). If the equation of state is used to obtain the ratio of density and temperature as: where R is the gas constant. It is known that p = p e , so ρT = ρ e T e . The final system of equations is: where T e T = ρ ρ e . At this point, a similarity parameter can be introduced to the system to obtain a similarity solution [46]. The similarity parameter, η, can be defined as: where s = µ e ρ e u e x. Let's assume that the stream function is The u and v velocities can be calculated from the stream function as: (20) In this step, the variables in Equation (15), u, v, ∂u ∂x , ∂u ∂y , and ∂ ∂y µ ∂u ∂y can be calculated. The first derivative of η with respect to y and the first derivative of s with respect to x will be required for the chain rule.
It is better to note that, in ∂η ∂y calculation, T e T = ρ ρ e relation is used. The u velocity can be calculated as: The same procedure can be applied for v velocity as: Once u velocity is obtained, the derivatives with respect to x and y can be calculated as: ∂u ∂x = ∂u ∂η ∂η ∂x (38) = ∂(u e f ) ∂η ∂η ∂x (39) =u e f ∂η ∂x . (40) All terms in Equation (15) are known. If the above terms are substituted into Equation (15) and the necessary simplifications are done, the final equation will be: It has to be noted that if ρ = ρ e and µ = µ e , in other words, if the flow is incompressible, Equation (41) becomes an incompressible Blasius equation ( f + f f = 0). Equation (41) can be further simplified as: = T e u 2 e ρ 2s When these terms are substituted into Equation (17), the new equation will be: Equation (51) where c p = γ γ−1 R, M = u e a e , and a e = √ γRT e . In the final system of equations, theμ can be calculated from Sutherland Viscosity Law [47]. The dimensional viscosity function is: where The derivative of the viscosity is also required. The derivative terms can be calculated as: The final system of equations is: It has to be emphasized thatρ is a function of η and the final system of equations is coupled, so they have to be solved together. The boundary conditions of the system for an adiabatic system are: The boundary condition for the isothermal wall depends on the wall temperature. For example, if the wall temperature equals the boundary-layer edge temperature, it will beρ = 1, and it will be replaced with the last boundary condition of the system. In the adiabatic boundary condition, the derivative of the temperature with respect to wall-normal direction will be 0. During the numerical procedures, the difference will be emphasized one more time.

Numerical Procedure
In this section, the compressible Blasius equation will be solved with the fourth-order Runge-Kutta method [48] and Newton's iteration method [49]. Different methods can be used for this problem; however, we used Runge-Kutta and Newton's method because of their extensive usage in the literature and accuracy. To start the numerical procedure, high-order differential equations can be reduced to the first-order differential equations as: if Equations (64)-(68) are substituted into Equations (58) and (59), the final version of these equations can be written as: The final system of equations can be written in the matrix form as: The adiabatic boundary conditions for the system are: The isothermal boundary conditions for the system are: The functions can be introduced in Julia as shown in Listing 1, where cµ is the second coefficient of the Sutherland Viscosity Law, T∞ is the temperature at the boundary-layer edge, M∞ is the Mach number at the boundary-layer edge, γ is the specific heat ratio, Pr is the Prandtl number and y 1 , y 2 , y 3 , y 4 , and y 5 are the terms given in Equations (64)-(66), Equation (67), and Equation (68). In the functions given in Listing 1, only 2 parameters are dimensional, which are cµ and T∞. In this tutorial paper, Kelvin is the unit of both parameters. If the temperature unit is required to be different, such as Fahrenheit or Rankine, the units of cµ and T∞ must be transformed into the new unit accordingly. function Y3 (y 1 , y 3 , y 4 , y 5 , cµ, T∞) 10 return −y 3 * ((y 5 /(2 * (y 4 ))) − (y 5 /(y 4 + cµ/T∞))) − y 1 * y 3 * ((y 4 + cµ/T∞)/(sqrt(y 4 ) * (1 + cµ/T∞))) function Y5 (y 1 , y 3 , y 4 , y 5 , cµ, T∞, M∞, Pr, γ) 18 return −y 5^2 *((0.5/y 4 ) − (1/(y 4 + cµ/T∞))) − Pr * y 1 * y 5 /sqrt(y 4 ) * 19 (y 4 + cµ/T∞) In this paper, implementation of the Runge-Kutta method will be provided. The derivation of the Runge-Kutta method and how it calculates the function value at the next step can be checked from Reference [49]. The implementation of the Runge-Kutta method for the compressible Blasius problem can be seen in Listing 2, where N is the number of elements. It has to be emphasized that the number of node points is N + 1, which means that terms must be calculated until (N + 1) th node. The first point is the boundary condition, so there will be N number of calculations.
The initialization and the boundary conditions can be introduced as shown in Listing 3, where adi is a flag for the adiabatic or isothermal condition selection and Tw is the dimensionless wall temperature. It is nondimensionalized with T e , so if the temperature at the boundarylayer edge, T e , is 300 K and the wall temperature is required to be 150 K, Tw must be entered as 0.5. Another important point about Listing 3 is the indices. In the derived formulations, indices start from 0. However, in both Julia and MATLAB, indices start from 1. This is the reason why indices are starting from 1 in Listing 3 boundary conditions part. Listing 2. Implementation of Runge-Kutta method in Julia environment. It requires four slope calculation to estimate the function value in the next node value. § ¤ 1 function RK ( N, ∆η, y 1 , y 2 , y 3 , y 4 , y 5 , cµ, T∞, Pr, γ, M∞) , cµ, T∞, M∞, Pr, γ) 9 10 # Second Step , cµ, T∞, M∞, Pr, γ) 16

# Third
Step return y 1 , y 2 , y 3 , y 4 , y 5 39 end ¦ ¥ In the system of equations, there are five equations and five boundary conditions; however, two boundary conditions are located at the end of the domain. In order to start the calculation, all values at the η = 0 should be given. α 0 and β 0 in the Listing 3 are the initial guesses for the missing boundary conditions. They can be any value. Once they are introduced to the system, compressible Blasius equations can be solved. When the equations are solved with guessed initial conditions, the solution vector must satisfy the boundary conditions at the end of the domain. However, it will not converge at the first try because the guessed boundary conditions are not correct. To overcome this problem, different methods can be used, such as the shooting method, bisection method, or Newton's iteration method. In this paper, Newton's iteration method will be used because it is fast and it is not hard to implement. In order to use it, the algorithm needs to run with the initial guesses one time. Once the y 2 (corresponds to u) and y 4 (corresponds to T) at the end of the domain are obtained, an arbitrary small number can be added to one of the initial guesses. The algorithm can be run one more time with the new boundary condition guesses. After that, the same small number can be added to the other initial guess and the algorithm can be run one more time. After running the algorithm 3 times, there will be 3 different y 2 -y 4 pairs. It has to be noted that when the small number is added to the second boundary condition (in the third run), other boundary conditions should be equal to the value in the first run. In other words, after adding a small value in the second run, it should be subtracted in the third run. The main purpose of running three times is to determine the more accurate boundary condition guess. The new boundary conditions can be calculated with: where α and β are the initially guessed boundary conditions. dα and dβ are required for the new boundary conditions. These values can be approximated from the Taylor series expansion of the y 2 and y 4 , which can be shown as: The partial differentials can be approximated with the finite difference as: where y 2,old and y 4,old are the values obtained from the first run, y 2,new,1 and y 4,new,1 are the values obtained from the second run, and y 2,new,2 and y 4,new,2 are the values obtained from the third run. Once everything is calculated, the system of equations in Equation (86) can be used to calculate dα and dβ. The implementation of the explained method in Julia can be seen in Listing 4.

¦ ¥
The same procedure will run until y 2 , and y 4 at the end of the domain will be 1. It is important to decide the upper limit of the domain. If it is small, it will force the value at that point to be 1 where it should not be. It is also important to choose the small number, ∆, smaller than convergence criteria which will finalize the simulation. If ∆ is higher than the convergence criteria, the simulation might run until it reaches the maximum iteration number. In the code provided in GitHub, convergence criteria is taken as 1 × 10 −9 and the small number is taken as 1 × 10 −10 . The results of the code for M = 4.5 and M = 2.8 are illustrated in Figure 3, where total temperatures are 311 K for both. The freestream temperature is calculated from isentropic relation and it is 61.584 K for M = 4.5 and 121.11 K for M = 2.8. The results are compared with the Iyer's [20] BL2D boundary-layer solver, which is used in NASA's well-known compressible boundary-layer stability solver LASTRAC [21].

Comparison of Julia and MATLAB
The design process requires lots of simulations in order to obtain the final and optimized design. It is highly beneficial to have a fast CFD solver. One of the crucial factors that affects the speed of the solver is the language. The same script may lead to different central processing unit (CPU) times with different coding languages. Additionally, similar simulations will be required multiple times. Eventually, the total time spent on simulations might be drastic with a slow solver.
MATLAB is one of the languages that is widely used. It is one of the favorite coding language for most of the students because of its user-friendly syntax, easy debugging feature, and built-in functions. One of the most important drawbacks of this language is that it is not free. It is also slower than high-performance languages, such as Fortran and C/C++. Julia is a user-friendly, open-source language that can increase productivity drastically [13]. Another great feature of Julia is that it is completely free. Julia can call C, Fortran, and Python libraries. It is great for experienced engineers who think that their previous code in other coding languages will be useless.
One of the great concerns about language selection is the speed of the code. For large-scaled projects, most of the time, a fast solver is the most important point. In this section, the same problem will be solved with Julia and MATLAB codes. The solution times will be compared with each other. The purpose of this comparison is to provide a rough estimation about code execution speeds. Three different test cases will be applied for both languages. The cases are unsteady inviscid Burgers' equation with the first-order backward finite difference scheme, heat equation with second-order central finite difference scheme, and compressible Blasius equation with fourth-order Runge-Kutta and Newton's iteration method. Burgers' and heat equations will be tested with f or loops with file operations, vectorized operations with file operations, f or loops without file operations, and vectorized operations without file operations. For the compressible Blasius equation solver, the code developed for this paper will be used. The test cases will simulate real-life problems by solving the problem and exporting the solution vector to a text file when the file operations are included. In real-life problems, most of the time, post-processing is required after the simulation. In order to do that, saving data into a file is required. If it is a steady problem, exporting can be done at the end of the simulation, if there are not any other limitations or additional requirements. On the other hand, if the problem is unsteady, exporting the data in different time steps during the simulations is required. This is the reason why there will be two different simulations where data will be exported and will not be exported.
The first case is unsteady, inviscid, Burgers' equation in one dimension, which can be represented in the conservative form: The equation is solved with a first-order backward finite difference scheme. The details of the scheme will not be provided because the purpose of the test case is to measure the speed difference of two similarly developed codes. However, codes that are used in this paper are available on GitHub. Interested readers can check the implementation details from the codes. The number of elements in the problem is taken as 2500, 5000, and 10,000. The same simulation will be run with an increasing number of elements to show the solution time change trend. The execution time will be calculated by BenchmarkTools in Julia and tic/toc functions in MATLAB. The standard deviation will be calculated manually by using 10 data points obtained from the runs. The time step is taken as half of the grid spacing. The domain is limited with [0, π] and the initial conditions for velocity, u(x i , t), are taken as: The solution vector is written to a ".txt" file for every hundredth iteration. The mean execution times with the standard deviation of the data obtained by Julia and MATLAB solvers are given in Tables 1 and 2. Table 1 provides the execution times with file operations, and Table 2 excludes file operations in the calculations. The results show that MATLAB is slow with file operations. There is approximately 15 times' difference between Julia and MATLAB mean execution times with file operations and f or loops, but the speed-up difference is decreasing to 8 with vectorization. Without file operations, Julia is 3 times faster than MATLAB with f or loops. However, MATLAB vectorization is faster than Julia. One interesting point of this test case is that MATLAB execution times are reduced by vectorization, except for the N = 10,000 case with file operations excluded. Detailed investigations about this trend indicate that there is a relation between the L1, L2, L3 cache size of the CPU and the vectorization performance. Tests are completed in 3 different computers with varying cache sizes. After a certain number of elements, vectorization starts to increase the mean execution time. The limiting number of elements is related to cache size. In computers with higher cache sizes, the negative effect of vectorization started after N = 5000. In computers with lower cache sizes, the negative effect started after N = 2500. This trend is not observed in Julia language. In Julia, f or loops are highly optimized and manual vectorization leads to an increase in the mean execution times because manual vectorization creates temporary arrays during the calculations. Creating and deleting these temporary arrays require more time than calculation with f or loops. In MATLAB, results indicate that temporary array usage is faster than f or loops up to certain array size.
where α is a constant which is taken as 0.25∆x. The time step is taken as ∆x. This assures that the coefficient of the second derivative will satisfy the stability condition. The boundary conditions of the system are 1 for each side and the initial conditions for the remaining nodes are 0. The heat equation is solved with the second-order central finite difference with 250 × 250, 500 × 500, and 1000 × 1000 elements. The domain is limited with [0, π] 2 . The execution times of the two codes are given in Table 3 with file operations and in Table 4    Lastly, the derived compressible Blasius equations for the present study will be solved in both MATLAB and Julia. The difference of that case is to test the function calls because sometimes dividing the solver into smaller functions may lead to longer solution times. The problem will be solved with 50,000, 100,000, and 200,000 elements. Table 5 gives the solution times of two codes developed in MATLAB and Julia. In this problem, Julia is drastically faster than MATLAB, and the time differences are increasing with the problem size. With 50,000 elements, Julia is approximately 15 times faster than MATLAB, with 100,000 elements, it is 32 times faster, and with 200,000 elements, it is 120 times faster. Although time differences are varying with problems, Julia with f or loops exhibited better performance than MATLAB in every problem. On the other hand, MATLAB showed better performance when both of the codes are developed in vectorized form. In general, MATLAB file operations are slower than Julia. It has to be noted that MATLAB has special data exporting options which might be faster, such as .mat extensions. In order to conduct an exact comparison, regular .txt extension with conventional exporting commands are used. The main purpose of these time comparisons is to provide an approximate performance differences between Julia and MATLAB under different conditions. In this paper, Julia is compared with MATLAB. Interested readers can check Lubin and Dunning's paper [50] for other coding language comparisons.

Conclusions
Compressible Blasius equation, which comes from boundary-layer theory, is extensively used by researchers to validate the CFD simulation results. One can estimate the number of elements required to capture the boundary-layer by using the solution of the compressible Blasius equation. Although it is crucial to understand the boundary-layer theory, undergraduate-or graduate-level boundary-layer classes may not be adequate for a student to fully understand it due to time limitations. A step-by-step tutorial may help students to understand the theory better. Both compressible and incompressible boundary-layer problems require numerical solution. Deriving the equations from scratch and implementing the numerical methods may shorten the learning curve for a student or an engineer.
In this paper, compressible Blasius equation and energy equation are derived from scratch. The final system of equations is solved in the Julia environment. For the numerical implementation, the fourth-order Runge-Kutta and Newton's iteration methods are employed. It has to be noted that other methods such as Runge-Kutta-Fehlberg, compact finite difference, a high-order finite-difference can also be used to solve the final system of equations. However, authors preferred the Runge-Kutta and Newton's iteration method due to their accuracy and wide usage in the literature. Moreover, the authors compared the Julia and MATLAB solver speed to give an initial impression about performance of Julia. The results showed that MATLAB is slower than Julia in file operations. Additionally, Julia is faster than MATLAB with f or loops. On the other hand, Julia vectorization affects the solution times negatively. However, MATLAB vectorization decreases the solution time for small-sized problems. When the problem size increases, MATLAB vectorization also has a negative effect on the solution time. It has to be noted that these test cases are relatively less demanding cases. In real-life problems, simulations require longer codes with more complex operations. In longer runs, the time difference in between these two languages may increase.