考慮內壓影響的加筋柱殼高效屈曲優化

楊木森，張 可，吳 浩，郝 鵬，王 博，田 闊

(1.大連理工大學工程力學系工業裝備結構分析國家重點實驗室，大連 116024;2.北京宇航系統工程研究所，北京 100076)

0 引言

網格加筋柱殼具有高比強度和剛度，被廣泛應用于大直徑運載火箭中，在大直徑運載火箭中網格加筋殼的占重比也越來越高[1-2]。雖然網格加筋殼具有較好的承載性能，但結構的復雜也使網格加筋殼的失效模式變得復雜[2-3]。大直徑運載火箭在飛行過程中，不斷變化的軸內壓耦合工況是火箭貯箱的典型服役工況[4]，整體失穩、蒙皮失穩、筋條失穩等均為其典型失效模式。隨著箭體結構尺寸跨越式增大，采用有限元法進行精細模型屈曲分析的計算成本激增，不再適合大直徑運載火箭的初始設計。因此，考慮內壓影響的加筋柱殼高效屈曲優化對大直徑運載火箭貯箱初始設計具有重要意義[3]。

運載火箭在發射、飛行過程中，其推進劑貯箱承受著發動機推進的軸壓和箱內燃料液體的內壓耦合作用，對這種軸內壓耦合的工況，學者們研究了不同參數對柱殼結構失穩模式的影響，包括內壓大小、材料強度及徑厚比等。Abdi等[5]開展了3種不同加筋結構的失穩模式研究，認為加筋密度與結構的可靠性有直接關系。Rotter[6-7]研究了內壓對柱殼貯箱底部的影響，認為內壓會導致柱殼基部嚴重彎曲，發生前期彈塑性屈曲失效，從而形成“象腳型”屈曲。Khelil[8]研究了內壓與非均勻軸壓耦合工況下各向同性殼的臨界應力。Franzoni等[9]提出了一種以振動相關技術為基礎的非破壞性試驗流程來評估含內壓的正置正交加筋柱殼軸向屈曲載荷。因此，內壓對貯箱加筋柱殼結構承載力的影響規律尚不明確，傳統均勻加筋筒段難以適用于復雜的載荷工況，亟需發展適用內壓工況的高承載效率的加筋策略。

針對加筋柱殼結構，等效剛度法(Smeared Stiffener Methond)是最常用的快速計算方法，該方法具有較高的分析效率[10]，其思路基于解析法對筋條進行剛度等效，再將筋條的等效剛度系數與蒙皮的剛度系數疊加，代入瑞利-里茲(Rayleigh-Ritz)公式計算得出整體屈曲載荷。傳統的等效剛度法往往會導致預測屈曲結果普遍偏大，同時，對于不同的加筋構型，等效剛度法需要進行不同的解析推導，推導過程復雜，增加了選型與優化的難度，不適用于創新構型[11]。筆者先前工作基于漸近均勻化方法和瑞利-里茲法，建立了基于NIAH(Numerical Implementation Asymptotic Homogenization)的一種快速屈曲分析方法，該方法基于加筋單胞有限元模型進行等效，無須進行解析推導，相較于等效剛度法具有更高的等效精度[3,12]，但該方法僅適用于純軸壓下的加筋柱殼。目前,尚未有學者提出針對軸內壓耦合工況的加筋柱殼承載力快速計算方法。

本文首先設計了非等厚筋條的加筋構型，并基于漸近均勻化法和瑞利-里茲法，建立了考慮內壓影響的網格加筋柱殼整體失穩載荷快速計算方法，并與傳統等效剛度法對比了計算效率和計算精度；然后，基于自適應加點的代理模型開展給定質量約束的大直徑網格加筋柱殼極限承載力最大化設計；最后，分析了內壓對同等網格參數(最優設計)的加筋柱殼整體失穩載荷的影響規律，并從失穩波形的角度解釋了該規律的機理。

1 大直徑加筋柱殼承載力快速計算方法

1.1 網格加筋柱殼整體失穩載荷的快速計算方法

針對軸內壓下的網格加筋柱殼結構整體屈曲失穩載荷，求解流程如圖1所示，具體過程如下[1]：

第一步，選擇網格加筋柱殼的代表性單胞，并建立有限元模型；

