长江科学院院报

长江科学院院报 >

2026 , Vol. 43 >Issue 2: 181 - 191

DOI: https://doi.org/10.11988/ckyyb.20241284

水工结构与材料

非线性稳态热传导问题的数值流形法

  • 张丽美 , 1, 2 ,
  • 殷跃平 1 ,
  • 郑宏 3 ,
  • 朱赛楠 1 ,
  • 魏云杰 1 ,
  • 张楠 1 ,
  • 杨龙 1, 2
展开
  • 1 中国地质环境监测院, 北京 100081
  • 2 中国地质大学 工程技术学院, 北京 100083
  • 3 北京工业大学 城市建设学部, 北京 100124

张丽美(1992-),女,河南安阳人,工程师,博士研究生,主要从事数值流形法研究。E-mail:

Copy editor: 刘运飞

收稿日期: 2024-12-19

  修回日期: 2025-03-10

  录用日期: 2025-03-10

  网络出版日期: 2025-07-11

基金资助

四川省科技厅项目(N5100012024000900)

云南省重点研发计划项目(202403AA08001)

收起

Numerical Manifold Method for Nonlinear Steady-state Heat Conduction Problems

  • ZHANG Li-mei , 1, 2 ,
  • YIN Yue-ping 1 ,
  • ZHENG Hong 3 ,
  • ZHU Sai-nan 1 ,
  • WEI Yun-jie 1 ,
  • ZHANG Nan 1 ,
  • YANG Long 1, 2
Expand
  • 1 China Institute of Geological Environment Monitoring, Beijing 100081, China
  • 2 College of Engineering and Technology, China University of Geosciences (Beijing), Beijing 100083, China
  • 3 Department of Urban Construction, Beijing University of Technology, Beijing 100124, China

Received date: 2024-12-19

  Revised date: 2025-03-10

  Accepted date: 2025-03-10

  Online published: 2025-07-11

Fold

摘要

利用数值流形法具有统一解决连续和非连续问题的优势,进行了非线性稳态热传导问题分析。首次构建了基于NMM的非线性稳态热传导离散格式,并采用Newton-Raphson迭代算法对非线性方程组进行求解,从而实现了问题域内温度场和热流量场的精确计算。与有限元法相比,NMM在处理复杂边界条件和材料界面问题时展现出更强的适应性,其前处理过程也更为灵活便捷。为验证NMM的有效性,针对二维复杂边界和双层材料等问题域的热传导算例进行模拟,结果表明该方法具有较高的计算精度和良好的鲁棒性。

关键词： 非线性稳态热传导; 数值流形法; Newton-Raphson迭代算法; 二维; 温度场

本文引用格式

张丽美 , 殷跃平 , 郑宏 , 朱赛楠 , 魏云杰 , 张楠 , 杨龙 . 非线性稳态热传导问题的数值流形法[J]. 长江科学院院报, 2026 , 43(2) : 181 -191 . DOI: 10.11988/ckyyb.20241284

Abstract

[Objective] This study addresses the numerical modeling of nonlinear steady-state heat conduction processes where the thermal conductivity varies with temperature. The governing equation for such problems is a second-order quasi-linear partial differential equation, whose nonlinear nature makes analytical solutions extremely challenging, thus necessitating efficient numerical approaches. This work employs the Numerical Manifold Method (NMM) based on quadrilateral mesh covers to analyze and solve two-dimensional nonlinear steady-state heat conduction problems. [Methods] Within the NMM over traditional methods framework, a discrete formulation suitable for nonlinear steady-state heat conduction was established by incorporating three typical boundary conditions: Dirichlet, Neumann, and Robin. The classical Newton-Raphson iterative algorithm was adopted to solve the resulting nonlinear system of equations. A complete numerical solution procedure was implemented on the MATLAB platform to ensure algorithm stability and computational efficiency. To systematically verify the accuracy and robustness of the proposed NMM in handling nonlinear heat conduction, a series of representative numerical examples were designed and conducted. These examples covered various scenarios, including continuous homogeneous materials, discontinuous media containing circular holes, and heterogeneous materials. The simulation results were compared against analytical solutions, existing literature data, or Finite Element Method (FEM) solutions. [Results] 1) Compared to the traditional Finite Element Method (FEM), NMM demonstrates significant theoretical and practical advantages when simulating problem domains with complex geometries or internal discontinuities. This advantage primarily stems from its distinctive numerical characteristics: in NMM, the interpolation subdomains are independent of the subdomains used for numerical integration, whereas in FEM they coincide entirely on the same mesh. Furthermore, FEM is prone to mesh distortion when handling complex boundaries, which can degrade accuracy and impair computational efficiency. Leveraging its physical cover system, NMM can accurately describe complex geometric boundaries. At material interfaces, the different heat conduction behaviors across materials are naturally captured through physical covers and local functions without introducing additional interface conditions, thereby simplifying the computational process and enhancing efficiency. 2) The proposed NMM not only achieves high accuracy in temperature field and heat flux distribution across all examples but also exhibits excellent stability and convergence when dealing with discontinuous interfaces and complex geometries, fully validating the method’s effectiveness and reliability for nonlinear steady-state heat conduction problems. [Conclusion] This study successfully applies NMM to solve two-dimensional nonlinear steady-state heat conduction problems. Through comprehensive comparative analysis and numerical validation, the unique advantages of this method in handling complex engineering thermal problems are highlighted. It provides a novel solution for the numerical simulation of nonlinear heat conduction problems and extends the application scope of NMM in the field of computational thermal physics.

Key words： nonlinear steady-state heat conduction; numerical manifold method; Newton-Raphson iterative algorithm; two-dimensional; temperature field

0 引言

