γ－Reθ轉(zhuǎn)捩模型的標(biāo)定與應(yīng)用

牟 斌，江 雄，肖中云，陳作斌

（中國空氣動(dòng)力研究與發(fā)展中心，四川 綿陽 621000）

0 引 言

轉(zhuǎn)捩邊界層流動(dòng)研究在CFD工程應(yīng)用的許多方面非常重要，如層流翼型、風(fēng)輪機(jī)、輪船船體、渦輪機(jī)葉柵等流動(dòng)特性研究方面。這些流動(dòng)的雷諾數(shù)范圍確定了其層流區(qū)域和湍流區(qū)域往往具有相同量級(jí)，全湍流計(jì)算在阻力等方面帶來的偏差不可忽視，必須考慮轉(zhuǎn)捩的影響。

在過去的幾十年中，全湍流的計(jì)算取得了輝煌的成就，并發(fā)展起以并行計(jì)算、非結(jié)構(gòu)網(wǎng)格為代表的現(xiàn)代CFD技術(shù)；而轉(zhuǎn)捩的數(shù)值模擬則嚴(yán)重滯后。Menter［1］將其歸結(jié)為兩方面的原因：一、轉(zhuǎn)捩機(jī)制復(fù)雜，包含了自然轉(zhuǎn)捩、bypass轉(zhuǎn)捩、分離誘導(dǎo)轉(zhuǎn)捩，同時(shí)湍流邊界層在順壓梯度下可能再層流化，導(dǎo)致很難用一個(gè)簡單的數(shù)學(xué)模型描述如此多的效應(yīng)。二、轉(zhuǎn)捩流動(dòng)中線性和非線性效應(yīng)相關(guān)，RANS平均消除了線性擾動(dòng)的增長，因此傳統(tǒng)的雷諾平均方程不能描述轉(zhuǎn)捩流動(dòng)。

在工程應(yīng)用中，考慮轉(zhuǎn)捩效應(yīng)的經(jīng)典方法是經(jīng)驗(yàn)關(guān)聯(lián)，就是把自由流的湍流度和當(dāng)?shù)氐膲毫μ荻鹊韧ㄟ^實(shí)驗(yàn)數(shù)據(jù)關(guān)聯(lián)到轉(zhuǎn)捩動(dòng)量厚度雷諾數(shù)，如Abu－Ghannam和Shaw的關(guān)聯(lián)等。遺憾的是以前的關(guān)聯(lián)往往是非局部變量的函數(shù)，需要作諸如邊界層厚度判定，駐點(diǎn)、分離點(diǎn)查找，沿流線積分等，這些算法受到外形、網(wǎng)格拓?fù)浣Y(jié)構(gòu)的制約，應(yīng)用的范圍很窄。

為了與現(xiàn)代 CFD 兼容，Menter與 Langtry［1－5］設(shè)計(jì)了γ－Reθ轉(zhuǎn)捩模型。該模型采用經(jīng)驗(yàn)關(guān)聯(lián)的方法，并不模擬轉(zhuǎn)捩的物理機(jī)理。其成功之處在于用應(yīng)變率雷諾數(shù)觸發(fā)轉(zhuǎn)捩，并通過求解額外的兩個(gè)輸運(yùn)方程避免了非局部變量的采用，從而形成與現(xiàn)代CFD相容的框架。在 Menter［2－5］早期的文章中，由于商業(yè)保密的原因，兩個(gè)重要的經(jīng)驗(yàn)關(guān)聯(lián)函數(shù)沒有給出，國外許多學(xué)者為了應(yīng)用，依據(jù)平板實(shí)驗(yàn)結(jié)果提出了各種in－h(huán)ouse公式［6－8］，并在各種 CFD 工程應(yīng)用中獲得成功；國內(nèi)空氣動(dòng)力中心張玉倫［13－14］，西工大張曉東［15］等也提出了各自標(biāo)定的經(jīng)驗(yàn)關(guān)系式。雖然2009年Menter及Langtry［1］給出了一組關(guān)聯(lián)函數(shù)，然而，轉(zhuǎn)捩模型對不同的代碼及標(biāo)定過程中的微小差異非常敏感，因此這一組關(guān)聯(lián)函數(shù)在CFX中能夠得到很好的結(jié)果，在另一個(gè)解算器中的表現(xiàn)可能就會(huì)差一些［9］。為了較好模擬轉(zhuǎn)捩效應(yīng)，有必要在代碼中實(shí)現(xiàn)轉(zhuǎn)捩模型并標(biāo)定關(guān)聯(lián)函數(shù)。

