particleforce类内ErgunWenYuDrag的表达式感觉有错误
-
最近在研究DPMFoam求解器时，看到了ErgunWenYuDrag力的表达式源程序是：
src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/ErgunWenYuDrag/ErgunWenYuDragForce.C +97 template<class CloudType> Foam::forceSuSp Foam::ErgunWenYuDragForce<CloudType>::calcCoupled ( const typename CloudType::parcelType& p, const scalar dt, const scalar mass, const scalar Re, const scalar muc ） const { scalar alphac(alphac_[p.cell()]); if (alphac < 0.8) { return forceSuSp ( Zero, (mass/p.rho()) *(150.0*(1.0 - alphac)/alphac + 1.75*Re)*muc/(alphac*sqr(p.d()))//感觉和原始文献里面有很大的区别 ); } else { return forceSuSp ( Zero, (mass/p.rho()) *0.75*CdRe(alphac*Re)*muc*pow(alphac, -2.65)/(alphac*sqr(p.d()))//这个也是 ); } }公式显示貌似有些问题。
翻译一下公式是:
\begin{equation}
\frac{m_{d}}{\rho _{d}}\frac{\left ( \frac{150\left ( 1-\alpha _{c} \right )}{\alpha _{c}} +1.75Re\right)\mu _{c}}{\alpha _{c}d^{2}}\left( \alpha {c}< 0.8\right )
\end{equation}
原始文献里面的计算公式:
\begin{equation}
\frac{150\left(1-\alpha{c} \right )^{2}\mu {c}}{\alpha{c}d^{2}}+1.75\frac{\rho _{c}\alpha {d}\left | v{r} \right |}{d}
\end{equation}
两个差了很多。不知道有木有大神解答一下，纠结了一整天了
-
经过仔细比对，发现OF里面ErgunWenYuDrag模型的实现确实有些许问题，下面是我的详细分析。
首先，参照文献Ding and Gidaspow, 1990所述,欧拉-欧拉模型的连续相方程为：
\begin{equation}
\label{E1}
\frac{\partial \left( \alpha_{c}\rho_{c}\mathbf{U}_ {c}\right)}{\partial b}+\nabla \cdot \left( \alpha_{c}\rho_{c}\mathbf{U}_ {c}\mathbf{U}_ {c}\right )-\nabla \cdot\left( \alpha_{c}\rho_{c}\tau \right )=-\alpha_{c}\nabla p+\alpha_{c}\rho_{c}\mathbf{g}-\beta \left( \mathbf{U}_ {c}-\mathbf{U}_ {d}\right )
\end{equation}
当曳力模型为ErgunWenYuDrag时，$\beta$定义为：
\begin{equation}
\label{E2}
\beta=\begin{Bmatrix}
Re< 1000 & \frac{24}{\alpha_cRe} \left[1+0.15\left(\alpha_cRe\right)^{0.687} \right] \\
Re \geq 1000 & 0.44
\end{Bmatrix}
\end{equation}
上面公式\eqref{E2}中$C_D$为阻力系数，与雷诺数$Re$相关，阻力系数和雷诺数的计算式为：
\begin{equation}
\label{E3}
\beta=\begin{Bmatrix}
Re< 1000 &\frac{24}{\alpha_cRe}\left[1+0.15\left(\alpha_cRe\right)^{0.687} \right ] \\
Re\geq 1000& 0.44
\end{Bmatrix}
\end{equation}
\begin{equation}
\label{E4}
Re=\frac{\rho_c\left|\mathbf{U}_ {c}-\mathbf{U}_ {d}\right|d}{\mu_c}
\end{equation}
然后参照我研读DPMFoam代码以及东岳老师的多相流数学模型一文，总结出了DPMFoam中欧拉-拉格朗日模型中的连续相的方程为（这个方程暂时不知道对不对）：
\begin{equation}
\label{E5}
\frac{\partial \left( \alpha_{c}\rho_c\mathbf{U}_ {c}\right)}{\partial b}+\nabla \cdot \left( \alpha_{c}\rho_c\mathbf{U}_ {c}\mathbf{U}_ {c}\right )-\nabla \cdot\left( \alpha_{c}\rho_{c}\tau \right )=-\alpha_{c}\nabla p+\alpha_{c}\rho_{c}\mathbf{g}-\sum_{i=1}^{n_c}\frac{m_d}{\rho_d}Kd\left(\mathbf{\bar{U}}_ {c}-\mathbf{U}_ {di} \right )/V_c
\end{equation}
公式\eqref{E5}中$n_c$为单个流体网格中的颗粒数，$V_c$为网格体积。
OpenFOAM中把颗粒受力$\mathbf{F}$写为$\mathbf{F}=Sp\left( \mathbf{U}-\mathbf{U}_ d\right )+\mathbf{Su}$,然后通过particleforce类中的calcCoupled或calcNonCoupled函数来计算不同颗粒所受作用力的Sp和$\mathbf{Su}$，根据解读的代码，ErgunWenYuDrag模型力中的$\mathbf{Su}=\mathbf{0}$,Sp计算式为：
\begin{equation}
\label{E6}
Sp=\begin{Bmatrix}
\alpha_c< 0.8 &\frac{m_d}{\rho_d}\frac{\left(\frac{150\left(1-\alpha_c \right )}{\alpha_c}+1.75Re \right)\mu_c}{\alpha_cd^2} \\
\alpha_c\geq 0.8& \frac{m_d}{\rho_d}\frac{3}{4}\frac{CdRe\left(\alpha_cRe\right)\mu_c\alpha_{c}^{-2.65}}{\alpha_cd^2}
\end{Bmatrix}
\end{equation}
公式\eqref{E6}中CdRe函数为：
\begin{equation}
\label{E7}
\beta=\begin{Bmatrix}
Re< 1000 &24\left[1+0.15\left(Re\right)^{0.687} \right ] \\
Re\geq 1000& 0.44Re
\end{Bmatrix}
\end{equation}
把\eqref{E6}写成与\eqref{E2}相同的形式，得到：
\begin{equation}
\label{E8}
Sp=\begin{Bmatrix}
\alpha_c< 0.8 &\frac{m_d}{\rho_d}\left( \frac{150\left(1-\alpha_c \right )\mu_c}{\alpha_cd^2}+1.75\frac{\rho_c\left|\mathbf{U}_ {c}-\mathbf{U}_ {d}\right|}{\alpha_cd}\right )\\
\alpha_c\geq 0.8& \frac{m_d}{\rho_d}\frac{3}{4}C_D\frac{\rho_c\left|\mathbf{U}_ {c}-\mathbf{U}_ {d}\right|}{d}\alpha_{c}^{-2.65}
\end{Bmatrix}
\end{equation}
比较\eqref{E8}和\eqref{E2}发现，\eqref{E8}中当$\alpha_c<0.8$时，后面一项中的分母上多了一个$\alpha_c$,当$\alpha_c\geqslant 0.8$时，分子上少了一个$\alpha_c$，同时都少一个$1-\alpha_c$是因为要把网格内颗粒相加再除以网格体积，得到$1-\alpha_c$。
最后总结一下，第一次提出问题时我说两个模型公式相差很大有些表述不准确，事实上是有细微区别，现在暂时不知道一个$\alpha_c$的影响有多大，也有可能是我哪里推导错了，不过我检查了一遍，暂时没发现错误，现在发送出来讨论一下。
此外，最近一直研究DPMFoam，不知道上面写出的欧拉-朗格朗日方程对不对？还有，如果颗粒比网格大了，用单个网格内的颗粒相加除以网格体积是不是就没有意义了，也就是说模型方程是不是就不能再这么写了？ -
9月CFD算法编程课： http://dyfluid.com/class.html
需要帮助debug算例的看这个 https://cfd-china.com/topic/8018
-
@东岳 李老师，关于您这张图中的公式（16）我有几个疑问，如果说得不对望见谅。
1.对于连续相动量方程的拉力梯度项前是否需要乘以一个相分数主要决定于$\mathbf{F}$中的压力梯度力项是否提取出来了，如果提取出来了和前面的压力梯度项组合就会出现乘以相分数的情况，如果不提取出来，把颗粒的压力梯度力放在$\mathbf{F}$中，前面就不用再乘以相分数。不知道我这么解释是否正确？这个解释我是参考的文章Zhou et al. 2010
2.对于您式子（16）中最后一项为何$\mathbf{F}$的前面要乘以一个颗粒相分数，并且为何是除以颗粒的体积，按道理不应该是除以流场网格体积吗？
希望能得到李老师的解答。 -
。总体比较来看是有细微区别。