鄒 宏，伍 劍，徐志純，楊 立

(成都理工大學地質災害防治與地質環(huán)境保護國家重點實驗室，四川 成都 610059)

結構的安全性、適用性與耐久性對結構的可靠性提出了較高的要求，保證結構在規(guī)定的使用期內能夠承受設計考慮的各種作用，滿足設計要求的各項使用功能，這是工程結構可靠性的基本內容。為了對工程結構的可靠性進行定量分析，將結構在規(guī)定的條件下和規(guī)定時間內完成規(guī)定功能的概率定義為結構的可靠度[1]。在結構可靠度的計算方法當中，Monte Carlo方法是最基本的方法，當計算數據量足夠時，也是相對比較準確的方法。如果通過真實試驗進行Monte Carlo模擬計算可靠度，相關試驗可能會因為條件不足或代價高昂而難以實施。可使用有限元模型替代試驗模型，即采用計算機模擬或預測一個結構關于結構參數或設計變量的性能或響應[2]，使用Monte Carlo方法計算其可靠度。

在有限元模型的可靠度計算方面，近些年已開展了一些研究工作。徐軍等[3]將可靠度計算與響應面有限元直接耦合，提出了可靠度響應面有限元法。蔡陽等[4]建立了重力式擋土墻可靠度分析的數學模型，并敘述了其Monte Carlo模擬的計算思路。伍國軍等[5]編制了Matlab-Abaqus聯合計算可靠度程序，對圓形隧洞開挖錨固承載力進行了可靠度分析。任斌斌等[6]使用Python語言獲取離散隨機場，與ABAQUS有限元模型相結合，計算了邊坡的可靠度。蔡德詠等[7]通過ABAQUS二次開發(fā)實現了可靠性分析與有限元程序相結合，并將其應用于復合材料定向管優(yōu)化設計方面。上述研究在有限元模型的可靠度計算方面取得了一定的進展，但仍然存在不足，首先，上述程序設計往往針對單一結構類型進行，不具有普適性，其他人員使用時需重新編寫；其次，由于有限元腳本語言與編寫程序所用語言不一致，相結合后可能會降低計算效率。

本文運用ABAQUS二次開發(fā)技術，采用Python語言將Monte Carlo模擬所需的重復模擬過程程序化，根據設計參數的概率分布自動生成隨機參數表，程序讀取隨機參數表自動完成模型建立、作業(yè)提交、結果處理等工作。對于工程設計人員，該方法簡化了重復建模過程的復雜程度，縮短了計算周期，提高了可靠度計算的效率和自動化程度；對于研究人員，模塊化結構使其可便捷地替換可靠度計算對象，添加其他可靠度計算方法，具有很強的可移植性和可擴展性。

1 ABAQUS二次開發(fā)與Python腳本

ABAQUS在各個仿真領域都日益得到越來越廣泛的應用，特別是其開放的二次開發(fā)功能。ABAQUS二次開發(fā)可分為子程序開發(fā)和用戶圖形界面程序開發(fā)兩類[8]。子程序開發(fā)使用Fortran語言，主要用于材料本構關系、自定義單元等子程序的編寫。用戶圖形界面程序開發(fā)基于Python語言，主要用于對原有ABAQUS/CAE界面的繼承和擴展，開發(fā)專用的前后處理模塊以及GUI工具等。

作為ABAQUS二次開發(fā)的工具語言，Python具有強大的功能。Python程序不僅可以實現ABAQUS/CAE中的所有前后處理操作，還可以實現許多超出ABAQUS基本功能的操作，同時也能減少很多重復性工作，大大提高計算效率。Python具有面向對象、適應性強、可擴展性強等特點[9]，加上其豐富且強大的擴展程序庫，無論是數據處理還是科學計算都具有顯著的優(yōu)勢。

2 可靠度計算程序的設計

ABAQUS模型可靠度計算程序主要包括參數生成、模型建立、結果處理、循環(huán)執(zhí)行、可靠度計算5個模塊，具體的程序流程見圖1。

圖1 有限元法可靠度計算流程

2.1 參數生成模塊

