基于FEM-SPH耦合算法的磨粒仿真研究

程蔚

瀏覽：6943 評論：14 收藏：57

基于FEM-SPH耦合算法的磨粒仿真研究

1.工程背景

隨著半導體行業的蓬勃發展，單晶碳化硅作為典型的第三代半導體材料被廣泛應用于集成電路生產、光學襯底材料制備等加工過程中，晶片表面質量的好壞直接決定了半導體器件的使用性能及工作壽命^[1-2]，這也就對以單晶碳化硅為代表的光學材料的加工質量提出了更高的要求，其中在研磨拋光過程中產生的表面/亞表面損傷缺陷（見圖1-1亞表面損傷示意圖）是影響晶片性能的重要因素，而對損傷層（SSD）中的裂紋產生、擴展表征與控制更是不容忽視。通常，通過研磨拋光單晶碳化硅晶片實驗來觀測裂紋是行之有效的手段，但裂紋產生、擴展是一個復雜的動態過程，通過實驗后觀測的手段顯然與晶片表面裂紋擴展的真實情況相差甚遠，其次單晶碳化硅屬于典型的硬脆材料，其在加工過程中的裂紋產生、擴展往往在微秒間完成，這也是難以用實驗的手段對裂紋擴展進行觀測的另一原因。因此本文采用ANSYS/LSDYNA軟件建立了單顆磨粒劃擦單晶碳化硅模型，通過仿*真手段來分析單晶碳化硅在研磨加工過程中的裂紋產生及擴展情況。

2.研磨加工模型分析及模型簡化

2.1研磨原理

固結磨磨料研磨加工系統主要由研磨盤、載物盤、研磨液供給系統組成。研磨時，向下壓力P使緊貼上盤面作自轉運動的工件與下盤面接觸作公轉，依靠相對運動實現磨粒對工件的研磨加工^[3]。研磨加工系統及示意圖如圖2.1所示。本文所采用的圖案磨盤為課題組自行設計的雙圖案排布磨盤，在實驗中表現出比傳統圖案磨盤更高的效率和加工精度，其微粉金剛石磨粒通過釬焊技術排布在磨盤上，釬焊效果圖如圖2.2所示。此外因研磨相對磨削具備更低的轉速，在加工過程中的散熱及振動很小，故忽略加工系統的振動及研磨液對磨粒加工的熱作用。

2.2模型簡化

為了使仿真模型更加準確反映實際磨粒加工導致的損傷及裂紋擴展情況，做出如下三點模型簡化，最終模型簡化示意如圖2.3所示。

（1）研磨過程簡化成單顆磨粒對工件的變切深刻劃過程，磨粒設為圓錐狀，工件設為長方體。理由：研磨是通過磨盤上的微粉金剛石磨粒與工件之間的相互沖擊、劃擦等力學物理作用來去除材料（見圖2.1），因此，將研磨盤加工等效為單顆磨粒加工；在極短時間內研磨盤的復雜運動可以等效為直線劃擦運動；此外根據實際金剛石磨粒壓頭形狀，將磨粒設為圓錐形狀，工件設為長方體。

（2）忽略加工系統的振動及工件材料的熱傳導。理由：研磨相對磨削具備更低的轉速，在加工過程中的散熱及振動很小。

（3）磨粒設為剛體。理由：磨粒為金剛石材質，為自然界最硬的物質，其硬度、強度遠遠大于單晶碳化硅。

2.3模型幾何尺寸

刻劃實驗中采用圓錐角為120°，圓錐倒角半徑0.2mm的金剛石壓頭，故仿真中設置磨粒特征尺寸與之相同。工件尺寸與磨粒特征尺寸配合到同一尺度下即可，均為mm級，設為2×1×0.3mm³。

長方體工件（長×寬×高）：2.0mm×1.0mm×0.3mm

圓錐磨粒尺寸：圓錐角120°，圓錐倒角半徑R=0.2mm

3.仿真方法

3.1仿真算法選擇

FEM-SPH（Finite Element Method-Smoothed Particle Hydrodynamics）耦合算法，即有限元-光滑質點流體動力學耦合的方法。

本文所用的加工材料單晶碳化硅屬于典型的硬脆材料，對于硬脆材料的研磨加工仿真，單純用有限元法（FE）來模擬往往會因為網格畸變而導致計算精度低甚至計算強行終止的情況。這是因為在磨屑的形成過程中，材料的去除會經歷彈性變形階段、彈塑性轉變階段、塑性變形階段、塑脆性轉變階段及脆性斷裂，且這一過程中工件材料會產生大量的不連續裂紋。因此，網格劃分的精度及材料失效判據的參數選取直接影響了模型計算的精度與效率，而采用SPH方法則避免了上述問題^[4]。

