CFD中文网

出口A和入口B尺寸一样，可以做成一一对应的网格

没有太多用openfoam做实际应用的经验，所以想请教一下各位。

有一个如图的模拟，计算域只有高度不等的A腔和B腔，因为不关注管道内的流动，所以蓝色线表示的连通管道就省略了。现在想把出口A的压强（减去恒定的高度压力差）和速度（旋转90度）直接映射到入口B上，并且每个时间步都执行这一映射。请问有相应的边界条件和tutorial可以参考吗？

或者这种情况，还有什么简化的实现方法嘛？

@李东岳 谢谢李老师。

这应该是TJunctionFan这个tutorial吧，这个确实能得到局部的速度增加，但是如果管道足够长的话，管道的入口和出口处仍然是相同的速度和流量。

回复您之前的问题，是inlet，outlet这两个地方。我想实现的结果就是在outlet的速度比inlet速度高，这确实是一个不符合物理的问题。。流量不守恒。没有外力的情况下是不可能实现的。硬要做的话，好像是可以通过添加momentum source实现，但是确实是一个非物理的问题，打扰了。。

对的，其实就是中间有加速装置，简化成流量越来越高，所以想着fan是一种实现方式，不知道还有没有其他方式

如图所示，一个长圆管中，从左到右的流速越来越高。请问这样的流速（或者流量）随距离增加的管流如何实现？最好是在RANS模拟中。

我目前想到的办法是使用createBaffles在圆管内部每隔一小段距离建立若干个patch面，每个内部patch都使用fan边界条件（如下图所示，压强场非均匀增加，其他场为cyclic）。结果得到的管中压强场虽然是非线性降低，但是管中左右两边的速度仍然相同，就跟普通的管流一样，请问问题出在哪里？

或者各位前辈有没有其他更好的办法？

@bestucan

@dingcy 在 颗粒碰撞聚合 中说：

(1)MPPICFoam能不能算连续介质中固体颗粒的碰撞聚合问题；

在您回帖的链接中好像没有看到关于这个问题的解答，请问颗粒的碰撞聚合问题在MPPICFoam中能实现嘛？

好像欧拉颗粒的话可以通过population Balance Model实现，不知道拉格朗日颗粒的碰撞聚合问题在OpenFOAM中怎么实现？

比如在injectionModel里的颗粒速度，UMag，想设置为随模拟时间变化的函数。请问应该怎么调取这个时间变量runTime.value()？是否应该在哪里include关于时间的H文件？

@李东岳 谢谢李老师解答。

那是不是可以理解为，在比如下图中（打开颗粒对流体的影响），颗粒的中心点在网格1，但是颗粒的实际大小有如图显示那么大。现在就只有网格1受到颗粒的影响，而2,3,4号网格都完全感受不到颗粒的存在？

不知道这对总体结果会造成多大的影响？

ErgunWenYuDragForce.C中，计算颗粒曳力的代码如下：

     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()))
         );
     }

这里第一行的alphac(alphac_[p.cell()])，是插值计算网格的alphac值(cellPoint方式)，还是只读取网格中心值(cell方式)？因为没有查到关于这个[p.cell()]的解释，所以想请教一下论坛里的前辈

在WallSpringSliderDashpot.C中，有这样几句：

    scalar etaN = alpha_*sqrt(p.mass()*kN)*pow025(normalOverlapMag);
    vector fN_PW =
        rHat_PW
       *(kN*pow(normalOverlapMag, b_) - etaN*(U_PW & rHat_PW));
    ...

    scalar etaT = etaN;
    ...
    (if条件：颗粒不滑移的条件下)
    fT_PW = - kT*tangentialOverlap_PW - etaT*USlip_PW;

这里的颗粒碰撞中的切向阻尼系数，etaN，等于法向值，然后分别应用到求碰撞力的切向和法向分量，fN, fT，中。是不是意味着，我在颗粒配置文件中通过alpha设置的回弹系数值e, 就是同时设定了切向和法向的回弹系数，eT, eN, 都等于e？

请问这个理解对不对？

@李东岳 谢谢李老师，这个确实是我想找的功能。

不好意思才看到您的回复，不知道为啥论坛里的消息有时候没有提示，也没有邮件提示。。

比如用patchInjection，想让inlet patch在第1，3，5...秒内入射100个颗粒，在第2，4，6...秒内入射10000个颗粒。

请问可以通过在kinematicCloudProperties文件中调用time()变量实现吗？

@星星星星晴 谢谢解答，下面是calcVelocity的部分，但是因为在那些acp, ancp, bcp的存在，让我越来越混乱了

