• <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
    老汉色av国产亚洲站长工具| 人妻一区二区av| 日本欧美视频一区| 宅男免费午夜| 欧美日本中文国产一区发布| 午夜福利在线免费观看网站| 高清黄色对白视频在线免费看| 午夜免费男女啪啪视频观看| 亚洲欧美日韩另类电影网站| 亚洲欧洲日产国产| 欧美最新免费一区二区三区| 美女国产高潮福利片在线看| 欧美另类一区| 大香蕉久久成人网| 国产午夜精品一二区理论片| 国产亚洲一区二区精品| 亚洲情色 制服丝袜| 最黄视频免费看| 制服丝袜香蕉在线| 亚洲精品久久久久久婷婷小说| av在线观看视频网站免费| 午夜av观看不卡| 99久久精品国产亚洲精品| 少妇的丰满在线观看| 色综合欧美亚洲国产小说| 日本猛色少妇xxxxx猛交久久| 久久午夜综合久久蜜桃| 国产一区二区 视频在线| 在线观看免费高清a一片| 亚洲综合精品二区| 中文字幕av电影在线播放| 国产精品女同一区二区软件| 精品第一国产精品| 最近中文字幕高清免费大全6| av在线播放精品| 国语对白做爰xxxⅹ性视频网站| 在线 av 中文字幕| 免费黄网站久久成人精品| 日本欧美视频一区| 久久久久精品人妻al黑| 肉色欧美久久久久久久蜜桃| 美女脱内裤让男人舔精品视频| 国产黄色免费在线视频| 最近最新中文字幕免费大全7| 日本猛色少妇xxxxx猛交久久| 国产成人av激情在线播放| 高清不卡的av网站| 色播在线永久视频| 一二三四在线观看免费中文在| 女人被躁到高潮嗷嗷叫费观| 男人添女人高潮全过程视频| 国产在线视频一区二区| 51午夜福利影视在线观看| 亚洲av中文av极速乱| 成人国产av品久久久| 国产精品 国内视频| 91国产中文字幕| 极品人妻少妇av视频| 夫妻午夜视频| 99久久综合免费| 老司机在亚洲福利影院| 97在线人人人人妻| 亚洲久久久国产精品| 美女高潮到喷水免费观看| 亚洲国产日韩一区二区| 女性生殖器流出的白浆| 伊人久久国产一区二区| 久久久久久久国产电影| 天堂俺去俺来也www色官网| 国产无遮挡羞羞视频在线观看| 日本欧美国产在线视频| 国产精品国产三级国产专区5o| 飞空精品影院首页| 老司机靠b影院| 美女午夜性视频免费| 捣出白浆h1v1| 电影成人av| 亚洲国产av新网站| 日本av手机在线免费观看| 日本黄色日本黄色录像| 欧美日韩福利视频一区二区| av天堂久久9| 黄色视频在线播放观看不卡| 999久久久国产精品视频| a级片在线免费高清观看视频| 精品视频人人做人人爽| 久久99一区二区三区| bbb黄色大片| svipshipincom国产片| 亚洲五月色婷婷综合| 国产乱人偷精品视频| 午夜免费观看性视频| 美女大奶头黄色视频| 亚洲激情五月婷婷啪啪| 亚洲一级一片aⅴ在线观看| 久久人人爽av亚洲精品天堂| 婷婷色综合www| 国产一区有黄有色的免费视频| 亚洲在久久综合| 国产精品久久久久久人妻精品电影 | 国产精品蜜桃在线观看| 亚洲精品久久成人aⅴ小说| 成人免费观看视频高清| 国产男女超爽视频在线观看| 亚洲成人av在线免费| 精品国产一区二区三区四区第35| 黄色 视频免费看| 久久人人爽人人片av| 777米奇影视久久| 少妇人妻 视频| 国产精品久久久久久人妻精品电影 | 不卡视频在线观看欧美| 国产色婷婷99| 日韩 欧美 亚洲 中文字幕| 大香蕉久久网| 十八禁人妻一区二区| 日韩大片免费观看网站| 高清不卡的av网站| 免费看av在线观看网站| 久久精品熟女亚洲av麻豆精品| 黄色一级大片看看| 丝袜人妻中文字幕| 久久久久国产一级毛片高清牌| 国产精品三级大全| www日本在线高清视频| 国产精品蜜桃在线观看| 人人妻人人爽人人添夜夜欢视频| 超色免费av| 大码成人一级视频| 一本大道久久a久久精品| 亚洲专区中文字幕在线 | 精品少妇内射三级| 国产成人啪精品午夜网站| h视频一区二区三区| 久久毛片免费看一区二区三区| 秋霞伦理黄片| 国产乱人偷精品视频| 999久久久国产精品视频| 99久国产av精品国产电影| 美女福利国产在线| 嫩草影院入口| 天堂8中文在线网| 成人国产麻豆网| 日韩制服骚丝袜av| 亚洲精品久久久久久婷婷小说| 精品国产一区二区久久| 久久久久久久精品精品| 亚洲欧美一区二区三区国产| 久久精品aⅴ一区二区三区四区| 老司机影院毛片| 大香蕉久久成人网| 亚洲国产精品一区三区| 少妇人妻精品综合一区二区| 在线观看一区二区三区激情| 99国产综合亚洲精品| 婷婷色综合大香蕉| 日韩成人av中文字幕在线观看| 狂野欧美激情性bbbbbb| 午夜精品国产一区二区电影| 国产又色又爽无遮挡免| 亚洲七黄色美女视频| 日韩不卡一区二区三区视频在线| 中文字幕最新亚洲高清| 一边摸一边做爽爽视频免费| 考比视频在线观看| 电影成人av| 国产精品人妻久久久影院| 午夜日韩欧美国产| 亚洲av在线观看美女高潮| a级片在线免费高清观看视频| 一级毛片黄色毛片免费观看视频| 欧美日韩视频精品一区| 一级毛片电影观看| 国产亚洲av高清不卡| 国产野战对白在线观看| 黄色视频不卡| 这个男人来自地球电影免费观看 | 制服诱惑二区| 国产一区二区 视频在线| 777久久人妻少妇嫩草av网站| 国产精品麻豆人妻色哟哟久久| 久久久国产精品麻豆| av免费观看日本| 免费av中文字幕在线| 女性生殖器流出的白浆| 天堂中文最新版在线下载| 操美女的视频在线观看| 日日摸夜夜添夜夜爱| 亚洲精品第二区| 久久久精品区二区三区| kizo精华| 亚洲成色77777| 美女主播在线视频| 国精品久久久久久国模美| 热re99久久精品国产66热6| 久久婷婷青草| 在线看a的网站| 狂野欧美激情性xxxx| av国产久精品久网站免费入址| 狠狠精品人妻久久久久久综合| 中文字幕人妻丝袜一区二区 | 欧美日韩视频高清一区二区三区二| 水蜜桃什么品种好| 日韩精品免费视频一区二区三区| 肉色欧美久久久久久久蜜桃| 777米奇影视久久| 我要看黄色一级片免费的| av在线app专区| 看十八女毛片水多多多| 成人国产av品久久久| 国产伦人伦偷精品视频| 亚洲成色77777| 肉色欧美久久久久久久蜜桃| 国产精品久久久久久人妻精品电影 | 婷婷成人精品国产| 在线精品无人区一区二区三| 亚洲成人av在线免费| 99久久精品国产亚洲精品| 亚洲国产欧美网| 老司机在亚洲福利影院| 秋霞伦理黄片| 黄频高清免费视频| 一本—道久久a久久精品蜜桃钙片| 少妇猛男粗大的猛烈进出视频| www.自偷自拍.com| 久久久久久久久久久免费av| 日日撸夜夜添| 国语对白做爰xxxⅹ性视频网站| 欧美日韩精品网址| 久久亚洲国产成人精品v| 街头女战士在线观看网站| 91老司机精品| 纯流量卡能插随身wifi吗| 免费在线观看完整版高清| 免费高清在线观看视频在线观看| 一级毛片 在线播放| 大片免费播放器 马上看| av线在线观看网站| 最近中文字幕2019免费版| 国产麻豆69| 日日摸夜夜添夜夜爱| 国产深夜福利视频在线观看| 精品国产乱码久久久久久男人| 亚洲一级一片aⅴ在线观看| 欧美国产精品va在线观看不卡| 免费黄频网站在线观看国产| 亚洲av日韩精品久久久久久密 | 青春草国产在线视频| 欧美人与善性xxx| 中文欧美无线码| 国产人伦9x9x在线观看| 国产爽快片一区二区三区| 成年人免费黄色播放视频| 黄片无遮挡物在线观看| 美女福利国产在线| 国产野战对白在线观看| 国产女主播在线喷水免费视频网站| 一级片'在线观看视频| 最近的中文字幕免费完整| 国产一级毛片在线| 啦啦啦 在线观看视频| 国产欧美日韩综合在线一区二区| 久久久精品94久久精品| 久久97久久精品| 激情视频va一区二区三区| 中文字幕人妻丝袜一区二区 | xxxhd国产人妻xxx| 国产精品蜜桃在线观看| 中文字幕人妻丝袜制服| 久久韩国三级中文字幕| 深夜精品福利| 久久久精品94久久精品| 男女无遮挡免费网站观看| 在线观看免费午夜福利视频| 国产日韩欧美视频二区| 人成视频在线观看免费观看| 女人爽到高潮嗷嗷叫在线视频| 亚洲精华国产精华液的使用体验| 一边亲一边摸免费视频| 日韩制服丝袜自拍偷拍| 精品国产一区二区三区久久久樱花| 亚洲av福利一区| a级毛片黄视频| 男女午夜视频在线观看| 午夜免费观看性视频| 亚洲精品久久久久久婷婷小说| 伦理电影免费视频| 丝袜脚勾引网站| 欧美成人精品欧美一级黄| 岛国毛片在线播放| 一区在线观看完整版| 啦啦啦在线免费观看视频4| 在线亚洲精品国产二区图片欧美| 人体艺术视频欧美日本| 国产极品粉嫩免费观看在线| av不卡在线播放| 亚洲欧美色中文字幕在线| 最近最新中文字幕大全免费视频 | 日本午夜av视频| av在线播放精品| 免费黄色在线免费观看| 国产欧美亚洲国产| 最近最新中文字幕免费大全7| 中文乱码字字幕精品一区二区三区| 在线观看免费视频网站a站| 国产黄频视频在线观看| 国产淫语在线视频| 老鸭窝网址在线观看| 亚洲国产中文字幕在线视频| 男人舔女人的私密视频| 亚洲熟女精品中文字幕| 欧美国产精品va在线观看不卡| 精品一品国产午夜福利视频| 亚洲色图 男人天堂 中文字幕| 高清欧美精品videossex| kizo精华| 18禁裸乳无遮挡动漫免费视频| a 毛片基地| 久久av网站| 成人国产麻豆网| 2021少妇久久久久久久久久久| 考比视频在线观看| 日本wwww免费看| 日韩 亚洲 欧美在线| 亚洲成av片中文字幕在线观看| 免费不卡黄色视频| 中文字幕精品免费在线观看视频| 久热这里只有精品99| 一级毛片黄色毛片免费观看视频| 99九九在线精品视频| 国产淫语在线视频| 免费不卡黄色视频| 妹子高潮喷水视频| 满18在线观看网站| 两个人看的免费小视频| 久久精品国产a三级三级三级| 人成视频在线观看免费观看| 欧美日韩亚洲国产一区二区在线观看 | 日日啪夜夜爽| 丁香六月天网| 一区二区日韩欧美中文字幕| 亚洲av电影在线进入| 久久99精品国语久久久| 女人精品久久久久毛片| 国产成人啪精品午夜网站| 丝袜美足系列| 国产精品一国产av| 日本猛色少妇xxxxx猛交久久| 亚洲精品国产一区二区精华液| 精品少妇久久久久久888优播| 一级黄片播放器| 夫妻性生交免费视频一级片| 国产高清国产精品国产三级| www.自偷自拍.com| 午夜福利网站1000一区二区三区| 丁香六月欧美| 久久精品aⅴ一区二区三区四区| 69精品国产乱码久久久| 老司机亚洲免费影院| 亚洲精品美女久久久久99蜜臀 | 99热国产这里只有精品6| 欧美日韩精品网址| 亚洲人成77777在线视频| 亚洲欧洲精品一区二区精品久久久 | 中文字幕av电影在线播放| 国产精品久久久av美女十八| av在线观看视频网站免费| 婷婷色综合www| 亚洲久久久国产精品| 自拍欧美九色日韩亚洲蝌蚪91| 久久久久久免费高清国产稀缺| 亚洲人成77777在线视频| 国产一区亚洲一区在线观看| 国产视频首页在线观看| 国产亚洲av片在线观看秒播厂| 国产精品偷伦视频观看了| 一本大道久久a久久精品| 精品少妇黑人巨大在线播放| 自拍欧美九色日韩亚洲蝌蚪91| 欧美 日韩 精品 国产| 纵有疾风起免费观看全集完整版| 黑人猛操日本美女一级片| 精品国产一区二区三区四区第35| bbb黄色大片| 亚洲精品第二区| 香蕉丝袜av| 亚洲av电影在线观看一区二区三区| 亚洲精品国产av蜜桃| 国产男女超爽视频在线观看| 一级片'在线观看视频| 久久 成人 亚洲| 建设人人有责人人尽责人人享有的| 美女脱内裤让男人舔精品视频| 午夜激情久久久久久久| www.自偷自拍.com| 婷婷色麻豆天堂久久| 欧美日韩国产mv在线观看视频| 亚洲精品美女久久久久99蜜臀 | 十八禁人妻一区二区| 久久精品国产综合久久久| 久久狼人影院| 欧美日韩福利视频一区二区| 建设人人有责人人尽责人人享有的| 成年美女黄网站色视频大全免费| 热re99久久国产66热| 在线亚洲精品国产二区图片欧美| 美女大奶头黄色视频| 日本猛色少妇xxxxx猛交久久| 丝袜在线中文字幕| 丝袜人妻中文字幕| 亚洲精品成人av观看孕妇| 黑丝袜美女国产一区| 精品一区二区三区av网在线观看 | www日本在线高清视频| 啦啦啦 在线观看视频| 各种免费的搞黄视频| 如何舔出高潮| 日韩av在线免费看完整版不卡| 午夜91福利影院| 久久久久精品人妻al黑| 亚洲欧洲精品一区二区精品久久久 | 精品一区在线观看国产| 亚洲成人免费av在线播放| 日本色播在线视频| 黄网站色视频无遮挡免费观看| 久久人妻熟女aⅴ| 免费观看性生交大片5| 国产又色又爽无遮挡免| 777米奇影视久久| av线在线观看网站| 美女主播在线视频| 国产精品国产三级国产专区5o| 嫩草影视91久久| 国产一区二区激情短视频 | 性少妇av在线| tube8黄色片| 熟女少妇亚洲综合色aaa.| 国产午夜精品一二区理论片| 人人妻人人添人人爽欧美一区卜| www.av在线官网国产| 午夜福利视频精品| 在线天堂最新版资源| 欧美 亚洲 国产 日韩一| 亚洲av欧美aⅴ国产| 亚洲欧美一区二区三区国产| 制服诱惑二区| 久久精品国产a三级三级三级| 十分钟在线观看高清视频www| 国产免费又黄又爽又色| 国产片内射在线| 久久亚洲国产成人精品v| 久久精品久久精品一区二区三区| 大陆偷拍与自拍| 国产免费视频播放在线视频| 国产av码专区亚洲av| 男人爽女人下面视频在线观看| 男的添女的下面高潮视频| 亚洲婷婷狠狠爱综合网| 精品人妻熟女毛片av久久网站| 免费久久久久久久精品成人欧美视频| 建设人人有责人人尽责人人享有的| 国产免费福利视频在线观看| 热re99久久国产66热| 亚洲,欧美精品.| 精品午夜福利在线看| 亚洲图色成人| 日本av免费视频播放| 午夜激情av网站| 国产伦人伦偷精品视频| 精品少妇内射三级| 捣出白浆h1v1| 在线观看国产h片| 国产精品香港三级国产av潘金莲 | 一本色道久久久久久精品综合| 国产精品蜜桃在线观看| av在线观看视频网站免费| 97人妻天天添夜夜摸| 丝袜美腿诱惑在线| av卡一久久| 黄色一级大片看看| 日日爽夜夜爽网站| av.在线天堂| 一二三四在线观看免费中文在| 国产精品无大码| 老熟女久久久| www.精华液| 欧美人与善性xxx| 亚洲国产欧美网| 午夜福利网站1000一区二区三区| 黄网站色视频无遮挡免费观看| 国产极品粉嫩免费观看在线| 欧美在线黄色| 丰满迷人的少妇在线观看| 看十八女毛片水多多多| 在线 av 中文字幕| 国产精品无大码| 女人高潮潮喷娇喘18禁视频| 亚洲欧洲日产国产| 国产一区二区三区综合在线观看| 99九九在线精品视频| 欧美精品一区二区大全| 蜜桃国产av成人99| 日韩一卡2卡3卡4卡2021年| 日韩精品免费视频一区二区三区| 色视频在线一区二区三区| 亚洲av日韩精品久久久久久密 | 黄网站色视频无遮挡免费观看| 日韩精品有码人妻一区| 搡老岳熟女国产| 制服人妻中文乱码| 看十八女毛片水多多多| 久久久久精品性色| 欧美精品一区二区免费开放| 韩国精品一区二区三区| 午夜久久久在线观看| 国产欧美亚洲国产| 五月天丁香电影| 男人操女人黄网站| 国产亚洲精品第一综合不卡| 亚洲欧洲日产国产| 19禁男女啪啪无遮挡网站| 美国免费a级毛片| 亚洲激情五月婷婷啪啪| 自线自在国产av| 大码成人一级视频| 中文字幕人妻丝袜一区二区 | 午夜日本视频在线| 亚洲av成人不卡在线观看播放网 | 日本一区二区免费在线视频| 999精品在线视频| 亚洲精品乱久久久久久| 黄片小视频在线播放| 久久韩国三级中文字幕| 国产精品久久久久久人妻精品电影 | 少妇精品久久久久久久| 精品第一国产精品| 国产高清不卡午夜福利| 一级黄片播放器| 国产精品秋霞免费鲁丝片| 秋霞伦理黄片| 国产免费视频播放在线视频| 成人毛片60女人毛片免费| 欧美 日韩 精品 国产| xxx大片免费视频| 91精品三级在线观看| 国产精品久久久久成人av| 免费高清在线观看视频在线观看| 极品少妇高潮喷水抽搐| av一本久久久久| 女性生殖器流出的白浆| 国产精品一区二区精品视频观看| 人人妻人人澡人人看| 久热这里只有精品99| 男男h啪啪无遮挡| 国产女主播在线喷水免费视频网站| 如何舔出高潮| 国产亚洲欧美精品永久| 老汉色av国产亚洲站长工具| 在线 av 中文字幕| 看非洲黑人一级黄片| 成年动漫av网址| 亚洲色图综合在线观看| 纵有疾风起免费观看全集完整版| 如何舔出高潮| 熟女少妇亚洲综合色aaa.| 岛国毛片在线播放| 国产精品一区二区精品视频观看| 一区二区av电影网| av视频免费观看在线观看| 午夜日韩欧美国产| 国产精品国产三级专区第一集| videos熟女内射| 91老司机精品| 天堂中文最新版在线下载| 国产成人欧美| 国产男人的电影天堂91| 黄色一级大片看看| 狂野欧美激情性bbbbbb| 99热网站在线观看| 高清视频免费观看一区二区| 超碰成人久久| av福利片在线| 美女中出高潮动态图| 色视频在线一区二区三区| 满18在线观看网站| 天天躁夜夜躁狠狠躁躁| 夫妻午夜视频| 人体艺术视频欧美日本| 性色av一级| 色综合欧美亚洲国产小说| 啦啦啦 在线观看视频| 午夜激情av网站| 亚洲国产av影院在线观看| h视频一区二区三区| 好男人视频免费观看在线| 汤姆久久久久久久影院中文字幕| 丁香六月天网| 国产av国产精品国产| 亚洲av中文av极速乱| 天天躁狠狠躁夜夜躁狠狠躁| 校园人妻丝袜中文字幕| a级毛片在线看网站| 亚洲成人免费av在线播放| 久久精品久久精品一区二区三区| 日韩av不卡免费在线播放| 午夜福利视频在线观看免费| 色视频在线一区二区三区| 欧美精品人与动牲交sv欧美| 日日啪夜夜爽| 最近中文字幕高清免费大全6| a级毛片在线看网站| 丝袜喷水一区| 久久 成人 亚洲| 亚洲精品,欧美精品| 国产精品香港三级国产av潘金莲 |