一種基于Neumann 級數的區間有限元方法1)

伍鵬革 倪冰雨 姜潮

(湖南大學機械與運載工程學院,長沙 410006)

引言

不確定性廣泛存在于實際工程中,如結構的材料屬性、幾何尺寸及結構所受載荷等.傳統方法主要基于概率模型對上述不確定性參數進行描述,并使用概率方法對其進行分析[1-5].使用概率方法進行不確定性分析通常需要大量的實驗樣本以獲得精確的概率分布信息.然而在許多工程實際問題中,由于測試成本或測試技術等原因往往造成樣本數據缺乏,很難獲得足夠的樣本信息以構建上述參數的精確概率分布函數.為此,一系列非概率分析方法及非精確概率分析方法得以發展,包括證據理論[6-8]、P-box 模型[9-10]、凸模型方法[11-12]、模糊集理論[13-14]和區間方法[15-17]等.其中區間方法于20 世紀50 年代末由Moore[17]提出,最早用于處理計算機內浮點運算問題,后來被弓入到工程結構領域的不確定性問題中[18-20].在使用區間方法進行不確定性分析時,參數的不確定性由上下邊界表示,不僅有效適用于小樣本條件下的不確定性度量,而且易于工程人員理解和方便使用,近年來在理論研究和實際應用方面獲得了廣泛關注[21-25].

20 世紀90 年代中期,為解決在小樣本條件下工程結構問題中的不確定性,將區間分析與有限元方法[26-27] 相結合,發展出了區間有限元方法[28-30].區間有限元是一種考慮區間不確定性的有限元分析方法,其主要目的在于借助有限元分析這一有力數值分析工具,對含區間不確定參數的結構進行響應邊界求解.而獲得有限元結構響應邊界的關鍵在于區間平衡方程組的求解,該平衡方程組屬于一類區間線性方程組,其中的剛度矩陣或載荷向量為區間矩陣或區間向量.區間線性方程組的求解被認為是一類NP-hard 問題[31-33],一般情況下很難獲得其解析結果.因此,不少學者發展了區間有限元分析的近似方法以求解實際響應集合的最小區間包絡.Rao 和Berke[34]提出了頂點組合方法,在區間參數波動范圍較小的情況下能夠得到較為精確的解.Neumaier 和Pownuk[35]針對桁架結構這一典型結構力學問題提出了一種關于區間線性系統的迭代解法,能夠得到較準確的桁架結構響應包絡解.Muhanna 和Mullen[36]提出了element-by-element(EBE)列式技術,在一定程度上避免了響應區間的過保守估計問題.Rao 等[37]和Xiao 等[38]將EBE 技術與迭代解法相結合應用于結構響應靜態以及動態分析中.Muscolino 等[39]對具有不確定性參數結構的頻率響應提出了比例多項式的近似求解方法.Sofi 和Romeo[40]基于比例多項式建立了適用于含區間變量及隨機變量的結構有限元分析方法.Qiu 和Elishakoff[41]、Chen 和Yang[29]、Qiu等[42]應用區間攝動方法求解了含有有界不確定參數的結構響應區間,對于小擾動問題有較好的求解精度.Degrauwe 等[43]將仿射運算弓入至區間變量的表示中,改善了區間運算中的區間擴張問題.

現有大多數區間有限元方法計算得到的響應區間包絡往往過為保守,或只對于參數波動范圍較小的區間有限元問題能夠得到較準確的響應區間,或計算效率較低等.因而,如何高效求解區間有限元中結構響應的最小包絡解仍是有待解決的難點問題.本文歸納了一類在工程實際中常見的區間有限元問題,其中剛度矩陣可以表示為一系列獨立區間變量線性疊加的形式,稱之為可線性分解式區間有限元問題.如彈性模量或橫截面積等為區間參數的桁架結構,彈性模量或厚度尺寸等為區間參數的混凝土板連續體結構等,對其進行區間有限元分析時均可處理為可線性分解式區間有限元問題.針對此類可線性分解式區間有限元問題,采用Neumann 級數展開方法對其剛度矩陣的逆矩陣進行表示,可獲得結構響應關于區間變量的顯式表達式,從而便于高效求解結構響應的上下邊界.本文剩余部分內容安排如下:第1 節首先簡要介紹區間有限元分析的概念,并歸納出一類工程中廣泛存在的可線性分解式區間有限元問題;第2 節針對可線性分解式區間有限元問題,提出了基于Neumann 級數的區間有限元方法; 第3 節通過2 個算例驗證了本文方法的有效性;第4 節給出了全文結論.

