基于滑动Kriging插值的EFG-SBM求解含侧边界的稳态热传导问题

doi:10.16183/j.cnki.jsjtu.2020.215

基于滑动Kriging插值的EFG-SBM求解含侧边界的稳态热传导问题

王峰¹^,², 陈佳莉¹^,², 陈灯红^,³, 范勇¹^,², 李志远⁴, 何卫平¹^,²

1.三峡大学水利与环境学院,湖北宜昌 443002

2.三峡大学湖北省水电工程施工与管理重点实验室,湖北宜昌 443002

3.三峡大学土木与建筑学院,湖北宜昌 443002

4.大连理工大学建设工程学部,辽宁大连 116024

Element-Free Galerkin Scaled Boundary Method Based on Moving Kriging Interpolation for Steady Heat Conduction Analysis with Temperatures on Side-Faces

WANG Feng¹^,², CHEN Jiali¹^,², CHEN Denghong^,³, FAN Yong¹^,², LI Zhiyuan⁴, HE Weiping¹^,²

1. College of Hydraulic and Environmental Engineering, China Three Gorges University, Yichang 443002, Hubei, China

2. Hubei Key Laboratory of Construction and Management in Hydropower Engineering, China Three Gorges University, Yichang 443002, Hubei, China

3. College of Civil Engineering and Architecture, China Three Gorges University, Yichang 443002, Hubei, China

4. Faculty of Infrastructure Engineering, Dalian University of Technology, Dalian 116024, Liaoning, China

通讯作者: 陈灯红,男,副教授;E-mail:d.chen@ctgu.edu.cn.

责任编辑: 陈晓燕

收稿日期: 2020-07-8

基金资助:

国家自然科学基金(51809154)
湖北省教育厅科学技术研究项目(Q20181207)
湖北省高等学校优秀中青年科技创新团队计划(T2020005)
湖北省水电工程施工与管理重点实验室开放基金(2019KSD04)

Received: 2020-07-8

作者简介 About authors

王峰(1987-),男,山东省莱阳市人,副教授,从事无网格法及比例边界有限元法研究.

摘要

采用基于滑动Kriging插值的无单元伽辽金比例边界法(EFG-SBM)求解侧边界有温度载荷的稳态热传导问题,该方法通过无单元伽辽金法(EFG)和滑动Kriging插值离散环向边界.由于滑动Kriging插值形函数具备Kronecker delta函数插值特性,克服了移动最小二乘逼近难以直接准确施加本质边界条件的不足.作为一种新型的边界型无网格法,EFG-SBM兼有EFG法和比例边界有限元法(SBFEM)的优点.该方法继承了SBFEM的半解析特性,通过引入比例边界坐标系,可将偏微分控制方程环向离散,径向上解析求解.与传统的SBFEM相比,环向边界通过节点进行离散,前处理和后处理简便.通过数值算例可以看出,相比基于拉格朗日多项式的SBFEM,基于滑动Kriging插值的EFG-SBM计算精度更高.相比有限元法(FEM),该方法能更好地反映尖角处热奇异性以及无限域温度分布状态.

关键词： 无单元伽辽金比例边界法; 滑动Kriging插值; 热传导; 比例边界有限元法

Abstract

The element-free Galerkin scaled boundary method (EFG-SBM) based on moving Kriging (MK) interpolation is used to solve steady heat conduction problems with temperature loads on side-faces, in which the circumferential boundary is discretized based on MK interpolation and the element-free Galerkin (EFG) method. As the shape functions constructed from the MK interpolation possess the Kronecker delta interpolation property, the MK shape functions overcome the shortcomings of moving least squares (MLS) approximation which is difficult to impose essential boundary conditions directly and accurately. As a new boundary-type meshless method, EFG-SBM has advantages of the EFG and scaled boundary finite element method (SBFEM). This method inherits the semi-analytical property of SBFEM by introducing the scaled boundary coordinate system, in which the governing differential equations are weakened in the circumferential direction and can be solved analytically in the radial direction. Unlike the traditional SBFEM, the preprocessing and postprocessing processes of EFG-SBM are simplified since only the nodal data structure is required in the circumferential direction. Numerical examples show that the EFG-SBM based on MK interpolation can obtain a higher accuracy than the SBFEM based on Lagrange polynomials. Compared with the finite element method (FEM), this method can better characterize the thermal singularity at the sharp corner and the temperature distribution of the infinite region.

