可視化水平區融法生長非線性光學晶體之應用研究(1/2)

行政院國家科學委員會專題研究計畫 期中進度報告

可視化水平區融法生長非線性光學晶體之應用研究(1/2)

計畫類別： 個別型計畫 計畫編號： NSC92-2214-E-002-038- 執行期間： 92 年 08 月 01 日至 93 年 07 月 31 日 執行單位： 國立臺灣大學化學工程學系暨研究所 計畫主持人： 藍崇文 報告類型： 精簡報告 處理方式： 本計畫可公開查詢 中 華 民 國 93 年 5 月 28 日

92年國科會研究期中報告書

可視化水平區融法生長非線性光學晶體之應用研究

執行單位：國立台灣大學化工系

計畫編號：NSC 92-2214-E-002 -038 –

執行時間：92年8月1日至93年7月31日

計畫主持人：藍崇文 教授

目錄

中文摘要 1 英文摘要 2 第一章 緒論 3 第二章 系統架設及研究方法 8 2-1 系統架設 8 2-2 研究方法-系統架設 .13 2-3 研究方法-模擬方面 .18 第三章 實驗結果與討論 ..34 3-1 長晶系統之測試 34 3-2 模擬部分 37 第四章 結論及建議 39 4-1 系統之改進與建議 .39 4-2 預計進行之研究 .40 參考文獻 42

中文摘要

高品質單晶的生長及摻雜控制乃至於製成的最佳化，對於非線性 光學晶體在工業上應用扮有極為重要的角色。然而在自然融熔態生長 的晶體，或多或少都會偏離化學計量組成，形成對位缺陷。此對位缺 陷對晶體的性質，尤其在元件的加工及應用上，有極大的不良影響。 此外，在摻雜的應用上，更因為摻雜離子在固體及液體中的溶解度不 同，造成偏析及晶體的不均勻性。因此本計畫將使用可視化的水平區 熔系統，以等計量比鈮酸鋰單晶作為生長對象，來進行摻雜均勻的等 計量比晶體生長，以期望能得到低耗能且摻雜均勻的等計量比單晶體 生長技術，並輔以電腦模擬計算，發展雙層的系統熱流模式，並分析 最佳化晶體生長程序，以期開發高效率低耗能的等計量且均勻性高的 單晶技術。

Abstract

Growing of high quality single crystal, controlling of dopant segregation, and optimizing the growing process play an important role in industrial application of nonlinear optical crystals. However, crystals growing from the melting ones always distract from its stoichiometric composition and have anti-site defect. This anti-site defect has bad effect on the properties of crystals, especially on its application. In addition, because the dopant has different solubility in solid and liquid, it will make segregation and non-uniform of crystals when doping. Through the process by using visual horizontal zone melting method to get stoichiometric lithium niobate single crystal, we expect to get stoichiometric single crystal growing technique with low energy cost, high efficiency, and uniform dopant concentration. Besides doing the experiment, we also develop two-level numerical method to simulate the thermal and flow field of the crystal growth and analyze the optimum crystal growing process. We expect to get the single crystal growing technique with low energy cost, high efficiency, and uniform dopant concentration.

第一章 緒論

非線性光學單晶在光電元件中是屬於重要之零組件，而高品質的 非線性光學單晶之生長及製備技術更是光電產業中原料供給的關鍵 技術。光折變晶體可應用於相位共軛光學元件、光學記憶元件、圖形 辨識元件等光電元件【1-3】，近年來，由於光資訊處理、計算、影像 儲存等應用的快速發展及需求，光折變晶體之開發及生長技術更是倍 受重視。 晶體的生長方式有柴氏法（Czochralski Method，簡稱 Cz 法）

【4-10】、頂晶熔鹽法（Top Seed Solution Growth，簡稱 TSSG 法） 【11】、布氏法（Bridgman Method）【12】和懸浮熔區法（Floating Zone Method，簡稱 FZ 法）【13】。一般而言，使用 Cz 法生長晶體可得到 較大直徑的晶體，但 Cz 法之生長需較長的時間，且常伴隨著晶體內 雜質濃度分布不均的問題。至於 TSSG 生長雖然應用於許多塊晶的生 長上，但其生長速度較 Cz 法更慢，且存在著晶體直徑控制不易的問 題，在工業上要運用於大量生產較為困難。而布氏法可分為垂直布氏 法（Vertical Bridgman Method，簡稱 VB 法）及水平布氏法（Horizontal Bridgman Method，簡稱 HB 法），以往布氏法大多用於 II-VI 及 III-V 族半導體之生長，由於其操作簡單及成本低廉，因此近來使用十分普

遍，其中 VB 法之低溫度梯度更大大降低生長晶體時因熱應力而產生 的缺陷，不過其晶體在軸向偏析的問題仍待解決。FZ 法雖然無法在 生長直徑上與 Cz 法及 VB 法比擬，而其在生長過程中具有之大溫度 梯度，亦使得晶體受熱應力之影響產生較多缺陷，且其不適用於工業 大量生產，但期能以較高之長速生長，且軸向雜質分布之控制較容 易，因此在複雜摻雜及高計量控制之晶體製備較為容易。相似於 FZ 法，還有一種長晶法稱為區熔法（Zone Refining Method），其於 FZ 法之差別在於區熔法有坩堝壁可支撐熔區而 FZ 法之熔區是靠熔液之 表面張力支撐，其在生長晶體過程之特性皆相似於 FZ 法。 應用於工業上高品質的非線性光學晶體必須考慮到兩項最重要 的條件，分別是晶體品質及晶體中之摻雜分布，然而在現今之生長技 術中，晶體之化學計量及摻雜分布不均卻是大塊晶體普遍存在的問 題。於是本計畫擬以一非傳統的生長技術，將區熔（zone refining） 及其連續的製程的概念，來生長目前相當熱門的非線性光學等計量比 單晶 LiNbO3（SLN），稱為水平式區熔法（Horizontal Zone-melting

Method，簡稱 HZM 法）。傳統的單晶生長幾乎為一批次程序（batch process），若用於生長晶體的原始材料不是單一成分，例如 LiNbO3

摻雜 MgO 或 Fe；或是材料可以非化學計量（non-stoichiometry）存 在，在生長過程中，由於偏析（segregation）的影響，摻雜於晶體中

