基于徑向點插值方法的柔性螺旋槳氣動彈性模擬

張 宇，王曉亮

(上海交通大學 航空航天學院，上海 200240)

流固耦合問題廣泛存在于航空航天、船艦和能源等領域[1-3].氣動彈性問題的本質與流固耦合問題一致，常見于飛行器設計領域，主要涉及機翼顫振和螺旋槳變形等[4-5].作為常規飛行器(推進器)中的重要分支——螺旋槳動力飛行器(推進器)在航空(航海)領域始終扮演著關鍵角色.螺旋槳作為早期的動力系統，因其在低速時具有高推力、低能耗等優點，被廣泛使用在平流層飛艇、運輸機和傾轉旋翼機上[6-7].典型的飛行器有A400M運輸機、C-130運輸機、運12、AG600水陸兩用飛機、V-22傾轉旋翼機等.普通螺旋槳一般由合金材料制造，故在設計時把螺旋槳槳葉視為剛體而忽略在運行過程中可能發生的變形.但隨著飛行器的載荷能力和機動性能的不斷提高，普通的合金材料對高強度氣動載荷和周期性疲勞振動的承受能力會大大降低.因此，新型的碳纖維材料逐漸進入了人們的視野[8].相較于傳統的金屬螺旋槳，碳纖維材質的螺旋槳具有更好的抗疲勞性，能在一定程度上降低整機重量，提高有效荷載能力.但由于粘接和復合工藝的限制會造成彈性模量較低，在運轉時易發生變形，此時的螺旋槳推進性能將會發生改變.如果忽略氣動彈性，往往又會造成設計點的偏移，使真實情況與預期有較大的偏差.因此，在對柔性碳纖維螺旋槳進行推進性能預測時，有必要考慮材料的氣動彈性效應.

在針對螺旋槳氣動彈性問題的研究方法上，Sodja等[9]結合葉素動量理論和非線性梁理論，建立了柔性螺旋槳的氣動彈性簡化模型，并預測了前掠、后掠和平直型螺旋槳的推進性能結果.Hsiao等[10]耦合了可壓/不可壓流體求解代碼和有限元固體求解代碼，對多層復合材料螺旋槳的推進性能進行了研究，最后給出了合適的鋪陳方式以提高推力.Zhang等[11]依托ANSYS Workbench軟件中的System Coupling平臺，研究了在流固耦合作用下夾層復合材料螺旋槳的推進效率和結構響應.Das等[12]使用一階切應變理論和雷諾平均Navier-Stokes(RANS)方程研究了全尺寸多層粘合彎扭螺旋槳的推進特性.王建等[13]采用面元法與有限元法相結合的流固耦合分析方法對碳纖維螺旋槳進行了水動力性能分析，并借助iSIGHT軟件構建起基于響應面近似模型碳纖維螺旋槳的多目標優化策略.婁本強等[14]分別使用直接耦合數值計算和瞬態激勵實驗方法，研究了某型螺旋槳葉片在空氣中和浸入靜水域中的固有流固耦合振動特性.

從上述文獻可見，雖然現有研究已經建立起了關于螺旋槳氣動彈性的分析技術，但已有的研究方法一般是通過簡化的理論方法將流場和位移場方程耦合并求解，這種方式在應用對象上有一定的局限性，對于復雜的工況不能完全反映其流場情況.此外，通過封裝好的平臺實現耦合過程中的數據傳遞，容易導致數據不透明，降低了信息提取的自由度.

本文借助成熟的計算流體力學(CFD)/計算固體力學(CSD)求解器分別求解流場和位移場，CFD/CSD模塊彼此獨立，可操作性較強，便于數據信息的提取.應用無網格的徑向點插值方法(RPIM)完成固體網格節點朝流體網格節點的位移傳遞[15].結合位移傳遞矩陣和虛位移原理，生成節點力傳遞矩陣，該方法能保證數據在傳遞中力、力矩和能量的守恒.在動網格方面，采取Delaunay映射法進行流場網格更新，相較于傳統的彈簧光順等方法，Delaunay映射法具有更高的計算效率和穩健性，適用于任意拓撲類型的網格[16].

1 幾何模型

1.1 螺旋槳幾何模型

螺旋槳模型的基礎剖面翼型為Clark Y，通過片條理論獲得螺旋槳葉片的氣動響應，再由遺傳算法優化獲得最佳的螺旋槳設計參數.該螺旋槳的直徑Dp為4.6 m，槳轂直徑Dg為0.92 m.由于槳轂對葉片的變形以及對螺旋槳周圍的流場結構影響較小，為簡化問題去除槳轂結構，最終獲得的外形如圖1所示.

