動(dòng)導(dǎo)數(shù)數(shù)值預(yù)測(cè)中的相關(guān)問題

袁先旭， 陳琦， 謝昱飛， 陳堅(jiān)強(qiáng)

中國(guó)空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所， 綿陽(yáng)　621000

?

動(dòng)導(dǎo)數(shù)數(shù)值預(yù)測(cè)中的相關(guān)問題

袁先旭， 陳琦， 謝昱飛*， 陳堅(jiān)強(qiáng)

中國(guó)空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所， 綿陽(yáng)621000

動(dòng)導(dǎo)數(shù)是飛行器動(dòng)態(tài)穩(wěn)定性分析、彈道設(shè)計(jì)和控制系統(tǒng)設(shè)計(jì)的重要參數(shù)，其預(yù)測(cè)方法主要有工程近似方法、數(shù)值模擬和風(fēng)洞試驗(yàn)?；趧?dòng)導(dǎo)數(shù)、交叉導(dǎo)數(shù)的概念，介紹了強(qiáng)迫振蕩法、自由振蕩法等預(yù)測(cè)動(dòng)導(dǎo)數(shù)的數(shù)值模擬方法，重點(diǎn)就超聲速、高超聲速動(dòng)導(dǎo)數(shù)數(shù)值預(yù)測(cè)中存在的問題開展討論，包括時(shí)間步長(zhǎng)、子迭代步數(shù)的選取；振蕩頻率對(duì)辨識(shí)結(jié)果的影響；強(qiáng)迫振蕩法和自由振蕩法辨識(shí)動(dòng)導(dǎo)數(shù)可能存在差異的原因分析；交叉導(dǎo)數(shù)和交叉耦合導(dǎo)數(shù)的辨識(shí)問題以及飛行器構(gòu)型對(duì)動(dòng)導(dǎo)數(shù)預(yù)測(cè)的影響等。并結(jié)合算例進(jìn)行了具體分析，在總結(jié)現(xiàn)有研究經(jīng)驗(yàn)的同時(shí)，針對(duì)當(dāng)前研究中存在的困惑和難題，提出了相應(yīng)的改進(jìn)建議。

動(dòng)導(dǎo)數(shù)； 交叉導(dǎo)數(shù)； 強(qiáng)迫振蕩； 自由振蕩； 振蕩頻率

動(dòng)導(dǎo)數(shù)是飛行器控制系統(tǒng)設(shè)計(jì)中的重要參數(shù)之一，對(duì)飛行器的動(dòng)態(tài)穩(wěn)定性和飛行品質(zhì)有重要影響。隨著現(xiàn)代先進(jìn)飛行器對(duì)機(jī)動(dòng)性和敏捷性的要求越來越高，特別是大迎角機(jī)動(dòng)飛行時(shí)，傳統(tǒng)上視為小量而被忽略的交叉導(dǎo)數(shù)和交叉耦合導(dǎo)數(shù)也逐漸引起人們的重視。總體來看，隨著飛行器機(jī)動(dòng)性能的提升，動(dòng)導(dǎo)數(shù)的重要性也日益凸顯。

在20世紀(jì)90年代以前，動(dòng)導(dǎo)數(shù)的預(yù)測(cè)方法主要是經(jīng)驗(yàn)和半經(jīng)驗(yàn)方法，典型的有修正牛頓理論、內(nèi)伏牛頓流理論、牛頓-玻爾茲曼理論以及修正激波-膨脹波理論等[1]。90年代之后，隨著計(jì)算流體力學(xué)(CFD)技術(shù)和計(jì)算機(jī)硬件水平的迅速發(fā)展，數(shù)值求解Euler/Navier-Stokes方程獲取動(dòng)導(dǎo)數(shù)的方法迅速普及。通過數(shù)值模擬強(qiáng)迫振蕩過程或自由振蕩過程獲取動(dòng)導(dǎo)數(shù)的方法，國(guó)內(nèi)外都開展了廣泛的研究工作，國(guó)外如Weinacht[2]和Qin等[3]；國(guó)內(nèi)如劉偉[4-6]，袁先旭[7-9]等較早開展了動(dòng)導(dǎo)數(shù)的數(shù)值模擬預(yù)測(cè)研究工作。

2010年左右，采用數(shù)值模擬方法預(yù)測(cè)動(dòng)導(dǎo)數(shù)的文獻(xiàn)在國(guó)內(nèi)開始大量涌現(xiàn)，研究外形也不再是單一的鈍頭體，開始向各種復(fù)雜外形發(fā)展，研究?jī)?nèi)容也不斷豐富，促進(jìn)了數(shù)值預(yù)測(cè)動(dòng)導(dǎo)數(shù)方法的發(fā)展。但時(shí)至今日，動(dòng)導(dǎo)數(shù)數(shù)值預(yù)測(cè)中仍然存在很多困難和疑問，需要深入開展相關(guān)的研究工作。本文以討論的形式，在總結(jié)現(xiàn)有工作的同時(shí)，對(duì)當(dāng)前數(shù)值預(yù)測(cè)動(dòng)導(dǎo)數(shù)中存在的問題進(jìn)行剖析。期望通過對(duì)存在問題的討論和研究，達(dá)到共同促進(jìn)數(shù)值預(yù)測(cè)動(dòng)導(dǎo)數(shù)技術(shù)發(fā)展的目的。

1　動(dòng)導(dǎo)數(shù)概念

