談尚炯，程 凱，顧璐寅，劉 巖

(1.上海汽輪機廠有限公司，上海 200240;2.上海電氣電站設備有限公司，上海 200240)

在使用商業(yè)有限元軟件進行長葉片應力計算的過程中，網(wǎng)格劃分是重要的組成部分，網(wǎng)格的質量與單元類型直接影響到計算結果的精確性。在長葉片應力計算的前處理中，已經(jīng)形成了較為規(guī)范化的網(wǎng)格劃分方法，即盡量采用線性八節(jié)點六面體單元對葉片大部分規(guī)則區(qū)域劃分結構化網(wǎng)格，在非規(guī)則區(qū)域輔以少量四面體與五面體非結構網(wǎng)格，以保證網(wǎng)格的合理過渡與整體質量。但是，在實際的長葉片有限元計算中發(fā)現(xiàn)，不同的網(wǎng)格單元類型也會導致計算結果的差異，尤其是分別使用完全積分單元、減縮積分單元與非協(xié)調單元時，應力差異最大可達30%以上，甚至跨越屈服極限。這說明，這三種單元類型在計算所得的應力中，存在不符合實際情況的結果。所以，對這三種網(wǎng)格類型進行深入分析，找出究竟何種網(wǎng)格類型適用于長葉片的應力計算，這對提高長葉片應力結果的精確性有較大的工程意義。

1 研究問題描述

完全積分單元，是指當單元具有規(guī)則形狀時，所用的高斯積分點的數(shù)目足以對單元剛度矩陣中的多項式進行精確積分［1］。但是，當承受彎曲載荷時，線性完全積分單元會出現(xiàn)剪切閉鎖問題，造成單元過于剛硬。剪切閉鎖的原理如圖1 所示。發(fā)生剪切閉鎖時，即使網(wǎng)格劃分很細，計算精度仍然很差。

圖1 完全積分單元在彎矩作用下的變形示意圖［3］

相反，減縮積分單元可以避免剪切閉鎖問題的發(fā)生。減縮積分單元是在普通的完全積分單元的基礎上在每個方向減少一個積分點，這種單元形式大大減小了剛度矩陣的規(guī)模，節(jié)約了計算時間。但是由于減縮積分單元的積分點位于單元的中心，單元的變形不會引起積分點處應變能的產生，實際出現(xiàn)的單元變形模式像沙漏的形狀，稱為沙漏模式［1］。沙漏模式的原理見圖2，沙漏模式下的結構表現(xiàn)得非常柔軟，在應力結果上與真實情況相差很大，并且會導致不真實低頻模態(tài)的出現(xiàn)。

圖2 減縮積分單元在彎矩作用下的變形示意圖［3］

非協(xié)調單元也可以糾正完全積分單元出現(xiàn)的剪切閉鎖問題。非協(xié)調單元將增強單元變形梯度的附加自由度引入到線性單元中，以避免單元交界處的位移場出現(xiàn)重疊或裂痕。但是如果單元形狀較差，非協(xié)調單元的分析精度就會下降。如果所關心部位的單元扭曲較大，尤其是出現(xiàn)交錯扭曲，分析精度就會受到嚴重影響。

綜上所述，剪切閉鎖與沙漏模式是一對矛盾，想同時消除剪切閉鎖與沙漏模式是不可能的。上述三種單元類型沒有絕對的優(yōu)劣之分，我們要做的只是針對不同的應力狀態(tài)選擇合適的單元類型，既要避免大面積的剪切閉鎖問題，也要限制沙漏模式的發(fā)生。

文獻［3-7］中緩解與克服剪切閉鎖問題的方法有:

1)使用減縮積分單元，但是減縮積分單元會導致沙漏模式的產生。

2)使用二階完全積分單元，但是使用這種單元無法進行接觸分析。對于長葉片結構，計算模型涉及葉根與葉根槽的接觸，以及圍帶的接觸等。而且，使用二階單元，會大大增加計算量，從而提高葉片設計開發(fā)的成本。

3)使用細化的網(wǎng)格。一方面，這只能緩解剪切閉鎖問題的發(fā)生，而不能從根本上解決問題;另一方面，在長葉片的有限元計算過程中，為了控制計算規(guī)模，網(wǎng)格數(shù)量會受到限制，因而單元尺寸不可能任意縮小。

限制沙漏模式的方法有:

1)使用完全積分單元，但是完全積分網(wǎng)格會導致剪切閉鎖問題的產生。

