热动力学问题的数值计算

第17卷第1期高压物理学报Vo.17,No.12003年3月CHINESE JOURNAL OF HIGH PRESSURE PHYSICS文章编号:10005773(2003)01-0035-10热动力学问题的数值计算林华令,黄风雷(北京理工大学,北京100081)摘要:对含热传导的流体动力学方程组,用有限元方法进行数值求解。采用傅里叶热传导计算热流、用热流连续条件计算单元间接触面的温度、用三角形传热法计算体单元表面的热流,考虑各向同性弹塑性流体材料模型、三项式物态方程和导热系数与状态的相关性,给出了傅里叶热传导、接触传热、热应力应变效应、以及混合物冲击压缩特性等算例。对混合物冲击温度的数值模拟表明,小颗粒混合物在冲击压缩过程中,颗粒间的温度有差别、稍有波动,并随时间趋向于一致,以至热平衡。关词:热动力学;有限元法;数值模拟;热效应;混合物中图分类号:O521;O241.82文献标识码:A1引言由受热或受冲击引起的热动力过程涉及物质的运动变形、热传导和状态变化等数值模拟中通常根据问题所侧重的方面将物质输运和热输运分离进行数值分析,因而研制的计算程序可能需要两套,套用于数值模拟动力学过程另一套用于数值模拟热传导过程,两套程序通过外部数据文件交换初始现场和计算结果,其复杂性和近似性显而易见,当然也限定了其适用性和通用性。也有含热传导的动力学计算程序,但计算方法中往往在某些问题方面进行了简化,通常采用近似的能量与温度关系、不考虑温度对热力学状态的影响3,当然对数值模拟的能力也就有了一定的限制。对于物质运动变形与传热时间尺度相差悬殊的问题,采用物质输运和热输运分离的方法是可取的;但在一些特定的问题中热传导可能对物质的运动变形、以至状态变化有较大的影响,分离的数值模拟方法对结果的准确性影响就不可小视,得到的结果甚至可能不真实。在航天飞行器受气流冲击、炸药起爆、强激光对物质热冲击5、混合物的冲击压缩、冲击压缩后的热传导和热辐射1等研究方面,可能涉及高温高压、而又需要考虑物质间传热的问题。研制一个含热传导适用于高温高压、可处理流体弹塑性及失效、处理模型能力较强的计算程序是很有必要的。2计算方法2.1坐标描述在固定的直角坐标系中,质点的初始(t=0)构形为x(j=1,2,3),在t时刻的现时构形为x(i1,2,3),表示成(X;,t)(1)初始条件x;(X1,0)=X(2)收稿日期:200204-18;修回日期:200206-30中国煤化工作者简介:林华令(1963-),男,硕士,高级工程师CNMHG高压物理学报第17卷x(X,0)=v(X)其中v(X)为初速,上标·号表示对拉格朗日时间t的导数,x,表示为坐标i方向的质点速度。2.2控制方程运动方程为9式中:o、P、f和x分别为柯西( Cauchy)应力当前密度、单位质量体积力和质点加速度。(4)式中出现的重复下标按“约定求和法则”对其指标从1到3求和,下同。求解的运动方程还要满足边界应力条件位移条件和接触间断条件质量守恒方程为pu Po其中:v为相对初始体积的比容,A为常态(常温、常压)密度。材料的初始密度记为m0,则初始相对比容v=P/P,v可以不为1。当w>1时,为硫松材料;当v<1时,为密实材料能量守恒采用内能守恒方程形式给出9:1e=一(p+q)v+w甲/3x;+ah(6)式中:e为相对初始体积的比内能,s为偏应力,p为压强,q为人为体粘性Pl(colEm-C1a Es)(7)EM≥0式中:C和C1为无量纲常数(在DYNA程序中缺省值分别为1.5和0.06),l为长度特征量(在二维计算中取l=√A,A为二维单元的面积;在三维中取l=√,V为三维单元的体积),a为当地声速。(6)式中的E,为应变率张量(1/2)(ax;/0x;+x/7x;)内能守恒方程(6)中考虑了热流和热源的影响,其中的平为热流密度向量,h为单位质量介质在单位时间内热源的供热率。按傅里叶( Fourie)传热定律0,热流密度向量q表示为P,=-k(aT/3x;)式中:T为绝对温度为导热系数。我们仅考虑傅里叶热传导,并且假定传热是各向同性的,但考虑k与物质状态的相关性2.3材料模型和物态方程材料模型用以确定应力、应变关系,可采用多种材料模型。我们以各向同性弹塑性流体模型为例考虑了屈服的硬化效应[10,=0o+Epe+(a+a2 p)max(o, p)(10)式中:,和分别为当前屈服强度和初始屈服强度,E为塑性硬化模量,eP为有效塑性应变,a1和a2为压力硬化参数。弹塑性流体即在流体背景上加上弹塑性σ=5-(p+q)6式中:为克罗内克符号( Kronecker Delta:当t=j时,=1;否则b=0)。流体压强p由当前比容v、能量守恒方程(6)以及物态方程计算得到弹性区的本构关系用虎克( Hooke)定律的增量形式54=2G(E4-EM6/3)+5424+ss(12)式中:G为剪切模量,为旋率张量2=(1/2)(ax;/7x(13)塑性区采用屈服条件判别中国煤化工/2-CNMHG(14)林华令等:热动力学问题的数值计算由(12)式得到的s如果满足∫≤0,说明材料处于弹性变形状态;如果出现∫>0,则说明材料超出弹性状态,产生了塑性变形,此时要对应力作修正,使应力松弛回到屈服面,以确定实际应力和弹塑性状态具体如下:由(12)式得到的称为偏应力试探值5,引入修正因子ma,/√(2/3)356(15)实际偏应力5通过修正因子m与偏应力试探值相乘得到5y=m5(16)由实际偏应力s得到的实际偏应力率5。=ds;/dt,可得弹性应变率偏量sc=(l/2G)(s6;-s.2-smk)(17)而塑性应变率可由偏应变率ε,与其弹性分应变率偏量E;相减得到/3有效塑性应变e为e=regie)vade(20)可以采用各种形式的物态方程,本文中以含温度的三项式物态方程为例61p=p(p)+R了3xe++1Aar(21)e-Ae()+32++1mor其中R为摩尔气体常数H为摩尔质量,A为常态密度A为常态电子热容系数,p(p)为冷压,e(p)为单位质量的冷能,()为 Gruneisen参数,o=p/A,r=a[T/Tmn(p)]为晶格的热贡献,Tm(p)为熔化温度,a和为常数(可近似取a=0.0409,7=1.606)。根据冷压p冷能e、 Gruneisen参数γ和熔化温度Tm等与密度p的关系,可用三项式物态方程(21)描述压强和温度变化较宽的固液状态,材料物态参数的确定可参考有关资料。冷能、冷压和 Gruneisen参数采用28/Le(1-2)-84/32 d(p u2)/du其中g为材料参数。2.4数值离散方法对控制方程进行数值离散采用了DYNA所给出的有限元方法,其中运动方程与质量守恒方程的离散求解方法无需变动不再重述。对内能守恒方程则要补充热流和热源两项计算将(6)式对时间作中心差分,可得到内能守恒的差分方程e"t=e"+△x"l/2-「+p)+q141(1-矿)+△Q+Ahx12△r(23)式中:上标n表示!时刻的量,n+1表示增加一个时间步长△后在t“+1=p+△时的量,n+1/2为时间半点r”时的量对应的各时刻量定义在同一单元或网格区域内。△z和△Q分别为单位初始体积的畸变能增量和传热增量△z1=(v”1+)(s1+(24)中国煤化工(g,)"S1CNMHG高物理学报第17卷式中:V。为单元初始体积,S为单元第i个表面积,9)"为单元沿第i个表面外法向的热流密度。计算传热增量时,在(22)式中采用了上一时间整点P的量,传热增量△Q采用了面积分代替差分的方法,而表面外法向的热流密度可表示为(9)"=-k"(T-T)/d式中:T和k”分别为单元温度和导热系数,d为单元中心到第i个表面的距离,T为单元第i个表面的温度表面温度应根据单元与其相邻单元或边界特性确定。对有相邻单元的表面温度,假定热量不在相邻单元的表面间累积则过表面的热流连续,表面温度T为T=(T+a”T")/(1+a")(27)√kd/(k"a(28)式中:T和k为与第i个面相邻的单元温度和导热系数,d为该单元中心到这个表面的距离对给定边界条件的情况,根据给定的条件,计算表面温度T或热流密度(9)1。对绝热边界,取(qn)"=0;对处于恒温T。环境中的单元取T=T;对给定温度边界条件或热流密度边界条件,根据当前时刻的边界温度或热流密度,确定表面温度T或热流密度();。对于接触面上的单元表面温度,应根据相邻单元温度、相邻表面积等,按表面处热流连续条件计算表面温度T7,热流连续条件为T-ts=2kd(29)式中:J为与第i个面相邻的单元数,T和k为与第i个面相邻的第j个单元的温度和导热系数,d和S"为该第j个单元中心到这个表面的距离和相邻的面积。从(29)式可得Tn-(T+)/(…+)(30)Nk",d对于三维体单元,采用三角形传热法。如果单元表面是四边形,可能出现四点不在一个平面上的情况,为此将四边形表面划分成四个三角形,这四个三角形对应的一个边分别取为表面四边形的一个边它们的对角点取为四边形的中心,中心坐标按四边形的四个角点取平均得到。按(27)式分别计算各个三角形表面的温度和热流密度,累加这四个三角形的传热增量即为通过这四边形表面的传热增量。根据运动方程与质量守恒方程的解,结合各向同性弹塑性流体模型可以更新应力、应变,具体方法可参阅有关资料2。根据内能守恒的差分方程(23),结合含温度的三项式物态方程,可以采用迭代求解法,计算状态量e+、p1、T+1等以更新热力学状态量实际导热系数是与状态有关的。对于金属导热系数与德拜( Debye)温度的关系为式中:6b为 Debye特征温度,C为常数。而6b与 Gruneisen参数存在关系y=-dInep/(dnv)(33)结合(32)式和(33)式,可得y(udu(34)式中k为常态导热系数。在实际计算中采用 Griineiseny=2/3+(y其中y。为常态 Gruneisen参数。由此,得到导热系数与THE中国煤化工(35)CNMHG第1期林华令等:热动力学问题的数值计算k= ko vexp[-2(%-2/3)(-1)](36)对于非金属也有导热系数与状态量的关系3),因热传导进程比冲击波运动慢,所以时间步长仍可按DYNA中的方法确定;对不考虑物质运动的纯热传导问题,为节省计算时间,可适当增大时间步3算例基于以上计算方法,分别研制了二维有限元热动力学计算程序和三维有限元热动力学计算程序程序具有对热动力学耦合问题进行数值模拟的能力,可用于分析热应力效应冲击压缩过程中的热传导等问题;程序也具有对纯动力学和纯热传导问题进行数值模拟的能力;也可以修改、扩充材料模型和物态方程,以使其具有更强的数值模拟能力。3.1数值模拟傅里叶热传导问题仅考虑材料的傅里叶热传导,不考虑物质的运动和形变以检验数值模拟热传导的准确性。在数值计算中,可以由跳过运动方程与质量守恒方程的求解实现。不考虑应力应变时,内能守恒方程(6)转化成傅里叶热传导方程。以二维(x,y)为例,傅里叶热传导方程可表示成aT0式中:为定容比热容,k为导热系数。对金属按(32)~(36)式,不考虑形变时,k近似为常数作为算例,用边长a=1cm的正方形高温铝作二维热传导计算铝的初始温度T=2000K,边界取恒温T1=293K环境条件。如不考虑晶格和电子热运动对物态方程的贡献,此时c=(3R/21)P和为常数,从(37)式可见,该式为线性热传导方程。对高温铝的二维线性热传导问题由分离变量法可得解析解r(x,y,D)=T1+16(7)sin(2m+1)(2)x)n(2n+1)(2)yD(38)式中D=k/c为导温系数。在傅里叶热传导数值模拟计算中,将a和阳取为零,初始密度A取为常态密度A,即p=阳、k=kn,其它材料参数见表1。铝中网格按20×20划分,相当于每个网格宽度为0.05cm图1给出了高温铝经历了20ms热传导时的数值模拟温度分布。从图1可见,等温线具有较好的对称性,由于环境温度低于铝的温度故铝边缘的等温线上的温度低于内部温度。0图2给出了y=0.5cm对称线上的温度分布0从图2可见,三种网格划分10×10、20×20和30×30所得的温度分布,与解析解(38)式所得的结果具有较好的一致性;但相对来说,较粗的网格划分10×10所得的结果比细分网格与解析解比较要差些,图1高温铝热传导等温线20的网格划分数值模拟结果已与解析解接近两端附近,即在边缘处一个单元的范围内数值模拟的温度与解析解有差别,这种差别来源于数值模拟EYnf temperature in the aluminium中国煤化工:293KCNMHG.=20高压物理学报第17卷结果后处理方法的不足因数值模拟保存的温度是单元中心的温度,温度曲线是根据单元中心温度内插或外推得到的,属于正常现象;网格分得越细,数值模拟越接近解析解,细分网格会增大计算机时,所以在数值模拟中要把握网格的划分。如要考虑铝的晶格和电子热运动对物态方程的贡献,则c,与温度相关,此时,(37)式成为非线性热传导方程。要给出非线性热传导问题的解析解很困难。对髙温铝的非线性热传导进行了数值模拟,得到了有关结果,铝的材料参数见表1,铝中网格按20×20划分数值模拟结果表明:是否考虑晶格和电子热运动对物态方程的贡献对热传导结果的影响很小。高温铝中心处在20ms时的温度相对偏差为0.05%这表明在2000K以内,铝的晶格和电子热对c的影响较小,因而在一定范围内,线性热传导方程是非线性热传导的一种很好的近似实际材料在常压、高温时的密度低于常态密度。按照物态方程(21)式和(22)式,可得2000K的铝在常压(P=0)下的初始相对比容v=1.11323,对应的初始密度p=2.502g/cm3。数值模拟结果表明:修正初始密度后,是否考虑晶格和电子热对物态方程的贡献,对热传导的结果影响很小;但是否考虑修正初始密度,对热传导的结果有一定的影响初始密度对热传导的影响可从图3中看出。图3中的v=1曲线为未修正初始密度的热传导数值模拟结果,Po=0曲线为按常压修正初始密度的热传导数值模拟结果。高温铝中心处在20ms时的温度相差56K,修正初始密度的温度要低些,这符合修正初始密度后的定容比热容cn=(3R/2)p变小导温系数D=k/c变大、使得热传导更快的道理综上对傅里叶热传导的计算和比较说明我们研制的程序可用于纯热传导计算,也具有计算线性和非线性热传导问题的能力。ND-Consider effects uf thermal0.204图2沿对称线的温度分布图3不同数值模拟方式的比较Fig 2 Distributions of temperature along the symmetrical Fig 3 Comparison among the results from the differentline y=0. 5 cm in the aluminium(t=20 msAmerical simulation modes( Along y=0. 5 cm.= 20 ms3.2热效应的数值模拟计算和分析热的应力应变效应在工程技术方面有重要的应用,解析求解往往变得无能为力,数值模拟方法是可实现的途径。采用弹塑性流体模型的有限元热动力学数值模拟计算方法可以用来计算热的应力应变效应同样以边长a=1cm的正方形20K高温铝处于293K环境中的二维热传导模型为例,进行热效应的数值模拟。铝中网格按20×20划分,铝的有关材料参数见表1。衰1材料参数Table1 Parameters of materialsMaterials中国煤化工/(g/cm2)O)/CW/(cmK)]Aluminium 2.78327.642.0097.91962.83530.0CNMHG 2.37林华令等:热动力学问题的数值计算y=0.5cm对称线上20ms时的温度分布见图3的ND曲线。将考虑热效应的温度分布曲线(ND线)与未考虑热效应的温度分布曲线(p0=0线)进行比较可见:由于有热效应使得高温铝在热传导的过程中因冷却而收缩,热效应的温度曲线低于不考虑热效应的温度曲线图4给出了高温铝经历20ms热传导时的等位移线。从图4可见:整体位移方向向中心;四个顶点是位移最大的地方,最大位移为0.0197cm。图5和图6分别给出了高温铝经历20ms热传导时的等温线和x方向主应力等值线由数值模拟后处理得到的有关物理、力学量等值线和随时间、空间的分布情况,可用以分析受热和冷却过程物体内部的热效应。下面给出两物体热接触时的热效应。热接触物体分别为圆柱形的铜和铝,其中0.4cm×0.2cm钢的初始温度为1500K,0.4cm×0.1cm铝的初始温度为293K,铜的下底面与铝的上底面接触,接Min(-)=3.035Max(+)}=970×10tA=5640x10F2128-82Max(+)=1.865×1图4高温铝热效应等位移线图5高温铝热效应等温线Fig4 Contours of displacement in the aluminiumFig 5 Contours of temperature in the aluminiumCThe intial aluminum: 2000 K, the enviomment: 293 K, t=20 ms) (The initial aluminium: 200 K, the enviomment: 293 K, t =20 m)触面间热流连续、并可自由滑移,铝的下底面为固定Min(-=-0.3395不动、绝热边界,其它材料界面取自由边界、恒温Max(+)=0.1864293K环境条件。按轴对称进行二维(r,z)热动力学有限元数值模拟每个网格宽度为0.01cm,材料参数见表1图7给出了铜和铝接触热传导6ms时的温度分布。从图?可见:热从高温的铜通过与铝接触的界面传入铝中,铜的温度下降,铝的温度升高,6ms时最高温度已位于原来低温的铝内,从初始温度、接触条件和环境条件等分析,计算结果是合理的,界面间连续的等温线也表明了对接触热传导问题数值模拟有较好的准确性;材料的热效应比较明显,初始温度较高的铜由于降温,边缘收缩,而初始处于常温的铝,由于升温,则表现了膨胀;材料的收缩与影胀,Fi图6x方向主应力等值线使得接触界面滑移但热流仍存在,也表明了对接触中国煤化工 the aluminium问题的热传导、热效应数值模拟是合理的。CNMHG高压物理学报第178给出了这两个柱体侧面上、下4个顶点的径向位移随时间的变化曲线,其中A、B为铜的上、下顶点C、D为铝的上、下顶点。从图8可见各质点位移随时间的变化不一致,并且有起伏,说明材料边界膨胀和收缩不同步膨胀和收缩量也有差别。铝膨胀到一定程度要停止膨胀,而转为收缩,这符合铝初期因吸热而膨胀、后期因放热而收缩的规律;铜的上端初期因放热,到一定程度转为膨胀,这应该是铜的轴向也要收缩的二维效应(见图7)引起的;铜的下端在6ms时间内一直收缩但收缩速率随时间变以预计下端收缩到一定程度也会停止。2.5xl0-2.0x10-2x101-4x105x1012.5×M0Mms)图7高温铜与铝接触热传导温度分布图8质点位移随时间的变化Fig 7 Contours of temperature inFig.8 Curves of node displacement vs timethe copper and al(A and B: upper and bottom top nodes on the copper(The initial copper: 1500 K, aluminium: 293 Kthe enviornment: 293 K,t=6 ms)on the aluminium side face33混合物冲击压缩过程的数值模拟混合物冲击压缩特性的有限元数值模拟属于三维问题。将混合物组元材料颗粒按原定的配比随机分布在预先划定的三维网格中,用三维动力学有限元程序模拟其冲击压缩过程,以研究混合物的冲击压缩特性,可得到混合物的形变过程组元间相互作用过程、趋于平衡过程和冲击压缩数据等,同时可以验证“叠加原理”进一步采用混合物冲击压缩和热传导分离的数值模拟方法研究了混合物的冲击温度,分析了冲击压缩后的热传导对颗粒间温度分布的影响印,采用热动力学计算程序数值模拟混合物的冲击压缩特性,其优点是在模拟冲击压缩过程中考虑热传导的影响对于颗粒度较小的混合物,如合金,尤为合适。采用铜和铝两种材料组成的混合物,以方形(尺寸为0.006cm×0.02cm×0.002cm)三维网格单元按60×20×20形式划分,各单元随机分布铜和铝两种材料,体积配比为1:1,对应的混合物颗粒度为1μm。冲击压编采用速度加载M()+M()方法(加载速度方向后半部质点赋速度w,前半部质点赋速度0),引入质点运动周期性边界条件,以消除边缘稀疏的影响。对周期性边界,考虑了界面传热的周期性,即界面单元的外界面温度按单元温度和对应周期性边界上单元的温度,采用热流连续条件计算。材料模型中不考虑弹塑性效应,采用三项式物态方程,物态参数见表1。图9给出了w=4km/s时计算的铜/铝混合物靶内沿加载速度方向同一层上不同单元温度随时的变化曲线。从图9可见,在冲击压缩时间内各个单元温中国煤化工变化稍有波动,到一定时间趋于一致,以至达到热平衡,经历初始压缩到热传导,计算结果体现不出单元间趋于热平衡的过程。CNMHG考虑单元间第1期林华令等:热动力学问题的数值计算图10给出了靶内相邻9层上各9个单元平均温度随时间的变化曲线。从图10可见,各层平均温度基本一致,说明对于小颗粒混合物,在冲击压缩过程中颗粒间的热传导使混合物各层之间达到比较致的温度、趋于热平衡,其平衡温度可以作为混合物的冲击温度。25002.1x101.5x115009x10P7×10图9温度随时间的变化图10各层温度随时间的变化Fig 9 Temperature vs. time( Following impact ofFig 10 Temperature of layer vs. timeCu-Al mixture;w=4 km/s; temperature of(Following impact of Cu-Al mixture;nine elements on one layer)w=4 km/s)发展了热动力学问题的数值模拟方法和二维、三维有限元计算程序,可处理滑移接触传热、体单元接触传热、各种热边界条件等问题有多种导热系数与状态的相关性、多种材料模型和物态方程可供选择利用。程序具有模拟傅里叶热传导、热应力应变效应、冲击压缩过程中的热传导等问题的能力适用围宽。对混合物热动力学数值模拟表明,小颗粒混合物在冲击压缩过程中,能在较短时间内达到热平References:[1] Wilkins M L, Blum R E, Cronshagen E, et al. A Method for Computer Simulation of Problems in Solid Mechanicsnd Gas Dynamics in Three Dimensions and Time [R]. UCRI-51574, 1974.1-7[2J Hallquist JO. ISDYNA Theoretical Manual [R]. Livermore Software Technology Corporation, Copyright 19911998.29.1-29.6.[3] Huang C G, Duan Z P. Numerical Investigation on the Dynamic Responses of Structures Including the Effects ofThermal Conductivity [JJ. Explosion and Shock Waves, 2001, 21(4): 253-258. (in Chinese)黄展光段祝平含热传导的冲击动力学有限元程序的研究与应用[爆炸与冲击,2001,21(4):253-258.[4] Zhang G R,Chen D N. Shock Initiation of Condensed Explosive [M]. Beijing: National Defense Industry Press1991.100-101.( in Chinese)章冠人,陈大年凝聚炸药起爆动力学[M]北京国防工业出版社,1991.100-10L5] Jiang S E, Yang J M, Ding Y N, et al. Investigation on Shack and Heat Waves in Aluminum Irradiated by laser[]. Chinese Journal of High Pressure Physics, 2000, 14(4)中国煤化工江少恩扬家敏,丁耀南,等激光坚动铝靶产生冲击波和热波的实CNMHGO,144)64-268[6] Lin H L, Huang F L, Yu W R. Numerical Simulation of Shocs remperature ot MIxtures During Shock Loading44理学报第17卷[J]. Chinese Journal of High Pressure Physics, 2002, 16(1):46-56(in Chinese)林华令,黄风雷,于万瑞混合物冲击温度的数值模拟[J].高压物理学报,2002,16(1):46-56.[7] Zhang G R Problem of Heat Conduction Behind a Shock Wave Front [J]. Chinese Journal of High PressurePhysics, 1994.8(1): 9-13. (in Chines章冠人,冲击波阵面后温度的热传导问题[J]高压物理学报,1994,8(1):9-13[8] Hong Y J, Tan H,Jin X. The Effect of Initial Density on Radiation Characteristics of the Wave Front [].ChineseJournal of High Pressure Physics, 2000, 14(2): 133-138.(in Chinese)洪延姬,谭华金星.介质初始密度对波阵面辐射特性的影响[]高压物理学报,200,14(2):13-138J Hallquist J O. Theoretical Manual for DYNA3D [R]. UCID-19401, 1983. 5-44.[10] Du X. The Elements of Continuum Mechanics [M]. Beijing: Tsinghua University Press. 1985. 87-89:119 (in Chinese)杜珣.连续介质力学引论[M].北京:清华大学出版社,1985[11] Yun SR, Tu H J, Liang D S, et al. Calculation Method of Detonation Mechanics [M]. Beijing: Beijing Institute ofTechnology Press, 1995. 238-239. (in Chinese)恽寿榕涂侯杰梁德寿,等.爆炸力学计算方法[M].北京:北京理工大学出版社,1995.238-23912] Zhang C B Li M S, Zhang S Z. A Three Pahse Equation of State for 2024Al [J]. Chinese Journal of High PressurePhysics, 1989, 3(4):279-283. (in Chines张春斌李茂生,张世泽2024A的三相物态方程[]高压物理学报,1989,3(4)1279-283.[13] Tang W H. The Pressure and Temperature Dependence of Thermal Conductivity for Nonmetal Crystals [JChinese Journal of High Pressure Physics, 1994, 8(2): 125-132. (in Chinese汤文辉非金属晶体导热系数与压力和温度相关性[冂高压物理学报,1994,8(2):125-132.[14] Gu C H, Li D Q, Chen SX, et al. Equations of Mathematical Physics [M]. Beijing: Higher Education Press, 19876-79.(in Chinese)谷超豪李大潜陈恕行,等数学物理方程[M]北京:高等教育出版社,1984.76-79[15] L in H L Simulation of Shock Compression Behavior of Mixture by Using the Finite Element Method [J].ChineseJournal of High Pressure Physics, 1998, 12(1):40-46. ( in Chinese)林华令.有限元法模拟混合物的冲击压缩特性[高压物理学报,1998,12(1):40-46Numerical Simulation of Thermodynamics ProblemLIN Hua-Ling, HUANG Feng-I(Beijing Institute of Technology, Beijing 100081, China)Abstract: The equations of hydrodynamics including conduction of heat are solved numerically byusing the finite element method. The heat flow is calculated according to the Fourier's heat conduetion. The temperature on the contact side of elements is calculated by the continuity condition of heatflow. According to triangular partitions divided from a quadrangle side of body element, convectionheat flow is calculated. Considering the isotropic-elastic-plastic-hydrodynamic material model, thequation of state with three contributions, and the correlation between thermal conductivity and temperature,some computational examples, including the Fourier's heat conduction, contact heat transferdeformation effects of heat, and shock compression of mixture, are given by the finite element code.The numerical simulation of shock temperature of mixture indicates that the temperatures amongute grains fluctuate, and tend to reach a thermal equiKey words: thermodynamics; finite element method, numerYHa中国煤化工t: mixtureCNMHG