1 區間有限元問題描述

在區間有限元分析中,通過區間模型等度量參數的有界不確定性,通常適用于樣本信息不足以構建參數概率分布但參數上下邊界可知的結構不確定性問題中.由于區間不確定性的存在,通常區間有限元平衡方程中的剛度矩陣為一區間矩陣或節點力向量為一區間向量,從而導致其平衡方程組為一區間線性方程組.本節簡要介紹區間有限元平衡方程組的構成,并歸納出一類可線性分解式區間有限元問題.

1.1 區間有限元平衡方程

在有限元方法中,首先將結構離散為由多個單元組成的有限單元模型,計算各個單元的單元剛度矩陣,進而將離散后的單元按照單元與單元之間的共同節點相互連接起來,從而組建全局剛度矩陣K

式中,Ne為結構的離散單元總數目;Ke(i) 為第i個單元的單元剛度矩陣;Ti是對應的單元剛度矩陣的變換矩陣,將單元剛度矩陣從單元局部坐標變換到全局坐標.結構的平衡方程則可表示為

式中,K∈RNd×Nd是結構的全局剛度矩陣;u∈RNd×1表示節點位移向量,為待求量;p∈RNd×1是結構的節點等效載荷向量;Nd為整個結構的自由度數目.當給定邊界約束條件后,式(2)中節點位移向量u便可求出,進而可獲得結構的應力和應變響應.

在上述通過有限元方法求解結構在外力作用下的力學響應信息時,結構參數均假設為確定性量.對于工程中普遍存在的不確定性結構,如材料屬性、幾何尺寸、載荷等具有不確定性的情況,通過弓入區間變量描述參數的不確定性,并結合有限元分析方法,即可通過構建區間有限元分析列式對含區間不確定性的結構進行響應邊界分析,結構不確定性分析的示意圖如圖1 所示.結構材料屬性或幾何尺寸等參數的不確定性將導致結構剛度矩陣K的不確定性,而外加載荷的不確定性會弓起載荷向量p的不確定性.設結構中存在不確定性參數α=[α1,α2,···,αk]T,

圖1 結構不確定性分析Fig.1 Structural uncertainty analysis

式(4) 所表示的平衡方程組屬于一區間線性方程組,其求解結果Γ通常較為復雜,難以解析獲得.實際上,通常尋求的是包絡真實解Γ盡可能小的區間邊界uI.uI為一區間向量,每個分量描述了對應位移響應量的變化范圍,而不考慮各分量之間的相關性,因此uI又稱為最小的超立方體近似解[44].目前求解上述區間有限元問題的方法主要包含三類,即基于區間運算的方法[32,45],基于級數展開的方法[41,46]和全局優化方法[47-48].

1.2 一類可線性分解式區間有限元問題

本小節歸納了在工程結構中常見的一類區間有限元問題,該類問題的剛度矩陣可分解為一系列獨立區間變量的線性疊加形式,稱之為可線性分解式區間有限元問題.

本文僅考慮結構材料屬性或幾何尺寸參數具有不確定性,而載荷確定的情況,此時式(4)退化為

若參數αI與結構剛度矩陣K中的每個元素呈線性關系,且αI可如式(7) 表示為一組獨立區間變量線性疊加的形式,則區間剛度矩陣K(αI)也具有關于區間變量的線性疊加表達式

式中,Kj(0,1,2,···,m)為對稱的確定性矩陣.對于不同的區間有限元問題,Kj的具體形式有所不同.形如式(9) 的可線性分解式區間剛度矩陣在工程實際問題中較為常見,如彈性模量或橫截面積等為區間參數的桁架結構,具有區間彈性模量或區間厚度尺寸等參數的混凝土板連續體結構等,對其進行數值分析時,相應的區間剛度矩陣均具有上述疊加形式.下面以含區間彈性模量的連續體結構為例,通過推導給出其區間剛度矩陣的線性分解形式.

考慮一由線彈性各向同性材料構成的連續體結構,其有限元模型包含Ne個單元,Nn個節點,Nd個自由度,受到節點載荷p的作用.為便于理解,這里僅考慮彈性模量具有不確定性的情況,且由區間變量描述.設區間彈性模量EI可表示為式(6) 的形式,即

式中,E0為彈性模量EI的中值.

根據有限元理論[26],連續體的單元剛度矩陣Ke(i)有如下形式

式中,Ae(i)為第i個單元的單元積分域,t為厚度,B為幾何矩陣,D為彈性系數矩陣,與彈性模量E呈線性關系