本文首先在內(nèi)部軟件mbns3d上實(shí)現(xiàn)了轉(zhuǎn)捩模型，依據(jù)零壓力梯度平板轉(zhuǎn)捩實(shí)驗(yàn)數(shù)據(jù)確定了一組關(guān)聯(lián)函數(shù)，然后針對A翼型做了驗(yàn)證計(jì)算。隨后數(shù)值模擬了NLR7301兩段翼型，研究了轉(zhuǎn)捩、遠(yuǎn)場邊界條件及遠(yuǎn)場范圍對精確計(jì)算阻力的影響，獲得了有益結(jié)論。

1 γ－Reθ 轉(zhuǎn)捩模型簡介

（1）控制方程

無量綱的γ－Reθ輸運(yùn)方程的守恒形式為：

式中各項(xiàng)的詳細(xì)解釋可見文獻(xiàn)［5］。

（2）關(guān)聯(lián)函數(shù)

模型中存在三個(gè)關(guān)聯(lián)函數(shù)：

對于關(guān)聯(lián)函數(shù)Reθt，不同文獻(xiàn)中的形式略有不同，本文采用Langtry［5］建議的形式。下面的工作是確定其余兩個(gè)關(guān)聯(lián)函數(shù)。

2 γ－Reθ 關(guān)聯(lián)函數(shù)的確定

標(biāo)定工作在內(nèi)部軟件mbns3d軟件上進(jìn)行，mbns3d采用多塊結(jié)構(gòu)網(wǎng)格（含拼接、對接、重疊），有限體積求解，具備多重網(wǎng)格、預(yù)處理、并行計(jì)算功能。

標(biāo)定計(jì)算中，對流項(xiàng)離散選用三階精度的Roe's FDS格式，粘性項(xiàng)二階中心格式離散，AFADI方法求解離散方程。為了適應(yīng)低速情況和加速收斂，采用了預(yù)處理技術(shù)和多重網(wǎng)格技術(shù)。

標(biāo)定依據(jù)文獻(xiàn)［5］提供的 T3A，T3B，T3A－和Schubauer＆Klebanof等平板試驗(yàn)數(shù)據(jù)，試驗(yàn)條件見表1。

表1 平板轉(zhuǎn)捩試驗(yàn)的進(jìn)口條件Table 1 Inlet condition for the flat plate test cases

平板長度選為2．0，來流遠(yuǎn)場邊界距平板前緣0．22，法向遠(yuǎn)場邊界距壁面0．22，網(wǎng)格為H型，壁面第一層網(wǎng)格法向距離1．0×10－5，平板前緣第一網(wǎng)格的x向距離1．5×10－3。網(wǎng)格規(guī)模為273×97，其中壁面上有241網(wǎng)格點(diǎn)，滿足四重多重網(wǎng)格計(jì)算要求。

表1中的FSTI為是平板前緣的湍流強(qiáng)度，由于網(wǎng)格外邊界距平板前緣0．22m，因此外邊界的湍流強(qiáng)度實(shí)際上是依據(jù)下式估計(jì)［5］：

其中β＝0．09，β＊＝0．0828，x是下游距進(jìn)口的流向距離，V是平均速度。

關(guān)聯(lián)函數(shù)中Reθc和Flength是耦合性很強(qiáng)的函數(shù)，一個(gè)確定，另一個(gè)也隨之確定。為了簡便，這里先選定Reθc函數(shù)。其選取有一定的原則，保證轉(zhuǎn)捩能夠發(fā)生；在低雷諾數(shù)轉(zhuǎn)捩區(qū)盡量大，在高雷諾數(shù)轉(zhuǎn)捩區(qū)適當(dāng)小。這里采用以下形式：

下面的主要任務(wù)是確定Flength函數(shù)，先選取一個(gè)平板算例，按以下幾個(gè)步驟進(jìn)行迭代：

（1）以常數(shù)Flength進(jìn)行計(jì)算，并與試驗(yàn)結(jié)果比較，記下與試驗(yàn)吻合較好的Flength值；

