李廣凱，劉傳東，肖仁軍，馬懷發(fā)

（1.山東泰山抽水蓄能電站有限責(zé)任公司，山東 泰安 271000；2.中國水利水電科學(xué)研究院 流域水循環(huán)模擬與調(diào)控國家重點(diǎn)實(shí)驗(yàn)室，北京 100038；3.中國水利水電科學(xué)研究院 工程抗震研究中心，北京 100048）

1 研究背景

接觸問題是工程中普遍存在的力學(xué)問題，接觸問題的特點(diǎn)和難點(diǎn)是接觸邊界和接觸力的未知性。由于接觸界面的區(qū)域大小、位置以及接觸狀態(tài)都是未知的、隨時(shí)間變化的，因此接觸問題表現(xiàn)出顯著的非線性特征。接觸問題的非線性決定了接觸分析過程中需要經(jīng)常插入接觸界面的搜尋判定。有限元法為分析接觸問題提供了有效的工具，在進(jìn)行有限元求解接觸問題時(shí)，接觸體被離散為空間單元集，接觸界面的相互作用被轉(zhuǎn)化為離散的單元表面之間的作用或者離散的節(jié)點(diǎn)與單元表面的作用。接觸問題的動(dòng)力分析主要包括選取時(shí)間積分方案和接觸算法，其中接觸算法包含接觸搜索算法和接觸力算法。

接觸搜索的目的是確定整個(gè)系統(tǒng)中有哪些部位發(fā)生了接觸或者有哪些原已接觸的部分發(fā)生了滑移或脫離。在接觸動(dòng)力分析中，接觸搜索占較大的計(jì)算量，因此接觸搜索算法的計(jì)算效率至關(guān)重要。由于接觸搜索是接觸力計(jì)算的基礎(chǔ)，所以其計(jì)算精度也非常關(guān)鍵。在一般情況下接觸搜索需要先進(jìn)行全局搜索，再進(jìn)行局部搜索。通過全局搜索先粗略找到圍繞特殊點(diǎn)所有潛在可能的接觸單元面。全局算法有主從面算法［1］、級(jí)域算法［2-3］、單曲面算法［4］、位碼算法［5］和LC-Grid法［6］等。主從面算法將相互接觸的兩個(gè)邊界面分別指定為主面和從面，主面上的節(jié)點(diǎn)定義為主節(jié)點(diǎn)，從面上的節(jié)點(diǎn)定義為從節(jié)點(diǎn)，但在單一曲面發(fā)生嚴(yán)重形變時(shí)，不能處理自接觸搜索；級(jí)域法［2-3］涉及不同級(jí)別接觸元素的定義和管理，在接觸搜索過程中需要逐級(jí)進(jìn)行，數(shù)據(jù)管理復(fù)雜，不便于并行計(jì)算，并且只在一定的條件下有效。單曲面算法［4］使用了子域排序算法來尋找接觸面，通過三層嵌套排序來完成三維空間的搜索過程，優(yōu)點(diǎn)是能夠處理自接觸問題，但該算法編程復(fù)雜，運(yùn)算量大，可并行性差；位碼算法采用了位碼構(gòu)造思想，將三維空間中分類和搜索的問題轉(zhuǎn)化為一維數(shù)組分類排序，進(jìn)行搜索處理，既不需要事先指定接觸區(qū)域，又避免了三維嵌套搜索，便于進(jìn)行大規(guī)模計(jì)算，是計(jì)算效率較好的接觸搜索算法；LC-Grid類似于位碼算法，其計(jì)算時(shí)間和存儲(chǔ)要求與接觸面數(shù)量成線性比例，但計(jì)算效率依賴于接觸空間分解尺度。