2)使用二階減縮積分單元，同樣，這種單元類型也無法用于網(wǎng)格接觸分析，顯然不適用于長葉片的分析過程。

3)使用細化的網(wǎng)格。同樣，這無法從根本上解決問題，且為了控制計算規(guī)模，這種方法可行性不強。

另外，文獻中還提到可以通過添加自由度或添加人工剛度等方式來抑制沙漏模式的出現(xiàn)。但是這不是我們目前階段的研究重點。我們的目標是針對長葉片的結構與實際受力狀態(tài)，找出適合長葉片有限元計算的單元類型。

為了研究上述問題，我們選用變截面扭轉的某汽輪機長葉片為原型，并進行了簡化，以排除網(wǎng)格質量對計算結果的影響。在同一套網(wǎng)格下使用不同的單元類型，計算該葉片在離心力載荷與汽流力作用下的應力分布。

2 有限元模型

本文采用某產品葉片為母型，建立簡化模型進行分析。為了保留葉片型線的扭轉與漸縮對葉身內部應力分布的作用，并排除網(wǎng)格質量對計算結果的影響，將葉片簡化為如圖3(a)所示的扭轉減縮體。圖中，葉型部分高度為2 000 mm，僅取母型葉片的底部截面與頂部截面的外切矩形，其余部分由幾何線性插值獲得。簡化后，頂部截面與底部截面的面積比為0.165，主慣性軸夾角為58.5°。

葉根部分由中間體截面向下拉伸而成。簡化模型忽略了原有葉型的圍帶部分，在后續(xù)的計算中，只考慮單個自由葉片在運行過程中的應力分布。網(wǎng)格劃分時，對整個葉片全部使用六面體網(wǎng)格，節(jié)點數(shù)為11 915，單元數(shù)為9 528，如圖3(b)所示。網(wǎng)格劃分完成后，對網(wǎng)格單元分別賦予一階完全積分單元、一階減縮積分單元與一階非協(xié)調單元的單元類型。

本文采用非線性有限元軟件ABAQUS 對該葉片的穩(wěn)態(tài)應力進行分析求解。為了更精確地模擬葉片在工作過程中受到的彎力，在有限元模型中，除1 500 r/min 的離心力載荷外，還將CFD 中計算得到的汽流力等效加載到葉身表面。

圖3 葉片簡化模型幾何與網(wǎng)格示意圖

3 結果分析

首先考察葉片整體的應力分布，圖4 至圖6分別給出完全積分單元、減縮積分單元和非協(xié)調單元三種情況下的葉片應力分布云圖。

圖4 完全積分單元的應力云圖

圖5 減縮積分單元的應力云圖

圖6 非協(xié)調單元的應力云圖

由圖4 至圖6 可以看出:

1)三種網(wǎng)格單元類型下計算所得的應力場分布趨勢基本相似，無論選擇哪一種單元進行計算，所得到的葉片受力定性規(guī)律是相同的;

2)應力在葉型的根部表現(xiàn)出相對較大的應力集中，三種網(wǎng)格單元類型的計算結果都顯示應力最大值出現(xiàn)在葉型底部出汽邊位置;

3)完全積分單元所得的應力極值最大，達到749 MPa;減縮積分單元所得的應力極值最小，為644 MPa，僅為完全積分單元計算結果的86%;非協(xié)調單元所得的最大應力值為672 MPa，介于兩者之間。

沿該葉片的出汽邊取52 個點(參見圖7)，輸出其應力值與位移值，如圖8 與圖9 所示。首先，從葉片出汽邊沿葉高方向的應力分布來看，沿著葉高應力呈逐漸減小的趨勢，這是因為沿著葉高的每個截面承受的離心力載荷逐漸減小;其次，變形量的主要成分是離心力導致的扭轉恢復角，變形量沿著葉高方向逐漸增大的趨勢也符合葉身扭轉的真實情況，當然由于該簡化模型是對原有模型的定性描述，與真實葉片的計算結果會有一定差距。

對比三種網(wǎng)格單元類型下的應力計算結果，可以發(fā)現(xiàn):

1)三種單元類型的計算結果表現(xiàn)出一致的趨勢，總體來說，完全積分單元的應力結果較大，而減縮積分單元與非協(xié)調單元的計算結果較小。

圖7 葉片出汽邊的52 個應力和位移量考察點

圖8 沿葉片出汽邊的應力分布圖

圖9 沿葉片出汽邊的位移量分布圖