根據隨機變量的概率分布，產生足夠多的樣本值即隨機數，這一過程稱為對該隨機變量的隨機抽樣[2]。參數生成的過程就是利用隨機抽樣原理，將重要的設計參數按其各自的分布類型生成許多組隨機數，參數生成的主要步驟如下。

a)篩選并列舉出重要的設計參數。設計參數主要分為兩大類：一類是結構的基本屬性，如幾何尺寸、材料屬性、接觸條件等；另一類是施加在結構上的直接作用或間接作用，如重力荷載、溫度作用、地震作用等。

b)定義隨機數生成函數。Python語言的擴展程序庫NumPy提供了常見的隨機數生成函數，包括但不限于均勻分布、高斯分布、對數正態(tài)分布等常見分布類型，如果提供的類型不能滿足要求，可根據相關公式自行編寫。

c)通過提前對實驗數據或統(tǒng)計資料的分析，確定設計參數的分布類型、均值、標準差等統(tǒng)計數據，根據參數的分布類型調用對應的隨機數生成函數。

d)保存抽樣結果為文本文件，方便后續(xù)調用。新建一個名稱為input.txt的文本文件，將第三步生成的隨機數按列寫入此文件，每列之間用空格間隔。生成的文本文件的列數等于設計參數的個數，行數等于預先設定的隨機數抽樣次數。參數生成模塊GrowthParameter.py的部分程序如下。

L =100

#定義參數

def GAUS(Meanvalue,Standarddeviation,size):

#定義分布函數

L = np.random.normal(loc=Meanvalue,

scale=Standarddeviation,size=size)

return(L)

LR=UNIF1(L,0.1*L,10)

#調用分布函數

file1 = open(′input.txt′,′w′)

#保存為文本文件

for i in range(len(LR)):

file1.write(str(LR[i])+′ ′)

file1.close()

2.2 模型建立模塊

在結構的可靠度計算中，對結果影響最大的是模型的準確性，準確的有限元模型將會得到優(yōu)良的可靠度計算結果。同時，模型的建立也是可靠度計算過程中最復雜的步驟。在有限元模型的建立和分析中，使用參數化有限元分析方法，利用參數來描述結構特征，通過參數來表征分析過程，從而實現可變結構參數的有限元分析[10]。

模型建立模塊的主要工作是按照ABAQUS分析流程編寫可變參數的有限元分析的命令流文件，提交命令流文件并生成結果數據庫文件?？芍苯泳庉嬌擅盍魑募?，但對使用者的能力要求較高，也可通過界面輸入方式完整地完成一次有限元分析流程，在此流程中獲取命令流文件。提交命令流文件的方法有2種：①生成并修改后綴名為inp的命令流文件，將其提交到ABAQUS Command中完成計算并得到結果數據庫文件；②生成并修改后綴名為py的命令流文件，將其提交到ABAQUS GUI中完成計算并得到結果數據庫文件。

上面2種方法都能獲得相同的結果數據庫文件，但計算速度存在差別。通過第2種方式提交計算時，是先將Python命令流運行后生成后綴名inp的輸入文件，再將此文件提交計算，因此第2種方式的計算速度較慢。此外，第1種方法可同時批量計算多個有限元模型，相較于第2種更具有優(yōu)勢。

2.3 結果處理模塊

可靠度計算模塊需要從有限元模擬結果中獲得結構的響應，這個響應是一個或多個具體的值，例如最大應力、最大應變等。ABAQUS輸出的結果數據來自于整個模型或者模型的大部分區(qū)域，不但包含的數據信息量非常大，而且分為2種數據保存類型，包括以分析步劃分的場輸出數據和以幀劃分的歷史輸出數據[11]。因此，為了提高結果提取效率，需編寫專門的模塊讀取結果數據庫文件，根據結構響應的數據類型及相關要求從數據庫當中輸出特定的數據。

ABAQUS的結果數據保存在工作目錄下后綴名為odb的數據庫文件當中，可通過Python腳本讀取結果數據。例如，獲取特定區(qū)域中的最大應力(應變)值，首先獲取特定區(qū)域的節(jié)點(單元)號并創(chuàng)建節(jié)點(單元)集，按照節(jié)點(單元)集中的每一個節(jié)點(單元)的編號到數據庫中讀取其對應的應變(應力)，通過比較大小獲得最大應力(應變)值，最后刪除ABAQUS工作目錄的所有文件，避免后續(xù)新的模型文件因權限問題而無法創(chuàng)建。

