基于半數(shù)值方法的限彎器剪彎剛度分析

丁樂聲，陳金龍，陳 瀟，郭 豪，張 聰

(1. 南方海洋科學(xué)與工程廣東省實(shí)驗(yàn)室(湛江)，廣東 湛江 524006; 2. 大連理工大學(xué) 寧波研究院，浙江 寧波 315700; 3. 江蘇亨通海洋光網(wǎng)系統(tǒng)有限公司，江蘇 蘇州 215000; 4. 華中科技大學(xué) 船舶與海洋工程學(xué)院， 湖北 武漢 430074)

柔性管纜是海洋漁業(yè)養(yǎng)殖[1-2]、油氣和風(fēng)電[3]等海洋資源開發(fā)的關(guān)鍵裝備之一。其在多種載荷共同作用下易產(chǎn)生過大曲率，從而導(dǎo)致彎曲失效[4]。因此美國石油協(xié)會(API)發(fā)表過相關(guān)規(guī)范[5]，其中指出柔性管纜易發(fā)生彎曲失效的位置需要安裝限彎器。并且需對不同工況的限彎器進(jìn)行整體分析[6]，以確定其抵抗變形能力符合使用要求。

限彎器是由多組同尺寸變徑半圓筒式子結(jié)構(gòu)相互嵌套而成，當(dāng)其子結(jié)構(gòu)受力將相互鎖合，從而限制彎曲變形進(jìn)一步發(fā)生[7]。在限彎器鎖合過程中，其鎖合結(jié)構(gòu)具有接觸非線性，進(jìn)而導(dǎo)致剪彎剛度分析困難。針對該類問題，部分研究人員進(jìn)行了相關(guān)研究。其中理論方面，Souza和Ramos[8]建立了剪彎荷載下連續(xù)異形梁結(jié)構(gòu)受彎的理論模型，但無法應(yīng)用于限彎器這一包含接觸非線性的梁結(jié)構(gòu)。張聰[9]建立了限彎器剪彎剛度的整體數(shù)值分析方法。Noh和Ha[10]運(yùn)用ABAQUS的接觸分析模塊[11]進(jìn)行了限彎器與纜的多體接觸分析，并進(jìn)行限彎器結(jié)構(gòu)拓?fù)鋬?yōu)化設(shè)計。

綜上所述，目前限彎器這一包含接觸非線性梁結(jié)構(gòu)的剪彎剛度分析計算都基于整體數(shù)值方法，無法直接給出限彎器截面尺寸對剪彎剛度的影響，導(dǎo)致限彎器剪彎剛度設(shè)計多基于經(jīng)驗(yàn)。為此提出剪彎剛度半數(shù)值分析方法以便通過公式直接分析部分參數(shù)與剪彎剛度關(guān)系，并滿足限彎器剪彎剛度迭代設(shè)計的快速分析需求。先基于懸臂梁方程，建立包含鎖合結(jié)構(gòu)扭轉(zhuǎn)剛度的剪彎剛度理論模型；然后通過構(gòu)建考慮接觸邊界的鎖合結(jié)構(gòu)數(shù)值模型，得到其扭轉(zhuǎn)剛度曲線；最后結(jié)合理論分析和數(shù)值計算求得考慮接觸非線性的限彎器剪彎剛度曲線，并同傳統(tǒng)的整體數(shù)值方法相比較，驗(yàn)證該半數(shù)值分析方法的有效性。

1 剪彎剛度理論分析

如上所述，限彎器是典型的非線性短梁結(jié)構(gòu)。傳統(tǒng)的限彎器剪彎剛度采用整體數(shù)值分析，先建立三維懸臂梁有限元模型計算該結(jié)構(gòu)在剪力作用下的端部撓度；然后使用懸臂梁撓度方程計算出剪彎剛度。如圖1所示，限彎器的整體結(jié)構(gòu)可分為子梁結(jié)構(gòu)和鎖合結(jié)構(gòu)；在彎曲過程中子梁結(jié)構(gòu)主要發(fā)生撓曲變形，鎖合結(jié)構(gòu)主要發(fā)生扭轉(zhuǎn)變形[12]，其中鎖合結(jié)構(gòu)發(fā)生的扭轉(zhuǎn)會導(dǎo)致梁結(jié)構(gòu)端部發(fā)生垂直于軸線方向的位移進(jìn)而使結(jié)構(gòu)撓度增加；因此限彎器的剪彎剛度受子梁結(jié)構(gòu)彎曲剛度和鎖合結(jié)構(gòu)扭轉(zhuǎn)剛度共同影響。

