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

    渦輪燃?xì)庵髁?盤(pán)腔二次流的分區(qū)耦合仿真I:算法及平臺(tái)實(shí)現(xiàn)*

    2015-05-31 11:14:16中國(guó)航空工業(yè)集團(tuán)公司沈陽(yáng)發(fā)動(dòng)機(jī)設(shè)計(jì)研究所周曉菲康曉恩
    航空制造技術(shù) 2015年12期
    關(guān)鍵詞:輪緣湍流分區(qū)

    中國(guó)航空工業(yè)集團(tuán)公司沈陽(yáng)發(fā)動(dòng)機(jī)設(shè)計(jì)研究所 于 霄 周曉菲 康曉恩

    重慶大學(xué)動(dòng)力工程學(xué)院 葉 建 武芏茳

    為滿足軍用和民用飛行器對(duì)其動(dòng)力裝置提出的持續(xù)攀升的性能指標(biāo)要求,過(guò)去幾十年間,航空發(fā)動(dòng)機(jī)渦輪部件的進(jìn)口溫度不斷提高[1-2],高溫燃?xì)馊肭洲D(zhuǎn)靜盤(pán)腔可能帶來(lái)的負(fù)面效應(yīng)也日趨嚴(yán)重,由此導(dǎo)致燃?xì)庵髁骱捅P(pán)腔二次流的耦合作用問(wèn)題受到越來(lái)越多研究者的關(guān)注。在研究渦輪部件氣動(dòng)熱力學(xué)的3種基本方法中,數(shù)值模擬手段的重要性不言而喻,而已有的研究表明[3-4]:輪緣間隙附近冷氣和燃?xì)獾膿交爝^(guò)程非常復(fù)雜,RANS模擬難以準(zhǔn)確捕捉,URANS結(jié)果與試驗(yàn)更為接近,但仍有很大差距,LES則是準(zhǔn)確預(yù)測(cè)該類(lèi)流動(dòng)的可能選項(xiàng)之一。雖然RANS模擬的預(yù)測(cè)精度有限,但從實(shí)踐的角度看,該方法幾乎仍是工程設(shè)計(jì)中可選用的唯一方案,提高RANS對(duì)燃?xì)庵髁?盤(pán)腔二次流耦合問(wèn)題的預(yù)測(cè)能力仍然有著重要的工程應(yīng)用價(jià)值。盤(pán)腔內(nèi)部的二次流動(dòng)與渦輪主流是2種性質(zhì)截然不同的流動(dòng)[5],適用于每種流動(dòng)的“最佳”湍流模型顯然是不同的,但目前已有的耦合模擬幾乎都只用了一個(gè)模型。如果可以針對(duì)這類(lèi)流動(dòng)發(fā)展一個(gè)分區(qū)耦合算法,在主流和二次流區(qū)域分別選擇各自“最佳”的湍流模型,有可能會(huì)提高RANS方法的預(yù)測(cè)精度。

    構(gòu)建針對(duì)燃?xì)庵髁?盤(pán)腔二次流的分區(qū)耦合平臺(tái),可以歸結(jié)到復(fù)雜多物理場(chǎng)耦合這一范疇。隨著計(jì)算機(jī)運(yùn)算能力的不斷增強(qiáng)以及人們對(duì)工程和科學(xué)問(wèn)題預(yù)測(cè)精度提出的更高要求,越來(lái)越多的研究人員認(rèn)識(shí)到:大多數(shù)的科學(xué)問(wèn)題和工程實(shí)際問(wèn)題往往處于多物理場(chǎng)的復(fù)雜耦合環(huán)境中,此時(shí),單憑某一軟件實(shí)現(xiàn)這類(lèi)問(wèn)題的高精度預(yù)測(cè)變得越來(lái)越困難,解決方案之一即是進(jìn)行多物理場(chǎng)的耦合仿真。但是,目前大多數(shù)的軟件產(chǎn)品還只能處理某一類(lèi)問(wèn)題或者只對(duì)該類(lèi)問(wèn)題具有最佳的預(yù)測(cè)精度,要實(shí)現(xiàn)多場(chǎng)耦合,其現(xiàn)實(shí)途徑就是通過(guò)某種技術(shù)平臺(tái)將多個(gè)軟件產(chǎn)品集成到一起,以期完成復(fù)雜多物理場(chǎng)的耦合模擬。

    工業(yè)界的迫切需求,以及該研究方向取得突破可能帶來(lái)的巨大收益使得國(guó)內(nèi)外研究者,尤其是國(guó)外的CAE巨頭開(kāi)展了大量的相關(guān)工作,并已經(jīng)取得了較為豐碩的成果。歸納起來(lái),現(xiàn)有的耦合平臺(tái)大體上可分為3類(lèi):(1)商業(yè)化的耦合平臺(tái),如德國(guó)SCAI研究中心的MpCCI[6]和美國(guó)Ansys公司的Workbench[7]; (2)開(kāi)源的耦合平臺(tái),如法國(guó)研究機(jī)構(gòu)CERFACS 開(kāi)發(fā)的 OpenPALM[8],美國(guó)Sandia國(guó)家實(shí)驗(yàn)室開(kāi)發(fā)的LIME[9]等;(3)基于Python等腳本語(yǔ)言的自編程耦合平臺(tái)[10-11]。

    要實(shí)現(xiàn)燃?xì)庵髁?盤(pán)腔二次流的分區(qū)耦合模擬,上述3類(lèi)平臺(tái)均可以使用,但考慮到編程的靈活性、平臺(tái)調(diào)試和修改的方便程度以及計(jì)算的精度和速度等因素,本文采用第三種方案——即基于Python語(yǔ)言構(gòu)建耦合平臺(tái)。下文首先從分區(qū)界面的定義、界面的數(shù)據(jù)交換2個(gè)方面進(jìn)行分區(qū)耦合的算法設(shè)計(jì),而后從耦合計(jì)算的執(zhí)行流程、數(shù)據(jù)交換的實(shí)現(xiàn)、Fluent執(zhí)行的自動(dòng)化和整體計(jì)算流程的控制4個(gè)方面介紹了耦合平臺(tái)的構(gòu)建。

    圖1 高壓渦輪主流和盤(pán)腔二次流系統(tǒng)示意圖

    分區(qū)耦合算法設(shè)計(jì)

    對(duì)于本文所研究的燃?xì)庵髁?盤(pán)腔二次流分區(qū)耦合仿真問(wèn)題,在構(gòu)建耦合平臺(tái)之前,有以下關(guān)鍵問(wèn)題需要解決:(1)燃?xì)庵髁骱捅P(pán)腔二次流作為一個(gè)在空間連續(xù)分布、且有著復(fù)雜相互作用的流體系統(tǒng),分區(qū)的界面應(yīng)如何選取; (2)確定分區(qū)界面后,在進(jìn)行耦合計(jì)算時(shí),主流區(qū)和二次流區(qū)域如何耦合在一起,這些問(wèn)題所涉及的內(nèi)容就是分區(qū)耦合的算法設(shè)計(jì),它既是構(gòu)建耦合平臺(tái)的前提,也是核心的技術(shù)難點(diǎn)。本文通過(guò)分析燃?xì)庵髁骱捅P(pán)腔二次流所具有的主要耦合特性,從分區(qū)界面的定義出發(fā),探討合理的界面選取方案;而后討論當(dāng)2個(gè)區(qū)域分別使用不同的湍流模型時(shí),分區(qū)界面應(yīng)交換哪些數(shù)據(jù)、如何進(jìn)行交換,以此實(shí)現(xiàn)高效的耦合。

    1 分區(qū)界面的定義

    高壓渦輪主流和盤(pán)腔二次流所組成的流動(dòng)系統(tǒng)如圖1所示,其復(fù)雜性主要表現(xiàn)在3個(gè)方面:(1)主流流道和盤(pán)腔的幾何結(jié)構(gòu)非常復(fù)雜,導(dǎo)致高質(zhì)量計(jì)算網(wǎng)格的生成難度較大;(2)導(dǎo)葉葉排和動(dòng)葉葉排包括葉盤(pán)之間的相對(duì)運(yùn)動(dòng),對(duì)計(jì)算域的劃分以及算法的選取都提出了特殊的要求;(3)多種因素綜合作用導(dǎo)致流場(chǎng)的流動(dòng)結(jié)構(gòu)非常復(fù)雜。

    針對(duì)單獨(dú)的主流或盤(pán)腔流動(dòng)問(wèn)題,目前的仿真方案已經(jīng)相對(duì)成熟。本文中只考慮穩(wěn)態(tài)定常模擬時(shí),對(duì)于單級(jí)渦輪,許多商業(yè)軟件均能生成高質(zhì)量的計(jì)算網(wǎng)格,靜子和轉(zhuǎn)子的界面通過(guò)摻混面方法實(shí)現(xiàn)流場(chǎng)數(shù)據(jù)的傳遞;對(duì)于轉(zhuǎn)靜盤(pán)腔,其網(wǎng)格大多數(shù)情況下需要手動(dòng)生成,但計(jì)算中通過(guò)設(shè)置恰當(dāng)?shù)谋诿孢吔鐥l件可以模擬轉(zhuǎn)盤(pán)的影響。如若將兩者耦合在一起,則需要解決一系列的難題。

    首先,耦合計(jì)算域中間存在靜子/轉(zhuǎn)子的摻混面以及主流區(qū)/二次流區(qū)的分界面,這兩個(gè)界面的關(guān)系是什么,其相對(duì)位置的選取如何。從摻混面的概念出發(fā),兩個(gè)界面應(yīng)該分開(kāi)、避免相交,也就是說(shuō),盤(pán)腔的計(jì)算域只是與轉(zhuǎn)子或者靜子的計(jì)算域存在交接面,而不會(huì)同時(shí)與兩者相交,這意味著盤(pán)腔的網(wǎng)格或者與靜子網(wǎng)格固結(jié)在一起不動(dòng),或者和轉(zhuǎn)子網(wǎng)格固結(jié)在一起轉(zhuǎn)動(dòng)。從已有的結(jié)果看,不同研究者選取的方案并不一致,2種方案都有一定的合理性,其選擇也可能與輪緣間隙的軸向相對(duì)位置有關(guān)。

    其次,燃?xì)庵髁骱捅P(pán)腔二次流的分區(qū)界面應(yīng)如何選取?從上一段的分析可知,盤(pán)腔計(jì)算域只與靜子或者轉(zhuǎn)子計(jì)算域中的一個(gè)交接,這一問(wèn)題就簡(jiǎn)化為:對(duì)于靜子/盤(pán)腔網(wǎng)格或者轉(zhuǎn)子/盤(pán)腔網(wǎng)格,二者的分界面取在何處。從主流和二次流耦合系統(tǒng)的幾何結(jié)構(gòu)看,2個(gè)大的幾何空間通過(guò)一個(gè)狹小的通道——輪緣間隙連接在一起,理論上,分區(qū)界面只能設(shè)定在該區(qū)域,如若不然,界面將可能穿過(guò)存在復(fù)雜流動(dòng)結(jié)構(gòu)的主流或二次流區(qū)域,這將給界面數(shù)據(jù)交換帶來(lái)大的困難。實(shí)際上,即便在輪緣間隙區(qū),也可能存在著較為復(fù)雜的流動(dòng),但考慮到定常模擬時(shí),該區(qū)域的時(shí)均流動(dòng)相對(duì)簡(jiǎn)單,因而將分界面取在該處是可行的;至于其具體位置,可取在輪緣間隙區(qū)中部,界面方向則和兩個(gè)壁面盡量垂直。

    再次,耦合系統(tǒng)的計(jì)算網(wǎng)格的生成??紤]到通過(guò)分區(qū)界面,兩個(gè)區(qū)域需要頻繁地進(jìn)行數(shù)據(jù)交換,為保證數(shù)據(jù)傳遞的精度,最好在界面兩側(cè)使用完全匹配的計(jì)算網(wǎng)格。由此可采用如下的計(jì)算域劃分策略:在主流區(qū)創(chuàng)建一個(gè)和盤(pán)腔輪緣間隙匹配的區(qū)域,而后確定分區(qū)界面,在其兩側(cè)生成完全匹配的計(jì)算網(wǎng)格。

    基于上述分析,燃?xì)庵髁髋c盤(pán)腔二次流分區(qū)界面的定義應(yīng)遵循如下原則:(1)盤(pán)腔的網(wǎng)格或者與靜子網(wǎng)格、或者與轉(zhuǎn)子網(wǎng)格固結(jié)在一起,二者的交接面即是主流和二次流的分區(qū)界面,該分區(qū)界面不與靜子/轉(zhuǎn)子的摻混面相交。(2)主流和二次流的分區(qū)界面應(yīng)位于具有最小流通面積的輪緣間隙區(qū)域,界面的朝向與兩側(cè)壁面盡量垂直,保證最小的截面積,通過(guò)該界面,流動(dòng)的方向應(yīng)盡量保持一致。(3)分區(qū)界面兩側(cè)的計(jì)算網(wǎng)格應(yīng)保證完全匹配。

    2 分區(qū)界面的數(shù)據(jù)交換

    從已有的研究可知:在定常雷諾平均模擬中,對(duì)于渦輪燃?xì)庵髁?,使用S-A模型的綜合效果最優(yōu),對(duì)于盤(pán)腔二次流,則推薦SSTk-ω模型,據(jù)此,我們嘗試在主流區(qū)選用S-A模型、二次流區(qū)選用SSTk-ω模型進(jìn)行分區(qū)耦合。

    根據(jù)上一小節(jié)的原則定義分區(qū)界面后,就可以考慮界面間的數(shù)據(jù)交換了。理想情況下,為實(shí)現(xiàn)2個(gè)區(qū)域間的強(qiáng)耦合,可以將分區(qū)界面作為一個(gè)“內(nèi)邊界”,2個(gè)區(qū)域使用不同的控制方程分別推進(jìn),每計(jì)算一步,沿著內(nèi)邊界進(jìn)行一次數(shù)據(jù)的交換,交換的數(shù)據(jù)既包括u、v、w、p、ρ、T等基本的流場(chǎng)變量,也包括湍流粘性系數(shù)μt。由于兩側(cè)求解的湍流輸運(yùn)方程不同,還應(yīng)該通過(guò)界面重構(gòu)湍流輸運(yùn)量:對(duì)于主流區(qū)而言,已知二次流邊界上的μt、k、ω?cái)?shù)值,容易得到S-A模型的輸運(yùn)量,反過(guò)來(lái)對(duì)于二次流區(qū)域而言,已知主流邊界上的μt和,要得到k和ω則較為困難(見(jiàn)圖2)。

    圖2 主流和盤(pán)腔二次流之間的數(shù)據(jù)交換

    實(shí)際的情形與此不同,由于本文選用商業(yè)軟件Fluent作為求解器,我們對(duì)其內(nèi)核缺少了解,更無(wú)法修改源代碼,兩個(gè)區(qū)域的耦合只能通過(guò)基于Python的耦合平臺(tái)并結(jié)合Fluent提供的用戶定義函數(shù)UDF實(shí)現(xiàn),顯然,這是一種“松耦合”:兩個(gè)Fluent程序被同時(shí)執(zhí)行,分別計(jì)算主流流動(dòng)和二次流動(dòng),分區(qū)界面是其“外邊界”而非“內(nèi)邊界”,分區(qū)耦合平臺(tái)利用UDF重設(shè)邊界條件實(shí)現(xiàn)區(qū)域間的數(shù)據(jù)交換。這樣帶來(lái)的問(wèn)題是:我們無(wú)法像強(qiáng)耦合中的內(nèi)邊界一樣在界面交換所有流場(chǎng)數(shù)據(jù),而只能根據(jù)界面處的流動(dòng)情況設(shè)定恰當(dāng)?shù)倪吔鐥l件,并根據(jù)邊條的要求提供部分流場(chǎng)變量,以此實(shí)現(xiàn)界面的數(shù)據(jù)傳遞。

    (1)界面處流動(dòng)方向的設(shè)定。

    對(duì)于在實(shí)際工況運(yùn)行的主流/二次流系統(tǒng)而言,輪緣間隙附近的流動(dòng)異常復(fù)雜,受間隙兩側(cè)當(dāng)?shù)鼐植苛鲃?dòng)參數(shù)的影響,周向不同位置處間隙徑向流動(dòng)的方向很可能是變化的:有的地方從盤(pán)腔進(jìn)入主流,有的地方則從主流進(jìn)入盤(pán)腔。但在Fluent中,分區(qū)界面作為一個(gè)確定的外部邊界,必須預(yù)先給定邊條的類(lèi)型,或者是進(jìn)口,或者是出口,其前提就是要明確界面處的流動(dòng)方向。

    考慮下述3個(gè)因素,我們將分區(qū)耦合界面的流動(dòng)方向設(shè)定為從盤(pán)腔流向主流:盤(pán)腔冷氣的主要作用是封住輪緣間隙,防止外部高溫燃?xì)馇秩?,在正常工作狀態(tài)下,盡管沿不同周向位置可能有局部瞬態(tài)的燃?xì)馇秩?,但從平均效果看,盤(pán)腔冷氣還是往外流入主流區(qū)的;本文只對(duì)耦合系統(tǒng)進(jìn)行定常仿真,由于未計(jì)入非定常效應(yīng),周向局部可能存在的燃?xì)馊肭脂F(xiàn)象被大大弱化,當(dāng)盤(pán)腔二次流入口給定的總壓相對(duì)較高時(shí),即便存在回流,分區(qū)界面處流動(dòng)的主要方向還是指向燃?xì)庵髁饕粋?cè)的;設(shè)定分區(qū)界面的流動(dòng)方向以后,F(xiàn)luent的計(jì)算依然允許在界面處存在與設(shè)定方向相反的流動(dòng),因而回流的影響可以得到考慮。

    (2)界面處流動(dòng)參數(shù)的傳遞。

    要到達(dá)主流二次流耦合的效果,必須保證沿著分區(qū)界面流場(chǎng)的質(zhì)量、動(dòng)量和能量守恒,在Fluent中,分區(qū)界面作為一個(gè)“外邊界”,并不能直接交換基本流場(chǎng)變量,而只能進(jìn)行邊界條件的設(shè)定,借助UDF,我們能夠?qū)崿F(xiàn)邊界上每個(gè)網(wǎng)格單元一一對(duì)應(yīng)的邊條信息傳遞,待流場(chǎng)迭代收斂后達(dá)到和交換基本流場(chǎng)變量類(lèi)似的效果,最終保證界面處三大守恒律的滿足。

    確定分區(qū)耦合界面的流動(dòng)方向是從盤(pán)腔流向主流后,對(duì)于盤(pán)腔流動(dòng),分區(qū)界面為計(jì)算域出口,設(shè)定壓力出口邊條;對(duì)于燃?xì)庵髁鳎謪^(qū)界面為計(jì)算域入口,設(shè)定速度進(jìn)口邊條,具體的設(shè)置如下:對(duì)于盤(pán)腔流動(dòng),分區(qū)界面處設(shè)定壓力出口邊條,需要主流區(qū)提供的參數(shù)包括背壓、回流總溫和回流的方向向量;對(duì)于主流燃?xì)?,分區(qū)界面設(shè)定為速度進(jìn)口邊條,需要二次流區(qū)提供的參數(shù)包括3個(gè)速度分量,溫度和回流時(shí)的背壓。盡管沿著分區(qū)界面?zhèn)鬟f的邊條參數(shù)是有限的,但這些參數(shù)已經(jīng)保證了計(jì)算結(jié)果的適定性,通過(guò)分區(qū)界面?zhèn)鬟f數(shù)據(jù),2個(gè)計(jì)算域反復(fù)迭代直至流場(chǎng)收斂,此時(shí)界面兩側(cè)的流場(chǎng)變量是連續(xù)分布的,也就保證了分區(qū)界面處三大守恒律的滿足。

    (3)界面處湍流參數(shù)的傳遞。

    通過(guò)界面處流動(dòng)參數(shù)的傳遞,保證了沿著分區(qū)界面兩側(cè)流場(chǎng)變量是連續(xù)分布的,分區(qū)界面也滿足質(zhì)量、動(dòng)量、能量三大守恒律,但這并不能保證耦合仿真結(jié)果的預(yù)測(cè)精度,其原因在于:對(duì)基于雷諾平均的仿真方法而言,湍流模型是決定預(yù)測(cè)成敗的關(guān)鍵,僅僅在分區(qū)界面兩側(cè)保證流場(chǎng)變量連續(xù)是遠(yuǎn)遠(yuǎn)不夠的。當(dāng)流動(dòng)從一個(gè)區(qū)域穿過(guò)分區(qū)界面流入另一個(gè)區(qū)域時(shí),湍流量也會(huì)隨流進(jìn)入該區(qū)域,即在計(jì)算域入口(或計(jì)算域出口的回流區(qū)),除了需要流動(dòng)參數(shù)的邊界信息外,還應(yīng)該提供湍流參數(shù)的信息。

    主流區(qū)和二次流區(qū)域計(jì)算時(shí)使用的湍流模型不同,分別是S-A和SSTk-ω模型,S-A模型求解關(guān)于變量 的一個(gè)輸運(yùn)方程,SSTk-ω模型則求解關(guān)于湍動(dòng)能k和湍流比耗散率ω的兩個(gè)輸運(yùn)方程。無(wú)論S-A還是SSTk-ω模型,均是基于Boussinesq的渦粘性假設(shè),顯然,界面湍流參數(shù)的傳遞應(yīng)該滿足這樣的原則:傳遞后沿界面兩側(cè)湍流粘性系數(shù)μt的分布連續(xù)。

    基于上述原則,當(dāng)流動(dòng)從盤(pán)腔流入主流時(shí),對(duì)于主流區(qū)(求解S-A模型)而言,分區(qū)界面的速度進(jìn)口邊條要求二次流區(qū)域(求解SSTk-ω模型)提供一個(gè)湍流參數(shù),這個(gè)參數(shù)可以是修改后的湍流粘性,也可以是湍流粘性比μt/μ。從兩方程模型降為一方程模型,這是很容易實(shí)現(xiàn)的,考慮到UDF中可以直接讀出界面單元的湍流粘性μt和分子粘性μ,這里選用湍流粘性比μt/μ作為傳遞的參數(shù),具體執(zhí)行時(shí),F(xiàn)luent根據(jù)邊界UDF給定的μt/μ反算出湍流粘性,用于輸運(yùn)方程的求解。

    當(dāng)流動(dòng)從主流區(qū)回流盤(pán)腔,對(duì)于二次流區(qū)域(求解SSTk-ω模型)而言,分區(qū)界面的壓力出口邊條要求主流區(qū)(求解S-A模型)提供2個(gè)湍流參數(shù),分別是湍動(dòng)能k和湍流比耗散率ω,這意味著要從一方程模型重構(gòu)二方程模型的參數(shù),其實(shí)現(xiàn)難度很大。該問(wèn)題可描述為:以保證界面處湍流粘性系數(shù)連續(xù)分布為原則,將湍流粘性系數(shù)μt或湍流粘性比μt/μ作為參數(shù)傳遞(事實(shí)上這也是主流區(qū)的唯一選擇),實(shí)現(xiàn)對(duì)界面處湍動(dòng)能k和湍流比耗散率ω的重構(gòu)。事實(shí)上,對(duì)于k、ω、ε和μt,我們知道兩個(gè)關(guān)系式:

    其中,ρ是密度,經(jīng)驗(yàn)系數(shù),給定μt后,3個(gè)未知量k、ω、ε,2個(gè)方程,無(wú)法求解,必須補(bǔ)充1個(gè)關(guān)系式或者給定3個(gè)未知量中的1個(gè)。

    由于并不存在一般的準(zhǔn)則,要解決這一問(wèn)題,必須對(duì)界面所在位置——輪緣間隙處的流動(dòng)進(jìn)行具體分析,該流動(dòng)具有2個(gè)主要的特點(diǎn):沿著整個(gè)圓周,輪緣間隙是一個(gè)很窄的狹縫,子午面內(nèi)狹縫的長(zhǎng)度有限;輪緣間隙轉(zhuǎn)盤(pán)側(cè)相對(duì)于間隙靜盤(pán)側(cè)以很高的速度運(yùn)動(dòng),該速度遠(yuǎn)遠(yuǎn)大于間隙在子午面內(nèi)的流動(dòng)速度。

    據(jù)此,如圖3所示,可以將間隙內(nèi)的流動(dòng)分解為2部分的矢量和,其一是周向的Couette流,一側(cè)平面的速度很大;其二是子午面內(nèi)的Poiseuille流,該流動(dòng)的速度很小。通過(guò)輪緣間隙,無(wú)論進(jìn)入還是離開(kāi)盤(pán)腔,大部分流體質(zhì)點(diǎn)在子午面內(nèi)運(yùn)動(dòng)的距離很短,但該過(guò)程中均沿周向運(yùn)動(dòng)了較長(zhǎng)的距離,那么可以認(rèn)為,流體質(zhì)點(diǎn)穿過(guò)輪緣間隙的流動(dòng),近似為一個(gè)平面Couette流動(dòng)。

    圖3 Couette流和Poiseuille流示意圖

    我們將耦合界面設(shè)置在輪緣間隙中間位置,顯然該位置也處于Couette流動(dòng)中部,Couette流是一種相對(duì)簡(jiǎn)單的剪切流,若流動(dòng)是湍流態(tài),則湍動(dòng)能的生成和耗散近似平衡,更近一步,可以假設(shè)界面處的湍流流動(dòng)處于平衡態(tài),即有如下關(guān)系:

    其中,湍動(dòng)能的生成項(xiàng)為:

    其中,是速度張量, 是梯度張量。引入Boussinesq的渦粘性模型有:

    其中,Sij是張量。由上面的等式可以知道,耗散率ε計(jì)算公式為:

    知道ε以后,利用公式(1)、(2)可以算出k、ω。顯然,這些參數(shù)應(yīng)該在主流區(qū)進(jìn)行計(jì)算,而后將k、ω傳遞到二次流區(qū)域。

    分區(qū)耦合平臺(tái)構(gòu)建

    通過(guò)上一節(jié)對(duì)分區(qū)耦合算法的研究,我們對(duì)如何定義分區(qū)界面,耦合過(guò)程中沿分區(qū)界面?zhèn)鬟f哪些數(shù)據(jù)等問(wèn)題有了較為深刻的理解和認(rèn)識(shí),接下來(lái)需要考慮的問(wèn)題是:如何構(gòu)建分區(qū)耦合平臺(tái)以成功實(shí)現(xiàn)燃?xì)庵髁髋c盤(pán)腔二次流的耦合仿真。本節(jié)從分析耦合計(jì)算的執(zhí)行流程入手,概括了耦合平臺(tái)的2大主要工作,而后分析界面數(shù)據(jù)交換的實(shí)現(xiàn)方法以及Fluent程序的自動(dòng)化執(zhí)行,在此基礎(chǔ)上完成了耦合平臺(tái)的流程框圖并介紹了平臺(tái)的Python實(shí)現(xiàn)。

    1 耦合計(jì)算的執(zhí)行流程

    設(shè)計(jì)分區(qū)耦合平臺(tái)的重要前提之一是要搞清楚分區(qū)耦合過(guò)程的執(zhí)行流程,現(xiàn)在設(shè)想人工執(zhí)行這一流程的步驟如下(假設(shè)計(jì)算網(wǎng)格已準(zhǔn)備好):

    (1)運(yùn)行Fluent程序A,讀入燃?xì)庵髁鞯挠?jì)算網(wǎng)格,進(jìn)行相關(guān)設(shè)置(湍流模型、邊界條件等等),完成主流流場(chǎng)初始化;再運(yùn)行一個(gè)Fluent程序B,讀入盤(pán)腔二次流的計(jì)算網(wǎng)格,進(jìn)行相關(guān)設(shè)置(湍流模型、邊界條件等等),完成二次流流場(chǎng)初始化。

    (2)設(shè)定計(jì)算步數(shù),程序A(主流計(jì)算)進(jìn)行指定步數(shù)的迭代求解;設(shè)定計(jì)算步數(shù),程序B(二次流計(jì)算)也同樣進(jìn)行指定步數(shù)的迭代求解。

    (3)程序A(主流計(jì)算)輸出分區(qū)界面壓力p,總溫Tt,湍流粘性μt等數(shù)據(jù);程序B(二次流計(jì)算)輸出分區(qū)界面速度分量u、v、w,溫度T以及和μt等數(shù)據(jù)。

    (4)處理程序B的輸出數(shù)據(jù),更新程序A(主流計(jì)算)分區(qū)界面的邊界條件;處理程序A的輸出數(shù)據(jù),更新程序B(二次流計(jì)算)分區(qū)界面的邊界條件。

    (5)跳轉(zhuǎn)到第(2)步,(2)、(3)、(4)、(5)步順序執(zhí)行。

    (6)若流場(chǎng)滿足收斂條件,程序A、程序B分別輸出計(jì)算結(jié)果,而后分別結(jié)束。

    從上述步驟容易看出,要通過(guò)程序代碼將耦合計(jì)算的執(zhí)行流程自動(dòng)化,耦合平臺(tái)的主要工作包括2大部分:

    第一是流程的控制,即控制2個(gè)Fluent程序的啟動(dòng),初始化(讀入網(wǎng)格、進(jìn)行相關(guān)的設(shè)置),迭代求解,收斂判定,以及計(jì)算結(jié)果的輸出,程序結(jié)束等等,這一部分的重要問(wèn)題是保證兩個(gè)軟件的同步,即是說(shuō)迭代求解快的程序等待迭代求解慢的程序,實(shí)現(xiàn)界面數(shù)據(jù)傳遞的同步。

    第二是完成分區(qū)界面的數(shù)據(jù)交換,兩個(gè)Fluent程序分別迭代計(jì)算一定步數(shù)后,沿著分區(qū)界面分別輸出對(duì)方程序需要的數(shù)據(jù),對(duì)輸出的數(shù)據(jù)進(jìn)行恰當(dāng)?shù)奶幚恚?個(gè)Fluent程序再分別讀入處理后的數(shù)據(jù),更新分區(qū)界面的邊界條件。

    2 數(shù)據(jù)交換的實(shí)現(xiàn)

    通過(guò)上一小節(jié)搞清楚了分區(qū)耦合計(jì)算的執(zhí)行流程后,我們首先考慮如何實(shí)現(xiàn)數(shù)據(jù)的交換,因?yàn)橄鄬?duì)于整體的流程控制而言,它屬于“局部”的問(wèn)題,對(duì)流程依賴不多,可以“獨(dú)立”的解決。

    分析發(fā)現(xiàn),無(wú)論對(duì)于燃?xì)庵髁骰蛘弑P(pán)腔二次流動(dòng),耦合計(jì)算過(guò)程中,除分區(qū)界面處邊界條件的數(shù)值需要不斷更新外,其他的邊界條件均無(wú)需做任何改動(dòng);而界面處邊條的類(lèi)型并不發(fā)生變化,只是需要更新具體的數(shù)值,并且這些數(shù)值不是單個(gè)的數(shù)據(jù),而是界面上的一組數(shù)據(jù)分布。

    利用Fluent的用戶自定義函數(shù)UDF,可以輕松實(shí)現(xiàn)上述更新。圖4給出了耦合邊界數(shù)據(jù)交換的示意圖,燃?xì)庵髁鱂luent程序A以及盤(pán)腔二次流Fluent程序B完成一定步數(shù)的迭代后,開(kāi)始交換數(shù)據(jù),利用 UDF,程序A將界面上的相關(guān)信息輸出到aout.dat文件,程序B則將界面信息輸出到bout.dat文件,基于Python的主控程序讀入aout.dat、bout.dat,按照要求對(duì)其中的數(shù)據(jù)進(jìn)行處理,而后分別生成程序B和程序A的輸入文件 bin.dat、ain.dat,接下來(lái)同樣利用UDF,程序A和B分別讀入ain.dat和bin.dat,完成界面邊界條件數(shù)據(jù)的更新。這樣就實(shí)現(xiàn)了一次分區(qū)界面的數(shù)據(jù)交換。

    圖4 耦合邊界的數(shù)據(jù)交換

    3 Fluent執(zhí)行的自動(dòng)化

    一般而言,我們?cè)谑褂肍luent進(jìn)行仿真工作時(shí),通常會(huì)啟動(dòng)Fluent的圖形用戶界面(簡(jiǎn)稱(chēng)GUI),以便實(shí)現(xiàn)人與程序的實(shí)時(shí)交互。但是,在某些情況下,我們可能對(duì)交互性要求不高,只是希望程序能自動(dòng)完成一系列的計(jì)算工作,F(xiàn)luent提供了命令行的運(yùn)行方式,合理利用它很容易實(shí)現(xiàn)程序執(zhí)行的自動(dòng)化。

    該過(guò)程主要是通過(guò)Fluent 的一個(gè)命令集合即進(jìn)程文件(Journal File)實(shí)現(xiàn)的,Journal文件可以重播用戶曾經(jīng)進(jìn)行的所有操作,其創(chuàng)建途徑有2個(gè):一是在用戶進(jìn)入圖形用戶界面后,系統(tǒng)自動(dòng)記錄用戶的操作和命令輸入,自動(dòng)生成進(jìn)程文件;另一個(gè)則是用戶使用文本編輯器直接用Scheme 語(yǔ)言創(chuàng)建進(jìn)程文件。

    對(duì)于本項(xiàng)目而言,主控程序需要“同時(shí)”反復(fù)調(diào)用2個(gè)Fluent程序,分別進(jìn)行燃?xì)庵髁骱捅P(pán)腔二次流的迭代計(jì)算,這一過(guò)程中,需要修改的參數(shù)很少,沒(méi)有必要也基本上無(wú)法調(diào)用Fluent的GUI進(jìn)行相關(guān)設(shè)置。在主控程序執(zhí)行前,我們可以利用Scheme語(yǔ)言分別創(chuàng)建好主流計(jì)算和二次流計(jì)算的Journal文件,然后主控程序分別調(diào)用這2個(gè)文件,進(jìn)行迭代求解,如果有相關(guān)計(jì)算參數(shù)需要修改,則主控程序直接讀入Journal文件進(jìn)行修改再保存即可,這樣就完全實(shí)現(xiàn)了Fluent迭代過(guò)程的自動(dòng)化。

    4 整體計(jì)算流程的控制及實(shí)現(xiàn)

    通過(guò)上述3個(gè)小節(jié)的描述將相關(guān)細(xì)節(jié)討論的比較清楚之后,接下來(lái)考慮整體計(jì)算流程的控制,設(shè)想我們用某種工具或語(yǔ)言編寫(xiě)耦合計(jì)算程序,參考圖5,按照計(jì)算流程該平臺(tái)需要完成的主要工作如下:

    圖5 耦合平臺(tái)的流程框圖

    (1)計(jì)算啟動(dòng),平臺(tái)初始化,這一階段需要考慮多種情況,是進(jìn)行全新的計(jì)算,還是讀入已有的結(jié)果進(jìn)行續(xù)算?如有需要,主控程序分別讀入主流和二次流計(jì)算的Journal文件,進(jìn)行相關(guān)部分的修改然后存盤(pán)。

    (2)主控程序調(diào)用Fluent分別執(zhí)行兩個(gè)Journal,進(jìn)行耦合計(jì)算第n次的迭代求解。

    (3)主控程序判斷主流和二次流的第n次迭代求解是否均已完成,若都已完成,則分別讀入兩個(gè)程序的輸出文件aout.dat和bout.dat,按照要求對(duì)其中的數(shù)據(jù)進(jìn)行處理,再分別輸出兩個(gè)程序的輸入文件bin.dat和ain.dat。

    (4)主控程序判斷計(jì)算是否收斂,若已收斂,跳到第(5)步執(zhí)行,否則對(duì)兩個(gè)Journal文件進(jìn)行適當(dāng)?shù)男薷?,返回第?)步執(zhí)行。

    (5)耦合計(jì)算已收斂,將相關(guān)結(jié)果存盤(pán),主控程序結(jié)束。

    本文的耦合平臺(tái)選用開(kāi)源的面向?qū)ο蟮腜ython 語(yǔ)言編寫(xiě),由于該語(yǔ)言具有良好的可移植性,在避免使用依賴于系統(tǒng)特性的前提下,耦合平臺(tái)幾乎無(wú)需修改就可以在任何系統(tǒng)操作系統(tǒng)上運(yùn)行。作為一種近似解釋型語(yǔ)言,Python相對(duì)于編譯型語(yǔ)言運(yùn)行速度較慢,但耦合平臺(tái)要實(shí)現(xiàn)的只是接口和主控功能,在耦合計(jì)算過(guò)程中,絕大部分時(shí)間是“被耦合”的Fluent軟件在進(jìn)行計(jì)算,耦合程序真正執(zhí)行的時(shí)間是很短的,因而它執(zhí)行的快慢幾乎不會(huì)影響整體的耦合計(jì)算效率。

    結(jié) 論

    為實(shí)現(xiàn)渦輪燃?xì)庵髁骱捅P(pán)腔二次流的分區(qū)耦合仿真,本文從分析該耦合系統(tǒng)所具有的主要流動(dòng)特征和耦合特性入手,探討了分區(qū)耦合算法的設(shè)計(jì)及分區(qū)耦合平臺(tái)的構(gòu)建,主要結(jié)論如下:

    (1)燃?xì)庵髁髋c盤(pán)腔二次流分區(qū)界面的定義應(yīng)遵循3個(gè)原則:盤(pán)腔的網(wǎng)格與靜子或轉(zhuǎn)子網(wǎng)格固結(jié)在一起,二者的界面不與摻混面相交;分區(qū)界面應(yīng)位于輪緣間隙區(qū)域,通過(guò)該界面的流動(dòng)方向盡量保持一致;分區(qū)界面兩側(cè)的計(jì)算網(wǎng)格應(yīng)保證完全匹配。

    (2)分區(qū)界面處的流動(dòng)方向設(shè)定為從盤(pán)腔指向主流,界面參數(shù)的傳遞結(jié)合UDF通過(guò)設(shè)定邊界條件實(shí)現(xiàn),既應(yīng)滿足界面兩側(cè)流場(chǎng)的質(zhì)量、動(dòng)量和能量守恒,也應(yīng)保證湍流粘性系數(shù)的連續(xù)分布。

    (3)耦合平臺(tái)的2大主要功能包括實(shí)現(xiàn)對(duì)計(jì)算流程的控制以及完成分區(qū)界面的數(shù)據(jù)交換,用Python語(yǔ)言搭建平臺(tái),結(jié)合Fluent的UDF功能和命令行執(zhí)行,可以實(shí)現(xiàn)燃?xì)庵髁骱捅P(pán)腔二次流分區(qū)耦合仿真提出的全部要求。

    [1]Bogard D G, Thole K A. Gas turbine film cooling. Journal of Propulsion and Power,2006, 22: 249-270.

    [2]劉大響, 金捷. 21世紀(jì)世界航空動(dòng)力技術(shù)發(fā)展趨勢(shì)與展望. 中國(guó)工程科學(xué),2004, 6(9): 1-8.

    [3]Valencia A G, Dixon J A, Soghe R D,et al. An investigation into numerical analysis alternatives for predicting re-Ingestion in turbine disc rim cavities. ASME Paper GT2012-68592,2012.

    [4]Mahoney O, Hills T S D, Chew N J, et al. Large-Eddy Simulation of Rim Seal Ingestion.ASME GT2010-22962, 2010.

    [5]Chew J W, Hills N J. Computational fluid dynamics for turbomachinery internal air systems. Phil. Trans. R. Soc. A, 2007, 365: 2587-2611.

    [6]Wolf K. MpCCI-The General Code Coupling Interface, LS-DYNA Forum,Frankenthal, 2007.

    [7]ANSYS Inc, Workbench User's Guide,Canonsburg:AHSYS Inc, 2013.

    [8]Piacentini A, Morel T, The'venin A, et al. Open-PALM: An open source dynamic parallel coupler. Coupled problems, Greece, 2011.

    [9]Pawlowski R, Bartlett R, Belcourt N,et al. A Theory Manual for Multi-Physics Code Coupling in LIME//Sandia Technical Report SAND2011-2195. Albuquerque:Sandia National Laboratories, 2011.

    [10]Schluter J U, Wu X, Weide E, et al.A python approach to multi-code simulations CHIMPS//Center for Turbulence Research,Annual Research Briefs, 2005:97-110..

    [11]Wang P, Zheng Y, Zou Z, et al. A novel multi-fidelity coupled simulation method for flow systems. Chinese Journal of Aeronautics,2013,26(4): 868-875.

    猜你喜歡
    輪緣湍流分區(qū)
    上海實(shí)施“分區(qū)封控”
    淺談液態(tài)和固態(tài)輪緣潤(rùn)滑裝置的差異性
    地鐵車(chē)輛輪緣厚度偏磨問(wèn)題研究
    重氣瞬時(shí)泄漏擴(kuò)散的湍流模型驗(yàn)證
    關(guān)于優(yōu)化四方平臺(tái)動(dòng)車(chē)組輪對(duì)踏面旋修的研究
    浪莎 分區(qū)而治
    干式輪緣潤(rùn)滑器對(duì)地鐵車(chē)輛車(chē)輪保護(hù)效果的研究
    基于SAGA聚類(lèi)分析的無(wú)功電壓控制分區(qū)
    基于多種群遺傳改進(jìn)FCM的無(wú)功/電壓控制分區(qū)
    “青春期”湍流中的智慧引渡(三)
    男男h啪啪无遮挡| 精品一区二区三区视频在线| 99热6这里只有精品| 国产深夜福利视频在线观看| 日日摸夜夜添夜夜添av毛片| 欧美日韩综合久久久久久| 国产午夜精品久久久久久一区二区三区| 亚洲欧美一区二区三区黑人 | 18禁在线无遮挡免费观看视频| 国产精品99久久99久久久不卡 | 国产亚洲欧美精品永久| 夜夜骑夜夜射夜夜干| 亚洲国产成人一精品久久久| 高清视频免费观看一区二区| 久久久午夜欧美精品| 高清av免费在线| 我要看黄色一级片免费的| 国产欧美日韩一区二区三区在线 | 久久精品久久精品一区二区三区| 男人爽女人下面视频在线观看| 高清毛片免费看| 久久鲁丝午夜福利片| 欧美日韩成人在线一区二区| 男女高潮啪啪啪动态图| 精品人妻熟女av久视频| 午夜影院在线不卡| 夫妻性生交免费视频一级片| 日韩精品有码人妻一区| 黄色怎么调成土黄色| 建设人人有责人人尽责人人享有的| 日本欧美视频一区| 亚洲第一区二区三区不卡| 亚洲伊人久久精品综合| 少妇猛男粗大的猛烈进出视频| 久久99热6这里只有精品| 两个人免费观看高清视频| 我的老师免费观看完整版| 国产黄色视频一区二区在线观看| 午夜精品国产一区二区电影| a级片在线免费高清观看视频| 三级国产精品欧美在线观看| 爱豆传媒免费全集在线观看| 插阴视频在线观看视频| 国产日韩欧美亚洲二区| 波野结衣二区三区在线| 国产伦精品一区二区三区视频9| 成人漫画全彩无遮挡| 观看美女的网站| 久久精品国产亚洲网站| 日韩一区二区视频免费看| 成人毛片a级毛片在线播放| 中文精品一卡2卡3卡4更新| 欧美 亚洲 国产 日韩一| 老女人水多毛片| .国产精品久久| 搡女人真爽免费视频火全软件| xxxhd国产人妻xxx| 美女视频免费永久观看网站| 国产成人freesex在线| 成人黄色视频免费在线看| 黄片播放在线免费| 啦啦啦啦在线视频资源| 久久av网站| 亚洲精品,欧美精品| 免费大片18禁| 蜜桃国产av成人99| 99re6热这里在线精品视频| 哪个播放器可以免费观看大片| 国产高清国产精品国产三级| 韩国av在线不卡| 在线播放无遮挡| 国产精品国产三级专区第一集| 一级毛片电影观看| 一级黄片播放器| 国产不卡av网站在线观看| 婷婷色麻豆天堂久久| 狂野欧美激情性xxxx在线观看| 国产成人aa在线观看| 久热这里只有精品99| 久久热精品热| 久久久精品94久久精品| 国产精品国产三级国产专区5o| 久热这里只有精品99| 国产黄色免费在线视频| 成人18禁高潮啪啪吃奶动态图 | 久久久久久久久久人人人人人人| 最新中文字幕久久久久| 少妇猛男粗大的猛烈进出视频| 久久久久久久大尺度免费视频| 亚洲无线观看免费| 国产精品无大码| av国产精品久久久久影院| 纯流量卡能插随身wifi吗| 久久久久久久久久久免费av| 女人精品久久久久毛片| 美女cb高潮喷水在线观看| 国产免费又黄又爽又色| 成人毛片a级毛片在线播放| 精品一区在线观看国产| 韩国高清视频一区二区三区| 国产探花极品一区二区| 最新中文字幕久久久久| 亚洲综合精品二区| 久久国产精品男人的天堂亚洲 | 九色成人免费人妻av| 熟女人妻精品中文字幕| 国产精品不卡视频一区二区| 一区在线观看完整版| 日日摸夜夜添夜夜爱| 久久精品国产亚洲av天美| 精品国产露脸久久av麻豆| 99久久中文字幕三级久久日本| 一区二区三区乱码不卡18| 国产成人精品一,二区| 成年美女黄网站色视频大全免费 | 亚洲精品久久成人aⅴ小说 | 高清欧美精品videossex| 男人添女人高潮全过程视频| a级毛片在线看网站| 色哟哟·www| 97精品久久久久久久久久精品| 免费看不卡的av| 中文欧美无线码| 美女福利国产在线| 欧美精品高潮呻吟av久久| 交换朋友夫妻互换小说| 在线观看国产h片| 国产精品久久久久久精品古装| 欧美人与善性xxx| 日韩精品有码人妻一区| 2021少妇久久久久久久久久久| 国产精品麻豆人妻色哟哟久久| 狠狠婷婷综合久久久久久88av| 国产精品国产三级国产av玫瑰| 夜夜骑夜夜射夜夜干| 亚洲丝袜综合中文字幕| 一区二区三区精品91| 免费观看无遮挡的男女| 少妇精品久久久久久久| 日韩av不卡免费在线播放| 美女xxoo啪啪120秒动态图| 如何舔出高潮| 国产熟女欧美一区二区| 免费久久久久久久精品成人欧美视频 | 一级片'在线观看视频| 精品一区二区三卡| 午夜免费男女啪啪视频观看| 日韩强制内射视频| 精品一区在线观看国产| 99热这里只有是精品在线观看| 久久久a久久爽久久v久久| 日本-黄色视频高清免费观看| av免费观看日本| 日本爱情动作片www.在线观看| 久久精品夜色国产| 少妇被粗大的猛进出69影院 | 亚洲三级黄色毛片| 国产免费福利视频在线观看| 91国产中文字幕| 亚洲国产精品一区三区| 在线免费观看不下载黄p国产| 国产一区有黄有色的免费视频| 另类亚洲欧美激情| 国产av精品麻豆| 欧美人与性动交α欧美精品济南到 | 国产在线视频一区二区| 我的女老师完整版在线观看| 少妇 在线观看| 看免费成人av毛片| 免费不卡的大黄色大毛片视频在线观看| 老司机影院成人| 人人妻人人爽人人添夜夜欢视频| 哪个播放器可以免费观看大片| av卡一久久| 欧美一级a爱片免费观看看| 2022亚洲国产成人精品| 精品卡一卡二卡四卡免费| 男人添女人高潮全过程视频| 日本色播在线视频| 久久国内精品自在自线图片| 国产伦理片在线播放av一区| 国产一区二区在线观看av| av女优亚洲男人天堂| 欧美激情国产日韩精品一区| 各种免费的搞黄视频| 人妻系列 视频| 日韩熟女老妇一区二区性免费视频| 精品久久久久久电影网| 老司机亚洲免费影院| 欧美老熟妇乱子伦牲交| 国产一区二区在线观看av| 欧美精品亚洲一区二区| 久久久国产精品麻豆| 少妇的逼好多水| 欧美亚洲日本最大视频资源| 在线看a的网站| 黄色怎么调成土黄色| 国产精品麻豆人妻色哟哟久久| 亚洲精品自拍成人| 精品少妇内射三级| 成人午夜精彩视频在线观看| 精品国产乱码久久久久久小说| 人妻人人澡人人爽人人| 亚洲精华国产精华液的使用体验| 国产国语露脸激情在线看| 久久久久精品性色| 成人免费观看视频高清| 男女边吃奶边做爰视频| 蜜桃久久精品国产亚洲av| 国产免费一级a男人的天堂| 黑人欧美特级aaaaaa片| 久久久a久久爽久久v久久| av播播在线观看一区| 欧美 亚洲 国产 日韩一| 亚洲久久久国产精品| 爱豆传媒免费全集在线观看| 婷婷成人精品国产| 久久人妻熟女aⅴ| 国产成人av激情在线播放 | 3wmmmm亚洲av在线观看| 国产老妇伦熟女老妇高清| 国产探花极品一区二区| 一边亲一边摸免费视频| 99热这里只有精品一区| a级毛片在线看网站| 嫩草影院入口| 色5月婷婷丁香| 各种免费的搞黄视频| 国产一级毛片在线| 纵有疾风起免费观看全集完整版| 免费观看性生交大片5| 夜夜爽夜夜爽视频| 视频区图区小说| 狂野欧美激情性bbbbbb| 国产有黄有色有爽视频| av在线播放精品| 国产高清三级在线| 欧美日韩在线观看h| 热99国产精品久久久久久7| 久久精品国产亚洲av涩爱| 日韩中文字幕视频在线看片| 三上悠亚av全集在线观看| 亚洲精品乱码久久久久久按摩| 少妇猛男粗大的猛烈进出视频| 精品一区二区三卡| 亚洲国产av影院在线观看| 亚洲精品成人av观看孕妇| 大陆偷拍与自拍| 中国国产av一级| 丝袜喷水一区| 日韩伦理黄色片| 久久久久国产精品人妻一区二区| 久久精品人人爽人人爽视色| 丰满迷人的少妇在线观看| 国产精品麻豆人妻色哟哟久久| 只有这里有精品99| 老司机影院成人| 久久久国产精品麻豆| 亚洲内射少妇av| 精品卡一卡二卡四卡免费| 国产伦精品一区二区三区视频9| 亚洲国产欧美在线一区| 美女福利国产在线| 天天影视国产精品| 大香蕉久久成人网| 成人无遮挡网站| 国产午夜精品一二区理论片| 国产欧美日韩综合在线一区二区| 丰满少妇做爰视频| 久久这里有精品视频免费| 69精品国产乱码久久久| 国产精品久久久久成人av| 午夜福利,免费看| 99热网站在线观看| 精品国产国语对白av| 80岁老熟妇乱子伦牲交| 亚洲色图 男人天堂 中文字幕 | 日韩伦理黄色片| a级毛片黄视频| 欧美日本中文国产一区发布| 桃花免费在线播放| 亚洲四区av| 另类亚洲欧美激情| 国产极品粉嫩免费观看在线 | 精品亚洲乱码少妇综合久久| 久久久久久久久久久久大奶| av国产精品久久久久影院| 国产日韩欧美在线精品| 黄片无遮挡物在线观看| 午夜福利影视在线免费观看| 国产成人午夜福利电影在线观看| 日本色播在线视频| 一区二区三区乱码不卡18| 午夜视频国产福利| 日韩制服骚丝袜av| 久久久精品区二区三区| 久久人妻熟女aⅴ| 啦啦啦视频在线资源免费观看| 亚洲第一av免费看| 国产一区二区在线观看日韩| 中文乱码字字幕精品一区二区三区| 伊人久久精品亚洲午夜| 国产一区亚洲一区在线观看| 蜜桃国产av成人99| 日韩一本色道免费dvd| 一区二区三区免费毛片| 亚洲久久久国产精品| 男女边吃奶边做爰视频| av黄色大香蕉| 如日韩欧美国产精品一区二区三区 | 中文乱码字字幕精品一区二区三区| 亚洲少妇的诱惑av| 日产精品乱码卡一卡2卡三| 插阴视频在线观看视频| 99国产精品免费福利视频| 免费久久久久久久精品成人欧美视频 | 91久久精品国产一区二区三区| 亚洲精品乱码久久久v下载方式| 久久久久久久久久久久大奶| 亚洲国产最新在线播放| 国产精品国产av在线观看| 美女国产视频在线观看| 亚洲av综合色区一区| 国产永久视频网站| 午夜老司机福利剧场| 欧美激情国产日韩精品一区| 国产精品三级大全| 十八禁高潮呻吟视频| 黄片无遮挡物在线观看| 国产精品99久久99久久久不卡 | 欧美另类一区| 欧美成人午夜免费资源| 亚洲精品日韩在线中文字幕| 久久99热这里只频精品6学生| 丁香六月天网| 成人手机av| 最近手机中文字幕大全| 一本一本久久a久久精品综合妖精 国产伦在线观看视频一区 | 超碰97精品在线观看| 精品熟女少妇av免费看| 久久99热这里只频精品6学生| 91成人精品电影| 久久99精品国语久久久| 夜夜骑夜夜射夜夜干| 欧美激情国产日韩精品一区| 自拍欧美九色日韩亚洲蝌蚪91| 性色av一级| 久久99一区二区三区| 日日摸夜夜添夜夜爱| 成人午夜精彩视频在线观看| 国产一区二区在线观看av| 久久久久国产精品人妻一区二区| 欧美人与善性xxx| 一区二区av电影网| 在线免费观看不下载黄p国产| 午夜福利,免费看| 777米奇影视久久| 久久精品国产a三级三级三级| 亚洲精品自拍成人| 日韩欧美精品免费久久| 久久久精品免费免费高清| 久久99蜜桃精品久久| 中文字幕av电影在线播放| 亚洲av综合色区一区| 日韩熟女老妇一区二区性免费视频| 欧美日韩视频精品一区| 卡戴珊不雅视频在线播放| 极品少妇高潮喷水抽搐| 高清黄色对白视频在线免费看| xxx大片免费视频| 最后的刺客免费高清国语| 又黄又爽又刺激的免费视频.| 高清午夜精品一区二区三区| 97超视频在线观看视频| 国产精品熟女久久久久浪| h视频一区二区三区| 交换朋友夫妻互换小说| 91久久精品国产一区二区三区| 亚洲欧美日韩另类电影网站| 久久av网站| 亚洲欧美成人综合另类久久久| 精品少妇内射三级| 国产高清三级在线| 在线看a的网站| 国产视频首页在线观看| 亚洲精品一区蜜桃| 久久免费观看电影| 亚洲欧美日韩卡通动漫| 亚洲三级黄色毛片| 3wmmmm亚洲av在线观看| 美女国产高潮福利片在线看| 国产毛片在线视频| 97超视频在线观看视频| 麻豆精品久久久久久蜜桃| 亚洲中文av在线| 精品午夜福利在线看| 国产精品一二三区在线看| 大香蕉久久网| 欧美精品国产亚洲| 丁香六月天网| 免费播放大片免费观看视频在线观看| 女性被躁到高潮视频| 日韩,欧美,国产一区二区三区| 人人妻人人澡人人看| 建设人人有责人人尽责人人享有的| 色94色欧美一区二区| 亚洲国产精品一区三区| 一级爰片在线观看| 制服诱惑二区| 在线 av 中文字幕| 亚洲av二区三区四区| 国产精品一区二区在线不卡| 久久国产精品男人的天堂亚洲 | 久久国产精品大桥未久av| 成人18禁高潮啪啪吃奶动态图 | 日韩,欧美,国产一区二区三区| 人人妻人人澡人人看| 色婷婷久久久亚洲欧美| 少妇的逼水好多| 国产成人免费观看mmmm| 亚洲成人一二三区av| 欧美精品人与动牲交sv欧美| 亚洲精品国产av成人精品| 国产成人精品福利久久| 亚洲国产精品一区三区| 大又大粗又爽又黄少妇毛片口| 最黄视频免费看| 一边摸一边做爽爽视频免费| 天堂中文最新版在线下载| 亚洲成人手机| 日韩,欧美,国产一区二区三区| 午夜久久久在线观看| 少妇精品久久久久久久| 亚洲国产av新网站| 天天影视国产精品| 国产亚洲最大av| 午夜福利,免费看| 日韩一本色道免费dvd| 一级黄片播放器| 一区在线观看完整版| av线在线观看网站| 国产亚洲欧美精品永久| 亚洲性久久影院| 国产免费视频播放在线视频| 欧美xxⅹ黑人| 人妻系列 视频| 久久精品国产自在天天线| 精品人妻一区二区三区麻豆| 亚洲欧洲日产国产| 国产在线视频一区二区| 国产在视频线精品| 亚洲av国产av综合av卡| 18禁观看日本| 欧美精品国产亚洲| 伦精品一区二区三区| 国产成人a∨麻豆精品| 欧美日韩一区二区视频在线观看视频在线| 在线 av 中文字幕| 在线观看免费日韩欧美大片 | 中文字幕av电影在线播放| 国产成人精品在线电影| 亚洲国产精品成人久久小说| 中文欧美无线码| 丝瓜视频免费看黄片| 亚洲精品色激情综合| 日韩欧美一区视频在线观看| 最近中文字幕高清免费大全6| 伦理电影大哥的女人| 免费观看的影片在线观看| 亚洲激情五月婷婷啪啪| 国产爽快片一区二区三区| 亚洲精品色激情综合| 久久久欧美国产精品| av在线播放精品| 高清午夜精品一区二区三区| 99热全是精品| 亚洲国产av新网站| 午夜免费观看性视频| 日本爱情动作片www.在线观看| 国产白丝娇喘喷水9色精品| 我要看黄色一级片免费的| 午夜日本视频在线| 成年av动漫网址| 久久久久视频综合| 国产伦理片在线播放av一区| 成人无遮挡网站| 国产精品 国内视频| 久热久热在线精品观看| 亚洲精品色激情综合| 亚洲精品一二三| av卡一久久| 边亲边吃奶的免费视频| 香蕉精品网在线| 18禁裸乳无遮挡动漫免费视频| 久久精品熟女亚洲av麻豆精品| 久久ye,这里只有精品| 女性生殖器流出的白浆| 亚洲av免费高清在线观看| 免费高清在线观看视频在线观看| 18禁观看日本| 有码 亚洲区| 欧美+日韩+精品| 欧美日韩国产mv在线观看视频| 我的老师免费观看完整版| 欧美3d第一页| 我的女老师完整版在线观看| 蜜桃在线观看..| 少妇猛男粗大的猛烈进出视频| 夜夜看夜夜爽夜夜摸| 国产精品国产三级国产专区5o| 最黄视频免费看| 欧美97在线视频| 久久久久久久精品精品| 寂寞人妻少妇视频99o| 久久久精品区二区三区| 成人国语在线视频| 91精品一卡2卡3卡4卡| 久久午夜福利片| 国产永久视频网站| 高清毛片免费看| 国产精品欧美亚洲77777| av有码第一页| 国产成人aa在线观看| 精品国产乱码久久久久久小说| 国产色爽女视频免费观看| 狂野欧美激情性xxxx在线观看| 九草在线视频观看| 简卡轻食公司| 日韩精品有码人妻一区| 秋霞伦理黄片| 卡戴珊不雅视频在线播放| 美女大奶头黄色视频| 中文字幕亚洲精品专区| 2018国产大陆天天弄谢| 伦理电影大哥的女人| 国产精品一区www在线观看| 久久精品久久久久久噜噜老黄| 在线观看一区二区三区激情| 在线观看美女被高潮喷水网站| 免费高清在线观看日韩| 麻豆乱淫一区二区| 最近最新中文字幕免费大全7| 亚洲人成网站在线观看播放| 国模一区二区三区四区视频| 日产精品乱码卡一卡2卡三| 乱码一卡2卡4卡精品| 美女xxoo啪啪120秒动态图| 人体艺术视频欧美日本| 韩国高清视频一区二区三区| 国产欧美另类精品又又久久亚洲欧美| 日日摸夜夜添夜夜添av毛片| 一本久久精品| 成人黄色视频免费在线看| 国产色婷婷99| 亚洲色图综合在线观看| 精品人妻偷拍中文字幕| 乱人伦中国视频| 制服人妻中文乱码| 免费观看av网站的网址| 飞空精品影院首页| 国产男女内射视频| 国产成人精品在线电影| a级毛色黄片| 黑人高潮一二区| 看免费成人av毛片| 最近2019中文字幕mv第一页| 极品人妻少妇av视频| 亚洲无线观看免费| 五月天丁香电影| 99热国产这里只有精品6| 在线观看免费视频网站a站| 丝袜美足系列| 汤姆久久久久久久影院中文字幕| 99久久人妻综合| 中文字幕最新亚洲高清| 狠狠婷婷综合久久久久久88av| 99热这里只有是精品在线观看| 色婷婷av一区二区三区视频| 久久亚洲国产成人精品v| av不卡在线播放| 精品国产乱码久久久久久小说| 日韩视频在线欧美| 人人澡人人妻人| 日本黄大片高清| 狠狠精品人妻久久久久久综合| 国产亚洲欧美精品永久| 人人妻人人澡人人爽人人夜夜| 美女国产高潮福利片在线看| 国产精品秋霞免费鲁丝片| 国产不卡av网站在线观看| 久久国内精品自在自线图片| 婷婷色综合大香蕉| 少妇的逼水好多| 春色校园在线视频观看| 欧美人与性动交α欧美精品济南到 | 国产亚洲午夜精品一区二区久久| 少妇高潮的动态图| 少妇熟女欧美另类| 黑人高潮一二区| 成人免费观看视频高清| 国产视频内射| 在线观看三级黄色| 久久精品人人爽人人爽视色| 久久久精品区二区三区| 日韩av免费高清视频| 久久精品国产亚洲av涩爱| 一本久久精品| 亚洲精品视频女| 黑人巨大精品欧美一区二区蜜桃 | 校园人妻丝袜中文字幕| 国产高清不卡午夜福利| 黄色视频在线播放观看不卡| 婷婷色综合大香蕉| 欧美97在线视频| xxxhd国产人妻xxx|