第二步，基于漸近均勻化方法計算典型單胞的等效剛度系數；

第三步，將單胞等效剛度系數代入至Rayleigh-Ritz公式，計算得出整體屈曲載荷值和屈曲模態波數。

該方法不需要解析推導代表性單胞剛度陣，并適用于任意周期性的網格加筋類型，有限元分析僅計算單胞模型，較精細模型分析將大幅提高計算效率。

圖1 基于NIAH的快速計算方法的流程圖

針對圓筒殼在均勻軸壓下的穩定性問題，學者們已做了大量的理論研究和試驗研。針對正交各向異性筒殼結構，根據圓筒殼穩定性基本方程組[13]

(1)

K11=A11α2+A66β2

K22=A22β2+A66α2

K12=(A12+A66)αβ

K66=D11α4+2(D12+2D66)α2β2+

(2)

對于任意的正整數m和n，數值求解N的最小值Nmin，即為單位長度軸壓臨界軸力

Nij=Nmin

(3)

N=

(4)

同理可求得初始內壓工況下的單位長度軸壓臨界軸力

Nij=Nmin

(5)

1.2 網格加筋柱殼局部載荷的快速計算方法

在工程實際中，網格加筋柱殼失效模式多種多樣，蒙皮局部失穩和筋條局部失穩是兩種典型的局部失效模式。本節基于板殼屈曲理論，發展了非等厚筋條加筋柱殼的蒙皮局部穩定性和筋條局部穩定性的承載力理論計算公式。

1.2.1 蒙皮局部穩定性

當加筋類型為正置正交時，根據彈性力學四邊簡支板軸壓屈曲理論[14]，蒙皮局部失穩時的單位長度軸力為

(6)

(7)

當加筋類型為等邊三角形時，由文獻[15]可知蒙皮單位長度臨界軸力為

(8)

式中，C1=9.72。

考慮到豎筋和斜筋寬度不相同，如圖2所示，由文獻[15]可知斜筋的等效面積為B=Asin3θ，θ=30°，那么，加筋柱殼的單位長度軸力為蒙皮受力和筋條受力的總和

(9)

Tlj=2πRNlj

(10)

圖2 等邊三角單胞示意圖

1.2.2 筋條局部穩定性

簡化筋條軸壓局部屈曲問題，將豎筋條局部失穩看作三邊簡支一邊自由的矩形板軸壓屈曲問題，其筋條局部失穩載荷為

(11)

當加筋類型為正置正交時，豎筋屈曲失穩的單位長度軸力為

(12)

Tlj=2πRNlj

(13)

當加筋類型為等邊三角形時，考慮到豎筋和斜筋寬度不相同，那么

(14)

因此，加筋柱殼筋條局部失穩時的整體臨界軸力為

Tlj=2πRNlj

(15)

1.3 快速計算方法精度驗證

為了驗證1.1節中提出的軸內壓下的加筋柱殼快速計算方法預測精度，本節采用典型等邊三角網格加筋殼模型，模型幾何參數：直徑D=4000mm，長度L=6000mm，蒙皮厚度ts=5.5mm，筋條高度H=44.5mm，筋條寬度tw=13mm，筋條間距bs=135mm。加筋柱殼的整體數值模型采用有限元仿真軟件ABAQUS建立，蒙皮網格單元大小為50mm，筋條高度方向離散兩個單元。設置底部簡支、頂部簡支且放松軸向位移的邊界條件，計算線性屈曲載荷。

表1給出了軸內壓下的加筋柱殼快速計算方法、線性有限元方法和等效剛度法的承載力計算結果的對比。圖3為加筋柱殼整體失穩示意圖。

圖3 加筋柱殼整體失穩示意圖

表1 失穩模式及計算時間對比

Tab.1 Comparison of buckling mode and calculation time

有限元方法快速計算方法等效剛度法屈曲載荷/kN146834152572157029軸向屈曲半波數442環向屈曲半波數 9910計算時間/s960.3721.923.18

