本案例利用OpenFOAM计算水银通过具有均匀热通量壁面的圆形管道内的层流流动现象。

1 案例描述

案例模型如下图所示，水银在管道入口处为充分发展湍流速度分布，平均速度为0.005 m/s。利用OpenFOAM计算管道压降及出口中心点位置温度，并与解析解进行比较。案例采用二维轴对称模型计算。管道半径为0.0025 m，管道长度0.1 m。本案例利用ICEM CFD生成全四边形计算网格。

2 介质参数

流体介质为水银，其材料介质参数为：

• 密度： kg/m3
• 粘度：  kg/m-s
• 比热： J/kg-K
• 导热率： W/m-K

3 边界条件

4 解析计算方式

1、压降计算

理论值来自 ：

Theodore L. Bergman, Adrienne S. Lavine, Frank P. Incropera, David P. DeWitt. Fundamentals of Heat and Mass Transfer [ 7ed. ] . P523。

”

5 OpenFOAM计算

注意的问题：

• 计算模型。采用导入Fluent cas的形式，注意本算例计算的是轴对称问题，应用到wedge边界。
• 充分发展流动。这里使用codeFixedValue来完成。
• 求解器。本算例可以选择使用buoyantSimpleFoam求解器进行计算，也可以改造icoFoam求解器使其具备温度求解功能。这里改造求解器。
• 边界条件。本算例多出一个热流边界，此边界在OpenFOAM中使用fixedGradient来解决。

5.1 改造求解器icoTempFoam

icoFoam求解器本身不具备温度求解功能，这里通过改造求解器源代码使其能够求解温度场。

注：此部分内容参阅http://openfoamwiki.net/index.php/How_to_add_temperature_to_icoFoam。

”

• 利用下面的命令拷贝icoFoam源代码

cd $FOAM_SOLVERS/incompressible
mkdir -p $WM_PROJECT_USER_DIR/applications/solvers
cp -r icoFoam $WM_PROJECT_USER_DIR/applications/solvers/icoTempFoam
cd $WM_PROJECT_USER_DIR/applications/solvers/icoTempFoam
mv icoFoam.C icoTempFoam.C

处理完毕后，icoTempFoam文件夹中的文件组织结构如下图所示。

• 修改Make/files文件内容。该文件指定了要编译的源文件名称及编译后的目标文件路径。

icoTempFoam.C
EXE = $(FOAM_USER_APPBIN)/icoTempFoam

此时可以尝试利用wmake命令进行编译，若能够正常编译，则表示文件准备没有问题。后面为icoTempFoam求解器添加温度求解功能。

• 修改creatFields.H文件，在其中添加热扩散系数。

dimensionedScalar DT
{
    "DT",
    dimViscosity,
    transportProperties.lookup("DT")
};

添加完毕后如下图所示。

• 在createFields.H文件中添加温度物理量T

Info<< "Reading field T\n" << endl;
volScalarField T
(
    IOobject
    (
        "T",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);

添加完毕后如下图所示。

• 在源文件icoTempFoam.C中添加温度场控制方程

fvScalarMatrix TEqn
{
    fvm::ddt(T) fvm::div(phi,T) - fvm::laplacian(DT,T)
};
TEqn.solve();

添加完毕后如下图所示。

• 文件准备完毕后，运行命令wmake编译求解器

wmake

如下图所示表示编译成功。

5.2 准备网格

案例以OpenFOAM官方算例cavity作为模板。本案例采用导入Fluent Cas文件方式处理网格。

• 准备文件

run
cp -r $FOAM_TUTORIALS/incompressible/icoFoam/cavity/cavity .
mv cavity VM02
cd VM02

• 将网格文件WM02.cas放入VM02文件夹中

算例文件结构如下图所示。

注：为了将二维网格转换为轴对称网格，需要使用非官方工具makeAxialMesh，见https://github.com/cfdcoach/MakeAxialMesh.git。但这工具我一直没有编译成功。

”

这里采用ICEM CFD中的Extrude Mesh功能对2D网格进行处理，将其转化为三维结构，如下图所示(OUTLET边界在图中没有标识出来)。

ICEM CFD处理完毕后导出计算网格。

注：在ICEM CFD中旋转2D网格时，确保两个wedge面的对称面为轴平面。换句话说，旋转时可以先往一个方向旋转-1°，再在同一方向旋转2°，这样网格就关于XY平面两侧各偏转1°。

”

运行命令转换网格：

fluentMeshToFoam VM02.msh
checkMesh

网格转换完毕且检查没有问题后，打开文件constant/polyMesh/boundary修改边界，指定边界SIDE1与SIDE2的边界类型为wedge。修改后的boundary文件内容如下所示。

FoamFile
{
    version     2.0;
    format      ascii;
    class       polyBoundaryMesh;
    location    "constant/polyMesh";
    object      boundary;
}
// * * * * * * * * * * * * * * * //
5
(
    INLET
    {
        type            patch;
        nFaces          20;
        startFace       31180;
    }
    // 修改SIDE1边界类型为wedge
    SIDE1
    {
        type            wedge;
        inGroups        List<word> 1(wall);
        nFaces          16000;
        startFace       31200;
    }
    OUTLET
    {
        type            patch;
        nFaces          20;
        startFace       47200;
    }
    WALL
    {
        type            wall;
        inGroups        List<word> 1(wall);
        nFaces          800;
        startFace       47220;
    }
    // 修改SIDE2边界类型为wedge
    SIDE2
    {
        type            wedge;
        inGroups        List<word> 1(wall);
        nFaces          16000;
        startFace       48020;
    }
)

5.3 设置材料属性

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      transportProperties;
}
// * * * * * * * * * * * * * * * * * //
nu              [0 2 -1 0 0 0 0] 1.12573e-7;
DT              [0 2 -1 0 0 0 0] 4.51557e-6;

