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的版本里还存在