陳　斌，張延慶

(北京工業(yè)大學(xué)建筑工程學(xué)院，北京　100124)

?

廣義協(xié)調(diào)COONS曲面厚薄板通用單元振動分析

陳斌，張延慶

(北京工業(yè)大學(xué)建筑工程學(xué)院，北京100124)

摘要：厚薄板通用單元一直是板殼有限元研究中的重點和難點。根據(jù)單向厚板的基本方程推導(dǎo)得到內(nèi)含剪切變形因子的插值基函數(shù)，基于Mindlin-Reissner厚板理論并結(jié)合廣義協(xié)調(diào)思想，通過第2類有限元COONS曲面構(gòu)造法及撓度和轉(zhuǎn)角導(dǎo)數(shù)相關(guān)法構(gòu)造出一個新的有效的厚薄板通用矩形單元，采用集中質(zhì)量矩陣法得到單元的質(zhì)量矩陣。通過編制MATLAB程序，采用子空間迭代法計算結(jié)構(gòu)在不同板厚時的自振頻率。將計算結(jié)果與理論解和大型通用有限元軟件ABAQUS解進行比較，該單元具有良好的數(shù)值計算結(jié)果，能夠?qū)崿F(xiàn)動力學(xué)范圍內(nèi)薄板、中厚板以及厚板完全通用，具有工程應(yīng)用價值。

關(guān)鍵詞：結(jié)構(gòu)力學(xué)；自振頻率；厚薄板通用單元；COONS曲面；廣義協(xié)調(diào)；MATLAB

對于結(jié)構(gòu)計算而言，靜力問題分析雖是首要的，但是往往動力作用引起的破壞才是致命的，這是引起結(jié)構(gòu)嚴(yán)重破壞的主因。例如：地震作用引起的建筑物的坍塌；風(fēng)振作用引起的高塔、橋梁的振動破壞；碰撞或者是爆炸對結(jié)構(gòu)引起的沖擊破壞等等。因此，在工程結(jié)構(gòu)的研究、設(shè)計和安全性評估中，進行結(jié)構(gòu)的振動特性分析是相當(dāng)重要的，其中結(jié)構(gòu)振動頻率是非常重要和基礎(chǔ)的動力特性。當(dāng)前，板結(jié)構(gòu)是建筑結(jié)構(gòu)中最為常用也是最重要的結(jié)構(gòu)形式，因此對板結(jié)構(gòu)進行自由振動分析非常有必要。

由于板的特殊性，薄板與厚板的力學(xué)性能有所不同，對于薄板和厚板通常需分開單獨計算，實際應(yīng)用中對于薄板與厚板的劃分無明顯界限。為了提高計算效率和精度，人們開始探索厚薄板通用單元，而要實現(xiàn)厚薄板通用，則必須克服薄板情況下的剪切閉鎖問題。此外，對于厚薄板殼的振動問題，采用常規(guī)的解析法求解是很困難的，而采用有限元法求解則能得到很好的效果[1-7]。有限元COONS曲面法作為有限單元法與計算幾何學(xué)相結(jié)合的產(chǎn)物，構(gòu)造單元簡便高效，具有可以直接得到位移場函數(shù)顯式的優(yōu)點。

然而，到目前為止關(guān)于厚薄板通用的COONS曲面單元的振動分析研究尚少。文獻(xiàn)[2—3]通過Timoshenko深梁理論推導(dǎo)剪切變形因子并從薄板理論出發(fā)構(gòu)造出廣義協(xié)調(diào)COONS曲面混合單元，能夠?qū)崿F(xiàn)動力下的中厚薄板通用。由于基于薄板理論構(gòu)造，所以該單元缺乏轉(zhuǎn)角位移場和剪切應(yīng)變矩陣，對于厚跨比大的厚板(h/L>0.2時)，計算結(jié)果并不理想。因此，本文從Mindlin-Reissner厚板理論出發(fā)并結(jié)合廣義協(xié)調(diào)思想，采用第2類COONS曲面構(gòu)造法及撓度和轉(zhuǎn)角導(dǎo)數(shù)相關(guān)法構(gòu)造出一個新的厚薄板通用單元，通過編制MATLAB程序進行振動分析，為運用COONS曲面法研究板殼振動提供新的思路。

1單元構(gòu)造

