應用晶格波茲曼兩相流模擬法於射出成型與覆晶封裝底膠填充流動研究

科技部補助專題研究計畫成果報告 期末報告

應用晶格波茲曼兩相流模擬法於射出成型與覆晶封裝底膠填充 流動研究

計 畫 類 別 ： 個別型計畫

計 畫 編 號 ： MOST 105-2221-E-006-094- 執 行 期 間 ： 105年08月01日至106年10月31日

執 行 單 位 ： 國立成功大學航空太空工程學系（所）

計 畫 主 持 人 ： 楊文彬 共 同 主 持 人 ： 陳介力

計畫參與人員： 碩士班研究生-兼任助理：楊為凱 碩士班研究生-兼任助理：陳螢萱

報 告 附 件 ： 出席國際學術會議心得報告

中　華　民　國　106　年　10　月　02　日

中 文 摘 要 ： 多組分和多相流的問題在自然界和工業中很常見，隨著科技與製程 進步，許多電子或生醫等元件也逐漸朝向微奈米化的方向製作，例 如微射出成型的應用。傳統的計算流體力學方法是基於求解巨觀統 御方程式來描述物理現象，因此在某些微奈米尺度的問題中較不適 用。晶格波茲曼法是一種基於分子動力論與統計物理學的數值計算 方法，已被發展用來研究熱流、電磁、聲學等多重物理問題。晶格 波茲曼法雖然也被用來研究多組分與多相流的問題，但在射出成型 流動的研究上較為稀少。由於此兩研究主題在工業製程中相當重要

，特別是在產品越來越微小化的趨勢下，有必要發展適合微奈米尺 度的數值模擬方法。

本計畫將發展晶格波茲曼兩相流模擬法，應用在射出成型模流 分析。第一年計畫主要推導晶格波茲曼的高分子非牛頓流體計算法

。之後則建構兩相流計算模型，結合高分子模型來進行射出成型的 模流分析，並著重於微尺度特徵元件的製程研究。並和其他數值方 法或商用套裝軟體結果作比較，以期建立一個適合應用在模流模擬 的數值方法。

中 文 關 鍵 詞 ： 晶格波茲曼法、微射出成型、單向自由表面

英 文 摘 要 ： Multicomponent and multiphase flow problems are common in nature and industry. As advances in technology and

manufacturing process, e.g. the micro injection molding, the size of electronic components or biomedical equipment are gradually reduced with the trend toward micro- and nano- scale. Based on solving macroscopic governing

equations for physical phenomena, the traditional methods of computational fluid dynamics may not applicable to some issues in micro- and nano- scale, in which the microscopic nature of matter dominates the physical transport

phenomena. In recent years, the lattice Boltzmann method (LBM), based on gas kinetic theory and statistical physics, is developed for the researches of fluid flows, heat

convection, electromagnetic, acoustic, and problems

coupling with multi-physics. The LBM has also been used to study multicomponent and multiphase flow problems, but there are few LBM researches about the injection molding process. Because these two topics are very important in the micro- and nano- scale manufacturing processes in industry and biomedical applications, it is necessary to develop a useful scheme, as the Lattice Boltzmann Method, for the numerical simulation of these problems. Furthermore, the LBM is also a promising method for multi-physics modeling and multiscale analysis.

In this project, the two-phase and non-Newtonian LBM models are constructed for the moldflow analysis in injection molding process. The work of the first year of this project is to develop the non-Newtonian LBM model for simulation of polymer fluids. In the followed study, the

calculation schemes of LBM on the two-phase flows will be developed. By the implementation of this project, the

lattice Boltzmann methods for multicomponent and multiphase flow are developed as an alternative scheme of numerical analysis in the micro injection molding process.

英 文 關 鍵 詞 ： lattice Boltzmann method, micro injection molding, free surface flow

科技部補助專題研究計畫成果報告

（□期中進度報告/■期末報告）

應用晶格波茲曼兩相流模擬法於射出成型與覆晶封 裝底膠填充流動研究

計畫類別：■個別型計畫 □整合型計畫 計畫編號：MOST-105-2221-E-006-094

執行期間： 105 年 8 月 1 日至 106 年 7 月 31 日 執行機構及系所：國立成功大學

計畫主持人：楊文彬 共同主持人：

計畫參與人員：

本計畫除繳交成果報告外，另含下列出國報告，共 _1__ 份：

□赴國外移地研究心得報告

□赴大陸地區移地研究心得報告

■出席國際學術會議心得報告及發表之論文

□國際合作研究計畫國外研究報告

中 華 民 國 105 年 8 月 19 日

I

目錄

中文摘要... II 英文摘要... III

一、 前言 ... 1

二、 研究目的 ... 1

三、 研究方法 ... 3

四、 結果與討論 ... 6

五、 結論 ... 9

參考文獻... 11 附件一國科會補助專題研究計畫成果報告自評表...

附件二國科會補助專題研究計畫項下出席國際學術會議心得報告...

II

中文摘要

多組分和多相流的問題在自然界和工業中很常見，隨著科技與製程進步，許多電子或生醫等元件 也逐漸朝向微奈米化的方向製作，例如微射出成型的應用。傳統的計算流體力學方法是基於求解巨觀 統御方程式來描述物理現象，因此在某些微奈米尺度的問題中較不適用。晶格波茲曼法是一種基於分 子動力論與統計物理學的數值計算方法，已被發展用來研究熱流、電磁、聲學等多重物理問題。晶格 波茲曼法雖然也被用來研究多組分與多相流的問題，但在射出成型流動的研究上較為稀少。由於此兩 研究主題在工業製程中相當重要，特別是在產品越來越微小化的趨勢下，有必要發展適合微奈米尺度 的數值模擬方法。

