曾衛(wèi)剛，李維東，王發(fā)民

（中國(guó)科學(xué)院力學(xué)研究所，北京 100190）

0 引 言

傳統(tǒng)的再入飛行器大都采用大鈍頭低升阻比構(gòu)型，高溫氣體效應(yīng)會(huì)對(duì)其氣動(dòng)特性產(chǎn)生較大的影響，國(guó)內(nèi)外學(xué)者對(duì)此做了較多的分析研究［1－5］。對(duì)于新型的近空間飛行器（如：HTV1，HTV2），其采用的是細(xì)長(zhǎng)體高升阻比的乘波體布局，高空高馬赫數(shù)飛行條件下高溫氣體效應(yīng)對(duì)其氣動(dòng)特性的影響是十分重要的工程科學(xué)問(wèn)題。但是以往高溫氣體效應(yīng)的研究對(duì)象多集中在再入飛行器，而關(guān)于高溫氣體效應(yīng)對(duì)乘波體氣動(dòng)特性影響的研究并不多見(jiàn)。1995年Inger G R［6－7］就化學(xué)非平衡條件下，壁面催化對(duì)乘波體駐點(diǎn)區(qū)及上下表面中軸線的氣動(dòng)熱特性影響做了較為詳細(xì)的理論分析。2005年Starkey R P［8］以乘波體作為再入飛行器構(gòu)型，對(duì)乘波體在壁面完全催化條件下的壁面熱流及輻射平衡溫度進(jìn)行了初步研究。上述研究探討的是高溫氣體效應(yīng)對(duì)乘波體氣動(dòng)特性的影響，而關(guān)于考慮高溫氣體效應(yīng)的乘波體設(shè)計(jì)，Anderson等人［9］在1992年做了初步研究，其研究結(jié)果表明在馬赫數(shù)不超過(guò)50的情形下，無(wú)粘流場(chǎng)中的平衡氣體效應(yīng)對(duì)乘波體氣動(dòng)優(yōu)化設(shè)計(jì)影響不大。

本文的工作就是：針對(duì)乘波體在高空以高馬赫飛行的特點(diǎn)，給出了一種能考慮高溫平衡氣體效應(yīng)的乘波體設(shè)計(jì)方法；并通過(guò)數(shù)值模擬的方法研究了高溫平衡氣體效應(yīng)對(duì)乘波體氣動(dòng)力特性的影響及高溫平衡氣體效應(yīng)對(duì)乘波體關(guān)鍵熱防護(hù)部位的氣動(dòng)熱特性和輻射平衡溫度的影響。

1 考慮平衡氣體效應(yīng)的乘波體設(shè)計(jì)方法

1.1 平衡氣體錐形繞流的乘波體生成方法

對(duì)于新型的在高空以高馬赫數(shù)飛行的乘波體構(gòu)型飛行器，其飛行過(guò)程中要面臨著高溫真實(shí)氣體效應(yīng)等氣動(dòng)問(wèn)題，本文針對(duì)乘波體飛行器的這一特點(diǎn)給出了考慮平衡氣體效應(yīng)的錐導(dǎo)乘波體設(shè)計(jì)方法。

平衡氣體錐形繞流激波后的流動(dòng)亦屬于軸對(duì)稱(chēng)等熵流動(dòng)，其控制方程可以描述為：

求得平衡氣體錐形繞流的流場(chǎng)后，我們便可通過(guò)流線追蹤得到如圖1所示長(zhǎng)度固定為10m的考慮平衡氣體效應(yīng)的錐導(dǎo)乘波體。

圖1 平衡氣體錐導(dǎo)乘波體示意圖Fig.1 Waverider derived from chemical equilibrium conical flow

1.2 乘波體鈍化方式

這里我們采用加材料的方式用圓弧對(duì)錐導(dǎo)乘波體進(jìn)行鈍化，圖2給出的本文所設(shè)計(jì)的鈍前緣錐導(dǎo)乘波體構(gòu)型示意圖。

圖2 鈍前緣乘波體示意圖Fig.2 Configuration of Waverider with blunt leading edge

2 數(shù)值計(jì)算方法及計(jì)算程序驗(yàn)證

2.1 數(shù)值計(jì)算方法

在笛卡爾坐標(biāo)系下三維無(wú)源非定常N－S方程組的守恒積分形式為：

