唐敬哲，汪 偉，鄭延豐,2，羅堯治,2

(1. 浙江大學(xué)空間結(jié)構(gòu)研究中心，浙江，杭州 310058；2. 浙江省空間結(jié)構(gòu)重點實驗室，浙江，杭州 310058)

基于向量式力學(xué)所提出的有限質(zhì)點法(Finite Particle Method，F(xiàn)PM)是一種面向工程應(yīng)用、面向結(jié)構(gòu)復(fù)雜行為分析的新型數(shù)值計算方法[1]。該方法的本質(zhì)是對質(zhì)點的動力平衡狀態(tài)的描述，控制方程質(zhì)點間相互獨立，并以離散運動路徑的方式處理幾何非線性，在結(jié)構(gòu)復(fù)雜非線性問題分析中具有優(yōu)勢和靈活性。目前該方法已被應(yīng)用于機構(gòu)運動[2?3]、結(jié)構(gòu)倒塌破壞[4?5]、固體彈塑性大變形[6]、薄膜結(jié)構(gòu)形態(tài)分析[7?8]、結(jié)構(gòu)多尺度精細(xì)化[9]等結(jié)構(gòu)復(fù)雜行為的分析。目前，有限質(zhì)點法基于GPU并行加速的相關(guān)研究也已經(jīng)初見成效，這種并行策略也已經(jīng)應(yīng)用在有限質(zhì)點法通用計算平臺的研發(fā)中[10]。借此，有限質(zhì)點法的研究也逐漸向大型工程問題的精細(xì)化數(shù)值模擬的方向所發(fā)展。

界面斷裂、接觸等強非線性復(fù)雜行為廣泛存在于結(jié)構(gòu)破壞、連續(xù)倒塌、上下部結(jié)構(gòu)協(xié)同等實際工程的精細(xì)化分析中。有限質(zhì)點法對于強非線性行為模擬的適應(yīng)性以及并行計算能力的優(yōu)勢使其以較低的成本完成界面行為的細(xì)觀模擬。界面力學(xué)行為的模擬主要面臨兩個問題：一是界面作用對偵察；二是作用力計算。其中界面作用對偵查一向是最消耗時間和計算資源的步驟。目前的大多數(shù)的偵查算法都依照兩個步驟進行，即先通過某種空間排序的手段將搜索區(qū)域局部化，再在局部范圍進行精細(xì)的作用關(guān)系判斷[11]。較有代表性的易于并行實現(xiàn)的偵查算法就是盒圍法[12?13]。這種方法將空間點按某種排序算法進行空間排序，置于空間網(wǎng)格中，將偵察范圍逐漸局部化。特別在將結(jié)構(gòu)進行質(zhì)點化離散的諸多數(shù)值方法中，如物質(zhì)點法、離散元法和光滑粒子流體動力學(xué)等，大多使用這種方法的變體。如Chen等[14]使用桶排序的盒圍法對有限元法與物質(zhì)點法耦合的點-面接觸進行了研究。Gan等[15]使用基數(shù)排序?qū)崿F(xiàn)了離散元的顆粒系統(tǒng)的并行化接觸偵察。Xia和Liang[16]基于四叉樹鄰域搜索對基于光滑例子流體動力學(xué)的盒圍法進行了優(yōu)化。Hu等[17]拓展出了條帶盒圍法進一步提高了搜索效率，成功應(yīng)用在有限元與光滑例子流體動力學(xué)耦合的模擬中。盡管這些偵察算法具有易于并行實現(xiàn)的優(yōu)點，也有不少效果顯著的并行化研究[15? 16,18]，但具體的并行實現(xiàn)仍然是界面作用偵察算法計算效率問題的研究熱點與難點。

在界面作用力計算方面，界面的接觸與斷裂行為采用不同的作用力模型。界面接觸作為常見的界面行為，其常用的力學(xué)模型主要有兩種，即適用于隱式求解的拉格朗日乘子法和通用的罰函數(shù)法，這二者都是通過對接觸區(qū)域施加位移限制條件來達到接觸處理的作用。顯然，罰函數(shù)法適用于有限質(zhì)點法的顯式時間積分框架，因此有過一些成功的應(yīng)用，例如喻瑩和羅堯治[19]在離散化的梁桿系結(jié)構(gòu)中，提出了針對空間點-線、線-線接觸的偵測算法和彈性接觸力模型。除了罰函數(shù)法，張鵬飛等[20]引入了顯式的防御節(jié)點法，確保點-面接觸過程中不會發(fā)生穿透，應(yīng)用于薄殼結(jié)構(gòu)的接觸分析中。鄭延豐[21]同樣使用了防御節(jié)點法，并使用圖形處理器(Graphics Processing Unit,GPU)對精細(xì)化接觸行為的模擬進行了加速。對于界面斷裂，學(xué)者們通常使用基于內(nèi)聚力模型的各種變體[22? 24]來描述斷裂過程中的裂尖應(yīng)力。特別地，張鵬飛[25]基于此模型開發(fā)了可動態(tài)插入的有限質(zhì)點法斷裂單元模擬裂紋展開。這些實現(xiàn)充分利用了有限質(zhì)點法將結(jié)構(gòu)進行質(zhì)點化離散后在強非連續(xù)問題分析中的優(yōu)勢，但也都缺乏普遍的適用性。

