1問(wèn)題的提出

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

2 FEM-SPH耦合模型算例

2.1模型建立

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

由于在磨削加工中，實(shí)際是金剛石磨粒的刀尖圓弧半徑劃過(guò)工件表面實(shí)現(xiàn)的材料去除，因此在介觀尺度下，不規(guī)則形狀的磨粒可以簡(jiǎn)化成球體，工件簡(jiǎn)化成與磨粒尺度相匹配的長(zhǎng)方體，工件在7.5μm的切深范圍內(nèi)采用SPH算法建模，剩下部.分采用FEM算法建立有限元網(wǎng)格，SPH粒子總數(shù)為144000個(gè)，粒子間隔為0.25μm，SPH粒子下的FEM網(wǎng)格工件網(wǎng)格大小并不影響計(jì)算結(jié)果，為提高計(jì)算時(shí)間，可適當(dāng)取大網(wǎng)格間距，本文中取1μm，即4個(gè)SPH粒子與1個(gè)有限元網(wǎng)格匹配。磨粒仿真模型如圖2-1所示。幾何模型的具體參數(shù)如表2-1所示。因?yàn)槟チ榻饎偸馁|(zhì)，其硬度和彈性模量遠(yuǎn)遠(yuǎn)大于單晶碳化硅工件，因此在研磨過(guò)程中，磨粒幾乎不會(huì)發(fā)生變形，因此將磨粒（密度3560kg/m3、泊松比0.2、楊氏模量1000GPa）設(shè)為剛體。單晶碳化硅是典型的的各向異性材料，本文仿真選用6H-SiC，單晶碳化硅（6H-SiC）工件的本構(gòu)參數(shù)如表2-2所示。

表2-1幾何模型具體參數(shù)

參數(shù)	數(shù)值
磨粒半徑 切深	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本構(gòu)參數(shù)[1,3]

參數(shù)	取值	參數(shù)	取值
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算例實(shí)現(xiàn)

2.2.1選用的軟件

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

2.2.2聯(lián)合建模的單位制

全程建模均采用統(tǒng)一的μg-μm-μs單位制。

2.2.3建模步驟

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

2.2.4求解算例

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

3 FEM-SPH、SPH對(duì)比

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

圖3-1電腦配置

圖3-2不同算法計(jì)算時(shí)間對(duì)比

4分析與討論

4.1接觸反力分析

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

圖4-1三向接觸力分析

4.2能量分析

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

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

4.3最大等效應(yīng)力

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

圖4-3最大等效應(yīng)力

   4.4材料去除狀態(tài)分析

球形磨粒在研磨加工硬脆工件時(shí)，一般是以脆性斷裂的形式去除材料。圖4-4分別給出了在相應(yīng)時(shí)刻的塑性應(yīng)變、等效應(yīng)力圖?？梢钥闯瞿チＵ胺?5°范圍內(nèi)所受應(yīng)力最為嚴(yán)重，也就意味著工件在這一范圍內(nèi)最容易產(chǎn)生損傷堆積，同時(shí)磨粒以粉末狀飛濺去除與實(shí)際加工單晶碳化硅工件相一致，也初步驗(yàn)證了仿真結(jié)果的正確性。而從磨粒對(duì)SPH工件的三向壓力來(lái)看工件損傷區(qū)域的方向更為直觀，如圖4-5（a-c）所示，通過(guò)FCOMP→Misc→history#2可以看裂紋損傷情況如圖4-5（d）所示。

圖4-4應(yīng)力應(yīng)變圖

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

5 總結(jié)

經(jīng)多次調(diào)試與仿真實(shí)驗(yàn)，得出結(jié)論：

（1）計(jì)算效率方面：發(fā)現(xiàn)其計(jì)算時(shí)間比單純用SPH快了16倍（SPH算法耗時(shí)約16小時(shí)，而采用FEM-SPH耦合算法仿真耗時(shí)1小時(shí)13分鐘左右），仿真結(jié)果文件D3PLOT文件平均十秒一個(gè)，仿真過(guò)程極度絲滑。

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

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

參考文獻(xiàn)

[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