用twoPhaseEulerFoam求解水砂输送的问题

Zhy2022

各位前辈好，我想在一个长方体几何模型中从左往右输送水（85%）砂（15%）混合物，观察砂的堆积形态（如下图所示）。

我的速度边界为：

object      U.particles;
internalField   uniform (0 0 0);

boundaryField
{
    INLET
    {
        type               fixedValue;
        value              uniform (0.5 0 0);
    }

    OUTLET
    {
        type               inletOutlet;
        inletValue         uniform (0 0 0);
        value              uniform (0 0 0);
    }

    WALL
    {
        type               zeroGradient;
    }

    frontAndBackPlanes
    {
        type               empty;
    }
}

object      U.water;
internalField   uniform (0.5 0 0);

boundaryField
{
    INLET
    {
        type               fixedValue;
        value              uniform (0.5 0 0);
        /*type               interstitialInletVelocity;
        inletVelocity      uniform (0.8 0 0);
        alpha              alpha.water;
        value              $internalField;*/
    }

    OUTLET
    {
         type               inletOutlet;
        inletValue         uniform (0 0 0);
        value              uniform (0 0 0);
    }

    WALL
    {
        type               noSlip;
    }

压力边界为：

internalField       uniform 0;

boundaryField
{
    
    INLET
    {
        type            fixedFluxPressure;
        value           $internalField;
    }
    
    OUTLET
    {
         type            fixedValue;
         value           uniform 0;
    }

    WALL
    {
        type            zeroGradient;
    }

    frontAndBackPlanes
    {
        type            empty;
    }
}

但不知道为何第一个迭代步中，e.water的迭代次数就超了，导致压力的迭代次数也很大，算了几个时间步后就发散了：

实时证明第一个时间步算出来的速度和压力结果就不对，还请各位前辈告知e.water是什么量，这样设置边界条件对不对。谢谢！@东岳

李东岳

把你的能量方程关掉。在phaseProperties里面选pureIsothermalPhase，具体你找找，手头没有openfoam

Zhy2022

@东岳谢谢东岳老师，也就是e.water是跟传热有关的参数？我在of8里面又跑了一下（之前是of5）能跑动，只是压力迭代次数任然比较高（70次）左右，初始残差0.7左右，也老是降不下来

不知道压力这样设置有没有问题

李东岳

最后一个图说明不了什么问题可以
对 e跟温度有关

Zhy2022

@东岳好的！感谢东岳老师的答疑解惑！计算过程中除了压力初始残差下不去之外，maxCo也是随计算时间逐渐增大的，照这个进度下去感觉算久了就会发散。后续还会继续跟进这个例子！

Zhy2022

@东岳东岳老师，我感觉压力边界设置得还是有些问题，计算结果表明prgh在整个过程中都没什么变化，p的变化也很奇怪：
（prgh）
（p）

李东岳

两个都很合理啊。p_rgh不应该变化太多。p正好是存在泥沙的时候，在下面压力比较大。应为重力的因素。

Zhy2022

@东岳那么东岳老师，为何p_rgh是上方大，下方几乎为0呢，而且p的量级也没有p_rgh大

Zhy2022

@Zhy2022 明白了。。p_rgh只是一个参考压力系

Jacob Z

出口流速zerogradient会不会不一样？

Zhy2022

@Jacob-Z 回流很明显，容易发散

Zhy2022

@东岳东岳老师您好，再向您请教一个问题，为什么我的两相速度能够计算出来：
（颗粒速度）
（液体速度）
但是残差监测里面显示N/A呢？

李东岳

你的residual里面应该没有监控U.particle，默认是U，你需要指定fields (U.particle)

Zhy2022

@东岳非常感谢东岳老师这么晚还在回复问题！！！希望您也能休息好！！！目前有三个问题：
（1）我在residuals字典里面写明了U.particles（U.particle监控不了），才出现的N/A这个问题；
（2）模拟结果表明砂堤没有沉积（注入的固相百分数就为0.15，而且停止注入固相后，固相整体被往前推动），不知道是不是颗粒动理论模型设置的问题

（3）液相、固相注入速度都是（0.5 0 0），但计算过程中底部达到4m/s+，这是否合理（液相的internalFiled也给的0.5，感觉液-固混合液比纯液跑得还快了）：

（边界给的固定值（0 0 0），插值看不到结果）

希望您有空的时候能帮忙解答，不胜感谢！

李东岳

液相、固相注入速度都是（0.5 0 0），但计算过程中底部达到4m/s+，

这是用层流模拟的么？好像速度确实有点大。不过也不好判断。没准是重力的作用加速了？

Zhy2022

@东岳 RAS，但是重力也不该水平方向这么大的速度哇，主要是没堆积

Zhy2022

@东岳东岳老师您好，还想请问一下：我按照of自带算例将固相模型改为从kineticTheory改为phasePressure，以实现层流模拟，报错：
2020-10-28 23-28-43屏幕截图.png
我检查了初始场文件发现没有问题，不知这是怎么回事。。还请您看一看这个gs0.particles是什么东西。。

李东岳

把你fvScheme发一下

Zhy2022

@东岳老师，您看：

ddtSchemes
{
    default     Euler;
}

gradSchemes
{
    default     Gauss linear;
}

divSchemes
{
    default                         none;

    "div\(phi,alpha.*\)"            Gauss vanLeer;
    "div\(phir,alpha.*\)"           Gauss vanLeer;

    "div\(alphaRhoPhi.*,U.*\)"      Gauss limitedLinearV 1;
    "div\(phi.*,U.*\)"              Gauss limitedLinearV 1;

    "div\(alphaRhoPhi.*,(h|e).*\)"  Gauss limitedLinear 1;
    "div\(alphaRhoPhi.*,K.*\)"      Gauss limitedLinear 1;
    "div\(alphaPhi.*,p\)"           Gauss limitedLinear 1;

    div(alphaRhoPhi.particles,Theta.particles) Gauss limitedLinear 1;

    "div\(alphaRhoPhi.*,(k|epsilon).*\)"  Gauss limitedLinear 1;

    div((((alpha.water*thermo:rho.water)*nuEff.water)*dev2(T(grad(U.water))))) Gauss linear;

    div((((thermo:rho.particles*nut.particles)*dev2(T(grad(U.particles))))+(((thermo:rho.particles*lambda.particles)*div(phi.particles))*I)))  Gauss linear;

    div((((alpha.particles*thermo:rho.particles)*nuEff.particles)*dev2(T(grad(U.particles))))) Gauss linear;//改层流后多出
}

laplacianSchemes
{
    default     Gauss linear uncorrected;
    bounded     Gauss linear uncorrected;
}

interpolationSchemes
{
    default     linear;
}

Zhy2022

@东岳老师，问题解决了：
固相的速度边界和湍流模型必须配套使用才行！

Zhy2022

@东岳东岳老师，有个问题不太想得通，不知道您对此看法如何：
采用TFM模拟固-液两相流，固相密度设置为常数（fluidisedBed案例中也是如此），视为不可压缩，但固相湍流模型（kineticTheory或phasePressure）都是属于可压缩模型，它们是否矛盾？

Zhy2022

@李东岳东岳老师，在kineticTheoryModel.C中有这样一段代码：

Foam::tmp<Foam::volScalarField>
Foam::RASModels::kineticTheoryModel::k() const
{
    NotImplemented;
    return nut_;
}


Foam::tmp<Foam::volScalarField>
Foam::RASModels::kineticTheoryModel::epsilon() const
{
    NotImplemented;
    return nut_;
}

但是感觉固相湍流开关打开的话，其湍流黏度完全是按照如下公式计算的：
（其中a为固相，b为液相）
也就是固相被看作连续相后不存在真实的k和epsilon哇？

Zhy2022

同样工况下，fluent模拟结果和OF模拟结果存在差异：
fluent的砂堤快速堆积，而OF的砂堤会先跑满整个空间，再堆积回缩，堆积也达不到最大堆积程度0.63，
而且OF中固相Y方向速度有局部异常向下的特点，不知道是不是固相计算有什么问题，不知东岳老师@李东岳或者其他前辈能否给一些建议？

李东岳

你这个是否有参考的sci或者参考数据可以对比？有的话可以把openfoam算例清理干净上传给我我给你调一下，如果不涉密的话

Zhy2022

@李东岳谢谢东岳老师！暂无可参考的sci，目前仅能通过实验结果进行粗略对比。
最大问题还是固相分数回流、沉降达不到最大堆积分数0.62和固相Y速度不太对劲。
case如何上传给您？实验视频有10个G不太方便发送，我以截图的方式一并发给您

Zhy2022

@李东岳东岳老师，有一篇Mfix的模拟结果，出口设置和我的模型略有不同，影响不大，一并发给您

veen

@Zhy2022
之前在openfoam里边算过空气中沙物质沉降，算出来最大体积堆积分数和fluent是可以对上的，会不会是alphaMax 这个参数在设置的问题呢？
替代文字

Zhy2022

@veen 谢谢前辈，我打的最大相分数给的0.62，但是堆积只有0.59左右，这也是个很奇怪的问题

Zhy2022

@veen 前辈，您的设置能否分享一下（不用网格）？

李东岳

@Zhy2022 不好意思算例太大了，你把算例清理好只保留必要文件发给我就好，不需要结果文件

veen

@Zhy2022 刚刚看了我的alphamax也没改，链接在这：https://github.com/Veenxz/Veen_OpenFOAM_Cases/tree/main/2D_Plane_EE_Multiphase

veen

@Zhy2022 你直接blockmesh改一下，然后改一下流体参数应该就好了，曳力函数也需要改一下。

Zhy2022

@李东岳好的东岳老师，其实里面没有结果文件，只是我实验做的一点结果和文章我再给您发一个只有算例的

Zhy2022

@veen 谢谢前辈！目前打算先用gidaspow算算后面会修正曳力模型

李东岳

我大略看了下你这个2D网格但是前后不是empty？why

Zhy2022

@李东岳李老师，并非2D网格，宽度方向划分了5个网格（很窄）；前后目前看作壁面，后期会看作具有一定渗透性的出口

Zhy2022

@veen 前辈，参考您的设置，我将碰撞回复系数和镜面系数都进行了调整，但是不知道这两个参数分别代表什么，之前也没考虑过。
调整后的计算结果跟预想一样，沉降更迅速，但仍会出现颗粒堆积达不到alphaMax、先铺开再回缩（到一定高度后再往后堆积）以及颗粒Y方向局部有向下速度的情况。

（X方向）
（Y方向）

不知道是不是压力出口设置的问题（prgh给的固定值0）；另外，我看您相分数设置了回流，之前我设置回流的时候会报错，不知道这是不是原因之一。

@李东岳不知东岳老师您有没有其他看法。

veen

@Zhy2022 我那个也就是个Test case，从空气的计算来看没什么问题，你可以先计算一个纯气流的稳定场，然后map到多相流里边再开始算沉积，随着网格分辨率的提高，结果应该会越来越契合实际。
因为没看到你的case，我觉得可能时间步、网格、求解方法可以再调整调整。最大相分数不对我也不清楚为什么，你可以先设置为0.72看看堆积的相分数会不会再增加，更重要的还是看你想说明什么问题再去调整。
另外，水沙输运的求解也可以试试sedFoam。

Zhy2022

@veen 感谢指点！目前的网格已经是加密后的了，对入口射流细节的捕捉更加到位。时间步采用自适应，也能最大程度保证计算准确性。
现在就想明确两个问题：为什么最大堆积分数不像fluent那样，（1）在砂堤沉降后就达到alphaMax（明明已经没有速度）；（2）为什么砂堤铺开后会回缩一段距离后再往前堆积。
我会根据后续模拟结果再做反馈！

veen

@Zhy2022 之前你圈起来的两个参数应该是固相之间的相互作用力或者双流体之间的相互作用力系数，你可以查一下文献。

你的入口边界是均一的固相体积分数？因为你的固相体积分数比较大，不知道你的实验是什么样的设定，但是实际情况中水流携沙也是随着离壁面距离增加而减小的，我的建议可以按实际情况去设定。

在沙堤沉降就达到alphaMax这个是没有问题的吧，或者我没太理解你的意思。
我建议你可以在模拟中用一下网格自适应，提高一下局部的分辨率。

关于沙堤堆积后再推进和水沙运动理论有关，在水流携沙运动中，开始的时候水流中沙物质分数较大，在非饱和路径中以堆积为主，随着沙物质的堆积，计算域内非堆积搬运出现概率逐渐增加，沙丘开始逐步运动，你可以求一下沙丘的几何中心，看一看沙丘的运动情况，最开始应该是随着堆积几何中心移动最快，然后出现一段比较稳定的区域，最后如果计算区域和时间足够长，它会以比较衡定的速度运动。（以上是我自己的理解，也有可能有些问题）

Zhy2022

@veen 多谢前辈的指点！
入口给的均一固相分数（根据实际工况设置）；
砂堤沉降就达到alphaMax是fluent的模拟结果，我的这个模拟工况就是没能达到，所以感觉有问题；
在我的工况下，砂堤运动规律大概是清楚的，因此推测固相分数在模拟中回缩可能是有问题的
不太明白砂堤几何中心该如何获得。
下面附上实验结果和模拟结果（没法传动图）：

Jacob Z

随便一问，是西南石油大学的设备吗？

Zhy2022

@Jacob-Z 犀利，欢迎指导交流

veen

@Zhy2022 加油！openfoam里边双流体模型我没有很深入去研究过，几何中心的求取直接根据定义写代码也行https://en.wikipedia.org/wiki/Centroid，用opencv做图像处理之后提取也行。
你算例中的问题可能与设置有关，可以再调试调试算理。

Zhy2022

@veen 好的，谢谢您！

veen

有没有什么进展或者文章见刊啊？

Zhy2022

@veen 还在整理研究结果，但是也有了一点新认识：
1）沉降效果与曳力模型和颗粒粒径有关；
2）砂堤回缩与边界条件和固相应力模型（猜测，用JJ比JJShaffer要好，且正在植入其他模型）有关；
3）达不到alphaMax我猜想可能也和固相应力模型有关。
争取早日把成果研究整理分享出来

veen

@Zhy2022