滾塑成型工藝加熱階段的數(shù)值研究

劉學(xué)軍

（北京工商大學(xué)材料與機(jī)械工程學(xué)院，北京100048）

0 前言

滾塑工藝的成型周期主要受外部環(huán)境與模具的對(duì)流換熱、內(nèi)部空氣和塑料粉末的混合物與模具的對(duì)流換熱、塑料層的熔融與固化相變等多方面的影響。為此，對(duì)滾塑成型工藝建立一個(gè)準(zhǔn)確的流動(dòng)與傳熱的理論模型，然后通過(guò)數(shù)值計(jì)算對(duì)其進(jìn)行模擬就顯得尤為重要。這種數(shù)值模擬可以研究多個(gè)滾塑參數(shù)對(duì)其成型周期的影響，從而為滾塑成型工藝的最優(yōu)化提供理論依據(jù)。

郭同凱［1］和高鐿等［2］分別為滾塑成型工藝建立了傳熱模型。但都沒(méi)有給出任何理論分析解或數(shù)值計(jì)算結(jié)果。高桂天等［3］為卡車翼子板的滾塑成型工藝建立了一個(gè)簡(jiǎn)單的一維熱傳導(dǎo)及熔融相變的傳熱模型。他們根據(jù)這個(gè)模型給出了加熱時(shí)間的理論計(jì)算結(jié)果，并將其與實(shí)驗(yàn)結(jié)果相比較。何躍龍等［4］對(duì)一個(gè)圓柱形塑料制品的滾塑工藝建立了一維的傳熱模型。根據(jù)該模型給出了模具加熱時(shí)間、塑料粉末熔融燒結(jié)時(shí)間以及溫度場(chǎng)的理論分析解。

Gogos等［5］對(duì)一個(gè)球形中空制品的滾塑成型工藝的加熱階段建立了一維瞬態(tài)的傳熱模型，利用有限差分法對(duì)其進(jìn)行了真正的數(shù)值計(jì)算。在該模型里，他們假設(shè)塑料粉末是均勻混合的。在加熱階段，塑料粉末只有一個(gè)作為集總參數(shù)的溫度在隨時(shí)間變化。計(jì)算區(qū)域只包含模具本身，不包括塑料粉末和內(nèi)部空氣，但在模具內(nèi)壁面處建立了模具與塑料粉末之間的熱平衡方程。筆者認(rèn)為該傳熱模型有以下幾個(gè)方面需要改進(jìn)：（1）塑料粉末的導(dǎo)熱率比較小，因此忽略其內(nèi)部熱阻對(duì)其采用集總參數(shù)法會(huì)產(chǎn)生較大的誤差。（2）由于計(jì)算區(qū)域不包括內(nèi)部空氣，因此該模型無(wú)法計(jì)算出加熱階段內(nèi)部空氣的溫度。（3）該模型需要估計(jì)模具內(nèi)部的對(duì)流換熱系數(shù)，這帶有較大的不確定性。

本文將對(duì)球形中空塑料制品的滾塑成型工藝的加熱階段建立一個(gè)更接近于實(shí)際的傳熱模型。在這個(gè)模型里，計(jì)算區(qū)域?qū)U(kuò)大到包括模具、黏附在模具內(nèi)壁面的塑料粉末層以及內(nèi)部空氣，這樣內(nèi)部空氣的溫度隨時(shí)間的變化規(guī)律也就可以計(jì)算出來(lái)。由于這是耦合換熱，因此模型中不再出現(xiàn)內(nèi)部的熱邊界條件，內(nèi)部的對(duì)流換熱系數(shù)也就不再需要估計(jì)。本文將利用有限容積法通過(guò)Fluent軟件對(duì)文獻(xiàn)［6］中的三種實(shí)驗(yàn)條件進(jìn)行數(shù)值模擬，并將計(jì)算結(jié)果與該文獻(xiàn)中的實(shí)驗(yàn)結(jié)果進(jìn)行比較，以驗(yàn)證本模型的準(zhǔn)確性。最后再應(yīng)用此模型分別對(duì)不同外部加熱溫度、外部對(duì)流換熱系數(shù)以及塑料粉末層厚度的情形進(jìn)行數(shù)值計(jì)算，以研究這些參數(shù)對(duì)滾塑成型工藝加熱時(shí)間的影響，并得出了一些有益的結(jié)論。

