樊一達，毛玉明，舒忠平，王吉飛，張洋洋，于哲峰,*

1.上海交通大學 航空航天學院，上海 200240 2.上海宇航系統(tǒng)工程研究所，上海 201109 3.上海飛機設計研究院，上海 201210

運載火箭在飛行中受到的外載荷包括發(fā)動機的推力、飛行氣動力、陣風載荷等。在計算運載火箭的動力學響應時，往往將火箭簡化為一維梁模型，并將各種外載荷轉(zhuǎn)換為各橫截面上的軸力、彎矩等內(nèi)力，由于這種轉(zhuǎn)化結果比較粗略，在結構設計時需要預留較大的安全裕度。而為了應對不斷增長火箭結構輕量化設計需求，應采用更為精細的三維結構有限元模型計算部件載荷，這對于結構優(yōu)化、減重具有很大意義。此時就需要把氣動載荷數(shù)據(jù)轉(zhuǎn)換三維結構模型上，火箭氣動模型與實際結構外形相同，其表面有許多凸起整流罩，而結構模型往往只建立主要承力結構，其表面形狀較氣動模型而言更為簡單，并且氣動模型與結構模型網(wǎng)格疏密程度也不同。這就產(chǎn)生了氣動載荷轉(zhuǎn)換為結構載荷的問題，即流-固載荷傳遞問題。

在過去的半個世紀中，國內(nèi)外學者針對這種不匹配網(wǎng)格之間的數(shù)據(jù)傳遞問題提出了多種插值方法來解決，如反距離插值法(Inverse Distance Weighting Interpolation, IDW)、Kriging方法和徑向基函數(shù)法(Radial Basis Function，RBS)等。Frank對多種插值方法進行了精度和計算效率的比較，結果表示徑向基函數(shù)法的綜合表現(xiàn)較好。吳宗敏和Buhmann從數(shù)學的角度分析了不同基函數(shù)的徑向基函數(shù)特性。Harder和Duchon等先后提出了無限板樣條(Infinite Plate Spline, IPS)、薄板樣條法(Thin Plate Spline, TPS)，也都是徑向基函數(shù)中的一種，其中TPS法是航空宇航工程中最普遍使用的一種插值方法。

力等效是將一個集中載荷等效分配到若干點上的方法。王專利和徐建新等證明了該方法的可靠性。高尚君等提出了一種基于彈簧-懸臂梁模型最小變形能的氣動載荷分配方法，解決了氣動壓力點與部分結構網(wǎng)格節(jié)點重合而導致的難以求解的問題。

壓力插值法具有局部載荷等效程度高的優(yōu)點，但要求結構模型外形與氣動模型外形保持一致。而力等效法可以保證載荷轉(zhuǎn)換區(qū)域總的合力和合力矩等效，并且當結構模型外形與氣動模型不一致時仍然適用，但局部載荷未必等效。因此，本文針對簡化后、無次要承力部件的結構模型外形與氣動模型外形不一致的問題，提出一種壓力插值與力等效混合使用的流-固載荷轉(zhuǎn)換方法，即結構模型外形與氣動模型外形一致的區(qū)域使用壓力插值法，不一致區(qū)域用力等效法。本文先簡要介紹樣條插值法和力等效法的原理，然后闡述混合法的原理，最后給出某算例驗證的結果。

1 載荷轉(zhuǎn)換原理

1.1 TPS法原理

TPS法是一種三維插值方法，適用于曲面插值，由薄板受分布載荷引起彎曲小變形的模型推導得出，其微分方程為

(1)

式中：為板的抗彎剛度，為板的變形。引入球坐標=sincos,=sinsin,=cos，式(1)的基本解設為

(2)

式中：和是待定系數(shù);為集中載荷。對于作用在(，，)處的個集中載荷，=1，2，…，，板的總變形為基本解式(2)的疊加，可表示為

(3)

(4)

此時，式(4)的+4個未知量(,,,,,,…,)可由力系的平衡方程以及個已知插值點的位移 的方程求出，即

(5)

=1,2,…,

(6)

將式(5)和式(6)寫為矩陣形式，得到含+4個未知數(shù)的線性方程為

(7)

