王超房，黃　明，石憲章，申長(zhǎng)雨，趙振峰

?

注塑成型Hele-Shaw流動(dòng)模擬中熱對(duì)流的異步長(zhǎng)求解

王超房，黃明，石憲章，申長(zhǎng)雨，趙振峰

（鄭州大學(xué)橡塑模具國(guó)家工程研究中心，河南 鄭州 450002）

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

關(guān)鍵詞:熱對(duì)流；傳熱；Hele-Shaw；注塑成型；數(shù)值模擬；聚合物加工

引 言

目前，注塑成型 CAE分析已經(jīng)成為模具制造前不可或缺的環(huán)節(jié)，其將模具設(shè)計(jì)和工藝參數(shù)設(shè)置建立在科學(xué)的分析基礎(chǔ)之上，對(duì)縮短生產(chǎn)周期、提高產(chǎn)品質(zhì)量具有重要指導(dǎo)作用[1-5]。然而，相關(guān)的CAE計(jì)算方法并非完美無(wú)缺，許多問(wèn)題仍有待于深入研究[6]，其中 Hele-Shaw流動(dòng)中的熱對(duì)流求解就是其中之一。

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

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

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

1　Hele-Shaw流動(dòng)控制方程

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

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

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

2　熱對(duì)流項(xiàng)（溫度場(chǎng)）求解

2.1 溫度場(chǎng)分步求解方法

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

即在給定的時(shí)間步內(nèi):① 不考慮熱對(duì)流和熱傳導(dǎo)，假定熔體處于靜止的條件下，求解方程（6）得到剪切熱引起的溫度場(chǎng)變化；② 不考慮熱傳導(dǎo)，只計(jì)算溫度場(chǎng)隨熔體運(yùn)動(dòng)產(chǎn)生的變化，即求解對(duì)流方程（7）；③ 然后再假定熔體處于靜止?fàn)顟B(tài)，只計(jì)算熱傳導(dǎo)，即只對(duì)每個(gè)節(jié)點(diǎn)（或單元）分別求解一個(gè)一維熱傳導(dǎo)方程。上述溫度場(chǎng)的計(jì)算采用與前面壓力場(chǎng)計(jì)算相同的有限元網(wǎng)格，該網(wǎng)格沿中面劃分，即中面網(wǎng)格。同時(shí)，對(duì)于每一個(gè)單元還要按一定的準(zhǔn)則在厚度方向分層，如圖1所示。同一單元的不同層以及不同單元的相同層可以有不同的層厚，并假定熔體除了前沿之外，只能在同一層內(nèi)流動(dòng)。溫度定義在單元的節(jié)點(diǎn)和各層的交界處，且速度、比熱容、剪切速率和黏度由壓力場(chǎng)計(jì)算得到。

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

2.2 對(duì)流方程異步長(zhǎng)求解

對(duì)于給定的時(shí)間步長(zhǎng)Δt，很容易得到式（6）的解

式（9）適用于每一層。

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

本節(jié)重點(diǎn)討論式（7）的求解。在一個(gè)速度為v的局部流場(chǎng)中，在時(shí)間步 tΔ內(nèi)，一個(gè)微粒由A點(diǎn)移動(dòng)到B點(diǎn)，如圖2所示，這樣在計(jì)算B點(diǎn)的當(dāng)前溫度時(shí)，如果不考慮其他因素，就需要找到該點(diǎn)的物質(zhì)在上一時(shí)刻所在的位置A，而A點(diǎn)在上一時(shí)刻的溫度就是B點(diǎn)當(dāng)前的溫度，即

在Hele-Shaw模型中，壓力是線性插值，通過(guò)式（4）只能計(jì)算得到單元中心速度。考慮到對(duì)于中面網(wǎng)格，某些轉(zhuǎn)角處的節(jié)點(diǎn)速度無(wú)法定義，因此本文算法將溫度定義在單元中心，再通過(guò)將單元溫度加權(quán)平均到節(jié)點(diǎn)，然后對(duì)節(jié)點(diǎn)溫度進(jìn)行單元插值形成溫度場(chǎng)?？紤]到流體運(yùn)動(dòng)是層流，求解熱對(duì)流過(guò)程主要就是在同一層內(nèi)沿速度的逆方向追蹤當(dāng)前相關(guān)點(diǎn)B在上一時(shí)間步的位置A，如圖3所示，然而在算法上該過(guò)程有兩個(gè)缺點(diǎn):一是在速度很大的區(qū)域，逆向搜索會(huì)穿越很多單元，計(jì)算很繁復(fù)；二是當(dāng)搜索穿越單元遇到筋板或熔接線等地方時(shí)，搜索方向會(huì)遇到歧路。

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

