CFD中文网

@Vortex 前辈您好，请问可以向您请教一下LESGO温度分层算例的相关设置吗？

@xpc 在 如何在动量源项中用一个关于时间的速度函数去驱动流体运动？ 中说：

各位老师好，不知道大家是否有人在fvOptions动量源项meanVelocityForce中的某个方向加入过一个函数去驱动整个流体的运动？
感觉这个方法不太行，自带算例的Ubar只能是常数，大家有没有人做过类似的工作？或者这种问题一般如何解决？谢谢大家！

momentumSource
{
    type            meanVelocityForce;
    active          yes;
    selectionMode   all;
    fields          (U);
    Ubar    ((cos(2*pi/(12*3600)*time())+1)    0     0  );

用vectorSemiImplicitSource应该可以实现

momentumSource 
{   
    type            vectorSemiImplicitSource;
    

    selectionMode   cellSet;
    cellSet         UpperFluid;

    volumeMode      specific;

    injectionRateSuSp
    {
        U       ((0.0016 0 0) 0);
    }

}

@Vortex 在 中性大气环境湍流动能的自保持 中说：

@李东岳 介绍个专门做ABL的code，叫LESGO (https://lesgo.me.jhu.edu/)，Meneveau大佬以前是拿它研究LES模型的。不稳定、中性、稳定的ABL都可以做，还能考虑科氏力。

想问一下LESGO有现成的考虑浮力驱动的求解器吗？

@落花风 抱歉我暂时没用python画过这么复杂几何+非结构网格的云图....我的算例是基于均匀网格的所以相对比较好处理

https://fluidfoam.readthedocs.io/en/latest/
我一直都是用这个包来读OpenFOAM cell的数据来做后处理的，非常好用

你有热通量输入到计算域里，但是没有边界让输入的热量离开，系统的温度无法达到统计稳态那就自然随着时间不断提高了

如果你求解的是瞬态问题且不存在稳定解（不是湍流那种统计稳态），那求解结果肯定不满足稳态NS方程。

瞬态求解你可以理解为下一时间步的流场$U_{n+1}$是未知的，你只能通过$U_{n}$和离散NS方程结合来求解出$U_{n+1}$。假如你$U_{n}$的流场已经满足稳态NS方程了，那必有$U_{n+1}-U_{n} = 0$，也就是下一步的流场和本步完全一样，流动不随时间变化，这就是前面所说的稳定解。层流倒是可能会出现这种情况。所以是否能满足还是针对你考虑的物理问题。

不懂就问，隐式大涡模拟和在粗网格下算DNS有区别吗？

建筑环境湍流那边很多人都这么算的，比如

Anderson W. Amplitude modulation of streamwise velocity fluctuations in the roughness sublayer: evidence from large-eddy simulations[J]. Journal of Fluid Mechanics, 2016, 789: 567-588.

@李东岳 非常感谢！我也稍微想通一点了，壁面函数的作用只会出现在边界处的nut上，时间步第n步通过壁面函数的到的patch上的nut，会在第n+1步的动量方程粘性项求解中起作用 所以第一层网格中心上的nut只要通过k求解得到就好了，不需要被壁面函数更新。希望这次没有理解错了

本人最近在学习OpenFOAM中的壁面函数。因为自己涉及到更多的是大涡模拟的计算，所以看的基本上是关于湍流运动粘度nut相关的内容。大概把理论部分过了一遍然后好奇OpenFOAM具体是如何执行的，就在几乎零C++认识的基础上看了下OpenFOAM的代码，发现有特别多不懂的地方。

一开始自己的疑问是，在大涡模拟中并没有与nut相关的方程，那么经过壁面函数计算得到的nut是如何进入到方程时间推进的计算中的？
看了下自己用的最多的 one-equation sgs tke 模型 kEqn的代码correct()

void kEqn<BasicTurbulenceModel>::correct()
{
    if (!this->turbulence_)
    {
        return;
    }

    // Local references
    const alphaField& alpha = this->alpha_;
    const rhoField& rho = this->rho_;
    const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
    const volVectorField& U = this->U_;
    volScalarField& nut = this->nut_;         //nut在这里定义
    fv::options& fvOptions(fv::options::New(this->mesh_));

    LESeddyViscosity<BasicTurbulenceModel>::correct();

    volScalarField divU(fvc::div(fvc::absolute(this->phi(), U)));

    tmp<volTensorField> tgradU(fvc::grad(U));
    volScalarField G(this->GName(), nut*(tgradU() && dev(twoSymm(tgradU())))); //在这里用到了nut，并作为计算源项G的一个输入参数
    tgradU.clear();

    tmp<fvScalarMatrix> kEqn
    (
        fvm::ddt(alpha, rho, k_)
      + fvm::div(alphaRhoPhi, k_)
      - fvm::laplacian(alpha*rho*DkEff(), k_)
     ==
        alpha*rho*G  //生成项G在这里被使用
      - fvm::SuSp((2.0/3.0)*alpha*rho*divU, k_)
      - fvm::Sp(this->Ce_*alpha*rho*sqrt(k_)/this->delta(), k_)
      + kSource()
      + fvOptions(alpha, rho, k_)
    );

    kEqn.ref().relax();
    fvOptions.constrain(kEqn.ref());
    solve(kEqn);
    fvOptions.correct(k_);
    bound(k_, this->kMin_);

    correctNut();  // 求解sgs k 方程结束后，nut被更新，函数结束
}



void kEqn<BasicTurbulenceModel>::correctNut()
{
    this->nut_ = Ck_*sqrt(k_)*this->delta(); // 新一步得到的sgs k用来计算cell center上的nut
    this->nut_.correctBoundaryConditions();  // 基于wall function更新边界上nut的数值
    fv::options::New(this->mesh_).correct(this->nut_);

    BasicTurbulenceModel::correctNut();
}

由上面的代码看来，nut在边界处的数值似乎并没有进入到矩阵求解的系数中（个人的理解，如有不对烦请指正），但从下面这段代码可以看到，calcNut()是在更新壁面上的nut数值，而并没有修改第一层网格中心上nut的数值

tmp<scalarField> nutkWallFunctionFvPatchScalarField::calcNut() const
41 {
42 const label patchi = patch().index();
43
44 const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
45 (
46 IOobject::groupName
47 (
48 turbulenceModel::propertiesName,
49 internalField().group()
50 )
51 );
52
53 const scalarField& y = turbModel.y()[patchi];
54 const tmp<volScalarField> tk = turbModel.k();
55 const volScalarField& k = tk();
56 const tmp<scalarField> tnuw = turbModel.nu(patchi);
57 const scalarField& nuw = tnuw();
58
59 const scalar Cmu25 = pow025(Cmu_);
60
61 tmp<scalarField> tnutw(new scalarField(patch().size(), 0.0));
62 scalarField& nutw = tnutw.ref();
63
64 forAll(nutw, facei)
65 {
66 label faceCelli = patch().faceCells()[facei];
67
68 scalar yPlus = Cmu25*y[facei]*sqrt(k[faceCelli])/nuw[facei];
69
70 if (yPlus > yPlusLam_)
71 {
72 nutw[facei] = nuw[facei]*(yPlus*kappa_/log(E_*yPlus) - 1.0);
73 }
74 }
75
76 return tnutw;
77 }

我个人的理解是，可能在经过壁面函数更新后，boundary上的nut经过某段代码赋值给第一层网格中心点上的nut，使得经过壁面函数修正的nut得以进入迭代计算中。可惜的是本人暂时还没找到执行上述猜想的代码，那表明OpenFOAM实际的执行方式与我理解的不同？

知乎上一位前辈似乎也提到了相同的疑问
https://zhuanlan.zhihu.com/p/32520364

这么说我的猜测的对的咯？那为什么我copy其中一个值，在internal field里搜不到？难道数据结构跟我想的不一样？那internal field里对应点上的nut值又是什么？当计算到最靠近壁面的那个cell时，到底应该用internal field里的nut还是boundary上的nut？

我目前的猜测，对于最靠近壁面那个点，其实有两个值，一个值是根据正常turbulence model算出来的，储存在该变量的internal field里；另一个是根据wall function算出来的，储存在该变量的boundary field里。然后真正计算矩阵的时候，会检查一下有没有用wall function，假如用了，就忽略internal field里的值，把wall function算出来的值覆盖上去。这里没有直接填充是为了数据处理安全，反正数据多了可以不用，少了就麻烦了。

但好像问题还是没有被解决，毕竟没有贴出来相应模块的代码......想请各位对这方面有深刻理解OFer指教指教，非常感谢！

我自己用python做后处理比较多，你可以先用下面这个python的库

https://github.com/fluiddyn/fluidfoam

把OpenFOAM数据读进程序然后写为别的格式，如二进制或者hdf5之类的，再用Matlab读新生成的数据文件

想请问一下各位OFer，通过命令

postProcess -func "grad(U)"

得到的速度梯度张量共有九个元素，打开数据文件后发现这九个元素是以一维数列的形式保存的，那这九个元素的排列顺序是什么呢？Google了大半天都没找出来.....
之前自己默认顺序是dudx,dudy,dudz,dvdx,dvdy,dvdz,dwdx,dwdy,dwdz，但算出来的结果有点怪，所以怀疑自己是不是在排列顺序上理解错了

buoyantBoussinesqPimpleFoam在OpenFOAM-v2006的版本里还存在

@dzw05 你遇到的极值局部点距离给定压力条件的边界接近吗？我今天下午算的时候发现时间步过大时，距离fixedValue p_rgh边界最近的cell中心也会出现温度极值....调小了时间步后就没有了。我算的这段时间一直感觉温度和压强的联系太密切了...

https://cfd-china.com/topic/4194/求助-buoyantboussinesqpimplefoam-遇到的非物理温度边界/6?_=1605935257119

我之前用buoyantBoussinesqPimpleFoam遇到了完全一样的情况，后来发现是在fvSolution里面pimple的部分给了一个pRefCell和pRefValue，每次都是这个点温度莫名增高和降低。

现在我的解决办法是找到一个边界直接给定一个压力边界值，这样就不需要调用这两个给出莫名其妙结果的设置了。

如果有什么更好的解决办法可以一起交流一下。

@benqing 谢谢你的解答！我再想想其他解决办法

各位OFer：

请问是否可能在一个算例中对速度的inlet/outlet采用cyclic的边界条件，但对其他变量如pressure采用zeroGradient的边界条件呢？Google了一圈好像都没有发现理想的解决办法...

@东岳 是的，一开始我也以为是网格的问题，甚至试过用blockmesh生成一个长方体计算域（四周是cyclic，上下分别为slip和noslip），里面都是立方体的网格，但还是会出现相同的情况。

但是我翻回自己学这个solver时测试的一个所谓"RB对流"算例，长方体的计算域每个面都是noSlip，通过上下表面温度差来驱动流动。这个算例里的温度场就非常正常，全场的温度最大最小值分别是下壁面和上壁面的温度....所以我怀疑会不会时cylic边界的问题？但是cyclic和cyclicAMI我都尝试了，似乎问题还是在。