張龍飛，胡全星

（中國航天科技集團公司第四研究院四○一所，陜西西安710025）

固體火箭發(fā)動機試車架中板簧彈阻力有限元計算方法

張龍飛，胡全星

（中國航天科技集團公司第四研究院四○一所，陜西西安710025）

以固體火箭發(fā)動機試車架中重要的零部件板簧為研究對象，為了實現(xiàn)對多段板簧彈阻力的分析計算統(tǒng)一處理，編寫Matlab-m有限元計算程序?qū)υ囼炦^程中板簧的變形產(chǎn)生的彈阻力Q進行分析計算，計算考慮了由軸向力引起的幾何剛度矩陣，得到了彈阻力Q與位移Δ的解析關(guān)系，進而可以確定板簧的最佳工作點。設計的算法、計算思路和結(jié)果為下一步進行板簧設計、結(jié)構(gòu)優(yōu)化、失穩(wěn)分析等提供一定的參考作用。

板簧；Matlab；有限元；剛度矩陣

固體火箭發(fā)動機試車架是地面點火試驗的主要設備，其動架和定架之間靠板簧連接，可使整個結(jié)構(gòu)動態(tài)響應好，對測力傳感器校定影響小，提供沿發(fā)動機軸向運動的小位移自由度，使發(fā)動機推力全部作用到推力傳感器上。但板簧是彈性體，在推力測量過程中會產(chǎn)生阻力，所以有必要對板簧彈阻力進行計算與估計。對此，王德厚等［1］基于力學的基本理論建立微分方程組計算得到了單工作段板簧彈阻力的理論解；文獻［2］建立了雙工作段板簧的有限元模型，通過ANSYS分析得到了具體彈阻力值的大小。為了實現(xiàn)對多段板簧進行更精確的結(jié)構(gòu)分析，計算彈阻力的解析表達式至關(guān)重要，建立微分方程的方法過于復雜，計算過程容易出錯；采用有限元分析軟件（如ANSYS等）只能輸入具體參數(shù)值得到數(shù)值結(jié)果。

筆者建立了3種結(jié)構(gòu)形式板簧的有限元模型，基于有限元的基本理論，通過在Matlab中編寫m文件程序的方法對單段到多段板簧彈阻力的計算實現(xiàn)了統(tǒng)一處理，均可以得到含有板簧幾何參數(shù)變量的解析表達式［3］。首先應用單段板簧的計算結(jié)果和航天標準中的理論解進行對比，用以驗證算法的有效性。隨后分別給定單段到三段板簧的具體參數(shù)進行算例計算，并對3種板簧的彈阻力變化進行對比分析。基于本文的計算結(jié)果，通過結(jié)構(gòu)優(yōu)化可以設計板簧合適的工作點，保證結(jié)構(gòu)穩(wěn)定的條件下提供最小的彈阻力，在試驗中測試到更加精確的推力值。

1 板簧有限元模型

1.1 有限單元法靜力分析理論

結(jié)構(gòu)靜力分析的平衡方程可以表示為

式中：K為剛度矩陣；δ為位移向量；R為載荷向量。

如果為梁單元，考慮軸向力時，則應包含幾何剛度矩陣，即

對于由一系列單元組成的結(jié)構(gòu)，整體剛度矩陣（K+Kg）將總體坐標的結(jié)點位移δ和整個結(jié)構(gòu)的總體載荷R聯(lián)系在一起，總體剛度K、Kg分別是由單元剛度矩陣ke和幾何單元剛度矩陣kg經(jīng)過總裝而成［4］。

1.2 板簧結(jié)構(gòu)以及單元劃分

板簧作為固體火箭發(fā)動機試車架中動架和靜架的連接構(gòu)件，具有重復性好、彈阻力小和結(jié)構(gòu)簡單的優(yōu)點，廣泛應用在各種試車架中［5］。板簧的結(jié)構(gòu)形式如圖1所示。

