楊石扣，艾華東

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

1 研究背景

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

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

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

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

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

(1)

(2)

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

(3)

(4)

CT(x)=

(5)

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

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

(7)

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

(9)

(10)

其中di的形式為

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

3 位移邊界處理方法

(12)

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

(13)

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

Kd=R。

(14)

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

(15)

(16)

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

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

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

(17)

其中，

(18)

(19)

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

(20)

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

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

公式，編寫了相應的C++程序。為了驗證程序和文中所述方法的可行性和正確性，選取2個典型算例進行討論。

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

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

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

(21)

計算時取薄板彈性模量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)測點，計算荷載q1和q2一步加載完成，且取罰函數(shù)s=1 000E，并將計算結果和解析解繪于圖5。比較后可知，文中方法的計算結果與解析解相吻合，誤差較小，均控制在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

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

圖6 不同加載步數(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為一個寬度a=5 m、高度b=5 m和厚度W=1 m的薄板，如圖7所示，左右兩邊及下邊均被固定，上邊的位移給定為u=0，v=-η(1-x2/a2)，不計體力。當a=b時，其位移場v的解析解[26]為

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

(22)

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

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

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

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

5 結 論

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

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