結果處理模塊ExtractionData.py的部分程序如下。

Odb=openOdb(r′Job-1.odb′)

#打開ODB文件

inX=Odb.rootAssembly.instances[′PART-1-1′]

endNode = inX.nodeSets[′SET-1′]

#獲取節(jié)點集的編號

Var= Odb.steps[′Step-1′].frames[-1].fieldOutputs[′U′]

#提取位移

nset_val =Var.getSubset(region=endNode).values

stress_data = map(lambda x:[

x.nodeLabel,x.data[1]],nset_val)

Res=stress_data[0][1]

del mdb.jobs[′Job-1′]

# 刪除作業(yè)

2.4 循環(huán)執(zhí)行模塊

在使用Monte Carlo方法時，為批量完成模型建立、任務提交、結果輸出等過程，需要設置專門的循環(huán)執(zhí)行模塊。該模塊的主要任務是讀取參數生成模塊所生成的隨機參數表，將其中的變量作為參數輸入模型建立模塊中，調用結果處理模塊讀取模型建立模塊生成的結果數據庫文件，保存結構的響應。

循環(huán)執(zhí)行模塊的具體流程如下：讀取參數生成模塊輸出的input.txt文件，通過按行讀取以實現循環(huán)，通過while函數判斷是否結束循環(huán);將讀取的單行數據處理后，分別對模型建立模塊生成的命令流文件當中的參數符號進行賦值，將修改過后的命令流文件提交到ABAQUS中，計算得到后綴名為odb的結果數據庫文件，調用結果處理模塊處理數據庫文件，獲取結構的響應，并將文件保存為output.txt文件。

循環(huán)執(zhí)行模塊LoopExecution.py的部分程序如下。

from ExtractionData import output

#導入結果處理模塊

from ModelBuilding import create

#導入模型建立模塊

path = r″input.txt″

#讀取參數數值表

file = open(path,″r″)

mystr = file.readline()

#一次讀取一行

L = float(mystr.split()[0])

#為參數賦值

create(L)

#創(chuàng)建模型

value=output()

#數據輸出

file2= open(′output.txt′,′w′)

#結果保存

file2.write(value+′ ′)

2.5 可靠度計算模塊

使用Monte Carlo模擬法計算可靠度的基本思路是：當重要設計參數x1,x2,…,xn(n為參數數量)的概率分布類型、均值、標準差等數據已知時，利用算法產生符合相應重要參數概率分布的隨機數矩陣，矩陣的列數等于參數的數量，矩陣的行數等于每個參數生成的隨機數的數量，每次從中抽取一行隨機數X=(xn1,xn2,…,xnm)T(m為參數的抽樣數)組成隨機樣本輸入結構的參數化有限元模型，提交計算后提取得到一組隨機抽樣值S=(sn1,sn2,…,snk)T(k為響應的數量)，設極限狀態(tài)函數：

Znk=sn0-snk

(1)

式中sn0——失效狀態(tài)值。

當snk超過sn0，即Znk<0,則認為結構失效。由大數定律中的伯努利定理可知，當抽樣次數足夠大時，隨機事件出現的頻率近似于它的概率[12]。因此，將結構失效的次數n與總模擬次數m之比n/m近似為結構的失效概率pf，查表可得可靠度指標β。

3 可靠度計算程序的應用

3.1 模型介紹

某簡支鋼桁架橋二維模型[13]見圖2，該鋼桁架橋模型由23根桿件組成，所有水平桿都具有完全相同的彈性模量和橫截面積，斜桿也是如此。該鋼桁架橋的長度為24 m，高度為2 m，支座A為固定鉸支座，支座B為滑動鉸支座。荷載為6個相互獨立的集中荷載，大小為50 kN。

圖2 鋼桁架橋模型

3.2 參數統(tǒng)計