1 流動(dòng)與傳熱模型

本模型假設(shè)在整個(gè)滾塑成型過(guò)程中，塑料粉末總是緊密地貼附在模具內(nèi)壁面上形成一層均勻的塑料導(dǎo)熱層。即不考慮塑料粉末以及隨后熔融塑料的流動(dòng)，但考慮內(nèi)部空氣的流動(dòng)。這個(gè)假設(shè)與實(shí)際情況之間具有一定的偏差，本文后面將此模型的計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果相比較可以證明該假設(shè)所帶來(lái)的誤差在工程上仍是可以接受的。如圖1所示，本模型的計(jì)算區(qū)域包括模具、塑料粉末層以及內(nèi)部空氣，將模具的導(dǎo)熱、塑料粉末層的導(dǎo)熱和隨后的熔融以及內(nèi)部空氣的自然對(duì)流換熱耦合起來(lái)求解?？紤]到此球形計(jì)算區(qū)域沿經(jīng)度方向的對(duì)稱性，因此這是一個(gè)二維瞬態(tài)的傳熱學(xué)問(wèn)題，即所有的計(jì)算物理量都是時(shí)間坐標(biāo)t和空間坐標(biāo)（r，θ）的函數(shù)。

圖1 滾塑成型工藝加熱階段的流動(dòng)與傳熱模型Fig.1 Heat transfer and fluid flow model of heating phase of rotational molding

對(duì)于模具，其控制方程是二維瞬態(tài)的導(dǎo)熱方程。外壁面（ro）的邊界條件為第三類熱邊界條件，內(nèi)壁面（rm）的邊界條件為熱流密度連續(xù)條件。對(duì)于塑料層，其控制方程也是二維瞬態(tài)的導(dǎo)熱方程。內(nèi)、外壁面（rp，rm）的邊界條件均為熱流密度連續(xù)條件。塑料層的熔融相變將采用Enthalpy-porosity方法來(lái)計(jì)算。這種方法不直接跟蹤熔融界面的移動(dòng)，而是根據(jù)熵平衡在每一次迭代中計(jì)算液態(tài)相的百分比（介于0和1之間）。詳細(xì)內(nèi)容可參見(jiàn)文獻(xiàn)［7］。

對(duì)于內(nèi)部空氣，其流動(dòng)和傳熱模式為自然對(duì)流換熱。控制方程為二維瞬態(tài)層流的連續(xù)性方程、動(dòng)量方程和能量方程。其壁面（rp）上采用無(wú)滑移速度邊界條件和熱流密度連續(xù)條件。在以上所有方程中除內(nèi)部空氣的密度之外，均采用常物性假設(shè)。考慮到在整個(gè)加熱過(guò)程中，內(nèi)部空氣的溫度變化比較大，為保證計(jì)算結(jié)果的準(zhǔn)確性，本模型不采用自然對(duì)流的Boussinesq假設(shè)，而是把內(nèi)部空氣當(dāng)成是密度隨溫度變化的理想氣體。

2 計(jì)算參數(shù)的確定

筆者將采用本文的理論模型去數(shù)值模擬文獻(xiàn)［6］中的3種實(shí)驗(yàn)條件，如表1所示。

表1 計(jì)算的3個(gè)實(shí)驗(yàn)條件Tab.1 Three experimental conditions calculated

在對(duì)以上3個(gè)實(shí)驗(yàn)條件進(jìn)行數(shù)值模擬時(shí)，本文將采用與文獻(xiàn)［6］完全相同的物性參數(shù)、結(jié)構(gòu)參數(shù)以及工藝參數(shù)（如表2所示）以保證計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果的可比性。

3 Fluent軟件的設(shè)置

