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

    一種基于聚類分析的二維激波模式識別算法

    2020-09-10 03:25:34常思源白曉征劉君
    航空學(xué)報 2020年8期
    關(guān)鍵詞:波點激波壁面

    常思源,白曉征,劉君

    大連理工大學(xué) 航空航天學(xué)院,大連 116024

    激波的探測與可視化是計算流體力學(xué)(CFD)領(lǐng)域中一個非常富有挑戰(zhàn)的問題。在某些流場參數(shù)的等值線云圖中,通常認(rèn)為激波間斷位于大量等值線聚集處,但對于存在接觸間斷、膨脹波和旋渦等復(fù)雜波系干擾的流場,該方法很容易對激波產(chǎn)生誤判,因此有必要發(fā)展更加精準(zhǔn)的激波探測(Shock Wave Detection)技術(shù)。

    至今,很多學(xué)者[1-14]基于激波捕捉(Shock-Capturing)法計算出的流場提出了各種各樣激波探測的方法。1985年,Buning和Steger[1]首次提出將沿激波法向(近似用壓強梯度方向代替)馬赫數(shù)等于1的等值面視為激波面;該方法隨后被Darmofal[2]稱為基于正則馬赫數(shù)(Normal Mach Number)的激波探測法,并逐漸在流場可視化中被廣泛應(yīng)用;Liou等[3]在此基礎(chǔ)上介紹了3種濾波技術(shù),以剔除辨識出的偽激波;隨后,Lovely和Haimes[4]又成功將其推廣到非定常流動中用以探測運動激波。1993年,Pagendarm和Seitz[5]詳細(xì)地介紹了一種基于流場密度方向?qū)?shù)(Directional Derivative)來辨識激波的方法,即將密度二階方向?qū)?shù)為0的等值面視為激波陣面,并結(jié)合密度的一階方向?qū)?shù)(密度梯度在當(dāng)?shù)厮俣仁噶可系耐队?來過濾噪點;其后,Ma等[6]的數(shù)值算例表明基于正則馬赫數(shù)過濾噪點可能會獲得更好的結(jié)果。1994年,Van Rosendale[7]提出了一種針對二維非結(jié)構(gòu)網(wǎng)格的激波擬合算法,通過比較網(wǎng)格節(jié)點之間的密度梯度來移動網(wǎng)格使其邊緣與激波線對齊;Ma等[6]認(rèn)為該算法中局部加權(quán)密度梯度的策略可以用于識別出激波附近的網(wǎng)格節(jié)點,但可能需要結(jié)合其他流場特征來實現(xiàn)更準(zhǔn)確的辨識。2011年,Kanamori和Suzuki[8]巧妙地將特征線理論(Theory of Characteristics)運用于二維激波探測中,并成功推廣到非定常[8]和三維流場[9]計算中,相比前述方法,該算法雖然較為復(fù)雜,但不需要人為設(shè)置一些經(jīng)驗性的閾值參數(shù)進行濾波就可以獲得更加良好的結(jié)果。2017年,Akhlaghi等[10]通過分析數(shù)值紋影圖內(nèi)某些流動參數(shù)的高斯分布(Gaussian Distribution),提供了一種新的激波辨識方法,對連續(xù)流和稀薄流都比較適用。近年來還有學(xué)者[11-12]嘗試將深度學(xué)習(xí)(Deep Learning)融入大規(guī)模流場的激波識別中,取得了不錯的結(jié)果。

    國內(nèi)針對激波探測的研究報道相對較少,馬千里等[13]基于壓強梯度計算正則馬赫數(shù),并利用激波物理特征,結(jié)合光線投射法提出了一種基于兩級采樣的多激波提取方法。2013年,Wu等[14]在其綜述中總結(jié)了以往激波辨識方法,并指出了對運動激波和弱激波的精確識別是非常有挑戰(zhàn)的研究方向。

    整體來看,以上這些激波探測法具有一個共同的特征,即它們都是基于“當(dāng)?shù)貥?biāo)準(zhǔn)(Local Criteria)”進行辨識的;也就是說,激波位置的識別僅通過分析局部一個網(wǎng)格節(jié)點/單元(最多考慮一層相鄰網(wǎng)格)上的流動參數(shù)來實現(xiàn)。Paciorri和Bonfiglioli[15]認(rèn)為這些探測技術(shù)的當(dāng)?shù)匦再|(zhì)正是其局限性和缺陷的根源。一方面,這些方法往往將探測到的一系列散點[6]或小線段[8-9]視為激波線/面,但由于數(shù)值耗散和振蕩,這些激波線/面處經(jīng)常會出現(xiàn)噪聲、偽激波分支或空隙[14]。另一方面,所有這些算法都不能判斷出流場中是否存在激波干擾點,即不能清晰地給出流場中激波干擾的模式,如λ激波、異側(cè)激波相交、激波-壁面反射等。

    近年來,有一批學(xué)者[16-19]重新開啟了激波裝配(Shock-Fitting)法的研究與應(yīng)用,相比當(dāng)前主流的激波捕捉法,該方法必須顯式地處理激波,因此在開始計算前,通常都需要從激波捕捉流場獲取相關(guān)信息,如合理的初始流動參數(shù)、激波的近似位置、激波的干擾結(jié)構(gòu)等。然而,由于上述激波探測法的缺陷,其很難直接用于激波裝配的初始化中[15],最終造成這些激波裝配算法難免需要依據(jù)先驗知識進行人工干預(yù)來設(shè)置初始激波及激波干擾點的位置。

    針對上述問題,Paciorri和Bonfiglioli[15]于2012年介紹了一種二維激波模式識別技術(shù),首先根據(jù)傳統(tǒng)基于當(dāng)?shù)孛芏榷A方向?qū)?shù)的方法[6]提取出激波點云,隨后利用霍夫變換(Hough Transform)搜尋這些散點中的特征直線并去除噪點,再采用最小二乘(Least-Squares)法擬合得到各條激波曲線,最后通過一套模糊邏輯(Fuzzy Logic)算法識別出激波干擾的模式。然而,該方法實施起來比較繁瑣,對于流場中長度較短的激波分支容易產(chǎn)生誤判,且很難適用于非定常流場。

    本文提出了一種面向全局(Global)的二維激波模式識別算法,對探測出的激波帶網(wǎng)格單元進一步處理,通過聚類分析(Cluster Analysis)方法將激波單元劃分成許多簇;分析各個簇的相鄰關(guān)系就可以清晰地判斷出激波干擾的模式,最終基于Bézier曲線擬合算法優(yōu)化的激波線和數(shù)值捕捉的激波相比在形狀和位置上都較為吻合。該算法不僅邏輯簡單,且不受網(wǎng)格類型的限制,同時在非定常流動中也具有較高的可靠性。

    1 二維激波模式識別算法

    結(jié)合一個定常無黏激波反射算例(幾何模型和流場如圖1所示,其中p表示壓強,下標(biāo)∞表示來流狀態(tài)),從3個主要方面介紹了該二維激波模式識別算法的步驟和特點。在該算例中,來流馬赫數(shù)Ma∞=1.5,楔角θ為10°,入口和出口高度分別為2.17L和2L。入射斜激波S1在上壁面發(fā)生馬赫反射產(chǎn)生了強度較強的馬赫干擾S2,而反射激波S3經(jīng)膨脹波異側(cè)干擾后強度下降,最后在下壁面產(chǎn)生了正規(guī)反射,次反射激波S4與出口邊界相交。此外,雖然該流場是基于均勻的非結(jié)構(gòu)三角形網(wǎng)格(網(wǎng)格尺寸Δ=0.015L)計算的,但是該識別算法同樣適用于結(jié)構(gòu)網(wǎng)格、混合網(wǎng)格和非均勻網(wǎng)格等,在2.2節(jié)有相關(guān)描述。

    1.1 辨識激波單元

    本文算法是針對網(wǎng)格單元而非網(wǎng)格散點實現(xiàn)的,因此首先需要從圖1所示流場提取出表征激波位置的網(wǎng)格單元,即圖2所示激波單元。考慮到引言所述若干激波探測方法的適用性和魯棒性,本文基于流場單元格心的壓強梯度,設(shè)計了一套針對激波的辨識準(zhǔn)則,具體步驟為

    圖2 定常激波反射:初始辨識的激波單元

    (1)

    由于考慮了流向,δp的正、負(fù)號分別對應(yīng)流體處于壓縮和膨脹的狀態(tài),而激波處通常表現(xiàn)出很強的壓縮性,因此首先要過濾掉δp<0的單元。

    2) 由于數(shù)值計算難免會存在耗散和誤差,會造成一些非激波區(qū)域的δp值為一個較小的正量,故過濾掉滿足式(2)的單元:

    δp<ξ1δpmax

    (2)

    式中:δpmax為流向壓強梯度的全場最大值;ξ1為全局濾波系數(shù),當(dāng)流場中激波的強度比較懸殊時,該系數(shù)過大可能會濾掉某些弱激波而導(dǎo)致部分激波帶缺失,因此經(jīng)大量數(shù)值試驗的調(diào)試,取ξ1=0.01。流場中除了激波間斷外可能還存在接觸間斷,由于數(shù)值捕捉出的接觸間斷兩側(cè)壓強近似相等,因此基于式(2)流向壓強梯度閾值法可以有效地過濾掉接觸間斷。

    3) 然而幾乎不可能只通過設(shè)置滿意的ξ1,既能完整地辨識多激波又能去除全部噪聲。此外,因全局濾波強度不大,經(jīng)第2步辨識的激波帶可能會出現(xiàn)一些偽激波單元,為了在一定程度上提高辨識方法的普適性,考慮局部特性對滿足式(3)的單元進行二次去噪:

    (3)

    定常激波反射算例經(jīng)上述方法辨識出的具有一定厚度的“激波帶”如圖3所示,受激波強弱影響,激波帶的厚度并不均勻,通常有2~5層激波單元。由于數(shù)值誤差的存在,激波帶往往存在很多“毛刺”和“空穴”,為了使后續(xù)聚類分析更加準(zhǔn)確,從兩方面對原始辨識的激波帶進行了光順優(yōu)化。首先一方面標(biāo)記“空穴單元”為激波單元,另一方面剔除“毛刺單元”。所謂空穴單元是指其所有節(jié)點均落在激波帶上的初始非激波單元(圖3中綠色單元),而毛刺單元是指共面相鄰激波單元的個數(shù)小于2的初始激波單元(圖3中藍色單元)。

    圖3 定常激波反射:激波帶的優(yōu)化

    1.2 激波單元的聚類分析

    聚類分析是一種十分重要的統(tǒng)計數(shù)據(jù)分析方法,在數(shù)據(jù)挖掘、模式識別、圖像分析和機器學(xué)習(xí)等許多領(lǐng)域受到廣泛應(yīng)用;它的目標(biāo)是將數(shù)據(jù)集分成許多簇,使得同一簇內(nèi)的數(shù)據(jù)點相似度盡可能大,而不同簇間的數(shù)據(jù)點相似度盡可能小[20]。本文采用經(jīng)典的K-means聚類(K-Means Clustering)[21]算法對優(yōu)化后的激波帶進行聚類處理,下面舉例簡要介紹該算法對圖4(a)數(shù)據(jù)點的聚類流程:

    圖4 K-means聚類算法示意圖

    1) 選擇K個初始“聚類中心”。在圖4(b)中,初始化選擇了兩個聚類中心,即圖中紅、藍“×”形點。

    2) 對任意一個數(shù)據(jù)點,求其到K個聚類中心的“距離”(取歐氏距離),將該數(shù)據(jù)點歸到距離其最近的聚類中心所在的“簇”。如圖4(c)所示,所有數(shù)據(jù)點在經(jīng)過第一輪遍歷后被分為紅、藍標(biāo)識的兩個簇。

    3) 利用均值等方法更新K個簇的中心值。此時對當(dāng)前標(biāo)記為紅色和藍色的數(shù)據(jù)點分別求出其新的質(zhì)心,作為新的聚類中心,如圖4(d)所示,紅色和藍色聚類中心的位置發(fā)生了變動。

    4) 重復(fù)步驟2)~3),若所有聚類中心在迭代更新后位置都保持不變,則迭代結(jié)束;否則繼續(xù)迭代。最終得到的兩個簇及其聚類中心如圖4(f)所示。

    K-means聚類算法操作簡單,易于實現(xiàn),且十分高效,一般僅需5~10次迭代即可收斂。但是研究表明,該算法對初始聚類中心的位置選取比較敏感,不同的隨機初始聚類中心得到的聚類結(jié)果可能差異顯著。因此,在對激波單元聚類時,合理的初始聚類中心位置顯得至關(guān)重要,其既不能偏離激波帶過遠,也不能分布的過密或過稀。下面介紹下針對激波帶初始聚類中心的選取方法。

    首先任意選取一激波單元作為搜索起點,以該單元格心為圓心,n倍當(dāng)?shù)鼐W(wǎng)格尺寸為半徑r作“搜索圓”(經(jīng)大量算例評估,一般n取6~7即可獲得合理的結(jié)果),如圖5所示;隨后將格心點落在該圓區(qū)域內(nèi)且未被搜尋到的激波單元歸為一簇,計算該簇中所有激波單元格心坐標(biāo)的平均值,作為新的搜索圓圓心坐標(biāo);反復(fù)迭代直到所有激波單元均被搜尋到,則所有搜索圓的圓心即為初始聚類中心。

    圖5 定常激波反射:確定初始聚類中心

    取所有激波單元的格心作為數(shù)據(jù)點集,對所有激波單元進行K-means聚類;經(jīng)過若干次迭代后最終所有激波單元被劃分到不同的簇,如圖6所示,為方便顯示,緊鄰的簇用不同顏色加以區(qū)分。

    圖6 定常激波反射:K-means聚類結(jié)果

    為了后續(xù)便于分析激波拓?fù)浣Y(jié)構(gòu)及其干擾模式,根據(jù)每個簇的相鄰特征將簇進一步分成4類:

    1) 終止簇。若流場內(nèi)部某簇僅有1個相鄰簇,則該簇為終止簇;其一般位于流場內(nèi)部激波起始/終止處,注意該簇不與邊界相鄰。

    2) 普通簇。若流場內(nèi)部某簇僅有2個相鄰簇,則該簇為普通簇;其一般呈“串狀”普遍分布于激波帶中,是數(shù)量最多的簇。

    3) 分叉簇。若流場內(nèi)部某簇有3個及以上相鄰簇,則該簇為分叉簇;其一般分布于厚度較大的激波帶處,如三波點、四波點和弱激波處等。

    4) 邊界簇。若某簇直接與流場邊界相鄰,則該簇為邊界簇;其通常存在于激波-壁面反射處、壁面拐角處以及流場進出口邊界處等。

    其中,終止簇、分叉簇和邊界簇統(tǒng)稱關(guān)鍵簇,其在流場中個數(shù)雖少,卻往往位于激波拓?fù)浠蛐螤畎l(fā)生變化的關(guān)鍵位置處。

    1.3 激波的模式識別與擬合

    如何自動準(zhǔn)確地辨識出流場中激波干擾點的位置是一個難題。下面通過對激波帶K-means聚類得到的各個簇進行“近鄰分析”,給出了一種激波模式識別和激波線擬合的策略:

    1) 判斷是否存在需要進行合并的簇。具體分兩方面,若某分叉簇與其他分叉簇相鄰,則將這些分叉簇合并;若某邊界簇與其他邊界簇(注意它們必須對應(yīng)同一個邊界)相鄰,則將這些邊界簇合并。

    2) 若進行了簇合并的操作,需要重新對簇進行分類,并計算新簇的聚類中心。如圖6(a)無需進行簇合并,而圖6(b)經(jīng)簇合并后變成圖7,即原4個相鄰分叉簇合并后形成1個新的普通簇,一同緊靠下壁面的2個相鄰邊界簇合并后形成1個新的邊界簇。

    圖7 定常激波反射:簇合并后

    3) 識別關(guān)鍵激波點。規(guī)定若某分叉簇與周邊3個簇(見圖6(a))或4個簇(見圖23)相鄰,則該分叉簇對應(yīng)的聚類中心即為三波點或四波點;若某分叉簇僅與1個簇相鄰,則其聚類中心即為激波終止/起始點(見圖25);若某邊界簇與2個簇相鄰(圖7下邊界處),則取該邊界簇對應(yīng)的邊界中心點作為激波-壁面反射點;同理,若某邊界簇僅與1個簇相鄰(圖7右邊界處),則邊界中心點即為激波-邊界干擾點。通過以上近鄰分析,便可識別出激波-激波干擾或激波-邊界干擾等模式。

    4) 記錄各條激波所包含的簇。從第3)步得到的某一關(guān)鍵簇開始搜索,依次記錄下該分支上兩兩相鄰的所有普通簇,直至搜索到另一關(guān)鍵簇停止,此時該一系列簇即對應(yīng)一條激波;繼續(xù)搜索確定下一激波分支上的簇,直到遍歷完流場中所有的簇。

    5) 擬合激波線。依次連接激波包含的一系列簇所對應(yīng)的聚類中心,便可得到各條激波線;然而此時得到的激波線可能比較曲折(圖7中的白色實線),為此,本文進一步采用著名的Bézier曲線擬合算法[22]獲取更加光滑的激波線,下面對其作簡要介紹。

    當(dāng)給定一系列有序散點P0、P1、…、Pn時,則其n階Bézier曲線的一般參數(shù)化形式為

    (4)

    如4個點對應(yīng)的三階Bézier曲線表達式為

    B3(λ)=P0(1-λ)3+3P1λ(1-λ)2+

    3P2λ2(1-λ)+P3λ3λ∈[0,1]

    (5)

    式中:Pi稱為Bézier曲線的控制點,且該曲線開始于點P0(對應(yīng)λ=0)并結(jié)束于點Pn(對應(yīng)λ=1)。值得注意的是,Bézier曲線是有方向性的,控制點的次序不同會擬合得到不同的曲線;且Bézier曲線是直線的充分必要條件是所有控制點共線。

    因此,根據(jù)Bézier曲線的特點,若第4)步記錄的某條激波包含n+1個簇,則將兩個關(guān)鍵簇的聚類中心分別作為起始點P0和終止點Pn,而將其余依次緊鄰的普通簇的聚類中心作為控制點P1~Pn-1,進行n階Bézier曲線擬合便得到更加光滑的激波線,如圖8所示。

    圖8 定常激波反射:擬合的激波線

    圖9對比了最終擬合的激波線與流場壓強云圖,4條激波線、三波點、激波-壁面反射點和激波-邊界干擾點的位置都比較準(zhǔn)確,顯示出該算法具有良好的有效性和準(zhǔn)確性。

    圖9 定常激波反射:激波線擬合結(jié)果與云圖對比

    1.4 算法小結(jié)

    該二維激波模式識別算法的總體流程如圖10所示,主要包含激波辨識、聚類分析和識別擬合3方面內(nèi)容;此三者層層遞進,激波辨識是基礎(chǔ),聚類分析是核心,而識別擬合是關(guān)鍵。

    圖10 算法流程圖

    從激波捕捉流場出發(fā),通過流向壓強梯度閾值法進行二級濾波辨識提取出具有一定厚度的激波帶,到最終準(zhǔn)確識別出激波-激波干擾點和激波-邊界干擾點等關(guān)鍵位置,且擬合出的各條激波線的形狀和位置都具有較高的精度。總之,該算法的最大優(yōu)勢是基本不需要人工參與,且算法簡單高效,可以準(zhǔn)確實現(xiàn)從“激波帶”到“激波線”的擬合識別。

    2 算例驗證

    通過4個超聲速流動算例從多方面考察了該算法的準(zhǔn)確性、適用性和可靠性。以下激波捕捉流場均基于格心型有限體積法[23]計算,空間對流項通量采用van Leer格式,時間推進采用顯式4步Runge-Kutta格式,保證在光滑流場中時空均具有二階精度,全場均不考慮黏性。

    2.1 定常斜激波

    首先,采用兩種條件下的定常斜激波算例,驗證本算法最終擬合出的激波線的準(zhǔn)確性。計算區(qū)域為x∈[0,L],y∈[0,0.9L],在x=0.1L處存在楔角,均勻超聲速來流經(jīng)過楔面會產(chǎn)生一道筆直的斜激波,上邊界為超聲速出口,激波不會反射。為了產(chǎn)生相對較弱和較強的斜激波,楔面角θw分別取3°和20°,對應(yīng)的來流馬赫數(shù)Ma∞分別為1.469和1.959。計算網(wǎng)格尺寸均為0.02L,全場分別有5 257個和4 472個三角形網(wǎng)格。

    兩種狀態(tài)下的斜激波經(jīng)辨識和聚類分析的結(jié)果如圖11所示,弱激波辨識出的激波帶較寬,最終一共得到25個聚類中心;而強激波辨識出的激波帶更加銳利,最終一共得到了22個聚類中心。將這些聚類中心作為控制點按順序代入式(4)擬合得到Bézier曲線,并與相應(yīng)的理論值進行對比(圖12),可見擬合值具有較高的位置精度。表1進一步對比了激波角,可見無論是弱激波還是強激波,最終擬合出的激波線的激波角度和理論值的相對偏差都不到1%。

    圖11 定常斜激波:聚類結(jié)果

    圖12 定常斜激波:激波線擬合值與理論值對比

    表1 定常斜激波:激波角擬合值和理論值對比

    2.2 定常激波反射

    仍考慮圖1所示的定常激波反射流動,測試本算法對不同類型網(wǎng)格的可靠性。首先分別基于全場均勻的四邊形結(jié)構(gòu)網(wǎng)格和三角形/四邊形混合網(wǎng)格重新進行數(shù)值計算,平均網(wǎng)格尺寸均為0.015L。 圖13給出了這兩種類型網(wǎng)格下局部流場壓強等值線云圖(圖例與圖1一致)。

    圖13 定常激波反射:不同種類網(wǎng)格下的壓強云圖

    根據(jù)1.1節(jié)所述的激波探測方法,采用相同的濾波系數(shù),分別獲得滿足條件的激波帶。在對激波單元進行K-means聚類分析時,結(jié)構(gòu)網(wǎng)格算例全場一共確定了64個初始聚類中心,而混合網(wǎng)格算例共搜尋到了82個,經(jīng)若干步后得到收斂的簇及其聚類中心。此時,根據(jù)簇的類別和相鄰簇的個數(shù)判斷某些簇是否需要進行合并,圖14和圖15分別給出了兩種網(wǎng)格下壁面馬赫反射處和正規(guī)反射處聚類合并后的結(jié)果。

    圖14(a)和圖15(a)中各存在一個分叉簇,因其與周圍3個激波帶分支中的3個普通簇相鄰,可以判定該分叉簇的聚類中心位于流場三波點附近。此外,反射激波S3在受到膨脹波干擾后,強度顯著減弱;觀察圖14(b)和圖15(b)壁面激波反射處的激波帶,其厚度要明顯大于三波點附近激波帶的厚度,相應(yīng)地,各簇包含的單元個數(shù)也明顯增加,但并沒有存在錯誤的分叉簇,且下壁面和出口邊界處的兩個邊界簇都被準(zhǔn)確地識別出來。值得注意的是,由于激波-壁面反射點比較靠近出口邊界,圖14(b)中存在兩個邊界簇直接相鄰的情況,但由于它們鄰近的邊界不同,因此在簇合并操作中并未錯誤地將它們合并。

    圖14 定常激波反射:基于結(jié)構(gòu)網(wǎng)格的聚類結(jié)果

    圖15 定常激波反射:基于混合網(wǎng)格的聚類結(jié)果

    以上采用的都是均勻網(wǎng)格,下面進一步考察了本算法在非均勻網(wǎng)格中的適用性。如圖16所示,在非結(jié)構(gòu)網(wǎng)格的基礎(chǔ)上,沿入射激波的法向進行加密,第一層網(wǎng)格高度為0.002L,向激波兩側(cè)各推進5層四邊形網(wǎng)格,網(wǎng)格高度增長率取1.2,隨后將激波右側(cè)四邊形網(wǎng)格沿對角線剖分成三角形網(wǎng)格。圖17顯示了對辨識出的激波帶進行聚類分析后的結(jié)果,雖然激波帶厚度差異明顯,但并不影響分叉簇的位置識別。

    圖16 定常激波反射:非均勻網(wǎng)格下的壓強云圖

    圖17 定常激波反射:基于非均勻網(wǎng)格的聚類結(jié)果

    最后對比了以上3種網(wǎng)格下最終的擬合結(jié)果(圖18),4條擬合激波線的形狀和位置都與數(shù)值捕捉的激波吻合得很好;從局部放大圖看出,三者擬合出的三波點和邊界干擾點位置稍有不同,這是由于流場數(shù)值誤差和初始激波帶辨識差異造成的,難以完全根除,但偏差都在一個網(wǎng)格尺度左右,仍顯示出較好的一致性。

    2.3 同側(cè)和異側(cè)激波干擾

    對波系較為復(fù)雜的多激波干擾問題進行了數(shù)值驗證,幾何模型和流場壓強等值線如圖19所示,全場采用均勻的非結(jié)構(gòu)Delaunay三角形網(wǎng)格離散,平均網(wǎng)格尺寸Δ1=0.015L。Ma∞=3.0的高超聲速均勻氣流在經(jīng)過收縮管道時,經(jīng)上壁面20°楔角壓縮產(chǎn)生了一道斜激波S1,與此同時,先后經(jīng)下壁面10°和20°楔角二次壓縮產(chǎn)生了兩道斜激波S2和S3,隨后S2和S3發(fā)生了同側(cè)正規(guī)激波相交并匯聚成了一道新激波S4,最終激波S1與S4發(fā)生了異側(cè)正規(guī)相交并產(chǎn)生了兩條透射激波S5和S6。全場一共產(chǎn)生了六道激波,并存在一個三波點和一個四波點結(jié)構(gòu)。

    圖19 激波干擾:幾何模型與壓強云圖

    根據(jù)流向壓強梯度過濾得到的激波帶分布如圖20所示,在兩個激波干擾處聚集了大量的激波單元。通過進一步增補原激波帶外緣兩側(cè)的空穴單元,優(yōu)化后的激波帶更加光順。隨后對優(yōu)化的激波帶進行K-means聚類,結(jié)果如圖21(a)所示,在每條激波帶分支中都分布著大量呈串狀依次排列的普通簇,而在分支交匯處卻聚集著若干分叉簇,這些分叉簇都至少與周圍3個簇直接相鄰,因此需要進行簇合并操作。

    圖20 激波干擾:辨識出的激波帶

    在將某些簇合并后(圖21(b)),清晰地存在兩個較大的分叉簇,分別與周邊3和4個普通簇相鄰,因此其聚類中心可以被視為三波點和四波點位置。最終分別將各個分支中的聚類中心進行擬合得到了所有6條激波線,結(jié)果如圖22所示。擬合的激波線、激波相交點均落在數(shù)值捕捉激波的過渡區(qū)內(nèi),顯示出該算法良好的精度。

    圖22 激波干擾:擬合的激波線

    實際上,激波(尤其是激波-激波干擾點)的探測精度在很大程度上取決于流場的分辨率。也就是說,捕捉得到的激波過渡區(qū)寬度越小,越有利于精確探測到激波及其干擾點的位置。下面考察了網(wǎng)格的疏密程度對激波探測的影響。

    在其他相同計算條件下,僅僅改變流場網(wǎng)格的尺寸,分別設(shè)計疏(Δ2=0.02L)和密(Δ3=0.01L)兩套非結(jié)構(gòu)網(wǎng)格,重新進行流場計算。經(jīng)激波辨識、聚類分析和簇合并后的結(jié)果如圖23所示;網(wǎng)格越稀,劃分的簇越少,且聚類中心兩兩相距越遠。圖24對比了3種網(wǎng)格尺寸下最終擬合出的激波線,可以發(fā)現(xiàn),三者大部分基本完全吻合,但激波干擾點處差異較大。對于同側(cè)激波干擾結(jié)構(gòu),由于三道激波的斜率都比較接近,因此當(dāng)網(wǎng)格較稀時,數(shù)值捕捉激波的過渡區(qū)較寬,會造成兩條入射激波帶提前交匯(對比圖23中兩圖可見),分叉簇的位置會發(fā)生較大改變,最終造成了辨識出的三波點位置出現(xiàn)較大差異。另一方面,對于異側(cè)激波干擾產(chǎn)生的四波點結(jié)構(gòu),雖然稀網(wǎng)格對應(yīng)的分叉簇也較大,但其聚類中心的位置和密網(wǎng)格的結(jié)果偏差并不顯著。

    圖23 激波干擾:不同網(wǎng)格尺寸下的聚類結(jié)果

    圖24 激波干擾:3種網(wǎng)格尺寸下擬合的激波線對比

    該算例表明,本文算法對同側(cè)激波干擾點的識別精度網(wǎng)格依賴性相對較強,即網(wǎng)格疏密對結(jié)果的影響較大;而對于異側(cè)激波干擾點,包括上文激波反射算例中的三波點結(jié)構(gòu),其識別精度受網(wǎng)格尺度的影響較??;而對于干擾點以外的激波線,擬合精度受網(wǎng)格尺度的影響最小,始終吻合得比較好。

    2.4 前向臺階非定常流動

    采用著名的前向臺階超聲速繞流算例[24]測試該算法在復(fù)雜非定常流動中的有效性。均勻自由來流馬赫數(shù)Ma∞=3.0,且初始流場取來流參數(shù),通道入口和出口處的高度分別為L和0.8L,臺階高度及其上表面長度分別為0.2L和2.4L。全場采用均勻非結(jié)構(gòu)三角形網(wǎng)格劃分,網(wǎng)格尺度為0.01L。本算例中的時間均為無量綱量。

    在流場的發(fā)展歷程中,臺階前首先出現(xiàn)一道向左運動的弓形激波,因其強度較大,當(dāng)受到上壁面干擾時會發(fā)生反射,并從運動正規(guī)反射逐漸過渡到馬赫反射;而反射激波又從下壁面反射回上壁面,最終在管道壁面形成了4次激波-壁面干擾(Shock-Wall Interaction,SWI)。此外,氣流在臺階拐角處發(fā)生過度膨脹,過膨脹氣流會與下壁面相撞形成一道相對較弱的孤立斜激波,該激波最終與第一道反射激波發(fā)生異側(cè)正規(guī)相交,相交后的弱激波又與第二道反射激波匯聚形成同側(cè)激波-激波干擾(Shock-Shock Interaction,SSI)結(jié)構(gòu)。

    圖25依次給出了5個關(guān)鍵時刻的辨識激波帶的聚類結(jié)果,通過二次濾波,在辨識出運動強激波的同時,那道較弱的孤立斜激波也被很好地提取出來;此外,圖25(d)和25(e)中并沒有錯誤地將三波點處的接觸間斷提取出來。隨后,通過聚類分析和簇合并操作,可以準(zhǔn)確地定位分叉簇與邊界簇的位置,從而清晰地識別出各個時刻流場激波的總個數(shù)及相互干擾的模式(如SWI或SSI)。

    圖25 前向臺階繞流:不同時刻激波帶的聚類結(jié)果

    以兩時刻為例,圖26定性比較了最終擬合出的激波線與流場壓強云圖??梢园l(fā)現(xiàn),不僅各條近似直線的反射激波位置準(zhǔn)確,而且臺階前的大曲率弓形激波也與捕捉激波吻合得很好,反映了該算法對各種形狀激波的擬合都具有較好的可靠性。在得到不同時刻的擬合激波線后,可以方便地將其在同一張圖中顯示出來,進而有利于分析激波的演化過程,深化對運動激波干擾的認(rèn)識。從圖27可以看出,相比激波-上壁面干擾點,臺階后的激波-下壁面干擾點向左移動的速度更快;且弓形激波在上壁面反射形成的馬赫干擾的長度在迅速增大;而在t=1.5時刻后,弓形激波的位置變化比較緩慢,三波點近似沿弓形激波向下滑動。

    圖26 前向臺階繞流:擬合的激波線與壓強云圖對比

    圖27 前向臺階繞流:4個時刻擬合的激波線

    3 結(jié) 論

    1) 該算法在傳統(tǒng)基于流場當(dāng)?shù)貐?shù)進行激波探測的方法基礎(chǔ)上,對提取的激波單元進一步處理;先后結(jié)合K-means聚類算法和近鄰分析成功實現(xiàn)了激波干擾點的準(zhǔn)確定位;最后采用Bézier曲線算法擬合得到各條質(zhì)量較高的激波線,數(shù)值試驗表明激波線的大部分位置與形狀受流場網(wǎng)格疏密影響不大。

    2) 對于同側(cè)激波干擾,干擾點的位置受流場本身分辨率的影響較大,密網(wǎng)格下一般會獲得更為精準(zhǔn)的結(jié)果;而對于異側(cè)激波干擾,其干擾點的識別結(jié)果通常比較準(zhǔn)確,且網(wǎng)格依賴性不強。

    3) 該算法自動化程度較高,基本不需要人工干預(yù),且適用于任意網(wǎng)格類型,對復(fù)雜非定常流場中的運動激波辨識也具有良好的可靠性和準(zhǔn)確性。

    總之,該面向全局的激波探測算法可以有效識別出更為精確的激波拓?fù)浣Y(jié)構(gòu),未來可以用于CFD中某些對激波位置有特殊需求的技術(shù)中,例如激波裝配[16-19]、網(wǎng)格自適應(yīng)[25]、流場特征可視化[26]等。由于三維激波相交干擾結(jié)構(gòu)比較復(fù)雜,將該方法推廣到三維是今后努力的方向。

    猜你喜歡
    波點激波壁面
    二維有限長度柔性壁面上T-S波演化的數(shù)值研究
    血液動力學(xué)中血管流激波與駐波的相互作用
    基于HIFiRE-2超燃發(fā)動機內(nèi)流道的激波邊界層干擾分析
    斜激波入射V形鈍前緣溢流口激波干擾研究
    讓注意力到你身上來 波點的世界怎能錯過
    波點之美
    適于可壓縮多尺度流動的緊致型激波捕捉格式
    波點女孩
    壁面溫度對微型內(nèi)燃機燃燒特性的影響
    頑趣波點
    Coco薇(2016年3期)2016-04-06 02:40:48
    99久久综合免费| 美女cb高潮喷水在线观看| 成人无遮挡网站| 免费在线观看成人毛片| 亚洲精品日本国产第一区| 久久久久久久久久久丰满| 午夜视频国产福利| 日本免费在线观看一区| 草草在线视频免费看| 国产高清国产精品国产三级| 免费在线观看成人毛片| 伦精品一区二区三区| av天堂中文字幕网| 99热国产这里只有精品6| 精品一区在线观看国产| 亚洲精品日本国产第一区| 日本91视频免费播放| 亚洲国产日韩一区二区| 精品国产乱码久久久久久小说| 18禁在线播放成人免费| 久久精品夜色国产| 久久亚洲国产成人精品v| 精品熟女少妇av免费看| 国产精品人妻久久久影院| 好男人视频免费观看在线| 晚上一个人看的免费电影| 欧美老熟妇乱子伦牲交| 国产一区二区三区av在线| 国产精品一区二区在线不卡| 丰满饥渴人妻一区二区三| 赤兔流量卡办理| 国产女主播在线喷水免费视频网站| 色吧在线观看| 国内精品宾馆在线| 亚洲丝袜综合中文字幕| 永久网站在线| 日韩中文字幕视频在线看片| 在线 av 中文字幕| av视频免费观看在线观看| 亚洲av欧美aⅴ国产| 亚州av有码| 性色avwww在线观看| 久久久久久人妻| 九草在线视频观看| 精品久久久久久电影网| 国产亚洲精品久久久com| 久热久热在线精品观看| 亚洲性久久影院| 亚洲精品亚洲一区二区| 成人漫画全彩无遮挡| 乱系列少妇在线播放| 久久99热6这里只有精品| 国产午夜精品一二区理论片| 这个男人来自地球电影免费观看 | 美女内射精品一级片tv| 精品亚洲成国产av| 人体艺术视频欧美日本| 一级毛片 在线播放| av国产久精品久网站免费入址| 亚洲av不卡在线观看| 久久亚洲国产成人精品v| 久久热精品热| 国产精品欧美亚洲77777| 欧美 亚洲 国产 日韩一| av在线app专区| 国产在线视频一区二区| 免费播放大片免费观看视频在线观看| 午夜91福利影院| 熟妇人妻不卡中文字幕| 看十八女毛片水多多多| 日韩中字成人| 日本爱情动作片www.在线观看| 啦啦啦啦在线视频资源| 精品一品国产午夜福利视频| 久久久a久久爽久久v久久| 少妇人妻久久综合中文| 美女主播在线视频| 日本av免费视频播放| 人人妻人人澡人人看| 国产一区亚洲一区在线观看| 欧美国产精品一级二级三级 | 在线观看一区二区三区激情| 欧美97在线视频| 99久久精品一区二区三区| 人人妻人人看人人澡| 高清视频免费观看一区二区| 久久久久久久久久久免费av| 成人免费观看视频高清| 黄色日韩在线| 精品人妻偷拍中文字幕| 免费久久久久久久精品成人欧美视频 | 两个人免费观看高清视频 | 丰满饥渴人妻一区二区三| 亚洲无线观看免费| 日韩精品有码人妻一区| 我要看黄色一级片免费的| 在线观看人妻少妇| 国产永久视频网站| 9色porny在线观看| 亚洲国产成人一精品久久久| 亚洲成人av在线免费| 亚洲伊人久久精品综合| 亚洲国产最新在线播放| 大又大粗又爽又黄少妇毛片口| 亚洲欧美一区二区三区黑人 | 久久久欧美国产精品| 亚洲av国产av综合av卡| 中文字幕av电影在线播放| 丝袜在线中文字幕| 国产精品一区二区三区四区免费观看| 一本—道久久a久久精品蜜桃钙片| 一区二区三区乱码不卡18| av在线老鸭窝| 精品人妻偷拍中文字幕| 国产精品女同一区二区软件| 久久精品国产亚洲av天美| 99精国产麻豆久久婷婷| 一边亲一边摸免费视频| 国产亚洲一区二区精品| 成人免费观看视频高清| 国产男女内射视频| 免费久久久久久久精品成人欧美视频 | 国产午夜精品久久久久久一区二区三区| 我要看日韩黄色一级片| 如何舔出高潮| 亚洲第一区二区三区不卡| 午夜激情久久久久久久| 国产日韩欧美视频二区| 内射极品少妇av片p| 精品国产乱码久久久久久小说| 18禁在线无遮挡免费观看视频| 伊人久久国产一区二区| 男女啪啪激烈高潮av片| 啦啦啦中文免费视频观看日本| 国产乱人偷精品视频| 伦理电影大哥的女人| 日日撸夜夜添| 久久精品国产亚洲网站| 精品少妇黑人巨大在线播放| 亚洲伊人久久精品综合| 麻豆乱淫一区二区| 夫妻性生交免费视频一级片| 黑人猛操日本美女一级片| 免费在线观看成人毛片| 老司机影院毛片| 亚洲自偷自拍三级| 超碰97精品在线观看| 久久久久国产网址| 亚洲精品乱久久久久久| 日韩一区二区三区影片| 国产熟女欧美一区二区| 国产亚洲最大av| 在线观看www视频免费| www.av在线官网国产| 一区二区av电影网| 最近的中文字幕免费完整| 欧美人与善性xxx| 人人妻人人澡人人看| 国产精品嫩草影院av在线观看| 在线观看免费视频网站a站| 婷婷色麻豆天堂久久| 日日摸夜夜添夜夜爱| 伊人久久国产一区二区| 2022亚洲国产成人精品| 能在线免费看毛片的网站| 国产淫语在线视频| 99热国产这里只有精品6| 热99国产精品久久久久久7| 赤兔流量卡办理| 成人黄色视频免费在线看| 日韩欧美精品免费久久| 国产精品人妻久久久久久| 久久精品国产自在天天线| 亚洲国产最新在线播放| 国产精品国产av在线观看| 亚洲美女黄色视频免费看| av卡一久久| 日韩视频在线欧美| 搡女人真爽免费视频火全软件| 午夜福利影视在线免费观看| 国产又色又爽无遮挡免| 国产一区二区在线观看av| 欧美精品亚洲一区二区| 免费看光身美女| 成人国产麻豆网| 中文字幕精品免费在线观看视频 | av在线app专区| 精品人妻偷拍中文字幕| av国产久精品久网站免费入址| av视频免费观看在线观看| 精品久久久久久久久亚洲| videossex国产| 久久精品熟女亚洲av麻豆精品| 免费高清在线观看视频在线观看| 一二三四中文在线观看免费高清| 在线观看免费视频网站a站| 亚洲熟女精品中文字幕| 在现免费观看毛片| av有码第一页| 国产高清国产精品国产三级| 99热这里只有是精品50| 欧美日韩一区二区视频在线观看视频在线| 亚洲精品久久久久久婷婷小说| 日韩av不卡免费在线播放| 桃花免费在线播放| 久久99热这里只频精品6学生| 美女cb高潮喷水在线观看| 亚洲第一av免费看| 成人亚洲精品一区在线观看| 制服丝袜香蕉在线| 97精品久久久久久久久久精品| videos熟女内射| 欧美人与善性xxx| 国产欧美日韩综合在线一区二区 | 美女福利国产在线| 晚上一个人看的免费电影| 亚洲久久久国产精品| 久久99一区二区三区| 亚洲av男天堂| 最后的刺客免费高清国语| 久久韩国三级中文字幕| 免费在线观看成人毛片| 狠狠精品人妻久久久久久综合| 赤兔流量卡办理| 久久热精品热| 免费av中文字幕在线| 精品一区二区三区视频在线| 国产色爽女视频免费观看| 九九在线视频观看精品| 七月丁香在线播放| 丰满人妻一区二区三区视频av| 免费大片黄手机在线观看| 美女大奶头黄色视频| 欧美精品国产亚洲| 中文字幕亚洲精品专区| 永久网站在线| av卡一久久| 精品亚洲成国产av| 国产精品秋霞免费鲁丝片| 国产黄色视频一区二区在线观看| 色网站视频免费| 久久久久人妻精品一区果冻| 中文字幕久久专区| 亚洲综合色惰| 亚洲美女搞黄在线观看| 人体艺术视频欧美日本| 国产成人免费无遮挡视频| 青春草视频在线免费观看| 伊人亚洲综合成人网| 一二三四中文在线观看免费高清| 在线精品无人区一区二区三| 久久久久久伊人网av| 人妻系列 视频| 老司机影院毛片| 伦理电影免费视频| 国产色婷婷99| 色吧在线观看| 国产精品久久久久久久久免| 一级,二级,三级黄色视频| 国产免费一区二区三区四区乱码| 国产精品99久久久久久久久| 国产视频首页在线观看| 在线观看人妻少妇| av福利片在线| 免费在线观看成人毛片| 男女边摸边吃奶| 性色avwww在线观看| 伦精品一区二区三区| 欧美精品国产亚洲| 91成人精品电影| av免费观看日本| 国产美女午夜福利| 你懂的网址亚洲精品在线观看| 一区二区三区四区激情视频| 国产精品麻豆人妻色哟哟久久| 黄片无遮挡物在线观看| 国产精品99久久99久久久不卡 | 成人亚洲欧美一区二区av| 国产亚洲午夜精品一区二区久久| 国产深夜福利视频在线观看| 精品久久久久久久久亚洲| 下体分泌物呈黄色| 全区人妻精品视频| 男男h啪啪无遮挡| 99热这里只有是精品50| 日日撸夜夜添| 涩涩av久久男人的天堂| 少妇丰满av| 亚洲国产精品成人久久小说| 久久av网站| 亚洲av成人精品一二三区| 观看美女的网站| 亚洲欧美精品自产自拍| 草草在线视频免费看| 精品一区二区三区视频在线| 国产探花极品一区二区| 久久久久久久久久久久大奶| 青青草视频在线视频观看| 国产中年淑女户外野战色| 精品一区二区三卡| 伊人亚洲综合成人网| 少妇裸体淫交视频免费看高清| 亚洲国产成人一精品久久久| 久久久久网色| 久久精品熟女亚洲av麻豆精品| 在线观看免费日韩欧美大片 | 日韩制服骚丝袜av| 国产免费一区二区三区四区乱码| 在线观看人妻少妇| 人妻夜夜爽99麻豆av| 又粗又硬又长又爽又黄的视频| a 毛片基地| 久久久精品免费免费高清| 国产黄色免费在线视频| 一级二级三级毛片免费看| 晚上一个人看的免费电影| 国产精品一区二区在线观看99| 一级毛片 在线播放| 女性被躁到高潮视频| 国产欧美日韩精品一区二区| 国产精品无大码| 久久国产亚洲av麻豆专区| 日本91视频免费播放| 永久免费av网站大全| 全区人妻精品视频| 精品国产乱码久久久久久小说| 精品国产国语对白av| 一级毛片aaaaaa免费看小| 日本色播在线视频| 99久久综合免费| 精品一区二区三卡| 如何舔出高潮| 国产男女内射视频| 亚洲成人手机| 熟妇人妻不卡中文字幕| 日本黄色日本黄色录像| 久久久久久久久久成人| 99热国产这里只有精品6| 最后的刺客免费高清国语| 成年女人在线观看亚洲视频| 亚洲精品日韩在线中文字幕| 国产欧美另类精品又又久久亚洲欧美| 国产乱来视频区| 黑丝袜美女国产一区| 高清av免费在线| 免费黄色在线免费观看| 少妇高潮的动态图| 99久久中文字幕三级久久日本| 三级国产精品片| 少妇高潮的动态图| 日韩av免费高清视频| 国产欧美日韩一区二区三区在线 | 黑丝袜美女国产一区| 各种免费的搞黄视频| 国产无遮挡羞羞视频在线观看| 国产精品一区二区在线观看99| √禁漫天堂资源中文www| 免费黄网站久久成人精品| 人妻夜夜爽99麻豆av| 一个人看视频在线观看www免费| 69精品国产乱码久久久| 亚洲精品久久午夜乱码| 啦啦啦啦在线视频资源| 精品国产一区二区三区久久久樱花| 亚洲自偷自拍三级| 热re99久久精品国产66热6| 国产在线免费精品| 日韩三级伦理在线观看| 人妻系列 视频| 免费大片18禁| 精品一区二区三区视频在线| 亚洲av二区三区四区| 国产成人精品婷婷| 色视频在线一区二区三区| 另类精品久久| 波野结衣二区三区在线| 欧美xxxx性猛交bbbb| 性高湖久久久久久久久免费观看| 国产成人91sexporn| 一区二区三区免费毛片| 午夜免费观看性视频| 欧美日韩精品成人综合77777| 最近手机中文字幕大全| 国产精品人妻久久久久久| 欧美精品一区二区免费开放| 久久久久人妻精品一区果冻| 久久久久久久国产电影| 亚洲精品aⅴ在线观看| 久久综合国产亚洲精品| 久久精品国产亚洲av天美| 久久久久久久久久久久大奶| av国产久精品久网站免费入址| 日本欧美视频一区| 亚洲精品视频女| 能在线免费看毛片的网站| 男女无遮挡免费网站观看| 亚洲精品aⅴ在线观看| 97超视频在线观看视频| 久久人人爽av亚洲精品天堂| 十八禁高潮呻吟视频 | 91久久精品国产一区二区三区| 韩国高清视频一区二区三区| 亚洲欧美日韩东京热| 免费播放大片免费观看视频在线观看| 久久国内精品自在自线图片| 日韩一区二区视频免费看| 免费人妻精品一区二区三区视频| 黑人猛操日本美女一级片| 亚洲国产欧美在线一区| 五月玫瑰六月丁香| 免费看av在线观看网站| 国内少妇人妻偷人精品xxx网站| 精品少妇久久久久久888优播| 在线观看人妻少妇| 日韩欧美 国产精品| 97在线视频观看| 久久久a久久爽久久v久久| 久久国产乱子免费精品| 熟女人妻精品中文字幕| 插逼视频在线观看| 亚洲精品国产成人久久av| 国产日韩欧美亚洲二区| 少妇人妻一区二区三区视频| 亚洲精品第二区| 三级国产精品欧美在线观看| 天美传媒精品一区二区| 久久久国产欧美日韩av| a级毛片在线看网站| 中文字幕精品免费在线观看视频 | 人人妻人人添人人爽欧美一区卜| 久久久久久久久久久丰满| 久久久久久久久久久久大奶| 99九九线精品视频在线观看视频| 我要看日韩黄色一级片| av黄色大香蕉| 婷婷色综合大香蕉| 欧美丝袜亚洲另类| 亚洲国产精品一区二区三区在线| 久久人人爽人人爽人人片va| 丝袜在线中文字幕| av.在线天堂| 精品久久国产蜜桃| 精品人妻一区二区三区麻豆| 国产精品国产三级专区第一集| 日本vs欧美在线观看视频 | 精品久久久噜噜| 深夜a级毛片| 视频中文字幕在线观看| 三上悠亚av全集在线观看 | 国产精品一区www在线观看| av有码第一页| 成人18禁高潮啪啪吃奶动态图 | 亚洲精品乱码久久久久久按摩| 色吧在线观看| 免费在线观看成人毛片| 老司机影院毛片| 欧美少妇被猛烈插入视频| 久久精品国产亚洲av天美| 亚洲一区二区三区欧美精品| 欧美日韩国产mv在线观看视频| av视频免费观看在线观看| 亚洲精品日韩av片在线观看| 一区二区三区乱码不卡18| 国产成人freesex在线| 国产国拍精品亚洲av在线观看| av线在线观看网站| 免费高清在线观看视频在线观看| 精品久久久精品久久久| 久久久久久久久久成人| 国产亚洲午夜精品一区二区久久| 精品亚洲成国产av| 精品国产乱码久久久久久小说| 免费在线观看成人毛片| 人妻夜夜爽99麻豆av| 国产成人a∨麻豆精品| 亚洲久久久国产精品| 免费av中文字幕在线| 男人添女人高潮全过程视频| 久热这里只有精品99| 久久久久久伊人网av| 免费黄网站久久成人精品| 国产 一区精品| 一区二区三区免费毛片| 欧美日本中文国产一区发布| 久久午夜综合久久蜜桃| av黄色大香蕉| 日本午夜av视频| 亚洲精品成人av观看孕妇| 在线天堂最新版资源| 国产免费视频播放在线视频| 久久av网站| 一边亲一边摸免费视频| 国产乱人偷精品视频| 男女啪啪激烈高潮av片| 亚洲四区av| 不卡视频在线观看欧美| 91久久精品国产一区二区三区| 国产精品麻豆人妻色哟哟久久| 日本午夜av视频| 免费在线观看成人毛片| 人妻夜夜爽99麻豆av| 一本久久精品| av天堂中文字幕网| 国产成人免费无遮挡视频| 欧美人与善性xxx| 国产日韩一区二区三区精品不卡 | 亚洲无线观看免费| 好男人视频免费观看在线| 中文字幕人妻熟人妻熟丝袜美| 日日啪夜夜撸| 亚洲精品乱码久久久v下载方式| 一本色道久久久久久精品综合| 国产精品不卡视频一区二区| 精品卡一卡二卡四卡免费| a级毛片在线看网站| 女的被弄到高潮叫床怎么办| 国产精品国产三级国产av玫瑰| 亚洲欧美成人精品一区二区| 热re99久久国产66热| 国产女主播在线喷水免费视频网站| 免费看日本二区| 中国三级夫妇交换| 国产精品国产av在线观看| 精品一区在线观看国产| 国产精品国产av在线观看| 亚州av有码| 色视频在线一区二区三区| 国产精品国产av在线观看| 亚洲精品国产成人久久av| 久久精品熟女亚洲av麻豆精品| 亚洲自偷自拍三级| 国产色爽女视频免费观看| 精品久久久久久电影网| 三上悠亚av全集在线观看 | 少妇猛男粗大的猛烈进出视频| 在线精品无人区一区二区三| 一级毛片久久久久久久久女| 欧美高清成人免费视频www| 成人影院久久| 尾随美女入室| 亚洲av电影在线观看一区二区三区| 赤兔流量卡办理| 亚洲欧洲精品一区二区精品久久久 | 久久99蜜桃精品久久| 久久精品国产亚洲网站| 蜜臀久久99精品久久宅男| 秋霞在线观看毛片| 亚洲一区二区三区欧美精品| 99热6这里只有精品| 国产黄色免费在线视频| 黑人巨大精品欧美一区二区蜜桃 | 久久久久视频综合| 欧美+日韩+精品| 性色avwww在线观看| 免费高清在线观看视频在线观看| videos熟女内射| 青春草亚洲视频在线观看| 免费av不卡在线播放| 国产 精品1| 最黄视频免费看| 日韩精品免费视频一区二区三区 | 成年人免费黄色播放视频 | tube8黄色片| 国内少妇人妻偷人精品xxx网站| 成人18禁高潮啪啪吃奶动态图 | 毛片一级片免费看久久久久| 国产av精品麻豆| 男人和女人高潮做爰伦理| 高清午夜精品一区二区三区| 午夜影院在线不卡| 观看av在线不卡| 婷婷色av中文字幕| 国产精品秋霞免费鲁丝片| 精品视频人人做人人爽| 国产日韩一区二区三区精品不卡 | 欧美3d第一页| 色婷婷av一区二区三区视频| 国产免费又黄又爽又色| 美女中出高潮动态图| 成人国产av品久久久| 大香蕉久久网| av国产久精品久网站免费入址| 在线天堂最新版资源| av不卡在线播放| 国产熟女欧美一区二区| 久久国内精品自在自线图片| 两个人的视频大全免费| 水蜜桃什么品种好| 国产精品福利在线免费观看| 精品一区二区三区视频在线| 大陆偷拍与自拍| 午夜福利影视在线免费观看| 久久人妻熟女aⅴ| 欧美精品亚洲一区二区| 黄片无遮挡物在线观看| 成年av动漫网址| 99热全是精品| 尾随美女入室| 观看美女的网站| 久久精品国产鲁丝片午夜精品| 亚洲综合色惰| 午夜影院在线不卡| 亚洲精品国产成人久久av| 日日摸夜夜添夜夜添av毛片| 久久国产亚洲av麻豆专区| 国产熟女欧美一区二区| 一级毛片我不卡| 成人毛片a级毛片在线播放| 美女主播在线视频| 久久精品夜色国产| 中文欧美无线码| 欧美日韩亚洲高清精品| 欧美xxⅹ黑人| 国产精品一二三区在线看| 欧美亚洲 丝袜 人妻 在线|