Skip to content
Posted by curiosityFluidsAdmin13 Comments on Numerical Solution of the Compressible Laminar Boundary Layer Equations CFD from Scratch

Numerical Solution of the Compressible Laminar Boundary Layer Equations

In this post I go over the numerical solution to the compressible boundary layer equations. This is very useful when a quick estimate of shear stress, wall heat flux, or boundary layer height if necessary.

The sections of this post are:

  1. Introduction
  2. Compressibility transformation
  3. Using the general parabolic form
  4. Numerical solution using Crank-Nicolson
  5. Results and comparison to OpenFOAM
  6. Conclusions

Introduction

For a compressible flow, the boundary layer equations are:

r_o^j\frac{\partial \rho}{\partial t}+\frac{\rho u r_o^j}{\partial x}+\frac{\rho v r_o^j}{\partial y} = 0                 (Continuity)

\rho\left[\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}\right] = - \frac{\partial p}{\partial x}+\frac{\partial}{\partial y}\left(\mu\frac{\partial u}{\partial y}\right)                 (Momentum)

\rho c_p\left[\frac{\partial T}{\partial t}+u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}\right]=k \frac{\partial}{\partial y}\left(k \frac{\partial T}{\partial y}\right)+\mu \left(\frac{\partial u}{\partial y}\right)^2+u\frac{\partial p}{\partial x}+\frac{\partial p}{\partial t}               (Energy)

Since (in general) we are only concerned with the steady state boundary layer these equations reduce to:

\frac{\rho u r_o^j}{\partial x}+\frac{\rho v r_o^j}{\partial y} = 0                 (Continuity)

\rho\left[u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}\right] = - \frac{\partial p}{\partial x}+\frac{\partial}{\partial y}\left(\mu\frac{\partial u}{\partial y}\right)                 (Momentum)

\rho c_p\left[u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}\right]=k \frac{\partial}{\partial y}\left(k \frac{\partial T}{\partial y}\right)+\mu \left(\frac{\partial u}{\partial y}\right)^2+u\frac{\partial p}{\partial x}               (Energy)

These equations are a parabolic partial differential equation. Essentially what this means for us (the person trying to solve it) is that the boundary layer equations are an initial value problem. So if we know the initial boundary layer profile (for example on a flat plate) then the solution can be achieved by marching forward along the wall boundary. This is what we will do here.

In this post intend to solve the above equations allowing for flow with pressure gradient, heat transfer, and compressibility. 

Compressibility Transformation

Because of the complicated nature of the above equations, it is common practice to use a compressibility transformation. In general, compressibility transformations refer to a change of variables that allow the compressible flow equations to be solved using the same numerical methods and approaches as the incompressible boundary layer equations. In particular these transformation typically correct for the fact that density is now a variable that is directly couple to the pressure and temperature fields through the equation of state, and the temperature dependence of viscosity.

Here we will use the Levy-Lees coordinate transformation. This change of variables is:

d\bar{s}= \rho_eU_e\mu_er_o^{2j}dx       (For the direction tangential to the wall)

d\eta = \frac{\rho_eU_er^j_o}{\sqrt{2\bar{x}}}\frac{\rho}{\rho_e}dy (For the direction normal to the wall)

Additionally we define some new dependent variables:

A scaled velocity parameter: F=u/U_e

A scaled total enthalpy parameter: G=H/H_e

The Chapman-Rubeson constant: C=\frac{\rho\mu}{\rho_e\mu_e}

The Prandtl number: Pr=\frac{\mu c_p}{k}

A dimensionless pressure gradient parameter:\beta=\frac{2\bar{s}}{U_e}\frac{dU_e}{d\bar{s}}

When these variables and the transformations are applied to the boundary layer equations the following equations result:

Continuity:
2\bar{s}F_{\bar{s}}+V^\prime+F=0

Momentum:
2\bar{s}FF_{\bar{s}}+VF^\prime=\beta\left(\frac{\rho_e}{\rho}-F^2 \right)+\left(CF^\prime\right)^\prime