局部搜索目的是計(jì)算節(jié)點(diǎn)到單元面的精確距離，找出接觸投影點(diǎn)所在單元面的相對(duì)位置關(guān)系，從而確定接觸節(jié)點(diǎn)與單元面間的接觸狀態(tài)。局部搜索算法主要有點(diǎn)面算法［1，4］、小球算法［7］、基于光滑曲面的搜索算法［8-9］以及內(nèi)外算法［10］等。點(diǎn)面算法在計(jì)算“點(diǎn)”在單元面上的投影點(diǎn)時(shí)，需迭代求解點(diǎn)面投影的偏微分方程，有時(shí)會(huì)出現(xiàn)迭代穩(wěn)定性問題；小球算法是一種幾何近似方法，該方法將小球嵌入表面單元，僅對(duì)小球施加不可侵入條件，主要用于滑動(dòng)和摩擦力不是至關(guān)重要的問題；基于光滑曲線（曲面）的搜索算法采用樣條函數(shù)逼近接觸表面更加準(zhǔn)確地描述接觸邊界，從而使接觸搜索甚至接觸力的計(jì)算精度提高，但此技術(shù)在離散模型中引入了冗余信息，加大了數(shù)據(jù)處理的難度；內(nèi)外算法僅僅是邏輯操作過程，簡單穩(wěn)定，計(jì)算精度高。但這些方法幾乎都出現(xiàn)搜索盲區(qū)問題。另外，除內(nèi)外算法外，其它算法都需要通過點(diǎn)到投影面的非線性投影方程，確定節(jié)點(diǎn)和潛在接觸面的距離。盡管小球算法可以回避求解非線性方程，但會(huì)降低計(jì)算的準(zhǔn)確性。

在完成接觸搜索之后，接觸力求解通常采用拉格朗日乘子法或罰函數(shù)法［11］引入附加條件構(gòu)造修正基本方程泛函的方法。其中罰函數(shù)法需要根據(jù)經(jīng)驗(yàn)選擇接觸剛度參數(shù)即彈簧常數(shù)，而拉格朗日乘子法對(duì)接縫的處理更真實(shí)、更自然。接觸面的本構(gòu)模型按照接觸面的性質(zhì)，有光滑接觸模型和摩擦接觸模型，如果接觸面是絕對(duì)光滑的或摩擦可以忽略，則采用光滑接觸模型，但在處理壩體伸縮縫和地基夾層接觸時(shí)，一般采用摩擦接觸模型，并按庫倫摩爾準(zhǔn)則考慮摩擦條件。

強(qiáng)震作用下的高混凝土壩系統(tǒng)災(zāi)變過程伴隨著材料和接觸復(fù)合非線性問題，對(duì)其進(jìn)行全過程精細(xì)化數(shù)值模擬和全面深入的抗震安全評(píng)價(jià)，需要求解未知量高達(dá)百萬甚至千萬級(jí)，盡管各種商業(yè)軟件在求解一般常見問題時(shí)可以顯示出所謂的強(qiáng)大功能，但在計(jì)算方法和計(jì)算實(shí)施方案上很難適用于解決高混凝土壩復(fù)雜工況下地震響應(yīng)分析所遇到的復(fù)雜而又特殊的問題，因此，高性能并行計(jì)算是解決該問題的必然選擇。

本文擬在研究現(xiàn)有接觸算法的基礎(chǔ)上，提出多體接觸問題快捷高效的面面接觸算法。這種算法既可解決搜索盲區(qū)問題，又能回避求解非線性投影方程。本文算法采用拉格朗日乘子法動(dòng)力接觸方程的增量形式求解接觸力，另外，為了便于并行計(jì)算及其程序編制，本文還根據(jù)接觸面分布特點(diǎn)進(jìn)行預(yù)分區(qū)處理，并基于提出的多體接觸問題面面接觸算法，開發(fā)實(shí)現(xiàn)其算法的FORTRAN 源代碼程序，為后續(xù)實(shí)現(xiàn)高混凝土壩和巖體高邊坡靜動(dòng)力分析的并行計(jì)算提供支持。

2 搜索算法分析

2.1 全局搜索算法本文將主從面算法與位碼算法相結(jié)合，實(shí)現(xiàn)接觸面的全局搜索。利用主從面算法［1］在建立幾何模型時(shí)，根據(jù)可能接觸邊界的分布特點(diǎn)將可能接觸邊界面分別指定為主面和從面，并將其定義為“接觸面對(duì)”，再利用類似于位碼算法［5］將包含模型接觸表面的三維空間劃分成若干個(gè)立方格，并對(duì)這些立方格在三個(gè)緯度方向進(jìn)行整數(shù)編碼，通過這些位置號(hào)碼從三維映射到一維，基于一維排序和搜索的算法來實(shí)現(xiàn)對(duì)區(qū)域內(nèi)接觸節(jié)點(diǎn)的檢測判別。在有限元剖分網(wǎng)格時(shí)主面單元上的節(jié)點(diǎn)定義為主節(jié)點(diǎn)，從面單元上的節(jié)點(diǎn)定義為從節(jié)點(diǎn)。通過搜索與主面接觸的從節(jié)點(diǎn)，即找出由從節(jié)點(diǎn)與主面上的單元面所構(gòu)成的“接觸點(diǎn)面對(duì)”。

