本案例利用OpenFOAM计算管道中湍流流动压降并对计算结果进行验证。

参考文献：F.M. White. Fluid Mechanics. 3rd Edition. McGraw-Hill Co., New York, NY. 1994.

”

1 案例描述

本案例计算空气流经光滑的水平管道，管道长度2m，半径0.002 m，空气密度1.225 kg/m3，粘度1.7894e-5 kg/m-s，管道入口速度50 m/s，出口压力 0 Pa。计算管道的压降。

采用轴对称模型。入口采用速度边界，50m/s；出口为压力出口，静压 0Pa；采用稳态求解。

湍流模型采用SST k-omega模型， 为保证Y+=1，采用Y+计算器可得第一层网格高度约为5e-6m，径向采用32层网格节点，轴向网格尺寸0.0005m。

2 案例验证

流经管道的雷诺数：

为湍流流动，采用布拉休斯公式计算摩擦系数（F. M. White. Fluid Mechanics(7nd ed). P366页）

则压力降：

3 OpenFOAM设置

本算例涉及到稳态湍流计算，可以使用simpleFoam求解器进行计算。

3.1 文件及网格准备

利用官方算例pitzDaily作为模板。

run
cp -r $FOAM_TUTORIALS/incompressible/simpleFoam/pitzDaily .
mv pitzDaily/ VM03

将网格文件VM03.msh拷贝到VM03文件夹中。

网格文件在ICEM CFD中进行处理。需要注意轴对称模型需要确保对称面为XY面，否则在网格转换时会报错。

• 利用命令转换计算网格

fluentMeshToFoam VM03.msh

• 转换完毕后需要修改constant/polyMesh/boundary文件，将SIDE1与SIDE2边界类型修改为wedge，文件内容如下所示

FoamFile
{
    version     2.0;
    format      ascii;
    class       polyBoundaryMesh;
    location    "constant/polyMesh";
    object      boundary;
}
// * * * * * * * * * * * * * * * //
 
5
(
    SIDE1
    {
        type            wedge;
        inGroups        List 1(wall);
        nFaces          2784;
        startFace       5324;
    }
    TOPWALL
    {
        type            wall;
        inGroups        List 1(wall);
        nFaces          232;
        startFace       8108;
    }
    OUTLET
    {
        type            patch;
        nFaces          12;
        startFace       8340;
    }
    INLET
    {
        type            patch;
        nFaces          12;
        startFace       8352;
    }
    SIDE2
    {
        type            wedge;
        inGroups        List 1(wall);
        nFaces          2784;
        startFace       8364;
    }
)

• 利用命令checkMesh检查网格

checkMesh

网格检查结果如下图所示。

网格质量不怎么好，长宽比有点儿高。注意对于轴对称模型，一定是在修改了wedge边界后再进行网格质量检查。

3.2 修改物性参数

本案例需要修改momentumTransport及transportProperties文件。

1、momentumProperties文件

此文件中需要指定流体介质的运动粘度。

修改文件内容为：

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      transportProperties;
}
// * * * * * * * * * * * * * * * * * * * //
 
transportModel  Newtonian;
 
nu              [0 2 -1 0 0 0 0] 1.46073e-05;

2、momentumTransport文件

本算例采用sst k-omega模型进行计算。文件内容为：

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      momentumTransport;
}
// * * * * * * * * * * * * * * * //
simulationType RAS;
RAS
{
    model           kOmegaSST;
    turbulence      on;
    printCoeffs     on;
}

3.3 设置边界条件

删除0文件夹中的多余文件：

cd 0
rm epsilon f nuTilda v2
cd ..

此时0文件夹中的文件如下图所示。

湍流参数采用公式进行估算。

1、U文件

修改U文件，内容如下所示。

FoamFile
{
    version     2.0;
    format      ascii;
    class       volVectorField;
    object      U;
}
// * * * * * * * * * * * * * * * * * * //
 
dimensions      [0 1 -1 0 0 0 0];
internalField   uniform (50 0 0);
 
boundaryField
{
    INLET
    {
        type            fixedValue;
        value           uniform (50 0 0);
    }
 
    OUTLET
    {
        type            zeroGradient;
    }
 
    TOPWALL
    {
        type            noSlip;
    }
 
    "(SIDE1|SIDE2)"
    {
        type            wedge;
    }
}

