FLAC3DV3.0邊坡穩(wěn)定剪切強(qiáng)度折減法的多線程并行算法

瞿秀祥，林　杭，王　敏，張道勇

(1. 中南大學(xué) 資源與安全工程學(xué)院，湖南 長沙 410083；2. 中機(jī)國際工程設(shè)計研究院有限公司，湖南 長沙 410007)

運(yùn)用FLAC3D編制強(qiáng)度折減法計算程序在邊坡穩(wěn)定性分析中應(yīng)用廣泛，研究邊坡強(qiáng)度折減法計算時間具有重要經(jīng)濟(jì)效益和實(shí)際意義。隨著計算機(jī)的迅速發(fā)展，有限元強(qiáng)度折減法在邊坡工程中受到關(guān)注[1?7]。其中較多學(xué)者選擇FLAC3D軟件進(jìn)行邊坡強(qiáng)度折減法數(shù)值模擬試驗(yàn)研究。薛雷等[8]以 FLAC3D應(yīng)用軟件為計算平臺，基于計算收斂性準(zhǔn)則利用內(nèi)嵌FISH語言二次開發(fā)出能夠自動搜索安全系數(shù)的整體強(qiáng)度折減代碼和局部折減強(qiáng)度代碼。周元輔等[9]依托珍珠壩邊坡采用 FLAC3D對位移突變判據(jù)及塑性區(qū)判據(jù)進(jìn)行分析，并與ANSYS計算結(jié)果進(jìn)行比較得出：在滯滑帶明確的情況下采用塑性區(qū)貫通率增量突變判據(jù)，在滯滑帶未知的情況下采用位移增量突變判據(jù)。李躍等[10]以攀鋼蘭尖鐵礦山第7排土場邊坡為實(shí)例，采用FLAC3D軟件進(jìn)行三維強(qiáng)度折減法，數(shù)值模擬試驗(yàn)結(jié)果與現(xiàn)場調(diào)查及平面計算結(jié)果較為一致。張昊等[11]通過理論分析和FLAC3D強(qiáng)度折減法計算探討了強(qiáng)度折減法邊坡穩(wěn)定性過程中安全系數(shù)和滑動面之間的關(guān)系。敬靜等[12]利用FLAC3D軟件建立節(jié)理巖質(zhì)邊坡樁基加固的數(shù)值模型，通過強(qiáng)度折減法計算邊坡安全系數(shù)。采用FLAC3D軟件進(jìn)行邊坡強(qiáng)度折減法穩(wěn)定性分析時，判斷邊坡是否穩(wěn)定是基于模型計算最終不平衡比率滿足10?5的要求[13]。強(qiáng)度折減法對參數(shù)的折減一般采用普通的二分法進(jìn)行迭代計算，需要較多的迭代次數(shù)及計算時間。因此，研究如何減少FLAC3D邊坡穩(wěn)定強(qiáng)度折減法過程中迭代次數(shù)及計算時間具有實(shí)際意義。本文引入多線程并行運(yùn)算強(qiáng)度折減法的計算原理，并與傳統(tǒng)的二分法強(qiáng)度折減法進(jìn)行比較，分析2種方法迭代次數(shù)的差異。推導(dǎo)出二分法和多線程并行運(yùn)算迭代次數(shù)計算公式。將FLAC3D命令流內(nèi)嵌于Python腳本中，通過Python腳本同時調(diào)用多個FLAC3D應(yīng)用程序進(jìn)行多線程數(shù)值模擬計算。通過具體算例分析多線程并行運(yùn)算強(qiáng)度折減法與傳統(tǒng)二分法強(qiáng)度折減法迭代次數(shù)及計算時間上的差別。

1　強(qiáng)度折減法原理

強(qiáng)度折減法將邊坡的安全系數(shù)定義為使邊坡剛好達(dá)到臨界破壞狀態(tài)時，對其強(qiáng)度參數(shù)進(jìn)行折減的程度。邊坡采用Mohr-Coulomb準(zhǔn)則判定，影響其穩(wěn)定的強(qiáng)度參數(shù)是內(nèi)聚力c和內(nèi)摩擦角φ。強(qiáng)度折減法是先確定邊坡的原始黏結(jié)力ci和原始內(nèi)摩擦角 φi，然后將 ci，φi正切值同時除以折減系數(shù)K[2?3, 14?15]：