热传导分析在土木工程[1]、水利工程[2]和地热开采[3]等诸多领域都具有重要意义。以岩土工程为例,温度场的改变会显著影响黏土材料的物理力学性质[4]。考虑热传导系数随温度变化时,控制方程呈现非线性特征,此类问题的解析解通常难以获得,因此发展相应的数值方法尤为重要。目前,求解热传导问题的数值方法主要包括有限差分法[5]、有限元法[6]、边界元法[7]、无网格方法[8-9]和比例边界有限元[10]以及数值流形法[1,11]((Numerical Manifold Method, NMM)等。
近几十年来,国内外研究人员致力于解决非线性热传导问题。Annasabi等[6]使用B样条函数结合有限元法分析了非线性热传导问题。李腊梅等[12]将计算模型离散为块体和接触边界,提出一种非连续介质中热传导过程的数值模拟方法。刘承论等[13]利用边界元法计算了三维非稳态热传导问题。梁钰等[14]提出一种快速算法,并应用到导热系数随温度变化的瞬态热传导问题中。Sivin等[15]根据几何结构和问题的复杂程度,采用非标准网格技术对区域进行离散化,结合不动点迭代的比例边界有限元法求解了非线性瞬态热传导问题。Mierzwiczak等[16]基于Kirchhoff变换,将非线性热传导问题转化为线性,利用奇异边界元法成功求解了不规则区域中的稳态温度场。传统边界元法能够高效求解线性热传导问题。然而,对于非线性热传导问题,由于控制方程中通常包含与温度相关的非线性项,这些非线性项无法直接通过边界积分方程表示,导致传统边界元法难以直接应用。因此,必须采用特殊的方法来处理非线性问题,例如引入径向积分法、域分解技术或与其他数值方法耦合。Yang等[17]开发了新的径向积分边界元法(Radial Integration Boundary Element Method,RIBEM),并应用Newton-Raphson方法求解离散代数方程,计算了二维和三维非线性热传导问题。吴泽艳等[9]提出一种混合的无网格方法求解非线性瞬态热传导问题。Khosravifard等[18]采用了改进的无网格径向点插值方法,分析了涉及非均匀热源功能梯度材料的非线性瞬态传热问题。无网格法在计算非线性热传导中,其计算精度和效率与节点分布密切相关,节点布置不当容易出现计算效率低和不收敛的情况。
1991年,Shi[19]提出了NMM,包括数学覆盖和物理覆盖2套系统。通常以有限元的三角形或者四边形网格作为数学覆盖,共用一个节点的所有单元组成一个数学片,数学片被几何边界、材料界面或裂纹切割形成物理片,所有的物理片联结形成物理覆盖,物理覆盖需要精确的匹配问题域。物理片之间相互切割形成流形单元,流形单元是最小的积分单元。目前,NMM已经成功应用于热传导[20]、渗流[21]、裂纹扩展[22]、无界域问题[23]和边坡稳定性[24]等多个领域中。其中,在热传导计算中,胡国栋等[25]利用NMM研究了二维稳态功能梯度材料中的热传导问题。Zhang等[11]发展了无网格数值流形法求解功能梯度材料瞬态非线性热传导问题。Tan等[26]基于NMM模拟求解了稳态和瞬态热传导问题。Wang等[27]对NMM的局部网格进行加密,同时引入热力耦合模型研究岩石热破裂问题。Li等[28]创新性地将NMM应用于复杂工况下的压电材料分析,成功解决了机电荷载作用下含有任意形状孔洞缺陷的压电板的力学响应问题。Zhang[29]构建了三维土石混合体的数学模型,并将三维NMM应用于路基渗透性分析,为复杂地质条件下的路基渗透特性研究提供了新的方法。传统的有限元法在计算复杂边界问题时,网格容易产生畸变,导致计算精度下降并影响计算效率。NMM通过其独特的物理覆盖系统,能够精确描述复杂几何边界。此外,在材料界面处,NMM通过物理覆盖和局部函数,能够自然地描述不同材料间的热传导行为,无需引入额外的界面条件,从而简化了计算过程并提高了计算效率。本文首次构建了基于NMM的非线性稳态热传导离散格式,采用Newton-Raphson迭代算法对非线性方程组进行求解,并基于MATLAB软件实现数值求解过程。通过5个典型的算例验证了本文方法的正确性和鲁棒性,为非线性热传导问题的求解提供了一种新方法。

1 控制方程

如图1所示,二维问题域Ω由Dirichlet、Neumann和Robin三类边界条件所围成。
图1 非线性热传导的问题域

Fig.1 Problem domain of nonlinear heat conduction

非线性稳态热传导的控制方程可以描述为二阶拟线性偏微分方程[20],即
$\frac{\partial }{\partial {x}_{i}}\left(k\left(T\right(x\left)\right)\frac{\partial T\left(x\right)}{\partial {x}_{i}}\right)+Q\left(x\right)=0,　x\in \Omega 　。$
式中:重复索引i在二维问题域中表示1, 2。T(x)表示问题域的温度;k(T)为随温度变化的导热系数;Q(x)是热源函数。
为了确定温度场,附加Dirichlet边界条件,即
$T\left(x\right)=\overline{T}\left(x\right)\mathrm{ },　x\in {\Gamma }_{\mathrm{D}}　。$
Neumann边界条件为
$-\mathrm{k}\left(\mathrm{T}\right(\mathrm{x}\left)\right)\frac{\partial T\left(x\right)}{\partial {x}_{i}}\mathrm{n}=\stackrel{-}{q}\left(x\right)\mathrm{ },　x\in {\Gamma }_{\mathrm{N}}　。$
Robin边界条件为
$-k\left(T\right(x\left)\right)\frac{\partial T(x,t)}{\partial {x}_{i}}n=h(T-{T}_{\mathrm{\infty }})\mathrm{ },　x\in {\Gamma }_{\mathrm{R}}　。$
式中:ГD、 ГN、ГR分别为Dirichlet、Neumann和Robin边界;$\overline{T}$和$\stackrel{-}{q}$分别为ГD和ГN上给定的温度和热流量;h为换热系数;T∞是环境温度;n=(n1,n2)是单位外法向向量。

2 NMM求解非线性热传导

2.1 NMM简介

NMM引入2种类型的覆盖系统,分别为用于建立单位分解函数的数学覆盖和构造局部逼近函数的物理覆盖,它们可以统一地解决连续和不连续问题。本节只对NMM进行简单的介绍,具体细节可以参考文献[30]—文献[32]。
一般来讲,组成数学覆盖的数学片可以是任何形状,可以是有限元网格中共用1个节点的所有单元所组成,也可以是移动最小二乘法中节点的影响域。如图2所示,ABCD围成问题域Ω,EF为材料界面。以红色圆点所连接的四个有限单元构成1个数学片,用MPi表示,i=1,2,…,n;其中n表示数学片的个数。所有的数学片联结形成的数学覆盖须完全覆盖住问题域Ω。
图2 NMM中的数学片、物理片和流形单元

Fig.2 Mathematical patches, physical patches, and manifold elements in NMM

问题域边界和材料界面切割数学片形成物理片,1个数学片可以生成1个或者多个物理片。PPi-j是从数学片MPi中生成的第j个物理片,j=1,2,…,mi;其中,mi是指同一个数学片中生成物理片的个数。图2显示了3种典型的物理片:第1类物理片如PP1、PP5、PP6、PP7和PP8,它们由数学片MP1、MP5、MP6、MP7和MP8生成,这类物理片没有被切割,因此,第1类物理片就是数学片本身;第2类物理片如PP2、PP3,它是由数学片被边界切割并丢弃域外部分所形成;第3类物理片如PP4-1和PP4-2,是数学片MP4被材料界面EF切割生成的。一共生成np个物理片,所有的物理片组成了物理覆盖,精确覆盖问题域。物理片PP5、PP6、PP7和PP8的相互切割,最后形成了流形单元E1。当有限元网格用作数学覆盖时,流形单元是基本的积分单元。
对于每1个数学片,都有一个充分光滑的权函数wi(x)满足:
$\left\{\begin{array}{l}{w}_{i}\left(\mathrm{x}\right)\le 1\mathrm{ },　\mathrm{x}\in M{\mathrm{P}}_{i}　;\\ {w}_{i}\left(\mathrm{x}\right)=0\mathrm{ },　\mathrm{x}\notin M{\mathrm{P}}_{i}　;\\ \stackrel{n}{\sum _{i=1}}{w}_{i}\left(x\right)=1\mathrm{ },　x\in \Omega 　。\end{array}\right.$
每一个物理片都有一个对应的权函数$\left\{{w}_{i-j}\right\}$,它取自数学片上的权函数在物理片上的限制,其表达式为
${w}_{i-j}\left(x\right)=\left\{\begin{array}{l}{w}_{i}\left(x\right)\mathrm{ },　x\in \mathrm{P}{\mathrm{P}}_{i-j}　;\\ 0\mathrm{ },　x\notin \mathrm{P}{\mathrm{P}}_{i-j}　。\end{array}\right.$
至此,构造出NMM的逼近u(x)为
$\mathrm{u}\left(\mathrm{x}\right)=\stackrel{{n}_{\mathrm{p}}}{\sum _{i,j=1}}{w}_{i-j}\left(x\right){u}_{i-j}\left(x\right)　。$
式中:ui-j(x)代表物理片PPi-j的局部逼近,本文中ui-j(x)取未知常数。为简化表示,所有物理片都按单个指标排序,方程简化为
$u\left(x\right)=N\left(x\right)a　。$
式中:N(x)是具有权函数的行向量,wi(x)为其中的元素;a=(ai),并且ai=u(xi)。

2.2 弱形式

Dirichlet边界条件是按照罚函数方法施加,带罚因子的弱形式为
$\begin{array}{l}{\int }_{\Omega }\left[-\frac{\partial }{\partial {x}_{i}}\left(k\left(T\right)\frac{\partial T}{\partial {x}_{i}}\right)-Q\right]\delta Td\mathrm{\Omega }+{\int }_{{\Gamma }_{\mathrm{D}}}\mathrm{\lambda }(\mathrm{T}-\stackrel{-}{T})\mathrm{\delta }T\mathrm{d}\Gamma +\\ {\int }_{{\Gamma }_{\mathrm{N}}}\stackrel{-}{q}\mathrm{\delta }T\mathrm{d}\Gamma +{\int }_{{\Gamma }_{\mathrm{R}}}h(T-{T}_{\mathrm{\infty }})\mathrm{\delta }T\mathrm{d}\Gamma =0　。\end{array}$
式中:δT(x)是T(x)的变分; λ是罚因子,本文采用罚函数方法施加Dirichlet边界条件是因为其简单方便。另外,拉格朗日乘子法[33]和本质边界条件精确施加[34]也可以使用。
温度T(x)采用分离变量法离散,
$T\left(x\right)=N\left(x\right)T　。$
式中:T是所有物理片上未知温度的列向量。变分δT(x)可以近似为
$\mathrm{\delta }T\left(x\right)=N\left(x\right)\delta T　。$
通过将式(10)代入式(11),可以得到
$K\left(k\left(T\right)\right)T=F　。$
式中:K(k(T))和F分别表示热传导矩阵和热荷载列阵,可以通过单元组装形成。在任意一个单元内Ke和Fe的形式分别为:
${K}^{e}={\int }_{{\Omega }^{\mathrm{e}}}k\left(T\right){B}^{T}B\mathrm{d}\Omega +{\int }_{{\mathrm{\Gamma }}_{{}_{\mathrm{D}}}^{\mathrm{e}}}\lambda {N}^{T}N\mathrm{d}\Omega +{\int }_{{\Gamma }_{R}^{e}}h{N}^{T}N\mathrm{d}\Omega 　,$
${F}^{e}={\int }_{{\mathrm{\Omega }}^{\mathrm{e}}}{N}^{T}\mathrm{Q}\mathrm{d}\mathrm{\Omega }+{\int }_{{\Gamma }_{D}^{e}}\mathrm{\lambda }{N}^{\mathrm{T}}\stackrel{-}{T}\mathrm{d}\Omega -{\int }_{{\Gamma }_{N}^{e}}{N}^{\mathrm{T}}\stackrel{-}{q}\mathrm{d}\Omega +{\int }_{{\Gamma }_{R}^{e}}h{T}_{\mathrm{\infty }}{N}^{T}\mathrm{d}\Omega 　。$
方程中的矩阵形式可以写为
$B=[{B}_{1},{B}_{2},\dots,{B}_{{\mathrm{n}}_{\mathrm{p}}}]　。$
其中,Bi=${\left(\begin{array}{ll}\frac{\partial {w}_{i}}{\partial {x}_{1}}& \frac{\partial {w}_{i}}{\partial {x}_{2}}\end{array}\right)}^{\mathrm{T}}$。

2.3 Newton-Raphson法

本文采用New-Raphson迭代方法,用来求解非线性方程组,残差可以定义为
R=F-KT　。
使用泰勒公式,可以得到
${R}^{n+1}={R}^{n}+\frac{\partial R}{\partial T}\mathrm{\Delta }T+\frac{{\partial }^{2}R}{\partial {T}^{2}}{\left(\mathrm{\Delta }T\right)}^{2}+\dots 　。$
式中:Rn和Rn+1分别为第n和n+1次的残差。
取式(17)的线性项,可以得到下面的方程
${R}^{n+1}={R}^{n}+\frac{\partial R}{\partial T}\mathrm{\Delta }T　。$
式(16)对T的一阶偏导数为
$\frac{\partial R}{\partial T}=-\frac{\partial K}{\partial T}T-K　。$
式(19)右端第1项相对较小,可以忽略,只保留第2项,式(19)可以等价为
$\frac{\partial R}{\partial T}=-K　。$
将式(20)代入式(18),并假定Rn+1为0,可以得出以下公式
$　　\mathrm{\Delta }T={K}^{-1}{R}^{n}　。$
然后,向量T可以更新为
${T}^{\mathrm{n}+1}={T}^{n}+\theta \mathrm{\Delta }T　。$
式中θ被称为松弛因子,它的取值范围是0~1,旨在提高收敛性。在本文计算过程中,θ=1。

3 数值算例

本文提出了5个算例,用来验证NMM在非线性热传导问题中的计算精度。其中,第1个和第4个算例分别为规则四边形和矩形区域的热传导问题,均可以简化为有解析解的一维热传导模型。其余3个算例分别为矩形板中间含圆孔、二维不规则板、由2种材料组成的机械构件的热传导。为了评估数值结果,定义相对误差Re为
${R}_{\mathrm{e}}=\frac{\left|T-{T}^{\mathrm{N}\mathrm{M}\mathrm{M}}\right|}{T}\times 100\mathrm{\%}　。$
式中:T代表解析解或参考解;TNMM代表由NMM计算的数值解。

3.1 方板上的热传导

如图3所示的第1个算例,方板1.0 m×1.0 m上随温度变化的导热系数为k(T)=15+0.01T2,左、右2条边界上施加的温度分别为Tl=100 ℃, Tr=200 ℃。上、下2条边界为绝热边界。
图3 方板及边界条件

Fig.3 Square plate and boundary conditions

本算例采用的罚因子λ=1×1010。此算例可简化为一维热传导问题,推导出该算例的解析解[35]为
$\begin{array}{l}T\left({x}_{1}\right)=\frac{50\times {5}^{\frac{1}{3}}\times {6}^{\frac{2}{3}}}{{\left(-a{x}_{1}-b+\sqrt{150\mathrm{ }000+(a{x}_{1}{)}^{2}+2a{x}_{1}b+{b}^{2}}\right)}^{\frac{1}{3}}}-\\ {5}^{\frac{2}{3}}\times {6}^{\frac{1}{3}}{\left(-a{x}_{1}-b+\sqrt{150\mathrm{ }000+(a{x}_{1}{)}^{2}+2a{x}_{1}b+{b}^{2}}\right)}^{\frac{1}{3}}。\end{array}$
式中:a=24 833.3367; b=4 833.33。
首先,为了验证基于Newton-Raphson法的NMM的收敛性。如图4,采用了4种密度的物理片覆盖问题域,物理片的个数分别为121、441、961、1 681。表1列出了采用4种不同物理片,计算沿着x1方向的温度值。由结果可以看出,随着物理片数目增加,NMM计算的温度值越来越接近解析解,即使采用121个物理片覆盖问题域,最大的误差仅为0.106 5%,可以得到相对精确的解。
图4 问题域的离散

Fig.4 Discretization of problem domain

表1 物理片沿着x1方向的温度

Table 1 Temperature distribution along x1 direction with different numbers of physical patches

x1/m 不同数量物理片对应的温度/℃ 解析解
温度/℃
121 441 961 1 681
0.1 118.325 5 118.388 7 118.412 4 118.424 8 118.451 6
0.2 132.746 6 132.776 4 132.787 3 132.792 9 132.806 3
0.3 144.733 1 144.746 9 144.751 8 144.754 4 144.761 0
0.4 155.102 1 155.107 2 155.109 0 155.109 9 155.112 6
0.5 164.305 3 164.305 3 164.305 3 164.305 3 164.305 3
0.6 172.607 6 172.611 5 172.612 8 172.613 4 172.615 1
0.7 180.214 1 180.220 4 180.222 4 180.223 4 180.225 9
0.8 187.253 1 187.260 9 187.263 4 187.264 6 187.267 4
0.9 193.819 0 193.827 7 193.830 6 193.831 9 193.834 8
采用441个物理片覆盖问题域,针对罚因子λ在104~1010范围内的取值进行了测试。图5展示了x1=0.0~1.0 m,x2=0.5 m时,解析解与不同罚因子下NMM计算结果的相对误差分布。结果表明,当罚因子λ=104时,NMM计算结果和解析解存在显著偏差;当λ>105时,误差逐渐收敛并趋于稳定。基于此,本算例最终取值λ=1010作为罚因子。文献[36]同样对NMM渗流场计算中罚因子的取值进行了讨论。
图5 x2=0.5 m时不同罚因子下的误差对比

Fig.5 Error comparison with different penalty factors when x2=0.5 m

同时,图6绘制了解析解和采用4种物理片的NMM计算的温度值。图7和图8分别展示了NMM计算结果和解析解在x1方向的温度相对误差值和热流量分布。很明显随着物理片密度的增大,NMM计算结果与解析解之间的误差越来越小,进一步表明所提NMM有很好的收敛性。
图6 不同物理片计算的温度值

Fig.6 Computed temperatures with different physical patches

图7 解析解和不同物理片NMM计算的温度相对误差

Fig.7 Relative errors of temperature between NMM results and analytic solution with different physical patches

图8 x2=0.5 m时沿着x1=0~1 m不同物理片的热流量

Fig.8 Heat flux along x1=0~1 m with different physical patches at x2=0.5 m

3.2 含圆孔矩形板的热传导

如图9所示,研究一个含有圆孔的矩形板,矩形板的尺寸为4.0 m×2.0 m,分析圆孔半径分别为r=0.4 m和r=0.75 m的2种情况。随温度变化的导热系数为k(T)=100+0.2T。左、右2条边界施加的温度分别为100 ℃和400 ℃,其他的边界都是绝热边界。本算例罚因子λ=1×1012。
图9 带圆孔的矩形板及边界条件

Fig.9 Rectangular plate with a circular hole and boundary conditions

图10显示了2种半径情况下,覆盖问题域的物理片数量分别为840个和744个。本研究计算了顶部边界DC,x1=-2.0~2.0 m, x2=1.0 m的温度分布,并与文献[35]中采用96个边界节点和144个内部节点的RIBEM计算结果进行了对比。图11的结果表明,NMM和RIBEM两种方法计算的温度分布具有良好的一致性。此外,图12进一步展示了整个计算域内半径r=0.4 m和r=0.75 m两种情况的温度等值线分布,其与RIBEM方法绘制的等值线表现出高度吻合的特征,这充分验证了NMM的可靠性和准确性。
图10 不同半径的问题域离散情况

Fig.10 Discretization of problem domain with different radii

图11 沿着上边界CD的温度分布

Fig.11 Temperature distribution along upper boundary CD

图12 不同半径的温度分布等值线

Fig.12 Temperature contours with different radii

3.3 二维不规则板的热传导

计算1个内部含有圆孔洞且半径r=1 cm的不规则板的算例,模型尺寸和边界条件如图13(a)、图13(b)所示。热传导系数随温度变化,其函数表达式为k(T)=50+T。BC和CD边界的温度值是400 ℃,AE边界的温度是100 ℃,其余边界皆为绝热边界。本算例罚因子λ=1×107。
图13 几何尺寸和边界条件

Fig.13 Geometric dimensions and boundary conditions

如图14所示,采用NMM计算时,使用540个物理片覆盖问题域。考虑到几何模型和边界条件均具有对称特性,为提高计算效率,研究选取半圆形区域作为分析对象。计算过程中,迭代了8次,最后2次的误差为4.66×10-5。
图14 540个物理片离散问题域

Fig.14 Discretization of problem domain with 540 physical patches

鉴于该算例具有不规则几何特征,解析解难以直接获取。为验证NMM的计算精度,本研究采用COMSOL Multiphysics软件建立有限元对比模型。如图15所示,有限元模型利用5 182个三角形单元进行空间离散,形成包含10 606个自由度的计算网格。表2详细列出了半圆形区域MN边界处温度场的数值对比结果。计算数据显示,NMM与有限元解的最大相对误差仅为0.270 5%,且即使采用540个物理片的简化离散方案,NMM仍能保持较高计算精度。这一对比验证表明,NMM在复杂几何问题中具有优异的计算性能,即使在较少的物理片下亦可获得可靠的工程预测结果。
图15 5 182个三角形单元的有限元网格

Fig.15 Finite element meshes of 5 182 triangular elements

表2 半圆MN上的温度对比

Table 2 Comparison of data along semicircle MN

角度/(°) 2种数值方法模拟的温度/℃ 相对误差/%
FEM NMM
0 206.212 206.189 0.010 8
15 213.797 214.249 0.211 4
30 234.151 233.911 0.102 6
45 261.848 261.402 0.170 1
60 291.467 290.781 0.235 5
75 319.142 318.279 0.270 5
90 342.639 342.345 0.085 8
105 361.005 360.606 0.110 5
120 374.235 373.799 0.116 7
135 382.982 382.663 0.083 4
150 388.234 387.971 0.067 8
165 390.946 390.748 0.050 7
180 391.773 391.636 0.035 1
图16(a)和图16(b)分别展示了FEM和NMM对不规则带孔板温度场的模拟结果。通过对比温度场等值线云图可以发现,2种数值方法呈现出高度一致的温度分布特征。
图16 FEM和NMM计算的温度等值线

Fig.16 Temperature contours obtained by FEM and NMM

3.4 含两层材料矩形板的热传导

图17为分层的2种材料Ω=Ω1∪Ω2,几何尺寸和边界条件[16]已给出。左侧正方形区域的热传导系数为k1(T1)=k01(1+k11T1)=1+0.1T,右侧矩形的热传导系数为k2(T2)=k02(1+k12T2)=2(1+0.25T2),左边边界施加的温度Tl=750 K,右边界的对流换热系数和环境温度分别为h=20 W/(m2K)和T∞=300 K,上下边界为绝热边界。
图17 含2层材料矩形板的几何尺寸和边界条件

Fig.17 Geometric dimensions and boundary conditions of two-layer rectangular plate

温度场可以简化为一维,只沿着x1方向发生变化,该问题的解析解[16]为:
${T}_{1\mathrm{e}\mathrm{x}\mathrm{a}\mathrm{c}\mathrm{t}}=\frac{-{k}_{01}+\sqrt{{k}_{01}^{2}-2{k}_{01}{k}_{11}\left(q\right({x}_{1}+1\left)\right)-{k}_{01}{T}_{\mathrm{l}}-0.5{k}_{01}{k}_{11}{T}_{l}^{2}}}{{k}_{01}{k}_{11}},\mathrm{ }-1\le {x}_{1}\le 0　;$
${T}_{2\mathrm{e}\mathrm{x}\mathrm{a}\mathrm{c}\mathrm{t}}=\frac{-{k}_{02}+\sqrt{{k}_{02}^{2}-2{k}_{02}{k}_{12}(q{x}_{1}-{k}_{02}c-0.5{k}_{02}{k}_{12}{c}^{2})}}{{k}_{02}{k}_{12}},\mathrm{ }0\le {x}_{1}\le 2　。$
式中:c=660.354 644 5 K; q=6 411.232 529 W/m2。
如图18所示,采用616个物理片离散此问题域,罚因子λ=1×1010。计算过程中,迭代了5次,且最后2次的误差为2.237 2×10-5。表3列出了由NMM和解析解计算沿着x方向的温度值以及误差值。可以看出本文方法和解析解间最大的相对误差仅仅为1.144×10-5。图19绘制了数值解和解析解的温度值。图20绘制了在问题域内的温度云图分布。本算例再次验证本文方法的精度和可行性。
图18 616个物理片离散问题域

Fig.18 Problem domain discretized with 616 physical patches

表3 沿着x1方向的温度

Table 3 Temperatures along x1 direction

x1/m 2种数值方法模拟的温度/K 相对误差/10-5
NMM 解析解
-0.9 741.510 5 741.516 8 0.850
-0.8 732.932 1 732.936 8 0.643
-0.7 724.253 2 724.256 5 0.457
-0.6 715.465 1 715.472 4 1.011
-0.5 706.577 0 706.580 5 0.497
-0.4 697.571 5 697.577 0 0.790
-0.3 688.449 6 688.457 4 1.137
-0.2 679.215 7 679.217 1 0.214
-0.1 669.843 7 669.851 3 1.144
0 660.354 6 660.354 6 0.000
0.2 656.482 9 656.483 2 0.057
0.4 652.588 8 652.589 0 0.029
0.6 648.671 3 648.671 5 0.029
0.8 644.730 0 644.730 4 0.061
1.0 640.764 8 640.765 2 0.056
1.2 636.775 4 636.775 5 0.011
1.4 632.760 4 632.760 7 0.046
1.6 628.720 1 628.720 5 0.062
1.8 624.654 1 624.654 3 0.028
2 620.561 4 620.561 6 0.039
图19 沿着x1方向的温度分布

Fig.19 Temperature distribution along x1 direction

图20 温度分布云图

Fig.20 Temperature contours

3.5 机械构件上的对流换热

如图21所示,考虑一个内部含2个不同半径孔洞的机械构件[37],构件的几何尺寸和边界条件在图中已经标明。构件由2种材料组成,其中区域Ω1的导热系数为k1(T)=550-0.1T,区域Ω2和Ω3是同一种材料,热传导系数为k2(3)(T)=80+0.5T。左、右2个内部圆孔的温度分别为100 ℃和500 ℃。整个构件处于热对流环境中,对流换热系数和环境温度分别为:h=100 W/(m2·℃),T∞=50 ℃。
图21 构件的几何尺寸和边界条件

Fig.21 Geometric dimensions and boundary conditions of component

如图22所示,此构件采用四边形网格作为数学覆盖,经过材料边界和几何边界切割后,生成792个物理片。由于此算例几何形状较为复杂,解析解不易求得,因此将采用有限元软件ANSYS计算相应结果作为对比,ANSYS软件采用八节点二次单元对问题域进行离散。有限元模型所用单元数总计476个,节点总数为1 149个。
图22 792个物理片离散问题域

Fig.22 Discretization of problem domain with 792 physical patches

为了验证本方法的正确性,将NMM计算所得上半圆的结果与有限元计算结果进行对比。鉴于该几何形状上下对称,仅取上半圆进行分析。对上半圆选取18个节点进行数据对比,结果如表4所示。有限元和本文所提方法最大相对误差仅仅为0.442%,表明本文方法具有较高的计算精度。图23展示了该模型的温度分布等值线图,从中可以看出NMM与FEM的计算结果高度吻合。
表4 界面圆环上的温度计算结果及对比

Table 4 Calculated temperature results and comparison along interface ring

角度/
(°)		 左边的界面圆环 右边的界面圆环
温度/℃ 相对误
差/%		 温度/℃ 相对误
差/%
ANSYS NMM ANSYS NMM
0 256.32 255.19 0.442 451.69 451.87 0.040
22.5 248.49 247.77 0.288 450.65 450.74 0.019
45.0 227.69 227.25 0.192 447.11 447.32 0.046
67.5 200.23 199.99 0.117 439.59 439.72 0.030
90.0 173.61 173.29 0.186 424.92 424.93 0.002
112.5 153.16 152.86 0.193 400.40 400.65 0.062
135.0 139.73 139.491 0.175 373.54 373.84 0.080
157.5 132.23 131.93 0.230 354.18 354.02 0.046
180.0 129.83 129.27 0.432 347.26 347.50 0.069
图23 温度分布等值线

Fig.23 Temperature contours

4 结束语

本文采用四边形网格覆盖的NMM求解二维非线性稳态热传导问题,通过连续材料、含圆孔的非连续材料以及非均值材料等典型算例,系统验证了NMM方法的可行性和计算精度。与传统FEM相比,NMM在模拟不规则问题域热传导方面具有显著优势。这种优势主要源于NMM独特的数值特性:在NMM中,插值子域与积分子域相互独立,而FEM中2种子域是相同的。此外,NMM能够精确地描述复杂几何边界,并充分利用局部域中解的行为特征,这使得其在裂纹扩展模拟方面展现出比FEM更强的灵活性。基于当前研究成果,后续研究致力于将NMM方法拓展至二维热裂纹问题中。

参考文献

原文顺序 | 文献年度倒序 | 文中引用次数倒序
[1]
林绍忠, 明峥嵘, 祁勇峰. 用数值流形法分析温度场及温度应力[J]. 长江科学院院报, 2007, 24(5): 72-75.

(LIN Shao-zhong, MING Zheng-rong, QI Yong-feng. Thermal Field and Thermal Stress Analysis Based on Numerical Manifold Method[J]. Journal of Yangtze River Scientific Research Institute, 2007, 24(5): 72-75.(in Chinese))

[2]
任继勋, 佳琳, 阳建新, 等. 渗流条件下地下水含盐量对基坑坑底冻结温度场影响的数值模拟[J]. 长江科学院院报, 2024, 41(1): 151-158, 166.

(REN Ji-xun, JIA Lin, YANG Jian-xin, et al. Numerical Simulation on the Effect of Groundwater Salinity on Freezing Temperature Field at the Bottom of Foundation Pit under Seepage Conditions[J]. Journal of Changjiang River Scientific Research Institute, 2024, 41(1): 151-158, 166.(in Chinese))

[3]
郭平业, 卜墨华, 张鹏, 等. 高地温隧道灾变机制与灾害防控研究进展[J]. 岩石力学与工程学报, 2023, 42(7): 1561-1581.

(GUO Ping-ye, BU Mo-hua, ZHANG Peng, et al. Review on Catastrophe Mechanism and Disaster Countermeasure of High Geotemperature Tunnels[J]. Chinese Journal of Rock Mechanics and Engineering, 2023, 42(7): 1561-1581.(in Chinese))

[4]
江文豪, 冯晨, 李江山. 饱和黏土一维非线性固结与热传导耦合模型[J]. 岩石力学与工程学报, 2023, 42(10): 2588-2600.

(JIANG Wen-hao, FENG Chen, LI Jiang-shan. Coupled Model for One-dimensional Nonlinear Consolidation and Heat Conduction in Saturated Clay[J]. Chinese Journal of Rock Mechanics and Engineering, 2023, 42(10): 2588-2600.(in Chinese))

[5]
史策. 热传导方程有限差分法的MATLAB实现[J]. 咸阳师范学院学报, 2009, 24(4): 27-29, 36.

(SHI Ce. Heat Conduction Equation Finite Difference Method to Achieve the MATLAB[J]. Journal of Xianyang Normal University, 2009, 24(4): 27-29, 36.(in Chinese))

[6]
ANNASABI Z, ERCHIQUI F. Robust Kirchhoff Transformation Using B-spline for Finite Element Analysis of the Non-linear Heat Conduction[J]. International Communications in Heat and Mass Transfer, 2021, 120: 104985.

DOI

[7]
FENG W Z, GAO X W. An Interface Integral Equation Method for Solving Transient Heat Conduction in Multi-medium Materials with Variable Thermal Properties[J]. International Journal of Heat and Mass Transfer, 2016, 98: 227-239.

DOI

[8]
王峰, 林皋, 郑保敬, 等. 非线性热传导问题的基于滑动Kriging插值的MLPG法[J]. 大连理工大学学报, 2014, 54(3): 339-344.

(WANG Feng, LIN Gao, ZHENG Bao-jing, et al. MLPG Method Based on Moving Kriging Interpolation for Solving Nonlinear Heat Conduction Problems[J]. Journal of Dalian University of Technology, 2014, 54(3): 339-344.(in Chinese))

[9]
吴泽艳, 郑保敬, 叶永, 等. 非线性热传导方程MLPG/RBF-FD无网格数值模拟[J]. 工程热物理学报, 2022, 43(1): 164-172.

(WU Ze-yan, ZHENG Bao-jing, YE Yong, et al. Numerical Simulation for the Nonlinear Heat Conduction Equations Based on MLPG/RBF-FD Meshless Method[J]. Journal of Engineering Thermophysics, 2022, 43(1): 164-172.(in Chinese))

[10]
李庆华, 冯子超, 陈莘莘, 等. 稳态非线性热传导问题的比例边界有限元法[J]. 华东交通大学学报, 2023, 40(6): 110-114.

(LI Qing-hua, FENG Zi-chao, CHEN Shen-shen, et al. Scaled Boundary Finite Element Method for Steady-state Nonlinear Heat Conduction Problem[J]. Journal of East China Jiaotong University, 2023, 40(6): 110-114.(in Chinese))

[11]
ZHANG L, GUO F, ZHENG H. The MLS-based Numerical Manifold Method for Nonlinear Transient Heat Conduction Problems in Functionally Graded Materials[J]. International Communications in Heat and Mass Transfer, 2022, 139: 106428.

DOI

[12]
李腊梅, 冯春. 一种非连续介质中热传导过程的数值模拟方法[J]. 工程力学, 2016, 33(1):25-31,46.

(LI La-mei, FENG Chun. A Numerical Simulation Method for Heat Conduction in Discontinuous Media[J]. Engineering Mechanics, 2016, 33(1): 25-31, 46.(in Chinese))

[13]
刘承论, 秦忠诚. 三维非稳态热传导问题的边界元法[J]. 岩石力学与工程学报, 2004, 23(18):3168-3173.

(LIU Cheng-lun, QIN Zhong-cheng. Boundary Element Method for 3D non-steady Heat Conduction[J]. Chinese Journal of Rock Mechanics and Engineering, 2004, 23(18): 3168-3173.(in Chinese))

[14]
梁钰, 郑保敬, 高效伟, 等. 基于POD模型降阶法的非线性瞬态热传导分析[J]. 中国科学:物理学力学天文学, 2018, 48(12):36-45.

(LIANG Yu, ZHENG Bao-jing, GAO Xiao-wei, et al. Reduced Order Model Analysis Method via Proper Orthogonal Decomposition for Nonlinear Transient Heat Conduction Problems[J]. Scientia Sinica (Physica, Mechanica & Astronomica), 2018, 48(12): 36-45.(in Chinese))

[15]
SUVIN V, OOI E T, SONG C, et al. Temperature-dependent Nonlinear Transient Heat Conduction Using the Scaled Boundary Finite Element Method[J]. International Journal of Heat and Mass Transfer, 2025, 243: 126780.

DOI

[16]
MIERZWICZAK M, CHEN W, FU Z J. The Singular Boundary Method for Steady-state Nonlinear Heat Conduction Problem with Temperature-dependent Thermal Conductivity[J]. International Journal of Heat and Mass Transfer, 2015, 91: 205-217.

DOI

[17]
YANG K, FENG W Z, WANG J, et al. RIBEM for 2D and 3D Nonlinear Heat Conduction with Temperature Dependent Conductivity[J]. Engineering Analysis with Boundary Elements, 2018, 87: 1-8.

DOI

[18]
KHOSRAVIFARD A, HEMATIYAN M R, MARIN L. Nonlinear Transient Heat Conduction Analysis of Functionally Graded Materials in the Presence of Heat Sources Using an Improved Meshless Radial Point Interpolation Method[J]. Applied Mathematical Modelling, 2011, 35(9): 4157-4174.

DOI

[19]
SHI G H. Manifold Method of Material Analysis[C]// U.S. Army Research Office. Transactions of the Ninth Army Conference on Applied Mathematics and Computing. Minneapolis, Minnesota, U.S.A, June 18-21,1991: 51-76.

[20]
ZHANG L, KONG H, ZHENG H. Numerical Manifold Method for Steady-state Nonlinear Heat Conduction Using Kirchhoff Transformation[J]. Science China Technological Sciences, 2024, 67(4): 992-1006.

DOI

[21]
贾真, 杨冬梅, 郑宏. 数值流形法渗流分析中处理弱不连续界面的新方法[J]. 长江科学院院报, 2024, 41(12): 133-137, 154.

DOI

(JIA Zhen, YANG Dong-mei, ZHENG Hong. A New Method to Deal with Weak Discontinuous Interface in Seepage Analysis Based on Numerical Manifold Method[J]. Journal of Changjiang River Scientific Research Institute, 2024, 41(12): 133-137, 154.(in Chinese))

[22]
李伟, 郑宏, 王海龙, 等. 求解断裂问题的新型无网格数值流形法[J]. 岩石力学与工程学报, 2020, 39(增刊1): 2655-2664.

(LI Wei, ZHENG Hong, WANG Hai-long, et al. A New Meshless Numerical Manifold Method for Solving Fracture Problems[J]. Chinese Journal of Rock Mechanics and Engineering, 2020, 39(Supp.1): 2655-2664.(in Chinese))

[23]
王方义, 郑宏. 无界域问题的数值流形法[J]. 长江科学院院报, 2023, 40(7): 110-117.

DOI

(WANG Fang-yi, ZHENG Hong. Numerical Manifold Method for Unbounded Domain Problems[J]. Journal of Changjiang River Scientific Research Institute, 2023, 40(7): 110-117.(in Chinese))

[24]
陈远强, 郑宏, 陈涛. 基于数值流形法的重力坝抗滑稳定性分析[J]. 长江科学院院报, 2016, 33(9): 133-137.

DOI

(CHEN Yuan-qiang, ZHENG Hong, CHEN Tao. Analysis of Anti-sliding Stability of Gravity Dam Using Numerical Manifold Method[J]. Journal of Yangtze River Scientific Research Institute, 2016, 33(9): 133-137.(in Chinese))

DOI

[25]
胡国栋, 张慧华, 谭育新. 功能梯度材料稳态热传导问题的数值流形方法研究[J]. 应用力学学报, 2017, 34(2): 311-317, 406.

(HU Guo-dong, ZHANG Hui-hua, TAN Yu-xin. Numerical Manifold Study of Steady Heat Conduction Problems in Functionally Graded Materials[J]. Chinese Journal of Applied Mechanics, 2017, 34(2): 311-317, 406.(in Chinese))

[26]
TAN F, TONG D, LIANG J, et al. Two-dimensional Numerical Manifold Method for Heat Conduction Problems[J]. Engineering Analysis with Boundary Elements, 2022,137:119-138.

[27]
WANG K, TANG C, QIAN X, et al. Numerical Manifold Method with Local Mesh Refinement for Thermo-mechanical Coupling Analysis in Rocks[J]. Computers and Geotechnics, 2025, 179: 107009.

DOI

[28]
LI C L, GUO D L, ZHANG H H. The Numerical Manifold Method for Piezoelectric Materials with Hole Flaws under Electro-mechanical Loadings[J]. Engineering Analysis with Boundary Elements, 2025, 173: 106149.

DOI

[29]
ZHANG Y, ZHENG H, LIN S. Random Structure Modeling of Soil and Rock Mixture and Evaluation of Its Permeability Using Three-dimensional Numerical Manifold Method[J]. Computers and Geotechnics, 2025, 180: 107089.

DOI

[30]
苏海东, 颉志强, 龚亚琦, 等. 基于独立覆盖的流形法的收敛性及覆盖网格特性[J]. 长江科学院院报, 2016, 33(2):131-136.

DOI

(SU Hai-dong, XIE Zhi-qiang, GONG Ya-qi, et al. Characteristics of Convergence and Cover Mesh in Numerical Manifold Method Based on Independent Covers[J]. Journal of Yangtze River Scientific Research Institute, 2016, 33(2): 131-136.(in Chinese))

DOI

[31]
WU W, WAN T, YANG Y, et al. Three-dimensional Numerical Manifold Formulation with Continuous Nodal Gradients for Dynamics of Elasto-plastic Porous Media[J]. Computer Methods in Applied Mechanics and Engineering, 2022, 388: 114203.

DOI

[32]
YANG Y, LI J, WU W. Modeling Wave Propagation Across Rock Masses Using an Enriched 3D Numerical Manifold Method[J]. Science China Technological Sciences, 2024, 67(3):835-852.

DOI

[33]
HU M, WANG Y, RUTQVIST J. On Continuous and Discontinuous Approaches for Modeling Groundwater Flow in Heterogeneous Media Using the Numerical Manifold Method: Model Development and Comparison[J]. Advances in Water Resources, 2015, 80: 17-29.

DOI

[34]
ZHENG H, LI W, DU X. Exact Imposition of Essential Boundary Condition and Material Interface Continuity in Galerkin-based Meshless Methods[J]. International Journal for Numerical Methods in Engineering, 2017, 110(7): 637-660.

DOI

[35]
YANG K, WANG J, DU J M, et al. Radial Integration Boundary Element Method for Nonlinear Heat Conduction Problems with Temperature-dependent Conductivity[J]. International Journal of Heat and Mass Transfer, 2017, 104: 1145-1151.

DOI

[36]
ZHANG L, YIN Y, ZHENG H, et al. Singularity Treatments in Transient Confined Seepage Using Numerical Manifold Method[J]. Engineering Analysis with Boundary Elements, 2025, 171: 106100.

DOI

[37]
YANG K, LI H Y, PENG H F, et al. New Interface Integration BEM for Solving Multi-medium Nonlinear Heat Transfer Problems[J]. Engineering Analysis with Boundary Elements, 2020, 117: 66-75.

DOI

Options
文章导航

/

〈 〉

*严禁使用本系统存储、处理、传输涉密信息*
版权所有 © 《长江科学院院报》编辑部
技术支持：北京玛格泰克科技发展有限公司