的濃度分佈或晶體的組成沿晶體之分佈將會不同，這對於晶體生長及 其應用上將增加許多限制及增加材料的浪費。我們希望將我們以往在 區熔法上的成果，運用在水平式的單晶生長，將使用可視化的水平區 熔法此連續式（continuous）操作，以等化學計量比的鈮酸鋰單晶為 對象，來達到晶體等化學計量比及均勻摻雜，以及對晶體生長過程中 做晶體接種的觀察，並且可以得到直徑均勻的單晶。另外，在熔區加 熱的方法，本實驗室採用一種相當新穎的方法，稱為浸入式加熱法 【14】，此方法的特色是可以減少加熱材料所需的能量，以達到利用 最少的能量得到所需的熔區。在本計畫中，我們將使用水平式區熔法 配合浸入式加熱法來生長晶體，以達到晶體中摻雜均勻分佈及減少能 量的使用，並進行單晶生長的可視化研究。 另外，將對水平式區熔法進行電腦模擬之研究，主要是模擬長晶 系統整體的流場與溫場。對於生長高品質單晶而言，熱場的設計是很 重要的，譬如：尋找適當的熱場分佈、保溫材的設置與瞭解熔液流動 與熱傳之間的相互影響等。本研究模擬的方法主要建立雙層模式。關 於上層模式部分，有一些發展很完全的商業軟體，例如Fluent【32】、 CFDRC【34】、ANSYS【33】、FIDAP【31】與STHAMAS【29】等， 可以幫我們處理很複雜的幾何形狀，也可以幫我們計算三維的熱輻 射，但是唯一的缺點是這些商業軟體在處理自由凝固界面的問題能力

比較弱。文獻中報導的整體模式大多是模擬柴氏法為主【19-20】，其 他相關的研究如柴氏法生長矽單晶熱場設計最佳化等【21-23】，均是 使用二維軸對稱系統，只有少數利用超級電腦模擬三維柴氏法的整體 模式，如Derby團隊【24-25】使用三維整體模式模擬整體長晶系統， Lan【26】與Müller【27】有做相關的研究回顧。一般在做整體模式 的研究時，大多選用商業軟體模擬複雜的長晶系統，相關的研究包括 Evstratov【28】等人模擬柴氏法生長單晶，使用商業軟體CFD-ACE+ 計算包括整體熱傳、氣體流動、熱輻射與界面形狀；Müller【29】利 用他自己實驗室發展出來的軟體-STHAMAS，模擬晶體生長程序的熱 場與熱場設計等。Lukanin【30】等人結合二維與三維模式模擬柴氏 法長晶，局部模式為三維系統，包括矽熔液、坩堝、晶體與所延伸包 含的氣體；整體模式為二維系統，所得到的熱場可帶入三維模式當作 邊界條件。他們計算所的計算商業軟體為CGSim。以上均為使用商業 軟體計算整體模式的例子。我們的上層模式為長晶系統的熱流模式， 我們選用商用CFD軟體Gambit/Fluent來處理。此層將包含三維的熱輻 射，氣體對流，及各部份材料的熱傳導的整體計算，計算所得到的表 面熱通量可帶入下層模式當作邊界條件；而下層模式需處理自由凝固 界面，則是以我們實驗室建立的多重網格有限體積法（MG/FVM）的 熱流模式為主，此部份將處理晶體凝固時的界面變形，液相熱流，與

氣液界面熱毛細流的問題。此三維長晶系統之熱流模擬程序的建立可 以幫助我們模擬整體長晶的溫場與熱場，並分析最佳晶體生長程序， 將來可以應用於其他長晶系統三維熱流模擬的建立。

第二章 系統架設及研究方法

2-1 系統架設

2-1-1 水平式區熔系統 ( Horizontal Zone-melting method )

在本計畫中，我們所採用的水平式區熔法(HZM)，實驗設備如圖 一所示，由圖二的概念圖可知，此系統為利用一個半圓柱狀的的白金 舟（圖三），其大小為直徑 2 公分長 15 公分厚度 0.03 公分，裡面填 充晶體生長所需材料及適量之摻雜，將白金舟置於系統中，在熔區部 份加入一長為 8 公分、寬 0.3 公分、厚 0.05 公分的白金片（圖四）， 使白金片與材料直接接觸，之後給予一約 80 安培的高電流通過白金 片，使白金片產生約 1200℃高溫，讓白金舟內產生一個穩定之熔區。 一開始所得的熔區與船中的子晶相接處，之後慢慢移動白金舟，向保 溫區靠近，可得與白金船相同直徑大小的晶體。由於晶體之生長過程 中，會使得熔區中之熔液之組成不斷改變，影響了生長出來晶體之品 質，再加上由於摻雜在固液相中溶解度的不同，偏析現象造成晶體生 長過程中軸向的雜質分佈不均。因此，在生長過程中，浸入式加熱器 固定不動，白金舟持續往保溫區移動，相當於以連續的方式進料，以 替代晶體生長時熔區中所失去的熔液，並且能使熔液中之雜質能保持 在一定的濃度，達到生長出軸向雜質分佈均勻之晶體。此外，區熔段

也可事先摻入較多的摻雜，以減短達到 Zone melting 的時間。 目前我們已經完成長晶系統中的保溫區段與預熱區段，而熔區加 熱部分，則尚在測試，以找出能有效達到加熱效果的白金片大小與所 需的電流。 圖一 水平式區熔法(HZM)儀器圖 圖二 水平式區熔法概念圖

圖三(a) 白金舟全長 15 公分、直徑 2 公分 圖四 浸入式加熱白金片，大小為長 8 公分、寬 0.3 公分、 厚 0.05 公分 2-1-2 預熱裝置(Pre-heater) 由於使用潛入式加熱器加熱所產生的功率較小，因此其加熱範圍 也甚小。為了減輕白金片所產生的功率以及有效的控制熔區範圍，避 免白金片因過負載而斷裂，因此在進入生長區之前，先將進料預熱， 一方面可以減輕白金片所需產生的功率；另一方面也可以減少生長區

的熱散失。圖五為本系統所使用的預熱裝置。此裝置由本實驗室自行 組裝，此構造為一個氧化鋁管，其長度約為 38 公分，管外纏繞加熱 絲，並以高溫膠固定，作為加熱氧化鋁管之用，而管外包覆 6 層高溫 棉，做為絕熱用。此裝置外接一變壓器(Sin Chyuan Hsin Co.)，可輸 入的電壓為 110V 與 220V 兩種，輸出的電壓範圍 0-260V。此預熱裝 置的溫度由輸出的電壓所控制。