為了提高界面力學(xué)模型的廣泛適用性，學(xué)者們在界面區(qū)域構(gòu)造了特殊的單元，如接觸單元、界面單元等，在局部范圍內(nèi)同時完成作用對偵查和作用力計算的過程，早期有Goodman等[26]的零厚度接觸面單元、Desai等[27]的薄層單元等。在此基礎(chǔ)上，諸多商用軟件也出現(xiàn)了通用的單元種類，較有代表性的有LS-DYNA中的接觸單元[28]、FLAC3D中的界面單元[29]以及ABAQUS中的內(nèi)聚單元[30]等，但也存在若干局限性。例如，多以有限元和有限差分法為基礎(chǔ)的商業(yè)軟件，對于存在強非線性的局部界面行為的處理較為復(fù)雜，經(jīng)常出現(xiàn)難以收斂的情況，且計算效率較低。另外，也缺乏統(tǒng)一的單元形式來模擬界面粘結(jié)、接觸、斷裂的耦合情況。

為了充分發(fā)揮有限質(zhì)點法在處理強非線行為中的優(yōu)勢以及并行優(yōu)勢，本文提出了一種適用于有限質(zhì)點法并行計算體系下的三角形界面單元，能夠通用化地模擬界面之間的粘結(jié)、斷裂、接觸等界面作用行為以及之間的耦合。此界面單元基于GPU并行模式進行高效的界面作用偵察，使用內(nèi)聚力模型和罰函數(shù)法作為界面斷裂和接觸力計算的力學(xué)模型，具有廣泛的適用性。本文首先詳述了其并行化的作用對偵察，并介紹了界面在粘結(jié)、斷裂、以及接觸這三種作用下的狀態(tài)判據(jù)和作用力計算方法。最后，依托于有限質(zhì)點法通用計算平臺FPM[10]，通過數(shù)值算例的形式對此界面單元的計算效果進行了正確性和有效性驗證。

1 考慮斷裂-接觸的界面單元

1.1 三角形界面單元

本文對于界面相互作用的模擬融合了粘結(jié)、斷裂和接觸這三種行為，為了區(qū)別于僅考慮接觸的單元類型，本文采用了界面單元這一描述方式。本文將發(fā)生相互作用的邊界面離散化為一系列的三質(zhì)點三角形界面單元，根據(jù)可能發(fā)生的作用情況設(shè)置為目標(biāo)界面單元與從界面單元。每個界面單元將其面積加權(quán)均分給其質(zhì)點，即成為每個界面單元質(zhì)點的代表面積。圖1展示了目標(biāo)界面單元、從界面單元以及界面質(zhì)點的幾何關(guān)系。

圖1 有限質(zhì)點法界面單元Fig. 1 FPM interface element

區(qū)別于其他有限質(zhì)點法單元，界面單元質(zhì)點并不具有獨立的質(zhì)量，而是通過與殼、四面體實體或六面體實體等其他有限質(zhì)點法單元共用質(zhì)點而為離散模型的邊界面賦予界面作用模擬的能力。界面作用限定于從界面單元的質(zhì)點(后文中將簡稱為從界面質(zhì)點)與目標(biāo)界面單元之間。每個時刻針對每一個從界面質(zhì)點搜索與其實際發(fā)生作用的目標(biāo)界面單元，生成一個個點-面形式的界面作用對，作用對的點-單元之間產(chǎn)生法向或切向的作用力，實現(xiàn)粘結(jié)、開裂、接觸等行為的模擬。

