CFD中文网

Haining LUO

@李东岳
是的，但是我没有在ICEM import Model里面找到相应的*.ansa格式呀，应该用ansa导出什么格式呢？

Haining LUO

工程上遇到的问题通常几何复杂，觉得ANSA这个软件对几何的处理挺强大的，个人感觉比ICEM强不少！所以想要用ANSA来清理几何、划分面网格，然后导入ICEM划分体网格（prism + tetra)

发帖的缘由：CATIA数模很多时候有缝或者有些尖角，不擅长用CATIA处理，所以想用ANSA来绕过这个坎

我是ANSA小白。当然，如果ANSA修补后的几何能够保存下来（也就是不到画面网格那一步，到处格式*.stp *.igs ?），也就满足了我的述求！望各位不吝赐教。

Haining LUO

@veen 从源码编译of是个挺痛苦的事情，建议不把of安装在root，可以考虑试一试easybuild，它是基于module的环境管理，应该比较成熟了。我从0开始，最开始在PC上面编of-3.0，到后来在实验室的机群上home编译of-4.0, foam-extend等等，多版本的管理用module还是方便些。印象中，gcc（如果你用gcc的话）和mpi的版本很重要。

Haining LUO

Thinkpad X13

Haining LUO

哈哈哈哈，好吧～

Haining LUO

@东岳麻烦请教一下，可以告知吗？

Haining LUO

有没有试过改project type？

我只做过把domain投影到database上面，在这个选项里面可以修改投影方式（版本：Pointwise V18.2 R1）

Haining LUO

想知道上图是用什么做的形状优化，是的……我截的@东岳的盆友圈

还有，这个不像是paraview的可视化的感觉，有高光！感觉不错呀

Haining LUO

是OpenFOAM吗？
如果是的话参考
https://www.cfd-online.com/Forums/openfoam-meshing/61596-cyclic-bcs-pointwiseopenfoam-export.html

Haining LUO

这个情况我倒是没有遇到过，不过我编了个synthetic的湍流生成条件，核比较少的时候patch像是一个完整的剖面，而decomposePar之后（没有遇到你这个问题）而且核数特别高（~500）之后，发现patch也被分割了的感觉。。我想说的是核数少的时候(~30)根本没发现patch上面会被不同的processor分割。

有一个计算圆形（是的，我的BC仅仅适用于圆形）中心的程序，我想不同processor算出来的ctr不同！

    // Get range and orientation
    boundBox bb(patch().patch().localPoints(), true);

    vector ctr = 0.5*(bb.max() + bb.min());

看来得定义一个常量才稳妥

Haining LUO

@xpqiu
cavity 里面是因为没有进出口所以才可以设置成 p zeroGradient ? cavity里面对p的BC都是Neuman所以必须设定pRef。（U全部是fixedValue）

我的模拟里面有入口和出口，出口可能离上游复杂流动不够远但觉得 p fixedValue uniform有点太强，换成zeroGradient，(所以p的BC全部变成了Neuman) 即使我设定了pRef =>
continuity error cannot be removed by adjusting the outflow.
please check the velocity boundary condition and/or potentialFoam to initial list the outflow.
数学上有什么解释吗？有没有其他的边界条件可选？

跟cavity里面不同的地方在我的算例里面有进口1，2，U都设定fixedValue类的BC，出口U zeroGradient。但求解U的过程没有问题（求解U用的是fvc::grad(p)），icoFoam走到压力方程才报错。按道理p方程都给了Neumann，给了pRef，为啥还是adjuctPhi error?

参考贴：http://www.cfd-china.com/topic/305/adjustphi的作用是检查边界条件

这里是第一个时间步：

Time = 1e-05

