張火明，陳國慶，陸萍藍，余潤梁，管衛(wèi)兵

(1.中國計量大學 計量測試工程學院，浙江 杭州 310018；2.中國計量大學 工程訓練中心，浙江 杭州310018；3.自然資源部第二海洋研究所，浙江 杭州 310012)

圖1 無接觸拉管工藝Figure 1 Contactless tube drawing process

近年來，國內(nèi)外的很多學者對玻璃管的研究主要集中在加工制造方法、理論分析和材料特性參數(shù)等方面。李海蓉等[3]對現(xiàn)有的聚合物擠出過程熔體溫度在線計量技術進行了綜述，為聚合物擠出工業(yè)熔體溫度在線計量技術的開發(fā)和應用提供了參考。苗培培、饒堯、張偉等[4-6]利用光束傳播法對耦合器的拉制過程進行了數(shù)值模擬，得到耦合器輸出光功率隨拉伸長度的變化規(guī)律以及光纖拉錐后的內(nèi)應力變化。陳珊等[7]在2014年設計并制造了一種基于激光加熱的玻璃基毛細管拉伸裝置，對玻璃管進行拉伸試驗，探究實際實驗中玻璃管拉伸規(guī)律。王珀琥和吳必成[8-9]等通過仿真和實驗分析得到了不同溫度、拉伸速度、拉力大小、加熱時間等因素對玻璃管拉制成形的影響。黃建斌等[10]通過有限元方法模擬與仿真，探討了光纖熔融拉錐后其內(nèi)應力大小與光纖被加熱區(qū)域的降溫時間的關系。Quentin等[11]通過數(shù)值模擬研究了單根光纖成型的基本物理過程和靈敏度，分析研究了關鍵參數(shù)和物理機制對內(nèi)應力的影響。Xue等[12]首先建立了共軛多相模型，該模型中包含了光纖加熱爐的傳熱，和使用凈輻射法評估外部輻射熱通量，將蒙特卡羅示蹤法與Polyflow仿真軟件相結(jié)合，以獲得各種爐壁和變形預制件之間的所有影響因素。Florians等[13]設計了一種用于預測冷管拉制過程中管徑大小的模型，利用此模型還能研究相對厚度對管變形的影響，并通過實驗驗證該模型的準確性。本文以毛細管為例，考慮石英材料的粘彈性行為，基于順序耦合熱應力分析法，對管坯在高溫無接觸拉制情況下的成形過程進行數(shù)值模擬分析，研究溫度和拉伸力等工藝參數(shù)對毛細管管成品的幾何參數(shù)的影響規(guī)律。

1 有限元模型建立

1.1 幾何模型的建立及管坯拉制控制方程

本文進行模擬的毛細管管坯為內(nèi)外徑尺寸分別為4 mm和6 mm，高度120 mm的薄管模型，加熱爐的內(nèi)外徑分別為40 mm、95 mm，高度為100 mm，模型如圖2(a)所示。使用ABAQUS軟件自帶的建模版塊建立數(shù)值模擬模型，網(wǎng)格均采用楔形單元形狀，劃分技術為掃掠網(wǎng)格劃分，單元類型為C3D6T，共計11 466個節(jié)點、15 249個網(wǎng)格，網(wǎng)格如圖2(b)所示。加熱爐和管坯之間還存在空氣，在軟件中將加熱爐與管坯之間的模型賦予空氣的特性參數(shù)，材料特性參數(shù)見下表1和表2。

圖2 毛細管模型及網(wǎng)格劃分示意圖Figure 2 The Capillary model and the grid diagram

表1 數(shù)值模擬模型環(huán)境參數(shù)

表2 管坯材料特性參數(shù)

可將拉制過程看作是非穩(wěn)態(tài)或穩(wěn)態(tài)的不可壓縮流動，用Navier-Stokes(N-S)方程來描述，即

(1)

其中，ρ為流體物質(zhì)密度，v為流體速度，t為時間，p為流體各向壓力，F(xiàn)為質(zhì)量力，μ為流體的動態(tài)粘度系數(shù)。