動(dòng)導(dǎo)數(shù)的概念來自工程設(shè)計(jì)，長(zhǎng)期以來，對(duì)動(dòng)導(dǎo)數(shù)的確切含義仍存在不同看法和爭(zhēng)議，導(dǎo)致在非定常氣動(dòng)力建模方面也不盡一致。

一般認(rèn)為，最早的數(shù)學(xué)模型由Bryan和Williams[10]首先提出，Bryan將氣動(dòng)力、力矩視為擾動(dòng)速度、控制角度和其速率的瞬時(shí)值的函數(shù)，以俯仰為例，即

(1)

式中：Cm為俯仰力矩系數(shù)；α和β分別為迎角和側(cè)滑角；p、q和r分別為偏航、俯仰和滾轉(zhuǎn)角速度。Tobak和Schiff[11]采用指示函數(shù)作為氣動(dòng)力的泛函，建立了非線性指示泛函理論。仍以俯仰為例，其非定常動(dòng)態(tài)俯仰力矩的依賴關(guān)系為一微分-積分系統(tǒng)：

Cm(t)=Cm(t0)+

(2)

式中：L和V∞分別為參考長(zhǎng)度和來流速度；Cm α和Cm q分別為俯仰力矩系數(shù)對(duì)迎角和俯仰角速度的導(dǎo)數(shù)；ξ∈[t0,t]。雖然式(2)在數(shù)學(xué)上是完備的，但非線性指示函數(shù)的確定非常困難，不可能直接求解這一微分-積分系統(tǒng)，必須對(duì)其進(jìn)行簡(jiǎn)化處理。

Etkin和Reid[12]認(rèn)為，非定常氣動(dòng)力、力矩是狀態(tài)變量的泛函。在任何給定的時(shí)刻決定氣動(dòng)力的流場(chǎng)實(shí)際上不僅取決于瞬時(shí)的姿態(tài)，嚴(yán)格來說還和它的整個(gè)過去的歷史有關(guān)，這個(gè)泛函關(guān)系式為

L(t)=L(α(τ))∞≤τ≤t

(3)

任玉新和劉秋生[13]利用Tobak和Schiff的非線性指示函數(shù)的方法，發(fā)展了改進(jìn)的Etkin模型，在基準(zhǔn)狀態(tài)參數(shù)中包含了時(shí)間變量，擴(kuò)展了動(dòng)導(dǎo)數(shù)的應(yīng)用范圍。

這里根據(jù)Etkin模型，給出動(dòng)導(dǎo)數(shù)的概念。對(duì)式(3)，當(dāng)α(τ)可在t附近展開為收斂的泰勒級(jí)數(shù)時(shí)，有

(4)

(5)

(6)

據(jù)此，以俯仰運(yùn)動(dòng)為例，可給出動(dòng)態(tài)俯仰力矩系數(shù)的表達(dá)式：

(7)

取α0為平衡迎角，θ為俯仰角，定義為θ=α-α0。將式(7)在平衡迎角α0處進(jìn)行泰勒展開，可得到

(8)

2　動(dòng)導(dǎo)數(shù)數(shù)值預(yù)測(cè)方法

數(shù)值預(yù)測(cè)動(dòng)導(dǎo)數(shù)，目前常用的方法主要有強(qiáng)迫振蕩法和自由振蕩法兩種。以俯仰方向?yàn)槔?，?qiáng)迫振蕩法一般是給定一個(gè)如下形式的簡(jiǎn)諧振蕩：

α=α1+θ=α1+αmsin(2πft)

(9)

式中：α1為振蕩起始迎角；αm為振蕩幅值;f為振蕩頻率。通過數(shù)值求解飛行器作強(qiáng)迫俯仰振蕩的時(shí)間歷程，即可給出動(dòng)態(tài)俯仰力矩系數(shù)Cm的時(shí)間歷程曲線和Cm-α遲滯圈。根據(jù)動(dòng)態(tài)氣動(dòng)力模型式(8)，略去二階項(xiàng)和高階項(xiàng)，對(duì)上述遲滯圈在一個(gè)周期內(nèi)進(jìn)行積分，即可辨識(shí)得到迎角α1處的動(dòng)導(dǎo)數(shù)，即

(10)

式中：Cm0為起始迎角處的俯仰力矩系數(shù)；T為振蕩周期。自由振蕩法則是放開飛行器的轉(zhuǎn)動(dòng)自由度，飛行器在自身氣動(dòng)力矩的作用下自由振蕩。仍以俯仰為例，單自由度無機(jī)械阻尼俯仰自由振蕩方程可寫為

(11)

式中：I為無量綱化的轉(zhuǎn)動(dòng)慣量，通過數(shù)值計(jì)算得到自由俯仰振蕩的時(shí)間歷程曲線θ=θ(t)和Cm=Cm(t)后，可確定自由振蕩的頻率f=1/T，并選取同一個(gè)周期的兩個(gè)峰值解θ1、θ2，可進(jìn)而得到平衡迎角處的動(dòng)導(dǎo)數(shù)：

(12)

上述方法只能辨識(shí)出平衡迎角處的動(dòng)導(dǎo)數(shù)，應(yīng)用受到一定的限制。袁先旭[14]提出了一種新型振蕩法，可辨識(shí)出自由振蕩經(jīng)過的任一迎角處的動(dòng)導(dǎo)數(shù)，拓寬了該方法的應(yīng)用范圍。