Courant Number mean: 1.07669e-05 max: 0.00554049
U_tau ajust for bulk velocity 0.3
max of yPlus = 84.6308
min of yPlus = 0.733485
U_tau ajust for bulk velocity 0.3
max of yPlus = 84.6308
min of yPlus = 0.733485
DILUPBiCCCG:  Solving for Ux, Initial residual = 1, Final residual = 2.61292e-07, No Iterations 1
DILUPBiCCCG:  Solving for Uy, Initial residual = 1, Final residual = 2.73053e-07, No Iterations 1
DILUPBiCCCG:  Solving for Uz, Initial residual = 1, Final residual = 6.39136e-07, No Iterations 1


--> FOAM FATAL ERROR: 
Continuity error cannot be removed by adjusting the outflow.
Please check the velocity boundary conditions and/or run potentialFoam to initialise the outflow.
Total flux              : 1.58773e-07
Specified mass inflow   : 2.54763e-05
Specified mass outflow  : 0
Adjustable mass outflow : 0

Haining LUO

@小葱小虫
 @Izumi 应该说的就是nNonOrthogonalCorrectors

另外你两个网格之间差别是什么，我也是用pointwise画网格，出现aspect ratio那么大可能是手误也有可能是中间有breakpoint，你用examine aspect ratio看看那些畸形的网格，用connector->distribution->delete breakpoint去掉。因为如果一个方向的网格意外的小，虽然你的求解器是simpleFoam没有CFL数，但也对稳定性有影响吧？

Haining LUO

@李东岳

我用了上面的设置，正在算，同样的网格（Hexahedra: 23万，对应上面“稀疏”的case）coupled的solver收敛迅速！虽然如你所说震荡仍然很大，但迭代次数U[x,y,z]Iters_0=1000 => U[x,y,z]Iters_0={1,2,3}。看来对于我的算例symGaussSeidel明显不如PBiCCCG，一查结果发现后者解非对称而symGaussSeidel是解对称矩阵的。请问有办法知道自己求解的问题是否是对称的吗，比如检查A矩阵UEqn.A()？

图对应的网格很稀疏（Hexahedra: 23万），实验是PIV，是中间截面的数据

导师想让我做DNS，icoFoam应该可以解DNS，不过不知道PISO的DNS能力相比pimple怎么样？

我把Re2400_mesh2_coupledSolver放在网盘了：（这是稀疏的网格，用pointwise后来发现这个网格有有些手误的地方，应该做镜像然后merge，这样就完全对称；你可以也试一试，两个solver就是两重天，我这个白痴:crying: ）
https://pan.baidu.com/s/1smtVW1f
密码：w39p

Haining LUO

@李东岳对的，实验上无量纲化了，我算的无量纲后的速度是4，实验在那个位置小于2，贴近壁面的地方<0.25

Haining LUO

@haining-luo
忘了加上：圆管截面上是O-H型的网格，能够贴体加密，算圆管流动没有问题

还有fvSolution

solvers
{
    p
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-06;
        relTol          0;
    }

    U
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-05;
        relTol          0;
    }
}

PISO
{
    nCorrectors     2;
    nNonOrthogonalCorrectors 1;
    pRefCell        0;
    pRefValue       0;
}

