含裂紋結(jié)構(gòu)的區(qū)間B樣條小波模糊有限元分析

彭惠芬，孟廣偉，周立明，李 鋒

（1.吉林大學(xué)機(jī)械科學(xué)與工程學(xué)院，130022 長(zhǎng)春，phfdaqing@163.com;2.東北石油大學(xué)機(jī)械科學(xué)與工程學(xué)院，163318 大慶黑龍江）

含裂紋結(jié)構(gòu)的區(qū)間B樣條小波模糊有限元分析

彭惠芬1，2，孟廣偉1，周立明1，李 鋒1

（1.吉林大學(xué)機(jī)械科學(xué)與工程學(xué)院，130022 長(zhǎng)春，phfdaqing@163.com;2.東北石油大學(xué)機(jī)械科學(xué)與工程學(xué)院，163318 大慶黑龍江）

為了解決工程問題中材料和載荷的不確定性給含裂紋結(jié)構(gòu)數(shù)值分析帶來的困難，提高數(shù)值計(jì)算的精度和效率，將模糊理論與小波有限元相結(jié)合，提出了基于區(qū)間B樣條小波(B-spline wavelet on the interval，BSWI)模糊有限元分析法.該方法將啞節(jié)點(diǎn)斷裂單元鑲嵌到含裂紋結(jié)構(gòu)的小波有限元模型中，建立了含裂紋結(jié)構(gòu)的小波有限元模型，推導(dǎo)了小波有限元模糊平衡方程，并采用λ水平截集及分解定理求解模糊平衡方程;在此基礎(chǔ)上，利用虛擬裂紋閉合法計(jì)算了不同裂紋長(zhǎng)度下應(yīng)力強(qiáng)度因子隸屬函數(shù)值，并將計(jì)算結(jié)果與解析解進(jìn)行了比較.結(jié)果分析表明該方法可用較少單元更為真實(shí)、準(zhǔn)確地反映結(jié)構(gòu)響應(yīng)的變化情況，為工程實(shí)際中實(shí)現(xiàn)含不確定參數(shù)的斷裂數(shù)值分析提供了一種新途徑.

區(qū)間B樣條小波;裂紋;模糊數(shù);應(yīng)力強(qiáng)度因子;有限元

現(xiàn)有的斷裂力學(xué)數(shù)值分析方法大多都是基于傳統(tǒng)有限元模型的.傳統(tǒng)有限元理論對(duì)裂紋等奇異性問題，存在計(jì)算精度低和效率下降的弊端.近年來，人們提出了一些新的有限元法，如有限差分法［1］、邊界元法［2－3］和無網(wǎng)格法［4－5］等，但由于缺少相應(yīng)軟件的支持，使這些數(shù)值方法的發(fā)展受到了限制.

小波有限元法(Wavelet Finite Element Method，WFEM)是近幾年發(fā)展起來的一種新的數(shù)值分析方法.以小波函數(shù)或小波尺度函數(shù)作為插值函數(shù)，充分利用小波尺度函數(shù)的多分辨特性，可根據(jù)實(shí)際需要任意改變分析尺度，具有算法穩(wěn)定性好、計(jì)算精度和效率高的優(yōu)點(diǎn).因此，在裂紋等奇異性問題數(shù)值計(jì)算方面具有誘人的優(yōu)越性.目前，小波有限元法研究對(duì)象基本上是確定性問題.然而，在實(shí)際工程結(jié)構(gòu)中，由于各種因素的影響，使得結(jié)構(gòu)的物理參數(shù)、幾何參數(shù)及載荷等具有不確定性，從而導(dǎo)致結(jié)構(gòu)的響應(yīng)也具有不確定性，這種不確定性需要用模糊理論來進(jìn)行研究［6－7］.本文針對(duì)工程實(shí)際中幾何參數(shù)及載荷等不確定性裂紋結(jié)構(gòu)的有限元分析問題，提出了基于BSWI模糊有限元分析法，建立了BSWI有限元模糊斷裂分析模型，采用區(qū)間數(shù)分解方法求解模糊平衡方程，基于輸入材料和載荷模糊數(shù)的隸屬函數(shù)，利用虛擬裂紋閉合法，計(jì)算并分析了含裂紋結(jié)構(gòu)響應(yīng)量的可能性分布.數(shù)值算例證明了該方法切實(shí)可行，并具有較高的計(jì)算精度.

1 BSWI單元構(gòu)造

采用BSWI尺度函數(shù)的張量積插值構(gòu)造單元，位移函數(shù)表示為

式中:Φ1，Φ2分別為m階j尺度空間的n(n=2j+m －1)個(gè)一維BSWI尺度函數(shù)［8］;p、q分別為待求小波系數(shù)列向量;ε，η分別為單元局部坐標(biāo).局部坐標(biāo)與整體坐標(biāo)關(guān)系為