圖1 螺旋槳幾何外形示意圖Fig.1 Schematic diagram of propeller geometry

1.2 CFD計算域模型

CFD計算域模型如圖2所示，整個模型包括一個旋轉域和一個靜止域，FLUID_PROP為包含螺旋槳PROP的旋轉域，FLUID為靜止域.計算域尺寸如圖3所示.計算域入口Inlet到旋轉域FLUID_PROP的距離為10Dp，旋轉域FLUID_PROP到出口Outlet的距離為30Dp，螺旋槳中心與遠場壁面的距離為8Dp，旋轉域FLUID_PROP的直徑為1.5Dp，旋轉域FLUID_PROP沿來流方向的尺寸為0.5Dp.旋轉域FLUID_PROP與靜止域FLUID之間的交界面類型為Interface，通過插值進行數據交換，各邊界名稱與類型統計如表1所示.

圖2 流體計算域Fig.2 Fluid computational domain

圖3 計算域尺寸Fig.3 Computational domain size

表1 邊界條件分類Tab.1 Classification of boundary conditions

2 計算方法與可靠性驗證

2.1 Delaunay映射方法

流場網格更新通過Delaunay方法實現，該方法具有很高的穩健性，避免了求解大型矩陣的高耗時缺點.Delaunay四面體如圖4所示.流體域內任意一個網格節點P在Delaunay四面體ABCD的位置由體積坐標(e1,e2,e3,e4)唯一確定：

圖4 Delaunay四面體Fig.4 Delaunay tetrahedron

(1)

(2)

(3)

(4)

(5)

(6)

顯然對于體積坐標有：

(7)

因此，當物面網格節點發生位移u后，流體域內網格節點的位移可由下式獲得：

(8)

使用Delaunay映射方法對NACA0012二維網格進行變形的結果如圖5所示.翼型表面的運動方式為ΔY=-0.1sin(2πX)，X、Y為無量綱長度.變形前流場網格的最小正交質量為0.445，變形后流場網格的最小正交質量為0.124.從變形前后的翼型周圍網格分布來看，Delaunay映射方法能較好地完成貼體網格及流場網格的更新.

圖5 NACA0012二維翼型網格變形Fig.5 NACA0012 2D airfoil grid deformation

2.2 可靠性驗證

采取成熟的商業軟件ANSYS Fluent?和Abaqus?分別求解氣動載荷和結構位移.流場網格劃分工具為ICEM CFD，在靜止域FLUID采用結構網格剖分，在旋轉域FLUID_PROP內使用非結構網格填充，在螺旋槳周圍進行網格加密，所獲得的最終網格如圖6所示.網格單元總數為 2 134 807，旋轉域內節點數為 304 975，單片槳葉表面節點數為 13 535.考慮到螺旋槳周圍流場存在強烈的流動分離情況，湍流模型選取k-ω剪切應力輸運(SST)模型.同時鑒于螺旋槳運行時的動平衡特點，采取多重參考系(MRF)方法求解旋轉流場，該方法通過求解準靜態的RANS方程獲取旋轉域內的流場信息.MRF方法的主要思想是將螺旋槳相對于靜止流場的運動轉變為旋轉流體相對于靜止螺旋槳的運動，通過隱式交界面進行旋轉域與靜止域之間的數據交換，以保證各通量守恒[17].盡管該方法是一種近似處理，與較真實的瞬態方法有一定的差異，但對螺旋槳的流動預測依然是有效的，且所需計算的資源較瞬態方法低很多[18].

圖6 螺旋槳周圍流體域網格Fig.6 Fluid grids around propeller

為驗證本文所提計算方法CFD_MRF的可靠性，在此選取某三葉螺旋槳縮比模型進行風洞試驗分析，該三葉螺旋槳的實物模型如圖7所示.實驗時的當地大氣壓為9.31×105Pa，實驗風洞溫度為298 K(25 ℃)，來流風速為22 m/s，轉速設置為 3 186、3 536、3 886、4 236、4 586、4 936、5 086、5 236 r/min.上述工況下的推力和扭矩的變化曲線如圖8和9所示.其中：F為螺旋槳推力；M為螺旋槳扭矩；n為螺旋槳轉速；EXP表示風洞實驗.由圖8和9可知，CFD_MRF方法與實驗獲得的結果相近，推力和扭矩的平均誤差分別為5.94%和5.90%，且變化趨勢保持一致.該算例表明所提數值方法能夠較好地模擬螺旋槳的真實氣動效應.