在圖1中，進行了單元以及幾何參數(shù)的標記，①～⑤表示單元號。三工作段板簧（或多段）為新提出的板簧形態(tài)，相比于單段和雙段板簧，具有可修改參數(shù)多，在大噸位高精度試車架中，提高側(cè)向剛度的同時，為設計彈阻力值及其變化提供更多的選擇范圍。

結(jié)點單元劃分如下：

1）單工作段板簧2個結(jié)點1個單元。

2）雙工作段板簧4個結(jié)點3個單元。

3）三工作段板簧6個結(jié)點5個單元。

1.3 板簧幾何模型簡化

運用Matlab對板簧進行分析計算前，需建立其簡化模型［6］。為了保證計算的準確性，建立板簧簡化模型采取以下規(guī)則：

1）板簧與動架、靜架的連接采用剛體焊接結(jié)構(gòu)，分析中認為焊接質(zhì)量可靠，不考慮焊縫對傳力效果的影響，把板簧當作連續(xù)體來處理。

2）刪除板簧一些微小特征，如螺栓連接的螺孔、螺紋等。其對計算結(jié)果產(chǎn)生的影響可以忽略，完全能保證結(jié)果精確度。

3）應對板簧受到的載荷進行合理的分配，使其作用在適當?shù)奈恢谩?/p>

筆者簡化板簧為梁模型，假設為二維梁單元，結(jié)點具有3個自由度（軸向位移u，橫向位移v，轉(zhuǎn)角θ），根據(jù)板簧實際的安裝狀態(tài)，其約束設置為：

下端結(jié)點固定（u1=0，v1=0，θ1=0）；上端結(jié)點轉(zhuǎn)角約束（θ頂端=0）。

筆者分析所用到的各變量定義如表1所示。

表1 板簧分析中所用到的符號

1.4 板簧受力分析及提出問題

板簧的受力如下圖2所示。

當板簧一端相對另一端產(chǎn)生一個位移Δ時，計算引起位移彈阻力Q的大小，分析Q隨著板簧的軸向壓力P的變化趨勢，進而確定板簧合適的工作狀態(tài)。通過編寫程序，實現(xiàn)對單段到多段板簧的統(tǒng)一計算，可得到含有板簧幾何參數(shù)變量的彈阻力Q的解析表達式。

2 計算結(jié)果及其驗證

2.1 計算結(jié)果

Matlab軟件已成為當今工業(yè)設計計算的主要工具，主要特點是計算元素為矩陣，在處理參數(shù)多、計算量大的問題時快速而精確［7］。程序結(jié)構(gòu)為建立3個函數(shù)。

1）主函數(shù)：該函數(shù)主要負責變量、結(jié)點、單元的定義，原始數(shù)據(jù)的輸入、函數(shù)的調(diào)用、結(jié)果的輸出。

2）單元剛度矩陣計算子函數(shù)：該函數(shù)主要負責接收相關(guān)變量，返回值為單元剛度矩陣和單元幾何剛度矩陣。

3）單元剛度矩陣組裝子函數(shù)：該函數(shù)主要是接收相關(guān)變量，返回值為單元剛度矩陣組裝的參數(shù)值。

幾何單元剛度矩陣和單元剛度矩陣［8］分別具有以下組成形式：

運行程序分別對單段、雙段以及三段板簧進行計算，得到含有板簧結(jié)構(gòu)幾何尺寸、軸向力、橫向力的表達式如下：

式中：P1=25b3E3h16h23；P2=18l3P3r2（3l+2r）+ 10b2E2h13P［3h13r2+h23l（13l+15r）］；P3= 3bEl P2［15h23l2（l+4r）+2h13r2（26l+5r）］；Q1= 3l3P2r2（32l2+39lr+12r2）；Q2=25b2E2h13［h13r3+ 2h23l（4l2+6lr+3r2）］；Q3=10bEl P［h13r2（24l2+ 19lr+3r2）+2h23l2（4l2+18lr+9r2）］

3）三工作段板簧

2.2 算法正確性驗證