最近幾年，求解周期性非定常流場(chǎng)的諧波法等方法得到了快速發(fā)展，并被應(yīng)用于動(dòng)導(dǎo)數(shù)的預(yù)測(cè)[15-17]。由于不需要求解一個(gè)完整的振蕩周期，計(jì)算效率有大幅提高，但其本質(zhì)仍可歸類為強(qiáng)迫振蕩法，這里不再展開敘述。

上述方法均以俯仰為例，給出俯仰動(dòng)導(dǎo)數(shù)的辨識(shí)過程，偏航和滾轉(zhuǎn)的動(dòng)導(dǎo)數(shù)可類似給出。另外，上述方法得到的是組合導(dǎo)數(shù)，要將組合項(xiàng)分離開，則需加入平動(dòng)運(yùn)動(dòng)的模擬，具體方法可參閱文獻(xiàn)[18]。綜合以上方法不難看出，數(shù)值模擬預(yù)測(cè)動(dòng)導(dǎo)數(shù)的方法，其過程中存在一些共同的問題：首先，需要求解非定常Navier-Stokes方程，一般采用雙時(shí)間步方法，時(shí)間步長(zhǎng)的選取、子迭代步數(shù)以及子迭代收斂判據(jù)等都會(huì)影響辨識(shí)結(jié)果；其次，如果是自由振蕩，則需要給定飛行器的轉(zhuǎn)動(dòng)慣量，而如果是強(qiáng)迫振蕩，則需要給定飛行器的振蕩頻率，不同的振蕩頻率和振蕩幅值也會(huì)影響動(dòng)導(dǎo)數(shù)的辨識(shí)結(jié)果；最后，復(fù)雜構(gòu)型飛行器存在流動(dòng)干擾，也會(huì)影響動(dòng)導(dǎo)數(shù)的預(yù)測(cè)精準(zhǔn)度，而交叉導(dǎo)數(shù)和交叉耦合導(dǎo)數(shù)的辨識(shí)方法目前還存在爭(zhēng)議。對(duì)以上這些問題，本文將結(jié)合算例進(jìn)行具體分析。

3　動(dòng)導(dǎo)數(shù)預(yù)測(cè)中存在的問題與分析

3.1時(shí)間步長(zhǎng)、子迭代步數(shù)的選取

對(duì)飛行力學(xué)方程和Navier-Stokes方程之間的耦合求解問題，通常有松耦合和緊耦合之分。松耦合方法將飛行力學(xué)方程和Navier-Stokes方程分別獨(dú)立求解，在時(shí)間域上交錯(cuò)推進(jìn)，從而獲得耦合系統(tǒng)的響應(yīng)。已經(jīng)證明，該方法在時(shí)間推進(jìn)上只有一階精度，過大的時(shí)間步長(zhǎng)將顯著降低非定常流場(chǎng)求解精度，進(jìn)而降低動(dòng)導(dǎo)數(shù)的預(yù)測(cè)精度。緊耦合方法[19]在亞迭代中實(shí)現(xiàn)兩套系統(tǒng)之間的信息交換，耦合推進(jìn)的時(shí)間精度能達(dá)到二階，因而對(duì)時(shí)間步長(zhǎng)的選取可比松耦合方法稍大。

以松耦合方法為例，分析時(shí)間步長(zhǎng)、子迭代步數(shù)對(duì)動(dòng)導(dǎo)數(shù)預(yù)測(cè)結(jié)果的影響。選取的計(jì)算模型為鈍度比為0.3的鈍錐，其外形和網(wǎng)格如圖1所示。來流馬赫數(shù)Ma∞=6.85，以底部直徑為參考長(zhǎng)度的雷諾數(shù)Re=1.45×106，強(qiáng)迫振蕩的振幅為1°。

圖1　鈍錐外形和計(jì)算網(wǎng)格Fig.1　Geometry and computational grid of blunt cone

表1給出了模擬結(jié)果，無量綱時(shí)間步長(zhǎng)分別取為0.10，0.05，0.20；子迭代步數(shù)分別取為6，15，25。時(shí)間步長(zhǎng)和子迭代步數(shù)的選取對(duì)靜導(dǎo)數(shù)影響較小，但對(duì)動(dòng)導(dǎo)數(shù)的影響很大。同樣的時(shí)間步長(zhǎng)，隨著子迭代步數(shù)的增加，動(dòng)導(dǎo)數(shù)預(yù)測(cè)結(jié)果有向試驗(yàn)數(shù)據(jù)逼近的趨勢(shì)；同樣的子迭代步數(shù)，隨著時(shí)間步長(zhǎng)的降低，動(dòng)導(dǎo)數(shù)預(yù)測(cè)結(jié)果同樣有向試驗(yàn)數(shù)據(jù)逼近的趨勢(shì)。

表1的計(jì)算結(jié)果似乎表明，時(shí)間步長(zhǎng)越小越好、子迭代步數(shù)越多越好；而當(dāng)時(shí)間步長(zhǎng)、子迭代步數(shù)滿足了一定的收斂性要求后，動(dòng)導(dǎo)數(shù)的預(yù)測(cè)結(jié)果與兩者的選取是無關(guān)的。因此，認(rèn)為要保證動(dòng)導(dǎo)數(shù)的預(yù)測(cè)精準(zhǔn)度，時(shí)間步長(zhǎng)和子迭代步數(shù)必須滿足亞迭代收斂性的要求，以降低數(shù)值誤差。建議在計(jì)算過程中，檢驗(yàn)并給出亞迭代收斂曲線。

