張鈞奕， 蘇新兵， 趙希瑋， 王 振

(空軍工程大學(xué)航空工程學(xué)院， 西安， 710038)

前掠翼具有高升阻比、失速迎角大等優(yōu)勢(shì)，但前掠翼的靜氣彈發(fā)散問題一直制約著其發(fā)展和應(yīng)用[1-3]。國(guó)內(nèi)外學(xué)者針對(duì)前掠翼氣動(dòng)彈性現(xiàn)象與抑制靜氣彈發(fā)散進(jìn)行了許多研究[4-5]，主要手段有氣動(dòng)彈性剪裁[6]、操縱面偏轉(zhuǎn)[7]與應(yīng)用復(fù)合材料[8]。近年來，流動(dòng)控制技術(shù)發(fā)展迅速，給前掠翼靜氣彈抑制提供了新的思路。流動(dòng)控制主要包括主動(dòng)流動(dòng)控制與被動(dòng)流動(dòng)控制，主動(dòng)流動(dòng)控制的優(yōu)點(diǎn)在于可以根據(jù)不同的飛行條件與流動(dòng)特點(diǎn)進(jìn)行相應(yīng)的改變以適應(yīng)不同工作環(huán)境。

等離子體流動(dòng)控制是基于“等離子體氣動(dòng)激勵(lì)”的新概念主動(dòng)流動(dòng)控制技術(shù)[9-10]，它的優(yōu)點(diǎn)是無需運(yùn)動(dòng)部件，對(duì)于機(jī)翼流場(chǎng)主動(dòng)控制具有重要意義。國(guó)內(nèi)外學(xué)者針對(duì)等離子體流動(dòng)控制技術(shù)進(jìn)行了大量的試驗(yàn)和仿真研究，其作用原理是通過設(shè)置于物體表面的電極將附近的空氣電離，產(chǎn)生“離子層”，并使等離子體在電場(chǎng)力作用下運(yùn)動(dòng)，通過帶電粒子與中性氣體分子的碰撞，誘導(dǎo)物體表面空氣定向流動(dòng)，產(chǎn)生射流[11-12]。目前一種重要手段是通過表面介質(zhì)阻擋放電[13](DBD)使等離子體激勵(lì)器產(chǎn)生定向射流，主要有兩大作用：一是增大附面層的流動(dòng)速度；二是向附面層注入能量，提高低能流和高能流的摻混，增強(qiáng)其抵抗逆壓梯度的能力[14]。Shyy Wei等提出了對(duì)介質(zhì)阻擋放電進(jìn)行簡(jiǎn)化的唯象學(xué)模型[15-16]，為等離子體流動(dòng)控制數(shù)值仿真提供了新的思路。

隨著計(jì)算機(jī)性能的提高，流固耦合(CFD/CSD耦合)技術(shù)成為研究氣動(dòng)彈性問題的一種重要數(shù)值仿真方法[17]，劉萬鋼等對(duì)自由飛行狀態(tài)下的細(xì)長(zhǎng)彈箭進(jìn)行了靜氣彈計(jì)算[18]，雷帥等基于AnsysWorkbench對(duì)HIRENASD機(jī)翼靜氣彈特性進(jìn)行了仿真研究[19]。本文借鑒上述研究思路與方法，將等離子體流動(dòng)控制技術(shù)與流固靜力耦合技術(shù)相結(jié)合，通過數(shù)值模擬，初步驗(yàn)證了在前掠翼表面布置等離子體激勵(lì)器進(jìn)行流動(dòng)控制以抑制靜氣彈變形的可行性，其中等離子體流動(dòng)控制的仿真計(jì)算采用唯象學(xué)模型，流固靜力耦合采用Fluent與Mechanical進(jìn)行聯(lián)合仿真，唯象學(xué)模型通過UDF在動(dòng)量方程中添加體積力源項(xiàng)實(shí)現(xiàn)。

1 計(jì)算方法與模型

1.1 流場(chǎng)求解技術(shù)

采用基于有限體積法的Fluent模塊求解三維定常Navier-Stokes方程，在笛卡爾坐標(biāo)系(x1,x2,x3)中，定義速度分量為(u1,u2,u3)，根據(jù)Einstein求和約定，無熱源三維Navier-Stokes方程的守恒形式為：

