VOF添加斥力模型
-
- 遍历所有相分数在0-1直接的网格，并做标记；
- 在标记的网格上，通过方程10判断是要遍历的网格，取第i个网格，算$\bfF_{rep,i}=\sum^N \frac{K}{h_{ij}^3}\bfn$；
- 对i遍历，有做标记网格上的所有的$\bfF_{rep,i}$，加到动量方程或压力方程
你这样debug一下：
用你的2D算例，网格搞到几千，然后算一下，
- 把你的相分数、斥力的云图弄出来，
- 显示斥力矢量图，类似下面这个图（显示的是网格上的表面张力）
首先要把气泡往外推开的过程复现
他们挺有意思，之前还用多相VOF来算这个，一个气泡一哥相，暴力！
Direct numerical simulation of deformable rising bubbles at low Reynolds numbers -
斥力矢量图我晚些再搞一下。
气泡往外推开是复现了的。
现在发现压力参考点附近的异常是结果，原因还是斥力有些突兀。不使用压力参考点也会出现这种压力振荡。
下面是不使用压力参考点的一次模拟。压力振荡是存在的，而且是和斥力的出现是同步的。
最开始的时候气泡弹开了，如果不加斥力，上面两个气泡是会融合的。
3bbsNoReferencevideo1.avi压力振荡最厉害的一段时间：3-4s。
3bbsNoReferencevideo2.avi -
现在大致上解决了，确实是添加的斥力太大的问题，之前代码写得不太对，写成了$K/{h^3}\nabla \alpha $，而正确的表达式应该是$\frac{K}{{{h^3}}}\frac{{\nabla \alpha }}{{\left| {\nabla \alpha } \right|}}$，导致某些情况下两个气泡接近时候，由于界面上某些网格的梯度特别大，导致实际添加的斥力会非常非常大。
现在修改了代码，之前的问题很大程度上解决了，还有一点问题是还是会出现添加斥力比较大的时候，因为总会有气泡不断接近到比较近的时候，从而产生一个比较大的斥力，也会导致全局的压力振荡。就比如下图中这些上下跳的曲线。
表面张力是${\sigma \kappa \nabla \alpha }$，和斥力的表达式还不一样，我再都算一下volScalarField，不往面上插值什么的，比较一下它俩。
-
还有一点问题是还是会出现添加斥力比较大的时候，因为总会有气泡不断接近到比较近的时候，从而产生一个比较大的斥力
是的，在最开始最开始你在讨论这个问题的时候，我就想到了可能会有这个问题。所以就跟你讨论了一下。因为我最近从OpenFOAM学习到一个方法来处理这个问题。
你这个添加斥力的模型，从另外一种角度，可以近似的看做添加一种分散力。同理的一个现象，在气固流动那面，存在一个相压力模型，在相分数超过设定值的时候，添加一个分散力，将alpha打散。正常的算法就是在压力方程中，比如你得处理方式，把这个分散力加进去。但是OpenFOAM那面把这个问题在alpha方程中也处理了。我测试过，如果仅仅是压力方程处理，可能会导致发散（alpha越界）。如果压力方程+相方程同时处理，就非常稳定。
我觉得你可以借鉴这个方法。
http://dyfluid.com/reactingTwoPhaseEulerFoam.html
PS 你这个斥力是$1/h^3$，那面的力是exp指数级别的函数（最终是一个非常陡暴增，力会突然变得高几十个数量级）
-
谢谢李老师提供的思路。我的理解是有两种方法：
- 处理相方程，粗略看了一下，感觉这个方法挺复杂。
- 可以改一改斥力的表达式，让它不那么大。
我在输出两个volVectorField，用来比较斥力和表面张力的大小，具体编程上也想请教一下，控制体体心的alpha的单位梯度怎么写。我目前是这样的，单位梯度必须加一个1e-300防止发散。
// Cell gradient of alpha const volVectorField gradAlpha(fvc::grad(alpha1)); repulsion_Compare = repulsion_Sum*gradAlpha/(mag(gradAlpha)+dimensionedScalar(dimensionSet(0,-1,0,0,0,0,0), 1e-300)); surfaceTensionMy=mixture.sigmaK()*gradAlpha;OpenFOAM在算表面张力的时候，算曲率 $\kappa$ 的时候用到了面心的alpha单位梯度，也分享一下代码。
const surfaceVectorField& Sf = mesh.Sf(); // Cell gradient of alpha const volVectorField gradAlpha(fvc::grad(alpha1_, "nHat")); // Interpolated face-gradient of alpha surfaceVectorField gradAlphaf(fvc::interpolate(gradAlpha)); // Face unit interface normal surfaceVectorField nHatfv(gradAlphaf/(mag(gradAlphaf) + deltaN_)); correctContactAngle(nHatfv.boundaryFieldRef(), gradAlphaf.boundaryField()); // Face unit interface normal flux nHatf_ = nHatfv & Sf; // Simple expression for curvature K_ = -fvc::div(nHatf_); -
我没有写全，就在function里面测试了一下这个力的大小（很多参数都是拍脑袋硬植入）：
functions { setDeltaT { type coded; libs ("libutilityFunctionObjects.so"); writeControl writeTime; executeControl timeStep; code #{ #}; codeExecute #{ const volScalarField& alpha ( mesh().lookupObject<volScalarField>("alpha.air") ); const surfaceVectorField& Sf = mesh().Sf(); const volVectorField gradAlpha("gradAlpha", fvc::grad(alpha)); surfaceVectorField gradAlphaf(fvc::interpolate(gradAlpha)); const dimensionedScalar deltaN_ ( "deltaN", 1e-8/pow(average(mesh().V()), 1.0/3.0) ); surfaceVectorField nHatfv(gradAlphaf/(mag(gradAlphaf) + deltaN_)); volVectorField nHatv("nHatv", gradAlpha/(mag(gradAlpha) + deltaN_)); surfaceScalarField nHatf = nHatfv & Sf; volScalarField K("K", -fvc::div(nHatf)); volScalarField interface("interface", alpha); volVectorField force("force", nHatv); List<label> markedCell; forAll(mesh().C(), i) { force[i] = Zero; if (alpha[i] < 0.99 && alpha[i] > 0.01) { interface[i] = 1.0; markedCell.append(i); } else { interface[i] = 0.0; } } forAll(markedCell, i) { label celli = markedCell[i]; forAll(markedCell, j) { label cellj = markedCell[j]; scalar dij = mag(mesh().C()[celli] - mesh().C()[cellj]); scalar posneg = nHatv[celli] & nHatv[cellj]; if (posneg < 0.0) { if (dij < 0.15) { scalar K = 1e-2; force[celli] += K*nHatv[celli]/pow3(dij); } } } } markedCell.clear(); K.write(); nHatv.write(); gradAlpha.write(); force.write(); interface.write(); #}; } } -
李老师，想请教一下源项的添加的问题。
比如我现在在有一个repulsion源项，要想植入，应该是有三个地方要添加。
一，速度预测方程，即UEqn.H中的fvMatrix UEqn。但我不清楚这里是否要加上源项。
加的理由是：fvOptions这种人工添加的源项都在。
不加的理由是：写的那个icoFoam解析，以及这里实际上只保留了对流项和粘度项，压力、表面张力和重力都没有加，这个斥力在形式上是和表面张力有关的，就像表面张力那样处理。fvVectorMatrix UEqn ( fvm::ddt(rho, U) + fvm::div(rhoPhi, U) + MRF.DDt(rho, U) + turbulence->divDevRhoReff(rho, U) == fvOptions(rho, U) + repulsion_Vector // add );
二，真实速度方程，一般的想法都是加上这个源项，就是最后的速度修正项是
但是在reactingTwoPhaseEulerFoam中，我看到这些源项都是直接加到最后的速度上的，不知道这里是您省略它们的具体形式了，还是这些$M$确实就不除以$ A_{d,p} $
就是这两种方式应该是哪一种？还是两种都可以？
三，上面两种不同方式的速度方程，对应不同的压力泊松方程，也就差个$ A_{d,p} $
-
速度预测方程，即UEqn.H中的fvMatrix UEqn。但我不清楚这里是否要加上源项。
不需要加
源项你直接1）加到phi上面组建压力方程，2）最后的速度更新要加上源项
我那个M都忘记除掉A了，是要除掉，我有遗漏，我周末更新过去
http://dyfluid.com/reactingTwoPhaseEulerFoam.html 这个里面41-46之前写的一直没更新 29之前是我最近更新的 没问题 -
