卿光輝，張春潮

（中國(guó)民航大學(xué)航空工程學(xué)院，天津 300300）

一種簡(jiǎn)化的區(qū)間B樣條小波雜交應(yīng)力元

卿光輝，張春潮

（中國(guó)民航大學(xué)航空工程學(xué)院，天津 300300）

結(jié)合雜交應(yīng)力元理論，分別以2階、4階0尺度小波函數(shù)為基礎(chǔ)，構(gòu)造了一種對(duì)應(yīng)于雜交應(yīng)力元的簡(jiǎn)單形函數(shù)。這種形函數(shù)不僅克服了區(qū)間網(wǎng)格數(shù)成2的指數(shù)級(jí)增加而導(dǎo)致未知量增加的困難，同時(shí)簡(jiǎn)化了計(jì)算過(guò)程。為了驗(yàn)證理論的正確性，分析了二維彈性板在承受均布載荷下中線位移和應(yīng)力的情況。在相同網(wǎng)格數(shù)的情況下，4階小波尺度函數(shù)構(gòu)造的雜交應(yīng)力元計(jì)算結(jié)果比2階的計(jì)算結(jié)果更加準(zhǔn)確。

有限元；雜交應(yīng)力元；小波函數(shù)；形函數(shù)

小波有限元繼承了傳統(tǒng)有限元法離散逼近的優(yōu)點(diǎn)，能有效處理復(fù)雜的邊界條件[1]。同時(shí)，由于小波函數(shù)擁有多分辨、多尺度特性，所以得到一種提高精度的細(xì)化算法，即小波有限元在不改變網(wǎng)格劃分的條件下提高其分辨率，使其可以在問(wèn)題大梯度處、小梯度處分別采用高階單元和低階單元，從而提高分析精度和分析效率[2]。

雜交單元應(yīng)力場(chǎng)和應(yīng)變場(chǎng)都是獨(dú)立假設(shè)的，所以可作為獨(dú)立的變量來(lái)求解，而不像位移法那樣求解單元的應(yīng)力依賴于求解位移的導(dǎo)數(shù)。因此，基于雜交單元方法，不但應(yīng)力精度更高，而且還具有數(shù)值穩(wěn)定性好的優(yōu)點(diǎn)[3]。傳統(tǒng)位移法等參元所遇到解的自鎖問(wèn)題對(duì)雜交單元來(lái)說(shuō)是不存在的，或容易解決。雜交單元靈活地構(gòu)造和選配單元應(yīng)力場(chǎng)，對(duì)于確保裂紋前沿等奇異性問(wèn)題和復(fù)合材料界面問(wèn)題的可靠性非常重要[4]。

由于小波函數(shù)在數(shù)學(xué)計(jì)算上的復(fù)雜性，目前只在科研和工程上廣泛應(yīng)用，但可以預(yù)見(jiàn)，隨著研究的不斷深入，小波雜交元將會(huì)成為有限元一個(gè)強(qiáng)有力的分支，未來(lái)有可能具有與普通的位移元方法一樣廣泛的工程應(yīng)用前景。

文獻(xiàn)[5]中采用小波單元時(shí)，在區(qū)域內(nèi)以n×n（n= 2j+m-2，j＞jO，jO為保證至少有一個(gè)內(nèi)部小波的最小尺度）的形式劃分網(wǎng)格，這樣的網(wǎng)格劃分具有一定的局限性。特別是所用的樣條小波函數(shù)的階數(shù)較大時(shí)，網(wǎng)格將會(huì)成2的指數(shù)級(jí)增加，毫無(wú)疑問(wèn)，對(duì)于一般的問(wèn)題，這樣處理會(huì)給數(shù)值計(jì)算帶來(lái)困難。

