海洋柔性立管防彎器結構非線性剛度拓撲優化1)

范志瑞 楊志勛 許 琦 蘇 琦 牛 斌 趙國忠

* (大連理工大學工程力學系,工業裝備結構分析國家重點實驗室,大連 116024)

? (哈爾濱工程大學機電工程學院,哈爾濱 150001)

** (大連理工大學機械工程學院,大連 116024)

引言

海洋柔性立管是連接海上浮體及海底井口的關鍵裝備,被譽為海洋浮式生產系統的“血管”[1-2].如圖1 所示,在服役過程中,由于風浪流及浮體運動的作用,立管與浮體間的連接部位會因彎曲變形過大而造成結構的失效破壞.因此,防彎器在防止立管的過彎失效,以及提高立管結構安全性等方面具有重要的意義[3].

圖1 柔性立管與防彎器結構示意圖Fig.1 Schematic diagram about the flexible riser and bend stiffener

國內外學者對防彎器結構的分析方法做了大量研究,并提出了不同的結構分析模型.如Deruntz[4]所建立的細長梁分析模型,Boet 和Out[5]提出的變截面梁理論模型,Tong等[6]提出的斜拉力懸臂梁模型,以及Drobyshevski[7]建立的異形梁模型等.上述研究中多采用梁模型對防彎器結構開展性能分析,可有效降低防彎器結構建模和分析的成本.然而,梁結構的分析精度依賴于實際防彎器結構的長細比,特別對于深水應用的防彎器,當長細比較小時,采用梁結構可能會造成較大的分析誤差.為了能夠更加準確地獲得防彎器性能,本文研究中將采用能夠反映防彎器截面/整體真實構型的平面/實體有限元模型.

結構優化是提高防彎器結構性能的有效設計方法.Tanaka等[8]以柔性立管的最大曲率,以及防彎器尺寸、最大應變和端部最大彎矩為約束,采用遺傳算法對防彎器進行了輕量化設計.Yang 和Wang[9]采用最優拉丁超立方方法建立了考慮材料、幾何和載荷不確定性的柔性立管疲勞分析代理模型,通過優化提高了防彎器結構的疲勞可靠性.Tang等[10]以質量和疲勞強度設置為目標,以柔性立管的最大曲率為約束,對防彎器的結構尺寸參數進行了多目標優化.然而,上述基于尺寸優化的結構設計所提供的設計自由度有限,難以充分發揮防彎器結構的性能.拓撲優化方法[11-14]能夠在給定的目標函數和約束條件下,獲得設計域內最優的材料分布,發現結構的創新構型,從而相比傳統的尺寸優化在更大程度上提高結構的性能.因此,本文研究中采用拓撲優化方法對防彎器進行結構優化.

實際的海洋工程中,防彎器由合成樹脂材料構造而成,同時在位運行中,由于惡劣的海洋環境載荷往往發生幾何大變形.因此在防彎器結構分析和優化設計中需要考慮幾何和材料的非線性影響[10].與線彈性結構的拓撲優化不同,在結構的非線性分析中,外載荷與位移響應之間表現為非線性映射關系.這使得結構剛度的描述既可以采用切線剛度陣,也可以采用割線剛度陣,且兩者所描述的結構剛度性能往往具有較大差異,如圖2 所示.Wallin等[15]的研究表明以結構割線剛度最大化和切線剛度最大化為目標時,優化結果具有很大的不同.Kemmler等[16]對非線性剛度優化中的目標函數進行了總結,其包括割線剛度(end-compliance)[15,17]、切線剛度(endstiffness)[15]、應變能[18]以及余能[19]等.實際工程中應根據實際需要選取恰當的目標函數.

圖2 切線剛度 Ktan 與割線剛度 Kcot 的定義,其中 F 和 u 分別為結構的載荷和位移Fig.2 The definition of tangent stiffness Ktan and secant stiffness Kcot,where F and u are the force and displacement of a structure,respectively

綜上所述,本文在考慮材料和幾何非線性情況下,以結構剛度最大化為目標,對柔性立管的防彎器結構進行拓撲優化.并在此基礎上,對不同設計假設下所得優化結果進行對比,探尋不同載荷方向對防彎器承載能力的影響.

1 海洋柔性立管防彎器拓撲優化問題描述

