设置vectorCodedSource类型源项，源项量纲和程序量纲不匹配！

lizhisongsjtu

用interFoam计算两相流，在计算域内选择一块区域作为源区域，施加动量源项。
在 system/fvOptions 中添加了 vectorCodedSource类型的源项，并参考了网站链接的方法设置了源项。

根据报错信息，初步确定是添加源项的量纲出了问题。请高手指点，如何解决这个问题

最核心的一句报错信息是

--> FOAM FATAL ERROR: (openfoam-2206)
Not implemented

    From virtual void Foam::fv::codedMomentumSourceFvOptionvectorSource::addSup(const volScalarField&, Foam::fvMatrix<Foam::Vector<double> >&, Foam::label)

附fvOptions代码

codedSource
{
    type            vectorCodedSource;
    active          true;
    name            codedMomentumSource;

    // selectionMode   all;
    selectionMode   cellSet;
    cellSet         myUSourceCellsSet;  // 可以用 topoSet 定义的 cellSet
    fields          
    (
        U
    );

    codeAddSup
    #{  
        // https://www.cfd-china.com/topic/3255/%E5%88%86%E5%8C%BA%E7%BD%91%E6%A0%BC%E7%9A%84%E4%B8%80%E5%B0%8F%E6%AE%B5%E4%BB%A3%E7%A0%81?_=1731583465684
        Info << "++++From codedSource in fvOptions" << nl;

        const label cellZoneID = mesh().cellZones().findZoneID("myUSourceCellsSet"); // find ID for the cellZone "BANANA"
        
		const cellZone& zone = mesh().cellZones()[cellZoneID];
		const cellZoneMesh& zoneMesh = zone.zoneMesh();
		const labelList& cellsZone = zoneMesh[cellZoneID]; //list of all cellZone cell ID's  
        
        const Time& time = mesh().time();
        const scalarField& V = mesh_.V(); // Cell volume
            
        Foam::Field<Foam::Vector<double> >& USource = eqn.source(); // get source from fvMatrix
		//scalarField & Udiag = eqn.diag(); // get diagonal of fvMatrix
                
        const scalarField& rho = mesh().lookupObject<scalarField>("rho"); // get density field
        
        const vector g (0, 0, -9.81);      // gravitational vector


        forAll(cellsZone,i)
        {
            const label celli = cellsZone[i];
        
            // USource[celli] -= rho[celli]*g/V[celli] * sin(time.value());
            USource[celli] = sin(time.value());
        }
    #};

    //  chatGPT 的解决方案 + CFD-online解决方案
    // codeAddSup
    // #{
    //     Info << "++++From codedSource in fvOptions" << nl;

    //     const scalarField& V = mesh_.V();
    //     vectorField& USource = eqn.source();
    //     const scalarField& cellx = mesh_.C().component(0) ; 
    //     const scalarField& celly = mesh_.C().component(1) ;
        
    //     forAll(USource, cellI)
    //     {
    //         Info << "cellI = " << cellI << nl;
    //         // scalar v1 = 512*( 2*pow(cellx[i],6) - 6*pow(cellx[i],5) + 7*pow(cellx[i],4)
    //         //                  - 4*pow(cellx[i],3) + pow(cellx[i],2))*pow(celly[i],7) ;
    //         // // .... and the other 9 scalars
            
    //         // // Please ADD to the source, don't overwrite it.
    //         // // Overwriting it means you don't care about source terms in 
    //         // // DESCRETISATION OPERATORS (laplacian, div ...)
    //         USource[cellI][0] += 1.0 * USource[cellI][1];
    //     }
    // #};

    codeInclude
    #{
        #include "fvCFD.H"
    #};

    codeCorrect
    #{
        // Use Correct to make adjustments to the U values after they get calculated
        Info << "++++codeCorrect" << nl;
    #};
    codeConstrain
    #{
        // Use Constrain to modify the U equation before it gets solved
        Info << "++++codeConstrain" << nl;
    #};
}

lizhisongsjtu

addSup的构造函数

//- Explicit/implicit matrix contributions
        virtual void addSup
        (
            fvMatrix<Type>& eqn,
            const label fieldi
        );

        //- Explicit/implicit matrix contributions to compressible equation
        virtual void addSup
        (
            const volScalarField& rho,
            fvMatrix<Type>& eqn,
            const label fieldi
        );

学流体的小明

加源项的时候，这个源项的表达式好像有点问题，比如下面添加vectorSemiImplicitSource的时候源项是( (299.4600 0 0) 0)而不是(299.4600 0 0)。原因是加入的源项表达式为：
$$
S=S_C+S_Pϕ_P
$$
标量源项的代码也是这么写的：

        injectionRateSuSp
        {
            variable_name     (Sc Sp);
        }

momentumSource
{
   type vectorSemiImplicitSource;
   active on;
   

   vectorSemiImplicitSourceCoeffs
   {
      selectionMode all;
      //volumeMode      absolute; // specific
      volumeMode        specific;
      injectionRateSuSp
      {
         U           ( (299.4600 0 0) 0); //partial p / partial x
      }
   }
}

两个参考：
https://caefn.com/openfoam/fvoptions-acousticdampingsource
https://xiaopingqiu.github.io/2016/03/20/fvOptions2/

