* Lab 3 - 2D BL,
Laminar and Turbulent*

**Objectives**

This lab exercise involves solutions to the 2-D, steady, incompressible, boundary layer ordered, Navier-Stokes equations via the variable-implicit Taylor Series modified Galerkin Weak Statement.

First, the case of laminar flow is examined via the following objectives:

- Verification of the input aPSE template code.
- Corroboration of solution convergence results with theory via a uniform mesh refinement study for Reynolds numbers of 10
^{4}, 10^{5}, and 10^{6}. - Comparison of thermal boundary layer solutions for both upwards and downwards flow with the isothermal case for Re = 10
^{4}. - Confirmation of Galerkin Weak Statement (GWS
^{h}) solution optimality in comparison with a problem execution in the Finite Volume Statement (GWS^{h}) form.

In addition, the case of turbulent flow is examined via the following objectives:

- Verification of the input aPSE templates for the Prandtl Mixing Length (MLT) and Turbulent Kinetic Energy (TKE) turblence closure models.
- Comparison of turbulent solutions with experimental data from the Stanford 1968 Validation Symposium.

**Problem Statement **

The first part of the lab considers the case of external laminar fluid flow over a flat plate. Buoyancy effects are also considered. If the additional assumptions of steady, incompressible flow with constant properties is applied, the well known Navier-Stokes equations reduce to a more tractable form. The non-dimensionalized continuity, momentum, and energy equations are shown below.

It should be noted that the irreversible work term shown in the energy equation is left out of the numerical analysis. Since the Eckert number is proportional to the Mach number squared, the assumption of incompressible flow (M < 0.3) renders this term negligible.

The resulting equations are essentially as the unsteady, convection-diffusion problems from last semester. Here, however, the time term is replaced with the "x" distance along the wall. An implicit Taylor series expansion of the Galerkin Weak Statement then allows for a numerically tractable expression:

The residuals used in the definition of the problem statement are shown below in their aPSE pseudo-template form. The blue lines indicate the code that would be changed to implement a Finite Volume Statement modification.

RESIDUALSINITIAL VALUE TERM IN F(U) ------------------------------------- ()()(UBAR)(0;1)(A3000)(-U) # ()()(UBAR)(0;1)(A3000F)(-U) RESIDUAL TERMS IN F(U) ------------------------------------- ()()(V)(0;0)(A3001)(U) +(REI)()()(0;-1)(A211)(U) +(DPDX)()()(0;1)(A200)() +(GDOTI,GR,RE2I)()()(0;1)(A200)(TEMP) CONTINUITY RESIDUAL ------------------------------------- (VMUL)()()(0;0)(A201V)(V) +(VMUL,HALF)()()(0;1)(A201H)(UPR) INITIAL VALUE TERM IN F(TEMP) ------------------------------------- ()()(UBAR)(0;1)(A3000)(-TEMP) # ()()(UBAR)(0;1)(A3000F)(-TEMP) RESIDUAL TERMS IN F(TEMP) ------------------------------------- ()()(V)(0;0)(A3001)(TEMP) +(PEI)()()(0;-1)(A211)(TEMP)

With the problem statement in a Taylor series-modified Galerkin Weak Statement form, Newton's iterative solution proceedure is implimented to converge upon a boundary layer flow distribution at each position across the wall:

where

where

where

As shown above, the residuals are functions of the two velcity dimensions as well as the temperature. The Jacobian then must be computed for each of the 9 derivative combinations. The resulting aPSE pseudo-code is given below. Again, the commented blue lines indicate the sole code modifications necessary to convert the problem to a finite volume implimentation.

JACOBIANSINITIAL VALUE TERM IN DF(U)/D(U) ------------------------------------- ()()(U)(0;1)(A3000)() # ()()(U)(0;1)(A3000F)() RESIDUAL TERMS IN DF(U)/D(U) ------------------------------------- ()()(V)(0;0)(A3001)() +(REI)()()(0;-1)(A211)() RESIDUAL TERM IN DF(U)/D(V) ------------------------------------- ()()(U)(0;0)(A3100)() RESIDUAL TERM IN DF(U)/D(TEMP) ------------------------------------- (GDOTI,GR,RE2I)()()(0;1)(A200)() RESIDUAL TERM IN DF(V)/D(U) ------------------------------------- (VMUL,HALF,ATERM)()()(0;1)(A201H)() RESIDUAL TERM IN DF(V)/D(V) ------------------------------------- (VMUL)()()(0;0)(A201V)() RESIDUAL TERM IN DF(TEMP)/D(U) ------------------------------------- (HALF)()(TEMP)(0;1)(A3000)() (-,HALF)()(TEML)(0;1)(A3000)() # (HALF)()(TEMP)(0;1)(A3000F)() # (-,HALF)()(TEML)(0;1)(A3000F)() RESIDUAL TERM IN DF(TEMP)/D(V) ------------------------------------- ()()(TEMP)(0;0)(A3100)() RESIDUAL TERM IN DF(TEMP)/D(TEMP) ------------------------------------- ()()(UBAR)(0;1)(A3000)() RESIDUAL TERMS IN DF(TEMP)/D(TEMP) ------------------------------------- +()()(V)(0;0)(A3001)() +(PEI)()()(0;-1)(A211)()

