chtMultiRegionFoam求解器及算例分析

最新推荐文章于 2024-06-19 21:06:31 发布

CFD_Tyro

最新推荐文章于 2024-06-19 21:06:31 发布

阅读量3k

点赞数 10

分类专栏：理论基础文章标签：算法

本文链接：https://blog.csdn.net/CFD_Tyro/article/details/126347830

版权

本文详细分析了chtMultiRegionFoam求解器的算例结构与运行步骤，包括blockMeshDict、snappyHexMeshDict和splitMeshRegions工具的使用。此外，介绍了求解器的工作原理，解析了固体和流体的求解过程，展示了流体方程的求解策略。最后，探讨了求解器中涉及的关键算法，如压力修正和能量方程的解算。

摘要由CSDN通过智能技术生成

1. 算例分析

1.1. 算例结构

算例目录heatTransfer/chtMultiRegionFoam/heatedDuct

0
- fluid
  - p
  - p_rgh
  - T
  - U
- heater
  - T
- metal
  - T
constant
- fluid
  - g
  - momentumTransport
  - physicalProperties
- heater
  - fvModels
  - physicalProperties
- metal
  - physicalProperties
- geometry
  - fluidToMetal.stl
  - heatedDuct.stl
  - metalToHeater.stl
- regionProperties
system
- fluid
  - fvSchemes
  - fvSolution
- heater
  - fvSchemes
  - fvSolution
- metal
  - fvSchemes
  - fvSolution
- blockMeshDict
- controlDict
- decomposeParDict
- fvSchemes
- fvSolution
- meshQualityDict
- snappyHexMeshDict

1.2. 算例运行

通过Allrun了解其执行步骤：

blockMesh
snappyHexMesh -overwrite
splitMeshRegions -cellZones -overwrite
decomposePar -allRegions
chtMultiRegionFoam
reconstructPar -allRegions

1.2.1. blockMeshDict

该字典文件代码段。注意到blocks段落中出现了fluid关键词，这是icoFoam算例中没有的。通过观察blockMesh的运行输出文件，该关键词的效果是生成了名为fluid的cellZone，其包含 $20 * 24 * 60$ 个网格

convertToMeters 0.001;

vertices
(
    (-0.001 -10.001  -0.001) // (x0, y0, z0)
    (50.001 -10.001  -0.001) // (x1, y0, z0)
    (50.001  50.001  -0.001) // (x1, y1, z0)
    (-0.001  50.001  -0.001) // (x0, y1, z0)
    (-0.001 -10.001 150.001) // (x0, y0, z1)
    (50.001 -10.001 150.001) // (x1, y0, z1)
    (50.001  50.001 150.001) // (x1, y1, z1)
    (-0.001  50.001 150.001) // (x0, y1, z1)
);

blocks
(
    hex (0 1 2 3 4 5 6 7) fluid (20 24 60) simpleGrading (1 1 1)
);

defaultPatch
{
   
    type patch; // 定义默认面的类型为patch
}

运行该文件会在constant/下生成polyMesh/，用于记录网格信息。

1.2.2. snappyHexMeshDict

该工具可以自动的从STL和OBJ文件生成六面体和多面体网格，网格依靠迭代将一个初始网格细化，并将细化后的网格变形以依附于表面。运行该工具需提供

STL文件格式的几何文件（binary或者ascii格式）存储在constant/triSurface子字典中
一个背景六面体网格（必须是纯六面体），定义计算域范围和背景网格，一般由 blockMesh 生成。stl文件必须和背景网格产生交叉，snappyHexMesh就是将这些交叉处的网格细化
snappyHexMeshDict 字典，它提供了生成网格的必要信息，位于 system 子文件夹下

#includeEtc "caseDicts/mesh/generation/snappyHexMeshDict.cfg"

castellatedMesh true; // 是否切割网格
snap            true; // 是否进行网格贴合或对齐
addLayers       false; // 是否添加边界层

geometry // 所用的几何文件
{
   
    heatedDuct // 整个构体几何外型
    {
   
        type triSurfaceMesh; // 类型，除此之外还有searchableBox
        file "heatedDuct.stl"; // 文件名

        regions // 当stl文件内包含多个solid时需指定
        {
   
            metalInlet     {
    name metalInlet;     } // stl文件内的域名
            heaterInlet    {
    name heaterInlet;    }
            fluidInlet     {
    name fluidInlet;     }
            metalOutlet    {
    name metalOutlet;    }
            heaterOutlet   {
    name heaterOutlet;   }
            fluidOutlet    {
    name fluidOutlet;    }
            metalExternal  {
    name metalExternal;  }
            heaterExternal {
    name heaterExternal; }
        }
    }

    metalToHeater // 金属套管和加热器的接触面
    {
   
        type triSurfaceMesh;
        file "metalToHeater.stl";
    }

    fluidToMetal // 流体和金属套管接触面
    {
   
        type triSurfaceMesh;
        file "fluidToMetal.stl";
    }
};

