LES直流槽道边界层模拟，如何得到正则化速度u+以及正则化坐标y+？

学流体的小明

各位老师同学好，我正在用大涡模拟做一个直流槽道边界层模拟，雷诺数
$${\rm Re} _{\tau } = u _\tau h / \nu = 1000$$
其中$h$为槽道半高，$ u _\tau$为摩擦速度，$\nu$为运动粘度。根据Pope在Turbulent Flows一书中提到的公式
$${{\mathop {\rm Re} \nolimits} _\tau } \approx 0.09 {{\mathop {\rm Re} \nolimits} ^ {0.88}}$$
估计槽道的雷诺数为
$${\mathop{\rm Re}\nolimits} = {U_m}\left( {2h} \right)/\nu =3.9578e4\approx4e4$$
其中${U_m}$为平均流速。

所以我将OpenFOAM自带的槽道流算例channel395拿来，对于它的半高$h=1m$，运动粘度$\nu=2e-5m^2/s$的槽道，${U_m=0.4m/s}$。我修改了一下平均速度和网格。

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      fvOptions;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

momentumSource
{
    type            meanVelocityForce;

    selectionMode   all;

    fields          (U);
    Ubar            (0.4 0 0);
}


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

网格划分，粘性长度$\delta_\nu=h/{Re_{\tau}}=1e-3$。LES的网格尺寸建议为
$$50 \le \Delta {x^ + } \le 150,\Delta y_{wall}^ + < 1,15 \le \Delta {z^ + } \le 40$$
我使用的网格$\Delta {x^ + } =100$，$\Delta {z^ + } =50$，$\Delta {z^ + }$超了点，不过影响应该不会太大。
关键在于$\Delta {y}$，在于槽道半高范围内60个网格，最后一个网格的长度是第一个网格长度的25倍，算下来第一个网格的尺寸0.001103718，就是1个粘性长度多一点。
网格是符合要求的。

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

scale   1;

vertices
(
    (0 0 0)
    (4 0 0)
    (0 1 0)
    (4 1 0)
    (0 2 0)
    (4 2 0)
    (0 0 2)
    (4 0 2)
    (0 1 2)
    (4 1 2)
    (0 2 2)
    (4 2 2)
);

blocks
(
    hex (0 1 3 2 6 7 9 8)   (40 60 30) simpleGrading (1  25 1)
    hex (2 3 5 4 8 9 11 10) (40 60 30) simpleGrading (1 -25 1)
);

edges
(
);

boundary
(
    bottomWall
    {
        type            wall;
        faces           ((0 1 7 6));
    }
    topWall
    {
        type            wall;
        faces           ((4 10 11 5));
    }

    sides1_half0
    {
        type            cyclic;
        neighbourPatch  sides1_half1;
        faces           ((0 2 3 1));
    }
    sides1_half1
    {
        type            cyclic;
        neighbourPatch  sides1_half0;
        faces           ((6 7 9 8));
    }

    sides2_half0
    {
        type            cyclic;
        neighbourPatch  sides2_half1;
        faces           ((2 4 5 3));
    }
    sides2_half1
    {
        type            cyclic;
        neighbourPatch  sides2_half0;
        faces           ((8 9 11 10));
    }

    inout1_half0
    {
        type            cyclic;
        neighbourPatch  inout1_half1;
        faces           ((1 3 9 7));
    }
    inout1_half1
    {
        type            cyclic;
        neighbourPatch  inout1_half0;
        faces           ((0 6 8 2));
    }

    inout2_half0
    {
        type            cyclic;
        neighbourPatch  inout2_half1;
        faces           ((3 5 11 9));
    }
    inout2_half1
    {
        type            cyclic;
        neighbourPatch  inout2_half0;
        faces           ((2 8 10 4));
    }
);

mergePatchPairs
(
);

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

接下来是边界条件。
边界条件k，后面都是周期边界的cyclic

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "1";
    object      k;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    bottomWall
    {
        type            fixedValue;
        value           uniform 0;
    }
    topWall
    {
        type            fixedValue;
        value           uniform 0;
    }

边界条件nut

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "1";
    object      nut;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 2 -1 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    bottomWall
    {
        type            zeroGradient;
    }
    topWall
    {
        type            zeroGradient;
    }