5.4 边界条件与初始条件

边界条件中需要加入温度场文件T，同时还需要修改p文件与U文件。

1、U文件

算例中入口采用充分发展流动，这里采用codeFixedValue进行指定。U文件内容如下所示。

FoamFile
{
    version     2.0;
    format      ascii;
    class       volVectorField;
    object      U;
}
// * * * * * * * * * * * * * * * * * * * //
dimensions      [0 1 -1 0 0 0 0];
internalField   uniform (0 0 0);
boundaryField
{
    INLET
    {
        // 指定类型为codeFixedValue
        type            codedFixedValue;
        // 指定边界初始值
        value           uniform (0 0 0);
        // 指定的名称标识符
        name            parabolicVelocity;
        // 编译时所需的信息，按实际需求给
        codeOptions
        #{
            -I$(LIB_SRC)/finiteVolume/lnInclude \
            -I$(LIB_SRC)/meshTools/lnInclude
        #};
        codeInclude
        #{
            #include "fvCFD.H"
            #include <cmath>
            #include <iostream>
        #};
        code
        #{
            // 下面三行为标准写法，一般不用修改
            const fvPatch& boundaryPatch = patch();
            const vectorField& Cf = boundaryPatch.Cf();
            vectorField& field = *this;
            // U_0为2倍的平均速度；p_ctr为中心点偏移量；p_r为半径
            scalar U_0 = 0.01, p_ctr = 0, p_r = 0.0025;
            forAll(Cf, faceI)
            {
                field[faceI] = vector(U_0*(1-(pow(Cf[faceI].y()-p_ctr,2))/(p_r*p_r)),0,0);
            }
        #};
    }
    "(SIDE1|SIDE2)"
    {
        type            wedge;
    }
    OUTLET
    {
        type            zeroGradient;
    }
    WALL
    {
        type            zeroGradient;
    }
}

2、p文件

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;
    }
    "(SIDE1|SIDE2)"
    {
        type            wedge;
    }
    OUTLET
    {
        type            fixedValue;
        value           uniform 0;
    }
    WALL
    {
        type            zeroGradient;
    }
}

3、T文件

将p文件拷贝一份，命名为T文件，修改其内容如下所示。

这里用fixedGradient指定壁面热流，需要利用傅里叶定律进行转换：

换算得到：

因此T文件修改为：

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      T;
}
// * * * * * * * * * * * * * * //
dimensions      [0 0 0 1 0 0 0];
internalField   uniform 300;
boundaryField
{
    INLET
    {
        type            fixedValue;
        value           uniform 300;
    }
    "(SIDE1|SIDE2)"
    {
        type            wedge;
    }
    OUTLET
    {
        type            zeroGradient;
    }
    WALL
    {
        type            fixedGradient;
        gradient        uniform 585.48;
    }
}

5.5 求解控制

1、controlDict设置

文件内容如下所示。

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      controlDict;
}
// * * * * * * * * * * * * * * * * * * //
application     icoTempFoam;
startFrom       startTime;
startTime       0;
stopAt          endTime;
// 计算的是瞬态，这个时间应当足够长
// 30s基本已经达到稳定
endTime         30;
// 时间步长0.01s能够确保库朗数小于1
deltaT          0.01;
writeControl    timeStep;
writeInterval   300;
purgeWrite      0;
writeFormat     ascii;
writePrecision  6;
writeCompression off;
timeFormat      general;
timePrecision   6;
runTimeModifiable true;

2、fvSchemes文件

文件内容修改为：

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * //
ddtSchemes
{
    default         Euler;
}
gradSchemes
{
    default         Gauss linear;
    grad(p)         Gauss linear;
}
divSchemes
{
    default         none;
    div(phi,U)      Gauss linear;
    div(phi,T)      Gauss upwind; //增加
}
laplacianSchemes
{
    default         Gauss linear orthogonal;
    laplacian       Gauss linear corrected; //增加
}
interpolationSchemes
{
    default         linear;
}
snGradSchemes
{
    default         orthogonal;
}

3、fvSolution文件

文件内容为：

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * //
solvers
{
    p
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-06;
        relTol          0.05;
    }
    pFinal
    {
        $p;
        relTol          0;
    }
    U
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-05;
        relTol          0;
    }
    // 增加T方程控制
    T
    {
        solver          PBiCGStab;
        preconditioner  DILU;
        tolerance       1e-7;
        relTol          0;
    }
}
PISO
{
    nCorrectors     2;
    nNonOrthogonalCorrectors 0;
    pRefCell        0;
    pRefValue       0;
}

6 计算结果

• 30s时刻速度分布

• 30 s时刻温度分布

• 30 s时刻出口位置速度分布

• 统计出口轴心温度，如下图所示为341.43K，文献值为341K

• 创建一个Slice

• 在Slice1上创建一个IntegrateVariable

• 在IntegrateVariables1上创建一个Calculator1，获取面上的平均压力

• 得到入口面上的平均压力，如下图所示，其压力值为7.20312e-5 m2/s2

出口静压在边界条件中被设置为0，因此无需通过计算器获取。

7 验证结果

通过前面的计算，可得到计算结果与理论值的比较，如下表所示。

这结果已经相当可以了。

声明：原创文章，欢迎留言与我讨论，如需转载留言