Keywords： element-free Galerkin scaled boundary method (EFG-SBM); moving Kriging (MK) interpolation; heat conduction; scaled boundary finite element method (SBFEM)

PDF (1618KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

王峰, 陈佳莉, 陈灯红, 范勇, 李志远, 何卫平. 基于滑动Kriging插值的EFG-SBM求解含侧边界的稳态热传导问题[J]. 上海交通大学学报（自然版）, 2021, 55(11): 1483-1492 doi:10.16183/j.cnki.jsjtu.2020.215

WANG Feng, CHEN Jiali, CHEN Denghong, FAN Yong, LI Zhiyuan, HE Weiping. Element-Free Galerkin Scaled Boundary Method Based on Moving Kriging Interpolation for Steady Heat Conduction Analysis with Temperatures on Side-Faces[J]. Journal of shanghai Jiaotong University, 2021, 55(11): 1483-1492 doi:10.16183/j.cnki.jsjtu.2020.215

稳态传热是指系统各点温度仅随空间位置变化,不随时间发生改变^[1].在水利工程中,拱坝封拱温度计算下游以年平均气温、上游以年平均水温作为边界条件,此时求出的坝体温度场为稳定温度场.

目前求解稳态传热问题主要依靠数值计算方法,常见的有:有限体积法、有限元法、边界元法、无网格法^[2,3,4,5,6]等.其中有限元法是应用最为广泛的数值方法,然而该方法在求解热奇异性(热流奇异或温度奇异)问题时,裂纹附近需要加密网格^[7].此外,对于涉及无限域的稳态热传导问题,有限元法需要将无限域截断为一系列有限区域,并施加一定的人工边界条件^[8],这带来了计算时间的增加.

比例边界有限元法(Scaled Boundary Finite Element Method, SBFEM)是1997年由Song等^[9]基于坐标变换提出的一种半解析半数值的方法,能精确高效地求解热奇异性、无限域热传导问题^[10,11].该方法在计算域引入相似中心,建立含有径向和环向的比例边界坐标系.经过相似变换,物理变量转化为域边界点和相似中心所连径向坐标变化表示,偏微分控制方程转化为一阶常微分方程组.只需通过有限元法离散物理域边界,空间求解维数降低一维.该方法集成了有限元法和边界元法的优点,在环向可以获得有限元法精度解,相比边界元法,无需基本解.然而,比例边界有限元法的计算精度受环向离散网格影响,传统的拉格朗日多项式基函数可能导致相邻节点之间的低阶连续性.为此,Lobatto 多项式^[12]、非均匀有理B样条^[13]已经被用于构造环向形函数,计算精度和效率明显提高.

为了避免网格影响,学者们提出了无网格法^[14,15],该方法非常适合求解涉及网格重构和畸变等问题.其中具有代表性的是无单元Galerkin法(Element-Free Galerkin Method, EFG)^[16]、无网格局部Petrov-Galerkin (Meshless Local Petrov-Galerkin, MLPG) 法^[17]、径向基点插值法(Radial Point Interpolation Method, RPIM)^[18].EFG法是借助于移动最小二乘 (Moving Least Squares, MLS) 法与Galerkin弱形式,构造近似场函数的一种离散方法.相比有限元法,计算精度和收敛速度明显提高,无需后处理便可获得光滑梯度场.将EFG法引入到比例边界坐标系环向边界,可获得新的边界型无网格法—无单元伽辽金比例边界法(Element-Free Galerkin Scaled Boundary Method, EFG-SBM)^[19].该方法只用节点离散环向边界,前处理简单.相比基于拉格朗日多项式的SBFEM,该方法形函数更加光滑连续,具有更高的计算精度和更快的收敛速度.

然而,与EFG法类似,EFG-SBM中的形函数通常是由移动最小二乘近似构造的,不具有Kronecker delta插值特性,不能直接施加本质边界条件,故形函数转化矩阵法^[19]被提出,该方法将虚节点值转换为实际节点值,方便了本质边界条件的施加.另外一种处理本质边界条件的方法是选取具有Kronecker delta插值特性的形函数,例如插值型移动最小二乘法^[20]、滑动Kriging插值^[21].滑动Kriging插值(Moving Kriging, MK)是基于变异函数理论和结构分析的线性无偏插值,直接通过节点进行插值,建立的形函数具有Kronecker delta插值性质,能直接施加本质边界条件.目前已被广泛应用于无网格法中,能显著提高计算效率^{[22,23,24,25,26]}.由于MK插值模型参数取值直接影响计算精度,Zheng等^[27]提出了模型相关参数与节点间距之间关系表达式.

本文将MK插值的EFG-SBM法扩展用于二维含侧边界的稳态传热问题求解,文中所有算例均无量纲.环向边界用节点进行离散,通过MK插值来构造形函数,计算域边界依据端点分段为不同的光滑线段.由于MK插值构建的形函数具有Kronecker delta插值性质,能直接施加本质边界条件.

1 MK插值

在EFG-SBM中,为了不增加角点试函数的光滑性,MK插值是沿着每条线段而不是整个域边界实现的.设问题域边界Γ由M条光滑边组成,N个离散节点s_i(i=1,2,…,N)分布在一条边上.设点s插值域内有n个节点s_k(k=1,2,…,n),则点s处温度近似函数T^h(s)通过MK插值逼近为

(1)T^h(s)=

\overset{n}{\sum_{k = 1}}

ϕ_k(s)T(s_k)=Φ(s)T

式中:给定的n个节点处的温度值向量为T= ${[T (s_{1}) T (s_{2}) \dots T (s_{n})]}^{T}$ .

形函数Φ(s)= $[ϕ_{1} (s) ϕ_{2} (s) \dots ϕ_{n} (s)]$ 为

(2)Φ(s)=p^T(s)A+r^T(s)B

p(s)为m次完备的单项式基,即

(3)p^T(s)=

[\begin{array}{l} p_{1} (s) & p_{2} (s) & \dots & p_{m} (s) \end{array}]

对于线性基,

(4)p^T(s)=[1 s], m=2

对于二次基,

(5)p^T(s)=[1 s s²], m=3

(6)A=(P^TR^-1P)^-1P^TR^-1

(7)B=R^-1(I-PA)

P为节点处基函数值矩阵;R为对角线为1的对称相关矩阵;I为n阶单位矩阵.

(8)P=

[\begin{array}{l} p_{1} (s_{1}) & p_{2} (s_{1}) & \dots & p_{m} (s_{1}) \\ p_{1} (s_{2}) & p_{2} (s_{2}) & \dots & p_{m} (s_{2}) \\ ︙ & ︙ & ⋱ & ︙ \\ p_{1} (s_{n}) & p_{2} (s_{n}) & \dots & p_{m} (s_{n}) \end{array}]

(9)R=

[\begin{array}{l} 1 & R (s_{1}, s_{2}) & \dots & R (s_{1}, s_{n}) \\ R (s_{2}, s_{1}) & 1 & \dots & R (s_{2}, s_{n}) \\ ︙ & ︙ & ⋱ & ︙ \\ R (s_{n}, s_{1}) & R (s_{n}, s_{2}) & \dots & 1 \end{array}]

矢量r(s)表达式为

(10)r^T(s)=

[\begin{array}{l} R (s_{1}, s) & R (s_{2}, s) & \dots & R (s_{n}, s) \end{array}]

R(s_i,s_j)为任意两节点s_i和s_j之间的对称相关矩阵,一般采用Gaussian型函数

(11)R(s_i,s_j)=

e^{- θ r_{ij}^{2}}

式中:θ>0,为函数相关参数;r_ij=‖s_i-s_j‖.

考虑[0,1]内的一组离散节点(0,0.25,0.5,0.75,1),取线性基,当θ取不同值时,s=0.5处所对应的形函数ϕ及形函数导数ϕ_,_s分别如图1和2所示,其中θ取值范围为0.1~500.由图1和2可以看出,参数θ对形函数及形函数导数影响很大.

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 MK形函数随参数θ的变化

Fig.1 ϕ versus θ

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 MK形函数导数随参数θ变化

Fig.2 ϕ_,_s versus θ

一维空间[0,1]通过6个节点(0,0.2,0.4,0.6,0.8,1)进行离散,图3所示为采用二次基时一维MK形函数,很明显,MK形函数具有Kronecker delta插值性质,即

(12)ϕ_k(s_j)=δ_kj=

\{\begin{array}{l} 1 & (k = j; k, j = 1,2, \dots, n) \\ 0 & (k \neq j; k, j = 1,2, \dots, n) \end{array}

同时也满足单位分解特性,即

(13)

\overset{n}{\sum_{k = 1}}

ϕ_k(s)=1

(14)

\overset{n}{\sum_{k = 1}}

ϕ_k(s)s_k=s

对于温度函数T(s),在曲线域进行插值,插值域半径r_i为

(15)r_i=βds

式中:β为系数;ds为点s附近节点距离.

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 采用二次基时一维MK形函数

Fig.3 One-dimensional MK shape functions with quadratic basis

2 侧边界含温度荷载的EFG-SBM方程

2.1 二维稳态热传导控制方程

设问题域边界Γ由温度边界Γ_D、热流密度边界Γ_N及对流换热边界Γ_R组成,当不考虑内热源,二维稳态热传导控制方程为

(16)Δ²T=0, 在Ω中

相应的边界条件为

(17)T=

\bar{T}

, 在Γ_D上

(18)-n·kΔT=

\bar{q}

, 在Γ_N上

(19)n·kΔT=h(T_a-T), 在Γ_R上

式中:Ω为计算域; $\bar{T}$ 和 $\bar{q}$ 分别为给定的温度和热流密度;n为边界Γ的单位外法线向量;k为热传导系数;h为表面换热系数;T_a为周围环境温度.简便起见,这里并没有考虑对流换热边界Γ_R作用.

2.2 侧边界施加温度荷载求解

如图4所示,引入相似中心为(x₀,y₀)的比例边界坐标系oξs.ξ为径向坐标,从相似中心连至边界点,相似中心ξ值为0,边界上ξ值为1.s为环向坐标,表示沿着边界Γ从s₀开始逆时针旋转的弧长,终点为s₁.由式(16)可得

(20)E⁰ξ²T(ξ)_,_ξξ+(E⁰-E¹+(E¹)^T)ξT(ξ)_,_ξ-E²T(ξ)=0

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 比例边界坐标系

Fig.4 Scaled boundary coordinate system

式中:

(21)E⁰=∫_SB¹(s)^TB¹(s)

|J|

(22)E¹=∫_SB²(s)^TB¹(s)

|J|

(23)E²=∫_SB²(s)^TB²(s)

|J|

(24)

\begin{array}{l} B^{1} (s) = b^{1} (s) Φ (s) \\ B^{2} (s) = b^{2} (s) Φ {(s)}_{, s} \end{array}\}

(25)b¹(s)=

\frac{1}{|J|} [\begin{array}{l} y_{s} {(s)}_{, s} \\ - x_{s} {(s)}_{, s} \end{array}]

(26)b²(s)=

\frac{1}{|J|} [\begin{array}{l} - (y_{s} (s) - y_{0}) \\ x_{s} (s) - x_{0} \end{array}]

边界(ξ=1)上Jacobian矩阵为

(27)J=

[\begin{array}{l} x_{s} (s) - x_{0} & y_{s} (s) - y_{0} \\ x_{s} {(s)}_{, s} & y_{s} {(s)}_{, s} \end{array}]

Jacobian行列式为

(28)

|J|

(x_{s} (s) - x_{0})

y_s(s)_,_s-

(y_{s} (s) - y_{0})

x_s(s)_,_s

(x_s(s),y_s(s))为边界ξ=1上点的笛卡尔坐标.

相似中心和起点s₀、终点s₁的连线为侧边界,当侧边界s=s₀或s=s₁上温度不为0时,此时温度函数T(ξ)分解为约束温度函数 ${\bar{T}}_{s}$ (ξ)(自由度数为n_s)和未约束温度函数T_u(ξ)(自由度数为n_u)^[28^-^29].此时比例边界方程(20)分解为

(29)

[\begin{array}{l} E_{uu}^{0} & E_{us}^{0} \\ E_{su}^{0} & E_{ss}^{0} \end{array}]

ξ²

{[\begin{array}{l} T_{u} (ξ) \\ {\bar{T}}_{s} (ξ) \end{array}]}_{, ξξ}

[\begin{array}{l} E_{uu}^{0} - E_{uu}^{1} + (E_{uu}^{1})^{T} & E_{us}^{0} - E_{us}^{1} + (E_{us}^{1})^{T} \\ E_{su}^{0} - E_{su}^{1} + (E_{su}^{1})^{T} & E_{ss}^{0} - E_{ss}^{1} + (E_{ss}^{1})^{T} \end{array}]

{[\begin{array}{l} T_{u} (ξ) \\ {\bar{T}}_{s} (ξ) \end{array}]}_{, ξ}

[\begin{array}{l} E_{uu}^{2} & E_{us}^{2} \\ E_{su}^{2} & E_{ss}^{2} \end{array}] [\begin{array}{l} T_{u} (ξ) \\ {\bar{T}}_{s} (ξ) \end{array}]

式(29)的第一行方程可得

(30)

E_{uu}^{0}

ξ²T_u(ξ)_,_ξξ+

(E_{uu}^{0} - E_{uu}^{1} + (E_{uu}^{1})^{T})

ξT_u(ξ)_,_ξ-

E_{uu}^{2}

T_u(ξ)=-

E_{us}^{0}

ξ²

{\bar{T}}_{s}

(ξ)_,_ξξ-(

E_{us}^{0}

E_{us}^{1}

E_{us}^{1}

)^T)ξ

