Darcy-Brinkman-Forchheimer方程的无量纲化

根据文献[chen2008],使用Darcy-Brinkman-Forchheimer扩展模型描述孔隙流动的方程如下:

不可压方程:\(\nabla\cdot\vec{u}=0\)

动量方程:\(\rho\frac{\partial \vec{u}}{\partial t}+\nabla\cdot\left(\frac{\rho\vec{u}\vec{u}}{\epsilon}\right)=-\nabla(\epsilon p^*)+\mu\nabla^2\vec{u}-\frac{\mu\epsilon}{K}\vec{u}-\frac{\rho\epsilon C_F|\vec{u}|}{\sqrt{K}}\vec{u}\)

在这里的\(\vec{u}\)是Darcy velocity,或local average velocity,或superficial velocity。其意义是,单位截面积孔隙介质中流出的流体通量。这里的“单位截面积”是固体和液体一起算,所以这个Darcy velocity会比孔隙中流体质点的“真实速度”小一些。二者满足关系:Darcy velocity=真实速度*孔隙度(孔隙度一般记作\(\epsilon\))。而p*这个星号的含义是它是intrinsic也就是流体质点上的真实压力。另外,Forchheimer系数由\(C_F=1.75/\sqrt{150\epsilon^3}\)定义,它和孔隙度一样是一个无量纲数。

我们定义如下特征量:特征速度U(定义为远处来流速度),特征尺度D(依问题不同而异,比如说管流的D可取管径,圆柱绕流的D可取圆柱直径等等)。

在这样的流动中,比较关键的无量纲参数有雷诺数\(Re=U\rho D/\mu\)和达西数(Darcy number)\(Da=K/D^2\)。这里我们认为流体密度\(\rho\)是一个常量。

定义无量纲速度\(\vec{u}’=\vec{u}/U\),无量纲压力\(p’=p^*/(\rho U^2)\).

不可压方程转换为无量纲形式之后就是乘了一个常数,形式不变:

\(\nabla\cdot\vec{u}’=0\)

动量方程比较麻烦。首先无脑转换:

\(U\rho\frac{\partial \vec{u}’}{\partial t}+U^2\nabla\cdot\left(\frac{\rho\vec{u}’\vec{u}’}{\epsilon}\right)=-(\rho U^2)\nabla(\epsilon p’)+\mu U\nabla^2\vec{u}’-\frac{U\mu\epsilon}{K}\vec{u}’-\frac{\rho U^2\epsilon C_F|\vec{u}’|}{\sqrt{K}}\vec{u}’\)

然后把每一项除以\(\rho U^2\) :

\(\frac{1}{U}\frac{\partial \vec{u}’}{\partial t}+\nabla\cdot\left(\frac{\vec{u}’\vec{u}’}{\epsilon}\right)=-\nabla(\epsilon p’)+\frac{\mu}{\rho U}\nabla^2\vec{u}’-\frac{\mu\epsilon}{\rho UK}\vec{u}’-\frac{\epsilon C_F|\vec{u}’|}{\sqrt{K}}\vec{u}’\)

为了凑雷诺数和达西数,把每一项乘以D:

\(\frac{D}{U}\frac{\partial \vec{u}’}{\partial t}+D\nabla\cdot\left(\frac{\vec{u}’\vec{u}’}{\epsilon}\right)=-D\nabla(\epsilon p’)+\frac{D^2\mu}{\rho UD}\nabla^2\vec{u}’-\frac{\mu\epsilon D^2}{\rho UDK}\vec{u}’-\frac{D\epsilon C_F|\vec{u}’|}{\sqrt{K}}\vec{u}’\)

也就是:

\(\frac{D}{U}\frac{\partial \vec{u}’}{\partial t}+D\nabla\cdot\left(\frac{\vec{u}’\vec{u}’}{\epsilon}\right)=-D\nabla(\epsilon p’)+\frac{D^2}{Re}\nabla^2\vec{u}’-\frac{1}{Re}\frac{\epsilon}{Da}\vec{u}’-\frac{\epsilon C_F|\vec{u}’|}{\sqrt{Da}}\vec{u}’\)

至此无量纲化还未完成,因为时间、空间、U、D都是带量纲的。为了把它们干掉,按照[nithiarasu2002]文中方法,把坐标和时间也无量纲化为:\(x’=x/D\)(y,z无量纲化方法相同,略去),\(t’=tU/D\)。这里实际上是认为特征尺度/特征速度=特征时间(=D/U)。

如此一来,空间偏导和时间偏导也就相应地无量纲化了。特别地,无量纲的nabla算子满足\(\nabla’=D\nabla\).如此一来我们就得到了最终的无量纲动量方程:

\(\frac{\partial \vec{u}’}{\partial t’}+\nabla’\cdot\left(\frac{\vec{u}’\vec{u}’}{\epsilon}\right)=-\nabla'(\epsilon p’)+\frac{1}{Re}\nabla’^2\vec{u}’-\frac{1}{Re}\frac{\epsilon}{Da}\vec{u}’-\frac{\epsilon C_F|\vec{u}’|}{\sqrt{Da}}\vec{u}’\)

至于连续方程,在无量纲nabla算子下也是差个常数,仍然为:

\(\nabla’\cdot\vec{u}’=0\)

无量纲坐标的特点是,在该坐标系下,那个特征尺度为1.不是网格雷诺数。

最后是 Courant 数。这是一个无量纲参数。其定义为\(C=U\Delta t/\Delta x\),也就是\(C=\Delta t’/\Delta x’\).

也有一种姿势是\(C=u\Delta t/\Delta x\),u是网格中最大速度什么的,那么就是\(C=u’\Delta t’/\Delta x’\)

参考文献

[chen2008] Chen X, Yu P, Winoto S H, et al. Numerical analysis for the flow past a porous square cylinder based on the stress-jump interfacial-conditions[J]. International Journal of Numerical Methods for Heat & Fluid Flow, 2008, 18(5): 635-655.

[nithiarasu2002] Nithiarasu P, Seetharamu K N, Sundararajan T. Finite element modelling of flow, heat and mass transfer in fluid saturated porous media[J]. Archives of Computational Methods in Engineering, 2002, 9(1): 3.

发表评论

电子邮件地址不会被公开。 必填项已用*标注