想问问这么理解对不对：
（1）这个kinematicParcel中的calcVelocity部分，主要是为了计算颗粒相带给流体相多少动量转移，用于计算流体相的速度？
（2）而如果只计算单个颗粒相本身的速度，其实就是可以通过F = Sp(U - Up) + Su这个合力能计算的，Sp(U - Up)就是单个颗粒所受的曳力，Su就是升力，重力等等的合力？

    // Calculate the integration coefficients
    const vector acp = (Fcp.Sp()*td.Uc() + Fcp.Su())/massEff;
    const vector ancp = (Fncp.Su() + Su)/massEff;
    const scalar bcp = Fcp.Sp()/massEff;

    // Integrate to find the new parcel velocity
    const vector deltaU = cloud.UIntegrator().delta(U_, dt, acp + ancp, bcp);
    const vector deltaUncp = ancp*dt;
    const vector deltaUcp = deltaU - deltaUncp;

    // Calculate the new velocity and the momentum transfer terms
    vector Unew = U_ + deltaU;

    dUTrans -= massEff*deltaUcp;

    Spu = dt*Fcp.Sp();

    // Apply correction to velocity and dUTrans for reduced-D cases
    const polyMesh& mesh = cloud.pMesh();
    meshTools::constrainDirection(mesh, mesh.solutionD(), Unew);
    meshTools::constrainDirection(mesh, mesh.solutionD(), dUTrans);

    return Unew;

在LiftForce的计算中，表达式为：

value.Su() = mass/p.rho()*td.rhoc()*Cl*((td.Uc() - p.U())^curlUc);

为什么这里会出现一个异或运算符“^”？看文献中的Lift的公式应该是

$F_{lift}=rho_f * V * Cl * (U_f-U_p) * |dUf / dy|$

1. 难道这个"^"就是简单的理解为乘号？

2. 这个curlUc等于是速度梯度|dUf / dy|吗？

3.另外想问问，我的理解是在SaffmanMeiLift.C或者TomiyamaLift.C中的运算，都只是为了求参数Cl的值，然后不管用哪个模型，这个Cl最终都要在LiftForce.C计算，以 value.Su() 的形式添加给粒子，这个理解对不对?

@bestucan 感谢回复，这么一说确实是，就没法达到稳定状态。

不过如果按照您说的，用拉格朗日粒子，如果我设置粒子为passive scalar，用icoUncoupledKinematicParcelFoam对一个不稳定的场来求解粒子，这样粒子速度也会不断增加，最终导致发散吗？

现在我得到了一个不均匀分布的外力场，想查看流体在这个外力场持续的作用下的流向和流速变化情况。

我有尝试过将外力场转化为压力场p，并保持速度场U不变，然后进行模拟。但是因为这样的p和U并不是稳态，压力场p很快就变化掉了，不能持续。所以这个办法应该不可行。

想问问各位前辈，在fvOption里面能不能添加自定义分布的，并且持续的外力场？或者有其他什么思路可以实现？

@马乔 想借这里请问您一下，如果我只想输出颗粒的Uc的话（写入到每个时间步的/lagrangian/kinematicCloud中），应该怎么操作？

我代码只有入门水平，实在看不太明白您的一系列操作修改，望您有空能指教

@浪迹天大 谢谢您的回复，我还是想问问这个position()是不是所说的颗粒的位置值？但是好像参数的格式与我最终要传到UTurb里的格式不太一样，所以导致无法编译成功。。想问问应该怎么一层一层传导？

@李东岳 谢谢李老师关注，我在做的是颗粒流射入腔体内的一个模拟，确实stochastic模型本身是只适合于各向同性流动的，算出来的结果很不好，所以我现在想将UTurb的计算过程进行修改，所以还是针对颗粒相的。

之前有试过在UTurb的方向值dir中添加与流体速度Uc的相关性（修改如下），结果能稍微有些改进。现在想将与颗粒位置的相关性也添加进来，但是position和Uc不一样，没有在这个C文件中声明，所以不知道该怎么调用它？

        {
            tTurb = 0;

            const scalar sigma = sqrt(2*k/3.0);

            // Calculate a random direction dir distributed uniformly
            // in spherical coordinates

            const scalar theta = rnd.sample01<scalar>()*twoPi;
            const scalar u = 2*rnd.sample01<scalar>() - 1;

            const scalar a = sqrt(1 - sqr(u));

            // 修改部分开始
            const scalar Uc_x = Uc.component(0);
            const scalar Uc_y = Uc.component(1);
            const scalar Uc_z = Uc.component(2);
        
            scalar f_x = mag(sqrt(sqr(Uc_y)+sqr(Uc_z)))/mag(Uc);
            scalar f_y = mag(sqrt(sqr(Uc_x)+sqr(Uc_z)))/mag(Uc);
            scalar f_z = mag(sqrt(sqr(Uc_x)+sqr(Uc_y)))/mag(Uc);

            scalar dir_x = f_x*a*cos(theta); 
            scalar dir_y = f_y*a*sin(theta);
            scalar dir_z = f_z*u;
            const vector dir(dir_x, dir_y, dir_z);
            // 修改部分结束
            (源文件是只有一句：const vector dir(a*cos(theta), a*sin(theta), u);)

            UTurb = sigma*mag(rnd.GaussNormal<scalar>())*dir;
        }

@浪迹天大 请问您说的改接口是指怎么操作？

我看到插值读取Uc的代码里，确实有颗粒位置值position()带入，您提到的改接口的意思是把这个position()赋值给一个新变量嘛？

Uc_ = td.UInterp().interpolate(this->position(), tetIs);

现想通过修改StochasticDispersionRAS的C文件，使颗粒的随机扰动UTurb，和颗粒所处的位置（距离原点的位置）相关。请问在颗粒的C文件中应该怎样正确地调用颗粒位置的值呢？

我有试过“const vector & position”，和“position()”之类的，但是并不能成功编译。