CFD中文网

@李东岳 静态网格也是有这个问题的。按理说不应该啊...我再跑跑其他类似的原生算例

p：

#include "../varsSettings" //Include global variables

dimensions	[0 2 -2 0 0 0 0];

internalField	uniform	0;

boundaryField
{
    inlet
    {
        type            zeroGradient;
    }

    outlet
    {
        type            fixedValue;
        value           $internalField;
    }

    frontAndBack
    {
        type            empty;//zeroGradient;
    }

    topAndBottom
    {
        type            zeroGradient;
    }
    
    frame
    {
        type            zeroGradient;//fixedFluxPressure;
    }
}

U：

#include "../varsSettings" //Include global variables

dimensions	[0 1 -1 0 0 0 0];

internalField	uniform	$flowVelocity;

boundaryField
{
    inlet
    {
        type		fixedValue;
        value		$internalField;
    }

    outlet
    {
        type		inletOutlet;
        inletValue	uniform (0 0 0);
        value		$internalField;
    }

    frontAndBack
    {
        type		empty;//slip;
        //value		uniform (0 0 0);	
    }

    topAndBottom
    {
        type		slip;
        value		uniform (0 0 0);
    }
    
    frame
    {
        type            noSlip;//movingWallVelocity;
        //value           uniform (0 0 0);
    }
}

算例文件 OpenfOAM2012 网格数4w
functionalTest_NoPZRefine2_2D_staticMesh.zip

@李东岳 李老师好，这里我请教一个问题->链接

为了简化算例，我把网格三维转化为二维，简化了结构形式，以便排错。

【算例云图】

【简化版算例下载】
functionalTest_NoPZRefine2_2D.zip
生成网格createMesh，
跑算例runCaseOnly，
全流程runCase

感谢~

@李东岳 是的老师，这里我为了排错，先把动网格用constraints六自由度都固定住。
以下附图：

各位老师好，

【问题描述】
本人在使用pimpleFoam模拟固定式网框在自由来流（u=0.5m/s）中的响应。整个流场的计算结果符合预期，可是网框表面处的流速竟然高达1e2量级，明显与现实不符，结果云图如下：（侧视图）

【分析】
网框截面原本比本帖中展示的要更小，起初正方形边长仅分布四个网格，为了排除网格数量对结果的影响，特意扩大了截面，增加边长网格数量。但是现在依然存在这个问题。

【初始条件设置】
以下摘抄0文件夹设置：
p：

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
	version     2.0;
	format      ascii;
	class       volScalarField;
	location    "0";
	object      p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#include "../varsSettings" //Include global variables

dimensions	[0 2 -2 0 0 0 0];

internalField	uniform	0;

boundaryField
{
    inlet
    {
        type            zeroGradient;
    }

    outlet
    {
        type            fixedValue;
        value           $internalField;
    }

    frontAndBack
    {
        type            zeroGradient;
    }

    topAndBottom
    {
        type            zeroGradient;
    }
    
    frame
    {
        type            fixedFluxPressure;
    }
}


// ************************************************************************* //

U:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volVectorField;
    location    "0";
    object      U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#include "../varsSettings" //Include global variables

dimensions	[0 1 -1 0 0 0 0];

internalField	uniform	$flowVelocity;

boundaryField
{
    inlet
    {
        type		fixedValue;
        value		$internalField;
    }

    outlet
    {
        type		inletOutlet;
        inletValue	uniform (0 0 0);
        value		$internalField;
    }

    frontAndBack
    {
        type		slip;
        value		uniform (0 0 0);	
    }

    topAndBottom
    {
        type		slip;
        value		uniform (0 0 0);
    }
    
    frame
    {
/*        type            noSlip; */
        type            movingWallVelocity;
        value           uniform (0 0 0);
    }
}


// ************************************************************************* //

pointDisplacement:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  5                                     |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       pointVectorField;
    object      pointDisplacement;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 1 0 0 0 0 0];

internalField   uniform (0 0 0);

boundaryField
{
	inlet
	{
		type            fixedValue;
		value           uniform (0 0 0);
	}

	outlet
	{
		type            fixedValue;
		value           uniform (0 0 0);
	}
	
	frontAndBack
	{
		type            fixedValue;
		value           uniform (0 0 0);
	}
	
	topAndBottom
	{
		type            fixedValue;
		value           uniform (0 0 0);
	}

	frame
	{
		type            calculated;
		value           uniform (0 0 0);
	}
}

// ************************************************************************* //

【算例设置】
functionalTest_NoPZRefine1.zip
生成网格createMesh，
跑算例runCaseOnly，
全流程runCase

希望各位老师能给出建议~@李东岳

@liujm
问题已解决，添加libs("libfvMotionSolvers.so")，或者参考这里

@Joker
您好，请问您使用的OF是什么版本呢？我在用displacementLaplacian后，使用moveDynamicMesh预览网格运动的时候，系统报错，似乎没有displacementLaplacian这个求解器。请问是我的字典设置有误吗？

Create time

Create mesh for time = 0

Selecting dynamicFvMesh dynamicMotionSolverFvMesh
Selecting motion solver: displacementLaplacian


--> FOAM FATAL IO ERROR: (openfoam-2012)
Unknown solver type displacementLaplacian

Valid solver types :
7(coded displacementInterpolation displacementLayeredMotion multiSolidBodyMotionSolver sixDoFRigidBodyMotion solidBody velocityDisplacement)


file: /OpenFOAM-v2012/myProject/tutorials/floatingPorousPimpleFoam/Condition1_Uin=1.2/constant/dynamicMeshDict at line 18 to 26.

    From static Foam::autoPtr<Foam::motionSolver> Foam::motionSolver::New(const Foam::polyMesh&, const Foam::IOdictionary&)
    in file motionSolvers/motionSolver/motionSolver.C at line 137.