為了使N－S方程組封閉，我們需要補(bǔ)充相應(yīng)的氣態(tài)方程、熱力學(xué)參數(shù)計(jì)算公式及輸運(yùn)參數(shù)的計(jì)算公式。

對(duì)于完全氣體，只需補(bǔ)充相應(yīng)的顯示表達(dá)式即可。對(duì)于平衡氣體，其熱力學(xué)參數(shù)及輸運(yùn)參數(shù)之間沒(méi)有明顯的顯示表達(dá)式，對(duì)于平衡氣體熱力學(xué)特性參數(shù)，本文采用Tannehill等人［10－12］給出的曲線擬合方法進(jìn)行計(jì)算；對(duì)于平衡氣體輸運(yùn)參數(shù)，本文采用Srinivasan等人［13］給出的多項(xiàng)式擬合方法進(jìn)行計(jì)算。

計(jì)算過(guò)程中空間離散采用AUSM格式，并配合Minimod限制器來(lái)抑制數(shù)值振蕩，時(shí)間項(xiàng)用LUSGS隱式迭代格式進(jìn)行推進(jìn)。另外計(jì)算過(guò)程中遠(yuǎn)場(chǎng)采用自由來(lái)流條件；對(duì)稱(chēng)面采用對(duì)稱(chēng)邊界條件；出口采用數(shù)值外插得到；壁面采用無(wú)滑移條件。

2.2 計(jì)算程序驗(yàn)證

為了檢驗(yàn)計(jì)算程序的可靠性，本文對(duì)文獻(xiàn)［14］給出算例1進(jìn)行程序驗(yàn)證。該算例是一高超聲速球頭繞流，球頭半徑R＝0.25m，來(lái)流馬赫數(shù) M∞＝11.26，來(lái)流密度ρ＝3.99×10－5kg／m3，來(lái)流靜溫T∞＝182.33K，壁溫Tw＝1000K。

圖3給出的是表面壓力系數(shù)分布，圖4給出的是壁面熱流系數(shù)分布。從圖中我們可以看出，與文獻(xiàn)［14］的計(jì)算結(jié)果相比，本文計(jì)算程序所得到的壁面壓力系數(shù)和壁面熱流系數(shù)都較為接近，并且二者在捕捉平衡氣體與完全氣體的的差別時(shí)具有較好的一致性。

圖3 球面壓力分布對(duì)比Fig.3 Comparison of surface pressure distributions over sphere

圖4 球面熱流分布對(duì)比Fig.4 Comparison of surface heating rates distributions over sphere

3 平衡氣體對(duì)乘波體氣動(dòng)力特性的影響

乘波體作為具有典型細(xì)長(zhǎng)體特征的高超聲速飛行器，其鈍頭體繞流與尖頭體繞流具有很大的不同，并且這將在一定程度上改變平衡氣體效應(yīng)對(duì)乘波體氣動(dòng)性能的影響。

圖5～圖6分別給出的是60km高空，來(lái)流馬赫數(shù)為20時(shí)，用兩種氣體模型計(jì)算所得的鈍化半徑為0.025m的乘波體的熵層云圖，從圖中我們可以看出，平衡氣體效應(yīng)大大地降低了由鈍化效應(yīng)引起的熵增，這主要是因?yàn)槠胶鈿怏w使得激波更貼近物面，從而導(dǎo)致熵層影響減弱，這與Singh［15－16］在研究鈍體繞流時(shí)的結(jié)論是一致的。另外無(wú)論是完全氣體還是平衡氣體，乘波體上表面的熵層都很厚，基本上都處在熵層影響范圍內(nèi)，而下表面的熵層的影響范圍只局限在乘波體前沿線附近，這主要是因?yàn)槌瞬w上下表面對(duì)應(yīng)的流線與鈍化圓弧可構(gòu)成二維類(lèi)鈍錐，并且這種類(lèi)鈍錐上表面的等效錐角較小，下表面的等效錐角較大，并且由Singh［15－16］的研究結(jié)果我們可知在其他條件一致的條件下，錐角越小則熵層影響范圍就會(huì)越大。

圖5 完全氣體條件下流場(chǎng)熵層分布云圖Fig.5 Entropy contour of flow field under the perfect gas condition