(1)

式中：φ為狀態(tài)矢量；Fi為對(duì)流通量矢量；Fvi為黏性通量矢量，S為源項(xiàng)矢量，用列向量形式表示如下：

式中：ρ、p、f、E、H、T分別表示密度、壓力、體積力、總能量、總焓和溫度；黏性應(yīng)力張量為：

(2)

式中：k、μ、λ分別表示流體熱傳導(dǎo)系數(shù)、分子黏性系數(shù)和第二黏性系數(shù)；δij為克羅尼柯爾符號(hào)，當(dāng)i=j時(shí)τij為正應(yīng)力，i≠j時(shí)τij為切應(yīng)力。

1.2 結(jié)構(gòu)求解技術(shù)

采用基于有限元法的Static Structural模塊求解結(jié)構(gòu)靜力平衡方程：

(3)

式中：δ(t)為節(jié)點(diǎn)位移分量矩陣；R(t)為節(jié)點(diǎn)載荷矩陣；M、C、K分別為結(jié)構(gòu)的質(zhì)量矩陣、阻尼矩陣和剛度矩陣。

Kδ(t)=R(t)

(4)

采用基于Newton-Raphson方法[20]的迭代求解進(jìn)行結(jié)構(gòu)計(jì)算，逐步增加特定載荷求解結(jié)果，考慮結(jié)構(gòu)的幾何非線性，在每一次迭代時(shí)根據(jù)結(jié)構(gòu)新的幾何位置坐標(biāo)形成新的剛度矩陣。

1.3 前掠翼模型

本文采用的前掠翼數(shù)值計(jì)算模型如圖1所示，選用NACA0015翼型，翼型后緣為半徑0.003 m的圓弧(此模型也被用作進(jìn)行基于Coanda效應(yīng)對(duì)鈍后緣機(jī)翼流動(dòng)控制的相關(guān)研究)；機(jī)翼前、后緣掠角分別為41.19°和54.07°；機(jī)翼半展長(zhǎng)為0.4 m；平均空氣動(dòng)力弦為0.204 m；展弦比為4.3；根梢比為3.4。

圖1 前掠翼模型尺寸

1.4 等離子體唯象學(xué)模型

介質(zhì)阻擋放電等離子體激勵(lì)器的仿真模擬有多種方法，唯象學(xué)模型是對(duì)復(fù)雜的等離子體激勵(lì)流動(dòng)的一種簡(jiǎn)化，以對(duì)流體施加一定的電場(chǎng)力作為等離子體激勵(lì)形式，通過在動(dòng)量方程中添加體積力源項(xiàng)實(shí)現(xiàn)功能，有利于嵌入流體仿真軟件進(jìn)行計(jì)算。本文采用的等離子體唯象學(xué)模型的平面圖如圖2所示。

圖2 等離子體唯象學(xué)模型

在△OAB區(qū)域內(nèi)電場(chǎng)力沿線段AB呈線性分布，方向?yàn)橛葾指向B，其具體分布為：

E(x,y)=E0-k1x-k2y

(5)

Ftave=?αρcecΔtE(x,y)δ

(6)

式中：?為電源電壓頻率；α為電荷碰撞效率因子；ρc為電荷數(shù)密度；ec為元電荷電量常數(shù)；Δt為一個(gè)周期內(nèi)的放電時(shí)間；δ為狄拉克函數(shù)，在基于有限體積法的N-S方程求解中，狄拉克函數(shù)在定義域內(nèi)的積分值為1。

等離子體激勵(lì)區(qū)域如圖3所示，等離子體激勵(lì)器布置在機(jī)翼前緣外側(cè)上表面處，激勵(lì)區(qū)域呈三棱柱形，其在翼型平面內(nèi)的投影均為圖2所示的三角形激勵(lì)區(qū)域。由于布置在機(jī)翼表面，等離子體激勵(lì)區(qū)域與機(jī)翼表面的重合面并非嚴(yán)格意義上的平面，但由于其彎度較小，在計(jì)算時(shí)可以看作平面處理。本文所用的流固耦合方法與等離子體唯象學(xué)模型在文獻(xiàn)[13]、[19]中已有驗(yàn)證，可以保證仿真計(jì)算的正確性和可信度。