式中:x1，y1分別為整體坐標(biāo)下單元起始坐標(biāo);Lex為單元長(zhǎng)度;Ley為寬度.

單元節(jié)點(diǎn)位移列陣為

將各節(jié)點(diǎn)坐標(biāo)和位移代入式(1)和式(2)，即得:

由彈性力學(xué)知，平面應(yīng)力問題的勢(shì)能泛函為

式中:Ωe為單元求解域;t為單元厚度;F={fx，fy}T為體力向量;p={px，py}T為面力向量;R={u，v}T為位移場(chǎng)向量;D為彈性矩陣;ε為應(yīng)變矩陣，其中

式中:E，μ分別為材料的彈性模量和泊松比.

將式(3)、式(5)和式(6)聯(lián)立，由變分原理，令δΠp=0，可得單元?jiǎng)偠染仃嚍?/p>

2 BSWI斷裂單元構(gòu)造

BSWI啞節(jié)點(diǎn)斷裂單元(構(gòu)造如圖1所示)，是用來從小波有限元分析結(jié)果中提取相關(guān)信息，進(jìn)而計(jì)算裂紋尖端后面的張開位移和裂紋尖端前面的虛擬裂紋擴(kuò)展量［9－10］.

沿裂紋方向，啞節(jié)點(diǎn)斷裂單元從兩側(cè)的小波單元中共提取5個(gè)節(jié)點(diǎn):節(jié)點(diǎn)1和節(jié)點(diǎn)2對(duì)應(yīng)于裂紋尖端，并用特殊剛度的彈簧連接;節(jié)點(diǎn)3和節(jié)點(diǎn)4在裂紋尖端的后面;節(jié)點(diǎn)5在裂紋尖端的前面.在實(shí)際應(yīng)用時(shí)，節(jié)點(diǎn)1和節(jié)點(diǎn)2，節(jié)點(diǎn)3和節(jié)點(diǎn)4的坐標(biāo)分別是重合的.

圖1 啞節(jié)點(diǎn)斷裂單元及節(jié)點(diǎn)位移矢量

設(shè)Kx和Ky分別為x和y方向的彈簧剛度，由圖1可知，裂紋尖端垂直節(jié)點(diǎn)力為

裂紋尖端后面的張開垂直位移為

裂紋前面的虛擬擴(kuò)展為

3 BSWI模糊有限元分析

式中:Em為的主值，其值與不考慮彈性模量模糊性時(shí)的E值相同;EL，ER分別為的左、右展形.

根據(jù)式(7)，相應(yīng)地，小波模糊單元?jiǎng)偠染仃嚍?/p>

將式(14)代入式(12)得:

將式(15)代入式(13)得:

式中［Km］e與式(7)普通 BSWI單元?jiǎng)偠染仃囅嗤?

式中:Fm為的主值，F(xiàn)L，F(xiàn)R分別為的左、右展形.

假設(shè)彈簧連接裂尖節(jié)點(diǎn)i和i+1，組集啞節(jié)點(diǎn)斷裂單元和各BSWI模糊單元的剛度矩陣可得含裂紋結(jié)構(gòu)總模糊剛度矩陣為

L－R型模糊結(jié)構(gòu)平衡方程為

對(duì)模糊結(jié)構(gòu)小波有限元平衡方程作λ水平截集，可得區(qū)間方程:

根據(jù)區(qū)間數(shù)運(yùn)算法則和區(qū)間數(shù)分解形式的唯一性，λ水平截集下結(jié)構(gòu)區(qū)間方程的解［11］為

4 BSWI虛擬裂紋閉合法

虛擬裂紋閉合法是用于二維斷裂參數(shù)問題計(jì)算的一步分析法，具有精度高、方法簡(jiǎn)單的優(yōu)點(diǎn).本文針對(duì)I型裂紋斷裂問題進(jìn)行討論，但思路和方法適用于其他斷裂模式.如圖2所示，將裂紋從α擴(kuò)展到α+Δα所需要的功與將裂紋從α+Δα閉合到α所需要的功是相等的.

圖2 虛擬裂紋閉合法

在虛擬裂紋線上裂紋尖端節(jié)點(diǎn)力在節(jié)點(diǎn)位移上做功為

假設(shè)虛擬裂紋尖端后面的張開位移和初始裂紋尖端后面的張開位移近似相等，式(22)改寫為

數(shù)值實(shí)驗(yàn)表明虛擬裂紋閉合法對(duì)有限元網(wǎng)格尺寸并不敏感，因此，I型裂紋能量釋放率G1可近似表達(dá)為