Energy:
2\bar{s}FG_{\bar{s}}+VG^\prime=\frac{C}{Pr}G^{\prime \prime}+\left(\frac{C}{Pr}\right)^\prime G^\prime+\frac{U_e^2}{H_e}\left(\left(1-\frac{1}{Pr}\right)\left(C^\prime F F^\prime+C\left(F^\prime\right)^2+CFF^{\prime\prime}\right)+CFF^\prime\left(1-\frac{1}{Pr}\right)^\prime\right)

In the above equations the primes (‘) denote a derivative in the \eta direction, and the subscript \bar{s} denote a derivative in the transformed streamwise (\bar{s}) direction.

These equations can be found in the EXCELLENT textbook on boundary layer analysis from Schetz. You can find the details of the textbook in the references section of this post.

Using the General Parabolic Form

The transformed equations above are non-linear and quite annoying to solve. However, there is a simple trick that can be performed that makes numerical solution much easier.

The general parabolic PDE for two independent variables is:

A\frac{\partial^2 u}{\partial x^2}+B\frac{\partial^2 u}{\partial x \partial y}+C\frac{\partial^2 u}{\partial y^2}+E\frac{\partial u}{\partial y}+Fu=G

It obvious from inspection that our new continuity, momentum and energy equations ARE NOT in the general parabolic form. However, we can make an approximation by lumping some of the terms into the coefficients of the terms of the general parabolic equations.

So for example the momentum equation becomes:

F\prime\prime+A_1F\prime+B_1F +C_1+D_1F_{\bar{s}}=0

Where:

A_1=-\left(\frac{V}{C}-\frac{1}{C}\frac{\partial C}{\partial \eta}\right)

B_1=0

C_1=\frac{\beta}{C}\left(\frac{\rho_e}{\rho}-F^2\right)

D_1=-\frac{2\bar{s}F}{c}

Similarly for the energy equation:

G\prime\prime+A_2G\prime+B_2G+C_2+D_2G_{\bar{s}}=0

Where:

A_2=-\frac{Pr}{C}\left(V-\frac{\partial}{\partial \eta}\left(\frac{C}{Pr}\right)\right)

B_2=0

C_2=\frac{Pr}{C}\frac{U_e^2}{H_e}\left[\left(1-\frac{1}{Pr}\right)\left(\frac{\partial C}{\partial \eta}F\frac{\partial F}{\partial \eta}+C\left(\frac{\partial F}{\partial \eta}\right)^2+CF\frac{\partial^2F}{\partial \eta^2}\right)+CF\frac{\partial F}{\partial \eta}\frac{\partial }{\partial \eta}\left(1-\frac{1}{Pr}\right)\right]

D_2=-\frac{2\bar{s}F Pr}{c}

We now have two equations written in the general parabolic form which can be easily solved numerically (explained in the next section). Supplementing these two equations is the continuity equation which can simply be integrated at each streamwise location:

dV=-\left(2 \bar{s} F_{\bar{s}}+F\right)d\eta

We now have 3 equations that can be solved for the 3 unknowns: F, G, and V. The boundary conditions for these equations are:

F(\eta=1)=1, G(\eta=1)=1, F(\eta=0)=0 ,G(\eta=0)=G_w and  V(\eta=0)=V_w

Solution of the Equations using Crank-Nicolson

The equations above are now in a form that can be solved in many different ways depending on your preferred flavor of PDE solution. In this post I am going to use my preferred method for this problem which is the Crank-Nicolson method.

Recall that in fact the boundary layer problem is an initial value problem.  To solve this we start at a known station on the boundary and march the solution downstream. This can be done explicitly, or implicitly.  When we solve explicitly, this means that all of the derivatives and coefficients are calculated at the present location and used to predict the next step downstream. An implicit solution means that all of the derivatives are evaluated both at the current step (therefore as knowns) and at the next downstream step (therefore remain unknowns). This creates a matrix which must then be solved. The method we will use here (Crank-Nicolson) is an implicit method.

Additionally, we are going to use an uncoupled approach. This means that the momentum, energy and continuity equations are solved in series, not simultaneously. In a coupled approach these equations are solved simultaneously in one large solution at each step.