在圓柱坐標系下，設管坯熔體的流動速度為v=uez+wer，其中，u，w分別為沿z軸(軸向)和沿r軸(徑向)方向的速度分量，ez為兩個方向的單位矢量，在管坯連續(xù)穩(wěn)定拉制的過程中，我們可以將重力的影響忽略，即ρF=0，將N-S方程在z、r方向展開。

z方向上的動量守恒方程為

計算機實驗教學中心是地方高校學生計算機實驗教學重要基地，在普及計算機知識、全面培養(yǎng)學生計算機應用能力和創(chuàng)新能力中發(fā)揮越來越重要作用。地方高校在轉(zhuǎn)型背景下，要不斷深化教育教學改革，積極探索計算機實驗教學中心的建設與發(fā)展，利用有限的資金投入，優(yōu)化組合計算機資源，實現(xiàn)計算機資源高效、共享利用，建設具有“規(guī)模化，功能化、網(wǎng)絡化、高效化、安全化、規(guī)范化、節(jié)約化”的計算機實驗教學中心，為未來申請“省級示范中心”打好堅實基礎。

(2)

r方向上的動量守恒方程為

(3)

管坯拉制過程中遵守質(zhì)量守恒，出入加熱爐中的管坯質(zhì)量保持不變，用方程表示如下:

(4)

對于不可壓縮物質(zhì)，其密度不會隨著溫度、時間等參量的變化而改變，所以式(4)可改寫成[11]

(5)

上述方程構(gòu)成了熔融管坯流動的二維控制方程組。

根據(jù)實際生產(chǎn)情況，將管坯溫度設置為25 ℃，空氣上方溫度為128 ℃，下方溫度為35 ℃，加熱爐外表面溫度為25 ℃，同時通過設置空氣與加熱爐之間接觸面的溫度來調(diào)整加熱爐的加熱溫度，兩個面之間引入間隙熱傳導來模擬面間的熱交換，其他面之間接觸均采用罰函數(shù)算法；在軟件的相互作用屬性模塊將管坯的六個自由度與管坯下表面的圓心約束在一起，便于拉伸力的施加；溫度場選用Coupled Temperature-displacement分析步計算，共計兩個分析步。第一個分析步為將坯料送至加熱爐中，第二個分析步用來穩(wěn)定管坯與加熱爐內(nèi)的溫度場。分析步設定界面中選擇瞬態(tài)求解，每個分析步為600 s，步長為固定步長類型，增量步長為10。而在應力應變場的計算則選擇Visco分析步，分析步步數(shù)和時間的設置均與溫度場一致。第一個分析步用于將管坯送入加熱爐內(nèi)，第二個分析步施加拉伸力；輸出選擇S、U、NT11等場變量，研究溫度和拉伸力對成品的外徑和壁厚的影響。

1.2 數(shù)值模擬方法的選擇

在毛細管拉制數(shù)值模擬過程中，涉及到溫度和變形，需要對應力應變場和溫度場兩個物理場進行耦合分析。ABAQUS中提供了順序耦合和完全耦合兩種求解耦合熱應力的方法。成型過程中，毛細管的應力應變場對溫度場的影響較小，故采用順序耦合熱應力分析，先求解加熱爐內(nèi)和毛細管的溫度場，再將溫度場的計算結(jié)果文件作為預定義場導入到應力應變場的計算中，計算出此溫度場下毛細管的應變。

數(shù)值模擬所采用的毛細管材料為石英，是一種熱流變材料，其粘度隨著溫度的升高而降低，在數(shù)值模擬過程中應對材料的粘彈性行為加以考慮[11]。熱流變材料滿足時溫等效性，即在較高的溫度和較短的應力松弛時間情況下獲得的應力松弛變量，與較低的溫度和較長的應力松弛時間情況下獲得的結(jié)果可認為是相同的。在ABAQUS軟件中的時溫等效模型為WLF方程，該方程適用于高分子聚合物的粘彈性分析過程，而Arrehenius公式對玻璃等非晶體材料的時溫現(xiàn)象的描述比WLF方程更準確，故采用Arrehenius定律來描述毛細管粘度隨著溫度的變化規(guī)律

(6)

式(6)中,T為溫度，η(T)為T時刻毛細管的粘度，Tref為參考溫度，取293 K，ηref為毛細管在Tref時的參考粘度，取9.565×108N·s/m2，α為與玻璃性質(zhì)本身有關的參數(shù)，取20 270[11]。

根據(jù)王炯[14]改進的Maxwell松弛模量表達式，有

(7)

式(7)中,τ為KWW松弛時間，τ=η/G，β為指數(shù)擴展因子，這里根據(jù)毛細玻璃管的特性，取β=0.5，G為剪切模量或體積模量。

綜合式(6)和式(7)，利用Fortran語言編寫相應的用戶子程序，完成變溫松弛模量的計算，并參考文獻[12]使用ABAQUS自帶的Prony級數(shù)對應力松弛模量實驗數(shù)據(jù)擬合得到粘彈性參數(shù)。

2 模擬結(jié)果處理及分析

在實際生產(chǎn)過程中，通過調(diào)節(jié)加熱爐溫度、管坯內(nèi)的氣體壓力和拉伸力來控制成品的幾何尺寸，本文根據(jù)實際生產(chǎn)情況，將加熱溫度T設定在1 100～1 400 ℃范圍內(nèi)，每隔100 ℃進行取值；管內(nèi)氣體壓力P保持在0 kPa；拉伸力F范圍在450～850 mg，每隔100 mg取一個值?；陧樞蝰詈蠠釕Ψ治龇ǖ奶攸c，先對加熱爐內(nèi)和管坯的溫度場進行數(shù)值模擬，再對毛細管拉伸過程中的應力應變進行分析。

2.1 溫度場模擬分析

數(shù)值模擬中將加熱爐與空氣相交的內(nèi)圓柱表面定義為加熱面，該面的溫度為加熱溫度，熱量從加熱爐內(nèi)圓柱表面逐漸向外擴散。下面以加熱溫度為1 200 ℃情況為例，對加熱爐爐內(nèi)的溫度場進行分析。

從圖3中可以看出，加熱爐內(nèi)溫度由內(nèi)到外開始上升，大約從30 s開始達到穩(wěn)定狀態(tài)，且達到最大值1 200 ℃。模型中將外表面設置為25 ℃，熱量由內(nèi)表面向外擴散，所以外表面溫度要低于內(nèi)表面。

圖3 加熱爐內(nèi)溫度分布圖Figure 3 Furnace temperature distribution map

為了更加直觀地對加熱爐內(nèi)不同位置處不同時間點的溫度進行分析，選取爐內(nèi)上中下三點(Node3 715、Node1 971、Node227)進行比較。結(jié)合圖3和圖4，可以看出，三點的溫度大約在30 s左右時，都達到了最大值。在大約100 s時，上部節(jié)點溫度開始急劇下降，這是因為管坯從上往下送進爐內(nèi)，對爐內(nèi)溫度有一定的影響；隨著時間的增加，管坯進入到中部，中部節(jié)點的溫度開始下降，大約到了600 s左右，第二個分析步生效，管坯不再向前行進，上中兩節(jié)點的溫度開始上升。管坯的運動位置并未達到加熱爐下部，故下節(jié)點溫度受管坯的影響較小，變化較小。

圖4 三節(jié)點溫度變化曲線Figure 4 Three nodes temperature curve

2.2 內(nèi)外徑分析

管坯進入加熱爐之后，對其施加不同的拉伸力和加熱溫度，以分析管坯的外徑(R1)、內(nèi)徑(R2)和應力狀態(tài)變化情況。為了更好地對管坯在t=1 200 s時內(nèi)外徑以及應力狀態(tài)進行分析，在ABAQUS后處理模塊中以成品外表面軸向設置路徑Path-1，該路徑各點與中心軸線的距離為管坯外徑R1；內(nèi)表面為Path-2，該路徑各點與中心軸線的距離為管坯內(nèi)徑R2。圖5為管坯熔融拉制之后的流變形狀圖。

