OpenFOAM中重力的植入方式
考虑重力的NS方程可以写为：
ρ ∂ u ⃗ ∂ t + ρ ∇ ( u ⃗ u ⃗ ) = ∇ ( μ ∇ u ⃗ ) − ∇ P + ρ g ⃗ (1) \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} ρ∂t∂u+ρ∇(uu)=∇(μ∇u)−∇P+ρg(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中的重力的表达式为： G ⃗ = − ( g ⃗ ⋅ s ⃗ ) ∇ ρ (2) \vec G=-(\vec g \cdot \vec s)\nabla \rho\tag {2} G=−(g⋅s)∇ρ(2)其中 s ⃗ \vec s s是待求点的位置矢量，这显然与方程1使用的重力表达式相差很大，接下来我们就一步一步推导这个结果。
首先明确，在涉及到重力的求解器中，其求解的压力p并不是总压P，而是动压，即： P = p + ρ g ⃗ ⋅ s ⃗ (3) P=p+\rho \vec g \cdot \vec s\tag {3} P=p+ρg⋅s(3)之所以要这样处理是为了方便用户设置初场，大家可以设想最简单的溃坝算例，如果直接对总压设置初场，那我们要设置的则是一个梯形的压力分布，计算域的y坐标越高其净水压强就越大，而如果我们拆分成动压和静压的形式，就可以直接设置初场为0了。可能有小伙伴对 g ⃗ ⋅ s ⃗ \vec g \cdot \vec s g⋅s项不太理解，其实我们只要对其展开就很好说明问题了： g ⃗ ⋅ s ⃗ = ( 0 , g , 0 ) ⋅ ( x , y , z ) = g y (4) \vec g \cdot \vec s=(0,g,0)\cdot(x,y,z)=gy\tag {4} g⋅s=(0,g,0)⋅(x,y,z)=gy(4)在动量方程中，我们存在对总压的压力梯度项，将总压拆成静压和动压即可得到: ∇ P = ∇ ( p + ρ g ⃗ ⋅ s ⃗ ) = ∇ p + ∇ ( ρ g ⃗ ⋅ s ⃗ ) = ∇ p + ( g ⃗ ⋅ s ⃗ ) ∇ ρ + ρ ∇ ( g ⃗ ⋅ s ⃗ ) (5) \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} ∇P=∇(p+ρg⋅s)=∇p+∇(ρg⋅s)=∇p+(g⋅s)∇ρ+ρ∇(g⋅s)(5)其中关键点在于处理 ∇ ( g ⃗ ⋅ s ⃗ ) \nabla(\vec g \cdot \vec s) ∇(g⋅s)项,根据式4有: ∇ ( g ⃗ ⋅ s ⃗ ) = ∇ ( g y ) = ( ∂ ( g y ) ∂ x , ∂ ( g y ) ∂ y , ∂ ( g y ) ∂ z ) = ( 0 , g , 0 ) = g ⃗ (6) \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} ∇(g⋅s)=∇(gy)=(∂x∂(gy),∂y∂(gy),∂z∂(gy))=(0,g,0)=g(6)因此，式5最终写为： ∇ P = ∇ p + ( g ⃗ ⋅ s ⃗ ) ∇ ρ + ρ ∇ ( g ⃗ ⋅ s ⃗ ) = ∇ p + ( g ⃗ ⋅ s ⃗ ) ∇ ρ + ρ g ⃗ (7) \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} ∇P=∇p+(g⋅s)∇ρ+ρ∇(g⋅s)=∇p+(g⋅s)∇ρ+ρg(7)将式7带回方程1，即可得到OpenFOAM代码中使用的含重力的动量方程了：
ρ ∂ u ⃗ ∂ t + ρ ∇ ( u ⃗ u ⃗ ) = ∇ ( μ ∇ u ⃗ ) − ∇ p − ( g ⃗ ⋅ s ⃗ ) ∇ ρ (8) \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} ρ∂t∂u+ρ∇(uu)=∇(μ∇u)−∇p−(g⋅s)∇ρ(8)
