请教一个问题,关于Foam::fvMatrix<Type>::A() 函数的



  • template<class Type>
    Foam::tmp<Foam::volScalarField> Foam::fvMatrix<Type>::A() const
    {
        tmp<volScalarField> tAphi
        (
            new volScalarField
            (
                IOobject
                (
                    "A("+psi_.name()+')',
                    psi_.instance(),
                    psi_.mesh(),
                    IOobject::NO_READ,
                    IOobject::NO_WRITE
                ),
                psi_.mesh(),
                dimensions_/psi_.dimensions()/dimVol,
                zeroGradientFvPatchScalarField::typeName
            )
        );
    
        tAphi().internalField() = D()/psi_.mesh().V();
        tAphi().correctBoundaryConditions();
    
        return tAphi;
    }
    

    为什么要除以V()???感觉跟李东岳的icoFoam解析对不上了。
    tAphi().internalField() = D()/psi_.mesh().V();



  • 为了跟动网格的处理在代码形式上统一吗?还是为了与边界处理的一致?



  • 试着回答一下。本身A()这个成员函数的目的就是返回代数矩阵的格心系数。而实际当中fvmatrix在构造时都已经包含了网格体积,即完成了体积分或者面积分运算。因此在取回格心系数(central coeffcient)时需要格外剔除网格体积项。否则在内网格部分就和D()这个成员函数重复了。





  • @youmengtian 你还记得在哪些文件或者代码里描述了这种处理的吗?如果可以的话,指点一下。谢谢



  • simpleFoam



  • @yuan_neu

    在icoFoam中插入以下代码

    fvVectorMatrix UEqn
    (
        fvm::ddt(U)
        +fvm::div(phi,U)
        -fvm::laplacian(nu,U)
    )
    Info<<"fvc::laplacian(nu,U).dimensions="<<fvc::laplacian(nu,U)().dimensions()<<endl;
    Info<<"fvm::laplacian(nu,U).dimensions="<<fvm::laplacian(nu,U)().dimensions()<<endl;
    Info<<"phi.dimensions="<<phi.dimensions()<<endl;
    Info<<"U.dimensions="<<U.dimensions()<<endl;
    Info<<"fvc::div(phi,U).dimensions="<<fvc::div(phi,U)().dimensions()<<endl;
    Info<<"fvm::div(phi,U).dimensions="<<fvm::div(phi,U)().dimensions()<<endl;
    Info << "Dimension of UEqn is "<<UEqn.dimensions()<<endl;
    // ...
    

    输出大概是

    fvc::laplacian(nu,U):[0 1 -2 0 0 0 0]
    fvm::laplacian(nu,U):[0 4 -2 0 0 0 0]
    phi:[0 3 -1 0 0 0 0]
    U:[0 1 -1 0 0 0 0]
    fvc::div(phi,U):[0 1 -2 0 0 0 0]
    fvm::div(phi,U):[0 4 -2 0 0 0 0]
    UEqn:[0 4 -2 0 0 0 0]
    UEqn.A():[0 0 -1 0 0 0 0]
    UEqn.H():[0 1 -2 0 0 0 0]
    

    可以看出:

    • 实际的情况是fvm项是乘以了体积的积分方程离散,输出类型是fvMatrix,fvc项是没有乘以体积的微分方程离散,输出类型是vol<Type>Field

    • 从后面的代码UEqn == -fvc::grad(p)来看,二者是可以用operator==连接的。

    • 根据fvMatrix的operator==源代码可以发现,operator==将GeometricField加到fvMatrix中的时候乘以了体积。

      /*
      /src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
      */
      //line: 1458
      template<class Type>
      Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
      (
          const fvMatrix<Type>& A,
          const DimensionedField<Type, volMesh>& su
      )
      {
          checkMethod(A, su, "==");
          tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
          tC.ref().source() += su.mesh().V()*su.field(); //!!!
          return tC;
      }
      
    • fvMatrix的部分函数比较诡异,UEqn.A()的量纲和UEqn相差[0 4 -1 0 0 0 0],而UEqn.A()和H()相差待求变量速度的量纲[0 1 -1 0 0 0 0],合理的解释是:

      • UEqn是积分方程的离散,对于icoFoam,待求变量是速度U,量纲是[0 1 -1 0 0 0 0];
      • UEqn方程系数和phi的量纲一致,是[0 3 -1 0 0 0 0],而方程右端源项是[0 4 -2 0 0 0 0]
      • UEqn.A()是对角线系数除以体积,所以是[0 0 -1 0 0 0 0]
      • UEqn.H()是非对角线部分作用于待求变量并被方程右端源项减去,再除以体积得到的,所以是[0 1 -2 0 0 0 0]
    • 如果你详细看源代码,你会发现UEqn.D()是无量纲的,而UEqn.A()是强行把量纲附上去的,我觉得OF这些地方是有偷懒或者不一致的地方的。



  • @程迪 :kiss_big: :kiss_big: :kiss_big:



  • @程迪 如果你详细看源代码,你会发现UEqn.D()是无量纲的,而UEqn.A()是强行把量纲附上去的,我觉得OF这些地方是有偷懒或者不一致的地方的。

    这句话我也有所感觉。


登录后回复
 

与 CFD中文网 的连接断开,我们正在尝试重连,请耐心等待