本文構(gòu)建了描述“接觸面對(duì)”的主面單元和從面單元的共享實(shí)常數(shù)組：主面單元用正整數(shù)表示，從面單元用其對(duì)應(yīng)的負(fù)整數(shù)表示，當(dāng)絕對(duì)值相同時(shí)主面單元和從面單元構(gòu)成“接觸面對(duì)”；不同的“接觸面對(duì)”通過不同的實(shí)常數(shù)定義；一組實(shí)常數(shù)可以對(duì)應(yīng)多個(gè)邊界面。

以下是對(duì)一組“接觸面對(duì)”的全局搜索步驟：（1）搜索預(yù)處理：統(tǒng)計(jì)“接觸面對(duì)”單元面上的主節(jié)點(diǎn)和從節(jié)點(diǎn)個(gè)數(shù)，計(jì)算各單元面外法向量及各節(jié)點(diǎn)外法向量。其中節(jié)點(diǎn)外法向量由包含該節(jié)點(diǎn)的單元法向量取平均得到。（2）根據(jù)“接觸面對(duì)”的位置關(guān)系，構(gòu)造包含接觸單元面的最小長方體。在接觸問題的有限元計(jì)算過程中，全局搜索是相對(duì)耗時(shí)較大的計(jì)算環(huán)節(jié)，為了避免在每個(gè)時(shí)間步都進(jìn)行全局搜索，適當(dāng)擴(kuò)大該長方體的區(qū)域，確定接觸搜索范圍。（3）確定立方格的網(wǎng)格尺寸，立方格尺寸與接觸面中平均的單元尺寸相接近或可以稍大一些。（4）建立典型的立方格結(jié)構(gòu)，其三個(gè)坐標(biāo)（x，y，z）方向的網(wǎng)格尺寸一致。按先x方向，再y方向，最后z方向的順序?qū)α⒎礁襁M(jìn)行編號(hào)。（5）對(duì)當(dāng)前“接觸面對(duì)”的所有主從節(jié)點(diǎn)循環(huán)，根據(jù)節(jié)點(diǎn)坐標(biāo)確定其所在的立方格的編號(hào)［14］。（6）對(duì)立方格循環(huán)，在當(dāng)前立方格中，對(duì)每一個(gè)主節(jié)點(diǎn)，找距離與其最近的從節(jié)點(diǎn)，并記錄該從節(jié)點(diǎn)所對(duì)應(yīng)的接觸單元編號(hào)，建立點(diǎn)面接觸測試對(duì)。

基于以上全局搜索流程編制的程序框圖如圖1所示。

圖1 全局搜索程序框圖

2.2 局部搜索算法本文結(jié)合點(diǎn)面算法與內(nèi)外算法完成接觸局部搜索。由全局搜索獲得“接觸面對(duì)”的點(diǎn)面關(guān)系，采用內(nèi)外算法判斷“接觸面對(duì)”主面上的節(jié)點(diǎn)落入哪些從面接觸面單元內(nèi)，為了避免搜索盲區(qū)，本文采用了相關(guān)單元節(jié)點(diǎn)法向“平均向量”。如圖2所示，在接觸單元中，節(jié)點(diǎn)B的法向量n取單元AB的法向量na和單元BC的法向量nb平均值，即判斷主節(jié)點(diǎn)沿平均法線方向的投影點(diǎn)是在目標(biāo)面單元的內(nèi)部或外部。

下面給出一組“點(diǎn)面接觸測試對(duì)”的局部搜索步驟：（1）對(duì)測試對(duì)中的面單元數(shù)進(jìn)行循環(huán)，根據(jù)“點(diǎn)”法向量與單元面法向量，進(jìn)一步判斷是否為潛在接觸關(guān)系。（2）針對(duì)潛在接觸單元，計(jì)算“點(diǎn)”到單元面投影點(diǎn)位置及其貫入量。（3）由內(nèi)外算法［10］，判別“點(diǎn)”是落入到目標(biāo)單元面的有效區(qū)域的內(nèi)部還是外部。（4）如圖3所示，如“點(diǎn)”在目標(biāo)單元面內(nèi)部，則由三角形面積坐標(biāo)法，計(jì)算投影點(diǎn)在單元內(nèi)的局部坐標(biāo)，從而得到投影點(diǎn)對(duì)應(yīng)的單元插值形函數(shù)，進(jìn)而得出接觸關(guān)聯(lián)矩陣。

