B樣條函數(shù)法似計(jì)算方法

張友華,袁 波*,徐子杰,鄭世宇

(1.貴州大學(xué) 空間結(jié)構(gòu)研究中心,貴州 貴陽 550025;2.貴州大學(xué) 貴州省結(jié)構(gòu)工程重點(diǎn)實(shí)驗(yàn)室,貴州 貴陽 550025)

B樣條函數(shù)法是基于變分原理和B樣條為基礎(chǔ)建立的數(shù)值計(jì)算方法,作為一種良好的逼近工具,被應(yīng)用于分析工程中的各類結(jié)構(gòu)。樣條函數(shù)方法不僅擁有通用有限元法的許多優(yōu)點(diǎn),包含計(jì)算精度高、系數(shù)矩陣具有對稱性、適應(yīng)不同邊界條件等;基于B樣條的優(yōu)異性能,該方法還具有計(jì)算量更小、更易于計(jì)算機(jī)程序編制等優(yōu)點(diǎn)。隨著科學(xué)技術(shù)的進(jìn)步,工程對于數(shù)值計(jì)算方法要求更精確,實(shí)際情況中結(jié)構(gòu)承受的荷載并非線性分布,因此,基于樣條函數(shù)法研究結(jié)構(gòu)受任意分布荷載的數(shù)值計(jì)算方法具有現(xiàn)實(shí)意義。

1946年,美國數(shù)學(xué)家I.J.Schoenoerg首先提出B樣條函數(shù)。1976年,國內(nèi)外學(xué)者開始將B樣條函數(shù)用于結(jié)構(gòu)分析,并提出了多種樣條函數(shù)方法。石鐘慈[1]基于三次B樣條及變分法推導(dǎo)出樣條有限元法,分析各種邊界條件的板梁組合彎曲問題。秦榮[2]基于B樣條、三角函數(shù)及能量法提出樣條有限點(diǎn)法,分析彈性薄板受均布荷載彎曲問題。此后,秦榮[3]在樣條函數(shù)的基礎(chǔ)上提出樣條邊界元法、樣條子域法、板殼分析等新理論新方法,并用于規(guī)則結(jié)構(gòu)及非規(guī)則結(jié)構(gòu)的分析求解。秦榮在文獻(xiàn)[4]中闡述樣條函數(shù)法中結(jié)構(gòu)受集中力、均布荷載、線荷載、任意荷載的數(shù)值積分方法。此后,多位學(xué)者在文獻(xiàn)[2]中對樣條函數(shù)法及任意荷載的計(jì)算方法進(jìn)行了大量的研究。

此前,基于最小勢能原理,運(yùn)用B樣條函數(shù)的分部積分公式計(jì)算任意荷載函數(shù)與B樣條位移函數(shù)乘積的積分,該方法存在計(jì)算高次B樣條函數(shù)構(gòu)成的基函數(shù)運(yùn)算量大、處理復(fù)雜荷載不便、編程工作量大等問題。在本文中,提出使用n次Lagrange插值多項(xiàng)式構(gòu)造需求的荷載函數(shù),發(fā)揮多項(xiàng)式與樣條函數(shù)的乘積便于數(shù)值積分的特性,簡化荷載列陣的計(jì)算,本方法的具體步驟和數(shù)值積分公式均在算例中體現(xiàn)。參考文獻(xiàn)[5-10]中的算例,基于樣條函數(shù)法解圓柱殼問題時(shí),由于坐標(biāo)軸建立的特殊性,結(jié)構(gòu)自重荷載與環(huán)向坐標(biāo)位置為三角函數(shù)關(guān)系,且可同文獻(xiàn)中的計(jì)算結(jié)果進(jìn)行直接比較,驗(yàn)證本方法的誤差大小,故取文獻(xiàn)[6]相同尺寸的圓柱殼算例。

1 樣條函數(shù)方法解圓柱殼結(jié)構(gòu)

1.1 位移函數(shù)

圖1為圓柱殼結(jié)構(gòu)模型,坐標(biāo)軸x和y布置在板的中性平面內(nèi),坐標(biāo)原點(diǎn)布置在圓柱殼的頂端角點(diǎn)上,取x軸沿殼體中曲面的環(huán)向,y軸沿殼體中曲面的縱向,z軸沿殼體中曲面的法向。結(jié)構(gòu)在x與y方向上的長度分別取為a和b,則取值范圍可表示為0≤x≤a及0≤y≤b。沿x方向及y方向分別作均勻劃分為N段、M段,本文選擇以三次B樣條函數(shù)乘積的線性組合作為結(jié)構(gòu)的位移函數(shù)。表示如下:

(1)

圖1 圓柱殼結(jié)構(gòu)模型Fig.1 Structural model of cylindrical shell

現(xiàn)將(1)式寫成矩陣形式:

(2)

(3)