從表1中可以看出，針對該模型，采用精細有限元模型計算方法耗時960.37s，所得屈曲載荷FFEM=146834kN。從圖3加筋柱殼的失穩波形圖中看出，環向半波數為9，縱向半波數為4。采用本文提出的快速計算方法耗時21.92s，所得屈曲載荷FRBAM=152572kN，相較于精細有限元模型所得承載力的誤差為3.76%，環向半波數為9，縱向半波數為4；采用等效剛度分析方法耗時3.18s，所得屈曲載荷FSSM=157029kN，相較于精細有限元模型所得承載力的誤差為6.94%，環向半波數為10，縱向半波數為2。對比這3種算法可以看出，快速計算方法相較于有限元方法，屈曲載荷計算結果僅差3.76%，計算效率提高了44倍，同時捕捉到軸向屈曲半波數與環向屈曲半波數與有限元結果相同；快速計算方法相較于等效剛度法，計算精度更高，同時捕捉到的軸向屈曲半波數與環向屈曲半波數更精準。

表2給出了快速計算方法與傳統的等效剛度法的網格加筋筒等效剛度陣差異對比，發現兩種方法計算的等效剛度陣中的A陣基本一樣，傳統的等效剛度法的B陣和D陣都比快速計算方法的大。由于傳統的等效剛度法對筋條進行等效，難以準確描述筋條與蒙皮之間的耦合關系，且基于梁假設的解析方法對筋條的彎曲剛度估計過高，導致預測的屈曲載荷普遍偏大。

表2 快速計算方法與傳統的等效剛度法的ABD陣差異分析

為驗證本文提出的整體失穩載荷快速計算方法的適用性，針對加筋柱殼的筋條參數，基于最優拉丁超立方法抽取100個樣本點進行誤差分析。圖4給出了本文提出的快速計算方法與有限元方法的相對誤差，可以看出快速計算方法較有限元方法的相對誤差均在10%以內，其中81個樣本點的相對誤差在5%以內。圖5給出了等效剛度法與有限元方法的相對誤差，可以看出等效剛度法較有限元方法有6個采樣點的相對誤差超過10%，最大誤差可達40.23%，所有樣本點相對誤差均超過5%。上述誤差分析說明，本文提出的快速計算方法預測誤差魯棒性更強，表明了本文提出的快速計算方法的可信性。

圖4 快速計算方法與有限元方法的相對誤差

圖5 等效剛度法與有限元方法的相對誤差

2 大直徑網格加筋柱殼承載力快速優化設計

2.1 基于自適應代理模型的快速優化方法

動態代理模型以其高效的優化效率被廣泛研究發展，其中應用最廣泛的是Jones于1998年提出的高效全局優化方法[17-18](Efficient Global Optimization，EGO)，國內外許多學者發展并改進該算法應用于眾多行業中[19-24]。該方法通過構造一種期望改進(Expected Improvement，EI)的加點準則，考慮了預測模型的均值和方差，尋找該準則的最大期望提高點作為新的樣本點來更新預測模型，反復迭代直至結束。由于EI加點準則同時考慮了模型的預測最優值和誤差性，該算法具有優異的全局和局部搜索能力。面向運載火箭貯箱的初步設計的高效率需求，本文基于加筋柱殼承載力快速計算方法，搭建了基于EGO代理模型的加筋柱殼高效優化設計框架，如圖6所示。

圖6 基于EGO代理模型的快速優化流程

針對兩種網格加筋類型，以筋條高度H、筋條寬度tw、筋條間距bs、蒙皮厚度ts為設計變量，以整體承載力和局部承載力都滿足相應的設計要求為約束條件，開展內壓為0.1MPa與軸壓耦合工況下網格加筋結構輕量化設計，優化列式如下

(16)

式中，Tlj為整體和局部承載的臨界軸力，Tsj為整體承載和局部承載的設計載荷。

2.2 加筋柱殼優化算例

本節采用直徑D為9500mm，長度L為6000mm的兩種加筋典型的柱殼算例，如圖7所示。加筋柱殼承受0.1MPa內壓和軸壓耦合作用。材料選用鋁合金，彈性模量E為68264MPa，泊松比μ為0.33。加筋柱殼可設計的變量分別為蒙皮厚度ts、筋條高度H、橫筋寬度tw1、豎筋寬度tw2、筋條間距bs，優化質量約束為4299.00kg。加筋柱殼極限承載力采用1.1節中提出的快速計算方法，同時考慮局部失穩的情況。兩種加筋類型的優化區間和最優設計參數如表3所示。

