计算作用于cellZone上的压力和剪切力

[liujm]

此处附上改写的代码，烦请各位老师指出错误所在，不胜感激：

    if (porosity_)
    {   
        const volVectorField& U  = lookupObject<volVectorField>(UName_);
        const volScalarField& p  = lookupObject<volScalarField>(pName_);   // Add code.
        const volScalarField rho(this->rho());
        const volScalarField mu(this->mu());

        const HashTable<const porosityModel*> models = obr_.lookupClass<porosityModel>();

        if (models.empty())
        {
            WarningInFunction
                << "Porosity effects requested, but no porosity models found "
                << "in the database"
                << endl;
        }

        forAllConstIters(models, iter)
        {
              // Non-const access required if mesh is changing
            porosityModel& pm = const_cast<porosityModel&>(*iter());

            vectorField fPTot(pm.porousForce(U, rho, mu));  // Modified code.
            
            const labelList& cellZoneIDs = pm.cellZoneIDs();
            
            for (const label zonei : cellZoneIDs)
            {
                const cellZone& cZone = mesh_.cellZones()[zonei];  // Get cells in zonei.
                const vectorField d(mesh_.C(), cZone);
                const vectorField fP(fPTot, cZone);
                
                const vectorField Md(d - coordSys_.origin());
                const labelList& cellsID = cZone;  // Add code: Get cellsID(index) list in zonei.
                vectorField fN(cZone.size(), Zero);
                vectorField fT(cZone.size(), Zero);
                scalar                  pRef                            = pRef_/rho(p);   // Add code: Scale pRef by density for incompressible simulations
                tmp<volSymmTensorField> tdevRhoReff                     = devRhoReff();
                const                   volSymmTensorField& devRhoReff1 = tdevRhoReff();
                
                Info << "cellZoneIDs zonei = " << zonei << "\n"<< endl;
                Info << "cells in zonei: cZone.size() = " << cZone.size() << "\n"<< endl;
                Info << "fN.size() = " << fN.size() << "\n"<< endl;
                Info << "fT.size() = " << fT.size() << "\n"<< endl;
                Info << "fP.size() = " << fP.size() << "\n"<< endl;

                forAll (cellsID, celli)
                {
                    label cellID               = cellsID[celli];
                    const labelList& cellFaces = mesh_.cells()[cellID];
                    forAll (cellFaces, facei)
                    {
                        label faceID = cellFaces[facei];
                          // fN.
                        vector cellFacesFN  = rho(p) * mesh_.Sf()[faceID] * (p[cellID] - pRef);
                        fN    [celli]      += cellFacesFN;
                          // fT.
                        vector cellsFacesfT  = mesh_.Sf()[faceID] & devRhoReff1[faceID];
                        fT    [celli]       += cellFacesFT;
                        
                        Info << "No." << celli << "cell's ID(index) is " << cellID << "\n" << endl;
                        Info << "fN = " << fN[celli] << "\n" << endl;
                        Info << "fT = " << fT[celli] << "\n" << endl;
                        Info << "fP = " << fP[celli] << "\n" << endl;
                    }
                }                     
                addToFields(cZone, Md, fN, fT, fP);

                applyBins(Md, fN, fT, fP, d);
            }
        }
    }

[liujm]

总体而言，就是要计算cellZone所有网格上的法向力和切向力，而不是仅仅计算作用在固壁边界上的力。
这里通过图示表达目的：

[李东岳]

你的力定义在面上，如何对你的cellZone里面的所有cells取力？

[liujm]

@李东岳 老师，可以直接对区域内的所有cell压力求积分得出吗？其实这篇帖子是这个帖子的延申 。
一方面，要保证动网格能够识别到patch名，另一方面也要计算cell上的力，请问李老师有没有更好的实现方法？

[李东岳]

一般力都定义在面上，压力、表面张力
当然存在一些体积力，比如重力，但不会产生剪切
但是你好像需要一个产生剪切的体积力，这是什么情况？

[liujm]

@李东岳 老师，我想要模拟浮式多孔介质模型在流体中的运动，因此需要获得物体的力（矩），这里的剪切力其实对应一楼【代码片段1】的54行和56行。

[李东岳]

你这个图里面，第一个图是可以理解的，因为force定义在面上，因此可以做面积分。

第二个myForce，你要对cell做操作，如何做面积分？

[liujm]

@李东岳 李老师，这也正是我遇到的问题 。

退一步我做一个简化，在图示的myForces.H的结构最外层cell上套一层patch边界，可不可以仅用patch计算cellZone最外层网格的力和力矩[1]，而不阻碍水体穿越被patch包裹的cellZone区域？

这个想法可行吗，如果可行的话老师方便指点一下，如何实现水体不穿越patch的功能？

[1] 只计算外层网格表面的受力是因为网衣通常不会太厚，有文献采用2m厚的多孔介质模型，网格尺寸0.5m来模拟。

在OF2206中，
对于patch，有addToPatchFields函数：

void Foam::functionObjects::forces::addToPatchFields
(
    const label patchi,
    const vectorField& Md,
    const vectorField& fP,
    const vectorField& fV
)
{
    sumPatchForcesP_ += sum(fP);
    sumPatchForcesV_ += sum(fV);
    force().boundaryFieldRef()[patchi] += fP + fV;

    const vectorField mP(Md^fP);
    const vectorField mV(Md^fV);

    sumPatchMomentsP_ += sum(mP);
    sumPatchMomentsV_ += sum(mV);
    moment().boundaryFieldRef()[patchi] += mP + mV;
}

对于多孔介质的cellZone，有addToInternalField函数：

void Foam::functionObjects::forces::addToInternalField
(
    const labelList& cellIDs,
    const vectorField& Md,
    const vectorField& f
)
{
    auto& force = this->force();
    auto& moment = this->moment();

    forAll(cellIDs, i)
    {
        const label celli = cellIDs[i];

        sumInternalForces_ += f[i];
        force[celli] += f[i];

        const vector m(Md[i]^f[i]);
        sumInternalMoments_ += m;
        moment[celli] = m;
    }
}

[liujm]

已经想到合适的解决方案啦！ 封楼

[Linzd]

@liujm 您好，您最终采用了什么方案呢

[liujm]

@Linzd 多孔介质区域在toposet的时候本身只是选中了cell集合，只有在壁面上才能通过原始代码计算法向力和切向力。