本文在不影響計(jì)算結(jié)果精度的情況下，對(duì)小波插值函數(shù)進(jìn)行了簡(jiǎn)化。即通過(guò)直接使用某一區(qū)間段上的樣條小波尺度函數(shù)，結(jié)合雜交應(yīng)力元理論的推導(dǎo)方法，按照形函數(shù)的性質(zhì)構(gòu)造了一種簡(jiǎn)單的形函數(shù)。這種處理方法克服了網(wǎng)格數(shù)成2的指數(shù)級(jí)增加而導(dǎo)致未知量增加的困難，并使得計(jì)算更加方便、快捷。

1 雜交應(yīng)力元

首先，在單元內(nèi)部假設(shè)獨(dú)立的應(yīng)力場(chǎng)

其中：P={σ1…σM}為應(yīng)力矩陣；σi為假設(shè)的應(yīng)力模式；βi為對(duì)應(yīng)的應(yīng)力參數(shù)。位移場(chǎng)

根據(jù)（H-R）廣義變分原理[6]，可獲得應(yīng)力參數(shù)和節(jié)點(diǎn)位移之間的關(guān)系為

因此，矩陣形式的控制方程為

其中

為單元?jiǎng)偠染仃?；fe為等效節(jié)點(diǎn)載荷；G為杠桿矩陣；H為單元柔度矩陣，其相應(yīng)的表達(dá)式為

其中

為柔度矩陣；E為彈性模量；μ為泊松比。

根據(jù)一般有限元的方法，可以求得未知節(jié)點(diǎn)位移ae。然后由式（2）求得單元位移場(chǎng)，并由式（1）和式（3）求得單元應(yīng)力場(chǎng)[6]為

2 小波單元的構(gòu)造

以下由B樣條小波函數(shù)的尺度函數(shù)構(gòu)造滿足形函數(shù)條件的小波插值函數(shù)，并在此基礎(chǔ)上構(gòu)造出求解平面問(wèn)題的小波雜交應(yīng)力單元。

2.1 插值函數(shù)的選擇

0尺度2階和4階樣條小波尺度函數(shù)的表達(dá)式分別如下[2]

2.2 二維四節(jié)點(diǎn)2階樣條小波單元

自然坐標(biāo)系中二維四節(jié)點(diǎn)單元如圖1所示。

圖1 自然坐標(biāo)系中二維四節(jié)點(diǎn)單元Fig.1 2-D four-node element of natural coordinates

在文獻(xiàn)[6-7]中，是以矩陣的形式假設(shè)位移模式。但在這里，采用以下形式假設(shè)單元中任一點(diǎn)的橫向位移模式[8]為

其中，a1、a2、a3、a4為待定系數(shù)。當(dāng)ξ=-1，η=-1時(shí)，u= u1=a4；當(dāng)ξ=1，η=-1時(shí)，u=u2=a2；當(dāng)ξ=1，η=1時(shí)，u=u3=a1；當(dāng)ξ=-1，η=1時(shí)，u=u4=a3。

于是可得

因此，式（12）可寫成

同理可得單元內(nèi)任一點(diǎn)縱向位移v，其形式與式（12）相同。因此，小波形函數(shù)為

將式（13）代入式（2），可求得二維四節(jié)點(diǎn)小波單元的位移場(chǎng)

其中：ue為單元位移場(chǎng)

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

為小波插值函數(shù)。

2.3 二維四節(jié)點(diǎn)4階樣條小波單元

采用與2階樣條小波相似的方法，假設(shè)單元中任一點(diǎn)的橫向位移為

同理，可以求出式（17）中的系數(shù)

所以可得到其形函數(shù)為

3 應(yīng)用實(shí)例

算例1 如圖2和圖3所示，二維空間中有一兩端固支正方形彈性平面體，厚度t=0.1 m，長(zhǎng)度L=1 m，彈性模量E=3×104Pa（本文的彈性模量與參考文獻(xiàn)相比做了變化，但只會(huì)影響位移的數(shù)量級(jí)，計(jì)算誤差時(shí)沒(méi)有影響，即E的階數(shù)變小使變形圖的變形更明顯），泊松比μ=0.16，在其面內(nèi)承受y方向的均布載荷q=1×105N/m[2]。