SPH方法是一種純Lagrange方法，其基本思想是:將連續的流體（或固體）用相互作用的質點組來描述，各個物質點上承載各種物理量，包括質量、速度等，通過求解質點組的動力學方程并跟蹤每個質點的運動軌道，求得整個系統的力學行為。SPH用質點代替了網格，因此不會出現網格大變形問題，同時采用SPH方法進行建模不受尺度的限制，對于微納尺度的動態仿真（本文工件模型尺寸為μm級）結果仍然具有較高可信度。但FE方法在計算效率上相比SPH方法仍然具有極大的優勢，因此本文采用FEM-SPH耦合方法建立單顆磨粒劃擦單晶碳化硅模型。

3.2仿真流程

本文選用ANSYS 19.0版本進行模型建模相關操作，全程建模均采用統一的g-mm-ms單位制。首先，在Workbench中選用Workbench LSDYNA模塊，完成磨粒和工件的有限元模型工作。然后忽略其他信息，直接將該模型保存為x_t格式導入到ANSYS 經典版中，選用ANSYS/LSDYNA模塊進行單元、材料參數、網格劃分工作；單元選用顯示3D SOLID164實體單元，另外選擇顯示Thin Shell163單元用于圓錐模型的映射網格（mapped）劃分，之后設置金剛石磨粒的材料本構（rigid），包括：密度、楊氏模量、泊松比；工件的材料本構選用JH-2模型，由于ANSYS/LSDYNA中并沒有此種本構模型，故本文采用替代法設置：即先賦予工件一種彈性材料本構方便后續網格劃分，模型檢查等工作，之后通過在LSDYNA的二次開發軟件LS-PREPOST中修改導出的K文件，將工件彈性材料本構再次替換成JH-2本構，完成工件材料的材料賦予。完成單元、材料定義后，在對模型進行網格劃分，磨粒采用mapped劃分，工件采用六面體網格劃分，網格劃分完成后可以通過Part options查看工件Part的總網格數對模型網格大小（精度）進一步調整；之后對仿真總時長、磨粒初始刻劃速度、邊界條件及自由度約束進行設置（這一步驟也可以在LS-PREPOST/UE軟件中通過編輯關鍵字進行設置），之后對沙漏、動能、內能進行控制（一般設為默認），然后在ASCII OUTPUT中設置需要輸出的物理量（能量、接觸力、損傷等），最后寫出K文件，然后在LSDYNA的LS-PREPOST前后處理模塊打開保存的K文件，進行網格化工件模型到SPH工件模型的轉化替換工作及關鍵字的刪除、添加編輯，之后設定SECTION-SPH關鍵字，選用材料模型MAT 110（JH-2），設定材料模型所需參數，并將材料關聯到工件Part；之后設定DEFINE-CURVE 關鍵字定義磨粒軌跡與速度并關聯到rigid Part部分。對磨粒的約束通過PRESCRIBED_MOTIOM_RIGID關鍵字定義，較為簡單，而對SPH模型的約束不同于網格化模型，需定義BOUNDARY及 SPH_SYMMETRY_PLANE 關鍵字，對 SPH 質點采用虛粒子約束法，即在 SPH 質點周圍假設仍有其他粒子對模型包圍從而對其形成約束。對磨粒與SPH工件的耦合接觸也不同于以往有限元網格模型，本文仿真應用了AUTOMATIC_NODES_TO_SURFACE 關鍵字進行定義。對于工件損傷裂紋的查看需要借用UE編輯器修改EXTENT_BINARY關鍵字完成，將所有已經定義的關鍵字在 Part 中進行關聯后，并用UE編輯器檢查K文件后，整個求解模型及模型檢查完成。

將上述結果仍保存為K文件格式，添加到 ANSYS 中的 LS-DYNA Solver 中進行最終求解計算，求解得出的 D3PLOT 文件用 LS-PREPOST 打開來獲取應力應變云圖、損傷圖、接觸力曲線等仿真結果。至此，磨粒劃擦模型分析基本結束，為了后期方便對模型修改，本文中的關鍵字修改可采用UE編輯器替代LS-PREPOST對K文件的編輯工作，其優點是避開繁瑣的多次軟件的打開、關閉及保存操作，直接在軟件底層對關鍵字命令編輯，省時高效。建模流程如圖 3.1所示。

3.3仿真難點

（1）網格化工件到SPH工件的粒子化轉化與粒子間隔的控制（前后模型尺寸不變、材料相同（包括本構替換和密度替換）、粒子間隔可控制）

