热载荷作用下颗粒增强复合材料温升分布数值模拟

杨万友， 王家序， 黄彦彦， 周青华， 杨 勇

(四川大学 空天科学与工程学院，成都 610065)

卫星、宇宙飞船等航天器常常服役于复杂空间环境中，受到高低温、强辐射、微重力、真空等严酷空间环境的共同作用，材料性能被严重影响.尤其是宽范围的高低温交变，将导致材料发生热变形和产生内部热应力，造成材料性能退化和航天器可靠性降低.此外，针对航天器使用的高性能颗粒增强复合材料，热量在传递的过程中还会受到增强体的影响，使材料表面及内部温度场产生扰动，从而造成表面热变形不协调、内部热应力集中等不利后果，进一步影响航天器的可靠性[1-2].

颗粒增强复合材料作为一种典型的工程材料，其目的是通过往基体材料中添加不同材料属性的增强体以提升基体材料性能.由于热传导属性的差异，增强体的存在将会影响热载荷作用下颗粒增强复合材料的热平衡.对于材料热平衡问题的理论分析方法，国内外学者开展了许多相关研究.早期的理论研究多基于均质材料这一假设，因此可用解析方法直接进行求解.如Carslaw等[3]给出了移动热源作用下半空间表面温升解析解.其后圆形热源[4]、正方形和椭圆形热源[5]作用下半空间表面温升分布的解析公式被其他学者们推导出来.Liu等[1,6]提出了点力或者点热源作用下求解半空间3种状态(暂瞬态、恒瞬态和稳态)下温升、热应力[1]和表面热位移[6]的频率响应函数，这避免了采用Green函数计算温度场时函数存在奇点的问题[7-8].

然而颗粒增强复合材料是由多种性质不同的材料复合而成，呈现出典型的非均质特性.由于其复杂的内部微观结构，难以找到解析解求解相关的热平衡等问题，所以，有限元法(FEM)和半解析法(SAM)等数值方法逐渐被广大学者采用来求解非均质材料问题.如Bazyar等[9]利用比例边界FEM求解了各向异性非均质材料的二维热传导问题；Mortazavi等[10]建立了有限元模型来预测各向异性随机双相复合材料的有效导热系数和弹性模量，并将之与其他模型进行对比分析研究.此外，余天堂等[11]提出了以导热系数为基本参数的热传导扩展FEM来研究非均质材料中的热传导问题；陈康等[12]给出了一种适用于梯度复合材料热传导分析的梯度单元法，结果表明，梯度单元和均匀单元得到的温度场基本一致.类似研究还包括文献[13-15].尽管FEM在各种商业软件支持下广泛应用于颗粒增强复合材料热传导分析，但是该方法具有材料微观结构建模困难、计算耗时长、结果过度依赖网格尺寸和类型等缺点.

为了克服上述FEM的缺点，学者们提出利用SAM来求解与颗粒增强复合材料相关的科学问题.在Eshelby[16]提出等效夹杂方法(EIM)求解非均质材料弹性场这一开创性工作之后，SAM被逐渐发展、丰富和完善[17-22]，可求解含任意形状、分布增强体颗粒增强复合材料的粗糙面弹性接触问题.然而对于颗粒增强复合材料热平衡分析，目前鲜有文献报道.Hatta等[23]类比Eshelby[16]提出的弹性场等效夹杂方法，将应力对应热流密度、应变对应温度梯度、弹性模量对应导热系数，提出了适用于求解含圆/球形、椭圆/球形增强体颗粒增强复合材料的稳态热传导全空间等效夹杂方法.然而，全空间解只能处理远端热载荷问题，真实热载荷作用下颗粒增强复合材料的热平衡问题则需要在半空间条件下求解.

本文将颗粒增强复合材料中增强体设定为椭球，利用镜像法[24]完成热传导等效夹杂方法从全空间到半空间的转换，用于研究热载荷作用下颗粒增强复合材料内部和表面温升分布规律.分析了不同导热系数、深度和位置朝向的增强体对颗粒增强复合材料热平衡行为的影响，并通过对含任意分布球形增强体颗粒增强复合材料温升分布进行求解，验证了算法处理真实颗粒增强复合材料热平衡分析的能力.