出于罰函數(shù)法對顯式數(shù)值算法的普遍通用性，本文選用罰函數(shù)法作為有限質(zhì)點法界面模擬力學(xué)描述的基本形式，并在此基礎(chǔ)上添加基于點-面自由度耦合的界面粘結(jié)模擬與基于內(nèi)聚力模型的斷裂模擬。界面作用對中的從界面質(zhì)點p與目標(biāo)界面單元et之間的界面力學(xué)模型如圖2所示。如果考慮界面粘結(jié)的影響，質(zhì)點p與粘結(jié)點pc保持共同運動，以加速度相同來計算粘結(jié)應(yīng)力。當(dāng)由于粘結(jié)所產(chǎn)生的界面等效應(yīng)力超出了給定的臨界值σcr時，則判定該界面作用對的粘結(jié)失效，轉(zhuǎn)為內(nèi)聚力開裂模型。裂縫完全展開后，繼而轉(zhuǎn)為接觸模型。該接觸模型等效于在作用發(fā)生點pc添加了沿法向界面彈簧Kn和切向非線性彈簧Ks，如圖2所示。

圖2 有限質(zhì)點法界面單元力學(xué)模型Fig. 2 Mechanical model of the FPM interface element

本節(jié)將詳細(xì)介紹界面作用力的并行化計算過程，包括作用對偵察、界面狀態(tài)判斷以及作用力計算這幾個步驟。

1.2 界面作用對偵察

本文假定每一個參與作用對偵察的從界面質(zhì)點在某一時刻僅與一個目標(biāo)界面單元發(fā)生作用，因此作用對偵查的目標(biāo)即為針對每一個從界面質(zhì)點確定唯一的點-面作用對。在所有可能發(fā)生相互作用的點-面對中找到真正發(fā)生作用的點-面對是任何偵察算法中最重要且最耗時的一步。為了進一步加速偵察的效率，本文在經(jīng)由GPU并行的有限質(zhì)點法架構(gòu)[10]中，進行并行化的作用對偵察。

在每一個時間步，首先對偵察范圍進行區(qū)域分割，進行一次全局搜索。區(qū)域分割的目的是為了將龐大的偵察范圍拆解為大量的可同時并獨立地進行偵察的較小空間，方便并行執(zhí)行以提高效率。參考在散粒系統(tǒng)中常用的“桶排序(Bucket Sort)”算法的思想，整個接觸偵察區(qū)域根據(jù)坐標(biāo)范圍被劃分為若干立方體網(wǎng)格，即若干個“桶”容器。若給定單元格邊長l，可以計算在全局坐標(biāo)系三個方向下的單元格個數(shù)Ni,j,k。以x方向上的網(wǎng)格數(shù)Ni為例：

式中：xmax和xmin分別為整個界面單元分布區(qū)域的包圍盒x坐標(biāo)的最大值和最小值，此包圍盒尺寸在邊界位置宜適當(dāng)放大； f loor()為向下取整函數(shù)。這樣，整個偵察區(qū)域即被分割為Ni×Nj×Nk個包圍盒單元格，被分割區(qū)域內(nèi)的任何一點均隸屬于唯一的包圍盒單元格中。

界面作用對偵察的過程分為兩步。首先，將參與接觸搜索的所有界面質(zhì)點按照當(dāng)前時刻的位置坐標(biāo)裝入劃分好的包圍盒單元格中，如圖3所示。在某時刻的搜索開始階段，對于任意界面質(zhì)點p，根據(jù)其在上一時刻結(jié)束后的位置坐標(biāo)可以確定其當(dāng)前時刻的包圍盒單元格歸屬Bp(i,j,k)，(i,j,k)表示該單元格在所有單元格中的位置索引。以x方向上的索引i為例：

圖3 界面作用對偵察：全局搜索Fig. 3 Global search of interaction pairs

式中，xp為質(zhì)點p在上一時刻結(jié)束后坐標(biāo)。此過程以參與偵察的質(zhì)點數(shù)為并行度，每一個GPU線程負(fù)責(zé)確定一個質(zhì)點的包圍盒單元格歸屬。

下一步，對于任意目標(biāo)界面單元e，包絡(luò)其三角形幾何形狀的包圍盒單元格集合Be={B(i,j,k)}即為該單元的包圍盒單元格集，被放入這個集合中的所有從界面質(zhì)點即為潛在作用質(zhì)點。這一過程以參與偵察的目標(biāo)界面單元數(shù)為并行度，即每一個GPU線程負(fù)責(zé)確定一個目標(biāo)界面單元的包圍盒單元格集并將這些網(wǎng)格中的質(zhì)點索引放入各個單元獨有的數(shù)據(jù)結(jié)構(gòu)中進行存儲。至此，界面作用對的偵察從全局范圍縮小至由各個目標(biāo)界面單元的包圍盒單元格集所控制的局部范圍內(nèi)，如圖4所示。

圖4 界面作用對偵察：局部搜索Fig. 4 Local search of interaction pairs

