楊石扣，艾華東

(江西理工大學(xué) 土木與測(cè)繪工程學(xué)院，江西 贛州 341000)

1 研究背景

數(shù)值流形法(Numerical Manifold Method,NMM)是Shi[1]在關(guān)鍵塊體理論(Key Block Theory,KBT)和非連續(xù)變形分析(Discontinuous Deformation Analysis,DDA)的基礎(chǔ)上提出的一種具有雙重覆蓋系統(tǒng)的數(shù)值計(jì)算方法，可以統(tǒng)一處理連續(xù)和非連續(xù)變形問題。有限元法處理位移邊界條件較為方便，然而對(duì)于數(shù)值流形法[2]、廣義有限元法[3]和無網(wǎng)格法[4]等數(shù)值分析方法而言，其處理方法往往存在各自缺陷[5]。目前，數(shù)值流形法和無網(wǎng)格法處理位移邊界條件常采用罰函數(shù)法和Lagrange乘子法[4,6-7]。采用罰函數(shù)法得到的整體剛度矩陣具有對(duì)稱、正定且實(shí)施方便等特點(diǎn)，且數(shù)學(xué)網(wǎng)格一般與邊界不匹配，故而在數(shù)值流形法中應(yīng)用較為廣泛[5]，但罰函數(shù)法選擇合適的罰值困難，取值過大則影響方程組工作性態(tài)，取值過小則計(jì)算誤差較大。Zhang等[8]提出采用位移約束方程法施加本質(zhì)邊界條件，并應(yīng)用于無網(wǎng)格法，Wu等[9]則提出采用BFCM(Boundary Flux Collocation Method)處理無網(wǎng)格法中的本質(zhì)邊界問題。蔡永昌等[4,10]在無網(wǎng)格法高斯數(shù)值積分點(diǎn)上對(duì)剛度矩陣和荷載向量進(jìn)行修正，提出了一種準(zhǔn)確施加位移邊界和材料不連續(xù)條件的方法，具有實(shí)施簡(jiǎn)單、穩(wěn)定和求解精度高等優(yōu)點(diǎn)。Mihara等[11]將罰函數(shù)法應(yīng)用于混合虛功原理，提出了一種混合型罰函數(shù)法(Hybrid-Type Penalty Method,HPM)，可方便地進(jìn)行非線性分析。Kourepinis等[12]針對(duì)高階數(shù)值流形法的邊界協(xié)調(diào)問題，使用投影矩陣提出了一種邊界處理方法，可方便地處理DDA、NMM和XFEM中的邊界問題。Zheng等[13]和李偉等[14]提出了一種精確的邊界施加方法，可方便地處理本質(zhì)邊界、材料接觸邊界和滲流邊界。蘇海東等[5]分析了獨(dú)立覆蓋流形法中本質(zhì)邊界條件不易施加的原因，并提出了改進(jìn)邊界覆蓋函數(shù)和直接設(shè)定獨(dú)立覆蓋函數(shù)2種方法，可嚴(yán)格滿足邊界條件。對(duì)于地震響應(yīng)和固體振動(dòng)分析，Yang等[15]、Wei等[16]和Wang等[17]提出了相應(yīng)的邊界設(shè)置方法和控制方程。

目前大多數(shù)關(guān)于三維數(shù)值流形法的研究[18-25]，主要采用罰函數(shù)法施加位移邊界條件，但對(duì)于沿著某一方向施加邊界位移的情況未作具體討論，且對(duì)于多步加載，計(jì)算誤差累積效應(yīng)將會(huì)更加明顯。本文通過修正三維數(shù)值流形法控制方程中的位移邊界部分，推導(dǎo)出含沿著某一方向施加邊界位移條件和相應(yīng)的分步加載情況的計(jì)算公式，以期擴(kuò)大控制方程在位移邊界處理上的適用范圍和減少多步加載誤差累積效應(yīng)。通過2個(gè)典型算例的對(duì)比分析，驗(yàn)證了該方法的準(zhǔn)確性。