對圖1 中柔性立管及防彎器的工作過程進行抽象,可得圖3 所示簡化的力學模型.其中,模型的上端固定,在立管下端部施加波浪載荷.本文研究中涉及2D 和3D 兩種防彎器設計方案,為了避免兩種方案中復雜的載荷轉換關系,采用了Dirichlet 邊界條件來模擬波浪載荷的作用,即在y-z平面內,立管端部施加位移ap,且ap與y軸之間的夾角記為 θ .

圖3 海洋柔性立管、防彎器力學模型及邊界條件Fig.3 The mechanical model and boundary conditions of ocean flexible riser and the bend stiffener

關于Dirichlet 邊界條件下的線彈性結構拓撲優化問題,Niu等[20]采用了以結構支反力Rp與節點位移ap之間的點積進行定義的目標函數c,即

然而,如前所述,當結構呈現非線性特征時,結構剛度的描述比較復雜.以下將對上式的力學意義進行進一步分析.

圖4 給出了不同結構剛度下支反力隨位移變化示意圖.由圖可知,當位移加載結束時,曲線2 的支反力大于曲線1,即Rp2＞Rp1.同時,曲線2 所對應的結構割線剛度也大于曲線1 的對應值,即 α2＞α1.因此,就割線剛度角度而言,曲線2 所對應的結構剛度大于曲線1.與此同時,由于位移ap已經給定,式(1)所示目標函數的值與具體的加載過程無關,而僅與加載結束時的支反力水平相關.此時,式(1)與結構割線剛度所反映的結構性能一致,即支反力越大,結構的割線剛度越大.

圖4 式(1)所定義目標函數的力學含義Fig.4 Nature of the objective defined by Eq.(1)

此外,在圖4 中,曲線2 始終位于曲線1 的上方,因此曲線2 與坐標軸圍成的面積大于曲線1,即曲線2 所對應的結構應變能大于曲線1.并且,在加載結束點處,曲線2 的切線斜率大于曲線1.但是,上述情況并不表示式(1)所示目標函數可以反映結構的應變能和切線剛度水平.這是因為應變能水平與整體的加載過程相關,而切線剛度則與加載結束點附近的鄰域相關.如圖5 和圖6 所示,曲線2 的割線剛度大于曲線1,但是前者的應變能和切線剛度均小于后者.

圖5 應變能與式(1)所定義目標函數的不同Fig.5 The Difference between the strain energy and the objective defined by Eq.(1)

圖6 切線剛度與式(1)所定義目標函數的不同Fig.6 The difference between the tangent stiffness and the objective defined by Eq.(1)

綜上分析,本研究中采用Dirichlet 邊界條件下割線剛度作為優化目標,立管防彎器結構的拓撲優化列式可描述為

式中,z為所有設計域內單元密度設計變量集合,a為節點位移場,R(a,z) 為不平衡載荷向量,n為結構內單元的總數,g為材料用量體分比約束,ξ0為給定的材料用量體分比上限,ξ=V/V0為優化后的材料體分比.其中,V0為結構的總體積,V為優化后材料所占的體積.

2 Dirichlet 邊界條件下結構非線性分析

2.1 Dirichlet 邊界條件下的結構非線性分析方法

彈性體的虛功平衡方程可表示為

式中,u為彈性體的位移場,E為Green-Lagrange 應變張量,S為二階Piola-Kirchhoff 應力張量,t為作用在彈性體外邊界的牽引力.

對式(3)進行泰勒展開,并取一階截斷,可得

與此同時,引入Galerkin 方法對上式進行有限元離散化,即

式中,N為形函數矩陣,a為單元位移場u的節點值.

對式(4)右端第一項進行推導可得非線性迭代的平衡條件

式中,Fint和Fext分別為內外力載荷向量,其表達式如下

對式(4)右端第二項進行推導可得切線剛度陣的表達式為

式中,D為本構矩陣,R為結構應力相關矩陣,B和H均為單元的應變-位移矩陣.關于非線性分析的具體理論詳見文獻[21].

在本文研究中采用Dirichlet 邊界條件來模擬海洋波浪載荷的作用.如圖7 所示,結構的邊界 ?Ω 上施加位移條件的邊界包含兩種,即固定邊界 ?Ωz和非零位移邊界 ?Ωp,并且 ?Ωu=?Ωz∪?Ωp,?Ωz∩?Ωp=? .除Dirichlet 邊界外的其它自由邊界記作 ?Ωf,即?Ωf=?Ω-?Ωu.為了后續推導方便,在此約定下標“ u ”,“ f ”,“ z ”和“ p ”分別表示與邊界 ?Ωu,?Ωf,?Ωz和?Ωp相關的物理量.