圖4　單元內(nèi)逆向追蹤Fig.4 Reverse tracking inside an element

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

其中，Li、Lj、Lk為 a點(diǎn)的面積坐標(biāo)，和由式（12）計(jì)算得到

但上述做法對(duì)注塑成型而言效率很低，因?yàn)樵诔尚瓦^(guò)程中不同區(qū)域的速度差異極大，有的點(diǎn)可能需要?jiǎng)澐稚习賯€(gè)子時(shí)間步，而有的點(diǎn)可能根本不需要。都采用統(tǒng)一的子時(shí)間步，會(huì)導(dǎo)致計(jì)算時(shí)間不必要的浪費(fèi)。為解決這個(gè)問(wèn)題，本文建議對(duì)不同的點(diǎn)采用不同的子時(shí)間步。眾所周知，在注塑流動(dòng)模擬中，總的時(shí)間步長(zhǎng)Δt并不是由溫度計(jì)算自己規(guī)定的，而是取決于流動(dòng)前沿節(jié)點(diǎn)的充填時(shí)間。如果對(duì)流計(jì)算中某個(gè)點(diǎn)能夠允許的最大時(shí)間步小于Δt，則必須設(shè)定自己的子時(shí)間步

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

上述計(jì)算必須要遵守次序，即必須先處理時(shí)間滯后的點(diǎn)。每計(jì)算完一個(gè)點(diǎn)的一個(gè)子時(shí)間步，就搜索下一個(gè)時(shí)間滯后的點(diǎn)，雖然算法不復(fù)雜，但會(huì)降低計(jì)算效率。為解決此問(wèn)題，本文采用二分法確定子時(shí)間步，即子時(shí)間步數(shù)滿足 2N0，相應(yīng)的子時(shí)間步長(zhǎng)為Δt/2N0，然后根據(jù)N0的大小對(duì)節(jié)點(diǎn)分組。實(shí)踐表明，N0通常不會(huì)超過(guò)10。這樣計(jì)算次序就變成先處理時(shí)間滯后的組，而組數(shù)是十分有限的。

通過(guò)上述算法得到某個(gè)時(shí)間步的熱對(duì)流解后，還要求解式（8），即對(duì)各個(gè)節(jié)點(diǎn)在厚度方向求解一次時(shí)域?yàn)?t,t0+Δt)的一維熱傳導(dǎo)方程。在厚度方向，將各層作為一個(gè)一維單元并進(jìn)行有限元離散，對(duì)時(shí)間后差分，則得到標(biāo)準(zhǔn)的三對(duì)角有限元方程

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

3　算 例

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

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

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

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

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

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

圖7　算例2的網(wǎng)格Fig.7 Mesh of Example 2

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

（1）如果將N在整個(gè)計(jì)算中取為定值1～8，計(jì)算會(huì)產(chǎn)生溢出。取為16時(shí)，計(jì)算可以完成。這表明如果在少量的時(shí)間段不滿足式(13)，可以進(jìn)行計(jì)算。但是，這取決于具體問(wèn)題和計(jì)算機(jī)系統(tǒng)，是一個(gè)很難定量給出的準(zhǔn)則，為了計(jì)算的穩(wěn)定可靠，對(duì)流計(jì)算的時(shí)間步長(zhǎng)應(yīng)該滿足式(13)。

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

（3）采用本文的異步長(zhǎng)算法，充填完成時(shí)體平均溫度和壓力分布如圖9所示。截止充填結(jié)束，總用時(shí)206 s，對(duì)流計(jì)算用時(shí)65 s。

圖8　使用統(tǒng)一子時(shí)間步長(zhǎng)充填完成時(shí)計(jì)算結(jié)果Fig.8 Computing results of Example 2 at filling end time by using uniform sub time step

圖9　使用異步長(zhǎng)充填完成時(shí)計(jì)算結(jié)果Fig.9 Computing results of Example 2 at filling end time by using variable time step

由于存在高剪切，熔體前沿溫度要遠(yuǎn)遠(yuǎn)高于熔體入口溫度，計(jì)算結(jié)果充分反映了這一點(diǎn)，與文獻(xiàn)[5]的論述是一致的。采用異步長(zhǎng)計(jì)算結(jié)果與采用統(tǒng)一的最小時(shí)間步長(zhǎng)結(jié)果相仿，但計(jì)算時(shí)間大大縮減，說(shuō)明本方法在保證計(jì)算精度和求解穩(wěn)定性的同時(shí)，可以明顯提高計(jì)算效率。

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

圖10　面窗中面網(wǎng)格Fig.10 Mid-plane mesh of transparent window