圖6 平衡氣體條件下流場(chǎng)熵層分布云圖Fig.6 Entropy contour of flow field under the equilibrium gas condition

圖7給出的是該乘波體中軸線上的壓力分布，從圖在我們可以看出，對(duì)于乘波體下表面，無(wú)論是平衡氣體還是完全氣體，其壓力分布都會(huì)產(chǎn)生一個(gè)先下降后上升的過(guò)程，這與張涵信［17］在分析鈍體繞流的壓力分布特點(diǎn)時(shí)的結(jié)果也是相吻合的。另外對(duì)于乘波體上表面，平衡氣體的壓力在前緣線后的區(qū)域都相對(duì)偏低，而對(duì)于下表面平衡氣體的壓力在前緣線稍后的地方偏低，但過(guò)后又會(huì)慢慢超過(guò)完全氣體的壓力。按照文獻(xiàn)［18］對(duì)高超聲速鈍體繞流壁面壓力分布的分析我們可知高超聲速鈍體繞流的物面壓力只與來(lái)流條件和熵層外邊界的當(dāng)?shù)匦甭视嘘P(guān)，并且來(lái)流條件一定時(shí)當(dāng)?shù)匦甭试酱髣t物面壓力越大。由Singh［15－16］的研究結(jié)果我們可知，由于平衡氣體的激波更靠近物面，這會(huì)減弱熵層影響范圍，降低熵層外邊界的當(dāng)?shù)匦甭?。另外由前面分析我們可知乘波體上表面基本上處于熵層影響范圍內(nèi)，故對(duì)于乘波體上表面，平衡氣體的壓力在前緣線后都相對(duì)偏低；而乘波體下表面熵層只影響靠近前緣線附近的區(qū)域，故對(duì)于下表面，平衡氣體的壓力在前緣線稍后的地方會(huì)偏低，對(duì)于乘波體下表面遠(yuǎn)離前緣線的地方，由于熵層影響較弱，加之化學(xué)反應(yīng)會(huì)使得壓力相對(duì)升高，故乘波體下表面離前緣線較遠(yuǎn)的地方平衡氣體的壓力會(huì)比完全氣體的壓力稍大一點(diǎn)。

圖7 乘波體上下表面中軸線壓力分布Fig.7 Pressure distributions along the central line of waverider surface

圖8 升阻比隨鈍化半徑的變化關(guān)系Fig.8 Variation of lift to drag ratio with radio of bluntness

圖8～圖9分別給出的是升阻比和俯仰力矩系數(shù)隨設(shè)計(jì)半徑的變化關(guān)系曲線。圖8說(shuō)明鈍化會(huì)使得乘波體下表面高壓氣體泄露到上表面，從而導(dǎo)致升阻比迅速下降。圖9說(shuō)明平衡氣體會(huì)對(duì)乘波體產(chǎn)生一個(gè)附加的低頭力矩，這主要是由平衡氣體使得乘體上表面壓力有所降低的同時(shí)，又會(huì)使得乘波體下表面遠(yuǎn)離前緣線區(qū)域的壓力升高所致?？傮w看來(lái)乘波體前緣激波較弱，不足以強(qiáng)烈壓縮空氣引起激烈的化學(xué)反應(yīng)，因而平衡氣體效應(yīng)對(duì)乘波體的氣動(dòng)力影響不是太大。

圖9 俯仰力矩系數(shù)隨鈍化半徑的變化關(guān)系Fig.9 Variation of pitching moment coefficients with radius of bluntness

4 平衡氣體效應(yīng)對(duì)氣動(dòng)熱特性影響

乘波體飛行時(shí)以長(zhǎng)時(shí)間中低熱流進(jìn)行氣動(dòng)加熱，其壁面防熱應(yīng)為疏導(dǎo)式、主動(dòng)防熱。輻射平衡條件下壁面熱流應(yīng)滿足：

其中在式（8）中（q（Tweq））Convect為對(duì)流引起的氣動(dòng)加熱；（q（Tweq））Radiation為壁面輻射所帶走的熱量，其具體計(jì)算式為：

式（9）中ε為輻射系數(shù)，本為取ε＝0.87。

本文通過(guò)數(shù)值模擬的方法與工程算法相結(jié)合的方式來(lái)研究平衡氣體對(duì)乘波體在輻射平衡條件下的駐點(diǎn)熱流和輻射平衡溫度的影響。

