基于SPH 算法的機(jī)翼油箱晃動(dòng)分析

趙博偉，賴(lài) 輝

（成都飛機(jī)工業(yè)（集團(tuán)）有限責(zé)任公司技術(shù)中心，四川 成都 610092）

液體的自由表面由于受到激勵(lì)或外加擾動(dòng)而產(chǎn)生的運(yùn)動(dòng)稱(chēng)為液體晃動(dòng)。航空航天、船舶、石油化工及核動(dòng)力等領(lǐng)域廣泛存在液體晃動(dòng)問(wèn)題。在航空領(lǐng)域，飛機(jī)起飛、飛行及著陸過(guò)程中由于氣動(dòng)紊亂、發(fā)動(dòng)機(jī)轉(zhuǎn)子的不平衡、飛機(jī)發(fā)射等原因?qū)е碌挠拖鋬?nèi)燃油的晃動(dòng)，會(huì)造成不良的影響。第一，晃動(dòng)會(huì)使燃油循環(huán)往復(fù)地沖擊油箱結(jié)構(gòu)，產(chǎn)生沖擊載荷，造成油箱結(jié)構(gòu)的疲勞破壞[1]；第二，飛機(jī)全機(jī)的重心分布會(huì)受到燃油重心變化的影響，機(jī)翼上的重心分布也是影響飛機(jī)姿態(tài)穩(wěn)定性的因素之一[2-3]。

因此，研究燃油晃動(dòng)產(chǎn)生的沖擊力、晃動(dòng)力矩所引發(fā)的油箱結(jié)構(gòu)強(qiáng)度問(wèn)題，對(duì)于飛行器機(jī)翼油箱（尤其是某型無(wú)人機(jī)大展弦比機(jī)翼油箱）結(jié)構(gòu)和強(qiáng)度設(shè)計(jì)可以提供一定的工程技術(shù)支持和指導(dǎo)。

本項(xiàng)目基于有限元分析軟件ABAQUS，建立機(jī)翼油箱的有限元模型，并進(jìn)行流固耦合動(dòng)力分析，選用光滑粒子流體動(dòng)力學(xué)（Smoothed Particle Hydrodgnamics，SPH）算法對(duì)機(jī)翼油箱進(jìn)行燃油晃動(dòng)下的仿真分析，判斷其是否漏油，是否滿足強(qiáng)度剛度要求。

1 計(jì)算方法

1.1 流固耦合

流固耦合力學(xué)是研究變形固體在流場(chǎng)作用下的各種行為以及固體位形對(duì)流場(chǎng)影響這兩者相互作用的一門(mén)科學(xué)。按其耦合機(jī)理可分為兩大類(lèi)：一類(lèi)為耦合作用只發(fā)生在兩相交界面上，方程上的耦合是由兩相耦合面上的平衡來(lái)引入的；二類(lèi)為兩域部分或全部重疊，建立方程需要針對(duì)具體物理現(xiàn)象，通過(guò)描述問(wèn)題的微分方程來(lái)反映耦合效應(yīng)。

1.2 SPH 方法

本文采用光滑粒子流體動(dòng)力學(xué)方法，是一種無(wú)網(wǎng)格化的Lagrange（拉格朗日）算法，具有自適應(yīng)、無(wú)網(wǎng)格、粒子形式以及拉格朗日單元的特征，避免了傳統(tǒng)拉格朗日解法中的網(wǎng)格變形問(wèn)題，適合處理大變形與沖擊問(wèn)題[4]。

拉格朗日形式的流體動(dòng)力學(xué)N-S 方程可表示為：

式中，F(xiàn)為單位質(zhì)量的力；σ為總應(yīng)力張量，由各項(xiàng)同性壓力和黏性應(yīng)力兩部分組成。應(yīng)用SPH 粒子近似法對(duì)動(dòng)量方程進(jìn)行變換可得如下表達(dá)式：

根據(jù)單位梯度SPH 核粒子近似式，對(duì)動(dòng)量方程進(jìn)行變換整理，并把總應(yīng)力張量分解為各向同性壓力和黏性應(yīng)力，再用SPH 粒子近似法對(duì)能量方程進(jìn)行變換整理，可得到能量方程的最終粒子近似式。

1.3 關(guān)鍵問(wèn)題

為了準(zhǔn)確模擬流體動(dòng)力學(xué)問(wèn)題，防止在沖擊域內(nèi)的求解結(jié)果產(chǎn)生非物理震蕩，對(duì)SPH 算法進(jìn)行特別處理，采用Monaghan 型人工黏度[5]。表達(dá)式如下所示。

式中，αΠ和βΠ為常數(shù)，αΠ取值為1，βΠ取值為1 或者2，這2 個(gè)取值能有效地防止粒子的非物理穿透；c表示聲速；v為粒子的速度矢量；為防止粒子相互靠近時(shí)發(fā)生數(shù)值發(fā)散，計(jì)算中引入因子φ=0.1hij。