式中B為裂紋體的厚度.

對(duì)于各向同性均勻線彈性材料而言，應(yīng)力強(qiáng)度因子KI和能量釋放率GI的關(guān)系為

當(dāng)λ取遍［0，1］中的一切值后，根據(jù)模糊分解定理可得結(jié)構(gòu)模糊應(yīng)力強(qiáng)度因子為

若只取有限個(gè)λ值進(jìn)行計(jì)算，則:

5 算例分析

如圖3所示，共線雙邊裂紋單向拉伸平板.其中:L=2 m，W=2 m，裂紋長(zhǎng)度a，板厚t=0.01 m，彈性模量 E=(210，0.04，0.04)LRGPa，泊松比μ =0.3，σ =(0.2，0.01，0.01)LRMPa，計(jì)算時(shí)，彈性模量、載荷均為L(zhǎng)－R型模糊數(shù)，其隸屬函數(shù)取為線性函數(shù)為

假設(shè)為平面應(yīng)力問題，試用基于BSWI小波有限元的模糊虛擬裂紋閉合法計(jì)算其應(yīng)力強(qiáng)度因子.

圖3 裂紋板單向拉伸

為更好地利用BSWI有限元虛擬裂紋閉合法，考慮結(jié)構(gòu)幾何和載荷的對(duì)稱性，取該結(jié)構(gòu)1/2作為計(jì)算模型，在其底部施加對(duì)稱邊界條件，采用2個(gè)二階三尺度區(qū)間B樣條小波(BSWI23)和1個(gè)啞節(jié)點(diǎn)斷裂單元求解.

表1給出了當(dāng)隸屬度為1.0時(shí)，采用2個(gè)BSWI23單元的模糊虛擬裂紋閉合法計(jì)算不同裂紋長(zhǎng)度下應(yīng)力強(qiáng)度因子與解析解的相對(duì)誤差.從表1中可看出最大相對(duì)誤差值為1.583%，說明本文提出的計(jì)算方法切實(shí)可行的，并可用較少單元獲得較高計(jì)算精度.

表1 隸屬度為1.0時(shí)模糊算法與解析解的相對(duì)誤差

表2給出了不同裂紋長(zhǎng)度α下的應(yīng)力強(qiáng)度因子KI隸屬函數(shù)值.由表2看出不同隸屬度下應(yīng)力強(qiáng)度因子隨裂紋長(zhǎng)度變化的關(guān)系.

圖4，5分別為裂紋面上距裂尖r=0.125，0.250 m處垂直位移v的隸屬函數(shù)圖.由于選取左、右基準(zhǔn)函數(shù)L(x)，R(x)相同，垂直位移v分布均為中心對(duì)稱分布;由圖4，5可見，對(duì)于不同r值，當(dāng)隸屬度較低時(shí)，垂直位移v的變化區(qū)間較大，當(dāng)隸屬度較高時(shí)，垂直位移v變化區(qū)間較小，當(dāng)隸屬度為1.0時(shí)，左右端點(diǎn)值均相等，分別為v0.125=4.569 6×10－5m和v0.250=6.446 6 ×10－5m，結(jié)果與確定性小波有限元結(jié)果一致.由此看出，利用本文建立的區(qū)間B樣條小波有限元模糊虛擬裂紋閉合法對(duì)含裂紋結(jié)構(gòu)進(jìn)行斷裂分析比確定性小波有限元分析更能全面反映結(jié)構(gòu)響應(yīng)的真實(shí)變化情況.

表2 不同裂紋長(zhǎng)度下應(yīng)力強(qiáng)度因子隸屬函數(shù)值 105Pa·m1/2

圖4 裂紋面上距裂尖r=0.125 m垂直位移v隸屬函數(shù)

圖5 裂紋面上距裂尖r=0.250 m垂直位移v隸屬函數(shù)

6 結(jié) 論

1)將本文方法的計(jì)算結(jié)果與解析解進(jìn)行了比較，結(jié)果表明:當(dāng)隸屬度為1.0時(shí)，不同裂紋長(zhǎng)度下，應(yīng)力強(qiáng)度因子隸屬函數(shù)值最大相對(duì)誤差為1.583%.驗(yàn)證了該方法在解決含有不確定性參數(shù)的裂紋結(jié)構(gòu)有限元分析方面的可靠性.

2)利用本方法計(jì)算了不同裂紋長(zhǎng)度和隸屬度下的應(yīng)力強(qiáng)度因子KI隸屬函數(shù)值.結(jié)果表明:對(duì)于不同r值，當(dāng)隸屬度較低時(shí)，垂直位移v的變化區(qū)間較大，當(dāng)隸屬度較高時(shí)，垂直位移v變化區(qū)間較小，當(dāng)隸屬度為1.0時(shí)，左右端點(diǎn)值均相等.由此看出，本文建立的計(jì)算方法對(duì)含裂紋結(jié)構(gòu)進(jìn)行斷裂分析比確定性有限元分析法更能全面反映結(jié)構(gòu)響應(yīng)的真實(shí)變化情況.

