基于Ansys子模型法的肘板結構優化

2014-12-05 02:18:50史戰新

艦船科學技術 2014年8期

史戰新

(武漢第二船舶設計研究所，湖北武漢430064)

0 引言

肘板結構是容器結構節點的重要連接形式，外壓容器平面艙壁與耐壓殼體連接處的肘板趾端由于變形不協調，應力集中嚴重，容易產生疲勞裂紋。提高疲勞強度的根本方法是降低熱點處的拉應力，而結構形式是降低應力集中提高疲勞壽命的主導因素。本文基于Ansys 子模型分析法對弧形肘板進行參數化建模，獲得彈塑性應力應變歷程，利用局部應力應變疲勞估算法進行弧形肘板疲勞特性分析，最后采用Matlab 與Ansys 聯合仿真進行遺傳算法優化設計。優化后弧形肘板較初始假定模型重量減輕，疲勞壽命大大提高。

1 Ansys 子模型法

1.1 子模型原理

子模型法基于圣維南原理，即如果實際分布載荷被等效載荷代替后，應力和應變只在載荷施加的位置附近有改變。因此只要子模型切割邊界避開載荷集中及應力集中位置，子模型內部就可以得到較精確的解［1］。Ansys 子模型法先建立局部模型，再切割邊界位移插值，將粗的整體模型的求解結果如位移作為局部子模型的邊界條件進行細化網格求解。

1.2 模型優點［2］如下:

1)減少甚至取消有限元實體模型中所需的復雜傳遞區域;

2)方便用戶在感興趣的關鍵區域就不同的設計進行快速優化;

3)通過粗細網格對比幫助用戶驗證網格劃分是否夠細。

1.3 采用子模型的限制

1)全模型切割邊界插值只能經過實體或殼單元，子模型中可以有其他單元如梁單元作為加強筋或由殼體單元變為實體單元;

2)子模型邊界必須遠離應力集中區域，在子模型計算后應就邊界處的應力水平與整體全模型進行比對。

1.4 采用子模型法的步驟

1)生成并分析較粗糙的整體模型，保留模型db 文件及結果rst 文件;

2)生成子模型并劃分好網格，生成節點;

3)提取保存子模型邊界節點為后綴為.node 文件，默認文件名與子模型一致，命令為“nwrite’;

4)在子模型工作環境下恢復整體模型db 文件，在后處理模塊讀入整體模型結果文件，利用命令“cbdof”形成子模型邊界插值文件，后綴為.cbdo;

5)恢復子模型db 文件，在求解模塊讀入上一步生成的cbdo 文件，并根據子模型范圍內實際位移邊界條件及外載荷加載求解;

6)驗證子模型切割邊界應力分布與整體全模型切割邊界應力分布是否基本一致。

2 假定整體模型

2.1 結構模型

為了真實模擬耐壓平面艙壁的受力狀態，根據圣維南原理，當離開不連續體的距離大于時(R 為耐壓殼體內徑，δ 為殼體厚度)，邊緣及應力集中處的附加應力的影響可以忽略不計。因此建立內徑為9 200 mm，長度為8 300 mm 的假定耐壓圓柱殼體模型，耐壓平面艙壁距離左右邊界分別為2 900 mm和5 400 mm。

由于本文研究肘板結構的局部應力應變響應，故在建立幾何模型時考慮肋骨、各類肘板等詳細的邊角特征，并進行有限元模型板厚偏置。

2.2 邊界條件

在計算壓力下，加筋板殼將會產生大變形，局部結構不連續處會產生大的應力集中，導致構件局部進入塑性狀態，計算壓力下肘板結構力學模型已不適用線彈性克?；舴?勒夫薄板理論，從而應該應用彈塑性力學進行有限元求解。考慮到減少橫艙壁構架承受沿構架長度方向的壓力，形成復雜彎曲梁，同時為減小耐壓殼板與橫艙壁構架的應力集中程度，橫艙壁構架布置于耐壓艙室外。耐壓艙室殼板及艙壁板施加計算壓力4.9 MPa。耐壓殼體耐壓端口邊界施加全約束，采用mm，N，MPa，t 單位制，如圖2所示。

