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

    一種用于間斷裝配的流場(chǎng)結(jié)構(gòu)辨識(shí)算法

    2019-08-22 01:07:22陳澤棟鄒東陽(yáng)
    關(guān)鍵詞:方法

    劉 君,陳澤棟,陳 潔,鄒東陽(yáng)

    (1.大連理工大學(xué) 航空航天學(xué)院,大連 116024;2.大連理工大學(xué) 工程力學(xué)系,大連 116024)

    0 引 言

    流體微團(tuán)加速到超聲速以后無(wú)法感知到達(dá)之前遇到的空間物體,只能在物體前通過(guò)急劇壓縮的激波才能把速度降下來(lái)。因此,激波是可壓縮流動(dòng)的基本特征。早期計(jì)算流體力學(xué)(CFD)的算法無(wú)法模擬這一現(xiàn)象,直到1954年Lax提出處理包含有間斷流動(dòng)問(wèn)題的弱解理論后[1-2],國(guó)內(nèi)外眾多科學(xué)家開(kāi)始研究在計(jì)算過(guò)程中能自動(dòng)分辨出激波的算法,提出TVD、NND、ENO、WENO等具有里程碑意義的激波捕捉格式[3],解決了航天航空等領(lǐng)域很多工程模擬問(wèn)題,極大提高了CFD地位,使其發(fā)展為一個(gè)重要的學(xué)科分支。盡管建立這些格式時(shí)考慮了激波間斷,但是本質(zhì)上還是1950年von Neumann提出的把間斷描述成流動(dòng)參數(shù)連續(xù)變化過(guò)渡區(qū)的思路[4],構(gòu)造出來(lái)的格式存在以下難以克服的理論困境。

    數(shù)值計(jì)算的離散格式是對(duì)偏微分方程的近似,按照CFD收斂性理論,數(shù)值解隨著網(wǎng)格尺度減小趨向理論解,CFD追求的極限是方程的理論解,原則上說(shuō)對(duì)于方程不能描述的流動(dòng)現(xiàn)象,數(shù)值解和實(shí)驗(yàn)值不具有可比性。海平面環(huán)境下馬赫數(shù)小于10的激波厚度大約1×10-8m,流體在這么小的空間內(nèi)受到劇烈壓縮,暫且不說(shuō)準(zhǔn)確模擬所需要的網(wǎng)格尺度及其帶來(lái)的計(jì)算量,從連續(xù)介質(zhì)假設(shè)出發(fā)建立的N-S方程是否適用于激波模擬也存在爭(zhēng)議,因?yàn)镾toke假設(shè)(3λ+2μ=0)不再成立[5],計(jì)算時(shí)需要對(duì)體積膨脹黏性系數(shù)λ進(jìn)行修正才能得到與實(shí)驗(yàn)相符的激波內(nèi)部溫度[6]。從事稀薄氣體動(dòng)力學(xué)研究的學(xué)者認(rèn)為激波內(nèi)部結(jié)構(gòu)涉及到熱力學(xué)非平衡,連續(xù)介質(zhì)假設(shè)本身就不成立,“從理論上說(shuō),Navier-Stokes方程不能用來(lái)描述激波結(jié)構(gòu)的[7]”。無(wú)黏流體力學(xué)把激波表示成間斷,兩側(cè)參數(shù)和運(yùn)動(dòng)速度滿足蘭金-許貢紐(Rankine-Hugoniot,R-H)關(guān)系式,這是迄今為止幾乎很少有爭(zhēng)議的理論。所謂間斷在數(shù)學(xué)上就是在同一空間位置對(duì)應(yīng)至少兩組(激波前后)流體變量,在多維情況下可能更多,例如二維馬赫反射三岔點(diǎn)有四組流體變量。很多基于弱解理論構(gòu)造的激波捕捉格式考慮了間斷,但是用參數(shù)連續(xù)變化的思路描述激波注定不可能計(jì)算出真正的數(shù)學(xué)間斷,捕捉到的激波主要是現(xiàn)象模擬。

    按照上述觀點(diǎn),常規(guī)N-S方程及其理論解不能用來(lái)描述激波內(nèi)部物理現(xiàn)象,那么基于離散格式的數(shù)值解至少在過(guò)渡區(qū)內(nèi)沒(méi)有物理意義,采用Euler方程產(chǎn)生的后果是無(wú)法評(píng)價(jià)過(guò)渡區(qū)內(nèi)數(shù)值解的精度(激波前后參數(shù)相差很大),激波捕捉格式對(duì)于激波相交等現(xiàn)象更是無(wú)能為力。

    20世紀(jì)60年代Moretti等人將激波裝配方法和Lax-Wendroff格式結(jié)合,在CDC 600型計(jì)算機(jī)上進(jìn)行,總耗時(shí)為360 s,給出了超聲速鈍頭體問(wèn)題的數(shù)值解,2002年Moretti回顧從業(yè)36年來(lái)CFD領(lǐng)域模擬激波的進(jìn)展,針對(duì)超聲速鈍頭體繞流比較了捕捉法和裝配法計(jì)算結(jié)果,認(rèn)為準(zhǔn)確模擬激波還需要發(fā)展裝配法[8-9]。對(duì)于50年前提出的裝配法,文獻(xiàn)[10]的評(píng)價(jià)具有代表性——“激波裝配法具有計(jì)算精度高、物理概念清晰等優(yōu)點(diǎn),但是計(jì)算過(guò)程復(fù)雜、計(jì)算量大,并且必須事先知道激波的大概位置。一般來(lái)說(shuō),激波裝配法僅適用于那些激波比較簡(jiǎn)單而事先可以估計(jì)激波形狀和位置的定常流動(dòng)問(wèn)題?!?/p>

    造成計(jì)算復(fù)雜的主要原因是求解過(guò)程中激波運(yùn)動(dòng)和變形。裝配法減少了網(wǎng)格量,造成計(jì)算量大的原因是激波及其相交點(diǎn)等空間網(wǎng)格點(diǎn)需要判別、標(biāo)記或特殊處理,增加代碼編寫(xiě)復(fù)雜性和并行計(jì)算難度。隨著CFD網(wǎng)格技術(shù)發(fā)展,這些難點(diǎn)也得到緩解。國(guó)外文獻(xiàn)把最早Moretti提出的稱為邊界激波裝配法 ,僅能處理形狀較為規(guī)則和位移不大的脫體激波。后來(lái)Richtmyer等人[11]提出了浮動(dòng)激波裝配方法,在靜止網(wǎng)格基礎(chǔ)上嵌入R-H關(guān)系式,允許激波在背景網(wǎng)格上滑動(dòng),解決了內(nèi)部激波裝配難題。近年來(lái)Zhong等人將這種方法和高精度有限差分格式相結(jié)合有效提高內(nèi)部光滑流場(chǎng)的計(jì)算精度,成功應(yīng)用于許多高超聲速流動(dòng)的求解中[12-13]?;诮Y(jié)構(gòu)網(wǎng)格的邊界激波裝配法和浮動(dòng)激波裝配法在模擬激波相交和馬赫反射等復(fù)雜波系結(jié)構(gòu)時(shí)面臨一定困難。Paciorri和Bonfigolioli等人[14-16]采用非結(jié)構(gòu)網(wǎng)格技術(shù)和局部網(wǎng)格重構(gòu)提出了激波陣面追蹤法,記錄激波在背景網(wǎng)格上的位置,對(duì)激波附近的背景網(wǎng)格進(jìn)行挖洞,然后再利用Delanay方法在洞內(nèi)重構(gòu)以激波陣面分界的新網(wǎng)格,時(shí)間推進(jìn)過(guò)程中根據(jù)R-H關(guān)系式得到新的激波位置,重新判別、挖洞和填補(bǔ)背景網(wǎng)格。這種算法不需要考慮網(wǎng)格節(jié)點(diǎn)之間的連接關(guān)系,也就減少了對(duì)于特殊點(diǎn)進(jìn)行處理的麻煩,增加了代碼的普適性,但是在激波邊界節(jié)點(diǎn)對(duì)應(yīng)的網(wǎng)格單元只能采用一階插值,像常用的多塊重疊網(wǎng)格技術(shù)一樣存在插值引起的精度降低問(wèn)題。其次,不允許激波在一個(gè)時(shí)間步內(nèi)跨越多個(gè)網(wǎng)格節(jié)點(diǎn),通常選擇足夠小的時(shí)間步長(zhǎng),因而對(duì)計(jì)算效率有所影響。上述這3類裝配法時(shí)間推進(jìn)過(guò)程中需要判別激波位置,增加了計(jì)算量,而且激波前后網(wǎng)格數(shù)目是變化的,不利于并行計(jì)算。2015年劉君等人[17]提出了基于非結(jié)構(gòu)動(dòng)網(wǎng)格的激波(間斷)裝配方法,從能夠描述網(wǎng)格變形的任意拉格朗日-歐拉坐標(biāo)系(Arbitrary Lagrangian-Eulerian,ALE)下的控制方程出發(fā),得到激波位置及其兩側(cè)的參數(shù)后,裝配的激波就變成了一個(gè)移動(dòng)邊界,對(duì)于存在內(nèi)部激波的計(jì)算域則變成類似于拼接網(wǎng)格的多個(gè)分區(qū),激波運(yùn)動(dòng)由上下游參數(shù)驅(qū)動(dòng),采用成熟的動(dòng)網(wǎng)格技術(shù)[18]進(jìn)行移動(dòng)和變形。與Paciorri及Bonfigolioli的裝配法相比,在保留非結(jié)構(gòu)網(wǎng)格的優(yōu)點(diǎn)外,新的裝配法不需要在每一步對(duì)計(jì)算網(wǎng)格進(jìn)行挖洞重構(gòu),避免了插值誤差,易于并行計(jì)算,并且時(shí)間步長(zhǎng)不受限制。采用這種算法模擬了內(nèi)部存在正規(guī)反射和馬赫反射等復(fù)雜激波結(jié)構(gòu)的算例[19-21],從效果看,較好地解決了傳統(tǒng)裝配法遇到的“計(jì)算過(guò)程復(fù)雜、計(jì)算量大”問(wèn)題。

    為了突破“激波裝配法僅適用于那些激波比較簡(jiǎn)單而事先可以估計(jì)激波形狀和位置的定常流動(dòng)問(wèn)題”的限制,劉君等人提出了在采用捕捉法計(jì)算的流場(chǎng)基礎(chǔ)上辨識(shí)并裝配主要激波結(jié)構(gòu)、計(jì)算過(guò)程中自動(dòng)裝配內(nèi)部細(xì)致結(jié)構(gòu)和新生激波的解決思路。在技術(shù)途徑上提出了新的數(shù)據(jù)結(jié)構(gòu),在激波邊界通量計(jì)算時(shí),如果采用Roe、van Leer、AUSM等格式就是捕捉法,如果采用R-H關(guān)系式就是裝配法,實(shí)現(xiàn)兩種方法之間任意轉(zhuǎn)換。按照上述思路,為了解決裝配法的應(yīng)用難題,亟需發(fā)展激波、接觸間斷等流場(chǎng)特征結(jié)構(gòu)的辨識(shí)算法。

    目前,常用的激波辨識(shí)算法主要有三種[22]:第一種根據(jù)最大變量梯度確定激波面,該方法實(shí)現(xiàn)簡(jiǎn)單,但需要設(shè)置合適的過(guò)濾器才能消除錯(cuò)誤的結(jié)果;第二種為法向馬赫數(shù)法[24],該方法認(rèn)為在激波面上流動(dòng)的法向馬赫數(shù)等于1,該方法與第一種方法相比適用性更強(qiáng),但也需要設(shè)置適當(dāng)?shù)倪^(guò)濾器才能獲得良好結(jié)果,且難以用于探測(cè)接觸間斷;第三種為基于特征線理論的方法[25-26],該方法與前兩種方法相比更精確、更有效,但實(shí)施過(guò)程復(fù)雜且計(jì)算量大。

    綜上所述,這些激波辨識(shí)算法難以在復(fù)雜三維問(wèn)題中廣泛應(yīng)用,且上述算法辨識(shí)出的激波通常為具有一定寬度的激波帶,難以直接用于裝配法。因此,需要發(fā)展新的辨識(shí)算法來(lái)推動(dòng)激波裝配法實(shí)用化進(jìn)程。

    本文提出了一種光源射線法,該方法可將捕捉法中探測(cè)到的具有一定寬度的激波帶擬合成一系列離散線段(二維)或三角形面片(三維),利用曲線/曲面擬合算法進(jìn)一步將其擬合成連續(xù)的曲線/曲面。該曲線/曲面可作為激波裝配法中的激波面,并且可根據(jù)計(jì)算需要對(duì)其進(jìn)行任意形式的網(wǎng)格劃分。為便于表述,先用二維模型的實(shí)際操作過(guò)程進(jìn)行原理說(shuō)明,然后簡(jiǎn)單介紹三維推廣過(guò)程,并展示實(shí)際工程問(wèn)題的辨識(shí)結(jié)果。

    1 流場(chǎng)結(jié)構(gòu)辨識(shí)算法

    1.1 二維辨識(shí)算法流程

    (a)捕捉法流場(chǎng) (b)激波單元 (c)激波曲線

    (d)分片擬合

    (e)確定特征點(diǎn)圖1 光源射線法示意圖Fig.1 Sketch of source-rays method

    1.2 脫體激波辨識(shí)

    針對(duì)二維超聲速鈍頭體繞流問(wèn)題,采用1.1節(jié)所述方法進(jìn)行激波面辨識(shí)。算例中,空氣來(lái)流馬赫數(shù)Ma=5.0、無(wú)量綱密度ρ=1、無(wú)量綱溫度T=1,來(lái)流攻角為0°。鈍頭體中心位置坐標(biāo)為(0,0),鈍頭體半徑為0.0254。網(wǎng)格為160×100的四邊形網(wǎng)格,其中沿鈍頭體法向的節(jié)點(diǎn)數(shù)為100。計(jì)算網(wǎng)格和捕捉法流場(chǎng)壓力云圖分別如圖2(a)和圖2(b)所示。

    (a)網(wǎng)格 (b)壓力云圖 (c)激波點(diǎn)集 (d)激波面圖2 基于壓力梯度的激波辨識(shí)結(jié)果Fig.2 Shock wave surface extraction based on pressure gradient

    考慮到很多高精度格式的限制器選擇壓力梯度來(lái)進(jìn)行調(diào)節(jié),在此也將其作為激波特征參數(shù)來(lái)探測(cè)脫體激波。對(duì)于捕捉法得到的壓力場(chǎng),用格林公式計(jì)算出每個(gè)網(wǎng)格單元的壓力梯度矢量。作為間斷面,激波厚度理論上處于分子自由程大小這一數(shù)量級(jí),其兩側(cè)壓力的梯度值為無(wú)窮大。但是,在捕捉法數(shù)值計(jì)算中,激波表現(xiàn)為帶狀分布的過(guò)渡區(qū)且具有較大壓力梯度。對(duì)于不同的問(wèn)題以及不同的計(jì)算網(wǎng)格,激波帶上的壓力梯度值都有所不同。因此,難以通過(guò)一個(gè)固定的壓力梯度值將激波帶提取出來(lái)。本文為了準(zhǔn)確獲得激波點(diǎn)集的分布,首先計(jì)算出流場(chǎng)中的最大壓力梯度值|p|max,然后針對(duì)不同問(wèn)題設(shè)置相應(yīng)的閾值系數(shù)k,壓力梯度閾值|p|thre=k|p|max,所有滿足條件|p|≥|p|thre的單元格心點(diǎn)構(gòu)成激波點(diǎn)集T。Case1取閾值系數(shù)k=0.1,獲得離散點(diǎn)集,如圖2(c)所示。針對(duì)上述激波離散點(diǎn)集,采用光源射線法進(jìn)行曲線擬合。光源S位置固定在(x,y)=(0.1,0)處,相鄰兩條射線間的夾角取常值φi=4.5°。辨識(shí)出的結(jié)果如圖2(d)所示,從圖中可以看出辨識(shí)出的激波曲線與捕捉法流場(chǎng)中的激波點(diǎn)集T吻合度很高,能準(zhǔn)確描述激波形狀,對(duì)激波曲線重新進(jìn)行網(wǎng)格離散即可作為激波邊界直接用于激波裝配法。

    離散點(diǎn)集的分布取決于所選擇參數(shù)和閾值,為了比較不同壓力梯度閾值下所辨識(shí)出的激波曲線,在捕捉法流場(chǎng)的基礎(chǔ)上另外補(bǔ)充兩個(gè)辨識(shí)狀態(tài):Case2取k=0.2,Case3取k=0.3。得到的點(diǎn)集及其最終辨識(shí)出來(lái)的脫體激波曲線如圖3所示。從3個(gè)狀態(tài)的辨識(shí)結(jié)果可以看出,不同的壓力梯度條件下所獲取的離散點(diǎn)集有所不同,導(dǎo)致最終辨識(shí)出的激波曲線長(zhǎng)度有所區(qū)別,但三條激波曲線都能與捕捉法流場(chǎng)中激波的位置和形狀保持一致,將其作為激波裝配法中的激波邊界可快速得到收斂解。得到初始激波位置以后的裝配法計(jì)算方法和過(guò)程參見(jiàn)文獻(xiàn)[17,19-21],這里不再贅述。

    (a)Case1 (b)Case2 (c)Case3

    圖3 不同壓力梯度閾值下的脫體激波辨識(shí)結(jié)果比較
    Fig.3 Comparison of shock surface extraction results under different pressure gradient thresholds

    1.3 內(nèi)部激波及接觸間斷辨識(shí)

    對(duì)波系更為復(fù)雜的馬赫反射問(wèn)題進(jìn)行數(shù)值模擬,計(jì)算區(qū)域及邊界條件如圖4所示。邊長(zhǎng)為L(zhǎng)AF=0.4、LAB=0.1、LDE=0.312、LEF=0.85,斜楔BC傾角為14°;邊界AF為超聲速入口、DE為超聲速出口,其余為滑移壁。計(jì)算網(wǎng)格為300×218的四邊形網(wǎng)格,其中入口和出口邊界上的網(wǎng)格節(jié)點(diǎn)數(shù)為300。來(lái)流馬赫數(shù)Ma=1.9,攻角為0°。

    圖4 計(jì)算域Fig.4 Calculation domain

    由于斜楔BC的偏轉(zhuǎn)作用,自由來(lái)流在點(diǎn)B處產(chǎn)生一道入射激波(Shock1);激波達(dá)到壁面FE處時(shí)發(fā)生馬赫反射,產(chǎn)生一道反射激波(Shock2)和一條馬赫桿(Shock3);反射激波在壁面CD處再次發(fā)生反射,產(chǎn)生第二道反射激波(Shock4)。流場(chǎng)穩(wěn)定后的密度云圖以及波系結(jié)構(gòu)如圖5所示。采用激波及接觸間斷探測(cè)方法可從流場(chǎng)中獲得主要的四條激波點(diǎn)集和一條接觸間斷點(diǎn)集分布。為了檢驗(yàn)不同的探測(cè)方法對(duì)最終結(jié)果的影響,激波點(diǎn)集采用密度梯度作為參考量進(jìn)行探測(cè),接觸間斷點(diǎn)集由密度梯度和壓力梯度共同決定。同樣,先計(jì)算出全場(chǎng)最大密度梯度值|ρ|max和最大壓力梯度|p|max。激波探測(cè)的密度梯度閾值所有滿足的單元節(jié)點(diǎn)構(gòu)成激波點(diǎn)集,這里取接觸間斷單元由密度梯度和壓力梯度共同決定,滿足密度梯度較大且壓力梯度較小,即其中所有滿足該條件的單元節(jié)點(diǎn)構(gòu)成曲線擬合所需的接觸間斷點(diǎn)集。這里取其中在不同密度梯度閾值下辨識(shí)出的接觸間斷長(zhǎng)度有一定區(qū)別,但位置和形狀基本一致。激波和接觸間斷離散點(diǎn)集分布如圖5所示。

    獲得離散點(diǎn)集后,通過(guò)在區(qū)域內(nèi)合理布置光源可擬合出每條激波和接觸間斷曲線。首先在Ⅱ區(qū)設(shè)置一個(gè)光源,光源發(fā)出等間距射線,采用上文光源射線法可辨識(shí)出表征激波1和激波2的曲線,如圖6(a)所示。同理,在Ⅲ區(qū)布置光源可辨識(shí)出激波4和接觸間斷曲線,如圖6(b)所示。在I區(qū)設(shè)置光源可辨識(shí)出馬赫桿曲線,如圖6(c)所示。所有曲線的最終辨識(shí)結(jié)果與流場(chǎng)密度云圖的比較如圖7所示??梢钥闯觯孀R(shí)出的曲線與采用捕捉法獲得的波系結(jié)構(gòu)高度吻合,直接用于裝配法有助于快速得到收斂解。得到初始激波位置以后的裝配法計(jì)算方法和過(guò)程參見(jiàn)文獻(xiàn)[17,19-21],后文也會(huì)有簡(jiǎn)單介紹,這里不再贅述。

    圖5 捕捉法密度云圖及間斷點(diǎn)集分布Fig.5 Density contours and discontinuity point set distribution under shock-capturing schemes

    (a)激波1-2 (b)激波4和接觸間斷 (c)激波3

    圖6 激波和接觸間斷辨識(shí)過(guò)程
    Fig.6 Shock wave and contact discontinuity surface extraction process

    圖7 激波及接觸間斷辨識(shí)結(jié)果與流場(chǎng)比較Fig.7 Comparison between feature surface extraction result and flow field structure

    值得注意的是,這里光源的位置分布并不唯一,可采用不同的分布方式達(dá)到相似的效果,比如可先在Ⅰ區(qū)設(shè)置光源同時(shí)擬合出激波1和激波3曲線,再在Ⅲ區(qū)設(shè)置光源得到激波2、激波4和接觸間斷曲線。

    2 三維實(shí)現(xiàn)

    2.1 三維辨識(shí)算法流程

    本文通過(guò)引入單位球可確定各射線之間的連接關(guān)系,如圖9所示。對(duì)半徑為1的球面進(jìn)行三角形離散,獲得一系列三角形面片,各頂點(diǎn)連接關(guān)系確定,且三角形面片滿足連續(xù)且不相交條件。

    圖8 射線方向示意圖

    Fig.8 Ray direction

    圖9 單位球

    Fig.9 Unit sphere

    將單位球的球心作為光源S,球心與單位球上三角形面片的三個(gè)頂點(diǎn)ai、aj、ak的連線li、lj、lk構(gòu)成一個(gè)子區(qū)域Ap,如圖10(a)所示。若Ap內(nèi)離散點(diǎn)的個(gè)數(shù)不小于3,則對(duì)Ap內(nèi)所有離散點(diǎn)進(jìn)行平面擬合,擬合出的平面與三條射線的交點(diǎn)分別為bi,p、bj,p和bk,p;若該子區(qū)域內(nèi)離散點(diǎn)的個(gè)數(shù)小于3,則跳過(guò)該區(qū)域。所有區(qū)域擬合完成后,每條射線li與不同子區(qū)域內(nèi)的擬合平面形成若干交點(diǎn)。將各射線li上的交點(diǎn)取平均值獲得特征點(diǎn)bi,則各子區(qū)域Ap的三條射線邊界上的特征點(diǎn)bi、bj和bk構(gòu)成目標(biāo)三角形面片,如圖10(b)所示。最終,整個(gè)離散點(diǎn)集T被擬合成一系列連續(xù)的三角形面片,三角形面片頂點(diǎn)的連接關(guān)系由單位球確定??蓪?duì)上述三角形面片構(gòu)成的空間曲面重新進(jìn)行網(wǎng)格劃分,獲得可用于裝配法的激波面網(wǎng)格。

    (a)分區(qū)擬合 (b)三角面片圖10 子區(qū)域內(nèi)三角形面片擬合過(guò)程Fig.10 Triangular fitting process in sub-region

    2.2 三維測(cè)試算例

    通過(guò)以下兩個(gè)算例驗(yàn)證本文辨識(shí)算法在三維問(wèn)題中的有效性。

    在橢球面附近隨機(jī)生成100 000個(gè)點(diǎn)作為離散點(diǎn)集,隨機(jī)點(diǎn)坐標(biāo)B(x′,y′,z′)滿足:

    其中,

    式中,i、j、k、l、m為區(qū)間[0,1]內(nèi)的隨機(jī)數(shù)。

    獲得的離散點(diǎn)集如圖11所示。單位球共有736個(gè)節(jié)點(diǎn)、1468個(gè)三角形單元(如圖9),球心位置坐標(biāo)為(0,0,0)。單位球位置及橢球面辨識(shí)結(jié)果如圖12所示。從結(jié)果可以看出,辨識(shí)出的曲面光滑完整且形狀與離散點(diǎn)分布形狀一致。

    圖11 橢球隨機(jī)點(diǎn)集圖12 單位球及擬合結(jié)果Fig.11 Random point set of ellipsoidFig.12 Unit sphere and the extracted surface of ellipsoid

    圖13 雙橢球隨機(jī)點(diǎn)集圖14 雙橢球面擬合結(jié)果Fig.13 Random point set of double-ellipsoidFig.14 Extracted surface of double-ellipsoid

    2.3 三維繞流激波辨識(shí)

    對(duì)三維超聲速球頭繞流進(jìn)行模擬,來(lái)流馬赫數(shù)Ma=5.0,無(wú)量綱密度ρ=1.0,無(wú)量綱溫度T=1.0,來(lái)流攻角和側(cè)滑角均為0°。球頭的球心位置坐標(biāo)為(0,0,0),球頭半徑為0.0254。計(jì)算網(wǎng)格共有508 365個(gè)的六面體單元,球面網(wǎng)格如圖15(a)所示。圖15(b)為流場(chǎng)壓力云圖,圖15(c)為根據(jù)流場(chǎng)壓力梯度獲得的激波帶離散點(diǎn)集。算例中,壓力梯度閾值系數(shù)k=0.05離散點(diǎn)為滿足|p|≥|p|thre的所有單元節(jié)點(diǎn),共有23 640個(gè)。單位球共包含1734個(gè)節(jié)點(diǎn),3464個(gè)三角形面片,球心坐標(biāo)為(2,0,0)。最終辨識(shí)出的激波面如圖15(d)所示,除了邊緣附近外的其他區(qū)域辨識(shí)結(jié)果與離散點(diǎn)集吻合較好,將邊緣位置稍作處理即可作為初始激波面用于激波裝配法。對(duì)于不同的初始位置的激波,最終收斂解相同,故本例不再展示裝配法的計(jì)算過(guò)程,計(jì)算方法參見(jiàn)文獻(xiàn)[17,19-21]。

    (a)球頭表面網(wǎng)格 (b)捕捉法壓力云圖

    (c)激波點(diǎn)集 (d)激波面圖15 三維超聲速球頭繞流激波辨識(shí)過(guò)程Fig.15 Shock wave surface extraction process of 3D blunt body flow

    在球頭繞流算例的基礎(chǔ)上,對(duì)波系較為復(fù)雜的球柱錐組合體外形進(jìn)行模擬,檢驗(yàn)特征面提取方法處理三維復(fù)雜問(wèn)題的能力。計(jì)算模型如圖16(a)所示,來(lái)流馬赫數(shù)Ma=4.04,來(lái)流攻角為20°。波系結(jié)構(gòu)如圖16(b)所示,在最外層形成一道弓形激波,在錐裙與柱體的結(jié)合處形成一個(gè)內(nèi)嵌激波。探測(cè)出的激波點(diǎn)集和提取出的激波面分別如圖16(c)和圖16(d)所示。

    (a)計(jì)算模型 (b)捕捉法壓力云圖

    (c)激波點(diǎn)集 (d)激波面圖16 球柱錐組合體繞流激波辨識(shí)過(guò)程Fig.16 Shock wave surface extraction process of cylinder with a hemispherical nose and a conical flow

    3 裝配法實(shí)例

    利用上述方法獲得初始激波面后,將其作為激波邊界用于激波裝配法,進(jìn)行流動(dòng)數(shù)值模擬,考查所述激波辨識(shí)方法在實(shí)際使用中的效果。

    對(duì)于定常激波裝配法,初始激波位置及形狀會(huì)影響計(jì)算的收斂過(guò)程,但最終的收斂結(jié)果都相同。因此,本節(jié)直接采用文獻(xiàn)[19]算例的部分結(jié)果說(shuō)明激波裝配法與捕捉法的區(qū)別。該算例為具有復(fù)雜波系結(jié)構(gòu)的馬赫反射問(wèn)題,來(lái)流馬赫數(shù)2.0,偏轉(zhuǎn)角δ=14°。計(jì)算區(qū)域?yàn)長(zhǎng)×L的正方形,邊界上網(wǎng)格節(jié)點(diǎn)均勻分布,節(jié)點(diǎn)間距為0.01L,內(nèi)部區(qū)域?yàn)槿切尉W(wǎng)格,單元數(shù)為22 429。流場(chǎng)波系結(jié)構(gòu)及邊界條件如圖17所示,入射激波IS和反射激波RS、馬赫桿MS交于三波點(diǎn)TP,SS為接觸間斷。

    首先采用激波捕捉法獲得初始流場(chǎng),在此基礎(chǔ)上采用本文辨識(shí)算法辨識(shí)出激波及接觸間斷的位置和形狀,將其作為裝配法的初始間斷面,結(jié)合R-H關(guān)系式進(jìn)行裝配法流場(chǎng)計(jì)算。得到初始激波位置以后的計(jì)算過(guò)程參見(jiàn)文獻(xiàn)[19],計(jì)算結(jié)果完全一致。下面給出部分結(jié)果,如圖18所示,其中圖18(a)為采用捕捉法的計(jì)算結(jié)果,圖18(b)為采用裝配法的計(jì)算結(jié)果。從圖中可以看出,捕捉法流場(chǎng)中激波附近出現(xiàn)等值線的扭曲和震蕩,且激波呈帶狀分布,具有很寬的過(guò)渡區(qū)。與此相比,采用裝配法所獲得的流場(chǎng)中各區(qū)域之間界線分明,激波是一個(gè)沒(méi)有寬度的界面,激波兩側(cè)變量值嚴(yán)格滿足R-H關(guān)系式。

    圖17 邊界條件及流場(chǎng)結(jié)構(gòu)示意圖[19]Fig.17 Boundary conditions and sketch of flow field structure[19]

    (a)捕捉法計(jì)算結(jié)果

    (b)裝配法計(jì)算結(jié)果圖18 裝配法與捕捉法計(jì)算結(jié)果對(duì)比(馬赫數(shù)云圖)[19]Fig.18 Comparison between shock-capturing solution and shock-fitting solution[19]

    4 結(jié) 論

    通常對(duì)于具有一定“厚度”的空間點(diǎn)集,采用傳統(tǒng)的擬合方法難以將其擬合成形狀一致的光滑曲面,本文提出的光源射線法很好地解決了此類問(wèn)題。

    對(duì)于從捕捉法流場(chǎng)中獲取的具有一定“厚度”的激波帶和接觸間斷帶,采用該方法進(jìn)行空間曲面擬合后可作為激波裝配法中的初始間斷面。尤其對(duì)于具有復(fù)雜波系結(jié)構(gòu)的流場(chǎng),該方法提供了一種裝配初始流場(chǎng)的途徑。本文詳細(xì)介紹了該算法在二維問(wèn)題中的實(shí)現(xiàn)過(guò)程,并將其成功推廣應(yīng)用于三維問(wèn)題。算例結(jié)果表明,辨識(shí)出的激波曲面整體上與捕捉法流場(chǎng)中激波帶分布相吻合,但是三維問(wèn)題中激波曲面邊緣附近的辨識(shí)效果有待進(jìn)一步提高。對(duì)于復(fù)雜工程問(wèn)題,往往存在較為復(fù)雜的三維波系結(jié)構(gòu),間斷面的提取過(guò)程更為困難和繁瑣,該方法還需進(jìn)一步完善。

    對(duì)于傳統(tǒng)激波裝配法,在初始激波面給定后需要經(jīng)過(guò)若干步迭代才能收斂到準(zhǔn)確的形狀和位置。尤其對(duì)于三維問(wèn)題,收斂過(guò)程更加漫長(zhǎng),甚至出現(xiàn)難以收斂的情況。采用本文方法可提供較為準(zhǔn)確的初始激波,有助于縮短收斂過(guò)程。

    本文所提方法除了用于激波裝配法外,未來(lái)還可應(yīng)用于自適應(yīng)網(wǎng)格算法。根據(jù)不同的流動(dòng)特征可擬合出相應(yīng)的空間曲面,進(jìn)而可以任意調(diào)整該曲面附近區(qū)域的網(wǎng)格分布,達(dá)到網(wǎng)格根據(jù)流動(dòng)特征自動(dòng)調(diào)整的目的。

    猜你喜歡
    方法
    中醫(yī)特有的急救方法
    中老年保健(2021年9期)2021-08-24 03:52:04
    高中數(shù)學(xué)教學(xué)改革的方法
    化學(xué)反應(yīng)多變幻 “虛擬”方法幫大忙
    變快的方法
    兒童繪本(2020年5期)2020-04-07 17:46:30
    學(xué)習(xí)方法
    可能是方法不對(duì)
    用對(duì)方法才能瘦
    Coco薇(2016年2期)2016-03-22 02:42:52
    最有效的簡(jiǎn)單方法
    山東青年(2016年1期)2016-02-28 14:25:23
    四大方法 教你不再“坐以待病”!
    Coco薇(2015年1期)2015-08-13 02:47:34
    賺錢(qián)方法
    欧美大码av| 国产综合懂色| 少妇熟女aⅴ在线视频| 中文字幕av成人在线电影| 国内揄拍国产精品人妻在线| 男女之事视频高清在线观看| 国内精品久久久久久久电影| 午夜视频国产福利| 久久久久亚洲av毛片大全| 男女床上黄色一级片免费看| 国产精华一区二区三区| 久久久精品大字幕| 精品一区二区三区人妻视频| 一进一出好大好爽视频| 男人的好看免费观看在线视频| 一本综合久久免费| 亚洲精品一区av在线观看| 国产亚洲欧美在线一区二区| 国产成人影院久久av| 欧美中文综合在线视频| 免费av不卡在线播放| 久久精品国产自在天天线| 少妇裸体淫交视频免费看高清| 黄色片一级片一级黄色片| 国产探花在线观看一区二区| 国产极品精品免费视频能看的| 他把我摸到了高潮在线观看| 久久性视频一级片| 久久人妻av系列| av天堂中文字幕网| 久久久久久久久久黄片| av女优亚洲男人天堂| 高潮久久久久久久久久久不卡| 午夜两性在线视频| 中文字幕人妻丝袜一区二区| 成人性生交大片免费视频hd| 久久人妻av系列| 九九在线视频观看精品| 亚洲精品色激情综合| 国产成人影院久久av| 国产在视频线在精品| 变态另类成人亚洲欧美熟女| 久久久久久久久久黄片| 亚洲人成电影免费在线| 欧美日本视频| 3wmmmm亚洲av在线观看| 两个人视频免费观看高清| 成年女人毛片免费观看观看9| 一本综合久久免费| 国产午夜精品久久久久久一区二区三区 | 国产极品精品免费视频能看的| 国产私拍福利视频在线观看| av女优亚洲男人天堂| 精品国内亚洲2022精品成人| 中文字幕av成人在线电影| 国产亚洲欧美在线一区二区| 午夜激情欧美在线| 国产免费男女视频| 国产真人三级小视频在线观看| 日本a在线网址| 中出人妻视频一区二区| 欧美日韩综合久久久久久 | 欧美最新免费一区二区三区 | 亚洲精品国产精品久久久不卡| 别揉我奶头~嗯~啊~动态视频| 高潮久久久久久久久久久不卡| 三级国产精品欧美在线观看| 乱人视频在线观看| 精品午夜福利视频在线观看一区| 变态另类成人亚洲欧美熟女| 床上黄色一级片| 99精品久久久久人妻精品| 国产精品电影一区二区三区| 2021天堂中文幕一二区在线观| 欧美性感艳星| svipshipincom国产片| 欧美日韩综合久久久久久 | 欧美日韩黄片免| 波多野结衣高清作品| 亚洲av免费高清在线观看| 午夜精品一区二区三区免费看| 三级国产精品欧美在线观看| 超碰av人人做人人爽久久 | 男人舔奶头视频| 90打野战视频偷拍视频| 亚洲av不卡在线观看| 国产精品一区二区三区四区免费观看 | 别揉我奶头~嗯~啊~动态视频| 99久久精品一区二区三区| 久久精品人妻少妇| 日本a在线网址| 校园春色视频在线观看| 日本黄色视频三级网站网址| 亚洲国产欧洲综合997久久,| 99久久精品国产亚洲精品| 国产成+人综合+亚洲专区| 91av网一区二区| 99国产精品一区二区蜜桃av| 国产久久久一区二区三区| 亚洲专区国产一区二区| 亚洲熟妇熟女久久| 国产成人aa在线观看| 免费av观看视频| 午夜影院日韩av| 欧美激情在线99| 法律面前人人平等表现在哪些方面| h日本视频在线播放| 久久久久久久久大av| 男插女下体视频免费在线播放| 欧美一级a爱片免费观看看| 一区二区三区高清视频在线| 日韩欧美免费精品| 亚洲第一电影网av| 日韩人妻高清精品专区| av视频在线观看入口| 丰满人妻熟妇乱又伦精品不卡| 久久香蕉精品热| 99精品欧美一区二区三区四区| 日韩有码中文字幕| 日韩精品中文字幕看吧| 岛国视频午夜一区免费看| 亚洲国产色片| 亚洲欧美日韩东京热| 国产主播在线观看一区二区| 欧美色视频一区免费| 波多野结衣巨乳人妻| 国产午夜精品久久久久久一区二区三区 | 成年女人毛片免费观看观看9| 国语自产精品视频在线第100页| 国产私拍福利视频在线观看| xxxwww97欧美| 久久国产精品影院| 亚洲熟妇中文字幕五十中出| 亚洲一区高清亚洲精品| 亚洲电影在线观看av| 国产成人福利小说| 亚洲av不卡在线观看| 99视频精品全部免费 在线| 天堂动漫精品| 精品午夜福利视频在线观看一区| 成年女人看的毛片在线观看| 我要搜黄色片| 天堂网av新在线| 九九在线视频观看精品| 乱人视频在线观看| 最近最新中文字幕大全电影3| 国内少妇人妻偷人精品xxx网站| 国产亚洲欧美在线一区二区| 国产精品永久免费网站| 此物有八面人人有两片| 一区二区三区免费毛片| 午夜免费男女啪啪视频观看 | 亚洲在线自拍视频| 亚洲欧美日韩无卡精品| 给我免费播放毛片高清在线观看| 香蕉久久夜色| 国产精品,欧美在线| 国产成人aa在线观看| 国产探花极品一区二区| 国产精品一区二区三区四区久久| 亚洲第一欧美日韩一区二区三区| 熟女电影av网| 欧美一区二区亚洲| 两个人看的免费小视频| 18禁黄网站禁片免费观看直播| 18美女黄网站色大片免费观看| 脱女人内裤的视频| 久久精品夜夜夜夜夜久久蜜豆| 免费观看的影片在线观看| 天堂网av新在线| 2021天堂中文幕一二区在线观| 亚洲精品影视一区二区三区av| 中文字幕久久专区| av欧美777| 日韩亚洲欧美综合| 又黄又爽又免费观看的视频| 麻豆国产97在线/欧美| 国产精品亚洲av一区麻豆| 又紧又爽又黄一区二区| 俄罗斯特黄特色一大片| 久久久精品大字幕| 丁香欧美五月| 久久九九热精品免费| 午夜免费男女啪啪视频观看 | 99久久精品一区二区三区| 亚洲成人中文字幕在线播放| 久久久国产成人免费| 一本综合久久免费| 变态另类丝袜制服| 色尼玛亚洲综合影院| 一进一出抽搐gif免费好疼| 超碰av人人做人人爽久久 | 日本 av在线| 久久久国产成人免费| 国产精品乱码一区二三区的特点| 美女大奶头视频| 啪啪无遮挡十八禁网站| 波多野结衣巨乳人妻| 日本 av在线| 久99久视频精品免费| 国内精品一区二区在线观看| 成人特级av手机在线观看| 国产精品久久电影中文字幕| 亚洲精品在线观看二区| 久久久色成人| 老司机在亚洲福利影院| 18禁在线播放成人免费| 成人欧美大片| 免费av毛片视频| 非洲黑人性xxxx精品又粗又长| 手机成人av网站| 女生性感内裤真人,穿戴方法视频| 99久国产av精品| 老汉色∧v一级毛片| 午夜亚洲福利在线播放| 亚洲av中文字字幕乱码综合| 久久精品国产清高在天天线| 免费看日本二区| 国产一区二区三区视频了| 免费观看人在逋| 夜夜躁狠狠躁天天躁| 18禁黄网站禁片免费观看直播| 日韩亚洲欧美综合| 欧美激情久久久久久爽电影| 90打野战视频偷拍视频| 中文字幕久久专区| 久久精品国产清高在天天线| 中出人妻视频一区二区| 99国产极品粉嫩在线观看| 欧美xxxx黑人xx丫x性爽| 在线国产一区二区在线| 一二三四社区在线视频社区8| 欧美成人免费av一区二区三区| 一进一出抽搐gif免费好疼| 一个人观看的视频www高清免费观看| 亚洲av中文字字幕乱码综合| 国产私拍福利视频在线观看| 国产成人啪精品午夜网站| 午夜福利高清视频| 男插女下体视频免费在线播放| 国产成+人综合+亚洲专区| 波多野结衣高清作品| 午夜福利成人在线免费观看| 亚洲乱码一区二区免费版| 久久精品91无色码中文字幕| 伊人久久精品亚洲午夜| 亚洲国产欧美人成| 淫妇啪啪啪对白视频| 中文字幕人成人乱码亚洲影| 桃色一区二区三区在线观看| 亚洲 国产 在线| 十八禁人妻一区二区| 又黄又粗又硬又大视频| 乱人视频在线观看| 小蜜桃在线观看免费完整版高清| 国产高清视频在线观看网站| 日韩欧美国产在线观看| 日本 av在线| 99久久精品国产亚洲精品| 亚洲中文字幕一区二区三区有码在线看| 最近最新中文字幕大全免费视频| 国产成+人综合+亚洲专区| 久久精品影院6| 亚洲第一电影网av| 少妇丰满av| 亚洲人成伊人成综合网2020| 悠悠久久av| 欧美成人一区二区免费高清观看| 久久久久国产精品人妻aⅴ院| 久久久久久久午夜电影| 国产探花在线观看一区二区| 亚洲激情在线av| 欧美3d第一页| 精品人妻1区二区| 亚洲av免费高清在线观看| 精品福利观看| 亚洲欧美日韩卡通动漫| 天天一区二区日本电影三级| 欧美性感艳星| 欧美日本视频| 国产一区二区三区在线臀色熟女| 波多野结衣巨乳人妻| 亚洲av免费在线观看| 一卡2卡三卡四卡精品乱码亚洲| 亚洲精品亚洲一区二区| 在线观看免费午夜福利视频| 不卡一级毛片| 国产精华一区二区三区| 在线天堂最新版资源| 两个人视频免费观看高清| 亚洲成av人片在线播放无| 久久精品影院6| 一本精品99久久精品77| 99久久精品热视频| 国产精品久久久久久久久免 | 午夜福利欧美成人| 特大巨黑吊av在线直播| 欧美日韩在线观看h| 日韩 亚洲 欧美在线| 亚洲国产高清在线一区二区三| 99热这里只有是精品50| 少妇被粗大猛烈的视频| 国产精品久久久久久久久免| 免费观看a级毛片全部| 麻豆精品久久久久久蜜桃| 国产爱豆传媒在线观看| 在现免费观看毛片| 夜夜看夜夜爽夜夜摸| 亚洲aⅴ乱码一区二区在线播放| 亚洲无线观看免费| 99九九线精品视频在线观看视频| 国产精品人妻久久久久久| 精品午夜福利在线看| av黄色大香蕉| 99热网站在线观看| 欧美日韩亚洲高清精品| 国产一区二区在线观看日韩| 久久久久久久国产电影| 久久久久精品久久久久真实原创| 亚洲av电影不卡..在线观看| 日韩伦理黄色片| 国产单亲对白刺激| 精品人妻视频免费看| 久久久久性生活片| 神马国产精品三级电影在线观看| 99热6这里只有精品| 1000部很黄的大片| 简卡轻食公司| 三级国产精品欧美在线观看| 免费观看的影片在线观看| 久久久久久久久久人人人人人人| 嫩草影院入口| av天堂中文字幕网| 一本久久精品| 永久网站在线| 五月天丁香电影| av天堂中文字幕网| 亚洲va在线va天堂va国产| av天堂中文字幕网| 一区二区三区乱码不卡18| 黄色日韩在线| 日韩欧美精品免费久久| 在线天堂最新版资源| 午夜福利在线观看吧| 精品酒店卫生间| 两个人视频免费观看高清| 亚洲在久久综合| 成人漫画全彩无遮挡| 赤兔流量卡办理| 精品午夜福利在线看| 国产精品麻豆人妻色哟哟久久 | 成人午夜精彩视频在线观看| 丰满人妻一区二区三区视频av| 亚洲国产精品国产精品| 国产 亚洲一区二区三区 | 内地一区二区视频在线| 日韩不卡一区二区三区视频在线| 久久精品国产亚洲av天美| 久久久久精品久久久久真实原创| 国产在线一区二区三区精| 亚洲国产av新网站| 乱码一卡2卡4卡精品| 国产精品一及| 十八禁国产超污无遮挡网站| 高清av免费在线| 联通29元200g的流量卡| 国国产精品蜜臀av免费| 欧美日本视频| 日本av手机在线免费观看| 国产精品99久久久久久久久| 综合色丁香网| 国产精品.久久久| 卡戴珊不雅视频在线播放| 白带黄色成豆腐渣| 欧美成人一区二区免费高清观看| 欧美日韩一区二区视频在线观看视频在线 | 国产综合精华液| 欧美+日韩+精品| 国产黄色小视频在线观看| 久久精品国产亚洲av涩爱| 超碰av人人做人人爽久久| 欧美激情国产日韩精品一区| 一二三四中文在线观看免费高清| 人人妻人人澡人人爽人人夜夜 | 内射极品少妇av片p| 少妇裸体淫交视频免费看高清| 一个人看视频在线观看www免费| 超碰av人人做人人爽久久| 久久久久久久久久成人| 少妇猛男粗大的猛烈进出视频 | 免费av毛片视频| 乱码一卡2卡4卡精品| 免费不卡的大黄色大毛片视频在线观看 | av免费观看日本| 91精品国产九色| 少妇丰满av| 丝瓜视频免费看黄片| 久久99热这里只有精品18| 精品亚洲乱码少妇综合久久| 亚洲精品色激情综合| 欧美97在线视频| 一级a做视频免费观看| 中文字幕av在线有码专区| 综合色av麻豆| 精品少妇黑人巨大在线播放| 美女内射精品一级片tv| 免费观看在线日韩| 熟女人妻精品中文字幕| 黄片无遮挡物在线观看| 九九爱精品视频在线观看| 久久久久久伊人网av| 国产亚洲av片在线观看秒播厂 | 内地一区二区视频在线| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 可以在线观看毛片的网站| 97在线视频观看| 亚洲国产日韩欧美精品在线观看| 久久99热这里只频精品6学生| 色综合色国产| 亚洲欧美清纯卡通| 中文字幕av成人在线电影| a级毛色黄片| 日本午夜av视频| 直男gayav资源| 丰满人妻一区二区三区视频av| 黄色欧美视频在线观看| 嫩草影院精品99| 免费看光身美女| 一个人观看的视频www高清免费观看| 最近中文字幕2019免费版| 亚洲最大成人手机在线| 亚洲av中文av极速乱| 亚洲成人一二三区av| 日韩精品有码人妻一区| 国产日韩欧美在线精品| or卡值多少钱| 色综合色国产| 午夜精品在线福利| 国产黄色免费在线视频| 欧美日韩亚洲高清精品| 日韩大片免费观看网站| 天堂av国产一区二区熟女人妻| 国产高清国产精品国产三级 | 18+在线观看网站| 观看美女的网站| 人人妻人人澡人人爽人人夜夜 | 人人妻人人澡人人爽人人夜夜 | 夫妻性生交免费视频一级片| 精品一区二区三区视频在线| 精品国内亚洲2022精品成人| 丰满人妻一区二区三区视频av| 啦啦啦韩国在线观看视频| 日韩强制内射视频| 九草在线视频观看| 高清欧美精品videossex| 中文字幕av成人在线电影| 超碰97精品在线观看| 九九久久精品国产亚洲av麻豆| 国产精品麻豆人妻色哟哟久久 | 欧美不卡视频在线免费观看| 禁无遮挡网站| freevideosex欧美| 亚洲久久久久久中文字幕| 久久久久久久久久黄片| 免费av观看视频| 亚洲欧美日韩卡通动漫| 国产午夜福利久久久久久| 免费黄色在线免费观看| 日韩中字成人| 亚洲在线观看片| 春色校园在线视频观看| 91午夜精品亚洲一区二区三区| 又爽又黄无遮挡网站| 免费av毛片视频| av专区在线播放| 一个人免费在线观看电影| av福利片在线观看| 国内精品宾馆在线| 国产三级在线视频| 直男gayav资源| 亚洲内射少妇av| 高清视频免费观看一区二区 | 在线观看av片永久免费下载| 国产单亲对白刺激| 日韩欧美三级三区| 日韩成人伦理影院| 国产精品久久视频播放| 日韩av不卡免费在线播放| 欧美日韩亚洲高清精品| av又黄又爽大尺度在线免费看| 18+在线观看网站| 少妇猛男粗大的猛烈进出视频 | 日本猛色少妇xxxxx猛交久久| 欧美+日韩+精品| 最后的刺客免费高清国语| 午夜福利高清视频| 狠狠精品人妻久久久久久综合| 久久久久久伊人网av| 亚洲av中文字字幕乱码综合| 一二三四中文在线观看免费高清| 纵有疾风起免费观看全集完整版 | 久久精品人妻少妇| 精品99又大又爽又粗少妇毛片| 国产精品人妻久久久影院| 日韩成人伦理影院| 一区二区三区乱码不卡18| 男人舔奶头视频| 国产精品蜜桃在线观看| 国产成人精品婷婷| 亚洲最大成人中文| 亚洲最大成人手机在线| 精品人妻视频免费看| 国产色爽女视频免费观看| 国产亚洲午夜精品一区二区久久 | 黄片无遮挡物在线观看| 一级av片app| 亚洲自拍偷在线| 乱码一卡2卡4卡精品| 国产黄a三级三级三级人| 国产伦精品一区二区三区四那| 国产av不卡久久| 日韩国内少妇激情av| 亚洲四区av| 黄片无遮挡物在线观看| 国产一区有黄有色的免费视频 | 久久精品熟女亚洲av麻豆精品 | 中文资源天堂在线| 能在线免费看毛片的网站| freevideosex欧美| 精品人妻一区二区三区麻豆| 淫秽高清视频在线观看| 久久久久久久久久人人人人人人| 久久精品国产亚洲网站| 国产精品一及| 高清视频免费观看一区二区 | 自拍偷自拍亚洲精品老妇| 精品人妻一区二区三区麻豆| 听说在线观看完整版免费高清| 国产老妇女一区| 国产乱来视频区| 99久久中文字幕三级久久日本| 亚洲欧美成人精品一区二区| 一级毛片我不卡| 亚洲乱码一区二区免费版| 国产高清三级在线| 美女xxoo啪啪120秒动态图| av.在线天堂| 久久精品久久久久久久性| 丰满人妻一区二区三区视频av| 国产伦一二天堂av在线观看| 午夜精品在线福利| 深夜a级毛片| 久久久久久久亚洲中文字幕| 久久久精品欧美日韩精品| 免费大片黄手机在线观看| 久久国内精品自在自线图片| 亚洲精品自拍成人| 成年av动漫网址| 亚洲精品色激情综合| 91精品一卡2卡3卡4卡| 久久韩国三级中文字幕| 欧美精品一区二区大全| 欧美 日韩 精品 国产| 国产精品一及| 日韩在线高清观看一区二区三区| 亚洲精品影视一区二区三区av| 边亲边吃奶的免费视频| 国产精品久久久久久久久免| 国产精品人妻久久久影院| 精品少妇黑人巨大在线播放| 高清欧美精品videossex| 免费高清在线观看视频在线观看| 亚洲不卡免费看| 色哟哟·www| 国产成人freesex在线| 婷婷色综合www| 人人妻人人澡欧美一区二区| 亚洲国产精品国产精品| a级毛片免费高清观看在线播放| www.av在线官网国产| 在线 av 中文字幕| 国产 一区精品| 男女那种视频在线观看| 免费看美女性在线毛片视频| 美女脱内裤让男人舔精品视频| 秋霞在线观看毛片| 精品一区二区免费观看| 国产亚洲精品av在线| 欧美xxxx黑人xx丫x性爽| 亚洲激情五月婷婷啪啪| 又爽又黄a免费视频| 亚洲精品国产成人久久av| 亚洲成色77777| 51国产日韩欧美| 久久这里只有精品中国| 欧美变态另类bdsm刘玥| 国产成人a∨麻豆精品| 婷婷色麻豆天堂久久| 尤物成人国产欧美一区二区三区| 91久久精品国产一区二区三区| 亚洲经典国产精华液单| 成人性生交大片免费视频hd| 男女下面进入的视频免费午夜| av专区在线播放| 午夜福利视频精品| 亚洲精品乱码久久久久久按摩| 99热这里只有精品一区| 97在线视频观看| 欧美最新免费一区二区三区| 国产中年淑女户外野战色| 美女国产视频在线观看| 亚洲国产日韩欧美精品在线观看| 久久久久九九精品影院| 中文欧美无线码| 成人美女网站在线观看视频|