圖7 三葉螺旋槳實物圖Fig.7 Picture of three-bladed propeller

圖8 CFD方法與風洞實驗所獲得的推力對比Fig.8 Comparison of thrusts between CFD and wind tunnel experiment

圖9 CFD方法與風洞實驗所獲得的扭矩對比Fig.9 Comparison of torques between CFD and wind tunnel experiment

螺旋槳由碳纖維復合材料構成，蒙皮鋪層采用3621碳布預浸料(平紋布0°/90°)、M40JB(單向帶)混合鋪層形式.腹板鋪層采用3621碳布預浸料鋪層形式，其密度為 1 600 kg/m3，彈性模量近似取為70 GPa，泊松比為0.3，采取C3D8R單元進行固體網格剖分.為了消除網格因素帶來的誤差，通過改變螺旋槳翼型弦向網格份數和展向網格份數形成稀疏、粗糙、中等和精細4種密度等級的固體網格，4種網格的信息如表2所示.在螺旋槳葉片表面施加均布壓力，并比較其變形結果.Abaqus網格無關性分析如圖10所示.其中：ξp為螺旋槳葉片無量綱徑向位置；ux為螺旋槳葉片尾緣在x方向上的位移.由圖10可知，中等密度的網格已經能很好地滿足計算精度的要求.

表2 固體模型網格參數(單個螺旋槳)Tab.2 Parameters of meshes in solid model (single propeller)

圖10 Abaqus網格無關性分析Fig.10 Grid independence analysis in Abaqus

3 RPIM

RPIM結合了徑向基函數的對稱正定特點和無網格方法的優點，可以生成非奇異的插值矩陣，保證了數值穩定性，適用于任意分布的節點.而傳統的徑向基函數(RBF)，只能對氣動壓力分布進行較好的映射，而無法保證積分載荷(力與力矩)的守恒.而在RPIM中引入了線性基，可對積分載荷進行守恒插值，插值效果更接近真實情況.

3.1 位移傳遞矩陣

假設已知S個點的位移，待求點xh的位移u(xh)可由如下插值格式獲得：

RT(xh)a+PT(xh)b

(9)

式中：ai為徑向基函數Ri的未知權重系數；bj為多項式基函數Pj的未知系數；M為Pj中基底的個數；R為徑向基函數矩陣；P為多項式基函數矩陣；a為權重系數向量；b為未知系數向量.

未知系數ai和bj通過求解類似于式(1)的方程獲得，這里需要用到已知的S個點的位移，第t個點的位移ut為

(10)

t=1,2,…,S

將式(10)寫成矩陣的形式為

u=Ra+Pb

(11)

式中：u為所有已知點的位移向量.

為保證方程封閉，需在式(10)的基礎上施加如下約束：

(12)

j=1,2,…,M

為保證數據傳遞中力和力矩的平衡，要求多項式基函數至少取一階，故M可取為4，則P(x)=[1xyz].聯立式(10)和(12)可得：

(13)

或

(14)

其中：

此處的RBF選取Wendland二階光順格式[19]

Ri(xk,yk,zk)=(1-?)4(4?+1)

(15)

式中：?=li,k/rs，rs為徑向基函數的支撐半徑，li,k為點(xi,yi,zi)與點(xk,yk,zk)之間的距離.由于距離沒有方向之分，所以矩陣R是對稱正定的，則矩陣G也是對稱的.

假設G可逆，則系數向量a為

a=R-1u-R-1Pb

(16)

再由式(13)可得：

b=S2u

(17)

(18)

將式(17)代入式(16)可得：

a=S1u

(19)

S1=R-1-R-1PS2

(20)

將式(17)～(20)代入式(9)可得待求點的位移為

u(xh)=[RT(xh)S1+PT(xh)S2]u

(21)

因此位移傳遞矩陣H可寫為

式中：N為待求點個數.

3.2 節點力傳遞矩陣

為防止造成能量損失，可通過虛位移原理獲得節點力傳遞矩陣.假設Ff和Fs分別為流體網格節點和固體網格節點上的力向量，uf和us為對應節點上的節點位移，由虛位移原理有：

(22)

已知位移傳遞矩陣H，即

uf=Hus

(23)

代入式(22)則有：

Fs=HTFf

(24)

3.3 CFD/CSD迭代過程