本計畫將發展晶格波茲曼兩相流模擬法，應用在射出成型模流分析。第一年計畫主要推導晶格波 茲曼的高分子非牛頓流體計算法。之後則建構兩相流計算模型，結合高分子模型來進行射出成型的模 流分析，並著重於微尺度特徵元件的製程研究。並和其他數值方法或商用套裝軟體結果作比較，以期 建立一個適合應用在模流模擬的數值方法。

關鍵詞：晶格波茲曼法、微射出成型、單向自由表面

III

英文摘要

Multicomponent and multiphase flow problems are common in nature and industry. As advances in technology and manufacturing process, e.g. the micro injection molding, the size of electronic components or biomedical equipment are gradually reduced with the trend toward micro- and nano- scale. Based on solving macroscopic governing equations for physical phenomena, the traditional methods of computational fluid dynamics may not applicable to some issues in micro- and nano- scale, in which the microscopic nature of matter dominates the physical transport phenomena. In recent years, the lattice Boltzmann method (LBM), based on gas kinetic theory and statistical physics, is developed for the researches of fluid flows, heat convection, electromagnetic, acoustic, and problems coupling with multi-physics. The LBM has also been used to study multicomponent and multiphase flow problems, but there are few LBM researches about the injection molding process. Because these two topics are very important in the micro- and nano- scale manufacturing processes in industry and biomedical applications, it is necessary to develop a useful scheme, as the Lattice Boltzmann Method, for the numerical simulation of these problems. Furthermore, the LBM is also a promising method for multi-physics modeling and multiscale analysis.

In this project, the two-phase and non-Newtonian LBM models are constructed for the moldflow analysis in injection molding process. The work of the first year of this project is to develop the non-Newtonian LBM model for simulation of polymer fluids. In the followed study, the calculation schemes of LBM on the two-phase flows will be developed. By the implementation of this project, the lattice Boltzmann methods for multicomponent and multiphase flow are developed as an alternative scheme of numerical analysis in the micro injection molding process.

Keywords : lattice Boltzmann method, micro injection molding, free surface flow

1

一、前言

在自然界跟工業應用中，兩相流是相當普遍存在的現象，如金屬鑄造模腔的填充過程(filling process in casting)與射出成型過程(injection molding process)。在這些流動過程，由於液體和空氣的密度差距很 大，因此也可將空氣對液體流動的影響忽略，亦即視為具有自由表面的流動。兩相流的流場現象對於 製程產品的優劣有很大的影響，例如在射出成型過程中，高分子聚合物的注射壓力、注射點的位置、

排氣通道的位置等參數，都會影響最終產品的品質。

在射出成型或底填封膠流動的過程中，流動的可視化相當困難，且模具的製作以及實驗的成本昂 貴，因此數值模擬在這些主題的研究中，扮演相當重要的角色。常見的數值方法如有限元素法[1-3]、

有限體積法[4-5]、以及有限差分法[6]等，都被用來進行底填封膠過程的研究。而商用軟體如 Moldflow 和ANSYS，也提供研究人員在模穴流動方面一個有用的模擬分析工具，例如 ANSYS 軟體就被一些文 獻使用來研究影響空氣空孔形成的因素[7]。

二、研究目的

由於科技進展迅速，許多工業或醫學元件的特徵尺寸都在朝微/奈米尺度邁進，因此對於微奈米尺 加工製程的精度和品質的掌握也更加重要，尤其是一些微小特徵的射出成型製程。傳統的計算流體力 學方法都是基於離散巨觀的統馭方程式來求解問題，對於基於分子尺度影響所造成的物理現象則未必 能精確描述。而晶格波茲曼法是基於分子運動力學的理論建構出來的數值模擬方法，在微奈米流道的 數值計算上，可扮演相當重要的角色。本專題計畫將研究晶格波茲曼法在射出成型模流分析方面的應 用。計畫目標主要推導發展LBM 的高分子非牛頓流體計算方法，並研究各種流體模型在複雜幾何邊界 或多孔介質中的流動現象。之後則建構LBM 的兩相流計算模型，並結合高分子流體模型來進行射出成 型的模流分析，研究的主題主要著重於微尺度特徵元件的製程研究。

國內外有關本計畫之研究情況

在具有移動界面的不混溶兩相流動中，由於介面上壓力和速度具有許多不連續性，使得求解其流 體統御方程式較為困難，此外，界面的運動行為也必須被準確的描述。處理介面移動問題的方法可分 為兩類：屬於拉格朗日法(Lagrangian)的表面追蹤法(surface tracking)，以及歐拉法(Eulerian)的表面捕捉 法(surface capturing) [8]。前者將界面視為具有不連續性，通常透過計算網格來追蹤界面的變化。後者 則先將統御方程式離散後，再將界面的變化納入考慮。在處理具自由表面的問題方面，流體體積法 (volume of fluid, VOF)[9]是一種廣泛被使用的歐拉法，可用來模擬模流充填過程[10-11]。而 Fyfe 等人[12]

使用的拉格朗日法以及Cleary 等人[13]使用的 smoothed particle hydrodynamics 法(SPH)則屬於追蹤流動 前緣界面的方法。

傳統的計算流體力學，在計算前必須將統御方程式進行離散，對於複雜邊界之處理更是困難，導 致前處理時間大大的增加。相較於傳統的數值方法，晶格波茲曼法(lattice Boltzmann method, LBM)不直 接求解Navier-Stokes 方程式中的巨觀(macroscopic)物理量，而是以流體的微觀(microscopic)模型為出發 點，依據統計熱力學與分子動力學之理論為基礎的波茲曼方程式，建構出一簡化的動力模型，藉著虛 擬粒子在此簡化模型中運作並且滿足特殊的守恆定律，反應出流體運動的微觀物理本質，並遵守流體 的巨觀物理現象，為介於巨觀與微觀之間的介觀(mesocopic)模型，並兼具二者的優點。晶格波茲曼法 的運作機制可分為碰撞(collision)及傳遞(propagation)，此時運動的是粒子的分佈函數。在所模擬的流體 與熱傳計算區域中，區域內離散為若干格點，每個格點依速度向量的多寡擁有個別分佈函數，該分佈

