OpenFOAM计算圆球绕流过程中，如何输出切向粘性阻力系数和法向粘性阻力系数？或是如何从输出的结果中计算得到？

Prometheus10

各位老师好！

最近计算圆球绕流，通过在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文件还不太看得明白，希望能有老师指点，谢谢！

李东岳

阻力系数计算在这里面有 http://www.dyfluid.com/theory.pdf 在OpenFOAM的计算里面，没看到关于pi的相关操作

Prometheus10

@李东岳好的，谢谢东岳老师。
我也自己看了下[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等：

东岳老师，各位大佬，能否指点一下这几个变量通常要去如何理解它的物理意义？
谢谢！

lwjetmann

@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进行修改，输出投影后的结果

Prometheus10

@lwjetmann 感谢老师，组里同学大家都不太会OpenFOAM，我马上也看一下。再次感谢。

Prometheus10

@lwjetmann 老师您好，再请教一下您。

1. 用wallshearstress后处理算壁面粘性剪切力，然后在paraview中计算出在壁面各点法向和切向的投影；
2. 基于wallshearstress functionObject进行修改，输出投影后的结果

按照您提到的上述方案1，对于圆球绕流（球体界面是“noSlip”边界条件），我顺利的输出了粘性剪切力，结果也验证的上。

想请您再帮忙看看另外一种情况：绕流气泡，即将边界条件改为Symmetry：
因为wallshearstress只对wall类型的边界计算粘性剪应力，因此不能再通过【方案1】输出界面上的粘性剪应力结果。
关于设置Symmetry边界条件的原因是源于这里：Stokes流中的圆球绕流的阻力系数，怎么设置才能算准。

对此，针对绕流气泡的结果（边界条件为Symmetry）：

我尝试手动在结果中将Symmetry边界条件改为Wall，然后执行一次wallshearstress；
通过对比文献中气泡的压差阻力系数、法向粘性应力和切向粘性应力结果，输出的结果并不对

除此以外，是否只有像您提到的【方案二】基于wallshearstress functionObject进行修改，才能输出“绕流气泡”的粘性剪切力结果。
如果您方便的话，能否请您指点一些修改思路，谢谢老师。

lwjetmann

@Prometheus10 在 OpenFOAM计算圆球绕流过程中，如何输出切向粘性阻力系数和法向粘性阻力系数？或是如何从输出的结果中计算得到？中说：

我尝试手动在结果中将Symmetry边界条件改为Wall，然后执行一次wallshearstress；
通过对比文献中气泡的压差阻力系数、法向粘性应力和切向粘性应力结果，输出的结果并不对

是的，因为wallshearstress代码会自动过滤非wall边界，可以手动改为wall。既然计算剪切力的公式都是一样的，这么改一下不会影响计算正确性的。和文献对不上是不是symmetry边界条件的问题呢？（你们这个帖子挺长，我看大意是不是：楔形网格算球noslip正确，但算slip和symmetry和文献对不上，且symmetry效果更好些？）
用整球网格算会如何呢？做一些其它情况的验证呢，比如壁面滑移的方腔流动？

我之前做过低雷诺数下的一阶maxwell速度滑移的圆柱绕流，算出来的受力以及粘性正应力占比和文献比得上。不过你的自由滑移情况我没考虑过，不知道为啥出问题。

@Prometheus10 在 OpenFOAM计算圆球绕流过程中，如何输出切向粘性阻力系数和法向粘性阻力系数？或是如何从输出的结果中计算得到？中说：

除此以外，是否只有像您提到的【方案二】基于wallshearstress functionObject进行修改，才能输出“绕流气泡”的粘性剪切力结果。

这个只是方便些，直接输出壁面切向和法向粘性力场，不用在paraview中写公式，但是需要改OpenFOAM代码。得先看明白wallshearstress中咋算的，主要是改calcShearStress函数，不难，慢慢看慢慢改。

CFD中文网

OpenFOAM计算圆球绕流过程中，如何输出切向粘性阻力系数和法向粘性阻力系数？或是如何从输出的结果中计算得到？