1 基本理论和方程

(1)

图1 稳态热传导等效夹杂方法Fig.1 EIM for steady state heat conduction

图2 镜像法基本原理Fig.2 Basic principle of method of images

采用镜像法[24]，实现稳态热传导EIM从全空间到半空间的转换.基于镜像法原理(见图2)，半空间夹杂引起的温度梯度量是两个全空间夹杂引起的温度梯度量和对称面法向热流引起的温度梯度量之和，其求解方法详见文献[1,23].

进一步引入类似弹性场中的Eshelby张量Sij，文献[25]已给出其推导过程及具体形式，此处不再赘述.利用Eshelby张量Sij关联扰动温度梯度与本征温度梯度：

(2)

将上式代入式(1)整理得

(3)

(4)

式中：Φ,i为牛顿势函数Φ的一阶导数，详情参见文献[25]；xi为目标温度坐标点到增强体单元中心坐标的向量.颗粒增强复合材料总温升Tt为均质材料温升T与增强体扰动温升Td之和[26]：

Tt=T+Td

(5)

(6)

式中：Sij为关联第j个等效夹杂相应温度梯度与第i个本征温度梯度之间相互影响的系数；n为增强体的数目.

2 颗粒增强复合材料热平衡算法验证

利用热载荷作用下含单椭球增强体颗粒增强复合材料模型验证本文所提算法的有效性.模型示意图如图1(a)所示，其中椭球形增强体尺寸设置为a=2b=2c=r/2，并令其中心距离表面h=r.外部加载热流密度可表示为

(7)

模型计算区域大小设置为6r×6r×3r，等距离散为128×128×64个网格点，利用上述算法对不同导热系数(K*=4Km或K*=Km/4，其中，基体材料导热系数Km=50.2 W/(m·K)[1])增强体进行计算.同时，利用Abaqus软件中轴对称方式建立同样条件的有限元验证模型，如图3所示.另外值得注意的是，第1节中已直接或间接说明了本文所提模型/算法与相关的有限元分析模型基于下列3点假设：① 颗粒增强复合材料是半无限体；② 热传导分析是在稳态时进行；③ 颗粒增强体之间不重叠.

图4所示为利用上述两种方法计算得到的颗粒增强复合材料下表面温升分布云图，NT11表示节点温度输出变量.当增强体与基体导热系数关系为K*=4Km时，增强体内部温度等值线分布相较基体稀疏，温度变化率较小；而K*=Km/4时，增强体内部温度等值线分布较基体紧凑，温度变化率较大.图5所示为利用两种方法计算得到的模型结果沿z轴的温升分布.从图5中可以看出，增强体单元对基体温度场有一处明显的扰动，并且K*=4Km时会在增强体距离表面近端(z/r=0.5)以及远端(z/r=1.5)附近分别造成相对均质材料的温度降低和升高，而增强体材料属性为K*=Km/4时造成的影响与前者相反且相对较弱.此外，利用FEM计算的结果与本文算法模拟的结果在数值上能够较好地吻合，验证了所提算法的有效性.

图3 含单一椭球增强体二维有限元对称稳态热传导模型Fig.3 2D axisymmetric FEM including an ellipsoidal reinforcement for steady state heat conduction

图4 含单一椭球增强体下表面温升分布Fig.4 Lower surface temperature rise distribution including an ellipsoidal reinforcement

图5 含单一椭球增强体沿z轴的温升分布Fig.5 Temperature rise distribution along the z axis including an ellipsoidal reinforcement

3 计算结果和讨论

采用图1(a)的模型，除了特别说明外，其余参数设置参见第2节，以椭球形增强体导热系数、位置、朝向、数目为研究参数，分析其对热载荷作用下颗粒增强复合材料温升分布影响.

3.1 单个增强体热平衡分析

设置椭球形增强体尺寸参数为a=1.25b=1.25c=r，增强体中心至表面的距离h=0.75r.由于椭球本身各轴长短不一，在空间中有着复杂的位置关系.如图6所示，可采用欧拉角α、β和γ来描述椭球增强体所处的空间位置.目标点从全局坐标到局部坐标的转换可表示为

Rxyz(m)=Rxyz(m)+

