CFD中文网

@cccrrryyy
老师，实现这个的方法是，将我们算出来的粘性值，通过DEFINE_PROPERTY宏，加入到物性的分子粘性中，这样，在计算动量方程时用到的粘性，就是有效粘性（分子粘性+湍流粘性）了
后面要和燃烧模型耦合会有点麻烦，因为燃烧模型会用到湍流粘性，我估计fluent会直接去找他们自己模型算的湍流粘性值，不知道这里会不会需要做一些额外的处理

@cccrrryyy
谢谢老师，目前算是解决了。我后来不把UDF挂载在k-epsilon模型旁边，而是直接把这个UDF作为湍流模型计算，fluent选择laminar模式，解决了几个bug之后，现在算的结果还可以，只是和OpenFOAM算的结果有一定数量上的差异。
目前看来fluent植入自己的新湍流模型是可行的，但是不建议这么做，调试起来很麻烦，不是项目要求我也不会做这个

最近在开发DES模型，看到李老师的《无痛苦N-S方程笔记》中好像有一处笔误，前来po个帖子。
kOmegaSSTDES模型中，k方程源项有一个系数$F_{DES}$
《笔记》写的表达式是：
$F_{DES}=\mathop{\max}(\frac{L_t}{C_{DES}\Delta}(1-F_{DES}),1)$
这个系数的表达式应该是：
$F_{DES}=\mathop{\max}(\frac{L_t}{C_{DES}\Delta}(1-F_{SST}),1)$
其中，$F_{SST}$有三个选择，0、$F_{1}$、$F_{2}$。这里的$F_{1}$和$F_{2}$是SST模型中的函数，可以在SST模型解释中找到，原文和OpenFOAM SSTDES模型程序如下


论文截图

OpenFOAM-9程序截图

——240430 孙逸凡

我补充一句，最近也在写湍流模型UDF，需要用到壁面距离，发现C_WALL_DIST不能用之后不信邪，测试了很多（一个多小时测了十几二十个case），终于测试出来发现：
C_WALL_DIST只有在选择了fluent已有湍流模型（我选了k-omega）的时候才能用，而当选择了laminar模式，用自己的湍流模型UDF的时候，如果涉及到C_WALL_DIST，就会报错。我的fluent版本为2022R2

各位老师，我在进行fluent UDS计算时遇到了很奇怪的问题。算的是一个对称半域的平面射流，从左下角射入，右侧和上边是压力出口，中间对称边界，速度和湍动能分布如下


然后目前植入的UDS方程是这样的

很常规的输运方程，通过UDF自定义扩散系数加入了湍流扩散
然后源项系数S的分布是这样的

和湍动能很像，上游剪切层的值很大
但是算出来的UDS结果却没有这种明显特征（射流入口的值为0.2）

出口处的诡异分布就先不管了，但是上游剪切层源项S那么大的区域，UDS算出来就好像啥也没有似的，另外的混合层算例也遇到过这样的情况，就很迷，源项明明很大，10的4、5次方，但是phi的分布平平无奇，一点也没有遵循源项的分布规律
所以想请问下各位老师，是我的理解有问题么，源项很大不一定导致结果很大？我有一个推测是fluent会不会因为UDS源项过大，自己做一些用户看不到的调整，不会把这么大的源项全给UDS方程

我之前一直在OpenFOAM做这个模型的开发，植入的方程和这个差不多，因为UDS的一些限制做了一些阉割，不过导出过源项系数S看过和fluent算的差不多，phi也是和湍动能很类似的分布，而不是fluent UDS算出来的这种很平均、上游不够大的情况，这也是我很怀疑是fluent有自己操作的原因，但是还没有实锤，因为网上搜不到类似的问题，我在想有没有办法能验证这种猜测

@李东岳 谢谢老师
老师说的pEqn.H里的更新，是指phi = phiHbyA - pEqn.flux()对吗，这个是SIMPLE算法里压力修正完了更新一下速度，我理解的没错吧老师（您的simpleFoam解析我可是学了很多遍 XD ）
那下一次迭代的时候，UEqn里用的phi用的就是上一次pEqn里更新过的，是吗老师