在Fluent應(yīng)用程序中選擇瞬態(tài)的二維平面層流模型，并選擇melting／solidification模塊。在設(shè)置環(huán)境條件時(shí)，需要將豎直方向上的加速度設(shè)為-9.8m／s2以考慮內(nèi)部空氣自然對(duì)流的影響。在設(shè)置空氣物性時(shí)，要選擇理想氣體。計(jì)算采用SIMPLE算法［8］。為保證瞬態(tài)問(wèn)題數(shù)值解的穩(wěn)定性，在沿時(shí)間方向上的差分格式采用一階全隱格式。關(guān)于松弛因子，壓力項(xiàng)采用0.3，密度項(xiàng)、體積力項(xiàng)和能量項(xiàng)都采用1，動(dòng)量項(xiàng)采用0.5，液態(tài)相百分比采用0.8。時(shí)間步長(zhǎng)設(shè)為5s，每一個(gè)時(shí)刻的迭代次數(shù)設(shè)為90次。計(jì)算結(jié)果表明90次的迭代可保證每一個(gè)時(shí)刻的速度殘差不超過(guò)10-3，溫度殘差不超過(guò)10-9，且熔融開(kāi)始后，塑料層中液態(tài)相百分比的相對(duì)變化不會(huì)超過(guò)千分之一。

表2 數(shù)值模擬所采用的參數(shù)Tab.2 Parameters for numerical simulation

在迭代計(jì)算過(guò)程中，當(dāng)模具內(nèi)壁面溫度升高到塑料粉末的熔融溫度時(shí)，液態(tài)相的百分比開(kāi)始大于0，這表明塑料粉末開(kāi)始熔融。隨后液態(tài)相的百分比隨時(shí)間越來(lái)越大。當(dāng)其值達(dá)到1時(shí)，說(shuō)明全部塑料粉末都已經(jīng)熔融完畢，加熱過(guò)程結(jié)束，此時(shí)即可終止迭代。

4 數(shù)值計(jì)算結(jié)果及分析

4.1 計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果的比較

首先應(yīng)用本文的理論模型對(duì)表1所列的三種實(shí)驗(yàn)條件進(jìn)行數(shù)值模擬。表3中給出了本文模型計(jì)算所得的加熱時(shí)間與文獻(xiàn)［6］中實(shí)驗(yàn)所得的加熱時(shí)間相比較的結(jié)果。加熱時(shí)間是指從裝有塑料粉末的模具在加熱室里開(kāi)始受熱直到所有塑料粉末都已熔融沉積在模具內(nèi)壁面為止所花費(fèi)的時(shí)間。

從表3中可以看出對(duì)于加熱時(shí)間來(lái)說(shuō)，本模型的計(jì)算結(jié)果與文獻(xiàn)［6］中的實(shí)驗(yàn)結(jié)果吻合得較好。在3種實(shí)驗(yàn)條件下，最大的相對(duì)誤差不超過(guò)9%，在文獻(xiàn)［3］中，理論加熱時(shí)間與實(shí)際加熱時(shí)間的相對(duì)誤差為13.8%，仍認(rèn)為在工程允許的范圍內(nèi)，本文計(jì)算的相對(duì)誤差更小，從而證明了本文模型的準(zhǔn)確性。從該表還可以看出計(jì)算所得的加熱時(shí)間比實(shí)測(cè)的加熱時(shí)間都要短，本文認(rèn)為造成該誤差的主要有兩個(gè)原因：（1）本文假設(shè)自始至終塑料粉末都均勻緊密地黏貼在模具的內(nèi)壁面上，這樣塑料粉末與模具內(nèi)壁面之間始終是以純導(dǎo)熱方式來(lái)進(jìn)行傳熱。但實(shí)際的情況是塑料粉末和內(nèi)部空氣的混合物與模具內(nèi)壁面進(jìn)行對(duì)流換熱。因?yàn)樵搶?duì)流換熱是在一個(gè)封閉空腔內(nèi)進(jìn)行的，因此其換熱系數(shù)很小，不會(huì)超過(guò)5W／（m2·K）［5-6］。也就是說(shuō)，本文的模型在一定程度上過(guò)分地強(qiáng)化了模具內(nèi)壁面向塑料粉末的傳熱，因此使得計(jì)算所得的加熱過(guò)程提前結(jié)束。（2）文獻(xiàn)［6］通過(guò)一套測(cè)量系統(tǒng)來(lái)連續(xù)地監(jiān)測(cè)模具內(nèi)部的空氣溫度。當(dāng)該空氣溫度開(kāi)始急劇增大，即認(rèn)為所有的塑料粉末已經(jīng)完全熔融沉積，此時(shí)所測(cè)得的時(shí)間即為整個(gè)加熱過(guò)程所花費(fèi)的時(shí)間。本文認(rèn)為用這樣的方法來(lái)判斷塑料粉末完成熔融在時(shí)間上可能有一定的滯后，也就是說(shuō)文獻(xiàn)［6］所測(cè)得的加熱時(shí)間可能比實(shí)際值要稍微長(zhǎng)一些。

