三維20 節(jié)點(diǎn)六面體和10 節(jié)點(diǎn)四面體單元的高精度中節(jié)點(diǎn)集中質(zhì)量矩陣1)

侯松陽 王東東 吳振宇 林智煒

(廈門大學(xué)土木工程系,福建省濱海土木工程數(shù)字仿真重點(diǎn)實(shí)驗(yàn)室,福建廈門 361005)

引言

在結(jié)構(gòu)動(dòng)力分析[1-6]中,質(zhì)量矩陣的構(gòu)造形式對(duì)有限元計(jì)算精度具有顯著的影響.其中,常用的有限元質(zhì)量矩陣構(gòu)造形式主要包括一致質(zhì)量矩陣[1,7-8]、高階質(zhì)量矩陣[1,9-11]和集中質(zhì)量矩陣[12-19].其中,集中質(zhì)量矩陣具有構(gòu)造形式簡單、存儲(chǔ)空間小、計(jì)算效率高和可與中心差分法等顯式積分方案結(jié)合使用等優(yōu)點(diǎn).而常用的集中質(zhì)量矩陣構(gòu)造方法主要包括行求和法[1,13]、節(jié)點(diǎn)積分法[14-16]和一致質(zhì)量矩陣主對(duì)角元素放大法(HRZ 法)[17]等.在集中質(zhì)量矩陣的精度分析方面,Ainsworth 等[16]給出拉格朗日單元采用節(jié)點(diǎn)積分方案求解標(biāo)量波特征值問題時(shí)的精度度量理論.Yang 等[18]針對(duì)8 節(jié)點(diǎn)平面單元提出一種基于數(shù)值微分流形的集中質(zhì)量構(gòu)造方法.Duczek等[19]研究了Lobatto 單元行求和、節(jié)點(diǎn)積分和HRZ 集中質(zhì)量矩陣之間的等價(jià)性.最近,Li 等[20]研究了采用行求和集中質(zhì)量矩陣時(shí),拉格朗日單元節(jié)點(diǎn)布置對(duì)頻率計(jì)算精度和收斂階次的影響,詳細(xì)分析了節(jié)點(diǎn)均勻分布單元集中質(zhì)量矩陣的頻率精度受限性.

在有限元分析領(lǐng)域,20 節(jié)點(diǎn)六面體單元和10 節(jié)點(diǎn)四面體單元應(yīng)用非常廣泛[21-25].與27 節(jié)點(diǎn)六面體拉格朗日單元相比,20 節(jié)點(diǎn)六面體單元不含內(nèi)部節(jié)點(diǎn),充分利用了單元間的節(jié)點(diǎn)共享特性,可顯著縮減計(jì)算模型的帶寬,提高計(jì)算效率.與此同時(shí),10 節(jié)點(diǎn)四面體單元在復(fù)雜形狀網(wǎng)格剖分方面具有明顯優(yōu)勢(shì).然而,當(dāng)采用行求和方法構(gòu)造這兩類單元的集中質(zhì)量矩陣時(shí),主對(duì)角線上均出現(xiàn)負(fù)質(zhì)量元素,難以直接用于結(jié)構(gòu)動(dòng)力分析.目前,最常用的消除集中矩陣負(fù)質(zhì)量元素的有效方法是HRZ 法[17],但是三維情況下該方法缺乏理論層面的精度分析.針對(duì)8 節(jié)點(diǎn)平面單元的HRZ 集中質(zhì)量矩陣,Hou 等[26]建立了相應(yīng)的頻率精度表達(dá)式,通過分析證明HRZ 集中質(zhì)量矩陣并不能提供最優(yōu)的頻率計(jì)算精度,進(jìn)一步發(fā)展了具有更優(yōu)精度的8 節(jié)點(diǎn)平面單元中節(jié)點(diǎn)集中質(zhì)量矩陣.但是,該研究尚未考慮更為復(fù)雜的三維單元.另外,由于中節(jié)點(diǎn)集中質(zhì)量矩陣的角節(jié)點(diǎn)質(zhì)量為0[26],其對(duì)時(shí)程動(dòng)力分析的影響也需要進(jìn)一步厘清.