{\bar{T}}_{s}

(ξ)_,_ξ+

E_{us}^{2} {\bar{T}}_{s}

(ξ)

假设 ${\bar{T}}_{s}$ (ξ)可以展成幂级数的形式,即

(31)

{\bar{T}}_{s}

(ξ)=

ξ^{p_{i}} {\bar{T}}_{sp}

式(30)为非齐次常微分方程,该方程的通解是对应齐次线性方程的通解加非齐次线性方程的特解.式(30)的特解也可以设置成幂级数的形式,即

(32)

T_{u}^{*}

(ξ)=

ξ^{p_{i}} T_{up}^{*}

代入式(30)可得

(33)

T_{up}^{*}

{[p_{i}^{2} E_{uu}^{0} - p_{i} ((E_{uu}^{1})^{T} - E_{uu}^{1}) - E_{uu}^{2}]}^{- 1} [- p_{i}^{2} E_{us}^{0} + p_{i} ((E_{us}^{1})^{T} - E_{us}^{1}) + E_{us}^{2}] {\bar{T}}_{sp}

对于有限域,式(33)的解表示为

(34)T(ξ)=

[\begin{array}{l} T_{u} (ξ) \\ {\bar{T}}_{s} (ξ) \end{array}]

[\begin{array}{l} \overset{n_{u}}{\sum_{i = 1}} ξ^{- λ_{i}} T_{u} c_{1 i} \\ 0 \end{array}]