對(duì)于完全氣體件下乘波體駐點(diǎn)熱流的工程算法，本文采用Vanmol D O［19］給出的駐點(diǎn)熱流公式：

對(duì)于平衡氣體駐點(diǎn)熱流的工程算法，本文采用Andrushchenko V A［20］等人給出平衡氣體駐點(diǎn)熱流公式的公式：

（10）與（11）兩式中V∞為來(lái)流速度，hw為壁面焓值，h0為滯止焓，ρ∞為來(lái)流密度，R為鈍化半徑，Г為鈍化度，其定義為：

其中Rn為乘波體前緣線駐點(diǎn)處的曲率半徑。

圖10～圖11分別給出的是駐點(diǎn)熱流及輻射平衡溫度隨設(shè)計(jì)鈍化半徑的變化關(guān)系曲線，從圖中我們可以看出，工程算法與數(shù)值方法的計(jì)算結(jié)果都很接近；平衡氣體效應(yīng)會(huì)使得駐點(diǎn)熱流和輻射平衡溫度都有所增大，并且隨著設(shè)計(jì)鈍化半徑的增大，兩種氣體模型計(jì)算所得到的駐點(diǎn)熱流和駐點(diǎn)輻射平衡溫度都會(huì)變小。

圖10 駐點(diǎn)熱流隨鈍化半徑的的變化Fig.10 Variation of stagnation point heat flux with radius of bluntness

圖11 駐點(diǎn)輻射平衡溫度隨鈍化半徑的的變化Fig.11 Variation of stagnation point radiative equili－brium temperature with radius of bluntness

圖12～圖13分別給出的是60km高空，來(lái)流馬赫數(shù)為20時(shí)鈍化半徑為0.05m的乘波體前緣區(qū)熱流分布云圖及輻射平衡溫度云圖。從圖中我們同樣可以明顯看出對(duì)于前緣駐點(diǎn)區(qū)，平衡氣體較完全氣體的熱流及輻射平衡溫度都要有所增大。

圖14～圖15分別給出的相同條件下乘波體上下表面中軸線熱流分布及輻射平衡溫度分布，我們可以看出平衡氣體使得乘波體下表面的熱流及輻射平衡溫度有所升高，但對(duì)乘波體上表面的熱流及輻射平衡溫度溫度幾乎沒(méi)有影響，這主要是因?yàn)槌瞬w下表面化學(xué)反應(yīng)較上表面的化學(xué)反應(yīng)更為激烈。

圖12 乘波體駐點(diǎn)區(qū)熱流分布云圖Fig.12 Heat flux contour near the nose region of waverider

圖13 乘波體駐點(diǎn)區(qū)輻射平衡溫度分布云圖Fig.13 Radiative equilibrium temperature contour near the nose region of waverider

圖14 上下表面中軸線上的熱流分布Fig.14 Distribution of heat flux along the central line

圖15 上下表面中軸線上的輻射平衡溫度分布Fig.15 Distribution of radiative equilibrium temperature along the central line

5 非設(shè)計(jì)攻角下平衡氣體效應(yīng)影響

對(duì)于高超聲速飛行器，其非設(shè)計(jì)工況下的氣動(dòng)性能是相當(dāng)重要的，這里我們就平衡氣體效應(yīng)對(duì)鈍前緣乘波體在較為適宜的飛行攻角下的氣動(dòng)性能影響進(jìn)行分析。出于氣動(dòng)防熱要求，乘波體在使用時(shí)必須進(jìn)行鈍化，為了使得乘波體駐點(diǎn)區(qū)防熱都滿足要求，并兼顧乘波體需有較高的升阻比，我們選擇乘波體的設(shè)計(jì)鈍化半徑為0.05m，并在此基礎(chǔ)上考慮非設(shè)計(jì)攻角下平衡氣體效應(yīng)對(duì)其氣動(dòng)特性的影響。

