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

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

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

1 研究問題描述

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

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

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

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

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

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

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

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

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

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

限制沙漏模式的方法有:

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

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

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

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

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

2 有限元模型

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

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

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

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

3 結(jié)果分析

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

圖4 完全積分單元的應(yīng)力云圖

圖5 減縮積分單元的應(yīng)力云圖

圖6 非協(xié)調(diào)單元的應(yīng)力云圖

由圖4 至圖6 可以看出:

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

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

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

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

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

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

圖7 葉片出汽邊的52 個應(yīng)力和位移量考察點(diǎn)

圖8 沿葉片出汽邊的應(yīng)力分布圖

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

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

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

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

［6］郭建生.有限元方法中單元類型和密度對計算結(jié)果的影響［J］.湖北工業(yè)大學(xué)學(xué)報，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.