ξ^{p_{i}} [\begin{array}{l} T_{up}^{*} \\ {\bar{T}}_{sp} \end{array}]

式中:λ_i为特征值分解所得特征值;c₁_i为积分常数.

对于有限域,对偶变量q(ξ)可以分为

(35)

[\begin{array}{l} q_{u} (ξ) \\ q_{s} (ξ) \end{array}]

[\begin{array}{l} E_{uu}^{0} & E_{us}^{0} \\ (E_{us}^{0})^{T} & E_{ss}^{0} \end{array}]

{[\begin{array}{l} T_{u} (ξ) \\ {\bar{T}}_{s} (ξ) \end{array}]}_{ξ}

[\begin{array}{l} (E_{uu}^{1})^{T} & (E_{su}^{1})^{T} \\ (E_{us}^{1})^{T} & (E_{ss}^{1})^{T} \end{array}] [\begin{array}{l} T_{u} (ξ) \\ {\bar{T}}_{s} (ξ) \end{array}]

在边界ξ=1处,将式(34)代入式(35)可得

(36)

[\begin{array}{l} q_{u} (ξ) \\ q_{s} (ξ) \end{array}]

\overset{n_{u}}{\sum_{i = 1}}

c₁_i

[\begin{array}{l} q_{u i} \\ q_{s i} \end{array}]