（2）Flength值不變，代入常數(shù)Reθc進(jìn)行計(jì)算，當(dāng)獲得與（1）相同結(jié)果時(shí)記錄此時(shí)Reθc值；

（3）對比Reθc的表達(dá)式獲得此時(shí)的eθt值；

最后以擬合公式對四個(gè)算例進(jìn)行驗(yàn)證計(jì)算，如有必要，可對Reθc表達(dá)式進(jìn)行適當(dāng)調(diào)整。平板算例中，F(xiàn)length函數(shù)只有在當(dāng)?shù)刂迪虏艑D(zhuǎn)捩進(jìn)行作用，因此，擬合公式幾乎不需要調(diào)整。

表2 平板算例關(guān)聯(lián)函數(shù)值Table 2 Correlations for Reθcand Flengthas function ofeθt

表2 平板算例關(guān)聯(lián)函數(shù)值Table 2 Correlations for Reθcand Flengthas function ofeθt

Case ~Reθt ~Reθc Flength T3B 107．136 102．2 37．01 T3A 205．816 178 31．74 T3A － 714．516 508．6 0．46 S＆K 927．066 685．1 0．40

圖1 壁面摩擦阻力系數(shù)與試驗(yàn)的比較Fig．1 Comparisons of skin friction coefficient distribution between CFD calculations and wind tunnel tests

3 驗(yàn)證算例與討論

3．1 A－aerofoil

A－aerofoil是驗(yàn)證轉(zhuǎn)捩模型的標(biāo)準(zhǔn)算例之一。目前存在兩套試驗(yàn)數(shù)據(jù)，ONERA F1風(fēng)洞迎角13．1°和F2風(fēng)洞迎角13．3°的試驗(yàn)數(shù)據(jù)，這里采取了與文獻(xiàn)［5］相同的比較方式。計(jì)算參數(shù)：M＝0．15，Re＝2．07×106，α＝13．1°，13．3°。湍流度 FSTI＝0．2%，粘性比RT＝10。計(jì)算網(wǎng)格（圖4）481×97（流向×法向），遠(yuǎn)場延伸100倍弦長，網(wǎng)格第一層1．0×10－5C，滿足y＋≤1的要求。

數(shù)值模擬顯示，翼型上表面在0．12C處出現(xiàn)層流分離，且導(dǎo)致分離誘導(dǎo)轉(zhuǎn)捩，向下游發(fā)展為湍流，并在翼型后緣出現(xiàn)湍流分離。轉(zhuǎn)捩模型計(jì)算得到升力為1．5741，阻力 0．0180，試驗(yàn)值分別 為 1．562、0．0208。全湍流1．4470、0．0282，升力偏小，阻力偏大。與文獻(xiàn)［5］中全湍流阻力差別較大。全湍流結(jié)果在后緣與試驗(yàn)結(jié)果吻合更好些。

湍流項(xiàng)和轉(zhuǎn)捩方程離散均采用二階精度。圖5比較了壓力分布，湍流計(jì)算壓力系數(shù)比試驗(yàn)明顯偏低。圖6中比較了湍流及轉(zhuǎn)捩變量用一階精度離散的摩阻分布，采用轉(zhuǎn)捩模型大大提高了摩阻計(jì)算精度，同時(shí)摩阻分布顯示一階精度離散導(dǎo)致轉(zhuǎn)捩位置略微提前。

圖5 壁面壓力系數(shù)曲線比較Fig．5 Comparison of wall pressure coefficient curve with expriment

圖6 摩阻系數(shù)曲線比較Fig．6 Comparison of wall friction coefficient curve with experiment

3．2 NLR7301兩段翼型

NLR7301兩段翼型的風(fēng)洞試驗(yàn)是在70年代末期在NLR3m×2m低速風(fēng)洞中完成的，試驗(yàn)結(jié)果包括了總體氣動(dòng)特性、壓力分布、典型站位的速度型分布等多種數(shù)據(jù)［9］。該翼型的襟翼偏角為20°、來流馬赫數(shù)0．185，試驗(yàn)雷諾數(shù)2．51×106、主翼／襟翼重疊區(qū)均為5．3%c；試驗(yàn)構(gòu)型包含了兩種不同的縫隙寬度，一種為2．6%c，另一種為1．3%c。本文計(jì)算了主翼／襟翼的縫隙寬度2．6%c的構(gòu)型，其中c為弦長。