本文針對(duì)三維20 節(jié)點(diǎn)六面體和10 節(jié)點(diǎn)四面體單元行求和集中質(zhì)量矩陣出現(xiàn)負(fù)質(zhì)量的問題,提出了三維廣義中節(jié)點(diǎn)集中質(zhì)量矩陣的構(gòu)造形式,并將HRZ 集中質(zhì)量矩陣作為特例涵蓋其中.然后,以20 節(jié)點(diǎn)六面體單元為例,構(gòu)建了三維廣義中節(jié)點(diǎn)集中質(zhì)量矩陣的頻率計(jì)算精度表達(dá)式,進(jìn)而分析了三維HRZ 集中質(zhì)量矩陣的頻率精度.根據(jù)三維廣義中節(jié)點(diǎn)集中質(zhì)量矩陣的頻率精度表達(dá)式,通過對(duì)其中的待定參數(shù)進(jìn)行優(yōu)化,構(gòu)造了20 節(jié)點(diǎn)六面體單元的高精度集中質(zhì)量矩陣.接著,利用中節(jié)點(diǎn)集中質(zhì)量矩陣構(gòu)造形式簡單的特點(diǎn),將其推廣到10 節(jié)點(diǎn)四面體單元.最后,對(duì)于時(shí)程動(dòng)力分析,通過靜力凝聚消除了中節(jié)點(diǎn)集中質(zhì)量矩陣的角節(jié)點(diǎn)零質(zhì)量影響,構(gòu)造了三維有限元中節(jié)點(diǎn)集中質(zhì)量矩陣的降階動(dòng)力計(jì)算模型,在保證動(dòng)力計(jì)算精度的同時(shí)提高了計(jì)算效率.

1 結(jié)構(gòu)動(dòng)力分析的有限元離散方程

不失一般性,考慮如下的經(jīng)典波動(dòng)方程等效積分弱形式[1]

其中,δ 表示變分符號(hào),上標(biāo)T 代表轉(zhuǎn)置.對(duì)于膜或聲場(chǎng)問題,場(chǎng)變量u退化為標(biāo)量u,表示膜的橫向位移或聲壓,外力或源項(xiàng)b退化為標(biāo)量b.與之對(duì)應(yīng)的梯度矩陣和材料矩陣D分別為

其中,c為波速.對(duì)于彈性力學(xué)問題,矢量場(chǎng)u={ux,uy,uz}T表示彈性體內(nèi)一點(diǎn)的位移,b={bx,by,bz}T表示體力,相應(yīng)的梯度矩陣和材料彈性矩陣D分別為

其中,cp和cs表示壓力波和剪力波的波速[11]

其中,λ 和 μ 為拉梅常數(shù),ρ 為材料密度.

其中,NI(ξ) 表示單元第I個(gè)節(jié)點(diǎn)的形函數(shù),nen為單元節(jié)點(diǎn)個(gè)數(shù).將式(6)代入式(1)中,可得結(jié)構(gòu)動(dòng)力分析的有限元離散方程

其中,M和K分別為整體質(zhì)量矩陣和整體剛度矩陣,f為力向量,可分別由單元質(zhì)量矩陣Me、剛度矩陣Ke和力向量fe通過組裝算子 A 組裝得到.Me,Ke和fe的定義為

其中,1 是單位矩陣,對(duì)于聲場(chǎng)問題,1=1;對(duì)于彈性連續(xù)體問題,1=diag{1,1,1}.

對(duì)于自由振動(dòng)分析,在式(7)中忽略外力項(xiàng),并引入如下的簡諧波假定

其中,ωh為數(shù)值頻率,? 為與之對(duì)應(yīng)的振動(dòng)模態(tài).將式(12)代入式(7),可得自由振動(dòng)分析的有限元離散方程

2 三維20 節(jié)點(diǎn)六面體和10 節(jié)點(diǎn)四面體單元的高精度集中質(zhì)量矩陣構(gòu)造方法

為表述簡潔起見,在下面的討論中,分別用H20和T10 表示三維20 節(jié)點(diǎn)六面體和10 節(jié)點(diǎn)四面體單元.

