求助:孔隙介质中流体被加热上浮过程的压力方程求解问题



  • 孔隙介质中的流体被加热上浮过程可以与自然界中的温泉对应起来,温泉的流体循环过程其实就是冷的自然水下渗进入岩石孔隙中受地热加热之后密度变小而上涌,最后喷出地表的过程。求解过程应该与buoyantPimpleFoam类似,只是动量方程和压力方程不同。一般在地球科学领域内,这类温泉的流体循环流动速度被认为是符合达西定律(Darcy’s Law)的。主要控制方程如下:

    \begin{equation}
    \vec u = - \frac{k}{\mu }\left( {\nabla P - \rho \vec g} \right)
    \end{equation}

    \begin{equation}
    \phi \rho \beta \frac{{\partial P}}{{\partial t}} = \nabla \cdot \left( {\frac{{k\rho }}{\mu }\left( {\nabla P - \rho \vec g} \right)} \right) - \phi \rho \alpha \frac{{\partial T}}{{\partial t}}
    \end{equation}

    \begin{equation}
    \left( {\phi {\rho_f}{c_{pf}} + (1 - \phi ){\rho_r}{c_{pr}}} \right)\frac{{\partial T}}{{\partial t}} = \nabla \cdot {k_m}\nabla T - \vec u\rho {c_{pf}}\nabla T
    \end{equation}

    方程中需要求解的未知量为温度T,压力P和速度v。密度暂且近似为温度的线性函数,其他的量均为常数(是温度和压力的函数,由流体也就是水的热力学参数公式给出)。
    压力的初始条件为静水压。

    方程1,3我都会求解,方程1其实是一个显示求解梯度的问题,求出达西速度;方程3其实就是一个标量运输问题,涉及到diffusion和advection,求出温度场。并且做了benchmark都没问题。可是方程2里面即有压力瞬态项,又有温度瞬态项,不知道如何耦合求解。我现在正在同步学习OpenFoam和有限体积法的理论,参照buoyantPimpleFoam求解器的思路,但是还是不确定如何下手。希望论坛里的各位高人看到之后能指点一二,感激不尽!


  • CORE 网格教授 OpenFOAM教授 管理员

    $\phi$ $\beta$ 都是常量?公式11好像能化简成:?
    \begin{equation}
    \phi \rho \beta \frac{{\partial P}}{{\partial t}} = \nabla \cdot (\rho \vec{u}) - \phi \rho \alpha \frac{{\partial T}}{{\partial t}}
    \end{equation}

    你的三个方程是整套的方程么?不需要连续性方程啥的么



  • 首先非常感谢东岳老师回复答疑!估计李老师您对海底黑烟囱也有了解的,这个问题其实就是模拟seafloor hydrothermal circulation(海水通过洋壳的孔隙和裂隙下渗至深部被岩浆房等热源加热,密度变小,在上浮作用力下最终喷出至海底)。

    (一)$\phi$表示孔隙度是常数;$\alpha$和$\beta$分别表示热膨胀系数和压缩系数,虽然是温度和压力的函数,可以用相应的理论和公式计算出来,这里可以简化认为是个常数。

    (二)确实还有一个连续性方程的(如下公式4),这几个方程和这个nature文章的method部分的公式3,4,5,6对应的。

    \begin{equation}
    \frac{{\partial \phi \rho }}{{\partial t}} = - \nabla \cdot \left( {\vec u\rho } \right)
    \end{equation}

    方程(2)也就是压力方程,是将达西速度带入这个连续性方程,然后对密度$\rho$求了对温度$T$和压力$P$的全导数得来的。



  • @东岳
    我查了一些资料和东岳流体网站上的一些solver的讲解,这个问题是不是得用这个PIMPLE算法求解,但是不确定努力的方向是否正确


  • CORE 网格教授 OpenFOAM教授 管理员

    从Nature的文藏来看:

    1. 求解压力方程:
      \begin{equation}
      \nabla \cdot \left( {\frac{{k\rho }}{\mu }\left( {\nabla P - \rho \vec g} \right)} \right) = \phi \frac{\partial \rho}{\partial t}
      \end{equation}

    2. 求解速度方程:
      \begin{equation}
      \vec u = - \frac{k}{\mu }\left( {\nabla P - \rho \vec g} \right)
      \end{equation}

    3. 求解温度方程

    明天再详细看看,
    你在研究Nature么 :wocao: 有前途..



  • @东岳
    太感谢了!我博士论文就是做这个流体循环并且耦合矿物反应的模拟,目前我们用的是有限元方法, 我从您的网站和论坛上了解到有限体积法和openfoam,感觉这个更有意思。所以我想用openfoam尝试做后面的模拟任务。

    您的意思是说直接求解方程(8)表示的这个压力方程对吧?可否推荐一个相似的solver的求解过程或者其他相似的资料给我参考一下。看了buoyantPimpleFoam的求解跟我这个问题最相近,但是它的求解方程式N-S方程,我这个问题应该是NS方程的一个简化形式。我尝试根据您的指导把代码写出来,再请教您


  • CORE 网格教授 OpenFOAM教授 管理员

    有几个数值问题:

    • 密度变化么?可压缩?
    • 如果恒定密度,压力方程非常好求,求完压力解速度就行,可以参考laplacianFoam,里面求解了个压力播送方程,跟你的压力方程差不多


  • @东岳
    我这里的流体是水(water),不可压缩。密度是变化的,$\rho=\rho(T,P)$,可以通过IAPWS方程求解得到,简化认为是温度的线性函数。


  • CORE 网格教授 OpenFOAM教授 管理员

    如果不可压

    1. 求解压力方程:
      \begin{equation}
      \nabla \cdot \left( {\frac{{k\rho }}{\mu } {\nabla P } } \right) = 0
      \end{equation}

    2. 求解速度方程:
      \begin{equation}
      \vec u = - \frac{k}{\mu } {\nabla P }
      \end{equation}

    3. 求解温度方程

    Pretty simple :mianmo:


 

Forest
Mountains