圖3 等離子體激勵(lì)區(qū)域

1.5 網(wǎng)格劃分

采用結(jié)構(gòu)網(wǎng)格分別對(duì)流體域與固體域進(jìn)行網(wǎng)格劃分，計(jì)算采用半模，遠(yuǎn)場(chǎng)尺度約為20倍機(jī)翼弦長(zhǎng)，無量綱化壁面距離y+≈1，流體與固體網(wǎng)格數(shù)量分別為2 815 826與116 812，等離子體激勵(lì)流體域網(wǎng)格數(shù)為31 407，流場(chǎng)網(wǎng)格如圖4所示。

圖4 流場(chǎng)網(wǎng)格

較高的網(wǎng)格質(zhì)量是流固耦合仿真運(yùn)算的關(guān)鍵，關(guān)系網(wǎng)格運(yùn)動(dòng)是否合理和耦合運(yùn)算的成敗，由于等離子體激勵(lì)區(qū)域呈三棱柱形狀，所以在劃分網(wǎng)格時(shí)采用Y形網(wǎng)格進(jìn)行處理，保證了即使在尖角處網(wǎng)格縱橫比仍大于0.8。

流體域面網(wǎng)格與結(jié)構(gòu)域面網(wǎng)格采用節(jié)點(diǎn)一致化方案，優(yōu)點(diǎn)是可以消除因界面網(wǎng)格交錯(cuò)產(chǎn)生的插值運(yùn)算誤差，提高計(jì)算速度和精度。

等離子體唯象學(xué)模型通過在Fluent中加載UDF的形式在動(dòng)量方程中添加體積力源項(xiàng)實(shí)現(xiàn)，傳統(tǒng)的UDF以網(wǎng)格坐標(biāo)的不等式來限定等離子體激勵(lì)區(qū)域，這種方法多用于網(wǎng)格數(shù)量較小時(shí)的流動(dòng)計(jì)算，如對(duì)平板流動(dòng)進(jìn)行激勵(lì)。但對(duì)于網(wǎng)格數(shù)量較大的機(jī)翼流固耦合仿真運(yùn)算，本身運(yùn)算量很大，以不等式形式在外流場(chǎng)通過遍歷網(wǎng)格坐標(biāo)來定位激勵(lì)區(qū)域?qū)⒑馁M(fèi)大量的時(shí)間，因此采用劃分多流體域的形式將等離子體激勵(lì)區(qū)域的流體獨(dú)立劃分，劃分時(shí)可采用自頂向下與自底向上相結(jié)合的方法與布爾運(yùn)算分割流體域。等離子體激勵(lì)流體域與外流體域的交界面類型為interior，優(yōu)點(diǎn)是減小插值運(yùn)算誤差，提高計(jì)算精度。

2 計(jì)算結(jié)果與分析

流體計(jì)算采用壓力遠(yuǎn)場(chǎng)邊界，其中p=1.01×105Pa，T=288.15 K，Ma=0.4，α=3°，湍流模型選擇k-ωSST模型，黏度采用Sutherland動(dòng)力學(xué)黏度關(guān)聯(lián)式進(jìn)行計(jì)算，動(dòng)態(tài)網(wǎng)格采用平滑、分層方法，收斂殘差為10-4，采用二階迎風(fēng)格式，壓力速度耦合方式采用Coupled算法，機(jī)翼表面為無滑移壁面。

在唯象模型中，a=5.773 mm、b=10 mm、d=0.83 mm、?=3 kHz、ρc=1×1017m-3、ec=1.602×10-19、U0=20 kV、Eb=3×103kV/m、t=67 μs，由此可以得到E0=24 MV/m、k1=21 MV/m、k2=36.4 MV/m，代入式(6)中即可求得等離子體激勵(lì)區(qū)域內(nèi)各點(diǎn)的體積力。以上面給定參數(shù)為基礎(chǔ)，通過調(diào)節(jié)式(6)中的變量，可以實(shí)現(xiàn)不同倍數(shù)的激勵(lì)強(qiáng)度。

