基于SPH 算法的機翼油箱晃動分析

趙博偉，賴 輝

（成都飛機工業（集團）有限責任公司技術中心，四川 成都 610092）

液體的自由表面由于受到激勵或外加擾動而產生的運動稱為液體晃動。航空航天、船舶、石油化工及核動力等領域廣泛存在液體晃動問題。在航空領域，飛機起飛、飛行及著陸過程中由于氣動紊亂、發動機轉子的不平衡、飛機發射等原因導致的油箱內燃油的晃動，會造成不良的影響。第一，晃動會使燃油循環往復地沖擊油箱結構，產生沖擊載荷，造成油箱結構的疲勞破壞[1]；第二，飛機全機的重心分布會受到燃油重心變化的影響，機翼上的重心分布也是影響飛機姿態穩定性的因素之一[2-3]。

因此，研究燃油晃動產生的沖擊力、晃動力矩所引發的油箱結構強度問題，對于飛行器機翼油箱（尤其是某型無人機大展弦比機翼油箱）結構和強度設計可以提供一定的工程技術支持和指導。

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

1 計算方法

1.1 流固耦合

流固耦合力學是研究變形固體在流場作用下的各種行為以及固體位形對流場影響這兩者相互作用的一門科學。按其耦合機理可分為兩大類：一類為耦合作用只發生在兩相交界面上，方程上的耦合是由兩相耦合面上的平衡來引入的；二類為兩域部分或全部重疊，建立方程需要針對具體物理現象，通過描述問題的微分方程來反映耦合效應。

1.2 SPH 方法

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

拉格朗日形式的流體動力學N-S 方程可表示為：

式中，F為單位質量的力；σ為總應力張量，由各項同性壓力和黏性應力兩部分組成。應用SPH 粒子近似法對動量方程進行變換可得如下表達式：

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

1.3 關鍵問題

為了準確模擬流體動力學問題，防止在沖擊域內的求解結果產生非物理震蕩，對SPH 算法進行特別處理，采用Monaghan 型人工黏度[5]。表達式如下所示。

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

使用SPH 方法進行流體模擬時，選取初始光滑長度會對計算結果有重要影響。對于油箱晃動中復雜流體大變形這種過程，必須考慮粒子間距的變化來重新計算適合的光滑長度。故使用按照平均密度來修正光滑長度的方法。

對不可壓縮流，通過采用人工壓縮率將其視為可壓縮流，從而達到求解壓力對時間導數的要求[4]。

2 晃動分析

2.1 計算模型

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

表1 材料參數

表2 機翼油箱部件材料分配和板厚

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

表3 機翼油箱充油量具體參數

為減少機翼油箱模型的單元數及單元尺寸，對油箱采用四節點的縮減積分殼單元（S4R）和三節點殼單元（S3）進行劃分，燃油采用PC3D 粒子單元（由C3D8R實體單元轉換）。其中機翼油箱單元數為22 655，燃油單元數根據多油及半油分別為63 466 和45 694。機翼油箱燃油及單元劃分如圖1 所示。

圖1 機翼油箱燃油及單元劃分

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

2.2 仿真計算

機翼油箱晃動采用SPH 算法，使用ABAQUS/Explicit求解器進行求解。前文提到過機翼油箱的材料為復合材料及鋁材。

2.2.1 半油晃動應力和應變

在半油晃動過程中，鋁材的最大Mises 應力為66 MPa，纖維材料的最大縱向應力為45 MPa，最大橫向應力為17 MPa，均滿足強度要求。最大橫向應變為2 054 με。最值點時，部分部件的應力和應變分布如圖2-圖4所示。各部件最大應力應變如表4 所示。

表4 各部件最大應力應變

圖2 半油晃動，機翼下蒙皮的應力應變在最值點時的分布

圖3 半油晃動，機翼框架鋁材的應力應變在最值點時的分布

圖4 半油晃動，機翼框架復合材料的應力和應變在最值點時的分布

2.2.2 半油晃動速度與位移

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

2.2.3 半油晃動壓力

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

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

2.2.4 多油晃動應力和應變

在多油晃動過程中，鋁材的最大Mises 應力為89 MPa，纖維材料的最大縱向應力為85 MPa，最大橫向應力為18 MPa，均滿足強度要求。最大橫向應變為2 075 με。最值點時，各部件的應力和應變分布如圖5-圖7 所示。各部件最大應力應變如表6 所示。

表6 各部件最大應力應變

圖5 多油晃動，機翼下蒙皮的應力應變在最值點時的分布

圖6 多油晃動，機翼框架鋁材的應力應變在最值點時的分布

圖7 多油晃動，機翼框架復合材料的應力和應變在最值點時的分布

2.2.5 多油晃動速度與位移

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

2.2.6 多油壓力

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

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

2.3 小結

本文中對機翼油箱晃動考慮了多油（68%）和半油（50%）2 種充油量。油箱的應力和應變峰值均小于材料的拉伸強度，滿足強度要求。其中多油較半油的應力和應變較大。連接單元處的最大相對位移均小于0.250 mm，不會發生漏油。晃動過程中，剛體運動占主導，蒙皮和框架的變形很小，兩油艙各壁面的最大平均壓力均小于機翼油箱的沖壓載荷，多油較半油的平均壓力峰值較小。

3 結論

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