2

函數在該格點的離散晶格範圍內作碰撞機制後，透過格點與格點間的速度向量經傳遞機制傳遞到下一 格點，而透過計算格點的分布函數可得到巨觀物理量。晶格波茲曼法屬於顯性法則(explicit numerical method)，對於程式的邏輯與撰寫上，也較傳統計算流體力學來的容易許多。由於晶格波茲曼法具有明 顯的計算優勢，目前已在各領域得到廣泛的應用，如熱對流、熱波傳遞、多孔介質、磁流體、生物流 體等。獲得相當準確的計算結果。

射出成型是高分子流體的流動行為，屬於非牛頓流體。雖然晶格波茲曼模型針對非牛頓流體做相 關研究的起步較晚，但由於計算局部剪應力時，不需透過差分的方法求得速度場的空間導數，對於處 理複雜邊界條件也較傳統計算流體力學有利，所以近年來LBM 的非牛頓流體發展相當迅速。Aharonov 和 Rothman [14]用晶格波茲曼法去模擬局部鬆弛係數的非牛頓流體引起了極大的關注，而鬆弛係數可 由黏性係數計算得到，但計算過程為隱性計算法則必須使用迭代才可求出，造成計算上的負擔。Chen 等人[15]將分佈函數透過多重尺度展開後，利用差分的方式簡化後，剪切速率可利用平衡態分佈函數求 得，再透過剪切速率算出局部鬆弛係數。此法為顯性計算法則，避免了迭代過程中所耗費的計算時間。

近年來 LBM 方法在黏彈流體的數值模型上有比較多的發展，Su 等人[16]建構一個基於 Oldroyd-B 模型的LBM 方法來模擬黏彈流體在渠道流以及具轉角的緊縮流道內的流動現象，其計算結果與解析解 和其他文獻的實驗結果都很吻合，顯示其晶格波茲曼模型可用於計算較複雜的流體流變行為。Zou 等 人[17-18] 也 針 對 相 同 的 流 道 進 行 LBM 的 計 算 模 擬 ， 但 其 方 法 是 同 時 結 合 有 限 體 積 法 來 求 解 Navier–Stokes 方程式以及流體本構方程式，其黏彈流體模型則採用 Oldroyd-B 和 linear PTT 兩種模型。

其計算結果也是與解析解和其他文獻的數值結果相吻合。Zou 等人指出，這個混合 LBM 和 FVM 的方 法具有LBM 的高效率和穩定的優點，同時也保留 FVM 的精確度和廣泛適用性。Gupta 等人[19]則使用 Shan-Chen LBM 模型來模擬多重組分的流體現象，並結合有限差分法來計算流體內的高分子材料運動 行為。其方法可用來研究高分子流體的剪切、延伸、以及震盪等流動現象，將LBM 方法的適用性更推 廣一步。

晶格波茲曼法應用於多相流的模擬方面，主要可分為 color-fluid model、interparticle-potential model、

the free-energy model、以及 mean-field theory model 等四種模型[20]。 Shan 和 Chen 提出的單組分多相 流體模型(single component multiphase model, SCMP)以及多組分多相流體模型(multi-component and multiphase model, MCMP) [21-22]屬於 interparticle-potential model，為近年來較被廣泛使用來研究多重 組分或多相流的晶格波茲曼方法，特別是在包含多孔介質的流動研究方面。Huang 等人[23]使用 Shan–Chen (SC model)多相流晶格波茲曼模型來研究多孔介質內不互溶的液/氣兩相流動，得到良好的 結果，並指出SC model 具有簡單、容易處理多孔介質邊界等優點，並可以模擬具有高密度比的兩相流 體流動。

雖然近年來晶格波茲曼法在高分子流體流動以及多相流等方面的研究都有持續的發展，但關於應 用在射出成型方面的研究卻很少。Ginzburg 和 Steiner [24-25]較早使用 LBM 來研究金屬材料鑄造的模 穴流動，其模型結合Bingham 模型來描述黏塑性流體。該研究包含二維和三維的計算，考慮單一流體 在真空中的流動，亦即忽略空氣對金屬熔融液體的影響。其計算結果可準確模擬出流體與真空介面的 移動。在高分子材料的製程中，螺桿擠出機扮演相當重要的角色，Buick 和 Cosgrove [26]使用 LBM 來 研究螺桿擠出過程中，桿內的流體流動受剪切應力的影響行為，其研究結果雖然具有相當準確性，但 其模擬的模型僅是將螺桿擠出機內的流體混合腔室簡化為一個二維的頂蓋驅動流，並未考量三維的真 實螺桿擠出流動過程。不過，隨後Buick[27]更進一步的將此研究推展到非牛頓流體的研究，使用 LBM 方法結合power-law 模型來模擬擠出機內的非牛頓流體流動現象，使的 LBM 方法應用在高分子材料製 程的研究上更進一步。Latt 等人[28]使用 LBM 來模擬射出成型的流動過程，其研究採用 SC model 來模 擬高分子流體/空氣之兩相流，並與 Moldflow 軟體的模擬結果作比較，獲得不錯的研究成果。

3

三、研究方法

塑膠射出成型是一個相當複雜的高分子熱力與流動的變化過程，熔融塑料在經過充填、保壓、冷 卻過程中，需經歷溫度壓力、液相與固相變化等過程，因此使用 Hele-Shaw 流動將此複雜行為簡化，

從三維流動簡化成二維流動，其中假設熔融塑料於充填過程中為不可壓縮流，密度、比熱、熱傳導係 數為定值。模擬不可壓縮流的巨觀方程式為

0

⋅ =u (1)

  ^{1 2}

u uu p u

t 



        



   (2)

   

T Tu T

t 

     



 (3)

其中u

為熔融塑料的速度、t為時間、 為密度、p為壓力、 為黏度、T為能量和 為熱擴散係數。