轉(zhuǎn)捩模型的應(yīng)用就是為了提高阻力計(jì)算精度，文獻(xiàn)［10］指出，高升力翼型的遠(yuǎn)場邊界距離對氣動(dòng)特性、尤其是阻力特性影響顯著，這要求對遠(yuǎn)場邊界進(jìn)行修正，或者擴(kuò)大遠(yuǎn)場邊界。然而對于本文使用的代碼，遠(yuǎn)場邊界取多大合適，這里首先采用環(huán)量修正遠(yuǎn)場邊界條件在不同大小遠(yuǎn)場的網(wǎng)格下進(jìn)行驗(yàn)證計(jì)算。

環(huán)量修正遠(yuǎn)場邊界條件其實(shí)質(zhì)為：由于翼型的存在，遠(yuǎn)場的自由來流受到了擾動(dòng)，不再是均勻流動(dòng)。遠(yuǎn)場的擾動(dòng)可以近似用位于翼型壓心點(diǎn)渦代替，由于阻力一般比升力小2個(gè)量級(jí)，因此阻力代表的點(diǎn)源不予考慮。點(diǎn)渦誘導(dǎo)速度為：

ξ、η為氣流坐標(biāo)系下坐標(biāo)。從上式可以看到，誘導(dǎo)速度與升力成正比，與網(wǎng)格遠(yuǎn)場范圍成反比，因此在升力較大（如高升力翼型），而遠(yuǎn)場不夠遠(yuǎn)的范圍下，誘導(dǎo)速度必須考慮。ξ、η表達(dá)式為：

代碼中實(shí)際應(yīng)用的遠(yuǎn)場條件通過變化得到：

壓力、密度可以由等熵條件求得。

計(jì)算基準(zhǔn)網(wǎng)格如圖7所示，為2004年中國空氣動(dòng)力研究與發(fā)展中心與中國航空計(jì)算技術(shù)研究所聯(lián)合舉辦的“CFD統(tǒng)一算例研討”活動(dòng)發(fā)布的標(biāo)準(zhǔn)網(wǎng)格。采用多塊對接結(jié)構(gòu)網(wǎng)格，網(wǎng)格分7塊，網(wǎng)格單元數(shù)為15萬，主翼布置457個(gè)點(diǎn)，襟翼布置249個(gè)點(diǎn)，尾跡區(qū)89個(gè)，壁面第一層網(wǎng)格距離1．0×10－5c，遠(yuǎn)場邊界延伸12c。為研究網(wǎng)格外邊界的影響，這里在基準(zhǔn)網(wǎng)格基礎(chǔ)上又做了7套網(wǎng)格進(jìn)行對比計(jì)算。其他網(wǎng)格均是在基準(zhǔn)網(wǎng)格基礎(chǔ)上應(yīng)用雙曲型生成方法推進(jìn)得到，網(wǎng)格徑向間距放大因子取為1．1，網(wǎng)格生成過程始終保證了大網(wǎng)格內(nèi)部網(wǎng)格與小網(wǎng)格完全一致，計(jì)算迎角為6°。

計(jì)算中采用Roe格式，應(yīng)用4重多重網(wǎng)格及預(yù)處理，隱式離散用AFADI求解。翼型湍流度取 ，粘性比RT＝5。在網(wǎng)格外邊界較大時(shí)，湍流強(qiáng)度可能衰減厲害，導(dǎo)致計(jì)算中轉(zhuǎn)捩位置偏后，甚至可能出現(xiàn)氣動(dòng)力振蕩。因此本算例應(yīng)用了Spalart［11］建議的環(huán)境源項(xiàng)法來控制湍流度衰減，本文具體做法是使用物面距離來判斷，即y＞10，控制源項(xiàng)作用，反之則使其自由衰減。也可以使用x坐標(biāo)來控制。計(jì)算中應(yīng)用了環(huán)量修正［12］來作對比計(jì)算。

圖7 NLR7301計(jì)算基準(zhǔn)網(wǎng)格及拓展網(wǎng)格Fig．7 Base grid and big grid of NLR7301

