OpenFOAM计算圆球绕流过程中,如何输出切向粘性阻力系数和法向粘性阻力系数?或是如何从输出的结果中计算得到?
-
各位老师好!
最近计算圆球绕流,通过在controlDict中添加functions,能在每步输出周围流体对界面施加的总阻力系数。
由于总阻力系数是“压差阻力系数”和“粘性阻力系数”两个分量的总和,进一步开启log选项后,能分别显示出这两项,如下:
也就是说:
总阻力系数$C_D=C_{D,p}+C_{D,\tau}$其中压差阻力系数$C_{D,p}=-4\int_0^\pi p\cos\theta\sin\theta d\theta $
粘性阻力系数$\begin{aligned}
C_{D,\tau}& =C_{D,\tau1}+C_{D_{\tau2}} \
&=4\int_0^\pi\tau_{r\theta}\sin\theta\sin\theta d\theta+4\int_0^\pi\tau_{rr}\cos\theta\sin\theta d\theta
\end{aligned}$我想运行的时候输出切向粘性阻力系数$C_{D,\tau1}$ 和法向粘性阻力系数$C_{D,\tau2}$
或是从运行完之后的流场结果中计算出来。forces.C和forces.H,forceCoeffs.C以及forceCoeffs.H文件还不太看得明白,希望能有老师指点,谢谢!
-
@李东岳 好的,谢谢东岳老师。
我也自己看了下[forces.C、forces.H,forceCoeffs.C、forceCoeffs.H]四个文件,正如您得教材里所说的,在forceCoeffs.C中,总阻力的计算是:
$Cd=\frac{F_{drag}}{(\sum|\mathbf{S}_f|)\frac{1}{2}\rho|\mathbf{U}|^2}$
顺着这个我在forces.C中找了下各个力的计算公式:
其关系到的几个变量fN、fT、fP以及Md,在forces.H最先出现,
但是整个文档都看了下,好像都没有给出这几个量的定义,类似的还有pf、pm等:
东岳老师,各位大佬,能否指点一下这几个变量通常要去如何理解它的物理意义?
谢谢! -
@Prometheus10 在 OpenFOAM计算圆球绕流过程中,如何输出切向粘性阻力系数和法向粘性阻力系数?或是如何从输出的结果中计算得到? 中说:
但是整个文档都看了下,好像都没有给出这几个量的定义,类似的还有pf、pm等:
forces.C中都有定义吧,在calcForcesMoments函数里(下面这个是新版的定义,我记得以前fP是表示porosity的)
void Foam::functionObjects::forces::calcForcesMoments() { initialise(); reset(); const point& origin = coordSysPtr_->origin(); 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); } } 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); } } } reduce(sumPatchForcesP_, sumOp<vector>()); reduce(sumPatchForcesV_, sumOp<vector>()); reduce(sumPatchMomentsP_, sumOp<vector>()); reduce(sumPatchMomentsV_, sumOp<vector>()); reduce(sumInternalForces_, sumOp<vector>()); reduce(sumInternalMoments_, sumOp<vector>()); }
@Prometheus10 在 OpenFOAM计算圆球绕流过程中,如何输出切向粘性阻力系数和法向粘性阻力系数?或是如何从输出的结果中计算得到? 中说:
我想运行的时候输出切向粘性阻力系数和法向粘性阻力系数
这个用不着改forces functionObject ,可以:
- 用wallshearstress后处理算壁面粘性剪切力,然后在paraview中算在壁面各点法向和切向的投影;
- 基于wallshearstress functionObject进行修改,输出投影后的结果
-
@lwjetmann 感谢老师,组里同学大家都不太会OpenFOAM,我马上也看一下。再次感谢。
-
@lwjetmann 老师您好,再请教一下您。
1. 用wallshearstress后处理算壁面粘性剪切力,然后在paraview中计算出在壁面各点法向和切向的投影; 2. 基于wallshearstress functionObject进行修改,输出投影后的结果
按照您提到的上述方案1,对于圆球绕流(球体界面是“noSlip”边界条件),我顺利的输出了粘性剪切力,结果也验证的上。
想请您再帮忙看看另外一种情况:绕流气泡,即将边界条件改为Symmetry:
因为wallshearstress只对wall类型的边界计算粘性剪应力,因此不能再通过【方案1】输出界面上的粘性剪应力结果。
关于设置Symmetry边界条件的原因是源于这里:Stokes流中的圆球绕流的阻力系数,怎么设置才能算准。对此,针对绕流气泡的结果(边界条件为Symmetry):
我尝试手动在结果中将Symmetry边界条件改为Wall,然后执行一次wallshearstress; 通过对比文献中气泡的压差阻力系数、法向粘性应力和切向粘性应力结果,输出的结果并不对
除此以外,是否只有像您提到的【方案二】基于wallshearstress functionObject进行修改,才能输出“绕流气泡”的粘性剪切力结果。
如果您方便的话,能否请您指点一些修改思路,谢谢老师。 -
@Prometheus10 在 OpenFOAM计算圆球绕流过程中,如何输出切向粘性阻力系数和法向粘性阻力系数?或是如何从输出的结果中计算得到? 中说:
我尝试手动在结果中将Symmetry边界条件改为Wall,然后执行一次wallshearstress;
通过对比文献中气泡的压差阻力系数、法向粘性应力和切向粘性应力结果,输出的结果并不对是的,因为wallshearstress代码会自动过滤非wall边界,可以手动改为wall。既然计算剪切力的公式都是一样的,这么改一下不会影响计算正确性的。和文献对不上是不是symmetry边界条件的问题呢?(你们这个帖子挺长,我看大意是不是:楔形网格算球noslip正确,但算slip和symmetry和文献对不上,且symmetry效果更好些?)
用整球网格算会如何呢?做一些其它情况的验证呢,比如壁面滑移的方腔流动?我之前做过低雷诺数下的一阶maxwell速度滑移的圆柱绕流,算出来的受力以及粘性正应力占比和文献比得上。不过你的自由滑移情况我没考虑过,不知道为啥出问题。
@Prometheus10 在 OpenFOAM计算圆球绕流过程中,如何输出切向粘性阻力系数和法向粘性阻力系数?或是如何从输出的结果中计算得到? 中说:
除此以外,是否只有像您提到的【方案二】基于wallshearstress functionObject进行修改,才能输出“绕流气泡”的粘性剪切力结果。
这个只是方便些,直接输出壁面切向和法向粘性力场,不用在paraview中写公式,但是需要改OpenFOAM代码。得先看明白wallshearstress中咋算的,主要是改calcShearStress函数,不难,慢慢看慢慢改。