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

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

李东岳

Komen那个是的

但是我看EugeneDeVilliers博士论文那个用的6万网格计算的u+y+还可以（图5.5）

我觉得这个是很重要的内容，如果openfoam官方的算例，出现了网格分辨率不够导致的错误，这个bug必须应该处理

hongjiewang

@学流体的小明先前我做过这个工作，用于大涡模拟入口的前置算例，我当时画u+--y+使用的你说的方法二得到的，效果还行

hongjiewang

@学流体的小明

对于此，我想请问，Um即平均速度如何给定，雷诺数的计算，假如我是在半宽通道中跑的，即下表面为wall，上表面为symmetry,此时，雷诺数公式里面应该把2h换成水力半径吗

coolhhh

@李东岳在 LES直流槽道边界层模拟，如何得到正则化速度u+以及正则化坐标y+？中说：

Komen那个是的

但是我看EugeneDeVilliers博士论文那个用的6万网格计算的u+y+还可以（图5.5）

我觉得这个是很重要的内容，如果openfoam官方的算例，出现了网格分辨率不够导致的错误，这个bug必须应该处理

李老师，EugeneDeVilliers博士论文图5.5，图名写normalised by the DNS shear velocity，是否有可能正则化用的$U_{\tau }$不一样？

学流体的小明

@hongjiewang
Um是通过在动量方程中直接添加源项实现的，在v2012版本中是：