3)小波有限元與傳統(tǒng)有限元相比，在計(jì)算裂紋等奇異性問題方面可用較少單元和自由度數(shù)，獲得較高計(jì)算精度和效率，因此，在數(shù)值計(jì)算方法上為求解具有材料和載荷不確定性的裂尖奇異性問題提供了一種可行的途徑.

［1］ZHANG Yang.A finite difference method for fractional partial differential equation［J］.Applied Mathematics and Computation，2009，215(2):524－529.

［2］程長(zhǎng)征，牛忠榮，葉建喬.邊界元法計(jì)算淺表面裂紋應(yīng)力強(qiáng)度因子［J］.合肥工業(yè)大學(xué)學(xué)報(bào)(自然科學(xué)版)，2009，32(4):503－507.

［3］SUN Yuzhou，WANG Jinyan，XU Junfeng.Boundary element method for the plane elastic problem with cracks［J］.Journal of Zhong yuan Institute of Technology，2004，15(5):46－54.

［4］張雄，劉巖，馬上.無網(wǎng)格法的理論及應(yīng)用［J］.力學(xué)進(jìn)展，2009，39(1):1－36.

［5］RAO B N，RAHMAN S.Probabilistic fracture mechanics by Galerkin meshless methods-Part I:Rates of stress intensity factors［J］.Computational Mechanics，2002，28(5):351－364.

［6］蘇靜波，卲國(guó)建.基于區(qū)間分析的工程結(jié)構(gòu)不確定性研究現(xiàn)狀與展望［J］.力學(xué)進(jìn)展，2005，35(3):338－344.

［7］RASHID M K，AL-ARAIMI S A.Fuzzy algorithm and structural stiffness in error attenuation of intelligent toolpost［J］.Journal of Intelligent Manufacturing，2005(16):277－286.

［8］何正嘉，陳雪峰，李兵，等.小波有限元理論及其工程應(yīng)用［M］.北京:科學(xué)出版社，2006:40－42.

［9］XIE D.Damage progression in tailored laminated panels with a cutout and delamination growth in sandwich panels with tailored face sheets［D］.Clemson:Clemson U-niversity，2002.

［10］XIE D，Jr BIGGERS S B.Progressive crack growth analysis using interface element based on the virtual crack closure technique［J］.Finite Elements in Analysis and Design，2006，42(11):977－984.

［11］肖盛燮，王平義，呂恩琳.模糊數(shù)學(xué)在土木與水利工程中的應(yīng)用［M］.北京:人民交通出版社，2004:349－368.

Fuzzy finite element analysis of B-spline wavelet on the interval of structure with cracks

PENG Hui-fen1，2，MENG Guang-wei1，ZHOU Li-ming1，LI Feng1

(1.College of Mechanical Science and Engineering，Jilin University，130022 Changchun，China，phfdaqing@163.com;2.College of Mechanical Science and Engineering，Northeast Petroleum University，163318 Daqing，China)

To overcome the difficults caused by uncertainties of material and load in numerical analysis of structure with cracks and to improve the accuracy and efficiency of numerical calculation，a fuzzy finite element analysis method of B-spline wavelet on the interval(BSWI)is put forward by combing the fuzzy theories with wavelet finite element.Fracture elements of dummy nodes are embedded into the wavelet finite element model with cracks，fuzzy equilibrium equations on the basis of finite element of BSWI are established，and the fuzzy equilibrium equations are solved by using λ level set and decomposition theorem.Subsequently，the membership function values of stress intensity factor with different crack lengths are calculated by using virtual crack closure technique，and the calculated results are compared with analytical solution.The analysis results show that the proposed method can accurately reflect changes in structural response with fewer elements and can be a new way for engineering fracture analysis of complex structures with uncertainties.

BSWI;crack;fuzzy numbers;stress intensive factor;finite element method

O242;TB115

A

0367－6234(2011)09－0134－05

2010－12－20.

高等學(xué)校博士學(xué)科專項(xiàng)科研基金資助項(xiàng)目(20060183063);吉林省科學(xué)技術(shù)廳基金資助項(xiàng)目(20090540);吉林大學(xué)“985工程”資助項(xiàng)目.

彭惠芬(1969—)，女，博士研究生，副教授;

孟廣偉(1959—)，男，教授，博士生導(dǎo)師.

(編輯 張 紅)