鄭 炎， 肖永建， 許震宇

(西藏職業(yè)技術(shù)學(xué)院 建筑系， 西藏 拉薩 850030)

在很多工程問題中，常常要求求解廣義代數(shù)特征值問題。而目前的有限單元法(FEM)、有限差分法(FDM)以及邊界單元法(BEM)等數(shù)值方法在處理復(fù)雜幾何形狀的工程的本征值問題時(shí)，網(wǎng)格很難劃分，并且劃分網(wǎng)格很花費(fèi)時(shí)間。近年來，無網(wǎng)格法越來越受到學(xué)者們的重視，其最初的思想可以追溯到用光滑粒子法(SPH)模擬天體物理學(xué)現(xiàn)象[1]，許多學(xué)者致力于用無網(wǎng)格方法求解偏微分方程(PDEs)的數(shù)值解。目前，無網(wǎng)格方法比較多，例如，無網(wǎng)格伽遼金法(EFGM)[2]，重構(gòu)核點(diǎn)法(RKM)[3]，基本解方法(MFS)。

本文采用一種簡單的由Chen和Tanaka[4，5]提出的無網(wǎng)格方法——邊界點(diǎn)法(BKM)求解各種形狀的板在各種邊界條件下的固有頻率。BKM是基于徑向基函數(shù)(RBF)方法的一種，由于徑向基函數(shù)方法不采用移動(dòng)最小二乘法，用一維的距離變量來構(gòu)造多元空間變量，不用劃分網(wǎng)格。而BKM繼承徑向基函數(shù)的無網(wǎng)格特性，與維數(shù)和幾何形狀無關(guān)，是完全的無網(wǎng)格方法，不需積分，容易學(xué)習(xí)和編程。Chen和Hon[6]從數(shù)學(xué)上說明了BKM是譜收斂的，譜收斂就是說近似程度只和選擇的物理解的光滑程度有關(guān)系。

Kang[7，8]和Chen[9]分別用無網(wǎng)格方法給出了夾支板和膜的自由振動(dòng)的固有頻率的數(shù)值解，然而對(duì)于計(jì)算復(fù)雜的高階邊界的板的自有振動(dòng)的固有頻率的數(shù)值解仍然有一定的難度。Kang[7]將勁度矩陣有2N×2N維降到N×N，這種方法能夠加快計(jì)算速度，節(jié)約計(jì)算機(jī)內(nèi)存空間，但降維的過程中有矩陣逆的操作，使矩陣的條件數(shù)變差，容易使勁度矩陣奇異，最終無法求解。另外，在處理多邊形的角點(diǎn)的法向量時(shí)， Kang簡單的將相鄰兩邊的該點(diǎn)的法向量的矢量和代替該點(diǎn)的法向量，其法向量大于單位1。在計(jì)算夾支板和膜的固有頻率時(shí)，是可行的，但是計(jì)算簡支邊和自由邊時(shí)，誤差很大，無法滿足精度要求。本文的角點(diǎn)法向量也采用兩邊法向量的矢量和的方向，但是法向量的大小采用單位1，這樣更合理，數(shù)值解也更精確。

1 板自由振動(dòng)的控制方程

各向同性的薄板的自由振動(dòng)方程為

w(r,t)=eiωtW(r)

(2)

其中，ω表示圓頻率,W(r)為板的振型函數(shù), 將方程(2)代入方程(1)有

(3)

(4)

其中k為板的頻率參數(shù)。方程(3)的通解[10]為

W(r)=AjJ0(k|r-rj|)+BjI0(k|r-rj|)

(5)

其中，J0(r),I0(r)分別為第一類和第一類修正的貝塞爾函數(shù)，|r-rj|為邊界點(diǎn)之間的距離(圖1)，r為邊界點(diǎn)的矢量,AJ,Bj為未知系數(shù)。

nk，τs分別表示點(diǎn)k，s處的法向量和切向量圖1 任意形狀的薄板的邊界點(diǎn)的分布

如圖1所示，在薄板的邊界取N個(gè)點(diǎn)，任意一點(diǎn)的模態(tài)可以表示為

(6)

2 邊界條件

薄板邊界上的位移為

u=W(r)

(7)

薄板邊界沿法向量方向的轉(zhuǎn)角為

(8)

邊界上的彎矩為

(9)

邊界上的剪力為

(10)

相鄰自由邊的交點(diǎn)處的集中反力

(11)

