1問題的提出

單純用FEM算法建立有限元網格模型在模擬大變形問題經常會出現網格畸變，且FEM算法在模擬不連續的問題，如斷裂等問題并不具有優勢，SPH算法由于不用依賴網格算法，可以很好解決這一問題，但隨之帶來的邊界難以處理，計算效率低的問題也一直難以很好解決。為此本文嘗試用FEM-SPH耦合算法來耦合兩者優點，以期獲得理想的仿真結果。本文以單顆球形磨粒等切深劃擦碳化硅工件的FEM-SPH耦合模型為例，驗證這一耦合算法的高效性、正確性。

2 FEM-SPH耦合模型算例

2.1模型建立

圖2-1磨粒仿FEM-SPH模型

由于在磨削加工中，實際是金剛石磨粒的刀尖圓弧半徑劃過工件表面實現的材料去除，因此在介觀尺度下，不規則形狀的磨?？梢院喕汕蝮w，工件簡化成與磨粒尺度相匹配的長方體，工件在7.5μm的切深范圍內采用SPH算法建模，剩下部.分采用FEM算法建立有限元網格，SPH粒子總數為144000個，粒子間隔為0.25μm，SPH粒子下的FEM網格工件網格大小并不影響計算結果，為提高計算時間，可適當取大網格間距，本文中取1μm，即4個SPH粒子與1個有限元網格匹配。磨粒仿真模型如圖2-1所示。幾何模型的具體參數如表2-1所示。因為磨粒為金剛石材質，其硬度和彈性模量遠遠大于單晶碳化硅工件，因此在研磨過程中，磨粒幾乎不會發生變形，因此將磨粒（密度3560kg/m3、泊松比0.2、楊氏模量1000GPa）設為剛體。單晶碳化硅是典型的的各向異性材料，本文仿真選用6H-SiC，單晶碳化硅（6H-SiC）工件的本構參數如表2-2所示。

表2-1幾何模型具體參數

參數	數值
磨粒半徑 切深	5μm 4μm
工件整體尺寸	15×15×10μm3
SPH工件尺寸	15×15×7.5μm3
FEM工件尺寸	15×15×7.5μm3

表2-2 6H-SiC的JH-2本構參數[1,3]

參數	取值	參數	取值
R（kg/m3）	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	PHE（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

2.2算例實現

2.2.1選用的軟件

本文為方便建模過程，采用聯合建模。選用的軟件有：ANSYS 19.0、WORKBENCH協同仿真平臺（主要用于幾何建模）、LSDYNA（用于K文件求解）、LS-PREPOST（用于模型前后處理、K文件修改、模型檢查）、UE編輯器（K文件修改）。

2.2.2聯合建模的單位制

全程建模均采用統一的μg-μm-μs單位制。

2.2.3建模步驟

