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 运行求解器
主要命令解释：
$ cd OneDimUnsteadyFlowCase/ # 从当前目录切换到算例目录
$ blockMesh > log.blockMesh # 运行 blockMesh 画网格，并将标准输出重定向到 log.blockMesh
$ ../build/OneDimUnsteadyFlow # 运行求解器，注意求解器的相对路径
另外，我们也可以看到求解器打印的线性方程组与数值解法中所描述的是一致的。
6. 后处理
我们对比第 /(40 / /mathrm{s}/) 时数值结果与解析结果
云图：
曲线图：
参考文献
[1] H. Versteeg , W. Malalasekera. Introduction to Computational Fluid Dynamics, An: The Finite Volume Method 2nd Edition[M]. Pearson. 2007
原创文章，作者：306829225，如若转载，请注明出处：https://blog.ytso.com/tech/pnotes/277958.html