2.1 H20 和T10 單元的HRZ 集中質(zhì)量矩陣

由于滿足變分一致性,式(9)定義的質(zhì)量矩陣通常被稱為一致質(zhì)量矩陣.當(dāng)采用行求和法構(gòu)造集中質(zhì)量矩陣時(shí),僅需將一致質(zhì)量矩陣每行的元素累加到主對(duì)角元素.因此,對(duì)于標(biāo)量波動(dòng)問題,考慮規(guī)則的H20 和T10 單元構(gòu)形,例如長方體和直邊四面體,其行求和集中質(zhì)量矩陣分別為

其中,下標(biāo) R S 表示行求和法;me為單元的總質(zhì)量.為清晰起見,圖1 和圖2 給出了H20 和T10 單元的節(jié)點(diǎn)質(zhì)量分布,其中MNLM (mid-node lumped mass matrix)表示中節(jié)點(diǎn)集中質(zhì)量矩陣.

圖1 H20 單元3 種集中質(zhì)量矩陣的節(jié)點(diǎn)質(zhì)量分布示意圖Fig.1 Schematic illustration of the nodal mass distribution for three mass lumping schemes of H20 element

圖2 T10 單元3 種集中質(zhì)量矩陣的節(jié)點(diǎn)質(zhì)量分布示意圖Fig.2 Schematic illustration of the nodal mass distribution for three mass lumping schemes of T10 element

由式(14)和式(15)可得,對(duì)于H20 和T10 單元,行求和集中質(zhì)量矩陣中的角節(jié)點(diǎn)質(zhì)量均為負(fù)值.根據(jù)HRZ 方法,將式(9)中的一致質(zhì)量矩陣的主對(duì)角元素在滿足質(zhì)量守恒的條件下進(jìn)行比例縮放,即可得到非負(fù)的HRZ 集中質(zhì)量矩陣

其中,r為縮放系數(shù)

由式(16) 和式(17),可得H20 和T10 單元的HRZ 集中質(zhì)量矩陣

圖1 和圖2 也給出了這兩種單元的HRZ 集中質(zhì)量矩陣節(jié)點(diǎn)質(zhì)量分布情況.

2.2 H20 和T10 單元的非負(fù)中節(jié)點(diǎn)集中質(zhì)量矩陣

如前所述,H20 和T10 單元的行求和集中質(zhì)量矩陣,負(fù)質(zhì)量都出現(xiàn)在角節(jié)點(diǎn)上,而中節(jié)點(diǎn)則沒有負(fù)質(zhì)量問題.因此,可以基于行求和集中質(zhì)量矩陣中的非負(fù)質(zhì)量元素,通過單元質(zhì)量守恒將其進(jìn)行比例放縮,同時(shí)將行求和集中質(zhì)量矩陣中的負(fù)對(duì)角元素置零,構(gòu)造一種新型非負(fù)集中質(zhì)量矩陣.該集中質(zhì)量矩陣的構(gòu)造流程為

基于式(20)和式(21),H20 和T10 單元的新型集中質(zhì)量矩陣分別為

由式(22)和式(23)可見,對(duì)于H20 和T10 單元,角節(jié)點(diǎn)的質(zhì)量均為零,質(zhì)量只被分配在中節(jié)點(diǎn)上,所以把該集中質(zhì)量矩陣稱之為中節(jié)點(diǎn)集中質(zhì)量矩陣,簡稱MNLM 集中質(zhì)量矩陣.同樣,圖1 和圖2 中也描述了H20 和T10 單元的MNLM 集中質(zhì)量矩陣.

3 H20 單元集中質(zhì)量矩陣的頻率精度分析

本節(jié)以H20 單元為例,研究集中質(zhì)量矩陣構(gòu)造形式對(duì)頻率計(jì)算精度的影響.為了更全面衡量各類集中質(zhì)量矩陣的精度,構(gòu)造如下H20 單元廣義集中質(zhì)量矩陣