對于磨粒仿真，一般采用分部建模法：在ANSYS軟件中建立磨粒的有限元模型，在LS-PREPOST軟件中建立工件的SPH模型，但此種建模方法不利于后期模型尺寸的對比修改，故本文采用集中替換法，即直接在ANSYS WORKBENCH中建立磨粒和工件的有限元模型，之后在LS-PREPOST軟件中進行SPH轉化與替換。這樣處理的優點是：便于在同一建模環境下（ANSYS 建模環境）對幾何模型尺寸的對照修改與集中編輯。在ANSYS中建立好磨粒、工件系統模型后在LS-PREPOST中進行FE-SPH轉化，如圖3.2所示。為了仿真計算的準確性，在FE-SPH轉化過程中，應滿足：前后模型尺寸不變、材料相同（包括本構替換和密度替換）、粒子間隔可控制3個條件。對應軟件操作：在LS-PREPOST中編輯SPHGEN關鍵字創建SPH模型，并選擇Solid Nodes創建方法；之后在密度（Den）一欄填入單晶碳化硅的真實密度（g-mm-ms單位制下）完成材料密度的替換；最后通過UE編輯器編輯保存的K文件，將預設的彈塑性材料修改為*MAT_110號材料（JH-2本構）完成本構替換。另外，為保證SPH粒子間隔的可控性，SPH粒子間隔h₀與替換前的有限單元網格大小應具有確定的對應關系，圖3.2給出了一種對應方法（即上文提及的Solid Nodes創建SPH法）：有限網格單元的4個節點（分別表示為1、2、3、4）對應轉化為4個SPH粒子，同時網格大小等于粒子間隔長度h0。那么在ANSYS進行網格劃分時就可以通過設置不同尺度的網格來控制SPH粒子間隔h0（l=h0），這極大方便了本文仿真模型中對于切深的設置（見本文表1所示）。

（2）單晶碳化硅的損傷參數如何輸出、裂紋的云圖后處理如何顯示

JH-2本構模型包含有關于硬脆材料的損傷失效模型算法，因此無需另外定義單元失效準則就用模擬出單晶碳化硅的損傷及裂紋擴展現象。這也是JH-2本構能夠模擬單晶碳化硅加工裂紋擴展的本質原因。

為了更好地輸出損傷參數，在仿真計算中首先識別非損傷區與損傷區，將斷裂應變滿足D₁[P*+T*]D₂
≥E_fmin的SPH粒子標紅，以示區別，這部分區域就表示加工過程中的損傷。之后將輸出裂紋的損傷參數定義附加寫入D3PLOT文件的時間歷程變量數目（history #2）。軟件操作步驟：用UE編輯器編輯關鍵字*DATABASE_EXTENT_BINARY，將第2個時間歷程變量修改為2（即NEIPH=2），該參數就表示陶瓷材料的損傷參數。對應的關鍵字編輯如下：

*DATABASE_EXTENT_BINARY

$NEIPH,NEIPS,

2,0

（3）SPH工件邊界的處理（“虛粒子約束”法）

SPH算法是模擬大變形、防止網格畸變的重要手段，但相對FEM算法，其計算效率較低，同時SPH粒子的邊界條件較難處理。對于此種問題，文獻中^[5]通常采用“虛粒子約束”法來定義SPH模型的工件，本文亦采用“虛粒子約束”法處理單晶碳化硅SPH工件的邊界問題。圖3.3給出了“虛粒子約束”算法的基本原理：“虛粒子約束”法是對工件可能移動方向進行約束，靠近SPH粒子邊界處2h（h為粒子光滑長度）范圍內設置出虛粒子。對于靠近邊界的SPH粒子，通過對自身的映射，自動創建具有相同質量、壓力、絕對速度的虛粒子，使得真粒子能正常進行鄰域搜索，以達到約束邊界的目的。與之對應的關鍵字是*BOUNDARY_SPC_SYMMETRY_PLANE。

（4）有限網格磨粒與SPH粒子化工件的耦合接觸

不同于FE算法，本文中SPH模型以粒子代替網格，相當于有限網格的磨粒與粒子化工件的不連續加工過程。因此，FE的面面接觸算法已經不再適用。本文對于有限元單元與光滑粒子接觸界面的相互作用（磨粒與工件），則是通過罰函數算法來定義，耦合接觸算法采用自動點面接觸算法，主面設為磨粒，從面設置為SPH工件（MSTYP=3，SSTYP=4），其對應的關鍵字為*CONTACT_AUTOMATIC_NODES_TO_SURFACE。本文針對接觸算法經過多次仿真實驗驗證，得出結論：自動點面接觸算法、侵蝕點面接觸算法（*CONTACT_ERODING_NODES_TO_SURFACE）均可應用于FEM-SPH耦合模型的接觸設置中，但自動點面接觸算法的計算效率更高且計算不易報錯。故本文最終選用自動點面接觸耦合算法。

