计算作用于cellZone上的压力和剪切力
-
各位老师和前辈们好，
我想要通过划定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); } } }祝好~
-
-
此处附上改写的代码，烦请各位老师指出错误所在，不胜感激：
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); } } } -
-
@李东岳 李老师，这也正是我遇到的问题
。退一步我做一个简化，在图示的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; } }