經(jīng)過了全局和局部搜索的篩查，對于每個目標(biāo)界面單元，需對其包圍盒單元格集中的從界面質(zhì)點進行逐一判斷以最終確定實際發(fā)生界面作用的點-面作用對。如圖5所示，為了保證界面作用方向的統(tǒng)一性，對任意一個界面單元e，其外法向ne定義三個質(zhì)點逆時針圍繞為正向，而任意一個界面質(zhì)點p的外法向np取為與之相連的所有界面單元外法向的加權(quán)平均值：

圖5 Inside-Outside算法示意圖Fig. 5 Diagram of the Inside-Outside algorithm

式中：v12和v23代表單元e三個質(zhì)點按照逆時針方向相連的邊向量；Ae表示單元e的面積。利用式(3)可以唯一確定質(zhì)點的法線方向，界面作用只發(fā)生在外法線方向相向的從界面質(zhì)點與目標(biāo)界面單元之間。

常規(guī)的點-面投影算法可以判斷點的投影是否落在三角形內(nèi)部以進一步判斷點面是否發(fā)生相互作用，但在目標(biāo)界面外凸或內(nèi)凹時可能會發(fā)生重判或漏判。本文參考“Inside-Outside”判斷法[31]來唯一確定每個從界面質(zhì)點的目標(biāo)界面單元歸屬。

如圖5所示的質(zhì)點p，對于單元e的三邊分別判斷其內(nèi)外狀態(tài)，計算：

式中，np為質(zhì)點p的外法向。若dij≤0則認(rèn)為質(zhì)點p對于邊ij處于In狀態(tài)，反之為Out狀態(tài)。如果質(zhì)點p對于單元e的三邊同時處于In或Out狀態(tài)，則判定質(zhì)點p歸屬于單元e。

在確定了質(zhì)點的單元歸屬后，采用點-面投影算法計算質(zhì)點在單元內(nèi)的投影位置?？紤]到三角形單元的形函數(shù)即為面積坐標(biāo)，投影點的位置使用面積坐標(biāo)進行表示：

式中，x1、x2、x3、xp分別為單元e的三個質(zhì)點以及質(zhì)點p的位置坐標(biāo)。由式(5) 即可求得投影點pc在單元e中的面積坐標(biāo)值Li。

得到了單元內(nèi)投影點位置后，可以計算質(zhì)點p對單元e的法向嵌入深度(或間隙)gn，gn取發(fā)生嵌入為負(fù)，產(chǎn)生間隙為正：

式中，npc為單元e內(nèi)投影點(即界面作用發(fā)生點)的外法向，可以取單元e的外法向。為了將目標(biāo)界面平滑化，本文將npc按照單元e三個質(zhì)點的質(zhì)點外法向按形函數(shù)(即面積坐標(biāo)Li)進行單元內(nèi)插：

這樣處理可以使得界面單元內(nèi)任意一點法線在整個目標(biāo)界面內(nèi)保持光滑連續(xù)，有效減小相鄰界面單元作用力的突變。

上述過程以并行方式進行，每個線程負(fù)責(zé)對一個目標(biāo)界面單元包圍盒單元格集中的所有從質(zhì)點依次進行循環(huán)，計算此從質(zhì)點是否歸屬于該目標(biāo)界面單元，如果確定歸屬則進而計算二者之間的法向嵌入深度gn。如果gn≤0，則發(fā)生相互作用，從而使得每個從質(zhì)點的點-面作用對得以唯一確定。

下面討論特殊情況。一個從質(zhì)點可能同時處于若干目標(biāo)界面單元的局部偵察范圍。由于并行實現(xiàn)時各個線程之間具有同時性與獨立性，不同線程針對不同單元的局部偵察過程可能同時對同一個從質(zhì)點進行歸屬計算。如果某從質(zhì)點的外法向恰好指向目標(biāo)界面單元的質(zhì)點或兩個相鄰目標(biāo)界面單元的公共邊，這種情況在界面初始粘結(jié)時十分常見，此時“Inside-Outside”算法會出現(xiàn)重判。這種情況可以借由并行算法中的原子操作進行處理。

當(dāng)某個目標(biāo)界面單元判定某從質(zhì)點與之發(fā)生接觸后，會將嵌入深度、單元索引、投影點面積坐標(biāo)、投影點外法向?qū)懭朐搹馁|(zhì)點獨有的數(shù)據(jù)結(jié)構(gòu)中。如果發(fā)生了重判，則會出現(xiàn)同一時刻不同線程(即不同目標(biāo)界面單元)向同一個內(nèi)存地址(即某從質(zhì)點記錄接觸信息的數(shù)據(jù)結(jié)構(gòu))進行寫入的操作。本文使用原子交換操作(AtomicExchange)[32]保證同一時刻只有一組接觸信息成功寫入該質(zhì)點的數(shù)據(jù)結(jié)構(gòu)中。此過程如圖6所示。