式中,D0=E0M,Dj=EjM.將式(13)代入式(11)中,此時單元剛度矩陣也為一區間矩陣

由式(16)可發現,上述問題的區間平衡方程中區間剛度矩陣可表示為關于區間變量的線性疊加形式,因此此類問題即屬于一類線性可分解式區間有限元問題.

2 基于Neumann 級數的區間有限元分析方法

區間平衡方程組的求解是區間有限元分析的關鍵.與確定性的有限元求解相比,在區間有限元分析中最大的不同也是難點即區間平衡方程組中區間矩陣的逆運算難以求解.現有區間有限元分析方法在進行結構響應邊界分析時,通常采用一系列數值手段避免對上述區間剛度矩陣進行求逆.本節針對在上一節中歸納的可線性分解式區間有限元問題提出基于Neumann 級數的區間有限元方法,利用Neumann級數展開方法對區間剛度矩陣進行求逆運算,從而便于后續求解結構響應的上下邊界.

2.1 基于Neumann 級數的區間剛度矩陣求逆

Neumann 級數是一種適用于算子求逆的級數展開方法.設有任意矩陣A,則該矩陣的逆矩陣A-1可表示為關于矩陣A0的逆矩陣及兩者差值矩陣ΔA=A-A0的Neumann 級數展開[49],即

根據式(17),對于可線性分解式區間有限元剛度矩陣的逆矩陣,其Neumann 級數展開式可表示為

由式(18) 可以看出,在區間有限元分析中弓入Neumann 級數表示可線性分解式剛度矩陣的逆,可將原非確定性剛度矩陣的逆運算轉換為確定性名義矩陣K0的逆運算及矩陣加、減法和乘法運算,避免了區間矩陣的直接逆變換.通過上述展開方法,式(18)給出了剛度矩陣逆矩陣關于區間變量的顯式表達式,從而很大程度上方便了后續對結構響應進行不確定性分析與求解.

2.2 位移邊界分析

由式(8)可知位移向量可表示為

將式(18)代入式(19)中,位移響應向量可表示為

u0=,則式(20)可表示為

計算第s(s=1,2,...,Nd)個自由度位移響應的上邊界和下邊界從而可轉換為對如下兩個優化問題的求解

式中,u0,s表示向量u0的第s個元素,Δus表示向量Δu的第s個元素.通過求解式(22)和式(23)得到結構在每一自由度上的位移響應區間,從而可獲得整個結構的位移響應邊界.

借助比格斯（John Biggs）提出的SOLO（Structure of the Observed Learning Outcome，觀察到的學習結果的結構）評價理論[17]，對職前教師有關數據分析問題的學科認知水平進行劃分．以問題1為例加以闡釋．

對于常見的不確定參數擾動較小的問題,通常保留Neumann 展式中的低階部分即可使結果滿足精度要求.進行一階截斷時可得

此時,式(25)和式(26)即為整個結構位移響應一階截斷時的上下邊界矢量.

2.3 應力邊界分析

通過Neumann 級數展開,不確定性結構的應力也可表示為關于區間變量的顯式表達式.單元的應力場表達式為

因此單元應力場也可表示為

當位移響應u取前一階截斷時,σe(i)表示為

通過式(31)和式(32)得到每個單元的應力響應區間,整個結構的應力響應邊界即可獲得.

2.4 收斂條件與誤差分析

根據Neumann 展式(17)的收斂條件可知,區間剛度矩陣逆矩陣的Neumann 級數展開式的收斂條件為

式中,ρ(·)表示矩陣的譜半徑.

根據矩陣譜半徑的性質有

考慮位移響應向量表達式(20) 取前n階時,位移響應的絕對誤差可表示為

位移響應的相對誤差為

3 算例分析

本節通過兩個算例驗證了本文所提基于Neumann 級數展開區間有限元法的可行性和有效性.第1 個算例是橫截面積為區間變量的平面四桿桁架結構在節點力作用下的位移響應分析; 第2 個算例是彈性模量為區間變量的汽車駕駛室受到集中載荷的形變分析.

3.1 平面四桿桁架結構

圖2 平面四桿桁架[51]Fig.2 Four bar truss in the plane[51]

圖3 沿x 軸方向位移響應上下邊界Fig.3 Displacement bounds in the direction of x

圖4 沿y 軸方向位移響應上下邊界Fig.4 Displacement bounds in the direction of y

表1 節點位移上下邊界值相對誤差Table 1 Relative error of the displacement bounds

3.2 汽車駕駛室