表3 加熱時(shí)間的計(jì)算值與實(shí)驗(yàn)值的比較Tab.3 Comparison of calculated results of heating time with experimental results

圖2 模具外壁面溫度的計(jì)算值與實(shí)測(cè)值的比較Fig.2 Comparison of calculated values with experimental values of mold temperatures

圖3 內(nèi)部空氣溫度的計(jì)算值與實(shí)測(cè)值的比較Fig.3 Comparison of calculated values with experimetal values of internal air temperatures

圖2和圖3分別顯示了對(duì)于實(shí)驗(yàn)條件1，本模型計(jì)算所得的模具外壁面溫度和內(nèi)部空氣溫度與文獻(xiàn)［6］中相應(yīng)的實(shí)測(cè)值相比較的結(jié)果?？梢钥闯鲞@兩個(gè)溫度的數(shù)值計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果吻合得也比較好，這從另一個(gè)側(cè)面進(jìn)一步證實(shí)了本文理論模型的準(zhǔn)確性。從圖2可以看出，模具外壁面溫度的計(jì)算值在整個(gè)加熱過(guò)程都要低于實(shí)測(cè)值。造成這個(gè)誤差的主要原因可能還是如前所述。在本文模型中，模具的內(nèi)壁面與塑料粉末層總保持緊密接觸，即模具與塑料粉末層之間的熱阻比實(shí)際情況偏小。因此在加熱過(guò)程中就有更多的熱量傳給了塑料粉末層和內(nèi)部空氣，而留給模具本身的熱量相對(duì)較少，因此計(jì)算出來(lái)的模具溫度就比實(shí)際值偏低一些。而獲得較多熱量的塑料粉末層則能更快地完成熔融，從而使得理論上的加熱過(guò)程比實(shí)際提前結(jié)束。另外同樣也獲得較多熱量的內(nèi)部空氣理論上溫升就會(huì)較快，這也解釋了為什么在圖3中內(nèi)部空氣溫度的計(jì)算值隨時(shí)間增大的速率要高于實(shí)測(cè)值。但在加熱過(guò)程的最后大約300s，該計(jì)算值隨時(shí)間增大的很慢，無(wú)限地趨近于塑料粉末的熔融溫度128℃，這一點(diǎn)也與實(shí)測(cè)的內(nèi)部空氣溫度隨時(shí)間的變化規(guī)律有一些差別。對(duì)于實(shí)驗(yàn)條件2和3，從表3可以看出其加熱時(shí)間的計(jì)算值與實(shí)測(cè)值之間的相對(duì)誤差比實(shí)驗(yàn)條件1更小，因此可以推斷其模具外壁面溫度和內(nèi)部空氣溫度的計(jì)算值與實(shí)測(cè)值之間應(yīng)該都吻合得比實(shí)驗(yàn)條件1更好，對(duì)這兩個(gè)實(shí)驗(yàn)條件就沒(méi)必要再做上述的比較。

4.2 外部加熱溫度和對(duì)流換熱系數(shù)對(duì)加熱時(shí)間的影響

首先采用與表2中完全相同的計(jì)算參數(shù)，只是將加熱溫度從250℃增大到380℃，把數(shù)值模擬的結(jié)果畫(huà)在圖4中以表示外部加熱溫度對(duì)滾塑工藝加熱時(shí)間的影響。從圖4可以看出加熱時(shí)間隨加熱溫度的增大而減少，但在加熱溫度超過(guò)300℃后，加熱時(shí)間隨加熱溫度的增大而減少的速度變得越來(lái)越慢。也就是說(shuō)，盡管提高外部加熱溫度確實(shí)可以縮短滾塑工藝的加熱時(shí)間，但在加熱溫度超過(guò)300℃后，該措施的有效性就變得越來(lái)越差。實(shí)際上如果加熱溫度太高，塑料粉末就容易發(fā)生氧化，致使?jié)L塑制品的力學(xué)性能變差。因此從加熱的有效性以及力學(xué)性能的強(qiáng)化這兩方面來(lái)綜合考慮，合適的外部加熱溫度應(yīng)該是300℃左右。