圖6 并行化的界面作用對偵察Fig. 6 Parallel procedures for searching interaction pairs

經(jīng)過“Inside-Outside”算法以及并行原子寫操作的雙重保障，每個從質(zhì)點的唯一點-面作用對得以最終確定，可以進行后續(xù)的作用力計算。

1.3 界面粘結(jié)狀態(tài)判斷

經(jīng)過了界面作用對的偵察與最終確定，對于任意從界面質(zhì)點p，如果它對單元e發(fā)生了嵌入，則p?e作用對成立，此單元標(biāo)記為et。由于一個時刻界面作用對對于每個從界面質(zhì)點是唯一的，界面作用力的計算將以界面質(zhì)點為并行度執(zhí)行，即每個線程負(fù)責(zé)計算一個界面作用對的作用力。

每個從界面質(zhì)點均存在一個粘結(jié)狀態(tài)指示量sbond，并在時間循環(huán)開始前根據(jù)界面設(shè)置情況設(shè)置為sbond=1 (初始考慮粘結(jié))與sbond=0(初始不考慮粘結(jié))。對于考慮粘結(jié)的情況，只需要首個迭代步中進行界面作用對的偵察與確定即可，這是因為粘結(jié)的存在會使得點-面相對位置保持不變。

對于界面作用對中的從界面質(zhì)點p與界面作用發(fā)生點pc，結(jié)合有限質(zhì)點法的點公式，有：

式中：下標(biāo)p和pc分別為從界面質(zhì)點與界面作用發(fā)生點；力向量的上標(biāo) e xt 、 int 和 b ond分別為質(zhì)點的外力、單元變形所產(chǎn)生的單元內(nèi)力以及界面作用的粘結(jié)力。界面作用發(fā)生點pc雖然并非在網(wǎng)格離散時生成的質(zhì)點，但可將其視為et上為了抵御質(zhì)點p的作用所產(chǎn)生的防御點，其質(zhì)量、力等物理量可由三角形界面單元的形函數(shù)得到：

式 中：Li為pc在et中 的面積 坐 標(biāo)；mi和Fi分別 為et的三個質(zhì)點的質(zhì)量與合力。

粘結(jié)作用力的計算原理為從界面質(zhì)點p與界面作用發(fā)生點pc的加速度相同，輔以粘結(jié)作用力對于兩質(zhì)點互為相互作用力的關(guān)系：

1.4 界面開裂作用力計算

當(dāng)界面等效應(yīng)力 σeq超出了界面的臨界應(yīng)力σcr時，界面開始發(fā)生斷裂。本文使用外加內(nèi)聚力模型[33]對裂縫開展過程進行模擬，該模型中變形與應(yīng)力的本構(gòu)關(guān)系如圖7所示。

圖7 外加內(nèi)聚力模型Fig. 7 Extrinsic cohesive model

建立接觸點pc處的局部坐標(biāo)系 (ξ,η,npc)，則p與pc的相對剪切變形gs為：

式 中：xp和xpc分別 為p與pc的 位移；xpc同樣 在et內(nèi)根據(jù)形函數(shù)插值求得。

類似于界面等效應(yīng)力的定義，這里使用嵌入深度gn與相對剪切變形gs定義界面等效變形geq：

式中， max()函數(shù)是為了保證只有界面間隙計入等效變形而界面嵌入不計入。根據(jù)線性外加內(nèi)聚力模型，開裂所產(chǎn)生的內(nèi)聚應(yīng)力 σcrack與等效變形之間geq的關(guān)系滿足線性關(guān)系：

式中，gc為等效變形的臨界值，可由斷裂釋放能Ec計算，gc=2Ec/σcr。

式中，Ap為質(zhì)點p的代表區(qū)域面積。當(dāng)geq＜gm時則視為卸載，此時用gm替換式(16)中的geq即可。當(dāng)geq＞gc時，界面完全斷裂，記scrack=0，界面內(nèi)聚力為0，徹底轉(zhuǎn)為接觸模型。

1.5 界面接觸作用力計算

式中：gn由式 (6) 計算得來；Kn為界面法向接觸剛度；Ap為質(zhì)點p的代表區(qū)域面積。

