離子推力器柵極組件有限元的建模及熱變形研究

梁秀強，袁杰紅，周仕明

（國防科技大學 空天科學學院，長沙 410073）

0 引言

由于柵極組件在工作狀態(tài)下溫度較高[1]，當結(jié)構(gòu)熱變形嚴重時，會導致雙柵間短路或是柵極聚焦性能變差，從而影響推力器的整體工作性能及可靠性和壽命。

當前對離子推力器柵極熱變形的探究主要分為實驗法、理論法和數(shù)值仿真法。實驗方面主要是NASA開展了相關實驗，實驗系統(tǒng)從接觸式機械測量系統(tǒng)[2]發(fā)展到非接觸式光學測量系統(tǒng)[3]。理論法和數(shù)值仿真法目前所做的工作較少，主要是將柵極等效為無孔平板模型[4-6]，或者建立孔數(shù)較少的簡化模型進行分析[7]，計算其熱變形及熱應力。

由于柵極組件為多孔薄壁扁球殼結(jié)構(gòu)（30 cm柵極約有17 000多孔），為其建模及網(wǎng)格劃分帶來了巨大困難，而且全尺寸有限元建模及計算費時費力，不便于工程應用，因而前人的工作多集中于使用理論推導的柵極等效模型以及簡化的有限元模型來分析其熱變形。但是，如此簡化分析的方法弱化了柵極組件最重要的多孔特征，不符合柵極實際結(jié)構(gòu)，由此計算得到的熱變形結(jié)果可信度不高。

基于此，為高效、精確地計算柵極熱變形及熱應力，提出了柵極組件的全尺寸有限元建模方法，并分別使用殼單元、體單元和梁單元建立有限元模型，對比三種單元的計算結(jié)果并分析熱變形與柵極開孔結(jié)構(gòu)的關系，為離子推力器工程化研制提供技術支撐。

1 有限元建模

對于柵極殼單元及體單元模型，將其分為兩個區(qū)域：一是內(nèi)部開孔區(qū)域；二是邊緣未開孔區(qū)域。對于內(nèi)部開孔區(qū)域，將其視為含孔的六邊形胞元組合而成。在平面上建立含孔的正六邊形面，使用四邊形網(wǎng)格對其進行網(wǎng)格劃分，然后將其沿與水平夾角30°方向陣列，得到若干個網(wǎng)格模型，然后將這些網(wǎng)格模型分別沿水平方向以及與水平夾角60°方向陣列，控制陣列數(shù)量，使得開孔不會超出柵極邊界。對于邊緣未開孔區(qū)域，建立一條由外圍含孔正六邊形的外部邊組成的邊緣線，以此為工具，切割出邊緣未開孔區(qū)域的幾何面，并在其上開孔作為螺釘固定孔，對其劃分網(wǎng)格，便得到邊緣未開孔區(qū)域的平面網(wǎng)格模型，如圖1所示。

圖1 柵極組件有限元建模圖Fig.1 Method of creating grids assembly FEM model

將上述兩部分網(wǎng)格模型繞平面法線進行旋轉(zhuǎn)陣列，然后將所有節(jié)點向母線旋轉(zhuǎn)得到的柵極球殼面上垂直投影，并賦予殼單元屬性，得到柵極殼單元模型，進一步通過殼單元法向拉伸單元可建立柵極體單元模型。采用類似的方法，可以建立柵極梁單元模型，不同之處在于其內(nèi)部含孔區(qū)域是將六邊形邊劃分為分段等截面矩形梁單元。另外在完成手動建模的基礎上，利用MSC.Patran的二次開發(fā)工具PCL語言實現(xiàn)參數(shù)化建模。如圖2所示為三種單元類型模型示意圖。

2 熱變形計算與分析

2.1 溫度場和邊界約束條件

采用文獻[6]中的方法，由于溫度分布T(t,r)與徑向坐標r（方向為柵極中心溫度最高處至邊緣處，范圍為0～0.15 m）之間的對稱性，根據(jù)一維熱傳導問題的基本解來進行溫度預估，并在保證精度的基礎上，將其擬合為二次多項式T(r)=Ar2+Br+C形式，施加到有限元模型。

柵極組件是通過邊緣未開孔區(qū)域沿周向均布分布的24個螺釘固定孔，用螺釘固定在推力器上的，所以選取孔周圍的節(jié)點作為約束點，限制其6個自由度，以此作為邊界約束條件。

圖2 三種單元類型模型示意圖Fig.2 Three types of FEM models

2.2 三種單元模型計算結(jié)果對比分析

建立屏柵三種單元模型，施加溫度場，定義邊界約束條件，計算結(jié)果顯示三者具有相似的熱變形及熱應力分布規(guī)律，如圖3所示為殼單元模型Z方向（垂直于柵極投影平面方向，下同）位移及Mises應力分布云圖。屏柵三種單元模型Z方向位移沿徑向分布曲線如圖4所示。