castellatedMeshControls // 切割或细化网格的字典
{
   
    refinementSurfaces // 指定需要细化的面，相应还有refinementRegions
    {
   
        heatedDuct
        {
   
            level (0 0); // 该面默认的细化程度：细化等级的最小、最大值
            regions // 除了默认，指定以下region的细化方式以及面的类型
            {
   
                metalInlet     {
    level (0 0); patchInfo {
    type patch; } }
                heaterInlet    {
    level (0 0); patchInfo {
    type patch; } }
                fluidInlet     {
    level (0 0); patchInfo {
    type patch; } }
                metalOutlet    {
    level (0 0); patchInfo {
    type patch; } }
                heaterOutlet   {
    level (0 0); patchInfo {
    type patch; } }
                fluidOutlet    {
    level (0 0); patchInfo {
    type patch; } }
                metalExternal  {
    level (1 1); patchInfo {
    type wall;  } }
                heaterExternal {
    level (1 1); patchInfo {
    type wall;  } }
            }
        }

        fluidToMetal
        {
   
            level (1 1);
            faceZone fluidToMetal; // 生成faceZone的名称
            cellZone metal; // 生成cellZone的名称
            mode        insidePoint; // 如何生成cellZone
            insidePoint (0.025 0.0025 0.075);
        }

        metalToHeater
        {
   
            level (1 1);
            faceZone metalToHeater;
            cellZone heater;
            mode        insidePoint;
            insidePoint (0.025 -0.005 0.075);
        }
    }

    nCellsBetweenLevels 1;

    insidePoint (0.025 0.025 0.075);
}

addLayersControls
{
   
    relativeSizes       true; // 相对边界层厚度
    minThickness        1; // 边界层网格最小厚度
    finalLayerThickness 1; // 壁面最近的边界层厚度
    expansionRatio      1; // 膨胀因子
    layers
    {
   }
}

执行该命令后，

在polyMesh/faceZones中生成了fluidToMetal和metalToHeater，在polyMesh/cellZones中生成了fluid，metal和heater
metalExternal、heaterExternal、fluidToMetal、metalToHeater的网格被细化

1.2.3. splitMeshRegions

该工具将网格划分为多个区域，可以将已有的多个cellZones转化为一个cellZone

splitMeshRegions -cellZones

Each region is defined as a domain whose cells can all be reached by cell-face-cell walking without crossing

boundary faces
additional faces from faceset (-blockedFaces faceSet).
any face between differing cellZones (-cellZones)

可见，上述命令将根据已有的cellZone将创建region，该region中的网格彼此相连，直到接触到另一个cellZone为止。

该工具运行时会先检验system/xx/fvScheme文件是否存在，如果存在就会创建名为xx的region，生成如下文件

0
- cellToRegion
- heater
  - cellToRegion 该文件给出当前region的编号以及region的边界条件，还有region之间交界面的信息，如metal_to_fluid
- metal
  - cellToRegion
- fluid
  - cellToRegion
constant
- cellToRegion
- heater
  - polyMesh
- metal
  - polyMesh
- fluid
  - polyMesh

2. chtMultiRegionFoam求解器解析

该求解器采用了分离式的求解策略，即上一个方程的求解结果作为输入条件应用于下一个方程的求解。在共轭传热问题中，求解步骤为

将初始固体温度作为边界条件输入到流体温度方程中求解；
解出流体温度，作为边界条件输入到固体温度方程中；
求出新的固体温度，重复上面两步，直到收敛。

压力基求解器的释义：The solver used to solve the fluid equations is a pressure bases solver. That means that a pressure equation (similar to the pressure equation used in an incompressible solver) is used to establish the connection between the momentum and the continuity equation.

2.1. 程序结构

2.1.1. chtMultiRegionFoam.C