4.仿真參數

仿真參數的設置原則一般遵循三點^[6]：一是計算時間合適，二是結合實際工藝參數，三是適當放大來凸顯作用規律。

4.1加工參數的設置

本文中磨粒的變切深刻劃是通過磨粒的運動完成，工件底面設為全約束。根據實際研磨實驗中磨粒相對工件的速度、研磨盤直徑及加工深度^[7]，并適當合理放大，設定磨粒的初始速度為50m/s，最大切深設為30μm，磨粒從切入工件到離開工件切深的變化范圍為0-30μm。具體仿真參數的設置如下表1所示。

表1 仿真參數設置

參數	數值
工件尺寸	2×1×0.3mm
切深（ap）	0μm-30μm-0μm
模型粒子數量	371280
粒子間隔	0.012mm

4.2材料本構參數的設置

磨粒的本構設置較為簡單，磨粒為金剛石，設為剛體（密度3560kg/m³、泊松比0.2、楊氏模量1000GPa），工件為單晶碳化硅材料，在碳化硅材料的變切深刻劃過程中，材料的應力應變呈現非線性關系，這種情況下通常使用的彈塑性材料本構模型已經不再適用，目前針對陶瓷、碳化硅等硬脆性材料的研究，學者們普遍使用 Johnson 和 Holmquist 提出的 JHC（Johnson-Holmquist-Ceramic）本構模型（JH-2）來進行表征^[7]。此模型能夠較為準確反映碳化硅加工過程中的動態損傷及裂紋擴展情況，適合本文的仿真要求。單晶碳化硅的 JH-2 材料本構模型參數是在查閱文獻的基礎上，通過多次模擬與實驗驗證獲得的。表2給出了單晶碳化硅的本構參數。

表2 SiC的JH-2本構參數^[6][9][10]

參數	取值	參數	取值
R（kg/m³）	3215	EPSI	1
G（GPa）	193	T（GPa）	0.75
A	0.96	SFMAX	0.11
B	0.35	HEL（GPa）	11.7
C	0.009	PHEL（GPa）	5.13
M	1	BETA	1
N	0.65	K1（GPa）	220
D1	0.48	D2	0.48
K2（GPa）	361	K3（GPa）	0

5.結果分析

5.1材料去除狀態分析

磨粒從0μm切深位置切入工件到工件的最大深度30μm處，之后做返程運動切出工件，在整個變切深刻劃過程中（刻劃距離50×0.04=2mm），材料的去除形態為“檳榔”形溝槽，由圖5.1材料去除狀態圖可知：材料的去除方式主要為脆性斷裂，磨屑的形狀主要是粉末狀。這主要是由于加工材料為單晶碳化硅，其高硬度和低斷裂韌性導致材料很容易發生脆性斷裂，加工磨屑因此呈現粉末形態，這也是其不同于加工鑄鐵等彈塑性材料形成的條狀切屑的本質原因。在LS-PREPOST通過三點截面設置法（選取417139、380731、380600三點）觀察磨粒加工截面[如圖5.2（a）所示]發現：工件破碎深度大于最大切深30μm，磨粒兩側出現磨屑堆積。磨粒刻劃過的區域磨屑堆積嚴重，磨粒未刻劃區域也有少量磨屑堆積。這是因為工件以脆性斷裂方式去除材料時，在Z向不斷的壓力作用下，裂紋向下擴展；而X向磨粒兩側對工件的擠壓效應是導致材料產生橫向裂紋的原因，磨粒兩側的擠壓使得工件內外產生應力差，工件受拉后在靠近磨粒側面處產生橫向裂紋，橫向裂紋擴展到工件表面形成破碎，便產生磨屑堆積在磨粒兩側，具體如圖5.2（b）所示。

5.2刻劃力分析

在ANSYS對模型進行前處理時定義輸出 RCFORC，即可用 LS-PREPOST 提取出單顆磨粒在變切深刻劃過程中受到的切向力、法向力和軸向力。刻劃力如圖5.3所示。