圖7 Dirichlet 邊界條件結構分析示意圖Fig.7 Illustration of Dirichlet boundary conditions

Dirichlet 邊界條件作用下的Newton-Raphson 迭代過程中,邊界上無任何外載荷作用,即

嚴格意義上而言,Newton-Raphson 迭代過程中內外力的平衡條件R=0 是指所有位于邊界 ?Ωf上的節點所對應的不平衡載荷向量Rf取值為0,即

而對于施加位移約束的邊界 ?Ωu,內外力的差值則不為0.該差值即為結構的支反力

結合式(10),Dirichlet 邊界下的支反力可進一步表示為

此時,上述支反力又可進一步分解為位于邊界?Ωz上的Rz和邊界 ?Ωp上的Rp.相應地,不平衡載荷向量R可分解為以下3 部分

對R(a+da) 進行泰勒展開,并取一階截斷可得

式中,?R(a)/?a為結構的切線剛度陣K.且根據式(14),其可進一步寫為如下分塊矩陣格式

2.2 Neo-Hookean 材料分析

研究采用Neo-Hookean 超彈性材料來模擬防彎器和立管的材料非線性行為,其應變能密度函數表達式為

式中,J=det(F),F為應變梯度張量,C=FTF為右Cauchy 應變張量,K為體模量,G為剪切模量.

根據二階Piola-Kirchhoff 應力張量S和剛度張量D的定義可得兩者的表達式分別為

3 數值處理方法及靈敏度分析

3.1 對稱算子

在實際運行中,防彎器結構往往受到往復性波浪載荷的作用.為了保證防彎器具有較好的承載能力,通常要求防彎器為對稱結構.為了保證優化后的結構的對稱性,研究中引入了對稱算子M,此時初始設計變量z與變換后的設計變量之間滿足如下關系

以下引入3 自由度系統來簡要闡述對陣矩陣的構造方法.初始的設計變量所構成的向量z為

同時要求設計變量具有如下對稱關系

則對稱矩陣可以定義為如下兩種方式

在式(26)中,對稱算子的構造涉及所有設計變量,并且M=MT.在式(27)中,由于z3變量取值總與z1相等,因此在構造對稱算子M時,可略去z3.相應地,此時對稱算子M為非對稱矩陣.需注意的是,在實際優化中,M通常為稀疏矩陣.并且,對稱矩陣更容易實施轉置等運算.因此,研究中采用了式(26)所定義的對稱算子.

3.2 過濾及懲罰策略

在拓撲優化過程中,常會遇到棋盤格、網格依賴以及灰度單元等數值不穩定問題.對于棋盤格和網格依賴性問題,拓撲優化中常采用的方法包括靈敏度過濾[22-23]和密度過濾[24]等.密度過濾可有效克服優化過程中的棋盤格現象.在本文研究中,為了提高結構分析和優化的效率,基于MPI (massage passing interface)進程通信協議及PETSc庫[25-26]建立了并行分析和優化框架.然而,在并行計算過程中,單元信息及設計變量會在不同的進程中存儲.當采用傳統的密度過濾時,需要頻繁的進程通信來獲取臨近的單元信息,從而造成大量的進程通信成本.

為了克服傳統密度過濾方法在并行計算中的不足,本研究采用了基于Helmholtz-PDE 方程的過濾方法[27],以下簡稱PDE 過濾.PDE 過濾的控制方程可寫為如下格式

式中,為過濾前所有單元密度構成的向量,ρ 為過濾后所有節點密度構成的向量.Kf和Tf的定義如下

式中,Nf為形函數,Kd為過濾半徑相關的正定矩陣,d為過濾的維度,ri為維度i的過濾半徑;vi為過濾的方向矢量.當各個方向的過濾半徑相同,且vi為單位向量時,各向異性的過濾則退化為各向同性過濾.過濾后單元內各高斯積分點處的單元密度可通過下式求得

研究中,采用SIMP 方法對單元的材料屬性進行插值.在傳統線彈性拓撲優化中,SIMP 方法的插值格式通常表述為對材料彈性模量或者單元剛度矩陣的插值.在非線性優化問題中則略有區別,通常將其表述為對應變能密度函數的插值,即

為了改善優化過程中的灰度單元現象.在此,引入了Heaviside 懲罰策略[28-30],并且Heaviside 函數表達式為

圖8 Heaviside 函數的性質Fig.8 The characteristic of Heaviside function