[\begin{array}{l} q_{up}^{*} \\ q_{sp}^{*} \end{array}]

式中:

(37)q_u_i=

((E_{uu}^{1})^{T} - λ_{i} E_{uu}^{0}) T_{up}^{*}

(38)q_s_i=

((E_{us}^{1})^{T} - λ_{i} (E_{us}^{0})^{T}) T_{up}^{*}

(39)

q_{up}^{*}

((E_{uu}^{1})^{T} + p_{i} E_{uu}^{0}) T_{up}^{*}

((E_{su}^{1})^{T} + p_{i} E_{us}^{0}) {\bar{T}}_{sp}

(40)

q_{sp}^{*}

((E_{us}^{1})^{T} + p_{i} E_{us}^{0}) T_{up}^{*}

((E_{ss}^{1})^{T} + p_{i} E_{ss}^{0}) {\bar{T}}_{sp}

此时q_u(ξ)的通解为

(41)q_u(ξ)=

{\tilde{Φ}}_{21} ξ^{- {\tilde{Λ}}_{0}}

c₁+

ξ^{p_{i}} q_{up}^{*}

在边界ξ=1处,

(42)T_u=

{\tilde{Φ}}_{11}

c₁+

T_{up}^{*}

(43)q_u=

{\tilde{Φ}}_{21}