由于程序?qū)Χ喽伟寤傻挠嬎憬y(tǒng)一處理，所以只需要對單段進行驗證即可。航天標準通過建立板簧的梁模型，基于力學的基本理論得到板簧撓曲線的微分方程，通過計算給出單工作段板簧的彈阻力Q的理論解如下所示：

對有限元的計算解式（3）進行變量替換得：

對比式（6）、（7）可以看出，只要兩個式子前部分系數(shù)c1、c2相等，就可以說明式（3）計算的結(jié)果是正確的。

從圖3可以看出，兩個系數(shù)值基本吻合，在多段板簧的計算中，由于程序計算是重復迭代進行的，至此可以驗證所編寫程序的正確性。

3 算例

板簧的材料選用彈簧鋼，彈性模量E= 206 GPa。

1）單工作段板簧

參數(shù)取值為：L=240 mm；b=100 mm；h= 4.5 mm。

彈阻力Q的計算結(jié)果為

當P=13 800.85 N，Δ≤0.5 mm時，Q≤33.39 N。

2）雙工作段板簧

參數(shù)取值為：L=240 mm；l=24 mm；r= 192 mm；b=100 mm；h1=2.85 mm；h2=15 mm。

彈阻力Q的計算結(jié)果為

當P=13 800.85 N，Δ≤0.5 mm時，Q≤4.41 N。

3）三工作段板簧

參數(shù)取值為：L=240 mm；l=18 mm；c=12 mm；r=180 mm；h1=2.85 mm；h2=5 mm；h3=15 mm；b=100 mm。

彈阻力Q的計算結(jié)果為

當P=13 800.85 N，Δ≤0.5 mm時，Q≤10.27 N。

4 結(jié)果分析

從彈阻力表達式可以看出，無論是單段、雙段還是多段板簧，當軸向壓力P固定不變時，彈阻力Q均和位移Δ成比例關(guān)系。在試驗過程中，為了提高測試精度，有必要進一步分析Q的變化趨勢以及確定板簧的最佳工作點，根據(jù)Q的表達式，就可以得到當Δ固定不變時彈阻力Q的系數(shù)隨著P的變化趨勢，分別繪制其趨勢如圖4所示。

從圖4可以看出，無論是單段、雙段還是三段板簧，彈阻力的系數(shù)值隨著壓力P是呈直線減小的關(guān)系。針對算例中板簧幾何參數(shù)的取值，計算斜率得ksy=－0.005，kdy=－0.004 44，kty=－0.004 40，即彈阻力隨著壓力P的衰減率是相當?shù)?，對于確定的壓力P，Qsy＞Qty＞Qdy。對于目前廣泛使用的雙段板簧，更改參數(shù)l、r（板簧的總長2l+r=240 mm不變），l、r分別取以下3種模型：

得到彈阻力Q隨著壓力P變化如圖5所示?？梢钥闯?，從模型1到3，l長度減小，r長度增加，彈阻力的值呈增大的趨勢，并隨著壓力P的變化，彈阻力的下降速率相當。

從圖3可以看出單工作段板簧在kL≤1.0的一段，彈阻力變化平穩(wěn)，最有利于用在帶原位校準裝置的試車架中，這時因初始彈阻力已在原位校準中被消除，由彈阻力帶來的誤差主要是因發(fā)動機在工作過程中質(zhì)量逐漸減小所引起的彈阻力變化量所決定的。工作點選在這一段中，保證了彈阻力變化量和推力測量誤差最小。

當板簧受壓力時，P達到臨界載荷Pcr時，彈阻力Q=0，亦即達到板簧的失穩(wěn)狀態(tài)，從Q直線變化的趨勢可以得到，在發(fā)動機試驗過程中，對于特定的發(fā)動機質(zhì)量（壓力P），選定或者設計板簧的時候，在保證板簧不被破壞的基礎上選接近壓板簧Q=0的點為工作點，此時的彈阻力基本上接近于零，用在不帶原位校準裝置的試車架中可減小推力測量誤差。