圖16～圖17分別給出的是60m高空，來(lái)流馬赫數(shù)為20時(shí)，乘波體升阻比和俯仰力矩系數(shù)隨飛行攻角的變化關(guān)系曲線。從圖中我們可以看出，在乘波體較為適宜的飛行攻角范圍內(nèi)，平衡氣體效應(yīng)使得乘波體的升阻比略微降低；并且平衡氣體效應(yīng)對(duì)乘波體俯仰力矩的影響相對(duì)攻角的影響幾乎可以忽略。這與葉友達(dá)［4－5］在研究平衡氣體對(duì)再入飛行器氣動(dòng)力特性的影響時(shí)所得到規(guī)律有很大的區(qū)別，其主要原因是因?yàn)槌瞬w前緣激波較弱，不能像大鈍體再入飛行軌道器那樣強(qiáng)烈壓縮空氣引起激烈的化學(xué)反應(yīng)來(lái)改變飛行器氣動(dòng)力特性。

圖16 升阻比隨攻角的變化關(guān)系Fig.16 Variation of lift to drag ratio with angle of attack

圖18～圖19分別給出的是相同條件下兩種氣體模型計(jì)算所得到的乘波體駐點(diǎn)熱流及輻射平衡溫度隨攻角的變化關(guān)系曲線，從圖中我們可以看出隨著攻角的增大，乘波體駐點(diǎn)熱流及輻射平衡溫度的變化并不明顯。

圖20～圖21分別給出的是用平衡氣體模型計(jì)算所得到的不同攻角下乘波體下表面中軸線上的熱流及輻射平衡溫度相對(duì)于完全氣體模型的增值分布。從圖中我們可以看出隨著攻角的增大，兩種氣體模型計(jì)算所得到的乘波體下表面中軸線上更多區(qū)域的熱流差值和輻射平衡溫度差值會(huì)趨近于零，這說(shuō)明攻角的增大會(huì)減弱平衡氣體效應(yīng)對(duì)乘波體下表面氣動(dòng)熱及輻射平衡溫度的影響，會(huì)使得平衡氣體效應(yīng)對(duì)乘波體下表面的熱流及輻射平衡溫度影響區(qū)域減小。

圖17 俯仰力矩系數(shù)隨攻角的變化關(guān)系Fig.17 Variation of pitching moment coefficients with angle of attack

圖18 駐點(diǎn)熱流隨攻角的變化Fig.18 Variation of stagnation point heat flux with angle of attack

圖19 駐點(diǎn)輻射平衡溫度隨攻角的變化Fig.19 Variation of stagnation point radiative equilibrium temperature with angle of attack

圖20 下表面中軸線熱流增值分布Fig.20 Heat flux increment distributions along the central line of lower surface

圖21 下表面中軸線輻射平衡溫度增值分布Fig.21 Radiative equilibrium temperature increment distributions along the central line of lower surface

6 結(jié) 論

本文設(shè)計(jì)了考慮平衡氣體效應(yīng)的乘波體，并對(duì)平衡氣體效應(yīng)的影響進(jìn)行了計(jì)算和分析，得出了如下結(jié)論：

（1）平衡氣體效應(yīng)對(duì)鈍前緣乘波體氣動(dòng)力特性的影響是通過(guò)與熵層效應(yīng)相互作用來(lái)完成的；鈍化半徑的增大會(huì)使得乘波體升阻比急劇降低，并且乘波體前緣激波較弱，平衡氣體效應(yīng)對(duì)乘波體氣動(dòng)力特性的影響并不十分明顯。

（2）平衡氣體效應(yīng)對(duì)乘波體的局部氣動(dòng)熱特性和輻射平衡溫度有較大的影響；并且鈍化半徑的增大有助于降低乘波體駐點(diǎn)熱流及輻射平衡溫度，但是當(dāng)鈍化半徑增大到一定程度后，這種由鈍化所帶來(lái)的熱防護(hù)作用就不再明顯。

（3）在乘波體較為適宜的飛行攻角范圍內(nèi)，攻角的增大不會(huì)顯著改變平衡氣體效應(yīng)對(duì)乘波體氣動(dòng)力特性的影響，也不會(huì)顯著改變平衡氣體效應(yīng)對(duì)乘波體前緣駐點(diǎn)區(qū)的熱流與輻射平衡溫度的影響，但卻會(huì)在一定程度上減小平衡氣體效應(yīng)對(duì)乘波體下表面熱流及輻射平衡溫度的影響區(qū)域。

由于篇幅所限，本文只是初步地研究了高溫平衡氣體效應(yīng)對(duì)乘波體布局設(shè)計(jì)及氣動(dòng)特性的影響，而實(shí)際上采用求解高溫非平衡氣體N－S方程的方法研究高溫氣體效應(yīng)對(duì)高超聲速飛行器氣動(dòng)特性的影響更為準(zhǔn)確，對(duì)此作者在今后的工作中會(huì)進(jìn)行更深一步的研究。