進(jìn)行理論分析時將限彎器結(jié)構(gòu)簡化為由兩部分組成的端部受剪力P的周期性懸臂短梁結(jié)構(gòu)，其中直線為子梁結(jié)構(gòu)、圓點(diǎn)為鎖合結(jié)構(gòu)，如圖2所示；限彎器端部的總撓度w由所有子梁結(jié)構(gòu)撓度w1與所有鎖合結(jié)構(gòu)扭轉(zhuǎn)所對應(yīng)撓度w2組成，即w=w1+w2。子梁結(jié)構(gòu)是典型的直梁，可引入小變形與直線梁假設(shè)；同時考慮其截面剪切變形，由Timoshenko梁理論[13-14]進(jìn)行懸臂梁結(jié)構(gòu)撓度計算，如式(1)所示。值得注意的是限彎器結(jié)構(gòu)中的鎖合結(jié)構(gòu)會使整體產(chǎn)生轉(zhuǎn)角，導(dǎo)致在引入直線梁假設(shè)后基于長度計算力矩時產(chǎn)生誤差，但當(dāng)其圓心角較小時該影響可忽略[15]。

圖2 限彎器結(jié)構(gòu)簡化懸臂梁模型Fig. 2 Simplified cantilever beam model with bend restrictor

(1)

式中：w1為所有子梁結(jié)構(gòu)撓度，P為端部剪力，E為材料彈性模量，n為子梁結(jié)構(gòu)節(jié)數(shù)，l為限彎器有效節(jié)長，I為梁截面慣性矩，c為剪切系數(shù)，Mr為剪切模量，A為梁截面面積。若忽略梁截面剪切變形，則式(1)退化為基于Euler-Bernoulli梁理論的懸臂梁撓度方程(2)：

(2)

(3)

(4)

n節(jié)子梁結(jié)構(gòu)需要n-1節(jié)鎖合結(jié)構(gòu)連接，進(jìn)一步求得有鎖合結(jié)構(gòu)扭轉(zhuǎn)所對應(yīng)撓度w2。

(5)

(6)

同理，將式(2)、(5)聯(lián)立，并也代入總撓度公式和等效剪彎剛度方程，可得到基于Euler-Bernoulli梁理論的忽略梁截面剪切變形的剪彎剛度K計算方法。

(7)

對比可知，方程(6)僅分母比方程(7)多一項(xiàng)，該項(xiàng)為短梁截面的剪切應(yīng)變造成，并恒為正值；因此基于Timoshenko梁理論的計算方程(6)更適用于短梁結(jié)構(gòu)，并且結(jié)果更為保守。同時也需要注意，以上理論模型是引入了直線梁假設(shè)和小變形假設(shè)建立的，限彎器本身的曲率與變形將導(dǎo)致該模型存在誤差。

2 限彎器剪彎剛度計算

2.1 算例幾何與材料參數(shù)

限彎器在進(jìn)行測試和使用時常因節(jié)省成本而選擇較短的組裝長度，以5節(jié)典型限彎器子結(jié)構(gòu)為例，子結(jié)構(gòu)形狀如圖3所示。使用過程中，限彎器子結(jié)構(gòu)頭尾相互嵌套于海纜外部，以提高局部剪彎剛度。構(gòu)件的幾何參數(shù)以及材料力學(xué)性能參數(shù)在表1中給出，并以5節(jié)限彎器子結(jié)構(gòu)嵌套為例進(jìn)行計算。