其中,α 和 β 是待定參數(shù).為了保證集中質(zhì)量矩陣的非負(fù)性,角節(jié)點(diǎn)質(zhì)量 α 需滿足 α ∈[0,1].當(dāng) α=7/31 時(shí),式(24)中的廣義集中質(zhì)量矩陣即退化為式(18)的HRZ 集中質(zhì)量矩陣;而當(dāng) α=0 時(shí),即為式(22) 的MNLM 集中質(zhì)量矩陣.

為了進(jìn)行頻率精度分析,考慮圖3 所示的H20單元典型節(jié)點(diǎn)影響域.由H20 單元構(gòu)成的有限元典型節(jié)點(diǎn)包括一個(gè)角節(jié)點(diǎn)xabc和3 個(gè)中間節(jié)點(diǎn)xa(b+1)c,x(a-1)bc和xab(c+1).θ 和 φ 是與諧波相關(guān)的不同方向的入射角,hi為xi方向的單元長度,不同方向的波數(shù)ki與k之間的關(guān)系滿足

圖3 H20 單元典型節(jié)點(diǎn)影響域及波動(dòng)方向示意圖Fig.3 Illustration of the typical nodal influence domains and wave propagation angles for a mesh with H20 elements

假設(shè)離散節(jié)點(diǎn)的波動(dòng)形式為諧波,4 個(gè)典型節(jié)點(diǎn)系數(shù)可以表示為

方便表述起見,引入算子cos(±lkx±m(xù)ky±nkz)

當(dāng)式(28)僅包括兩個(gè)下標(biāo)時(shí),有

對(duì)于4 個(gè)典型節(jié)點(diǎn),將式(26)和式(27)代入式(7),并利用簡諧波假定化簡可得

其中,AJI=AIJ,KI,J=KIJ為剛度矩陣元素,為了不造成下標(biāo)混淆,下標(biāo)的行號(hào)和列號(hào)之間采用“,”隔開.

由于式(30)有非0 解,則矩陣A的行列式應(yīng)為0

值得指出的是,式(42)是關(guān)于數(shù)值頻率 ωh的非線性方程,但由于其復(fù)雜性,從中直接求解 ωh通常是非常困難的.另一方面,在理論分析中我們關(guān)注的是頻率誤差,而非頻率本身,因而可直接引入頻率誤差e[1]

結(jié)合解析頻率 ω=kc和式(43),有

將式(44)代入式(42),并忽略e的高次項(xiàng)[26],可得

進(jìn)一步在式(45)中引入余弦項(xiàng)關(guān)于kh的泰勒展開,可得H20 單元廣義集中質(zhì)量矩陣的頻率誤差表達(dá)式

其中h為特征單元長度

將 α=7/31 代入式(48),即得H20 單元HRZ 集中質(zhì)量矩陣的頻率誤差表達(dá)式

對(duì)比式(48)和式(50),可見H20 單元的HRZ集中質(zhì)量矩陣的精度并非最優(yōu).反之,由式(48)可知,α=0 對(duì)應(yīng)的H20 單元集中質(zhì)量矩陣具有更優(yōu)的頻率計(jì)算精度,此時(shí)的集中質(zhì)量矩陣即為MNLM 集中質(zhì)量矩陣.進(jìn)一步,將 α=0 代入式(48) 中,即得MNLM 集中質(zhì)量矩陣的頻率誤差表達(dá)式

基于式(50)和式(51),若忽略頻率誤差的高階項(xiàng),可得三維H20 單元MNLM 集中質(zhì)量矩陣與HRZ集中質(zhì)量矩陣兩者頻率誤差之間存在如下關(guān)系

4 基于中節(jié)點(diǎn)集中質(zhì)量矩陣的時(shí)程分析

對(duì)于H20 和T10 單元,式(22)和式(23)表明MNLM集中質(zhì)量矩陣的角節(jié)點(diǎn)質(zhì)量為0.利用這一特點(diǎn),可通過靜力凝聚建立相應(yīng)的結(jié)構(gòu)動(dòng)力分析降階模型[27],減小整體計(jì)算規(guī)模,提高計(jì)算效率.

