解析型幾何非線性圓拱單元

丁 敏，石家華，王斌泰，王宏志，鄧 婷,2，羅 雙,3，蔣秀根

(1. 中國農(nóng)業(yè)大學水利與土木工程學院，北京 100083；2. 唐山市路北區(qū)鳳凰新城管委會，河北，唐山 064400；3. 中國農(nóng)業(yè)科學院基本建設局，北京 100081)

圓拱作為一種壓彎結構，由于軸線曲率的影響，使得軸力與彎矩相互耦合，其截面內(nèi)力以軸力為主、同時存在彎矩及剪力。一方面在實際工程中，由于拱結構的長細比較大[1]，截面剛度較小，為保證拱結構內(nèi)力和位移計算結果的精確性，研究中考慮桿件幾何非線性影響是十分必要的；另一方面由于拱截面內(nèi)力以軸壓力為主，彎矩及剪力為輔，拱發(fā)生失穩(wěn)破壞是一種常見的破壞形式，而拱結構失穩(wěn)破壞是典型的非線性問題。

拱的失穩(wěn)破壞具有突然性，如極值點失穩(wěn)、分岔失穩(wěn)，拱結構穩(wěn)定性是影響其安全性的重要因素之一[2?10]。圓拱由于受面內(nèi)荷載易發(fā)生平面內(nèi)失穩(wěn)，當拱在平面內(nèi)荷載作用下，所承受的荷載達到一定的臨界值時，拱的平衡狀態(tài)會發(fā)生突然性的改變，拱軸線離開原來的變形狀態(tài)，向平面撓曲狀態(tài)轉化，從而喪失穩(wěn)定性[11?13]。平面內(nèi)失穩(wěn)可分為平面內(nèi)分岔失穩(wěn)與平面內(nèi)極值點失穩(wěn)，而結構只有在徑向均布荷載作用下且忽略軸向變形時，才會發(fā)生分岔失穩(wěn)[14]。因此在分析圓拱結構幾何非線性問題時應考慮是否忽略軸向變形：忽略軸向變形的無壓縮模型可進行分岔失穩(wěn)分析，考慮軸向變形的可壓縮模型可進行極值點失穩(wěn)分析。

對于圓拱的幾何非線性問題，國內(nèi)外學者采用靜力法、能量法和有限元法等建立了各類幾何非線性圓拱模型，求解相應問題。Timoshenko[15]考慮幾何非線性對結構的影響，建立了受壓圓環(huán)和曲桿的屈曲分析模型，從物理方程出發(fā)得到了二階彎曲撓度位移微分方程，結合圓拱邊界條件求出了精確的圓拱分岔失穩(wěn)特征方程及臨界荷載。Dinnik[16]在均布靜水壓力作用下的圓拱研究中，運用力法先求出內(nèi)力的表達式，再根據(jù)物理方程、幾何方程求出位移表達式，進而根據(jù)邊界條件得到了不同支座情況下的失穩(wěn)臨界荷載。項海帆[14]分析了圓拱的平面內(nèi)屈曲問題，詳細列出了無鉸拱、兩鉸拱以及三鉸拱在不同弧心角時的臨界荷載系數(shù)。鄧婷[11]基于靜力法基本方程求得了考慮軸力產(chǎn)生的二階效應的幾何非線性圓拱位移及內(nèi)力解析表達式，以矩陣形式的直接剛度法構造了有限單元，并用靜力法進行了分岔失穩(wěn)分析。申波等[17]利用小變形理論及ANSYS 研究了均勻徑向壓力作用下深圓拱的平面內(nèi)穩(wěn)定性，從微段平衡方程入手，計算得到了軸向位移通解。程鵬、童根樹[18]依據(jù)平截面假定，采用有限變形理論中的應變位移關系，完整地考慮了橫向應力和剪應力的二階效應，用虛功原理推導了拱平面內(nèi)非線性分析的平衡方程，給出了圓弧拱平面內(nèi)彎曲失穩(wěn)一般理論?？梢?，研究主要集中在特定工況下圓拱結構穩(wěn)定問題的解析解方面，對于任意工況下圓拱結構的內(nèi)力、位移和穩(wěn)定分析尚缺乏統(tǒng)一模型。本課題組提出了解析形函數(shù)法[19?22]，該方法綜合了靜力法、能量法及有限元法的優(yōu)點，是一種高精度高效率的研究方法，可用較少的單元精確分析復雜結構復雜荷載的問題。