3.1 晶格波茲曼法與熱晶格波茲曼法

本文使用的數值方法為LBM，主要探討二維問題，因此離散速度場使用 D2Q9 模型[29]，如圖 1。

圖1 D2Q9 模型

在不考慮外力作用的晶格 Boltzmann 方程式與忽略黏滯熱耗散與壓力功的熱晶格 Boltzmann 方程式分別如下

       

       

, , 1 , ,

, , 1 , ,

eq

i i i i i

eq

i i i i i

g

f x c t t t f x t f x t f x t

g x c t t t g x t g x t g x t





 

         

 

         

    

     (4)

其中9 個方向的平衡態分布函數fi^eq與gi^eq之定義分別為

2 2

2 4 2

2 2

2 4 2

( )

1 2 2

( )

1 2 2

eq i i

i i

s s s

eq i i

i i

s s s

c u c u u

f c c c

c u c u u

g

c c c

 

 

   

     

 

   

     

 

    

  

    

  

(5)

其中 ₀ ⁴

 9，當i1~ 4，_i取¹

9 ，當i5 ~ 8，_i取 ¹

36，透過求解晶格的分布函數 f_i與g_i的演化得到流

體巨觀的密度 _i

i

f ^{、動量 i i}

i

u c f

^^{ 與能量 i}

i

g ^{，且T ，而其中}c^i為晶格的速度向量，其數 值為

 0,0 0 ( 1) ( 1)

cos ,sin 1 ~ 4

2 2

( 5) ( 5)

2 cos ,sin 5 ~ 8

2 4 2 4

i

i

i

c i

i i

c i

i i

c i

 

   

 

 

 

  

      

        







(6)

4

晶格 Boltzmann 方程式與平衡態分布函數經過 Chapman-Enskog 展開可以推得巨觀方程式(1~3)，

推導過程中可得黏度的定義以及熱擴散係數與普朗特數Pr 的定義分別為

 

 

2 2

0.5 0.5 2 1

Pr 2 1

s

s g

g

c t

c t

 

 

 

 

  

  

  



(7)

3.2 邊界條件設定

根據邊界格式的處理方法，速度邊界可分為平直邊界條件及曲面邊界條件，但在此研究中不會使 用到曲面邊界，而平直邊界條件又可分為反彈邊界、修正反彈邊界、半反彈邊界。Ziegler [30]於 1993 年對邊界條件做討論。而本研究主要使用半反彈邊界條件（half bounce-back rule）半反彈邊界的原理為 流 向 邊 界 的 粒 子 在 下 一 步 時 ， 該 粒 子 會 沿 著 原 方 向 反 彈 至 計 算 的 網 格 中 間 之 節 點 ， 可 表 示 為

( ,_b ) _i( _f, )

fi x t   t f x t ，其中i 與i 呈相反方向(i  i)，x^_b是壁面格點，x_f 是流體格點，此邊界條件 具有二階精度。

在壓力驅動的無限管道中的流動問題，可將壓力轉換為外力，但是由於壓力梯度在許多情況下無 法轉換為外力項，因此Zou 跟 He [30]提出了 D2Q9 模型的速度邊界與壓力邊界，根據巨觀物理量與分 布函數的關係，在邊界格點上定義出滿足未知數的方程式。速度邊界與壓力邊界的方法大致相同，其 差異在於速度邊界的設定上需給定速度，而壓力邊界的設定上則是需給定密度。舉例來說，入口邊界 條件設定為速度邊界，如圖2，由於左邊無分布函數傳入，因此f₁、f₅、f₈、 為未知的分布函數及密 度，而其他的 f₀、 f₂、f₃、f₄、 f₆、 f₇與給定入口處的速度u_x、u_y皆為已知的分布函數，因此可推得 入口處的未知量為

   

1 3

2 4 9 3 6 7

2 3

1 2

1

x

x

f f u

f f f f f f

u





 

 

        

 

 

5 7 2 4

8 6 2 4

1 1 1

2 2 6

1 1 1

2 2 6

y x

y x

f f f f u u

f f f f u u

 

 

    

    

(8)

圖2 速度邊界的入口圖示

熱模型的邊界條件則採用Shan [32]提出的被動標量模型中的邊界設定方法，等溫邊界條件設定為

( , ) 2 ( , )

i i i

g x t   g x t (9)

由於在邊界上的速度為零，亦可定義為

( , ) ^{eq eq} ( , )

i i i i

g x t g g g x t (10)

其中i 與i呈相反方向(i  i)。

本文中的邊界條件設定方法皆利用Chang [33]於 2011 年提出一種未知分量指標算法，亦可稱為格 點判別法，此方法用來增進LBM 在求解流場邊界之未知量時的程式撰寫效率，與傳統的半反彈方法項 較之下，除了精確度更高以及合理的數解析理論依據，亦可維持晶格波茲曼法簡單處理的特性。

5

3.3 單向自由表面 LBM 模型

Korner 與Thurey [34]提出的單向自由表面 LBM 為一種模擬兩相流的方法，深受廣大學者的青睞，

其類似於VOF 方法，在此方法中，採用體積函數來表示一個格子中所占用的液體之比例。而在使用此 方法對兩相流的模擬時，必須假設氣相與液相之間沒有接觸，使用自由表面將兩者分開，確保氣液兩 相完全分開，自由表面必須是完全封閉的。在使用此方法時，為了滿足以上假設，需要將流體格點計 算域劃分成3 個區域，如圖 3。

圖3 自由表面區分格點的示意圖

其中分別是液體、氣體、介面區域分別標示成Fluid、Interface、Gas。依照介面格點計算的過程，

主要分成三個步驟，分別是重建分布函數、介面移動量的計算、格點標示重新更新。