從圖3、圖4可以看出，三種模型計算得到的熱變形分布規(guī)律基本一致，呈中心對稱分布，在邊緣未開孔區(qū)域吻合較好，并且最大位移均出現(xiàn)在徑向坐標約為100 mm處；熱變形在徑向坐標為100～170 mm區(qū)域內(nèi)變化較為劇烈，而徑向坐標小于100 mm區(qū)域內(nèi)熱變形相差較小；邊緣未開孔區(qū)域的Mises應力遠大于其他地方，在螺釘開孔處會產(chǎn)生應力集中。另外，從計算精度來看，殼單元和體單元模型的計算結(jié)果基本一致，梁單元模型結(jié)果偏大。從計算效率來看，梁單元模型計算時間較少，殼單元模型次之，體單元模型計算效率最低，甚至當其單元規(guī)模超過70萬時，將導致數(shù)據(jù)文件過大而無法計算。

對于梁單元模型計算結(jié)果偏大的原因，可通過理論分析加以解釋，由材料力學知識可知單位寬度矩形截面梁抗彎剛度的計算如式（1）：

圖3 殼單元模型Z方向位移及Mises應力分布云圖Fig.3 Z-directional defomation and stress distribution of shell finite element models

圖4 屏柵三種單元模型Z方向位移沿徑向分布曲線Fig.4 Z-directional defomation distribution curve of three finite element models along r axes

而根據(jù)薄板彎曲理論，單位寬度薄板抗彎剛度計算如式（2）：

比較兩式可知，由于計算因子1/1-μ2大于1，因此在厚度相同的情況下，梁單元模型的整體剛度偏小，同一溫度場下會產(chǎn)生較大熱變形。

2.3 開孔率及開孔方式對柵極熱變形的影響

柵極開孔率由孔徑和孔間距決定。首先保持孔間距不變，通過改變孔徑大小調(diào)節(jié)柵極開孔率，不同開孔率柵極殼單元模型Z方向位移沿徑向分布曲線以及中心Z方向位移如圖5所示。

圖5 不同開孔率柵極殼單元模型Z方向位移沿徑向分布曲線及中心Z方向位移曲線Fig.5 Z-directional defomation distribution curvealongr axes and center point deformation of different aperture ratio model

從圖5可以看出，在相同溫度條件下，隨著開孔率增大柵極熱變形呈現(xiàn)出減小的趨勢。不同開孔率下，柵極邊緣Mises應力云圖如圖6所示，可以看出，隨著開孔率增大，柵極邊緣約束處附近應力值以及應力梯度相應減小，從而導致柵極中心熱變形減小，此規(guī)律與材料力學中簡支梁中點相應位移變化下的支點附近應力分布變化情況一致。另一方面，結(jié)構(gòu)產(chǎn)生的熱應力與材料熱膨脹系數(shù)、彈性模量以及溫度差有關，溫度場不變的情況下，增大開孔率時，熱應力變化很小，而熱應力所合成的內(nèi)力則會減小，因而約束處附近的應力值及應力梯度也相應減小。

另一方面，在開孔率不變的情況下，可以選擇“開大孔”或是“開小孔”。“開大孔”即是選取較大的孔間距，由于柵極口徑（針對30 cm口徑柵極）是一定的，較大孔間距將會導致開孔數(shù)量減小，因此要保證開孔率不變，需要選取較大的孔徑，“開小孔”則與之相反。為研究“開大孔”及“開小孔”對柵極熱變形的影響，建立了如表1所列三種開孔方式下的柵極殼單元模型，對比計算結(jié)果。

圖6 不同開孔率下柵極邊緣Mises應力云圖Fig.6 Thermal stress on edge of different aperture ratio grids

表1 三種開孔方式Table1 Three aperture methods

三種開孔方式下的殼單元模型Z方向位移沿徑向分布曲線如圖7所示。可以看出，在開孔率幾乎不變的情況下，“開大孔”或是“開小孔”對柵極的熱變形影響非常小。在開孔率不變的情況下，不同開孔方式的柵極整體剛度幾乎不變，因而在相同溫度場作用下，熱變形基本一致。

圖7 三種開孔方式下殼單元模型Z方向位移沿徑向分布曲線Fig.7 Z-directional defomation distribution curve along r axes of three aperture methods

3 結(jié)論

建立全尺寸柵極殼單元、體單元和梁單元模型，在此基礎上，計算比較了三種單元模型熱變形及熱應力結(jié)果，研究了開孔率及開孔方式對柵極熱變形的影響。

由于目前實驗條件尚未成熟，無法得到柵極工作時的真實溫度場以及熱變形位移數(shù)據(jù)，后續(xù)將開展實驗，通過人為加熱的方式，測量柵極組件溫度場和熱變形數(shù)據(jù)。進一步工作便是在已有的柵極模型的基礎上，根據(jù)實驗條件開展熱傳導和熱變形仿真分析，對比實驗結(jié)果驗證模型準確性。