2-1-3 熔區加熱裝置

對於熔區加熱的部分，我們採用目前相當新穎且省能源的方式， 此方法稱潛入式加熱法(Immersed heating method)，此方法的特色是利 用白金片做為加熱體，將白金片與要加熱的材料直接接觸，利用電流 通過金屬會因為金屬有電阻的關係，將電能轉換成熱能，而達到將材 料熔化的目的。由於白金片與材料直接接觸，因此白金片所產生的熱 能可以完全的被材料所吸收，所散失到空氣的熱能極為少，可以減少 加熱材料所需的能量。其略圖如圖七所示，此加熱器由白金片與一變 壓器所組成，此變壓器的輸入電壓為 220V；電流為 4A，可轉換輸出 電壓為 7V；電流為 120A。由於高電流通過白金片後會產生很高的溫 度，為了連接白金片與變壓器，我們在兩者之間加入一長為 60 公分、 寬 5 公分、厚 5 公厘的銅片做為之間的連接器，避免高溫的白金片將

變壓器的接頭熔化。 2-1.4 保溫裝置 由於在晶體生長之後，若降溫太快，晶體以因為熱應力的關係容 易產生裂痕(crack)，因此在晶體生長之後的保溫動作極為重要，在本 系統中，我們採用 0.5mm 的鎳鉻絲作為加熱體，將加熱絲纏繞成螺 旋狀，固定在氧化鋁磚中，而管外包覆 6 層高溫棉，做為絕熱用，如 圖六所示。在溫度量測與控制方面，我們採用 K-type 的熱電偶，與 BTC-9090 溫控器 圖五 預熱裝置側視圖 圖六 保溫裝置側視圖

圖七 熔區加熱裝置略圖 2-2 研究方法-實驗方面 2-2-1 熔區材料之製備 對於鈮酸鋰晶體，我們將生長化學計量(stoichiometric)之鈮酸鋰晶 體。我們使用 99.99％氧化鈮(Nb2O5)、99.99％之碳酸鋰(Li2CO3)及 99.9 ％氧化鎂(MgO)為原料。為了生長摻雜氧化鎂的化學計量組成之鈮酸 鋰，使用一般的柴式拉晶法(Czochralski method)是無法達到的，必須 使用連續式操作(continuous process)才有可能達到。本計畫中所提出 之水平式區熔法即為一連續操作程序，針對熔區部分，我們將配製富 含鋰(Li-rich，Li2O 58-60 mol％)的鈮酸鋰摻雜 1mol％Mg 為原料；由

於在利用浸入式加熱器加熱時熔區上方的粉末會有粉末硬化(soaking) 的現象發生，因此我們使用含 1mol％Mg 的共熔組成鈮酸鋰晶體

（CLN）、Li2CO3、 Nb2O5以及 MgO 粉末混和成為富含鋰(Li-rich， Li2O 60 mol％)的鈮酸鋰摻雜 1mol％Mg 做為熔區的部分。 2-2-2 進料鈮酸鋰的配製 針對連續進料部分，我們將配製化學計量之鈮酸鋰粉末(Li:Nb=1:1) 做為原料。由相圖（圖五）可知，當鈮酸鋰熔液組成為 Li2O 58-60 mol ％時，其最初凝固所得晶體組成即為化學計量組成，若我們於此時補 充同為化學計量組成的鈮酸鋰進料，則熔液中之組成成分便維持不 變，於是我們可生長出組成維持均一化學計量組成的晶體。由於區熔 法其熔區必須維持一固定大小才能得到均勻組成的晶體，因此進料的 體積也必須維持固定，並且結構必須非常緻密，避免進料體積不均勻 影響晶體的生長，因此本實驗室才採用緻密的燒結棒作為進料棒，其 製作流程如下: (1) 將碳酸鋰(Li2CO3)、氧化鈮(Nb2O5)以及氧化鎂(MgO)粉末至入烘箱 中，烘箱溫度維持 100℃24 小時，使粉末乾燥。 (2) 將乾燥後的粉末，以莫爾比，Li2CO3：Nb2O5：MgO＝49.5：49.5： 1 的比例配製，並將配製好的粉末均勻混和。 (3) 把混和好的粉末放入烘箱中，以 1000℃的溫度維持 24 小時做粉末 的燒結，即可得到摻雜 1mol％MgO 的等計量比鈮酸鋰粉末(SLN)。

約 20 公分，直徑 2.2 公分的半圓柱，如圖九所示。 (5) 將裝填好的粉料棒以 700kg/cm2的壓力維持 1 小時將其壓錠，可得 長約 18 公分直徑 1.8 公分的半圓柱，如圖十所示。 (6) 將壓錠好的進料棒放入高溫爐中，以 1150℃的溫度維持 4 小時， 使的進料棒更為緻密，其大小為長 18 公分直徑 1.7 公分，緻密度 約為 0.94％的燒結棒，如圖十一(a)、(b)所示 圖八 鈮酸鋰相圖

圖九 未壓錠前的進料棒 圖十 壓錠後的進料棒 圖十一(a) 燒結後的進料棒 圖十一(b)燒結後的進料棒 2-2-3 晶體生長 當儀器設備架設完成後，必須先針對晶體生長所需的溫度範圍量 測系統中的溫度分布，並藉由系統之溫場量測，針對保溫部分及預熱 部分做適當之設計。 將子晶、熔區粉末及進料棒置入白金舟中，其中子晶大小為長 1.5 公分、直徑 1.8 公分，熔區約 1 公分，進料棒約 12.5 公分(如圖十 二(a)、(b))，再將白金舟置於生長系統中適當的位置後，我們將白金

片的電流，可得到一與白金舟半徑差不多大小的熔區，經過 2~3 小時 後，確定熔區大小穩定，便將熔區慢慢向子晶相接觸，待子晶晶種回 熔並穩定，觀察於此時並可觀察子晶熔化情形以判斷熔區大小與溫度 是否適當。 當系統穩定後，開始以設定之速度生長晶體（針對個別需要有不 同需求之長速）。晶體生長開始的同時，白金舟也慢慢向保溫區移動， 減少生長晶體脫離生長區域所產生的溫度梯度，以減少熱應力對晶體 所產生的影響，晶體生長完成後，將白金舟移動到保溫區高溫爐的中 心。將保溫區高溫爐的溫度慢慢下降至室溫後，方可將晶體取出並對 晶體進行其他後處理工作。 圖十二(a) 白金舟內粉末的擺設 圖十二(b) 長 1.5 公分直徑 1.8 公分 的子晶<001>