圖4 加熱時(shí)間隨外部加熱溫度的變化Fig.4 Variation of heating time with outer heating temperature

加熱時(shí)間（t）由兩部分組成，一部分是將模具從室溫加熱到塑料粉末熔融溫度所花費(fèi)的時(shí)間（t1），另一部分是塑料粉末從剛開(kāi)始熔融到全部完成熔融所花費(fèi)的時(shí)間（t2）。很顯然t＝t1＋t2。上述計(jì)算結(jié)果表明對(duì)于加熱 溫 度 250 ℃ 的 情 形，t1／t ＝ 48.3%，t2／t＝51.7%；對(duì)于加熱溫度380 ℃的情形，t1／t＝41.2%，t2／t＝58.8%。通過(guò)兩種情形的對(duì)比可以發(fā)現(xiàn)隨著加熱溫度的提高，塑料粉末的熔融時(shí)間在整個(gè)加熱時(shí)間中所占的百分比明顯增加。這說(shuō)明造成高溫加熱有效性變差的主要因素是塑料粉末的熔融過(guò)程拖累了整個(gè)加熱過(guò)程的結(jié)束。這可能是由于在較高的加熱溫度下，更多的熱量被模具吸收走，而留給塑料粉末用于熔融的熱量卻增加的不多，致使整個(gè)熔融過(guò)程不能很快結(jié)束?？偠灾荒軉渭兊赝ㄟ^(guò)提高外部加熱溫度來(lái)縮短滾塑工藝的加熱時(shí)間，而要采取一些其他的輔助措施。比如可以在塑料粉末開(kāi)始熔融的時(shí)候，采用模具內(nèi)部加熱的方法來(lái)加快塑料粉末的熔融過(guò)程。

接著本文采用與表2相同的計(jì)算參數(shù)，將外部加熱溫度定為300℃，將外部對(duì)流換熱系數(shù)從18W／（m2·K）依次增大到28W／（m2·K），把數(shù)值模擬的結(jié)果畫(huà)在圖5中以表示外部對(duì)流換熱系數(shù)對(duì)滾塑工藝加熱時(shí)間的影響。從該圖可以看出加熱時(shí)間隨外部對(duì)流換熱系數(shù)的增大而減少，但在對(duì)流換熱系數(shù)超過(guò)22 W／（m2·K）之后，加熱時(shí)間隨之減少的速度也開(kāi)始逐步放緩。但與圖4相比，可以看出圖5中曲線的彎曲程度并不嚴(yán)重，這說(shuō)明采用增大外部對(duì)流換熱系數(shù)的方法比采用提高外部加熱溫度的方法能更有效地縮短滾塑工藝的加熱時(shí)間?？梢酝ㄟ^(guò)改變外部加熱方式、改進(jìn)模具的表面結(jié)構(gòu)等措施來(lái)提高外部對(duì)流換熱系數(shù)。

4.3 塑料粉末層厚度對(duì)加熱時(shí)間的影響

圖5 加熱時(shí)間隨外部對(duì)流換熱系數(shù)的變化Fig.5 Variation of heating time with outer convective heat transfer coefficient

采用與表2完全相同的計(jì)算參數(shù)，并將加熱溫度保持為300℃。然后應(yīng)用本文計(jì)算模型分別對(duì)不同的塑料粉末層厚度的情形進(jìn)行數(shù)值模擬，將其結(jié)果畫(huà)在圖6中以表示滾塑成型工藝的加熱時(shí)間隨塑料粉末層厚度的變化規(guī)律。從該圖可以看出隨著塑料粉末層厚度的增大，需要的加熱時(shí)間也越來(lái)越長(zhǎng)。另外還可以看出圖中曲線的斜率越變?cè)酱?，這說(shuō)明加熱時(shí)間并不按線性關(guān)系隨塑料粉末層厚度增加，而是隨該厚度以越來(lái)越快的速度急劇增加。也就是說(shuō)對(duì)于厚度較大的塑料制品，其滾塑成型工藝的時(shí)間成本會(huì)增大很多，因此滾塑成型的塑料制品的厚度應(yīng)設(shè)計(jì)得盡可能的小。但另一方面，為提高塑料制品的強(qiáng)度，又希望其厚度能比較大。所以采用滾塑成型工藝的設(shè)計(jì)者應(yīng)綜合考慮塑料制品的成本及強(qiáng)度這兩方面的要求才能確定合適的塑料制品厚度。

