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

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