已有大佬做了基于cohesive+泰森多邊形(Voronoi)插件模擬晶粒之間的晶界開裂的視頻，這個視頻中晶粒模型主要是通過homtools插件建立的，cohesive單元的建立是通過Cohesive_generator_2D3D插件實現的。附上視頻鏈接：

   在Neper中建立晶粒模型劃分網格時也可以批量插入0厚度cohesive單元，我之前進行了嘗試，這里給出一個示例。在Neper中生成一個1×1的二維正方形晶體模型，平均晶粒尺寸為0.08，圓形度為0.1，共196個晶粒：

neper -T -n from_morpho -dim 2 -domain "square(1,1)" -morpho "diameq:lognormal(0.08,0.01),1-circularity:lognormal(0.100,0.03)" -oriformat plain -oridescriptor euler-bunge -reg 1 -o cp1

劃分網格并生成inp文件：   

neper -M cp1.tess -rcl 1.2 -elttype tri -pl 15 -interface cohesive -format inp -o abq_input

    注意這里的命令-pl是為了減輕網格的過度細化，避免局部網格過密，這一問題在neper的github討論中提到過，附上地址：https://github.com/neperfepx/neper/discussions/777

    命令-interface cohesive是在晶界處批量插入0厚度cohesive單元，遺憾的是，neper暫時不支持插入非0厚度的cohesive單元（https://github.com/neperfepx/neper/discussions/192）。此外，晶粒網格選擇了三角形單元，neper生成2維cohesive單元時似乎不支持四邊形單元劃分（https://github.com/neperfepx/neper/discussions/98），另外一個討論中（https://github.com/neperfepx/neper/discussions/481）作者說明了它不適用于六面體網格，但Neper2CAE似乎可以完成這項工作，感興趣的讀者可以自行查找。

    從neper導出inp文件后，需要對該inp文件進行一個簡單的前處理，為每一個晶粒建立單獨的材料并輸入材料參數，這部分工作可以參考Github上的工作Neper2Abaqus（https://github.com/smrg-uob/Neper2Abaqus?tab=readme-ov-file）。Neper2Abaqus是通過Bunge歐拉角和晶體坐標向量計算全局坐標向量并輸入到abaqus的材料屬性中，此外還需要為inp文件的每個晶界集合創建單獨的截面，以便指派材料屬性，并設置粘聚力模型的材料參數。

   在abaqus中導入inp文件之后，可以根據自身需要修改cohesive單元的控制屬性：

    這里對模型施加了一個拉伸荷載，采用黃永剛老師的晶體塑性子程序huang.for，然后提交作業就可以計算得到應力云圖：

這里附上計算文件：