/*--------------------------------*- 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);
}


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

程序会自动调整压力梯度以达到Ubar的预设值。
之后我是打算用给定的动量方程源项，相当于给定$\partial p / \partial x$，通过vectorSemiImplicitSource实现。也是在fvOptions中调用。4800是随便给的，和上面0.4m/s的平均速度没有联系哈。$\partial p / \partial x$和$\tau_w$的关系见方法四。
有一些参考：
http://xiaopingqiu.github.io/2016/03/20/fvOptions2/
https://caefn.com/openfoam/fvoptions-semiimplicitsource

momentumSource
{
   type vectorSemiImplicitSource;
   active on;
   

   vectorSemiImplicitSourceCoeffs
   {
      selectionMode all;
      volumeMode        specific; //absolute
      injectionRateSuSp
      {
         U           ( (4800 0 0) 0); //partial p / partial x
      }
   }
}

学流体的小明

那EugeneDeVilliers博士论文的结果也是很好的，你看我方法一算出来的$u_{\tau}$，也相当于DNS的结果吧，画出来之后粘性底层完全对不上。
而且就算真的有DNS结果，我画出来的曲线的形状，也无法和理论解对上。就关键是形状它就对不上。

学流体的小明

@hongjiewang
槽道的雷诺数公式都是用两倍半高来，不管实际模拟了多少的。
摩擦雷诺数是用的一倍半高。

其中h为槽道半高，uτ为摩擦速度，ν为运动粘度。根据Pope在Turbulent Flows一书中提到的公式

这个公式里面的Re是用的两倍半高，Re_tau是一倍半高。

李东岳

@coolhhh

normalised by the DNS shear velocity

我觉得这句话提供了很多信息，他们应该是用方法1计算的$u_\tau$。

如果这样的话，应该是个小bug，就是网格分辨率不够。

很难相信会留这么一个小bug在里面

李东岳

@学流体的小明在 LES直流槽道边界层模拟，如何得到正则化速度u+以及正则化坐标y+？中说：

通过vectorSemiImplicitSource实现

你可以试一下这个，我直觉感觉这个可能会导致流速持续的增加下不来了

学流体的小明

那不会的，我最开始就是用vectorSemiImplicitSource的，另外一个算例的正则化的结果和

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

类似。您说的流速持续增加很可能是压力梯度的具体数值没设置好，设置得大了。我之前设置的小了，发现平均速度、以及采样的某个点的速度总是在下降。就是驱动力不够，有耗散，速度就降下来了。

槽道的平均速度是会趋于稳定的。
$h=0.005m$，$U=4.0m/s$

李东岳

在动量方程添加一个固定的源项，直觉理解会一直往上增加速度。你这么一说，那就是要设置的刚刚好，这个设置的压力梯度与本身的压力梯度相抗衡。好像是需要这样。不能大、不能小。有道理啊！

coolhhh

@李东岳我用OpenFOAM-8也算了channel395，用tutorial算例计算，采用 @学流体的小明提供的MATLAB程序中method2算的$U_\tau$：

tutorial算例用的WALE模型，0文件夹其实不需设置k和nuTilda，只需要设置p、U和nut。
tutorial算例是给定了一个湍流初始场。一开始没注意删掉了，没有初始场到1000s计算不出湍流。初始湍流场的不同对结果是否影响会比较大？
一般认为开始计算的一段时间的湍流还没充分发展，计算平均场是都会剔除。tutorial算例提供的统计是从0时刻就开始统计。不过我尝试从1000s才开始统计，共计算5000s，结果与默认算例结果一致。
同时还提取了EugeneDeVilliers博士论文图5.5中6万网格的结果，其实跟我们的计算结果一致的，都是有偏差！推测他也是用method2算的。他的论文没有给出参考线，导致看起来结果很好。图中Spalding曲线是根据《无痛苦NS方程笔记》中提到的公式画出对比。

(1)tutorial算例默认设置结果：从0秒开始统计

u+ - y+.png

(2)从1000秒开始统计。1000s的结果其实是1000s这个时刻的瞬时结果，但对整个底面平均后，与后面时刻的结果基本一致。
u+ - y+.png

4. tutorial算例默认设置，method2计算结果都是在粘性子层结果较好，log区结果不吻合。尝试$U_\tau$按照method1取值，即$U_\tau=0.0079$，结果如下。说明了目前算例结果，改变$U_\tau$，只能满足粘性子层或者log区之一，是无法同时满足这两个区。想要同时满足，可能需要尝试加密网格之类方式
u+ - y+.png

coolhhh

@学流体的小明 @李东岳不好意思，刚再仔细看了下，我的getdata数据好像有点问题，我再核实一下

coolhhh

@学流体的小明 @李东岳
很抱歉撒上面的getdata数据不小心搞错了，更正如下

先说结论，EugeneDeVilliers博士论文图5.5中6万网格的结果，是按照method1计算的。
按照method1，$U_{\tau}=0.0079$，结果如下图所示。EugeneDeVilliers结果也是在粘性子层偏差较大，和method1计算结果基本一致
tutorial算例默认设置结果：从0秒开始统计，结果更正如下：
tutorial算例从1000秒开始统计，结果更正如下

学流体的小明

槽道湍流是无法一直流动下去的，它里面总会有耗散，所以需要施加驱动。
施加驱动的两种方式本质上都是给一个压力梯度，稳定之后，压力梯度就平衡了壁面摩擦力。
其实是这么来的：

1959afe8-18ce-45e1-8206-0f7c65a9f830-9c197790f8b2d16fe0b0eafc7d86828.jpg

也不是不能大、不能小，设置得大了，槽道流速增加，会达到新的平衡的。

不能大、不能小。

学流体的小明

做的真好，谢谢你的探索和分享
所以EugeneDeVilliers结果也不是那么好，只是他没给出参考线😂

李东岳

嗯，经过你们2个的验证，应该就是自带tutorials的问题了。

另外压力梯度那个，我最近在方程里面加了科氏力，也是一个规规矩矩的槽道，只不过尺度非常大（大气层）。目前算出来的结果就是速度持续的在旋转稳定不下来，与类似压力梯度的效果类似，我深入研究下咋回事，看看你那个图能不能解释一下

李东岳

http://dyfluid.com/boundaryFoam.html 我顺着这个里面的方程2，在y方向做积分，继续推了一下：

$$
\int\left(-\frac{\partial \tau_{xy}}{\partial y}\right)\rd y=-\int\left(\frac{\partial p}{ \partial x}\right)\rd y
$$

$$
\tau_{xy}^h -\tau_{w}=-\tau_{w}=\frac{\partial p}{ \partial x}h
$$

所以充分发展的管道，如果壁面剪切力知道的话，应该这样给一个压力梯度。

也不是不能大、不能小，设置得大了，槽道流速增加，会达到新的平衡的。

对，应该是。还得持续学习啊！

李东岳

我最近在方程里面加了科氏力，也是一个规规矩矩的槽道，只不过尺度非常大（大气层）。目前算出来的结果就是速度持续的在旋转稳定不下来

我发现添加压力梯度跟平均速度，在我这里还不太一样。添加速度的话，如果是x方向速度，因为科室力会产生旋转，x方向速度一直不变，但是y会增加，导致速度增加。压力梯度却没有。因为跟本帖关系也不大，我先研究一下，然后新开一个帖子讨论。

李东岳

@李东岳 https://www.cfd-china.com/topic/6259/添加动量方程源项导致速度持续增加

luofq-sysu

@学流体的小明感谢此贴，感谢几位的研究，解答了我的很多疑问。
并请教：槽道流LES计算的入口inflow速度场影响大吗？
因为我近期也开展了多个雷诺数Re180~5200的槽道流计算，但计算得到的壁面摩擦速度偏小（大约只有目标值的一半）。
对比了一下你的channel1000案例文件，发现一个主要的差别是初始场0/U，看起来你是采用了一种生成脉动速度的入口条件吗？

学流体的小明

@luofq-sysu

并请教：槽道流LES计算的入口inflow速度场影响大吗？

我的计算域没有入口出口之分，沿流动方向的两个边界都是周期性的边界条件，展向也是周期性的边界条件。所以无法回答你的问题。

看起来你是采用了一种生成脉动速度的入口条件吗？

是的，用的就是EugeneDeVilliers在他博士论文中提到的初始化槽道流流场的方法。

luofq-sysu

谢谢。@学流体的小明

我的计算域没有入口出口之分，沿流动方向的两个边界都是周期性的边界条件，展向也是周期性的边界条件。所以无法回答你的问题。

我的算例也是周期边界的槽道流，其实应该是internalField初始场的影响。
前面我采用的是uniform的初始速度场，计算湍流雷诺数远达不到目标值。最近按你算例的codeStream代码加入了初始脉动速度，计算雷诺数比较接近目标值了。
以Re_tau=1000为例，我现在计算得到的雷诺数数值是910

是的，用的就是EugeneDeVilliers在他博士论文中提到的初始化槽道流流场的方法。

EugeneDeVilliers大佬的论文提到了这种湍流需要加入初始扰动才能发展，这个初始化方法也有一个 perturbU的开源代码。不过codeStream的实现更方便

李东岳

最近按你算例的codeStream代码加入了初始脉动速度

这个codeStream写的扰动，可以保证$\nabla\cdot\bfU=0$么

luofq-sysu

这个codeStream写的扰动，可以保证$\nabla\cdot\bfU=0$么

我只看到codeStream跟perturbU的代码是一致的，具体的公式在EugeneDeVilliers博士论文5.1.2章也能找到