2 三維數(shù)值流形法的基本概念

圖1 數(shù)值流形法的覆蓋系統(tǒng)和流形單元Fig.1 Covers and manifold elements for NumericalManifold Method

三維數(shù)值流形法在數(shù)學(xué)覆蓋Mi上定義權(quán)函數(shù)φi(x)，權(quán)函數(shù)φi(x)滿足式(1)所示的條件。

(1)

(2)

傳遞系數(shù)相應(yīng)取值為

(3)

(4)

CT(x)=

(5)

式中n為多項(xiàng)式次數(shù)。

總位移函數(shù)u(x)可以寫為

(7)

在實(shí)際應(yīng)用中，數(shù)值計(jì)算均基于單元，故改寫式(7)可得單元e內(nèi)部任一點(diǎn)x的位移函數(shù)ue(x)為

(9)

(10)

其中di的形式為

由于本文采用六面體數(shù)學(xué)網(wǎng)格，故矩陣N取有限元六面體單元形函數(shù)，并采用高斯積分[25]。

3 位移邊界處理方法

(12)

式中Ω是以邊界Γt和Γu所形成的物理域。

(13)

將式(8)代入式(13)并考慮虛覆蓋系數(shù)δdT的任意性，可得

Kd=R。

(14)

式中K和R分別為剛度矩陣和荷載向量，其形式分別為

(15)

(16)

式中:B為應(yīng)變轉(zhuǎn)換矩陣；D為彈性矩陣。

考慮到計(jì)算中加載步對(duì)位移影響的修正，由第k加載步至第k+1加載步的位移增量(見圖2)為

圖2 沿某方向tn施加邊界位移示意圖Fig.2 Sketch of boundary displacement applied alonga direction tn

(17)

其中，

(18)

(19)

將式(19)代入式(16)可得

(20)

4 數(shù)值算例和結(jié)果

根據(jù)第3節(jié)的方法和

公式，編寫了相應(yīng)的C++程序。為了驗(yàn)證程序和文中所述方法的可行性和正確性，選取2個(gè)典型算例進(jìn)行討論。

4.1 矩形薄板受分布力荷載作用

算例1為一個(gè)寬度L=10 m、高H=20 m和厚W=1 m的薄板[26]，如圖3所示，左邊及下邊受連桿支承，右邊及上邊分別受均布?jí)毫1=20 kN/m2及q2=10 kN/m2作用，不計(jì)體力。其位移場(chǎng)u和v的解析解[26]分別為：

圖3 矩形薄板受分布力荷載作用Fig.3 Rectangular thin plate subjected todistributed load

(21)

計(jì)算時(shí)取薄板彈性模量E=300 MPa，泊松比μ=0.20。采用10×1×20網(wǎng)格數(shù)，則初始物理覆蓋數(shù)為462，流形單元數(shù)為200，如圖4所示。在x=0.001 m、沿z方向和z=10 m、沿x方向分別布置監(jiān)測(cè)點(diǎn)，計(jì)算荷載q1和q2一步加載完成，且取罰函數(shù)s=1 000E，并將計(jì)算結(jié)果和解析解繪于圖5。比較后可知，文中方法的計(jì)算結(jié)果與解析解相吻合，誤差較小，均控制在0.1%以內(nèi)。

圖4 數(shù)值模型網(wǎng)格劃分Fig.4 Meshing of numerical model

圖5 沿z和x方向的水平位移和豎直位移Fig.5 Horizontal and vertical displacements alongz direction and x direction