從算例計算中可以看出，當板簧的位移Δ≤0.5 mm時，彈阻力Q值很小，這樣數(shù)量級的彈阻力大小相比于發(fā)動機的推力值屬于可以接受的范圍。

5 結(jié)論

1）建立了3種結(jié)構(gòu)形式板簧的有限元模型，通過在Matlab中編寫m文件程序計算得到彈阻力Q與位移Δ的表達式。對板簧的靜力分析、失穩(wěn)分析等可以提供更精確的數(shù)學模型。

2）本設計程序具有通用性，將求解多段板簧的彈阻力與位移的關(guān)系統(tǒng)一處理。在將來板簧的設計計算中，對于不同規(guī)格的發(fā)動機選擇的多段板簧，只需要設置初始的輸入?yún)?shù)，不用對程序進行更改，就可以方便地進行計算。

3）板簧的彈阻力跟位移成比例關(guān)系，隨著壓力的增大呈減小趨勢。對于特定的板簧一般選取接近失穩(wěn)的點為工作點，此時彈阻力最小。

（Referenees）

［1］王德厚.單板簧的設計計算［R］.西安：航天動力技術(shù)研究院，1997.WANG Dehou.The design and calculation of the single plate spring［R］.Xi'an：The Research Institute of Space Power Technology，1997.（in Chinese）

［2］潘娜.固體火箭發(fā)動機試車架結(jié)構(gòu)有限元分析［D］.西安：航天動力技術(shù)研究院，2008.PAN Na.The finite element analysis of test frame structure in SRM［D］.Xi'an：The Research Institute of Space Power Technology，2008.（in Chinese）

［3］ROQUE C.Symbolic and numerical analysis of plates in bending using Matlab［J］.Joural of Symbolic Computation，2014，61：3- 11.

［4］RAH MAN T，VALDMAN J.Fast MATLAB assembly of FEM matrices in 2D and 3D：nodal elements［J］.Applied Mathematics and Computation，2013，219（13）：7151- 7158.

［5］薛群，徐向東.固體火箭發(fā)動機測試與試驗技術(shù)［M］.北京：中國宇航出版社，1994.XUE Qun，XU Xiangdong.The test technology of solid rocket motor［M］.Beijing：China Aerospace Press，1994.（in Chinese）

［6］李海濤.火箭發(fā)動機推力矢量測量理論、方法與自動測試技術(shù)研究［D］.長沙：國防科學技術(shù)大學，2005.LI Haitao.Research on the theory＆method of trust vector measurement and the automatic test technology in rocket engine ground test［D］.Changsha：National University of Defense Technology，2005.（in Chinese）

［7］ZANDER N，BOG T，ELHADDAD M，et al.A finite cell research toolbox for MATLAB［J］.Advances in Engineering Software，2014，74：49- 63.

［8］LOGAN D.A first course in the finite element method［M］.Beijing：Publishing House of Electronics Industry，2003.

The Finite Element Method for Caleulating Elastie Resistanee of Plate Spring in SRM Test Frame

ZH ANG Longfei，HU Quanxing
（The 401st Institute of the Fourth Institute of CASC，Xi'an 710025，Shaanxi，China）

plate spring；Matlab；finite element；stiffness matrix

V416.1

A

1673-6524（2015）04-0050-05

2015- 01- 13；

2015- 04- 09

張龍飛（1988－），男，碩士研究生，主要從事固體火箭發(fā)動機試驗測試研究。E-mail：704771327＠qq.com

Abstraet：The important parts（plate spring）in SRM test frame are taken as the research object.For the purpose of achieving a unified treatment in analysis and calculation of elastic resistance of multistage plate spring，the Matlab-m finite element calculation program is to be made for the analysis and calculation of elastic resistance Q produced from the deformation of the plate spring in the process of testing，which takes the geometric stiffness matrix of axial force into account.The analytical relationship of resistance Q and displacementΔis obtained，which makes it possible to pinpoint the optimal working point of plate spring.The designed algorithm，train of thought and results of calculation are preparations for plate spring design，structure optimization，and buckling analysis，etc.