鋼桁架橋在施工和運營過程中，構件的幾何誤差、材料的隨機性、荷載的不確定性都會對橋梁結構的可靠度產生影響，在設計中必須考慮這些因素以確保結構的安全性[14]。將水平桿彈性模量E1、斜桿彈性模量E2、水平桿截面面積A1、斜桿截面面積A2以及各集中力的大小P1～P6作為隨機參數。各隨機參數間相互獨立，除集中力的大小P1～P6的概率分布為Ⅰ型極值分布外，其他隨機參數均滿足對數正態(tài)分布，采用Latin超立方抽樣法[15]進行抽樣。采用表1所示的參數統(tǒng)計數據。

表1 隨機參數的統(tǒng)計特征

3.3 極限狀態(tài)函數

將結構的極限狀態(tài)定義為鋼桁架橋中點(圖2點C)的位移值不超過0.11 m。因此該鋼桁架橋的極限狀態(tài)函數為：

g(x)=0.11-S(x)

(2)

式中，S(x)為點C的位移值，利用結果處理模塊從有限元模型的結果文件中提取得到。

3.4 計算過程

步驟一根據鋼桁架橋主要參數的數量n、抽樣次數m、均值muX、標準差sigamaX以及分布類型，利用參數生成模塊生成隨機數矩陣，矩陣列數等于參數數量n，行數等于參數抽樣次數m，并將其保存為input.txt文件，主要流程見圖3。

圖3 步驟一流程

步驟二循環(huán)執(zhí)行模塊按行讀取上一步生成的input.txt文件，并將各參數的隨機值傳遞到模型建立模塊中。模型建立模塊由建立鋼桁架橋有限元模型的命令流文件參數化而成，完成模型建立、作業(yè)提交等任務，并生成后綴名為odb的結果文件，主要流程見圖4。

圖4 步驟二流程

步驟三循環(huán)執(zhí)行模塊按順序讀取上一步生成的結果文件，利用結果處理模塊從結果數據庫當中輸出特定數據。結果處理模塊由多個數據提取方法的命令流組成，可根據結構響應的類型選擇適合的提取方法，從結果文件中提取數據并保存為output.txt文件。鋼桁架橋模型選擇的是提取節(jié)點位移值的提取方法，主要流程見圖5。

圖5 步驟三流程

步驟四可靠度計算模塊根據按行讀取上一步生成的output.txt文件，將位移值代入功能函數g(x)。統(tǒng)計功能函數g(x)結果為負時的次數k4，用負結果數k4除以總抽樣次數m，得到鋼桁架橋的失效概率pf，查表可得對應的可靠度指標β，主要流程見圖6。

圖6 步驟四流程

3.5 計算結果

采用本文編寫的基于有限模型的可靠度計算程序對鋼桁架橋模型進行可靠度計算，不同抽樣次數下對應的鋼桁架橋失效概率結果見表2，可基本得到鋼桁架橋的失效概率為0.008 23。

表2 失效概率計算結果

當抽樣次數為105次時，隨機輸出變量(點C位移值)的頻率直方圖見圖7，位移值大多分布在0.05～0.12 m之間，位移值的均值為0.079 4 m，標準差為0.011 1 m。

圖7 位移值頻率直方圖

4 小結

a)該程序為ABAQUS建立的各類結構有限元模型的可靠度計算提供了一個便捷的工具，擴展了可靠度計算的途徑。提出了有限元模型可靠度通用計算程序開發(fā)的設計思路及主要功能，將其細化為5個模塊，便于在其他工程結構可靠度計算中修改使用。

b)通過對ABAQUS的二次開發(fā)，建立了針對有限元模型的可靠度計算程序。參數生成模塊方便了隨機數生成過程，模型建立模塊減少了繁瑣的ABAQUS界面操作，結果處理模塊簡化了隨機響應的提取。

c)使用設計的可靠度計算程序對鋼桁架橋有限元模型進行可靠度計算，驗證了該程序的可行性和便捷性。對于以下情況，此可靠度計算程序有較好的應用：①極限狀態(tài)函數為隱函數；②復雜的空間結構；③獲取結構可靠度的精確值；④建立有限元模型。