圖3 限彎器算例尺寸Fig. 3 Bend restrictor calculation example

表1 限彎器結(jié)構(gòu)參數(shù)Tab. 1 Bend restrictor structure parameter

2.2 半數(shù)值分析

根據(jù)理論模型，剪彎剛度求解需要鎖合結(jié)構(gòu)扭轉(zhuǎn)角，但限彎器中的每個鎖合結(jié)構(gòu)包含兩個接觸邊界，使其計算較為復(fù)雜。因此使用半數(shù)值方法進(jìn)行剪彎剛度分析，具體的流程[16]如圖4所示，先基于有限元軟件ABAQUS建立限彎器鎖合結(jié)構(gòu)局部三維模型，并計算提取不同力矩下鎖合結(jié)構(gòu)扭轉(zhuǎn)角度；同時基于理論分析中的式(1)或(2)使用不同梁理論計算子梁結(jié)構(gòu)彎曲撓度；再進(jìn)一步通過式(3)分析不同力矩下的鎖合結(jié)構(gòu)扭轉(zhuǎn)剛度，通過式(5)分析鎖合結(jié)構(gòu)所對應(yīng)的撓度，進(jìn)而算出限彎器端部撓度；最后再將其與各參數(shù)代入式(6)或(7)基于Matlab軟件計算限彎器剪彎剛度。

圖4 剪彎剛度半數(shù)值分析流程Fig. 4 Shear bend stiffness analysis process

限彎器鎖合結(jié)構(gòu)可簡化為平面對稱結(jié)構(gòu)，因此建立半體三維模型[2-3, 9]進(jìn)行數(shù)值分析以減少計算量。模型在鎖合結(jié)構(gòu)末端設(shè)置固定約束，并在半體模型軸對稱面設(shè)置對稱約束；在進(jìn)行鎖合結(jié)構(gòu)局部數(shù)值分析時使用扭矩來等效5節(jié)限彎器端部受剪力造成的力矩，分析中限彎器端部受到的最大剪力為20 kN，則等效的最大扭矩為21.6 kN·m；又由于采用半體三維模型，實(shí)際該模型端部z軸方向施加最大扭矩的一半；具體如圖5所示。

圖5 鎖合結(jié)構(gòu)相互作用及邊界條件設(shè)置Fig. 5 Interaction and boundary condition setting of locking structure

該模型有兩個接觸對，法向行為選擇硬接觸算法，切向行為則選擇罰函數(shù)算法。網(wǎng)格劃分單元選擇分析耗時短的C3D8R單元，共3 980個。該模型每迭代100次的計算時間為96 s。通過有限元分析方法不僅可以直接計算鎖合結(jié)構(gòu)扭轉(zhuǎn)變形進(jìn)而計算彎曲剛度，還可以直接讀取力矩荷載達(dá)到21.6 kN·m時鎖合結(jié)構(gòu)各接觸面的接觸應(yīng)力云圖，如圖6所示。從圖6中可知，鎖合結(jié)構(gòu)接觸面較小且接觸壓力分布不均勻[2]，這都表明該處有較強(qiáng)的接觸非線性。

圖6 鎖合結(jié)構(gòu)各部分接觸壓應(yīng)力云圖Fig. 6 Contact compressive stress cloud diagram of each part of locking structure

從分析結(jié)果中可提取不同力矩的鎖合結(jié)構(gòu)轉(zhuǎn)角，進(jìn)一步可基于數(shù)值分析得到力矩—轉(zhuǎn)角關(guān)系以及式(3)計算限彎器鎖合結(jié)構(gòu)非線性扭轉(zhuǎn)剛度，如圖7所示。圖7(a)力矩與轉(zhuǎn)角的關(guān)系曲線與切線相比，具有明顯的非線性；從幾個給定力矩條件下的母頭端部后側(cè)接觸壓力云圖可以看出在力矩較小時接觸面積隨著力矩快速增加；而力矩超過10 kN·m后接觸面趨于穩(wěn)定，隨著力矩的增加接觸非線性減弱，力的傳遞主要體現(xiàn)在接觸壓力的增加。圖7(b)中扭轉(zhuǎn)剛度的變化趨勢進(jìn)一步說明鎖合結(jié)構(gòu)接觸非線性在加載初期較為明顯，隨著荷載增加非線性逐漸減弱。