重建分布函數（reconstruction）主要由於在進行自由表面流動的計算時，只有液相格點和介面格點 有進行碰撞和傳遞的計算，而氣相格點的部分忽略，因此介面格點處從氣相傳遞過來的未知分布函數 需要透過重建分布函數的方式來求得，其定義為

 ^,  i^eq G^,  ^eq G^,  i ^,

i i

f x t   t f  u  f  u  f x t (11) 其中_G1，_G為氣相密度，亦為模擬中的參考密度。

介面移動量的計算是為了重新標示經過求解晶格摸茲曼方程所發生之碰撞傳遞後的格點類型，需 要計算先非氣相格點之間的質量變化。假設介面格點位置是x

，其相鄰的格點為x tc_i

，該介面格點與 相鄰格點之間的質量變化用分布函數來表示為

     

   

, ,

, 2

, ,

i i

i i

i

x tc t x t m x t t

f x tc t f x t

   

   

    

 

  



   (12)

其中m為質量、 為流體體積函數，兩個參數的關係式為

       , , , , m x t  x t  x t V x t 

(13)

其中^{V x t} ^{, 為體積且V 1}。液相格點被液體充滿，所以液相格點的液體體積函數 恆為1，且在體積為

1 的情況下，質量與密度相等，即m；氣相格點沒有液體，因此氣相格點的液體體積函數 恆為0；

對介面格點來說，其格點內部部分充滿液體，所以流體體積函數 介於0 到 1 之間，即0  1。 在求得介面格點的質量變化後，對應介面格點的質量可由下式求得

( , ) ( , ) _i( , )

i

m x t^   t m x t^ m x t^  t ⁽¹⁴⁾

其表示為介面格點隨時間演化的質量。

格點標示重新更新指的就是在每一個時間步階之分布函數傳遞後，要重新標示每個格點。液相格 點與氣相格點之間不會直接做變換，需透過介面格點做變換。在本研究中，設計由判定介面格點所含 之質量m來做格點標示的變換，判定的規則為

( , ) ( , ) for filled cell ( , ) 0 for empty cell

m x t t x t t

m x t t



    

  

 

 (15)

對於填滿的格點而言，當格點轉變成液相格點後，為了使介面保持封閉，即避免氣相與液相格點直接 接觸，須判斷液相格點周圍是否有氣相格點，若是有氣相格點則標示為新的介面格點，這樣就能確保

6

介面的封閉性。

3.4 自由表面的邊界處理

Xiu Qing Xing [35]等人提出在介面格點與牆面接觸時必須特別處理，當介面格點接觸到固體格點 的邊界，而邊界格點（boundary cells）是介面格點時，使用先前提到的重建分布函數（reconstruction）

的方式來計算從固體格點傳入介面格點的離散函數，其可滿足邊界條件的設定。介面格點的質量已填 滿時，格點經過重新標示後，變為液體格點，而當此格點接觸到固體邊界時，則進行半反彈邊界，利 用此模式來做邊界的計算。

熱模型中引入自由表面的文獻不多，主要多為利用模擬加熱粉末之相變化的應用。Ammer[36]此篇 亦是使用Shan 提出的被動標量模型，加入自由表面控制氣相、液相及固相之間複雜的交互關係，模擬 EBM 製程情形。文中提到，在介面格點的溫度邊界設定上，設定反射條件為熱邊界條件，即設定未知 的離散函數是由氣相格點傳入介面格點，而此未知的離散函數近似於液相格點傳入介面格點之離散函 數。在介面格點邊界設定上使用此設定可滿足流體於流道中之設定條件。

3.5 修改熱晶格波茲曼格式

由於熔融塑料的黏度很高，黏度與普朗特數（Prandtl number，Pr）成正比，因此 Pr 相對的也會很 高。本研究使用Xiang [37]於 2012 年提出的修改熱晶格波茲曼格式（modified lattice Boltzmann scheme）

來修改原本Shan 提出的熱晶格波茲曼模型，可降低因鬆弛時間_g所導致LBM 的不穩定性及模擬性能，

為了克服極限，文獻中引入一個新參數 於熱擴散方程式(7)中，公式變為

 

2 0.5

s g

c t

     (16)

由此公式可推得新的平衡態分布函數

 ^{2 2}  





2 4 2 2

1 2

1 2 2 2

i s

i

eq i

i i

s s s s

c c

c u c u u

g T

c c c c

 

     

 

      

 

 

 

    

    (17)

Xiang[6]比較使用原格式(  )與修改的格式(1   )做比較，最後求解的結果修改的格式可獲得較1 低的誤差。因此在本研究中，使用其能夠改善穩定性這點，來做為改善當Pr 提高時的不穩定性，並參 考文獻中提到定義新參數 值的方向，讓鬆弛時間g落在穩定範圍內，使結果不會發散。

四、結果與討論

本文將依照速度與溫度來說明，速度場的部分，以二維自由表面位移及質量守恆的公式，來驗證 程式的準確性，最後以單向自由表面LBM 模擬高黏度塑料填充，探討流體於模穴中的流動情形。

4.1 程式驗證

首先本文為了確保自由表面之程式，邊界條件設定，入口速度為均勻流速流入，並設定為速度邊 界條件，且入口速度已知；出口設定為壓力邊界條件格式；上下兩壁面設定為無滑動的半反彈邊界條 件格式。利用質量守恆來探討流體流入經過一段時間後，流道內所佔之體積是否符合流體力學的公式，

因此本文利用流體力學中質量守恆的公式與體積流率Q 來推算充填體積為

V  AVavg t (18)

其中V_avg表示為入口平均速度、A表示入口之截面積、V 為體積。但由於本文為二維模擬，因此將此 流場的深度設定為單位長度，上式將同除以深度得

covered avg

A HV t (19)

7

由此公式即可求得經過時間變化 t 後流體充填的面積A_covered。晶格單位與實際單位值相對應之比例參 數分別為，比例長度C h 4 10 m^5 、比例時間C t 8.83 10 s ^9 、比例速度Cu 4530.25 / sm 。驗證計 算200000 步時流體的充填面積，可由步數乘以比例時間求得。