為了說明文中方法的適應(yīng)性，將計(jì)算模型繞x、y和z軸分別旋轉(zhuǎn)α=30°、β=45°和γ=30°，新模型如圖4(b)所示?？紤]兩種罰函數(shù)s=1 000E和s=10E，在式(20)中按照考慮位移邊界誤差修正(M)和不考慮位移邊界誤差修正(NM)位移邊界誤差修正2種情況，將計(jì)算荷載q1和q2加載步分為1、2、5和10等情況分別進(jìn)行計(jì)算。在x=0.001 m、沿z方向和z=10 m、沿x方向分別布置監(jiān)測(cè)點(diǎn)，并將計(jì)算結(jié)果繪于圖6。比較后可知，罰函數(shù)較大的計(jì)算結(jié)果，位移誤差較??；考慮位移邊界誤差修正比不考慮位移邊界誤差修正的計(jì)算精度高；隨著分步加載數(shù)目的增多，累積計(jì)算誤差逐漸增大，對(duì)不考慮位移邊界誤差修正的誤差累積效應(yīng)更加明顯，對(duì)考慮位移邊界誤差修正的影響不大。

圖6 不同加載步數(shù)、不同罰函數(shù)和是否考慮位移邊界誤差修正時(shí)的位移誤差Fig.6 Displacement errors in cases of different loading steps and penalty functions and whether the modification of displacement boundary error is considered

4.2 矩形薄板受位移荷載作用

算例2為一個(gè)寬度a=5 m、高度b=5 m和厚度W=1 m的薄板，如圖7所示，左右兩邊及下邊均被固定，上邊的位移給定為u=0，v=-η(1-x2/a2)，不計(jì)體力。當(dāng)a=b時(shí)，其位移場(chǎng)v的解析解[26]為

圖7 矩形薄板受位移荷載作用Fig.7 Rectangular thin plate subjected to displacement load

(22)

計(jì)算時(shí)取薄板彈性模量E= 300 MPa，泊松比μ=0.20，罰函數(shù)s=1 000E。采用20×1×10(見圖8)、40×1×20和80×1×40共3種網(wǎng)格數(shù)分別進(jìn)行計(jì)算。

圖8 數(shù)值模型網(wǎng)格劃分Fig.8 Meshing of numerical model

監(jiān)測(cè)點(diǎn)坐標(biāo)為(-2.5，0，2.5)，將位移v的計(jì)算結(jié)果繪于圖9。比較后可知，文中方法的計(jì)算結(jié)果與解析解相吻合，誤差控制在1.21%以內(nèi)；隨著計(jì)算網(wǎng)格的加密，計(jì)算誤差逐漸減小，計(jì)算精度逐漸提高。

圖9 監(jiān)測(cè)點(diǎn)在不同網(wǎng)格密度下的位移誤差Fig.9 Displacement errors of a monitoring point under different grid densities

5 結(jié) 論

(1)本文提出了一種基于三維數(shù)值流形法的位移邊界處理方法，并編寫了相應(yīng)的C++程序。此法通過修正傳統(tǒng)三維數(shù)值流形法控制方程中的位移邊界部分，推導(dǎo)出含沿著某一方向施加邊界位移條件和相應(yīng)的分步加載情況的計(jì)算公式，擴(kuò)大了控制方程在位移邊界處理上的適用范圍和減少多步加載誤差累積效應(yīng)。

(2)2個(gè)算例的計(jì)算結(jié)果表明：文中方法的計(jì)算結(jié)果與解析解相吻合，修正公式適用于不同方向施加邊界位移的情況，具有較強(qiáng)的適應(yīng)性；罰函數(shù)較大的計(jì)算結(jié)果，位移誤差較?。豢紤]位移邊界誤差修正比不考慮位移邊界誤差修正的計(jì)算精度高；隨著分步加載數(shù)目的增多，累積計(jì)算誤差逐漸增大，對(duì)不考慮位移邊界誤差修正的誤差累積效應(yīng)更加明顯，對(duì)考慮位移邊界誤差修正的影響不大；隨著計(jì)算網(wǎng)格的加密，計(jì)算誤差逐漸減小，計(jì)算精度逐漸提高。