圖2 承受面內(nèi)均布載荷作用的兩邊固支板（30×30）Fig.2 Both sides of fixed plate under uniformly distributed loads（30×30）

在以下的數(shù)值對(duì)比中，取模型中有代表性的點(diǎn)進(jìn)行比較分析，即取中線（沿x方向）上節(jié)點(diǎn)的位移v與中線（沿y方向）上節(jié)點(diǎn)的應(yīng)力σx與對(duì)應(yīng)的理論解[9]比較，同時(shí)用ANSYS軟件分別建立30×30和20×20網(wǎng)格有限元模型分析計(jì)算。最后分別計(jì)算誤差，如表1和表2所示。其變形圖如4所示。

圖4 網(wǎng)格劃分Fig.4 Grid division

表1 用BSWI20插值的位移v和應(yīng)力σx計(jì)算結(jié)果及其與理論解的誤差Tab.1 Result and error of displacement v and stress σxfor BSWI20 interpolation

表2 用BSWI40插值的位移v和應(yīng)力σx計(jì)算結(jié)果及其與理論解的誤差Tab.2 Result and error of displacement v and stress σxfor BSWI40 interpolation

通過(guò)分析表1和表2可以注意到，位移和應(yīng)力的誤差均在允許范圍之內(nèi)。所以本文構(gòu)造的簡(jiǎn)化小波雜交應(yīng)力元是正確的，并且在劃分較少網(wǎng)格情況下本文的算法更加精確，特別是用4階樣條小波函數(shù)所得結(jié)果。觀察ANSYS所得應(yīng)力結(jié)果，發(fā)現(xiàn)在正負(fù)應(yīng)力交替的地方，誤差突然變大，說(shuō)明在此節(jié)點(diǎn)位置有一定的奇異性，而傳統(tǒng)有限元在網(wǎng)格劃分較少的情況下不能準(zhǔn)確地將這一變化刻畫出來(lái)，本文算法可以較為真實(shí)地反映此處應(yīng)力變化情況。同時(shí)還可看出，位移求解精度很高，而應(yīng)力求解結(jié)果與理論解相比誤差較大，其原因在于由位移結(jié)果所產(chǎn)生的應(yīng)力是通過(guò)式（1）和式（3）微分運(yùn)算得到的。

算例2 將算例1所建立的模型用不同的網(wǎng)格數(shù)進(jìn)行分析，比較在劃分不同網(wǎng)格的情況下，最大位移計(jì)算結(jié)果與理論結(jié)果的誤差百分比。本例沒(méi)有對(duì)最大應(yīng)力進(jìn)行比較，因?yàn)閼?yīng)力是采用位移結(jié)果且通過(guò)式（1）和式（3）微分運(yùn)算得到的，所以只用位移的精度來(lái)說(shuō)明此方法的正確性。最大位移v的比較結(jié)果如表3和表4所示（其中理論解：最大位移v=-2.281 4）。

表3 2階小波函數(shù)在不同網(wǎng)格數(shù)下的最大位移vTab.3 Maximum displacement v of 2-order wavelet function in different grid numbers

表4 4階小波函數(shù)在不同網(wǎng)格數(shù)下的最大位移vTab.4 Maximum displacement v of 4-order wavelet function in different grid numbers

從表3和表4可以看出，兩種小波尺度函數(shù)構(gòu)造的插值函數(shù)都能使計(jì)算結(jié)果非常精確，說(shuō)明此種簡(jiǎn)化的小波雜交應(yīng)力元的構(gòu)造是正確的。