通過不斷試算折減系數(shù)K，反復(fù)分析計算直至邊坡達(dá)到臨界破壞狀態(tài)，此時的折減系數(shù) Kcr就是邊坡的安全系數(shù)F。強(qiáng)度折減法中一般采用二分法對折減系數(shù)進(jìn)行迭代計算，其具體的計算流程如圖1。在FLAC3D軟件中計算收斂不平衡判據(jù)為最終不平衡力比率小于或等于10?5，導(dǎo)致需要花費(fèi)較多的計算時間。采用二分法對折減系數(shù)進(jìn)行迭代是不斷將折減系數(shù)的計算下限值Kmin和上限值Kmax之間的范圍縮小為上一次迭代時的 1/2，當(dāng)折減系數(shù)上限值Kmax與下限值Kmin小于給定計算精度η時則停止計算。根據(jù)上述折減系數(shù)二分法迭代算法，則折減系數(shù)迭代次數(shù)Imin可用下式表示：

式中：Kmax為折減系數(shù)上限值；Kmin為折減系數(shù)下限值，η為給定誤差。

由式(2)可知：采用二分法進(jìn)行折減系數(shù)迭代計算時，至少要進(jìn)行Imin次迭代，折減系數(shù)才能達(dá)到所需精度η：

式中：Ceil表示為向上取整。

圖1　二分法強(qiáng)度折減法安全系數(shù)求解流程Fig.1　Dichotomy strength reduction method of safety factor solving process

綜上可知：二分法是將折減系數(shù)的范圍(Kmin，Kmax)不斷以1/2倍縮小，即每迭代一次折減系數(shù)計算范圍就縮小 1/2，直至計算差值(Kmax? Kmin)小于給定計算精度η時，則停止計算。實(shí)際應(yīng)用中當(dāng)計算區(qū)域(Kmax? Kmin)較大、計算精度η較小時，根據(jù)式(3)可知相應(yīng)迭代次數(shù)Imin將顯著增大，其相應(yīng)的計算時間也隨之增大。

2　多線程并行運(yùn)算原理及其實(shí)現(xiàn)

2.1　多線程并行運(yùn)算原理

通過上文分析可知二分法是將折減系數(shù)的計算范圍(Kmin，Kmax)不斷以1/2倍縮小，因此，減少迭代次數(shù)最為行之有效的方法是將每次迭代折減系數(shù)的計算范圍(Kmin，Kmax)以更小的倍數(shù)縮小，多線程并行運(yùn)算強(qiáng)度折減法即采用該思想。

多線程并行運(yùn)算先將計算范圍(Kmin，Kmax)等分成T+1個計算范圍，形成(Kmin，K1，K2，……，KT，Kmax)的折減系數(shù)等差序列，每次迭代同時調(diào)用T個線程分別進(jìn)行折減系數(shù)為 K1，K2，……，KT的數(shù)值模擬試驗(yàn)，根據(jù)計算結(jié)果的收斂情況，確定下一步計算范圍的上下限值，周而復(fù)始，直至(Kmax?Kmin)達(dá)到計算精度η時則計算結(jié)束。當(dāng)T=1時，即采用單線程運(yùn)算，其計算范圍(Kmin，Kmax)以1/2縮小，因此T=1是多線程并行運(yùn)算中的一個特例，即傳統(tǒng)的二分法強(qiáng)度折減。以T線程并行運(yùn)算為例說明多線程并行運(yùn)算過程，如圖2。

圖2　多線程并行計算流程圖Fig.2　Multi-threaded parallel computing flow chart

根據(jù)多線程并行運(yùn)算原理，采用FLAC3D軟件進(jìn)行計算時，數(shù)值模擬計算共需要迭代次數(shù)IT滿足下式：

