钝体建筑扰流的大涡模拟

HITSC30

老师们朋友们大家好，请问有没有老师或者朋友一样是做高层建筑大涡模拟风荷载的。目前遇到一个问题，在运用fluent时我采用的网格边界层厚度足够小，满足y+<1，保证可以用近壁模型计算。但是这样的网格在FOAM中对应的是运行发散。排查后发现边界层网格在不设置，或者设置的非常厚的情况下运行才会满足。想请问一下各位老师是否碰见过类似的问题？
再者由于fixed value条件一般会导致脉动风压误差变大，所以想在FOAM中运用outflow出口的时候，在foam里面一设置就会出现错误。让我用potentialfoam初始化，但是该种情况下，最终的结果也会发散。有朋友清楚该种问题是什么情况吗？
谢谢大家。

李东岳

哪个求解器？用的什么湍流模型。网格是六面体么

HITSC30

@李东岳李老师您好，采用的是pisofoam算法，求解器p用的是GAMG，U用的是PBICG。湍流模型LES-WALE，使用的是fluentmeshing生成的多面体网格。该网格在fluent中已经运算过，得到的结果很不错。来到FOAM罢工了。个人估计是边界层的问题，模型尺寸为0.1m0.1m0.1m，整个计算域尺寸为8m4m1.6m。我的建筑面面网格的尺寸是0.002m，为了满足y+<1，建筑表面的边界层网格厚度是0.00005m。地面厚度的全局最大尺寸为0.05m，地面边界层网格的厚度是0.0005m。也曾在一篇帖子下发现您推荐的snappyHexMesh，但是生成质量不佳，有时候甚至边界层网格无法生成。目前用OpenFOAM做LES建筑风荷载的文献也没怎么查到，导致困在这了

李东岳

log贴一下，按照这种代码的方式

log

zhouxu

@李东岳速度出口是零梯度边界，压力出口是固定值，速度入口是固定值，压力入口是零梯度

gtian

@hitsc30 什么类型的建筑物？
https://link.springer.com/article/10.1007/s10546-020-00574-1
这是我之前做的，openfoam + 小方块。有什么需要可以邮件联系。tiang@sustech.edu.cn

HITSC30

@李东岳

\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volVectorField;
    location    "1";
    object      U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0 0 0);

boundaryField
{

    inlet
    {
        type            timeVaryingMappedFixedValue;
        offset          (0 0 0);
        setAverage      off;
    }


    outlet
    {
        type            zeroGradient;
    }

    upperwall
    {
        type            symmetry;
    }

    bottom
    {
        type            noSlip; 
    }
    
    symmetry
    {
        type            symmetry;
    }
    
    building
    {
        type            noSlip; 
    }   	
	
}

// ************************************************************************* //
/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "1";
    object      p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 0;

boundaryField
{
   
    inlet
    {
        type            zeroGradient;
    }

    outlet
    {
        type            fixedValue;
        value           uniform 0;
        //type            zeroGradient;
    }

    upperwall
    {
        type            symmetry; 
    }

    bottom
    {
        type            zeroGradient;
    }
    
    symmetry
    {
        type            symmetry;
    }

    building
    {
        type            zeroGradient;
    }

}

// ************************************************************************* //
/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "1";
    object      nut;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 2 -1 0 0 0 0];

internalField   uniform 0;

boundaryField
{
   
    inlet
    {
        type            calculated;
        value           uniform 0;
    }

    outlet
    {
        type            calculated;
        value           uniform 0;
    }

    upperwall
    {
        type            symmetry; 
    }

    bottom
    {
        type            nutUSpaldingWallFunction;
        value           $internalField;
        //type            fixedValue;
        //value           uniform 0;
    }
    
    symmetry
    {
        type            symmetry;
    }

    building
    {
        type            nutUSpaldingWallFunction;
        value           $internalField;
        //type            fixedValue;
        //value           uniform 0;
    }

}

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

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      turbulenceProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

simulationType LES;