計(jì)算結(jié)果見表3、表4，對于無環(huán)量修正情況，阻力隨著遠(yuǎn)場擴(kuò)大單調(diào)減少，升力單調(diào)增加，達(dá)到6000C后收斂。從120C的結(jié)果來看，阻力距離最后收斂值還差了27count。應(yīng)用標(biāo)準(zhǔn)網(wǎng)格計(jì)算的結(jié)果與最后收斂值相比大了接近100%。加入遠(yuǎn)場環(huán)量修正后，網(wǎng)格外邊界達(dá)到120C后氣動(dòng)力就不再變化，且最終收斂值與無環(huán)量修正時(shí)完全一致。表中CDp為壓力貢獻(xiàn)，CDv為摩擦貢獻(xiàn)?？梢钥吹侥ψ柝暙I(xiàn)在12C網(wǎng)格下就可以滿足需要，阻力的主要差別在于壓力的阻力分量。

因此，在高升力翼型設(shè)計(jì)中，計(jì)算網(wǎng)格遠(yuǎn)場邊界至少應(yīng)達(dá)到100C以上，盡量應(yīng)用環(huán)量修正遠(yuǎn)場邊界。

表3 無環(huán)量修正時(shí)計(jì)算的升力、阻力Table 3 Lift and drag without vortex correction

表4 有環(huán)量修正時(shí)計(jì)算的升力、阻力Table 4 Lift and drag with vortex correction

下面以300C的網(wǎng)格計(jì)算，加環(huán)量修正，研究轉(zhuǎn)捩對阻力的影響，計(jì)算方法及條件同上，分別計(jì)算了6°、10．1°、13．1°狀態(tài)。計(jì)算結(jié)果見表5。應(yīng)用轉(zhuǎn)捩模型后，阻力最大誤差4%，全湍流計(jì)算最大誤差達(dá)34．3%。本算例中，轉(zhuǎn)捩計(jì)算提高阻力計(jì)算精度的同時(shí)，也高估了升力，這是以后應(yīng)當(dāng)考慮的問題。

表5 計(jì)算升力、阻力與試驗(yàn)比較Table 5 Comparison of lift and drag with experiment

表6給出了NLR－7301機(jī)翼上翼面、下翼面和襟翼上表面的轉(zhuǎn)捩位置與計(jì)算的比較。本文計(jì)算中轉(zhuǎn)捩位置的判定是以摩阻開始增加為依據(jù)。從表中可以看到，計(jì)算主翼上表面和襟翼上表面轉(zhuǎn)捩位置與試驗(yàn)吻合較好，主翼下表面計(jì)算轉(zhuǎn)捩位置靠前。圖8給出了NLR7301不同迎角下的摩阻分布，從圖上可以看到，計(jì)算摩阻分布轉(zhuǎn)捩位置6°和13．1°的結(jié)果與試驗(yàn)轉(zhuǎn)捩位置很接近，但是在計(jì)算中用摩阻系數(shù)第一次開始增加來判斷不一定合理。以6°為例，摩阻系數(shù)從0．625C位置開始增加，但是增長極其緩慢，在0．68C處存在拐點(diǎn)，從這點(diǎn)開始才急劇增長。同時(shí)，圖8結(jié)果表明，應(yīng)用經(jīng)驗(yàn)關(guān)聯(lián)轉(zhuǎn)捩模型最重要的是關(guān)聯(lián)函數(shù)的標(biāo)定，由于標(biāo)定過程均在平板零壓力梯度流動(dòng)中進(jìn)行，因此在實(shí)際外形流動(dòng)的表現(xiàn)還有待改進(jìn)。

表6 NLR－7301兩段翼型轉(zhuǎn)捩位置（GA P2．6%）Table 6 Transition location of NLR－7301（GA P2．6%）

圖9的升力、阻力收斂歷程比較可以看到，轉(zhuǎn)捩與全湍流計(jì)算在多重網(wǎng)格迭代2000步后氣動(dòng)力就不再變化，轉(zhuǎn)捩對收斂速度影響很小，適合工程應(yīng)用。圖10給出NLR7301迎角6°湍流強(qiáng)度等值線，圖中白線部分為多塊網(wǎng)格中的間隙。由于轉(zhuǎn)捩模型采用了輸運(yùn)方程，因此湍流強(qiáng)度在網(wǎng)格間連續(xù)、光滑變化，相應(yīng)的間歇函數(shù)在網(wǎng)格塊間也光滑變化，這在一定程度上增強(qiáng)了代碼的健壯性。