边界条件nuTuilda

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "1";
    object      nuTilda;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 2 -1 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    bottomWall
    {
        type            fixedValue;
        value           uniform 0;
    }
    topWall
    {
        type            fixedValue;
        value           uniform 0;
    }

边界条件p

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "1";
    object      p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    bottomWall
    {
        type            zeroGradient;
    }
    topWall
    {
        type            zeroGradient;
    }

边界条件U

FoamFile
{
    version     2.0;
    format      ascii;
    class       volVectorField;
    location    "1";
    object      U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0.4 0 0);

boundaryField
{
    bottomWall
    {
        type            noSlip;
    }
    topWall
    {
        type            noSlip;
    }

算了很长时间，可以保证是完全发展了。channel395也就算到了1000s，后来我还往后延长算了一些，1000s是足够的。我这个channel1000算500000s肯定够了。
controlDict

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

application     pimpleFoam;

startFrom       latestTime;

startTime       0;

stopAt          endTime;

endTime         500000;

deltaT          0.1;

writeControl    timeStep;

writeInterval   100000;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;

functions
{
    #includeFunc yPlus
    #includeFunc wallShearStress

    fieldAverage1
    {
        type            fieldAverage;
        libs            (fieldFunctionObjects);
        writeControl    writeTime;

        fields
        (
            U
            {
                mean        on;
                prime2Mean  on;
                base        time;
            }

            p
            {
                mean        on;
                prime2Mean  on;
                base        time;
            }
        );
    }
}

算出来的结果使用postChannel工具进行空间平均，得到平均速度随y坐标变化的两组数，储存在 graphs/时间文件/Uf.xy 中。
例如130000s的时候

0.00110372  0.0133657
0.00337304  0.0407483
0.00576962  0.0693065
0.00830057  0.0985189
 0.0109734   0.127533
 0.0137962   0.155335
 0.0167772   0.181026
 0.0199254   0.204059
 0.0232501   0.224291
 0.0267613   0.241887
 0.0304693   0.257159
 0.0343852   0.270453
 0.0385207   0.282087
 0.0428881   0.292333
 0.0475004   0.301416
 0.0523713   0.309517
 0.0575154   0.316783
 0.0629479   0.323336
  0.068685   0.329275
 0.0747438   0.334681
 0.0811423   0.339623
 0.0878997    0.34416
 0.0950359   0.348342
  0.102572   0.352215
  0.110531   0.355819
  0.118936   0.359191
  0.127813   0.362361
  0.137187    0.36536
  0.147087   0.368214
  0.157542   0.370948
  0.168584   0.373584
  0.180244   0.376144
  0.192558   0.378647
  0.205563   0.381111
  0.219297   0.383551
  0.233801   0.385977
  0.249118     0.3884
  0.265294   0.390826
  0.282378   0.393261
  0.300419   0.395709
  0.319472   0.398174
  0.339593   0.400655
  0.360842   0.403151
  0.383283    0.40566
  0.406982    0.40818
   0.43201   0.410709
  0.458442   0.413247
  0.486356   0.415792
  0.515834    0.41834
  0.546966   0.420883
  0.579844   0.423407
  0.614565   0.425894
  0.651232   0.428324
  0.689956    0.43067
  0.730852   0.432905
   0.77404   0.434979
   0.81965   0.436819
  0.867818    0.43835
  0.918686   0.439484
  0.972407     0.4401

接下来就是正则化，画出速度分布了。别人的$ {\rm Re} _ {\tau } = 177.6 $结果是这样的。

正则化的关键在于计算摩擦速度$ u_{\tau} $，$ u^{+}=u/u_{\tau} $，$ y^{+}=yu_{\tau} / \nu $。
我这边使用了两种计算摩擦速度的方法。
第一种是认为此算例的摩擦雷诺数就是标称的1000，因此直接利用公式$u_{\tau }={\rm Re}_{\tau }\nu/h$计算摩擦速度。可以看到在50000s后，速度曲线在对数率部分与理论解吻合上了。但所有时间的速度曲线在粘性底层都无法和线性率的理论解对应。

第二种方法是，通过计算出来的速度-高度曲线，实际地进行计算，${u_\tau } = \sqrt {\nu {{\left. {\frac{{du}}{{dy}}} \right|}_{y = 0}}} $。这种方法画出的速度曲线在粘性底层（也就是线性率部分）和理论解是对应得上的。但是在对数率部分，数值结果明显过高，但可以看到50000s及之后的数值结果曲线和理论解的线基本平行的，就好像是理论解应该再高一些，数值结果就和理论解完美吻合了。

已经卡在这个问题上好久了。
应该不是网格的问题，现在的网格是符合要求的。
不是湍流没发展出来的问题，因为channel395计算的结果，我画出来也和上面这两张图差不多，所以感觉还是正则化的问题。
但感觉又不是正则化的问题，因为最终决定的因素就是摩擦速度，我也尝试过给一个具体的摩擦速度值，看看数值结果的曲线的形状能不能和理论解对得上，结果就像第一种方法画出来的图，对数率部分吻合了，粘性底层就差太多了。
看别人的图真真好看，比如那张绿底的。
求各位大佬指点迷津，我哪里做的不对？
也欢迎各位同学一起讨论，互相指点，少走弯路。

PS1：不知道上面的代码发出来会不会折叠，发帖预览的时候没有折叠感觉好长，有点担心实际效果。
PS2：插入在段落中的代码有时候不会被转写，至少是在预览界面。

李东岳

排版非常好。

你用的是层流模拟么，还是什么湍流模型？没看到你用什么湍流模型。

然后你这个槽道，u+y+是某个点上的还是底面平均之后的？

另外，你试一下用这个方式来计算u+与y+

{
    // Evaluate near-wall behaviour

    scalar nu = turbulence->nu()().boundaryField()[patchId][faceId];
    scalar nut = turbulence->nut()().boundaryField()[patchId][faceId];
    symmTensor R = turbulence->devSigma()().boundaryField()[patchId][faceId];
    scalar epsilon = turbulence->epsilon()()[cellId];
//    scalar omega = turbulence->omega()()[cellId];
    scalar k = turbulence->k()()[cellId];
    scalar magUp = mag(U[cellId] - U.boundaryField()[patchId][faceId]);

    scalar tauw = flowDirection & R & wallNormal;

    scalar uTau = ::sqrt(mag(tauw));

    scalar yPlus = uTau*y[cellId]/(nu + rootVSmall);

    scalar uPlus = magUp/(uTau + rootVSmall);

    scalar nutPlus = nut/nu;

    scalar kPlus = k/(sqr(uTau) + rootVSmall);

    scalar epsilonPlus = epsilon*nu/(pow4(uTau) + rootVSmall);

//    scalar omegaPlus = omega*nu/(sqr(uTau) + rootVSmall);

    scalar Rey = magUp*y[cellId]/nu;

    Info<< "Rey = " << Rey << ", uTau = " << uTau << ", nut+ = " << nutPlus
        << ", y+ = " << yPlus << ", u+ = " << uPlus
        << ", k+ = " << kPlus << ", epsilon+ = " << epsilonPlus
        << endl;
}

https://dyfluid.coding.net/p/OpenFOAM/d/OpenFOAM-8/git/tree/master/applications/solvers/incompressible/boundaryFoam/evaluateNearWall.H

学流体的小明

感谢李老师关注。

你用的是层流模拟么，还是什么湍流模型？没看到你用什么湍流模型。

我用的湍流模型和channel395的一样，是

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      turbulenceProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

simulationType LES;

LES
{
    LESModel        WALE;

    turbulence      on;

    printCoeffs     on;

    delta           cubeRootVol;

    cubeRootVolCoeffs
    {
        deltaCoeff      1;
    }

    PrandtlCoeffs
    {
        delta           cubeRootVol;
        cubeRootVolCoeffs
        {
            deltaCoeff      1;
        }

        smoothCoeffs
        {
            delta           cubeRootVol;
            cubeRootVolCoeffs
            {
                deltaCoeff      1;
            }

            maxDeltaRatio   1.1;
        }

        Cdelta          0.158;
    }

    vanDriestCoeffs
    {
        delta           cubeRootVol;
        cubeRootVolCoeffs
        {
            deltaCoeff      1;
        }

        smoothCoeffs
        {
            delta           cubeRootVol;
            cubeRootVolCoeffs
            {
                deltaCoeff      1;
            }

            maxDeltaRatio   1.1;
        }

        Aplus           26;
        Cdelta          0.158;
    }

    smoothCoeffs
    {
        delta           cubeRootVol;
        cubeRootVolCoeffs
        {
            deltaCoeff      1;
        }

        maxDeltaRatio   1.1;
    }
}

然后你这个槽道，u+y+是某个点上的还是底面平均之后的？

是底面平均之后的，postChannel工具会把槽道按照xz线（x和z坐标固定，y坐标变化）求平均，所以我使用postChannel输出的数据计算u+和y+，就是底面平均之后的。

学流体的小明

这个方法的重点是那个应力张量吧

symmTensor R = turbulence->devSigma()().boundaryField()[patchId][faceId];

我输出了wallShearStress，壁面切应力应该就是X方向的分量吧，和上面代码一个思路。能在paraview中看到上下壁面(patch)上的数据，话说做个平均应该就可以吧。顺便问下怎么做？

另外，wallShearStress计算出来的单位是（dimensions [0 2 -2 0 0 0 0];）也即m2/s2，而剪切应力的单位应该是（dimensions [1 -2 -2 0 0 0 0];）kg/m2/s2，刚好相差一个密度，所以OpenFOAM输出的wallShearStress还需要乘以密度，才是实际计算使用的剪切应力。

李东岳

paraview那面有个integrate滤镜，可以算出几何平均值，应该是一个单一的值。你看这种方式计算出来的$\tau_w$跟你的区别大么

李东岳

我按照你提供的信息，用RAS算了一下，结果还挺好。我下周测试下LES。LES要跑好久..

$\tau_w=0.00041199426,u_\tau=0.02029764, Re_\tau=1014$

捕获.PNG

kOmega模型的结果跟kEpsilon基本一样。这是算例文件，直接allrun，会自己处理一个图出来

wallFunction-kOmega.zip

因为channel395计算的结果，我画出来也和上面这两张图差不多，所以感觉还是正则化的问题。

自带的channel395也对不上么？

学流体的小明

以OpenFOAM自带的channel395算例为例，我这里有4种计算$u_{\tau}$的方法。

方法一：

$$ u_{\tau } = {\rm Re}_{\tau }\nu/h=395*2e-5/1 $$

得$u_{\tau } = 0.0079000$。

方法二：
过计算出来的速度-高度曲线，实际地进行计算。

$$ {u_\tau } = \sqrt {\nu {{\left. {\frac{{du}}{{dy}}} \right|}_{y = 0}}} = \sqrt {\nu {\frac{{u_1}}{{y_1}}}} $$

其中${u_1}$和${y_1}$分别是第一个输出点的流向速度和到壁面距离。
得到$u_{\tau } = 0.0069133$。

方法三：
通过paraview的filter - integrate variables - 查看 cellData，得到wallShearStress的和以及面积Area，计算固壁面上的wallShearStress平均值，再使用

$$ {u_\tau } = \sqrt {\frac{{{\tau _w}}}{\rho }} $$

计算摩擦速度$ {u_\tau }$。注意不可压缩求解器中没有密度，则认为$\rho=1$。
结果是$u_{\tau } = 0.0066952$。

方法四：
在完全发展的槽道流中

$$ \frac{{\partial p}}{{\partial x}} \times h = {\tau _w} $$

也就是压力梯度乘以槽道半高等于壁面摩擦应力。给定平均速度的情况下，OpenFOAM会调节压力梯度使得平均速度达到设定值。在channel395的算例中到计算快结束时候，选择了某一时刻的压力梯度$\frac{{\partial p}}{{\partial x}} = 4.44604e-05$。代入计算，
得到$u_{\tau } = 0.0066678$。

你看这种方式计算出来的$u_{\tau }$跟你的区别大么

可以看到方法一是一种结果，方法二、三、四的结果很接近。
现在的问题就是实际计算出来的摩擦速度和由${\rm Re}_{\tau}$给出的摩擦速度不一样，但方法一的结果才是正确的。

自带的channel395也对不上么？

channel395算例的结果
方法一的摩擦速度归一化结果是

方法二（方法三和四也是这样子的图）摩擦速度归一化结果是

channel1000的结果在一楼。
现在感觉是槽道流动还是没算对。

李东岳

234方法都应该趋向于1. 有可能自带的channel395网格需要调整，现在yplus是多少？fvScheme格式我看也可以改成liner格式。目前不确定影响多大。debug一次太费时间。

学流体的小明

channel395的平均yPlus是1.5910。channel1000的平均yPlus是0.8471。
channel395的fvScheme格式是下面这个，channel1000也没改fvScheme。

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
    default         backward;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    default         none;
    div(phi,U)      Gauss linear;
    div(phi,k)      Gauss limitedLinear 1;
    div(phi,B)      Gauss limitedLinear 1;
    div(B)          Gauss linear;
    div(phi,nuTilda) Gauss limitedLinear 1;
    div((nuEff*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

wallDist
{
    method meshWave;
}

学流体的小明

这个一维的算例也不能用LES，我试了一下，算出来的$\tau_w=5.41451763e-05$，和正确值0.00041199426差远了。
LES还得调三维算例。

李东岳

是的，LES就得算三维的，然后还需要时间平均，这一算起来debug起来就需要时间了。网格这个问题，我看EugeneDeVilliers说没问题。应该不是网格的事。channel395是EugeneDeVilliers做的。这个是他的博士论文。

https://www.jianguoyun.com/p/DT7jUyMQ9s3ZBhip1uYDIAA

学流体的小明

对，我最开始做槽道流算例的时候，用的他的初始化流场的方法，不然一直算的都是层流。不过chanel1000的这个不用特地初始化也可以。
下面这个是我的channel1000的算例，LES，就是在channel395的基础上改的，OpenFOAM版本是v2012。Allrun脚本可以直接运行，画图用matlab文件夹里面的boundary_layer_profile_2.m就行，改一下要画哪个时刻的就行。
各位老师同学有兴趣可以算一算，看看是哪里出问题了。帮帮我吧🙏
https://jbox.sjtu.edu.cn/l/L1Asg4 (提取码：1234)

学流体的小明

找到一个帖子，对channel395分析很多
https://www.cfd-online.com/Forums/openfoam-solving/155534-les-channel-flow-data-case-files-technical-report.html
我下载了他的数据，他是三套网格，我都用方法二画了速度分布图，M1的结果和我画出来的问题类似，M2和M3就都比较好了。
为什么$u_\tau$算不对？结论：还是网格不够细。

李东岳

结论：还是网格不够细。

还是要深入研究下。我看默认的算例6万网格，EugeneDeVilliers博士论文里面说6万网格没有问题也有数据。

不过，如果确定6万网格数据不对的话。可以找openfoam提bug让他们处理，这属于他们的bug了。

李东岳

Re_tau1000的k、epsilon这两个场的数据你有不，或者雷诺应力的

coolhhh

@学流体的小明你好，这个下载链接需要校内VPN，外校无法下载，想问下能提供个外校人员可下载链接吗？

对，我最开始做槽道流算例的时候，用的他的初始化流场的方法，不然一直算的都是层流。不过chanel1000的这个不用特地初始化也可以。
下面这个是我的channel1000的算例，LES，就是在channel395的基础上改的，OpenFOAM版本是v2012。Allrun脚本可以直接运行，画图用matlab文件夹里面的boundary_layer_profile_2.m就行，改一下要画哪个时刻的就行。
各位老师同学有兴趣可以算一算，看看是哪里出问题了。帮帮我吧🙏
https://jbox.sjtu.edu.cn/l/L1Asg4 (提取码：1234)

学流体的小明

找到一个数据库
http://turbulence.pha.jhu.edu/Channel_Flow.aspx

学流体的小明

@coolhhh 不好意思哈，以前也没怎么用过学校的网盘，不知道还需要VPN。放到百度网盘了。
这次画的网格更细了，应该会有一个比较好的结果，我也正在算。

链接：https://pan.baidu.com/s/1C89EfgcZAxFhADNkMEqSNA?pwd=tzkd
提取码：tzkd

学流体的小明

@coolhhh @李东岳
昨天仔细读了一下下面的这个文章，作者对OpenFOAM计算LES槽道进行了比较系统的分析。我发现他的图也是，只有最细的那套网格才算出了非常好的速度剖面，比较粗的两套网格都没算好，形状也和我的问题一样。
这篇文章还是在https://www.cfd-china.com/topic/2121/q-dns计算槽道流遇到了一些问题-求大神们指点看到的。之前读过，但是没仔细读😂。然后就浪费了好多好多时间。
确实是网格不够细的问题。

https://www.sciencedirect.com/science/article/pii/S0021999117304059?via%3Dihub

coolhhh

@学流体的小明收到，非常感谢

CFD中文网

LES直流槽道边界层模拟，如何得到正则化速度u+以及正则化坐标y+？