進行氣動彈性計算時，需要在每個迭代步交換數據，該過程如圖11所示.CFD模塊將由計算獲得的氣動力數據①插值傳遞給CSD模塊；然后CSD模塊以①為輸入載荷，通過計算過程②獲得位移場數據③，并插值傳遞給CFD模塊；CFD接收到位移后進行網格更新過程④，進而獲得網格更新后的氣動力數據⑤；通過插值將氣動力數據⑤傳遞給CSD模塊，如此完成一個迭代步的數據交換.隨著迭代次數的增加，迭代結果會趨于收斂，此時可終止計算.

圖11 CFD/CSD模塊迭代過程Fig.11 Iteration process of CFD/CSD modules

3.4 RPIM插值精度及其效率對比

為說明使用RPIM進行流固耦合數據傳遞的必要性，分別使用RPIM和RBF對某Clark Y機翼表面進行氣動力傳遞，徑向基函數均如式(15)所示.NACA0014機翼插值精度驗證如圖12所示.由圖12(a)可知，該機翼展向長度為10.5 m，機翼根部弦向長度為3.0 m，機翼后掠角為20.6°.該機翼的流場表面網格及固體網格如圖12(b)和(c)所示.其中，機翼表面網格節點數為 3 140，固體網格表面節點數為 1 878.設定來流馬赫數Ma=0.839 5，攻角α=3.06°，當地大氣壓為 100 312 Pa，當地溫度為300.9 K.

圖12 NACA0014機翼插值精度驗證Fig.12 Interpolation accuracy verification for NACA0014 wing

對固體而言，表面載荷最終會轉化為一致節點力，故此處直接選取網格節點力進行傳遞，避免再進行表面壓力轉換.表3統計了使用上述兩種方式計算獲得的各自傳遞矩陣并進行一次氣動力傳遞后，在固體模型3個方向的合力(Fx、Fy、Fz)和力矩(Mx、My、Mz)大小.其中：er為相對誤差；tc為計算耗時.由表3可知，在足夠的精度范圍內使用RPIM幾乎沒有能量損失，而由RBF獲得的結果卻與實際情況相差較遠，最大相對誤差為-348.8%.導致該結果的原因是RBF在數據傳遞過程中沒有引入線性基，無法保證力與力矩的守恒，使用RBF傳遞節點力無法保證數據的有效性.從計算時間上來看，應用兩種方式的總耗時分別為 20.03 s 和 32.13 s，RPIM的時間成本較RBF節省了約37.7%.對于相同的網格規模，雖然RPIM中的矩陣維度大于RBF，但由于RPIM只需求解一種傳遞矩陣(節點力傳遞矩陣或位移傳遞矩陣)，通過轉置即可獲得另一傳遞矩陣，而RBF需進行兩次獨立的矩陣求解，所以RBF在求解速度上較RPIM更慢.綜上，所提RPIM插值方法能以較高的精度保證數據的準確性，在相同規模的網格前提下，能以更快的速度進行數據傳遞.

表3 CFD/CSD模塊節點力傳遞結果Tab.3 Node force transfer results of CFD/CSD modules

4 螺旋槳氣動彈性分析

設定螺旋槳飛行高度Hf=0 km，轉速n=400～1 000 r/min，來流速度v=10 m/s，來流方向垂直于螺旋槳平面.當n=800 r/min時的氣彈耦合過程如圖13所示.其中：Is為計算迭代步數；ux,max為槳葉在x方向上的最大位移.由圖13可知，當迭代到第5步時，結果業已收斂.故在后續分析中可約定最大迭代步數為5步，以第0步和第5步的流場計算結果作為變形前后氣動彈性效應的對比依據.

圖13 當n=800 r/min時的氣彈耦合過程Fig.13 Aeroelastic coupling process at n=800 r/min

4.1 幾何外形變化

當n=400～1 000 r/min時，槳葉尾緣在x和y方向的變形量變化情況如圖14所示，其中uy為槳葉尾緣在y方向上的位移.由圖14可知，當ξp<0.3時，位移量與徑向位置呈現出近似二次關系；當超過該值后，位移量開始呈線性增加.說明槳轂的約束可以有效地擬制槳葉的變形能力，但作用范圍只有槳葉徑向30%的區域.槳葉在不同方向上的最大位移量與轉速之間的關系如圖15所示，其中umax為槳葉最大位移.由圖15可知，當螺旋槳受到旋轉氣動載荷時，其主要變形不僅僅發生在來流方向，在旋轉平面上依舊有相當可觀的變形發生，uy,max平均可以達到ux,max的52.1%.當n=800 r/min時，ux,max為130.8 mm，達到了槳葉半徑的5.7%.當n=1 000 r/min時，ux,max達到槳葉半徑的9.4%，該量級的變形量說明螺旋槳在運行過程中產生的氣彈效應不能被忽略.當n=1 000 r/min時，槳葉的扭轉角Δε隨槳葉徑向位置的變化情況如圖16所示.由圖16可知，附加于螺旋槳表面的氣動力會迫使槳葉的當地迎角增大，這將有利于推力的提高.