LES
{
    LESModel        WALE;

    turbulence      on;

    printCoeffs     on;

    delta           cubeRootVol;

    cubeRootVolCoeffs
    {
        deltaCoeff      1;
    }

    PrandtlCoeffs
    {
        delta           cubeRootVol;
        cubeRootVolCoeffs
        {
            deltaCoeff      1;
        }

        smoothCoeffs
        {
            delta           cubeRootVol;
            cubeRootVolCoeffs
            {
                deltaCoeff      1;
            }

            maxDeltaRatio   1.1;
        }

        Cdelta          0.158;
    }

    vanDriestCoeffs
    {
        delta           cubeRootVol;
        cubeRootVolCoeffs
        {
            deltaCoeff      1;
        }

        smoothCoeffs
        {
            delta           cubeRootVol;
            cubeRootVolCoeffs
            {
                deltaCoeff      1;
            }

            maxDeltaRatio   1.1;
        }

        Aplus           26;
        Cdelta          0.158;
    }

    smoothCoeffs
    {
        delta           cubeRootVol;
        cubeRootVolCoeffs
        {
            deltaCoeff      1;
        }

        maxDeltaRatio   1.1;
    }
}

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

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  4.x                                   |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
//shijianyijiedushuxiang
ddtSchemes  //
{
    default backward;
}
//tidu
gradSchemes
{
    default         Gauss linear;

    //grad(nuTilda)   cellLimited Gauss linear 1;
    //grad(U)         cellLimited Gauss linear 1;
    //grad(P)         Gauss linear;
}
//sandduxiang
divSchemes
/////////////////////////////////////////////////////
{
    default         none;

    div(phi,U)      bounded Gauss linearUpwind grad(U);
    div(phi,k)      bounded Gauss linearUpwind grad(k);
    div(phi,B)      bounded Gauss linearUpwind grad(B);
    div(phi,omega)      bounded Gauss linearUpwind grad(omega);
    div(phi,B)      bounded Gauss linear;
    //div(phi,nuTilda) Gauss limitedLinear 1;
    div((nuEff*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes  //
{
    default         Gauss linear limited corrected 0.33;
}

interpolationSchemes   //
{
    default         linear;
}

snGradSchemes  //
{
    default         limited corrected 0.33;
}

wallDist
{
    method meshWave;
}

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

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  4.x                                   |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
    p
    {
        solver          GAMG;
        tolerance       1e-6;
        relTol          0.1;
        smoother        GaussSeidel;
        //nCellsInCoarsestLevel 50;
    }

    pFinal
    {
        $p;
        tolerance       1e-6;
        relTol          0;
    }

    U
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-6;
        relTol          0.05;
    }
    UFinal
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-6;
        relTol          0;
    }
}

PISO
{
    nCorrectors     3;
    nNonOrthogonalCorrectors 1;
    //pRefCell        0;
    pRefPoint       (2 0 1.5);
    pRefValue       0;
}

relaxationFactors
{
    fields
    {
        p           0.3;

    }
    
    equations
    {
        U           1;
    }





    //"U.*"               0.5;
    //"nuTilda.*"         1;
}


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

李老师这是我的设置，感谢各位老师朋友的解答

HITSC30

@zhouxu 对的老师，对p用零梯度一直发散，potentialfoam后也存在这种情况，后面就用了固定值。不过看文献说p采用固定值会增大脉动风压的误差。

HITSC30

@gtian 谢谢！我第一个case也是规整的方形建筑物，采用的大涡模拟计算的。目前在看如何采用商软生成多面体网格，保证能有一个好的收敛结果。fluent用的网格在foam罢工了

gtian

@hitsc30

snappyHexMeshDict 中把snap 和 add 关掉，就会生成结构化的六面体网格，网格质量会号很多，只需要有stl文件即可.

