三維20 節點六面體和10 節點四面體單元的高精度中節點集中質量矩陣1)

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

(廈門大學土木工程系,福建省濱海土木工程數字仿真重點實驗室,福建廈門 361005)

引言

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

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

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

1 結構動力分析的有限元離散方程

不失一般性,考慮如下的經典波動方程等效積分弱形式[1]

其中,δ 表示變分符號,上標T 代表轉置.對于膜或聲場問題,場變量u退化為標量u,表示膜的橫向位移或聲壓,外力或源項b退化為標量b.與之對應的梯度矩陣和材料矩陣D分別為

其中,c為波速.對于彈性力學問題,矢量場u={ux,uy,uz}T表示彈性體內一點的位移,b={bx,by,bz}T表示體力,相應的梯度矩陣和材料彈性矩陣D分別為

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

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

其中,NI(ξ) 表示單元第I個節點的形函數,nen為單元節點個數.將式(6)代入式(1)中,可得結構動力分析的有限元離散方程

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

其中,1 是單位矩陣,對于聲場問題,1=1;對于彈性連續體問題,1=diag{1,1,1}.

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

其中,ωh為數值頻率,? 為與之對應的振動模態.將式(12)代入式(7),可得自由振動分析的有限元離散方程

2 三維20 節點六面體和10 節點四面體單元的高精度集中質量矩陣構造方法

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

2.1 H20 和T10 單元的HRZ 集中質量矩陣

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

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

圖1 H20 單元3 種集中質量矩陣的節點質量分布示意圖Fig.1 Schematic illustration of the nodal mass distribution for three mass lumping schemes of H20 element

圖2 T10 單元3 種集中質量矩陣的節點質量分布示意圖Fig.2 Schematic illustration of the nodal mass distribution for three mass lumping schemes of T10 element

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

其中,r為縮放系數

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

圖1 和圖2 也給出了這兩種單元的HRZ 集中質量矩陣節點質量分布情況.

2.2 H20 和T10 單元的非負中節點集中質量矩陣

如前所述,H20 和T10 單元的行求和集中質量矩陣,負質量都出現在角節點上,而中節點則沒有負質量問題.因此,可以基于行求和集中質量矩陣中的非負質量元素,通過單元質量守恒將其進行比例放縮,同時將行求和集中質量矩陣中的負對角元素置零,構造一種新型非負集中質量矩陣.該集中質量矩陣的構造流程為

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

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

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

本節以H20 單元為例,研究集中質量矩陣構造形式對頻率計算精度的影響.為了更全面衡量各類集中質量矩陣的精度,構造如下H20 單元廣義集中質量矩陣

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

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

圖3 H20 單元典型節點影響域及波動方向示意圖Fig.3 Illustration of the typical nodal influence domains and wave propagation angles for a mesh with H20 elements

假設離散節點的波動形式為諧波,4 個典型節點系數可以表示為

方便表述起見,引入算子cos(±lkx±mky±nkz)

當式(28)僅包括兩個下標時,有

對于4 個典型節點,將式(26)和式(27)代入式(7),并利用簡諧波假定化簡可得

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

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

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

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

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

進一步在式(45)中引入余弦項關于kh的泰勒展開,可得H20 單元廣義集中質量矩陣的頻率誤差表達式

其中h為特征單元長度

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

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

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

4 基于中節點集中質量矩陣的時程分析

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

將H20 或T10 單元的MNLM 集中質量矩陣代入式(7)的有限元離散運動方程,同時將零質量角節點和非零質量中節點進行分類,可得

其中,下標C和M分別表示單元的角節點和中節點,MMM,KCC和KMM均為正定矩陣.對式(53)中的零質量節點對應的元素,可通過靜力凝聚進行模型降階.將式(53)展開,有

根據式(54),角節點位移向量dC可表示為

將式(56)代入式(55)中,可得僅與中節點有關的有限元離散運動方程

對于結構時程分析,這里采用中心差分法進行時間離散.當采用等時間步長 Δt時,中心差分法第n步的速度和加速度分別為

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

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

5 數值算例

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

5.1 長方體空腔問題

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

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

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

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

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

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

5.2 彈性力學長方體問題

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

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

其中,對于壓力波頻率,下標需滿足l,m,n≥1;對于剪力波頻率,下標需滿足至少有兩個大于等于1.值得需要注意的是,下標均大于等于1 時,剪力波頻率為二重根.

在對該結構進行自由振動分析時,同樣采用圖4的有限元網格進行計算.前6 階頻率的收斂對比結果見圖7.可見,H20 和T10 單元的MNLM 集中質量矩陣的頻率計算精度均明顯優于HRZ,驗證了所提方法在彈性力學問題中同樣具有良好的適用性.

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

為了深入驗證所提集中質量矩陣的動力性能,進一步計算該彈性力學問題的時程響應.計算中考慮如下的位移解析解

其中,取t=0 即可得到有限元動力分析的初始位移和初始速度.

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

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

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

5.3 1/4 圓筒空腔問題

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

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

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

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

對于該問題,我們同樣進行時程分析,其中考慮的空腔內解析聲壓解為

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

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

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

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

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

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

6 結論

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