2、p文件

文件内容如下所示。

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      p;
}
// * * * * * * * * * * * * * * * //
dimensions      [0 2 -2 0 0 0 0];
internalField   uniform 0;
 
boundaryField
{
    INLET
    {
        type            zeroGradient;
    }
 
    OUTLET
    {
        type            fixedValue;
        value           uniform 0;
    }
 
    TOPWALL
    {
        type            zeroGradient;
    }
 
    "(SIDE1|SIDE2)"
    {
        type            wedge;
    }
}

3、k文件内容

k文件内容如下所示。

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "0";
    object      k;
}
// * * * * * * * * * * * * * * //
dimensions      [0 2 -2 0 0 0 0];
internalField   uniform 8.87476;
 
boundaryField
{
    INLET
    {
        type            fixedValue;
        value           uniform 8.87476;
    }
    OUTLET
    {
        type            zeroGradient;
    }
    TOPWALL
    {
        type            kqRWallFunction;
        value           uniform 8.87476;
    }
 
    "(SIDE1|SIDE2)"
    {
        type            wedge;
    }
}

4、nut文件

文件内容如下所示。

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "0";
    object      nut;
}
// * * * * * * * * * * * * * * * * * * //
dimensions      [0 2 -1 0 0 0 0];
internalField   uniform 0;
 
boundaryField
{
    INLET
    {
        type            calculated;
        value           uniform 0;
    }
    OUTLET
    {
        type            calculated;
        value           uniform 0;
    }
    TOPWALL
    {
        type            nutkWallFunction;
        value           uniform 0;
    }
 
    "(SIDE1|SIDE2)"
    {
        type            wedge;
    }
}

5、omega文件

文件内容如下。

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "0";
    object      omega;
}
// * * * * * * * * * * * * * * * * * * //
 
dimensions      [0 0 -1 0 0 0 0];
internalField   uniform 23795;
 
boundaryField
{
    INLET
    {
        type            fixedValue;
        value           $internalField;
    }
    OUTLET
    {
        type            zeroGradient;
    }
    TOPWALL
    {
        type            omegaWallFunction;
        value           $internalField;
    }
    "(SIDE1|SIDE2)"
    {
        type            wedge;
    }
}

3.4 设置计算控制文件

• 利用命令foamGet residuals添加残差字典文件

1、controlDict文件

文件内容如下所示。

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      controlDict;
}
// * * * * * * * * * * * * * * * * * * //
application     simpleFoam;
startFrom       startTime;
startTime       0;
stopAt          endTime;
endTime         2000;
deltaT          1;
writeControl    timeStep;
writeInterval   100;
// 只保留最后的5个文件
purgeWrite      5;
writeFormat     ascii;
writePrecision  6;
writeCompression off;
timeFormat      general;
timePrecision   6;
runTimeModifiable true;
functions
{  
    #includeFunc residuals
}

3.5 执行计算

本算例利用simpleFoam进行计算。

simpleFoam

也可以利用命令显示计算残差：

foamMonitor -l postProcessing/residuals/0/residuals.dat

残差显示如下图所示。

采用的计算网格质量太差了。可以通过改善计算网格来提高收敛性。

若有更好的计算机，可以尝试采用并行计算。

foamGet decomposeParDict

修改system/decomposeParDict文件

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      decomposeParDict;
}
// * * * * * * * * * * * * * * * //
numberOfSubdomains  40;
method              scotch;

利用命令进行计算：

decomposePar
mpirun -np 40 simpleFoam -parallel

计算完毕后可以利用命令组合结果：

reconstruct

或采用命令直接查看结果：

paraFoam -builtin

4 计算结果

在paraView中统计入口压力。

• 创建一个无限接近于入口的Slice，如下图所示

• 在Slice1上添加integrate Variables，采用默认参数

• 在IntegrateVariables1上添加Calculator

• 如下图所示设置参数

• 软件显示得到的平均压力为16985

注意到这里得到的是压力为运动压力，其单位为，其与密度的乘积为实际压力。由于出口压力被设置为0，因此管道内的压力降为：

利用经验公式计算得到的压力为22366 Pa，误差约为-6.97%。

误差相对有点大，改善计算网格以及利用充分发展入口条件都有利于提高计算精度。