使用SPH 方法進(jìn)行流體模擬時(shí)，選取初始光滑長(zhǎng)度會(huì)對(duì)計(jì)算結(jié)果有重要影響。對(duì)于油箱晃動(dòng)中復(fù)雜流體大變形這種過(guò)程，必須考慮粒子間距的變化來(lái)重新計(jì)算適合的光滑長(zhǎng)度。故使用按照平均密度來(lái)修正光滑長(zhǎng)度的方法。

對(duì)不可壓縮流，通過(guò)采用人工壓縮率將其視為可壓縮流，從而達(dá)到求解壓力對(duì)時(shí)間導(dǎo)數(shù)的要求[4]。

2 晃動(dòng)分析

2.1 計(jì)算模型

流固耦合計(jì)算中，將三維實(shí)體油箱模型簡(jiǎn)化為殼體模型，并且將多個(gè)部件合并為單一部件以提高計(jì)算效率。某型大展弦比無(wú)人機(jī)機(jī)翼油箱包括3 個(gè)獨(dú)立的油艙，分別為左艙、中艙和右艙，基于模型及載荷的對(duì)稱(chēng)性考慮，取一半機(jī)翼進(jìn)行分析（保留一半中艙及右艙），簡(jiǎn)化后，計(jì)算模型為上蒙皮，框架及下蒙皮共3 個(gè)。機(jī)翼油箱采用鋁板及復(fù)合材料板，具體參數(shù)如表1 所示。機(jī)翼油箱的復(fù)合材料采用橫向各向異性板建模，其纖維縱向（0°方向）沿著機(jī)翼后梁軸線方向。材料分配和板厚如表2 所示。

表1 材料參數(shù)

表2 機(jī)翼油箱部件材料分配和板厚

燃油為航空燃油，考慮2 種充油量：多油（68%）和半油（50%），如表3 所示。

表3 機(jī)翼油箱充油量具體參數(shù)

為減少機(jī)翼油箱模型的單元數(shù)及單元尺寸，對(duì)油箱采用四節(jié)點(diǎn)的縮減積分殼單元（S4R）和三節(jié)點(diǎn)殼單元（S3）進(jìn)行劃分，燃油采用PC3D 粒子單元（由C3D8R實(shí)體單元轉(zhuǎn)換）。其中機(jī)翼油箱單元數(shù)為22 655，燃油單元數(shù)根據(jù)多油及半油分別為63 466 和45 694。機(jī)翼油箱燃油及單元?jiǎng)澐秩鐖D1 所示。

圖1 機(jī)翼油箱燃油及單元?jiǎng)澐?/p>

機(jī)翼油箱晃動(dòng)周期為0.32 s，最大轉(zhuǎn)角0.262，初角速度3.50 m/s，最大角加速度為408.5 m/s2，變速歷時(shí)0.018 s。不考慮沖壓作用，機(jī)翼油箱的對(duì)稱(chēng)面，運(yùn)動(dòng)耦合在參考點(diǎn)上，對(duì)參考點(diǎn)施加繞著x軸的角加速度和初角速度。在靠近轉(zhuǎn)角峰值時(shí)，轉(zhuǎn)動(dòng)減速，隨后反向加速。油箱及燃油都施加重力載荷且相互之間切向無(wú)摩擦接觸，法向硬接觸。

2.2 仿真計(jì)算

機(jī)翼油箱晃動(dòng)采用SPH 算法，使用ABAQUS/Explicit求解器進(jìn)行求解。前文提到過(guò)機(jī)翼油箱的材料為復(fù)合材料及鋁材。

2.2.1 半油晃動(dòng)應(yīng)力和應(yīng)變

在半油晃動(dòng)過(guò)程中，鋁材的最大Mises 應(yīng)力為66 MPa，纖維材料的最大縱向應(yīng)力為45 MPa，最大橫向應(yīng)力為17 MPa，均滿足強(qiáng)度要求。最大橫向應(yīng)變?yōu)? 054 με。最值點(diǎn)時(shí)，部分部件的應(yīng)力和應(yīng)變分布如圖2-圖4所示。各部件最大應(yīng)力應(yīng)變?nèi)绫? 所示。

表4 各部件最大應(yīng)力應(yīng)變

圖2 半油晃動(dòng)，機(jī)翼下蒙皮的應(yīng)力應(yīng)變?cè)谧钪迭c(diǎn)時(shí)的分布

圖3 半油晃動(dòng)，機(jī)翼框架鋁材的應(yīng)力應(yīng)變?cè)谧钪迭c(diǎn)時(shí)的分布

圖4 半油晃動(dòng)，機(jī)翼框架復(fù)合材料的應(yīng)力和應(yīng)變?cè)谧钪迭c(diǎn)時(shí)的分布

2.2.2 半油晃動(dòng)速度與位移