求解式(7)，可求出,,,,,…,,同時把待插值的坐標一起代入式(4)可得待插值點的解。

1.2 力等效法原理

力等效法是將一個集中載荷等效分配到若干點上的方法。其基本規(guī)則是離氣動載荷點近的節(jié)點分配到的載荷多，遠的分配到的少，且保證分配前后的合力、合力矩不變。本文采用基于彈簧-梁模型的方法進行力的等效轉(zhuǎn)換，避免結構單元中心點位置與氣動壓力點位置相距過近時產(chǎn)生計算誤差過大的問題。

首先假設在-平面內(nèi)，載荷方向同升力方向。假設某點的氣動力為，則合力、合力矩等效方程為

(8)

(9)

(10)

式中：為要等效的結構單元中心點數(shù);、為結構單元中心點的坐標；、為氣動壓力點的坐標。

假設各個有限元節(jié)點和氣動壓力點之間有一根“懸臂梁”，在氣動載荷點用一個彈簧約束其滑動，其他方向位移固定，設其自由端上的有限元節(jié)點分配到的載荷為，此時其變形能為

(11)

式中：為假想梁的長度；為假想梁的抗彎剛度；為彈簧彈性系數(shù)。

此時假想梁和彈簧的總變形能為

(12)

為使系統(tǒng)的變形能最小，可采用拉格朗日乘子法建立拉格朗日函數(shù)

(13)

式中：,,為Lagrange乘子。

為方便計算，可令3=1，將式(13)分別對求偏導并令其為0，得

(14)

將式(14)代入式(8)～式(10)，并令=(3)得到簡化后的方程組為

(15)

解出、、后,回代式(14)中即可得到各有限元節(jié)點所分配到的載荷。

2 混合法原理

火箭表面的凸起處通常是氣動模型和結構模型外形不一致區(qū)域。某型運載火箭表面形狀的局部如圖1所示。因此考慮通過分析氣動壓力點距軸心的距離來判斷氣動壓力點是否位于表面凸起處。

對于具有復雜表面的火箭結構，分別從火箭氣動和結構模型中導出氣動壓力點信息與結構單元信息，對于每個氣動壓力點，沿軸向查找距其最近的2個結構單元中心點，分別對比這幾個點到火箭軸心的距離，從而判斷該氣動壓力點是否在結構模型表面上，實現(xiàn)氣動外形與結構外形一致和不一致區(qū)域的劃分。對于氣動外形與結構外形一致的區(qū)域采用壓力插值法，而對于外形不一致區(qū)域采用力等效法，最終結構模型的載荷為2種載荷的疊加。

圖1 火箭外形細節(jié)圖Fig.1 Detail drawing of rocket shape

2.1 區(qū)分氣動壓力點

首先計算出每個氣動壓力點、結構單元中心點到火箭軸心的距離。由于火箭沿軸向直徑會變化，因此查找每個氣動壓力點與沿軸向距該點最近的2個結構單元中心點，并且這2個結構單元中心點分別位于氣動壓力點左右兩側，采用點斜式建立上述2個結構單元中心點組成的直線方程，通過該方程求出氣動壓力點所在位置處的基準半徑，如圖2所示。比較氣動壓力點的基準半徑與實際半徑:

|-|>

(16)

式中：為判斷一點是否在面內(nèi)的臨界值，本文中取0.005 m，相當于氣動模型中的最小網(wǎng)格尺寸。

該研究參考AWL，采用語料庫技術，以國內(nèi)使用甚廣的大學英語教材《新視野》為研究對象，計算此套教材的AWL覆蓋率，試圖回答以下四個具體問題：

如若滿足式(16)，則認為該氣動壓力點位于凸出區(qū)域，即該點與結構模型不一致，反之，則認為該點與結構模型外表面一致。

圖2 識別不一致氣動壓力點的方法示意圖Fig.2 Schematic diagram of identification method for inconsistent aerodynamic pressure points

2.2 區(qū)分結構單元中心點

將第1步中區(qū)分出的不一致區(qū)域的氣動壓力點投影到基準半徑的圓環(huán)上，如圖3所示，其中

(17)

(18)

式中：、為面外氣動壓力點的坐標值；′、′為面外氣動壓力點投影到面內(nèi)后的坐標值。