Since the momentum and energy equations are in the identical parabolic form, I will just do a general proof for the solution procedure using a dummy variable, W.

Recall that the general parabolic equation is:

W\prime\prime+AW\prime+BW +C+DW_{\bar{s}}=0

Now, the Crank-Nicolson algorithm is as follows. Given an equation of the form (using our variables here):

\frac{\partial W}{\partial \bar{s}}=F\left(W,\eta,\bar{s},\frac{\partial W}{\partial \eta},\frac{\partial^2 W}{\partial \eta^2}\right)

The solution is obtained by:

\frac{W_i^{n+1}-W_i^{n}}{\Delta\bar{s}}=\frac{1}{2}\left(F\left(W,\eta,\bar{s},\frac{\partial W}{\partial \eta},\frac{\partial^2 W}{\partial \eta^2}\right)_i^{n+1}+F\left(W,\eta,\bar{s},\frac{\partial W}{\partial \eta},\frac{\partial^2 W}{\partial \eta^2}\right)_i^{n+1}\right)

Index n, refers to the streamwise location, and index i refers to the normal location. This is demonstrated in the figure below.

BoundaryLayerFD

In the equation, all derivatives are calculated using finite difference approximations. Now we will apply this to our general parabolic PDE. To start, we rearrange the equation to put the streamwise derivative on the LHS of the equation:

-DW_{\bar{s}}=W\prime\prime+AW\prime+BW +C

Insert the 2nd order central difference approximations for first and second normal derivatives, and use Crank-Nicolson for the streamwise derivative.

-D_i\frac{W_i^{n+1}-W_i^n}{\Delta \bar{s}}=0.5\left(\frac{W_{i+1}^n-2W_i^n+W_{i-1}^n}{\Delta \eta^2}+\frac{A_i}{2\Delta\eta}\left(W_{i+1}^n-W_{i-1}^n\right)+C_i+\frac{W_{i+1}^{n+1}-2W_i^{n+1}+W_{i-1}^{n+1}}{\Delta \eta^2}+\frac{A_i}{2\Delta\eta}\left(W_{i+1}^{n+1}-W_{i-1}^{n+1}\right)+C^{n}\right)

Next we place all of the knowns (index n) on the RHS, and all of the unkowns (index n+1) on the LHS, and collect like terms. This yields:

\left(\frac{A_i}{4\Delta\eta}-\frac{1}{2\Delta\eta^2}\right)W_{i-1}^{n+1}+\left(\frac{-D_i}{\Delta\bar{s}}+\frac{1}{\Delta\eta^2}\right)W_i^{n+1}+\left(\frac{-1}{2\Delta\eta^2}+\frac{-A}{4\Delta\eta}\right)W_{i-1}^{n+1} =\left(\frac{A_i}{4\Delta\eta}-\frac{1}{2\Delta\eta^2}\right)W_{i-1}^{n}+\left(\frac{-D_i}{\Delta\bar{s}}+\frac{1}{\Delta\eta^2}\right)W_i^{n}+\left(\frac{-1}{2\Delta\eta^2}+\frac{-A}{4\Delta\eta}\right)W_{i-1}^{n}+A_3^n

The above equation creates a tri-diagonal matrix that is solved at each step. Notice that in order to solve the equations this way we do not apply the Crank-Nicolson algorithm to the coefficients A, C, or D. This is because it would create a non-linear matrix solution. Because of this, the present approach is infact only semi-implicit. 

A figure showing the construction of the tridiagonal matrix formed from this equation is shown below:

tridiagonal

Solution Algorithm

Now that we understand the Crank-Nicolson method and how it can be used to formulate a tridiagonal matrix we can formulate our solution algorithm. For simplicity I have mapped out the algorithm in the figure below:

BLAlgorithm

As I said before, this algorithm is a segregated (or uncoupled), semi-implicit solution to the boundary layer equations.

I should point out that the algorithm shown in the figure above can easily be improved upon. For instance, the coefficients A thru D can be evaluated as an average between the current step at the next step. This is an application of the trapezoidal rule and the solution then becomes iterative. In these iterations as well, one can easily add an update step to account for changing thermophysical properties such as viscosity, thermal conductivity and specific heats that depend on the solution!