圖6 加熱時(shí)間隨塑料粉末層厚度的變化Fig.6 Variation of heating time with thickness of plastic powder layer

5 結(jié)論

（1）本文為滾塑成型工藝的加熱階段建立了一個(gè)新的理論模型。該模型將模具的導(dǎo)熱、塑料粉末層的導(dǎo)熱和隨后的熔融相變，以及內(nèi)部空氣的自然對(duì)流換熱耦合起來(lái)求解。在該模型中，采用有限容積法并通過(guò)Fluent軟件對(duì)滾塑成型工藝的加熱階段進(jìn)行了數(shù)值模擬。計(jì)算所得的加熱時(shí)間、模具溫度和內(nèi)部空氣溫度與文獻(xiàn)［6］中的實(shí)驗(yàn)結(jié)果都吻合得較好，從而驗(yàn)證了本文理論模型的正確性。

（2）提高外部加熱溫度可以縮短滾塑成型工藝的加熱時(shí)間。但在加熱溫度達(dá)到300℃時(shí)，繼續(xù)增大外部加熱溫度，加熱時(shí)間減少的速度明顯變慢，即通過(guò)提高加熱溫度來(lái)縮短加熱時(shí)間的有效性開(kāi)始變差。

（3）增大外部對(duì)流換熱系數(shù)也可以縮短滾塑成型工藝的加熱時(shí)間。但在對(duì)流換熱系數(shù)超過(guò)22W／（m2·K）之后，加熱時(shí)間隨之減少的速度開(kāi)始逐步放緩。與提高外部加熱溫度相比，增大外部對(duì)流換熱系數(shù)能更有效地縮短加熱時(shí)間。

（4）隨著塑料制品厚度的增大，滾塑成型工藝的加熱時(shí)間會(huì)急劇增加，即采用滾塑工藝的制造成本會(huì)大大增加。

［1］郭同凱.旋轉(zhuǎn)模塑過(guò)程的熱交換計(jì)算機(jī)模擬系統(tǒng)［J］.新技術(shù)新工藝，1997，（3）：7-9.Guo Tongkai.Heat Transfer Computer-simulating System of Rotational Molding Process［J］.New Technology and Process，1997，（3）：7-9.

［2］高 鐿，王福生，高幼銀.滾塑成型工藝中的傳熱模型［C］.2005中國(guó)國(guó)際滾塑論壇論文集，2005：62-66.

［3］高桂天，梁 策，史新星，等.滾塑成型加熱時(shí)間的計(jì)算［J］.吉林大學(xué)學(xué)報(bào)（工學(xué)版），2008，38（2）：340-343.Gao Guitian，Liang Ce，Shi Xinxing，et al.Heating Time Algorithm for Rotational Molding［J］.Journal of Jilin University（Engineering and Technology Edition），2008，38（2）：340-343.

［4］何躍龍，丁玉梅，李維維.單軸滾塑成型傳熱過(guò)程的數(shù)學(xué)建模［J］.塑料科技，2009，37（5）：35-39.He Yuelong，Ding Yumei，Li Weiwei.Mathematical Model of Heat Transfer During Rotational Molding of Single Axis［J］，Plastics Science and Technology，2009，37（5）：35-39.

［5］Gogos G，Olson L G，Liu Xuejun，et al.New Models for Rotational Molding of Plastics［J］.Polymer Engineering and Science，1998，38（9）：1387-1398.

［6］Olson L G，Crawford R，Kearns M，et al.Rotational Molding of Plastics：Comparison of Simulation and Experimental Results for an Axisymmetric Mold［J］.Polymer Engineering and Science，2000，40（8）：1758-1764.

［7］Voller V R，Prakash C.A Fixed-grid Numerical Modeling Methodology for Convection-diffusion Mushy Region Phase-change Problems［J］.Int J Heat Mass Transfer，1987，30（8）：1709-1719.

［8］陶文銓.數(shù)值傳熱學(xué)［M］.西安：西安交通大學(xué)出版社，2001：203-207.