a(m,n)[Pxyz(m)-Cxyz(m)]

(8)

m,n=1,2,3

图6 欧拉角和坐标转换关系Fig.6 Relationship of Euler angles and coordinate transformation

式中：m,n分别代表x,y,z；Rxyz(m)和Pxyz(m)分别为目标点的局部坐标和全局坐标；Cxyz(m)为椭球中心的全局坐标;a(m,n)为对应的旋转因子，具体可表示为

(9)

由上述欧拉角的定义，图1(a)模型中椭球增强体各欧拉角分别为α=0°、β=90° 和γ=90°.不论增强体处于什么位置，其对温度场或多或少都有扰动，且这种扰动直接表现为温度的变化.本节利用表面温度变化率η来分析各参数对材料热平衡行为的影响.具体表达式为

(10)

图7 不同导热系数增强体表面温升分布及温度变化率Fig.7 Distribution and change rate of surface temperature rise including an ellipsoidal reinforcement with different thermal conduction properties

图8 不同导热系数增强体下表面温升分布Fig.8 Lower surface temperature rise distribution including an ellipsoidal reinforcement with different thermal conduction properties

3.1.1导热系数 导热系数是影响材料热平衡的关键因素，本节将通过改变增强体与基体导热系数之比λ(λ=K*/Km或λ=Km/K*)使之从1～20之间变化，以探究颗粒增强复合材料内部热平衡规律，结果如图7和8所示.从图7(a)中可以看出，当λ=K*/Km时，随着λ的增大，表面温升分布趋于平缓，表面温度值逐渐减小；当λ=Km/K*时，则得到与前者相反的结论.图7(a)中还给出了K*=20Km、K*=Km/20和K*=Km3种情况下温度标尺相同的表面温升三维等值线图，可直观地看出表面温升情况.图7(b)结果表明，随着λ的增大，η先快速增长后增速逐渐放缓，且λ=K*/Km时相对λ=Km/K*增速更为显著.图8为图7(a)情况对应的下表面沿z轴的温升分布.图中曲线表明增强体距离表面近端和远端位置附近都会对温度场有较大扰动，且λ越大，扰动越明显.产生上述现象的原因在于增强体的导热系数较基体大时，则热量在增强体中集中，导致在增强体距表面近端处温度相对均质材料要低，而远端温度相对均质材料要高.增强体导热系数较基体小时则相反.

3.1.2深度 由于材料受摩擦热等外部热载荷作用时，是从表面往内部传递的，所以，颗粒增强复合材料内部热平衡往往与增强体的深度位置密切相关.设置椭球增强体中心到表面的距离h∈[0.6r, 2r]，研究不同深度增强体对颗粒增强复合材料热平衡的影响规律.图9所示为含不同深度增强体颗粒增强复合材料表面温升分布及温度变化率.图9(a)中表面温升分布曲线表明，不论增强体导热系数大于或小于基体，随着增强体深度的增加，表面温升分布都将越来越接近均质材料.在h/r=0.6时，K*=20Km和K*=Km/20两种情况下基体表面温升分布可从相应的三维等值线图直观看出.从图9(b)中可以看出，随着增强体深度的增加，无论其导热系数大于或小于基体，颗粒增强复合材料表面温度变化率都快速下降，逐渐趋近于0.上述结果表明，随着增强体深度的增加，表面热载荷作用对颗粒增强复合材料内部及表面温度场扰动逐渐减小.

图9 不同深度位置增强体表面温升分布及温度变化率Fig.9 Distribution and change rate of surface temperature rise including an ellipsoidal reinforcement at different depth locations

3.1.3朝向角度 颗粒增强复合材料中类纤维状增强体空间朝向具有随机性，且增强体的朝向角度已被证实对应力场会产生影响[17].本小节将改变其朝向角度，保持α=0° 和β=90° 不变，使欧拉角γ在0°～ 90° 之间变化，探究其对颗粒增强复合材料温度场的影响规律.图10给出了γ分别为0°、45° 和90° 时颗粒增强复合材料下表面温升分布云图，此时K*=4Km.由图可知，增强体不同空间朝向将对材料温升分布产生较大影响.从图11颗粒增强复合材料表面温度变化率曲线中可知，随着γ的逐渐增大，η先减小后增大，且γ=45° 时η最小.