首先，在Workbench中選用Workbench LSDYNA模塊，完成磨粒和工件的有限元模型工作，建立工件時為方便后面工件的SPH粒子化及部分工件的網格劃分，通過平面分割命令將工件分成兩個PART,球形磨粒可以直接調用sphere命令生成。在完成磨粒和工件的幾何建模后，忽略其他信息，將該模型保存為x_t格式導入到ANSYS 經典版中，繼續選用ANSYS/LSDYNA模塊進行單元、材料參數、網格劃分工作；單元選用顯示3D SOLID164實體單元，之后設置金剛石磨粒的材料本構（rigid），包括：密度、楊氏模量、泊松比；4H-SiC的材料本構選用任意一種彈性材料代替即可，后期在LSPP軟件中直接修改工件的本構參數，當然也可以直接通過UE編輯器修改工件本構的關鍵字。在完成單元定義、材料賦予工作后，開始對磨粒和工件模型網格劃分，磨粒采用四面體網格劃分，網格單元大小設為1μm，工件的網格劃分需要注意網格的大小。對于磨粒正下方的這部分工件，由于直接承擔著磨粒的劃擦作用，因此網格單元應當加密，對于下方的工件網格可以稀疏，其網格大小本質上并不會影響算例的準確性。因此，上方工件采用六面體網格劃分，網格大小為0.25μm，下方網格同樣采用六面體網格劃分，但單元大小取大網格劃分，設為1μm。網格劃分完成后，可以對仿真時間、磨粒劃擦速度、磨粒自由度約束、工件約束、能量控制等進行設置，也可以不設置，留在LSPP軟件中對這些參數進行集中處理，這也是聯合建模的最大優勢所在：在聯合建模中可以不必嚴格按照前處理→求解→后處理的CAE分析步驟進行，不必擔心參數設置不全的問題。之后設置需要輸出的三個物理量（能量、接觸力、損傷）寫出K文件完成ANSYS經典界面下的分析過程。然后在LSDYNA的LS-PREPOST前后處理模塊打開保存的K文件，進行網格化工件的SPH替換、工件材料的替換、接觸約束邊界條件的設置等操作。對于SPH工件的轉化一定要注意實際建模尺寸與轉化成SPH尺寸之間的對應關系，在選用Solid Nodes SPH創建方法時，實際建模尺寸與SPH尺寸之間是對應相等的（本文選用此種方法創建SPH粒子），而在選用Solid Center方法創建SPH粒子時，轉化后的SPH粒子總長度是要小于實際建模尺寸的，這是由于Solid Center是將每個網格的中心點轉化為一個SPH粒子，這可以理解成每一個網格都簡化成一個具有網格質量位于網格中心的理想質點，因此網格轉化成點，在尺寸上就減少了一個粒子間隔長度。因此，選用Solid Nodes方法創建上方SPH工件，下方網格工件保留不作處理。之后設定SECTION-SPH關鍵字，選用材料模型MAT 110（JH-2），設定材料模型所需參數，并將材料關聯到工件Part；之后設定DEFINE-CURVE 關鍵字定義磨粒軌跡與速度并關聯到rigid Part部分。對磨粒的約束通過PRESCRIBED_MOTIOM_RIGID關鍵字定義，較為簡單，而對FEM-SPH耦合工件的接觸設置用固連點面接觸（TIED_NODE_TO_SURFACE）定義，這就避免了單純SPH工件需要定義關鍵字*BOUNDARY及 SPH_SYMMETRY_PLANE來對邊界處粒子進行約束。本文的FEM-SPH耦合之處除了工件之間的耦合，還有磨粒與SPH工件的耦合，對對磨粒與SPH工件的耦合接觸采用AUTOMATIC_NODES_TO_SURFACE 關鍵字進行定義。最后，對仿真時間、沙漏能等進行最后的設置，另外對于工件損傷裂紋的查看需要借用UE編輯器修改EXTENT_BINARY關鍵字完成，將所有已經定義的關鍵字在 Part 中進行關聯后，并用UE編輯器最終檢查K文件后，完成球形磨粒劃擦6H-SiC工件算例。

2.2.4求解算例

最后在LSPP將文件保存為K文件格式，用LSDYNA Solver求解器求解K文件。在LSPP中可以打開binary文件（D3PLOT）查看云圖，繪制二維圖像等。

3 FEM-SPH、SPH對比

仿真效率的高低與電腦的配置有關，本次算例在如圖3-1所示的電腦配置下進行。為了較為客觀比較采用FEM-SPH算法與采用SPH算法建模的計算效率，本文進行了兩次仿真，仿真參數及幾何模型尺寸均設為相同，但采用SPH算法建模的計算效率低下，在設置相同仿真時間1μs時，計算時間實在太長（＞90h），因此經過多次仿真調試試驗后，在設置仿真時間為0.04μs時，采用SPH算法建模的計算時間約為14h37min，這比仿真1μs，采用FEM-SPH算法建模的計算時間（10h42min）還要長（如圖3-2所示），由此可見，采用FEM-SPH算法的計算效率要遠遠高于采用SPH算法建模的計算效率。

圖3-1電腦配置

圖3-2不同算法計算時間對比

4分析與討論

4.1接觸反力分析

