本算例演示利用OpenFOAM计算层流管道中的非牛顿流体流动，并对计算结果进行验证。

参考资料:W.F. Hughes and J.A. Brighton. Schaum's Outline of Theory and Problems of Fluid Dynamics. McGraw-Hill Book Co., Inc., New York, NY. 1991.

”

1 模型描述

本案例描述的模型如图所示。

其中管道长度0.1m，直径 0.0025 m，入口为平均速度2m/s的完全发展层流速度边界，管道出口为静压为0的压力出口。流经管道中的介质密度 1000 kg/m3，动力粘度为幂率流体。

其中，k=10，n=0.4。

2 OpenFOAM设置

本算例使用nonNewtonianIcoFoam求解器进行计算。

2.1 文件准备

打开终端进入工作路径下，运行下面的命令准备模板算例。

cp -r $FOAM_TUTORIALS/incompressible/nonNewtonianIcoFoam/offsetCylinder/ .
mv offsetCylinder VM07
cd VM07

删除多余的文件，最终算例文件结构如下所示。

2.2 网格准备

本算例采用导入ICEM CFD准备的msh文件的方式准备网格。

• 将文件VM07.cas拷贝至算例文件夹中，运行以下命令转换网格

fluentMeshToFoam VM07.msh

• 修改网格文件cosntant/polyMesh/boundary，将边界SIDE1与SIDE2的类型改为wedge

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

• 利用命令wedge检查网格，确保计算网格没有错误或警告信息，如下图所示

2.3 材料属性设置

算例中的流体介质为幂律流体。修改constant/transportProperties文件，详细内容如下所示。注意这里的粘度值要除以材料的密度1000 kg/m3。

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      transportProperties;
}
// * * * * * * * * * * * * * //
transportModel  powerLaw;
 
powerLawCoeffs
{
    nuMax [ 0 2 -1 0 0 0 0 ] 1;
    nuMin [ 0 2 -1 0 0 0 0 ] 1e-5;
    k [ 0 2 -1 0 0 0 0 ] 0.01;
    n [ 0 0 0 0 0 0 0 ] 0.4;
}

2.4 边界条件与初始条件

1、U文件

U文件内容如下所示。

FoamFile
{
    version     2.0;
    format      ascii;
    class       volVectorField;
    object      U;
}
// * * * * * * * * * * * //
dimensions      [0 1 -1 0 0 0 0];
internalField   uniform (2 0 0);
boundaryField
{
    INLET
    {
        type            codedFixedValue;     
        value           uniform (2 0 0);
        name            parabolicVelocity;
        code
        #{
            // 下面三行为标准写法，一般不用修改
            const fvPatch& boundaryPatch = patch();
            const vectorField& Cf = boundaryPatch.Cf();
            vectorField& field = *this;
            // U_0为2倍的平均速度；p_ctr为中心点偏移量；p_r为半径
            scalar U_0 = 4, p_ctr = 0, p_r = 0.00125;
 
            forAll(Cf, faceI)
            {
                field[faceI] = vector(U_0*(1-(pow(Cf[faceI].y()-p_ctr,2))/(p_r*p_r)),0,0);
            }
        #};
    }
    SIDE1
    {
        type            wedge;
    }
    OUTLET
    {
        type            zeroGradient;
    }
    WALL
    {
        type            noSlip;
 
    }
    SIDE2
    {
        type            wedge;
    }
}

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

2.5 计算控制参数设置

设置controlDict文件，内容如下所示。

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      controlDict;
}
// * * * * * * * * * * * * * * * //
application     nonNewtonianIcoFoam;
startFrom       startTime;
startTime       0;
stopAt          endTime;
endTime         0.1;
deltaT          0.00002;
writeControl    runTime;
writeInterval   0.001;
purgeWrite      5;
writeFormat     ascii;
writePrecision  6;
writeCompression off;
timeFormat      general;
timePrecision   6;
runTimeModifiable true;

非牛顿流体计算稳定性很差，建议先使用牛顿流体计算0.5s，再改换为非牛顿流体计算。

即先将transportProperties文件修改为下面的形式：

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

修改controlDict文件为：

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      controlDict;
}
// * * * * * * * * * * * * * * * * * //
application     nonNewtonianIcoFoam;
startFrom       startTime;
startTime       0.0;
stopAt          endTime;
endTime         0.05;
deltaT          0.00002;
writeControl    runTime;
writeInterval   0.001;
purgeWrite      5;
writeFormat     ascii;
writePrecision  6;
writeCompression off;
timeFormat      general;
timePrecision   6;
runTimeModifiable true;

执行命令nonNewtonianIcoFoam进行计算，计算完毕后修改controlDict文件。

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      controlDict;
}
// * * * * * * * * * * * * * * * * //
application     nonNewtonianIcoFoam;
startFrom       startTime;
startTime       0.05;
stopAt          endTime;
endTime         0.1;
deltaT          0.00002;
writeControl    runTime;
writeInterval   0.001;
purgeWrite      5;
writeFormat     ascii;
writePrecision  6;
writeCompression off;
timeFormat      general;
timePrecision   6;
runTimeModifiable true;
 

同时修改transportProperties文件。

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      transportProperties;
}
// * * * * * * * * * * * * * * * //
transportModel  powerLaw;
powerLawCoeffs
{
    nuMax [ 0 2 -1 0 0 0 0 ] 1;
    nuMin [ 0 2 -1 0 0 0 0 ] 1e-5;
    k [ 0 2 -1 0 0 0 0 ] 0.01;
    n [ 0 0 0 0 0 0 0 ] 0.4;
}

继续执行命令nonNewtonianIcoFoam进行计算计算直至结束。

3 计算结果

• 0.1s时刻压力分布

• 0.1s时刻速度分布

从速度分布来看，0.1s的计算时间还不够长。

• 0.1s时刻运动粘度

• 统计入口压力

统计得到的运动压力值为59.7596，流体密度值为1000 kg/m3，因此进出口压降为59.7596 kPa。

与文献中比较如下表所示。

误差稍微有点儿大。从计算结果云图上来看，需要继续往下计算，可以尝试计算到0.5s再统计结果。

计算文件：