CloudBird's Idea Space

考虑重力的NS方程可以写为：

$\rho \frac{\partial \vec u}{\partial t}+\rho \nabla(\vec u \vec u)=\nabla(\mu\nabla \vec u)-\nabla P +\rho \vec g \tag {1}$

在openFOAM中，许多求解器中都耦合了重力，以interFoam为例，其中的重力定义处的代码块为：

  fvc::reconstruct
            (
                (
                    mixture.surfaceTensionForce()
                  - ghf*fvc::snGrad(rho)
                  - fvc::snGrad(p_rgh)
                ) * mesh.magSf()
            )

将其中的重力抽离出来，其表达式为fvc::reconstruct( (ghf*fvc::snGrad(rho) ）* mesh.magSf() )。通过溯源代码，我们可以进一步得到其中ghf的表达式：

new surfaceScalarField
(
    "ghf",
    (gFluid[i] & fluidRegions[i].Cf()) - ghRef
)

其中gFluid[i]返回了通过字典(constant/g)读取得到的重力加速度项，fluidRegions[i].Cf()返回了当前网格表面的面心矢量,ghRef是给定的相对重力加速度，一般为0。综上所述，openFOAM中的重力的表达式为：

$\vec G=-(\vec g \cdot \vec s)\nabla \rho\tag {2}$

其中$\vec s$是待求点的位置矢量，这显然与方程1使用的重力表达式相差很大，接下来我们就一步一步推导这个结果。
首先明确，在涉及到重力的求解器中，其求解的压力p并不是总压P，而是动压，即：

$P=p+\rho \vec g \cdot \vec s\tag {3}$

之所以要这样处理是为了方便用户设置初场，大家可以设想最简单的溃坝算例，如果直接对总压设置初场，那我们要设置的则是一个梯形的压力分布，计算域的y坐标越高其净水压强就越大，而如果我们拆分成动压和静压的形式，就可以直接设置初场为0了。可能有小伙伴对$\vec g \cdot \vec s$项不太理解，其实我们只要对其展开就很好说明问题了：

$\vec g \cdot \vec s=(0,g,0)\cdot(x,y,z)=gy\tag {4}$

在动量方程中，我们存在对总压的压力梯度项，将总压拆成静压和动压即可得到:

$\nabla P=\nabla (p+\rho \vec g \cdot \vec s)=\nabla p+\nabla(\rho \vec g \cdot \vec s)=\nabla p+ ( \vec g \cdot \vec s)\nabla\rho+\rho\nabla( \vec g \cdot \vec s)\tag {5}$

其中关键点在于处理$\nabla(\vec g \cdot \vec s)$项,根据式4有:

$\nabla(\vec g \cdot \vec s)=\nabla (gy)=(\frac{\partial (gy)}{\partial x},\frac{\partial (gy)}{\partial y},\frac{\partial (gy)}{\partial z})=(0,g,0)=\vec g\tag {6}$

因此，式5最终写为：

$\nabla P=\nabla p+ ( \vec g \cdot \vec s)\nabla\rho+\rho\nabla( \vec g \cdot \vec s)=\nabla p+ ( \vec g \cdot \vec s)\nabla\rho+\rho \vec g\tag {7}$

将式7带回方程1，即可得到OpenFOAM代码中使用的含重力的动量方程了：

$\rho \frac{\partial \vec u}{\partial t}+\rho \nabla(\vec u \vec u)=\nabla(\mu\nabla \vec u)-\nabla p - ( \vec g \cdot \vec s)\nabla\rho \tag {8}$