圖4 速度場含自由表面之程式驗證的流道格點數

模擬的流道由兩塊平行板構成，計算格點數為25 250 （晶格單位的寬與長）且鬆弛時間設定為

1，但模擬之初始狀態流道內為氣相，入口為均勻流速流入，模型如圖 4。邊界條件設定為入口速 度為均勻流速流入，設定為速度邊界條件；出口設定為壓力邊界條件格式；上下兩壁面設定為無滑動 的半反彈邊界條件格式。透過單位轉換的計算推出，晶格單位下所對應之真實域的材料參數，如表1。

表1 速度場程式驗證的流道參數

在時間 t 後，式(20)可求得解析解之流體填充面積A_a  4 10^6。透過程式計算出流體格點數為2484 格，單一格點之面積計算為比例長度的平方C h² 1.6 10 m ^{9 2}，將兩數值相乘後可得模擬解之流體流

經面積A_s3.97 10 ^6。計算之誤差為0.64%，誤差量非常小。LBM 數值模擬的結果與理論結果非常接

近，說明利用本文所述之方法模擬射出成形的流動情形準確性很高。

溫度場的部分，本文將引入Shan 所提出的被動標量模型，並結合自由表面，同時分析不同普朗特數

下高黏度流體的熱傳問題。材料參數設定與邊界條件設定和前面塑料在流道內之填充的部份一樣，溫度

的部分設定入口為高溫 353K，其餘邊界皆為低溫 293K。模穴尺寸設計透過單位轉換的計算推出，晶 格單位下所對應之真實域的材料參數，如表2。

表2 溫度場程式的流道參數

如設定Pr 100 ，將原TLBM 及修改後 TLBM 兩結果放在一起做比較，繪出在相同時間步的結果，

圖5 中為每 30000 步間距的結果。先前的格式  的結果如圖 5 的左側，修正的格式1   的結果如圖1 5 的右側。

從結果來看兩者的溫度分布範圍皆介於最低溫與最高溫之間，而從溫度分布雲圖來看，每個時間 步的結果幾乎一模一樣，由此可見，改變 值計算出來的結果與原本結果一致，接著利用此數值的改變，

計算高Pr 的結果。

由於當Pr 值大於 500 時，先前的格式  的結果可能就會發散，因此本研究使用改變1  值，讓 鬆弛時間g落於穩定範圍內，而在此參數設定下，鬆弛時間g只要設定大於0.5001 且小於 2，結果較 不易發散。本研究使用在相同模型下的高Pr 數值結果整理分成 Pr 1000 、Pr 5000 、Pr 6000 三個 結果， 的設定分別為0.8、0.6、0.55，分別對應的鬆弛時間為 0.5006、0.5001、0.5001，結果如圖 6。

圖6 反映了流體在不同 Pr 下的流動範情形，顯示了不同時刻自由表面的位置與溫度分布情形，圖中可 以看到當 Pr 越大，流體中間的溫度分布越大，即中心紅色的區域越大。高 Pr 值即表示動黏滯係數遠 大於熱擴散係數，因此熱量主要以對流的方式傳遞，因此透過每個時刻的溫度分布雲圖來看，很明顯 Pr 越高傳播熱量的效果較好。利用此方法應用於計算高 Pr 值，成功模擬出高 Pr 值的結果。本研究不 僅成功地將自由表面結合溫度場的計算，應證了自由表面計算方法的能力，並且成功利用新參數 適 當調整Pr 值，使其能夠模擬到 Pr 為 6000。

8

圖5 比較新參數的影響

圖6 相同模型下，高 Pr 值的溫度分布雲圖

4.2 模擬微射出成形之填充過程

本研究的設計圖如圖7，實際尺寸（毫米 mm）與晶格單位的比例為 1：25，邊界條件設定和前面 驗證的部份一樣，材料參數設定Re 為 0.075、Pr 為 6000、 為 0.55。

圖7 流道幾何模型的尺寸

9

此部分主要觀察塑料流到微流道的情形，圖 8 左為溫度場分布情形、右為速度場分布情形。當塑 料流入微流道之後，溫度明顯逐漸下降，而速度值也明顯變慢，當當流體第二及第三個微流道時，第 一個微流道內的速度變為 0，而從溫度場來看，明顯的是此流道已被流體填滿，從時間步階為 700000 來看，流體已充滿於各個微流道中。

圖8 塑料流到微流道的情形

最後由於我們使用格點判別法[10]與自由表面結合，因此可以匯出各時間步階每個格點的格點標示形 式，由此可以得到個時間步階的介面移動位置，如圖9。其中橘色的格點為流體格點，標示為紅色的 格點即為介面格點。經由格點判別法與自由表面的結合，能夠明確的得到各時間步階下的介面移動位 置。

圖9 介面移動變化圖

五、結論

本研究採用晶格波茲曼法模擬射出成型的填充過程，不僅成功的利用自由表面方法模擬塑料填充 的過程，並利用自由表面方法結合被動標量模型，成功模擬出含有介面之速度場與溫度場的結果，最 後將其引入新參數 ，不僅能夠得到良好的數值結果，並能得到高普朗特數Pr 的模擬結果。

將自由表面引入LBM 中，依照射出成型的流動行為設計程式的邏輯運算，並使用格點判別法簡化 複雜的計算方式，不僅大大簡化計算時間，從驗證的結果來看，排除人為計算誤差可以得到極小的誤 差值。

10

溫度場的部分則是在被動標量模型中引入自由表面，其同樣使用格點判別法簡化複雜的程式，並 加入新參數 ，利用此數值與鬆弛時間的關係作為調整鬆弛時間，使其落於穩定範圍內，即

2_g0.5001，由此提升程式的穩定性及準確性，並利用此方法解決無法模擬高Pr 值的流動情形，最