彈性厚板理論(Mindlin-Reissner厚板理論)[8-9]假設(shè)原來垂直板中面的直線在變形后仍保持為直線，但由于考慮了剪切變形的影響，其不再垂直于變形后的板中面。因此，采用此理論的板單元撓度位移w和轉(zhuǎn)角位移θx，θy需各自獨立插值?；诤癜謇碚搧順?gòu)造厚薄板通用單元，關(guān)鍵在于如何在薄板情況下克服剪切閉鎖。為此，本文從單向厚板理論出發(fā)，推導(dǎo)出內(nèi)含考慮剪切變形影響的剪切變形因子λ的插值基函數(shù)。構(gòu)造單元時通過引入此插值基函數(shù)使板單元在薄板時的剪切剛度趨近于零，避免出現(xiàn)剪切閉鎖。

1.1剪切變形因子

如圖1所示，設(shè)單向厚板的撓度位移w為三次曲線，法線平均轉(zhuǎn)角θS為二次函數(shù)，則有式(1)：

(1)

受Timoshenko深梁理論推導(dǎo)的啟發(fā)，采用單向厚板理論推導(dǎo)得到[10-11]

(2)

圖1　單向厚板彎曲變形Fig.1　Bending deformation of one-way thick plate

將式(1)代入式(2)，積分得：

(3)

引入邊界條件w|S=0=wi,w|S=L=wj得：

C4=wi。

(4)

設(shè):

(5)

將式(4)和式(5)代入式(3)得單元撓度邊界

w=F0(t)wi+F1(t)wj+G0(t)θi+G1(t)θj。

(6)

將式(6)代入式(1)，得單元轉(zhuǎn)角邊界

θS=f0(t)wi+f1(t)wj+g0(t)θi+g1(t)θj。

(7)

式(6)、式(7)分別為單向厚板的撓度位移函數(shù)和轉(zhuǎn)角位移函數(shù)，均由物理方程推導(dǎo)得到。其中撓度位移插值基函數(shù)為

(8)

轉(zhuǎn)角位移插值基函數(shù)為

(9)

通過以上公式可以看出，當(dāng)板逐漸變薄時，即h/L→0時，λ→0，ρ→1時，

(10)

系數(shù)λ反映了剪切變形的影響，故稱為剪切變形因子，它是薄板COONS曲面單元與厚板COONS曲面單元相互聯(lián)系的橋梁。

綜上所述，當(dāng)h/L→0時，

θS→w′。

(11)

1.2單元撓度位移場

進行坐標(biāo)系轉(zhuǎn)換，無量綱自然坐標(biāo)系(ξ,η)與總體坐標(biāo)系(x,y)之間的關(guān)系為

ξ=x/a，η=y/b(0≤ξ,η≤1) ，

其中a和b為矩形單元的邊長。

擬合單元4條邊界的撓度位移場，并保證其在撓度和切向?qū)?shù)方向上保持協(xié)調(diào)：

[w(ξ,0)w(ξ,1)]=

(12)

(13)

通過計算表明：以線性插值基函數(shù)為混合函數(shù)的ACM元的計算精度較高，但收斂速度較慢，而以第1類Hermite插值基函數(shù)為混合函數(shù)的Ferguson曲面元收斂速度較快，但計算精度較差，據(jù)此本文對于單元邊界方向?qū)?shù)的擬合采用帶有剪切變形因子的線性組合函數(shù)為插值基函數(shù)進行擬合，這樣可以保證該單元在收斂速度和計算精度上有所提高。擬合結(jié)果如下：

[wη(ξ,0)wη(ξ,1)]=

(14)

(15)

其中L0=1-t，L1=t，t=ξ,η。

再引入廣義協(xié)調(diào)條件，使COONS曲面片邊界處的法線導(dǎo)數(shù)滿足廣義協(xié)調(diào)，即令

(16)

[Nw]=

[F0F00.5b(L0+F0)G0-0.5aG0(L0+F0)

F1F00.5b(L1+F1)G0-0.5aG1(L0+F0)

F1F10.5b(L1+F1)G1-0.5aG1(L1+F1)

F0F10.5b(L0+F0)G1-0.5aG0(L1+F1)]，

{δ}=[w1θx1θy1w2θx2θy2

w3θx3θy3w4θx4θy4]T。

[w]=[Nw]{δ}。

(17)

[Nw]的各項乘積因式中，前者是ξ的函數(shù)，后者是η的函數(shù)，前后順序不可交換(下同)。

1.3單元轉(zhuǎn)角位移場

