OpenFOAM 流体模拟入门

// blockMeshDict (简化)
vertices
(
    (0 0 0)   // 0-left bottom back
    (2 0 0)   // 1-right bottom back
    (2 1 0)   // 2-right top back
    (0 1 0)   // 3-left top back
    ...
);

blocks
(
    hex (0 1 2 3 4 5 6 7) (100 50 1) simpleGrading (1 1 1)
);

boundary
(
    inlet  { type patch; faces (...); }
    outlet { type patch; faces (...); }
    wall   { type wall;  faces (...); }
);

foamJob snappyHexMesh 2>&1 | tee log.snappy

// 0/U (稳态单位 m/s)
dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (10 0 0);

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform (10 0 0);
    }
    outlet
    {
        type            zeroGradient;
    }
    wall
    {
        type            noSlip;
    }
    symmetry
    {
        type            symmetry;
    }
}

# 1. 配置分区方案（scotch / simple / hierarchical）
# system/decomposeParDict
decompositionMethod scotch;
numberOfSubdomains  64;

# 2. 执行分区
decomposePar -force

# 3. 并行运行（推荐用 mpirun 或 foamJob）
mpirun -np 64 simpleFoam -parallel > log.solve 2>&1

# 4. 结果重构回单域
reconstructPar

# 生成 ParaView 可读的 .foam 文件
touch case.foam   # 空文件即可
paraFoam          # 或直接 paraview case.foam

# 批量切片采样
postProcess -func sampleDict

# 计算壁面 Y+ 值
yPlusRAS -latestTime

application     simpleFoam;
startFrom       startTime;
startTime       0;
stopAt          endTime;
endTime         2000;
deltaT          1;
writeControl    timeStep;
writeInterval   100;
functions       // 内嵌后处理函数
{
    forces      { type forces; patches (airfoil); ... }
}

# 计算壁面 Y+ 值（后处理）
yPlusRAS -latestTime

# 查看壁面附近流速分布
sample -latestTime

// myScalarTransportFoam.C
#include "fvCFD.H"

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"
    #include "createMesh.H"
    #include "createFields.H"

    Info << "\nStarting time loop\n" << endl;

    while (runTime.loop())
    {
        Info << "Time = " << runTime.timeName() << nl << endl;

        solve
        (
            fvm::ddt(T) + fvm::div(phi, T) - fvm::laplacian(DT, T)
            ==
            fvOptions(T)
        );

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

    Info << "End\n" << endl;
    return 0;
}

# 创建工作目录结构
mkdir -p $FOAM_RUN/mySolver
 foamNewSource application myScalarTransportFoam

# 编辑 Make/files 和 Make/options
# Make/files:
#     myScalarTransportFoam.C
#     EXE = $(FOAM_USER_APPBIN)/myScalarTransportFoam
# Make/options:
#     EXE_INC = -I$(LIB_SRC)/finiteVolume/lnInclude
#     EXE_LIBS = -lfiniteVolume

wmake  # 编译

from PyFoam.RunDictionary.SolutionDirectory import SolutionDirectory
from PyFoam.RunDictionary.ParsedParameterFile import ParsedParameterFile
import itertools

# 参数化研究：不同雷诺数下的翼型升阻力
case = SolutionDirectory("airfoil_case")
reynolds_numbers = [1e5, 5e5, 1e6, 2e6]
angles = [0, 5, 10, 15]

for re, alpha in itertools.product(reynolds_numbers, angles):
    # 克隆 case
    new_case = case.cloneCase(f"re_{re:.0e}_a{alpha}")

    # 修改控制参数
    control = ParsedParameterFile(new_case.controlDict())
    control["endTime"] = 5000
    control.writeFile()

    # 修改物性（运动粘度随雷诺数变化）
    transport = ParsedParameterFile(new_case.transportProperties())
    transport["nu"] = [0, 2, -1, 0, 0, 0, 0, 1.0 / re]
    transport.writeFile()

    # 提交运行
    import subprocess
    subprocess.run(["blockMesh"], cwd=new_case.name)
    subprocess.run(["simpleFoam"], cwd=new_case.name)

from paraview.simple import *

# 批量生成云图
reader = OpenFOAMReader(FileName='case.foam')
rep = Show(reader)
ColorBy(rep, ('POINTS', 'U', 'Magnitude'))
Render()

# 保存截图
SaveScreenshot('velocity_field.png', ImageResolution=[1920, 1080])

# OpenMPI 绑定策略
mpirun --bind-to core --map-by core -np 64 simpleFoam -parallel

# 跨节点通信优化（InfiniBand）
mpirun --mca btl openib,self -np 128 simpleFoam -parallel

# 查看各子域网格数与计算耗时
decomposePar -cellDist

# 后处理：检查并行效率
checkMesh -parallel