代码实现？

李东岳

在OpenFOAM里面替代文字
被省略了，如果你要植入的话：

volTensorField tauU("tauU", tau & U);
volVectorField unknown("unknown", fvc::div(tauU));

植入这一项打算对比收敛性或者准确性？

小龙

@李东岳东岳老师

volVectorField unknown("unknown", fvc::div(tauU));

这个是不是应该改成

volScalarField unknown("unknown", fvc::div(tauU));

这一项应该是个标量场？

李东岳

tauU是二阶张量，做div后是矢量

小龙

@李东岳这里的tau是μgrad(U)吗？（μ是动力粘度），还是可以在OpenFoam里调用的？

李东岳

tau是$\tau$，需要用Bousinessq假定来计算，参考 http://47.93.231.95/doku.php#boussinesq_approximation

Izumi

@李东岳请问tau如下定义可以吗？

volTensorField tau("tau", muEff*(dev2(gradU.T()) + gradU));

因为
\begin{equation}
dev2((\nabla u)^{T})=(\nabla u)^{T} - \frac{2}{3}tr((\nabla u)^{T})I = (\nabla u)^{T} - \frac{2}{3}(\nabla . u) I
\end{equation}
所以
\begin{equation}
\tau = \mu (dev2((\nabla u)^T) + \nabla u)
\end{equation}

小龙

@izumi tau应该是这样定义的

程迪

@izumi
不建议这么实现，参考rhoCentralFoam的实现，应该是把转置项显式离散，非转置项按laplacian隐式离散。

因为OF本身的限制，转置项会有各个分量间的耦合，所以必须显式离散，非转置项不用。

隐式的部分越多，时间步限制通常也越弱一些，所以OF一般实践上喜欢把粘性项拆拆拆之后尽可能多地采用隐式离散。实在不能隐式离散的源项、交叉项才用显式离散。

参考：https://chengdi123000.github.io/2018/01/05/rhoCentralFoam解析/#粘性应力项

Izumi

@程迪
谢谢！看rhoCentralFoam时没看懂，看了你的链接清楚一些了。
我主要是想在能量方程中添加$\nabla \cdot(\tau \cdot u)$，然后计算耗散函数$\tau:(\frac12(\nabla u + (\nabla u)^T))$。
根据CFD中的能量方程，$\nabla \cdot(\tau \cdot u)$好像没法展开成可以隐式离散的项。
看rhoCentralFoam中的能量方程，除时间项外，用的都是显式离散。所以，对于能量方程来说，是不是不必拆项并使用隐式格式？
至于动量方程中的$\nabla \cdot \tau$，对于一般的求解器，应该是包含在UEqn.H里的turbulence->divDevRhoReff(U)吧，不用我们自己定义。

rhoCentralFoam中的能量方程：

surfaceScalarField sigmaDotU
        (
            "sigmaDotU",
            (
                fvc::interpolate(muEff)*mesh.magSf()*fvc::snGrad(U)
              + fvc::dotInterpolate(mesh.Sf(), tauMC)
            )
          & (a_pos*U_pos + a_neg*U_neg)
        );

        solve
        (
            fvm::ddt(rhoE)
          + fvc::div(phiEp)
          - fvc::div(sigmaDotU)
        );

程迪

@izumi
我的理解是，拆项的话你得把phiEp和sigmaDotU表示成rhoE的函数才行，这似乎有点儿困难吧。。

Izumi

@程迪
rhoE应该是$\frac{\partial \rho e}{\partial t} + \frac{\partial \rho K}{\partial t}$
phiEp应该是$\nabla \cdot (\rho \mathbf{U} e)+\nabla \cdot (\rho \mathbf{U} K)+\nabla\cdot(p\mathbf{U})$
sigmaDotU应该是$\nabla \cdot(\tau \cdot \mathbf{U})$
感觉拆不了。

另外，请问对定常的求解器也有时间步限制吗？如果没有的话，是不是定常的求解器中，比如buoyanySimpleFoam，就不用管隐式还是显式了？

CFD中文网

代码实现？