在圖3中，Vij表示由節(jié)點(diǎn)i指向節(jié)點(diǎn)j的矢量；矢量n1I=V12×V1I；n為從節(jié)點(diǎn)I的法向矢量；投影點(diǎn)x將四邊形分割成4 個(gè)小區(qū)域，對(duì)應(yīng)的面積為Δ1、Δ2、Δ3和Δ4，由此可以求得x的局部坐標(biāo)。

圖2 兩個(gè)單元的平均向量

圖3 點(diǎn)與單元面的關(guān)系

圖4 局部搜索的程序

圖4給出了基于局部搜索流程編制的程序框架。

綜上所述，本文采用主從面算法與位碼算法相結(jié)合完成接觸全局搜索，點(diǎn)面算法與內(nèi)外算法相結(jié)合完成接觸局部搜索。同時(shí)在執(zhí)行接觸搜索前預(yù)先考慮了接觸區(qū)域分塊劃分，使得本文算法可以方便地進(jìn)行高混凝土壩和巖體高邊坡靜動(dòng)力穩(wěn)定分析的并行計(jì)算處理。

3 接觸力計(jì)算

根據(jù)可能接觸邊界的分布特點(diǎn)，將計(jì)算域分解成不同的子區(qū)域?？紤]由N個(gè)子區(qū)域Ωi(i=1，2，…，N)組成的多接觸體系統(tǒng)。對(duì)于子區(qū)域i，其動(dòng)力平衡方程可表示為：

式中：ρ為質(zhì)量密度；c為阻尼系數(shù)；為速度分量；為加速度分量；fj為體力分量；σjk為應(yīng)力張量分量。

不同子區(qū)域間的可能接觸邊界條件為分離狀態(tài)、黏著狀態(tài)或滑動(dòng)狀態(tài)。

接觸力求解采用拉格朗日乘子法。對(duì)每一個(gè)編號(hào)為i的子區(qū)域，獨(dú)立進(jìn)行網(wǎng)格剖分和有限元離散后，在引進(jìn)拉格朗日乘子λ后，動(dòng)力學(xué)方程式（1）的空間離散形式為：

式中：i為子區(qū)域編號(hào)；M為質(zhì)量矩陣；C為阻尼矩陣；K為剛度矩陣；F為荷載項(xiàng)；Bi=[BinB it]為位移與接觸力相關(guān)矩陣，Bin、Bit分別為接觸面法向和切向相關(guān)矩陣；λ為拉格朗日乘子力（接觸力）。

采用Newmark直接積分法，將式（2）寫成增量形式：

式中：為廣義剛度；為廣義力增量；

可能接觸的邊界條件的離散形式為：

式中：δ為接觸面間隙，會(huì)隨著時(shí)間步變化，應(yīng)根據(jù)可能接觸邊界的幾何位置自動(dòng)判定和修正，法向接觸間隙切向接觸間隙

將式（3）和式（4）寫成整體的形式：

由式（5）消去未知量dU，整理出關(guān)于乘子力的柔度矩陣D及接觸荷載向量Q：

其中：

將接觸力的計(jì)算轉(zhuǎn)化為不等式方程組的求解：

采用擬高斯迭代法［12-13］求解t+Δt時(shí)刻位移增量與接觸力。在處理位移約束邊界時(shí)，將在剛度矩陣對(duì)應(yīng)約束邊界點(diǎn)所在行的對(duì)角線元素置1，并將與這個(gè)對(duì)角線元素相應(yīng)行和列的其它元素都賦為0，位移方程右端賦予相應(yīng)約束位移值。這樣做可以避免含有邊界約束的接觸面形成接觸力柔度矩陣奇異的問題。

4 算法驗(yàn)證

下面將采用兩個(gè)算例驗(yàn)證本文提出的面面接觸算法。算例選取的一般原則是模型簡單典型，將計(jì)算結(jié)果與其理論解對(duì)比，再者借助于商業(yè)軟件相應(yīng)的功能對(duì)本文算法程序進(jìn)行補(bǔ)充驗(yàn)證。第一個(gè)為塊體接觸模型，主要驗(yàn)證算法的穩(wěn)定性；第二算例為經(jīng)典的赫茲問題，將數(shù)值計(jì)算的接觸力和接觸半徑結(jié)果與理論解比較，并借助商業(yè)軟件考查接觸面上的變形分布是否與本文算法結(jié)果是否一致，以驗(yàn)證本文算法的正確性和精確性。