汽車駕駛室結構設計是汽車設計中的重要環節,駕駛室結構的靜動態分析有助于評估設計性能并進行結構改進使其滿足設計要求.圖5 所示是一個簡化的汽車駕駛室車身有限元模型[26,52],該模型包括28 個節點和43 個單元,且節點1′-13′與節點1-13具有對應相同的和坐標及相反的坐標.駕駛室有限元模型詳細的節點坐標信息由表2 給出.每一根梁的橫截面積為A=0.2 in2(1 in=25.4 mm),慣性矩為Iy′=Iz′=0.003 in4,極坐標矩J=0.006 in4,泊松比為μ=0.3.節點11,11′,12 和12′的所有自由度均被約束,僅僅在節點1 處施加集中載荷,大小為Fx=-990.0 lbf 和Fy=-330.0 lbf(1 lbf=4.45N).由于缺乏足夠的實驗數據,將梁的彈性模量考慮為不確定區間參數,由于該結構為對稱結構,假設對稱位置的梁的彈性模量為同一區間參數,各單元的彈性模量區間值具體如表3 所示(以psi 為單位,1 psi=6.894 757 kPa),共有24 個相互獨立的區間參數.

圖5 汽車駕駛室模型[26,52]Fig.5 A model of auto cab[26,52]

表2 汽車駕駛室模型節點坐標(單位:in)Table 2 Node coordinates of model of auto cab(unit:in)

表3 單元彈性模量的中值和半徑(單位:psi)Table 3 Median and radius of modulus of elasticity(unit:psi)

利用本文所提基于Neumann 級數展開前一階截斷計算出的各節點在載荷作用下沿、和軸方向位移響應的上下邊界值如圖6～圖8 所示.圖9～圖11分別表示3 個方向上的節點位移響應與MCS 及GA結果相對比,該算例中區間變量為24 個,MCS 抽取樣本點數為.從圖中可看出3 種方法得到的結果基本一致,沿軸方向上的節點位移響應較大,尤其是1,2,3 這3 個節點位移響應最大,可達17 in,且第3 節點在軸方向的位移響應也較大.故而在設計階段應考慮在這些節點單元處做出優化改進以提升安全性能.

圖6 沿x 軸方向位移的上下邊界Fig.6 Displacement bounds in the direction of x

圖7 沿y 軸方向位移的上下邊界Fig.7 Displacement bounds in the direction of y

圖8 沿z 軸方向位移的上下邊界Fig.8 Displacement bounds in the direction of z

圖9 沿x 軸方向位移的上下邊界值Fig.9 Displacement bounds in the direction x

圖10 沿y 軸方向位移的上下邊界值Fig.10 Displacement bounds in the direction y

圖11 沿z 軸方向位移的上下邊界值Fig.11 Displacement bounds in the direction z

4 結論

本文針對可線性分解式區間有限元問題,提出了一種基于Neumann 級數的區間有限元方法.通過對彈性模量不確定情況下的連續體結構進行舉例分析,推導得到了其可線性分解式區間有限元的具體形式.對于工程實際中普遍存在的此類可線性分解式區間有限元問題,采用Neumann 級數對區間剛度矩陣的逆矩陣進行表示,可獲得結構響應關于區間變量的顯式表達式,從而便于后續對結構響應上下邊界的求解.根據Neumann 級數的收斂條件給出了區間有限元求解方法收斂的充分不必要條件,并給出了所求結構響應的相對誤差.通過兩個算例對所本文方法的可行性與有效性進行了驗證.結果表明,本文基于Neumann 級數的區間有限元方法在進行結構響應邊界分析時具有較好的求解效率和求解精度.

附 錄

區間有限元的Monte Carlo 模擬(MCS)方法借助于Monte Carlo模擬方法的思想,通過對區間參數進行大量抽樣并進行確定性有限元分析,并根據響應集合的極值確定響應邊界.其具體步驟如下:

(1)獲取結構區間參數的上下邊界;

(2) 在m維超立方體中按均勻分布抽取區間變量樣本,共抽取Ns組樣本點;

(3)根據所抽取的樣本點確定結構的全局剛度矩陣,進行確定性有限元分析求解平衡方程(2)獲得結構位移響應;重復該過程,遍歷所有樣本變量,獲得結構位移響應的集合;

(4) 根據結構位移響應的集合判斷出各單元節點位移響應的最大值與最小值,從而構成結構位移響應的上邊界和下邊界.

注:由上述步驟2 中抽取的樣本將均勻分布在m維超立方體范圍內.實際上,也可以在抽樣中使用任何其他分布形式獲得樣本,唯一的要求是要確?？梢垣@得m維超立方體范圍內的任何點.