OpenFOAM 编程 | One-Dimensional Transient Heat Conduction

2022/7/31 1:25:13

编程Tag： src 0.0 OpenFOAM dimensional 0.001 lnInclude Conduction FOAM

本文主要是介绍OpenFOAM 编程 | One-Dimensional Transient Heat Conduction，对大家解决编程问题具有一定的参考价值，需要的程序猿们随着小编来一起学习吧！

0. 写在前面

本文中将对一维瞬态热传导问题进行数值求解，并基于OpenFOAM类库编写求解器。该问题参考自教科书\(^{[1]}\)示例 8.1。

1. 问题描述

一维瞬态热传导问题控制方程如下

\[\rho c \frac{\partial T}{\partial t} = \nabla\cdot \left(k\nabla T\right) \]

其中，\(\rho c = 1.0\times10^{+7}\ \mathrm{J/m^3\cdot K}\)，\(k=10\ \mathrm{W/m\cdot K}\)。

假设等截面直杆长为 \(L=0.02\ \mathrm{m}\)，截面为边长 \(0.001\ \mathrm{m}\) 的正方形，全杆初始温度为 \(200 \mathrm{K}\) ，左侧边界条件为\(\nabla T = 0\)，右侧边界条件为\(T=0\)；杆长方向与 \(x\) 轴平行，此处一维问题不考虑 \(y\) 和 \(z\) 方向。

该问题存在解析解

\[\frac{T(x,t)}{200} = \frac{4}{\pi} \sum_{n=1}^\infty \frac{\left(-1\right)^{n+1}}{2n-1} \exp{\left(-\alpha\lambda_n^2t\right)} \cos \left(\lambda_nx\right) \]

其中，\(\lambda_n = \frac{\left(2n-1\right)\pi}{2L}\)，\(\alpha = \frac{k}{\rho c}\)。

2. 数值解法

对于该物理模型，采用均匀正六面体结构化网格，网格数量为10，相邻网格体心距离为 \(\Delta x = 0.002 \mathrm{m}\)，截面面积为\(S = 1\times 10^{-6} \mathrm{m}^2\)，网格体积为 \(V_P = 2\times10^{-9} \mathrm{m}^3\)，网格示意图如下。

对控制方程进行离散（时间项显式格式，固定时间步\(\Delta t = 0.001 \mathrm{s}\)（满足稳定条件）、界面插值采用中心差分格式），可以得到下面线性方程

\[\frac{\rho c V_P}{\Delta t} \left( T_P - T_P^0 \right) = k\sum_N \frac{ T_N - T_P }{\Delta x} S_{N,f} \]

其矩阵形式如下

\[\begin{bmatrix} 20.005 & -0.005 & 0 & \cdots & 0 & 0 & 0 \\ -0.005 & 20.005 & -0.005 & \cdots & 0 & 0 & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots\\ 0 & 0 & 0 & \cdots & -0.005 & 20.005 & -0.005 \\ 0 & 0 & 0 & \cdots & 0 & -0.005 & 20.015\\ \end{bmatrix} \begin{bmatrix} T_0 \\ T_1 \\ \vdots \\ T_8 \\ T_9 \end{bmatrix} = \begin{bmatrix} 20 T_0^0 \\ 20 T_1^0 \\ \vdots \\ 20 T_8^0 \\20 T_9^0 \end{bmatrix} \]

3. OpenFOAM求解

此处，我们把OpenFOAM作为类库使用，可以很快的完成一个求解器，不会涉及过多的底层工作。

3.1 求解器源码

#include "fvCFD.H"
#include <iostream>

int main(int argc, char* argv[])
{
#include "setRootCaseLists.H" 
#include "createTime.H" // 构造 runTime 对象
#include "createMesh.H" // 构造 mesh 对象
    // 密度 x 热容
    dimensionedScalar rhoC("rhoC", dimensionSet(1, -1, -2, 1, 0, 0, 0), scalar(1.e+7));
    // 热导率
    dimensionedScalar k("k", dimensionSet(1, 1, -3, 1, 0, 0, 0), scalar(10.0));
    // 温度场，需要从0文件夹中读取初始值
    volScalarField T(IOobject("T", "0", mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE), mesh);

    while ( runTime.loop() )
    {
        Info << "当前时间 : " << runTime.timeName() << " s" << endl << endl;
        // 构造线性方程组
        fvScalarMatrix TEqn(fvm::ddt(rhoC, T) == fvm::laplacian(k, T));
        // 求解
        TEqn.solve();
        // 更新边界值
        T.correctBoundaryConditions();

        if ( runTime.timeIndex() == 1 )
        { // 打印方程组；这段代码放在哪里无所谓，此代码没有在时间步内更新 TEqn
            Info << "#### UPPER\n" << TEqn.upper() << endl;
            Info << "#### DIAG \n" << TEqn.D() << endl;
            Info << "#### LOWER\n" << TEqn.lower() << endl;
            Info << "#### SOURCE\n" << TEqn.source() << endl; // Right Hand Side
            getchar(); // 此处暂停，按回车继续运行...
        }

        if ( runTime.writeTime() )
        {
            runTime.write();
        }
        runTime.printExecutionTime(Info);
    }
    return 0;
}