FOAM exiting

@coolhhh 厉害，测试代码copy了。不过，这里直接操作就又可以了

@coolhhh
感谢解答，我之所以直接操作mesh.C()，是因为我在源码中搜到了这种用法，OF2012

@卡洛 卡老师，能分享一下边界条件的修改思路吗

mesh_.C()
mesh_.C().primitiveField()
mesh_.C().internalField()
mesh_.C().boundaryField()
这四个返回值之间的区别是什么呢？

各位前辈大家好，

我在使用const vector cellCentre = mesh_.C()[celli]收集网格中心坐标时，编译器出现报错：

In file included from cfdTools/general/porosityModel/floatingDarcyForchheimer/floatingDarcyForchheimer.H:216,
                 from cfdTools/general/porosityModel/floatingDarcyForchheimer/floatingDarcyForchheimer.C:30:
cfdTools/general/porosityModel/floatingDarcyForchheimer/floatingDarcyForchheimerTemplates.C: In member function ‘void Foam::porosityModels::floatingDarcyForchheimer::applyInCalcForce(Foam::scalarField&, Foam::vectorField&, const scalarField&, const RhoFieldType&, const scalarField&, const vectorField&, Foam::vectorField&) const’:
cfdTools/general/porosityModel/floatingDarcyForchheimer/floatingDarcyForchheimerTemplates.C:119:48: error: no match for ‘operator[]’ (operand types are ‘const volVectorField’ {aka ‘const Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>’} and ‘const label’ {aka ‘const int’})
  119 |             const vector cellCentre = mesh_.C()[celli];
      |  

貌似mesh_.C()返回的数据类型不支持使用[]符号，可是为什么在其他代码中就可以这样使用 mesh_.C()[celli]呢？

后期，我通过使用const vector cellCentre = mesh_.C().primitiveField()[celli]通过编译。不过尚不清楚其与mesh_.C()返回的结果是否存在差异。

烦请各位老师解答~

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

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

退一步我做一个简化，在图示的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;
    }
}

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

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

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

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

    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);
            }
        }
    }

各位老师和前辈们好，

我想要通过划定cellZone区域来计算作用于该区域网格上的压力和剪切力，预期功能类似于forces.H中的void Foam::functionObjects::forces::calcForcesMoments()函数，然而，该函数是基于边界patch计算的【代码片段1】，我想要基于cellZone实现。而计算多孔介质部分的力【代码片段2】正是基于cellZone实现的，因此，

• i) 我想将基于patch的计算边界力的方法移植于多孔介质模块中，这样可行吗？

• ii) 评论区我放出改写的代码，烦请各位老师指出需要改进的地方

（如果有老师能够直接给出示例就更好了 ）

以下摘取了forces.H中void Foam::functionObjects::forces::calcForcesMoments()部分源码：
【代码片段1】

    if (directForceDensity_)
    {
        const auto& fD = lookupObject<volVectorField>(fDName_);

        const auto& Sfb = mesh_.Sf().boundaryField();

        for (const label patchi : patchSet_)
        {
            const vectorField& d = mesh_.C().boundaryField()[patchi];

            const vectorField Md(d - origin);

            const scalarField sA(mag(Sfb[patchi]));

            // Pressure force = surfaceUnitNormal*(surfaceNormal & forceDensity)
            const vectorField fP
            (
                Sfb[patchi]/sA
               *(
                    Sfb[patchi] & fD.boundaryField()[patchi]
                )
            );

            // Viscous force (total force minus pressure fP)
            const vectorField fV(sA*fD.boundaryField()[patchi] - fP);

            addToPatchFields(patchi, Md, fP, fV);
        }
    }
    else
    {
        const auto& p = lookupObject<volScalarField>(pName_);

        const auto& Sfb = mesh_.Sf().boundaryField();

        tmp<volSymmTensorField> tdevRhoReff = devRhoReff();
        const auto& devRhoReffb = tdevRhoReff().boundaryField();

        // Scale pRef by density for incompressible simulations
        const scalar rhoRef = rho(p);
        const scalar pRef = pRef_/rhoRef;

        for (const label patchi : patchSet_)
        {
            const vectorField& d = mesh_.C().boundaryField()[patchi];

            const vectorField Md(d - origin);

            const vectorField fP
            (
                rhoRef*Sfb[patchi]*(p.boundaryField()[patchi] - pRef)
            );

            const vectorField fV(Sfb[patchi] & devRhoReffb[patchi]);

            addToPatchFields(patchi, Md, fP, fV);
        }
    }

【代码片段2】

if (porosity_)
    {
        const auto& U = lookupObject<volVectorField>(UName_);
        const volScalarField rho(this->rho());
        const volScalarField mu(this->mu());

        const auto 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
            auto& pm = const_cast<porosityModel&>(*iter());

            const vectorField fPTot(pm.force(U, rho, mu));

            const labelList& cellZoneIDs = pm.cellZoneIDs();

            for (const label zonei : cellZoneIDs)
            {
                const cellZone& cZone = mesh_.cellZones()[zonei];

                const vectorField d(mesh_.C(), cZone);
                const vectorField fP(fPTot, cZone);
                const vectorField Md(d - origin);

                addToInternalField(cZone, Md, fP);
            }
        }
    }

祝好~

RT，请教大佬们，网上基本都是获取与当前网格相邻信息的代码，请问有没有类似mesh.cells()的方法获取网格的编号（索引）？

@HITSC30 感谢回答，两篇文章写得很详细！