本文以圓拱為研究對象，采用解析形函數(shù)法，基于小變形、大撓度假定，考慮軸力引起桿件二階彎矩的幾何非線性影響，在靜力法位移方程的基礎上，構造了圓拱解析形函數(shù)，結合能量法及勢能變分原理，構造解析型幾何非線性圓拱單元，采用本單元模型對圓拱的平面內(nèi)分岔失穩(wěn)和極值點失穩(wěn)進行對比分析，驗證本構造單元的精確性、可行性和適用性。

1 模型及假定

1.1 模型參數(shù)

建立隨動極坐標系，圓拱坐標系及內(nèi)力方向如圖1所示。按右手螺旋法則定義坐標系，坐標原點位于拱左端，軸向坐標x為圓拱軸線方向，向右為正，z軸為圓拱徑向方向，位于圓拱軸線平面內(nèi)，垂直于軸線，指向圓心為正。

圓拱參數(shù)：軸線為平面圓弧，其半徑為r，弧長為l，圓心角為θ ，且l=rθ。

圓拱荷載：軸向均布荷載qx、徑向均布荷載qz、均布力矩my，xi處軸向集中力Pxi、徑向集中力Pzi及集中力矩Myi，與坐標軸方向一致為正。

圓拱位移：軸向位移u、徑向位移w和轉角β，與坐標軸方向一致為正。

圖 1 圓拱坐標系及內(nèi)力Fig.1 Coordinate system and internal force of circular arch

圓拱內(nèi)力：軸力N、剪力V和彎矩M，當截面外法線方向與坐標軸方向一致時，定義為與坐標軸方向一致為正；當截面外法線方向與坐標軸正向相反時，內(nèi)力定義為與坐標軸方向相反為負。

δe={u1w1β1u2w2β2}T;

Fe={N1V1M1N2V2M2}T。

1.2 基本假定

本文推導過程基于以下假定：

1)圓拱軸線為等曲率平曲圓??；

2)將圓拱視為細長桿，忽略剪切變形，可采用Euler-Bernoulli彎曲理論計算截面的彎曲應力及彎曲變形；

3)圓拱截面按等截面考慮，對于變截面圓拱，可將其劃分或簡化為若干個等截面圓拱段；

4)圓拱為線彈性變形，圓拱的抗壓剛度EA及抗彎剛度EI均為常數(shù)，圓拱截面內(nèi)力與變形成正比；

5)符合大撓度假定，考慮軸力產(chǎn)生的二階效應，即P-Δ效應。

2 單元位移形函數(shù)

2.1 位移控制方程

綜合圓拱平衡方程、幾何方程及物理方程[11]，可得位移控制方程：

1)軸向位移控制方程

式中：E為材料彈性模量；A為截面面積；I為截面慣性矩。

2)徑向位移控制方程

2.2 位移方程

對位移控制方程式(1)～式(3)求解關于徑向位移w和軸向位移u的齊次微分方程，可得[11]：

2.3 位移基函數(shù)向量及導函數(shù)轉換矩陣

1)位移基函數(shù)向量

將式(4)和式(5)寫成關于位移系數(shù)的向量表達，可得位移向量表達式：

式中：d為位移系數(shù)向量，且d={d1d2d3d4d5d6}T；fw、fu分別為徑向和軸向位移基函數(shù)向量，且：

2)導函數(shù)轉換矩陣

使用導函數(shù)轉換矩陣Z來表述基函數(shù)導數(shù)f′與基函數(shù)f的關系，有：

結合式(7)可得，徑向與軸向位移基函數(shù)向量的導數(shù)轉換矩陣一致，且為：

2.4 單元位移形函數(shù)

根據(jù)兩端節(jié)點幾何邊界條件，可得：

δe={u(0)w(0)β(0)u(l)w(l)β(l)}T

寫成：

式中，A為節(jié)點位移矩陣，且：

根據(jù)式(10)，其位移方程式(6)可改寫為：

定義w=Nwδe，u=Nuδe，可得：

式中，Nw、Nu分別為徑向和軸向位移形函數(shù)。

同時根據(jù)幾何方程[11]，代入位移形函數(shù)式(12)，可得軸向應變ε、彎曲轉角 β和彎曲曲率κ表達如下：

3 可壓縮圓拱勢能及單元

3.1 單元勢能一般式