綜合上述PDE 過濾、SIMP 插值以及Heaviside 懲罰,最終單元應變能密度的插值格式為

3.3 靈敏度分析

采用伴隨法對優化靈敏度進行推導,則引入伴隨向量 λ 后的目標函數可寫為

同時,聯合式(16)可知,伴隨向量 λ 可通過下式求得

此時,目標函數對單元密度設計變量的導數可寫為

考慮PDE 過濾時,上式可進一步整理為

為了靈敏度分析的求解,對上式中的物理量進行如下變換

此時,式(39)可進一步整理為

根據式(7)中內力Fint的表達式,可得

式中,?=Bκ .

進一步地,考慮優化中的SIMP 插值以及Heaviside 懲罰,上式可進一步寫為如下離散格式

聯合式(44)、式(50)和式(53)即可求得目標函數對設計變量的導數.進一步地,根據式(23)可得目標函數對z的導數為

4 數值算例

采用上述所建立非線性分析和優化理論對防彎器結構進行優化.算例中,所有參數均采用無量綱化處理.防彎器的彈性模量E1=210,柔性立管的彈性模量為E2=2100,兩種材料的泊松比均為μ=0.3 .在Newton-Raphson 迭代中,采用均勻的載荷增量步.為了保證迭代的穩定性,所有算例均經過預先試算來確定恰當載荷增量步數目.一般而言,載荷增量步數目的設置要保證每個增量步只需要3 至4 次迭代即可收斂.迭代收斂的條件為不平衡載荷向量的2-范數小于10-6,即 ‖R‖≤10-6.

為了保證優化過程的穩定性,采用了連續化懲罰策略.在優化前期,采用較小的懲罰指數,以降低優化過程陷入局部最優的概率.后期則采用較大的懲罰指數,以使得優化結果逐漸趨于0/1 化.具體的懲罰策略如下:(1) 初始的SIMP 懲罰指數p和Heaviside 懲罰指數 β 均設置為1.0;(2) 在后續優化迭代中,SIMP 懲罰指數每5 步迭代增加0.05;(3) 當SIMP 懲罰指數達到3.0 后,其保持不變;Heaviside懲罰指數將會每5 步增加0.5,直至其達到12.0.

4.1 算例1 防彎器2D 拓撲優化設計

本算例假設防彎器在位運行中會受到沿各個方向的波浪載荷的作用(各方向概率相等),因此要求防彎器為旋轉對稱結構.此時,防彎器的設計問題可退化為如圖9 所示的2D 平面問題.防彎器的尺寸為R1=1800,R2=300,L1=2400,L2=900 .其左端固定,右端截面上施加方向豎直向上、大小為ap=200的指定位移.為了保證防彎器與柔性立管的配合,在防彎器與立管的接觸區域設置厚度為30 的不可設計域.

圖9 算例1 中懸臂梁結構的設計域及邊界條件Fig.9 Geometry of the admissible domain and boundary conditions about the cantilever beam considered in example 1

有限元計算中,該結構被離散為20 400 個4 節點平面單元,MPI 進程數為10.在結構優化中,約束材料用量體分比上限 ξ0=0.5,PDE 過濾半徑為30.此外,采用式(23)和式(26)構建了對稱算子M,以保證優化結果的對稱性.

優化后防彎器的拓撲優化構型如圖10 所示.對圖10 所示構型進行旋轉操作,即可得3D 的防彎器構型,如圖11 所示.由圖9 可知,施加指定位移場ap后,防彎器的變形主要以彎曲為主,同時伴隨面內的剪切.因此,在優化后的結果中,材料主要分布在防彎器的外側,用于提高截面的慣性矩,進而增強結構的彎曲剛度.同時,在防彎器的內部與立管相連接的區域出現了交叉的斜撐結構,其主要用于抵抗防彎器內部出現的剪切變形.

圖10 算例1 中防彎器2D 拓撲優化結果Fig.10 The 2D optimized result of bend stiffener in example 1

圖11 算例1 中旋轉操作后所得防彎器優化結果剖面圖Fig.11 Sectional view about the optimized results of bend stiffener in example 1 obtained by using rotation operation

4.2 算例2 防彎器3D 拓撲優化設計