各位大佬，请教一个代码问题：关于求解器在什么时候计算通量phi的值
以simpleFoam为例，在组装速度方程系数矩阵的时候，散度项会用到phi：

tmp<fvVectorMatrix> tUEqn
(
    fvm::div(phi, U)
  + MRF.DDt(U)
  + turbulence->divDevSigma(U)
 ==
    fvModels.source(U)
);

phi是在createPhi.H里初始化的，代码如下：

surfaceScalarField phi
(
    IOobject
    (
        "phi",
        runTime.timeName(),
        mesh,
        IOobject::READ_IF_PRESENT,
        IOobject::AUTO_WRITE
    ),
    fvc::flux(U)
);

fvc::flux(U)是确定了phi的一个计算方式？那，只要用到phi的地方就会更新phi的值吗，还是有什么代码会暗中计算一下，如果我在程序里写phi=fvc::flux(U)，也能正常计算是吗
主要是现在我需要在求解器里加入另一个速度的方程，需要计算这个速度的通量，甚至在主流U的方程里就会用到另一个速度的通量，但是通量的计算结果不太对，每隔一个数就会有一个0

tmp<fvVectorMatrix> tUEqn
(
    fvm::div(phi, U)
  + MRF.DDt(U)
  + turbulence->divDevSigma(U)
 ==
    fvModels.source(U)
  + fvc::div(phi2, U2)
);

类似上面这样的形式，直接在方程里加一个另一个速度的散度项，这样计算出来的phi2就不太对，请问这样直接引用phi2对吗各位老师，还是我需要在速度方程之前自己算一下这个phi2的值

各位老师，小生不才，在此请教一个问题，最近在算一个封闭长方体内的受限射流，15m/s的来流速度，射流孔径2.6mm，腔体大小70mm*70mm，网格用fluent meshing生成，形状如下：
网格最密的地方0.2mm，最粗0.8mm
使用的是Smagorinsky模型，pisoFoam瞬态求解器，dt=1e-7s，边界条件：【速度函数给的进口，pressureInletOutletVelocity出口】，【压力固定值0进口，totalPressure出口】，在之前算别的射流用的这一套
但是计算中总是会出现这样的情况：
使用RANS可以算，但是算的结果扩散不强，只是不会在入口附近发散。推测是网格不太好的原因（是入口太密，下游过于稀疏？），准备用ICEM画一个。
不知道各位大佬有没有在计算中遇到过类似的问题，或者有没有什么好的建议，在此先谢过大家花时间看我的帖子了

@evensun 哦哦如果从随机数来看，应该是白噪声类型的随机速度脉动是吗老师

老师您好，请问你这种入口脉动速度的给定是通过能量谱导出的吗，最近我在复现文献，但是对文献给出的能量谱怎么得到脉动速度不太懂

各位大佬大家好，我最近准备在rhoSimpleFoam的基础上加入自己的算法，算法中比较关键的点是有另一组速度、压力和温度场（暂且叫它们Uf，pf，Tf）。

之前在simpleFoam中顺利地植入了，但是在可压缩这里我发现：thermo.rho()这些函数在调用时，读取的物理量是求解器原来的U、p，如果我再另外定义Uf，pf的话，好像不能直接使用thermo.rho()函数进行另一组速度压力场的密度修正？

附上一段src/thermophysicalModels/basic/fluidThermo/fluidThermo.C中的代码，麻烦大佬们看下是不是这样的

Foam::fluidThermo::implementation::implementation
(
    const fvMesh& mesh,
    const word& phaseName
)
:
    p_(lookupOrConstruct(mesh, "p"))
{}


Foam::fluidThermo::implementation::implementation
(
    const fvMesh& mesh,
    const dictionary& dict,
    const word& phaseName
)
:
    p_(lookupOrConstruct(mesh, "p"))
{}

@李东岳

老师，我的求解器在计算时设置为层流，然后湍流的方程和粘性在求解器中解算。不知道是不是因为设置为层流的原因，wallshearstress在后处理的时候提示unable to find turbulence model in the database，所以没法用postProcess算

