最近涉及到一些化学反应，并不十分理解。在比如reactingTwoPhaseEulerFoam求解器的tutorial，bubbleColumnEvaporatingReacting中，展示了所有物质都在气相状态下的反应，CO^0.93 + H2O^0.24 = CO2^0.69 + H2^1。因为该反应被定义在reactions.gas文件中，反应物和产物的属性也在thermo.gas文件中，所有物质都为气态。
我想问问如果我的反应产物是有液态或者固态参与，比如A(气态) + B(气态) = C(气态) + D(液态)，应该如何设置这个reaction和thermo等文件。有没有相关的教程或者资料可以参考？
出口A和入口B尺寸一样，可以做成一一对应的网格
没有太多用openfoam做实际应用的经验，所以想请教一下各位。
有一个如图的模拟，计算域只有高度不等的A腔和B腔，因为不关注管道内的流动，所以蓝色线表示的连通管道就省略了。现在想把出口A的压强（减去恒定的高度压力差）和速度（旋转90度）直接映射到入口B上，并且每个时间步都执行这一映射。请问有相应的边界条件和tutorial可以参考吗？
或者这种情况，还有什么简化的实现方法嘛？
@李东岳 谢谢李老师。
这应该是TJunctionFan这个tutorial吧，这个确实能得到局部的速度增加，但是如果管道足够长的话，管道的入口和出口处仍然是相同的速度和流量。
回复您之前的问题，是inlet，outlet这两个地方。我想实现的结果就是在outlet的速度比inlet速度高，这确实是一个不符合物理的问题。。流量不守恒。没有外力的情况下是不可能实现的。硬要做的话，好像是可以通过添加momentum source实现，但是确实是一个非物理的问题，打扰了。。
对的，其实就是中间有加速装置，简化成流量越来越高，所以想着fan是一种实现方式，不知道还有没有其他方式
如图所示，一个长圆管中，从左到右的流速越来越高。请问这样的流速（或者流量）随距离增加的管流如何实现？最好是在RANS模拟中。
我目前想到的办法是使用createBaffles在圆管内部每隔一小段距离建立若干个patch面，每个内部patch都使用fan边界条件（如下图所示，压强场非均匀增加，其他场为cyclic）。结果得到的管中压强场虽然是非线性降低，但是管中左右两边的速度仍然相同，就跟普通的管流一样，请问问题出在哪里？
或者各位前辈有没有其他更好的办法？
比如在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|$
难道这个"^"就是简单的理解为乘号？
这个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);