平面圓拱結構承受荷載時，單元總勢能包括變形能和荷載勢能。變形能包括壓縮變形能和彎曲變形能；當考慮軸力產(chǎn)生的二階效應時，荷載勢能包括外荷載勢能和桿端軸力勢能。單元總勢能可寫成下列形式[23]：

式中：U為單元的變形能；Up為荷載勢能。

3.2 單元勢能變分

3.3 單元列式

1)等效平衡方程

將式(15)改寫，即可得以單元節(jié)點位移表達的圓拱等效平衡方程：

2)單元剛度矩陣

根據(jù)等效平衡方程式(15)與式(16)對應，可得單元剛度矩陣具體表達列式：

式中：KN為抗壓剛度矩陣；KM為抗彎剛度矩陣。且：

3)幾何剛度矩陣

幾何剛度矩陣可表達為：

4)等效節(jié)點力向量

等效節(jié)點向量可表達為：

4 無壓縮圓拱單元

在圓拱結構實際工程中，由于軸向壓力的存在，需考慮其穩(wěn)定問題，不考慮軸向變形的無壓縮圓拱模型可用于分析圓拱結構的分岔失穩(wěn)問題，因此本文在以上模型的基礎上，繼續(xù)構建無壓縮圓拱單元，為后續(xù)圓拱結構穩(wěn)定問題的研究做準備。

一般來說，構建無壓縮圓拱單元，可采用直接模型和退化模型兩種途徑。直接模型可基于本文模型的構建過程由 ε=0重新推導而得；退化模型可基于本文模型單元列式由EA=∞退化而得。如下是采用退化模型得到的無壓縮圓拱單元。

4.1 位移基函數(shù)及導數(shù)

1)位移參數(shù)

取EA=∞，代入位移參數(shù)

有：

2)位移基函數(shù)

取EA=∞，結合式(7)與式(21)，可得無壓縮圓拱的位移基函數(shù)為：

3)導數(shù)轉換矩陣

取EA=∞代入式(9)，可得無壓縮圓拱基函數(shù)的導數(shù)轉換矩陣為：

4.2 抗壓剛度矩陣

將式(18)抗壓剛度矩陣作如下改寫：

代入式(7)和式(9)計算，可得：

取EA=∞代入式(25)，可得到無壓縮圓拱抗壓剛度矩陣為：

4.3 單元列式

1)等效平衡方程

由式(16)，可得無壓縮圓拱等效平衡方程：

2)單元剛度矩陣

無壓縮模型單元剛度矩陣一般式形式上同式(17)，同時理論上取EA=∞，可得簡化式。

將式(22)、式(23)和式(26)代入式(17)和式(19)，可得無壓縮圓拱的單元剛度矩陣：

式中，KM表達式同式(19)。

3)等效節(jié)點力向量

將式(22)和式(23)代入式(20)，可得無壓縮圓拱等效節(jié)點力向量，其表達式同式(20)。

4.4 模型對比分析

經(jīng)按本文構建過程重新推導得到直接模型，通過對比可知，由本文模型理論退化得到的無壓縮圓拱單元模型與直接模型是一致的。在理論上，當EA=∞時，用本文模型計算無壓縮圓拱具有可行性和正確性。

退化模型與直接模型的結果一致，說明本文模型單元的準確性與適用性。

5 算例及對比

5.1 圓拱平面內(nèi)分岔失穩(wěn)分析

5.1.1 分岔失穩(wěn)分析方法

結構只有在徑向均布荷載qz作用下且忽略軸向變形時，才會發(fā)生分岔失穩(wěn)[14]。對于分岔失穩(wěn)，等效平衡方程式(27)有不確定解，其位移系數(shù)行列式為0，即：

根據(jù)式(29)可求得圓拱在各種約束作用下的平面內(nèi)分岔失穩(wěn)的最小特征值λcr。

5.1.2 分岔失穩(wěn)臨界荷載

由特征值λcr求解可得臨界軸力Ncr。對于圓拱分岔失穩(wěn)臨界狀態(tài)，有：Ncr=rqcr[11]，則分岔失穩(wěn)臨界荷載為：

定義臨界荷載系數(shù)k：

有臨界荷載：

5.1.3 不同結構形式的分岔失穩(wěn)分析

采用本單元模型分別對無鉸拱、單鉸拱、兩鉸拱、三鉸拱及左固右鉸拱五種結構形式的圓拱進行分岔失穩(wěn)分析，求得不同圓心角角度下的分岔失穩(wěn)最小特征值λcr和臨界荷載系數(shù)k，計算結果見表1。

