• <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ū)
    “青春期”湍流中的智慧引渡(三)
    国产欧美亚洲国产| 男女下面插进去视频免费观看| 久久久欧美国产精品| 狠狠精品人妻久久久久久综合| 一区二区三区乱码不卡18| 久久久久国产一级毛片高清牌| 一级,二级,三级黄色视频| 在线天堂中文资源库| 天天影视国产精品| 妹子高潮喷水视频| 亚洲国产最新在线播放| 中文字幕最新亚洲高清| 菩萨蛮人人尽说江南好唐韦庄| 欧美日韩福利视频一区二区| 欧美日韩综合久久久久久| 一区福利在线观看| 亚洲精品国产av蜜桃| 久久久久网色| 啦啦啦中文免费视频观看日本| 国产成人免费观看mmmm| 最新的欧美精品一区二区| 国产精品一区二区在线观看99| 亚洲国产看品久久| 男人舔女人的私密视频| 黄色视频在线播放观看不卡| 一区二区三区乱码不卡18| 日韩 亚洲 欧美在线| 国产精品一区二区精品视频观看| www日本在线高清视频| 国产黄色视频一区二区在线观看| 亚洲精品日本国产第一区| 夜夜骑夜夜射夜夜干| 一区二区三区四区激情视频| 国产视频首页在线观看| 99久久精品国产亚洲精品| 免费高清在线观看视频在线观看| 午夜福利在线免费观看网站| 亚洲成人av在线免费| 免费人妻精品一区二区三区视频| 久久久久精品久久久久真实原创| 久久久久人妻精品一区果冻| 亚洲婷婷狠狠爱综合网| 女性被躁到高潮视频| 日韩一区二区三区影片| 亚洲成人免费av在线播放| 欧美在线黄色| 国产成人午夜福利电影在线观看| 90打野战视频偷拍视频| 久久精品国产亚洲av高清一级| 不卡视频在线观看欧美| 亚洲人成网站在线观看播放| 日本欧美视频一区| 亚洲精华国产精华液的使用体验| 久久精品国产亚洲av涩爱| 人人妻人人澡人人爽人人夜夜| 在线观看人妻少妇| 热99国产精品久久久久久7| 91精品伊人久久大香线蕉| 久久久久人妻精品一区果冻| 久久久久国产一级毛片高清牌| 国产国语露脸激情在线看| 十八禁人妻一区二区| 成年美女黄网站色视频大全免费| 黄色视频不卡| 国产成人午夜福利电影在线观看| 国产精品久久久av美女十八| 丝瓜视频免费看黄片| 亚洲精品成人av观看孕妇| 日本欧美国产在线视频| 午夜日韩欧美国产| 久久青草综合色| 中国国产av一级| 大片电影免费在线观看免费| 黑人猛操日本美女一级片| 人人澡人人妻人| 色综合欧美亚洲国产小说| 超色免费av| 老司机影院毛片| 久久久欧美国产精品| 在线观看人妻少妇| 黄色视频在线播放观看不卡| 天美传媒精品一区二区| 国产一区二区在线观看av| 欧美在线一区亚洲| 亚洲欧美精品自产自拍| 伦理电影大哥的女人| 美女午夜性视频免费| 97人妻天天添夜夜摸| 中国国产av一级| 操出白浆在线播放| 你懂的网址亚洲精品在线观看| 国产精品香港三级国产av潘金莲 | 国产成人啪精品午夜网站| 91国产中文字幕| 男女下面插进去视频免费观看| 久久人妻熟女aⅴ| 免费高清在线观看日韩| 女人久久www免费人成看片| 十八禁网站网址无遮挡| 亚洲第一青青草原| 另类亚洲欧美激情| 男女无遮挡免费网站观看| 精品国产一区二区三区四区第35| 丁香六月天网| 一级毛片 在线播放| 国产欧美日韩综合在线一区二区| 如日韩欧美国产精品一区二区三区| 这个男人来自地球电影免费观看 | 久久精品亚洲熟妇少妇任你| 中文字幕精品免费在线观看视频| 又黄又粗又硬又大视频| 中文字幕精品免费在线观看视频| 亚洲综合精品二区| 美国免费a级毛片| 亚洲国产最新在线播放| 国产在视频线精品| 日韩欧美精品免费久久| 欧美乱码精品一区二区三区| 欧美乱码精品一区二区三区| 无限看片的www在线观看| 免费在线观看视频国产中文字幕亚洲 | 一区二区三区四区激情视频| 黄色怎么调成土黄色| 啦啦啦中文免费视频观看日本| av国产精品久久久久影院| 亚洲情色 制服丝袜| 精品第一国产精品| 成人漫画全彩无遮挡| 观看美女的网站| 国产成人av激情在线播放| 天美传媒精品一区二区| 日韩成人av中文字幕在线观看| 亚洲欧洲精品一区二区精品久久久 | 天天添夜夜摸| 亚洲久久久国产精品| 亚洲人成网站在线观看播放| 一级爰片在线观看| 自线自在国产av| 激情视频va一区二区三区| 亚洲精品日本国产第一区| 国产精品一区二区在线不卡| 亚洲国产av影院在线观看| 香蕉国产在线看| 90打野战视频偷拍视频| 91精品三级在线观看| 久久毛片免费看一区二区三区| 国产精品麻豆人妻色哟哟久久| 纯流量卡能插随身wifi吗| 国产乱来视频区| 老司机深夜福利视频在线观看 | 久久精品aⅴ一区二区三区四区| 97在线人人人人妻| 亚洲专区中文字幕在线 | 国产又爽黄色视频| 伊人久久大香线蕉亚洲五| 丰满饥渴人妻一区二区三| 午夜激情av网站| av女优亚洲男人天堂| 午夜福利影视在线免费观看| 高清不卡的av网站| 国产一区二区 视频在线| 免费av中文字幕在线| 少妇的丰满在线观看| av在线观看视频网站免费| 免费观看a级毛片全部| 超色免费av| 亚洲精品久久久久久婷婷小说| 天天躁夜夜躁狠狠久久av| 一级毛片黄色毛片免费观看视频| 九草在线视频观看| 高清视频免费观看一区二区| 国产不卡av网站在线观看| h视频一区二区三区| av在线老鸭窝| 在线观看人妻少妇| 精品亚洲成a人片在线观看| 99久久人妻综合| 黑人巨大精品欧美一区二区蜜桃| 国产无遮挡羞羞视频在线观看| 久久久久精品性色| 欧美国产精品va在线观看不卡| 欧美国产精品va在线观看不卡| 丰满少妇做爰视频| 99九九在线精品视频| 亚洲精品在线美女| 2021少妇久久久久久久久久久| 精品卡一卡二卡四卡免费| 国产成人免费无遮挡视频| 久久青草综合色| 男女边摸边吃奶| 亚洲第一区二区三区不卡| 亚洲少妇的诱惑av| 国产又爽黄色视频| h视频一区二区三区| 欧美97在线视频| 丰满饥渴人妻一区二区三| 一区二区三区激情视频| 亚洲欧美精品综合一区二区三区| 中文字幕精品免费在线观看视频| 亚洲精品久久午夜乱码| av网站在线播放免费| 操美女的视频在线观看| 亚洲欧美精品自产自拍| 国产精品国产三级国产专区5o| 亚洲精品第二区| 久热爱精品视频在线9| 一级a爱视频在线免费观看| 狂野欧美激情性bbbbbb| 免费久久久久久久精品成人欧美视频| 色综合欧美亚洲国产小说| 国产亚洲午夜精品一区二区久久| 国产精品欧美亚洲77777| 久久久精品94久久精品| 一区二区三区乱码不卡18| 亚洲精品日本国产第一区| 一二三四在线观看免费中文在| 精品国产超薄肉色丝袜足j| 九色亚洲精品在线播放| av不卡在线播放| 国产精品女同一区二区软件| 中文精品一卡2卡3卡4更新| 日韩不卡一区二区三区视频在线| 日韩av在线免费看完整版不卡| 五月开心婷婷网| 国产伦理片在线播放av一区| 美女午夜性视频免费| 97在线人人人人妻| 激情视频va一区二区三区| 精品国产超薄肉色丝袜足j| 国产成人欧美| 嫩草影视91久久| av在线老鸭窝| 日本欧美视频一区| 嫩草影视91久久| 97在线人人人人妻| 欧美97在线视频| 国产不卡av网站在线观看| 在线精品无人区一区二区三| 国产毛片在线视频| 涩涩av久久男人的天堂| 日日啪夜夜爽| 97人妻天天添夜夜摸| 超碰97精品在线观看| 性色av一级| 下体分泌物呈黄色| 亚洲精品日本国产第一区| 国产 精品1| 国产精品一二三区在线看| 晚上一个人看的免费电影| 丁香六月天网| 国产精品香港三级国产av潘金莲 | 国产精品嫩草影院av在线观看| 国产精品熟女久久久久浪| 亚洲欧美一区二区三区国产| 亚洲熟女毛片儿| 飞空精品影院首页| 亚洲少妇的诱惑av| 欧美 日韩 精品 国产| 欧美成人精品欧美一级黄| 久久久国产欧美日韩av| 美女午夜性视频免费| 午夜免费男女啪啪视频观看| 久久久久精品人妻al黑| 人体艺术视频欧美日本| 欧美在线一区亚洲| 中文欧美无线码| 国产一区二区三区av在线| 午夜激情久久久久久久| 久久热在线av| 丝袜美足系列| 日本欧美国产在线视频| 亚洲精品成人av观看孕妇| 欧美亚洲日本最大视频资源| 国产免费又黄又爽又色| 亚洲,欧美精品.| 成人免费观看视频高清| 亚洲国产精品一区二区三区在线| 欧美精品一区二区大全| 久久精品亚洲av国产电影网| 国产麻豆69| 国产成人啪精品午夜网站| 丰满少妇做爰视频| 在线免费观看不下载黄p国产| 久久久久久久大尺度免费视频| 在线观看人妻少妇| 久久天堂一区二区三区四区| 19禁男女啪啪无遮挡网站| 国产av码专区亚洲av| 亚洲精品自拍成人| 国产精品一国产av| 亚洲一级一片aⅴ在线观看| 一级毛片电影观看| av国产精品久久久久影院| 不卡av一区二区三区| 一区二区三区激情视频| 97在线人人人人妻| 午夜福利一区二区在线看| 天堂中文最新版在线下载| 久久性视频一级片| 成人国语在线视频| 久久热在线av| 最近手机中文字幕大全| 日韩中文字幕欧美一区二区 | 曰老女人黄片| 悠悠久久av| 国产成人欧美| 尾随美女入室| 中文字幕另类日韩欧美亚洲嫩草| 丝袜人妻中文字幕| 国产麻豆69| 精品一区在线观看国产| 午夜福利视频精品| 久久久久精品性色| 亚洲欧美色中文字幕在线| 午夜福利,免费看| 中文字幕另类日韩欧美亚洲嫩草| 亚洲精品久久久久久婷婷小说| 日韩一卡2卡3卡4卡2021年| 欧美少妇被猛烈插入视频| 黄色怎么调成土黄色| 久久人人爽人人片av| 男女午夜视频在线观看| 国产午夜精品一二区理论片| 90打野战视频偷拍视频| 亚洲国产日韩一区二区| 一级毛片黄色毛片免费观看视频| 亚洲精品,欧美精品| 欧美国产精品一级二级三级| 午夜日韩欧美国产| 久久人人爽人人片av| 亚洲成人免费av在线播放| 人人妻人人爽人人添夜夜欢视频| 色网站视频免费| 精品国产一区二区三区久久久樱花| 最黄视频免费看| 丰满乱子伦码专区| 男女国产视频网站| 男男h啪啪无遮挡| 99久久99久久久精品蜜桃| 水蜜桃什么品种好| 午夜91福利影院| 久久久精品区二区三区| 日韩 欧美 亚洲 中文字幕| 国产一区二区激情短视频 | 爱豆传媒免费全集在线观看| 天美传媒精品一区二区| 亚洲精品日本国产第一区| 一级黄片播放器| 涩涩av久久男人的天堂| h视频一区二区三区| 日韩欧美一区视频在线观看| 黑丝袜美女国产一区| 欧美精品人与动牲交sv欧美| 欧美亚洲 丝袜 人妻 在线| 国产极品粉嫩免费观看在线| 国产97色在线日韩免费| 亚洲欧洲精品一区二区精品久久久 | videosex国产| 亚洲精品日韩在线中文字幕| 亚洲国产av影院在线观看| 久久青草综合色| 丝袜脚勾引网站| 丝袜喷水一区| 免费观看人在逋| 日日撸夜夜添| 搡老乐熟女国产| 91精品伊人久久大香线蕉| 蜜桃国产av成人99| 亚洲欧美色中文字幕在线| 久久精品国产亚洲av涩爱| 男女下面插进去视频免费观看| 少妇被粗大猛烈的视频| www.熟女人妻精品国产| av卡一久久| 少妇 在线观看| 纯流量卡能插随身wifi吗| 久久久久久人妻| 中文字幕另类日韩欧美亚洲嫩草| 18在线观看网站| 最近最新中文字幕免费大全7| 一二三四中文在线观看免费高清| 少妇精品久久久久久久| 老汉色av国产亚洲站长工具| 永久免费av网站大全| 视频区图区小说| 中文字幕最新亚洲高清| 欧美国产精品va在线观看不卡| 一边摸一边做爽爽视频免费| 日韩精品免费视频一区二区三区| 国产高清不卡午夜福利| 国产精品偷伦视频观看了| 久久久久国产一级毛片高清牌| 国产av国产精品国产| 精品少妇一区二区三区视频日本电影 | 高清av免费在线| 大片免费播放器 马上看| 97精品久久久久久久久久精品| 国产日韩欧美亚洲二区| 欧美黄色片欧美黄色片| 亚洲精品,欧美精品| 大片电影免费在线观看免费| 久久精品国产综合久久久| 亚洲精品中文字幕在线视频| 在线观看一区二区三区激情| 亚洲国产日韩一区二区| 亚洲av中文av极速乱| 一本色道久久久久久精品综合| 久久午夜综合久久蜜桃| 97精品久久久久久久久久精品| 精品亚洲乱码少妇综合久久| 亚洲第一青青草原| 丝袜美腿诱惑在线| 亚洲一卡2卡3卡4卡5卡精品中文| 自拍欧美九色日韩亚洲蝌蚪91| 91精品国产国语对白视频| 亚洲三区欧美一区| 久久久久久人妻| 久久天堂一区二区三区四区| 你懂的网址亚洲精品在线观看| 亚洲欧洲国产日韩| 人人澡人人妻人| 久久久久国产精品人妻一区二区| 国产精品久久久久久久久免| 国产一区有黄有色的免费视频| 捣出白浆h1v1| 电影成人av| 各种免费的搞黄视频| 亚洲 欧美一区二区三区| 99热全是精品| 精品视频人人做人人爽| 捣出白浆h1v1| 国产黄色免费在线视频| 国产成人啪精品午夜网站| 国产高清国产精品国产三级| 老司机亚洲免费影院| 高清在线视频一区二区三区| 亚洲成av片中文字幕在线观看| 久久人人97超碰香蕉20202| 亚洲国产看品久久| 狠狠婷婷综合久久久久久88av| 久久这里只有精品19| 日日爽夜夜爽网站| 亚洲成人一二三区av| 蜜桃在线观看..| 日韩av不卡免费在线播放| 中文欧美无线码| 精品少妇久久久久久888优播| 丝袜脚勾引网站| 欧美在线黄色| 国产av码专区亚洲av| 日韩不卡一区二区三区视频在线| 久久久久久久久免费视频了| 悠悠久久av| 久久女婷五月综合色啪小说| 日韩一卡2卡3卡4卡2021年| 最近手机中文字幕大全| 99久久人妻综合| 99国产精品免费福利视频| 亚洲av电影在线进入| 国产麻豆69| 啦啦啦视频在线资源免费观看| a级毛片黄视频| 久久久国产欧美日韩av| 九草在线视频观看| 国产 一区精品| 欧美久久黑人一区二区| 最近2019中文字幕mv第一页| 亚洲,欧美精品.| 久久99精品国语久久久| 91老司机精品| 另类精品久久| 高清视频免费观看一区二区| 亚洲七黄色美女视频| 精品国产乱码久久久久久小说| 99热网站在线观看| 一区二区三区乱码不卡18| 成年av动漫网址| 欧美在线一区亚洲| 伦理电影免费视频| 秋霞在线观看毛片| 亚洲欧美成人精品一区二区| 大香蕉久久成人网| 日日摸夜夜添夜夜爱| 久久久久久人人人人人| 精品酒店卫生间| 街头女战士在线观看网站| 久久青草综合色| 亚洲精品在线美女| 黑丝袜美女国产一区| 搡老乐熟女国产| 秋霞在线观看毛片| 看十八女毛片水多多多| 国产成人免费观看mmmm| 欧美日韩国产mv在线观看视频| 午夜免费观看性视频| 日韩人妻精品一区2区三区| 日本猛色少妇xxxxx猛交久久| 黄色 视频免费看| 亚洲欧美精品自产自拍| 性少妇av在线| 国产免费一区二区三区四区乱码| 97在线人人人人妻| 免费女性裸体啪啪无遮挡网站| 亚洲一级一片aⅴ在线观看| 两个人看的免费小视频| 日日啪夜夜爽| 亚洲欧美一区二区三区黑人| 爱豆传媒免费全集在线观看| 五月天丁香电影| 热99久久久久精品小说推荐| 免费女性裸体啪啪无遮挡网站| 七月丁香在线播放| 最近手机中文字幕大全| 久久狼人影院| 18在线观看网站| 黑人欧美特级aaaaaa片| 国产一级毛片在线| 制服丝袜香蕉在线| 国产精品免费大片| 国产精品一二三区在线看| 亚洲,欧美,日韩| 国产熟女午夜一区二区三区| 91老司机精品| 亚洲图色成人| 91aial.com中文字幕在线观看| 久久久久国产精品人妻一区二区| 永久免费av网站大全| 一本久久精品| 捣出白浆h1v1| 人人妻人人添人人爽欧美一区卜| av在线老鸭窝| 精品国产一区二区久久| www.自偷自拍.com| 欧美日韩视频精品一区| 亚洲情色 制服丝袜| 国产精品久久久人人做人人爽| 高清av免费在线| 丝瓜视频免费看黄片| 哪个播放器可以免费观看大片| 女人爽到高潮嗷嗷叫在线视频| 精品一品国产午夜福利视频| 国产国语露脸激情在线看| 2021少妇久久久久久久久久久| 欧美激情极品国产一区二区三区| 一级毛片电影观看| 亚洲精品一区蜜桃| 天天躁夜夜躁狠狠久久av| 一二三四中文在线观看免费高清| 日韩一区二区视频免费看| 桃花免费在线播放| 久久人人爽人人片av| 丝袜美腿诱惑在线| 91精品国产国语对白视频| 成年动漫av网址| 国产片内射在线| 国产淫语在线视频| 午夜av观看不卡| 亚洲自偷自拍图片 自拍| 看免费成人av毛片| 少妇人妻久久综合中文| 999久久久国产精品视频| av网站免费在线观看视频| 欧美日韩国产mv在线观看视频| 欧美在线黄色| av又黄又爽大尺度在线免费看| 欧美日韩亚洲高清精品| 黑丝袜美女国产一区| 日本猛色少妇xxxxx猛交久久| av在线app专区| 大香蕉久久成人网| 超碰97精品在线观看| 18在线观看网站| 日本欧美国产在线视频| 纵有疾风起免费观看全集完整版| 亚洲国产最新在线播放| 亚洲,欧美精品.| 日本猛色少妇xxxxx猛交久久| 国产欧美日韩一区二区三区在线| 又大又爽又粗| a 毛片基地| 国产精品 国内视频| 久久久精品区二区三区| 十八禁人妻一区二区| 丰满饥渴人妻一区二区三| 亚洲国产精品一区二区三区在线| 国产在视频线精品| 欧美精品一区二区大全| 不卡av一区二区三区| 精品一区二区三卡| 一本—道久久a久久精品蜜桃钙片| 亚洲精品在线美女| 国产成人精品无人区| 久久人人爽av亚洲精品天堂| 精品亚洲成a人片在线观看| 久久精品亚洲av国产电影网| 我的亚洲天堂| 人人妻人人澡人人爽人人夜夜| 国产精品女同一区二区软件| 极品少妇高潮喷水抽搐| 精品国产一区二区久久| 麻豆精品久久久久久蜜桃| 亚洲欧美精品自产自拍| 精品一区二区三区av网在线观看 | 国产野战对白在线观看| 国产精品秋霞免费鲁丝片| 麻豆精品久久久久久蜜桃| 亚洲欧美精品自产自拍| 国产一区二区 视频在线| 男女床上黄色一级片免费看| 亚洲欧洲精品一区二区精品久久久 | 看免费av毛片| 男人添女人高潮全过程视频| 操出白浆在线播放| 久久人人爽人人片av| 欧美日韩成人在线一区二区|