式中：Kmax為折減系數(shù)的上限值；Kmin為折減系數(shù)的下限值，η為給定誤差。

則采用 T線程并行運(yùn)算時，至少要進(jìn)行 IT?min次迭代才能達(dá)到所需的計算精度：

式中：Ceil表示為向上取整。

通過比較式(3)和式(5)可知：采用二分法所需要的迭代次數(shù)是采用T線程并行運(yùn)算的一個特例。當(dāng)式(5)中的T=1時，即采用單線程運(yùn)算時所需要的迭代次數(shù)與采用二分法所需的迭代次數(shù)相等。圖3給出二分法 T分別等于 3，5和 7時計算所需迭代次數(shù)。

圖3　不同線程所需迭代次數(shù)Fig.3　Number of iterations required for different threads

多線程并行運(yùn)算時，應(yīng)先采用公式(5)進(jìn)行試算，從而確定合適的線程數(shù)。如圖3，當(dāng)=5 000時，采用5線程與7線程并行運(yùn)算所需迭代次數(shù)均為5次，實(shí)際應(yīng)用過程中，過多線程進(jìn)行應(yīng)用程序并行運(yùn)算可能造成計算機(jī)內(nèi)部進(jìn)程間計算資源相互競爭，導(dǎo)致計算時間延長，這種情況下建議采用5線程并行運(yùn)算。因此，采用多線程強(qiáng)度折減法計算時，應(yīng)根據(jù)迭代次數(shù)及計算機(jī)性能選擇并行運(yùn)算程序個數(shù)，不能盲目提高并行運(yùn)行線程個數(shù)以求快速得到計算結(jié)果。由于計算時間與計算機(jī)性能緊密相關(guān)，計算所需時間差異性較大，本文對計算時間不進(jìn)行定量理論計算與分析。

2.2　基于Python腳本的FLAC3D多線程并行運(yùn)算

Python語言設(shè)計明確、簡單，在科學(xué)計算領(lǐng)域應(yīng)用廣泛。本文將FLAC3D命令流內(nèi)嵌于Python腳本中，并由Python腳本同時調(diào)用多個FLAC3D應(yīng)用程序?qū)崿F(xiàn)并行運(yùn)算，分析計算模型收斂情況，決定下一步計算，實(shí)現(xiàn)圖2中并行運(yùn)算強(qiáng)度折減法。

為實(shí)現(xiàn)Python腳本調(diào)用FLAC3D應(yīng)用程序，先定義一個flac(reduction)函數(shù)，該函數(shù)主要實(shí)現(xiàn)折減系數(shù)為reduction的FLAC3D數(shù)值模擬計算，其中flac為函數(shù)名，括號內(nèi)參數(shù)為強(qiáng)度折減系數(shù)。通過預(yù)先安裝多線程包threading實(shí)現(xiàn)調(diào)用多個FLAC3D應(yīng)用程序。為避免多線程并行運(yùn)算過程中出現(xiàn)線程間計算資源相互競爭或死鎖現(xiàn)象，采用如下配套命令流：

Thread1=threading.Thread(target=flac,args=(K1, ))

Thread2=threading.Thread(target=flac,args=(K2, ))

……

ThreadT=threading.Thread(target=flac,args=(KT, ))

Thread1.start()

Thread2.start()

……

ThreadT.start()

Thread1.join()

Thread2.join()

……

ThreadT.join()

上述命令流中Thread1，Thread2，……ThreadT為線程名；K1，K2，……KT為線程名為Thread1，Thread2，……ThreadT所對應(yīng)的FLAC3D應(yīng)用程序計算的折減系數(shù)。由于Python腳本中對字母大小寫敏感，上述配套命令流中不得隨意更改其大小寫狀態(tài)，避免程序運(yùn)行過程中出現(xiàn)錯誤。

3　算例分析

3.1　數(shù)值模擬計算平臺設(shè)計

