nut壁面函数如何影响湍流模拟
-
本人最近在学习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指教指教，非常感谢！
-
@chszkc 在 nut壁面函数如何影响湍流模拟 中说：
在大涡模拟中并没有与nut相关的方程，那么经过壁面函数计算得到的nut是如何进入到方程时间推进的计算中的？
$\nu_t$在动量方程中起作用。在这一项里面$\nabla\cdot((\nu+\nu_t)\bfU\bfU)$，不在湍流模型里面。湍流模型只是求解。
nut在边界处的数值似乎并没有进入到矩阵求解的系数中（个人的理解，如有不对烦请指正）
$\nu_t$不存在PDE方程，因此你说的是对的。
boundary上的nut经过某段代码赋值给第一层网格中心点上的nut，使得经过壁面函数修正的nut得以进入迭代计算中。
不是。
网格第一层的$\nu_t$是通过$\varepsilon,k$来求解的，跟边界没有关系。外链网站的内容太多，没细看，暂不回答咯总体来讲，壁面函数的流程主要是手动计算壁面第一层的epsilon以及湍流产生项（在OpenFOAM里面一般称之为G），然后$\nu_t$通过代数公式求出即可。更详细的我更新在这里了 http://www.dyfluid.cn/theory.pdf http://dyfluid.com/docs/wallFunction.html
-
@李东岳 在 nut壁面函数如何影响湍流模拟 中说：
@chszkc 在 nut壁面函数如何影响湍流模拟 中说：
在大涡模拟中并没有与nut相关的方程，那么经过壁面函数计算得到的nut是如何进入到方程时间推进的计算中的？
$\nu_t$在动量方程中起作用。在这一项里面$\nabla\cdot((\nu+\nu_t)\bfU\bfU)$，不在湍流模型里面。湍流模型只是求解。
nut在边界处的数值似乎并没有进入到矩阵求解的系数中（个人的理解，如有不对烦请指正）
$\nu_t$不存在PDE方程，因此你说的是对的。
boundary上的nut经过某段代码赋值给第一层网格中心点上的nut，使得经过壁面函数修正的nut得以进入迭代计算中。
不是。
网格第一层的$\nu_t$是通过$\varepsilon,k$来求解的，跟边界没有关系。外链网站的内容太多，没细看，暂不回答咯总体来讲，壁面函数的流程主要是手动计算壁面第一层的epsilon以及湍流产生项（在OpenFOAM里面一般称之为G），然后$\nu_t$通过代数公式求出即可。更详细的我更新在这里了 http://www.dyfluid.cn/theory.pdf http://dyfluid.com/docs/wallFunction.html
李老师，会计算靠近壁面第一层网格的epsilon和湍流产生项吗？我在OpenFOAM中没有发现相关代码，希望老师点一下
-
@AppleKiller 附近壁面一层网格的生成方法我那个笔记里面有写
所以第一层网格中心上的nut只要通过k求解得到就好了，不需要被壁面函数更新。希望这次没有理解错了