c₁+

q_{up}^{*}

其中, $\tilde{Φ}$ 和 ${\tilde{Λ}}_{0}$ 为“划行划列”处理后的系数矩阵所形成的Hamiltonian矩阵 $\tilde{Z}$ 的特征系统.

由式(42)可得

(44)c₁=

{\tilde{Φ}}_{11}^{- 1} (T_{u} - T_{up}^{*})

代入式(43)可得

(45)q_u=

{\tilde{Φ}}_{21}

c₁+

q_{up}^{*}

{\tilde{Φ}}_{21} {\tilde{Φ}}_{11}^{- 1} (T_{u} - T_{up}^{*})

q_{up}^{*}

同样,对于无限域,可得

(46)q_u=-

{\tilde{Φ}}_{22} {\tilde{Φ}}_{12}^{- 1} (T_{u} - T_{up}^{*})

q_{up}^{*}

3 数值算例

3.1 矩形板稳态热传导问题

为了验证EFG-SBM求解侧边界含温度荷载的准确性,考虑如图5所示的矩形板^[30],图中:a为矩形长度.控制方程为

(47)Δ²T=0

该问题的温度解析解为

(48)T(x,y)=

\frac{100}{\sinh (0.6 π)}

sin

(\frac{π}{a} x)

sinh

(\frac{π}{a} y)

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 矩形板热传导边界条件

