关于MULES的疑问
-
Dear FOAMers:
大家好，本人最近在看MULES算法，有些疑问，希望大家能为我解惑，可能比较长，首先表示感谢
我使用的是interFOAM求解器，以OF6为例，为了简便，假设MULESCorr为false，在alphaEqn.H中的求解过程为：
for (int aCorr=0; aCorr<nAlphaCorr; aCorr++) { #include "alphaSuSp.H" //声明源项，全部为zeroField类型 surfaceScalarField phir(phic*mixture.nHatf()); //nHatf为自由表面的单位法向量 tmp<surfaceScalarField> talphaPhi1Un //组建alpha通量，采用时间项若采用Euler，则cnCoeff= 1.0 ( fvc::flux ( phiCN(), cnCoeff*alpha1 + (1.0 - cnCoeff)*alpha1.oldTime(), alphaScheme ) + fvc::flux //压缩项 ( -fvc::flux(-phir, alpha2, alpharScheme), alpha1, alpharScheme ) ); /** 省略一部分 */ else { alphaPhi10 = talphaPhi1Un; MULES::explicitSolve ( geometricOneField(), //此处不太理解，请看后文 alpha1, phiCN, alphaPhi10, Sp, (Su + divU*min(alpha1(), scalar(1)))(), oneField(), zeroField() ); } //#include "setAlphaMR.H" alpha2 = 1.0 - alpha1; mixture.correct(); }观察MULES::explicitSolve函数，其传入8个参数，查阅其代码，有：
void Foam::MULES::explicitSolve ( const RhoType& rho, //对应前面的geometricOneField() volScalarField& psi, //alpha1 const surfaceScalarField& phi, surfaceScalarField& phiPsi, const SpType& Sp, //显式计算的源项？ const SuType& Su, //隐式计算的源项？ const PsiMaxType& psiMax, const PsiMinType& psiMin ) { const fvMesh& mesh = psi.mesh(); psi.correctBoundaryConditions(); if (fv::localEulerDdt::enabled(mesh)) { const volScalarField& rDeltaT = fv::localEulerDdt::localRDeltaT(mesh); limit(rDeltaT, rho, psi, phi, phiPsi, Sp, Su, psiMax, psiMin, false); explicitSolve(rDeltaT, rho, psi, phiPsi, Sp, Su); } else { const scalar rDeltaT = 1.0/mesh.time().deltaTValue(); //时间步的倒数 limit(rDeltaT, rho, psi, phi, phiPsi, Sp, Su, psiMax, psiMin, false); //这个函数非常长，没仔细看不知道什么意思，但估计是保证alpha的结果在0-1之间？ explicitSolve(rDeltaT, rho, psi, phiPsi, Sp, Su); //最后更新值还是调用的这个六参数的函数 } }关于上述六参数的函数：
template<class RdeltaTType, class RhoType, class SpType, class SuType> void Foam::MULES::explicitSolve ( const RdeltaTType& rDeltaT, //时间步的倒数 const RhoType& rho, //这个地方为什么起了个密度的名字？ volScalarField& psi, const surfaceScalarField& phiPsi, const SpType& Sp, //两个源项 const SuType& Su ) { Info<< "MULES: Solving for " << psi.name() << endl; const fvMesh& mesh = psi.mesh(); scalarField& psiIf = psi; const scalarField& psi0 = psi.oldTime(); psiIf = 0.0; fvc::surfaceIntegrate(psiIf, phiPsi); //这里将通量除以了体积并幅值给psiIf //省略一小部分 else { psiIf = ( rho.oldTime().field()*psi0*rDeltaT + Su.field() - psiIf )/(rho.field()*rDeltaT - Sp.field()); } psi.correctBoundaryConditions(); }我试图把上述else中的代码变成公式：
此公式与离散相方程后的结果类似，只是多出了rho，查阅第一段代码可以发现，实际应用的时候传入的参数为geometricOneField()，即为1.0，因此这其实就是一个系数，但是我不明白这样定义的意义在哪里？是有其他的公式会用到这个参数吗？为什么起了密度这个名字？这个地方一开始把我搞晕了好久。
此外，在interFOAM中，Su和Sp均为0，我查看代码发现其一个表示显式源项，另一个表示隐式源项，如果我想在相方程中加入显式源项，请问该如何添加呢？直接调用函数的时候把我的源项表示的场传进去（给Su）就行了吗？
还有，上述代码中的注释都是我自己的理解，发帖时写上去的，请大神们看看有没有哪里理解错误的，请批评指正，谢谢！
不知道我说明白了没，希望大家为我解惑，谢谢！当然，如果有其他关于MULES的问题也可以在这里一起讨论！
-
因此这其实就是一个系数，但是我不明白这样定义的意义在哪里？是有其他的公式会用到这个参数吗？为什么起了密度这个名字？这个地方一开始把我搞晕了好久。
在某些情况下需要传入密度，这里只不过假定不可压缩，为
geometricOneField()Su和Sp均为0，我查看代码发现其一个表示显式源项，另一个表示隐式源项，如果我想在相方程中加入显式源项，请问该如何添加呢？直接调用函数的时候把我的源项表示的场传进去（给Su）就行了吗？
理论上是可以的，数值稳定性问题需要进一步研究
你的理解是正确的。图中公式也是正确的。另外，
limit函数是用来求$\lambda$的值，可以参考Boris and book 1973年关于Flux correct transport的文章，对应文章中的$\lambda$ -
@东岳 感谢东岳老师的回复，关于源项的问题，之前我是这么解决的：
if (MULESCorr) 6 { 7 fvScalarMatrix alpha1Eqn 8 ( 9 fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1) 10 + fv::gaussConvectionScheme<scalar> 11 ( 12 mesh, 13 phi, 14 upwind<scalar>(mesh, phi) 15 ).fvmDiv(phi, alpha1) //在这个位置添加显式的源项表达式 16 ); 17 18 solve(alpha1Eqn); 19 20 Info<< "Phase-1 volume fraction = " 21 << alpha1.weightedAverage(mesh.Vsc()).value() 22 << " Min(" << alpha1.name() << ") = " << min(alpha1).value() 23 << " Max(" << alpha1.name() << ") = " << max(alpha1).value() 24 << endl;之前在上述代码注释处添加了源项（那时候看不懂代码，就发现这里是相方程。。。），这个源项比较简单，就是一个时间的简单函数，且之前的结果看来稳定性还算可以，不过当时写程序的时候并没有在后面的MULES::correct中加入源项表达式，目前没有评估出这个差异会有多大。
此外，东岳老师，我还有个小疑问，这个bool变量MULESCorr决定的两种求解方式的区别大吗？您有做过类似的测试吗？谢谢！
