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

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

不可压方程:[latex]\nabla\cdot\vec{u}=0[/latex]

动量方程:[latex]\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}[/latex]

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

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

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

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

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

[latex]\nabla\cdot\vec{u}’=0[/latex]

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

[latex]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}'[/latex]

然后把每一项除以[latex]\rho U^2[/latex] :

[latex]\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}'[/latex]

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

[latex]\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}'[/latex]

也就是:

[latex]\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}'[/latex]

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

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

[latex]\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}'[/latex]

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

[latex]\nabla’\cdot\vec{u}’=0[/latex]

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

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

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

参考文献

[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.

发表评论

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