表1時(shí)間步長(zhǎng)、子迭代步數(shù)等對(duì)預(yù)測(cè)結(jié)果的影響

Table 1Effect of time step, step number of sub-iteration, etc. on prediction result

TimestepSubiterationstepnumberStaticderivativesDynamicderivatives0.106-0.058454-0.48250215-0.057769-0.27530325-0.057594-0.2009810.056-0.057845-0.30687815-0.057663-0.14117525-0.057607-0.1322060.026-0.057630-0.12471115-0.057576-0.11344825-0.057566-0.107140Experiment[18]-0.1056±0.0158

3.2振蕩頻率的影響

飛行器在實(shí)際飛行時(shí)，燃料的持續(xù)消耗將會(huì)導(dǎo)致飛行器的轉(zhuǎn)動(dòng)慣量發(fā)生改變；飛行高度的變化則會(huì)改變來流密度；在陣風(fēng)等擾動(dòng)的持續(xù)作用下，飛行器的受迫振蕩等，這些因素都會(huì)導(dǎo)致實(shí)際飛行時(shí)飛行器的振蕩頻率不斷地變化，因而需要考察動(dòng)導(dǎo)數(shù)隨頻率的變化關(guān)系。

針對(duì)超聲速帶翼導(dǎo)彈標(biāo)模外形，研究了動(dòng)導(dǎo)數(shù)預(yù)測(cè)結(jié)果隨頻率的變化問題。計(jì)算模型和網(wǎng)格如圖2所示。

圖2　帶翼導(dǎo)彈標(biāo)模外形和網(wǎng)格Fig.2　Basic finner geometry and grid

采用諧波平衡法和雙時(shí)間步方法預(yù)測(cè)了飛行器的動(dòng)導(dǎo)數(shù)，計(jì)算結(jié)果如圖3所示。其中“DTS”為雙時(shí)間步方法的計(jì)算結(jié)果，“Exp”為試驗(yàn)結(jié)果[20]，“1 Harmonic”表示諧波數(shù)取1時(shí)，諧波平衡點(diǎn)的動(dòng)導(dǎo)數(shù)預(yù)測(cè)結(jié)果。圖3的計(jì)算結(jié)果表明，一般情況下，兩種方法的預(yù)測(cè)結(jié)果均與試驗(yàn)值吻合良好。

圖3　動(dòng)導(dǎo)數(shù)預(yù)測(cè)結(jié)果與試驗(yàn)值比較Fig.3　Comparison between dynamic derivative prediction results and experimental data

圖4則是不同頻率時(shí)的動(dòng)導(dǎo)數(shù)預(yù)測(cè)結(jié)果，圖中k為無量綱減縮頻率。從預(yù)測(cè)結(jié)果來看，對(duì)超聲速飛行器，在頻率較大時(shí)，動(dòng)導(dǎo)數(shù)對(duì)頻率的變化不敏感，這符合以往的經(jīng)驗(yàn)認(rèn)識(shí)。在實(shí)際應(yīng)用時(shí)，也常常選取較大的振蕩頻率，以提高動(dòng)導(dǎo)數(shù)的預(yù)測(cè)效率。但在低頻時(shí)，動(dòng)導(dǎo)數(shù)隨振蕩頻率的降低迅速發(fā)生改變。若飛行器的固有振蕩頻率較低，而計(jì)算時(shí)選取的頻率較高，則預(yù)測(cè)結(jié)果可能會(huì)背離真實(shí)情況。

圖4　頻率對(duì)動(dòng)導(dǎo)數(shù)預(yù)測(cè)結(jié)果的影響　　Fig.4　Effect of frequency on dynamic derivativeprediction results

孫濤等[21]針對(duì)相同的外形，研究了振蕩頻率f的影響，也給出了相似的研究結(jié)論(見圖5)。

而對(duì)類航天飛機(jī)外形的飛行器，在研究其俯仰動(dòng)導(dǎo)數(shù)時(shí)發(fā)現(xiàn)(見圖6)：在低空低馬赫數(shù)時(shí)，不同頻率的預(yù)測(cè)結(jié)果差別不大；但高空高馬赫數(shù)時(shí)，不同頻率的結(jié)果差別明顯。

圖5　頻率對(duì)動(dòng)導(dǎo)數(shù)預(yù)測(cè)結(jié)果的影響[21]Fig.5　Effect of frequency on dynamic derivativeprediction results[21]

圖6　不同頻率時(shí)動(dòng)導(dǎo)數(shù)預(yù)測(cè)結(jié)果隨馬赫數(shù)的變化Fig.6　Dynamic derivative prediction results versusMach number at different frequencies

在以往的研究任務(wù)中，基于“動(dòng)導(dǎo)數(shù)隨頻率變化不大”的經(jīng)驗(yàn)認(rèn)識(shí)，或者在不清楚飛行器振蕩頻率時(shí)，通常會(huì)選取較大的振蕩頻率，一般給定f=2～10 Hz，以提高非定常計(jì)算效率。但理論分析[14]和實(shí)踐表明，振蕩頻率不能隨意給定。建議結(jié)合飛行器的典型飛行狀態(tài)、質(zhì)量特性和靜導(dǎo)數(shù)，估算給出典型振蕩頻率。

3.3自由振蕩與強(qiáng)迫振蕩的預(yù)測(cè)結(jié)果差異

