L
@李东岳 静态网格也是有这个问题的。按理说不应该啊...我再跑跑其他类似的原生算例
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
感谢~
各位老师好，
【问题描述】
本人在使用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
希望各位老师能给出建议~@李东岳
@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
各位前辈大家好，
我在使用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;
}
}
此处附上改写的代码，烦请各位老师指出错误所在，不胜感激：
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);
}
}
}
祝好~