圖7 限彎器鎖合結(jié)構(gòu)轉(zhuǎn)角與剛度Fig. 7 Torsion angle and rigidity of the bend restrictor locking structure

在半數(shù)值分析中為考慮鎖合結(jié)構(gòu)的扭轉(zhuǎn)剛度非線性，采用四階多項(xiàng)式對局部數(shù)值分析的結(jié)果進(jìn)行擬合，得到R平方系數(shù)[17](確定系數(shù))為0.997的扭轉(zhuǎn)剛度函數(shù)G(M)；其擬合結(jié)果在圖7(b)中給出。最后根據(jù)計算流程，將扭轉(zhuǎn)剛度函數(shù)G(M)與計算方程(1)、(5)或(2)、(5)分別聯(lián)立得到端部撓度，進(jìn)一步代入等效剪彎剛度方程即可獲得基于式(6)與(7)兩種理論的限彎器剪彎剛度。

2.3 整體數(shù)值分析

傳統(tǒng)方法中限彎器的剪彎剛度采用整體數(shù)值分析[18-19]；利用結(jié)構(gòu)對稱性建立上述算例的半體三維模型[5]，選擇與半數(shù)值分析方法相同的懸臂梁剪力加載，在根部施加完全固定約束，在對稱面設(shè)置對稱約束，并對半體模型端部施加10 kN的集中力，進(jìn)而模擬出實(shí)際結(jié)構(gòu)根部21.6 kN·m的力矩，如圖8所示。

圖8 限彎器相互作用及邊界條件設(shè)置Fig. 8 Interaction and boundary condition setting of bend restrictor

模型共有8個接觸對，接觸算法以及單元屬性與限彎器鎖合結(jié)構(gòu)數(shù)值分析模型相同；模型共包含28 390個C3D8R單元，是限彎器鎖合結(jié)構(gòu)模型的7倍，其中鎖合結(jié)構(gòu)處的單元尺寸和分布與鎖合結(jié)構(gòu)局部數(shù)值分析模型相同。整體數(shù)值分析可得到限彎器剪彎剛度用來驗(yàn)證半數(shù)值分析方法。并且整體數(shù)值分析能計算根部力矩和端部轉(zhuǎn)角關(guān)系，如圖9所示，進(jìn)一步說明了數(shù)值方法可以有效分析接觸非線性問題。

圖9 限彎器根部力矩—端部轉(zhuǎn)角曲線Fig. 9 Curve of bend restrictor root moment-end rotation angle

通過整體數(shù)值分析還可以進(jìn)一步分析小變形假設(shè)與直線梁假設(shè)對該結(jié)構(gòu)的影響。限彎器裝配體最大對數(shù)應(yīng)變和撓度分別如圖10、圖11所示。從圖10中可以看出該結(jié)構(gòu)整體應(yīng)變較低，應(yīng)變峰值出現(xiàn)在靠近懸臂梁末端的子梁結(jié)構(gòu)與鎖合結(jié)構(gòu)交界的應(yīng)力集中處，且均小于5%；從圖11中可以看出限彎器裝配體結(jié)構(gòu)最大撓度為140.5 mm，比其整體結(jié)構(gòu)長度1 174 mm低一個數(shù)量級。應(yīng)變所代表的局部變形與撓度代表的整體變形遠(yuǎn)小于結(jié)構(gòu)原始尺寸[20]，因此理論分析引入小變形假設(shè)所造成的誤差是有限的。限彎器裝配體中子梁結(jié)構(gòu)是典型的沿中性層對稱的直梁，直線梁假設(shè)導(dǎo)致的誤差主要由分析整體結(jié)構(gòu)力矩時產(chǎn)生，但該結(jié)構(gòu)對應(yīng)的圓心角為12.8°，而當(dāng)圓心角小于22°30′時曲線梁荷載可以近似的按直梁計算[15]，因此理論分析引入直線梁假設(shè)所造成的誤差也是有限的。

