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

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类型源项，源项量纲和程序量纲不匹配！