3.2 CMakeLists.txt

cmake_minimum_required (VERSION 3.8)
project(OneDimUnsteadyFlow)

# OpenFOAM 安装路径
set( FOAM_PREFIX "/opt/OpenFOAM-v2112" )
# 包含路径
set( FOAM_SRC ${FOAM_PREFIX}/OpenFOAM-v2112/src )
include_directories(
    ${FOAM_SRC}/atmosphericModels/lnInclude
    ${FOAM_SRC}/combustionModels/lnInclude
    ${FOAM_SRC}/conversion/lnInclude
    ${FOAM_SRC}/dummyThirdParty/lnInclude
    ${FOAM_SRC}/dynamicFaMesh/lnInclude
    ${FOAM_SRC}/dynamicFvMesh/lnInclude
    ${FOAM_SRC}/dynamicMesh/lnInclude
    ${FOAM_SRC}/engine/lnInclude
    ${FOAM_SRC}/faOptions/lnInclude
    ${FOAM_SRC}/fileFormats/lnInclude
    ${FOAM_SRC}/finiteArea/lnInclude
    ${FOAM_SRC}/finiteVolume/lnInclude
    ${FOAM_SRC}/functionObjects/lnInclude
    ${FOAM_SRC}/fvAgglomerationMethods/lnInclude
    ${FOAM_SRC}/fvMotionSolver/lnInclude
    ${FOAM_SRC}/genericPatchFields/lnInclude
    ${FOAM_SRC}/lagrangian/lnInclude
    ${FOAM_SRC}/lumpedPointMotion/lnInclude
    ${FOAM_SRC}/mesh/lnInclude
    ${FOAM_SRC}/meshTools/lnInclude
    ${FOAM_SRC}/ODE/lnInclude
    ${FOAM_SRC}/OpenFOAM/lnInclude
    ${FOAM_SRC}/optimisation/lnInclude
    ${FOAM_SRC}/OSspecific/POSIX/lnInclude
    ${FOAM_SRC}/overset/lnInclude
    ${FOAM_SRC}/parallel/lnInclude
    ${FOAM_SRC}/phaseSystemModels/lnInclude
    ${FOAM_SRC}/Pstream/lnInclude
    ${FOAM_SRC}/randomProcesses/lnInclude
    ${FOAM_SRC}/regionFaModels/lnInclude
    ${FOAM_SRC}/regionModels/lnInclude
    ${FOAM_SRC}/renumber/lnInclude
    ${FOAM_SRC}/rigidBodyDynamics/lnInclude
    ${FOAM_SRC}/rigidBodyMeshMotion/lnInclude
    ${FOAM_SRC}/sampling/lnInclude
    ${FOAM_SRC}/semiPermeableBaffle/lnInclude
    ${FOAM_SRC}/sixDoFRigidBodyMotion/lnInclude
    ${FOAM_SRC}/sixDoFRigidBodyState/lnInclude
    ${FOAM_SRC}/surfMesh/lnInclude
    ${FOAM_SRC}/thermophysicalModels/lnInclude
    ${FOAM_SRC}/topoChangerFvMesh/lnInclude
    ${FOAM_SRC}/transportModels/lnInclude
    ${FOAM_SRC}/TurbulenceModels/lnInclude
    ${FOAM_SRC}/waveModels/lnInclude
    .
)

link_directories(
    ${FOAM_PREFIX}/ThirdParty-v2112/platforms/linux64Gcc/boost_1_74_0/lib64
    ${FOAM_PREFIX}/ThirdParty-v2112/platforms/linux64Gcc/fftw-3.3.10/lib64
    ${FOAM_PREFIX}/ThirdParty-v2112/platforms/linux64Gcc/kahip-3.14/lib64
    ${FOAM_PREFIX}/ThirdParty-v2112/platforms/linux64GccDPInt32/lib
    ${FOAM_PREFIX}/ThirdParty-v2112/platforms/linux64GccDPInt32/lib/sys-openmpi

    ${FOAM_PREFIX}/OpenFOAM-v2112/platforms/linux64GccDPInt32Opt/lib
    ${FOAM_PREFIX}/OpenFOAM-v2112/platforms/linux64GccDPInt32Opt/lib/dummy
    ${FOAM_PREFIX}/OpenFOAM-v2112/platforms/linux64GccDPInt32Opt/lib/sys-openmpi
)
set(EXTRA_LIBS dl m)
set(LIBS
    Pstream
    OpenFOAM
    finiteVolume
    meshTools
    fileFormats
    ${EXTRA_LIBS}
)

