Darcy-Forchheimer model

OpenFOAM のDarcy-Forchheimer model で設定する D,F を整理。ついでに透水係数との関係も整理。

まずは、公式wikiより。なぜか太字が出ないのは御愛嬌。
https://openfoamwiki.net/index.php/DarcyForchheimer

The Darcy Forchheimer acts in the momentum equation as a sink term Sm. Considering the momentum equation, it follows:
\begin{align*}\frac{\partial\rho U}{\partial t}+\mathrm{\nabla}\left(\rho UU\right)=\mathrm{\nabla\sigma}+S_m
\end{align*}
The main important term is the source term Sm which is given as:
\begin{align*}S_m=-\left(\mu D+\frac{1}{2}\rho tr\left(U\cdot I\right)F\right)U\end{align*}
As the source term Sm can be expressed as a pressure gradient it follows Sm=∇p. Therefore, we can write:
\begin{align*}\nabla p=-\left(\mu D+\frac{1}{2}\rho tr\left(U\cdot I\right)F\right)U\end{align*}
While switching to the Cartesian coordinate system, we can achieve something similar to:
\begin{align*}\nabla p=-\left(\mu D_i-\frac{1}{2}\rho F_i\left|u_{kk}\right|\right)u_i\end{align*}
If we know the velocity dependent pressure we can calculate a polynomial function such as:
\begin{align*}-\mathrm{\nabla p}=Au+Bu^2\end{align*}
\begin{align*}A=\mu D\end{align*}
\begin{align*}B=\frac{1}{2}\rho F\end{align*}

ｐ：圧力損失P[N/m2], u：流速［m/s］, μ：動粘性係数[N/m2・s]=[kg*m/s2/m2*s]=[kg/m/s], ρ：流体密度[kg/m3]

CfdOF では、D、F を constant/fvOptions で指定しています。
F＝0 でダルシー則。単位は以下。

∇P [N/m3] ＝Au[N/m4*s*m/s]＋Bu２[N/m5*s2*m/s*m/s]
係数A[N/m4*s]=[kg*m/s2/m4*s]=[kg/m3/s]
係数B[N/m5*s2]=[ kg*m/s2/m5*s2]=[kg/m4]
D=A/μ[ N/m4*s*/N*m2/s]=[1/m2]
F=2B/ρ [kg/m4/kg*m3] =[1/m]

吉岡ほか（2010）「高透水性多孔質体中の非ダルシー流れに関する考察」では、ρg [kg/m3*m/s2=N/m3] で除し無次元化した式が紹介されています。係数 a、b に関して Ergun 等の近似式を検証しています。こちらの式の方が土木では‘なじみ’があります。
\begin{align*}-\mathrm{\nabla h}=au+bu^2\end{align*}
\begin{align*}a=\frac{150v\left(1-\emptyset\right)^2}{g\emptyset^3d^2}\end{align*}
\begin{align*}b=\frac{1.75\left(1-\emptyset\right)}{g\emptyset^3d}\end{align*}
h：水頭［m］, a：係数［s/m］, b：係数［s２/m２］, ｄ：粒径[m], φ：間隙率［－］, ν：動粘性係数［m２/s］, g：重力加速度［m/s２］

u=-1/a∇hより、透水係数 k＝1/a[m/s]

換算は以下の通り。
\begin{align*}D=A/\mu＝aρg/μ=\frac{150v\left(1-\emptyset\right)^2}{g\emptyset^3d^2}ρg/μ\end{align*} [s/m*N/m3/N*m2/s]=[ 1/m/m3*m2]=[1/m2]
\begin{align*}F=2B/\rho=2b\rho g/\rho=2\frac{1.75\left(1-\emptyset\right)}{g\emptyset^3d}\rho g/\rho\end{align*}    [1/m]
あるいは、
\begin{align*}D=\frac{\rho g}{k\mu}\end{align*}

***********************
20220611追記
見通しを良くしてみました。
https://phreeqc.blogspot.com/2022/06/blog-post_11.html