代码实现？

[小龙]

图中的式子，我想加在能量方程里，代码该怎么实现？（我还是一个编程菜鸟，不太会编程，希望有人可以帮我一下）

[小龙]

OpenFOAM中有考虑粘性耗散项的求解器吗？或者考虑粘耗散的算例？

[李东岳]

你这个不是粘性耗散吧？这个方程可以实现，首先对速度求梯度然后提取左下角和右下角分量处理（公式6：http://dyfluid.com/tensor.html ）

volTensorField gradU("gradU", fvc::grad(U));
scalarField parU1 = gradU.xz();
scalarField parU2 = gradU.zz();

类似这种代码。不过感觉不是很符合物理。这个方程怎么推导出来的？

[小龙]

@李东岳 这个公式里的最后一项应该是粘性耗散热吧？如果用这个编程 ，该怎么实现呢？是不是和这一项是一致的？ 我研究的模型中，速度在x和y方向的变化很小，忽略以后就是那样的了。（东岳老师给的链接页面不存在了）

[李东岳]

在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，就不用管隐式还是显式了？