通常情況下，自由振蕩方法和強(qiáng)迫振蕩方法的動(dòng)導(dǎo)數(shù)預(yù)測(cè)結(jié)果互為補(bǔ)充，可以相互驗(yàn)證。以構(gòu)型較為簡(jiǎn)單平頭雙錐帶翼飛行器外形(見圖7)為例，采用強(qiáng)迫振蕩法獲取了不同馬赫數(shù)時(shí)的動(dòng)導(dǎo)數(shù)，如圖8所示。

圖7　平頭雙錐外形示意圖　　Fig.7　Schematic diagram of flat-nose wingeddouble-cone body

圖8　平頭雙錐外形不同馬赫數(shù)時(shí)的動(dòng)導(dǎo)數(shù)Fig.8　Dynamic derivatives at different Mach numbers for flat-nose winged double-cone body

從模擬結(jié)果來看，隨著馬赫數(shù)的增加，俯仰動(dòng)導(dǎo)數(shù)的量值逐漸減小，表明飛行器受到的氣動(dòng)阻尼不斷降低；最后在Ma=7.0時(shí)，動(dòng)導(dǎo)數(shù)的符號(hào)也發(fā)生了變化(α=0°)，此時(shí)空氣對(duì)飛行器做正功，飛行器的自由振蕩將可能發(fā)散。

圖9則通過模擬自由振蕩的過程，得到了不同馬赫數(shù)時(shí)俯仰角的時(shí)間歷程曲線，起始迎角α1分別為1° 和4°?？梢钥吹剑婉R赫數(shù)時(shí)，俯仰角是收斂的，但隨著馬赫數(shù)增大，氣動(dòng)阻尼降低，俯仰角收斂的速度越來越慢；在Ma=7.0時(shí)，動(dòng)導(dǎo)數(shù)為正值，此時(shí)俯仰角的振蕩幅值不斷增加。自由振蕩結(jié)果驗(yàn)證了圖8的強(qiáng)迫振蕩動(dòng)導(dǎo)數(shù)預(yù)測(cè)結(jié)果。

圖9　自由振蕩模擬得到的俯仰角時(shí)間歷程曲線Fig.9　Time history curves of pitch angle from free oscillation simulation

但是對(duì)復(fù)雜構(gòu)型的升力體外形(見圖10)，圖11 給出了兩種方法預(yù)測(cè)的動(dòng)導(dǎo)數(shù)結(jié)果比較。同樣的振蕩頻率下，兩種方法得到的結(jié)果差別較大，在迎角大于10° 以后，兩種方法給出的動(dòng)導(dǎo)數(shù)符號(hào)都不一致。在另外一些飛行器的計(jì)算中，也發(fā)現(xiàn)相似的問題。

POD數(shù)據(jù)后處理應(yīng)用程序包由數(shù)據(jù)輸入模塊、POD數(shù)據(jù)處理核心模塊和數(shù)據(jù)輸出模塊等3個(gè)功能模塊組成(見圖1)。

從我們的計(jì)算經(jīng)驗(yàn)來看，自由振蕩法預(yù)測(cè)動(dòng)導(dǎo)數(shù)時(shí)，由于需要人工選取同一個(gè)周期的兩個(gè)峰值解θ1、θ2，選取的周期不同，特別是當(dāng)θ1、θ2與平衡迎角距離較遠(yuǎn)時(shí)，不完全滿足動(dòng)導(dǎo)數(shù)鄰域的概念，對(duì)動(dòng)導(dǎo)數(shù)的辨識(shí)結(jié)果有較大影響。因此，建議盡量采用小振幅強(qiáng)迫振蕩法預(yù)測(cè)動(dòng)導(dǎo)數(shù)，以典型狀態(tài)的自激振蕩結(jié)果作為驗(yàn)證。

圖10　升力體外形Fig.10　Geometry of lifting body

圖11　升力體外形強(qiáng)迫振蕩和自由振蕩方法辨識(shí)的動(dòng)導(dǎo)數(shù)的比較Fig.11　Comparison of dynamic derivatives for lifting body using forced and free oscillation methods

3.4交叉導(dǎo)數(shù)和交叉耦合導(dǎo)數(shù)的預(yù)測(cè)

以往，交叉導(dǎo)數(shù)和交叉耦合導(dǎo)數(shù)常常被認(rèn)為是小量而被忽略，沒有得到足夠的重視。近年來，隨著飛行器對(duì)機(jī)動(dòng)性和敏捷性的要求越來越高，特別是在大迎角機(jī)動(dòng)飛行時(shí)，交叉導(dǎo)數(shù)和交叉耦合導(dǎo)數(shù)的重要性日益凸顯。

圖12　交叉耦合導(dǎo)數(shù)隨迎角的變化情況[22]Fig.12　Variation of cross coupling derivatives with angle of attack[22]

目前，交叉導(dǎo)數(shù)和交叉耦合導(dǎo)數(shù)仍沿用單自由度強(qiáng)迫振蕩法，可能不能反應(yīng)飛行器運(yùn)動(dòng)之間的耦合效應(yīng)。建議結(jié)合控制系統(tǒng)設(shè)計(jì)要求設(shè)計(jì)典型耦合運(yùn)動(dòng)形式，并基于該運(yùn)動(dòng)形式的氣動(dòng)力矩時(shí)間歷程曲線辨識(shí)交叉導(dǎo)數(shù)和交叉耦合導(dǎo)數(shù)。