刚又找了一下之前的笔记，不是上面这个原因，因为你用的是vectorCodedSource。我之前也刚好写过scalarCodedSource，也发现量纲不对导致的问题。之前记录的笔记是：
有一个问题是感觉并不需要乘以cell的体积。官方给的文档里面是有乘以体积的，但这个网页就很迷惑，https://www.openfoam.com/documentation/guides/v2012/doc/guide-fvoptions-sources-coded.html ，上面的表达式是除以体积，下面的代码又是乘以体积了。找到一个人的注释 https://xiaopingqiu.github.io/2016/03/20/fvOptions2/

// fvMatrix<Type> 类中对“+=”操作符进行了重载，所以，eqn与Su的相加，相当于eqn+Su*mesh.V()，要不然eqn与Su的量纲不一致。
eqn += Su + fvm::SuSp(Sp, psi);

所以还是要乘以体积。

lizhisongsjtu

上传一下case文件

interFoam.sourceU.zip

李东岳

用下面这个验证过的代码模板，用来模拟树木冠层的。这个跟实验能对上，是要乘以体积的。

//fvModels如下指定
USource
{
    type            coded;
    selectionMode   all;
    field           U;

    codeAddSup
    #{
        const volVectorField& U = eqn.psi();
        vectorField& USource = eqn.source();

        const fvMesh& mesh = U.mesh();
        const DimensionedField<scalar, volMesh>& V = mesh.V();
        const volVectorField& C = mesh.C();
        scalar A = 0.919;
        scalar Cd = 0.15;

        forAll(U, i)
        {
            vector position = C[i];
            scalar pz = position.z();
            if (pz < 10.0)
            {
                USource[i] -= -Cd*A*mag(U[i])*U[i]*V[i];
            }
        }
    #};
}

http://dyfluid.com/code.html 最下面

我2021年那个代码主要是测试多区域网格分区赋值的。比如mesh分为2个zone，然后对某个zone处理，不过这部分也更新到下面了

const scalarField& V = mesh.V();                            // 网格体积

const surfaceVectorField& Sf = mesh.Sf();                   // 网格面矢量

const surfaceVectorField& Cf = mesh.Cf();                   // 网格面心位置

const surfaceScalarField& magSf = mesh.magSf();             // 网格面积的模

const surfaceVectorField& N = Sf/magSf;                     // 网格面法向量

const label& nCells = mesh.nCells();                        // 网格单元数量

const label& nInternalFaces = mesh.nInternalFaces();        // 网格内部面数量

const meshCellZones& cellZones = mesh.cellZones();          // 网格cellZone编号

wordList zoneNames = mesh.cellZones().names();              // 网格cellZone的名字

label zoneI =  mesh.cellZones().whichZone(celli);           // 判断celli在哪个cellZone里

label i = cellZones[0][i];                                  // 网格cellZone[0]的真实网格编号

const fvPatchList& patches = mesh.boundary();               // 边界网格，是一系列patch的集合

const polyBoundaryMesh& boundaryMesh = mesh.boundaryMesh(); // 边界网格，是一系列patch的集合

const label patchID = boundaryMesh.findPatchID("inlet");    // inlet边界的ID

volScalarField S(magSqr(symm(tgradU())));                   // 形变率的双点积
tmp<volTensorField> tgradU = fvc::grad(U);

李东岳

还有一个就是注意USource[i] -= -Cd*A*mag(U[i])*U[i]*V[i];这个是减号，不是加号。

lizhisongsjtu

感谢回复解惑。

1楼、2楼和4楼，即代码和算例，是在of2206上测试的，失败。

下面的一段代码，在of2206上仍然失败。在of8上暂时可以运行，源项效果有待检验。此上！

codedSource
{
    type            vectorCodedSource;
    name            sourceTime;

    selectionMode   cellSet;
    cellSet         myUSourceCellsZoneSet;
    
    fields          
    (
        U
    );
    
    codeAddSup
    #{
        const vectorField& C = mesh_.C();
        const scalarField& V = mesh_.V();
        vectorField& Usource = eqn.source();
        forAll(C, i)
        {
            // Info << "cellI = " << i << nl;
            const scalar y = C[i][0];
            const scalar vol = V[i];
            const scalar fx = vol * 1.0;

            Usource[i] -= vector(1.0, 0, 0);
        }
        Info << "***codeAddSup***" << nl;
    #};

    codeConstrain
    #{
    #};

    codeCorrect
    #{
    #};
}

李东岳

下面的一段代码，在of2206上仍然失败。在of8上暂时可以运行，源项效果有待检验。此上

我在2206上测试了一下，可以运行。

lizhisongsjtu

@李东岳我想确认一下，7楼的代码真的在2206上运行成功了吗？谢谢

附：of2206、of8和of12的日志。of2206失败。of8直接可以用7楼的代码，of12需放入fvModels中并略作修改。

log.interFoam.2206

log.interFoam.8

log.foamRun.12

李东岳

我想确认一下，7楼的代码真的在2206上运行成功了吗？谢谢

可以，不过我没用的interfoam，我用的什么求解器我忘了。可能是simplefoam、icoFoam之类。你试试这俩能不能运行？

lizhisongsjtu

解决方案：

在of8及同系列版本中，单相流和多相流均使用codeAddSup定义源项。
在of2206及同系列版本中，单相流使用codeAddSup定义源项，多相流使用codeAddSupRho定义源项。

李东岳

大佬这个流弊！屌爆了

学流体的小明

学到了学到了

CFD中文网

设置vectorCodedSource类型源项，源项量纲和程序量纲不匹配！