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