張友華,袁 波*,徐子杰,鄭世宇
(1.貴州大學 空間結構研究中心,貴州 貴陽 550025;2.貴州大學 貴州省結構工程重點實驗室,貴州 貴陽 550025)
B樣條函數法是基于變分原理和B樣條為基礎建立的數值計算方法,作為一種良好的逼近工具,被應用于分析工程中的各類結構。樣條函數方法不僅擁有通用有限元法的許多優點,包含計算精度高、系數矩陣具有對稱性、適應不同邊界條件等;基于B樣條的優異性能,該方法還具有計算量更小、更易于計算機程序編制等優點。隨著科學技術的進步,工程對于數值計算方法要求更精確,實際情況中結構承受的荷載并非線性分布,因此,基于樣條函數法研究結構受任意分布荷載的數值計算方法具有現實意義。
1946年,美國數學家I.J.Schoenoerg首先提出B樣條函數。1976年,國內外學者開始將B樣條函數用于結構分析,并提出了多種樣條函數方法。石鐘慈[1]基于三次B樣條及變分法推導出樣條有限元法,分析各種邊界條件的板梁組合彎曲問題。秦榮[2]基于B樣條、三角函數及能量法提出樣條有限點法,分析彈性薄板受均布荷載彎曲問題。此后,秦榮[3]在樣條函數的基礎上提出樣條邊界元法、樣條子域法、板殼分析等新理論新方法,并用于規則結構及非規則結構的分析求解。秦榮在文獻[4]中闡述樣條函數法中結構受集中力、均布荷載、線荷載、任意荷載的數值積分方法。此后,多位學者在文獻[2]中對樣條函數法及任意荷載的計算方法進行了大量的研究。
此前,基于最小勢能原理,運用B樣條函數的分部積分公式計算任意荷載函數與B樣條位移函數乘積的積分,該方法存在計算高次B樣條函數構成的基函數運算量大、處理復雜荷載不便、編程工作量大等問題。在本文中,提出使用n次Lagrange插值多項式構造需求的荷載函數,發揮多項式與樣條函數的乘積便于數值積分的特性,簡化荷載列陣的計算,本方法的具體步驟和數值積分公式均在算例中體現。參考文獻[5-10]中的算例,基于樣條函數法解圓柱殼問題時,由于坐標軸建立的特殊性,結構自重荷載與環向坐標位置為三角函數關系,且可同文獻中的計算結果進行直接比較,驗證本方法的誤差大小,故取文獻[6]相同尺寸的圓柱殼算例。
圖1為圓柱殼結構模型,坐標軸x和y布置在板的中性平面內,坐標原點布置在圓柱殼的頂端角點上,取x軸沿殼體中曲面的環向,y軸沿殼體中曲面的縱向,z軸沿殼體中曲面的法向。結構在x與y方向上的長度分別取為a和b,則取值范圍可表示為0≤x≤a及0≤y≤b。沿x方向及y方向分別作均勻劃分為N段、M段,本文選擇以三次B樣條函數乘積的線性組合作為結構的位移函數。表示如下:
(1)

圖1 圓柱殼結構模型Fig.1 Structural model of cylindrical shell
現將(1)式寫成矩陣形式:
(2)
(3)
上式(2)、(3)中,ψ?φ表示矩陣ψ與φ的Kronecker乘積,aij、bij和cij為樣條線點參數,φi(x)和ψj(y)都是與三次樣條函數有關的基函數。
其中,