圖10 限彎器對數(shù)應(yīng)變云圖Fig. 10 Logarithmic strain cloud diagram of bend restrictor

圖11 限彎器撓度云圖Fig. 11 Deflection cloud diagram of bend restrictor

3 各算法對比

3.1 計算效率對比

限彎器設(shè)計過程中常進(jìn)行迭代計算，以獲得較優(yōu)的設(shè)計[10, 21]。因此計算效率是限彎器分析方法的關(guān)鍵指標(biāo)。上文所述半數(shù)值分析與整體數(shù)值分析均在i3-8350k CPU中進(jìn)行，其耗時于表2給出?？芍狤uler-Bernoulli梁理論與Timoshenko梁理論2種計算方法效率極其接近，而半數(shù)值分析效率遠(yuǎn)高于整體數(shù)值分析。

表2 各算法計算效率對比Tab. 2 Calculation efficiency comparision of various algorithms

3.2 計算結(jié)果對比

剪彎剛度是限彎器的關(guān)鍵指標(biāo)，傳統(tǒng)計算采用整體數(shù)值方法；同時其也可基于上文所述的半數(shù)值分析流程計算。通過不同的計算方法得到了相應(yīng)的限彎器剪彎剛度曲線如圖12所示，可知基于Euler-Bernoulli梁理論與Timoshenko梁理論的兩組不同半數(shù)值計算結(jié)果與傳統(tǒng)數(shù)值方法的平均相對誤差均小于6%，并且兩種半數(shù)值方法均可以較好地反映出限彎器鎖合結(jié)構(gòu)的接觸非線性，Timoshenko梁理論對應(yīng)的算法由于考慮了短梁的剪切變形更加保守。

圖12 限彎器剪彎剛度Fig. 12 Shear bend stiffness of bend restrictor

半數(shù)值方法的誤差主要由理論模型中引入的小變形假設(shè)和直線梁假設(shè)引起，進(jìn)一步計算基于不同梁理論的兩種半數(shù)值方法與傳統(tǒng)數(shù)值方法的誤差，于表3中給出；證明了Timoshenko梁理論對應(yīng)的算法比Euler-Bernoulli梁理論準(zhǔn)確，且Timoshenko梁理論的平均相對誤差為3.398%。

表3 兩種理論算法誤差對比Tab. 3 Comparison of two theoretical calculation errors

4 結(jié) 語

基于Timoshenko梁理論和Euler-Bernoulli梁理論剪力限彎器剪彎剛度半數(shù)值分析方法，并采用該方法和傳統(tǒng)的整體數(shù)值方法分別進(jìn)行了剪彎剛度分析并比較，得到如下結(jié)論：

1) 在限彎器剪彎剛度分析時半數(shù)值方法的計算效率遠(yuǎn)高于傳統(tǒng)數(shù)值方法；并且基于不同梁理論半數(shù)值方法的平均相對誤差均小于6%。

2) 限彎器的鎖合結(jié)構(gòu)處受到接觸邊界的影響呈現(xiàn)接觸非線性，這對限彎器鎖合結(jié)構(gòu)剛度造成很大的影響；并且其非線性在加載初期較為明顯，隨著荷載增加逐漸減弱。

3) 限彎器的子梁結(jié)構(gòu)是典型的短梁結(jié)構(gòu)，因此基于Timoshenko梁理論的半數(shù)值分析更適用于限彎器剪彎剛度分析，結(jié)果誤差更小且更保守。

綜上所述，基于半數(shù)值方法可以高效地分析限彎器剪彎剛度，其結(jié)果可與數(shù)值方法相印證，為相關(guān)分析與設(shè)計提供參考。但是本文中所使用的梁理論模型需引入直線梁假設(shè)和小變形假設(shè)，仍對分析結(jié)果造成一定誤差，可以在未來進(jìn)一步研究。