高能夠模擬到Pr 6000 。在研究的過程中發現，在流體以均勻流速流進流道後，在流體受到黏滯力影響 逐漸變為完全發展流的過程中最容易使結果發散，所以只要此過程中流體沒有發散，就能夠得到穩定 的結果。

最後，將晶格波茲曼法之速度場與溫度場應用於微射出成型的填充過程中，不僅成功將單向自由 表面LBM 應用到微射出成型上，更模擬出高黏度的流場流動情形。由此說明了單向自由表面 LBM 在 模擬微射出成型的流動有很大的發展潛力，為數值模擬射出成型開闢了新道路。

11

參考文獻

[1] H. Wang, H.M. Zhou, Y. Zhang, D.Q. Li, K. Xu, Three-dimensional simulation of underfill process in flip-chip encapsulation, Comput. Fluids 44 (2011) 187–201.

[2] W.B. Young, Non-Newtonian flow formulation of the underfill process in flip-chip packaging, IEEE Trans. Compon. Packag. Manuf. 1 (2011) 2033–2037.

[3] J.W. Wan, W.J. Zhang, D.J. Bergstrom, Numerical modeling for the underfill flow in flip-chip packaging, IEEE Trans. Compon. Packag. Manuf. Technol. 32 (2009) 227–234.

[4] S.W. Moon, Z. Li, S. Gokhale, J. Wang, 3-D numerical simulation and validation of underfill flow of flip-chips, IEEE Trans. Compon. Packag. Manuf. 1 (2011) 1517–1522.

[5] C.Y. Khor, M.Z. Abdullah, M.A. Mujeebu, F.C. Ani, FVM based numerical study on the effect of solder bump arrangement on capillary driven flip chip underfill process, Int. Commun. Heat Mass Transfer 37 (2010) 281–286.

[6] T. Hashimoto, T. Shin-ichiro, K. Morinishi, N. Satofuka, Numerical simulation of conventional capillary flow and no-flow underfill in flip-chip packaging, Comput. Fluids 37 (2008) 520–523.

[7] Aizat Abas, M.H.H. Ishak, M.Z. Abdullah, F. Che Ani, Soon Fuat Khor, Lattice Boltzmann method study of bga bump arrangements on void formation, Microelectronics Reliability, in press, 2015

[8] Irina Ginzburg, Konrad Steiner, Lattice Boltzmann model for free-surface flow and its application to filling process in casting, Journal of Computational Physics 185 (2003) 61–99

[9] C.W. Hirt, B.D. Nicholls, Volume of fluid (VOF) method for the dynamics of free boundaries, J. Comput.

Phys. 39 (1981) 201.

[10] A.N. Alexandrou, E. Duc, V. Entov, Inertial, viscous and yield stress effects in Bingham fluid filling of a 2-D cavity, J. Non-Newtonian Fluid Mech. 96 (3) (2001) 383.

[11] D.M. Gao, A three-dimensional hybrid finite element-volume tracking model for mold filling in casting processes, Int. J. Numer. Methods Fluids 29 (1999) 877.

[12] D.E. Fyfe, E.S. Oran, M.J. Fritts, Surface tension and viscosity with Lagrangian hydrodynamics on a triangular mesh, J. Comput.Phys. 76 (1988) 394.

[13] P. W. Cleary, J. Hat, V. Ahuja, High pressure die casting simulation using Smoothed Particle Hydrodynamics, J. Cast Metals Res. 12 (2000) 335.

[14] E., Aharonov and D. H., Rothman, Non­Newtonian flow (through porous media): A lattice­Boltzmann method. Geophysical Research Letters, 20(8), 679-682,1993.

[15] Y.-L. Chen, X.-D. Cao, , and K.-Q. Zhu, A gray lattice Boltzmann model for power-law fluid and its application in the study of slip velocity at porous interface. Journal of Non-Newtonian Fluid Mechanics, 159(1), 130-136,2009.

[16] Jin Su, Jie Ouyang, Xiaodong Wang, and Binxin Yang, Lattice Boltzmann method coupled with the Oldroyd-B constitutive model for a viscoelastic fluid, PHYSICAL REVIEW E 88, 053304 (2013)

[17] Shun Zou, Xue-Feng Yuan, Xuejun Yang, Wei Yi, Xinhai Xu, An integrated lattice Boltzmann and finite volume method for the simulation of viscoelastic fluid flows, Journal of Non-Newtonian Fluid Mechanics 211 (2014) 99–113

[18] Shun Zou, Xinhai Xu, Juan Chen, Xiaowei Guo, and Qian Wang, Benchmark Numerical Simulations of Viscoelastic Fluid Flows with an Efficient Integrated Lattice Boltzmann and Finite Volume Scheme,

12

Advances in Mechanical Engineering February 2015 vol. 7 no. 2 805484

[19] A. Gupta, M. Sbragaglia, A. Scagliarini, Hybrid Lattice Boltzmann/Finite Difference simulations ofviscoelastic multicomponent flows in confined geometries, Journal of ComputationalPhysics291(2015)177–197

[20] Cyrus K. Aidun and Jonathan R. Clausen, Lattice-Boltzmann Method for Complex Flows, Annu. Rev.

Fluid Mech. 2010. 42:439–72

[21] X. Shan, H. Chen, Lattice Boltzmann model for simulating flows with multiple phases and components, Phys. Rev. E 47 (1993) 1815–1819.

[22] X. Shan, G. Doolen, Multicomponent Lattice-Boltzmann model with interparticle interaction, J. Stat.

Phys. 81 (1995) 379–393.

[23] Haibo Huang, Zhitao Li, Shuaishuai Liu and Xi-yun Lu, Shan-and-Chen-type multiphase lattice Boltzmann study of viscous coupling effects for two-phase flow in porous media, Int. J. Numer. Meth.

Fluids 2009; 61:341–354