In the second part of the lab, we "shift tracks" and briefly consider the more commonly experienced case of turbulent flow. The previous assumptions of 2-D, steady, incompressible flow with constant properties still applies. As detailed in the course lectures and problem experiences, the resulting time averaged Navier-Stokes equations become somewhat more difficult:

The important things to notice here is the approximation that the turbulent Prandtl number equals the Prandtl number. Also, it should be noticed that the pressure gradient is no longer just a function of x.

The required pseudo-code modifications for the residual and jacobian matrices are shown below. The red color highlights the additions from the laminar code. For turbulent flow either a Prandtl Mixing Length Theory (MTL) or a Turbulent Kinetic Energy (TKE) models are used to provide closure to system with its added unknowns.

RESIDUALSINITIAL VALUE TERM IN F(U) ------------------------------------- ()()(UBAR)(0;1)(A3000)(-U) RESIDUAL TERMS IN F(U) ------------------------------------- ()()(V)(0;0)(A3001)(U) +(REI)()()(0;-1)(A211)(U) +(REI)()(RET)(0;-1)(A3011)(U) +(DPDX)()()(0;1)(A200)() +(GDOTI,GR,RE2I)()()(0;1)(A200)(TEMP) CONTINUITY RESIDUAL ------------------------------------- (VMUL)()()(0;0)(A201V)(V) +(VMUL,HALF)()()(0;1)(A201H)(UPR) INITIAL VALUE TERM IN F(TEMP) ------------------------------------- ()()(UBAR)(0;1)(A3000)(-TEMP) RESIDUAL TERMS IN F(TEMP) ------------------------------------- ()()(V)(0;0)(A3001)(TEMP) +(PEI)()()(0;-1)(A211)(TEMP) +(PEI)()(RET)(0;-1)(A3011)(TEMP)

JACOBIANSINITIAL VALUE TERM IN DF(U)/D(U) ------------------------------------- ()()(U)(0;1)(A3000)() RESIDUAL TERMS IN DF(U)/D(U) ------------------------------------- ()()(V)(0;0)(A3001)() +(REI)()()(0;-1)(A211)() +(REI)()(RET)(0;-1)(A3011)() RESIDUAL TERM IN DF(U)/D(V) ------------------------------------- ()()(U)(0;0)(A3100)() RESIDUAL TERM IN DF(U)/D(TEMP) ------------------------------------- (GDOTI,GR,RE2I)()()(0;1)(A200)() RESIDUAL TERM IN DF(V)/D(U) ------------------------------------- (VMUL,HALF,ATERM)()()(0;1)(A201H)() RESIDUAL TERM IN DF(V)/D(V) ------------------------------------- (VMUL)()()(0;0)(A201V)() RESIDUAL TERM IN DF(TEMP)/D(U) ------------------------------------- (HALF)()(TEMP)(0;1)(A3000)() (-,HALF)()(TEML)(0;1)(A3000)() RESIDUAL TERM IN DF(TEMP)/D(V) ------------------------------------- ()()(TEMP)(0;0)(A3100)() RESIDUAL TERM IN DF(TEMP)/D(TEMP) ------------------------------------- ()()(UBAR)(0;1)(A3000)() RESIDUAL TERMS IN DF(TEMP)/D(TEMP) ------------------------------------- +()()(V)(0;0)(A3001)() +(PEI)()()(0;-1)(A211)() +(PEI)()(RET)(0;-1)(A3011)()

Figure 1: The given geometric flow profile for laminar flow over a flat plate in external flow.

Figure 1 above gives a pictoral representation of the problem statement for both the laminar and turbulent flow problems. The numerical model chosen requires an "initial condition" which, in this case, is an initial fluid velocity profile. As shown in the figure, the initial velocity profile is specified with a value of zero at the wall surface extending linearly over 1/8 of the boundary layer thickness to a uniform flow velocity. It should be noted that, as the mesh increases, the initial velocity profile input into aPSE must also change as shown below in Table 1.