All of these improvements are in fact included in the code used for comparison in the next section. This is because we are comparing the solver to a simulation in rhoCentralFOAM that also accounts for changing thermophysical properties.

Results and Comparison to rhoCentralFoam

As a test case, we are going to solve a laminar boundary layer behind an oblique shock on a flat.. I am not going to go into very much detail for the oblique shock case set-up since it is so simple, and I have done one very similar here.

A schematic of the simulation is shown here:

LaminarBLObliqueSetUp
Fig: Simulation Set-up

The above simulation set-up is a very simple way achieve an oblique shock solution. The boundary conditions for this simulation were:

  • P∞=101 Pa, T∞=298 K, U∞=[650.391, -236.723, 0] (Approximately Mach 2)
  • Wall temperature was: Tw=283.15
  • With a plate length (L) of 0.1 m these properties correspond to a free-stream Reynolds number of 4452.
  • Sutherland’s law was used for viscosity
  • NASA Janaf polynomials were used for specific heat capacity and enthalpy

Here’s a video of the simulation result:

animation
Fig: Simulation video (20 degree, Mach 2, Oblique Shock)

Now to use our compressible boundary layer solver we must recognize that the external properties of the boundary layer will be the “post-shock” conditions. So this can be calculated either using the known oblique shock relations, or by extracting them from the OpenFOAM simulation (I did the latter).

To compare the solutions I extracted the wall heat flux, and the wall shear stress from the simulation. This is done using the OpenFOAM utilities wallHeatFlux, and wallShearStress.

rhoCentralheatFluxComparison
Fig: Predicted Wall Heat Flux Comparison
rhoCentralShearStressComparison
Fig: Predicted Shear Stress Comparison

We can see that our solution is pretty good! Both the shear stress and the wall heat flux run right on top of each other! Obviously for this problem the solver outlined in this post is adequate.

Conclusions

In this post I outlined how one can write a numerical solver for the compressible laminar boundary layer equations. I then compared it to a solution achieved using OpenFOAM’s compressible flow solver rhoCentralFOAM and achieved very close results.

Please let me know if I made any mistakes! There were a lot of equations to transcribe from my notebook so I’m sure I messed some up.

If you want more detail for some parts I will gladly add to this post upon request! I skipped over some steps to keep the length of this post down. For instance I didn’t show how the algorithm changes when using averaged coefficients.

Some Useful References

[1] Schetz, J. and Bowersox, R. (2011). Boundary layer analysis. Reston, VA: American Institute of Aeronautics and Astronautics.

Categories

CFD from Scratch

Tags

, , , , , , , , ,

curiosityFluidsAdmin1 View All

3 thoughts on “Numerical Solution of the Compressible Laminar Boundary Layer Equations” Leave a comment ›

  1. Hello. This is very useful. Thank you. Would you be able to show how you could extend this method to three-dimensional boundary layers on an infinite swept flat plate, i.e. 3D compressible boundary layer equations but with d/dz derivative is zero? Thanks

    Reply

    • Hi Anthony, I’m not sure how you’d achieve that. It’s been so long since I wrote this code, i’m a little rusty. When I wrote it though I only ever used it to solve 2D or axisymmetric flow problems. At one point I briefly looked into using an inviscid code to get 3D surface streamlines and then use the streamline data as the boundary layer edge properties when you solve the BL equations. I didn’t get too far. You’d definitely need some surface curvature corrections in there too.

      My personal feeling, is that once you go 3D your time might be better spent just using CFD – haha.

      Reply

      • Hi and thanks for your answer. Yes I agree for full 3D. The case with d/dz = 0 is a simplified case however. This adds only one term in the x-momentum equation and one term in the energy equation + one equation (the z-momentum equation). Since you solve segregated, I guess it would be feasible just to add the solve z-momentum equation to your algorithm above. Would you be willing to share your 2-D code? I could test it for simple flows (flat plate, flate plate with pressure gradient,…) and then try to extend it to 3-D.
        Thanks,
        Anthony

        Reply

Leave a Reply

This site uses Akismet to reduce spam. Learn how your comment data is processed.

%d