alpha.water越界了?你看最小值都负零点零几了,可能是边界条件没设置好,其它物理量的边界条件也会影响。
学流体的小明
帖子
-
-
那我不知道咋办了,要不你用apt删了重装一下,或者编译源代码安装。
或者你电脑提前改变了apt的安装路径,安装到了其它地方?
下面这个里面安装的openfoam9就是在/opt目录下面。
https://blog.csdn.net/sagjhdj/article/details/123435344 -
apt这种安装方式应该会安装到/opt目录下面吧,你看看那里。
-
编程。
// 先创建一个场 volVectorField source_Vector ( IOobject ( "source_Vector", mesh.time().timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), mesh, dimensionedVector(dimensionSet(1,-2,-2,0,0,0,0), Zero) ); // 根据你的需要进行赋值 forAll(mesh,cellid) { source_Vector[cellid]=........; } // 添加方案一,添加到UEqn.H当中的动量方程中 fvVectorMatrix UEqn ( fvm::ddt(rho, U) + fvm::div(rhoPhi, U) + MRF.DDt(rho, U) + turbulence->divDevRhoReff(rho, U) == fvOptions(rho, U) + source_Vector // Add the source field ); // 添加方案二,添加到泊松方程里面,需要先插值为surfaceScalarField surfaceScalarField source_phig=-fvc::interpolate(source_Vector ); // 然后添加到pEqn的phig surfaceScalarField phig ( ( mixture.surfaceTensionForce() + source_phig // Add the source field - ghf*fvc::snGrad(rho) )*rAUf*mesh.magSf() );
我这个是改的interFoam,你类比一下就行,量纲应该不会一样。
-
以我对C++编程三脚猫程度的理解,应该是basicThermo.C没有被编译,因为在solidThermo.H/C中只包含了头文件basicThermo.H。
你可以试试在solidThermo的Make/files中这样写basicPorousThermo.C porousThermo.C LIB = $(FOAM_USER_LIBBIN)/libporousThermo
这样应该不需要你复制代码吧,也许(很没有信心)。
-
看报错信息是porousThermo没有找到basicPorousThermo,你在options中把basicPorousThermo链接上试试?
-
可以从边界上的face返回这个face属于的cell
const Foam::fvBoundaryMesh &procBoundary = mesh.boundary(); forAll(procBoundary, patch) { forAll(procBoundary[patch], facei) { const label bcell = boundaryMesh[patch].faceCells()[facei]; } }
你可以得到所有的边界单元编号,然后做你其它向做的事就应该比较方便了
-
@csj1246957849
因为提取等值面的时候进行了插值。
可以提取xz截面,显示网格线,然后你用我黄色画出来的看界面附近的网格,OpenFOAM用的有限体积法用到的是cell中心值,当你用setFields赋值的时候,z=0这条线和赋到的cell有各种关系。具体情况我这里表述很麻烦,你实际看一眼应该就能理解。
你上面图片显示的都是point的值,是paraview插值出来的,cell值前面是一个小方块,你可以看看cell值,就没有红色和蓝色的过渡。 -
interFoam好像不能用backward的ddt吧?应该是使用了Crank-Nicolson?
-
我刚好也用过window,window的时长还是要和你实际的流场有关,可以先用probe工具输出一下单点的监测历史,确定一个合适的window时间。
window 5;
意思是当前时间步往前的5s都是window的范围。举例:t = 1s,平均操作是对0-1s的数据;t = 6s,平均操作是1-6s的数据。
是不是平均的时间太长了? -
额,我那个帖子也是说的这个问题呀,你先串行自适应加密几次,然后把串行的网格直接替换constant/polyMesh就行了
-
继续往下计算的话,应该会算出来圆弧边的界面?
不过这个自适应加密对alpha.water的处理确实有点粗暴啊,我之前用的时候也没注意这一点。会不会是哪里没设置好?看看你的dynamicMeshDict。
还有就是fvSolution里面不要显式地把correctPhi改为no,自适应加密之后确实是需要进行通量修正的。
我的做法是先进行预加密,然后把在背景网格上已经有一定加密的网格文件直接替换掉constant/polyMesh。
https://www.cfd-china.com/topic/6177/paraview查看自适应加密网格出错 -
我也只是遇到过相同的自适应网格重组时候有网格问题,倒是没做过这种重叠网格的操作😂。感觉不太好搞,你得学习一下网格文件的结构。
-
感觉大部分动网格的类似问题都可以用这种方法处理:
先使用reconstructParMesh,再使用reconstructPar。
然后重组后的网格就在你重组的那个时间步里。 -
我前两天在解决我的问题的时候,也看了一下。
interfaceProperties类的对象mixture,自带可以输出sigmaK(),输出的是$ \sigma K $,表面张力乘以曲率。不知道你是单独把曲率算出来了还是输出的sigmaK()。
我看到的sigmaK()和你看到的是一样的,到处都有。那个图丢了,我后来也没再输出这个场。
但是我输出了表面张力的volVectorField,就是mixture.sigmaK()*fvc::grad(alpha1),这种表面张力就仅出现在界面上了。
内部曲率会对动量方程中表面张力以及相传输方程造成影响吗。
不会影响。
-
不需要加
源项你直接1)加到phi上面组建压力方程,2)最后的速度更新要加上源项
这种方法算得的结果就是那些有压力振荡不稳定的。
另一种方法:仅添加到
UEqn
当中,压力方程不做处理,斥力的贡献依靠UEqn.H
参与到后续的计算当中。结果就没有压力振荡了。
目前初步的结果是这样的,后续可能还有其它问题。
非常感谢李老师这么久这么多的指点,帮助太大了。
-
李老师,想请教一下源项的添加的问题。
比如我现在在有一个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} $
-
之前试算了很多,发现压力参考点不是最根本的原因,最根本的原因还是添加了这个斥力源项导致的流场压力不稳定。不用压力参考点也算出来了这种压力振荡。
应该是和具体的代码有关系,和SIMPLE算法里面预测方程有关系。
我没有往UEqn中加这个斥力项,仅加在了phig里面,也就是在压力泊松方程里面加了,后面速度修正那一行代码中,斥力也是体现在phig当中。
李老师您怎么写的代码? -
K调太小的话,气泡间斥力无法弹开气泡,两个气泡会融合。
现在还挺混乱下面这张图,左侧轴是压力测点的压力,右侧轴是volScalarField repulsion_Sum的最大值。- 有的时候没有斥力,也有全局的压力下降,比如下面的椭圆圈出来的。
- 有的时候斥力生效,但是压力没问题。方框中基本水平的蓝色和橙色线。
- 总体来说斥力生效会极大地增加压力出问题的风险。
红色箭头是斥力,白色箭头是表面张力,斥力还是比表面张力大一些。
-
谢谢李老师提供的思路。我的理解是有两种方法:
- 处理相方程,粗略看了一下,感觉这个方法挺复杂。
- 可以改一改斥力的表达式,让它不那么大。
我在输出两个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_);
-
现在大致上解决了,确实是添加的斥力太大的问题,之前代码写得不太对,写成了$K/{h^3}\nabla \alpha $,而正确的表达式应该是$\frac{K}{{{h^3}}}\frac{{\nabla \alpha }}{{\left| {\nabla \alpha } \right|}}$,导致某些情况下两个气泡接近时候,由于界面上某些网格的梯度特别大,导致实际添加的斥力会非常非常大。
现在修改了代码,之前的问题很大程度上解决了,还有一点问题是还是会出现添加斥力比较大的时候,因为总会有气泡不断接近到比较近的时候,从而产生一个比较大的斥力,也会导致全局的压力振荡。就比如下图中这些上下跳的曲线。
表面张力是${\sigma \kappa \nabla \alpha }$,和斥力的表达式还不一样,我再都算一下volScalarField,不往面上插值什么的,比较一下它俩。
-
https://zhuanlan.zhihu.com/p/66444729
之前在这里学习了一下,感觉讲得挺好的 -
网格上的表面张力是一个surfaceScalarField,要显示的话,看其他人的方法是使用foamToVTK先转化,再用paraview看,感觉有点麻烦。李老师您有什么快速便捷的方法吗?
-
也不算是空泡,其实第二个视频中也可以看到有斥力的那对泡,左边的泡变小了,反正又是一个很奇怪的问题,只是在有压力参考点的时候,缺少的那部分气体跑到了压力参考点附近。
压力不稳是最关键的问题,感觉解决它之后,第一个问题应该也能解决。 -
斥力矢量图我晚些再搞一下。
气泡往外推开是复现了的。
现在发现压力参考点附近的异常是结果,原因还是斥力有些突兀。不使用压力参考点也会出现这种压力振荡。
下面是不使用压力参考点的一次模拟。压力振荡是存在的,而且是和斥力的出现是同步的。
最开始的时候气泡弹开了,如果不加斥力,上面两个气泡是会融合的。
3bbsNoReferencevideo1.avi压力振荡最厉害的一段时间:3-4s。
3bbsNoReferencevideo2.avi -
更新一下,顺便就新的问题求助。
我现在全场加密之后算气泡,然后添加了气泡间斥力,是添加在了phig当中。pEqn.H surfaceScalarField phig ( ( mixture.surfaceTensionForce() + repulsion // 防止气泡融合的斥力 - ghf*fvc::snGrad(rho) )*rAUf*mesh.magSf() );
然后发现当斥力生效的时候,监测的压强出现了我上面帖子中同样的情况,流场中的压力同步地增大或减小,如下图3.5s、4.5s附近的压力曲线。
Fig1. 多个probe的p_rgh随时间的变化流场如下图,这是一个二维的算例(三维的一样有这个问题,所以拿二维的debug),三个泡在槽道中流动,无重力,上下边界是周期边界,左右是slip的wall,有压力参考点,位置就在左边图的左侧边界附近的“小气泡”的旁边。时间是3.5s,对应Fig1.压力曲线剧烈波动的那一阶段。
左边的图是alpha.water,中间图是p_rgh,右边图是我添加的斥力。
斥力仅在两个气泡接近到一定程度时,加在气泡边界上,就是右图中的两个小竖条。问题在于斥力生效的时候,压力参考点附近很奇怪的冒出了气体,这些气体似乎是从1号气泡上抽出来的,由于压力参考点附近有了气体,相当于有了一个气泡,所以由于表面张力,流场中的液体的压强就会下降,可以看到中间图中流场大部分p_rgh是负的几十帕,压力参考点附近有个红点。
这是Fig1.中流场压力同步地减小的直接原因。
Fig2. 流场图像。左,alpha.water;中,p_rgh;右,我添加的斥力源项。除了上述的“抽离气体”行为之外,还有可能凭空生成一些气体,表现为全场积分得到的气体体积超过初始生成的气体体积。
还可能是压力参考点附近出现alpha.water>1的网格,也会导致有grad(alpha.water),然后导致流场其它地方的压力下降。我之前采用自适应加密时候,看到的监测点压力曲线同步增大减小的直接原因就是这种alpha.water>1的现象。所以直接原因是压力参考点附近出现了气体,更深层次的原因我感觉是和流场的连续性有关。看计算的log文件,发现出现问题的这一个时间步,time step continuity errors明显大一些。下面是log文件,第一个时间步是没有斥力生效,第二个时间步时候斥力生效了,time step continuity errors也比第一个时间步大了两个数量级。
Courant Number mean: 0.062496661 max: 0.49675831 Interface Courant Number mean: 0.00012625052 max: 0.3551864 deltaT = 5.565631e-06 Time = 0.02073116103 PIMPLE: iteration 1 smoothSolver: Solving for alpha.water, Initial residual = 0.0006698861, Final residual = 2.4349078e-13, No Iterations 3 Phase-1 volume fraction = 0.99579311 Min(alpha.water) = 1.2926026e-11 Max(alpha.water) = 1.0000346 MULES: Correcting alpha.water Phase-1 volume fraction = 0.99579311 Min(alpha.water) = -0.0043784079 Max(alpha.water) = 1.0000346 smoothSolver: Solving for alpha.water, Initial residual = 0.00066896527, Final residual = 2.4380054e-13, No Iterations 3 Phase-1 volume fraction = 0.99579311 Min(alpha.water) = 1.2926035e-11 Max(alpha.water) = 1.0000345 MULES: Correcting alpha.water Phase-1 volume fraction = 0.99579311 Min(alpha.water) = -0.0043784079 Max(alpha.water) = 1.0000345 GAMG: Solving for p_rgh, Initial residual = 0.055540806, Final residual = 7.6812173e-05, No Iterations 3 Pressure gradient source: uncorrected Ubar = 1.0363001, pressure gradient = 333.76278 ~~time step continuity errors : sum local = 5.1237013e-08, global = -3.751065e-20, cumulative = -4.8212224e-09 The pressure at the reference cell is corrected from 0.13807379 to 0~~ GAMG: Solving for p_rgh, Initial residual = 0.0012770611, Final residual = 1.733568e-06, No Iterations 5 Pressure gradient source: uncorrected Ubar = 1.0363001, pressure gradient = 335.51285 time step continuity errors : sum local = 1.1477203e-09, global = -3.6709073e-20, cumulative = -4.8212224e-09 The pressure at the reference cell is corrected from 0.0099254384 to 0 GAMG: Solving for p_rgh, Initial residual = 0.00012244276, Final residual = 9.5463703e-09, No Iterations 43 Pressure gradient source: uncorrected Ubar = 1.0363001, pressure gradient = 334.6004 time step continuity errors : sum local = 6.3417109e-12, global = -3.7236763e-20, cumulative = -4.8212224e-09 The pressure at the reference cell is corrected from 0.018604574 to 0 ExecutionTime = 112420.57 s ClockTime = 112564 s fieldAverage fieldAverage write: Calculating averages volFieldValue gas_in_domain write: volIntegrate(region0) of alpha.air = 6.7310313e-08 Courant Number mean: 0.062903353 max: 0.50213272 Interface Courant Number mean: 0.00012710096 max: 0.37880418 deltaT = 5.5439467e-06 Time = 0.02073670498 PIMPLE: iteration 1 smoothSolver: Solving for alpha.water, Initial residual = 0.00066738165, Final residual = 2.5095471e-13, No Iterations 3 Phase-1 volume fraction = 0.99579311 Min(alpha.water) = 1.292604e-11 Max(alpha.water) = 1.0000345 MULES: Correcting alpha.water Phase-1 volume fraction = 0.99579311 Min(alpha.water) = -0.0043784079 Max(alpha.water) = 1.0000345 smoothSolver: Solving for alpha.water, Initial residual = 0.00066646251, Final residual = 2.5131431e-13, No Iterations 3 Phase-1 volume fraction = 0.99579311 Min(alpha.water) = 1.2926051e-11 Max(alpha.water) = 1.0000345 MULES: Correcting alpha.water Phase-1 volume fraction = 0.99579311 Min(alpha.water) = -0.0043784079 Max(alpha.water) = 1.0000345 GAMG: Solving for p_rgh, Initial residual = 0.75668049, Final residual = 0.00077440161, No Iterations 4 Pressure gradient source: uncorrected Ubar = 1.0362943, pressure gradient = 669.91266 ~~time step continuity errors : sum local = 1.9935763e-06, global = -6.2120888e-09, cumulative = -1.1033311e-08 The pressure at the reference cell is corrected from 5523.05 to 0~~ GAMG: Solving for p_rgh, Initial residual = 0.0010213942, Final residual = 2.0013854e-06, No Iterations 69 Pressure gradient source: uncorrected Ubar = 1.036315, pressure gradient = -548.48592 time step continuity errors : sum local = 1.5198456e-08, global = -6.2120888e-09, cumulative = -1.72454e-08 The pressure at the reference cell is corrected from 5341.8877 to 0 GAMG: Solving for p_rgh, Initial residual = 7.1997514e-05, Final residual = 8.5095615e-09, No Iterations 86 Pressure gradient source: uncorrected Ubar = 1.036315, pressure gradient = -549.38251 time step continuity errors : sum local = 6.2499269e-09, global = -6.2120888e-09, cumulative = -2.3457489e-08 The pressure at the reference cell is corrected from 5299.7144 to 0 ExecutionTime = 112459.83 s ClockTime = 112604 s fieldAverage fieldAverage write: Calculating averages
所以有大佬遇到过相同的问题吗?求指点,救救孩子吧,卡在压力这个问题上大半年了🥺🥺🥺🙏🙏🙏
-
Short-range repulsive force model for near-contact interactions of bubbles
Lingxin Zhang, Kai Peng, Xueming Shao, and Jian Deng
Phys. Rev. E 106, 045110 – Published 25 October 2022
https://doi.org/10.1103/PhysRevE.106.045110 -
在用二维算例debug了,应该是斥力加得太突兀的原因,也许
继续求指点 -
fvSolution里面的pFinal,relTol应该调到0吧。
矩阵最终求解的这些*Final一般都强制要求满足容差。
应该是这个原因导致你某一步没收敛,或者误差逐渐积累,导致发散。 -
学流体的小明 在 interFoam计算气泡槽道流时的压力问题 中说:
应该是 并行 + 自适应网格 的问题?
更正一下,是 并行+自适应网格+压力参考点 三个要素的共同作用,去掉任何一个都可以算出来好的结果。
一个新的发现是壁面上会有奇异点,如下图,着色是
p_rgh
,这些奇异点都处于并行分区的各个界面上。更奇怪的是,上壁面并没有这样的问题。Note:参考点放置在流场中间。
现在暂时放弃自适应加密了,直接全场加密算,压力就是正确的。
😂 -
单核跑算出来的结果非常好,图中所有测点的压强降到-140左右是因为气泡经过了参考点,参考点(此时在气泡内部)的压强为0,外部的水的压强就应该是-144。这个算例的参考点、测点的设置和我上传的那个不一样。
-
我搞了一个125000网格的小算例,也可以复现流场整体异常压力跳跃的问题,不用跑完,跑几分钟就可以用里面的gnuplot脚本画图看到流场的压力。
算例放上来了,李老师@李东岳 您有空帮我看看,非常感谢!
singleBubbleChannel_debug.zip -
学流体的小明 在 interFoam计算气泡槽道流时的压力问题 中说:
我按照您说的方法修改了代码,没有什么大的改变
甚至导致结果更差了。
不进行这样的修改的话,相同的设置下,相同测点的压力如下图,可以看到紫色线有明显的台阶上升,这是气泡经过了这个测点,测到了气泡内部的压强。
-
我按照您说的方法修改了代码,没有什么大的改变。虽然这行代码确实不太对。
还有就是不太明白setRefrence()
这个函数template<class Type> void Foam::fvMatrix<Type>::setReference ( const label celli, const Type& value, const bool forceReference ) { if ((forceReference || psi_.needReference()) && celli >= 0) { source()[celli] += diag()[celli]*value; diag()[celli] += diag()[celli]; } }
把方程中,RefCell位置的源项加上一个值,然后系数矩阵对角线上RefCell位置翻倍?
还发现一点是自适应加密的间隔可以调小,改成1。每次都运行自适应这个模块,但不是每一步都会导致网格变化,所以调小一些也没事。可能需要
nBufferLayers
调大到2。
自适应加密间隔调小之后,流场也算得比较好了。自适应加密间隔再从1调整到5,压力马上就出问题了。下面图中0.09s之后就是自适应加密间隔调大的后果。
-
@李东岳 东岳老师,最近试了很多参数,发现了点规律。
- 应该是自适应网格的问题。如果使用自适应网格而不显式地指定correctPhi,程序会自动调用correctPhi。并且自适应网格必须correctPhi,否则会导致alpha.water出问题,全场都加密。
bool correctPhi ( pimple.dict().getOrDefault("correctPhi", mesh.dynamic()) );
- 不使用自适应网格,但是correctPhi yes,流场的压力结果也正常。
- 单核跑,即使是自适应网格,自动调用correctPhi,流场也计算正常。
应该是 并行 + 自适应网格 的问题?是不是要在dynamicMeshDict中进行correctFluxes操作呀?但我看Tutorial里面都没correctFluxes。
-
学流体的小明 在 压力基求解器在OpenFOAM中的植入问题 中说:
CaseB中,同样的方法,压力p基本都是10000,或者有是9999.99,这和我限制了writePrecision 6有关系。
所以输出流场的时候,最好以binary二进制的格式输出。
后处理的文件,由于要使用其它软件读取,就用ASCII。 -
我也刚好遇到了这个问题,我的理解是在程序计算的过程中,数值精度已经足够高了,可以保证计算正确。但是在将计算结果写到文件的时候,由于输出有效数字的位数,经常会丢失信息。如果断点续算,读文件输入程序的信息就可能产生问题。
举个例子,我计算一个槽道流,压力参考点生效。
区分一下两个算例:CaseA,pRef=0;CaseB,pRef=10000。
这两个Case都可以把壁湍流的速度剖面算得很好,两个Case的结果完全一致,曲线可以重合。但是两个Case输出的压力结果相差很大。CaseA中,paraview查看某个瞬时的流场数据或者surfaces取样输出某个平面的结果,压力p都是
0.XXXX
小数点后很多位。
CaseB中,同样的方法,压力p基本都是10000,或者有是9999.99,这和我限制了
writePrecision 6
有关系。
-
有点意思,所以结论是“时间步长不能调太小?”
-
李东岳 在 interFoam计算气泡槽道流时的压力问题 中说:
为什么不同的网格点的压力是一样的?A点的值可能是100+0.1,B点的值可能是100+0.2,就是修正的量太大了,把它们本身的空间起伏都淹没了
李东岳 在 interFoam计算气泡槽道流时的压力问题 中说:
单气泡算例,p符合预期么(跟表面张力的关系)符合,多气泡算例的表面张力关系也是符合的。下面这个截面的图,
p_water = -289 Pa, p_gas = -165 Pa
,虽然气泡半径0.001 m,应该是差144 Pa,但考虑到变形,算出来的结果差124 Pa也还是符合预期的。
下面这张图,中间那个红点是压力参考点,
p = 0 Pa
,两侧的橘色的斑块是气泡。就可以看到在参考点之外,流场的液体的压力迅速就降低到-289 Pa
了。这算是没收敛么
-
@coolhhh 但这样会导致很多问题,入口速度设置好之后,流场还要向下游发展,如何保证湍流没有被耗散掉——这个时候应该不能用
meanVelocityForce
了。
我是要得到一个稳定的气泡输运过程,如果按照您的思路,气泡在上游初始化,之后得和湍流一起运动发展一段时间,到达计算域中游下游时候,我可以采集我需要的信息,这会造成我的计算域非常长。而且也很难保证最终是一个遍历态。
还是先想想从求解压力的方面入手吧😂 -
@coolhhh 感谢你的回复和讨论。
- 关于质量通量平衡这一点,确实就像你说的,我用周期边界应该不受这个影响。
coolhhh 在 interFoam计算气泡槽道流时的压力问题 中说:
需要设置在建筑物位置上方一定距离且要在建筑物影响范围外,此时建筑物位置附近就不受非物理压力脉动影响。参考点的位置确实会极大地影响到它附近的压力和流场周围的压力,比如下面这张图片,参考点附近的
pPrime2Mean
是很小的。
Fig.t=0.45s
,展向中间截面,pPrime2Mean
coolhhh 在 interFoam计算气泡槽道流时的压力问题 中说:
流场中是否有气泡完全不经过的地方?压力参考点的设置的不受气泡影响的位置处有的,从当前算出来的结果来看,气泡不会附着到壁面上,也很少靠近壁面,将参考点设置在壁面附近是一种选择,但这会带来另一个问题
——由于我的算例上下都是壁面,所以参考点设置在靠近壁面的地方会导致两个壁面上的结果不一致。下面图片是单相槽道流的壁面pPrime2Mean结果,参考点设置是pRefCell=0
,pRefValue=0
,可以看到下壁面的四个角附近的pPrime2Mean
是蓝色的很小的值。由于我是要对壁面上的压力进行分析,所以发现这个问题之后,我重新算了参考点设置在计算域中心的的算例。
coolhhh 在 interFoam计算气泡槽道流时的压力问题 中说:
尝试入口如果只加平均速度剖面,出口设置fixedValue=0,这个时候流场肯定没有压力脉动。再同样测试下气泡槽道中的同样位置的监测点压力,看看是否也是压力脉动过大。如果是的话,那就说明气泡产生的影响不可忽略,就类似于流场中放了个建筑绕流一样,建筑周边就有大的压力波动。感觉气泡产生的影响确实是不可忽略的。
coolhhh 在 interFoam计算气泡槽道流时的压力问题 中说:
原来的单向流槽道用meanVelocityForce驱动,因为边界没有指定压力值,应该也调用了压力参考点作用?查看算例是指定了pRefCell 1001; pRefValue 0;。整个流场应该是压力梯度驱动流动,但查看of自带的Channel395算例结果的压力均值,好像没有明显看到沿着X向的压力梯度?我计算的单相槽道流也是基于channel395的各种设置的,是调用了压力参考点,但是单相槽道流中,每次修正的值都非常小,压力数据是可以用的。下面的图是壁面上沿流向的一条线(
y=0,z=0.01,x
变化)的时间历程,pPrime=p-pMean
都是直接从OpenFOAM的后处理结果处理出来的,这张图就很nice。
因为是槽道流,设置的压力梯度是和壁面剪切力平衡了的,确实不会看到沿X向的压力梯度。meanVelocityForce
相当于加了源项,和重力的表现不一样。 -
fvSolution
/*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: v2012 | | \\ / A nd | Website: www.openfoam.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; location "system"; object fvSolution; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // solvers { "alpha.water.*" { nAlphaCorr 2; nAlphaSubCycles 3; cAlpha 1.0; MULESCorr yes; nLimiterIter 5; solver smoothSolver; smoother symGaussSeidel; tolerance 1e-8; relTol 0; } p_rgh { solver GAMG; tolerance 1e-07; relTol 0.01; smoother GaussSeidel; } "pcorr.*" { solver GAMG; smoother GaussSeidel; tolerance 1e-8; relTol 0; } ".*(rho|rhoFinal)" { solver diagonal; } p_rghFinal { $p_rgh; tolerance 1e-07; relTol 0; } U { solver smoothSolver; smoother GaussSeidel; tolerance 1e-08; relTol 0; nSweeps 1; } "(T|k|B|nuTilda).*" { solver PBiCGStab; preconditioner DILU; tolerance 1e-8; relTol 0; } } PIMPLE { momentumPredictor no; nOuterCorrectors 1; nCorrectors 3; nNonOrthogonalCorrectors 0; correctPhi yes; pRefPoint (0.0201 0.0101 0.0101); pRefValue 0; } // ************************************************************************* //
fvScheme
/*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: v2012 | | \\ / A nd | Website: www.openfoam.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; location "system"; object fvSchemes; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // ddtSchemes { default Euler; } gradSchemes { default Gauss linear; } divSchemes { default none; div(phi,alpha) Gauss vanLeer; div(phirb,alpha) Gauss linear; div(rhoPhi,U) Gauss limitedLinearV 1; div(phi,p) Gauss linearUpwind grad(p);//limitedLinear 1; div(phi,k) Gauss linear; div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear; } laplacianSchemes { default Gauss linear corrected; } interpolationSchemes { default linear; } snGradSchemes { default corrected; } // ************************************************************************* //
dynamicMeshDict
/*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: v2012 | | \\ / A nd | Website: www.openfoam.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; location "constant"; object dynamicMeshDict; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dynamicFvMesh dynamicRefineFvMesh; // How often to refine refineInterval 5; // Field to be refinement on field alpha.water; // Refine field inbetween lower..upper lowerRefineLevel 0.001; upperRefineLevel 0.999; // If value < unrefineLevel unrefine unrefineLevel 10; // Have slower than 2:1 refinement nBufferLayers 1; // Refine cells only up to maxRefinement levels maxRefinement 1; // Stop refinement if maxCells reached maxCells 36000000; // Flux field and corresponding velocity field. Fluxes on changed // faces get recalculated by interpolating the velocity. Use 'none' // on surfaceScalarFields that do not need to be reinterpolated. correctFluxes ( (phi none) (nHatf none) (rhoPhi none) (alphaPhi0.water none) (ghf none) (alphaPhiUn none) ); // Write the refinement level as a volScalarField dumpLevel true; // ************************************************************************* //
-
各位老师同学好,我现在在算一个槽道,里面有很多气泡。基于
interFoam
自己添加了一些需要的功能,所以水和气都是不可压的。现在的问题是流场的压力结果总是不好。
算例是这样的:x方向为流动方向,周期边界条件;y方向是壁面,上下都是壁面;z方向是展向,也是周期边界条件。流场通过
meanVelocityForce
驱动,流场的平均速度会基本保持在一个量上。没有重力。
之前先使用pimpleFoam
算过一个单相的湍流槽道流,达到稳定后,把湍流槽道流的速度场mapFields
到现在这个算例当中,然后用setFields
把气泡初始化,具体来说是alpha.water
、p_rgh
、p
都初始化。
使用了自适应网格加密,更好地捕捉界面,运行起来网格大约在7e6 ~ 8e6
。
由于压力边界条件没有在任何一个边界上给定一个压力值,所以interFoam
求解器(其它pimple
算法求解器同样也)会启用“压力参考点”这一套东西。就是指定一个pRefCell
,给一个pRefValue
,在每一次压力计算结束之后,流场会整体地加一个数值,使得这个pRefCell
的压力总等于pRefValue
。具体代码如下:p == p_rgh + rho*gh; if (p_rgh.needReference()) { p += dimensionedScalar ( "p", p.dimensions(), pRefValue - getRefCellValue(p, pRefCell) ); p_rgh = p - rho*gh; }
当槽道中仅有一个气泡时,每次修正的值并不大,还可以接受,流场中某点的压力随时间的变化(由于没有重力
p=p_rgh
)如下图,是可以接受的,由于有气泡经过,该点p会有几次峰。
但是在我计算下面这个很多气泡的流场时
流场中一些监测点的压力会变成下面这样,可以看到很多点的压力基本是一样的,而且振荡很厉害,完全不连续。这是因为每次修正的值,也就是上面公式中的
pRefValue - getRefCellValue(p, pRefCell)
很大,淹没了该点本身有意义的值。从上面的压力场中我无法进行进一步的分析,平均值和脉动平方均值都感觉没什么意义——因为修正操作造成的波动太大了。
自己分析有以下几种解决方案,也尝试了一下
- 是因为压力计算没有收敛吗?使用合适的
fvSolution
和fvScheme
可以解决?fvSolution
和fvScheme
我放在下一层。 - 关掉压力修正操作,流场的压力可能会一直增大或者减小,最后可能是
1e7Pa
这样很大的值,但是压力梯度本身还是有用的。 - 关键原因是求解方程,求解出来
pRefCell
的值可能就是500
,然后流场就会整体减去500
让它等于0
。
while (pimple.correctNonOrthogonal()) { fvScalarMatrix p_rghEqn ( fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA) ); p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell)); p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter()))); if (pimple.finalNonOrthogonalIter()) { phi = phiHbyA - p_rghEqn.flux(); p_rgh.relax(); U = HbyA + rAU()*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf); U.correctBoundaryConditions(); fvOptions.correct(U); } }
各位大佬有遇到过相似的问题吗?求指点🙏🙏🙏
- 是因为压力计算没有收敛吗?使用合适的
-
@SSSSK
我$ \Delta x^+=\Delta z^+=10 $,壁面法向$ \Delta y^+=2 - 10 $你可以看看这篇文章
https://www.sciencedirect.com/science/article/pii/S0021999117304059?via%3Dihub -
@SSSSK Uf是平均速度啊,都知道Uf(y)了,直接算剪切力 $ \mu\partial u/\partial y \big |_{y=0} $ 啊
-
@SSSSK 应该用空间平均后的数据,最后就是一个只随y坐标变化的函数,比如u(y)。建议把这个帖子好好翻一翻,好多东西应该都讨论过。
学流体的小明 在 LES直流槽道边界层模拟,如何得到正则化速度u+以及正则化坐标y+? 中说:
是底面平均之后的,postChannel工具会把槽道按照xz线(x和z坐标固定,y坐标变化)求平均
-
@SSSSK
我用给定槽道平均速度的方法算摩擦雷诺数550的槽道,根据这个摩擦雷诺数的变化曲线,感觉得到10s以后吧,你也可以画这样一条曲线看看。
-
simpleFoam为OpenFOAM中一个基于有限体积法的,用于求解稳态不可压缩牛顿流体N-S方程的求解器。
http://dyfluid.com/simpleFoam.html那就是提前算到稳态了呗,所以结束了。
interFoam吐核
openfoam9安装失败
如何在流场中加入体积力场
openfoam9安装失败
如何在流场中加入体积力场
关于编译动态库后在求解器中使用报错
关于编译动态库后在求解器中使用报错
在OpenFOAM中如何判断一个单元是否为边界单元?
overinterdymfoam的一些疑问
overinterdymfoam的一些疑问
openform的瞬时温度场
使用自适应网格细化并行计算结果在进行reconstructPar进行整合的时候会出现问题
openfoam对于气泡进行自适应网格细化出现棱角
openfoam如何中途调整更改网格?
openfoam如何中途调整更改网格?
关于InterFOAM曲率计算的问题
VOF添加斥力模型
VOF添加斥力模型
VOF添加斥力模型
VOF添加斥力模型
VOF添加斥力模型
VOF添加斥力模型
compressibleInterFoam与CalculiX求解流固耦合问题
VOF添加斥力模型
VOF添加斥力模型
VOF添加斥力模型
VOF添加斥力模型
VOF添加斥力模型
VOF添加斥力模型
pimplefoam在瞬态模拟过程中库朗数突然增大然后报错停止
interFoam计算气泡槽道流时的压力问题
interFoam计算气泡槽道流时的压力问题
interFoam计算气泡槽道流时的压力问题
interFoam计算气泡槽道流时的压力问题
interFoam计算气泡槽道流时的压力问题
interFoam计算气泡槽道流时的压力问题
压力基求解器在OpenFOAM中的植入问题
压力基求解器在OpenFOAM中的植入问题
时间步长对连续性方程的影响
interFoam计算气泡槽道流时的压力问题
interFoam计算气泡槽道流时的压力问题
interFoam计算气泡槽道流时的压力问题
interFoam计算气泡槽道流时的压力问题
interFoam计算气泡槽道流时的压力问题
interFoam计算气泡槽道流时的压力问题
LES直流槽道边界层模拟,如何得到正则化速度u+以及正则化坐标y+?
LES直流槽道边界层模拟,如何得到正则化速度u+以及正则化坐标y+?
LES直流槽道边界层模拟,如何得到正则化速度u+以及正则化坐标y+?
q-DNS计算槽道流遇到了一些问题,求大神们指点
请教一个openfoam边界问题:fixedprofile的使用