CN107563038B - 一种接触热阻有限元求解方法 - Google Patents

一种接触热阻有限元求解方法 Download PDF

Info

Publication number
CN107563038B
CN107563038B CN201710747272.1A CN201710747272A CN107563038B CN 107563038 B CN107563038 B CN 107563038B CN 201710747272 A CN201710747272 A CN 201710747272A CN 107563038 B CN107563038 B CN 107563038B
Authority
CN
China
Prior art keywords
contact
finite element
thermal resistance
regions
contact thermal
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201710747272.1A
Other languages
English (en)
Other versions
CN107563038A (zh
Inventor
徐立
谢鹏
李斌
杨中海
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
University of Electronic Science and Technology of China
Original Assignee
University of Electronic Science and Technology of China
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by University of Electronic Science and Technology of China filed Critical University of Electronic Science and Technology of China
Priority to CN201710747272.1A priority Critical patent/CN107563038B/zh
Publication of CN107563038A publication Critical patent/CN107563038A/zh
Application granted granted Critical
Publication of CN107563038B publication Critical patent/CN107563038B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Landscapes

  • Investigating Or Analyzing Materials Using Thermal Means (AREA)

Abstract

本发明属于三维热传导有限元数值求解技术领域，涉及一种新型的接触热阻有限元求解方法。本发明首先对存在接触问题的器件进行建模，然后开创性地将接触热阻问题作为一种边界条件引入热传导问题，并采用伽辽金残数加权的方法，得到接触热阻问题的有限元弱形式。接着采用四面体网格剖分模型，选择叠层基函数，离散有限元弱形式方程，得到有限元单元矩阵和右端向量，然后集成总的有限元线性方程组，经过直接法或者迭代法进行求解，从而快速准确地得到最终的数值计算结果。

Description

一种接触热阻有限元求解方法
技术领域
本发明属于三维热传导有限元数值求解技术领域，涉及一种接触热阻有限元求解方法。
背景技术
接触热阻(Thermal Contact Resistance)问题已经涉及到航天、机械制造、微电子、化工、低温超导、生物医学、核反应堆以及仪器仪表等众多科学与工程领域，其产生机理，广大学者普遍认为是由于粗糙表面间的不完全接触所造成的热流收缩而导致的。从理论上讲，完全接触的接触面应该保持同一温度，而在实际工程中，任何表面在微观上都是粗糙的，故任何固体表面之间都不可能完全接触，接触的地方直接导热，在不接触处存在空隙，产生热流收缩，存在传热阻力，即接触热阻。
实际中，由于接触热阻的影响，对物体进行热分析时会造成一定的误差，甚至得到完全错误的结果。近些年来随着计算机的不断发展，人们开始利用计算机进行数值分析研究，并且将会有越来越多的研究人员使用有限元方法进行接触热阻的数值模拟分析，数值分析可以为接触热阻实验提供很好的借鉴意义。到目前为止，接触热阻在使用有限元方法进行数值模拟分析方面发表的文献不多，主要是把界面微凸体等效成长方体、圆柱体、圆锥体等，或者是在两个物体的接触面上施加一个薄层来进行过渡。
以上的这些方法因为要改变原模型的结构，效率会很低下，误差也难以控制，甚至可能会得到完全错误的结果，在实际应用中会很难实施，尤其是对于形状不规则以及规模大的模型。
发明内容
针对上述存在问题或不足，为解决现有方法在实际应用中效率低下，误差难以控制，甚至得到完全错误的结果和难以实施的问题，本发明提供了一种接触热阻有限元求解方法，该方法将接触热阻问题转化为有限元边界条件进行处理，高效简便地解决了接触热阻难题。
其具体技术方案，包括以下步骤：
A.针对热分析的对象进行建模，建立对应的几何结构模型；
B.采用四面体网格对模型进行剖分；
C.将接触热阻转化为边界条件，用伽辽金残数加权法得到有限元的弱形式；
对于单位面积的交界面，接触热阻定义如下：
Figure GDA0002289098900000011
其中R表示接触热阻，uA、uB表示接触面两侧温度，q”表示平均热流密度，文字表述为：接触热阻等于两个接触面温度之差除以平均热流密度。
通过伽辽金方法，将接触热阻问题转化为边界条件之后，可以得到最终的有限元弱形式为
Figure GDA0002289098900000021
其中
Figure GDA0002289098900000022
为拉普拉斯算子，k为热传导系数，v、v1和v2为权函数，u、u1和u2为温度，δc为接触热导(接触热阻的倒数)，Ω为整个求解域，Γ12和Γ21为接触边界。
D.用叠层基函数进行目标离散，得到最终的有限元方程组。
E.求解步骤D中矩阵和右端项形成的线性方程组，得到最终的温度解。
步骤C中，将接触热阻问题转化为热传导问题的边界条件来进行处理是本发明的最优之处，通过考虑接触热阻的定义以及在接触面上的热通量对于接触面两侧来说是相同的(接触面两侧的温度值是不同的)这两点，开创性地将其转化为有限元边界条件，进而用有限元方法进行处理。
综上所述，本发明从接触热阻本质出发，虽然接触导致接触面两侧的温度不一致，但是接触面上的热通量是连续的，这里将接触热阻问题等效为热传导边值问题，用有限元方法进行求解。相比于以前的接触热阻处理方法，本发明最大的优势就是在不改变原模型的基础上，非常规范、高效地对接触热阻问题进行了求解，在实际工程领域非常易于实施。
附图说明
图1是本发明的流程图；
图2是接触热阻作为有限元边界条件进行处理的示意图；
图3是四节点四面体单元示意图；
图4是接触单元计算示意图。
具体实施方式
下面结合附图和具体实施例来详细描述本发明的技术方案。
一种接触热阻有限元计算方法，包括以下步骤：
A.针对热分析的对象进行建模，建立对应的几何结构模型。
B.采用四面体网格对模型进行剖分；剖分后的求解域被分割成三维四面体网格，从而将连续的几何结构空间转化为离散的网格空间。
C.将接触热阻转化为边界条件，用伽辽金残数加权法得到有限元的弱形式。
考虑一个域
Figure GDA0002289098900000031
(R3表示三维空间)，热传导的控制微分方程由下式给出:
Figure GDA0002289098900000032
其中，u代表在有限区域上随时间变化的温度分布，
Figure GDA0002289098900000033
为拉普拉斯算子，k是热传导系数，Q是内部产热量，ρ是密度，c是比热容，t是时间，Ω表示求解域。
为了说明问题，这里将求解域分为两个区域，进而考虑两个区域接触面上的接触热阻。就如图2中所示的那样：Ω:＝Ω1∪Ω2，将表示子域的接触面表示为
Figure GDA0002289098900000034
外部边界定义为
Figure GDA0002289098900000035
另外，
Figure GDA0002289098900000036
表示单位法向量，方向是在边界Γ12上指向区域Ωi外部。
为了更好地描述，这里定义体积分和面积分:
(u,v)Ω＝∫Ω(u,v)dV (4)
<u,v>Γ＝∫Γ(u,v)dS (5)
其中u、v表示任意两个函数，V表示体积，S表示面积。
由接触热阻的定义，对两个区域分别有：
Figure GDA0002289098900000037
Figure GDA0002289098900000038
其中δc是两个区域交界面上的接触热导值(接触热阻的倒数)，k1,k2为两个区域的热传导系数，u1,u2为接触面上的温度，n1,n2为法线。考虑稳态热分析过程，即去掉(3)中的时间项，也不讨论内部产热项Q。由伽辽金方法可以得到稳态热分析控制微分方程的残差表达式：
Figure GDA0002289098900000039
其中RΩ表示稳态热分析控制微分方程的残差。同样地，两个区域的残差表达式可以写为：
Figure GDA00022890989000000310
Figure GDA00022890989000000311
其中Rc1,Rc2表示两个区域上接触热导表达式的残差。由伽辽金残数加权法，根据(8)(9)(10)式，可以写为如下的表达式
Figure GDA0002289098900000041
其中c1,c2是任意实数，v,v1,v2为权函数。由格林公式，
Figure GDA0002289098900000042
可以展开为
Figure GDA0002289098900000043
同样地，在接触面上，(11)式后两项可以展开为
Figure GDA0002289098900000044
Figure GDA0002289098900000045
由于c1,c2,v1,v2的任意性，令c1＝c2＝-1，v＝v1＝v2，可以将(11)式写为如下的表达式，也就是最终得到的有限元伽辽金弱形式。
Figure GDA0002289098900000046
D.用叠层基函数进行离散，得到最终的有限元方程组。
如图3所示四面体单元中i,j,k,l代表四个顶点的编号，我们得到四个最基本的基函数，具体推导过程为一种公知过程，这里不再阐述：
Figure GDA0002289098900000047
Figure GDA0002289098900000048
Figure GDA0002289098900000049
Figure GDA00022890989000000410
式中
Figure GDA00022890989000000411
Figure GDA00022890989000000412
Figure GDA0002289098900000051
Figure GDA0002289098900000052
Figure GDA0002289098900000053
将(20)式、(21)式、(22)式、(23)式中的i,j,k,l轮换，得到aj,ak,al，bj,bk,bl，cj,ck,cl，dj,dk,dl。V为四面体的体积。为了使四面体的体积不为负值，单元节点编号i,j,k,l必须依照一定的顺序。在右手坐标系中，当i,j,k的方向转动时，右手螺旋应指向l的方向。
确定了基函数为N1,N2,N3,N4，对于有限元过程来说，把域Ω离散为M个单元之后，如同(15)所示的弱形式定积分，可以通过将每个单元的积分贡献简单相加，即
Figure GDA0002289098900000054
由(25)式，定义单元矩阵为Sij＝Kij+Cij，其中Kij为
Figure GDA0002289098900000055
矩阵K的计算在很多基础的热传导有限元书中都有介绍，本发明不再赘述，主要介绍接触热阻矩阵项C，以图4所示的两个四面体为例进行说明。两个四面体为abcd和efgh，三角面bcd和efg是同一个面，即其接触面。C矩阵可以写为如下的表达式
Figure GDA0002289098900000056
其中Γbcd、Γefg分别表示三角面bcd和三角面efg，
Figure GDA0002289098900000057
代表三角面bcd的基函数，
Figure GDA0002289098900000058
代表三角面efg的基函数。按照上面的公式就可以建立接触面两侧之间的关系，形成最终的矩阵。
E.求解步骤D中矩阵和右端项形成的线性方程组，得到最终的温度解。
在步骤D中已经求得有限元的矩阵和右端向量，可以形成最终的线性方程组，有很多方法可以进行方程组的求解，直接法和迭代法都在很多文献资料里有详细的论述，本发明不再进行具体描述。求解完有限元方程组，就得到了最终的温度解，可以发现接触面两侧温度值出现了跳跃，接触热阻问题得以解决。
综上所述，可见本发明针对先前接触热阻处理方法的弊端，将接触热阻问题等效为边值问题，用有限元方法进行数值求解，在不改变原模型的基础上，非常规范、高效地对接触热阻问题进行了求解，在实际工程领域非常易于实施。

Claims (1)

1.一种接触热阻有限元求解方法，包括以下步骤：
A.针对热分析的对象进行建模，建立对应的几何结构模型；
B.采用四面体网格对模型进行剖分；
C.将接触热阻转化为边界条件，用伽辽金残数加权法得到有限元的弱形式；
对于单位面积的交界面，接触热阻定义如下：
Figure FDA0002289098890000011
其中R表示接触热阻，uA、uB表示接触面两侧温度，q”表示平均热流密度，文字表述为：接触热阻等于两个接触面温度之差除以平均热流密度；
考虑一个域
Figure FDA0002289098890000012
热传导的控制微分方程由下式给出:
Figure FDA0002289098890000013
其中，u代表在有限区域上随时间变化的温度分布，
Figure FDA0002289098890000014
为拉普拉斯算子，k是热传导系数，Q是内部产热量，ρ是密度，c是比热容，t是时间，Ω表示求解域；
将求解域分为两个区域，进而考虑两个区域接触面上的接触热阻，首先定义体积分和面积分:
(u,v)Ω＝∫Ω(u,v)dV (3)
<u,v>Γ＝∫Γ(u,v)dS (4)
其中u、v表示任意两个函数，V表示体积，S表示面积；由接触热阻的定义，对两个区域分别有：
Figure FDA0002289098890000015
Figure FDA0002289098890000016
其中δc是两个区域交界面上的接触热导值，k1,k2为两个区域的热传导系数，u1,u2为接触面上的温度，n1,n2为法线；
考虑稳态热分析过程，即去掉(2)中的时间项，也不讨论内部产热项Q；由伽辽金方法可以得到稳态热分析控制微分方程的残差表达式：
Figure FDA0002289098890000021
其中RΩ表示稳态热分析控制微分方程的残差，同样地，两个区域的残差表达式可以写为：
Figure FDA0002289098890000022
Figure FDA0002289098890000023
其中Rc1,Rc2表示两个区域上接触热导表达式的残差，由伽辽金残数加权法，根据(7)(8)(9)式，可以写为如下的表达式
Figure FDA0002289098890000024
其中c1,c2是任意实数，v,v1,v2为权函数，Γ12和Γ21为接触边界；由格林公式，
Figure FDA0002289098890000025
可以展开为
Figure FDA0002289098890000026
同样地，在接触面上，(10)式后两项可以展开为
Figure FDA0002289098890000027
Figure FDA0002289098890000028
由于c1,c2,v1,v2的任意性，令c1＝c2＝-1，v＝v1＝v2，可以将(10)式写为如下的表达式，也就是最终得到的有限元伽辽金弱形式：
Figure FDA0002289098890000029
D.用叠层基函数进行目标离散，得到最终的有限元方程组；
E.求解步骤D中矩阵和右端项形成的线性方程组，得到最终的温度解。
CN201710747272.1A 2017-08-28 2017-08-28 一种接触热阻有限元求解方法 Active CN107563038B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710747272.1A CN107563038B (zh) 2017-08-28 2017-08-28 一种接触热阻有限元求解方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710747272.1A CN107563038B (zh) 2017-08-28 2017-08-28 一种接触热阻有限元求解方法

Publications (2)

Publication Number Publication Date
CN107563038A CN107563038A (zh) 2018-01-09
CN107563038B true CN107563038B (zh) 2020-05-12

Family

ID=60977094

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710747272.1A Active CN107563038B (zh) 2017-08-28 2017-08-28 一种接触热阻有限元求解方法

Country Status (1)

Country Link
CN (1) CN107563038B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108535313B (zh) * 2018-02-11 2021-01-29 中国矿业大学 一种用热线法测量两固体之间界面热阻的方法
CN108828003A (zh) * 2018-06-09 2018-11-16 安徽华兴车辆有限公司 一种固体材料接触热阻测量装置及测量方法
CN109657286A (zh) * 2018-11-28 2019-04-19 电子科技大学 一种有限元瞬态热分析区域分解求解方法
CN112836399B (zh) * 2020-11-20 2022-11-08 电子科技大学 一种基于有限元算法的非线性接触热阻热分析求解方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Thermal Analysis of High-Power Integrated Circuits and Packages Using Nonconformal Domain Decomposition Method;Yang Shao等;《IEEE TRANSACTIONS ON COMPONENTS, PACKAGING AND MANUFACTURING TECHNOLOGY》;20130831;第3卷(第8期);1321-1331 *
Thermal-aware DC IR-drop co-analysis using non-conformal domain decomposition methods;Yang Shao等;《Preceeding of the royal society》;20120215;1652-1675 *
热传导分析的加权残数方法;陈建军等;《西北电讯工程学院学报》;19870930;第14卷(第3期);31-37 *

Also Published As

Publication number Publication date
CN107563038A (zh) 2018-01-09

Similar Documents

Publication Publication Date Title
CN107563038B (zh) 一种接触热阻有限元求解方法
Gao et al. Free element collocation method: A new method combining advantages of finite element and mesh free methods
Gao et al. Element differential method for solving general heat conduction problems
Martínez et al. Orthotropic k-nearest foams for additive manufacturing
Zhang et al. A fast meshless method based on proper orthogonal decomposition for the transient heat conduction problems
Gao et al. Element differential method and its application in thermal‐mechanical problems
Sutradhar et al. Transient heat conduction in homogeneous and non-homogeneous materials by the Laplace transform Galerkin boundary element method
Chen et al. Fast finite difference approximation for identifying parameters in a two-dimensional space-fractional nonlocal model with variable diffusivity coefficients
Ooi et al. Construction of high‐order complete scaled boundary shape functions over arbitrary polygons with bubble functions
Zhang et al. Extract of stress intensity factors on honeycomb elements by the numerical manifold method
Zhang et al. Modeling 2D transient heat conduction problems by the numerical manifold method on Wachspress polygonal elements
Yang et al. A simple Galerkin meshless method, the fragile points method using point stiffness matrices, for 2D linear elastic problems in complex domains with crack and rupture propagation
Zhang et al. An improved meshless method with almost interpolation property for isotropic heat conduction problems
Ramos et al. An extension of the Hill–Mandel principle for transient heat conduction in heterogeneous media with heat generation incorporating finite RVE thermal inertia effects
Korpusov et al. Blow-up of solutions of a full non-linear equation of ion-sound waves in a plasma with non-coercive non-linearities
Gross et al. Spectral numerical exterior calculus methods for differential equations on radial manifolds
Han et al. Study on a BFC-Based POD-Galerkin reduced-order model for the unsteady-state variable-property heat transfer problem
Zang et al. Isogeometric boundary element for analyzing steady-state heat conduction problems under spatially varying conductivity and internal heat source
Vlasyuk et al. Numerical solution of three-dimensional problems of filtration consolidation with regard for the influence of technogenic factors by the method of radial basis functions.
Tamaki et al. Efficient dimension-by-dimension higher order finite-volume methods for a Cartesian grid with cell-based refinement
Zhu et al. Fast electrothermal coupling calculation method for supporting digital twin construction of electrical equipment
Li et al. A combination of isogeometric technique and scaled boundary method for the solution of the steady-state heat transfer problems in arbitrary plane domain with Robin boundary
CN105808508B (zh) 一种求解不确定热传导问题的随机正交展开方法
Carey et al. On penalty methods for interelement constraints
Li et al. A meshless model for transient heat conduction analyses of 3D axisymmetric functionally graded solids

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant