注塑成型Hele-Shaw流動模擬中熱對流的異步長求解

王超房，黃　明，石憲章，申長雨，趙振峰

?

注塑成型Hele-Shaw流動模擬中熱對流的異步長求解

王超房，黃明，石憲章，申長雨，趙振峰

（鄭州大學橡塑模具國家工程研究中心，河南 鄭州 450002）

摘要:在Hele-Shaw流動數值模擬中，速度是壓力的后處理結果。如果是點澆口，則澆口附近速度會隨單元尺寸縮小而趨于無窮大，導致能量方程作為一個整體求解時，時間步長必須非常小，否則會產生很大誤差；而根據熱對流物理意義分步求解，則需追蹤當前物質在上一時刻位置，當單元速度很高、逆向搜索需穿透多個單元時，搜索可能會失敗。鑒于此，基于分步求解法，研究提出一種變長度子時間步長方法處理對流項，確保搜索路徑局限在當前單元內，并采用二分法確定子時間步數量，使算法簡潔有效。算例表明，該方法在保證計算精度和求解穩定性的同時，可以明顯減少計算時間，提高計算效率。

關鍵詞:熱對流；傳熱；Hele-Shaw；注塑成型；數值模擬；聚合物加工

引 言

目前，注塑成型 CAE分析已經成為模具制造前不可或缺的環節，其將模具設計和工藝參數設置建立在科學的分析基礎之上，對縮短生產周期、提高產品質量具有重要指導作用[1-5]。然而，相關的CAE計算方法并非完美無缺，許多問題仍有待于深入研究[6]，其中 Hele-Shaw流動中的熱對流求解就是其中之一。

在注塑成型模擬中，由于熔體黏度嚴重依賴溫度，因此溫度場的計算至關重要。一般而言，溫度場的控制方程包含三要素:對流、傳導和能量轉換。所謂熱對流，就是承載一定溫度的物質由于自身運動而導致的溫度場變化，因此與速度場密不可分。在處理熱對流的數值方法中，大致有兩大類:一類是直接求解，即熱對流和熱傳導作為一個整體同時求解；另一類是分步求解，即熱對流從控制方程中分離出來單獨求解[7-8]。前者直接面對控制方程，可以認為是一個純數學問題，不用考慮其物理意義，基于該方法，文獻[9]開展了注塑成型Hele-Shaw流動溫度場的模擬，而文獻[10-13]則模擬了三維流動過程溫度場，為了抑制數值震蕩，上述研究均采用隱式格式離散方程，并使用上風類方法[14]處理熱對流項，但上風法僅考慮與節點相連的一層單元中的上游單元的貢獻，如計算中時間步長選擇不當，會導致節點速度與時間步長的乘積超越該節點周圍一層單元的空間尺度，造成求解誤差被無限放大，因此該時間步必須小于某個特定值，而特定值是由相關上游單元尺寸和節點速度共同決定的。針對上述問題，文獻[15]基于分步法，分別求解對流方程和熱傳導方程，實現了三維注塑充填溫度場數值模擬。

對于大部分注塑件，其幾何造型為薄殼結構，更適合用Hele-Shaw模型描述[16-18]，而該模型的直接求解對象是壓力，速度是通過節點壓力計算得到的。對于點澆口，由于澆口通常只能設置在一個有限元節點上，這導致在澆口附近，速度值會隨單元尺寸的縮小而趨于無窮大。這給溫度場計算帶來極大的困擾:除非采用極小的時間步，否則第一類求解方法必然產生很大誤差，而采用小時間步則會極大延長計算時間。第二類求解方法則是根據熱對流的物理意義，將其從控制方程中獨立出來，尋求當前時刻的物質在上一時間步的位置，從而得到熱對流的解，然后再計算熱傳導項。這種方法能有效避免第一類解法產生的問題，但算法相對復雜。如果某個單元的速度很高，則逆向搜索的路徑可能會穿透多個單元，使得計算復雜化，并且在筋板一類的地方逆向搜索可能失敗。

針對上述問題，本文提出一個基于第二類解法的異步長算法。在逆向搜索物質點運動軌跡時，對于速度高的質點，采用更多而更短的時間步，反之用較少同時較長的時間步，即將時間步分為若干子時間步，確保搜索路徑局限在一個單元內，采用二分法確定子時間步數量，使得此類解法更簡潔有效。并用算例驗證了該算法的有效性和可靠性。

1　Hele-Shaw流動控制方程

注塑充模過程屬于黏性不可壓縮流體的流動，基于Hele-Shaw模型假設[19-21]，熔體流動和傳熱控制方程可表述為

式中，v、p、η、ρ、cp、t、T、k和γ˙分別為速度、壓力、黏度、密度、比熱容、時間、溫度、熱傳導系數和剪切速率。對（2）式沿厚度方向(z向)積分，可以得到速度的解析式

式中，h+為型腔半厚度。將式（4）代入式（1），并在厚度方向積分，即可得到充模階段計算壓力的泊松方程

2　熱對流項（溫度場）求解

2.1 溫度場分步求解方法

利用熔體在注塑充模過程為層流且為對流占優的特點，Hele-Shaw流動模擬將熔體沿厚度方向分層，在不考慮熔體前沿的情況下，熔體只在同層內流動。與此相對應地，計算溫度時，認為對流只發生在同一層內，同時忽略同層內的熱傳導，而熱傳導只沿厚度方向進行。對于熱傳導，時間步長應該滿足步長與速度乘積在速度反方向上不超越該單元邊界的條件。否則，計算過程將是不穩定的（對于前差分），或產生較大誤差（對于后差分）。對于點澆口，這個步長限制會產生一個比較嚴重的問題，即當單元尺寸很小時，在點澆口周圍會產生很大的流動速度，從而時間步長必須非常小才能得到可靠的解。如果所有的點、所有的層都采用統一的最小時間步，將式（3）作為一個整體進行后差分求解，必然會非常耗時。針對以上問題，本文建議采用分步求解法，即根據能量方程的物理意義分別求解以下方程

即在給定的時間步內:① 不考慮熱對流和熱傳導，假定熔體處于靜止的條件下，求解方程（6）得到剪切熱引起的溫度場變化；② 不考慮熱傳導，只計算溫度場隨熔體運動產生的變化，即求解對流方程（7）；③ 然后再假定熔體處于靜止狀態，只計算熱傳導，即只對每個節點（或單元）分別求解一個一維熱傳導方程。上述溫度場的計算采用與前面壓力場計算相同的有限元網格，該網格沿中面劃分，即中面網格。同時，對于每一個單元還要按一定的準則在厚度方向分層，如圖1所示。同一單元的不同層以及不同單元的相同層可以有不同的層厚，并假定熔體除了前沿之外，只能在同一層內流動。溫度定義在單元的節點和各層的交界處，且速度、比熱容、剪切速率和黏度由壓力場計算得到。

圖1　三角形單元及分層Fig.1 Triangular element and layered

2.2 對流方程異步長求解

對于給定的時間步長Δt，很容易得到式（6）的解

式（9）適用于每一層。

圖2　熱對流的物理意義Fig.2 Physical meaning of thermal convection

本節重點討論式（7）的求解。在一個速度為v的局部流場中，在時間步 tΔ內，一個微粒由A點移動到B點，如圖2所示，這樣在計算B點的當前溫度時，如果不考慮其他因素，就需要找到該點的物質在上一時刻所在的位置A，而A點在上一時刻的溫度就是B點當前的溫度，即

在Hele-Shaw模型中，壓力是線性插值，通過式（4）只能計算得到單元中心速度。考慮到對于中面網格，某些轉角處的節點速度無法定義，因此本文算法將溫度定義在單元中心，再通過將單元溫度加權平均到節點，然后對節點溫度進行單元插值形成溫度場。考慮到流體運動是層流，求解熱對流過程主要就是在同一層內沿速度的逆方向追蹤當前相關點B在上一時間步的位置A，如圖3所示，然而在算法上該過程有兩個缺點:一是在速度很大的區域，逆向搜索會穿越很多單元，計算很繁復；二是當搜索穿越單元遇到筋板或熔接線等地方時，搜索方向會遇到歧路。

圖3　逆向追蹤Fig.3 Tracking in reverse direction of velocity

圖4　單元內逆向追蹤Fig.4 Reverse tracking inside an element

針對上述問題，本文提出一種變長度子時間步長求解方法，其基本思想為:將搜索路徑局限在當前單元內，如果給定的時間步長與速度之積要突破該單元，則將給定的時間步再劃分為若干子時間步，確保在子時間步內搜索路徑不會突破該單元。不失一般性，以圖4為例，在計算B點溫度時，根據上面準則將時間步劃分為N個等長的子時間步，而在每一個子時間步內的逆向搜索，都會從B點逆向行進到a點。在第n+1個子時間步，B點處的溫度為

其中，Li、Lj、Lk為 a點的面積坐標，和由式（12）計算得到

但上述做法對注塑成型而言效率很低，因為在成型過程中不同區域的速度差異極大，有的點可能需要劃分上百個子時間步，而有的點可能根本不需要。都采用統一的子時間步，會導致計算時間不必要的浪費。為解決這個問題，本文建議對不同的點采用不同的子時間步。眾所周知，在注塑流動模擬中，總的時間步長Δt并不是由溫度計算自己規定的，而是取決于流動前沿節點的充填時間。如果對流計算中某個點能夠允許的最大時間步小于Δt，則必須設定自己的子時間步

其中，d為待求溫度點到該點速度反方向與單元邊界交匯處的距離，N為一整數。仍以圖4為例，假定編號1的點（B點）需劃分的子時間步長為Δt/4，編號2、3、5、6需要的子時間步長為Δt/2，編號4的點不需劃分，即Δt。初始時刻t0的溫度在各節點已知，在經第一次計算后，得到了、和。當根據式（12）計算節點i溫度時，需知道2、 3、4、5、6各點在(t0+Δt/4)時刻溫度，該溫度通過對時間的插值求得

上述計算必須要遵守次序，即必須先處理時間滯后的點。每計算完一個點的一個子時間步，就搜索下一個時間滯后的點，雖然算法不復雜，但會降低計算效率。為解決此問題，本文采用二分法確定子時間步，即子時間步數滿足 2N0，相應的子時間步長為Δt/2N0，然后根據N0的大小對節點分組。實踐表明，N0通常不會超過10。這樣計算次序就變成先處理時間滯后的組，而組數是十分有限的。

通過上述算法得到某個時間步的熱對流解后，還要求解式（8），即對各個節點在厚度方向求解一次時域為(t,t0+Δt)的一維熱傳導方程。在厚度方向，將各層作為一個一維單元并進行有限元離散，對時間后差分，則得到標準的三對角有限元方程

其中，K為剛度陣，M為質量陣。處理完邊界項后，應用追趕法或三角分解很易得到方程的解，具體過程參見文獻[22]，這里不多贅述。

3　算 例

基于上述理論和算法，本文改進了自主研發“注射成型模擬仿真系統”（Z-Mold）的中面流動傳熱模塊，并以平板和某型號航天面窗為例開展模擬分析。需指出，前兩個算例只是為了突出本文算法的可行性，并不意味著一定符合工程實際設定（如澆口位置和注射時間）。此外，在極短時間內對模具內高速流動的熔體進行準確的溫度在線測量，目前還有一定困難，所以只定性地討論其合理性。

算例1 長、寬和厚度分別為200 mm、100 mm 和2 mm的平板，網格尺寸4 mm，點澆口注射，澆口設在長邊中點，如圖5所示。材料為PP（Adflex KS-221P），注射時間1 s，熔體溫度230℃，模具溫度50℃，由流率到壓力（V/P）控制設置為充填體積達到98%。

圖5　算例1的網格F i g . 5 M e s h o f E x a m p l e 1

圖6給出了本文方法（Z-Mold）模擬的體平均溫度。對于這樣的點澆口注塑，在開始一段時間熔體前沿總是呈半圓形狀，熱的熔體與模壁接觸時會在型腔表面形成冷凝層，熔體像“三明治”一樣被夾在冷凝層之間，熔體速度在厚度方向呈拋物線分布，即中心層速度最高，而且受模具傳熱的影響十分有限。同時由于剪切熱的存在，這樣會將一個略高于澆口溫度熔體帶到前沿，圖6（a）清楚顯示了這一過程。當熔體流動遇到邊界時，相應的部分流速會降低或停止，此時熱對流的作用降低或消失，熱傳導起主要作用，熔體溫度快速下降。熱對流在仍然流動的部分繼續以前的作用，熔體溫度在前沿保持最高，在前沿后方溫度下降也比靜止部分緩慢，圖6（b）清楚地展現了這一過程，說明本文算法給出的結果是合理的。

圖6　算例1體平均溫度Fig.6 Bulk temperature of Example 1

算例2 一長500 mm、寬200 mm、厚2 mm平板，點澆口注射，澆口設在短邊中點，網格剖分如圖7所示。材料為ABS，熔體溫度200℃，注射時間1 s，模具溫度40℃。本算例中，人為地構造一個高剪切場，用以考察本文算法對計算效率的改進。

圖7　算例2的網格Fig.7 Mesh of Example 2

如前所述，在注塑數值模擬中，每個時間步長是由充填程序決定的，如果此步長大于對流計算的步長要求[式(13)]，則需要在對流計算中選取合適的N減小步長。在本算例中，取決于每個步長的長短，式(13)中的N為2～32不等。

（1）如果將N在整個計算中取為定值1～8，計算會產生溢出。取為16時，計算可以完成。這表明如果在少量的時間段不滿足式(13)，可以進行計算。但是，這取決于具體問題和計算機系統，是一個很難定量給出的準則，為了計算的穩定可靠，對流計算的時間步長應該滿足式(13)。

（2）在每一個時間步中比較各個單元，選取滿足式(13)的最大的 N，即以最小的Δt/N作為本步計算的子時間步長。在充填完成時，體平均溫度和壓力分布如圖8所示。截止充填結束，總用時369 s，對流計算用時228 s。

（3）采用本文的異步長算法，充填完成時體平均溫度和壓力分布如圖9所示。截止充填結束，總用時206 s，對流計算用時65 s。

圖8　使用統一子時間步長充填完成時計算結果Fig.8 Computing results of Example 2 at filling end time by using uniform sub time step

圖9　使用異步長充填完成時計算結果Fig.9 Computing results of Example 2 at filling end time by using variable time step

由于存在高剪切，熔體前沿溫度要遠遠高于熔體入口溫度，計算結果充分反映了這一點，與文獻[5]的論述是一致的。采用異步長計算結果與采用統一的最小時間步長結果相仿，但計算時間大大縮減，說明本方法在保證計算精度和求解穩定性的同時，可以明顯提高計算效率。

算例 3 本文方法與其他計算模塊一起，成功分析、指導了某型號航天面窗的模具設計和工藝設定，其形狀類似于半徑145 mm的1/4球殼，網格尺寸10 mm，點澆口注射，如圖10所示。材料為光學級PC，注射時間4 s，熔體溫度268℃，模具溫度90℃，由流率到壓力（V/P）控制設置為充填體積達到98%。

圖10　面窗中面網格Fig.10 Mid-plane mesh of transparent window

圖11　充填結束時刻體平均溫度Fig.11 Bulk temperature at filling end time

圖11為模擬得到的體平均溫度分布。圖中流動前沿溫度最高，而澆口附近略低于熔體溫度。實際上，正如文獻[5]所論述的，注塑過程中流動前沿呈泉涌狀，其平均溫度在大多情況下是最高的，略高于入口溫度；而澆口附近，由于受模壁冷卻時間最長，同時不斷有新的熔體進入，該處平均溫度并不是最高。算例表明本文方法對于計算熱對流是有效、可靠的。

4　結 論

針對當前CAE在Hele-Shaw流動熱對流求解過程中存在的問題，本文開展了以下研究工作。

（1）采用分步法求解能量方程剪切項、對流項和熱傳導項，提出了采用變長度子時間步長方法處理對流項，確保搜索路徑局限在當前單元內，解決了逆向搜索路徑穿透多個單元時搜索可能失敗的問題，并采用二分法確定子時間步數量，使得算法簡潔有效。

（2）基于本文提出的算法，改進了自主研發“注射成型模擬仿真系統”（Z-Mold）的中面流動傳熱模塊。算例表明，Z-Mold模擬的體平均溫度在流動方向上有明顯的熱對流痕跡，剪切速率高、流速大的位置并未出現熱量聚集。在采用統一時間步計算對流時，如時間步大，不能保證計算的穩定（前差分）或精度（后差分），而時間步太小，又會耗時，本文方法有效地解決了該問題。這對于構造復雜、計算耗時的注塑模仿真，有著重要意義。

符 號 說 明

cp——比熱容，J·kg-1·K-1

d ——待求溫度點到該點速度反方向與單元邊界交匯處的距離，m

h——型腔半厚度，m

k——熱傳導系數，w·m-1·K-1

L ——單元面積坐標

N ——子時間步數

p ——壓力，Pa

T ——溫度，℃

t ——時間，s

Δt ——時間步長，s

v ——速度，m·s-1

x, y, z ——笛卡兒坐標，m

γ˙ ——剪切速率，s-1

δ ——權值

η ——黏度，Pa·s

ρ ——密度，kg·m-3

上角標

n——子時間步

下角標

i, j, k——單元節點編號

References

[1] MATIN I, HADZISTEVIC M, HODOLIC J, et al. A CAD/CAE-integrated injection mold design system for plastic products[J]. Int. J. Adv. Manuf. Technol., 2012, 63: 595-607.

[2] 黃明, 石憲章, 劉春太, 等. 基于統一網格的塑件成型與模具結構一體化分析[J]. 化工學報, 2012, 63(8): 2617-2622.

HUANG M, SHI X Z, LIU C T, et al. Integrated analysis of part molding and mold structural mechanics based on identical mesh[J]. CIESC Journal, 2012, 63(8): 2617-2622.

[3] YANG B, OUYANG J, LIU C T, et al. Simulation of non-isothermal injection molding for a non-Newtonian fluid by level set method[J]. Chin. J. Chem. Eng., 2010, 18 (4): 600-608.

[4] ZHOU H M, SHI S X, MA B. A virtual injection molding system based on numerical simulation[J]. Int. J. Adv. Manuf. Technol., 2009, 40: 297-306.

[5] 王利霞. 基于數值模擬的注塑成型工藝優化及制品質量控制研究[D].鄭州: 鄭州大學, 2004.

WANG L X. Process optimization and part quality control for plastic injection molding based on numerical simulation[D]. Zhengzhou: Zhengzhou University, 2004.

[6] WANG C F, HUANG M, SHEN C Y, et al. Warpage prediction of the injection-molded strip-like plastic parts[J]. Chin. J. Chem. Eng., 2015, doi: 10.1016/j.cjche.2015.12.012.

[7] DEMKOWICZ L, ODEN J T, RACHOWICZ W. A new finite element method for solving compressible navier-stokes equations based on an operator splitting method and h-p adaptivity[J]. Computer Methods in Applied Mechanics and Engineering, 1990, 84(3): 275-326.

[8] MINEV P, VABISHCHEVICH P N. An operator- splitting scheme for the stream function-vorticity formulation of the unsteady Navier–Stokes equations[J]. Journal of Computational and Applied Mathematics, 2016, 293: 147-163.

[9] 王利霞, 申長雨, 陳靜波, 等. 注塑成型充模過程的溫度場計算[J].計算力學學報, 2002, 19(2): 173-178.

WANG L X, SHEN C Y, CHEN J B, et al. Calculation of the temperature field in filling stage of injection molding[J]. Chinese Journal of Computational Mechanics, 2002, 19(2): 173-178.

[10] ZHOU H M, YAN B, ZHANG Y. 3D filling simulation of injection molding based on the PG method[J]. Journal of Materials Processing Technology, 2009, 204: 475-480.

[11] 李陽, 嚴波, 趙朋, 等. GLS/GGLS/SUPG 在三維注射成形充填模擬中的應用[J]. 化工學報, 2010, 61(2): 510-515.

LI Y, YAN B, ZHAO P, et al. GLS/GGLS/SUPG method in 3D numeical simulation of filling stage of plastic injection molding [J]. CIESC Journal, 2010, 61(2): 510-515.

[12] 曹偉, 王蕊, 申長雨. 注塑成型三維流動的數值模擬 [J]. 計算力學學報, 2005, 22(6): 792-795.

CAO W, WANG R, SHEN C Y. Numerical simulation for 3D flow in injection molding[J]. Chinese Journal of Computational Mechanics, 2005, 22(6): 792-795.

[13] 耿鐵, 李德群, 周華民, 等. 注塑充模過程中溫度場的全三維數值模擬[J]. 化工學報, 2010, 56(9): 1612-1617.

GENG T, LI D Q, ZHOU H M, et al. Simulation of temperature field in filling stage of injection molding using three-dimensional model[J]. CIESC Journal, 2010, 56(9): 1612-1617.

[14] BROOKS A N, HUGHES T J R. Streamline upwind/ Petrov-Galerkin formations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[J]. Computer Methods in Applied Mechanics and Engineering, 1981, 32(1/2/3): 199-256.

[15] 劉永志, 趙振峰, 申長雨. 基于分步法的三維注塑充填溫度場數值模擬[J].化工學報, 2011, 62(5): 1455-1459.

LIU Y Z, ZHAO Z F, SHEN C Y. 3D numerical simulation based on operator splitting method for temperature of plastic injection molding in filling stage[J]. CIESC Journal, 2011, 62(5): 1455-1459.

[16] 申長雨.注射成型模擬及模具優化設計理論與方法[M].北京: 科學出版社, 2009.

