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

(1.大連理工大學(xué)工程力學(xué)系工業(yè)裝備結(jié)構(gòu)分析國家重點(diǎn)實(shí)驗(yàn)室，大連 116024;2.北京宇航系統(tǒng)工程研究所，北京 100076)

0 引言

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

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

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

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

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

1.1 網(wǎng)格加筋柱殼整體失穩(wěn)載荷的快速計(jì)算方法

針對(duì)軸內(nèi)壓下的網(wǎng)格加筋柱殼結(jié)構(gòu)整體屈曲失穩(wěn)載荷，求解流程如圖1所示，具體過程如下[1]：

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

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

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

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

圖1 基于NIAH的快速計(jì)算方法的流程圖

針對(duì)圓筒殼在均勻軸壓下的穩(wěn)定性問題，學(xué)者們已做了大量的理論研究和試驗(yàn)研。針對(duì)正交各向異性筒殼結(jié)構(gòu)，根據(jù)圓筒殼穩(wěn)定性基本方程組[13]

(1)

K11=A11α2+A66β2

K22=A22β2+A66α2

K12=(A12+A66)αβ

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

(2)

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

Nij=Nmin

(3)

N=

(4)

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

Nij=Nmin

(5)

1.2 網(wǎng)格加筋柱殼局部載荷的快速計(jì)算方法

在工程實(shí)際中，網(wǎng)格加筋柱殼失效模式多種多樣，蒙皮局部失穩(wěn)和筋條局部失穩(wěn)是兩種典型的局部失效模式。本節(jié)基于板殼屈曲理論，發(fā)展了非等厚筋條加筋柱殼的蒙皮局部穩(wěn)定性和筋條局部穩(wěn)定性的承載力理論計(jì)算公式。

1.2.1 蒙皮局部穩(wěn)定性

當(dāng)加筋類型為正置正交時(shí)，根據(jù)彈性力學(xué)四邊簡支板軸壓屈曲理論[14]，蒙皮局部失穩(wěn)時(shí)的單位長度軸力為

(6)

(7)

當(dāng)加筋類型為等邊三角形時(shí)，由文獻(xiàn)[15]可知蒙皮單位長度臨界軸力為

(8)

式中，C1=9.72。

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

(9)

Tlj=2πRNlj

(10)

圖2 等邊三角單胞示意圖

1.2.2 筋條局部穩(wěn)定性

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

(11)

當(dāng)加筋類型為正置正交時(shí)，豎筋屈曲失穩(wěn)的單位長度軸力為

(12)

Tlj=2πRNlj

(13)

當(dāng)加筋類型為等邊三角形時(shí)，考慮到豎筋和斜筋寬度不相同，那么

(14)

因此，加筋柱殼筋條局部失穩(wěn)時(shí)的整體臨界軸力為

Tlj=2πRNlj

(15)

1.3 快速計(jì)算方法精度驗(yàn)證

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

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

圖3 加筋柱殼整體失穩(wěn)示意圖

表1 失穩(wěn)模式及計(jì)算時(shí)間對(duì)比

Tab.1 Comparison of buckling mode and calculation time

有限元方法快速計(jì)算方法等效剛度法屈曲載荷/kN146834152572157029軸向屈曲半波數(shù)442環(huán)向屈曲半波數(shù) 9910計(jì)算時(shí)間/s960.3721.923.18

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

表2給出了快速計(jì)算方法與傳統(tǒng)的等效剛度法的網(wǎng)格加筋筒等效剛度陣差異對(duì)比，發(fā)現(xiàn)兩種方法計(jì)算的等效剛度陣中的A陣基本一樣，傳統(tǒng)的等效剛度法的B陣和D陣都比快速計(jì)算方法的大。由于傳統(tǒng)的等效剛度法對(duì)筋條進(jìn)行等效，難以準(zhǔn)確描述筋條與蒙皮之間的耦合關(guān)系，且基于梁假設(shè)的解析方法對(duì)筋條的彎曲剛度估計(jì)過高，導(dǎo)致預(yù)測的屈曲載荷普遍偏大。

