• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    基于多場耦合的飛行器熱環(huán)境數(shù)值分析方法研究

    2014-05-04 07:34:07張勝濤陳方劉洪
    空氣動力學(xué)學(xué)報 2014年6期
    關(guān)鍵詞:熱流壁面步長

    張勝濤,陳方,劉洪

    (上海交通大學(xué) 航空航天學(xué)院,上海 200240)

    基于多場耦合的飛行器熱環(huán)境數(shù)值分析方法研究

    張勝濤,陳方,劉洪

    (上海交通大學(xué) 航空航天學(xué)院,上海 200240)

    新一代高超聲速飛行器的發(fā)展給防熱設(shè)計問題帶來嚴(yán)峻挑戰(zhàn)。根據(jù)飛行器熱環(huán)境多場耦合特性,提出了一種基于多場耦合的熱環(huán)境數(shù)值分析策略,并在此基礎(chǔ)上發(fā)展了基于Navier-Stokes方程的流場CED分析程序,通過有效的界面數(shù)據(jù)傳遞算法,實現(xiàn)了與結(jié)構(gòu)有限元熱分析軟件的耦合,形成了基于流場與結(jié)構(gòu)耦合傳熱的飛行器熱環(huán)境多場耦合數(shù)值分析方法。以典型圓管前緣為計算模型進(jìn)行了程序驗證,并對穩(wěn)態(tài)和非穩(wěn)態(tài)飛行環(huán)境下的流場與結(jié)構(gòu)耦合傳熱特征和規(guī)律進(jìn)行了數(shù)值分析研究。結(jié)果分析表明,該方法能夠有效地刻畫流場與結(jié)構(gòu)之間的耦合傳熱特征和規(guī)律,預(yù)測和分析飛行器熱環(huán)境的空間和時間分布特性,從而可為防熱設(shè)計的選材和優(yōu)化提供可靠的參考依據(jù)和分析手段。

    多場耦合;長時間氣動加熱;熱環(huán)境;防熱設(shè)計;高超聲速飛行器

    0 引 言

    隨著新一代高超聲速飛行器的發(fā)展,防熱設(shè)計遭遇了更為嚴(yán)峻的挑戰(zhàn)。一方面是飛行器要面臨高焓、中低熱流和長時間氣動加熱的惡劣熱環(huán)境;另一方面是要在滿足極為苛刻的質(zhì)量要求下實現(xiàn)非燒蝕防熱。為了能夠更合理的選材和優(yōu)化防熱,以便進(jìn)一步減小設(shè)計冗余,這就對飛行器熱環(huán)境的準(zhǔn)確預(yù)測提出了更高的要求。

    新一代高超聲速飛行器的熱環(huán)境特性與其自身氣動外形、所處飛行環(huán)境及其內(nèi)部結(jié)構(gòu)形式和材料物理屬性等密切相關(guān),外部流場的氣動加熱與內(nèi)部結(jié)構(gòu)熱響應(yīng)之間存在強(qiáng)烈的耦合關(guān)系。因此,必須研究和建立基于流場、熱、結(jié)構(gòu)多場耦合的飛行器熱環(huán)境數(shù)值分析方法和手段。

    關(guān)于流場、熱、結(jié)構(gòu)多場耦合基礎(chǔ)問題的研究,國外開展得較早。Thornton和Dechaumphai等人[1-3]于20世紀(jì)80年末首次提出流場、熱、結(jié)構(gòu)多場耦合問題,并對多場耦合特征和關(guān)鍵問題進(jìn)行了初步計算分析。Chen等人[4-6]通過松耦合的方式對流場和結(jié)構(gòu)進(jìn)行迭代求解,率先將多場耦合計算方法應(yīng)用于航天領(lǐng)域。20世紀(jì)初,國內(nèi)學(xué)者黃唐、毛國良等人[7]率先研究了二維流場、熱、結(jié)構(gòu)一體化數(shù)值模擬方法。此后國內(nèi)一些學(xué)者陸續(xù)開展了流場、熱、結(jié)構(gòu)多場耦合計算研究和分析[8-12]。但總體而言,國內(nèi)還主要側(cè)重于穩(wěn)態(tài)飛行環(huán)境下的流場、熱、結(jié)構(gòu)多場耦合基礎(chǔ)問題的理論方法研究。對于飛行器實際熱環(huán)境的準(zhǔn)確預(yù)測,還需要進(jìn)行非穩(wěn)態(tài)飛行環(huán)境下的流場、熱、結(jié)構(gòu)多場耦合分析,而這方面研究相對較少,對其多場耦合規(guī)律和特征的認(rèn)識還不夠深入,相關(guān)分析方法和研究手段還有待進(jìn)一步發(fā)展和完善。

    通過對飛行器熱環(huán)境多場耦合特性的分析,提出了一種基于多場耦合的熱環(huán)境數(shù)值分析策略,并在此基礎(chǔ)上建立了相應(yīng)的基于流場與結(jié)構(gòu)耦合傳熱的飛行器熱環(huán)境多場耦合數(shù)值分析方法。以典型圓管前緣為計算模型進(jìn)行了程序驗證,并對穩(wěn)態(tài)和非穩(wěn)態(tài)飛行環(huán)境下的流場與結(jié)構(gòu)耦合傳熱特征和規(guī)律進(jìn)行了數(shù)值分析研究。

    1 問題分析及耦合分析策略

    高超聲速飛行器在實際飛行過程中所經(jīng)歷的是一個沿特定飛行軌跡變化的持續(xù)、非瞬態(tài)的受熱過程,其物理實質(zhì)是一個在長時間持續(xù)非穩(wěn)態(tài)飛行環(huán)境下的流場、熱、結(jié)構(gòu)多場耦合問題。

    圖1 基于多場耦合的熱環(huán)境數(shù)值分析策略Fig.1 Multi-field coupling numerical analysis strategy for predicting aerothermal environment

    下面進(jìn)一步分析穩(wěn)態(tài)環(huán)境下的流場、熱、結(jié)構(gòu)多場耦合問題。對于新一代高超聲速飛行器防熱設(shè)計,由于要保持良好的氣動外形,不允許結(jié)構(gòu)有較大變形。本文暫不考慮結(jié)構(gòu)變形的影響,假定結(jié)構(gòu)為剛性體,這樣可簡化為流場與結(jié)構(gòu)耦合傳熱問題。

    對于流場與結(jié)構(gòu)耦合傳熱問題,通??刹捎脙深惙绞角蠼猓?)整體求解方式;2)分區(qū)迭代方式。其中,分區(qū)迭代方式是指各場分別采用各自的求解器在時域內(nèi)依次交替時間步進(jìn)行推進(jìn)求解,耦合作用通過界面?zhèn)鬟f耦合變量實現(xiàn)。這種方式可以充分利用各場特點(diǎn)和成熟算法,易于實現(xiàn)和模塊化,因而得到廣泛應(yīng)用。因此,本文采用分區(qū)迭代方式。

    對于高速可壓縮流動,流場被擾動后重新達(dá)到穩(wěn)定的特征時間一般為毫秒量級,而結(jié)構(gòu)傳熱的特征時間一般為秒量級。若采用分區(qū)迭代方式求解流場與結(jié)構(gòu)傳熱問題時,以流場計算時間步長作為耦合時間步長,將導(dǎo)致極大的計算量。研究表明,即使邊界處熱流很大,壁面溫度在較短時間內(nèi)變化也不是很大,其對流場的影響也僅限于附面層極小范圍內(nèi),故與結(jié)構(gòu)傳熱計算時間步長相比,可認(rèn)為流場是瞬時穩(wěn)定的[9]。基于上述實際,采用如圖2所示的流場與結(jié)構(gòu)耦合傳熱分區(qū)迭代求解策略。

    圖2 流場與結(jié)構(gòu)耦合傳熱分區(qū)迭代求解策略Fig.2 Partitioned iteration strategy of the flow-structural coupling heat transfer

    具體步驟為:1)根據(jù)初始時刻t0的結(jié)構(gòu)熱分析結(jié)果,通過耦合界面向流場提供壁面溫度邊界條件Tw;2)根據(jù)壁面溫度邊界條件Tw,計算定常流場得到壁面氣動熱流密度qw;3)通過耦合界面將壁面氣動熱流密度qw傳給結(jié)構(gòu),為熱分析提供熱流密度邊界條件;4)根據(jù)壁面熱流密度邊界條件qw,在t0至t0+Δtth時間段內(nèi)進(jìn)行瞬態(tài)傳熱分析,得到t0+Δtth時刻的結(jié)構(gòu)溫度場分布,至此完成一個耦合計算時間步長Δtth內(nèi)的迭代計算。

    2 控制方程及數(shù)值方法

    2.1 流場

    考慮基于量熱完全氣體模型的三維粘性可壓縮流動Navier-Stokes方程,其在笛卡爾坐標(biāo)系下的守恒積分形式為:

    式中:Q為流場守恒向量;Fc和Fν分別為對流矢通量和粘性矢通量。

    本文基于有限體積法對上述控制方程進(jìn)行離散求解。其中,對流通量采用M-AUSMPW+混合迎風(fēng)格式[13]進(jìn)行離散,并用多維限制器MLP[14]提高至三階精度;粘性通量采用二階中心差分格式進(jìn)行離散;時間推進(jìn)采用LU-SGS隱式方法[15]。為了加速收斂過程,采用當(dāng)?shù)貢r間步長和隱式殘值平均等加速收斂措施。

    2.2 結(jié)構(gòu)溫度場

    固體結(jié)構(gòu)內(nèi)部溫度場的控制方程為導(dǎo)熱微分方程。在笛卡爾直角坐標(biāo)系下,三維非穩(wěn)態(tài)導(dǎo)熱微分方程的表達(dá)式為:

    對于上述結(jié)構(gòu)溫度場控制方程,本文直接采用商用有限元熱分析軟件ANSYS作為求解器進(jìn)行離散求解。

    2.3 界面耦合關(guān)系及數(shù)據(jù)傳遞

    流場與結(jié)構(gòu)耦合傳熱的實質(zhì)是流場和結(jié)構(gòu)溫度場通過飛行器壁面發(fā)生耦合作用。如圖3所示,設(shè)Ωf和Ωs分別為流體域和固體域,Гfs=Ωf∩Ωs為耦合界面,則在耦合界面Гfs處滿足如下關(guān)系:

    對于上述耦合關(guān)系,采用Dirichlet-Neumann方法處理[16],即固體域向流體域提供壁面溫度條件,而流體域向固體域提供熱流密度邊界條件。

    圖3 流體域與固體域耦合關(guān)系示意圖Fig.3 Schematic illustration of thecoupling relationship between fluid and soliddomains

    由于流體域和固體域分別采用了不同的離散方法和數(shù)值格式,在耦合界面上存在兩套非匹配網(wǎng)格之間的數(shù)據(jù)傳遞,這在數(shù)學(xué)上是一個雙向插值問題,本文主要采用Shepard方法[17]來實現(xiàn)。該方法首先由氣象學(xué)家與地質(zhì)工作者提出,其基本思想是將待插值點(diǎn)的函數(shù)值定義為已知數(shù)據(jù)點(diǎn)的函數(shù)值的距離倒數(shù)加權(quán)之和:

    式中:M為待插值點(diǎn)的函數(shù)值;Mi為已知數(shù)據(jù)點(diǎn)的函數(shù)值;di為待插值點(diǎn)與已知數(shù)據(jù)點(diǎn)之間的歐氏距離;n為已知數(shù)據(jù)點(diǎn)的個數(shù);p為距離的冪次,其值越高插值結(jié)果越光滑。

    3 計算結(jié)果及分析

    3.1 計算模型及初場驗證

    高超聲速飛行器頭部通常采用鈍體前緣設(shè)計,故這里選取二維圓管前緣[10,18]作為計算模型。該圓管前緣的內(nèi)外半徑尺寸以及所采用的計算網(wǎng)格如圖4所示。其中,固體域采用有限元網(wǎng)格劃分,最小單元邊長為0.5mm;流體域采用結(jié)構(gòu)化網(wǎng)格劃分,網(wǎng)格大小為201×111,且在激波和近壁面附近進(jìn)行了加密處理。圓管內(nèi)部結(jié)構(gòu)材料為不銹鋼,初始壁面溫度及外部流場的來流參數(shù)如表1所示。

    表1 初始計算條件Table 1 Initial computational condition

    圖4 計算模型及計算網(wǎng)格Fig.4 Computational model and computational mesh

    圖5所示為初始時刻t=0s時外部流場溫度云圖以及沿駐點(diǎn)線溫度分布。可以看出,計算得到的溫度場分布宏觀上符合物理實際,其中最大溫度值與理論總溫一致,流場激波位置與理論位置相符,沿駐點(diǎn)線溫度分布符合正激波前后物理關(guān)系。

    為了考察壁面熱流計算的準(zhǔn)確性,以文獻(xiàn)[19]提出的準(zhǔn)則為依據(jù)確定壁面法向網(wǎng)格尺度。圖6給出了初始時刻t=0s時歸一化的壁面熱流密度分布與試驗結(jié)果的對比,可以看出計算結(jié)果與試驗結(jié)果符合得較好。

    圖5 初始時刻t=0s時溫度場分布Fig.5 Initial temperature fielddistribution at t=0s

    圖6 計算與試驗的壁面熱流分布對比Fig.6 Comparison of wall heat-fluxdistribution between the computation and the experiment

    3.2 穩(wěn)態(tài)飛行環(huán)境下熱環(huán)境特性分析

    首先考察了耦合計算時間步長對壁面溫度的影響情況,如圖7所示,總的持續(xù)時間為5s??梢钥闯觯?dāng)耦合計算時間步長為5s時,即熱載荷一次加載,不存在耦合作用,計算得到的壁面溫度最高;隨著耦合計算時間步長的減小,即耦合次數(shù)的增多,計算的壁面溫度越低。這表明,充分考慮了流場與結(jié)構(gòu)耦合傳熱效應(yīng)后,對結(jié)構(gòu)溫度的預(yù)測更接近實際物理過程,可以進(jìn)一步減小由于結(jié)構(gòu)溫度高估而引起的設(shè)計冗余,且這種效果會在氣動加熱越劇烈的地方越明顯,利于更合理地選材和進(jìn)行優(yōu)化設(shè)計。

    圖7 不同耦合時間步長對壁面溫度的影響Fig.7 Effects ofdifferent coupling computational time step on the predicted wall temperature

    從圖7還可以看出,隨著耦合計算時間步長的減小,其對壁面溫度的影響也逐漸減少,當(dāng)耦合計算時間步長減小至一定程度后,壁面溫度的變化已經(jīng)不再明顯。這表明,隨著耦合計算時間步長的減少,計算精度逐漸增加,流場與結(jié)構(gòu)之間的耦合傳熱逐漸接近實際的物理過程。

    圖8給出了采用耦合計算時間步長為Δt=0.5s時的壁面溫度分布隨時間的變化情況。隨著時間的持續(xù),壁面溫度從初始時刻的分布逐漸上升,其變化過程符合物理實際,最終t=5s時刻的壁面溫度分布量值上略低于試驗結(jié)果,但分布規(guī)律趨勢與試驗結(jié)果相一致。

    圖8 壁面溫度分布隨時間的變化Fig.8 Variation of wall temperaturedistribution with the time

    為了進(jìn)一步考察長時間持續(xù)飛行對熱環(huán)境特性的影響,這里將總的持續(xù)飛行時間延長至50s,采用耦合計算時間步長為0.5s,進(jìn)一步進(jìn)行了分析。

    圖9和圖10分別給出了壁面駐點(diǎn)處溫度、氣動加熱熱流及它們的時間導(dǎo)數(shù)隨長時間持續(xù)的變化情況。隨著飛行時間的持續(xù),結(jié)構(gòu)的溫度將會不斷地上升,其溫度變化率在起始時刻變化比較劇烈,但隨著時間的持續(xù),增幅將逐漸減小并漸趨于零;相對于結(jié)構(gòu)溫度的變化,壁面熱流密度隨著時間的持續(xù)會不斷地下降,其降幅也不斷減小。由此可以進(jìn)一步推知,隨著飛行時間的持續(xù),流場與結(jié)構(gòu)之間的耦合傳熱最終將會達(dá)到一個動態(tài)平衡狀態(tài)。通過上述分析可知,對于長時間持續(xù)飛行,即使馬赫數(shù)不是很高,結(jié)構(gòu)溫度也可能達(dá)到很高的值,飛行器防熱設(shè)計必須考慮飛行時間的累積效應(yīng)。

    圖9 駐點(diǎn)處溫度及其時間導(dǎo)數(shù)隨長時間持續(xù)變化Fig.9 Long-duration variation of the temperature and its timederivative at stagnation point

    圖10 駐點(diǎn)處熱流及其時間導(dǎo)數(shù)隨長時間持續(xù)變化Fig.10 Long-duration variation of the heat-flux and its timederivative at stagnation point

    圖11給出了外部流場和內(nèi)部結(jié)構(gòu)沿駐點(diǎn)線溫度分布隨長時間持續(xù)變化情況??梢灾庇^看出,外部流場的變化僅限在邊界層近壁處極小范圍內(nèi),流場溫度與結(jié)構(gòu)內(nèi)部溫度在壁面處連續(xù),外部流場氣動加熱通過壁面逐漸向結(jié)構(gòu)內(nèi)部傳遞,隨著時間的持續(xù),結(jié)構(gòu)內(nèi)部溫度逐漸增加。圖12給出了t=50s時外部流場和結(jié)構(gòu)內(nèi)部溫度分布云圖,可以直觀地看出結(jié)構(gòu)內(nèi)部傳熱具有明顯的多維特征,熱量除了主要沿壁面法向傳遞外,還會向其它方向傳遞。

    3.3 非穩(wěn)態(tài)飛行環(huán)境下熱環(huán)境特性分析

    高超聲速飛行器在實際飛行過程中一般會經(jīng)歷上升、巡航和下降等狀態(tài)。本文暫未考慮攻角的變化,根據(jù)氣動載荷限制條件給出了一條簡單的典型飛行軌跡,如圖13所示,總的持續(xù)飛行時間為410s,其中:上升時間為0~50s,巡航時間為50s~350s,下降時間為350s~410s。這里選取時間步長Δttr=5s將飛行軌跡離散劃分為82個離散飛行狀態(tài)點(diǎn)。對于本次計算,流場與結(jié)構(gòu)耦合傳熱的計算時間步長取為Δtcp=Δttr=5s;圓管前緣的材料采用C/C復(fù)合材料,其初始時刻的溫度為288.15K。

    圖14給出了來流總焓、圓管前緣駐點(diǎn)壁面溫度、氣動加熱熱流沿飛行軌跡的分布。可以看出:(1)從起始時刻t=0s至巡航狀態(tài)t=50s開始,結(jié)構(gòu)壁面溫度逐漸上升,但同時由于飛行速度的逐漸增大,流場的氣動加熱能力逐漸增強(qiáng),氣動加熱熱流密度并沒有降低,而是隨之逐漸增大。在這一過程中,通過耦合界面?zhèn)魅虢Y(jié)構(gòu)內(nèi)部的氣動加熱量越來越多,結(jié)構(gòu)處于不斷儲能的狀態(tài)。(2)在巡航狀態(tài)t=50s~350s,雖然飛行條件不再變化,流場的氣動加熱能力恒定,但由于流場與結(jié)構(gòu)之間的耦合傳熱作用,不斷有氣動加熱量持續(xù)地傳入結(jié)構(gòu)內(nèi)部,導(dǎo)致結(jié)構(gòu)壁面溫度會繼續(xù)上升,傳入結(jié)構(gòu)內(nèi)部的氣動加熱量越來越小,氣動加熱熱流密度開始逐漸下降,但結(jié)構(gòu)仍處于儲能狀態(tài)。如果巡航時間足夠長,隨著結(jié)構(gòu)內(nèi)部熱能的不斷儲存,最終流場與結(jié)構(gòu)之間的耦合傳熱將達(dá)到一個動態(tài)平衡狀態(tài)。(3)進(jìn)入下降狀態(tài)t=350s~410s,飛行速度開始降低,流場氣動加熱能力開始減弱,此時結(jié)構(gòu)溫度已經(jīng)很高,氣動熱流密度會繼續(xù)下降,外部氣動環(huán)境難以再繼續(xù)保持結(jié)構(gòu)的儲能狀態(tài),結(jié)構(gòu)開始轉(zhuǎn)為釋能狀態(tài),結(jié)構(gòu)內(nèi)部的熱能開始向外傳遞,從而結(jié)構(gòu)壁面溫度開始逐漸下降。

    圖15給出了結(jié)構(gòu)內(nèi)部溫度場沿飛行軌跡的分布,從中可以直觀看出結(jié)構(gòu)內(nèi)部溫度分布及熱能傳遞方向隨飛行軌跡的變化規(guī)律。由于流場與結(jié)構(gòu)之間的耦合傳熱效應(yīng),起初外部氣動環(huán)境使熱能不斷向結(jié)構(gòu)內(nèi)部傳遞,當(dāng)外部氣動環(huán)境難以維持結(jié)構(gòu)的儲能狀態(tài)時,結(jié)構(gòu)內(nèi)部熱能轉(zhuǎn)而向外傳遞。能夠準(zhǔn)確地預(yù)測和分析熱在結(jié)構(gòu)內(nèi)部的時空分布和傳遞規(guī)律,進(jìn)而加

    圖11 流場及結(jié)構(gòu)沿駐點(diǎn)線溫度分布隨長時間持續(xù)變化Fig.11 Long-duration variation of the temperature within the fluid and structure along stagnation line

    圖12 在t=50s時流場及結(jié)構(gòu)內(nèi)部溫度分布Fig.12 Temperaturedistribution within the fluid and structure at t=50s

    圖13 飛行軌跡及其離散劃分Fig.13 Flight trajectory anditsdiscretization

    圖14 來流總焓、駐點(diǎn)溫度和熱流密度沿飛行軌跡的分布Fig.14 Distribution of the fluid’s total enthalpy,temperature and heat-flux at the stagnation point along trajectory

    圖15 結(jié)構(gòu)內(nèi)部溫度場沿飛行軌跡的分布Fig.15 Distribution of the temperature field within the structure along trajectory

    以引導(dǎo)和控制,這對于新型防熱概念及方法的研究具有重要的意義。

    4 結(jié) 論

    通過對飛行器熱環(huán)境多場耦合特性的分析,提出了一種基于多場耦合的熱環(huán)境數(shù)值分析策略,建立了相應(yīng)的數(shù)值分析方法。以典型圓管前緣為計算模型,對穩(wěn)態(tài)和非穩(wěn)態(tài)飛行環(huán)境下的熱環(huán)境特性進(jìn)行了數(shù)值分析研究。結(jié)果分析表明:

    (1)充分考慮流場與結(jié)構(gòu)之間的耦合傳熱效應(yīng),對結(jié)構(gòu)溫度的預(yù)測更接近物理實際,可以進(jìn)一步減小設(shè)計冗余,而且這種效果會在氣動加熱越嚴(yán)重的地方越明顯;

    (2)對于長時間持續(xù)飛行,在穩(wěn)態(tài)飛行環(huán)境下,即使馬赫數(shù)不是很高,結(jié)構(gòu)溫度也可能達(dá)到很高的值,流場與結(jié)構(gòu)之間的耦合傳熱最終將達(dá)到一個動態(tài)平衡狀態(tài);在非穩(wěn)態(tài)飛行環(huán)境下,外部流場氣動加熱能力的變化,會改變結(jié)構(gòu)內(nèi)部熱能的分布和傳遞,引起熱流峰值或溫度峰值,飛行器熱環(huán)境隨飛行軌跡變化表現(xiàn)出動態(tài)的多場耦合時空特性。

    綜合以上分析,本文提出的基于多場耦合的飛行器熱環(huán)境數(shù)值分析策略以及所建立的數(shù)值分析方法,能夠有效地刻畫流場與結(jié)構(gòu)之間的耦合傳熱特征和規(guī)律,預(yù)測和分析飛行器熱環(huán)境空間和時間分布特性,可為防熱設(shè)計的選材和優(yōu)化提供可靠的參考依據(jù)和分析手段。同時值得指出的是,由于本文重點(diǎn)在于數(shù)值分析策略和方法的研究,暫未考慮壁面輻射特性的影響。實際上飛行器壁面輻射特性對氣動熱環(huán)境具有重要的影響,且隨著飛行時間的持續(xù)而不斷加強(qiáng),有待進(jìn)一步的分析和研究。

    [1]THORNTON E A,DECHAUMPHAI P.Coupled flow,thermal and structural analysis of aerodynamically heated panels[J].Journal of Aircraft,1988,25(11):1052-1059.doi:10.2514/3.45702

    [2]DECHAUMPHAI P,THORNTON E A,WIETING A R.Elowthermal-structural study of aerodynamically heated leading edges[J].Journal of Spacecraft,1989,26(4):201-209.doi:10.2514/3.26055

    [3]DECHAUMPH AI P,WIETING A R,PANDEY A K.Eluidthermal-structural interaction of aerodynamically heated leading edges[R].AIAA 89-1277,1989.doi:10.2514/6.1989-1227

    [4]CHEN Y K,HENLINE W D,TAUBER M E.Mars pathfinder trajectory based heating and ablation calculations[J].Journal of Spacecraft and Rockets,1995,32(2):225-230.doi:10.2514/3.26600

    [5]CHEN Y K,MILOS E S.Solution strategy for thermal response of non-ablating thermal protection systems at hypersonic speeds[R].AIAA 1996-0601.doi:10.2514/6.1996-615

    [6]CHEN Y K,MILOS E S.Two-dimensional implicit thermal response and ablation program for charring materials on hypersonic space vehicles[R].AIAA 2000-0206.doi:10.2514/6.2000-206

    [7]HUANG T,MAO G L,JANG G Q,et al.Twodimensional coupled flow-thermal-structural numerical simulation[J].ACTA Aerodynamica Sinica,2000,18(1):115-119.(in Chinese)

    黃唐,毛國良,姜貴慶,等.二維流場、熱、結(jié)構(gòu)一體化數(shù)值模擬[J].空氣動力學(xué)學(xué)報,2000,18(1):115-119.

    [8]GUI Y W,YUAN X J.Numerical simulation on the coupling phenomena of aerodynamic heating with thermal response in theregion of the leading edge[J].Journal of Engineering Thermophysics,2002,23(6):733-735.(in Chinese)

    桂業(yè)偉,袁湘江.類前緣防熱層流場與熱響應(yīng)耦合計算研究[J].工程熱物理學(xué)報,2002,23(6):733-735.

    [9]XIA G,LIU X J,CHENG W K,et al.Numerical simulation of coupled aeroheating and solid heat penetration for a hypersonic blunt body[J].Journal of National University of Defense Technology,2003,25(1):35-39.(in Chinese)

    夏剛,劉新建,程文科,等.鈍體高超聲速氣動加熱與結(jié)構(gòu)熱傳遞耦合的數(shù)值計算[J].國防科技大學(xué)學(xué)報,2003,25(1):35-39.

    [10]ZHAO X L,SUN Z X,TANG L S,et al.Coupledflow-thermal-structural analysis of hypersonic aerodynamically heated cylindrical leading edge[J].Engineering Applications of Computational Eluid Mechanics,2011,5(2):170-179.

    [11]NIE T,LIU W Q.Study of coupled fluid and solid for a hypersonic leading edge[J].ACTA Phys.Sin.,2012,61(18):184401-1-184401-7.(in Chinese)

    聶濤,劉偉強(qiáng).高超聲速飛行器前緣流固耦合計算方法研究[J].物理學(xué)報,2012,61(18):184401-1-184401-7.

    [12]LI X,ZHANG J E,HE Y L,et al.Numerical analysis of twodimensional fluid/thermal structure loosely-coupled simulation[J].Journal of Engineering Thermophysics,2012,33(1):87-90.(in Chinese)

    李欣,張劍飛,何雅玲,等.二維流場、熱結(jié)構(gòu)松耦合模擬研究[J].工程熱物理學(xué)報,2012,33(1):87-90.

    [13]KIM K H,KIM C.Accurate,efficient and monotonic numerical methods for multi-dimensional compressible flows part I:spatialdiscretization[J].Journal of Computational Physics, 2005,28(2):527-569.doi:10.1016/j.jcp.2005.02.021

    [14]YOON S H,KIM C,KIM K H.Multi-dimensional limiting process for three-dimensional flow physics analyses[J].Journal of Computational Physics,2008,227:6001-6043.doi:10.1016/j.jcp.2008.02.012

    [15]YOON S,JAMESON A.Lower-upper symmetry gauss-seidel method for the Euler and Navier-stokes equations[J].AIAA Journal,1988,26(9):1025-1026.doi:10.2514/3.10007

    [16]GILES M.Stability analysis of numerical interface conditions in fluid-structure thermal analysis[J].International Journal of Numerical Methods in Eluids,1997,25(4):421-436.doi:10.1002/(SICI)1097-0363(19970830)25:4<421:AID-ELD557>3.0.CO;2-J.

    [17]SHEPARD D.A twodimensional interpolation function for irregularly spaceddata[C].Proceedings of 23rd National Conference of the Association for Computing Machinery,New York,1968:517-524.

    [18]WIETING A R.Experimental study of shock wave interference heating on a cylindrical leading edge[R].NASA TM100484,1987.

    [19]CHENG X L,AI B C,WANG Q.A wall grid scale criterion based on the molecule mean free path for the wall heat flux computations by the Navier-Stokes equations[J].Chinese Journal of Theoretical and Applied Mechanics,2010,42(6):1083-1089.(in Chinese)

    程曉麗,艾邦成,王強(qiáng).基于分子平均自由程的熱流計算壁面網(wǎng)格準(zhǔn)則[J].力學(xué)學(xué)報,2010,42(6):1083-1089.

    Multi-field coupling numerical analysis approach for aerothermal environment of hypersonic vehicles

    ZHANG Shengtao,CHEN Eang,LIU Hong
    (School of Aeronautics and Astronautics,Shanghai Jiao Tong University,Shanghai 200240,China)

    Thedevelopment of new-generation hypersonic vehicles presents a major challenge in thedesign of thermal protection systems.Sustained hypersonic flight within the atmosphere can result in severe aerodynamic heating phenomena.It is a physical fact that significant interaction occurs between the external aerodynamic heating and the structural heat transfer within the vehicles.Through analyzing the multi-field coupling characteristics of aerothermal environment of hypersonic vehicles,a multi-field coupling numerical analysis approach for predicting aerothermal environment is proposed in this paper.This approach couples the computational fluiddynamics(CED)codes based on the Navier-Stokes equations with the general commercial finite element method(EEM)software by the reliable interfacialdata exchange method.Considering a cylindrical leading edge as test case,the fluid-structural thermal coupling characteristics along the steady and unsteady flight trajectory is numerically investigated.It is indicated that the proposed approach could accurately predict the fluid-structural thermal coupling characteristics,and achieve the analysis of spatial and temporaldistribution characteristics of aerothermal environment,thus providing the reliable analysis tool for the material selection and optimization in thedesign of thermal protection systems.

    multi-field coupling;long-duration aeroheating;aerothermal environment;thermal protectiondesign;hypersonic vehicles

    V211.22;V211.3

    Adoi:10.7638/kqdlxxb-2012.0220

    0258-1825(2014)06-0861-07

    2012-12-27;

    2013-05-26

    國家自然科學(xué)基金項目(11102111)

    張勝濤(1983-),男,山東聊城人,博士研究生,研究方向為高超聲速氣動熱數(shù)值模擬及熱防護(hù).E-mail:zhangst9656@163.com

    課題負(fù)責(zé)人:陳方(1977-),男,博士,副研究員,上海交通大學(xué)航空航天學(xué)院1326室.E-mail:fangchen@sjtu.edu.cn

    張勝濤,陳方,劉洪.基于多場耦合的飛行器熱環(huán)境數(shù)值分析方法研究[J].空氣動力學(xué)學(xué)報,2014,32(6):861-867.

    10.7638/kqdlxxb-2012.0220 ZH ANG S T,CHEN E,LIU H.Multi-field coupling numerical analysis approach for aerothermal environment of hypersonic vehicles[J].ACTA Aerodynamica Sinica,2014,32(6):861-867.

    猜你喜歡
    熱流壁面步長
    二維有限長度柔性壁面上T-S波演化的數(shù)值研究
    基于Armijo搜索步長的BFGS與DFP擬牛頓法的比較研究
    內(nèi)傾斜護(hù)幫結(jié)構(gòu)控釋注水漏斗熱流道注塑模具
    空調(diào)溫控器上蓋熱流道注塑模具設(shè)計
    聚合物微型零件的熱流固耦合變形特性
    中國塑料(2017年2期)2017-05-17 06:13:24
    壁面溫度對微型內(nèi)燃機(jī)燃燒特性的影響
    基于逐維改進(jìn)的自適應(yīng)步長布谷鳥搜索算法
    透明殼蓋側(cè)抽模熱流道系統(tǒng)的設(shè)計
    中國塑料(2014年5期)2014-10-17 03:02:17
    一種新型光伏系統(tǒng)MPPT變步長滯環(huán)比較P&O法
    電測與儀表(2014年2期)2014-04-04 09:04:00
    顆?!诿媾鲎步Ec數(shù)據(jù)處理
    一区二区三区高清视频在线| 可以在线观看的亚洲视频| 亚洲精品日韩av片在线观看| 亚洲成人久久爱视频| 午夜激情福利司机影院| 亚洲国产色片| 午夜亚洲福利在线播放| 97超视频在线观看视频| 久久久久久久久大av| 精品国产三级普通话版| 国产极品精品免费视频能看的| 免费av不卡在线播放| 国语自产精品视频在线第100页| 久久久久国产网址| 久久久久久久久久久免费av| 国产精品乱码一区二三区的特点| 国产伦理片在线播放av一区 | 久久草成人影院| 亚洲人成网站在线播| 亚洲国产精品久久男人天堂| 国产一级毛片在线| 尤物成人国产欧美一区二区三区| 黑人高潮一二区| 草草在线视频免费看| 国产成人91sexporn| 国产蜜桃级精品一区二区三区| 一区福利在线观看| 97超视频在线观看视频| 99热网站在线观看| 亚州av有码| 日韩成人av中文字幕在线观看| 99久久无色码亚洲精品果冻| 久久久色成人| 国产中年淑女户外野战色| 久久精品国产99精品国产亚洲性色| 99久久成人亚洲精品观看| 成人欧美大片| 欧美+日韩+精品| 不卡视频在线观看欧美| 你懂的网址亚洲精品在线观看 | 国产精品蜜桃在线观看 | 看免费成人av毛片| 成人特级黄色片久久久久久久| 日本-黄色视频高清免费观看| 悠悠久久av| 欧美日韩一区二区视频在线观看视频在线 | av在线播放精品| 性色avwww在线观看| 99久久精品国产国产毛片| 精品人妻熟女av久视频| 婷婷色av中文字幕| 亚洲国产精品久久男人天堂| 精品国内亚洲2022精品成人| 欧美潮喷喷水| 熟妇人妻久久中文字幕3abv| 三级男女做爰猛烈吃奶摸视频| 美女脱内裤让男人舔精品视频 | 亚洲成人精品中文字幕电影| 成人永久免费在线观看视频| 久久久国产成人免费| 男女啪啪激烈高潮av片| 搞女人的毛片| 99久久九九国产精品国产免费| 中国美白少妇内射xxxbb| 国产欧美日韩精品一区二区| 久久亚洲国产成人精品v| 久久久久久久久大av| 中文精品一卡2卡3卡4更新| 久久精品久久久久久噜噜老黄 | 久久久精品大字幕| 国产精品99久久久久久久久| 日本免费一区二区三区高清不卡| 国产不卡一卡二| 亚洲欧洲日产国产| 日韩亚洲欧美综合| 日本在线视频免费播放| 在线播放无遮挡| 身体一侧抽搐| 亚洲国产欧美在线一区| 日本一二三区视频观看| 亚洲av成人精品一区久久| 国产黄片美女视频| 久久热精品热| 亚洲18禁久久av| 校园人妻丝袜中文字幕| 欧美丝袜亚洲另类| 亚洲人成网站在线观看播放| 欧美日韩国产亚洲二区| 日韩,欧美,国产一区二区三区 | 天堂av国产一区二区熟女人妻| 中文字幕久久专区| 亚洲一区高清亚洲精品| 夜夜夜夜夜久久久久| 国产中年淑女户外野战色| 啦啦啦啦在线视频资源| 麻豆av噜噜一区二区三区| 免费在线观看成人毛片| 简卡轻食公司| 乱系列少妇在线播放| 我的女老师完整版在线观看| 日韩精品青青久久久久久| 亚洲精品乱码久久久久久按摩| 全区人妻精品视频| 九草在线视频观看| 日日摸夜夜添夜夜添av毛片| 日日摸夜夜添夜夜添av毛片| 少妇的逼水好多| 成年av动漫网址| 十八禁国产超污无遮挡网站| 欧洲精品卡2卡3卡4卡5卡区| 特大巨黑吊av在线直播| 欧洲精品卡2卡3卡4卡5卡区| 成人一区二区视频在线观看| 国内精品久久久久精免费| 少妇的逼水好多| 一个人看视频在线观看www免费| 伦理电影大哥的女人| 少妇熟女aⅴ在线视频| 亚洲色图av天堂| 精品久久久久久久人妻蜜臀av| 伦理电影大哥的女人| 国产精品无大码| 精品久久国产蜜桃| 国产成人a∨麻豆精品| 一本久久中文字幕| 久久综合国产亚洲精品| 欧美一区二区精品小视频在线| 变态另类成人亚洲欧美熟女| 国产成人freesex在线| 精品久久久噜噜| 伦精品一区二区三区| 国产精品乱码一区二三区的特点| 婷婷亚洲欧美| 欧美三级亚洲精品| 亚洲激情五月婷婷啪啪| 国产不卡一卡二| 欧美一区二区国产精品久久精品| 卡戴珊不雅视频在线播放| 国产一区亚洲一区在线观看| 国产日韩欧美在线精品| 十八禁国产超污无遮挡网站| 久久久欧美国产精品| 久久欧美精品欧美久久欧美| 久久久久久久久大av| 嫩草影院新地址| 欧美xxxx黑人xx丫x性爽| 欧美色欧美亚洲另类二区| 国产高清三级在线| 国产黄片美女视频| 国产乱人视频| 哪里可以看免费的av片| 一本久久精品| 国产精品久久久久久久久免| 啦啦啦观看免费观看视频高清| 国产成人午夜福利电影在线观看| 亚洲成人久久性| 3wmmmm亚洲av在线观看| 男女边吃奶边做爰视频| 一区二区三区免费毛片| 国产片特级美女逼逼视频| 尤物成人国产欧美一区二区三区| 麻豆成人av视频| 日本-黄色视频高清免费观看| av国产免费在线观看| 26uuu在线亚洲综合色| 国产精品久久久久久精品电影| 我要搜黄色片| 18禁裸乳无遮挡免费网站照片| 桃色一区二区三区在线观看| 久久欧美精品欧美久久欧美| 国产精品一区二区三区四区久久| 亚洲高清免费不卡视频| 日本黄色片子视频| 国产综合懂色| 少妇熟女aⅴ在线视频| 一进一出抽搐gif免费好疼| 2021天堂中文幕一二区在线观| 久久婷婷人人爽人人干人人爱| 免费看美女性在线毛片视频| 久久精品夜色国产| 99久久成人亚洲精品观看| 久久精品国产亚洲av涩爱 | 国产av不卡久久| 毛片一级片免费看久久久久| 亚洲欧美日韩高清专用| 亚洲不卡免费看| 伦精品一区二区三区| 日韩人妻高清精品专区| 日韩欧美三级三区| 毛片一级片免费看久久久久| 国产精品一二三区在线看| 99久久人妻综合| 干丝袜人妻中文字幕| 亚洲国产欧美在线一区| 超碰av人人做人人爽久久| 久久精品国产自在天天线| 深夜a级毛片| 精品久久国产蜜桃| 18禁裸乳无遮挡免费网站照片| 久久99热6这里只有精品| 精品国内亚洲2022精品成人| 一本精品99久久精品77| 国产精品乱码一区二三区的特点| 亚洲激情五月婷婷啪啪| 亚洲精品色激情综合| 99久久无色码亚洲精品果冻| 久久久国产成人精品二区| 亚洲人成网站在线观看播放| 亚洲中文字幕一区二区三区有码在线看| av免费观看日本| 美女大奶头视频| 我的女老师完整版在线观看| av视频在线观看入口| 干丝袜人妻中文字幕| 亚洲av中文字字幕乱码综合| 日本-黄色视频高清免费观看| 热99在线观看视频| 我要搜黄色片| 亚洲欧美日韩高清专用| 天堂影院成人在线观看| 久久久久免费精品人妻一区二区| 少妇熟女aⅴ在线视频| 国产精品福利在线免费观看| 国产精品日韩av在线免费观看| 国产亚洲精品久久久com| 国产av一区在线观看免费| 国产成人freesex在线| 欧美日韩一区二区视频在线观看视频在线 | 赤兔流量卡办理| 美女脱内裤让男人舔精品视频 | 国产免费一级a男人的天堂| 亚洲精品粉嫩美女一区| 99热6这里只有精品| 免费观看精品视频网站| 大型黄色视频在线免费观看| 青青草视频在线视频观看| 中出人妻视频一区二区| 亚洲,欧美,日韩| 男人舔女人下体高潮全视频| 在线观看午夜福利视频| av在线老鸭窝| 国产精品综合久久久久久久免费| 久久久精品欧美日韩精品| 如何舔出高潮| 亚洲经典国产精华液单| 久久久成人免费电影| 欧美色视频一区免费| 国产午夜精品一二区理论片| 亚洲国产精品sss在线观看| 日日摸夜夜添夜夜爱| 久久久a久久爽久久v久久| 两个人的视频大全免费| 婷婷六月久久综合丁香| 亚洲精品影视一区二区三区av| 午夜精品国产一区二区电影 | 欧美高清成人免费视频www| 22中文网久久字幕| 国国产精品蜜臀av免费| 久久久久网色| 午夜久久久久精精品| .国产精品久久| 久久亚洲精品不卡| 日日干狠狠操夜夜爽| 亚洲精品456在线播放app| 日韩欧美在线乱码| 天堂中文最新版在线下载 | 在线天堂最新版资源| 九九久久精品国产亚洲av麻豆| 欧美日韩在线观看h| 少妇人妻一区二区三区视频| 国产精品免费一区二区三区在线| 久久久精品欧美日韩精品| 国产成人aa在线观看| 尤物成人国产欧美一区二区三区| 国产高清不卡午夜福利| 久久精品久久久久久噜噜老黄 | 亚洲人成网站高清观看| 国产精品电影一区二区三区| 97超碰精品成人国产| 国产成人精品一,二区 | av在线观看视频网站免费| 成人三级黄色视频| 久久亚洲精品不卡| 国产精品一区二区三区四区久久| 久久久久久久午夜电影| 少妇裸体淫交视频免费看高清| 欧美激情在线99| 国模一区二区三区四区视频| 久久精品国产鲁丝片午夜精品| 国产精品1区2区在线观看.| 哪里可以看免费的av片| 国产麻豆成人av免费视频| 国产久久久一区二区三区| av天堂中文字幕网| 国产精品,欧美在线| 亚洲四区av| 精品国内亚洲2022精品成人| 亚洲电影在线观看av| 深夜精品福利| 国产亚洲av嫩草精品影院| a级毛色黄片| 少妇的逼水好多| 国产亚洲5aaaaa淫片| 国产精品99久久久久久久久| .国产精品久久| 日本色播在线视频| 一个人免费在线观看电影| 成人鲁丝片一二三区免费| 亚洲欧美成人精品一区二区| 午夜福利在线在线| 中国国产av一级| 精品一区二区三区视频在线| 欧美一级a爱片免费观看看| 美女脱内裤让男人舔精品视频 | 日产精品乱码卡一卡2卡三| 国产成年人精品一区二区| 日韩亚洲欧美综合| 老师上课跳d突然被开到最大视频| 一级黄色大片毛片| 亚洲精品日韩在线中文字幕 | 国产亚洲精品久久久com| 国产av一区在线观看免费| 夫妻性生交免费视频一级片| 99久久精品热视频| 亚洲婷婷狠狠爱综合网| 国产色婷婷99| 成人毛片60女人毛片免费| 性插视频无遮挡在线免费观看| kizo精华| 一级毛片久久久久久久久女| 成人午夜精彩视频在线观看| 日日撸夜夜添| 国产亚洲av片在线观看秒播厂 | 最近最新中文字幕大全电影3| 国产亚洲5aaaaa淫片| 91久久精品国产一区二区成人| 不卡视频在线观看欧美| 又粗又硬又长又爽又黄的视频 | 一进一出抽搐gif免费好疼| 国产亚洲精品久久久com| 欧美激情久久久久久爽电影| 在线观看美女被高潮喷水网站| 在线播放国产精品三级| 国产精品久久久久久亚洲av鲁大| 久久精品国产亚洲网站| 亚洲国产高清在线一区二区三| 亚洲精品久久国产高清桃花| 亚洲欧美日韩东京热| 欧美性猛交黑人性爽| 免费黄网站久久成人精品| 97人妻精品一区二区三区麻豆| 亚洲人与动物交配视频| 亚洲色图av天堂| 午夜亚洲福利在线播放| 国产伦理片在线播放av一区 | 久久午夜福利片| 国内精品一区二区在线观看| 深爱激情五月婷婷| 97超碰精品成人国产| a级毛片a级免费在线| 国产色爽女视频免费观看| 在线观看一区二区三区| 99热这里只有精品一区| 成人三级黄色视频| 麻豆av噜噜一区二区三区| 亚洲欧美精品专区久久| 一进一出抽搐动态| 91久久精品电影网| 国产精品一区www在线观看| 亚洲电影在线观看av| 久久韩国三级中文字幕| 人人妻人人澡人人爽人人夜夜 | 熟女电影av网| 97热精品久久久久久| 伦理电影大哥的女人| 观看美女的网站| 国产av麻豆久久久久久久| 男女边吃奶边做爰视频| 老女人水多毛片| 国内少妇人妻偷人精品xxx网站| 最后的刺客免费高清国语| 人妻制服诱惑在线中文字幕| 午夜亚洲福利在线播放| 国产免费男女视频| 超碰av人人做人人爽久久| 搞女人的毛片| 亚洲最大成人av| 精品国内亚洲2022精品成人| 五月玫瑰六月丁香| 观看免费一级毛片| 亚洲av.av天堂| 亚洲不卡免费看| 哪个播放器可以免费观看大片| 国产中年淑女户外野战色| 欧美三级亚洲精品| 久久韩国三级中文字幕| 亚洲最大成人中文| 69人妻影院| 久久久久久久久中文| 久久精品影院6| 一本精品99久久精品77| 久久久成人免费电影| 欧美成人免费av一区二区三区| 成人特级黄色片久久久久久久| 国产男人的电影天堂91| 热99re8久久精品国产| 国产伦理片在线播放av一区 | 一级黄片播放器| 国产女主播在线喷水免费视频网站 | 三级男女做爰猛烈吃奶摸视频| 一夜夜www| 久久亚洲国产成人精品v| h日本视频在线播放| 在线观看午夜福利视频| 美女脱内裤让男人舔精品视频 | 少妇人妻精品综合一区二区 | 久久人人爽人人爽人人片va| 天堂影院成人在线观看| 日韩制服骚丝袜av| 精品99又大又爽又粗少妇毛片| 久久婷婷人人爽人人干人人爱| 女的被弄到高潮叫床怎么办| 村上凉子中文字幕在线| 在线免费观看不下载黄p国产| 国产高潮美女av| 国产亚洲av嫩草精品影院| 欧美日韩综合久久久久久| 少妇猛男粗大的猛烈进出视频 | 中文精品一卡2卡3卡4更新| 国产大屁股一区二区在线视频| 成人国产麻豆网| 欧美丝袜亚洲另类| 国产精品一区二区三区四区免费观看| 美女国产视频在线观看| 亚洲最大成人av| 久久精品国产亚洲av香蕉五月| 国产蜜桃级精品一区二区三区| 夜夜爽天天搞| 国产精品麻豆人妻色哟哟久久 | 久久亚洲精品不卡| 午夜福利视频1000在线观看| 精品久久久噜噜| av在线蜜桃| 色播亚洲综合网| 日韩成人伦理影院| 久久这里只有精品中国| 国产精品美女特级片免费视频播放器| 日本黄大片高清| 一边摸一边抽搐一进一小说| 一级黄色大片毛片| 性色avwww在线观看| 亚洲在久久综合| 我要搜黄色片| 国产精品99久久久久久久久| 一级二级三级毛片免费看| 久久久a久久爽久久v久久| 一区福利在线观看| 亚洲不卡免费看| 26uuu在线亚洲综合色| 久久久国产成人精品二区| 一进一出抽搐gif免费好疼| 免费观看人在逋| 国产伦精品一区二区三区四那| 男人舔奶头视频| 免费观看在线日韩| 五月伊人婷婷丁香| 欧美成人一区二区免费高清观看| 波多野结衣巨乳人妻| 亚洲av男天堂| 亚洲人成网站在线播| 亚洲精品久久久久久婷婷小说 | 日韩人妻高清精品专区| 午夜福利在线观看吧| 99久久精品一区二区三区| 国产黄a三级三级三级人| 级片在线观看| 精品久久久久久久久亚洲| 99热全是精品| 我的老师免费观看完整版| 欧美在线一区亚洲| 青春草视频在线免费观看| 菩萨蛮人人尽说江南好唐韦庄 | 日产精品乱码卡一卡2卡三| 日韩欧美精品免费久久| 日本熟妇午夜| 在线免费观看的www视频| 欧美丝袜亚洲另类| 国产单亲对白刺激| 熟女电影av网| 亚洲国产精品成人综合色| 天堂av国产一区二区熟女人妻| 精品一区二区免费观看| 如何舔出高潮| 日本五十路高清| 亚洲国产色片| 亚洲经典国产精华液单| 一级毛片久久久久久久久女| 非洲黑人性xxxx精品又粗又长| 国产老妇女一区| 日本欧美国产在线视频| 1000部很黄的大片| 亚洲国产精品成人久久小说 | 1024手机看黄色片| 国产女主播在线喷水免费视频网站 | 国产在线男女| 久久久精品94久久精品| 日韩成人伦理影院| 中文字幕熟女人妻在线| 91精品一卡2卡3卡4卡| 26uuu在线亚洲综合色| 婷婷色av中文字幕| 成人三级黄色视频| 国产精品久久久久久久久免| 免费观看a级毛片全部| 别揉我奶头 嗯啊视频| 高清日韩中文字幕在线| 亚洲人成网站在线播| 久久国产乱子免费精品| 老司机福利观看| 日韩欧美 国产精品| 91在线精品国自产拍蜜月| 国产精品人妻久久久久久| 精品久久久久久久久亚洲| 97热精品久久久久久| 亚洲在线观看片| 五月伊人婷婷丁香| 99热全是精品| 亚洲成人精品中文字幕电影| 日本黄色片子视频| 日本色播在线视频| 午夜精品国产一区二区电影 | 精品久久久久久久久亚洲| 亚洲av电影不卡..在线观看| 草草在线视频免费看| 国产大屁股一区二区在线视频| 五月玫瑰六月丁香| 99精品在免费线老司机午夜| 天堂中文最新版在线下载 | 99热全是精品| 欧美日本亚洲视频在线播放| 久久国内精品自在自线图片| 免费av不卡在线播放| 亚洲七黄色美女视频| 成人一区二区视频在线观看| 久久久久久久久大av| 亚洲精品456在线播放app| 男人舔奶头视频| 国产精品一区www在线观看| 午夜福利在线观看免费完整高清在 | 亚洲av一区综合| 国产午夜精品一二区理论片| 久久人人爽人人爽人人片va| av在线老鸭窝| 久久中文看片网| 成人欧美大片| 免费电影在线观看免费观看| 日韩视频在线欧美| 婷婷色综合大香蕉| 床上黄色一级片| 国模一区二区三区四区视频| 波多野结衣巨乳人妻| 床上黄色一级片| 亚洲国产日韩欧美精品在线观看| 久久久久久国产a免费观看| 国产色婷婷99| 免费观看a级毛片全部| 久久精品国产亚洲网站| 精品久久久久久久人妻蜜臀av| 国产黄a三级三级三级人| 国产在视频线在精品| 久久99蜜桃精品久久| 99在线视频只有这里精品首页| 国产老妇女一区| 一个人免费在线观看电影| 亚洲美女搞黄在线观看| 日韩中字成人| 成人鲁丝片一二三区免费| 亚洲美女视频黄频| 久久精品国产亚洲网站| 人妻夜夜爽99麻豆av| 在现免费观看毛片| 在线天堂最新版资源| 六月丁香七月| 黄色配什么色好看| 观看美女的网站| 久久久久国产网址| 99久国产av精品国产电影| 在线天堂最新版资源| 69人妻影院| 少妇熟女aⅴ在线视频| 亚洲精品亚洲一区二区| av在线天堂中文字幕| 十八禁国产超污无遮挡网站| 国产 一区精品| 亚洲久久久久久中文字幕| 色综合亚洲欧美另类图片| 99久久成人亚洲精品观看| 国产成人91sexporn| 国产亚洲av片在线观看秒播厂 | 22中文网久久字幕| 亚洲欧美日韩高清专用| 国产女主播在线喷水免费视频网站 | 亚洲av熟女| 国产伦理片在线播放av一区 | 可以在线观看的亚洲视频| ponron亚洲| 最近的中文字幕免费完整| 在线观看66精品国产| 午夜老司机福利剧场| 人体艺术视频欧美日本| 淫秽高清视频在线观看| 成人午夜精彩视频在线观看| 国产成人福利小说| 亚洲人成网站在线观看播放| 国产精品日韩av在线免费观看|