不同外形特征的飛行器，其動(dòng)態(tài)特性相差較大，特別是對(duì)內(nèi)外流一體化飛行器外形，由于振蕩時(shí)存在內(nèi)外流干擾，其動(dòng)態(tài)特性比常規(guī)飛行器更為復(fù)雜。此外，為確保進(jìn)氣道的正常啟動(dòng)，內(nèi)外流一體化飛行器對(duì)姿態(tài)控制的精度要求更高，故其操縱性穩(wěn)定性分析評(píng)估也更加重要，需要準(zhǔn)確預(yù)測(cè)動(dòng)導(dǎo)數(shù)。

對(duì)研究的某環(huán)形進(jìn)氣道外形的飛行器(圖13是進(jìn)氣道入口附近的網(wǎng)格)，在研究其俯仰動(dòng)導(dǎo)數(shù)時(shí)發(fā)現(xiàn)：對(duì)常規(guī)純外流飛行器通常在計(jì)算啟動(dòng)后一個(gè)振蕩周期內(nèi)即可進(jìn)入遲滯圈，遲滯圈重復(fù)性較好；而內(nèi)外流一體化飛行器的相軌線一般經(jīng)過劇烈振蕩才能形成遲滯圈，且遲滯圈重復(fù)性較差，如圖14所示。

對(duì)另一型內(nèi)外流一體化飛行器，考慮湍流模型和不考慮湍流模型的計(jì)算，兩種動(dòng)導(dǎo)數(shù)的預(yù)測(cè)結(jié)果甚至出現(xiàn)了反號(hào)(見圖15)。這些現(xiàn)象都表明，內(nèi)外流一體化飛行器振蕩時(shí)，由于內(nèi)外流相互作用和相互干擾，導(dǎo)致其動(dòng)態(tài)特性也更為復(fù)雜，需要更深入地開展相關(guān)的研究工作。

圖13　進(jìn)氣道入口附近形狀和網(wǎng)格Fig.13　Geometry and grid near inlet

圖14　內(nèi)外流一體化飛行器遲滯圈Fig.14　Hysteresis loop for integrated internal-externalflow vehicle

圖15　考慮和不考慮湍流模型時(shí)計(jì)算的動(dòng)導(dǎo)數(shù)Fig.15　Comparison of dynamic derivatives with/without turbulence model

建議對(duì)于此類問題，應(yīng)采用多種模型、多種方法進(jìn)行對(duì)比分析，以保證動(dòng)導(dǎo)數(shù)預(yù)測(cè)的精準(zhǔn)度。

4　結(jié)　論

針對(duì)動(dòng)導(dǎo)數(shù)數(shù)值預(yù)測(cè)中存在的問題、困難開展討論，并結(jié)合具體算例進(jìn)行分析?？傮w來看，出現(xiàn)問題的原因可大致歸為兩類：

1) 方法問題。包括數(shù)值計(jì)算方法和動(dòng)導(dǎo)數(shù)辨識(shí)方法等。非定常計(jì)算的子迭代收斂判據(jù)、時(shí)間步長(zhǎng)的選取、計(jì)算模型與計(jì)算方法的選取和計(jì)算網(wǎng)格設(shè)計(jì)等，都會(huì)對(duì)動(dòng)導(dǎo)數(shù)預(yù)測(cè)結(jié)果產(chǎn)生影響。

2) 物理機(jī)制。包括動(dòng)導(dǎo)數(shù)、交叉導(dǎo)數(shù)的概念和適用范圍等。動(dòng)導(dǎo)數(shù)的概念本質(zhì)上屬于線化小擾動(dòng)理論范疇，在氣動(dòng)力非定常、非線性效應(yīng)很強(qiáng)的時(shí)候，動(dòng)導(dǎo)數(shù)的概念是否適用以及如何使用等問題仍是未解的難題。

3) 針對(duì)數(shù)值預(yù)測(cè)動(dòng)導(dǎo)數(shù)中存在的各種問題，結(jié)合我們自己的研究經(jīng)驗(yàn)，提出了一些意見和建議，供參考。

[1]劉緒. 高超聲速內(nèi)外流一體化飛行器動(dòng)態(tài)特性研究[D]. 長(zhǎng)沙: 國(guó)防科學(xué)技術(shù)大學(xué), 2011: 1-3.

LIU X. Investigation of dynamic characteristics of hypersonic airframe/propulsion integrative vehicle[D]. Changsha: National University of Defense Technology, 2011: 1-3 (in Chinese).

[2]WEINACHT P. Navier-Stokes predictions of pitch damping for a family of flared projectiles: AIAA-1991-3339[R]. Reston: AIAA, 1991.

[3]QIN N, LUDLOW D K, SHAW S T, et al. Calculation of pitch damping coefficients for projectiles: AIAA-1997-0405[R]. Reston: AIAA, 1997.