為避免剪切閉鎖，必須使厚板元的撓度位移場和轉(zhuǎn)角位移場在h/L→0時滿足導(dǎo)數(shù)關(guān)系，即滿足Kirchhoff薄板理論的直法線假設(shè)。式(10)和式(11)已指出，當(dāng)h/L→0時，撓度位移場w和轉(zhuǎn)角位移場θx，θy滿足導(dǎo)數(shù)關(guān)系。而h/L不趨于零時，撓度場和轉(zhuǎn)角場并不滿足導(dǎo)數(shù)關(guān)系，符合厚板理論的要求。因此，采用撓度位移場的“偏導(dǎo)數(shù)”作為轉(zhuǎn)角位移場，但求導(dǎo)時附加約束性求導(dǎo)條件[12]，

由此，就能得到厚薄板通用單元的轉(zhuǎn)角位移場函數(shù)：

具體表達(dá)式分別如式(18)和式(19)所示。

(18)

(19)

1.4單元剛度矩陣

應(yīng)變有5個分量，即

{χ}=[χxχy2χxy…γxγy]T=

單元剛度矩陣：

[k]e=[kb]e+[kS]e。

[kb]e=∫A[Bb]T[Db][Bb]dA為單元彎曲剛度矩陣，[kS]e=∫A[BS]T[DS][BS]dA為單元剪切剛度矩陣。

顯然，當(dāng)h/L→0時，[BS]→[0]，[kS]e→[0]，故不會出現(xiàn)剪切閉鎖。

1.5單元質(zhì)量矩陣

采用集中質(zhì)量矩陣：

其中ne為單元節(jié)點數(shù)，ρ為材料密度。

1.6編程求解

基于MATLAB[13]強大的矩陣運算能力，本文采用MATLAB編制了結(jié)構(gòu)計算程序，并采用子空間迭代法求解結(jié)構(gòu)的自振頻率。計算中只需修改板厚，即可實現(xiàn)厚薄板的自振分析，無須像ABAQUS分析時需根據(jù)厚跨比選用不同的單元進行計算。

2數(shù)值算例

數(shù)值算例描述：四邊簡支各向同性矩形彈性方板，邊長為10 m，材料密度ρ=2 400 kg/m3，彈性模量E=3.0×1010N/m2，泊松比v=1/6。采用MATLAB編制的程序，按12×12劃分單元，求解不同厚跨比(h/L)下前5階頻率值，并與理論解和ABAQUS解[14]進行對照，結(jié)果如表1和圖2—圖6所示。

圖2　不同厚跨比一階頻率比較Fig.2　Comparison of the first order frequency 　　　by different thickness span ratio

圖3　不同厚跨比二階頻率比較Fig.3　Comparison of the second order frequency　　　by different thickness span ratio

階次頻率/Hz薄板h/L=0.001薄板h/L=0.01中厚板h/L=0.1中厚板h/L=0.15中厚板h/L=0.2中厚板h/L=0.25厚板h/L=0.3厚板h/L=0.35本文解0.32493.247430.70143.75655.35065.61074.66182.626(1,1)理論解0.32523.252031.44845.56557.96368.64077.59785.192ABAQUS解0.32713.268330.03742.71053.76763.33371.60278.763本文解0.81158.110075.116104.530128.290147.120161.930171.560(1,2)理論解0.81308.130175.451104.520126.860143.750156.180165.440ABAQUS解0.83348.326173.473101.450123.620141.030154.720165.540本文解1.295012.9370114.83154.520184.160206.180222.620234.980(2,2)理論解1.300813.0080116.04155.350182.970201.740214.800223.770ABAQUS解1.330912.2890111.590149.320177.230198.060204.070207.700本文解1.622216.2020143.500191.480225.190248.610265.040273.790(1,3)理論解1.626116.260141.320185.910215.390234.720247.550256.970ABAQUS解1.738517.353140.200186.430207.980210.850229.580231.740本文解2.100620.968177.750230.750266.200290.070306.490318.080(2,3)理論解2.114821.147177.620227.780259.350278.790291.450300.110ABAQUS解2.220122.150173.260208.730219.160245.850265.140278.400

注：表1中薄板的理論解參考文獻(xiàn)[15]，中厚板及厚板的理論解參考文獻(xiàn)[10]；ABAQUS模擬計算中，薄板采用S4R5單元，中厚板和厚板采用S8R單元。

圖4　不同厚跨比三階頻率比較Fig.4　Comparison of the three order frequency 　　　by different thickness span ratio

圖5　不同厚跨比四階頻率比較Fig.5　Comparison of the four order frequency　　　by different thickness span ratio

圖6　不同厚跨比五階頻率比較Fig.6　Comparison of the five order frequency　　　by different thickness span ratio

