Stokes流中的圆球绕流的阻力系数，怎么设置才能算准。
-
东岳老师好，各位大佬、老师同学好：
最近需要做有关蠕动流(Re≤1)时的圆球绕流的阻力系数验证工作，但是曳力系数一直都算不准，算出来的阻力系数不满足Stokes定律Cd=24/Re；
因为我目前仅验证的是Re=1的情况，得到过的阻力系数基本上都是在30+，不管换ico, simple, pimpleFoam求解得到的结果都不对，下图是用pimpleFoam的求解结果，继续算下去仍然是降不下去的：
理论值(24/Re)应该接近24，我也去找了József Nagy演示圆球绕流的文件来看，发现他算的Re=1时阻力系数结果是37.4453，在这里附上我用的一个Case文件，希望各位指点迷津，感谢！
dyfluidAttachment.zip- 另外：ControlDict中functions内的forces_object和forceCoeffs_object我特地在圆柱绕流时验证过，我另外还拿Paraview自己测量过压差阻力，两个函数forces_object和forceCoeffs_object给出的数据应该是没错的，所以问题点应该不是在functions这里。*
-
@李东岳 我也不知道呀，东岳老师。这个是JózsefNagy.Case.zip 老师用的Case，他对此做了很大雷诺数范围的圆球绕流阻力系数验证，咋一看Re=1时结果挺好，但是切回到线性坐标看的话，Re=1时的阻力系数是37.4453（József Nagy老师在他的网页给的具体数据）
当然，József Nagy老师可能是以教学示例为主，没有过多讨论这个事。
看到其他论文用COMSOL算的结果，馋哭了，不得不说吻合得真好（目前使用openFoam是导师给的限定条件，怎么能向商软低头！）。
我也非常懵，Stokes流（Re≤1）的圆球绕流仿真还没有找到比较好的参考案例。
我最近在尝试看画个楔形网格来算。 -
@李东岳 东岳老师好， 我现在用二维楔形网格算出来的圆球绕流结果挺好，基本上符合了Schiller-Newman 公式(Cd=24/Re*(1+0.15*Re^0.687)). 算出来的Cd值大概是27.6(SN公式也是27.6)
我主要想验证的事是一些文献谈到的：
①颗粒-界面无滑移noSlip-则Cd值符合Stokes解析解Cd=24/Re(Schiller-Newman 公式是Stokes定律考虑了惯性项的修正公式）,
②纯净气泡-界面滑移slip-则Cd值符合H-R理论解Cd=16/Re。我计算域的边界设置就如同下图
我想的是只需要改变小孔处界面滑移性质(noSlip/slip)就能在球体颗粒/气泡之间切换。
基于此，颗粒算出来的Cd值大概是27.6，但是，在其他文件什么内容都不变的情况下，只将U,P文件的边界条件从noSlip改为了slip（也就是从颗粒转到气泡）
得到的阻力系数就有问题了，太低了，只有13+，理论来说应该16才对呀，我现在又不指导问题出在哪里了。
请东岳老师和各位大佬再度指点迷津，多谢！
我在这里上传了case文件，再次谢谢东岳老师。
Re.1.MB.WMesh.AllSlip.zip -
@李东岳 东岳老师好，目前明确提到这种处理方式的文献，我整理了一下，后续还有相关文献我再继续补充：
-
Bel Fdhila, R., Duineveld, P.C., 1996. The effect of surfactant on the rise of a spherical bubble at high Reynolds and Peclet numbers. Phys. Fluids 8, 310–321.
//一些学者认为是这篇文献的作者团队Duineveld是最先这么做处理的（还有的认为是改变界面滑移性质(noSlip/slip)就能在球体颗粒/气泡之间切换。） -
Kishore, N., Nalajala, V.S., Chhabra, R.P., 2013. Effects of Contamination and Shear-Thinning Fluid Viscosity on Drag Behavior of Spherical Bubbles. Ind. Eng. Chem. Res. 52, 6049–6056.
//这篇文章虽然是做非牛顿流体的，但是也对牛顿流体中的球体做了这种处理方式，并做了对比验证（改变界面滑移性质(noSlip/slip)就能在球体颗粒/气泡之间切换。） -
庞明军, 费洋, 陈小洪, 郭雨晨, 徐梦沁, 2019. 雷诺数和界面污染程度对气泡水动力学特性的影响. 农业工程学报 35, 99–105.
*费洋, 庞明军, 2017. 球形气泡界面变化对尾涡性质和尺寸的影响. 化工学报 68, 3409–3419.
Fei, Y., Pang, M., 2018. The influence of interface contaminated degree on the wake characteristics of a spherical bubble at moderate Reynolds number under the condition of isothermal flow. International Journal of Heat and Mass Transfer 121, 79–83.
常州大学的庞明军老师团队经常在Fluent中采用这种处理方式，得到的结果也非常好（改变界面滑移性质(noSlip/slip)就能在球体颗粒/气泡之间切换。）
庞明军老师团队⬇（两篇文献英汉都发了一下，中间那个CD的图在中文文献中有）
Duineveld团队⬇
以及其他学者的验证⬇
东岳老师，我最近在OpenFOAM的楔形网格，以及用SHM挖一个空心球体并铺设边界层的全三维网格里都继续算了一下，依然是：
①颗粒（无滑移界面noSlip）的阻力系数和Schiller-Newman 公式(Cd=24/Re*(1+0.15*Re^0.687))符合得挺好
②气泡（界面滑移slip）的时候，阻力系数值就很低，只有13~14左右，按照相关的经验公式或者理论公式来看，是不能低于16的，最好是在17.6~18（Mei公式-1994年和Taylor公式-1964年）。然后fvScheme也按照文献里说得格式调了调，并没有多大的区别。
我这段时间在想是不是自己对滑移条件的理解不正确，也就是说文献中说的边界滑移和OpenFOAM的slip边界条件会不会不是一回事，但是看下来感觉是一样的，没能进一步看出个所以然来 。再次有劳请大家指点迷津，多谢！
-
-
@李东岳 谢谢东岳老师！
之前上传的代码里用于gnuplot作图的AlldragPlot里面的这句代码要改一下 plot "./postProcessing/forceCoeffs/1/coefficient.dat" every ::10 using 1:3 w lp pt 5 axis x1y1 title "Curve of Drag Coefficient" 要把文件路径的1改为0； plot "./postProcessing/forceCoeffs/0/coefficient.dat" every ::10 using 1:3 w lp pt 5 axis x1y1 title "Curve of Drag Coefficient" -
@李东岳 东岳老师，我全部作无量纲处理了。
- 气泡直径d: 2m
- 速度U: 1m/s
- 粘度nu: 2 /m2s
因此Re=d*U/nu=1×2/2=1
同时我用量纲条件下计算过，真实气泡尺寸，超纯水粘度和微气泡浮升速度
- 气泡直径d: 1.2424e-4m ⏩124.24μm
- 速度U: 8.2745338e-3m/s ⏩8.2745338 mm/s
- 粘度nu: 0.000001028028079312 /m2s ⏩接近真实情况下超纯水粘度
因此Re=d*U/nu=1.2424e-4×8.2745338e-3/0.000001028028079312=1
上述微气泡直径124.24μm的速度取8.2745338 mm/s我想是合理的，因为更关注的是Re，同时大家粘度各不一样，所以速度稍微会有一点区别。
上图对比速度的数据引用于 Pawliszak, P., Ulaganathan, V., Bradshaw-Hajek, B.H., Manica, R., Beattie, D.A., Krasowska, M., 2019. Mobile or Immobile? Rise Velocity of Air Bubbles in High-Purity Water. J. Phys. Chem. C 123, 15131–15138.以下是有量纲和无量纲的两种情况下的阻力系数结果对比
以及采用实际量纲的case文件 ，请东岳老师再帮忙看看，谢谢东岳老师。
Re.1.MB.WMesh.AllNoSlip.Dimension.zip -
@李东岳 好的，东岳老师，这是调整好后的带量纲的气泡case文件（刚才我上传的是带量纲的颗粒case文件）
- 气泡直径d: 1.2424e-4m ⏩124.24μm
- 速度U: 8.2745338e-3m/s ⏩8.2745338 mm/s
- 粘度nu: 0.000001028028079312 /m2s ⏩接近真实情况下超纯水粘度
因此雷诺数等于1，用的是层流模型。Re=d*U/nu=1=(1.2424e-4)×(8.2745338e-3)/0.000001028028079312
- 网格已经用transformPoints将每单位长度缩放到气泡半径R:6.212e-5m⏩62.12μm，直径为124.24μm；
- deltaT用的0.0001，能保证计算过程中的最大Co一直小于0.8；
- 总时间换成了 0.75s；（d/U×50(50是之前无量纲case的总时间) = 0.750737s）
Re.1.MB.WMesh.Bubble.AllSlip.Dimension.zip
请东岳老师再帮忙看看，谢谢东岳老师。 -
@李东岳 东岳老师，稳态之前我用simpleFoam和icoFoam算过（但算的无量纲情况，得到的Cd值也是13.5左右）
刚才做了个带量纲的simpleFoam，结果显示如下，最终Cd值基本持平在13.496左右；
这里是带量纲的simpleFoam案例文件
Re.1.MB.WMesh.AllSlip.Dimension.simpleFoam.zip当时没有选择继续用simpleFoam的原因主要是：
-
①看到其他文献，比如庞明军老师团队的文献，用的瞬态求解器（他们的雷诺数稍大一些，Re最大到200，这个雷诺数范围内依然是符合轴对称流场结构）
-
②找了个wolfdynamics的培训视频中的圆柱绕流的simpleFoam案例来改，但算气泡的收敛判据我不知道怎么给才合适，只是看到simpleFoam一直跑，Cd值大概持平到13.5左右；主要责任是我自己也没往这里多想，因此之前主要都是用pimpleFoam在算。
再次感谢东岳老师
-
-
@李东岳 我在办公室差点哭出声。
依然十分感谢东岳老师，我想再转到Fluent算一下看看情况，有新的进展我再来补充上。
-
@Prometheus10 openfoam这面是没看出啥破绽。我捉摸能不能换个附加平均速度的周期性的网格试试。但是牵涉到画网格工作量比较大。不过你先试试fluent吧。我们再讨论 :-)
-
@Prometheus10 算了十秒还在33，我是不是应该放弃了，网格质量是不是不是特别好呀 感觉不应该呀，我之前用icoFoam算过雷诺数40的圆柱绕流，和文献对的还挺准
-
@Prometheus10 看到你后面解决啦，可能定常流算不准最可能的问题就是网格问题
-
@pengdi @李东岳 感谢pengdi大佬，感谢东岳老师。您们的讨论很大程度给了我帮助，已经很大程度解答了我的疑惑。
但是！
@pengdi @李东岳 东岳老师和pengdi大佬，我觉得OpenFOAM的slip条件这里，还不能直接否定它，还不该锤它。
我想再请两位老师看看我的这个理解过程，我还拿不准一些地方理解得对不对：
先上图：
（这个图是我自己做的，图和公式都是按照自己的理解分别拼一起的，特别是第三个切应力为0的那个表达式，我甚至觉得他写好像错了，但是我先表达一下我的逻辑和想法，请您们帮忙看看：）
因为OpenFOAM应该内置的是：
- Slip的底层代码执行的内容是：被称为perfect-slip或free-slip的理想滑移条件（不仅要约束法向速度为0，还要约束切应力为0），这一点，OpenFOAM没有错。
- noSlip的底层代码执行的内容是：作为无滑移noSlip（界面处速度为0）；这一点，OpenFOAM也没有错。
然而，我在这个帖子内提到的气泡界面“滑移”条件，虽然经常在各种文献中被大家称为“滑移slip”，但是它其实只需要一个约束条件（也就是中间这个，只约束法向速度为0，并不去约束切应力）。
不少文献将其误称为free-slip（比如常州的庞老师那篇文章....可能因为他们雷诺数比较大，然后在大雷诺数时这个边界条件带来的影响不是太明显，又或者是笔误...我很尊敬庞老师并想求签名...先在这里提前为我的comments抱歉。）但不管怎么说，这恐怕是一个错误的叫法。
因为，从下面的文献中的气泡界面法向速度分布的示意图可以看出，界面是有速度的，但是速度达不到主流速度那么大（我的同学和我讨论的结果，认为是Re=0.1时气泡界面速度大概约为主流速度的一半），因此气泡界面边界条件所指的“滑移”必然不是free-slip。
那么，符合这个帖子内提到的气泡界面“滑移”条件（也就是只约束法向速度为0，并不去约束切应力），反而是symmetry边界条件准确描述了这个问题。
（tomiyama这个图具体是需要看最下面的No Marangoni Stress这个地方的法向速度分布图）参考文献： Liu, B., Manica, R., Xu, Z., Liu, Q., 2020. The boundary condition at the air–liquid interface and its effect on film drainage between colliding bubbles. Curr. Opin. Colloid Interface Sci. 50, 101374. Hosokawa, S., Masukura, Y., Hayashi, K., Tomiyama, A., 2017. Experimental evaluation of Marangoni stress and surfactant concentration at interface of contaminated single spherical drop using spatiotemporal filter velocimetry. Int. J. Multiphase Flow 97, 157–167.----------------------------------------分割线：回复Pengdi老师------------------------
因此@pengdi ，PengDi老师如果您需要用slip如果是free-slip，那您的仿真其实并没错。
----------------------------------------分割线：新的一些问题和想法------------------------
然后，我现在还存在的疑惑是：
u·n=0是Partial-Slip吗？我看相关文档表示，OpenFOAM内置的Partial-Slip是继承于Slip（free-slip），然后给一个系数，感觉它主要是将界面速度执行为主流速度*系数，对于切应力方面，依然继承了切应力为0的约束条件。
因此我认为u·n=0并不是OpenFOAM中的Partial-Slip，这个观点是否正确？
再次谢谢东岳老师@李东岳 和Pengdi老师@pengdi
-
@Prometheus10 在 Stokes流中的圆球绕流的阻力系数，怎么设置才能算准。 中说：
，虽然经常在各种文献中被大家称为“滑移slip”，但是它其实只需要一个约束条件（也就是中间这个，只约束法向速度为0，并不去约束切应力）。
我认为主要是因为symmetry边界限制了法向速度，slip边界只能保证切向速度相同，但是不能限制法向速度，可能和你的观点差不多。今天我会做个算例来验证一下。是否应该限制法向速度我觉得也不太好说，可能需要结合具体的问题开展分析。
-
-
@xiezhuoyu 您好 @xiezhuoyu ，先尝试回答您的提问：
这个图带来了误解（这个图其实是partial-Slip的图，我放这个图主要是想请大家帮忙厘清下u·0和气泡“滑移界面”以及partial-Slip之间的关系）。那么，回到问题上来：
u·η=0（也就是symmetry边界的数学表达，COMSOL里可以看到对称边界的数学方程就是u·η=0）的作用是什么？ 我想象对于一个气泡界面处的流体而言，施加这个u·η=0的条件，并不是约束沿气泡界面法向的速度变化，而是要求流体沿壁面“贴体”进行运动，也就是说靠近壁面处的流体，不能沿法向运动，只能沿着切向运动。
但是！切向有没有运动，需要来看切应力的大小，也就是看有没有对切应力做约束。展开来说：如果此时在u·η=0的条件上继续添加切应力为0的约束条件，那么就形成free-slip边界条件；反之，如果添加的是切应力无限大的约束条件，那么就形成noSlip边界条件（u=0）。 -
妥了，实锤了，就是slip与symmetry的问题，估计slip没有做到完全的无剪切。我这几天详细看一下
针对我上面说的，这个不对，slip可以做到无剪切。
在测试算例，我尝试把球表面的值输出出来，第一个是slip，后面是symmetry：
U & n = -4.2351647e-22 max(U & n) = 4.3368087e-19 average(U & n) = -2.494512e-20 average(mag(U)) = 0.0045558926 p face value = 0.0002584113 p cell value = 0.0002584113 average(p face) = average(p) [0 2 -2 0 0 0 0] -1.1678857e-07U & n = 4.2351647e-22 max(U & n) = 4.3368087e-19 average(U & n) = -8.4957405e-21 average(mag(U)) = 0.0030605553 p face value = 0.00026689782 p cell value = 0.00026689782 average(p face) = average(p) [0 2 -2 0 0 0 0] -5.620298e-09简单嫩看出：
- 边界处slip跟symmetry的$\bfU_f\cdot\bfn_f\approx 0$，所以二者都可以把边界处法向速度限制为0。
- 压力严格的$\nabla p\cdot\bfn_f=0$，但是可以看出压力的值差距还是不小。看average(p face)这个值
- 因为slip是衍生来自symmetry，因此切向方向二者处理是一样的，都是零梯度。
所以我觉得slip与symmetry，二者都可以满足数学上的法向速度为0以及切向零梯度（在这里我不太清楚sci里面说的free slip是不是这个意思）。但是在求解的时候边界条件导致算法求解流程导致结果不一致，这体现在
average(p face)以及average(mag(U))都有差异。在求解的时候，边界条件导致算法求解流程导致结果不一致
这方面需要更详细看一下。
-
@pengdi Pengdi老师，@李东岳 东岳老师
我来补充一下低雷诺数条件下（Re=0.1时）对圆球型壁面设置为symmetry或者slip边界的切向速度的结果：
依然是Wedge楔形网格的圆球绕流文件，然后将仿真结果在Paraview中用 plotOverLine提取：首先是symmetry边界条件，第一个网格体心的速度大概是主流速度的1/2；
然后是slip边界条件，其第一个网格体心的速度就是主流速度；
沿用一下这个示意图，symmetry边界也就是u·n=0确实很好描述了气泡界面处的速度分布，和Multiphase flow handbook第25页给出的公式也很吻合。(我这里给出的结果，准确来说，应该要说成是紧挨界面处，而不是界面上的切向速度，因为我只会用Paraview显示网格体心速度，还没能掌握如何用Info速度界面上的切向速度)
Michaelides, E., Crowe, C.T., Schwarzkopf, J.D., 2016. Multiphase flow handbook. CRC Press.
参考文献： Liu, B., Manica, R., Xu, Z., Liu, Q., 2020. The boundary condition at the air–liquid interface and its effect on film drainage between colliding bubbles. Curr. Opin. Colloid Interface Sci. 50, 101374. -
@Prometheus10 您好，关于圆柱绕流的slip和symmetry边界计算我也计算对比过，我计算的雷诺数等于200，差别确实很大。我怀疑可能正方形绕流差别就没那么大了，今天比较忙，还没验证，我看明天能不能验证一下，稍后补充上来。
。CFD中文网2016-2023 | 京ICP备15017992号-2