將H20 或T10 單元的MNLM 集中質(zhì)量矩陣代入式(7)的有限元離散運(yùn)動(dòng)方程,同時(shí)將零質(zhì)量角節(jié)點(diǎn)和非零質(zhì)量中節(jié)點(diǎn)進(jìn)行分類,可得

其中,下標(biāo)C和M分別表示單元的角節(jié)點(diǎn)和中節(jié)點(diǎn),MMM,KCC和KMM均為正定矩陣.對(duì)式(53)中的零質(zhì)量節(jié)點(diǎn)對(duì)應(yīng)的元素,可通過靜力凝聚進(jìn)行模型降階.將式(53)展開,有

根據(jù)式(54),角節(jié)點(diǎn)位移向量dC可表示為

將式(56)代入式(55)中,可得僅與中節(jié)點(diǎn)有關(guān)的有限元離散運(yùn)動(dòng)方程

對(duì)于結(jié)構(gòu)時(shí)程分析,這里采用中心差分法進(jìn)行時(shí)間離散.當(dāng)采用等時(shí)間步長 Δt時(shí),中心差分法第n步的速度和加速度分別為

將式(61)代入式(57),可得第n+1 步的位移更新格式

與此同時(shí),將式(60) 代入式(61),可以得到n=-1步的起步位移向量

5 數(shù)值算例

由于本文主要研究不同方法的精度對(duì)比,簡潔起見且不失一般性,算例均采用無量綱單位.

5.1 長方體空腔問題

首先,考慮齊次邊界條件的長方體空腔頻率計(jì)算問題,長方體空腔的長度、寬度和高度分別為Lx=2,Ly=1 和Lz=1.該問題的頻率解析解[28]為

其中,c為波速,本算例中取為1.圖4 為該問題采用H20 和T10 兩種單元的網(wǎng)格劃分情況,其中,H20 單元模型包括216,512 和1000 個(gè)單元(1225,2673 和4961 個(gè)自由度),T10 單元模型包括1080,2560 和5000 個(gè)單元(1981,4401 和8261 個(gè)自由度).為了方便對(duì)比,圖4 有限元模型中,一個(gè)H20 單元對(duì)應(yīng)5 個(gè)T10 單元.

圖4 長方體空腔問題有限元網(wǎng)格Fig.4 Finite element meshes of the rectangular cavity problem

圖5 和表1 給出了采用圖4 系列有限元模型得到的前4 階頻率計(jì)算結(jié)果.從圖5 的頻率收斂結(jié)果可見,H20 和T10 單元的MNLM 集中質(zhì)量矩陣的頻率計(jì)算精度均明顯優(yōu)于HRZ 集中質(zhì)量矩陣.同時(shí),表1 的頻率誤差結(jié)果對(duì)比表明,對(duì)于H20 單元,MNLM集中質(zhì)量矩陣的前4 階頻率計(jì)算誤差與HRZ 之間比值約為0.66,與式(52)給出的理論結(jié)果吻合良好;對(duì)于T10 單元,MNLM 集中質(zhì)量矩陣與HRZ 集中質(zhì)量矩陣相比的頻率誤差比值小于0.40,精度優(yōu)勢(shì)更為顯著.另外,值得指出的是,無論是HRZ 還是MLNM 集中質(zhì)量矩陣,在相同自由度數(shù)條件下,T10單元的頻率計(jì)算精度均優(yōu)于H20 單元.

表1 長方體空腔問題前四階頻率計(jì)算精度對(duì)比Table 1 Accuracy comparison of the first four frequencies for the rectangular cavity problem

圖5 長方體空腔問題前四階頻率收斂特性對(duì)比Fig.5 Convergence comparison of the first four frequencies for the rectangular cavity problem

5.2 彈性力學(xué)長方體問題

第2 個(gè)算例是彈性力學(xué)長方體問題.該模型的幾何尺寸與長方體空腔問題一致,材料彈性拉梅常數(shù) λ=15/26,μ=5/13.邊界條件如圖6 所 示: 在x={0,Lx}處uy=uz=0,在y={0,Ly} 處ux=uz=0,在z={0,Lz} 處ux=uy=0.根據(jù)文獻(xiàn)[29],可導(dǎo)出該模型的壓力波和剪力波頻率解析解分別為

