根据文献[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.