對于切向接觸力，其計算通常基于摩爾庫侖準(zhǔn)則的基本形式，等效于彈塑性狀態(tài)下的切向彈簧，屈服函數(shù)為：

式中，fsmax為界面可以提供的最大切向接觸力，可由界面粘聚力c和摩擦角φ計算而出：

由于常用的雙線性塑性形式的摩爾庫侖準(zhǔn)則在時間步長不合理或剛度設(shè)置不合理時可能會在相對剪切變形接近轉(zhuǎn)向時難以收斂，本文使用反正切函數(shù)型的平滑摩擦力模型計算切向接觸力，在每一個途徑單元內(nèi)部按照如下兩個步驟進行：

1) 計算t時刻 (ξ,η) 平面內(nèi)p與pc的相對剪切變形增量 ?gs：

式中：vs=?gs/?t為相對剪切變形速度；C為無量綱的平滑參數(shù)，一般取0.1～0.001?？梢钥闯?，式 (21)實質(zhì)為變剛度的摩爾庫侖準(zhǔn)則。

將以上所描述的界面作用對的偵察、確定與作用力計算的計算過程嵌入并行有限質(zhì)點法的計算框架[10]中，即可到如圖8所示的界面單元計算流程。

2 數(shù)值算例

本課題組在先前的研究中設(shè)計研發(fā)了采用GPU加速的有限質(zhì)點法并行求解系統(tǒng)，并形成了有限質(zhì)點法通用計算平臺FPM。本節(jié)的算例均使用此平臺進行加速計算。關(guān)于有限質(zhì)點法的并行化以及FPM平臺的相關(guān)信息，可參考文獻[10]。

2.1 基礎(chǔ)-土體接觸

為了驗證本文的界面單元在處理接觸問題上的效果，本例考慮Qian等[34]使用的一個簡化的基礎(chǔ)-土體的接觸模型。如圖9所示，基礎(chǔ)為6 m×6 m×0.6 m的立方體；上部結(jié)構(gòu)的作用被簡化為作用在基礎(chǔ)上的法向均布荷載1 MPa；下部的土層簡化為10 m×10 m×4 m的立方體，四周側(cè)向約束，底部豎向約束。相互接觸的兩層三角形界面單元分別被設(shè)置在基礎(chǔ)下表面和土層上表面，土層上表面設(shè)置為目標(biāo)界面單元。

基礎(chǔ)與土層采用相同的彈性材料，其彈性模量為100 MPa，泊松比為0.2。接觸界面摩擦角為5.71°，即等效的摩擦系數(shù)為0.1。本例僅考慮接觸作用，不考慮界面粘結(jié)與斷裂。作為目標(biāo)界面的土層上表面的接觸單元網(wǎng)格尺寸為0.4 m，作為從界面的基礎(chǔ)下表面的網(wǎng)格為0.2 m，對應(yīng)的臨界步長約為0.2 ms。使用有限質(zhì)點法對本例進行計算時，計算步長取為0.1 ms，計算總時長1 s，均布荷載緩慢斜坡加載至基礎(chǔ)上表面。為了加速收斂速度，取虛擬的全局質(zhì)量阻尼系數(shù)100。作為對比，本例使用通用有限元軟件的罰函數(shù)法接觸模型，接觸剛度取558 MPa，其余計算參數(shù)均保持一致。為了驗證界面單元在不同網(wǎng)格形狀上的適應(yīng)性，在FPM平臺中分別使用四面體實體單元和六面體實體單元進行建模分析。

使用FPM平臺計算所得的結(jié)構(gòu)位移云圖如圖10所示，接觸面所在平面與Y=0平面交界線位置的界面單元的法向接觸應(yīng)力分布如圖11(a)和圖11(b)所示，并與Qian等[34]的結(jié)果進行了對比?？梢钥吹?，本文的有限質(zhì)點法界面單元的模擬結(jié)果呈現(xiàn)出了彈性力學(xué)解中的中心應(yīng)力均勻兩端應(yīng)力集中的狀態(tài)，且相較于有限元結(jié)果，本文中經(jīng)由法向量插值的曲面連續(xù)化處理使得法向接觸應(yīng)力的分布更加均勻。另外，Qian等[34]在模擬中使用了接觸面的曲面光滑算法RPIM，目的是讓接觸面上的接觸力分布更為平滑連續(xù)。而本文采用了基于形函數(shù)插值的法向量連續(xù)化處理，也同樣產(chǎn)生了一定的應(yīng)力平滑效果，雖然平滑效果有限，但計算代價極低。本例說明了本文實現(xiàn)的有限質(zhì)點法界面單元在處理接觸問題上的精度和有效性。