圖11　充填結(jié)束時(shí)刻體平均溫度Fig.11 Bulk temperature at filling end time

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

4　結(jié) 論

針對(duì)當(dāng)前CAE在Hele-Shaw流動(dòng)熱對(duì)流求解過(guò)程中存在的問(wèn)題，本文開(kāi)展了以下研究工作。

（1）采用分步法求解能量方程剪切項(xiàng)、對(duì)流項(xiàng)和熱傳導(dǎo)項(xiàng)，提出了采用變長(zhǎng)度子時(shí)間步長(zhǎng)方法處理對(duì)流項(xiàng)，確保搜索路徑局限在當(dāng)前單元內(nèi)，解決了逆向搜索路徑穿透多個(gè)單元時(shí)搜索可能失敗的問(wèn)題，并采用二分法確定子時(shí)間步數(shù)量，使得算法簡(jiǎn)潔有效。

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

符 號(hào) 說(shuō) 明

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

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

h——型腔半厚度，m

k——熱傳導(dǎo)系數(shù)，w·m-1·K-1

L ——單元面積坐標(biāo)

N ——子時(shí)間步數(shù)

p ——壓力，Pa

T ——溫度，℃

t ——時(shí)間，s

Δt ——時(shí)間步長(zhǎng)，s

v ——速度，m·s-1

x, y, z ——笛卡兒坐標(biāo)，m

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

δ ——權(quán)值

η ——黏度，Pa·s

ρ ——密度，kg·m-3

上角標(biāo)

n——子時(shí)間步

下角標(biāo)

i, j, k——單元節(jié)點(diǎn)編號(hào)

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] 黃明, 石憲章, 劉春太, 等. 基于統(tǒng)一網(wǎng)格的塑件成型與模具結(jié)構(gòu)一體化分析[J]. 化工學(xué)報(bào), 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] 王利霞. 基于數(shù)值模擬的注塑成型工藝優(yōu)化及制品質(zhì)量控制研究[D].鄭州: 鄭州大學(xué), 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] 王利霞, 申長(zhǎng)雨, 陳靜波, 等. 注塑成型充模過(guò)程的溫度場(chǎng)計(jì)算[J].計(jì)算力學(xué)學(xué)報(bào), 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] 李陽(yáng), 嚴(yán)波, 趙朋, 等. GLS/GGLS/SUPG 在三維注射成形充填模擬中的應(yīng)用[J]. 化工學(xué)報(bào), 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] 曹偉, 王蕊, 申長(zhǎng)雨. 注塑成型三維流動(dòng)的數(shù)值模擬 [J]. 計(jì)算力學(xué)學(xué)報(bào), 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] 耿鐵, 李德群, 周華民, 等. 注塑充模過(guò)程中溫度場(chǎng)的全三維數(shù)值模擬[J]. 化工學(xué)報(bào), 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] 劉永志, 趙振峰, 申長(zhǎng)雨. 基于分步法的三維注塑充填溫度場(chǎng)數(shù)值模擬[J].化工學(xué)報(bào), 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] 申長(zhǎng)雨.注射成型模擬及模具優(yōu)化設(shè)計(jì)理論與方法[M].北京: 科學(xué)出版社, 2009.

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

[17] 劉永志, 趙振峰, 馬蘭, 等. 注塑充填過(guò)程變厚度截面速度壓力場(chǎng)數(shù)值模擬[J].化工學(xué)報(bào), 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] 柳和生, 匡唐清, 周國(guó)發(fā), 等. 基于粘彈本構(gòu)模型的注射成型充填過(guò)程模擬研究[J]. 機(jī)械科學(xué)與技術(shù), 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] 林建忠, 阮曉東, 陳邦國(guó), 等. 流體力學(xué)[M]. 北京: 清華大學(xué)出版社, 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] 劉春太, 申長(zhǎng)雨, 陳靜波, 等. 注射模充模流動(dòng)和傳熱過(guò)程的理論與算法[J]. 高分子材料科學(xué)與工程, 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] 申長(zhǎng)雨, 郭恒亞, 趙振峰, 等. 模腔表面平均溫度邊界對(duì)注塑仿真的影響[J]. 鄭州大學(xué)學(xué)報(bào)(工學(xué)版), 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收到修改稿。

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

Received date: 2015-12-14.

中圖分類號(hào):TQ 320；TK 124

文獻(xiàn)標(biāo)志碼:A

文章編號(hào):0438—1157（2016）07—3047—08

DOI:10.11949/j.issn.0438-1157.20151900

基金項(xiàng)目:國(guó)家自然科學(xué)基金重點(diǎn)項(xiàng)目（11432003); 河南省高等學(xué)校重點(diǎn)科研項(xiàng)目(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