［1］MAUS J R，GRIFFITH B J，SZEMA K Y，et al.Hypersonic mach number and real gas effects on space shuttle orbiter［R］.AIAA 83－0343.

［2］MAUS J R，GRIFFITH B J，SZEMA K Y.Hypersonic mach number and real gas effects on space shuttle orbiter aerodynamics［J］.Journal of Spacecraft，1984，21（2）：136－141.

［3］MAUS J R，GRIFFITH B J，SZEMA K Y，et al.Understanding space shuttle flight data by use of wind tunnel and CFD results［R］.AIAA 83－2745.

［4］葉友達(dá).高空高速飛行器氣動(dòng)特性研究［J］.力學(xué)進(jìn)展，2009，39（4）：387－397.（YE You－da.Study on aerodynamic characteristics of high velocity vehicle flying at high altitude［J］.Advances In Mechanics，2009，39（4）：387－397.）

［5］葉友達(dá).近空間高速飛行器氣動(dòng)特性研究與布局設(shè)計(jì)優(yōu)化 ［J］.力學(xué)進(jìn)展，2009，39（6）：683－694.（YE You－da.Study on aerodynamic characteristics and design optimization for high speed near space vehicles［J］.Advances In Mechanics，2009，39（6）：683－694.）

［6］INGER G R.Nonequilibrium boundary－layer effects on the aerodynamic heating of hypersonic waverider vehicles［J］.Journal of Thermophysics and Heat Transfer，1995，9（4）：595－604.

［7］INGER G R.Nonequilibrium boundary layer effects on the aerodynamic heating of hypersonic vehicles［J］.Acta Astronautica，1995，36（4）：205－216.

［8］STARKEY R P.Design of waverider based re－entry vehicles［R］.AIAA 2005－3390.

［9］ANDERSON J D，CHANG J，MCLAUGHLIN T A.Hypersonic waveriders：effects of chemically reacting flow and viscous interaction［R］.AIAA－92－0302.

［10］SRINIVASAN S，TANNEHILL J C，WEILMUENSTER K J.Simplified curve fits for the thermodynamic properties of equilibrium air［R］.NASA CR－181245，1986.

［11］SRINIVASAN S，TANNEHILL J C，WEILMUENSTER K J.Simplified curve fits for the thermodynamic properties of equilibrium air［R］.NASA RP－1181，1987.

［12］TANNEHILL J C，MUGGE P H.Improved curve fits for the thermodynamic properties of equilibrium air for numerical computation using time－dependent or shockcaturing methods［R］.NASA CR－2470，1974.

［13］SRINIVASAN S，TANNEHILL J C.Simplified curve fits for the transport properties of equilibrium air［R］.New York：CR－178411，1987.

［14］KAZEM H，RAMIN K M，VAHID E.Dual－code solution procedure for efficient computing equilibrium hypersonic axisymmetric laminar flow［J］.Aerospace Science and Technology，2008，12（2）：135－149.

［15］SINGH D J，KUMAR A，TIWARI S N.Effect of nose bluntness on flow field over slender bodies in hypersonic flows［R］.AIAA 89－0270.

［16］SINGH D J，KUMAR A，TIWARI S N.Effect of nose bluntness on flowfield over slender bodies in hypersonic flows［J］.Journal of Thermophysics，1990，5（2）：166－171.

［17］張涵信.高超聲速運(yùn)動(dòng)中的熵層問(wèn)題［J］.航空學(xué)報(bào)，1965，1（2）：20－43.

［18］曾廣存.熵層分析與鈍體繞流［J］.航空學(xué)報(bào)，1965，1（2）：1－19.

［19］VANMOL D O，ANDERSON J D.Heat transfer characteristics of hypersonic waveriders with an emphasis on leading edge effect［R］.AIAA 92－2920.

［20］ANDRUSHCHENKO V A，SYZRANOVA N G，SHEVELEV Y D.An estimate of the heat fluxes to the surface of blunt bodies moving at hypersonic velocity in the atmosphere［J］.Journal of Applied Mathematics and Mechanics，2007，71（5）：747－754.