particleforce类内ErgunWenYuDrag的表达式感觉有错误
-
经过仔细比对，发现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，不知道上面写出的欧拉-朗格朗日方程对不对？还有，如果颗粒比网格大了，用单个网格内的颗粒相加除以网格体积是不是就没有意义了，也就是说模型方程是不是就不能再这么写了？ -
-
@东岳 李老师，关于您这张图中的公式（16）我有几个疑问，如果说得不对望见谅。
1.对于连续相动量方程的拉力梯度项前是否需要乘以一个相分数主要决定于$\mathbf{F}$中的压力梯度力项是否提取出来了，如果提取出来了和前面的压力梯度项组合就会出现乘以相分数的情况，如果不提取出来，把颗粒的压力梯度力放在$\mathbf{F}$中，前面就不用再乘以相分数。不知道我这么解释是否正确？这个解释我是参考的文章Zhou et al. 2010
2.对于您式子（16）中最后一项为何$\mathbf{F}$的前面要乘以一个颗粒相分数，并且为何是除以颗粒的体积，按道理不应该是除以流场网格体积吗？
希望能得到李老师的解答。 -
-
@上级
您好，请问您上面的公式 F=Sp(U−Ud)+Su用来计算颗粒受力的，这里面的U和Ud分别表示的是什么东西？和KinematicParcel中的Uc_以及U_是对应关系吗const parcelType& p = static_cast<const parcelType&>(*this); const forceSuSp Fcp = forces.calcCoupled(p, dt, mass, Re, mu); const forceSuSp Fncp = forces.calcNonCoupled(p, dt, mass, Re, mu); const forceSuSp Feff = Fcp + Fncp; const scalar massEff = forces.massEff(p, mass); // Update velocity - treat as 3-D const vector abp = (Feff.Sp()*Uc_ + (Feff.Su() + Su))/massEff; const scalar bp = Feff.Sp()/massEff; Spu = dt*Feff.Sp(); IntegrationScheme<vector>::integrationResult Ures = td.cloud().UIntegrator().integrate(U_, dt, abp, bp); vector Unew = Ures.value(); // note: Feff.Sp() and Fc.Sp() must be the same dUTrans += dt*(Feff.Sp()*(Ures.average() - Uc_) - Fcp.Su()); // Apply correction to velocity and dUTrans for reduced-D cases const polyMesh& mesh = td.cloud().pMesh(); meshTools::constrainDirection(mesh, mesh.solutionD(), Unew); meshTools::constrainDirection(mesh, mesh.solutionD(), dUTrans); return Unew; }如果我想要计算颗粒所受的ErgunWenYuDrag这个力的时候是不是就等于
(mass/p.rho() *(150.0*(1.0 - alphac)/alphac +1.75*Re)*muc/(alphac*sqr(p.d()))*（U-Ud）//这样就OK了吗？
。总体比较来看是有细微区别。