圖2 有限元模型Fig.2 The FEM model

2.3 模型單元選擇及計算設置

Shell281 是8 節點2 階可導形函數單元，具有彈塑性、大變形、大應變等非線性分析功能。

本文計算中打開大變形開關，采用稀疏矩陣求解器及完全NR(NROPT，FULL，ON)迭代，使用漸進式加載(KBC=0)，打開自動時間步長及線性搜索，最小子步為10，最大子步為100。收斂準則采用L2 范數即通過檢查所有自由度方向上的不平衡力(或力矩)的平方和的平方根指標進行收斂檢查，對于工程計算一般誤差小于5%就可以接受。本文設置力的L2 范數收斂準則為1%(默認為0.001)，以加快收斂。

2.4 整體模型粗網格計算結果

采用100 mm 四邊形網格劃分進行全模型進行大變形材料非線性分析，結果如圖3所示。

圖3 粗網格整體模型結果Fig.3 The coarse grid model results

3 子模型分析

3.1 三角形肘板子模型

由于子模型的切割邊界要遠離應立集中區域，邊界的選取也要經過嘗試驗證后才能確定，根據整體模型應力分布初步選取Z 方向范圍(1 300 ～8 300)、X 方向范圍(1 325 ～4 600)、Y 方向范圍(350 ～3 100)，其在全局坐標系下的位置與全模型中一致，如圖4所示。

圖4 子模型幾何模型Fig.4 Sub-model geometry

3.2 三角形肘板子模型網格劃分

采用25 mm 方形映射網格劃分，如圖5所示。

圖5 子模型10 mm 映射網格Fig.5 10 mm mapped meshing of sub-model

3.3 三角形肘板子模型邊界位移插值

子模型加載插值邊界及部分殼板、橫艙壁水壓力載荷邊界條件后，如圖6所示。

圖6 子模型邊界插值及加載Fig.6 Boundary interpolation and loading of sub-model

3.4 三角形肘板子模型邊界位移插值

求解后應力分布如圖7所示。

圖7 子模型未平均表面mises 應力Fig.7 Not average surface mises and first principal stress distribution of sub-model

驗證邊界插值的準確性，進行路徑插值分析比對，如圖8 ～圖9所示。根據對比，在邊界上子模型與全模型的應力分布與大小基本一致，可以認為子模型邊界選取基本合適。

圖8 切割邊界路徑1 ～3 應力比較Fig.8 Stress comparison of cutting edge path 1 ～3

圖9 切割邊界路徑4 ～8 應力比較Fig.9 Stress comparison of cutting edge path 4 ～8

盡管采用高階殼單元彈塑性分析一定程度上可以避免幾何奇異造成的應力集中，但應力集中無法消除。應力集中主要在3 個位置:①肘板端部焊趾及焊趾周圍殼板，應力集中最為嚴重;②各肋骨與橫梁相交處;③橫梁角隅。如圖10所示。

3 處應力集中是一個整體，相互影響，相互作用。其中位置1 主要由肘板與殼板剛度決定，位置2 主要由橫梁整體(含肘板)與殼板肋骨的剛度協調性決定，位置3 主要由橫梁角隅設計形式決定。一般解決位置3 處應力集中可以在橫梁角隅根部開應力釋放孔，圓孔半徑一般不超過橫梁高度的1/4。而通過改善肘板結構的形狀降低剛度劇烈差異來形成各結構位移協調，是改善位置1和位置2 處應力集中的關鍵。

圖10 應力集中位置未平均表面應力分布Fig.10 Not average surface stress distribution of stress concentration position

3.5 弧形肘板子模型

鑒于弧形肘板已成為各船級社共同推薦的結構形式，但在具體細節如外載環境、焊趾高度、弧度大小與橫梁關系等方面有差別。經上文計算可知子模型邊界劃分合理，肘板應力集中位置遠離子模型邊界，改變肘板局部形狀不影響子模型邊界的位移應力分布，故采用如下改進方案:

