interfoam中植入Fluent中溶质在单相扩散的功能

[麦迪文]

想在气液模型interfoam中实现Fluent里的液相溶质浓度扩散方程如下

暂不考虑源项
实现如下：

    fvScalarMatrix CEqn
    (
       fvm::ddt(alpha1, rho1,C)
      + fvm::div(rhophi1,C)
      - fvm::laplacian(alpha1*Ck, C)
    );

Ck是扩散系数，编译成功，运行出错，大家能不能给些建议

[麦迪文]

补充内容：
经测试把上述代码中时间项的相分数置为1，如下代码所示，程序能稳定运行

   fvScalarMatrix CEqn
    (
       fvm::ddt((alpha1+alpha2), rho1,C)
      + fvm::div(rhophi1,C)
      - fvm::laplacian(alpha1*Ck, C)
    );

感觉像是算法哪里出了问题？？？？？？

[麦迪文]

再次补充：
使用QT调试，查出问题

 void Foam::DILUPreconditioner::calcReciprocalD
 (
     scalarField& rD,
     const lduMatrix& matrix
 )
 {
     scalar* __restrict__ rDPtr = rD.begin();
 
     const label* const __restrict__ uPtr = matrix.lduAddr().upperAddr().begin();
     const label* const __restrict__ lPtr = matrix.lduAddr().lowerAddr().begin();
 
     const scalar* const __restrict__ upperPtr = matrix.upper().begin();
     const scalar* const __restrict__ lowerPtr = matrix.lower().begin();
 
     label nFaces = matrix.upper().size();

for (label face=0; face<nFaces; face++)
     {
         rDPtr[uPtr[face]] -= upperPtr[face]*lowerPtr[face]/rDPtr[lPtr[face]];//在这里报错。。。。。！！！！！！！！
     }

 
 

     label nCells = rD.size();
 
     for (label cell=0; cell<nCells; cell++)
     {
         rDPtr[cell] = 1.0/rDPtr[cell];
     }
 }

[麦迪文]

最后总结，给时间项中的相分数加上个极小值（1e-24），就能正常计算，是说明时间项中不能存在0系数？那两相欧拉模型中的时间项中的相分数是怎么处理的呢？

[李东岳]

你可以试试不要preconditioner，多相流里面也会出这个问题的，在你的fvSolution里面把preconditioner设置为none

另外好像你的对流项没包含相分数?

fvScalarMatrix CEqn
    (
       fvm::ddt(alpha1, rho1,C)
      + fvm::div(rhophi1,C)// + fvm::div(rhoAlphaPhi1,C)
      - fvm::laplacian(alpha1*Ck, C)
    );

[麦迪文]

对流项如下

rhoPhi1=fvc::interpolate(rho*alpha1)*phi;

在两相欧拉模型中的能量方程的矩阵结构在初始时刻和本文的差不多

   fvm::ddt(alpha1, rho1, he1) + fvm::div(alphaRhoPhi1, he1)
      - fvm::Sp(contErr1, he1)

不知其最后这项fvm::Sp(contErr1, he1)（隐式修正？），还有更后面有两项 heatTransferCoeff*he1/Cpv1
- fvm::Sp(heatTransferCoeff/Cpv1, he1是为了避免这类问题？

软件版本是5.x

[麦迪文]

谢谢老师
preconditioner设置为none后成功运行，并且暂时看上去和Fluent勉强趋势能勉强相近

[李东岳]

哦那你植入的可以

Sp是用来处理有界性的，如果越界可以用一下试试，没越界可以暂时不用管，