基于ABAQUS二次開發(fā)的可靠度計算程序的設(shè)計與應(yīng)用

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

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

結(jié)構(gòu)的安全性、適用性與耐久性對結(jié)構(gòu)的可靠性提出了較高的要求，保證結(jié)構(gòu)在規(guī)定的使用期內(nèi)能夠承受設(shè)計考慮的各種作用，滿足設(shè)計要求的各項使用功能，這是工程結(jié)構(gòu)可靠性的基本內(nèi)容。為了對工程結(jié)構(gòu)的可靠性進行定量分析，將結(jié)構(gòu)在規(guī)定的條件下和規(guī)定時間內(nèi)完成規(guī)定功能的概率定義為結(jié)構(gòu)的可靠度[1]。在結(jié)構(gòu)可靠度的計算方法當中，Monte Carlo方法是最基本的方法，當計算數(shù)據(jù)量足夠時，也是相對比較準確的方法。如果通過真實試驗進行Monte Carlo模擬計算可靠度，相關(guān)試驗可能會因為條件不足或代價高昂而難以實施?？墒褂糜邢拊Ｐ吞娲囼災Ｐ?，即采用計算機模擬或預測一個結(jié)構(gòu)關(guān)于結(jié)構(gòu)參數(shù)或設(shè)計變量的性能或響應(yīng)[2]，使用Monte Carlo方法計算其可靠度。

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

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

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

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

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

2 可靠度計算程序的設(shè)計

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

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

2.1 參數(shù)生成模塊

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

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

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

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

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

L =100

#定義參數(shù)

def GAUS(Meanvalue,Standarddeviation,size):

#定義分布函數(shù)

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

scale=Standarddeviation,size=size)

return(L)

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

#調(diào)用分布函數(shù)

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

#保存為文本文件

for i in range(len(LR)):

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

file1.close()

2.2 模型建立模塊

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

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

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

2.3 結(jié)果處理模塊

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

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

結(jié)果處理模塊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方法時，為批量完成模型建立、任務(wù)提交、結(jié)果輸出等過程，需要設(shè)置專門的循環(huán)執(zhí)行模塊。該模塊的主要任務(wù)是讀取參數(shù)生成模塊所生成的隨機參數(shù)表，將其中的變量作為參數(shù)輸入模型建立模塊中，調(diào)用結(jié)果處理模塊讀取模型建立模塊生成的結(jié)果數(shù)據(jù)庫文件，保存結(jié)構(gòu)的響應(yīng)。

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

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

from ExtractionData import output

#導入結(jié)果處理模塊

from ModelBuilding import create

#導入模型建立模塊

path = r″input.txt″

#讀取參數(shù)數(shù)值表

file = open(path,″r″)

mystr = file.readline()

#一次讀取一行

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

#為參數(shù)賦值

create(L)

#創(chuàng)建模型

value=output()

#數(shù)據(jù)輸出

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

#結(jié)果保存

file2.write(value+′ ′)

2.5 可靠度計算模塊

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

Znk=sn0-snk

(1)

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

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

3 可靠度計算程序的應(yīng)用

3.1 模型介紹

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

圖2 鋼桁架橋模型

3.2 參數(shù)統(tǒng)計

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

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

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

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

g(x)=0.11-S(x)

(2)

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

3.4 計算過程

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

圖3 步驟一流程

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

圖4 步驟二流程

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

圖5 步驟三流程

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

圖6 步驟四流程

3.5 計算結(jié)果

采用本文編寫的基于有限模型的可靠度計算程序?qū)︿撹旒軜蚰Ｐ瓦M行可靠度計算，不同抽樣次數(shù)下對應(yīng)的鋼桁架橋失效概率結(jié)果見表2，可基本得到鋼桁架橋的失效概率為0.008 23。

表2 失效概率計算結(jié)果

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

圖7 位移值頻率直方圖

4 小結(jié)

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

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

c)使用設(shè)計的可靠度計算程序?qū)︿撹旒軜蛴邢拊Ｐ瓦M行可靠度計算，驗證了該程序的可行性和便捷性。對于以下情況，此可靠度計算程序有較好的應(yīng)用：①極限狀態(tài)函數(shù)為隱函數(shù)；②復雜的空間結(jié)構(gòu)；③獲取結(jié)構(gòu)可靠度的精確值；④建立有限元模型。