另外还有一个问题，因为问题是强三维的，所以N-S方程icoFoam segerated的方式（即依次求解Ux, Uy, Uz）的方式是不是不够好？(https://openfoam.org/release/2-2-0/matrix-solvers/)

我这个例子里面是否可以考虑改求解器顺序(Uy, Ux, Uz)来验证问题是不是强三维coupled，但应该怎么改呢？

Haining LUO

目标算例：基于圆管直径的雷诺数2400 (0.3 m/s * 8e-3 m/ 1e-6 m^/s)，没有湍流模型(icoFoam)
网格：结构网格，第一层网格yPlus~1（基于Re=2400边界层近似公式得到uTau [~0.02m/s]，yPlus=uTau*deltaY/nu），圆管中心的yPlus～80，径向网格数～30（是不是太稀疏？）
边界条件：
速度：两个进口（垂直）均为turbulentInlet，0.3m/s的均匀速度加上随机脉动，出口法向零梯度，壁面无滑移
压力：法向梯度为零，出口为压力uniform参照值

主要遇到两个问题：

很贴近边界[第二第三个]上的网格有较大的速度，与进口速度一个量级；靠近mixing zone能看到较强的速度场不连续性
考虑到网格可能太稀疏，有加倍网格（三个方向都几乎加倍），但无论怎么调小时间步长Ux Uy Uz 迭代(UxIter_0)第一个时间步都飙到了1000，最终残差下不了1e-2

对于网格质量(pointwise做的)，这是较密的网格，fvScheme 里面 nonOrthorgonal Corrector =1：

Time = 0

Mesh stats
    points:           779916
    faces:            2299824
    internal faces:   2262096
    cells:            760320
    faces per cell:   6
    boundary patches: 4
    point zones:      0
    face zones:       0
    cell zones:       0

Overall number of cells of each type:
    hexahedra:     760320
    prisms:        0
    wedges:        0
    pyramids:      0
    tet wedges:    0
    tetrahedra:    0
    polyhedra:     0

Checking topology...
    Boundary definition OK.
    Cell to face addressing OK.
    Point usage OK.
    Upper triangular ordering OK.
    Face vertices OK.
    Number of regions: 1 (OK).

Checking patch topology for multiply connected surfaces...
    Patch               Faces    Points   Surface topology                  
    Port1               1056     1081     ok (non-closed singly connected)  
    Port2               1056     1081     ok (non-closed singly connected)  
    Port3               1056     1081     ok (non-closed singly connected)  
    wall                34560    34631    ok (non-closed singly connected)  

Checking geometry...
    Overall domain bounding box (-0.08 -0.08 -0.004) (0.08 0.004 0.004)
    Mesh (non-empty, non-wedge) directions (1 1 1)
    Mesh (non-empty) directions (1 1 1)
    Boundary openness (8.46492e-16 1.14034e-16 -8.71928e-17) OK.
    Max cell openness = 8.80197e-16 OK.
    **Max aspect ratio = 9.65639 OK.**
    Minimum face area = 8.44179e-09. Maximum face area = 4.11244e-07.  Face area magnitudes OK.
    Min volume = 2.00787e-12. Max volume = 4.4809e-11.  Total volume = 1.18509e-05.  Cell volumes OK.
    **Mesh non-orthogonality Max: 55.6104 average: 19.4413**
    Non-orthogonality check OK.
    Face pyramids OK.
    **Max skewness = 2.05024** OK.
    Coupled point location match (average 0) OK.

Mesh OK.

End

在较疏的case，我把nu变大十倍，nu=1e-5，UxIters_0 = {1,2,3}，最终残差可以收敛到1e-5，如下图（UMean在10*(D/U)时间的截面云图），进口是左边和下面，转角处有一个速度较大的网格点（离壁面第二个网格）。我觉得收敛得快是因为粘性项大了，数值稳定性强，所以我应该在雷诺数2400的算例改改fvScheme?还是得用其他的smooth solver?个人猜想速度大的那个网格产上的可能原因是我分块的时候网格随着分的块形变（下面看得很清楚），所以会有横向的flux，横向的flux通过什么输运到下游呢？通过那个斜着45度角的flux，但是仔细一看那个地方速度大的方向不是X而是Y方向（图中只有magnitude）。那会是什么原因呢？

附上实验数据，在那个角落速度不会那么大

所以问题是：

局部速度过大是个数值上的问题吗还是网格策略不对？
大雷诺数迭代不收敛应该是数值问题，是不是该试一试迎风格式什么的？还是换个U的求解器？

附上 fvScheme

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         Gauss linear;
    grad(p)         Gauss linear;
}

divSchemes
{
    default         none;
    div(phi,U)      Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear orthogonal;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         orthogonal;
}

fluxRequired
{
    default         no;
    p               ;
}

Haining LUO

@程迪