圖5 管坯拉制流變形狀Figure 5 Diagram of capillary drawing rheological shape

采用控制變量法研究管坯在熔融拉制成型過程中，各影響因素與其徑向位移之間的相互關系。保持管坯在熔融拉制過程中的拉伸力在750 mg，管內(nèi)氣體壓力設置為0 kPa，改變加熱爐的加熱溫度，Path-1的軸向距離(本文用A表示)與內(nèi)外徑的關系曲線圖見圖6(a)，圖6(b)為仿真軟件計算云圖。結(jié)合圖6(a)和圖6(b)可以看出，當加熱溫度在1 400 ℃的時候，坯料內(nèi)外徑變化量最大，拉伸長度最長；溫度為1 100 ℃時，大約在管坯58 mm左右軸向距離處，加熱爐內(nèi)的溫度達到管坯軟化溫度，管坯在拉伸力作用下，徑向位移急劇變大，此時位于加熱爐中的管坯形成了外徑最小的頸區(qū)。由圖6(a)可知，隨著溫度的升高，頸區(qū)出現(xiàn)的位置逐漸提前；隨著管坯軸向距離的增加，內(nèi)外徑變化量趨于平穩(wěn)，管坯表面越來越平坦。

圖6 毛細管在不同溫度下的軸向距離—內(nèi)外徑曲線，外徑云圖Figure 6 The axial distance inside and outsidediameter curve of capillary tube at different temperatures, and the outerdiameter nephogram

圖7為陳珊[7]通過實驗得出的不同溫度下的毛細管內(nèi)外徑變化。由圖可知，保持其它條件不變的情況下，隨著溫度的升高，毛細管的直徑逐漸變小，這與本實驗中得出的結(jié)論一致。

圖7 不同溫度下毛細管內(nèi)外徑變化Figure 7 Variation of capillary inner and outer diameters at different temperatures

將爐內(nèi)溫度控制在1 200 ℃，管坯管內(nèi)氣體壓力設為0 kPa，只改變拉伸力的大小。圖8(a)為450～850 mg拉伸力下，Path-1的軸向距離與徑向位移的關系曲線圖。由圖知內(nèi)外徑隨著拉力的增加呈現(xiàn)逐漸減小的趨勢，當拉伸力為850 mg時，頸區(qū)中部的壁厚最小。為了更加清晰直觀地觀察到該現(xiàn)象，本文以管坯在不同拉伸力下的外徑為例，截取了管坯在不同拉伸力下的外徑云圖和與之相對應的頸區(qū)局部放大圖，如圖8(b)。在拉伸力逐漸增大的情況下，管坯的頸縮現(xiàn)象提前出現(xiàn)，頸區(qū)處的外徑逐漸減小，軸向長度(即拉伸長度)隨之增加。

圖8 毛細管在不同拉伸力下軸向距離—內(nèi)外徑曲線，外徑云圖Figure 8 Axial distance inner and outer diameter curve and outer diameter nephogram of capillary tube under different stretching force

保持其它條件不變的情況下，增加拉伸力，拉伸速度也將逐漸增大，所以在研究拉伸力變化情況下對毛細管直徑的影響時，可以用拉伸速度與管徑變化的對應關系來近似代替拉伸力與管徑變化的關系，它們的變化趨勢大致相同。圖9為陳珊[7]通過實驗得出的不同速度下的毛細管內(nèi)外徑變化。

圖9 不同速度下毛細管內(nèi)外徑變化Figure 9 Variation of inner and outer diameter of capillary tube under different velocity

2.3 Mises應力分布規(guī)律