由表3和表4比較，隨著劃分網(wǎng)格數(shù)的增加，由兩種小波雜交應(yīng)力元所計(jì)算的最大位移結(jié)果與理論解的誤差呈逐漸減小趨勢(shì)。而且在相同網(wǎng)格數(shù)的情況下，4階小波尺度函數(shù)構(gòu)造的雜交應(yīng)力元計(jì)算結(jié)果比2階的計(jì)算結(jié)果更加準(zhǔn)確，說(shuō)明高階樣條小波的尺度函數(shù)作為插值函數(shù)與雜交應(yīng)力元結(jié)合具有優(yōu)越性。

4 結(jié)語(yǔ)

本文為區(qū)間B樣條小波有限元提出了一種簡(jiǎn)化的位移模式，即將矩陣形式換成了函數(shù)運(yùn)算，并與雜交應(yīng)力元結(jié)合，構(gòu)建了簡(jiǎn)化的小波雜交應(yīng)力元。最后，通過(guò)對(duì)二維彈性方形板應(yīng)力和位移的計(jì)算結(jié)果進(jìn)行分析，驗(yàn)證了該小波雜交應(yīng)力元的正確性和優(yōu)越性?？梢缘贸?，階次高的區(qū)間小波用較少的網(wǎng)格就能得到理想結(jié)果，節(jié)省單位剛度矩陣形成、方程求解、網(wǎng)格劃分的計(jì)算時(shí)間。網(wǎng)格數(shù)目越大，計(jì)算結(jié)果越精確。

[1]SWELDENSW.Wavelets：what next[J].IEEE，Proceeding of the IEEE，1996，4：680-685.

[2]何正嘉，陳雪峰.小波有限元理論及其工程應(yīng)用[M].北京：科學(xué)出版社.2006.

[3]卞學(xué)鐄.雜交應(yīng)力有限元法的研究進(jìn)展[J].力學(xué)進(jìn)展，2001（3）：344-349.

[4]吳長(zhǎng)春，黃茂光.雜交有限元進(jìn)展[J].自然科學(xué)進(jìn)展，1999（12）：3-10.

[5]劉艷紅，商中新.區(qū)間B樣條小波雜交應(yīng)力元分析及其應(yīng)用[J].中國(guó)民航大學(xué)學(xué)報(bào)，2013，31（2）：75-79.

[6]PIAN D P CHEN.On the suppression of zero energy deformation modes[J].Int J Numer Methods Engng，1973，19（12）：1741-1752.

[7]譚德坤，孫 輝，付雪峰.基于小波有限元的平面問(wèn)題求解方法[J].河北工程大學(xué)學(xué)報(bào)（自然科學(xué)版），2008（3）:106-109，112.

[8]韓建剛.小波有限元理論及其在結(jié)構(gòu)工程中的應(yīng)用[D].西安：西安建筑科技大學(xué)，2003.

[9]TIMOSHENKO S P，GOODIER J N.Theory of Elasticity[M].3rd ed.New York:McGraw-Hill Press，1970.

（責(zé)任編輯：楊媛媛）

Simplified interval B-spline wavelet on interval hybrid stress element

QING Guang-hui，ZHANG Chun-chao
（College of Aeronautical Engineering，CAUC，Tianjin 300300，China）

Combining with the hybrid stress finite element theory，based on 2-order/4-order 0 scale wavelet functions，a type of shape function is presented and the process of computer is reduced.The shape function can overcome the difficulty of unknown number increase exponentially with 2 to the power of n.In order to verify the correctness of current formulation，the median displacement and stress of 2-D elastic rectangular plate under uniformly distributed load is analyzed with different meshes.With the same mesh，the results of 4-order wavelet function are more accurate than that of 2-order wavelet function.

finite element；hybrid stress element；wavelet function；shape function

O343.1

：A

：1674－5590（2015）06－0050－05

2014-12-15；

：2015-01-04

：國(guó)家自然科學(xué)基金項(xiàng)目（60979001）

卿光輝（1968—），男，湖南新化人，教授，博士，研究方向?yàn)榻Y(jié)構(gòu)力學(xué).