圖6 長方體彈性連續(xù)體模型Fig.6 Description of the rectangular elastic continuum problem

其中,對(duì)于壓力波頻率,下標(biāo)需滿足l,m,n≥1;對(duì)于剪力波頻率,下標(biāo)需滿足至少有兩個(gè)大于等于1.值得需要注意的是,下標(biāo)均大于等于1 時(shí),剪力波頻率為二重根.

在對(duì)該結(jié)構(gòu)進(jìn)行自由振動(dòng)分析時(shí),同樣采用圖4的有限元網(wǎng)格進(jìn)行計(jì)算.前6 階頻率的收斂對(duì)比結(jié)果見圖7.可見,H20 和T10 單元的MNLM 集中質(zhì)量矩陣的頻率計(jì)算精度均明顯優(yōu)于HRZ,驗(yàn)證了所提方法在彈性力學(xué)問題中同樣具有良好的適用性.

圖7 長方體彈性體前六階頻率收斂特性對(duì)比Fig.7 Convergence comparison of the first six frequencies for the rectangular elastic continuum problem

為了深入驗(yàn)證所提集中質(zhì)量矩陣的動(dòng)力性能,進(jìn)一步計(jì)算該彈性力學(xué)問題的時(shí)程響應(yīng).計(jì)算中考慮如下的位移解析解

其中,取t=0 即可得到有限元?jiǎng)恿Ψ治龅某跏嘉灰坪统跏妓俣?

有限元分析中H20 和T10 單元分別采用圖4中的1000 個(gè)單元和5000 個(gè)單元進(jìn)行時(shí)程分析,時(shí)間步長統(tǒng)一取 Δt=0.01.圖8 給出了t=10 時(shí)刻的各方向位移誤差云圖,結(jié)果表明H20 和T10 單元的MNLM 集中質(zhì)量矩陣的位移計(jì)算精度均優(yōu)于對(duì)應(yīng)的HRZ 集中質(zhì)量矩陣.此外,圖9 給出了歸一化計(jì)算時(shí)間[30]對(duì)比結(jié)果.相比HRZ 集中質(zhì)量矩陣,H20單元的MNLM 集中質(zhì)量矩陣的計(jì)算效率提升了約4.3 倍,T10 單元的MNLM 集中質(zhì)量矩陣的計(jì)算效率提升了約3.0 倍.由此可得,MNLM 集中質(zhì)量矩陣不僅可以顯著提高計(jì)算精度,而且在進(jìn)行時(shí)程分析時(shí),與結(jié)構(gòu)動(dòng)力分析降階模型相結(jié)合可以提高計(jì)算效率.

圖8 長方體彈性體 t=10 時(shí)刻的位移誤差云圖Fig.8 Displacement error contour plots for the rectangular elastic continuum problem att=10

圖9 長方體彈性體時(shí)程分析效率對(duì)比Fig.9 Efficiency comparison of the transient analysis for the rectangular elastic continuum problem

5.3 1/4 圓筒空腔問題

最后一個(gè)算例是圖10 所示的1/4 圓筒空腔問題.該模型的內(nèi)徑ri=1,外徑ro=2,高度H=1,波速c=1.當(dāng)采用齊次邊界條件時(shí),該問題的頻率解析解為

圖10 1/4 圓筒空腔問題Fig.10 Description of the quarter-cylinder cavity problem

圖11 1/4 圓筒空腔問題有限元網(wǎng)格Fig.11 Finite element meshes of the quarter-cylinder cavity problem

圖12 1/4 圓筒空腔問題前四階頻率收斂特性對(duì)比Fig.12 Convergence comparison of the first four frequencies for the quarter-cylinder cavity problem

對(duì)于該問題,我們同樣進(jìn)行時(shí)程分析,其中考慮的空腔內(nèi)解析聲壓解為

有限元分析的初始條件可通過在式(70) 中令t=0 給出.同樣,式(11)中源項(xiàng)b也可由式(70)求得