圖10 基礎(chǔ)-土體接觸：四面體網(wǎng)格模型變形云圖Fig. 10 Interface contact analysis between foundation and soil: displacement contour of the model with tetrahedron mesh

在四面體模型中，實體單元表面為三角形，用一個界面單元進行覆蓋；六面體單元表面為四邊形，使用兩個三角形界面單元進行覆蓋。兩種網(wǎng)格劃分的模型使用FPM平臺計算所得的中心面接觸應(yīng)力分布如圖11(c)所示。除了四面體網(wǎng)格由于網(wǎng)格走向不對稱所造成的應(yīng)力分布不完全對稱的情況之外，兩種界面單元覆蓋方式下的接觸應(yīng)力分布具有極高的吻合度，證明了本文的三角形界面單元具有普遍的幾何適應(yīng)性，能夠在建模時以覆蓋在其他三維單元表面上的形式賦予其表面界面作用分析的能力。

圖11 基礎(chǔ)-土體接觸：Y=0平面界面單元接觸應(yīng)力分布Fig. 11 Interface contact analysis between foundation and soil: contact stress distribution on the Y=0 plane

2.2 磚塊粘合界面剪切

本例對采用文獻[22]所使用的Beyer等[35]進行的磚塊粘合界面剪切試驗進行數(shù)值模擬，用于測試本文界面單元在處理界面粘結(jié)后剪切開裂并耦合接觸時的模擬效果。原始試驗裝置如圖12(a)所示，三塊磚塊之間用砂漿粘結(jié)并列放置，兩側(cè)使用加載裝置進行水平向的施壓。中間磚塊頂端靠近界面的位置進行了豎向的位移約束，兩側(cè)磚塊在靠近接觸面的位置作用了豎向的移動支座，用于對界面施加剪切作用。

根據(jù)此試驗裝置，本文建立了簡化的有限質(zhì)點法數(shù)值模型，相關(guān)尺寸和邊界條件如圖12(b)所示。數(shù)值模型采用對稱的一半建立，磚塊使用四面體實體單元進行建模，中央半塊磚塊的頂部和左側(cè)施加對應(yīng)法向的固定約束，右側(cè)磚塊底部施加豎向的剪切位移約束，速度固定為25 mm/s，右側(cè)作用法向壓力P，分別取值為0.2 MPa、0.4 MPa和0.65 MPa。界面左右位置的兩個面覆蓋了界面單元，界面單元尺寸在3 mm～4 mm，磚塊的實體網(wǎng)格在界面位置與界面單元網(wǎng)格重合，邊界端網(wǎng)格尺寸適當(dāng)放大，如圖13所示，模型的臨界步長約為0.38 μs。

圖12 磚塊粘合界面剪切：試驗裝置圖與簡化的數(shù)值模型Fig. 12 Shear test of masonry wallettes: experiment setup and simplified numerical model

圖13 磚塊粘合界面剪切：FPM平臺中網(wǎng)格劃分Fig. 13 Shear test of masonry wallettes: the FPM mesh

根據(jù)材性試驗，磚塊彈性模量為14 GPa，泊松比為0.15，密度為940 kg/m3。界面斷裂釋放能為750 J/m2，界面粘結(jié)力的臨界值取為0.25 MPa，界面摩擦系數(shù)為0.77。根據(jù)文獻[24]所建議的接觸剛度的計算公式，可以得到近似的接觸剛度2590 GPa?？紤]到接觸位置過剛會導(dǎo)致臨界時間步長顯著減小，實際計算中減小兩個量級。計算中步長保持恒定為0.05 μs，通過保持位移加載速度的25 mm/s，模擬支座位移加大到10 mm的過程。在FPM平臺中計算時，需要設(shè)置兩個分析過程，首先僅對模型右側(cè)作用法向壓力P；待穩(wěn)定后再施加底部的位移荷載。由于本例為軸壓狀態(tài)下的剪切斷裂，求解過程中設(shè)置接觸面法向始終保持粘結(jié)，允許切向剪切斷裂。圖14展示了不同條件下FPM平臺所計算得出的界面剪切應(yīng)力-剪切變形曲線，其中應(yīng)力和變形均提取自作用位移荷載的第二個分析過程，應(yīng)力值取界面中心0.15 m高度范圍內(nèi)界面單元質(zhì)點應(yīng)力的平均值，此過程界面法向壓力保持恒定。

圖14 磚塊粘合界面剪切：界面剪切應(yīng)力-剪切位移曲線Fig. 14 Shear test of masonry wallettes: shear stress plotted against shear displacement