2-3 研究方法-模擬方面 在模擬方面為晶體生長系統之總體熱流研究，主要分為兩部分， 第一部份為整體模式，包括計算氣體流動、表面熱輻射與整體各材料 的熱傳計算，第二部份為處理坩堝與固液的熱傳與流動，包括固液生 長界面等。首先我們要先瞭解如何進行雙層模式的模擬，接下來再分 別介紹雙層模式所使用的主導方程式與邊界條件及其原理。 模擬步驟 圖十三為雙層模式的模擬步驟，首先我們必須先產生網格(Grid Generation)，可用的軟體包括 Gambit 與 TGrid 等網格產生軟體，可 建立複雜的長晶系統形狀，依據邊界特性建立網格，並且設定邊界型 態與定義流體的封閉區間與固體區塊，我們定義熔液為固體，也就是 先忽略熔液流動與界面變形的影響。所產生的網格可用計算流體力學 軟體(CFD)讀入，進入計算流體力學軟體之後，必須先定義系統尺寸 大小，然後選擇系統的物理模式與求解方法、操作條件、邊界條件後 與收斂條件後，計算部份分為兩個步驟，第一部份為計算觀察因子 (view factor)，第二部分開始計算氣體流動、表面熱輻射與整體各材 料熱傳之計算。皆下來我們要做後處理的部分(Postprocessing)，首先 我們將計算的結果輸入至 Tecplot，同時也將局部模式的網格輸入至 Tecplot，利用 Tecplot 將 Fluent 計算的熱通量內差成與我們系統的網

格數據，再將內差的結果輸入至局部模式，當作邊界條件帶入本實驗 室發展的程式，此部分計算的系統包括坩堝、進料、晶體與熔液部分， 這部分的計算包括熔液流動、各材料的間熱傳與固液生長界面的變形 模擬等。我們再將所變形的網格與局部模式所計算的表面溫度帶入整 體模式計算，如此進行上述步驟，直到局部模式與整體模式所得到的 溫場收斂為止。接下來我們會分別介紹整體模式與局部模式如何分別 運作。

圖十三 雙層模擬步驟流程圖

Local model

(in-house code) ¾ Melt flow ¾ Heat transfer

¾ Growth interface shape

Grid Generation

(Gambit, TGrid…)

¾ Set up complex geometry ¾ Mesh the model

¾ Set up the boundary type

CFD

(Fluent, CFDRC, ANSYS, CFX…) ¾ View factor calculation

¾ Gas flow ¾ S2S radiation

¾ Overall heat transfer

Postprocessing

(Tecplot, Plot3D…) ¾ Process results ¾ Linear interpolate ¾ Data output

Temperature _{Heat fluxes}

Local model

2-3-1 整體模式 對於整體晶體生長系統之熱流模擬，我們使用商業的計算流體力 學軟體 Fluent，此套軟體搭配建網格的 Gambit 軟體使用，它可以模 擬複雜幾何形狀的熱場與流場。首先我們要先建構網格，Fluent 可以 接受結構與非結構性網格，所以他可以適用各種複雜的形狀，可接受 網格的型態包括二維的三角形與四邊形網格，三維的網格包括四面 體、六面體、角錐形等網格。此外，Fluent 可以針對流動的狀況在某 些特定的邊界上加入邊界條件，所以可以適應的各種複雜的幾何形 狀。當網格讀入之後，首先必須要定義系統的尺寸，接下來設定系統 的物理模式、邊界條件、物性參數與操作條件之後，接下來計算分兩 部分，第一部份為根據系統幾何形狀計算視因子，第二部分計算包括 氣體流動、三維熱輻射、熱對流與總體材料熱傳導。Fluent 是以 c 語 言所寫的，所使用的離散方法為有限體積法（Finite volume）。 2-3-1.1 整體模式-主導方程式 針對此長晶系統，我們以笛卡兒座標（Cartesian coordinate）(x,y,z) 的方式來描述。對於流動、熱傳的主導方方程式描述如下：

( )

=0 ⋅ ∇ + ∂ ∂ υ g g t ρ ρ (1)

( )

g

(

g

)

P

( )

τ geG t ρ +∇⋅ ρ =−∇ +∇⋅ +ρ ∂ ∂ υυ υ (2)

τ：stress tensor，τ = _⎢⎣⎡

(

∇υ+∇υ

)

− ∇⋅υI_⎥⎦⎤ 3 2 T µ I：unit tensor T T t T i 2 ∇ = ∇ ⋅ + ∂ ∂ _α υ (3) 其中t、υ、P、T為時間、速度、壓力、溫度，μ為氣體的黏度， ρg為 氣體的密度，αi為系統內各材料的熱擴散係數，流動以層流為主，我 們對氣體的假設為理想氣體(ideal gas)， T M R P w op g × = ρ ，Pop為操作壓 力，R為氣體常數，Mw為氣體的分子量，除了氣體密度之外，所有物 種的物性都假設為常數。物質性質列於表一。 2-3-1.2 觀察因數(view factor)的計算與熱輻射 本 研 究 所 使 用 的 熱 輻 射 方 法 是 面 對 面 輻 射 (surface-to-surface radiation)，對於模擬封閉空間的輻射是一個相當好的方法。此方法假 設空間中的介質是透明的，不會參與輻射(例如：太空船熱排放系統， 太陽能收集器與摩托車的冷卻系統)，雖然此方法在計算觀察因數非 常耗時，尤其是在三維的系統或是需要參與輻射的面比較多的系統； 但是相對於其他假設介質會參與輻射的方法(如 Discrete Ordinate 與 Discrete Transfer)，此種方法在疊代時可以節省相當多的計算時間， 即使觀察因數的計算相當耗時，但是只要在一開始的時候計算一次之 後，之後程式所做的計算都可以用第一次計算的觀察因數，而不必再 重新計算，除非計算的過程中更改邊界型態的設定。