將機(jī)翼簡(jiǎn)化為具有各向同性材料的固體，激活大變形選項(xiàng)，在翼根施加固定約束，材料密度ρ=10.2 t/m3，楊氏彈性模量為E=200 MPa，泊松比為0.394，剪切模量為G=318.9 MPa。根據(jù)力與位移的L2范數(shù)小于缺省的收斂準(zhǔn)則作為收斂標(biāo)準(zhǔn)。在流固耦合面進(jìn)行力與位移的數(shù)據(jù)傳遞，耦合面變量收斂殘差為10-2。當(dāng)耦合迭代20次時(shí)，流場(chǎng)求解器與結(jié)構(gòu)求解器分別收斂且耦合面殘差收斂，可以認(rèn)為耦合計(jì)算收斂。

在上述條件下，分別對(duì)不施加激勵(lì)與施加不同強(qiáng)度激勵(lì)下的機(jī)翼外流場(chǎng)進(jìn)行計(jì)算。圖5給出了施加激勵(lì)前后機(jī)翼上表面流場(chǎng)的流線圖。

從圖5可以看出，未施加激勵(lì)時(shí)前掠翼分離區(qū)主要位于翼根后部；施加激勵(lì)后，加速的部分流體與未受加速的流體之間產(chǎn)生黏性剪切力，在激勵(lì)區(qū)域的兩側(cè)，受到剪切力作用的流體發(fā)生變形并誘導(dǎo)出具有一定強(qiáng)度的旋渦，由于上翼面流體具有沿展向的分速度，所以旋渦向翼根方向偏轉(zhuǎn)。

圖5 外流場(chǎng)流線圖

圖6給出了施加等離子體激勵(lì)前后機(jī)翼75%展向位置的壓力云圖。

圖6 y/b=75%位置壓力云圖

從圖6可以看出，未施加激勵(lì)時(shí)，機(jī)翼展向75%位置處下表面前緣對(duì)來流起阻礙作用導(dǎo)致壓力較大，上表面前緣由于流體的Coanda效應(yīng)導(dǎo)致壓力較小。施加激勵(lì)后，由于電場(chǎng)力做功對(duì)流體注入能量使總能量增大，壓力勢(shì)能和動(dòng)能均增大，上表面處于激勵(lì)位置的流體壓力增大，而下表面無明顯變化。

圖7給出了50%、75%、97.5%展向位置的速度云圖。

圖7 機(jī)翼外側(cè)速度云圖

從圖7可以看出，受激勵(lì)流體的動(dòng)能增加表現(xiàn)在受激勵(lì)位置的速度有不同程度的增大，其中在激勵(lì)位置內(nèi)側(cè)邊緣速度最大，這是由于此處流體受到電場(chǎng)力、壓差力和誘導(dǎo)的旋渦加速等多重作用。

圖8給出了施加激勵(lì)前后分別在2.5%、25%、50%、75%、97.5%展向位置的壓力系數(shù)變化情況，其中實(shí)線表示未進(jìn)行激勵(lì)，虛線表示進(jìn)行激勵(lì)，橫坐標(biāo)為空間絕對(duì)坐標(biāo)，由于機(jī)翼前掠，坐標(biāo)值從小到大分別表示翼尖到翼根、前緣到后緣的不同位置。

從圖8可以看出，與未施加激勵(lì)相比，施加激勵(lì)后流體受到電場(chǎng)力做功，能量增加，壓力增大，沿展向外側(cè)前緣吸力面的壓力系數(shù)絕對(duì)值先減小后增大，外側(cè)前緣上下表面壓差逐漸減小，且這種減小的程度隨激勵(lì)強(qiáng)度的增大而增大。激勵(lì)后機(jī)翼上表面外側(cè)前緣壓力增大，下表面壓力無明顯變化，由于壓差作用將產(chǎn)生低頭力矩，使翼尖的彈性變形減弱。從圖8中還可以看出，隨著激勵(lì)強(qiáng)度的增大，內(nèi)側(cè)上表面壓力系數(shù)的絕對(duì)值有一定程度的減小，而下表面變化不明顯，導(dǎo)致翼根的升力也有一定的減小。