式(11)中，u,θ,M,V分別表示邊界的位移、轉(zhuǎn)角、彎矩和剪力，R為相鄰兩自由邊交點(diǎn)處的集中反力；v為泊松比；n,τ分別表示邊界的法向量和切向量。對(duì)于各種邊界條件的板，只需滿足上述式(7)～(10)四個(gè)邊界條件中的兩個(gè)(有兩相鄰自由邊還需方程(11)等于零就可以。需要強(qiáng)調(diào)的是，對(duì)于多邊形的相鄰兩邊的交點(diǎn)的法向量的法向?yàn)樵擖c(diǎn)所屬的兩邊該點(diǎn)的法向量的矢量和的方向，但大小為單位1。

將式(6)代入式(11)，并使表達(dá)式等于零，可得邊界上任意一點(diǎn)的邊界條件方程，即

i=1,2,…,N

(12)

i=1,2,…,N

(13)

i=1,2,…,N

(14)

i=1,2,…,N

(15)

i=1,2,…,N

(16)

其中,rij=|ri-rj| 分別表示邊界上第i個(gè)點(diǎn)處的法向量和切向量，如圖1所示。

方程(12)～(16)可以用矩陣形式簡單的表示為

(17)

方程(17)左邊的勁度矩陣的部分公式推導(dǎo)如下

(18)

(19)

方程(17)左邊的勁度矩陣的其余公式的詳細(xì)推導(dǎo)結(jié)果略。式(18), (19)中，nix，niy和τix,τiy分別表示單位法向量ni=nix+niy和單位切向量τi=τix+τiy在第i個(gè)邊界點(diǎn)處的水平分量和豎直分量,θij為兩點(diǎn)矢量ri-rj方向與x軸的夾角。

(20)

(21)

=0

(22)

可根據(jù)方程(22)求出各種邊界的板的頻率參數(shù)k。

3 算 例

3.1 方形板

三邊簡支一邊自由的方形板(圖2(a))被分成4N-4(每邊N個(gè))個(gè)點(diǎn)(圖2(b))。部分點(diǎn)和

其相應(yīng)的法線方向如圖2(b)所示。方板的四個(gè)角的法線方向采用前述的處理方法處理。

圖2 (a)三邊簡支一邊自由的方形板(其中S-S表示簡支邊，F(xiàn)-F表示自由邊)(b)各邊分成N等分時(shí)方形板的離散邊界點(diǎn)和法線方向，ni表示第i節(jié)點(diǎn)的法線方向

從計(jì)算三邊簡支一邊自由的方形板的結(jié)果可以看出，其包含有四邊簡支方形板的頻率參數(shù)。因此，必須剔出這些不屬于三邊簡支一邊自由的方形板的頻率參數(shù)。為簡單起見，用SSSS(k)表示四邊簡支方板的的勁度矩陣，SSSF(k)、SF(k)分別表示三邊簡支一邊自由方板的勁度矩陣和其特有頻率的勁度矩陣。則三邊簡支一邊自由方板的特有頻率參數(shù)矩陣可以表示為

SK(k)=SSSF(k)[SSSS(k)]-1

(23)

在圖3中，三邊簡支一邊自由的方形板的頻率參數(shù)能從det[SF(k)]=0的對(duì)數(shù)曲線中得到。波峰表示四邊簡支方板的頻率參數(shù)，波谷表示三邊簡支一邊自由的方板的頻率參數(shù)，表1和表2分別給出了對(duì)應(yīng)的頻率參數(shù)值和與精確值的相對(duì)誤差。

圖3 N=6, ν=0.3時(shí)，矩形板的對(duì)數(shù)曲線det[SF(k)]，SF表示三邊簡支一邊自由板頻率參數(shù)，SS表示四邊簡支的頻率參數(shù)

參數(shù)計(jì)算值20(N=6)相對(duì)誤差24(N=7)相對(duì)誤差28(N=8)相對(duì)誤差精確解[11]SF14.44290%----4.4429SF27.02480%7.02480%7.02480%7.0248SF38.88580%8.88580%8.88580%8.8858SF49.9340-0.006%9.93460%9.93460%9.9346

表2 三邊簡支一邊自由方形板的頻率參數(shù)值

從表1和表2中可知，20個(gè)邊界點(diǎn)得到的頻率參數(shù)值已經(jīng)具有很高的精度，如果要得到更高階參數(shù)，須取更多的邊界點(diǎn)。采用同樣的方法，可以計(jì)算任意邊界條件下的方形板頻率參數(shù)值。

3.2 任意形狀的簡支板

圖4 (a) 任意形狀簡支板，S-S表示簡支邊，(b) 各邊分成N等分時(shí)任意形狀簡支板的離散邊界點(diǎn)的法線方向和切線方向，ni，τi表示第i節(jié)點(diǎn)的法線方向和切線方向

