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

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

(成都理工大學(xué)地質(zhì)災(zāi)害防治與地質(zhì)環(huán)境保護(hù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，四川 成都 610059)

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

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

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

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

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

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

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

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

圖1 有限元法可靠度計(jì)算流程

2.1 參數(shù)生成模塊

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

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

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

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

d)保存抽樣結(jié)果為文本文件，方便后續(xù)調(diào)用。新建一個(gè)名稱為input.txt的文本文件，將第三步生成的隨機(jī)數(shù)按列寫入此文件，每列之間用空格間隔。生成的文本文件的列數(shù)等于設(shè)計(jì)參數(shù)的個(gè)數(shù)，行數(shù)等于預(yù)先設(shè)定的隨機(jī)數(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)的可靠度計(jì)算中，對(duì)結(jié)果影響最大的是模型的準(zhǔn)確性，準(zhǔn)確的有限元模型將會(huì)得到優(yōu)良的可靠度計(jì)算結(jié)果。同時(shí)，模型的建立也是可靠度計(jì)算過程中最復(fù)雜的步驟。在有限元模型的建立和分析中，使用參數(shù)化有限元分析方法，利用參數(shù)來描述結(jié)構(gòu)特征，通過參數(shù)來表征分析過程，從而實(shí)現(xiàn)可變結(jié)構(gòu)參數(shù)的有限元分析[10]。

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

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

2.3 結(jié)果處理模塊

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

ABAQUS的結(jié)果數(shù)據(jù)保存在工作目錄下后綴名為odb的數(shù)據(jù)庫文件當(dāng)中，可通過Python腳本讀取結(jié)果數(shù)據(jù)。例如，獲取特定區(qū)域中的最大應(yīng)力(應(yīng)變)值，首先獲取特定區(qū)域的節(jié)點(diǎn)(單元)號(hào)并創(chuàng)建節(jié)點(diǎn)(單元)集，按照節(jié)點(diǎn)(單元)集中的每一個(gè)節(jié)點(diǎn)(單元)的編號(hào)到數(shù)據(jù)庫中讀取其對(duì)應(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é)點(diǎn)集的編號(hào)

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

循環(huán)執(zhí)行模塊的具體流程如下：讀取參數(shù)生成模塊輸出的input.txt文件，通過按行讀取以實(shí)現(xiàn)循環(huán)，通過while函數(shù)判斷是否結(jié)束循環(huán);將讀取的單行數(shù)據(jù)處理后，分別對(duì)模型建立模塊生成的命令流文件當(dāng)中的參數(shù)符號(hào)進(jìn)行賦值，將修改過后的命令流文件提交到ABAQUS中，計(jì)算得到后綴名為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

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

from ModelBuilding import create

#導(dǎo)入模型建立模塊

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 可靠度計(jì)算模塊

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

Znk=sn0-snk

(1)

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

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

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

3.1 模型介紹

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

圖2 鋼桁架橋模型

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

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

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

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

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

g(x)=0.11-S(x)

(2)

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

3.4 計(jì)算過程

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

圖3 步驟一流程

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

圖4 步驟二流程

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

圖5 步驟三流程

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

圖6 步驟四流程

3.5 計(jì)算結(jié)果

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

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

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

圖7 位移值頻率直方圖

4 小結(jié)

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

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

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