表2 快速計(jì)算方法與傳統(tǒng)的等效剛度法的ABD陣差異分析

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

圖4 快速計(jì)算方法與有限元方法的相對(duì)誤差

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

2 大直徑網(wǎng)格加筋柱殼承載力快速優(yōu)化設(shè)計(jì)

2.1 基于自適應(yīng)代理模型的快速優(yōu)化方法

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

圖6 基于EGO代理模型的快速優(yōu)化流程

針對(duì)兩種網(wǎng)格加筋類型，以筋條高度H、筋條寬度tw、筋條間距bs、蒙皮厚度ts為設(shè)計(jì)變量，以整體承載力和局部承載力都滿足相應(yīng)的設(shè)計(jì)要求為約束條件，開展內(nèi)壓為0.1MPa與軸壓耦合工況下網(wǎng)格加筋結(jié)構(gòu)輕量化設(shè)計(jì)，優(yōu)化列式如下

(16)

式中，Tlj為整體和局部承載的臨界軸力，Tsj為整體承載和局部承載的設(shè)計(jì)載荷。

2.2 加筋柱殼優(yōu)化算例

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

(a)正置正交

(b)等邊三角

表3 兩種加筋構(gòu)型優(yōu)化設(shè)計(jì)空間

由表3可知，最優(yōu)解對(duì)應(yīng)的結(jié)構(gòu)均到達(dá)質(zhì)量約束上限附近，材料得到充分利用。根據(jù)參數(shù)優(yōu)化結(jié)果可以得出：對(duì)于本文的大直徑筒殼模型，在0.1MPa內(nèi)壓與軸壓耦合工況下，在相同的質(zhì)量約束以及相同設(shè)計(jì)區(qū)間，等邊三角形網(wǎng)格加筋軸壓承載力優(yōu)于正置正交網(wǎng)格加筋筒。

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

(a)正置正交

(b)等邊三角

3 內(nèi)壓對(duì)大直徑網(wǎng)格加筋柱殼承載力的影響規(guī)律

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

本節(jié)以2.2節(jié)中的網(wǎng)格參數(shù)最優(yōu)設(shè)計(jì)為例，采用了4種不同直徑的加筋柱殼算例，直徑D分別為4000mm、5000mm、7500mm和9500mm，長度L均為6000mm?；诳焖儆?jì)算方法獲得的正置正交和等邊三角形網(wǎng)格加筋柱殼結(jié)構(gòu)整體失穩(wěn)載荷隨內(nèi)壓的變化規(guī)律如圖9所示。

(b)等邊三角

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

由式(6)可知

(17)

Plj=2πR×Nlj=

(18)

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

(a)正置正交

(b)等邊三角

4 結(jié)論

1)基于NIAH方法發(fā)展了考慮內(nèi)壓的非等厚筋條的加筋柱殼承載力快速計(jì)算方法，算例驗(yàn)證了快速計(jì)算方法預(yù)測誤差魯棒性更強(qiáng)；

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

3)內(nèi)壓對(duì)整體失穩(wěn)載荷的影響程度取決于筒殼縱橫失穩(wěn)波數(shù)的比值n2/m2，當(dāng)筒殼直徑為7500mm和9500mm時(shí)，對(duì)于最優(yōu)構(gòu)型的正置正交和等邊三角形加筋柱殼，內(nèi)壓增大導(dǎo)致筒殼失穩(wěn)波形轉(zhuǎn)變，進(jìn)而整體失穩(wěn)載荷呈現(xiàn)先快速增長后趨于不變的規(guī)律。

本文工作可為大直徑運(yùn)載火箭貯箱的初步設(shè)計(jì)提供較精確的分析和高效的優(yōu)化方法。