SHEN C Y. Injection Molding Simulation and Mold Optimization Design: Theory and Methods[M]. Beijing: Science Press, 2009.

[17] 劉永志, 趙振峰, 馬蘭, 等. 注塑充填過程變厚度截面速度壓力場數值模擬[J].化工學報, 2009, 60(7): 1818-1822.

LIU Y Z, ZHAO Z F, MA L, et al. Numerical simulation for velocity and pressure in filling stage of injection molded parts with uneven thickness[J]. CIESC Journal, 2009, 60(7): 1818-1822.

[18] 柳和生, 匡唐清, 周國發, 等. 基于粘彈本構模型的注射成型充填過程模擬研究[J]. 機械科學與技術, 2008, 27(5): 579-582.

LIU H S, KUANG T Q, ZHOU G F, et al. Numerical simulation of mold filling in injection molding of viscoelastic polymers[J]. Mechanical Science and Technology for Aerospace Engineering, 2008, 27(5): 579-582.

[19] 林建忠, 阮曉東, 陳邦國, 等. 流體力學[M]. 北京: 清華大學出版社, 2005: 85-86.

LIN J Z, RUAN X D, CHEN B G, et al. Fluid Mchanics[M]. Beijing: Tsinghua University Press, 2005: 85-86.

[20] HIEBER C A, SHEN S F. A finite element/finite difference simulation of the injection molding filling process[J]. J. Non-Newtonian Fluid Mech., 1980, 7: 1-32.

[21] 劉春太, 申長雨, 陳靜波, 等. 注射模充模流動和傳熱過程的理論與算法[J]. 高分子材料科學與工程, 2002, 18(6): 101-105.

LIU C T, SHEN C Y, CHEN J B, et al. Theory and algorithm of polymer melt flow and heat-transfer in injection-molding filling[J]. Polymer Materials Science and Engineering, 2002, 18(6): 101-105.

[22] 申長雨, 郭恒亞, 趙振峰, 等. 模腔表面平均溫度邊界對注塑仿真的影響[J]. 鄭州大學學報(工學版), 2005, 26(1): 9-12.

SHEN C Y, GUO H Y, ZHAO Z F, et al. Influence of the average temperature boundary condition on mold cavity surface to injection molding simulation[J]. Journal of Zhengzhou University (Engineering Science), 2005, 26(1): 9-12.

2015-12-14收到初稿，2016-03-19收到修改稿。

聯系人:黃明。第一作者:王超房(1982—)，男，博士研究生。

Received date: 2015-12-14.

中圖分類號:TQ 320；TK 124

文獻標志碼:A

文章編號:0438—1157（2016）07—3047—08

DOI:10.11949/j.issn.0438-1157.20151900

基金項目:國家自然科學基金重點項目（11432003); 河南省高等學校重點科研項目(15A430009, 15A430045)。

Corresponding author:HUANG Ming, huangming@zzu.edu.cn supported by the Key Program of National Natural Science Foundation of China (11432003) and the Key Research Project for Henan Universities (15A430009, 15A430045).

Thermal convection calculation with variable time step in Hele-Shaw flow simulation of injection molding

WANG Chaofang, HUANG Ming, SHI Xianzhang, SHEN Changyu, ZHAO Zhenfeng
(National Engineering Research Center for Advanced Polymer Processing Technology, Zhengzhou University, Zhengzhou 450002, Henan, China)

Abstract:In Hele-Shaw flow simulation, the directly solved variable is pressure and the velocity is only the post-treating result of pressure. Around the injection gate, the velocity may be very high along with reducing elemental size. This means that when the energy conservation equation is solved as a whole, the time step must be very short, otherwise, the error in the heat convection is unavoidable. The above problem can be overcome by using the operator-splitting method, in which the material at the current computing points needs to be tracked back to its position in the last time-step. However, this may lead to a new difficulty. If the elemental velocity is very high, the tracking needs to pass through a few elements and the reverse tracking may fail. To solve the problem, a new algorithm named sub time step with variable size was suggested to deal with the thermal convection in this paper, in which the sub time step using dichotomy was introduced that bounded the tracking path in a certain element. The new algorithm made the computation more simple and effective. The numerical examples showed that the new method had same accuracy as one using uniform small time step and high solving stability, but calculating time was dramatically reduced.

Key words:thermal convection; heat transfer; Hele-Shaw; injection molding; numerical simulation; polymer processing