圖14 槳葉尾緣的位移曲線Fig.14 Displacement curves of blade trailing edge

圖15 x和y軸方向上的最大位移變化曲線Fig.15 Maximum displacements along x and y axes

圖16 槳葉扭轉角沿徑向的變化Fig.16 Radial variations versus blade torsion angles

4.2 氣動特性變化

螺旋槳變形前后的槳葉在背風面上的靜壓分布如圖17所示，其中p為靜壓.由圖17可知，在背風面上氣動載荷幾乎不因槳葉的變形而發生變化.螺旋槳變形前后的槳葉在迎風面上的靜壓分布如圖18所示.由圖18可知，其靜壓分布與圖17有明顯的差異.綜合分析圖17和18可知，在靠近槳轂的區域，剛性與柔性槳葉的靜壓分布幾乎一致，說明變形未對該區域的流場造成影響；在遠離槳轂的區域，即如紅色虛線圈內所示，同一正靜壓在變形后的槳葉迎風面上占據更大的區域.這種靜壓分布結果會造成變形后的螺旋槳產生更大的推力和扭矩.結合圖16與文獻[20]中的翼型受力分析可知，造成柔性螺旋槳推力扭矩增大的原因在于槳葉的扭轉變形增大了當地翼型的安裝角，進而提高了當地迎角，使得迎風面正壓區域擴大.

圖17 螺旋槳背風面壓力分布云圖Fig.17 Distribution contour of propeller leeward pressure

圖18 螺旋槳迎風面壓力分布云圖Fig.18 Distribution contour of propeller windward pressure

4.3 推進性能變化

前進比J，推力系數KF，扭矩系數KM和推進效率η的定義分別為

(25)

(26)

(27)

(28)

式中：ρ為空氣密度.

變形前后推力系數隨轉速的變化曲線如圖19所示，其中：Cr為剛性與柔性螺旋槳物理量的相對改變量.由圖19可知，變形前后螺旋槳的推力變化趨勢是一致的，增長速度由快變慢，但變形后的螺旋槳較剛性螺旋槳能夠產生更大的推力，且差距隨著轉速的增加而增加，最大改變量達到7.2%.扭矩系數的變化情況如圖20所示.由圖20可知，當不考慮螺旋槳變形時，扭矩系數先增加繼而保持相對的穩定；而對于柔性螺旋槳而言，當n>500 r/min后，扭矩系數基本上保持線性增加，最大改變量可達9.9%.螺旋槳推進效率的變化情況如圖21所示.由圖21可知，隨著負載的增加，剛性和柔性螺旋槳的效率逐漸降低，柔性螺旋槳效率略低于剛性螺旋槳，兩者的變化趨勢整體上保持一致，最大改變量僅為-2.4%，可認為氣動彈性在本文工況下基本不會影響螺旋槳的推進效率.

圖19 氣動彈性對推力系數的影響曲線Fig.19 Influence of static aeroelasticity on thrust coefficient

圖20 氣動彈性對扭矩系數的影響曲線Fig.20 Influence of static aeroelasticity on torque coefficient

圖21 氣動彈性對效率的影響曲線Fig.21 Influence of static aeroelasticity on efficiency

5 結論

(1) 當螺旋槳運轉時，在旋轉平面上的變形約為來流方向上的52.1%，當轉速升至 1 000 r/min時，沿來流方向的最大變形量為槳葉半徑的9.4%，氣動彈性效應不能被忽略.

(2) 螺旋槳背風面上壓力分布幾乎不因槳葉的變形而發生變化，遠離槳轂時，槳葉的變形增大了當地翼型的迎角，使得同一正壓在變形后的槳葉迎風面上占據更大的區域.

(3) 相較于剛性螺旋槳，柔性螺旋槳能夠產生更大的推力，且差距隨轉速的增加而增加，最大改變量達7.2%.當轉速高于500 r/min后，扭矩系數基本保持線性增加，最大改變量達9.9%.隨著轉速的增加，柔性和剛性螺旋槳的效率略有降低，最大改變量為-2.4%，可認為柔性螺旋槳的推進效率基本不受氣動彈性的影響.

(4) 本文形成的柔性螺旋槳氣動彈性分析框架，對各類鋪層結構的螺旋槳均適用，實際應用中針對具體情況修改材料的力學屬性即可.