Blasius Flat-Plate Flow Benchmark

From OpenFOAMWiki
Revision as of 05:40, 25 February 2008 by Mike Jaworski (Talk | contribs)

1 Benchmarking OpenFOAM Solvers Against the Blasius Similarity Solution

Computational simulation is only as good as it accurately reflects the real world which makes benchmarking an important step in creating solvers. One option for benchmarking studies is to use analytical solutions, where available. One such solution exists that describes both velocity and temperature in a laminar boundary layer.

Two solvers will be benchmarked against the analytical solution available. The first is a modified version of icoFoam called icoHeatFoam (the steps to creating this in this tutorial: How to Add Temperature to IcoFoam). Later, the conjugateHeatFoam solver will be benchmarked as well.

First, a brief description of the Blasius flat-plate flow solution is given. The steps to creating the necessary case files and running the solution will be covered and then steps for post-processing.

2 Blasius Similarity Solution

The similarity solution describes the formation of a boundary layer. While solutions exist for stagnation flow, this case will look at flat-plate flow only. The full derivation of the similarity solution can be found in numerous fluid dynamics texts, such as "Viscous Fluid Flow" by Frank White, 2003. One can also see an abbreviated discussion Blasius Boundary Layer. For the time being, only the main points of the solution will be described.

The first item is the similarity variable itself: \eta = y \sqrt{\frac{U}{2 \nu x}}

where  y is the direction normal to the plate, x is the direction along its length with zero being the leading edge. /nu is the kinematic viscosity of the fluid and U is the free stream velocity. The main idea of the similarity solution is that the solution is a function of the similarity variable only, such that the solution will appear the same (i.e. similar) so long as the value of /eta is the same.