本文并行運(yùn)算采用的應(yīng)用軟件為 FLAC3DV3.0，Python版本為Python 2.7.5，所使用Python編譯器為Pycharm Community 4.5.3(開源，無需付費(fèi))。并預(yù)先安裝數(shù)值模擬試驗(yàn)中需要的 Python軟件包如表1。以上軟件包均遵守GNU協(xié)議，即軟件包開源無需付費(fèi)且用戶可隨意更改這些軟件包中的程序。計算機(jī)硬件對多線程并行運(yùn)算執(zhí)行效率影響較大，為此表2給出數(shù)值模擬實(shí)驗(yàn)的主要硬件設(shè)備及系統(tǒng)等參數(shù)指標(biāo)。

表1　預(yù)先安裝Python的軟件包Table1　Pre-install Python packages

表2　數(shù)值模擬實(shí)驗(yàn)平臺Table2　Numerical simulation experiment platform

3.2　數(shù)值模擬計算模型及參數(shù)

本文以均質(zhì)邊坡作為分析對象，根據(jù)平面應(yīng)變建立計算模型(模型厚度方向取 1 m)。邊坡的高度20 m，坡角45°，其余形態(tài)參數(shù)如圖4所示。模型邊界條件設(shè)置為：邊坡的左右和前后邊界約束法向位移，底部邊界為位移固定約束，其他面為自由邊界。視邊坡為理想彈塑性體，外荷載為巖土體自重，計算時先將巖土體視為彈性體，待自重應(yīng)力施加平衡后，再將模型調(diào)整為Mohr-Coulomb模型進(jìn)行彈塑性分析。模型網(wǎng)格采用均勻劃分，該邊坡模型共816個單元，1 850個節(jié)點(diǎn)，邊坡模型的主要力學(xué)參數(shù)如表3。數(shù)值模擬計算過程總運(yùn)行30 000步，其計算收斂準(zhǔn)則為最終不平衡比率小于或等于10?5。初始折減系數(shù)下限值 Kmin和折減系數(shù)上限值 Kmax分別為0和20，給定誤差η為0.001。

圖4　FLAC3D數(shù)值模擬計算模型Fig.4　FLAC3Dnumerical simulation model

表3　邊坡模型力學(xué)參數(shù)Table3　Mechanics parameters of slope model

3.3　二分法與4線程并行運(yùn)算效率比較

綜合考慮數(shù)值模擬試驗(yàn)計算平臺及硬件設(shè)施的實(shí)際情況，本文計算機(jī)的內(nèi)核是3.5 GHz，向上取整采用4線程并行運(yùn)算。采用4線程法與二分法數(shù)值模擬計算得到的邊坡安全系數(shù)F均為1.034。因此，采用多線程數(shù)值模擬試驗(yàn)結(jié)果與傳統(tǒng)的二分法數(shù)值模擬試驗(yàn)結(jié)果相同，不存在誤差。4線程并行運(yùn)算計算結(jié)果準(zhǔn)確可靠。由于本次數(shù)值模擬試驗(yàn)初始折減系數(shù)下限值 Kmin和折減系數(shù)上限值 Kmax分別為0和20，給定計算精度η為0.001，因此，=20 000。根據(jù)式(3)和式(5)，采用傳統(tǒng)二分法和4線程法所需迭代次數(shù)分別為15次和7次，這與數(shù)值模擬計算結(jié)果吻合。4線程強(qiáng)度折減法所需迭代次數(shù)是傳統(tǒng)二分法強(qiáng)度折減法迭代次數(shù)的0.46倍。因此，采用4線程并行運(yùn)算相比于傳統(tǒng)二分法能有效減少數(shù)值計算過程中的迭代次數(shù)。

圖5給出數(shù)值模擬計算過程中分別采用二分法和 4線程并行運(yùn)算迭代的 Kmax?Kmin變化趨勢。由圖5可知，4線程并行運(yùn)算在計算過程中的收斂速度明顯大于二分法的收斂速度，且所需迭代次數(shù)也較少。相比二分法將折減系數(shù)的計算范圍不斷劃分為2個區(qū)域，采用4線程并行運(yùn)算時，每進(jìn)行一次迭代都將折減系數(shù)的范圍劃分為5個相等的區(qū)域。每次迭代計算結(jié)束后，折減系數(shù)的范圍縮小到原來的 1/5。采用 4線程并行運(yùn)算，折減系數(shù)縮小范圍的速率要快于二分法。因此，采用4線程并行運(yùn)算能加快折減系數(shù)的收斂。

