CFD中文网

@RolandLent 源码在openfoam-adapter的FSI/Force.C和ForceBase.C中

十分感谢您的分享！

@Mina-Lee 在 大步长下变形网格发散 中说：

换句话说，是不是overset边界，根本不影响变形。fixed边界，才会影响变形。

确实是啊！displacementLaplacianFvMotionSolver类中组建cellDisplacement场时，只对pointDisplacement为fixedValue类型的patch进行了处理（改为cellMotion类型，并对face各点求平均），其它类型的没做任何处理；最后cellDisplacement再插值回pointDisplacement时，点位移场只有原本的fixedValue边界没被修改，其它的都被求解结果修改了。

我看代码时没深想到这一层，只在用的时候发现有时远场使用fixedValue和slip效果不同，我在这个贴子中的回复属于是瞎猫碰上死老鼠了，看来远场只要不是fixedValue，给任意都行。

最后我的问题还是用的单一区域网格以及求扩散方程变形成功实现的，为了验证还用StarCCM算了比过，吻合得很好。倒是在StarCCM中试过重叠网格，overset交界面网格难匹配，常常挖洞出错，算出来的结果锯齿也很厉害。至于rbf，有些懒得折腾了，以后遇到非用不可的情况再说。

这句话应该是您的“泣血之言”，指给后来的人看看

@Mina-Lee 在 大步长下变形网格发散 中说：

如果你要用重叠网格，记住，不要用逆距离，OF的这个方法存在BUG！

@Mina-Lee 感谢回复！不知道com版本有没有直接可用的rbf求解器呢？这个UP的视频我之前搜到过，当时还问了一下，他说dynamicOversetFvMesh就可以定义运动了。但是不知道overset的边界怎么在变形，而且他肯定是用的rbf

谢谢李老师提供的资讯，原来openfoam.org介绍slerp时说过“求解椭圆方程得位移场会导致旋转边界附近网格剪切变形”！

@李东岳 在 大步长下变形网格发散 中说：

是的，rbf对你说的正交性那个有很大的改善。这是算流固耦合？

是的，用的preCICE和calculiX算的流固耦合。

@李东岳 在 大步长下变形网格发散 中说：

你要不要玩玩extend那个，比solids4foam的简单很多，直接挂库就行

我看了一下，extend这个只针对给定旋转角度和角速度的固体运动情况，还不是普遍情况。

还有一点是，rbf插值算位移场在分核计算时要特别处理（对每一分区提取边界点，插值内部点），有些麻烦。针对我这个问题，非正交带来的误差应该还能接受，先这么算看看

如果原本不支持叠加的话，这种可能只有基于sixDoFRigidBodyMotion自定义一个动网格求解器了。

再次更新进展：
在displacementLaplacianFvMotionSolver基础上修改了一下，把边界（交界面）上固体域传来的位移分为刚性平移和局部变形，只针对后者求解laplacian方程得cellDisplacement场，插值到pointDisplacement场后再加上刚性平移场输出。这样改后，虽然时间步还是不能很大（否则端部还是出现“跟不上”的问题），但在和之前相同的时间步下，走了393个直径后的效果依旧很好：

还存在的一个问题是边界附近的网格变化没有考虑到正交（如右中图），可能是扩散方程求解变形的固有问题？也许对时间步长的限制也是因此导致的。了解到rbfMotionSolver这种求解方法，不知道是否可以解决这个问题？

@Elysiumor OpenFOAM有maxwellSlip边界条件

@bestucan 谢谢老师的手把手指导！原来是可以沿着基类往上找的，之前一看到继承了好几个基类就感觉太麻烦找不出来了。

使用grep在src中搜还可以搜到用法，也很不错，但是函数定义换行了就不好搜了。比如这个要是用 grep -rP "::average\n\(" . 啥也找不到。

编译本地doxygen也学到了，openfoam2406版本的方法整理一下

sudo apt-get install doxygen graphviz  
cd /usr/share/doc/openfoam2406-common/doc/Doxygen/
gzip -dk Doxyfile.gz  #解压doxygen配置文件
./Allwmake
firefox html/index.html

本地查看foam下的函数不会卡了，不过网页中没有的函数原型本地也没有

@bestucan 谢谢老师回复！我也知道这个搜索功能，只是我在里边实在找不到原型 ╯︿╰
再就是这种直接定义在FOAM名称空间下的函数，浏览器toolbox查看卡的不行，也很难搜索在哪儿定义的，不知道有没有好的解决方法

@Joker 不清楚算的是什么情况，但总该按实际情况给？
上下是壁面的话，就noslip（虽然两个角点有速度，就像cavity的两个角点一样）；
上下是自由流场的话，就symmetry或slip。
movingWallVelocity我理解是运动壁面的无穿透条件，上下边界没有给网格运动，只是网格节点变形，用不到这个BC吧

