地形网格cfmesh一键生成真的很方便,我前段时间做大尺度地形界面泥石流模拟,用这个真的省了不少事儿。
上级
帖子
-
-
据我所知目前openfoam都是pointparticle,没有resolve方法,dev版本里有重叠网格法,适合单个或者几个颗粒的计算。resolve法的话cfdem可以实现,当然找找的话也有其他的开源求解器,不过不是基于openfoam的。
-
看这个散点图,每0.5s一个,是不是0.5s保存一组数据,到云图里面测量的?
-
openfoam自带的颗粒处理方法是点源(point model),无法解析出颗粒的边界。
-
@杨英狄 目前是已经实现了,文章是前面那个题为"An optimized Eulerial-Lagrangian..."的文章。具体要请教什么呢?
-
@杨英狄 这个在肖恒老师的两篇文章详细讨论过吧,算是一种经典的处理粗颗粒的办法:
1:Diffusion-based coarse graining in hybrid continuum–discrete solvers: Theoretical formulation and a priori tests
2:Diffusion-based coarse graining in hybrid continuum–discrete solvers: Applications in CFD–DEM
然后也有很多这个方法的变种(基于OpenFOAM或者CFDEM)
1:An optimized Eulerian–Lagrangian method for two-phase flow with coarse particles: Implementation in open-source field operation and manipulation, verification, and validation
2: A semi-resolved CFD-DEM approach for particulate flows with kernel based approximation and Hilbert curve based searching strategy -
@李东岳 感觉和相变类似,但是我理解来看相变是相间的物质交换,这里是边界处的输入,单纯的是泥石流相的物质增加。我在想能不能实现一个边界条件,先判断坡面上哪里有泥石流(边界面相分数判断),然后把边界面有相分数的的网格面提取出来,对这些网格面改变$\alpha$的值。
-
目前正在用OpenFOAM做坡面泥石流相关的数值模拟,主要是基于InterMixingFoam求解器(interFoam的衍生),将泥石流看作两相流体的混合物(颗粒与水),然后气相用来捕捉泥石流混合物运动过程的自由面。这是大致的示意图。
前期修改求解器实现了泥石流内部两相的相互作用(考虑速度滑移,interFoam本身未考虑速度滑移),预测了现场泥石流的运动过程,这是计算的几何域和计算的泥石流运动的速度云图。
图二 计算几何域
图三 泥石流运动速度云图现在想要进一步延申。之前是直接将底部边界设置为noSlip边界,但是泥石流沿着坡面向下运动时会携带破面的物质向下走,导致泥石流相的体积会越来越大,文献中叫“坡面侵蚀”,想在数值上实现这个过程,我理解来看应该是下壁面上的一种边界条件,但是因为没有接触过类似的问题,现在不知道该用什么边界条件。
所以把问题发出来想咨询一下有没有研究类似问题或者类似过程的大佬可以给指个方向
-
@oitocfd 之前写过一个小函数,具体是在dpmfoam里用到的,是寻找一个颗粒所在的cell label以及这个cell周围的一圈cell label,不知道能不能帮到你。
template<class Type> const Foam::labelList& Foam::interpolationCellAround<Type>::cellpointCells ( const label celli, DynamicList<label>& storage ) const { const labelList& cPoints = this->psi_.mesh().cellPoints()[celli]; storage.clear(); forAll(cPoints, i) { const label pointi = cPoints[i]; const labelList& pointcellsi = this->psi_.mesh().pointCells()[pointi]; forAll(pointcellsi, j) { storage.append(pointcellsi[j]); } } // Filter duplicates if (storage.size() > 1) { sort(storage); label n = 1; for (label i = 1; i < storage.size(); i++) { if (storage[i-1] != storage[i]) { storage[n++] = storage[i]; } } // truncate addressed list storage.setSize(n); } return storage; } template<class Type> const Foam::labelList& Foam::interpolationCellAround<Type>::cellpointCells(const label celli) const { return cellpointCells(celli, labels_); }
interpolationCellAround是自定义的类,主要是函数的实现。
-
感谢两位大佬的帮助
-
计算滑块入水问题,初始想在计算域框选出一个三角形区域,不知道能否实现?
![Y}}FNTDQL8G3J4T)7OM9M2.png -
不知道这篇论文能否给点帮助(同学的论文):
"叶秉晟. 针对非定常空化流动的多相可压缩求解器开发及应用[D]. 中国科学院大学,2019." -
@李东岳 谢谢李老师,我再详细看一下推导过程
-
@mechanicsdog 谢谢您,万分感谢
-
@李东岳 读了东岳老师的这篇文章,不知道公式(22)的量纲是否有点问题,按照公式(21),弛豫时间的$\tau$量纲应该是$[M^{-1}L^{3}T^{1}]$,那么公式(22)右边第二项的两处量纲都是对不上的,大括号内$\mathbf{g}$那一串不是速度量纲,右边的$\tau$应该是少一个-1次方。另外,请问李老师,您文章onewaycoupledtestcase里面计算颗粒分离时,用的弛豫时间$10d^{2/3}$,请问一下这部分该怎么理解,2/3次方表示的是什么呢?另外,弛豫时间一般理解,量纲不是时间量纲么?谢谢李老师。
贴出您文章的公式(21)和(22)
-
@mechanicsdog 您好,我想再请教您一个问题,我按照您的方法实现颗粒上显示速度了,现在我在算一个带颗粒级配的流化床问题,颗粒直径是高斯分布,请问后处理的时候您知道如何显示颗粒的大小么,现在按照高斯点好像只能显示一样大小的颗粒,如果想这么显示您知道如何操作么?谢谢您
-
我也遇到过相同的问题,平均压力没问题,但是瞬时波动很大,暂时没有找出原因所在,难道是算法的问题么
-
想问作者一个问题,请问是如何在paraview里面把颗粒上色,显示不同的速度的大小的呢。
此外,是否压力本来就该取平均值,在controlDict里面添加平均库函数来获得平均场,在这里贴上我的用过的例子:functions { fieldAverage { type fieldAverage; libs ("libfieldFunctionObjects.so"); writeControl writeTime; timeStart 1; fields ( UField { mean on; prime2Mean on; base time; } alpha { mean on; prime2Mean on; base time; } alpha.air { mean on; prime2Mean on; base time; } U.air { mean on; prime2Mean on; base time; } ); } }
-
@BlookCFD 谢谢您的指点,我有看了一眼程序,我理解错了,OF里面是通过计算barycentric displacement来获得barycentric coordinates的,我一直被src/lagrangian/basic/particle/particle.C里面的一句代码给误导了,现贴出来:
particle.C 1062-1096 OpenFOAM6 void Foam::particle::correctAfterInteractionListReferral(const label celli) { // Get the position from the barycentric data const vector pos(coordinates_.b(), coordinates_.c(), coordinates_.d()); // Create some arbitrary topology for the supplied cell celli_ = celli; tetFacei_ = mesh_.cells()[celli_][0]; tetPti_ = 1; facei_ = -1; // Get the reverse transform and directly set the coordinates from the // position. This isn't likely to be correct; the particle is probably not // in this tet. It will, however, generate the correct vector when the // position method is called. A referred particle should never be tracked, // so this approximate topology is good enough. By using the nearby cell we // minimize the error associated with the incorrect topology. coordinates_ = barycentric(1, 0, 0, 0); if (mesh_.moving()) { Pair<vector> centre; FixedList<scalar, 4> detA; FixedList<barycentricTensor, 3> T; movingTetReverseTransform(0, centre, detA, T); coordinates_ += (pos - centre[0]) & T[0]/detA[0]; } else { vector centre; scalar detA; barycentricTensor T; stationaryTetReverseTransform(centre, detA, T); coordinates_ += (pos - centre) & T/detA; } }
其中的这句代码
coordinates_ += (pos - centre) & T/detA;
其实前面已经先定义了
coordinates_ = barycentric(1, 0, 0, 0);
后面计算单点的barycentric coordinates其实还是用的barycentric displacement的概念来转换的。
困扰了我两周的问题终于解决了,现在异常开心,在此对@东岳 和@BlookCFD 表示万分感谢。
-
上面那个公式引用是公式(8)
-
仔细读了代码,现在还剩最后一点小问题,其实说到底还是原来的问题,barycentric追踪是通过bary displacement元素的正负来判断颗粒运动轨迹与四面体网格面的交点的,例如这个图形
颗粒$\mathbf{p}$运动到$\mathbf{n}$,轨迹$\mathbf{pn}$与面$\underline{\rm{bcd}}$相交于点$\mathbf{s}$,那么矢量$\mathbf{pn}$的barycentric displacement的计算方法为:
\begin{equation}
\begin{aligned}
\label{eq:barycentric displacement}
\boldsymbol{\lambda}_{\mathbf{p}\rightarrow\mathbf{n}}
&=\boldsymbol{\lambda}\left(\mathbf{n}\right)-\boldsymbol{\lambda}\left(\mathbf{p}\right)=\frac{\mathbf{an}\cdot\mathbf{T}}{\rm{det}\mathbf{A}}-\frac{\mathbf{ap}\cdot\mathbf{T}}{\rm{det}\mathbf{A}}=\frac{\mathbf{pn}\cdot\mathbf{T}}{\rm{det}\mathbf{A}}\
&=\frac{\left[\mathbf{pn}\cdot\left(\mathbf{bd}\times\mathbf{bc}\right)\quad\mathbf{pn}\cdot\left(\mathbf{ac}\times\mathbf{ad}\right)\quad\mathbf{pn}\cdot\left(\mathbf{ad}\times\mathbf{ab}\right)\quad\mathbf{pn}\cdot\left(\mathbf{ab}\times\mathbf{ac}\right)\right]}{\mathbf{ab}\cdot\left(\mathbf{ac}\times\mathbf{ad}\right)};.
\end{aligned}
\end{equation}
此时barycentric dispalcement的符号为(-,+++),可以判断轨迹与面$\underline{\rm{bcd}}$相交,同根根据$\mathbf{p}$和$\mathbf{n}$的barycentric coordinates可以求出点$\mathbf{s}$的barycentric coordinates,再转换为点$\mathbf{s}的$gobal coordinates,然后进如下一个四面体网格追踪。
上面我只是简述了颗粒在一个四面体网格内的追踪过程,从计算barycentric dispalcement的角度来看,逆向变换矩阵$\mathbf{T}$是没有问题的,但是问题还是在于从barycentric coordinates to global coordinates的过程,矩阵$\mathbf{T}$的第一个元素有问题。如果把\ref{eq: the first element of T}变为$\boldsymbol\lambda_1=1-\boldsymbol\lambda_2-\boldsymbol\lambda_3-\boldsymbol\lambda_4$就好了,但是程序里没有这么写,说到底还是原来的问题 如果李老师有时间的话可不可以读一下particle类,我觉得您可以帮我释疑 -
@东岳 是的李老师,我的疑惑就在这儿,$\mathbf{bp}\neq\mathbf{ap}$,和代码对不上,细读代码发现不一样。
-
最近在研读OpenFOAM中particle类的barycentric coordinates 和global coordinates转换矩阵的程序,遇到一点问题,现总结出来和大家讨论一下。
OF5.0以后对颗粒追踪采用了barycentric方法,目的是更好的捕捉颗粒轨迹与tetrahedra网格面的交点,四面体内任意一点$\mathbf{p}$的global coordinates 可以表示为barycentric coordinates ($\lambda_i$)与顶点的乘积和:
\begin{equation}
\label{eq:barycentric coordinates}
\mathbf{p}=\sum_{i=1}^{4}\lambda_i\mathbf{v_i}
\end{equation}
其中$\mathbf{v_i}$是四面体顶点坐标。
示意图如下:
如上图,四面体中有一点$\mathbf{p}$,则$\mathbf{p}$点的重心坐标($\lambda_1,\lambda_2,\lambda_3,\lambda_4$)的计算方法为点$\mathbf{p}$将四面体分成的四个小四面体与大四面体的体积之比:
\begin{equation}
\label{eq:lambda1}
\lambda_{1}=\frac{V_{\rm{pbcd}}}{V_{\rm{abcd}}}=1-\lambda_{2}-\lambda_{3}-\lambda_{4},
\end{equation}
\begin{equation}
\label{eq:lambda2}
\lambda_{2}=\frac{V_{\rm{pacd}}}{V_{\rm{abcd}}}=\frac{6V_{\rm{pacd}}}{6V_{\rm{abcd}}}=\frac{\mathbf{ap}\cdot\left(\mathbf{ac}\times\mathbf{ad}\right)}{\mathbf{ab}\cdot\left(\mathbf{ac}\times\mathbf{ad}\right)};,
\end{equation}
\begin{equation}
\label{eq:lambda3}
\lambda_{3}=\frac{V_{\rm{pabd}}}{V_{\rm{abcd}}}=\frac{6V_{\rm{pabd}}}{6V_{\rm{abcd}}}=\frac{\mathbf{ap}\cdot\left(\mathbf{ad}\times\mathbf{ab}\right)}{\mathbf{ab}\cdot\left(\mathbf{ac}\times\mathbf{ad}\right)},
\end{equation}
\begin{equation}
\label{eq:lambda4}
\lambda_{4}=\frac{V_{\rm{pabc}}}{V_{\rm{abcd}}}=\frac{6V_{\rm{pabc}}}{6V_{\rm{abcd}}}=\frac{\mathbf{ap}\cdot\left(\mathbf{ab}\times\mathbf{ac}\right)}{\mathbf{ab}\cdot\left(\mathbf{ac}\times\mathbf{ad}\right)},
\end{equation}
其中 $V_{\rm{pbcd}}$, $V_{\rm{pacd}}$, $V_{\rm{pabd}}$, 和 $V_{\rm{pabc}}$ 是小四面体的体积, $\mathbf{ap}=\mathbf{p}-\mathbf{a}$, $\mathbf{ab}=\mathbf{b}-\mathbf{a}$, $\mathbf{ac}=\mathbf{c}-\mathbf{a}$, $\mathbf{ad}=\mathbf{d}-\mathbf{a}$.
定义从barycentric coordinates到global coordinates的转换矩阵$\mathbf{A}$,则
\begin{equation}
\begin{aligned}
\label{eq:barycentric to global}
\mathbf{p}=\mathbf{A}\boldsymbol{\lambda}\left(\mathbf{p}\right)\rightarrow
\begin{bmatrix}x_{\rm{p}}\\y_{\rm{p}}\\z_{\rm{p}}\end{bmatrix}=
\begin{bmatrix}x_{\rm{a}}&x_{\rm{b}}&x_{\rm{c}}&x_{\rm{d}}\\y_{\rm{a}}&y_{\rm{b}}&y_{\rm{c}}&y_{\rm{d}}\\z_{\rm{a}}&z_{\rm{b}}&z_{\rm{c}}&z_{\rm{d}}
\end{bmatrix}
\begin{bmatrix}\lambda_1\\
\lambda_2\\
\lambda_3\\
\lambda_4\end{bmatrix};.
\end{aligned}
\end{equation}
其中$\lambda_i$就是$\mathbf{p}$的重心坐标,矩阵$\mathbf{A}$就是由四面体各个顶点的global coordinates组成的坐标矩阵。
既然有了正向变换,同样需要逆向变换,定义一个逆向变换矩阵$\mathbf{T}=\left[\mathbf{bd}\times\mathbf{bc}\quad\mathbf{ac}\times\mathbf{ad}\quad\mathbf{ad}\times\mathbf{ab}\quad\mathbf{ab}\times\mathbf{ac}\right]$将global coordinates转换为barycentric coordinates:
\begin{equation}
\label{eq:global to barycentric}
\boldsymbol{\lambda}=\frac{\mathbf{ap}\cdot\mathbf{T}}{\rm{det}\mathbf{A}}=\frac{\left[\mathbf{ap}\cdot\left(\mathbf{bd}\times\mathbf{bc}\right)\quad\mathbf{ap}\cdot\left(\mathbf{ac}\times\mathbf{ad}\right)\quad\mathbf{ap}\cdot\left(\mathbf{ad}\times\mathbf{ab}\right)\quad\mathbf{ap}\cdot\left(\mathbf{ab}\times\mathbf{ac}\right)\right]}{\mathbf{ab}\cdot\left(\mathbf{ac}\times\mathbf{ad}\right)},
\end{equation}
公式\eqref{eq:global to barycentric}的四个组成成分其实就是上面的公式\eqref{eq:lambda1},\eqref{eq:lambda2},\eqref{eq:lambda3},\eqref{eq:lambda4}。$\rm{det}\mathbf{A}$是矩阵$\mathbf{A}$的行列式,表示为$\mathbf{ab}\cdot\left(\mathbf{ac}\times\mathbf{ad}\right)$。
至此,barycentric coordinates与global coordinates的正向变换和逆向变化都已经定义完成了,上面的分析是基于我看OF6的代码以及查阅相关的文献总结出来的。
但是有一个问题困扰了我很久了 ,那就是逆向变换矩阵的第一个元素是:
\begin{equation}
\label{eq:the first element of T}
\boldsymbol{\lambda}_1=\frac{\mathbf{ap}\cdot\left(\mathbf{bd}\times\mathbf{bc}\right)}{\mathbf{ab}\cdot\left(\mathbf{ac}\times\mathbf{ad}\right)},
\end{equation}
这个和公式\eqref{eq:lambda1}的定义不同,因为分子上算的不是图中小四面体pbcd的体积,就这儿我无法理解为什么OF里要这么定义,并且我看OF5-OF8里面都是这么定义的,我觉得有点问题,所有发出来和大家讨论一下,我觉得我自己是解决不了这个问题了 -
@闻久STU 谢谢哦,您这个方法我也记在小本本上了,以后肯定能用上
-
@闻久STU 您好,按照您的步骤,我遇到一点问题,我设置两个block,一个显示流场信息,一个显示颗粒位置信息,然后按照李老师所说,颗粒添加vector滤镜显示颗粒大小,您说的“把选中的数据放进去(左边copy active selection)”,我不太理解,这是我的设置截面,希望您能指导一下,谢谢您。
-
@东岳 谢谢李老师指导
-
我用DPMFoam计算颗粒的管道中的运动,因为对paraview使用不熟练,想后处理颗粒的时候处理成这样的:
请问该如何操作呢,或者有没有老师可以推荐好的颗粒后处理的软件。
这是我的模型:
还有就是想把颗粒随时间的数据做成视频,用paraview可以实现么?谢谢大家,期待能得到指导。 -
@upc_ngh
1.是的,$\mathbf{U}$代表的是KinematicParcel中的$\mathbf{U_c}$,也就是连续相的在颗粒所处位置的速度插值,具体的插值方式OF里面有cell,cellpoint那些,$\mathbf{U}_\rm{d}$表示的就是颗粒的速度。
2.是的,系数乘以速度差。 -
@东岳 李老师,您这个式子(16)的$\mathbf{g}$前面是不是还少一个$\alpha$?
-
@东岳 李老师,关于您这张图中的公式(16)我有几个疑问,如果说得不对望见谅。
1.对于连续相动量方程的拉力梯度项前是否需要乘以一个相分数主要决定于$\mathbf{F}$中的压力梯度力项是否提取出来了,如果提取出来了和前面的压力梯度项组合就会出现乘以相分数的情况,如果不提取出来,把颗粒的压力梯度力放在$\mathbf{F}$中,前面就不用再乘以相分数。不知道我这么解释是否正确?这个解释我是参考的文章Zhou et al. 2010
2.对于您式子(16)中最后一项为何$\mathbf{F}$的前面要乘以一个颗粒相分数,并且为何是除以颗粒的体积,按道理不应该是除以流场网格体积吗?
希望能得到李老师的解答。 -
@刘雄国 嗯,我去看一下OF7的代码,谢谢您
-
李老师@东岳 请问这篇文献出处是?我想研究一下
-
李老师您这篇文章的公式就是我上面式(4)的公式,但是我好像打错了 。总体比较来看是有细微区别。
-
经过仔细比对,发现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,不知道上面写出的欧拉-朗格朗日方程对不对?还有,如果颗粒比网格大了,用单个网格内的颗粒相加除以网格体积是不是就没有意义了,也就是说模型方程是不是就不能再这么写了? -
好的,谢谢李老师,我再研究研究
-
最近在研究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}
两个差了很多。不知道有木有大神解答一下,纠结了一整天了
求教边界层生成问题
OPENFoam里有没有颗粒的全解析求解器
流化床床层高度
OpenFOAM有对单个固体颗粒在流场中的全解析计算的求解器吗
基于扩散的方法解决粒子大于网格
基于扩散的方法解决粒子大于网格
壁面侵蚀边界的实现方法
壁面侵蚀边界的实现方法
OpenFOAM编程findCell的诡异问题,对同一坐标寻找cell结果不一致
topoSet初始能否框出一个三角形区域呢?
topoSet初始能否框出一个三角形区域呢?
compressibleInterFoam求解器怎么加入空化模型
DPMFoam算流化床压降波动过于剧烈
DPMFoam算流化床压降波动过于剧烈
DPMFoam算流化床压降波动过于剧烈
DPMFoam算流化床压降波动过于剧烈
DPMFoam算流化床压降波动过于剧烈
DPMFoam算流化床压降波动过于剧烈
OF5版本以后重心坐标与绝对坐标的相互转换的问题。
OF5版本以后重心坐标与绝对坐标的相互转换的问题。
OF5版本以后重心坐标与绝对坐标的相互转换的问题。
OF5版本以后重心坐标与绝对坐标的相互转换的问题。
OF5版本以后重心坐标与绝对坐标的相互转换的问题。
关于paraview后处理颗粒的问题
关于paraview后处理颗粒的问题
openFoam5.x中,MPPICFoam后处理时无法读入粒子文件
关于paraview后处理颗粒的问题
关于paraview后处理颗粒的问题
particleforce类内ErgunWenYuDrag的表达式感觉有错误
particleforce类内ErgunWenYuDrag的表达式感觉有错误
particleforce类内ErgunWenYuDrag的表达式感觉有错误
particleforce类内ErgunWenYuDrag的表达式感觉有错误
particleforce类内ErgunWenYuDrag的表达式感觉有错误
particleforce类内ErgunWenYuDrag的表达式感觉有错误
particleforce类内ErgunWenYuDrag的表达式感觉有错误
particleforce类内ErgunWenYuDrag的表达式感觉有错误
particleforce类内ErgunWenYuDrag的表达式感觉有错误