由以上表1和圖2—圖6可以看出：本單元在不同厚跨比下的前3階數(shù)值解基本與理論解吻合。隨著階數(shù)和板厚的增加，本文解的精度仍優(yōu)于ABAQUS解，能夠滿足工程計算需要。本單元能夠在動力學(xué)范圍內(nèi)實現(xiàn)薄板-中厚板-厚板的完全通用，能夠得到收斂速度快、計算精度較高的數(shù)值結(jié)果。

隨板厚和階數(shù)增加而出現(xiàn)的誤差與單元推導(dǎo)及程序編制時未考慮轉(zhuǎn)動慣量和擠壓變形的影響有關(guān)系。但相比于剪切變形的影響，以上的影響很小。為了提高精度，以后可在這方面進行相關(guān)研究。

3結(jié)論

基于Mindlin-Reissner厚板理論并結(jié)合廣義協(xié)調(diào)思想，采用第2類COONS曲面構(gòu)造法及撓度和轉(zhuǎn)角導(dǎo)數(shù)相關(guān)法構(gòu)造出一個新的厚薄板通用單元，通過自編MATLAB程序進行振動分析，數(shù)值算例表明該單元具有良好的收斂性和計算精度，能夠?qū)崿F(xiàn)動力學(xué)范圍內(nèi)薄板、中厚板及厚板的完全通用，具有工程實用性。

1)選用單向厚板基本方程推導(dǎo)出剪切變形因子，基于厚板理論構(gòu)造出的COONS曲面單元能夠?qū)崿F(xiàn)動力下的完全厚薄板通用，可作進一步的研究和推廣，以豐富大型有限元分析軟件的單元庫。

2)相比于第1類COONS曲面，第2類COONS曲面較復(fù)雜，自由度數(shù)也會增加，但收斂速度更快，計算精度更高。本文采用第2類COONS曲面構(gòu)造法，通過引入廣義協(xié)調(diào)條件消去自由度，從而得到消去節(jié)點扭矢的12參數(shù)矩形元，同時在單元邊界法向?qū)?shù)擬合中使用內(nèi)含剪切變形因子的線性組合函數(shù)作為插值基函數(shù)。數(shù)值算例結(jié)果表明，本單元具有良好的數(shù)值計算結(jié)果。

3)COONS曲面單元具有直接給出位移場顯式的優(yōu)點，已在靜力問題分析中取得了較好的收斂性和精度，然而在動力問題分析中缺少研究。本文通過編制的MATLAB程序，對推導(dǎo)出的基于厚板理論的厚薄板通用的COONS曲面單元進行振動特性分析。計算表明，COONS曲面單元在動力問題分析中同樣能夠?qū)崿F(xiàn)厚薄板通用，且具有良好的數(shù)值計算結(jié)果，為運用有限元COONS曲面法進行板殼的振動分析提供了借鑒。

參考文獻(xiàn)/References：

[1]張延慶.板殼有限元計算幾何理論[D]. 徐州:中國礦業(yè)大學(xué), 1990.

ZHANG Yanqing. Computational Geometric Theory of Finite Element of Plate and Shell[D]. Xuzhou: China University of Mining and Technology，1990.

[2]張延慶，龍馭球. 厚薄板通用的廣義協(xié)調(diào)單元COONS曲面構(gòu)造法[J]. 工程力學(xué), 1994, 11(2)：59-64.

ZHANG Yanqing, LONG Yuqiu. Formulation of a generalized conforming element for moderate thick plate based on COONS’ surface[J]. Engineering Mechanics，1994，11(2)：59-64.

[3]陳葉妮，張延慶. 用廣義協(xié)調(diào)COONS混合曲面元分析中厚板的振動[J]. 工業(yè)建筑, 2006(sup1)：612-614.

CHEN Yeni，ZHANG Yanqing. Vibration analysis of moderate thick plate by generalized conforming element and COONS surface approach[J]. Industrial Construction，2006(sup1)：612-614.

[4]曹金鳳，姜耀東，趙國景，等. 廣義協(xié)調(diào)厚薄板通用單元TMT的固有振動分析[J]. 中國礦業(yè)大學(xué)學(xué)報，2007, 36(1)：91-96.

CAO Jinfeng，JIANG Yaodong，ZHAO Guojing，et al. Free vibration analysis of plate using a generalized conforming thick/thin element TMT[J]. Journal of China University of Mining and Technology，2007, 36(1)：91-96.

[5]曹金鳳，姜耀東，韓志茹，等. 廣義協(xié)調(diào)厚薄板通用單元的固有振動分析[J]. 力學(xué)與實踐，2007，29(3)：30-35.