Applying this similarity variable to the boundary-layer equations, one can produce an ordinary differential equation describing the velocity profile:  f''' + f f'' = 0 This ODE needs to be numerically integrated.

A dimensionless temperature can be defined as follows:  \Theta = \frac{T(\eta) - T_e}{T_w - T_e} Where T_e and T_w are the free-stream and wall temperatures respectively. Further integration of the similarity solution yields the relation for \Theta.

3 Mesh Construction and Physical Properties

BlockMesh will be utilized to construct the case. It is suggested that one copy the case files from the icoFoam cavity tutorial to a new case which will be modified for this test.

First, the desired mesh will have a small section of entrance flow, the flat plate will be located at x=0 to ease later comparison with the analytical theory. The test will consist of a flat-plate held at constant temperature (275 K) and liquid water (Pr = 7) with an entrance temperature of 350 K. We will test a free-stream velocity of 2 cm/s which should produce a 1 cm thick boundary layer 3.5 cm downstream of the leading edge. The upper boundary of the computational domain should be "far" from the flat plate. We will take "far" to be \eta = 50 which is 10cm away.

Begin by editing the blockMeshDict file located in the <case>/constant/polymesh subdirectory. In this tutorial, the computational domain is broken into 4 regions:

  • 2 x-direction sections
    • An entrance length
    • The flat-plate region
  • 2 y-direction sections
    • A finely graded region next to the plate
    • A more coarsely graded region above

Edit blockMeshDict to reflect the following information:

// * * * * * * * * * //

convertToMeters 0.001;

vertices
(
    (0 0 0)    //0
    (35 0 0)   //1
    (35 10 0)  //2
    (0 10 0)   //3
    (0 0 1)    //4
    (35 0 1)   //5
    (35 10 1)  //6
    (0 10 1)   //7
    (0 100 0)  //8
    (35 100 0) //9
    (0 100 1)  //10
    (35 100 1) //11
    (-20 0 0)  //12
    (-20 0 1)  //13
    (-20 10 0) //14
    (-20 10 1) //15
    (-20 100 0) //16
    (-20 100 1)	//17
);

blocks          
(
    hex (0 1 2 3 4 5 6 7) (50 30 1) simpleGrading (3 5 1)
    hex (3 2 9 8 7 6 11 10) (50 45 1) simpleGrading (3 3 1)
    hex (12 0 3 14 13 4 7 15) (25 30 1) simpleGrading (.2 5 1)
    hex (14 3 8 16 15 7 10 17) (25 45 1) simpleGrading (.2 3 1)
);

edges           
(
);

patches         
(
    wall top
    (
        (8 9 11 10)
	(16 8 10 17)
    )
    wall inlet
    (
        (14 16 17 15)
	(12 14 15 13)
    )
    wall outlet
    (
        (9 2 6 11)
	(2 1 5 6)
    )
    wall plate
    (
       
	(0 4 5 1)

    )
    wall symmBound
    (
	(12 13 4 0)
    )
    empty frontAndBack 
    (
        (3 0 1 2)
        (7 6 5 4)
	(8 3 2 9)
	(10 11 6 7)
	(14 12 0 3)
	(15 7 4 13)
	(16 14 3 8)
	(17 10 7 15)
    )
);

mergePatchPairs 
(
);


// ************************************************************************* //

Aside: an important note regarding the use of the blockMesh utility is the order of defining blocks and patches. Blocks are defined such that the order of vertices create a local coordinate system. The number of cells and cell grading is determined according to this local coordinate system, which may not be the same for all blocks if the same ordering pattern is not used for all. In the case of patches, each must be defined such that the face normal (using the right-hand rule) points away, or out of, the computational domain. Failure to adhere to the guidelines for vertex ordering detailed in the user manual can result in negative volumes which will not appear until a solver attempts to use the mesh.

The transport properties of liquid water are: \nu = 1.0e-06 \alpha = DT = 1.438e-07

These items should be entered into the transportProperties file.

4 Boundary and Initial Conditions

The patches described above can be set to either wall or patch type in the boundary file because the laminar solvers do not distinguish between these patch types.

Appropriate boundary conditions for velocity are to set the entrance bottom plane and outlet as zeroGradient patches. The inlet and top patches are "fixedValue" patches with velocity of (0.02 0 0). The plate patch is "fixedValue" with a velocity of (0 0 0). Time can be saved by setting the internalField to a uniform (0.02 0 0).

For the pressure boundary conditions, the outlet is a fixedValue type with value 0. All other boundary patches are zeroGradient.

Temperature boundary conditions mimic the velocity field. The inlet and top patches are "fixedValue" with a value of 350. The symmBound and outlet patches receive a zeroGradient boundary. Finally, the plate itself is a fixedValue boundary of 275.

The following shows the information in the /0/T file:


dimensions      [0 0 0 1 0 0 0];

internalField   uniform 350;

boundaryField
{
    top
    {
        type            fixedValue;
	value           uniform 350;
    }
    inlet
    {
        type            fixedValue;
	value           uniform 350;
    }
    outlet
    {
        type            zeroGradient;
    }
    plate
    {
        type            fixedValue;
        value           uniform 275;
    }
    symmBound
    {
	type            zeroGradient;
    }
    frontAndBack
    {
        type            empty;
    }
}

5 Running the Case and post-processing

The last step before running the case is providing the correct settings to the control dictionary. Edit the controlDict file such that the case will run to 8 seconds, which should be enough time to reach steady state. Estimating time to run is complicated by having both velocity and temperature solutions evolving at the same time. The temperature equilibrates much slower than the velocity does (1/Pr as fast) and depends on the velocity solution.

A time step of 0.015 provides a maximum Courant number of approximately 0.9. An example controlDict file is shown here:

application     icoHeatFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         8;

deltaT          0.015;

writeControl    timeStep;

writeInterval   15;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression uncompressed;

timeFormat      general;

timePrecision   6;

runTimeModifiable yes;

One can run the case and examine the solution.

Magnitude of the velocity for the last time step.
Temperature solution for the last time step.

Case files are included in the following tarball: Media:icoHeatFoam-flatPlate2.tar.gz

The sample utility will be used to pull more quantitative data from the solution.

It is suggested that you copy an existing sampleDict file from another tutorial, such as the stressed plate tutorial in the user manual.

Several lines will be created to read the velocity and temperature of the fluid. These will be defined such that the starting point is 1cm above the plate and ends at y=0.0. This method is chosen to avoid a small bug in the sample utility which sometimes causes it to hang if the starting point is defined too close to the boundary of the domain.

The two items read will be the x component of the velocity and the temperature. Rather than us the Ucomponents utility as in the user manual tutorial, the velocity component we desire will be obtained by specifying U.component(0). The one drawback of this method is that the utility will create filenames with the special characters "(" and ")" in them.

The sampleDict file should resemble the following:

interpolationScheme cellPoint;

writeFormat     raw;

sampleSets
(
    uniform
    {
	name		x005;
	axis		y;
	start		(0.005 0.01 0.0005);
	end		(0.005 0.0   0.0005);
	nPoints		100;
    }
    uniform
    {
	name		x015;
	axis		y;
	start		(0.015 0.01 0.0005);
	end		(0.015 0.0   0.0005);
	nPoints		100;
    }
    uniform
    {
	name		x025;
	axis		y;
	start		(0.025 0.01 0.0005);
	end		(0.025 0.0   0.0005);
	nPoints		100;
    }
    uniform
    {
	name		x030;
	axis		y;
	start		(0.03 0.01 0.0005);
	end		(0.03 0.0   0.0005);
	nPoints		100;
    }
);

fields
(
    U.component(0)
    T
);

Run the sample utility which, as written, will proceed through every available time step available. Each subset in the sample utility produces its own file showing the <name> value at the start.

A closer examination of the data can be made from the last time step can be made. By renaming the data files to avoid the special characters, a gnuplot plotting script can be written which will compare the obtained data to the Blasius solution itself. This could look like the following:

reset
set term x11 enhanced font 'Helvetica, 18'

U0 = 0.2

eta(r, s) = r * sqrt(U0 / (2 * 1e-6 * s))

set key inside bottom right vertical Left
set pointsize 1.5

set xr [0:4]
set yr [0:0.025]
set xlabel '{/Symbol h} [-]'
set ylabel 'U_x [m/s]'

plot 'BlasiusF.dat' u ($1):( $3) t 'Blasius solution' w l,\
     'x030_Ux_T.xy' every 2 u ( eta($1, 0.03) ):($2/U0)  \
     		    t 'x = 3.0cm' w points pt 5,\
     'x025_Ux_T.xy' every 2 u ( eta($1, 0.025) ):($2/U0)  \
     		    t 'x = 2.5cm' w points pt 6,\
     'x015_Ux_T.xy' every 2 u ( eta($1, 0.015) ):($2/U0) \
     		    t 'x = 1.5cm' w points pt 9,\
     'x005_Ux_T.xy' every 2 u ( eta($1, 0.005) ):($2/U0) \
     		    t 'x = 0.5cm' w points pt 11

This and a similar script for comparing the \Theta variable produce the following graphs:

Comparison of velocity of OpenFOAM solution with Blasius similarity profile
Comparison of dimensionless temperature of OpenFOAM solution with Blasius similarity profile

All necessary files to produce these graphs are included in this tarball: Media:gnuplot_blasius_compare.tar.gz

6 Discussion of the results

The first main point is that this is not an optimal setup for comparison because of the closed boundary at the top of the solution domain. The growth of the boundary layer necessarily accelerates the fluid elsewhere to maintain fluid continuity. As a result, the "apparent" free-stream velocity is slightly higher than the 0.02m/s which was specified in the case files. By adjusting the free-stream velocity used in the Blasius solution, however, the OpenFOAM solution is in very good agreement with the analytical solution.

7 conjugateHeatFoam Benchmarking

This will be filled in when the solver is fixed.

--Mike Jaworski