2)靠近葉型根部與頂部位置，三種單元類型的計算結果表現(xiàn)出一定的差異，其中，在根部附近位置，完全積分單元計算得到的應力最高，而非協(xié)調單元計算得到的應力最低;在葉型頂部位置，完全積分單元計算得到的應力最高，而減縮積分單元計算得到的應力最低。

對比三種網(wǎng)格單元類型下的位移計算結果，可以發(fā)現(xiàn):

1)三種單元類型的計算結果表現(xiàn)了一致的趨勢，沿著葉高方向變形量逐漸增大;

2)沿著葉高方向，完全積分單元計算得到的位移量始終小于減縮積分;

3)非協(xié)調單元計算得到的位移量呈現(xiàn)出與其他兩種單元的計算結果比較大的差距。

根據(jù)三種網(wǎng)格單元類型的特性，對可能導致產生上述現(xiàn)象的原因作出如下推測:

1)減縮積分單元的計算結果出現(xiàn)了一定程度的沙漏現(xiàn)象。查看計算過程中偽應變能與內能，發(fā)現(xiàn)偽應變能僅占內能的0.23%。根據(jù)文獻［2］，當偽應變能占內能的比例低于1%時，沙漏模式對計算結果的影響不大。由此可見，本文算例中的沙漏模式對計算結果的影響很小。

2)完全積分單元的計算結果出現(xiàn)了一定程度的剪切閉鎖現(xiàn)象。由于我們在本文模型的計算過程中考慮了汽流彎力的作用，而完全積分單元的位移場不能模擬由于彎曲而引起的剪切變形與彎曲變形。另外，相比于減縮積分單元的位移計算結果，完全積分單元的結果偏大，這正是剪切閉鎖可能導致的現(xiàn)象。

3)非協(xié)調單元的位移計算結果與其他兩種單元類型的結果有較大差異。推測原因為:對于本文中的模型，葉型受到離心力后的扭轉恢復角較大，以至于產生了網(wǎng)格的平行扭曲，這會大大降低非協(xié)調單元的計算精度。

4 結論

對于汽輪機長葉片而言，影響有限元計算結果準確性的因素很多，本文從單元類型的角度進行分析，分別使用一階線性的完全積分單元、非協(xié)調單元以及減縮積分單元對簡化的葉片模型進行結構分析。分析結果顯示，減縮積分單元能較好地分析長葉片的受力狀態(tài)，在網(wǎng)格質量保證的情況下，沙漏模式的影響非常小;完全積分單元對于受彎力作用的長葉片可能出現(xiàn)比較明顯的剪切閉鎖現(xiàn)象;而非協(xié)調單元對網(wǎng)格質量比較敏感，在網(wǎng)格出現(xiàn)平行扭曲或交錯扭曲的情況下，計算精度會大幅下降。在對本文中的葉片簡化模型進行應力計算的過程中，綜合考慮計算規(guī)模和計算精度，減縮積分單元應當是可選用的適用性較好的單元類型。但對于造型復雜的真實葉片，其應力計算結果受到網(wǎng)格質量與單元類型的共同作用，難以單獨考察單元類型對計算結果的影響，因而如何選擇單元類型有待進一步研究。

［1］石亦平，周玉蓉.ABAQUS 有限元分析實例詳解［M］.北京:機械工業(yè)出版社，2006.

［2］曹金鳳，石亦平.ABAQUS 有限元分析常見問題解答［M］.北京:機械工業(yè)出版社，2009.

［3］包剛強，ERKE Wang，郝清亮，張國兵.對主流有限元軟件控制剪切自鎖和沙漏模式的比較和研究［C］//第八屆中國CAE 工程分析技術年會論文集.成都:第八屆中國CAE 工程分析技術年會，2012:517-524.

［4］張汶艷，祁頂春.淺談ABAQUS 單元選擇［C］//河南省汽車工程學會第五屆科研學術研討會論文集.鄭州:河南省汽車工程學會，2009:174-175.

［5］張利民，王克明，吳志廣.葉片模態(tài)分析的單元類型選擇［J］.沈陽航空航天大學學報，2011，28(2):21-24.

［6］郭建生.有限元方法中單元類型和密度對計算結果的影響［J］.湖北工業(yè)大學學報，2007，22(4):25-28.

［7］DUCHATEAUG，BERGEAUD V，GAIFFE S，SCHEURER B.Control of Hourglass Modes in Prisms with Pentagons or Hexagons As Bases［J］.Finite Element in Analysis and Design，2006(42):1199-1210.