CAO Jinfeng，JIANG Yaodong，HAN Zhiru，et al. Free vibration analysis for plates using a generalized conforming thick/thin triangular element[J]. Journal of Mechanics in Engineering，2007, 29(3)：30-35.

[6]韓志茹，姜耀東，王麗，等. 廣義協(xié)調(diào)厚薄通用板單元的動力特性分析[J]. 力學(xué)與實踐，2009, 31(3)：69-72.

HAN Zhiru，JIANG Yaodong，WANG Li，et al. The dynamic property of generalized conforming thick/thin plate elements[J]. Journal of Mechanics in Engineering，2009, 31(3)：69-72.

[7]黃修長，錢振華，李俊，等. 應(yīng)用基于解析式函數(shù)的廣義協(xié)調(diào)四邊形厚板元分析中厚板的自由振動[J]. 工程力學(xué)，2011, 28(9)：39-43.

HUANG Xiuchang，QIAN Zhenhua，LI Jun，et al. Vibration analysis of moderate thick plate by generalized conforming quadrilateral thick plate element based on analytical trial functions[J]. Engineering Mechanics，2011, 28(9)：39-43.

[8]MINDLIN R D. Influence of rotatory inertia and shear on flexural motions of isotropic elastic plates[J]. ASME Journal of Applied Mechanics，1951, 18(1)：31-38.

[9]REISSNER E. On bending of elastic plates[J]. Quarterly of Applied Mathematics，1947, 5(1)：55-68.

[10]曹志遠(yuǎn)，楊昇田. 厚板動力學(xué)理論及其應(yīng)用[M]. 北京:科學(xué)出版社，1983：57-124.

[11]鄭照北，王宏凱. 板殼有限元COONS曲面法[M]. 徐州:中國礦業(yè)大學(xué)出版社，1993：52-73.

[12]王宏凱，鄭照北. 一個有效的厚薄板通用COONS曲面矩形元[J]. 工程力學(xué)，1992, 9(1)：134-140.

WANG Hongkai，ZHENG Zhaobei. An efficacious COONS’ surface method for formulation of current finite element of thick and thin plates[J]. Engineering Mechanics，1992, 9(1)：134-140.

[13]劉衛(wèi)國. MATLAB程序設(shè)計教程[M].第2版. 北京：中國水利水電出版社，2010.

[14]莊茁，曲小川，廖劍暉, 等. 基于ABAQUS的有限元分析和應(yīng)用[M]. 北京：清華大學(xué)出版社，2009.

[15]曹志遠(yuǎn). 板殼振動理論[M]. 北京: 中國鐵道出版社, 1989：32-54.

Vibration analysis of generalized conforming and COONS’surface thick/thin plate element

CHEN Bin, ZHANG Yanqing

(College of Architecture and Civil Engineering, Beijing University of Technology, Beijing 100124, China)

Abstract:Thick/thin plate element has been the focus and difficulty in the study of shell finite element. According to the basic equation of one-way thick plate, the base function of interpolation containing the shear deformation factor is derived. Based on the Mindlin-Reissner thick plate theory and the idea of generalized coordination, the second kind of finite element COONS’ surface construction method and the derivative correlation method of deflection displacement and angular rotation displacement are used to derive a new effective thick/thin general rectangular plate element, then the lumped mass matrix method is used to get the mass matrix of the element. By using MATLAB programming, the natural frequencies of the plates with different thickness are acquired by a sub-space iteration method. Compared with the theoretical solution and the numerical value of ABAQUS, it is concluded that this element has good numerical calculation results, which can realize entirely common use for thin plate, cut deal and thick plate in dynamic range, and it has engineering application value.

Keywords:structural mechanics; natural frequency; thick/thin plate element; COONS’ surface; generalized conforming; MATLAB

文章編號：1008-1534(2016)03-0224-06

收稿日期：2016-04-04;修回日期：2016-04-24；責(zé)任編輯：馮民

作者簡介：陳斌(1990—)，男，浙江金華人，碩士研究生，主要從事工程結(jié)構(gòu)數(shù)值分析與優(yōu)化方面的研究。通訊作者：張延慶教授。E-mail:zhyq@bjut.edu.cn

中圖分類號：TU311.4

文獻(xiàn)標(biāo)志碼：A

doi:10.7535/hbgykj.2016yx03008

陳斌，張延慶.廣義協(xié)調(diào)COONS曲面厚薄板通用單元振動分析[J].河北工業(yè)科技,2016,33(3):224-229.

CHEN Bin, ZHANG Yanqing.Vibration analysis of generalized conforming and COONS’ surface thick/thin plate element[J].Hebei Journal of Industrial Science and Technology,2016,33(3):224-229.