第17卷第1期高壓物理學報Vo.17,No.12003年3月CHINESE JOURNAL OF HIGH PRESSURE PHYSICS文章編號:10005773(2003)01-0035-10熱動力學問題的數(shù)值計算林華令,黃風雷(北京理工大學,北京100081)摘要:對含熱傳導的流體動力學方程組,用有限元方法進行數(shù)值求解。采用傅里葉熱傳導計算熱流、用熱流連續(xù)條件計算單元間接觸面的溫度、用三角形傳熱法計算體單元表面的熱流,考慮各向同性彈塑性流體材料模型、三項式物態(tài)方程和導熱系數(shù)與狀態(tài)的相關性,給出了傅里葉熱傳導、接觸傳熱、熱應力應變效應、以及混合物沖擊壓縮特性等算例。對混合物沖擊溫度的數(shù)值模擬表明,小顆?；旌衔镌跊_擊壓縮過程中,顆粒間的溫度有差別、稍有波動,并隨時間趨向于一致,以至熱平衡。關詞:熱動力學;有限元法;數(shù)值模擬;熱效應;混合物中圖分類號:O521;O241.82文獻標識碼:A1引言由受熱或受沖擊引起的熱動力過程涉及物質(zhì)的運動變形、熱傳導和狀態(tài)變化等數(shù)值模擬中通常根據(jù)問題所側(cè)重的方面將物質(zhì)輸運和熱輸運分離進行數(shù)值分析,因而研制的計算程序可能需要兩套,套用于數(shù)值模擬動力學過程另一套用于數(shù)值模擬熱傳導過程,兩套程序通過外部數(shù)據(jù)文件交換初始現(xiàn)場和計算結(jié)果,其復雜性和近似性顯而易見,當然也限定了其適用性和通用性。也有含熱傳導的動力學計算程序,但計算方法中往往在某些問題方面進行了簡化,通常采用近似的能量與溫度關系、不考慮溫度對熱力學狀態(tài)的影響3,當然對數(shù)值模擬的能力也就有了一定的限制。對于物質(zhì)運動變形與傳熱時間尺度相差懸殊的問題,采用物質(zhì)輸運和熱輸運分離的方法是可取的;但在一些特定的問題中熱傳導可能對物質(zhì)的運動變形、以至狀態(tài)變化有較大的影響,分離的數(shù)值模擬方法對結(jié)果的準確性影響就不可小視,得到的結(jié)果甚至可能不真實。在航天飛行器受氣流沖擊、炸藥起爆、強激光對物質(zhì)熱沖擊5、混合物的沖擊壓縮、沖擊壓縮后的熱傳導和熱輻射1等研究方面,可能涉及高溫高壓、而又需要考慮物質(zhì)間傳熱的問題。研制一個含熱傳導適用于高溫高壓、可處理流體彈塑性及失效、處理模型能力較強的計算程序是很有必要的。2計算方法2.1坐標描述在固定的直角坐標系中,質(zhì)點的初始(t=0)構形為x(j=1,2,3),在t時刻的現(xiàn)時構形為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的導數(shù),x,表示為坐標i方向的質(zhì)點速度。2.2控制方程運動方程為9式中:o、P、f和x分別為柯西( Cauchy)應力當前密度、單位質(zhì)量體積力和質(zhì)點加速度。(4)式中出現(xiàn)的重復下標按“約定求和法則”對其指標從1到3求和,下同。求解的運動方程還要滿足邊界應力條件位移條件和接觸間斷條件質(zhì)量守恒方程為pu Po其中:v為相對初始體積的比容,A為常態(tài)(常溫、常壓)密度。材料的初始密度記為m0,則初始相對比容v=P/P,v可以不為1。當w>1時,為硫松材料;當v<1時,為密實材料能量守恒采用內(nèi)能守恒方程形式給出9:1e=一(p+q)v+w甲/3x;+ah(6)式中:e為相對初始體積的比內(nèi)能,s為偏應力,p為壓強,q為人為體粘性Pl(colEm-C1a Es)(7)EM≥0式中:C和C1為無量綱常數(shù)(在DYNA程序中缺省值分別為1.5和0.06),l為長度特征量(在二維計算中取l=√A,A為二維單元的面積;在三維中取l=√,V為三維單元的體積),a為當?shù)芈曀佟?6)式中的E,為應變率張量(1/2)(ax;/0x;+x/7x;)內(nèi)能守恒方程(6)中考慮了熱流和熱源的影響,其中的平為熱流密度向量,h為單位質(zhì)量介質(zhì)在單位時間內(nèi)熱源的供熱率。按傅里葉( Fourie)傳熱定律0,熱流密度向量q表示為P,=-k(aT/3x;)式中:T為絕對溫度為導熱系數(shù)。我們僅考慮傅里葉熱傳導,并且假定傳熱是各向同性的,但考慮k與物質(zhì)狀態(tài)的相關性2.3材料模型和物態(tài)方程材料模型用以確定應力、應變關系,可采用多種材料模型。我們以各向同性彈塑性流體模型為例考慮了屈服的硬化效應[10,=0o+Epe+(a+a2 p)max(o, p)(10)式中:,和分別為當前屈服強度和初始屈服強度,E為塑性硬化模量,eP為有效塑性應變,a1和a2為壓力硬化參數(shù)。彈塑性流體即在流體背景上加上彈塑性σ=5-(p+q)6式中:為克羅內(nèi)克符號( Kronecker Delta:當t=j時,=1;否則b=0)。流體壓強p由當前比容v、能量守恒方程(6)以及物態(tài)方程計算得到彈性區(qū)的本構關系用虎克( Hooke)定律的增量形式54=2G(E4-EM6/3)+5424+ss(12)式中:G為剪切模量,為旋率張量2=(1/2)(ax;/7x(13)塑性區(qū)采用屈服條件判別中國煤化工/2-CNMHG(14)林華令等:熱動力學問題的數(shù)值計算由(12)式得到的s如果滿足∫≤0,說明材料處于彈性變形狀態(tài);如果出現(xiàn)∫>0,則說明材料超出彈性狀態(tài),產(chǎn)生了塑性變形,此時要對應力作修正,使應力松弛回到屈服面,以確定實際應力和彈塑性狀態(tài)具體如下:由(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)可以采用各種形式的物態(tài)方程,本文中以含溫度的三項式物態(tài)方程為例61p=p(p)+R了3xe++1Aar(21)e-Ae()+32++1mor其中R為摩爾氣體常數(shù)H為摩爾質(zhì)量,A為常態(tài)密度A為常態(tài)電子熱容系數(shù),p(p)為冷壓,e(p)為單位質(zhì)量的冷能,()為 Gruneisen參數(shù),o=p/A,r=a[T/Tmn(p)]為晶格的熱貢獻,Tm(p)為熔化溫度,a和為常數(shù)(可近似取a=0.0409,7=1.606)。根據(jù)冷壓p冷能e、 Gruneisen參數(shù)γ和熔化溫度Tm等與密度p的關系,可用三項式物態(tài)方程(21)描述壓強和溫度變化較寬的固液狀態(tài),材料物態(tài)參數(shù)的確定可參考有關資料。冷能、冷壓和 Gruneisen參數(shù)采用28/Le(1-2)-84/32 d(p u2)/du其中g為材料參數(shù)。2.4數(shù)值離散方法對控制方程進行數(shù)值離散采用了DYNA所給出的有限元方法,其中運動方程與質(zhì)量守恒方程的離散求解方法無需變動不再重述。對內(nèi)能守恒方程則要補充熱流和熱源兩項計算將(6)式對時間作中心差分,可得到內(nèi)能守恒的差分方程e"t=e"+△x"l/2-「+p)+q141(1-礦)+△Q+Ahx12△r(23)式中:上標n表示!時刻的量,n+1表示增加一個時間步長△后在t“+1=p+△時的量,n+1/2為時間半點r”時的量對應的各時刻量定義在同一單元或網(wǎng)格區(qū)域內(nèi)?！鱶和△Q分別為單位初始體積的畸變能增量和傳熱增量△z1=(v”1+)(s1+(24)中國煤化工(g,)"S1CNMHG高物理學報第17卷式中:V。為單元初始體積,S為單元第i個表面積,9)"為單元沿第i個表面外法向的熱流密度。計算傳熱增量時,在(22)式中采用了上一時間整點P的量,傳熱增量△Q采用了面積分代替差分的方法,而表面外法向的熱流密度可表示為(9)"=-k"(T-T)/d式中:T和k”分別為單元溫度和導熱系數(shù),d為單元中心到第i個表面的距離,T為單元第i個表面的溫度表面溫度應根據(jù)單元與其相鄰單元或邊界特性確定。對有相鄰單元的表面溫度,假定熱量不在相鄰單元的表面間累積則過表面的熱流連續(xù),表面溫度T為T=(T+a”T")/(1+a")(27)√kd/(k"a(28)式中:T和k為與第i個面相鄰的單元溫度和導熱系數(shù),d為該單元中心到這個表面的距離對給定邊界條件的情況,根據(jù)給定的條件,計算表面溫度T或熱流密度(9)1。對絕熱邊界,取(qn)"=0;對處于恒溫T。環(huán)境中的單元取T=T;對給定溫度邊界條件或熱流密度邊界條件,根據(jù)當前時刻的邊界溫度或熱流密度,確定表面溫度T或熱流密度();。對于接觸面上的單元表面溫度,應根據(jù)相鄰單元溫度、相鄰表面積等,按表面處熱流連續(xù)條件計算表面溫度T7,熱流連續(xù)條件為T-ts=2kd(29)式中:J為與第i個面相鄰的單元數(shù),T和k為與第i個面相鄰的第j個單元的溫度和導熱系數(shù),d和S"為該第j個單元中心到這個表面的距離和相鄰的面積。從(29)式可得Tn-(T+)/(…+)(30)Nk",d對于三維體單元,采用三角形傳熱法。如果單元表面是四邊形,可能出現(xiàn)四點不在一個平面上的情況,為此將四邊形表面劃分成四個三角形,這四個三角形對應的一個邊分別取為表面四邊形的一個邊它們的對角點取為四邊形的中心,中心坐標按四邊形的四個角點取平均得到。按(27)式分別計算各個三角形表面的溫度和熱流密度,累加這四個三角形的傳熱增量即為通過這四邊形表面的傳熱增量。根據(jù)運動方程與質(zhì)量守恒方程的解,結(jié)合各向同性彈塑性流體模型可以更新應力、應變,具體方法可參閱有關資料2。根據(jù)內(nèi)能守恒的差分方程(23),結(jié)合含溫度的三項式物態(tài)方程,可以采用迭代求解法,計算狀態(tài)量e+、p1、T+1等以更新熱力學狀態(tài)量實際導熱系數(shù)是與狀態(tài)有關的。對于金屬導熱系數(shù)與德拜( Debye)溫度的關系為式中:6b為 Debye特征溫度,C為常數(shù)。而6b與 Gruneisen參數(shù)存在關系y=-dInep/(dnv)(33)結(jié)合(32)式和(33)式,可得y(udu(34)式中k為常態(tài)導熱系數(shù)。在實際計算中采用 Griineiseny=2/3+(y其中y。為常態(tài) Gruneisen參數(shù)。由此,得到導熱系數(shù)與THE中國煤化工(35)CNMHG第1期林華令等:熱動力學問題的數(shù)值計算k= ko vexp[-2(%-2/3)(-1)](36)對于非金屬也有導熱系數(shù)與狀態(tài)量的關系3),因熱傳導進程比沖擊波運動慢,所以時間步長仍可按DYNA中的方法確定;對不考慮物質(zhì)運動的純熱傳導問題,為節(jié)省計算時間,可適當增大時間步3算例基于以上計算方法,分別研制了二維有限元熱動力學計算程序和三維有限元熱動力學計算程序程序具有對熱動力學耦合問題進行數(shù)值模擬的能力,可用于分析熱應力效應沖擊壓縮過程中的熱傳導等問題;程序也具有對純動力學和純熱傳導問題進行數(shù)值模擬的能力;也可以修改、擴充材料模型和物態(tài)方程,以使其具有更強的數(shù)值模擬能力。3.1數(shù)值模擬傅里葉熱傳導問題僅考慮材料的傅里葉熱傳導,不考慮物質(zhì)的運動和形變以檢驗數(shù)值模擬熱傳導的準確性。在數(shù)值計算中,可以由跳過運動方程與質(zhì)量守恒方程的求解實現(xiàn)。不考慮應力應變時,內(nèi)能守恒方程(6)轉(zhuǎn)化成傅里葉熱傳導方程。以二維(x,y)為例,傅里葉熱傳導方程可表示成aT0式中:為定容比熱容,k為導熱系數(shù)。對金屬按(32)~(36)式,不考慮形變時,k近似為常數(shù)作為算例,用邊長a=1cm的正方形高溫鋁作二維熱傳導計算鋁的初始溫度T=2000K,邊界取恒溫T1=293K環(huán)境條件。如不考慮晶格和電子熱運動對物態(tài)方程的貢獻,此時c=(3R/21)P和為常數(shù),從(37)式可見,該式為線性熱傳導方程。對高溫鋁的二維線性熱傳導問題由分離變量法可得解析解r(x,y,D)=T1+16(7)sin(2m+1)(2)x)n(2n+1)(2)yD(38)式中D=k/c為導溫系數(shù)。在傅里葉熱傳導數(shù)值模擬計算中,將a和陽取為零,初始密度A取為常態(tài)密度A,即p=陽、k=kn,其它材料參數(shù)見表1。鋁中網(wǎng)格按20×20劃分,相當于每個網(wǎng)格寬度為0.05cm圖1給出了高溫鋁經(jīng)歷了20ms熱傳導時的數(shù)值模擬溫度分布。從圖1可見,等溫線具有較好的對稱性,由于環(huán)境溫度低于鋁的溫度故鋁邊緣的等溫線上的溫度低于內(nèi)部溫度。0圖2給出了y=0.5cm對稱線上的溫度分布0從圖2可見,三種網(wǎng)格劃分10×10、20×20和30×30所得的溫度分布,與解析解(38)式所得的結(jié)果具有較好的一致性;但相對來說,較粗的網(wǎng)格劃分10×10所得的結(jié)果比細分網(wǎng)格與解析解比較要差些,圖1高溫鋁熱傳導等溫線20的網(wǎng)格劃分數(shù)值模擬結(jié)果已與解析解接近兩端附近,即在邊緣處一個單元的范圍內(nèi)數(shù)值模擬的溫度與解析解有差別,這種差別來源于數(shù)值模擬EYnf temperature in the aluminium中國煤化工:293KCNMHG.=20高壓物理學報第17卷結(jié)果后處理方法的不足因數(shù)值模擬保存的溫度是單元中心的溫度,溫度曲線是根據(jù)單元中心溫度內(nèi)插或外推得到的,屬于正?，F(xiàn)象;網(wǎng)格分得越細,數(shù)值模擬越接近解析解,細分網(wǎng)格會增大計算機時,所以在數(shù)值模擬中要把握網(wǎng)格的劃分。如要考慮鋁的晶格和電子熱運動對物態(tài)方程的貢獻,則c,與溫度相關,此時,(37)式成為非線性熱傳導方程。要給出非線性熱傳導問題的解析解很困難。對髙溫鋁的非線性熱傳導進行了數(shù)值模擬,得到了有關結(jié)果,鋁的材料參數(shù)見表1,鋁中網(wǎng)格按20×20劃分數(shù)值模擬結(jié)果表明:是否考慮晶格和電子熱運動對物態(tài)方程的貢獻對熱傳導結(jié)果的影響很小。高溫鋁中心處在20ms時的溫度相對偏差為0.05%這表明在2000K以內(nèi),鋁的晶格和電子熱對c的影響較小,因而在一定范圍內(nèi),線性熱傳導方程是非線性熱傳導的一種很好的近似實際材料在常壓、高溫時的密度低于常態(tài)密度。按照物態(tài)方程(21)式和(22)式,可得2000K的鋁在常壓(P=0)下的初始相對比容v=1.11323,對應的初始密度p=2.502g/cm3。數(shù)值模擬結(jié)果表明:修正初始密度后,是否考慮晶格和電子熱對物態(tài)方程的貢獻,對熱傳導的結(jié)果影響很小;但是否考慮修正初始密度,對熱傳導的結(jié)果有一定的影響初始密度對熱傳導的影響可從圖3中看出。圖3中的v=1曲線為未修正初始密度的熱傳導數(shù)值模擬結(jié)果,Po=0曲線為按常壓修正初始密度的熱傳導數(shù)值模擬結(jié)果。高溫鋁中心處在20ms時的溫度相差56K,修正初始密度的溫度要低些,這符合修正初始密度后的定容比熱容cn=(3R/2)p變小導溫系數(shù)D=k/c變大、使得熱傳導更快的道理綜上對傅里葉熱傳導的計算和比較說明我們研制的程序可用于純熱傳導計算,也具有計算線性和非線性熱傳導問題的能力。ND-Consider effects uf thermal0.204圖2沿對稱線的溫度分布圖3不同數(shù)值模擬方式的比較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熱效應的數(shù)值模擬計算和分析熱的應力應變效應在工程技術方面有重要的應用,解析求解往往變得無能為力,數(shù)值模擬方法是可實現(xiàn)的途徑。采用彈塑性流體模型的有限元熱動力學數(shù)值模擬計算方法可以用來計算熱的應力應變效應同樣以邊長a=1cm的正方形20K高溫鋁處于293K環(huán)境中的二維熱傳導模型為例,進行熱效應的數(shù)值模擬。鋁中網(wǎng)格按20×20劃分,鋁的有關材料參數(shù)見表1。衰1材料參數(shù)Table1 Parameters of materialsMaterials中國煤化工/(g/cm2)O)/CW/(cmK)]Aluminium 2.78327.642.0097.91962.83530.0CNMHG 2.37林華令等:熱動力學問題的數(shù)值計算y=0.5cm對稱線上20ms時的溫度分布見圖3的ND曲線。將考慮熱效應的溫度分布曲線(ND線)與未考慮熱效應的溫度分布曲線(p0=0線)進行比較可見:由于有熱效應使得高溫鋁在熱傳導的過程中因冷卻而收縮,熱效應的溫度曲線低于不考慮熱效應的溫度曲線圖4給出了高溫鋁經(jīng)歷20ms熱傳導時的等位移線。從圖4可見:整體位移方向向中心;四個頂點是位移最大的地方,最大位移為0.0197cm。圖5和圖6分別給出了高溫鋁經(jīng)歷20ms熱傳導時的等溫線和x方向主應力等值線由數(shù)值模擬后處理得到的有關物理、力學量等值線和隨時間、空間的分布情況,可用以分析受熱和冷卻過程物體內(nèi)部的熱效應。下面給出兩物體熱接觸時的熱效應。熱接觸物體分別為圓柱形的銅和鋁,其中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)觸面間熱流連續(xù)、并可自由滑移,鋁的下底面為固定Min(-=-0.3395不動、絕熱邊界,其它材料界面取自由邊界、恒溫Max(+)=0.1864293K環(huán)境條件。按軸對稱進行二維(r,z)熱動力學有限元數(shù)值模擬每個網(wǎng)格寬度為0.01cm,材料參數(shù)見表1圖7給出了銅和鋁接觸熱傳導6ms時的溫度分布。從圖?可見:熱從高溫的銅通過與鋁接觸的界面?zhèn)魅脘X中,銅的溫度下降,鋁的溫度升高,6ms時最高溫度已位于原來低溫的鋁內(nèi),從初始溫度、接觸條件和環(huán)境條件等分析,計算結(jié)果是合理的,界面間連續(xù)的等溫線也表明了對接觸熱傳導問題數(shù)值模擬有較好的準確性;材料的熱效應比較明顯,初始溫度較高的銅由于降溫,邊緣收縮,而初始處于常溫的鋁,由于升溫,則表現(xiàn)了膨脹;材料的收縮與影脹,Fi圖6x方向主應力等值線使得接觸界面滑移但熱流仍存在,也表明了對接觸中國煤化工 the aluminium問題的熱傳導、熱效應數(shù)值模擬是合理的。CNMHG高壓物理學報第178給出了這兩個柱體側(cè)面上、下4個頂點的徑向位移隨時間的變化曲線,其中A、B為銅的上、下頂點C、D為鋁的上、下頂點。從圖8可見各質(zhì)點位移隨時間的變化不一致,并且有起伏,說明材料邊界膨脹和收縮不同步膨脹和收縮量也有差別。鋁膨脹到一定程度要停止膨脹,而轉(zhuǎn)為收縮,這符合鋁初期因吸熱而膨脹、后期因放熱而收縮的規(guī)律;銅的上端初期因放熱,到一定程度轉(zhuǎn)為膨脹,這應該是銅的軸向也要收縮的二維效應(見圖7)引起的;銅的下端在6ms時間內(nèi)一直收縮但收縮速率隨時間變以預計下端收縮到一定程度也會停止。2.5xl0-2.0x10-2x101-4x105x1012.5×M0Mms)圖7高溫銅與鋁接觸熱傳導溫度分布圖8質(zhì)點位移隨時間的變化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混合物沖擊壓縮過程的數(shù)值模擬混合物沖擊壓縮特性的有限元數(shù)值模擬屬于三維問題。將混合物組元材料顆粒按原定的配比隨機分布在預先劃定的三維網(wǎng)格中,用三維動力學有限元程序模擬其沖擊壓縮過程,以研究混合物的沖擊壓縮特性,可得到混合物的形變過程組元間相互作用過程、趨于平衡過程和沖擊壓縮數(shù)據(jù)等,同時可以驗證“疊加原理”進一步采用混合物沖擊壓縮和熱傳導分離的數(shù)值模擬方法研究了混合物的沖擊溫度,分析了沖擊壓縮后的熱傳導對顆粒間溫度分布的影響印,采用熱動力學計算程序數(shù)值模擬混合物的沖擊壓縮特性,其優(yōu)點是在模擬沖擊壓縮過程中考慮熱傳導的影響對于顆粒度較小的混合物,如合金,尤為合適。采用銅和鋁兩種材料組成的混合物,以方形(尺寸為0.006cm×0.02cm×0.002cm)三維網(wǎng)格單元按60×20×20形式劃分,各單元隨機分布銅和鋁兩種材料,體積配比為1:1,對應的混合物顆粒度為1μm。沖擊壓編采用速度加載M()+M()方法(加載速度方向后半部質(zhì)點賦速度w,前半部質(zhì)點賦速度0),引入質(zhì)點運動周期性邊界條件,以消除邊緣稀疏的影響。對周期性邊界,考慮了界面?zhèn)鳠岬闹芷谛?即界面單元的外界面溫度按單元溫度和對應周期性邊界上單元的溫度,采用熱流連續(xù)條件計算。材料模型中不考慮彈塑性效應,采用三項式物態(tài)方程,物態(tài)參數(shù)見表1。圖9給出了w=4km/s時計算的銅/鋁混合物靶內(nèi)沿加載速度方向同一層上不同單元溫度隨時的變化曲線。從圖9可見,在沖擊壓縮時間內(nèi)各個單元溫中國煤化工變化稍有波動,到一定時間趨于一致,以至達到熱平衡,經(jīng)歷初始壓縮到熱傳導,計算結(jié)果體現(xiàn)不出單元間趨于熱平衡的過程。CNMHG考慮單元間第1期林華令等:熱動力學問題的數(shù)值計算圖10給出了靶內(nèi)相鄰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)發(fā)展了熱動力學問題的數(shù)值模擬方法和二維、三維有限元計算程序,可處理滑移接觸傳熱、體單元接觸傳熱、各種熱邊界條件等問題有多種導熱系數(shù)與狀態(tài)的相關性、多種材料模型和物態(tài)方程可供選擇利用。程序具有模擬傅里葉熱傳導、熱應力應變效應、沖擊壓縮過程中的熱傳導等問題的能力適用圍寬。對混合物熱動力學數(shù)值模擬表明,小顆?；旌衔镌跊_擊壓縮過程中,能在較短時間內(nèi)達到熱平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]北京國防工業(yè)出版社,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)中國煤化工江少恩揚家敏,丁耀南,等激光堅動鋁靶產(chǎn)生沖擊波和熱波的實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)林華令,黃風雷,于萬瑞混合物沖擊溫度的數(shù)值模擬[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)洪延姬,譚華金星.介質(zhì)初始密度對波陣面輻射特性的影響[]高壓物理學報,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)杜珣.連續(xù)介質(zhì)力學引論[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的三相物態(tài)方程[]高壓物理學報,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湯文輝非金屬晶體導熱系數(shù)與壓力和溫度相關性[冂高壓物理學報,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)谷超豪李大潛陳恕行,等數(shù)學物理方程[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