上式(2)、(3)中,ψ?φ表示矩陣ψ與φ的Kronecker乘積,aij、bij和cij為樣條線點(diǎn)參數(shù),φi(x)和ψj(y)都是與三次樣條函數(shù)有關(guān)的基函數(shù)。

其中,

(4)

1.2 樣條基函數(shù)

本文選用如下φi(x)基函數(shù),(ψj(y) 基函數(shù)也選用如下相同格式,只需改hx為hy,改N為M),如下所示,式(2)的φ可以寫成

φ=ΦQ

(5)

其中,

(6)

(7)

式中,I為(N-3)行(N-3)列的單位矩陣,樣條基函數(shù)組合系數(shù)矩陣Φ=φ3(x/h-i)為三次樣條基函數(shù),Q為(N+3)行(N+3)列的樣條基函數(shù)組合系數(shù)矩陣。矩陣Ψ與矩陣φ的構(gòu)造方式相同,且選用相同樣條基函數(shù)組合系數(shù)矩陣Q(只需把hx改為hy)。

1.3 薄板的基本方程

論文假設(shè)薄板中平面應(yīng)力與彎曲應(yīng)力互不影響,薄板的基本方程則由平面問題的基本方程和彎曲問題的基本方程組合而成。

1.3.1幾何方程

位移函數(shù)采用(1)式所示的雙三次樣條函數(shù)的形式,則幾何方程為

(8)

其中,

(9)

(10)

(11)

1.3.2物理方程

(12)

式中,D和R都是彈性矩陣。對于正交異性板殼及各項(xiàng)同性板殼,D和R為

(13)

其中,

(14)

(15)

當(dāng)板殼各項(xiàng)同性時(shí),則式(14)(15)簡化為:

(16)

(17)

式中:Ex、Ey、Rx、Ry及μx、μy為x,y方向的彈性模量及泊松系數(shù),d為板殼厚度;Nx、Ny、Nxy為殼體的薄膜內(nèi)力;Mx、My、Mxy分別為殼體的彎矩和扭矩。

1.4 圓柱薄板的總勢能泛函

(18)

式中:Π為圓柱殼拱的能量泛函;U為位移向量,q為荷載向量。

將式(8)、(12)帶入圓形薄板總勢能泛函方程(18)中,可以將其改寫成如下形式:

(19)

運(yùn)用變分原理可以得到

Gδ=f

(20)

其中,

(21)

(22)

1.5 板塊的剛度方程

式(20)為圓柱殼的總剛度方程?，F(xiàn)將式(21)、(22)展開,得

(23)

(24)

(25)

本文對以上公式符號做出了如下定義

(26)

2 荷載列陣的數(shù)值積分

2.1 基于Lagrange插值的近似計(jì)算

由圖1中可知,圓形薄板任意位置處的自重荷載q均可以沿圓形薄板分解為切向和法相的荷載,則x,y,z方向的荷載可用函數(shù)表示,即

(27)

將上式(27)帶入式(24)中,得

(28)

在此處設(shè)Lagrange插值逼近的函數(shù)如下式:

(29)

式中:Pn(x) 表示正弦函數(shù)的n次Lagrange插值多項(xiàng)式;Ln(x)表示余弦函數(shù)的n次Lagrange插值多項(xiàng)式。

將式(29)帶入式(28)中,得

(30)

依據(jù)n次Lagrange中插值多項(xiàng)式函數(shù)的特點(diǎn),式(30)荷載列陣在x方向上的積分項(xiàng)可展開為

(31)

式中:an和bn分別為插值多項(xiàng)式Pn(x)及Ln(x)中冪函數(shù)xn的系數(shù)。

因此,只要求出每一項(xiàng)進(jìn)行累加即可得到目標(biāo)函數(shù)在區(qū)間上的定積分。

由式(5)及(31)可以得到下式:

(32)

2.2 數(shù)值積分公式

上式中矩陣Φ中的每一項(xiàng)是三次B樣條函數(shù),由B樣條函數(shù)的數(shù)值計(jì)算特性可知,單獨(dú)對B樣條函數(shù)進(jìn)行多重積分與多次求導(dǎo)是易于得出結(jié)果的,因此,可運(yùn)用分部積分法建立下列積分公式對式(32)中每一項(xiàng)進(jìn)行計(jì)算得到荷載列陣的值(依據(jù)本文需要,取到四次Lagrange中插值多項(xiàng)式函數(shù)),其中,每一項(xiàng)的數(shù)值積分方法如下所示。

在區(qū)間[0,a]上N等分,令t=(x/hx-i),通過換元法積分,則可得到下式:

(33)

(34)

再運(yùn)用分部積分法對式(34)中每一項(xiàng)進(jìn)行數(shù)值計(jì)算得到下式:

(35)

基于B樣條函數(shù)φn(x)的求積分公式便可計(jì)算式(35),積分公式如下所示:

(36)

式中:B樣條函數(shù)φn(x)右上角負(fù)號表示求積分次數(shù)。

2.3 構(gòu)造插值多項(xiàng)式

在考慮逼近曲線的光滑性、連續(xù)性、易于編程等特點(diǎn)后,本文選用二次、三次及四次等不同Lagrange插值多項(xiàng)式逼近正弦及余弦函數(shù)。由于結(jié)構(gòu)從坐標(biāo)原點(diǎn)處開始,沿x正負(fù)方向荷載變化規(guī)律一致,其逼近區(qū)間選為[0,π/2],其插值節(jié)點(diǎn)與插值多項(xiàng)式如下表1、表2所示,插值多項(xiàng)式函數(shù)圖如圖2所示。

表1 正弦函數(shù)插值節(jié)點(diǎn)與插值多項(xiàng)式Tab.1 Sinusoidal function interpolation nodes and interpolation polynomials

表2 余弦函數(shù)插值節(jié)點(diǎn)與插值多項(xiàng)式Tab.2 Cosine function interpolation nodes and interpolation polynomials

圖2 插值多項(xiàng)式函數(shù)圖Fig.2 Interpolation polynomial approximation rendering

正弦函數(shù)插值多項(xiàng)式的具體表達(dá)式如下式所示:

(37)

余弦函數(shù)插值多項(xiàng)式的具體表達(dá)式如下式所示:

(38)

3 數(shù)值算例

圓柱殼拱屋頂承受自重荷載。圓柱殼直線邊界自由。參考文獻(xiàn)中的算例,為同其他計(jì)算方法作比較,本文中算例曲線邊界約束設(shè)置為簡支,其結(jié)構(gòu)及尺寸如圖3所示。依據(jù)結(jié)構(gòu)對稱的受力特性,取對半結(jié)構(gòu)進(jìn)行計(jì)算,網(wǎng)格劃分為10×10的網(wǎng)格。本方法計(jì)算結(jié)果與相關(guān)文獻(xiàn)方法及精確解(扁殼解法)結(jié)果比較見表3和圖4,與精確解(扁殼解法)的相對誤差見表4。

表3 特殊點(diǎn)的位移和內(nèi)力Tab.3 Displacement and internal force of special points

表4 特殊點(diǎn)位移和內(nèi)力的相對誤差Tab.4 Relative error of displacement and internal force at special points

圖3 結(jié)構(gòu)尺寸Fig.3 Structural dimension diagram

圖4 內(nèi)力分布圖Fig.4 Internal force distribution diagram

綜合分析上述算例結(jié)果,如表4所示,二次、三次、四次插值方法在特殊點(diǎn)處的相對誤差均低于1.5%,均滿足工程要求的計(jì)算精度,且三次、四次插值方法的誤差低于二次插值方法,結(jié)合圖2中三次、四次插值函數(shù)逼近效果優(yōu)于二次插值函數(shù)分析,表明兩者之間存在著一致性。如圖4所示,三種插值方法的位移與內(nèi)力沿不同方向上的分布均與精確解(扁殼解法)一致。

4 結(jié)論

本文發(fā)展了一種插值逼近形式的數(shù)值計(jì)算方法求解結(jié)構(gòu)受任意荷載問題?；跇訔l函數(shù)法及Lagrange插值法,推導(dǎo)了近似計(jì)算方法的數(shù)值積分公式,并通過解圓柱殼結(jié)構(gòu)算例詳細(xì)展示了本方法的具體步驟。通過文中算例分析,得出以下結(jié)論:

1)本文使用三種不同插值多項(xiàng)式得到解與精確解(扁殼解法)及相關(guān)文獻(xiàn)解對比均具有較高計(jì)算精度,證明了本方法能夠有效分析結(jié)構(gòu)受任意荷載問題。

2)近似計(jì)算方法以n次Lagrange插值多項(xiàng)式作為荷載函數(shù),具有與原函數(shù)相同的光滑性、連續(xù)性、單調(diào)性等數(shù)學(xué)特性。本文基于Lagrange插值法的靈活性,取不同插值節(jié)點(diǎn)構(gòu)造了三種不同的插值多項(xiàng)式,基于其具有的廣泛適用性,對兩種不同分布荷載均構(gòu)造了插值函數(shù),說明本方法對于任意荷載均具有較強(qiáng)適用性。

3)本文通過插值逼近的方法及多項(xiàng)式與B樣條函數(shù)乘積易于數(shù)值積分的思路,用近似函數(shù)替代原荷載函數(shù),從而優(yōu)化計(jì)算流程,解決復(fù)雜荷載不便于積分運(yùn)算的問題。為未來研究結(jié)構(gòu)受任意荷載問題提供新思路,具有廣闊的發(fā)展和應(yīng)用前景。