圖8 壓力系數(shù)變化

圖9給出了不同激勵(lì)強(qiáng)度下機(jī)翼彈軸撓度的變化情況。

圖9 彈軸撓度

從圖9可以看出，在不施加激勵(lì)時(shí)，前掠翼由于氣動(dòng)彈性效應(yīng)使彈軸在翼尖處具有較大的撓度，進(jìn)行等離子體流動(dòng)控制后整個(gè)彈軸的撓度都有所降低，其中翼尖的撓度變化最為明顯，撓度降低了8.5%～23.7%，說明在進(jìn)行流動(dòng)控制后機(jī)翼受到的彎矩減小。

圖10給出了翼尖扭轉(zhuǎn)角與機(jī)翼結(jié)構(gòu)變形最大應(yīng)力隨激勵(lì)強(qiáng)度的變化情況。

圖10 翼尖扭轉(zhuǎn)角與最大應(yīng)力變化

從圖10可以看出，翼尖的扭轉(zhuǎn)角隨激勵(lì)強(qiáng)度的增加降低了8.6%～24.7%，最大應(yīng)力降低了9.8%～29.8%，說明在進(jìn)行流動(dòng)控制后機(jī)翼受到的扭矩減小。

由此可見，在前掠翼前緣外側(cè)上表面處布置等離子體激勵(lì)器后進(jìn)行流動(dòng)控制時(shí)，機(jī)翼彈軸的彎曲變形和翼尖的扭轉(zhuǎn)變形都得到了不同程度的抑制，說明機(jī)翼的受彎扭程度降低，抵抗靜氣彈發(fā)散的能力增強(qiáng)。

圖11給出了在進(jìn)行流動(dòng)控制過程中，升力系數(shù)、阻力系數(shù)、彎矩系數(shù)和扭矩系數(shù)隨激勵(lì)強(qiáng)度增加的變化情況。

圖11 氣動(dòng)特性參數(shù)變化

從圖11可以看出，隨著激勵(lì)強(qiáng)度增強(qiáng)，升力系數(shù)、彎矩、扭矩系數(shù)有不同程度的降低，阻力系數(shù)變化相對(duì)較小。說明將等離子體激勵(lì)器布置在機(jī)翼前緣外側(cè)上表面處，在一定程度上提升機(jī)翼抵抗彎扭變形的能力，對(duì)其靜氣彈發(fā)散問題能起到較好的抑制作用，但同時(shí)以犧牲部分升力等氣動(dòng)特性為代價(jià)，需要進(jìn)一步采取措施進(jìn)行氣動(dòng)補(bǔ)償。

3 結(jié)論

本文針對(duì)前掠翼靜氣彈發(fā)散的主動(dòng)抑制問題開展研究，基于等離子體流動(dòng)控制與流固雙向靜力耦合技術(shù)，通過仿真計(jì)算，得出如下結(jié)論：

1)利用等離子體流動(dòng)控制技術(shù)可以改變彈性前掠翼表面的氣動(dòng)力分布，通過等離子體的激勵(lì)作用，受激勵(lì)的流體從電場(chǎng)中吸收能量，壓力增大，流速加快。

2)在彈性前掠翼前緣外側(cè)上表面布置激勵(lì)器進(jìn)行等離子體流動(dòng)控制可以在一定程度上降低機(jī)翼彈軸的撓度，減小翼尖扭轉(zhuǎn)角，達(dá)到較好抑制靜氣彈發(fā)散的目的，且隨著激勵(lì)強(qiáng)度的增大，抑制效果逐漸增強(qiáng)。

3)將等離子體激勵(lì)器布置在機(jī)翼前緣外側(cè)上表面處，進(jìn)行等離子體流動(dòng)控制對(duì)前掠機(jī)翼彈性變形進(jìn)行抑制的過程中，機(jī)翼受彎曲和扭轉(zhuǎn)的強(qiáng)度減弱，阻力沒有明顯變化，但升力等氣動(dòng)特性有一定程度的損失，且隨著激勵(lì)強(qiáng)度的增大，升力損失逐漸增大，需要進(jìn)一步采取措施進(jìn)行氣動(dòng)補(bǔ)償。