表 1 圓拱分岔失穩(wěn)特征值及臨界荷載系數(shù)Table 1 Bifurcation instability eigenvalues and critical load coefficients of circular arches

5.1.4 對比與分析

為驗證構造單元的精確性，將本單元分岔失穩(wěn)分析結果與經(jīng)典模型(Timoshenko模型[15]、Dinnik模型[16])、靜力法模型[11]進行對比，給出與經(jīng)典模型臨界荷載系數(shù)的相對誤差，其中左固右鉸拱沒有經(jīng)典模型結果，給出與靜力法模型結果對比，對比結果見表2。

由表2看出，本文構造的非線性圓拱單元進行分岔失穩(wěn)分析得到的失穩(wěn)臨界荷載系數(shù)與經(jīng)典模型、靜力法模型計算結果一致，相對誤差為0%，這說明了本單元模型的精確性及運用本單元模型進行分岔失穩(wěn)分析的可行性。可見，采用靜力法獲得解析解的經(jīng)典模型，僅適用于特定工況下圓拱結構的平面內(nèi)分岔穩(wěn)定分析，而采用能量法建立的本文單元則適用于任意工況下圓拱結構的平面內(nèi)分岔穩(wěn)定分析。

同時可以看出各結構形式隨著圓心角角度增大，其失穩(wěn)臨界荷載系數(shù)變小；對于同一圓心角角度，上述五種結構形式，無鉸拱的失穩(wěn)臨界荷載系數(shù)最大，穩(wěn)定性最好，符合工程實際。

表 2 不同模型圓拱的分岔失穩(wěn)臨界荷載系數(shù)對比Table 2 Comparison of critical load coefficients for bifurcation instability of circular arches with different models

5.2 圓拱平面外極值點失穩(wěn)分析

對徑向均布荷載作用下的無鉸拱進行平面內(nèi)極值點失穩(wěn)分析。假定無鉸拱圓心角為120°，半徑為5 m，截面尺寸為400 mm×400 mm，彈性模量為206 GPa，受徑向均布荷載qz=10 kN/m。通過改變軸力特征值大小，求出跨中撓度，得到圓拱軸力特征值與跨中撓度變化曲線如圖2所示。

圖 2 軸力特征值-跨中撓度曲線Fig.2 Axial force eigenvalue vs. mid-span deflection curve

由圖2可見，當荷載達到臨界荷載時，圓拱出現(xiàn)極值點失穩(wěn)現(xiàn)象，且圓拱極值點失穩(wěn)臨界荷載與圓拱平面內(nèi)分岔失穩(wěn)臨界荷載值相同。

上述所有算例表明：本文構造單元既可用于圓拱幾何非線性的位移和內(nèi)力分析，也可用于圓拱的平面內(nèi)各類穩(wěn)定分析。

6 結論

本文采用解析法構造了幾何非線性圓拱單元，研究了圓拱結構幾何非線性內(nèi)力和位移計算以及平面內(nèi)穩(wěn)定分析問題，主要結論如下：

(1)基于圓拱幾何非線性靜力分析模型，構造了圓拱單元解析位移形函數(shù)。

(2)利用能量法，建立圓拱單元的勢能，采用勢能變分原理，得出等效平衡方程，給出單元剛度矩陣、幾何剛度矩陣及等效節(jié)點力向量表達，構造了解析型幾何非線性圓拱單元。

(3)通過理論上取EA=∞對本模型進行退化，得到不考慮軸向變形的無壓縮圓拱模型，進而得到了解析型無壓縮幾何非線性圓拱單元，可應用于平面內(nèi)分岔失穩(wěn)分析。

(4)本文構造的解析型非線性圓拱單元可用于圓拱的平面內(nèi)分岔失穩(wěn)和極值點失穩(wěn)分析。算例表明：本研究單元模型進行分岔失穩(wěn)分析得到的失穩(wěn)臨界荷載系數(shù)與經(jīng)典模型(Timoshenko模型、Dinnik模型)、靜力法模型計算結果一致，誤差為0%，這說明了本單元模型的精確性及運用本單元模型進行分岔失穩(wěn)分析的可行性；平面內(nèi)分岔失穩(wěn)的臨界荷載與平面內(nèi)極值點失穩(wěn)的臨界荷載一致；本研究模型能夠完整的呈現(xiàn)圓拱平面內(nèi)由位移非線性發(fā)展到極值點失穩(wěn)的過程。