采用肘板寬高等長構造弧形腹板，面板依腹板變化為矩形曲面，此類型肘板制造方便，肘板剛度較三角形肘板降低，橫梁端部變形與殼板協調性增加。計算結果表明肘板剛度的減小帶來了位置1 處由原來的應力993.75 MPa 降為918.37 MPa，降低7.5%，同時各肋骨線處應力集中明顯下降，如圖11所示。

從弧形肘板彈塑性分析來看，Full 模式顯示下，弧形肘板的最大應力點與最大應變點位置不同，最大應變點的應力較小，最大應力點的應變較小，同時它們的應力應變載荷歷程響應存在著差異，最大應變位置是未平均表面應力最大的位置，其應力在0.7 s 子步達到屈服極限，并開始塑性變形，而最大應力位置點則在0.95 s 子步進入塑性狀態，應力比前者稍大，如圖12 ～圖15所示。應力集中區域在進入塑性后，局部位置產生大應變，總體結構剛度矩陣隨載荷非線性變化。盡管不同位置點應力大小接近，但其應變歷程不同，而且在Power Graphics模式下，最大應變點與最大應力點重合，也表明弧形肘板結構應力集中位置先進入塑性，大應變產生后，結構應力分布根據新的剛度矩陣發生變化，應變最大點塑性區域擴展至內部及其他附近區域。

圖11 4.9 MPa 弧形腹板未平均表面應力分布圖Fig.11 Not average surface stress distribution of curved bracket at 4.9 MPa

圖12 最大應力點應力應變歷程曲線Fig.12 Stress-strain process curve of biggest stress points

圖13 最大應力點應力隨時間變化曲線Fig.13 Stress-time process curve of biggest stress points

圖14 最大應變點應力應變歷程曲線Fig.14 Stress-strain process curve of biggest strain points

圖15 最大應變點應力隨迭代時間變化曲線Fig.15 Stress-time process curve of biggest strain points

4 肘板結構疲勞分析

FE-SAFE 軟件疲勞分析中雙軸應變疲勞基于局部應力應變方法，主要用于低周疲勞壽命估算，其中Brown Miller 算法適合絕大多數金屬材料，對于韌性金屬給出比較精確的結果，對脆性金屬產生非保守值，使用彈塑性應力時，可提供多軸彈塑性修正。

Brown Miller 公式如下［3］:

平均應力修正采用適合低周疲勞修正的Morrow 法:

根據假定材料疲勞特性參數，利用曼森-柯芬應變疲勞壽命公式等，就可以獲得穩態應力-應變循環曲線、應變疲勞壽命曲線，同時考慮工程上鋼板表面一般粗糙度Ra=12.5，故軟件表面系數設為4 ～16μm 等級，材料缺陷、加載頻率、海水介質等影響因素已在經過試驗得出的材料應變疲勞曲線中予以考慮。

將4.9 MPa 下的結構彈塑性應力應變數據輸入Fesafe 軟件，采用Morrow 平均應力修正方法，考慮殘余應力［4］0.68σs=558 MPa，進行0 ～500 m 彈塑性循環加載如圖16所示，局部應力應變法獲得三角形肘板疲勞結果如圖17所示，圖中顯示疲勞壽命最小點肘板端部應力集中處，為101.948=88 次。

圖16 三角形肘板彈塑性疲勞壽命分析設置Fig.16 Elastic-plastic fatigue life analysis settings of triangle bracket

圖17 三角形肘板彈塑性疲勞壽命分析Fig.17 Elastic-plastic fatigue life analysis of triangle bracket

圖18 弧形肘板彈塑性疲勞壽命分析Fig.18 Elastic-plastic fatigue life analysis of curved bracket

同理計算3.5 節弧形肘板彈塑性應力應變響應，局部應力應變法獲得弧形肘板疲勞結果如圖18所示。圖18 顯示，疲勞壽命最小點位于肘板端部應力集中處，最小循環此數在肘板端部為102.142=138次，比三角形肘板壽命提高了56.8%。

5 肘板結構優化

5.1 肘板參數化子模型