int main(int argc, char *argv[])
{
   
    #define NO_CONTROL
    #define CREATE_MESH createMeshesPostProcess.H
    #include "postProcess.H"

    #include "setRootCaseLists.H"
    #include "createTime.H"
    #include "createMeshes.H"
    pimpleMultiRegionControl pimples(fluidRegions, solidRegions);
    #include "createFields.H"
    #include "initContinuityErrs.H"
    #include "createFluidPressureControls.H"
    #include "createTimeControls.H"
    #include "readSolidTimeControls.H"
    #include "compressibleMultiRegionCourantNo.H"
    #include "solidRegionDiffusionNo.H"
    #include "setInitialMultiRegionDeltaT.H"

    while (pimples.run(runTime))
    {
   
        #include "readTimeControls.H"
        #include "readSolidTimeControls.H"

        #include "compressibleMultiRegionCourantNo.H"
        #include "solidRegionDiffusionNo.H"
        #include "setMultiRegionDeltaT.H"

        runTime++;

        Info<< "Time = " << runTime.userTimeName() << nl << endl;

        // Optional number of energy correctors
        // 若没有在fvSolution中找到该项，则设为1，本算例中没有该项
        const int nEcorr = pimples.dict().lookupOrDefault<int>
        (
            "nEcorrectors",
            1
        );

        // --- PIMPLE loop
        while (pimples.loop())
        {
   
            List<tmp<fvVectorMatrix>> UEqns(fluidRegions.size()); // 建立一个列表，里面每个元素都是fvVectorMatrix

            for(int Ecorr=0; Ecorr<nEcorr; Ecorr++) // nEcorr=1
            {
   
                // solidRegions & fluidRegions 存在于constant/regionProperties
                forAll(solidRegions, i)
                {
   
                    Info<< "\nSolving for solid region "
                        << solidRegions[i].name() << endl;
                    #include "setRegionSolidFields.H" // volScalarField& e = thermo.he()
                    #include "solveSolid.H"
                }
                forAll(fluidRegions, i)
                {
   
                    Info<< "\nSolving for fluid region "
                        << fluidRegions[i].name() << endl;
                    #include "setRegionFluidFields.H"
                    #include "solveFluid.H"
                }
            }
        }

        runTime.write();

        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
    }

    Info<< "End\n" << endl;

    return 0;
}

2.1.2. createMeshes.H

regionProperties rp(runTime);

#include "createFluidMeshes.H"
#include "createSolidMeshes.H"

2.1.2.1. createFluidMeshes.H

const wordList fluidNames
(
    rp.found("fluid") ? rp["fluid"] : wordList(0)
);

PtrList<fvMesh> fluidRegions(fluidNames.size()); // 是一个指针列表，每一个指针所指向的元素都是fvMesh类型

forAll(fluidNames, i)
{
   
    Info<< "Create fluid mesh for region " << fluidNames[i]
        << " for time = " << runTime.timeName() << nl << endl;

    fluidRegions.set
    (
        i,
        new fvMesh
        (
            IOobject
            (
                fluidNames[i],
                runTime.timeName(),
                runTime,
                IOobject::MUST_READ
            )
        )
    );
}

2.1.2.2. createSolidMeshes.H

const wordList solidNames
(
    rp.found("solid") ? rp["solid"] : wordList(0)
);

PtrList<fvMesh> solidRegions(solidNames.size());

forAll(solidNames, i)
{
   
    Info<< "Create solid mesh for region " << solidNames[i]
        << " for time = " << runTime.timeName() << nl << endl;

    solidRegions.set
    (
        i,
        new fvMesh
        (
            IOobject
            (
                solidNames[i],
                runTime.timeName(),
                runTime,
                IOobject::MUST_READ
            )
        )
    );

    // Force calculation of geometric properties to prevent it being done
    // later in e.g. some boundary evaluation
    //(void)solidRegions[i].weights();
    //(void)solidRegions[i].deltaCoeffs();
}

2.1.3. createFields.H

#include "createFluidFields.H"
#include "createSolidFields.H"

从createFluidFields.H和log文件中可知创建以下场

Adding to thermoFluid // 读取自constant/fluid/physicalProperties
Adding to rhoFluid
Adding to UFluid
Adding to phiFluid // rhoFluid*UFluid
Adding to gFluid // 读取自constant/fluid/g
Adding to hRefFluid // 参考高度
Adding to pRefFluid // 未从constant中读取到
Adding to ghFluid // gFluid*cellCenter + gFluid*hRefFluid
Adding to ghfFluid // gFluid*faceCenter + gFluid*hRefFluid
Adding to turbulenceFluid // 读取自constant/fluid/momentumTransport
Adding to thermophysicalTransport // 同上
Adding to reactionFluid // 未读取到combustionProperties文件
Adding to KFluid // 0.5*(UFluid)^2
Adding to dpdtFluid
Adding to fieldsFluid

从createSolidFields.H创建了thermoSolid。对于heater，读取到了constant/heater/fvModels，定义了能量源项以实现加热的功能

energySource
{
   
    type            heatSource;
    selectionMode   all; // 选中heater cellZone
    q               1e7;
}

2.2. 固体求解

solveSolid.H求解

{
   
    while (pimple.correctNonOrthogonal())
    {
   
        fvScalarMatrix eEqn
        (
            fvm::ddt(rho, e)
          + thermo.divq(e) // 热传导，拉普拉斯项
          ==
            fvModels.source(rho, e) // 源项
        );

        eEqn.relax();

        fvConstraints.constrain(eEqn);

        eEqn.solve();

        fvConstraints.constrain(e);
    }
}