分析：由圖5.3（b）可知，X/Z始終小于1，法向力遠遠大于切向力；變切深加工過程中磨粒的切向力和法向力先隨著位移的增大而增大，當位移為1.0mm左右時，兩者力的大小開始減小，直至為0；而軸向力大小始終在0左右上下波動，基本不變。這主要是由于磨粒剛與工件接觸時刻劃力為零，隨著磨粒逐漸進入工件，切深不斷增大，磨粒與工件的接觸面積也隨之增大，從而導致了切向力和法向力的增大，當磨粒切深減小時，磨粒逐漸向上運動離開工件，切向力和法向力又開始對稱減小。劃刻劃力的波動則是由于磨屑的產生和飛濺所造成的。從圖中也可看出，當位移為0.3mm左右時，法向力急劇增大，這表明工件已經出現脆性斷裂。

圖5.3 刻劃力及力比（a）刻劃力（b）力比

5.3裂紋擴展分析

圖5.4為不同切深下單晶碳化硅襯底的裂紋產生及擴展情況。從沿刻劃方向（X向）的中位剖面圖中可以看出，在t=0.7μs，磨粒切入工件2.06μm，可以看到工件表面剛剛開始出現相連的紅色粒子，此時單晶碳化硅還處于彈塑性變形階段，并沒有出現裂紋；但t=1.1μs，磨粒切入工件3.29μm，紅色粒子繼續增加，而當t=3.2μs，磨粒切入工件9.52μm，紅色粒子開始飛濺，表明工件表面開始出現脆性斷裂，即裂紋開始產生；在這之后，工件進入脆性破碎加工階段，磨粒下方的 SPH 粒子存在大量致密的且相連的紅色粒子，說明出現明顯的裂紋擴展現象。裂紋不斷擴展、加深如圖5.4（d）所示。

為了計算在不同切深下（27μm、30μm）裂紋擴展的深度，自t=3.2μs后（裂紋產生），對裂紋深度進行統計分析，如圖5.5所示。

分析：從圖5.5（a）中可以看出，在27μm切深范圍內，可以發現大量相連的致密的紅色 SPH 粒子，說明工件加工面存在明顯的損傷。在 27μm切深下方約 33μm深的范圍內也存在大量的紅色 SPH 粒子，說明此時單晶碳化硅工件的亞表面損傷深度約為 33μm。雖然沿深度方向繼續往下也能發現少量紅色 SPH 粒子，但原本位置上的藍色SPH 粒子并沒有發生明顯的位移，所以這部分區域沒有發生明顯損傷，而此區域的少量紅色 SPH 粒子則是由于去除的SPH粒子飛濺溢散導致的。而在最大切深30μm范圍內，大量紅色SPH粒子致密相連更加明顯，且工件表面出現114μm左右的橫向裂紋，說明工件表面損傷愈加嚴重。此時損傷深度為26μm，有所減小，說明損傷深度并不是隨著切深增大而增大的。

6.仿真結果驗證

為了驗證本論文中所建立的變切深刻劃FEM-SPH 仿真模型準確性，應該進行仿真實驗與理論分析和實際實驗中的結果進行比較，當仿真實驗結果與實際實驗結果在規律和量級上保持一致時，認為本文中的仿真模型具備一定的準確性。

（1）從材料能量變化上來分析材料去除變形及損傷的大小是在仿真中獨有的優勢和特點^[10]。在LSPP中通過定義matsum可以查看SPH工件的能量變化。文獻^[11]指出，工件出現劃痕是首先磨粒對材料擠壓做功，使材料內部能量逐漸積累，當材料內部增加的能量難以維持平衡，材料開始對外做工，即出現材料變形和去除，而由于變形出現，對外做功，使材料內部能量下降。由圖5.6可以得出：磨粒在開始與工件接觸后，材工件內部能量迅速增長，形成尖峰，表明在研磨初期，磨粒能量完全被工件吸收且工件沒有對外做工，因此表現為塑性變形，沒有損傷和脆性斷裂出現，之后工件能量急劇下降，表明此時脆性斷裂出現，能量釋放轉化為磨屑動能、工件變形能、熱等形式。這與前文通過損傷云圖分析所得結果基本一致，同時驗證了仿真結果的前后一致性。

圖5.6 工件能量變化圖

（2）法向力、切向力對比驗證

對于磨粒的變切深刻劃仿真中的法向力，可以通過金剛石磨粒劃擦試驗輸出法向力、切向力與其對比來驗仿真力的準確性。文獻中^[12]實際實驗所得到的穩定狀態下法向力最大為20N左右，切向力最大為5N大小左右。圖5.7為本文仿真得出的力-時間曲線圖。對比可以發現仿真中的法向力最大為25N左右，切向力最大為10N左右，與文獻中的實驗結果數值處于同一量級，且數值上較為接近。由此驗證了本文仿真仿真模型的準確性。