在時(shí)程分析中,有限元模型采用圖11 中的1024 個(gè)H20 單元和5120 個(gè)T10 單元,時(shí)間步長為Δt=0.01.圖13 給出了x=(1.5,π/4,0.5) 處的聲壓數(shù)值解uh的時(shí)程曲線及誤差對(duì)比結(jié)果,圖14 給出了t=12時(shí)刻的位移誤差云圖.結(jié)果顯示,對(duì)于時(shí)程分析,H20 和T10 單元的MNLM 集中質(zhì)量矩陣的計(jì)算精度同樣優(yōu)于對(duì)應(yīng)的HRZ 集中質(zhì)量矩陣.

圖13 1/4 圓筒空腔問題聲壓時(shí)程對(duì)比Fig.13 Comparison of the acoustic pressure time history for the quarter-cylinder cavity problem

圖14 1/4 圓筒空腔問題 t=12 時(shí)刻的聲壓誤差云圖Fig.14 Acoustic pressure error contour plots for the quarter-cylinder cavity problem att=12

與此同時(shí),圖15 給出了MNLM 和HRZ 集中質(zhì)量矩陣的計(jì)算效率對(duì)比.該問題的計(jì)算效率結(jié)果再次表明,對(duì)于具有非均勻網(wǎng)格特征的1/4 圓筒空腔問題,相較于HRZ 集中質(zhì)量矩陣,H20 和T10 單元的MNLM 集中質(zhì)量矩陣的計(jì)算效率分別提升了約1.3 倍和3.0 倍,實(shí)現(xiàn)了計(jì)算效率的顯著提升.

圖15 1/4 圓筒空腔時(shí)程分析效率對(duì)比Fig.15 Efficiency comparison of the transient analysis for the quartercylinder cavity problem

6 結(jié)論

本文針對(duì)結(jié)構(gòu)分析中應(yīng)用廣泛的三維20 節(jié)點(diǎn)六面體單元和10 節(jié)點(diǎn)四面體單元,系統(tǒng)研究了其集中質(zhì)量矩陣構(gòu)造方法和精度.通過引用一種包含待定參數(shù)的三維20 節(jié)點(diǎn)六面體單元廣義集中質(zhì)量矩陣,建立了相應(yīng)的頻率精度表達(dá)式,進(jìn)而揭示了不同集中質(zhì)量矩陣的理論精度.研究表明,對(duì)角元素比例縮放法,即HRZ 方法,作為所提廣義集中質(zhì)量矩陣構(gòu)造方法的一個(gè)特例,雖然能夠保證集中質(zhì)量矩陣元素非負(fù)性,但其精度并非最優(yōu).隨后,通過對(duì)廣義集中質(zhì)量矩陣頻率精度進(jìn)行優(yōu)化,提出了一種精度更優(yōu)的三維20 節(jié)點(diǎn)六面體單元中節(jié)點(diǎn)集中質(zhì)量矩陣構(gòu)造方法,并將其推廣到10 節(jié)點(diǎn)四面體單元.中節(jié)點(diǎn)集中質(zhì)量矩陣同樣具有非負(fù)特性,但角節(jié)點(diǎn)的質(zhì)量為零,可以方便地進(jìn)行模型降階.典型算例的頻率和時(shí)程計(jì)算結(jié)果均表明,與HRZ 集中質(zhì)量矩陣相比,三維20 節(jié)點(diǎn)六面體單元和10 節(jié)點(diǎn)四面體單元的中節(jié)點(diǎn)集中質(zhì)量矩陣在計(jì)算精度和效率方面都得到了顯著提升.以20 節(jié)點(diǎn)六面體單元為例,就精度而言,如表1 所示,中節(jié)點(diǎn)集中質(zhì)量矩陣的頻率計(jì)算誤差比HRZ 集中質(zhì)量矩陣降低了約1/3;就效率而言,如圖9 和圖15 所示,對(duì)于彈性力學(xué)問題和勢(shì)問題,中節(jié)點(diǎn)集中質(zhì)量矩陣的時(shí)程分析計(jì)算效率約為HRZ 集中質(zhì)量矩陣的5 倍和2 倍.