線彈性靜力問題的Fourier譜元法初探

史玉潔

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

0 引言

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

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

本文利用離散傅里葉逼近的方法對未知函數(shù)逼近，借助有限元法中等參單元的思想，在坐標變換的時候同樣采用傅里葉譜逼近的方法，可推導出譜元法中單元剛度矩陣kei的形成過程，在此基礎上編制而成的譜元法計算程序，能夠較好的求解線彈性靜力問題。通過數(shù)值算例，對比譜元法與有限元方法的計算精度，能得到令人滿意的結果。

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

其中，

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

其中，

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

在實際應用中，式(2)是最常被采用的。

對三維函數(shù)u(ξ，η，ζ)，在標準的正方體單元內，分別定義ξ，η，ζ三個方向上的節(jié)點系{ξj}，{ηk}，{ζl}(見圖1)。

圖1 正方體單元

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

其中，

在單元內若節(jié)點(ξj，ηk，ζl)對應的單元節(jié)點編號為 i，則我們記 Njkl=hj(ξ)hk(η)hl(ζ)，則式(5)可以表示為:

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

2 分布函數(shù)

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

譜元法利用了有限元單元劃分以及等參化的思想，不同的是Fourier譜元法選取Fourier多項式的極值點作為配置點，形函數(shù)選用的是具有無窮階可微性質的譜展開式(h函數(shù))，如此可在有限的插值點上獲得指數(shù)型的收斂速度［2］。在求解微分方程的過程中，要計算函數(shù)對總體坐標系(x，y，z)的導數(shù)，因此需要導出總體坐標系下導數(shù)計算與局部坐標系下導數(shù)的關系。譜等參變換式為:

未知函數(shù)對整體坐標導數(shù)的表達式如下:

其中，

3 數(shù)值積分

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

作變換 ξ=πξ′，則:

將上述積分方法引入單元剛度矩陣的計算公式ke=即可得到剛度矩陣各元素的值。

對每一個單元進行循環(huán)，將單元體等參映射到標準的正方體單元，然后根據上述方法計算得到單元剛度矩陣ke，通過編碼法將所有單元的ke矩陣集成為對應的整體矩陣K，同時計算等效節(jié)點力向量，求解方程組即可得到彈性靜力問題的基本未知量——位移。

4 算例

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

圖2 狹長矩形截面懸臂梁

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

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

圖3 梁體示意圖

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

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

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

5 結語

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

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

圖4 內部節(jié)點劃分

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

［1］ 趙廷剛.若干發(fā)展方程的譜方法和譜元法［D］.上海:上海大學數(shù)學系博士學位論文，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］.聲學學報，2007(32):525-533.

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