圖5　4線程并行運(yùn)算與二分法Kmax?Kmin變化趨勢比較Fig.5　Comparison of 4 thread parallel computing and dichotomy Kmax?Kminchanges in trend

同時，在數(shù)值模擬計算過程中，通過 Python腳本中time模塊對二分法和4線程法數(shù)值模擬運(yùn)算時間進(jìn)行統(tǒng)計，所需計算時間分別為 3 123.9和1 735.6 s(不同硬件性能計算機(jī)相應(yīng)計算時間可能不同)，4線程并行運(yùn)算所需時間是傳統(tǒng)二分法強(qiáng)度折減法所需時間的0.55倍。采用4線程并行運(yùn)算能有效減少強(qiáng)度折減法迭代計算過程中所需要的計算時間。

綜上所述，相比于傳統(tǒng)的二分法折減系數(shù)迭代，4線程并行運(yùn)算能有效減少迭代次數(shù)及數(shù)值模擬分析計算所需要時間，其中迭代次數(shù)為傳統(tǒng)二分法的0.46倍，計算時間為二分法的0.55倍。從迭代次數(shù)及總的計算時間角度考慮，4線程并行運(yùn)算計算優(yōu)勢較為明顯。

4　結(jié)論

1) 將多線程并行運(yùn)算運(yùn)用于邊坡強(qiáng)度折減法，推導(dǎo)出多線程并行運(yùn)算和傳統(tǒng)二分法強(qiáng)度折減法所需迭代次數(shù)公式，并分析多線程并行運(yùn)算在迭代次數(shù)方面的優(yōu)勢。將FLAC3D命令流內(nèi)嵌入Python腳本中，開發(fā)出實(shí)現(xiàn) FLAC3D多線程并行運(yùn)算的Python腳本程序。

2) 通過具體算例進(jìn)行 4線程并行強(qiáng)度折減法運(yùn)算及傳統(tǒng)的二分法強(qiáng)度折減法運(yùn)算。計算結(jié)果表明：在本算例計算條件下，4線程并行運(yùn)算獲得最終計算結(jié)果需要迭代7次、耗時1 735.6 s；單線程運(yùn)算需要迭代15次、耗時3 123.9 s。4線程并行運(yùn)算迭代次數(shù)是傳統(tǒng)二分法的46%，而計算時間僅為55%。4線程并行運(yùn)算能有效減少迭代次數(shù)及計算時間。

3) 將 FLAC3D命令流內(nèi)嵌于 Python腳本中，然后通過Python腳本調(diào)用FLAC3D應(yīng)用程序，實(shí)現(xiàn)Python腳本與FLAC3D軟件的交互式運(yùn)算。Python腳本與FLAC3D相結(jié)合的交互式運(yùn)算實(shí)屬首創(chuàng)，旨在為FLAC3D數(shù)值模擬計算研究提供一種新的研究途徑。

4) 進(jìn)行多線程并行運(yùn)算試驗(yàn)時，應(yīng)根據(jù)計算機(jī)的性能而定，切不能盲目的提高多線程個數(shù)以求達(dá)到減少計算時間的目的。過度提高多線程應(yīng)用程序數(shù)量可能出現(xiàn)多線程運(yùn)算耗時長于傳統(tǒng)二分法的情況。本文尚未給出多線程數(shù)確定的具體方法，根據(jù)計算機(jī)的性能確定合適的并行運(yùn)算應(yīng)用程序數(shù)目將是下一步研究的重點(diǎn)。由于現(xiàn)今計算機(jī)性能普遍較高且支持多線程并行運(yùn)算，因此，基于Python腳本多線程并行運(yùn)算在FLAC3D軟件中的強(qiáng)度折減法具有較強(qiáng)實(shí)際意義和應(yīng)用價值。