在ANSYS對模型進行前處理時定義輸出 RCFORC，即可用 LS-PREPOST 提取出單顆磨粒在等切深刻劃過程中受到的切向力、法向力和軸向力。三向接觸力如圖4-1所示。分析：Y-force從磨粒進入工件開始一直在0上下波動。這是由于Y向與加工方向垂直；X、Z-force從磨粒進入工件后數值急劇上升到30mN、50mN，但Z向力更大，且在t=0.3μs后，Z、X向力出現波動，呈現先減小后緩慢增大再減小的波動，且X向力方向為負，但t=0.8μs后，三向力都可以減小。這是由于磨粒與工件接觸面積逐漸增大，受力面積增大，磨粒不斷壓迫（斜向下）工件破碎，法向力也就不斷增大，而切向力由受力分析可知其方向始終與磨?？虅澐较蛳喾?，即為負值。t=0.8μs后磨粒逐漸離開工件，所以力值逐漸變小。

圖4-1三向接觸力分析

4.2能量分析

在LSPP中通過定義matsum可以查看SPH工件（PART 4）的能量變化。SPH工件的內部能量隨時間的變化關系如圖4-2所示。由Griffith^[4]提出的斷裂能量分析可知：球形磨粒在開始與工件接觸后，材工件內部能量迅速增長，形成尖峰，表明在研磨初期，磨粒能量完全被工件吸收且工件沒有對外做工，因此表現為塑性變形，沒有損傷和脆性斷裂出現，之后工件能量急劇下降，表明此時脆性斷裂出現，能量釋放轉化為磨屑動能、工件變形能、工件表面自由能、熱的形式。

圖4-2 SPH工件能量變化圖

4.3最大等效應力

圖4-3表示劃擦過程中不同時刻的等效應力折線圖。分析：比較大的應力值主要集中在劃擦的初始階段?？梢园l現等效應力的趨勢與三向接觸力的趨勢有一定的對應性。當 t=0.076μs 時，也就是塑脆臨界共存態到脆性模式的轉變點，應力值急劇上升到12.2Gpa，然后持續上下震蕩，在t=0.4μs、t=0.7μs左右時，最大應力值出現尖峰，這與三向接觸力變化趨勢相對應，在此時法向力、切向力以及軸向力（軸向力在0左右波動）達到相應最大值。

圖4-3最大等效應力

   4.4材料去除狀態分析

球形磨粒在研磨加工硬脆工件時，一般是以脆性斷裂的形式去除材料。圖4-4分別給出了在相應時刻的塑性應變、等效應力圖。可以看出磨粒正前方45°范圍內所受應力最為嚴重，也就意味著工件在這一范圍內最容易產生損傷堆積，同時磨粒以粉末狀飛濺去除與實際加工單晶碳化硅工件相一致，也初步驗證了仿真結果的正確性。而從磨粒對SPH工件的三向壓力來看工件損傷區域的方向更為直觀，如圖4-5（a-c）所示，通過FCOMP→Misc→history#2可以看裂紋損傷情況如圖4-5（d）所示。

圖4-4應力應變圖

圖4-5X Y Z向磨粒對SPH工件的壓力顯示云圖（a-c），SPH工件損傷圖（d）

5 總結

經多次調試與仿真實驗，得出結論：

（1）計算效率方面：發現其計算時間比單純用SPH快了16倍（SPH算法耗時約16小時，而采用FEM-SPH耦合算法仿真耗時1小時13分鐘左右），仿真結果文件D3PLOT文件平均十秒一個，仿真過程極度絲滑。

（2）SPH邊界難以處理問題：FEM-SPH耦合模型不用考慮用虛粒子約束法定義SPH邊界，在耦合模型中直接將外層SPH粒子固連在FEM網格上即可。

因此，采用FEM-SPH耦合建模同時避免了單純S9PH算法的邊界難以處理、FEM網格算法難以處理大變形的兩個棘手問題。這為金屬切削、磨削加工、研磨拋光加工過程的仿真提供了一種更加高效的仿真手段，具有一定的理論與實踐意義。

參考文獻

[1]    Groenenboom P H L, Cartwright B K. Hydrodynamics and fluid-structure interaction by coupled SPH-FE method[J]. Journal of Hydraulic Research, 2010, 48(1): 61~73.

[2]    Xu J X, Liu X L. Analysis of structural response under blast loads using the coupled

[3]    SPH-FEM approach[J]. Journal of Zhejiang University SCIENCE A，2008, 9(9): 1184~1192.

[4]    Griffith, A. A.  The  phenomea  of rupture and flow  in  solids[D]. Phil. Trans. Roy. Soc. Lond. 1920, A221,163.

1.k