圖10(a)和圖10(b)分別是Path-1在不同溫度和拉伸力下的Mises應力(本文用M表示)曲線圖，從圖10中可以看出，Path-1應力峰值第一次出現(xiàn)在軸向距離為0的位置。因為在設置邊界條件時，對該位置施加了固定約束，所以在拉伸力的作用下，該位置的應力值稍大于其他位置。第二次峰值出現(xiàn)在軸向距離大約50 mm位置處。圖10(a)中，毛細管拉制過程中的Mises應力峰值隨著溫度的增加而變大；峰值過后，不同溫度下軸向各點的應力值波動范圍不大，趨于平穩(wěn)。而在圖10(b)中，不同拉伸力下的Mises應力峰值隨著拉伸力的增加而變大；峰值過后，不同拉伸力下軸向各點的應力值相差較大，但都趨于平穩(wěn)，說明管坯軸向受力均勻。

圖10 Path-1在不同溫度和不同拉伸力下的Mises應力曲線Figure 10 Mises stress curve of path-1 under different temperature and different tensile force

模擬得到的應力分布云圖如圖11和圖12，圖11為管坯不同溫度下的Mises應力分布云圖。

圖11 不同溫度下的Mises應力分布圖Figure 11 Mises stress distribution under different temperatures

從圖中可以清晰的看到，Mises應力最大值達到了875 GPa，應力最大值出現(xiàn)位置在管坯左端面，此區(qū)域因設定為固定約束，故在拉制過程中Mises應力較大；隨著溫度的增加，管坯拉伸頸區(qū)出現(xiàn)的位置也提前出現(xiàn)，并且在頸區(qū)上方大約5 mm左右出現(xiàn)第二次應力峰值，由于該區(qū)域的溫度未達到管坯軟化溫度點，而下方區(qū)域的溫度達到軟化點，在拉伸力的作用下，出現(xiàn)第二次應力峰值；拉伸頸區(qū)處于變形與未變形區(qū)域之間的過渡區(qū)域，應力值要遠小于其他區(qū)域。同樣，為了清晰直觀地觀察頸區(qū)應力變化情況，截取了管坯在不同拉伸力下的Mises應力圖及局部放大圖，見圖12。該圖中，管坯的管徑在拉伸力增大的情況下變得越來越小，而軸向距離越來越大，拉伸頸區(qū)提前出現(xiàn)，在頸區(qū)位置處的應力比其他區(qū)域的要小，與溫度升高情況下出現(xiàn)的現(xiàn)象一致。

圖12 不同拉伸力下的Mises應力分布圖Figure 12 Mises stress distribution under different tensile forces

3 結(jié) 論

本文首次將順序耦合熱應力分析方法運用到毛細管拉制成型過程中，在管坯拉制過程中設置不同的加熱溫度和施加不同的拉伸力，通過對管坯的徑向位移和Mises應力的變化分析，得到如下結(jié)論。

1)在毛細管拉制成型過程中，其內(nèi)外徑與溫度和拉伸力成反比:隨著加熱溫度升高，內(nèi)外徑逐漸減小，沿管坯的軸向距離隨之增大；隨著拉伸力增大，管坯的外徑逐漸減小，徑向位移同樣會增大;反之亦然。

2)Mises應力值與溫度和拉伸力成正比，并且出現(xiàn)了兩次應力峰值；第一次出現(xiàn)在軸向距離為0 mm的位置，第二次峰值出現(xiàn)在拉伸頸區(qū)上方5 mm位置處。

3)加熱溫度升高和拉伸力增大會將拉伸頸區(qū)出現(xiàn)的位置提前，并且該區(qū)域與其他區(qū)域相比，Mises應力值明顯較小。

4)在毛細玻璃管拉伸過程中，可以通過控制溫度、拉伸力大小等來控制其內(nèi)外徑和軸向距離。隨著溫度的升高，其拉伸頸區(qū)位置提前，應力峰值出現(xiàn)的位置同樣會提前，在拉伸過程中應該根據(jù)環(huán)境溫度來確定其應力峰值的位置，保證其拉伸成型的穩(wěn)定性。