圖10 NLR7301迎角6°湍流強(qiáng)度等值線Fig．10 NLR7301contours of turbulence intensity

4 小 結(jié)

（1）應(yīng)用本文方法確定了一組關(guān)聯(lián)函數(shù)，在平板測試中摩阻系數(shù)與試驗(yàn)吻合較好，說明了本文標(biāo)定方法的正確性；

（2）A翼型的計(jì)算中，轉(zhuǎn)捩結(jié)果的壁面壓力系數(shù)及摩阻系數(shù)與實(shí)驗(yàn)的吻合程度遠(yuǎn)遠(yuǎn)好于全湍流計(jì)算，說明本文確定的關(guān)聯(lián)函數(shù)在一定程度上能模擬轉(zhuǎn)捩效應(yīng)；

（3）NLR7301計(jì)算研究表明，要提高高升力翼型阻力計(jì)算精度，遠(yuǎn)場邊界必須足夠大（100C以上），盡量應(yīng)用環(huán)量修正遠(yuǎn)場邊界條件，同時(shí)應(yīng)用轉(zhuǎn)捩模型。

（4）轉(zhuǎn)捩模型還很不完善，標(biāo)定點(diǎn)太少，需要在今后應(yīng)用中不斷摸索。

致謝：感謝空氣動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室張玉倫、王光學(xué)、孟德虹同志的熱情幫助。

［1］ LANGTRY R B，MENTER F R．Correlation－based transition modeling for unstructured parallelized computational fluid dynamics codes［J］．AIAA Journal，2009，47（12）：2894－2906．

［2］ MENTER F R，LANGTRY R B，LIKKI S R，et al．A correlation based transition model using local variables：part I－model formulation ASME－GT2004－53452［A］．Proceedings of ASME Turbo Expo 2004［C］，Vienna，Austria，2004．

［3］ LANGTRY R B，MENTER F R，LIKKI S R，et al．A correlation based transition model using local variables：part II－model formulation ASME－GT2004－53454［A］．Proceedings of ASME Turbo Expo 2004［C］，Vienna，Austria，2004．

［4］ LANGTRY R B，MENTER F R．Transition modeling for general CFD Applications in aeronautics［R］．AIAA Paper 2005－0522．

［5］ LANGTRY R B．A correlation－based transition model using local variables for unstructured parallelized CFD codes［D］．［Ph．D．thesis］，Stuttgart，2006．

［6］ MISAKA T，OBAYASHI S．A correlation－based transition models to flows around wings［R］．AIAA Paper 2006－918．

［7］ KARL PETTERSSON，SIMONE CRIPPA．Implementation and verification of a correlation based transition prediction method［R］．AIAA 2008－4401．

［8］ PAUL MALAN，KEERATI SULUKSNA．Calibrating the gamma－Re＿theta transition model for commercial CFD［R］．AIAA 2009－1142．

［9］ VAN DEN BERG B．Boundary layer measurements on a two－dimensional wing with flap［R］．NLR TR 79009U，1979．

［10］CHRISTOPHER L RUMSEY，SUSAN X YING．Prediction of high lift：review of present CFD capability［J］．Progress in Aerospace Sciences，2002，38：145－180．

［11］PHILIPPE R SPALART，CHRISTOPHER L RUMSEY，Effective in flow conditions for turbulence models in aerodynamic calculations［J］．AIAA Journal，2007，45（10）：2544－2553．

［12］THOMAS J L，SALAS M D．Far field boundary conditions for transonic lifting solutions to the Euler equations［R］．AIAA－85－0020．

［13］張玉倫，王光學(xué)．γ－Reθ轉(zhuǎn)捩模型的標(biāo)定研究［A］．工程湍流轉(zhuǎn)捩學(xué)術(shù)研討會(huì)［C］．海口，2010．

［14］孟德虹，張玉倫．γ－Reθ轉(zhuǎn)捩模型在二維低速問題中的應(yīng)用研究［A］．工程湍流轉(zhuǎn)捩學(xué)術(shù)研討會(huì)［C］．?？?，2010．

［15］張曉東，高正紅．關(guān)于補(bǔ)充Langtry的轉(zhuǎn)捩模型經(jīng)驗(yàn)修正式的數(shù)值探討［J］．應(yīng)用數(shù)學(xué)和力學(xué)，2010，31（5）：544－552．