(4)
本文選用如下φi(x)基函數,(ψj(y) 基函數也選用如下相同格式,只需改hx為hy,改N為M),如下所示,式(2)的φ可以寫成
φ=ΦQ
(5)
其中,
(6)
(7)
式中,I為(N-3)行(N-3)列的單位矩陣,樣條基函數組合系數矩陣Φ=φ3(x/h-i)為三次樣條基函數,Q為(N+3)行(N+3)列的樣條基函數組合系數矩陣。矩陣Ψ與矩陣φ的構造方式相同,且選用相同樣條基函數組合系數矩陣Q(只需把hx改為hy)。
論文假設薄板中平面應力與彎曲應力互不影響,薄板的基本方程則由平面問題的基本方程和彎曲問題的基本方程組合而成。
1.3.1幾何方程
位移函數采用(1)式所示的雙三次樣條函數的形式,則幾何方程為
(8)
其中,
(9)
(10)
(11)
1.3.2物理方程
(12)
式中,D和R都是彈性矩陣。對于正交異性板殼及各項同性板殼,D和R為
(13)
其中,
(14)
(15)
當板殼各項同性時,則式(14)(15)簡化為:
(16)
(17)
式中:Ex、Ey、Rx、Ry及μx、μy為x,y方向的彈性模量及泊松系數,d為板殼厚度;Nx、Ny、Nxy為殼體的薄膜內力;Mx、My、Mxy分別為殼體的彎矩和扭矩。
(18)
式中:Π為圓柱殼拱的能量泛函;U為位移向量,q為荷載向量。
將式(8)、(12)帶入圓形薄板總勢能泛函方程(18)中,可以將其改寫成如下形式:
(19)
運用變分原理可以得到
Gδ=f
(20)
其中,
(21)
(22)
式(20)為圓柱殼的總剛度方程。現將式(21)、(22)展開,得
(23)
(24)
(25)
本文對以上公式符號做出了如下定義
(26)
由圖1中可知,圓形薄板任意位置處的自重荷載q均可以沿圓形薄板分解為切向和法相的荷載,則x,y,z方向的荷載可用函數表示,即
(27)
將上式(27)帶入式(24)中,得
(28)
在此處設Lagrange插值逼近的函數如下式:
(29)
式中:Pn(x) 表示正弦函數的n次Lagrange插值多項式;Ln(x)表示余弦函數的n次Lagrange插值多項式。
將式(29)帶入式(28)中,得
(30)
依據n次Lagrange中插值多項式函數的特點,式(30)荷載列陣在x方向上的積分項可展開為
(31)
式中:an和bn分別為插值多項式Pn(x)及Ln(x)中冪函數xn的系數。
因此,只要求出每一項進行累加即可得到目標函數在區間上的定積分。
由式(5)及(31)可以得到下式:
(32)
上式中矩陣Φ中的每一項是三次B樣條函數,由B樣條函數的數值計算特性可知,單獨對B樣條函數進行多重積分與多次求導是易于得出結果的,因此,可運用分部積分法建立下列積分公式對式(32)中每一項進行計算得到荷載列陣的值(依據本文需要,取到四次Lagrange中插值多項式函數),其中,每一項的數值積分方法如下所示。
在區間[0,a]上N等分,令t=(x/hx-i),通過換元法積分,則可得到下式:
(33)
(34)
再運用分部積分法對式(34)中每一項進行數值計算得到下式:
(35)
基于B樣條函數φn(x)的求積分公式便可計算式(35),積分公式如下所示:
(36)
式中:B樣條函數φn(x)右上角負號表示求積分次數。
在考慮逼近曲線的光滑性、連續性、易于編程等特點后,本文選用二次、三次及四次等不同Lagrange插值多項式逼近正弦及余弦函數。由于結構從坐標原點處開始,沿x正負方向荷載變化規律一致,其逼近區間選為[0,π/2],其插值節點與插值多項式如下表1、表2所示,插值多項式函數圖如圖2所示。

表1 正弦函數插值節點與插值多項式Tab.1 Sinusoidal function interpolation nodes and interpolation polynomials

表2 余弦函數插值節點與插值多項式Tab.2 Cosine function interpolation nodes and interpolation polynomials

圖2 插值多項式函數圖Fig.2 Interpolation polynomial approximation rendering
正弦函數插值多項式的具體表達式如下式所示:
(37)
余弦函數插值多項式的具體表達式如下式所示:
(38)
圓柱殼拱屋頂承受自重荷載。圓柱殼直線邊界自由。參考文獻中的算例,為同其他計算方法作比較,本文中算例曲線邊界約束設置為簡支,其結構及尺寸如圖3所示。依據結構對稱的受力特性,取對半結構進行計算,網格劃分為10×10的網格。本方法計算結果與相關文獻方法及精確解(扁殼解法)結果比較見表3和圖4,與精確解(扁殼解法)的相對誤差見表4。

表3 特殊點的位移和內力Tab.3 Displacement and internal force of special points

表4 特殊點位移和內力的相對誤差Tab.4 Relative error of displacement and internal force at special points

圖3 結構尺寸Fig.3 Structural dimension diagram

圖4 內力分布圖Fig.4 Internal force distribution diagram
綜合分析上述算例結果,如表4所示,二次、三次、四次插值方法在特殊點處的相對誤差均低于1.5%,均滿足工程要求的計算精度,且三次、四次插值方法的誤差低于二次插值方法,結合圖2中三次、四次插值函數逼近效果優于二次插值函數分析,表明兩者之間存在著一致性。如圖4所示,三種插值方法的位移與內力沿不同方向上的分布均與精確解(扁殼解法)一致。
本文發展了一種插值逼近形式的數值計算方法求解結構受任意荷載問題。基于樣條函數法及Lagrange插值法,推導了近似計算方法的數值積分公式,并通過解圓柱殼結構算例詳細展示了本方法的具體步驟。通過文中算例分析,得出以下結論:
1)本文使用三種不同插值多項式得到解與精確解(扁殼解法)及相關文獻解對比均具有較高計算精度,證明了本方法能夠有效分析結構受任意荷載問題。
2)近似計算方法以n次Lagrange插值多項式作為荷載函數,具有與原函數相同的光滑性、連續性、單調性等數學特性。本文基于Lagrange插值法的靈活性,取不同插值節點構造了三種不同的插值多項式,基于其具有的廣泛適用性,對兩種不同分布荷載均構造了插值函數,說明本方法對于任意荷載均具有較強適用性。
3)本文通過插值逼近的方法及多項式與B樣條函數乘積易于數值積分的思路,用近似函數替代原荷載函數,從而優化計算流程,解決復雜荷載不便于積分運算的問題。為未來研究結構受任意荷載問題提供新思路,具有廣闊的發展和應用前景。