phi是surfaceScalarField，在phi = mesh.Sf() & Uf;之前，phi还是定义在老网格的位置，通过Uf=fvc::interpolate(U)得到Uf后（我觉得是相对新网格的插值），这样就定义了一个新phi:phi = mesh.Sf() & Uf;，于是新的phi定义在新的网格面上，不再满足连续方程（旧phi是满足的），所以在进入pimple循环前通过pcorr来修正。

Haining LUO

@程迪明白你的意思了，在pimpleFoam和icoFoam里面都是在解p方程前调用了adjustPhi以保证相容性（通过修正流量保证连续方程的方式：其中涉及了两种边界条件fixedValue和inletOutletFvPatchField，修改的是inletOutletFvPatchField上的phi）。

但仔细看pimpleDyMFoam里面的correctPhi.H

if (mesh.changing()) // 如果是动网格，肯定会执行的部分
{
    forAll(U.boundaryField(), patchI)
    {
        if (U.boundaryField()[patchI].fixesValue())
        {
            U.boundaryField()[patchI].initEvaluate();
        }
    }

    forAll(U.boundaryField(), patchI)
    {
        if (U.boundaryField()[patchI].fixesValue())
        {
            U.boundaryField()[patchI].evaluate();

            phi.boundaryField()[patchI] =
                U.boundaryField()[patchI]
              & mesh.Sf().boundaryField()[patchI];
        }
    }
}

{  // 这里我想说妈蛋……这个大括号是什么鬼
    volScalarField pcorr
    (
        IOobject
        (
            "pcorr",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh,
        dimensionedScalar("pcorr", p.dimensions(), 0.0),
        pcorrTypes
    );

    dimensionedScalar rAUf("rAUf", dimTime, 1.0);  // 构建了一个单位且均匀的有时间量纲的场

    while (pimple.correctNonOrthogonal()) // 这里其实还没有进入真正的pimpleFoam循环（while (pimle.loop())），但corrNonOrtho_已经存在。但是这里没有adjustPhi !! 为什么 ?? 这里虽然求解的不是物理场p，但是pcorr也得有个边界条件吧? (于是我在createPcorrTypes.H里面找到pcorrTypes几乎是zeroGradient，如果p设定的条件中有fixedValue，对应的pcorr也是fixedValue，所以如果没有fixedValue的p边界这里也是需要adjustPhi的!)
    {
        fvScalarMatrix pcorrEqn
        (
            fvm::laplacian(rAUf, pcorr) == fvc::div(phi)  
        );

        pcorrEqn.setReference(pRefCell, pRefValue);
        pcorrEqn.solve();

        if (pimple.finalNonOrthogonalIter())
        {
            phi -= pcorrEqn.flux();
        }
    }

    #include "continuityErrs.H"
}

在这一段代码之后才是pimple.loop()，所以我觉得这是因为

        mesh.update();

        // Calculate absolute flux from the mapped surface velocity
        phi = mesh.Sf() & Uf;

用新网格算了新phi，而新的phi不是无散度的，因此用以上来修正之后再进入N-S方程求解。不过@李东岳指出的代码和我这里OpenFOAM/2.3.1-foss-2016a不大一样，我觉得这里求解的并不是Peric书方程12.13，你觉得呢？

Haining LUO

@程迪多谢分享，不过上面的帖子说的不是correctPhi.H吗？
你这里提到的是adjustPhi(phiHbyA, U, p)，我看icoFoam里面是在求解Poisson方程前，在pimpleoam和pimpleDyMFoam都是在pEqn.H里面，你提到的pRef也是，所以是为了求解Poisson方程作的准备。我觉得我们说的不是一个东西，这里是ajustPhi的帖子adjustPhi的作用是检查边界条件？

Haining LUO

@金石为开我觉得这里是因为网格动到了新的位置，而当前的phi相对于当前的网格就不再divergence free，网格更新了而phi场也需要更新，在当前网格下对absolute phi进行修正，而不是对应的方程12.13。假设了一个pcorrEqn场的拉普拉斯等于当前phi的不为零散度，然后用phi -= pcorrEqn.flux();修正。

CFD中文网

Haining LUO

帖子