收 率 (α， absorptivity) 、 反 射 率 (ρ， reflectivity) 與 穿 透 率 (τ， transmissivity)。這些性質會隨著入射角度，放射角度(θ,φ)與波長大 小(λ)改變，通常這些數據必須由實驗而來，或是藉由古典電磁波理 論(classical electromagnetic wave theory)計算得到。在 Fluent 的面對面 輻射方法假設所有的面都是擴散-灰表面 diffuse-gray。擴散(diffuse)的 意思即放射率與吸收率不因方向而改變，因此對於放射而言，黑體的 放射強度在每一個方向都是相同的。灰表面(gray)的意思對於某一表 面其放射率與吸收率不隨波長的改變而改變，吾人稱此為灰表面。所 以當面假設為擴散-灰表面時，意旨在固定溫度下，表面在任一方向 與任一角度，吸收與放射的熱輻射是固定的。 α(λ,θ,ϕ,TA)=α(TA) ε(λ,θ,ϕ,TA)=ε(TA) (4) 根據 gray-body 模式，表面輻射能 E，一部分會反射(ρE)，一部 分會吸收(αE)，另外一部分會透射(τE)。對於大部分的問題來說，表 面是不透明的，因此我們可以忽略透射的影響，由能量守衡來說，ρ +α=1。根據 Kirchoff’s Law【15】的原理，α=ε，因此ρ=1-ε。 離開表面的熱通量包含放射與反射，反射的熱通量與周圍環境所 給的熱通量有關，即可表示成離開其他所有表面的熱通量。所以由k 面的反射熱量即可表示成

qoutk k Tk kqin,k

4

. =ε σ +ρ (5)

其中 qout,k 為離開 k 面的熱通量，εk為放射率，σ 為波茲曼常數

(Boltzmann constant)，qin,k為周圍環境給k表面的熱通量。

觀察因數，Fjk表示有多少部分的能量由k面到 j面。所以熱通量 qin,k可以表示成與離開其他面熱通量有關的關係式：

∑

= = N j jk j out j out k Aq F q 1 , , (6) 其中Ak代表k面的表面積，Fjk表示k面與j面的觀察因數。對於N個 表面來說，視因子的關係式如下式： AjFjk=AkFkj for j=1,2,3…N (7) 所以

∑

= = N j j out kj k in F q q 1 , , (8) 因此

∑

= + = N j j out kj k k k k out T F q q 1 , 4 , ε σ ρ (9) 上式可以被寫成

∑

= + = N j j kj k k k E F J J 1 ρ (10) 其中Jk表示由k表面輻射熱，Ek表示 k表面的放射能(emissive power)， 因此有N條方程式，所以可以表示成N×N的矩陣： KJ=E (11)

vector。介於 j與k表面的觀察因數為： =

∫ ∫

k J A A kj k j j k k kj dAdA r A F δ π θ θ 2 cos cos 1 (12) 其中δkj定義為 dAj面對 dAk面看見的程度。面對面輻射對於有很多輻 射表面的例子是一個非常耗費計算資源的方法，為了要節省計算的時 間，我們必須將一些面聚集（clusters）起來，由一個面開始，加入其 他鄰界的面合在一起。對於計算觀察因數方面，Fluent 提供兩個方 法，第一個方法為”Adaptive”【16】，另外一個方法為”Hemicube” 【17】。Adaptive 的方法是利用多樣化的演算法(analytic or Gauss quadrature)計算觀察因數，Adaptive 的方法適合用於 2 維或是幾何形 狀 較 不 複 雜 的 系 統 ， Hemicube 的 方 法 使 用 微 分 面 對 面 的 方 法 (differential area-to-area method )計算觀察因數，而 Hemicube 較適用 於 3 維或是幾何形狀較複雜的系統。所以在本研究中，計算視因子均 使用 Hemicube 的方法。 所以在計算整體模式的時候，除了動量方程式與能量方程式之 外，另外還必須計算輻射非線性方程式。 2-3-1.3 邊界條件 對於整體模式而言，關於速度的邊界條件，在固體的邊界上以不

滑動（non-slip）的邊界條件設定，因為本系統固體都是靜止的，所 以氣體接觸固體部分的速度均設為零。另外，我們假設此次系統是封 閉的，所以氣體兩邊進出口的速度為零。 對於熱傳的邊界條件，在爐子外部設定一高斯溫度分佈 （Gaussian distribution）Ta(x)作為描述： ( ) ( )exp[ ( )2] a x x T T T x T_a = _L+ _H − _L − − p (13) 其中TH為最高溫，TL為最低溫，xp為最高溫的位置，a為溫度分佈的 寬度。固體與固體間為共軛熱傳，對於某一界面來說，進入的熱通量 等於出來的熱通量，即 qin =qout (14) 固體與氣體間的熱傳條件

∑

