张量二阶偏导数 偏(A_j B_i)/偏(x_i x_j) 如何编程

[学流体的小明]

想计算一个源项：
$$
\frac{\partial^2 (\overline{\rho} \overline{U}_j u_i^\prime) }{\partial x_i \partial x_j}
$$
其中上横线表示时间平均量，上标撇表示脉动量。
我把rhoMean、UMean和UPrime都计算好了，它们的类型分别是
volScalarField rhoMean; volVectorField UMean; volScalarField UPrime;
展开后是
$$
\begin{multline}
[{\frac{{{\partial ^2}}}{{\partial x_1^2}}(\bar \rho {{\bar U}_1}u_1^\prime ) + \frac{{{\partial ^2}}}{{\partial {x_1}\partial {x_2}}}(\bar \rho {{\bar U}_2}u_1^\prime ) + \frac{{{\partial ^2}}}{{\partial {x_1}\partial {x_3}}}(\bar \rho {{\bar U}_3}u_1^\prime ) + \frac{{{\partial ^2}}}{{\partial {x_2}\partial {x_1}}}(\bar \rho {{\bar U}_1}u_2^\prime ) +
\\
\frac{{{\partial ^2}}}{{\partial x_2^2}}(\bar \rho {{\bar U}_2}u_2^\prime ) + \frac{{{\partial ^2}}}{{\partial {x_2}\partial {x_3}}}(\bar \rho {{\bar U}_3}u_2^\prime ) + \frac{{{\partial ^2}}}{{\partial {x_3}\partial {x_1}}}(\bar \rho {{\bar U}_1}u_3^\prime ) + \frac{{{\partial ^2}}}{{\partial {x_3}\partial {x_2}}}(\bar \rho {{\bar U}_2}u_3^\prime ) + \frac{{{\partial ^2}}}{{\partial x_3^2}}(\bar \rho {{\bar U}_3}u_3^\prime )}]
\end{multline}
$$
然后该如何写这个偏导数呢？和AI战斗了好久它还是没懂，它一直想用fvc::laplacian()，但是fvc::laplacian()的含义是$\frac{\partial ^2}{\partial x_1^2}+\frac{\partial ^2}{\partial x_2^2}+\frac{\partial ^2}{\partial x_3^2}$。

[coolhhh]

@学流体的小明
问了AI，OpenFOAM需要手动做两次梯度计算

#include "fvCFD.H"      // OpenFOAM 的核心头文件，包含 mesh、field、fvc/fvm 工具
#include "fvcGrad.H"    // 显式声明 fvc::grad

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"
    #include "createMesh.H"
    // 创建速度场 U
    volVectorField U
    (
        IOobject
        (
            "U",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

    // 把 u1 分离出来
    volScalarField u1 = U.component(vector::X);  // 或者 U.x()

    // 1. 计算 u1 的一阶梯度 ∇u1
    volVectorField gradU1 = fvc::grad(u1);
    //    gradU1.x() = ∂u1/∂x1
    //    gradU1.y() = ∂u1/∂x2
    //    gradU1.z() = ∂u1/∂x3

    // 2. 对 ∂u1/∂x1 这个标量场，在 x2（Y）方向再求一次梯度
    //    fvc::grad(gradU1.x()) 是一个 volVectorField，
    //    其中 .y() 即 ∂(∂u1/∂x1)/∂x2 = ∂²u1/∂x1∂x2
    volScalarField d2u1dx1dx2 = fvc::grad( gradU1.x() ).y();

    // 写入结果
    d2u1dx1dx2.write();

    return 0;
}

[学流体的小明]

@coolhhh 感谢，本来也是想着OpenFOAM有没有简单的代码全部搞定，最后还是要手动一个一个算😂

[李东岳]

我推了一下，好像是这个$\nabla\cdot(\nabla\cdot(\bfU\bfU'))$，你验证下？