考慮到由于波浪、海流等環境載荷在特定海域、特定浮式平臺處的方向性,使得作用在防彎器上的載荷往往也呈現特定的方向性.此時,最優的立管防彎器不再具有軸對稱特征,必須按照3D 結構開展非線性拓撲優化設計.在本算例的結構分析和優化中,仍然采用如圖9 所示幾何尺寸.由于波浪載荷的往復性,此時要求優化后的結果具有x-z平面的對稱性,如圖12 所示.與算例1 相同,結構左端固定,右端立管截面施加沿y軸豎直向上的指定位移場ap,大小為200.由于載荷、邊界條件及結構變形也具有x-y平面的對稱性,僅選取防彎器和立管中z≥0 的部分進行有限元分析和優化.

圖12 算例2 中所采用的結構有限元模型Fig.12 The FEM model used in example 2

有限元計算中,結構整體被離散為494 080 個8 節點實體單元,此例求解時設置MPI 進程數為50.需注意的是,對圖10 所示2D 優化結果進行旋轉操作形成3D 結構時,由于各平面單元在旋轉路徑上掃略所產生的體積不同,會使得圖11 所示3D 結構的材料用量體分比與圖10 有所不同.本算例中,優化所采用的材料用量體分比約束的上限與圖11 所示3D 構型的實測體分比保持一致.此外,PDE 過濾半徑設為30.同時,與算例1 類似,為了保證防彎器與立管的配合,圖12中的區域(含柔性立管區域)設定為不可設計域.

考慮載荷方向的防彎器優化結果如圖13 和圖14 所示.其中,圖14 中所示構型為圖13 所得防彎器構型經x-y平面對稱操作后得到.由圖14 可知,優化后的防彎器呈現出明顯的“H”型構型.與算例1 不同,該構型主要用于抵抗沿y軸方向的載荷.其中,位于上下兩側的翼板主要用于提高結構的彎曲剛度;而中間的腹板則用于抵抗彎曲變形過程所產生的剪切變形.

圖13 算例2 中防彎器3D 拓撲優化結果Fig.13 The 3D optimized result of bend stiffener in example 2

圖14 算例2 中對稱操作后所得優化結果Fig.14 The optimized results obtained by symmetry operation in example 2

以下將對圖11 和圖14 所示優化結果在不同波浪載荷方向下的承載能力進行對比.如第3 節所述,當ap的大小一定時,位移施加端的支反力越大,則對應結構的剛度也越大.因此,選用支反力指標對兩種設計方案進行對比.在對比分析中,立管右端y-z平面內所施加位移ap的幅值固定為200.提取立管右端面所有節點的支反力,并且以所提取支反力的合力大小為指標,對結構的承載能力進行量化.不同θ取值情況下,兩種設計方法的支反力對比如圖15 所示.

由圖15 可知,由于圖11 所示優化結果為旋轉對稱結構(各方向載荷概率相等),因此其對各個方向的波浪載荷具有相同的承載能力.而對于圖14 所示優化結果,在不同波浪載荷方向下的承載能力則呈現出明顯的不同.具體而言,如圖15 所示在相當大的角度范圍內(除 θ 位于 [ 80,100] 和 [ 260,280]),3D拓撲優化結果的承載能力則顯著優于2D 的優化設計;而當載荷方向位于剛度薄弱區域,即θ ∈[80,100]∪[260,280]范圍時,其承載能力弱于圖11 所示優化結果.此時,對防彎器的安裝提出了更嚴格的要求,即在安裝過程中,需要對海域波浪載荷方向進行預先勘測,避免波浪載荷方向應避開防彎器的剛度薄弱區域.

圖15 不同載荷方向(θ)下,算例1 和算例2 優化結果的剛度對比Fig.15 Stiffness comparison between the optimized results of example 1 and 2 with different loading directions (θ)

5 結論

本文建立了Dirichlet 邊界條件下考慮材料和幾何非線性的防彎器結構拓撲優化數學模型及框架.為了保證優化結果的對稱性,研究引入了對稱算子,并給出了對稱算子的具體構建方法.采用伴隨法進行了靈敏度分析,并且引入了并行計算方法,提高了結構分析和優化的效率.數值算例部分,給出了基于2D 和3D 拓撲優化的防彎器創新設計方案.通過對此兩種方案的承載能力進行了對比,主要得到以下結論:

(1) 當不考慮特定的波浪載荷方向時,可采用旋轉對稱假設的防彎器結構創新設計,其能夠保證在任意載荷方向下,防彎器均具有較好的承載能力.

(2) 當海域內的波浪載荷具有特定的方向時,可采用平面對稱假設的防彎器設計方案.需注意的是,安裝過程中載荷方向應避開防彎器的剛度薄弱區域.