大家好，变形网格求解cellDisplacement时，需要用到cellMotion边界条件。我在看其updateCoeffs()函数时（见下），不理解这句：this->operator[](i) = pp[i].average(points, pointMotion); average函数咋有两个输入参数了呢？

    const fvPatch& p = this->patch();
    const polyPatch& pp = p.patch();
    const fvMesh& mesh = this->internalField().mesh();
    const pointField& points = mesh.points();
 
    word pfName = this->internalField().name();
    pfName.replace("cell", "point");
 
    const GeometricField<Type, pointPatchField, pointMesh>& pointMotion =
        this->db().objectRegistry::template
            lookupObject<GeometricField<Type, pointPatchField, pointMesh>>
            (pfName);
 
    forAll(p, i)
    {
        this->operator[](i) = pp[i].average(points, pointMotion);
    }

根据我的经验，/constant/polyMesh/boundary中左边界的patch要最后定义，不然与上下边界的角点的pointDisplacement BC会被覆盖。而且上下边界的pointDisplacement BC要用slip，再搭配合适的时间步长应该可以成功实现你要的运动。

给感兴趣的同学老师们更新一下进展：

1. directional diffusivity对我的三维情况，效果比inverDistance扩散模型更差，给的diffusivity directional (1,200,200)，只运动了不到两个直径网格便发散了
   

2. directional diffusivity好像只有给(a, b, 0)且a<<b这种模式才有效果，且仅针对单向平移的二维网格。用不同的参数试了试moving cone算例
   
   尝试把moving cone改为球，diffusivity directional (1,200,0)，效果也不错
   

3. 我的情况目前只有减小时间步长了，移动33个直径后还没发散
   

@李东岳 是扩散模型处理的：diffusivity quadratic inverseDistance (cylinder)

1楼的图给的是圆柱的对称平面（圆柱端部用半球过渡的），圆柱往x负方向动（水平往右）。到图中发散的时刻，圆柱已经走了16个直径。将时间步长减小为当前的十分之一，网格变形能跟上，但是算得太慢了。所以想能不能通过改扩散系数模型来让网格在大步长的情况下也能 “跟着动” 。

@李东岳 在 大步长下变形网格发散 中说：

要不然就是扩散模型边界条件没弄对

pointDisplacement边界条件很简单，远场是slip（或fixedValue也行，因为远场很远，10000个直径，网格扭曲也是局限在物面附近），圆柱物面是fixedValue（但实际应该是openfoam-adapter把固体域的位移map过来了：

            // Get the displacement on the patch
            fixedValuePointPatchVectorField& pointDisplacementFluidPatch(
                refCast<fixedValuePointPatchVectorField>(
                    pointDisplacement_->boundaryFieldRef()[patchID]));

            // Overwrite the nodes on the interface directly
            forAll(pointDisplacement_->boundaryFieldRef()[patchID], i)
            {
                for (unsigned int d = 0; d < dim; ++d)
                    pointDisplacementFluidPatch[i][d] = buffer[bufferIndex++];
            }

谢谢岳老师回复！刚去看了这个算例，确实使用的是directional diffusivity。顺着这个思路把其它用了directional的例程都看了一下，设置都是类似diffusivity directional (a,b,c)，我先照猫画虎试着算一个。

我用的是com版本的。重叠网格是好，可惜物面还需要根据固体域的变形来变形，而现在还没法实现不同的dynamicFvMesh同时用吧。即使能同时使用dynamicOversetFvMesh和dynamicMotionSolverFvMesh，固体也需要在刚性运动的坐标系下求解，这样得到的物面位移才是“纯变形”，感觉很难，preCICE那边也一堆想要用重叠网格算FSI的，也没见解决方案。商软（starccm）那面还是强呀，随意叠加刚性运动与变形。

至于原理，根据directionalDiffusivity.C中写的

void Foam::directionalDiffusivity::correct()
{
    const surfaceVectorField n(mesh().Sf()/mesh().magSf());
    faceDiffusivity_ == (n & cmptMultiply(diffusivityVector_, n));
}

给diffusivity directional (1,200,200)相当于给朝x方向的面1个单位的扩散系数，而给y和z向的面200个单位的扩散系数。Gpt解释说扩散系数越小代表网格越刚，变形越小。也许因此实现了在x方向跟随变形？

请教各位大神！OpenFOAM变形网格，步长大了之后，物面附近网格就跟不上了（如下图）。除了减小步长（实在是太耗费时间），还能通过使用不同diffusivity模型来改善吗，比如directional？这方面找不到啥帖子或讲解，自己测试也一头雾水，根本不知道原理上是怎么影响网格变形的。能给推荐几篇论文读读也十分感谢！

背景：
我在做一个流固耦合模拟，一根3D圆柱在横向平动时的变形，一端运动给定。因为preCICE耦合还不支持重叠网格，只能使用displacementLaplacian求解器，根据圆柱物面的平动位移+变形，将全场网格变形。试了几天后发现，OpenFOAM很难算大变形呐（虽然我这个问题主要是平动造成的大变形，可是又不能用solidBody）

或许根据对称反对称性质还可以进一步化简？OpenFOAM好像只有矢量与矢量的叉积

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

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

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

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

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

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

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

@李东岳 我用的v2006，Phi是速度势，用-grad(Phi)计算的速度和输出的U有些差异，尤其是在远场边界。

@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 ，可以：

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