Fig.5 Boundary conditions of heat conduction for rectangular plate

由于问题对称,只选取矩形板左半部分进行计算.计算过程中取a=10,β=4,相似中心选在左下角,如图6所示,此时左侧和底部是侧边界,无需通过节点进行离散.图7所示为相对均方根误差随θ变化趋势图,可以看出,θ=0.05~0.12可得到较理想的计算效果.图8所示为θ=0.1时温度等值线分布图,显然,温度等值线垂直于矩形板右侧边界,这跟矩形板右侧边界热流密度为0是相符的.矩形板左侧和底部边界为0,跟边界条件吻合.

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 矩形板相似中心和节点分布

Fig.6 Scaling center and nodes distribution for rectangular plate

图7

新窗口打开| 下载原图ZIP| 生成PPT

图7 相对均方根误差随θ的变化

Fig.7 Relative root mean square error versus θ

图8

新窗口打开| 下载原图ZIP| 生成PPT

图8 矩形板温度等值线图

Fig.8 Temperature isothermal line of rectangular plate

为了研究本文方法计算精度,给出了EFG-SBM和SBFEM数值计算结果比较(θ=0.1),见表1.从表1可以看出,EFG-SBM具有更高的计算精度,最大相对误差为0.005%,误差比SBFEM误差低一个量级.这是由于MK形函数高阶连续,插值域内可以共享更多节点信息.由此可以看出EFG-SBM的优越性,继承了SBFEM和EFG法的优点,只用节点离散环向边界,然后通过EFG法进行求解,形函数高阶连续以及贡献更多节点信息保证了计算精度更高,径向解析可求得域内点温度值.计算结果表明,EFG-SBM求解含侧边界的稳态传热问题是可行的.

表1 矩形板稳态热传导问题理论解与数值解比较

Tab.1 Comparison of theoretical and numerical solutions for steady heat conduction problem of rectangular plate

计算点 (x,y)	理论解	SBFEM解	EFG- SBM解	SBFEM 误差/%	EFG-SBM 误差/%
(5,1)	9.927	9.971	9.927	0.45	0
(5,2)	20.841	20.823	20.842	0.09	0.005
(5,3)	33.830	33.804	33.830	0.08	0
(5,4)	50.184	50.155	50.185	0.06	0.002
(5,5)	71.533	71.509	71.533	0.03	0

新窗口打开| 下载CSV

3.2 L型板稳态热传导问题

考虑如图9所示的L型板,其在尖角处具有热奇异性,左侧和底部绝热.EFG-SBM的节点分布如图10所示,相似中心选在坐标原点处,黑虚线是侧边界,选用线性基.为了比较计算结果,这里选用有限元方法,基于二维四节点拉格朗日矩形单元,共划分为 1200 个单元. $\bar{T}$ 考虑两种温度荷载形式:

(49)

{\bar{T}}_{1}

=50+50ξ

(50)

{\bar{T}}_{2}

=50+40ξ+15ξ²-5ξ³

图9

新窗口打开| 下载原图ZIP| 生成PPT

图9 L型区域热传导

Fig.9 Heat conduction in L-shaped domain

图10

新窗口打开| 下载原图ZIP| 生成PPT

图10 L型板相似中心和节点分布

Fig.10 Scaling center and nodes distribution for L-shaped plate

图11分别给出了 ${\bar{T}}_{1}$ 温度荷载下EFG-SBM和有限元方法计算得到的L型板温度场分布图,很明显,温度结果在整个计算域吻合很好.图12和13分别给出了L型板V型凹角温度和热流(q)分布图,可以看出,EFG-SBM求解得到的温度和热流等值线更加光滑平顺, 而有限元方法由于无法反映尖角处的热奇异性,温度和热流等值线不光滑甚至出现折线,温度和热流等值线存在较大误差.综合图11~13,EFG-SBM和有限元方法仅在V型凹角处存在细微差异,而EFG-SBM仅需通过节点离散模型边界,更加简便,计算更加高效.

