史玉潔

(山西省中部引黃工程建設(shè)管理局，山西太原 030000)

0 引言

譜方法是一種結(jié)合譜方法和有限元思想來(lái)求解偏微分方程的數(shù)值方法，它的最大優(yōu)點(diǎn)在于具有“無(wú)窮階”的收斂性，適當(dāng)?shù)淖V方法所求得的近似解以N－1次的收斂速度逼近精確解［1］，N代表選取的基函數(shù)的個(gè)數(shù)，這一點(diǎn)是有限元法無(wú)法比擬的;并且譜方法采用的插值函數(shù)(如Fourier級(jí)數(shù)、Legendre多項(xiàng)式、Chebyshev多項(xiàng)式)具有無(wú)限可微的性質(zhì)，相對(duì)于有限元法采用的有限次可導(dǎo)多項(xiàng)式作為插值函數(shù)的方法，具有顯著的精度優(yōu)勢(shì)［2］。

借助有限元中等參元的思想，Patera［3］結(jié)合譜方法的精度于1984年提出了譜元方法(SEM)。近年來(lái)，譜方法已被成功地應(yīng)用于求解各種實(shí)際問題。Usik Lee等［4］將Fourier譜方法運(yùn)用到一些梁結(jié)構(gòu)動(dòng)力分析中，通過譜元法與解析解以及有限元方法結(jié)果的對(duì)比證明了譜元法的高精度，M.Krawczuk等［5］討論了有縫的Timoshenko梁動(dòng)力分析的 Fourier譜元法，Dimitri Komatitsch等［6］在考慮大地模型的異質(zhì)性的情況下用譜方法模擬地震波的傳播，國(guó)內(nèi)林偉軍，秦國(guó)良等［7，8］長(zhǎng)期從事研究彈性波傳播模擬的Chebyshev譜元法。

本文利用離散傅里葉逼近的方法對(duì)未知函數(shù)逼近，借助有限元法中等參單元的思想，在坐標(biāo)變換的時(shí)候同樣采用傅里葉譜逼近的方法，可推導(dǎo)出譜元法中單元?jiǎng)偠染仃噆ei的形成過程，在此基礎(chǔ)上編制而成的譜元法計(jì)算程序，能夠較好的求解線彈性靜力問題。通過數(shù)值算例，對(duì)比譜元法與有限元方法的計(jì)算精度，能得到令人滿意的結(jié)果。

1 函數(shù)的Fourier譜近似

其中，

f(x)的傅里葉展開式可以由以下推導(dǎo):

其中，

在傅里葉譜方法中，節(jié)點(diǎn)的選取有以下選擇:

在實(shí)際應(yīng)用中，式(2)是最常被采用的。

對(duì)三維函數(shù)u(ξ，η，ζ)，在標(biāo)準(zhǔn)的正方體單元內(nèi)，分別定義ξ，η，ζ三個(gè)方向上的節(jié)點(diǎn)系{ξj}，{ηk}，{ζl}(見圖1)。

圖1 正方體單元

彈性靜力計(jì)算中的基本未知量位移u(ξ，η，ζ)的傅里葉譜展開式為:

其中，

在單元內(nèi)若節(jié)點(diǎn)(ξj，ηk，ζl)對(duì)應(yīng)的單元節(jié)點(diǎn)編號(hào)為 i，則我們記 Njkl=hj(ξ)hk(η)hl(ζ)，則式(5)可以表示為:

其中，N+1等于單元內(nèi)總節(jié)點(diǎn)數(shù)目;Ni為Fourier譜近似情況下，反映單元的位移形態(tài)，即單元位移的形函數(shù)。可以將它代入后，形成在傅里葉譜近似下的形函數(shù)矩陣N。需要注意的是，本節(jié)中的譜近似直接是在標(biāo)準(zhǔn)的正方形單元中進(jìn)行的，其ξ，η，ζ三個(gè)方向跨度區(qū)間均為(0，2π)。顯然在實(shí)際計(jì)算中是不可能全部得到這樣理想的單元的，因此同有限元方法一樣，譜元法也要對(duì)單元進(jìn)行等參化處理。

2 分布函數(shù)