[24] Irina Ginzburg and Konrad Steiner, A free-surface lattice Boltzmann method for modelling the filling of expanding cavities by Bingham fluids. Phil. Trans. R. Soc. Lond. A (2002) 360, 453–466

[25] Irina Ginzburg, Konrad Steiner, Lattice Boltzmann model for free-surface flow and its application to filling process in casting, Journal of Computational Physics 185 (2003) 61–99

[26] J.M. Buick, J.A. Cosgrove, Numerical simulation of the flow field in the mixing section of a screw extruder by the lattice Boltzmann model, Chemical Engineering Science 61 (2006) 3323 – 3326

[27] J.M. Buick, Lattice Boltzmann simulation of power-law fluid flow in the mixing section ofa single-screw extruder, Chemical EngineeringScience64(2009)52--58

[28] Jonas Latt, Guy Courbebaisse, Bastien Chopard, and Jean Luc Falcone, Lattice Boltzmann Modeling of Injection Moulding Process, Cellular Automata, P.M.A. Sloot, B. Chopard, and A.G. Hoekstra (Eds.): pp.

345–354, 2004.

[29] Y. H. Qian and S. A. Orszag, "Lattice BGK Models for the Navier-Stokes Equation: Nonlinear Deviation in Compressible Regimes," EPL (Europhysics Letters), vol. 21, p. 255, 1993.

[30] D. P. Ziegler, "Boundary conditions for lattice Boltzmann simulations," Journal of Statistical Physics, vol.

71, no. 5-6, pp. 1171-1177, 1993.

[31] Q. Zou and X. He, "On pressure and velocity boundary conditions for the lattice Boltzmann BGK model," Physics of Fluids, vol. 9, no. 6, pp. 1591-1598, 1997.

[32] X. Shan, "Simulation of Rayleigh-Be´nard convection using a lattice Boltzmann method," Phys.

Rev. E, vol. 55, no. 3, pp. 2780-2788, 1997.

[33] S.-C. Chang, Y.-S. Hsu, and C.-L. Chen, "Lattice Boltzmann simulation of fluid flows with fractal geometry an unknown-index algorithm," Chinese Society of Mechanical Engineers, vol.

32, pp. 523-531, 2011.

[34] C. Körner, M. Thies, T. Hofmann, N. Thürey, and U. Rüde, "Lattice Boltzmann Model for Free Surface Flow for Modeling Foaming," Journal of Statistical Physics, vol. 121, pp. 179-196, 2005.

[35] X. Q. Xing, D. L. Butler, and C. Yang, "Lattice Boltzmann-based single-phase method for free surface tracking of droplet motions," International Journal for Numerical Methods in Fluids, vol.

53, pp. 333-351, 2007.

13

[36] R. Ammer, M. Markl, U. Ljungblad, C. Körner, and U. Rüde, "Simulating fast electron beam melting with a parallel thermal free surface lattice Boltzmann method," Computers &

Mathematics with Applications, vol. 67, no. 2, pp. 318-330, 2014.

[37] X. Xiang, Z. Wang, and B. Shi, "Modified lattice Boltzmann scheme for nonlinear convection diffusion equations," Communications in Nonlinear Science and Numerical Simulation, vol. 17, pp. 2415-2425, 2012.

行政院國家科學委員會補助國內專家學者出席國際學術會議報告

2017年9月25 日

報告人姓名 陳介力 服務機構

及職稱 國立成功大學航太系 教授

時間 會議地點

4-5, Sep. 2017 London, UK

本會核定

補助文號 105-2221-E-006 -094 - 會議

名稱

(中文) 英國熱傳研討會

(英文) UK Heat Transfer Conference 2017 發表

論文題目

(中文) 晶格波茲曼法於雙攪拌器之熱傳與質傳效能增進模擬研究 (英文) Lattice Boltzmann Simulation of Enhanced Heat Transfer and 

Mixing Efficiency by Two Active Stirring Impellers 報告內容應

一、 參加會議經過及與會心得

第15屆英國熱傳研討會每兩年舉辦一次是國際重要之熱傳研討會之一。上次(2015) 在愛丁堡大學召開，本年度在Brunel University倫敦召開。本次共有120篇論文分二 天舉行研討。 

研討會主題包含很廣，主要場次大項有Micro  &  Nano  Scale,  Single  Phase  Heat  Transfer, Condensation, Heat Exchangers & Heat Pipes, Heat Transfer in Combustion & 

Thermal Management in Vehicles, Boiling & Evaporation, Computational Heat Transfer,  Heat  Transfer  for  Sustainable  Energy,  Energy  Recovery  &  Heat  Integration,  Heat  Transfer Fundamentals  II  &  Porous  Media, Cooling  of electronics  and  other  high  heat  flux devices, Heat Transfer Enhancement…..。本次發表之論文”  晶格波茲曼法於雙攪 拌器之熱傳與質傳效能增進模擬研究”就屬於熱傳增進之研究領域。 

雖然大會在前一天傍晚有會議reception， 但因前往曼徹斯特與利物浦拜會以往的 教授與同學，本人於9/4日上午8:30到達會場，8:40大會開始，其中曼徹斯特大學 國際著名教授Brian  E.  Launder  也獲邀發表對Prof  D.  Brian  Spalding一生貢獻熱流研 究之回顧。本次在Brunel University主辦，因中國學者在英國任教的學者較多，相對 中國的研究者也是僅次歐美學者的大團體。日本也只有4篇，台灣只有三篇文章參 與本次研討會。本人出席的場次中以Boiling  &  Evaporation,  Condensation,  Heat  Exchangers & Heat Pipes, Heat Transfer Enhancement 收穫最多。藉由本次研討與人 際交流對將來之研究主題之開發與成果展現均有非常大的助益。

Keynote speech  其中一場為新熱系統應用與設計，對強健設計的建構與可能新主題 都提出不錯的研究參考。 

  二、攜回資料名稱及內容

本次研討會攜回資料包含大會論文紙本與電子檔(記憶碟) 各一份。