= + + − = ∇ ⋅ − N j j out kj k k k k k T hT T T F q k 1 , 4 ) ( ε σ ρ n (15) 其中，n為 k表面上朝外的單位長度向量，h為熱傳係數。 2-3-2.1 局部模式 – 主導方程式 對於水平區熔法的生長系統的局部模式，包括坩堝(a)、晶體(c)、 進料(f)與熔液(m)。我們以笛卡兒座標（Cartesian coordinate）(x,y,z) 的方式來描述。對於熔液流動、熱傳的主導方方程式描述如下： 0 = ⋅ ∇ υ (16)

( )

(

)

(

m

)

G T T T p t + ⋅∇ =− ∇ + ∇⋅ ∇ + ∇ − − e ∂ ∂ _{ν ρβ} ρ υ υ υ υ υ 1 (17)

T T t T m 2 ∇ = ∇ ⋅ + ∂ ∂ _α υ (18) 其中t、υ、P、T為時間、速度、壓力、溫度，ν為熔液的動黏滯係數 (kinematic viscosity)，ρm為熔液的密度，αm為熔液的熱擴散係數，eG 為重力加速度，流動以層流為主，我們對液體的假設為 Boussinesq model，β為熱膨脹係數。在動量方程式中以熔液的熔點溫度 Tm為參 考溫度。 在固體的坩堝、晶體與進料部份，僅需要考慮熱傳方程式： U T 2T (i c,f,a) t T i x amp ⋅∇ = ∇ = + ∂ ∂ _α e (19) 其中，αi分別為晶體(c)、進料(f)與坩堝(a)之熱擴散係數。 2-3-2.2 邊界條件 關於速度的邊界條件，在固體邊界上是以不滑動(non-slip)的邊界 條件設定： cUampex r γ = υ (20) 其中，γc =ρc ρm。對於自由面上的速度邊界，即熔液與氣體接觸的邊 界條件是以剪切應力(shear stress)平衡： s T T ∂ ∂ ∂ ∂ = ⋅ ⋅τ s γ n (21) 其中，s與 n為在熔液與氣體接觸表面的切線方向與法線方向的單位 向 量 ； ∂γ ∂T 為 表 面 張 力 溫 度 係 數 (surface tension-temperature

對於熱傳的邊界條件，是以熔液/固體的熱通量相等作為設定， 如下所示： ⎟ ⋅ =0 ⎠ ⎞ ⎜ ⎝ ⎛ ₊ ∆ + − c e_x n amp c m dt dh U H q q (22) 0 = ⋅ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ₊ ∆ + − c e_x n amp f m dt dh U H q q (23) 其中，n為進料或生長界面朝向熔液的單位長度向量。qm，qf，qc為 在熔液、晶體與進料固體的熱通量，ΔH為熔化熱。 在做整體模式與局部模式結合之前，我們必須做一些測試的例子 來證明我們的結果可靠性，在測試的例子中，坩堝、晶體、進料與熔 液與加熱環境的熱交換主要以輻射與對流的方式，所以在設定坩堝、 進料、晶體與熔液表面的邊界條件上是以上述材料表面的熱通量與坩 堝的傳導熱通量相等作為設定，如式(24)(a)所示。而在與整體模式結 合方面，進料、晶體與熔液表面的熱通量就等於整體模式所計算在固 體邊的熱通量，如式(24)(b)所示： m c f a i q T T T T h T ki _i ( ) i( ) Global , , , 4 4− + = + − = ∇ ⋅ −n _∞ σε _∞ (24) (a) (b) 其中 n 為各材料表面朝外的單位長度向量，εi為各材料的表面放射 率，h為熱傳係數。在測試的例子中，外圍溫度T_∞是以高斯溫度分佈， （Gaussian distribution）T∞(x)作為描述： ( ) ( )exp[ ( )2] a x x T T T x T L H L p − − − + = _{∞ ∞ ∞} ∞ (25)

的寬度。 2-3-3 無因次化方程式 對於一個物理問題，我們經常以之前所述的主導方程式與邊界條 件轉換成無因次的方式來表示。這部份我們僅針對局部模式來作為探 討的系統。定義特性長度為直徑Dm，特性速度為αm/Dm，特性壓力為 2 2 / _m m mα D ρ ，特性時間Dm/αm 2 _{，無因次溫度 T}*_{以 T/T} m來定義，所以主導 方程式與邊界條件可以寫成： 在熔液部份： 0 * *⋅ = ∇ υ (26)

( )

(

)

G T RaT Pr P t e * * * * * * * * * * * * * + ∇ + ∇ ⋅ ∇ + −∇ = ⋅ ∇ ⋅ + ∂ ∂ υ υ υ υ υ _Pr (27) * * * *2 * * * T T t T ∇ = ∇ ⋅ + ∂ ∂ υ (28) 其中，Pr 為普蘭特數(Prandt number)，其定義為Pr=ν_m α_m ；Ra 為 Rayleigh number，其定義為 m m m TD g Ra ν α β∆ 3 = 。 在固相部份， T i c f a t T m i , , * 2 * * * = ∇ = ∂ ∂ α α (29) 對於自由面的剪切應力平衡的邊界條件，其無因次式為 s T Ma ∂∂ = ⋅ ⋅ * * s τ n (30) 其中，Ma為 Marangoni number，定義

(

)

m m m m TD T Ma α ν ρ γ ∂ ∆ ∂ = 。在固液界面 的熱通量平衡邊界條件的無因次式為：

0 * * * ⋅ = ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ + + − ex n m amp c m dt dh U St Q Q (31) 0 * * * ⋅ = ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ + + − ex n m amp f m dt dh U St Q Q (32) 其中Qm、Qc與Qf為在熔液、晶體與進料中的無因次熱通量，St為史 蒂芬常數(Stefan number)，St 其定義為 T Cp H m∆ ∆ 在系統表面與外部環境的熱交換的無因次式為 ki T_i Bi(T Ta) Radi(T Ta ) q i m, f,c,a 4 * 4 * * * *∇ = − + − + = ⋅ − ∗ Global n (33) 其中， * i k 為無因次的熱傳導係數，其定義為 m i i k k k*= 。Bi為 Biot number， 定 義 為 m m k hD

Bi= 。 Radi 為 輻 射 數 (Radiation number) ， 定 義 為

m m m i i k D T Rad 3 σε = ，εi為各材料的放射率。 * Global q 為由整體模式所得之無因 次表面熱通量。 2-3-4 數值方法 先前所述的主導方程式與相關連的邊界條件，在整體模式方面， Fluent 是以有限體積法的方式作方程式離散，我們選擇單重網格方 法 ； 局 部 模 式 也 是 利 用 有 限 體 積 法 作 離 散 ， 並 以 多 重 網 格 方 法 (multigrid method，簡稱 MG)加速計算效率。關於詳細的數值方法可 以參考 Liang[15]的論文。

2-3-5 模擬測試系統之描述 本次進度報告模擬部分主要以建立模擬程序為主，我們模擬所測 試的系統如圖十四。圖十四(a)整體模式方面包括石英外管、保溫材、 空氣、坩堝、晶體與熔液部分，圖十四(b)為圖十四(a)中 x=L/2 的截 面圖，其中 sample 包括熔液、晶體與進料。圖十四(c)局部模式主要 是包括坩堝、晶體與熔液部分。圖十五(a)與(b)為整體模式網格圖。 圖十五(c)與(d)為局部模式網格圖。

圖十四(a) 整體模式示意圖(c 晶體，m 熔液，f 進料，a 坩堝)

圖十四(b) 整體模式剖面圖(x=L/2)

圖十五(a)整體模式之截面網格圖(z=0)

圖十五(b)整體模式之截面網格圖(x=L/2)

第三章 結果及討論

3-1 長晶系統之測試 3-1-1 溫場量測 我們首先將熔區加熱，並藉由溫場的量測來了解所架設的系統對 晶體生長的可行程度。我們將預熱裝置通予電壓 110V，溫度約為 850 ℃，保溫區設定在 800℃，生長區的白金片通予 80A 的電流，其溫度 約為 1200℃，在把熱電偶固定在傳送系統上，在距離熔區上方 1 公 分的位置，以 1cm/min 的速度移動，量測每移動一公分所得到的溫 度，圖十六為所量測到系統內的溫度分佈。

3-1-2 其他因素之討論 在我們所使用的水平式區熔法中，測試其生長晶體時之儀器操作 上，發現在操作上有些許的設計需加以改進： (1) 在水平區熔中我們使用浸入式加熱法作為熔區的加熱裝置，但是 在使用白金片作為加熱體的過程中，白金片容易因為在加熱過程 中，因為高溫的原因，造成結構上的改變，因此白金片在通過高 電流達到熔區熔化的溫度後會有斷裂的情況發生造成熔區熱場的 不穩定，因此在白金片尺寸的選擇上需要在做改進，應選擇較厚 及較寬的白金片以增強白金片的強度及其加熱的區域。 (2) 在白金片與變壓器之間的連結銅片，在長期受到高溫的影響下，， 使得銅片溫度約 250℃，會嚴重產生氧化銅造成絕緣的作用(圖十 七（a）)，使得白金與銅片間的的導電並不順暢，造成熔區熱場的 不穩定，因此我們在銅片的邊緣加裝蛇管(coil)作為冷卻用，使得 銅片降溫至 60℃，可以有效的阻止氧化銅的產生(圖十七(b))。 (3) 晶體生長操作時，可發現熔區上方的熱對流強勁，在此過程中， 熔區中的鋰會與氧氣形成氧化鋰蒸汽附著在觀察區的石英蓋上， 使得石英蓋霧化。在此方面的雖然系統中會通以空氣將氧化鋰帶 離熔區上方但是若通過的空氣流量太大則會影響熔區的大小及其

的影響達到最小。 圖十七(a) 銅片表面產生氧化銅 圖十七(b) 加裝蛇管後，銅片表面 無氧化銅產生 在區熔法中，熔區的大小為控制生長晶體品質的重要關鍵， 在熔區大小的控制上，與白金片的大小及所通過的電流量有極為重要 的關係。我們採用長 8 公分寬 0.3 公分厚 0.05 公分的白金片，通予 80 安培的電流可以得到寬度約 0.7 公分大小的熔區，但是在達到 LN 熔點溫度的同時時，因為白金片在高溫下結構改變造成斷裂。 藉此經驗，我們可以得到在白金與電流的選擇上的經驗，關於加 熱片的選擇上，我們將選擇較厚的白金片，其厚度約為 0.1 公分、長 為 12 公分、寬為 1 公分，將可承受在高電流及高溫下白金結構的改 變，避免白金片在生長過程中斷裂；而在電流的選擇上，配合歐姆定 律我們可以計算出所需的電流大小以及可以得到的熔區尺寸。並藉由 此測試過程同時進行此長晶系統的試車工作，由此以進一步改良儀器 之設計。

3-2 模擬部分 關於局部模式方面，我們以最簡單的例子做測試，當我們不考慮 熔液的流動，只考慮熱傳導，固液界面因熔化熱的釋放而變形。圖十 一為固液界面的位置與固液界面的高度等高線，由圖十八可知，在相 同的加熱條件與物理參數之下，兩固液界面是完全一樣的。 凝固界面 熔化界面 圖十八 局部模式界面位置圖與高度等高線圖 關於整體模式方面，我們將石英管外溫度設為一高斯分佈，兩端 的溫度設為定值，此溫度分佈所加熱的區域，即 a 值，與局部模式所 加熱的區域大小是相同的，惟一不同的是最高溫的溫度。圖十九為我 們要輸入至局部模式的表面熱通量分佈圖。 熔化界面 凝固界面 cm x cm x 95 . 2 18 . 3 min max = = cm x cm x 81 . 3 05 . 4 min max = =

圖十九 表面熱通量分佈圖 此次報告主要討論雙層模式的物理系統與運作模式，我們分別討 論整體模式與局部模式，將來會將此二模式結合在一起，並且研究其 模式進行雙向之可能性。未來模擬的結果也將與實驗進行比較與討 論，找出兩者之差異。

第四章 結論及建議

4-1 系統之改進與建議 水平式區熔法晶體生長系統之測試，其熔區大小可以受到電流大 小的控制，然而由於熔區大小以及系統裝置仍然並非最為妥善，所以 在晶體生長過程中熔區仍然未完全穩定。此外，由系統之測試可歸納 出幾點改進要點 (1) 針對操作中之觀察方便，必須進一步修改儀器設計以減少氧化 蒸汽對觀察區的影響，方便實驗時之操作及觀察。 (2) 對於白金片的加熱過程中，其連接白金片與變壓器間的銅板會因 為高溫的關係而造成氧化，產生絕緣成造成電路的不穩定，因此 在銅片上加裝冷卻裝置可以減少氧化銅對系統之影響。 (3) 在白金片的選擇上，由於太薄的白金片在加熱的過程中會因會結 構的改變，以致於無法在長時間下承受高電流及高溫的環境而熔 化斷裂，因此需要選擇較厚及較長的白金片，以減少電流及高溫 對白金片結構的影響，我們將選擇長 12 公分、寬 1 公分、厚 0.1 公分的白金片作為浸入式加熱器的加熱體。

4-2 預計進行之研究 (1) 找出最符合經濟效益及加熱效果的白金片尺寸與所需的電流大 小，使得白金片的耗材及加熱的耗能達到最小值。 (2) 針對承載白金船的承載裝置做最佳的設計，以增加對白金舟的保 溫，減少加熱熔區所需的耗能 (3) 將生長後的晶體做晶體成分及品質分析 (4) 進一步改良 HZM 裝置，減少使用之白金坩堝的重量以及白金舟 對生長後晶體的影響，以減少耗費於坩堝上的成本。 (5) 數值模擬方面，將整體模式與局部模式結合在一起，並且研究其 模式進行雙向之可能性。 (6) 未來模擬的結果也將與實驗進行比較與討論，找出兩者之差異， 分析最佳長晶程序。 (7) 探討溫度梯度及長速對晶體品質及均勻性的影響，進而找到得到 平整界面應有的溫度梯度及可能的最高生長速度。

表一 物質物理性質與系統參數設定 Tube T ρ = 1750 kg m-3 T κ = 1.4 W m-1K-1 T Cp = 740 J kg-1k-1 ε_T = 0.35 Insulator i ρ = 3990 kg m-3 i κ = 0.4 W m-1K-1 i Cp = 740 J kg-1k-1 ε_i = 0.4 Crystal c ρ = 5500 kg m-3 c κ = 17 W m-1K-1 c Cp = 390 J kg-1k-1 ε_c = 0.3 Melt m ρ = 5500 kg m-3 m κ = 39 W m-1K-1 m Cp = 390 J kg-1K-1 ε_m = 0.3 β_m = 5x10-4 K-1 Feed f ρ = 5500 kg m-3 f κ = 17 W m-1K-1 f Cp = 390 J kg-1k-1 ε_f = 0.3 Air air ρ = 1.603 kg m-3(ideal gas) air κ = 0.05 W m-1K-1 air Cp = 523 J kg-1k-1 μair = 2.27x10-5 kg m-1s-1 Ampoule a ρ = 1800 kg m-3 a κ = 326W m-1K-1 a Cp = 1814 J kg-1k-1 ε_a = 0.4

Other physical parameter

ΔH = 460 J g-1 Tm = 937.4 oC a= 0.5 cm h = 4.67x10-1 W cm-2 K-1 g = 9.81 m s-2 Input parameters L = 7 cm L1 = 12 cm L2 = -5 cm Ri = 1.8 cm Rt = 2.0 cm Lc = 3 cm Lm = 1cm La = 3 cm Ra = 0.7cm Rc = 0.5 cm

參考文獻

1. G.C. Valley, M.B. Klein, R.A. Mullen, D. Rytz and B. Wechsler, Anm. Rev. Mater.

Sci. 18 (1998) 165

2. H. Abe, S. Toyoda, Y. Kakizaki, Y. Osugi and Y. Kozuka, J. Japan. Assoc. Crystal

Growth 17 (1990) 60

3. D.M. Sjhields and T.E. Luke, Opt. Commun. 55 (1985) 391

4. J.R. Carruthers, G.E. Prterson, M. Grasso, P.M. Bridenbaugh, J. Appl. Phys. 42 (1979) 1846

5. Y. Zhang, Y.H. Xu, M.H. Li and Y.Q Zhao, J. Crystal Growth 223 (2001) 537 6. I.Foldvari, K. Polgar, R. Voszka, R.N. Balasanyan, Crsy. Res. Technol. 19 (1984)

1659

7. G.I. Moalovichko, V.G. Grachev, E.P. Kokanyan, O.F. Schirmer, K. Betzler, B. Gather, F. Jermann, S. Klauer, U. Schlarb, M. Wohlecke, Appl. Phys. A 56 (1993) 103

8. K. Polgar, A. Peter, L. Kovacs, G.Corradi, Zs. Szaller, J. Crystal growth 177 (1997) 211

9. M. Wohlecke, G. Corradi, K. Betzler, Appl. Phys. B 63 (1996) 323

10. k. Kitamura, J.K. Yamamoto, N. Iyi, S. Kimura, T. Hayashi, J. Crystal Growth 116 (1992) 215

11. D. Serrano, V. Bermudez, L. Arizmendi and E. Dieguez, J. Crystal Growth 210 (2000) 670

12. H. Chen, H. Xia, J. Wang, J. Zhang, J. Xu and S. Fan, J. Crystal Growth 256 (2003) 219

13. J.C. Chen and C. Hu , J. Crystal Growth 149 (1995) 87

14. D.Q. Ni, W.Y. Wang, D.F. Zhang, X. Wu, X.L. Chen, K.Q. Lu, J. Crystal Growth 263 (2004) 421

15. R. Siegel, and J. R. Howell, “Thermal Radiation Heat Transfer”, 2nd edition, 1981 16. C.J. Saltiel, and J. Kolibal, Journal of Thermophysics and Heat Transfer 7 (1993)

175

17. A. Voigt, N. Hanssen, C. Weichmann, Mathematical And Computer Modeling 39 (2004) 145

18. M.C. Liang, P. h. D thesis (1999)

19. H. Zhang, L.L. Zheng, and V. Prasad, J. Heat Transfer 121 (1999) 148

20. F. Dupret, P. Nicodeme, Y. Ryckmans, P. Wouters, and M. J. Crochet, Int. J. Heat

21. E. Dornberger, E. Tomzig, A. Seidl, S. Schmitt, H. J. Leister, Ch. Schmitt, and G. Müller, J. Crystal Growth 180 (1997) 461

22. R. Assaker, N. van den Bogaert, and F. Dupret, J. Crystal Growth 180 (1997) 450 23. L.Y. Huang, P.C. Lee, C.K. Hsieh, W.C. Hsu, and C.W. Lan, J. Crystal growth 261

(2004)

24. A. Yeckel, Y. Zhou, M. Dennis, and J.J. Derby, J. Crystal Growth 191 (1998) 206 25. B. Vartak, Yong-Il Kwon, A. Yeckel, J.J. Derby, J. Crystal Growth 201 (2000) 704 26. C.W. Lan, Chemical Engineering Science 59 (2004) 1437

27. G. Müller, J. Friedrich, J. Crystal Growth 266 (2004) 1

28. I.Y. Evstratov, V.V. Kalaev, V.N. Nabokov, A.I. Zhmakin, Y.N. Makarov, A.G. Abramov, N.G. Ivanov, E.A. Rudinsky, E.M. Smirnov, S.A. Lowry, E. Dornberger, J. Virbulis, E. Tomzig, W.V. Ammon, Microelectronic Engineering 56 (2001) 139 29. G. Müller, J. Crystal Growth 237-239 (2002) 1628

30. D.P. Lukanin, V.V. Kalaev, Y.N. Makarov, T. Wetzel, J. Virbulis, W.V. Ammon, J.

of crystal growth 266 (2004) 20

31. C.H. Panzarella, M. Kassemi, American Society of Mechanical Engineers, Heat

Transfer Division, Micro/Macro Scale Phenomena in Solidification 218 (1992) 127

32. A.J. Nowak, R.A. Bialecki, A. Fic, G. Wecel, L.C. Wrobel, B. Sarler,

Computational Materials Science 25 (2002) 570

33. C. Gourbeyre, T. Chassagne, M. Le Berre, G. Ferro, E. Gautier, Y. Monteil, D. Barbier, Sensors and Actuators, A: Physical 99 (2002) 31

34. A.R. Qubbaj, J. Villarreal, S.R. Gollahalli, American Society of Mechanical