圖4 (a)中，一個(gè)由半徑為1 m的半圓和兩個(gè)邊長為1.5 m和2.5 m的邊組成的域?yàn)槿我庑螤畎宓男螤?，在其邊界?N-4點(diǎn)，每個(gè)點(diǎn)對(duì)應(yīng)的法線方向和切線方向如圖4 (b)所示。各邊的交點(diǎn)處的法向量和切向量處理如前文所述。

在計(jì)算任意形狀的簡支板頻率參數(shù)時(shí)，從計(jì)算結(jié)果的圖形中可以看出，其迭代出的結(jié)果含有對(duì)應(yīng)形狀的膜的頻率參數(shù)。因此，必須剔除這部分不屬于任意形狀簡支板的頻率參數(shù)。本文中，簡支板的頻率參數(shù)矩陣為(8N-8)×(8N-8)維，而對(duì)應(yīng)形狀的膜的頻率參數(shù)矩陣為(4N-4)×(4N-4)維，剔除這些非本征頻率參數(shù)必須降維。

(24)

方程組(24)中的第一個(gè)方程可以寫成

(25)

將方程(25)代入方程組(24)中第二個(gè)方程

(26)

膜的頻率參數(shù)[7]方程為

(27)

則任意形狀的簡支板的特征頻率參數(shù)方程為

(28)

方程(28)系數(shù)矩陣[A]有非零解，則其勁度矩陣的行列式為零

(29)

計(jì)算結(jié)果如圖5和表3所示。圖5為20個(gè)邊界點(diǎn)的計(jì)算結(jié)果，其中AS,M表示任意形狀簡支板和膜的頻率參數(shù)，6個(gè)波峰和7各波谷分別為det[AS(k)]的對(duì)數(shù)曲線對(duì)應(yīng)的任意形狀的膜和簡支板的頻率參數(shù)值。表3為20、30個(gè)邊界點(diǎn)的計(jì)算結(jié)果與ANSYS軟件計(jì)算結(jié)果的比較，由表中可以看出，計(jì)算精度很高。

表3 BKM法計(jì)算任意形狀簡支板的特征值

圖5 當(dāng) N=6, ν=0.3時(shí)，矩形板的對(duì)數(shù)曲線det[FSSS(k)]，SS表示三邊簡支一邊自由板頻率參數(shù)，M表示三邊簡支一邊自由膜的頻率參數(shù)

4 結(jié) 語

BKM法可視為一種Trefftz法，其實(shí)驗(yàn)功能要求測(cè)試函數(shù)滿足控制方程。“BKM法”可以說是一種結(jié)合DRM、RBF和非奇異一般解的邊界建模技術(shù)。

[1]Gingold R A, Maraghan J J. Smoothed particle hydrodynamics: theory and applications to non-spherical stars[J]. Monthly Notices of the Royal Astronomical Society, 1977, 181(2): 375-389.

[2]Lancaster P,Salkauskas K.Surfaces generated by moving least-squares methods[J]. Mathematics of Computation, 1981, 37(155): 141-158.

[3]Liu WK, Jun S. Multiple-scale Reproducing Kernel Particle Methods for large deformation problems [J]. International Journal for Numerical Methods in Engineering, 1998, 41(7): 1339-1362.

[4]Chen W, Tanaka M. A meshless, integration-free, and boundary-only RBF technique [J]. Computers and Mathematics with Applications , 2002, 43(3-5): 379-391.

[5]Chen W, Tanaka M. New insights into boundary-only and domain-type rbf methods[J]. International Journal of Nonlinear Sciences and Numerical Simulation, 2000, 1(3): 145-151.

[6]Chen W, Hon YC. Numerical convergence of boundary knot method in the analysis of helmholtz, modified helmholtz, and convection-diffusion problems[J]. Computer Methods in Applied Mechanics and Engineering, 2003, 192(15): 1859-1875.

[7]Kang S W, Lee J M, Kang Y J. Vibration analysis of arbitrarily shaped membranes using non-dimensional dynamic influence function [J]. Journal of Sound and Vibration, 1999, 221(1): 117-132.

[8]Kang S W, Lee J M. Free vibration analysis of arbitrarily shaped plates with clamped edges using wave-type functions [J]. Journal of Sound and Vibration, 2001, 241(1): 9-26.

[9]Chen J T, Chen I L, Chen K H, et al. A meshless method for free vibration analysis of circular and rectangular clamped plates using radial basis function [J]. Engineering Analysis with Boundary Elements, 2004, 28(5): 535-545.

[10]Leissa A W. The Free Vibration of Rectangular Plates[J]. Journal of Sound and Vibration, 1973, 31(3): 257-293.