图10 不同朝向角度椭球形增强体下表面温升分布Fig.10 Lower surface temperature rise distribution including an ellipsoidal reinforcement with different facing angles

图11 不同朝向角度椭球形增强体表面温度变化率Fig.11 Distribution and change rate of surface temperature including an ellipsoidal reinforcement with different facing angles

3.2 多增强体热平衡分析

增强体往往随机分布于颗粒增强复合材料中.本小节考虑多个增强体之间的相互作用，进一步分析其对颗粒增强复合材料温度场的影响.

图12(a)为含双球形增强体颗粒增强复合材料模型，增强体中心位于Oxz平面上，沿z轴对称分布，且其中心到表面的距离h=0.6r.球形增强体直径d=r，间距为l.增强体的导热系数为K*=4Km和K*=Km/4，其余参数保持不变.改变间距l，使其在r～3r之间变化，研究颗粒增强复合材料下表面温升分布规律.

图12(b)和12(c)所示为利用本文所提算法求出的含不同间距双球形增强体(K*=4Km)颗粒增强复合材料下表面温升分布云图.从图中可以直观看出不同间距的增强体对基体温度场的影响程度不同，而且当l增大时，会使基体温度场向周边扩散.

图13所示为不同间距双球形增强体颗粒增强复合材料表面温度变化率.图中曲线表明，随着l增大，η逐渐降低，且两种材料增强体造成的η差异也在逐渐减小.

图14(a)所示为中心位置坐标(x,y,z)服从正态分布的200个半径为0.2r球形增强体颗粒增强复合材料模型，即

x～N(0,r2)

y～N(0,r2)

z～N(r,r2)

图12 含双球形增强体模型及下表面温升分布Fig.12 Model including two spherical reinforcements and its lower surface temperature rise distribution

图13 不同间距双球形增强体表面温度变化率Fig.13 Change rate of surface temperature including two spherical reinforcements with different distances

图14 含正态分布球形增强体模型及下表面温升分布Fig.14 Model including spherical reinforcements in normal distribution and its lower surface temperature rise distribution

增强体导热系数为K*=20Km和K*=Km/20，其余参数参照前例，计算结果如图14(b)和14(c)所示.从图中可以看出，随机分布的增强体对材料下表面温升分布扰动较大，且位置接近表层的增强体对表面温度及周围温度影响较为显著，而较深的增强体对温度场的扰动则相对较小.对比K*=20Km和K*=Km/20两种情况，可发现后者温度变化范围更广，且具有较大的温升，说明即使是将增强体分散于基体中，其导热系数对基体温升分布也具有较大影响.

4 结论

(1) 本文提出了一种热载荷作用下颗粒增强复合材料热平衡数值建模方法.对于单椭球增强体模型，增强体导热系数与基体差异越大，越容易造成表面温度变化率增大.因此在颗粒增强复合材料制备过程中，增强体热属性不宜与基体差异太大，否则影响表面温度变化率，从而导致材料表面性能退化.

(2) 单椭球形增强体越靠近表面高温度区时，其对基体表面温度扰动越大.此外，单椭球形增强体朝向会对温度场产生较大影响，在本文设定情况下，欧拉角γ在45° 左右时对表面温度的影响最小.工程实际中，颗粒增强复合材料增强体的分布应尽量远离受热载荷表面区域，且椭球/纤维状增强体欧拉角γ应在45° 左右.

(3) 双球形增强体间距增大，表面温度变化率逐渐降低且由不同导热系数增强体造成的颗粒增强复合材料表面温度变化率差异也在逐渐减小.分布于颗粒增强复合材料中的增强体会对基体温度场造成扰动，且越靠近表面的增强体造成的扰动越显著.

本文所提材料温升分布求解模型可有效考虑材料微观结构对热载荷作用下材料内部热平衡的影响，为颗粒增强复合材料及其他材料散热、吸热等功能需求为导向的内部微观结构优化设计提供理论指导.耦合材料弹、塑、电、磁等不同场的求解模型，可进一步求解不同外载联合作用下非均质材料复杂工程力学问题.