Table 1: Required initial velocity profile variation for aPSE code.

============================================================================= Elements Initial Velocity Profile _____________________________________________________________________________ 16 [S 0. .50 15*1. S1.591377] $ U 32 [S 0. .25 .50 .75 29*1. S1.591377] $ U 64 [S 0. .125 .25 .375 .50 .625 .75 .875 57*1. S1.591377] $ U 128 [S 0. .0625 .125 .1875 .25 .3125 .375 .4375 .50 .5625 .625 .6875 .75 .8125 .875 .9375 113*1. S1.591377] $ U =============================================================================

We therefore need an estimation of the boundary layer thickness for the problem. To obtain an approximation that can serve as a guide, we turn to the approximation of laminar, steady, imcompressible, parallel flow over a flat plate with constant properties, negligible viscous dissipation, and with a negligible pressure gradient (dp/dx ~ 0). Under these approximations, we arrive at a solution for the boundary layer thickness:

For this analysis a plate length of 3 ft is chosen. Table 2 below shows the resulting boundary layer thicknesses and freestream flow velocities resulting from the specified Reynolds numbers. The data suggests a boundary layer thickness of 0.10 ft should suffice in the analysis. In actuality, a boundary layer thickness of 0.15 ft was chosen and used due to a miscalculation. This should not affect the results as long as proper meshes are used to ensure resolution across the actual boundary layers.

Table 2: Estimation of the boundary layer thickness for a 3 ft long flat plate.

=================================================== Reynolds U_inf BL Number Thickness Re_{L}(ft/s) (ft) (in) ___________________________________________________ 3x10^{4}1.591377 0.086603 1.039236 3x10^{5}15.91377 0.027386 0.328632 3x10^{6}159.1377 0.008660 0.103923 ===================================================

**Discussion of Results **

Table 3 lists the parameters used to drive the Newton iterative solution algorithm which was used to solve the laminar boundary layer problem. Table 4 lists the fluid transport properties used in the problem. Unless specified, all solutions below are isothermal at T = 70^{o}F).

Table 3: System constants pertaining to the Newton iterative solution algorithm and their numerical values.

============================================================= Variable Value Units _____________________________________________________________ Initial X 3 ft Final X 6 ft Initial X Step -0.0025 ft X Step Multiplier 1.5 Max. X Step 0.2 ft Convergence Criteria 1E-4 Iteration Criteria 0.001 ft/s Implicitness Factor 0.5 =============================================================

Table 4: System constants pertaining to fluid transport and their numerical values.

============================================================= Variable Value Units _____________________________________________________________ Viscosity 3.7044E-07 lbf-sec/ft^2 Kinematic Viscosity 1.5913E-04 ft^2/sec Density 2.3279E-03 slug/ft^3 =============================================================

Figures 2-4 detail the velocity profile change across the length of the flate plate for all three specified Reynolds numbers. Indeed, the choice of a 0.15 ft boundary layer was too high. The 128 element mesh, however, seems to capture the velocity profile and gives smooth variations as expected.

Figure 2: Laminar velocity profiles across the length of the flat plate at RE_{L} = 3x10^{4}.

Figure 3: Laminar velocity profiles across the length of the flat plate at RE_{L} = 3x10^{5}..

Figure 4: Laminar velocity profiles across the length of the flat plate at RE_{L} = 3x10^{6}..

Figure 5 shows a normalized comparison of the velocity profiles at the end of the flat plate for all three Reynolds numbers. It should be noted that the legend has accidentally been reversed in this figure. As expected, the higher Reynolds number results in a more compressed velocity profile.

Figure 5: Comparison of the laminar velocity profiles at the plate end for all three Reynolds numbers.

Figure 6 shows the element-wise x-component velocity norm. Since an "optimal" solution would spread the errors out equally across the elements, we can say that there is still room for improvement in the solution. Perhaps a non-uniform mesh would have done a better job. Such a mesh would have been smaller at the plate surface and might have minimized the error there, thus "evening-out" the profile.

Figure 6: A comparison of the x-velocity norms across the boundary layer for all three Reynolds numbers.

In an attempt to prove the underlying physics of the problem, Figure 7 gives a comparison of the x-component velocity norm between the Galerkin Weak Statement and Finite Volume Statement methods. As seen, the GWS solution is consistently lower which corresponds to theory. It is the most optimal solution!

Figure 7: A comparison of the optimality between the GWS^{h} and FVS^{h} finite element methods.

It is also desired to validate the asymptotic convergence theory for a uniform mesh refinement. Tables 5-7 show the results and indicate a convergence slope of 2 for all three Reynolds numbers. Since the aPSE program incorporated a linear element basis, a slope of 2 is correct!

Table 5: Convergence results for Re_{L} = 3x10^{4}.

=============================================================== Mesh Element Enorm Enorm_error Enorm_Slope Size Length (ft) _______________________________________________________________ 16 0.009375 0.00084199 0 0 32 0.0046875 0.000843322 4.441e-007 0 64 0.00234375 0.000843681 1.19733e-007 1.89106 128 0.00117188 0.000843772 3.02e-008 1.9872 ===============================================================

Table 6: Convergence results for Re_{L} = 3x10^{5}.

=============================================================== Mesh Element Enorm Enorm_error Enorm_Slope Size Length (ft) _______________________________________________________________ 16 0.009375 0.000192849 0 0 32 0.0046875 0.000192562 -9.55333e-8 0 64 0.00234375 0.000192498 -2.13333e-8 2.16289 128 0.00117188 0.000192483 -5e-9 2.09311 ===============================================================

Table 7: Convergence results for Re_{L} = 3x10^{6}.

================================================================ Mesh Element Enorm Enorm_error Enorm_Slope Size Length (ft) ________________________________________________________________ 16 0.009375 2.5244e-005 0 0 32 0.0046875 2.47772e-005 -1.556e-07 0 64 0.00234375 2.46475e-005 -4.32333e-08 1.84763 128 0.00117188 2.462e-005 -9.16667e-09 2.23767 ================================================================

Finally it is desired to explore the effects of buoyancy upon the laminar solution. As shown in Table 8, these solutions incorporate a plate surface temperature of 100^{o}F while the initial fluid remains at 70^{o}F. Figure 8, shown below, compares the exit velocity profiles at the minimum Reynolds number for three cases: (1) isothermal flow, (2) flow against the gravity gradient, and (3) flow with the gravity gradient. As expected, the fluid convection results in an increase in the velocity profile regardless of the flow direction. When the fluid flows along with the gravity gradient, however, the maximum velocities are seen.

Table 8: System constants pertaining to heat transfer and their numerical values.

============================================================= Variable Value Units _____________________________________________________________ Wall Temperature 100 F Flow Temperature 70 F Gravity 32.17398 ft/sec^2 Thermal Conductivity 4.1759E-06 btu/ft-sec-R Specific Heat (Cp) 9.4681E-06 btu/lbm-R Ratio Specific heat 1.4 Prandtl No. 6.7644E-01 =============================================================

Figure 8: A comparison of thermal boundary layer effects at different plate orientations with respect to the gravity vector.

Figures 9-12 show the results of the Bradshaw I-2400 experiment. Comparisons are made between the Mixing Length Theory (MLT) and Turbulent Kinetic Energy (TKE) closure models. Little variation is observed between the two in the velocity profile. The TKE closure model results in a much higher turbulent Reynolds number, however. In comparison to the laminar solutions, one can clearly see the impact of turbulent mixing in that the high viscosity at the wall is propagated deeper into the boundary layer resulting in a very steep velocity profile.

Finally, it should be noted that some disturbing discontinuities are observed in the turbulent Reynolds number for the MLT closure model. These are seen in Figures 10 and 12 and the source is unknown. Is this a problem with aPSE or a phenomenon seen with MLT codes in general?

Figure 9: The turbulent velocity profiles across a flat plate for Bradshaw's I-2400 experiment using the MLT closure model.

Figure 10: Turbulent viscotiy profiles across a flat plate for Bradshaw's I-2400 experiment using the MLT closure model.

Figure 11: A comparison of the MLT and TKE closure models for Bradshaw's I-2400 experiment as seen in the velocity profile.

Figure 12: A comparison of the MLT and TKE closure models for Bradshaw's I-2400 experiment as seen in the turbulent Reynolds number.

**Conclusion **

In conclusion, nearly all of the objectives for the lab have been met resulting in an opportunity to "get one's hands dirty" in a CFD analysis. While the author feels he has gained useful insight into the nuances involved in forming and modeling the non-dimensinoalized BL Navier-Stokes equations, he has also become supremely frustrated with the difficulty of aPSE and Tecplot operation. In short, he misses MATLAB :-)