機(jī)翼油艙的最大粒子速度為6.100 m/s（中油艙）以及3.400 m/s（右油艙），以及機(jī)翼油箱質(zhì)心的最大偏移為0.020 mm（x方向），0.550 mm（y方向），1.650 mm（z方向）。質(zhì)心偏移很小。而油箱整體保持穩(wěn)定，機(jī)翼油箱在晃動(dòng)過(guò)程中剛體運(yùn)動(dòng)占主導(dǎo)，局部變形（撓度）很小。機(jī)翼油箱的連接單元處，兩表面間相對(duì)位移均小于0.250 mm，可以預(yù)計(jì)不會(huì)發(fā)生燃油泄漏。

2.2.3 半油晃動(dòng)壓力

機(jī)翼兩油艙各壁面最大的瞬時(shí)平均壓力位于中油艙1 段的下壁面，為0.032 MPa，但仍小于機(jī)翼油箱的沖壓荷載0.043 MPa。各區(qū)段的平均壓力變化在位移峰值過(guò)后（0.100 s 和0.260 s 附近）出現(xiàn)壓力的峰值。前壁面和下壁面的平均壓力較大（表5）。

表5 機(jī)翼兩油艙各壁面的最大總壓力和平均壓力

2.2.4 多油晃動(dòng)應(yīng)力和應(yīng)變

在多油晃動(dòng)過(guò)程中，鋁材的最大Mises 應(yīng)力為89 MPa，纖維材料的最大縱向應(yīng)力為85 MPa，最大橫向應(yīng)力為18 MPa，均滿足強(qiáng)度要求。最大橫向應(yīng)變?yōu)? 075 με。最值點(diǎn)時(shí)，各部件的應(yīng)力和應(yīng)變分布如圖5-圖7 所示。各部件最大應(yīng)力應(yīng)變?nèi)绫? 所示。

表6 各部件最大應(yīng)力應(yīng)變

圖5 多油晃動(dòng)，機(jī)翼下蒙皮的應(yīng)力應(yīng)變?cè)谧钪迭c(diǎn)時(shí)的分布

圖6 多油晃動(dòng)，機(jī)翼框架鋁材的應(yīng)力應(yīng)變?cè)谧钪迭c(diǎn)時(shí)的分布

圖7 多油晃動(dòng)，機(jī)翼框架復(fù)合材料的應(yīng)力和應(yīng)變?cè)谧钪迭c(diǎn)時(shí)的分布

2.2.5 多油晃動(dòng)速度與位移

機(jī)翼油艙的最大粒子速度為6.600 m/s（中油艙）以及3.60 m/s（右油艙），以及機(jī)翼油箱質(zhì)心的最大偏移為0.020 mm（x方向），0.330 mm（y方向），1.420 mm（z方向）。質(zhì)心偏移很小。而油箱整體保持穩(wěn)定，機(jī)翼油箱在晃動(dòng)過(guò)程中剛體運(yùn)動(dòng)占主導(dǎo)，局部變形（撓度）很小。機(jī)翼油箱的連接單元處，兩表面間相對(duì)位移均小于0.250 mm，可以預(yù)計(jì)不會(huì)發(fā)生燃油泄漏。

2.2.6 多油壓力

機(jī)翼兩油艙各壁面最大的瞬時(shí)平均壓力位于中油艙1 段的前壁面，為0.028 MPa，但仍小于機(jī)翼油箱的沖壓荷載0.043 MPa。各區(qū)段的平均壓力變化在位移峰值過(guò)后（0.100 s 和0.260 s 附近）出現(xiàn)壓力的峰值。前壁面和下壁面的平均壓力較大（表7）。

表7 機(jī)翼兩油艙各壁面的最大總壓力和平均壓力

2.3 小結(jié)

本文中對(duì)機(jī)翼油箱晃動(dòng)考慮了多油（68%）和半油（50%）2 種充油量。油箱的應(yīng)力和應(yīng)變峰值均小于材料的拉伸強(qiáng)度，滿足強(qiáng)度要求。其中多油較半油的應(yīng)力和應(yīng)變較大。連接單元處的最大相對(duì)位移均小于0.250 mm，不會(huì)發(fā)生漏油?；蝿?dòng)過(guò)程中，剛體運(yùn)動(dòng)占主導(dǎo)，蒙皮和框架的變形很小，兩油艙各壁面的最大平均壓力均小于機(jī)翼油箱的沖壓載荷，多油較半油的平均壓力峰值較小。

3 結(jié)論

本文以某型大展弦比無(wú)人機(jī)機(jī)翼油箱晃動(dòng)分析為背景，選用高效的SPH 算法應(yīng)用通用有限元分析軟件ABAQUS 對(duì)機(jī)翼油箱的充液晃動(dòng)進(jìn)行仿真分析，通過(guò)數(shù)據(jù)研究表明，機(jī)翼在油箱晃動(dòng)過(guò)程中滿足強(qiáng)度要求，變形微小，整體穩(wěn)定。文中的數(shù)據(jù)論證對(duì)于飛機(jī)總體設(shè)計(jì)、航天器推進(jìn)貯箱的防晃設(shè)計(jì)等工作具有一定的理論價(jià)值和工程指導(dǎo)意義。