// Which of the steps to run
castellatedMesh true;
snap            false;
addLayers       false;
geometry
{
    flange11.stl
    {
        type triSurfaceMesh;
        name flange11;
    }

stl文件如下：

 solid flange
  facet normal 0 0 -1
    outer loop
      vertex 0.170 0.000 0.230
      vertex 0.170 0.020 0.230
      vertex 0.190 0.020 0.230
    endloop
  endfacet
  facet normal 0 0 -1
    outer loop
      vertex 0.170 0.000 0.230
      vertex 0.190 0.020 0.230
      vertex 0.190 0.000 0.230
    endloop
  endfacet

  facet normal -1 0 0
    outer loop
      vertex 0.170 0.000 0.230
      vertex 0.170 0.000 0.210
      vertex 0.170 0.020 0.210
    endloop
  endfacet
  facet normal -1 0 0
    outer loop
      vertex 0.170 0.000 0.230
      vertex 0.170 0.020 0.210
      vertex 0.170 0.020 0.230
    endloop
  endfacet
  facet normal 0 -1 0
    outer loop
      vertex 0.170 0.000 0.230
      vertex 0.190 0.000 0.210
      vertex 0.190 0.000 0.230
    endloop
  endfacet
  facet normal 0 -1 0
    outer loop
      vertex 0.170 0.000 0.230
      vertex 0.190 0.000 0.210
      vertex 0.170 0.000 0.210
    endloop
  endfacet
 
  facet normal 0 0 1
    outer loop
      vertex 0.170 0.000 0.210
      vertex 0.190 0.000 0.210
      vertex 0.190 0.020 0.210
    endloop
  endfacet
  facet normal 0 0 1
    outer loop
      vertex 0.170 0.000 0.210
      vertex 0.190 0.020 0.210
      vertex 0.170 0.020 0.210
    endloop
  endfacet
  facet normal 1 0 0
    outer loop
      vertex 0.190 0.000 0.230
      vertex 0.190 0.020 0.230
      vertex 0.190 0.020 0.210
    endloop
  endfacet  
  facet normal 1 0 0
    outer loop
      vertex 0.190 0.000 0.230
      vertex 0.190 0.020 0.210
      vertex 0.190 0.000 0.210
    endloop
  endfacet
  facet normal 0 1 0
    outer loop
      vertex 0.170 0.020 0.230
      vertex 0.170 0.020 0.210
      vertex 0.190 0.020 0.210
    endloop
  endfacet
  facet normal 0 1 0
    outer loop
      vertex 0.170 0.020 0.230
      vertex 0.190 0.020 0.230
      vertex 0.190 0.020 0.210
    endloop
 endfacet 
 endsolid flange

HITSC30

@gtian 感谢回复。我模拟的是建筑风压，感觉将边界层去掉的情况下，有点难搞
非常感谢您的帮助

gtian

@hitsc30 如果你提到的边界层指的是wall周围的对数区以下的粘性子层和buffer layer的话。可以有两种方式，一种是wall model LES，通过壁面模型来调nut，这时第一层网格在y+=30 即可。另一种是wall resolved LES，这种情况只能加密网格。blockmesh和shm都可以加密网格，blockmesh可以控制aspect ratio来加密，shm是4等分加密。这种情况下如果要resolve到粘性子层（viscous sublayer），需要在y+<5 以内至少有三个点，也就是第一点y+<1.

cccrrryyy

@hitsc30 你这个涉及的东西有点多，排查起来会比较复杂。首先建议抛开fluent里的设置，因为大涡模拟这块openfoam和fluent差的还有点多。

首先是网格，snappy如果不好生成的话，试试cfmesh，这个比snappy好用，比fluentmeshing我感觉也就是边界层生成差一些。考虑到你这个几何相对应该比较规整，cfmesh应该可以的。

其次是入口条件。openfoam大涡模拟的入口条件这块是不如fluent的vortex method的，fluent的很好用，不怎么发散，生成的湍流也很强。我看你用的是timeVaryingMappedFixedValue，这是用的实验数据或者测试数据吗？你fluent是怎么用的？

然后是nut的边界条件，我之前也是发现nutUSpaldingWallFunction很好用，但我不确定它应该和WALE一起。WALE一般是传热领域用的多，原因是这个模型是不需要壁面函数的，而传热主要就是要解析壁面，不能模化。WALE感觉和你加的nut边界条件可能会冲突？

再然后是数值格式，一个是pisoFoam，openfoam里面的pisoFoam外迭代只有一步，可能收敛不太好。试试pimpleFoam加上外迭代吧，这个其实你在fluent里面应该也设置过。或者把你目前计算的残差曲线贴出来看看。

最后，为了排查这些问题，建议先跑个非稳态的RANS试试，不然理不出头绪来。其实最好先跑个稳态RANS，没问题了再非稳态RANS，再没问题了再上LES。

HITSC30

@gtian 在钝体建筑扰流的大涡模拟中说：

@hitsc30 如果你提到的边界层指的是wall周围的对数区以下的粘性子层和buffer layer的话。可以有两种方式，一种是wall model LES，通过壁面模型来调nut，这时第一层网格在y+=30 即可。另一种是wall resolved LES，这种情况只能加密网格。blockmesh和shm都可以加密网格，blockmesh可以控制aspect ratio来加密，shm是4等分加密。这种情况下如果要resolve到粘性子层（viscous sublayer），需要在y+<5 以内至少有三个点，也就是第一点y+<1.
感谢回复，这个确实可以在后续重点突击一下。此前我曾用blockmesh做网格背景、snappyhexmesh加密做过一个网格。当时我得到的计算值也和我用fluentmeshing商软所做的商业网格差不太多，没太考虑网格的质量。您提出的这个是我当前需要考虑的一个方向。我当前觉得目前比较影响我接过的可能是离散格式，在这方面的知识储备非常的欠缺，只能选一些高阶离散格式，可能对结果有不小的影响。

HITSC30

@cccrrryyy 您好，感谢回复。目前接触到的在openfoam中做建筑风荷载模拟的论文主要是tamura教授在北京交通大学和重庆大学领导的学科组推出来的，都采用的是pisofoam算法。湍流风入口是由CDRFG方法生成的，在中文网有个帖子谈论过这种植入。也咨询过北交的师兄，他们植入的CDRFG算法也是类似的植入，该方法应该是没有问题的。
nut函数，我也和您的观点一样。该函数能自动切换成粘性层的线性表达式和对数表达式。再者是您提出的WALE，这个确实是需要改动一下.目前接触的论文中，foam中都是用的经典的smagorinksy。fluent中，国内硕士很多用的是WALE结果对的确实不是很好。您提到的WALE这个模型的具体作用，对代码认识比较浅，确实之前没考虑过，确实是我需要再次调整的方向。
当前的两个重点是离散格式和亚格子模型，然后是上面的师兄说的网格的质量。
ps：感觉foam计算y+很奇怪，我的边界层网格每次改变的不少，感觉y+变化趋势和我想的南辕北辙。

李东岳

你可以把你的网格截图一下看看能差到什么样子

gtian

@hitsc30
y+ 的理解不复杂，看code需要一点点的耐心就可以。我给你贴一下yPlusLES。

 if (isA<wallFvPatch>(currPatch))
            {
                yPlus.boundaryField()[patchi] =
                    d[patchi]
                   *sqrt
                    (
                        nuEff.boundaryField()[patchi]
                       *mag(U.boundaryField()[patchi].snGrad())
                    )
                   /nuLam.boundaryField()[patchi];
                const scalarField& Yp = yPlus.boundaryField()[patchi];

                Info<< "Patch " << patchi
                    << " named " << currPatch.name()
                    << " y+ : min: " << gMin(Yp) << " max: " << gMax(Yp)
                    << " average: " << gAverage(Yp) << nl << endl;
            }

![公式]( 1)
![公式2]( 2)
![替代文字]( 3)

我写的，你可以输出来看看每一个量的大小

            if (isA<wallFvPatch>(currPatch))
            {
                frictionvelocityLES.boundaryField()[patchi] = 
                   sqrt
                    (
                        nuEff.boundaryField()[patchi]
                       *mag(U.boundaryField()[patchi].snGrad())
                    );
                const scalarField& Yp = frictionvelocityLES.boundaryField()[patchi];

                Info<< "Patch " << patchi
                    << " named " << currPatch.name()
                    << " friction velocity u* : min: " << gMin(Yp) << " max: " << gMax(Yp)
                    << " average: " << gAverage(Yp) << nl << endl;

                 Info<< "Patch " << patchi
                    << " named " << currPatch.name()
                    << "nuLam : min: " << gMin(nuLam) << " max: " << gMax(nuLam)
                    << " average: " << gAverage(nuLam) << nl << endl;

                Info<< "Patch " << patchi
                    << " named " << currPatch.name()
                    << "d : min: " << gMin(d) << " max: " << gMax(d)
                    << " average: " << gAverage(d) << nl << endl;
            }

近壁面第一层网格的话，d= 0.5* 第一个cell的长度. du/dy = （u_infty - 0）/ (d-0)
nuefff = nu （层流 laminar） + nut（wall 上 nut =0 或者是你给的壁面模型的值）

我给的nut在壁面上就是0.

    object      nuSgs;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [ 0 2 -1 0 0 0 0 ];

internalField   uniform 0;

boundaryField
{
    bottomWall
    {
        type            fixedValue;
        value           uniform 0;
    }

HITSC30

@gtian 在钝体建筑扰流的大涡模拟中说：

@hitsc30
y+ 的理解不复杂，看code需要一点点的耐心就可以。我给你贴一下yPlusLES。

 if (isA<wallFvPatch>(currPatch))
            {
                yPlus.boundaryField()[patchi] =
                    d[patchi]
                   *sqrt
                    (
                        nuEff.boundaryField()[patchi]
                       *mag(U.boundaryField()[patchi].snGrad())
                    )
                   /nuLam.boundaryField()[patchi];
                const scalarField& Yp = yPlus.boundaryField()[patchi];

                Info<< "Patch " << patchi
                    << " named " << currPatch.name()
                    << " y+ : min: " << gMin(Yp) << " max: " << gMax(Yp)
                    << " average: " << gAverage(Yp) << nl << endl;
            }

![公式]( 1)
![公式2]( 2)
![替代文字]( 3)

我写的，你可以输出来看看每一个量的大小

            if (isA<wallFvPatch>(currPatch))
            {
                frictionvelocityLES.boundaryField()[patchi] = 
                   sqrt
                    (
                        nuEff.boundaryField()[patchi]
                       *mag(U.boundaryField()[patchi].snGrad())
                    );
                const scalarField& Yp = frictionvelocityLES.boundaryField()[patchi];

                Info<< "Patch " << patchi
                    << " named " << currPatch.name()
                    << " friction velocity u* : min: " << gMin(Yp) << " max: " << gMax(Yp)
                    << " average: " << gAverage(Yp) << nl << endl;

                 Info<< "Patch " << patchi
                    << " named " << currPatch.name()
                    << "nuLam : min: " << gMin(nuLam) << " max: " << gMax(nuLam)
                    << " average: " << gAverage(nuLam) << nl << endl;

                Info<< "Patch " << patchi
                    << " named " << currPatch.name()
                    << "d : min: " << gMin(d) << " max: " << gMax(d)
                    << " average: " << gAverage(d) << nl << endl;
            }

近壁面第一层网格的话，d= 0.5* 第一个cell的长度. du/dy = （u_infty - 0）/ (d-0)
nuefff = nu （层流 laminar） + nut（wall 上 nut =0 或者是你给的壁面模型的值）

我给的nut在壁面上就是0.

    object      nuSgs;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [ 0 2 -1 0 0 0 0 ];

internalField   uniform 0;

boundaryField
{
    bottomWall
    {
        type            fixedValue;
        value           uniform 0;
    }

非常感谢！

HITSC30

@李东岳 544c4ce2-63b4-47d3-bdbc-87861c0e7d5c-1651239158(1).png 抱歉，李博。这几天电脑卡的要死，开启fluentmeshing就卡死，一直想着要回复您。但是也没法打开软件。这是我再fluentmeshing中检测的问题。我当前的问题主要是，我画过粗略的多面体和精细的多面体网格是否对于结果优化都没太大的影响了。当前想着试一试snappyhexmesh画六面体拯救一下我，在使用中发现其对于边界层网格生成有点随缘化，网格在90度拐角处处理的不是很好。

李东岳

@hitsc30 你这个网格用openfoam跑没问题啊，为啥要用snappyHexMesh呢。你是说网格粗细对结果没啥影响？这是有可能的。你可以发线图看一下对比一下结果

CFD中文网

钝体建筑扰流的大涡模拟