熱動(dòng)力學(xué)問(wèn)題的數值計算

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