set( CMAKE_CXX_STANDARD 11 )
set( CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -Xlinker --no-as-needed -Xlinker --add-needed" )
add_definitions(-Dlinux64 -DWM_ARCH_OPTION=64 -DWM_DP -DWM_LABEL_SIZE=32 -DNoRepository -m64 -fPIC )

# 添加可执行文件
add_executable (${PROJECT_NAME} "main.cpp")

# 链接库
target_link_libraries(${PROJECT_NAME} ${LIBS})

4. 算例设置

4.1 system/blockMeshDict

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      blockMeshDict;
}

scale 1.0; // memter

length  0.02; // 长度 


nx  10; // x 方向 网格数量
ny  1;
nz  1;

vertices
(
    (0.0     0.0    0.0)
    ($length 0.0    0.0)
    ($length 0.001  0.0)
    (0.0     0.001  0.0)
    
    (0.0     0.0    0.001)
    ($length 0.0    0.001)
    ($length 0.001  0.001)
    (0.0     0.001  0.001)
);

edges
(
);

blocks
(
    hex (0 1 2 3 4 5 6 7)  ($nx $ny $nz) simpleGgrading (1 1 1) 
);

boundary
(
    left
    {
        type patch;
        faces
        (
            (0 4 7 3)
        );
    }
    right 
    {
        type patch;
        faces
        (
            (1 2 6 5)
        );
    } 
    other
    {
        type empty;
        faces
        (
            (2 3 7 6)
            (4 5 6 7) 
            (0 1 5 4)
            (1 0 3 2)
        );
    }
);

4.2 system/controDict

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      controlDict;
} 

application     OneDimUnsteadyFlow;
startFrom       startTime;
startTime       0;
stopAt          endTime;
endTime         125;
deltaT          0.001;
writeControl    adjustableRunTime;   
writeInterval   0.1;   // 0.1秒为间隔输出数据
purgeWrite      0;
writeFormat     ascii;
writePrecision  6;
writeCompression off;
timeFormat      general;
timePrecision   6;
runTimeModifiable no;

4.3 system/fvSchemes

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      fvSchemes;
} 

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         Gauss linear; 
}

divSchemes
{
    default         Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear uncorrected;
}

4.4 system/fvSolution

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      fvSolution;
} 

solvers
{
    T
    {
        solver          GAMG;
        smoother        GaussSeidel;
        tolerance       1e-08;
        relTol          0.0;
    } 
}

4.5 0/T

FoamFile
{
    version     2.0;
    format      ascii;
    arch        "LSB;label=32;scalar=64";
    class       volScalarField;
    location    "0";
    object      T;
} 

dimensions      [0 0 0 1 0 0 0];

internalField   uniform 200;

boundaryField
{
    left
    {
        type            zeroGradient;
    }
    right
    {
        type            fixedValue;
        value           uniform 0;
    }
    other
    {
        type            empty;
    }
}

5. 求解计算

.
├── build                                // build 目录，用于编译代码
├── CMakeLists.txt                       // 项目管理，内容见3.2节
├── main.cpp                             // 求解器源代码，内容见3.1节
└── OneDimUnsteadyFlowCase               // 算例所在目录 ***
    ├── 0                                    // 0 文件夹，保存初始条件
    │   └── T                                    // 本示例中只有温度场 T，故此处只有 T 文件，内容见 4.5 节
    ├── OneDimUnsteadyFlowCase.foam          // 算例目录名称+foam扩展名，空文件，仅作ParaView加载结果使用
    └── system                               // system 目录
        ├── blockMeshDict                        // blockMesh字典文件，内容见 4.1 节
        ├── controlDict                          // 求解器运行控制字典文件，内容见 4.2 节
        ├── fvSchemes                            // 有限体积数值格式字典文件，内容见 4.3 节
        └── fvSolution                           // 求解器参数设置字典文件，内容见 4.4 节

5.1 编译求解器

主要命令解释：

$ cd build/  # 从当前目录切换到路径 ./build
$ cmake ..   # 执行 CMake 生成构建文件（当前生成的是MakefIle）
$ make       # 执行 Make，编译代码

5.2 运行求解器