本文就弧形肘板利用Ansys 自動導圓功能在肘板腹板垂直邊之間構造弧線過渡，由于肘板兩邊V(I，5)，V(I，6)垂直，故隨著V(I，1)，V(I，2)，V(I，5)，V(I，6)，V(I，7)和V(I，8)變化，弧形是圓的一部分或者橢圓的一部分。其中肘板尺寸參數數組變量設置如圖19所示。

1)變量設置

變量范圍為:v(I，1)∈［400 1500］;v(I，2)∈［400 1800］;v(I，3)∈［50 80］;v(I，4)∈［50 80］;v(I，5)∈［5 100］;v(I，6)∈［5 300］;v(I，7)∈［10 30］;v(I，8)∈［10 30］;v(I，9)∈［60 80］;v(I，10)∈［80 90］;v(I，11)∈［100 180］。

圖19 子模型中局部肘板參數設置Fig.19 Local bracket parameter settings of sub-model

2)Ansys 子模型幾何參數篩選優化

為加快收斂及實現對種群個體的一個迭代子步內循環求解，需要對幾何參數變量進行篩選優化和迭代計算。若在Matlab 中對參數變量進行篩選，則會造成種群數組的不匹配，不利于進行交叉變異。因此在Ansys 中采用循環結構對Matlab 輸出參數進行求解，不合理參數的目標值直接賦予大數以便在Matlab 中排除，符合的繼續進行求解。

3)子模型目標參數輸出選擇

提高壓力容器疲勞強度的根本方法是降低熱點處的拉應力，即Anysys 應力結果的第一主應力S1，而S1 與Mises 應力正相關，且彈塑性材料的屈服判斷準則為Mises 準則，故本文采用Mises 應力反應疲勞壽命的高低。本文采用子模型質量和子模型最大平均應力點應力進行雙目標優化。

5.2 Matlab 優化控制程序

基于Sheffield 大學Gatbs 遺傳算法工具箱進行編程設計如下:

1)設計變量

將子模型中肘板的11 個尺寸參數作為設計變量，在Matlab 中設置變量范圍區間，初始二進制隨機種群，初始十進制隨機變量種群。

2)約束條件

為了保證優化后模型處于彈性應變狀態，假定材料屈服極限σs為820 MPa，對應應變為0.004，由于模型進入彈塑性后應力與應變不再是線性關系，因此子模型優化時約束最大應變小于0.004 較為合理，另外約束最大剪應力應小于0.57σs=467 MPa。

3)目標函數建立

將Ansys 輸出的肘板體積乘以密度轉換重量值作為目標函數1，將Ansys 中輸出的最大應力值作為目標函數2。根據約束條件加入懲罰因子對目標函數1 進行篩選，加速優化速度。

4)交叉率和變異率函數建立

根據遺傳算法的特點在迭代后期，種群個體趨向最優解，為了保持優良個體的穩定，應根據迭代次數采用較小的交叉率和變異率。log(pmmin/pmmax)為最大與最小變異率比值的對數解，保證在最后一次迭代變異率為pmin。

5)種群處理技巧

對于雙目標優化，種群個體要分為兩部分，分別輸入2 個目標函數，并根據適應度函數進行排序按照預定代溝進行優選。當種群個體為奇數時，分流會出現例如chrom(1:4.5)的情況，A 中個體數目其實等效于chrom(1:4)中數組數目，而在進行目標函數優選時，select 函數將變量個體數目乘以代溝后，采取四舍五入的取值規則，這樣就造成了目標值與變量值的數量不統一。為了實現無論種群個體奇數、偶數都能順利實現分流、優選、重組、交叉變異，采取種群個體數量先行取整分組。

6)遺傳迭代處理技巧

為了加快迭代效率，減少迭代數量，本文采用一次迭代循環內進行2 次(父代、子代)目標值計算，利用reins 函數將子代目標值計算的優選個體代替父代目標值計算選用種群中的不良個體，形成新的種群，如此循環得到最優解。

5.3 優化結果

初始種群包含10 個個體，迭代次數為30 次，初始種群變異率為0.1，最大變異率為0.1，最小變異率為0.01，初始種群交叉率為1，最大交叉率為0.9，最小交叉率0.5，遺傳算法迭代跟蹤如圖20 ～圖22所示，變量參數結果如表1所示。