將距投影點最近的結構單元中心點作為不一致區(qū)域的結構單元中心點，其余點作為一致區(qū)域的結構單元中心點。

圖3 面外氣動壓力點向面內(nèi)投影過程Fig.3 Projection process from aerodynamic pressure point out of plane to rocket surface

2.3 用兩種方法分別計算對應區(qū)域

首先使用TPS法對第1步中輸出的一致區(qū)域的氣動壓力點與第2步輸出的一致區(qū)域的結構單元中心點插值，得到一致區(qū)域的結構單元中心點的壓力。

然后使用力等效法轉(zhuǎn)換不一致區(qū)域的氣動力。先將不一致區(qū)域氣動壓力點的壓力乘單元面積轉(zhuǎn)換為氣動力，然后根據(jù)結構單元中心點的信息找到對應的結構網(wǎng)格節(jié)點，采用力等效法進行載荷轉(zhuǎn)換，得到結構網(wǎng)格節(jié)點上所受到的力。

3 算 例

本文的算例為某型運載火箭，氣動模型在FLUENT中建立，火箭外形及氣動載荷壓力分布如圖4所示，表面氣動壓力點共614 209個。結構模型在PATRAN中建立，結構有限元模型網(wǎng)格如圖5所示，表面結構單元共41 596個，結構網(wǎng)格節(jié)點共39 951個。

圖4 火箭外形及氣動壓力分布圖Fig.4 Launch vehicle shape and aerodynamic pressure distribution

圖6 氣動壓力點整體分布圖Fig.6 Overall distribution of aerodynamic pressure points

圖7 氣動壓力點局部細節(jié)圖(8 m

圖8 結構單元中心點整體分布圖Fig.8 Overall distribution of center points of structural elements

圖9 結構單元中心點局部細節(jié)圖(8 m

壓力插值是在氣動壓力點與結構單元中心點之間進行轉(zhuǎn)換，因此可直接使用TPS法轉(zhuǎn)換外形一致區(qū)域的氣動力，圖10為外形一致區(qū)域的結構單元壓力云圖。而力等效法是在氣動壓力點與結構網(wǎng)格節(jié)點之間進行轉(zhuǎn)換，所以要先獲取3 211個結構單元的節(jié)點，共4 201個，再使用力等效法將氣動壓力點載荷轉(zhuǎn)換到這些結構網(wǎng)格節(jié)點上，所得力的云圖如圖11所示。

最后對比載荷轉(zhuǎn)換前后整體的合力和合力矩，如表1所示,各誤差均小于3%，其中最嚴重的軸向載荷誤差小于1%。該算例的計算結果表明，該方法可實現(xiàn)表面有凸起物的火箭的三維氣動載荷直接向三維結構模型的等效轉(zhuǎn)換，并且能夠滿足工程上的精度要求。

圖10 一致區(qū)域施加的壓力載荷分布Fig.10 Distribution of pressure loads in consistent areas

圖11 不一致區(qū)域施加的節(jié)點力載荷分布Fig.11 Nodal force load distribution applied in inconsistent areas

表1 混合法所得合力、合力矩的誤差

4 結 論

針對在運載火箭結構載荷精細化計算過程中遇到的氣動模型與結構模型外形不匹配的問題，本文提出了一種基于壓力插值法和力等效法的混合法，該混合法能自動識別出流體模型和結構模型外表面是否一致，從而針對不同的區(qū)域采用不同的方法進行載荷轉(zhuǎn)換。

本方法適用于運載火箭或?qū)椚S有限元模型的流固載荷轉(zhuǎn)換，并且具有轉(zhuǎn)換精度高、轉(zhuǎn)換過程自動化、工程適用性好的特點。

本文中選取了TPS法作為壓力插值方法，實際上也可采用徑向基函數(shù)法等其他方法，選擇不同的方法可以得到不同的轉(zhuǎn)換精度。此外，判斷氣動壓力點是否在結構模型表面時所用臨界值目前是人為設置的，該值會影響一致和不一致區(qū)域的劃分，進而影響載荷轉(zhuǎn)換精度，其自動確定方法也有待進一步研究。