thermo.correct();

Info<< "Min/max T:" << min(thermo.T()).value() << ' '
    << max(thermo.T()).value() << endl;

固体能量方程本质上是热传导方程：
$\frac{\partial (\rho h)}{\partial t}=\frac{\partial}{\partial x_j}\left( \alpha \frac{\partial h}{\partial x_j} \right)$

$(\rho h)V=\int \alpha \frac{\partial h}{\partial x_j} dS \approx \sum(\alpha \nabla h)\cdot S$

$h$ is the specific enthalpy, kJ/kg, $\rho$ the density and $\alpha = \kappa / c_p$ is the thermal diffusivity which is defined as the ratio between the thermal conductivity $\kappa$ and the specific heat capacity $c_p$ . Note that $\kappa$ can be also anisotropic

对于不可压缩流体， $h$ 仅是温度和压力的函数， $h = h (T, P) = U + P V$

2.3. 流体求解

求解流程：

由连续性方程更新密度；
解动量方程，解出的速度不满足连续性方程；
解浓度输运方程，求出各物质的浓度以及化学反应；
解能量方程，求出温度；
解压力泊松方程使得速度满足连续性方程；
使用新的压力场，结合状态方程求出流体密度。

2.3.1. solveFluid.H

if (!pimple.flow())
{
   
    if (pimple.models())
    {
   
        fvModels.correct();
    }

    if (pimple.thermophysics())
    {
   
        tmp<fv::convectionScheme<scalar>> mvConvection(nullptr);

        if (Ecorr == 0)
        {
   
            #include "YEqn.H" // 求解组分方程
        }
        #include "EEqn.H"
    }
}
else
{
   
    if (Ecorr == 0)
    {
   
        if (!mesh.schemes().steady() && pimples.firstPimpleIter())
        {
   
            #include "rhoEqn.H"
        }

        if (pimple.models())
        {
   
            fvModels.correct();
        }

        #include "UEqn.H"
    }

    if (pimple.thermophysics())
    {
   
        tmp<fv::convectionScheme<scalar>> mvConvection(nullptr);

        if (Ecorr == 0)
        {
   
            #include "YEqn.H"
        }
        #include "EEqn.H"
    }

    if (Ecorr == nEcorr - 1)
    {
   
        tmp<fvVectorMatrix>& tUEqn = UEqns[i]; // 第 i个fluidRegion的U方程
        fvVectorMatrix& UEqn = tUEqn.ref();

        // --- PISO loop
        while (pimple.correct())
        {
   
            #include "../../buoyantFoam/pEqn.H"
        }

        if (pimples.pimpleTurbCorr(i))
        {
   
            turbulence.correct();
            thermophysicalTransport.correct();
        }

        if (!mesh.schemes().steady() && pimples.finalPimpleIter())
        {
   
            rho = thermo.rho();
        }
    }
}

2.3.2. src/finiteVolume/cfdTools/compressible/rhoEqn.H

连续性方程
$\frac{\partial \rho}{\partial t} + \frac{\partial \rho u_j}{\partial x_j} = \frac{\partial \rho}{\partial t}+\nabla\cdot (\rho U)=0$

{
   
    fvScalarMatrix rhoEqn
    (
        fvm::ddt(rho)
      + fvc::div(phi)
      ==      
        fvModels.source(rho)
    );

    fvConstraints.constrain(rhoEqn);

    rhoEqn.solve();

    fvConstraints.constrain(rho);
}

2.3.3. UEqn.H

可压缩流体动量方程
$\frac{\partial \rho U}{\partial t}+\nabla \cdot (\rho U U ) - \nabla \cdot \tau =-\nabla p_{rgh}- g\cdot h \nabla \rho+S$
上式与代码并不对应，下文的代码能够考虑旋转体系中的科里奥利力，施加了MRF，由于本算例并不涉及，所以使用上文的动量方程。关于turbulence.divDevTau(U)的解析可参考文献

// Solve the Momentum equation

MRF.correctBoundaryVelocity(U);

UEqns[i] =
(
    fvm::ddt(rho, U) + fvm::div(phi, U) // 左端前两项
    + MRF.DDt(rho, U)
    + turbulence.divDevTau(U) // 左端第三项
    ==
    fvModels.source(rho, U) // 右端第三项
);
fvVectorMatrix& UEqn = UEqns[i].ref();

UEqn.relax();

fvConstraints.constrain(UEqn);

if (pimple.momentumPredictor())
{
   
    solve
    (
        UEqn
     ==
        fvc::reconstruct
        (
            (
                - ghf*fvc::snGrad(rho) // 右端第二项
                - fvc::snGrad(p_rgh) // 右端第一项
            )*mesh.magSf()
        )
    );

    fvConstraints.constrain(U);
    K = 0.5*magSqr(U