圖14(a)～圖14(c)展示了當(dāng)法向壓應(yīng)力分別保持為0.2 MPa、0.4 MPa和0.65 MPa時，隨著剪切位移增大，界面剪切應(yīng)力的變化曲線，并與Beyer等[35]的原始試驗數(shù)據(jù)進行了比較。有限質(zhì)點法模擬結(jié)果與試驗結(jié)果呈現(xiàn)了相同的變化趨勢。從模擬結(jié)果可以看出，在法向壓應(yīng)力保持不變的情況下，界面剪切變形從0開始增長，界面間始終存在摩擦并很快達到最大摩擦力，保持滑動，切向內(nèi)聚力開始發(fā)揮作用；隨著剪切變形繼續(xù)增長，切向的內(nèi)聚力逐漸減小，直至內(nèi)聚力減小為零，界面間僅存在接觸造成的滑動摩擦的作用，剪切應(yīng)力大小逐漸穩(wěn)定。因此，本文所實現(xiàn)的界面單元在受壓剪切的情況下的本構(gòu)曲線形式實質(zhì)上是外加內(nèi)聚力曲線和切向摩擦力曲線這二者的疊加。對于模擬中限定的斷裂釋放能以及臨界應(yīng)力值，可以計算得到界面內(nèi)聚力在剪切位移為6 mm左右的位置失效，也與模擬結(jié)果相吻合。試驗結(jié)果曲線更為平滑，這主要是由于實際情況中摩擦力和界面內(nèi)聚力的耦合更加連續(xù)；而且本文采用了線性外加內(nèi)聚力模型，較真實情況有一定的簡化。

諸多文獻同樣使用此試驗進行了數(shù)值模擬，這里選擇了Snoozi和Molinari[36]、Spring和Paulino[24]以及Baek和Park[22]三者的模擬結(jié)果與本文的結(jié)果進行比較，如圖14(d)所示。Snoozi和Molinari[36]采用了與本文相同的線性外加內(nèi)聚力模型，但Snoozi和Molinari[36]采用了指數(shù)形式的摩擦力規(guī)則化函數(shù)以引入界面開裂損傷對摩擦力的影響。Spring和Paulino[24]以及Baek和Park[22]使用了基于勢能的Park-Paulino-Roesler內(nèi)聚力模型的變化形式，更符合真實物理過程，但計算相對復(fù)雜。需要說明的是，為了匹配試驗結(jié)果，Snoozi和Molinari[36]、Spring和Paulino[24]并未采取試驗得出的界面粘聚力臨界值0.25 MPa，Snoozi和Molinari[36]取0.3 MPa，而Spring和Paulino[24]取0.45 MPa。

圖14(e)～圖14(f)進一步展示了剪切應(yīng)力-剪切位移曲線隨著界面臨界應(yīng)力 σcr和斷裂能Ec取值不同而變化的情況。當(dāng)界面臨界應(yīng)力保持不變時，斷裂能越大，界面開裂過程的內(nèi)聚力作用距離越長，界面完全斷裂所需要變形也就越大；當(dāng)斷裂能保持不變時，界面臨界應(yīng)力越大，剪切應(yīng)力峰值越大，界面完全開裂所需要的剪切變形就越小，越快進入完全接觸摩擦的平臺段。

總地來說，本文基于有限質(zhì)點法所開發(fā)的界面單元在界面粘結(jié)-斷裂-接觸耦合的復(fù)雜行為時具有足夠的精度以及有效性。

3 結(jié)論

基于有限質(zhì)點法基本理論并依托有限質(zhì)點法通用計算平臺FPM，本文開發(fā)了受用于并行化的有限質(zhì)點法計算體系下的三角形界面單元，用于模擬界面之間的粘結(jié)、斷裂、接觸等復(fù)雜的界面作用行為。該單元采用并行化的作用對偵察算法，實現(xiàn)了基于等加速度原理的粘結(jié)作用、基于內(nèi)聚力模型的斷裂作用，以及基于罰函數(shù)法的接觸作用這三種作用下的狀態(tài)判據(jù)和作用力計算。數(shù)值算例證明本文所開發(fā)的界面單元在處理普通接觸問題以及粘結(jié)-斷裂-接觸耦合的復(fù)雜界面作用行為方面是正確并有效的。利用本文提出的界面單元，可對巖土和復(fù)合結(jié)構(gòu)等領(lǐng)域的實際工程常見的復(fù)雜界面行為進行較為高效的模擬，也為實際工程中復(fù)雜界面行為的精細(xì)化大規(guī)模數(shù)值分析奠定了基礎(chǔ)。