有限元方法通過坐標(biāo)變化將原來(lái)的不規(guī)則邊界單元轉(zhuǎn)化為標(biāo)準(zhǔn)單元，通過形函數(shù)將原有坐標(biāo)系(x，y，z)和新的坐標(biāo)系(ξ，η，ζ)建立一一對(duì)應(yīng)的關(guān)系，將原來(lái)整體坐標(biāo)系變換成新的各個(gè)單元下的局部坐標(biāo)系，在標(biāo)準(zhǔn)單元內(nèi)的計(jì)算提高計(jì)算的速度。如果坐標(biāo)變換和函數(shù)插值采用相同的節(jié)點(diǎn)和相同的插值函數(shù)，則稱此為等參單元。

譜元法利用了有限元單元?jiǎng)澐忠约暗葏⒒乃枷耄煌氖荈ourier譜元法選取Fourier多項(xiàng)式的極值點(diǎn)作為配置點(diǎn)，形函數(shù)選用的是具有無(wú)窮階可微性質(zhì)的譜展開式(h函數(shù))，如此可在有限的插值點(diǎn)上獲得指數(shù)型的收斂速度［2］。在求解微分方程的過程中，要計(jì)算函數(shù)對(duì)總體坐標(biāo)系(x，y，z)的導(dǎo)數(shù)，因此需要導(dǎo)出總體坐標(biāo)系下導(dǎo)數(shù)計(jì)算與局部坐標(biāo)系下導(dǎo)數(shù)的關(guān)系。譜等參變換式為:

未知函數(shù)對(duì)整體坐標(biāo)導(dǎo)數(shù)的表達(dá)式如下:

其中，

3 數(shù)值積分

由于傅里葉多項(xiàng)式的周期性，在本文的計(jì)算過程中積分區(qū)間均為［0，2π］，因此需要將其變換到積分區(qū)間［－1，1］上。首先，注意到被積函數(shù)的周期性，它們都是以2π為周期的函數(shù)，則有:

作變換 ξ=πξ′，則:

將上述積分方法引入單元?jiǎng)偠染仃嚨挠?jì)算公式ke=即可得到剛度矩陣各元素的值。

對(duì)每一個(gè)單元進(jìn)行循環(huán)，將單元體等參映射到標(biāo)準(zhǔn)的正方體單元，然后根據(jù)上述方法計(jì)算得到單元?jiǎng)偠染仃噆e，通過編碼法將所有單元的ke矩陣集成為對(duì)應(yīng)的整體矩陣K，同時(shí)計(jì)算等效節(jié)點(diǎn)力向量，求解方程組即可得到彈性靜力問題的基本未知量——位移。

4 算例

考慮狹長(zhǎng)矩形截面懸臂梁(見圖2)，梁長(zhǎng)為l，梁高為h，梁厚為一個(gè)單位，左端面上受合力為P的切向分布力作用。

圖2 狹長(zhǎng)矩形截面懸臂梁

上述問題屬于平面應(yīng)力問題，經(jīng)過分析可以得到該問題的應(yīng)力和位移的精確解，其位移表達(dá)式如下:

其中，I為矩形截面對(duì)z軸的慣性矩。我們?nèi)【匦谓孛娓遠(yuǎn)=2 m，梁長(zhǎng)l=10 m，合力P=105N，彈性模量 E=27 GPa，泊松比υ=0.2。將梁體按圖3劃分網(wǎng)格，然后分別用譜元法程序和ANSYS軟件對(duì)該網(wǎng)格進(jìn)行分析，比較它們算出的位移值與精確值之間的誤差，其結(jié)果如表1，表2所示。ANSYS單元采用4節(jié)點(diǎn)Solid42，位移單位為m。

圖3 梁體示意圖

表1 節(jié)點(diǎn)x方向位移(U)

表2 節(jié)點(diǎn)y方向位移(V)