4.1 塊體接觸分析考慮由兩個(gè)塊體組成的系統(tǒng)，如圖5所示，塊體1和塊體2在x方向、y方向和z方向的長度均為100 m，塊體1在3個(gè)尺寸方向均為10等分；塊體2在豎向10等分，在橫向兩個(gè)方向20 等分。塊體1 和塊體2 的密度均取10 kg/m3，阻尼比取200；兩塊體彈性模量為10 GPa，泊松比取0，摩擦系數(shù)取0.6。塊體1底面法向約束，塊體2突然施加100 N/m3體積力，由本文算法程序計(jì)算其動(dòng)力響應(yīng)，時(shí)間步長分別取0.01、0.04、0.06和0.10 s。由本例給出的彈性參數(shù)可以得到精確解，即塊體上表面的豎向位移為-1.5×10-4m，兩塊體接觸面的豎向位移為-1.0×10-4m。

圖6給出了兩塊體在突加重力作用下穩(wěn)定后的變形云圖。圖6顯示，最大變形位于塊體2上的上表面，其豎向位移值為-1.5×10-4m。

圖7和圖8分別給出了位于兩塊體接觸面塊體2底面中點(diǎn)和上表面中點(diǎn)的位移時(shí)程曲線。圖7和圖8顯示，時(shí)間步長為0.01 s 時(shí)的變形很快穩(wěn)定下來，其次是步長為0.04 s，然后為0.06 s 和0.10 s。由于本文算法動(dòng)力方程的時(shí)間積分采用了無條件穩(wěn)定的Newmark 積分，因此，時(shí)間步長的大小僅會(huì)影響計(jì)算精度而不會(huì)影響其數(shù)值計(jì)算的穩(wěn)定性。由圖7和圖8可以看出，在計(jì)算穩(wěn)定后，接觸塊體2底面和上表面豎向位移值分別為-1.0×10-4和-1.5×10-4m，與精確解完全相同。

圖5 塊體接觸模型

圖6 接觸塊體豎向變形云圖（單位：m）

圖7 塊體2底面中點(diǎn)位移時(shí)程

圖8 塊體2上表面中點(diǎn)位移時(shí)程

4.2 圓柱體接觸分析兩平行接觸的圓柱體接觸體［15］，如圖9（a）所示，考慮其對(duì)稱性，取1/4圓柱進(jìn)行分析。設(shè)兩圓柱體的半徑均為5 mm，彈性模量為210 GPa，泊松比為0.3，認(rèn)為兩圓柱為絕對(duì)光滑接觸，即不計(jì)接觸面的摩擦力。計(jì)算模型取圓柱體長1 mm。按圖9（b）所示的在底部界面施加法向位移約束，上部施加的豎向位移為-0.03 mm。模型網(wǎng)格劃分為兩種區(qū)域，計(jì)算模型網(wǎng)格如圖9（c）所示，即接觸部分的網(wǎng)格密度較大，網(wǎng)格尺寸取為0.1 mm，大約為五分之一圓柱半徑區(qū)域，而其余區(qū)域的網(wǎng)格密度較小。密度取7.85×10-6kg/mm3，阻尼系數(shù)取200。時(shí)間步長取0.01 s，計(jì)算至?xí)r長3 s，即得到靜態(tài)穩(wěn)定解。

圖9 兩平行圓柱體的接觸

在本文中，圓柱彈性模量E=210GPa，圓柱半徑R1=R2=5mm，圓柱體長L=1mm，當(dāng)δ=0.03mm時(shí)，由接觸相對(duì)位移反求得沿圓柱軸向分布力P=1234.42N/m，b=0.185mm，σmax=4256.4MPa。

由本文算法程序得到的沿Y向的應(yīng)力分量云圖如圖10所示，最大接觸應(yīng)力發(fā)生在接觸中心點(diǎn)，其值為4136.6 MPa，與赫茲理論值為4256.4 MPa 相比，相對(duì)誤差僅為0.9%；計(jì)算得到接觸半徑為0.188 mm，赫茲理論解值為0.185 mm，相對(duì)誤差為0.6%，本文算法的計(jì)算結(jié)果與理論解幾乎完全吻合。該算例結(jié)果表明，本文的接觸算法是正確的和有效的。

另外，采用ANSYS軟件對(duì)兩接觸圓柱體變形進(jìn)行補(bǔ)充計(jì)算，得到如圖11所示的兩平行接觸的圓柱體接觸體豎向變形云圖。從圖11可以看出，本文算法程序與ANSYS 軟件計(jì)算得到的變形分布一致。對(duì)比表1列出的沿兩圓柱上下接觸面由左到右的豎向位移，二者上接觸面位移的最大誤差為0.71%，計(jì)算結(jié)果幾乎完全一致。