图11

新窗口打开| 下载原图ZIP| 生成PPT

图11 ${\bar{T}}_{1}$ 温度荷载下L型板温度等值线图

Fig.11 Temperature isothermal line of L-shaped plate at ${\bar{T}}_{1}$

图12

新窗口打开| 下载原图ZIP| 生成PPT

图12 ${\bar{T}}_{1}$ 温度荷载下L型板V型凹角温度等值线图

Fig.12 Temperature isothermal line of L-shaped plate with V-shaped reentrant corner at ${\bar{T}}_{1}$

图13

新窗口打开| 下载原图ZIP| 生成PPT

图13 ${\bar{T}}_{1}$ 温度荷载下L型板V型凹角热流等值线图

Fig.13 Flux isothermal line of L-shaped plate with V-shaped reentrant corner at ${\bar{T}}_{1}$

图14所示为 ${\bar{T}}_{2}$ 温度荷载下EFG-SBM求解得到的L型板温度场分布图,可以看出,从 ${\bar{T}}_{1}$ 变化到 ${\bar{T}}_{2}$ ,温度变化趋势总体一致,温度荷载变大,影响到的区域范围越大,甚至能影响到L型板左侧和底部.

图14

新窗口打开| 下载原图ZIP| 生成PPT

图14 ${\bar{T}}_{2}$ 温度荷载下L型板温度等值线图

Fig.14 Temperature isothermal line of L-shaped plate at ${\bar{T}}_{2}$

3.3 半无限域稳态热传导问题

考虑如图15所示的半无限域,侧边界温度为50,边界温度为

(51)

\bar{T}

=100-5ξ

图15

新窗口打开| 下载原图ZIP| 生成PPT

图15 半无限域热传导

Fig.15 Heat conduction in a semi-infinite domain

EFG-SBM的节点分布如图16所示,相似中心选在坐标原点处,选用线性基.为了对计算结果进行比较,这里选用有限元方法离散一边长为20的区域,基于二维四节点拉格朗日矩形单元,单位长度上布置5个单元,此时节点总数共计 7701,共划分为 7500 个单元.图17所示为通过EFG-SBM和有限元方法计算得到的温度场分布图.很明显,有限元方法由于无法刻画无限域需截断边界,受截断区域影响,温度计算结果会产生不同.越接近截断边界,有限元结果跟EFG-SBM结果差别越大,EFG-SBM计算得到的温度等值线更加光滑.

图16

新窗口打开| 下载原图ZIP| 生成PPT

图16 半无限域相似中心和节点分布

Fig.16 Scaling center and nodes distribution for semi-infinite domain

图17

新窗口打开| 下载原图ZIP| 生成PPT

图17 截断区域温度等值线图

Fig.17 Temperature isothermal line of truncated domain

4 结论

(1) 传统移动最小二乘法形函数缺乏插值特性,将滑动Kriging插值法引入到EFG-SBM进行改进.由于滑动Kriging插值形函数具有Kronecker delta函数插值性质,在本质边界条件施加上无需其他处理.同时相比SBFEM中常用的拉格朗日多项式,滑动Kriging插值建立的形函数更加连续光滑.

(2) 将基于滑动Kriging插值的EFG-SBM用于求解二维含侧边界的稳态传热问题,计算域边界根据端点分段为不同光滑线段,每条线段通过滑动Kriging插值来构造弧长以及温度场函数.通过3个典型算例来说明EFG-SBM的特点,首先,在同等自由度离散情况下,相比SBFEM,EFG-SBM计算精度更高.其次,对于V型凹角结构,相比有限元法,EFG-SBM能更好地反映尖角处热奇异性(温度或热流奇异).最后,通过对半无限域稳态传热问题进行求解,有限元法截断边界对温度分布状态具有较大影响,实际问题求解应予以重视,而EFG-SBM可以很方便地刻画无限域.

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

杨世铭, 陶文铨. 传热学[M]. 北京: 高等教育出版社, 2006.