[4]劉偉, 張魯民. 鈍體俯仰阻尼導(dǎo)數(shù)數(shù)值計(jì)算[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 1997, 15(4): 427-435.

LIU W, ZHANG L M. Numerical calculation of damping in pitch derivatives of blunt cone[J]. Acta Aerodynamics Sinica, 1997, 15(4): 427-435 (in Chinese).

[5]劉偉, 牟斌. 高超聲速滾轉(zhuǎn)阻尼導(dǎo)數(shù)數(shù)值模擬[J]. 飛行力學(xué), 2000, 18(2): 27-29.

LIU W, MOU B. Numerical simulation of damping-in-roll derivatives of blunt cone for hypersonic flow[J]. Flight Dynamics, 2000, 18(2): 27-29 (in Chinese).

[6]劉偉, 劉君, 柳軍. 平衡氣體效應(yīng)對(duì)飛行器動(dòng)態(tài)特性的影響研究[J]. 飛行力學(xué), 2004, 22(4): 65-68.

LIU W, LIU J, LIU J. Investigation of equilibrium gas effect on dynamic characteristic of aerocraft[J]. Flight Dynamics, 2004, 22(4): 65-68 (in Chinese).

[7]袁先旭, 張涵信, 謝昱飛. 基于CFD 方法的俯仰靜、動(dòng)導(dǎo)數(shù)數(shù)值計(jì)算[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 2005, 23(4): 458-463.

YUAN X X, ZHANG H X, XIE Y F. The pitching static/dynamic derivatives computation based on CFD methods[J]. Acta Aerodynamics Sinica, 2005, 23(4): 458-463 (in Chinese).

[8]袁先旭, 陳堅(jiān)強(qiáng), 王文正. 平頭增阻再入體俯仰動(dòng)態(tài)特性計(jì)算與流動(dòng)機(jī)理分析[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 2007, 25(3): 300-305.

YUAN X X, CHEN J Q, WANG W Z. Pitching dynamic stability computation for plane nose reentry vehicle and flow mechanism analysis[J]. Acta Aerodynamics Sinica, 2007, 25(3): 300-305 (in Chinese).

[9]袁先旭, 張涵信, 謝昱飛. 飛船返回艙再入俯仰動(dòng)穩(wěn)定吸引子數(shù)值仿真[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 2007, 25(4): 431-436.

YUAN X X, ZHANG H X, XIE Y F. Numerical simulation for dynamic stability in pitching of unfinned reentry capsule and bifurcation with mach number prediction[J]. Acta Aerodynamics Sinica, 2007, 25(4): 431-436 (in Chinese).

[10]BRYAN G, WILLIAMS W. The longitudinal stability of aerial gliders[J]. Proceedings of the Royal Society of London, 1904, 73(488-496): 100-116.

[11]TOBAK M, SCHIFF L B. On the formulation of the aerodynamic characteristics in aircraft dynamics: NASA technical report, N76-15082/OSL[R]. Washington, D.C.: NASA, 1976.

[12]ETKIN B, TEICHMANN T. Dynamics of flight: Stability and control[J]. Journal of Guidance Control & Dynamics, 2009, 20(4): 839-840.

[13]任玉新, 劉秋生. 飛行器動(dòng)態(tài)穩(wěn)定性參數(shù)的數(shù)值計(jì)算方法[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 1996, 14(2): 117-126.

REN Y X, LIU Q S. A numerical method for evaluating aerodynamic stability parameters of vehicles[J]. Acta Aerodynamica Sinica, 1996, 14(2): 117-126 (in Chinese).

[14]袁先旭. 非定常流動(dòng)數(shù)值模擬及飛行器動(dòng)態(tài)特性分析研究[D]. 綿陽(yáng)： 中國(guó)空氣動(dòng)力研究與發(fā)展中心, 2002: 111-115.

YUAN X X. Numerical simulation for unsteady flows and research on dynamic characteristics of vehicle[D]. Mianyang: China Aerodynamics Research and Development Center, 2002: 111-115 (in Chinese).

[15]陳堅(jiān)強(qiáng), 陳琦, 袁先旭, 等. 舵面操縱動(dòng)態(tài)響應(yīng)的數(shù)值模擬研究[J]. 力學(xué)學(xué)報(bào), 2013, 45(2): 302-306.

CHEN J Q, CHEN Q, YUAN X X, et al. Numerical simulation study on dynamics response under rudder control[J]. Chinese Journal of Theoretical and Applied Mechanics, 2013, 45(2): 302-306 (in Chinese).

[16]陳琦, 陳堅(jiān)強(qiáng), 袁先旭, 等. 諧波平衡法在動(dòng)導(dǎo)數(shù)快速預(yù)測(cè)中的應(yīng)用研究[J]. 力學(xué)學(xué)報(bào), 2014, 46(2): 183-190.

CHEN Q, CHEN J Q, YUAN X X， et al. Application of a harmonic balance method in rapid predictions of dynamic stability derivatives[J]. Chinese Journal of Theoretical and Applied Mechanics, 2014, 46(2): 183-190 (in Chinese).

[17]楊小權(quán), 程蘇堃, 楊愛明, 等. 基于時(shí)間譜方法的振蕩翼型和機(jī)翼非定常黏性繞流數(shù)值模擬[J]. 航空學(xué)報(bào), 2013, 34(4): 787-797.

YANG X Q, CHENG S K, YANG A M, et al. Time spectral method for numerical simulation of unsteady viscous flow over oscillating airfoil and wing[J]. Acta Aeronautica et Astronautica Sinica, 2013, 34(4): 787-797 (in Chinese).

[18]劉偉, 楊小亮, 趙云飛. 高超聲速飛行器加速度導(dǎo)數(shù)數(shù)值模擬[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 2010, 28(4): 426-429.

LIU W, YANG X L, ZHAO Y F. Numerical simulation of acceleration derivative of hypersonic aircraft[J]. Acta Aerodynamics Sinica, 2010, 28(4): 426-429 (in Chinese).

[19]楊云軍, 崔爾杰, 周偉江. 細(xì)長(zhǎng)三角翼滾轉(zhuǎn)/側(cè)滑耦合運(yùn)動(dòng)的數(shù)值研究[J]. 航空學(xué)報(bào), 2007, 28(1): 14-19.

YANG Y J, CUI E J, ZHOU W J. Numerical research on roll and sideslip coupling motions about a slender delta-wing[J]. Acta Aeronautica et Astronautica Sinica, 2007, 28(1): 14-19 (in Chinese).

[20]EAST R A, HUTT G R. Comparison of predictions and experimental data for hypersonic pitching motion stability[J]. Journal of Spacecraft and Rockets, 1988, 25(3): 225-233.

[21]孫濤, 高正紅, 黃江濤. 基于CFD的動(dòng)導(dǎo)數(shù)計(jì)算與減縮頻率影響分析[J]. 飛行力學(xué), 2011, 29(4): 15-18.

SUN T, GAO Z H, HUANG J T. Identify of aircraft dynamic derivatives based on CFD technology and analysis of reduce frequency[J]. Flight Dynamics, 2011, 29(4):15-18 (in Chinese).

[22]COULTER S M, MARQUART E J. Cross and cross-coupling derivative measurements on the standard dynamics model at AEDC: AIAA-1982-0596[R]. Reston: AIAA, 1982.

[23]BUTLER R W. Aircraft motion sensitivity to cross and cross-coupling damping derivatives: NASA technical report, AD-A032 654/6SL[R]. Washington, D.C.: NASA, 1977.

[24]BUTLER R W, LANGHAM T F. Sensitivity of aircraft spinning motion to dynamic cross-coupling and acceleration derivatives: NASA technical report, AD-A060 516/2SL[R]. Washington, D.C.: NASA, 1979.

[25]LANGHAM T F. Aircraft motion sensitivity to dynamic stability derivatives: AIAA-1979-1621[R]. Reston: AIAA, 1979.

[26]何植岱. 飛機(jī)大迎角穩(wěn)定性對(duì)于交叉導(dǎo)數(shù)和加速度導(dǎo)數(shù)的靈敏度分析[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 1987, 5(4): 334-344.

HE Z D. Sensitivity of aircraft stability to cross-coupling derivatives and angular acceleration derivatives at high angles of attack[J]. Acta Aerodynamics Sinica, 1987, 5(4): 334-344 (in Chinese).

袁先旭男， 博士， 研究員。主要研究方向： 非定常流動(dòng)數(shù)值模擬、 飛行器動(dòng)態(tài)特性分析、 復(fù)雜流動(dòng)數(shù)值計(jì)算方法等。

Tel: 0816-2463168

E-mail: yuanxianxu@cardc.cn

陳琦男， 博士， 助理研究員。主要研究方向： 非定常流動(dòng)數(shù)值模擬、 飛行動(dòng)態(tài)特性分析、 計(jì)算軟件開發(fā)等。

Tel: 0816-2463304

E-mail: chenqi@mail.ustc.edu.cn

謝昱飛男， 博士， 副研究員。主要研究方向： 非定常流動(dòng)數(shù)值模擬、 飛行動(dòng)態(tài)特性分析等。

Tel: 0816-2463090

E-mail: xyf_5843@qq.com

Problems in numerical prediction of dynamic stability derivatives

YUAN Xianxu, CHEN Qi, XIE Yufei*, CHEN Jianqiang

Computational Aerodynamics Institute, China Aerodynamics Research and Development Center,Mianyang621000, China

Dynamic stability derivative is one of the important parameters in the design of control system, the design of flight orbit, and the analysis of dynamic stability of vehicles. The prediction methods of the dynamic stability derivatives mainly include the engineering approximations, numerical simulations and wind tunnel tests. The concepts of dynamic stability derivative and cross derivative are reviewed, and the forced oscillation and free oscillation methods for numerically predicting the dynamic stability derivatives are introduced. The problems existing in the numerical prediction are specially discussed, including the choice of time-marching step and step number of sub-iteration, the effects of the oscillating frequency on the identification results, the analysis of the probable reasons which lead to the difference between the forced oscillation and free oscillation methods in dynamic derivative identification, the identification of the cross and cross coupling derivatives, and the identification of dynamic stability derivatives of the integrated internal-external flow vehicles. The discussions are made with computational cases and detailed analysis, in which both experiences and puzzles exist, and the corresponding suggestions for improvement are presented.

dynamic derivative; cross derivative; forced oscillation; free oscillation; oscillating frequency

2016-04-21； Revised： 2016-04-29； Accepted： 2016-06-04； Published online： 2016-06-1209:01

s： National Natural Science Foundation of China (11172315, 11372341, 11532016)

. Tel.： 0816-2463090E-mail: xyf_5843@qq.com

2016-04-21； 退修日期： 2016-04-29； 錄用日期： 2016-06-04；

時(shí)間： 2016-06-1209:01

www.cnki.net/kcms/detail/11.1929.V.20160612.0901.002.html

國(guó)家自然科學(xué)基金 (11172315, 11372341, 11532016)

.Tel.： 0816-2463090E-mail: xyf_5843@qq.com

10.7527/S1000-6893.2016.0180

V211.3

A

1000-6893(2016)08-2385-10

引用格式： 袁先旭， 陳琦， 謝昱飛， 等． 動(dòng)導(dǎo)數(shù)數(shù)值預(yù)測(cè)中的相關(guān)問題[J]． 航空學(xué)報(bào), 2016, 37(8): 2385-2394. YUAN X X, CHEN Q, XIE Y F, et al. Problems in numerical prediction of dynamic stability derivatives[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(8): 2385-2394.

http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

URL： www.cnki.net/kcms/detail/11.1929.V.20160612.0901.002.html