圖10 本文算法計(jì)算的豎向（Y）應(yīng)力分量 （單位：MPa）

圖11 兩平行圓柱體接觸體豎向（Y）變形云圖（單位: mm）

表1 圓柱接觸面位置位移計(jì)算結(jié)果對(duì)比

5 討論及結(jié)論

通常位移邊界條件的引入有3種方法：（1）合成總剛矩陣時(shí)劃去邊界位移約束對(duì)應(yīng)的各行各列元素，緊縮總剛度及荷載列陣，將約束的影響作用轉(zhuǎn)移到荷載列陣，這個(gè)方法降低了位移方程的階數(shù)。但是在面面接觸算法中，若接觸面上含有邊界約束條件，形成接觸力柔度矩陣時(shí)，容易出現(xiàn)矩陣奇異。（2）對(duì)角元素充大數(shù)法，就是把邊界位移為零的那一行對(duì)角元素充一個(gè)大數(shù)，使得非對(duì)角元素相對(duì)地較小，獲得最后求解得到的位移值趨近于零的效果。盡管這種方法較為簡便，但取值不當(dāng)往往會(huì)影響計(jì)算結(jié)果的精度，有時(shí)還會(huì)使方程變成病態(tài)而得不到解；（3）“充0置1”，即將邊界位移約束的這一行的對(duì)角線元素置1，與這個(gè)對(duì)角線元素相應(yīng)行和列的其它元素都充0，位移方程右端的相應(yīng)行為約束位移值，這樣就保證了邊界行的位移等于約束值。通過本文反復(fù)的數(shù)值試驗(yàn)，認(rèn)為第三種方法的位移約束邊界處理方法最適宜本文面面接觸模型隱式接觸力的求解模式。

面面接觸模型的優(yōu)點(diǎn)可用于任意形狀多體接觸面，接觸面可以具有不同的網(wǎng)格。但這樣在接觸分析過程中需要進(jìn)行接觸狀態(tài)判定，以及確定接觸點(diǎn)位置。面面接觸算法通過全局搜索粗略判斷所有可能潛在接觸測試對(duì)，再由局部搜索找出接觸投影點(diǎn)所在單元面的相對(duì)位置關(guān)系，從而確定接觸節(jié)點(diǎn)與單元面間的接觸狀態(tài)。本文全局搜索算法融合了主從面算法與位碼算法，局部搜索算法繼承了點(diǎn)面算法與內(nèi)外算法的優(yōu)點(diǎn)。本文研究在現(xiàn)有算法作了以下工作：（1）在全局搜索中構(gòu)建了描述“接觸面對(duì)”的主面單元和從面單元的共享實(shí)常數(shù)，可以將主面單元和從面單元相關(guān)聯(lián)實(shí)現(xiàn)快速搜索；（2）在局部搜索中采用了相關(guān)單元節(jié)點(diǎn)法向“平均向量”，從而解決搜索盲區(qū)問題，同時(shí)避開求解了點(diǎn)到面投影的非線性方程；（3）在接觸力求解方面，采用拉格朗日乘子法的動(dòng)力方程的增量法，以便解決材料非線性問題；（4）在接觸搜索預(yù)處理過程中已考慮了接觸區(qū)域的分解，便于后續(xù)的并行程序編制和并行計(jì)算；（5）按照所提出的接觸算法策略編制了FORTRAN源代碼程序。

由本文算法程序得到了計(jì)算結(jié)果幾乎與赫茲接觸問題的理論解完全吻合，同時(shí)借助于ANSYS商業(yè)軟件提供的小球算法和拉格朗日乘子法的功能得到的計(jì)算數(shù)據(jù)也與本文算法計(jì)算結(jié)果幾乎完全一致。更重要是本文算法在接觸搜索預(yù)處理過程中已考慮了接觸區(qū)域的分解，便于實(shí)現(xiàn)后續(xù)對(duì)高混凝土壩和巖體高邊坡靜動(dòng)力穩(wěn)定分析的并行計(jì)算。

致謝：本文研究得到了中國鐵道科學(xué)研究院劉金朝教授的指導(dǎo)和幫助，北京希格瑪仿真技術(shù)有限公司完成了程序的調(diào)試和測試，作者在此謹(jǐn)表示衷心的感謝！最后，作者要特別感謝匿名評(píng)閱人對(duì)本文提出的非常中肯且具建設(shè)性的修改意見和建議。