(a)正置正交

(b)等邊三角

表3 兩種加筋構型優化設計空間

由表3可知，最優解對應的結構均到達質量約束上限附近，材料得到充分利用。根據參數優化結果可以得出：對于本文的大直徑筒殼模型，在0.1MPa內壓與軸壓耦合工況下，在相同的質量約束以及相同設計區間，等邊三角形網格加筋軸壓承載力優于正置正交網格加筋筒。

兩種加筋類型的加筋柱殼迭代過程如圖8所示，從圖8中可以看出，隨著樣本點的不斷加入，優化最終收斂。

(a)正置正交

(b)等邊三角

3 內壓對大直徑網格加筋柱殼承載力的影響規律

大直徑運載火箭在飛行過程中，貯箱內燃料不斷消耗隨之內壓逐漸降低，工況發生改變，因此需要探究內壓對大直徑網格加筋柱殼承載力的影響規律。由式(6)～式(15)可知，網格加筋柱殼的蒙皮局部失穩載荷和筋條局部失穩載荷與內壓大小無關，因此，本節僅探究內壓對網格加筋柱殼整體失穩載荷的影響規律。

本節以2.2節中的網格參數最優設計為例，采用了4種不同直徑的加筋柱殼算例，直徑D分別為4000mm、5000mm、7500mm和9500mm，長度L均為6000mm。基于快速計算方法獲得的正置正交和等邊三角形網格加筋柱殼結構整體失穩載荷隨內壓的變化規律如圖9所示。

(b)等邊三角

由圖9可知，對于正置正交和等邊三角兩種加筋形式最優結構的加筋筒，內壓為零時，隨筒殼直徑的增大整體失穩載荷基本不變；內壓大于0.1MPa時，筒殼極限承載力隨直徑的增大，基本呈線性增長的趨勢。對于正置正交最優結構的加筋筒殼，直徑為9500mm時，加筋柱殼的整體失穩載荷隨內壓增大，增長的速度越來越慢。對于等邊三角最優結構的加筋筒殼，直徑為7500mm和9500mm時，內壓大于0.7MPa后，隨內壓增大，筒殼整體失穩載荷基本不變。

由式(6)可知

(17)

Plj=2πR×Nlj=

(18)

從式(17)和(18)中可以看出，縱橫失穩波數比n2/m2是影響內壓項的重要參數。本文從失穩模式出發，統計了每個直徑和內壓對應的失穩模式和失穩波數，失穩波形n2/m2變化規律如圖10所示。從圖10中可以看出，對于正置正交和等邊三角加筋形式最優構型，內壓越小，筒殼的n2/m2越大，表明隨著內壓的增加導致加筋柱殼的失穩波形發生轉變，內壓較小時內壓隨筒殼整體失穩載荷的貢獻大，內壓較大時整體失穩載荷的內壓項大幅減小，最終呈現加筋柱殼整體失穩載荷先快速增加后基本不變。對于正置正交最優結構的加筋筒殼，n2/m2均大于1，因此極限承載力隨內壓增加一直增大。對于等邊三角最優構型的加筋筒，當直徑為7500mm和9500mm，內壓大于0.7MPa時，n2/m2是一個小量，整體失穩載荷隨內壓的增加而基本保持不變。

(a)正置正交

(b)等邊三角

4 結論

1)基于NIAH方法發展了考慮內壓的非等厚筋條的加筋柱殼承載力快速計算方法，算例驗證了快速計算方法預測誤差魯棒性更強；

2)考慮加筋筒殼整體失穩、蒙皮局部失穩和筋條局部失穩多種失效模式，開展了給定質量約束的加筋柱殼承載力最大化設計，優化結果表明：在0.1MPa內壓和軸壓耦合荷載作用下，等邊三角形網格加筋承載效率明顯優于正置正交網格加筋筒；

3)內壓對整體失穩載荷的影響程度取決于筒殼縱橫失穩波數的比值n2/m2，當筒殼直徑為7500mm和9500mm時，對于最優構型的正置正交和等邊三角形加筋柱殼，內壓增大導致筒殼失穩波形轉變，進而整體失穩載荷呈現先快速增長后趨于不變的規律。

本文工作可為大直徑運載火箭貯箱的初步設計提供較精確的分析和高效的優化方法。