各位老师打扰了，本人才疏学浅学艺不精，最近在后处理，花了好久终于看懂了网上解释壁面剪切应力的帖子。
由于使用的湍流模型得到的湍流粘性系数为张量，在求解器中编写后处理程序时 ，遇到了不知道该如何输出应力在边界上的值的问题：

先前看到，在OpenFOAM源码的wallShearStress.C中，壁面剪切应力由以下方法求得：

wallShearStressBf[patchi] = (-Sfp/magSfp) & tau.boundaryField()[patchi];

其中的tau，其实就是有效雷诺应力：

tau = model.devSigma();

但由于我的后处理程序植入在solver中（直接写bool不太会），因此自己在createFields里定义了一个应力张量场stress，然后把层流的粘性应力：

dimensionedScalar("nu", dimViscosity, 1.388e-5) * dev(twosymm(fvc::grad(U)))

和湍流粘性应力相加得到，经检验可以输出内部场（通过stress[cellID]）的应力值。

正准备和法向单位向量相乘，发现输不出来边界上的stress.boundaryField()[faceID]

我知道我这样的处理肯定是过于简单和理想化了（简直就是瞎搞，瞎搞都算不上）只能盼望各位老师赐教了，该如何得到边界上的值呢

@stardust
请问楼主写完了嘛，如果可以的话可不可以贴上来让大家学习一下 （最近正好在做一个类似的东西）

@cccrrryyy
当然老师您可能会问为什么非定常问题还要用稳态求解器，这个的话可能主要是这几点原因：

一是原有的经典模型如k-epsilon等对这些问题可以用稳态算法进行计算，通过模型理论中的平均方法得到稳态解（相关问题没有做足够的调研，应该也有使用瞬态求解器的）；

二是我现在使用的模型的提出者基本上都使用稳态算法进行求解的；

三是我初步使用了瞬态求解器进行计算，但不稳定，会发散（当然也有可能是算法和设置有问题导致的）

@cccrrryyy
老师说的是，我没有说太清楚。

毕竟所有的实际问题在严格的意义上都不是定常的。如果广义上这个流场的脉动比较小，那么理论上我们用一个稳态的求解器的话，是可以算出来一个残差很小的收敛解的，老师您看我这样理解是对的吗。

然而同时对于一些问题，虽然实际上并不能认为是一个严格的定常问题，但是我们现在使用的很多湍流模型用稳态求解器计算却能得到一个稳态的解，比如平面射流。

那么这个时候，对这样的非定常问题，如果一个湍流模型在稳态求解器的基础上，算得的结果呈现出了涡的发展变化过程，老师您是怎样看这个结果的呢

想请教各位大佬一个有关稳态求解器与湍流的问题

就是理论上，稳态的求解器应该算出来一个最终残差很小的稳定解对吗

但是如果使用的湍流模型说，自己会解算出湍流的拟序结构（according to 它的提出者），所以使用稳态求解器算出来的解，会呈现出一个脉动的特征（比如涡的变化发展等过程），以至于没有一个收敛的解...

请问各位大佬是怎么看这样的结果的呢，这样是可以理解或接受的吗，又或请问有没有相关的资料可以学习一下这方面内容的呢（本人菜鸡，跪求指教）

@李东岳 老师，公式中的p是坐标，然后是共顶点的两条边构成的向量叉乘吗？

如果是这样，叉乘不是两边的模长相乘，乘以两边夹角的正弦么，就相当于底乘高了，再乘上0.5就是面积了？（基于首句的假设）

$0.5\ast \vec{l_{AB}} \times \vec{l_{AC}} = 0.5 \ast l_{AB} \ast l_{AC} \ast \mathbf{sinA} = 0.5 \ast l_{AB} \ast H = S$

不知道是不是这样，求各位大佬轻喷

@李东岳
老师，请问对于定义在体心上的物理量，这里correct的边界场，是指边界上的一层网格吗

@wangfei9088
请问老师，有没有可能，是因为我们将p在pressureReference类下的对象命名为pressureReference后，再写pf时：

pressureReference pressureReference(pf, simple.dict());

编译器会认为这行第一个pressureReference就是一个对象了，不再识别为一个类的名称进而创建新的对象呢