通過表1，表2可以看到與精確值相比，譜元法計(jì)算誤差普遍小于有限元計(jì)算誤差，且在x方向位移計(jì)算上精度明顯較高，但在部分節(jié)點(diǎn)(如13，14點(diǎn))y方向位移譜元法計(jì)算誤差略高，其余計(jì)算值與精確值之間的相對(duì)誤差均控制在10%以內(nèi)。需要指出的是，在上述計(jì)算中采用的插值階數(shù)均為1，即Nξ=Nη=1，也就是插值均仍停留在單元邊界上，并沒有對(duì)單元內(nèi)部點(diǎn)進(jìn)行插值，難以充分發(fā)揮譜元法的優(yōu)勢(shì)。若提高插值階數(shù)，實(shí)則在單元上增加插值節(jié)點(diǎn)，從理論上來(lái)說(shuō)譜元法計(jì)算將體現(xiàn)出明顯的優(yōu)勢(shì)，但由于目前前處理軟件無(wú)法對(duì)單元內(nèi)部進(jìn)行節(jié)點(diǎn)編號(hào)、提取等相關(guān)操作，因此提高插值階數(shù)的問題有待解決。

5 結(jié)語(yǔ)

本文將傅里葉譜元法引入到動(dòng)力問題的矩陣形成過程當(dāng)中。從傅里葉多項(xiàng)式的基本性質(zhì)出發(fā)，推導(dǎo)其在空間內(nèi)對(duì)離散未知函數(shù)的表達(dá)式，進(jìn)而得到三維情況下形函數(shù)表達(dá)式，然后對(duì)單元進(jìn)行等參變換，得到整個(gè)問題從有限元方法的基本理論轉(zhuǎn)換到傅里葉譜元法當(dāng)中。

通過對(duì)實(shí)例的計(jì)算，我們發(fā)現(xiàn)要充分實(shí)現(xiàn)譜元法計(jì)算的譜精度，很重要一點(diǎn)在于提高插值的階次，即式(5)中Nξ，Nη，Nζ的大小。提高這些值的大小相當(dāng)于實(shí)現(xiàn)對(duì)單元內(nèi)部節(jié)點(diǎn)的插值(全區(qū)域插值)，這一點(diǎn)也是譜元法與有限元的重要區(qū)別之一。然而提高插值階次會(huì)遇到如下問題:對(duì)單元內(nèi)部點(diǎn)的劃分(見圖4)，并且提取其坐標(biāo)信息并對(duì)它們編號(hào)，目前我們并沒有找到能夠?qū)崿F(xiàn)這一功能的前處理軟件，這也成為阻礙譜元法深入發(fā)揮其自身優(yōu)勢(shì)的因素。

圖4 內(nèi)部節(jié)點(diǎn)劃分

就目前理論方面分析來(lái)看，其所謂的譜精度為各領(lǐng)域的實(shí)際工程都將會(huì)帶來(lái)諸多效益，相信未來(lái)譜元法將會(huì)受到越來(lái)越多研究工作者的青睞。

［1］ 趙廷剛.若干發(fā)展方程的譜方法和譜元法［D］.上海:上海大學(xué)數(shù)學(xué)系博士學(xué)位論文，2007.

［2］ Wang Xiuming，SERIANI Geza，Lin Weijun.“Some theoretical aspects of elastic wave modeling with a recently developed spectral element method，”Science in China Series G:Physics，Mechanics ＆ Astronomy，2007，2(50):185-207.

［3］ Patera A T.A spectral element method for fluid dynamics:Lanminar flow in z channel expansion［J］.J Fluid Mech，1984(9):67-68.

［4］ Usik Lee，Joohong Kim，Andrew Y.T.Leung.The Spectral Element Method in Structural Dynamics.The Shock and Vibration Digest，2000(32):415-465.

［5］ M.Krawczuk，M.Palaczb，W.Ostachowiczb.The dynamic analysis of a cracked Timoshenko beam by the spectral element method.Journal of Sound and Vibration，2003:1139-1153.

［6］ Komatitsch D，Ritsema J.Tromp J.The spectral element simulations of globe seismic wave propagation.Geophys J int，2002:390-412.

［7］ 林偉軍.彈性波傳播模擬的Chebyshev譜元法［J］.聲學(xué)學(xué)報(bào)，2007(32):525-533.

［8］ 張榮欣，秦國(guó)良.用切比雪夫譜元法求解均勻流場(chǎng)中的聲傳播問題［J］.西安交通大學(xué)學(xué)報(bào)，2009(7):120-124.