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

    一種分裂形式CPR 格式在欠解析流動(dòng)中的穩(wěn)定性研究

    2024-01-09 13:19:36賈斐然朱華君燕振國(guó)石京昶
    關(guān)鍵詞:散度激波高階

    賈斐然,朱華君,嚴(yán) 紅,燕振國(guó),*,石京昶

    (1.西北工業(yè)大學(xué) 動(dòng)力與能源學(xué)院,西安 710000;2.空氣動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室,綿陽(yáng) 621000)

    0 引言

    計(jì)算流體力學(xué)(Computational Fluid Dynamics,CFD)是一種通過數(shù)值手段研究流體力學(xué)問題的方法,是對(duì)理論研究和實(shí)驗(yàn)研究的補(bǔ)充。近些年來,為了對(duì)流動(dòng)問題進(jìn)行精細(xì)模擬,高階格式因其數(shù)值耗散誤差和色散誤差小的特性,逐漸成為CFD 中的一個(gè)重要研究方向[1]。

    重構(gòu)修正過程(correction procedure via reconstruction,CPR)方法的思想最早由Huynh[2]提出,當(dāng)時(shí)稱之為通量重構(gòu)(flux reconstruction,FR)方法。Huynh 將其用于求解結(jié)構(gòu)網(wǎng)格上的雙曲守恒律方程,王志堅(jiān)等[3-4]將此方法進(jìn)行了推廣,提出適用于非結(jié)構(gòu)網(wǎng)格的提升配點(diǎn)懲罰法(lifting collocation penalty formulation,LCP),這兩種方法聯(lián)系緊密,被統(tǒng)稱為CPR 方法[5]。CPR 方法的優(yōu)勢(shì)主要在于三個(gè)方面:第一,在求解線性標(biāo)量方程時(shí),通過調(diào)整修正函數(shù)類型,CPR 方法可以恢復(fù)成間斷伽遼金(discontinuous Galerkin,DG)[6]、譜差分(spectral difference,SD)等方法[7],但不同于DG 方法在求解過程中涉及非線性項(xiàng)的積分,CPR 方法是差分型方法,計(jì)算量小且計(jì)算復(fù)雜度低;第二,CPR 方法對(duì)復(fù)雜邊界有很好的適應(yīng)性,便于推廣到非結(jié)構(gòu)網(wǎng)格[8];第三,該方法具有良好的緊致性,便于進(jìn)行并行計(jì)算[9]。

    CPR 方法或DG 方法在求解非線性問題時(shí),即使流動(dòng)中沒有出現(xiàn)像激波這樣的不連續(xù)性情況,也可能會(huì)由于離散非線性對(duì)流項(xiàng)引起的混淆誤差的累積而引發(fā)穩(wěn)定性問題[10]。這類問題在欠解析流動(dòng)中尤為突出。為了增強(qiáng)欠解析流動(dòng)模擬的穩(wěn)定性,目前主要存在的穩(wěn)定機(jī)制有濾波[6]、過積分[11]以及分裂形式[12]。

    Gassner 等[13]研究了高階DG 離散在欠解析湍流模擬中的精度,發(fā)現(xiàn)低階近似顯示出難以接受的數(shù)值離散誤差,而高階離散受混淆誤差的影響,往往是不穩(wěn)定的。對(duì)流項(xiàng)中的非線性項(xiàng)使得相對(duì)較低波數(shù)的模態(tài)相互非線性疊加,產(chǎn)生較高波數(shù)的模態(tài),甚至是當(dāng)前基函數(shù)無(wú)法描述的模態(tài)。高于Nyquist 波數(shù)的未分辨波數(shù)被“混淆”到分辨波數(shù)范圍。Blaisdell 等[14]在譜方法中提出一種分裂形式,發(fā)現(xiàn)其減少了高波數(shù)范圍內(nèi)的混淆誤差,Gassner 和Abe 等[10,12,15]將這種思路推廣到DG 和CPR 方法中。Gassner 等[12]構(gòu)造了一種高階分裂形式間斷伽遼金譜元法(discontinuous Galerkin spectral element method,DGSEM)框架,并對(duì)無(wú)黏Taylor-Green 渦(TGV)算例進(jìn)行模擬,結(jié)果表明,與標(biāo)準(zhǔn)格式相比,這種新的分裂形式在求解欠解析湍流渦主導(dǎo)的流動(dòng)時(shí)增強(qiáng)了模擬的穩(wěn)定性。Winters 等[15]研究高階DG 方法在欠解析湍流計(jì)算中的精度和穩(wěn)定性,考慮無(wú)黏TGV 來分析DG 方法在極高雷諾數(shù)下進(jìn)行隱式大渦模擬(implicit large eddy simulation,ILES)的能力。為了抑制混淆誤差,采用過積分[13]和分裂形式兩種方法對(duì)控制方程進(jìn)行離散,并比較了這兩種方法在無(wú)黏TGV 算例中的表現(xiàn),結(jié)果表明,分裂形式具有更好的穩(wěn)定性。Chan 等[16]將高階熵穩(wěn)定DG 格式用于求解欠解析流動(dòng),這種采用LG(Legendre-Gauss)點(diǎn)并經(jīng)過熵投影的格式具有很好的魯棒性,與采用LGL(Legendre-Gauss-Lobatto)點(diǎn)的高階熵穩(wěn)定DGSEM 格式[12]相比,能夠分辨尺度更小的流動(dòng)結(jié)構(gòu)。

    CPR 方法在欠解析流動(dòng)中的研究相比于DG 方法較少。Ball 和Jameson 研究了兩種新的高階FR 格式[17-18]即優(yōu)化能量穩(wěn)定通量重構(gòu)(optimized energy stable flux reconstruction,OESFR)和優(yōu)化通量重構(gòu)(optimized flux reconstruction,OFR)在粗網(wǎng)格上求解黏性TGV 的能力,并與能量穩(wěn)定通量重構(gòu)(energy stable flux reconstruction,ESFR)方法進(jìn)行了比較。結(jié)果表明,OFR 格式精度最高,所計(jì)算的能譜與參考能譜吻合得最好,且在高波數(shù)下優(yōu)于ESFR 方法的能譜。Abe 等[10]針對(duì)可壓縮Euler 和Navier-Stokes(NS)方程,構(gòu)造并證明了一種穩(wěn)定且守恒的FR 格式。這種格式采用LGL 解點(diǎn),對(duì)對(duì)流項(xiàng)采用分裂形式,基于多項(xiàng)式分析嚴(yán)格推導(dǎo)了滿足離散守恒和動(dòng)能保持性質(zhì)的充分條件,通過黏性TGV 算例證明:基于這種分裂形式的FR 框架,使用無(wú)耗散動(dòng)能保持(kinetic energy preservation,KEP)通量[19],在非常高階(p15)和相對(duì)粗糙網(wǎng)格下的模擬仍然是穩(wěn)定的。Abe 提出了一種基于LG 點(diǎn)的滿足離散守恒律的分裂形式CPR 格式,但并未研究其求解欠解析流動(dòng)的特性。

    Ramírez 等[20]提出了一種通用的子單元限制策略來構(gòu)造魯棒的高精度節(jié)點(diǎn)DG 格式,主要思路是構(gòu)造兼容的低階有限體積(FV)離散,并與高階DG 進(jìn)行凸混合。這種策略可以有效處理強(qiáng)激波,在KPP 問題[21]、湍流和高超聲速Euler 模擬,及以激波和湍流為特征的MHD 問題上有很好表現(xiàn)。Zhu 等[22–24]針對(duì)基于LG 點(diǎn)的CPR 方法,提出了一類先驗(yàn)子單元限制格式。首先,利用激波偵測(cè)器檢測(cè)存在間斷的問題單元;然后,將問題單元分解為非均勻子單元,每個(gè)子單元包含一個(gè)解點(diǎn);最后,對(duì)光滑單元使用CPR 方法計(jì)算,對(duì)問題單元使用緊致非均勻非線性加權(quán)(compact nonuniform nonlinear weighted,CNNW)插值格式計(jì)算。這種新策略在分辨率和激波捕捉魯棒性之間具有良好的平衡。

    分裂形式CPR 格式在欠解析流動(dòng)中的研究較少,且主要以LGL 點(diǎn)作為解點(diǎn)。本文研究了以LG 點(diǎn)為解點(diǎn)的分裂形式CPR 格式在欠解析流動(dòng)中的應(yīng)用。首先,對(duì)比了基于LG 點(diǎn)的分裂形式和散度形式在求解欠解析流動(dòng)時(shí)的穩(wěn)定性,展現(xiàn)了分裂形式在增強(qiáng)格式穩(wěn)定性方面的優(yōu)勢(shì);其次,從分辨率和穩(wěn)定性的角度比較了LG 點(diǎn)和LGL 點(diǎn)分裂形式CPR 格式;最后,首次將基于LG 點(diǎn)的分裂形式CPR 格式與Zhu 等提出的CNNW 子單元限制格式相結(jié)合,求解含激波的Kelvin-Helmholtz 不穩(wěn)定性問題。

    1 控制方程與數(shù)值方法

    1.1 控制方程

    考慮守恒形式的三維N-S 方程:

    式中,U為守恒變量,F(xiàn)、G、H分別為x、y、z方向的無(wú)黏通量,F(xiàn)v、Gv、Hv分別為x、y、z方向的黏性通量,具體形式如下:

    式中,ρ 為密度,u、v、w分別為x、y、z方向的速度,p為壓力,E為單位質(zhì)量流體的總內(nèi)能。本文只涉及理想氣體,比熱比 γ=1.4。黏性應(yīng)力滿足:

    式中,μ為動(dòng)力黏度。

    1.2 高階CPR 格式

    本小節(jié)以一維守恒形式Euler 方程為例,介紹CPR 格式的構(gòu)造方法。守恒方程形式為:

    計(jì)算域 [a,b],首先將其剖分為N個(gè)區(qū)間,將其中第j個(gè)區(qū)間設(shè)為Ij=[xj-1/2,xj+1/2],j=1,·,N。為 了便于討論,將每個(gè)區(qū)間映射到標(biāo)準(zhǔn)單元I=[-1,1]上。假設(shè)每個(gè)區(qū)間Ij上的點(diǎn)x到標(biāo)準(zhǔn)單元I的點(diǎn) ξ存在一個(gè)線性映射關(guān)系:

    式中,xj=(xj-1/2+xj+1/2)/2 是區(qū)間Ij的中心,hj=xj+1/2-xj-1/2為區(qū)間Ij的長(zhǎng)度。

    在區(qū)間Ij上定義K個(gè)解點(diǎn)xj,k,k=1,·,K,解點(diǎn)上對(duì)應(yīng)的守恒變量為Uj,k,k=1,·,K。為了便于分析,本文對(duì)所有區(qū)間均采用相同類型的解點(diǎn),即LG 點(diǎn)或LGL 點(diǎn)。設(shè)標(biāo)準(zhǔn)單元上的解點(diǎn)位置為 ξk,k=1,·,K,則區(qū)間Ij上的解點(diǎn)位置由下式計(jì)算得到:

    由解點(diǎn)處的守恒變量Uj,k,k=1,·,K,通過構(gòu)造K-1 次Lagrange 插值多項(xiàng)式來逼近Uj:

    其中l(wèi)k(ξ)是Lagrange 插值基函數(shù),形式如下:

    同理,區(qū)間Ij上的通量函數(shù)由下式逼近:

    通量多項(xiàng)式Fj(ξ)的構(gòu)造只依賴于區(qū)間內(nèi)部的解點(diǎn)通量值,在相鄰區(qū)間交界面處的通量值不連續(xù),即Fj(1)≠Fj+1(-1)。因此,需要構(gòu)造一個(gè)連續(xù)通量多項(xiàng)式表達(dá)式為:

    式中,gL、gR為邊界的左右修正函數(shù),滿足:

    修正函數(shù)的形式一般有以下三種:

    式中,PK是K階Legendre 多項(xiàng)式。gDG、gGa和g2分別是可將CPR 格式恢復(fù)成DG 格式、SD 格式和Huynh 型格式的修正函數(shù)[2]。

    為了保證相鄰區(qū)間交界面處的連續(xù)性,需要引入Riemann 通量常用的Riemann 通量有局部Lax-Friedrichs(LLF)、Roe 通量[25]等。

    最后,得到Euler 方程的半離散形式:

    針對(duì)上式,本文使用三階TVD Runge-Kutta 格式[26-27]進(jìn)行時(shí)間離散。以上為一維CPR 格式的實(shí)現(xiàn)過程,二維及三維CPR 格式可以直接由一維形式擴(kuò)展而來[10,22]。

    1.3 分裂形式CPR 格式

    同樣以一維Euler 方程為例(計(jì)算坐標(biāo)下),分裂形式[10]是將方程中的對(duì)流項(xiàng)分裂成守恒形式和非守恒形式的組合,其一般形式為:

    式中,P為壓力項(xiàng),將單獨(dú)處理,不會(huì)被分裂。

    常用的分裂形式類型有Fe 分裂[28]、Bl 分裂[14]、KG1 分裂[29]和KG2 分裂[29]等。本文只考慮第一種分裂形式,表達(dá)式為:

    其中,?={1,u,H}T。式(21)中的 (Split)C可由式(22)代替。

    (Split)C表示分裂形式下重構(gòu)通量(通量散度)的ξ方向?qū)?shù)。首先,使用解點(diǎn)處的守恒變量值計(jì)算解點(diǎn)處分裂形式的通量散度項(xiàng)。在第n個(gè)單元上,假設(shè)符號(hào)函數(shù) {fn,gn,hn} 代表 {ρ,u,?},則式(22)中等號(hào)右端項(xiàng)在解點(diǎn)處可由以下形式表示:

    式中,引入符號(hào)ISP;K[ψn] 表示一個(gè)任意的函數(shù) ψn的K階多項(xiàng)式近似:

    式(23)的第一項(xiàng)可寫為:

    因此,解點(diǎn)處分裂形式的重構(gòu)通量散度項(xiàng)計(jì)算如下:

    Abe[10]證明了以LGL 點(diǎn)作為解點(diǎn)時(shí),F(xiàn)e 分裂形式滿足離散守恒律。同時(shí)也指出:如果以LG 點(diǎn)作為解點(diǎn),要使分裂形式滿足離散守恒律,需要對(duì)邊界通量項(xiàng)進(jìn)行修正,形式如下:

    式中I[?] 即為ISP;K[?]。

    1.4 子單元限制和激波偵測(cè)器

    本文所使用的子單元限制策略的主要思想是,利用激波偵測(cè)器偵測(cè)出流場(chǎng)中存在的大梯度或間斷的單元,并將其標(biāo)記為問題單元,然后將這些問題單元?jiǎng)澐譃榉堑染嘧訂卧?,每個(gè)子單元包含一個(gè)解點(diǎn),最后在子單元上添加限制手段,具體限制方式在文獻(xiàn)[24]中給出。光滑單元用分裂形式CPR 格式計(jì)算,而問題單元?jiǎng)t使用二階CNNW 格式計(jì)算以捕捉激波。值得注意的是,基于LG 點(diǎn)的分裂形式與子單元限制結(jié)合時(shí)并不直接滿足離散守恒律,仍然需要對(duì)分裂形式使用邊界通量修正,在2.1.2 小節(jié)對(duì)此進(jìn)行了數(shù)值驗(yàn)證。

    本文采用的激波偵測(cè)器為最高模態(tài)衰減(highest modal decay,MDH)。該方法利用解多項(xiàng)式的最高模態(tài)能量系數(shù)在光滑區(qū)域衰減較快、非光滑區(qū)域衰減較慢的特點(diǎn)[30],通過設(shè)置閾值來判斷單元內(nèi)是否存在大梯度或間斷。為了避免奇偶效應(yīng),Hennemann 等[31]使用最高和次高模態(tài)系數(shù)改進(jìn)了此方法,本文采用這一改進(jìn)方法。

    2 數(shù)值測(cè)試

    使用一維Sod 激波管、二維對(duì)流等熵渦[22]問題來測(cè)試格式的離散守恒律和數(shù)值精度;使用二維欠解析旋渦流動(dòng)[32-33]和三維黏性TGV[10,18]算例來測(cè)試格式的穩(wěn)定性;使用二維無(wú)黏Kelvin-Helmholtz 不穩(wěn)定性算例來測(cè)試分裂形式結(jié)合子單元限制技術(shù)的激波捕捉能力。為了便于描述,所測(cè)試的計(jì)算方案按以下方式命名:對(duì)流項(xiàng)-解點(diǎn)類型-修正函數(shù)-無(wú)黏通量。例如,Div-LG-gDG-Roe 表示對(duì)流項(xiàng)采用散度形式、解點(diǎn)選擇LG 點(diǎn)、修正函數(shù)選擇gDG、無(wú)黏通量選擇Roe 通量。表1 列出了本文所采用的算例及其對(duì)應(yīng)的計(jì)算條件設(shè)置。此外,對(duì)于黏性通量,均使用BR2 格式[34]。

    表1 算例類型及對(duì)應(yīng)的計(jì)算條件設(shè)置Table 1 Simulation cases and the corresponding calculation condition settings

    2.1 離散守恒律及精度測(cè)試

    2.1.1 Sod 激波管

    Sod 激波管問題計(jì)算域?yàn)?0 ≤x≤1,被劃分成等距的200 個(gè)網(wǎng)格。初始條件如下:

    該算例存在精確解[35]。左右邊界施加Dirichlet邊界條件,計(jì)算時(shí)間T=0.2,CFL=0.1。解多項(xiàng)式階為p2(即三階空間精度),總自由度為600。為了便于評(píng)估數(shù)值穩(wěn)定性,不采用任何激波捕捉方法。

    圖1 是使用LG 點(diǎn)和修正函數(shù)g2的計(jì)算結(jié)果。圖1(a)表明:基于LG 點(diǎn)的分裂形式CPR 格式直接求解Sod 激波管問題時(shí),并不能正確計(jì)算激波位置,即不滿足離散守恒律。圖1(b)是經(jīng)過邊界通量修正[10]后的計(jì)算結(jié)果,此時(shí)散度形式和Fe 分裂形式均能滿足離散守恒律。

    圖1 T=0.2時(shí)Sod 激波管問題的密度分布.修正采用LG 點(diǎn)、g2 修正函數(shù)和Roe 通量,總自由度為600(p2)Fig.1 Density profiles of the Sod shock tube problem at T=0.2:(a) without boundary flux fix;(b) with boundary flux fix.LG points,g2 correction function and Roe flux are used,and the total number of DoFs is 600(p2)

    2.1.2 對(duì)流等熵渦

    對(duì)流等熵渦問題除了可以用來驗(yàn)證格式的離散守恒律[22]外,還可以用于測(cè)試精度。該算例是在一個(gè)均勻流動(dòng)的基礎(chǔ)上添加一個(gè)等熵旋渦擾動(dòng)。均勻流動(dòng)設(shè)置如下:

    計(jì)算域 [-10,10]×[-10,10],邊界均為周期邊界,CFL=0.1,計(jì)算時(shí)間T=0.1。

    首先用該算例測(cè)試離散守恒律。全局離散守恒律誤差表達(dá)式如下:

    其中,q可以用于表示 ρ、ρu等守恒變量,下標(biāo)“0”表示初始時(shí)刻的守恒變量。

    表2 給出了對(duì)流項(xiàng)的不同形式在使用LG 點(diǎn)時(shí)的離散守恒律誤差,其中Div-subcell 和SplitFe-subcell分別表示結(jié)合子單元限制的散度形式和分裂形式CPR 格式。當(dāng)離散守恒律誤差接近機(jī)器零時(shí),認(rèn)為格式滿足離散守恒律。由表2 可知,當(dāng)分裂形式與子單元限制相結(jié)合時(shí),也需要經(jīng)過邊界通量修正才能滿足離散守恒律。

    表2 對(duì)流等熵渦的全局離散守恒律誤差Table 2 Errors of the global discrete conservation law in the computation of convecting isentropic vortex

    表3 和表4 分別是基于LG 點(diǎn)的散度形式和Fe 分裂形式的精度測(cè)試(p4)結(jié)果,這兩種對(duì)流項(xiàng)形式計(jì)算的誤差與收斂階基本相同。表5 和表6 則是基于LGL 點(diǎn)的散度形式和Fe 分裂形式的精度測(cè)試(p4)結(jié)果,F(xiàn)e 分裂形式的L1誤差小于散度形式,而L∞誤差大于散度形式。最后對(duì)比表4 和表6,發(fā)現(xiàn)基于LG 點(diǎn)的格式誤差明顯小于LGL 點(diǎn)(接近一個(gè)量級(jí))。

    表3 基于LG 點(diǎn)的散度形式精度測(cè)試(p4)Table 3 Accuracy test of divergence form based on LG points (p4)

    表4 基于LG 點(diǎn)的Fe 分裂形式精度測(cè)試(p4)Table 4 Accuracy test of the Fe split form based on LG points (p4)

    表5 基于LGL 點(diǎn)的散度形式精度測(cè)試(p4)Table 5 Accuracy test of the divergence form based on LGL points (p4)

    表6 基于LGL 點(diǎn)的Fe 分裂形式精度測(cè)試(p4)Table 6 Accuracy test of the Fe split form based on LGL points (p4)

    2.2 穩(wěn)定性測(cè)試

    2.2.1 欠解析旋渦流動(dòng)

    為了評(píng)估格式對(duì)實(shí)際流動(dòng)的欠解析模擬的影響,使用二維可壓縮Euler 方程和非定常流入邊界條件,模擬渦流沿下游傳播的被動(dòng)發(fā)生器。該問題是無(wú)黏的,使用Euler 方程來研究數(shù)值方法在黏性消失極限(極高雷諾數(shù))下的性質(zhì)和行為[15,36]。算例設(shè)置如下:

    計(jì)算域?yàn)?[0,20π]×[-π,π],計(jì)算網(wǎng)格數(shù) 120×12,CFL=0.5,T=150.0 。在y=±π位置施加壁面邊界條件,x=0位置施加流入邊界條件:

    式中,ρ∞=1、u∞=1是自由流密度和均勻流速度,是自由流靜壓,用于通過聲速c∞=u∞Ma-1定義流動(dòng)的參考馬赫數(shù)。此外,入口處擾動(dòng)參數(shù)定義為A=1/2、K=5、?=1,出口邊界為初始條件與出口邊界相同。T=20π左右時(shí),流場(chǎng)將達(dá)到漸近狀態(tài)。

    本算例使用的多項(xiàng)式階為p5。主要考慮三種CPR格式中常見的Riemann 求解器,即HLL 通量[37-38]、Roe 通量和LLF 通量。在本算例中,在多項(xiàng)式高階和欠解析條件下,Riemann 求解器會(huì)有不同的特性。本算例中選擇馬赫數(shù)為0.03,因?yàn)樵诘婉R赫數(shù)下,會(huì)放大HLL 和LLF 通量產(chǎn)生的偽反射和數(shù)值噪聲,從而更好地展現(xiàn)數(shù)值不穩(wěn)定性帶來的影響[32]。除了使用多種Riemann 求解器外,還使用兩種不同對(duì)流項(xiàng)形式的CPR 格式(散度形式和Fe 分裂)和兩種解點(diǎn)類型(LG 點(diǎn)和LGL 點(diǎn))。為了評(píng)估這些方面的影響,考慮流場(chǎng)中的渦量分布

    圖2 給出了在Roe 通量下,采用不同對(duì)流項(xiàng)形式和不同解點(diǎn)的渦量計(jì)算結(jié)果。所有情況均能穩(wěn)定計(jì)算至結(jié)束。比較圖2 中的第一和第三小圖,發(fā)現(xiàn)采用LG 點(diǎn)對(duì)流場(chǎng)中的微小結(jié)構(gòu)捕捉得更清楚。

    圖2 基于Roe 通量,不同對(duì)流形式和解點(diǎn)下的渦量場(chǎng)比較:LGL 點(diǎn),散度形式;LGL 點(diǎn),F(xiàn)e 分裂形式;LG 點(diǎn),散度形式;LG 點(diǎn),F(xiàn)e 分裂形式Fig.2 Comparison of the vorticity fields computed with different convention from and solution point combination based on the Roe flux:LGL points,Div form;LGL points,SplitFe form;LG points,Div form;LG points,SplitFe form

    圖3 是在LLF 通量下的渦量計(jì)算結(jié)果。值得注意的是,散度形式下LGL 點(diǎn)和LG 點(diǎn)分別在T=15.66和T=56.84時(shí)崩潰,因此沒有畫出其渦量分布。這是由于LLF 通量在靠近出口邊界處會(huì)產(chǎn)生偽反射(圖中被紅圈圈出的),這些偽反射與傳入的旋渦非線性相互作用,導(dǎo)致小尺度噪聲,使得計(jì)算不穩(wěn)定。此外,即使存在明顯的小尺度噪聲,分裂形式均能穩(wěn)定計(jì)算到T=150。

    圖3 基于LLF 通量和Fe 分裂形式、LGL 和LG 解點(diǎn)下的渦量場(chǎng)Fig.3 Comparison of the vorticity fields computed with the LGL and LG solution points based on the LLF flux and Fe split form

    圖4 是在HLL 通量下的渦量計(jì)算結(jié)果。值得注意的是,散度形式下LGL 點(diǎn)和LG 點(diǎn)分別在T=15.56和T=50.23時(shí)崩潰,因此沒有畫出其渦量分布。兩種分裂形式均能穩(wěn)定計(jì)算到T=150。

    圖4 基于HLL 通量和Fe 分裂形式、LGL 和LG 解點(diǎn)下的渦量場(chǎng)Fig.4 Comparison of the vorticity fields computed with the LGL and LG solution points based on the HLL flux and Fe split form

    2.2.2 黏性Taylor-Green 渦

    在本小節(jié)中,我們模擬了黏性TGV 問題。該問題通常用于研究旋渦動(dòng)力學(xué)和湍流轉(zhuǎn)捩與衰減[39]。該問題由最初包含光滑渦量分布的立方體流體組成,能夠反映數(shù)值格式的魯棒性以及對(duì)多尺度流動(dòng)結(jié)構(gòu)的模擬能力。在美國(guó)航空航天學(xué)會(huì)航空航天科學(xué)會(huì)議[40]上舉行的計(jì)算流體動(dòng)力學(xué)高階方法國(guó)際研討會(huì)上,該算例用于評(píng)估湍流模擬方法。許多作者已經(jīng)成功使用高階模擬方法預(yù)測(cè)該流場(chǎng)[13,18,41-42]。本文將使用TGV 來比較散度形式和分裂形式在湍流欠解析模擬中的精度和穩(wěn)定性。初始條件設(shè)置如下:

    其中,馬赫數(shù)和雷諾數(shù)分別設(shè)為0.1 和1 600,L為參考長(zhǎng)度。計(jì)算域是 -πL≤x,y,z≤πL,被劃分為互不重疊的均勻六面體單元。在本測(cè)試中,包括內(nèi)部解點(diǎn)在內(nèi)的自由度總數(shù)固定在 643左右,以便在不同階的解多項(xiàng)式之間進(jìn)行比較。表7 列出了網(wǎng)格數(shù)與解多項(xiàng)式階的對(duì)應(yīng)關(guān)系。計(jì)算時(shí)間為 0 ≤T≤20,使用固定時(shí)間步長(zhǎng) ?t=3.14×10-4。本算例只考慮了Fe 分裂形式。為了更好地描述該問題隨時(shí)間的演化過程,需要定義以下幾個(gè)參數(shù):

    表7 總單元數(shù)和解多項(xiàng)式階數(shù)Table 7 Total number of cells and order of the solution polynomial

    式中,u為速度矢量,ω為渦量矢量。

    圖5 分別給出了不同時(shí)刻下,多項(xiàng)式階為p3 的方案SplitFe-LG-g2-Roe 的z方向渦量等值面圖。從T=0開始,初始流場(chǎng)存在明顯的大渦結(jié)構(gòu),隨后不斷地拉伸和旋轉(zhuǎn),逐漸破裂成較小的渦結(jié)構(gòu)[43]。

    圖5 不同時(shí)刻黏性Taylor-Green 渦 z 方向渦量等值面圖Fig.5 Iso-surfaces of vorticity in z-direction for the viscous Taylor-Green vortex at different time instances

    2.2.2.1 多項(xiàng)式階數(shù)的影響

    為了節(jié)省計(jì)算資源和比較格式的穩(wěn)定性,本算例的計(jì)算均是在欠解析即網(wǎng)格數(shù)量不足條件下進(jìn)行的。我們從p1 開始,逐漸提高多項(xiàng)式階數(shù),直到出現(xiàn)不穩(wěn)定解。如圖6 所示,可以明顯看出:在保證總自由度不變時(shí),提高多項(xiàng)式階數(shù)使得數(shù)值解逼近參考解,這與COX 之前的結(jié)論一致[44]。本文主要關(guān)注的是欠解析情況下格式的非線性穩(wěn)定性問題。當(dāng)多項(xiàng)式階數(shù)提高時(shí),在T=5 時(shí)流場(chǎng)很容易由于混淆誤差而崩潰,這時(shí)就需要具有去混淆效果的分裂形式來增強(qiáng)格式的穩(wěn)定性。

    圖6 總自由度約為 643時(shí),采用Div-LG-g2-Roe 方案的動(dòng)能、動(dòng)能耗散率、擬渦能計(jì)算結(jié)果Fig.6 Computational results of the Div-LG-g2-Roe scheme with the total degree of freedom of 643: (a) kinetic energy;(b) kinetic energy dissipation rate;(c) enstrophy

    2.2.2.2 對(duì)流項(xiàng)形式的影響

    圖6 和圖7 對(duì)比了方案Div-LG-g2-Roe 和方案SplitFe-LG-g2-Roe。在圖6 中,多項(xiàng)式階從p1 到p4 的方案均能穩(wěn)定計(jì)算至T=20,而當(dāng)階數(shù)提高到很高階(p7)時(shí),方案在T=5左右崩潰。而圖7 中,所有多項(xiàng)式階的方案均穩(wěn)定。這一結(jié)果表明,即使是很高階(p8)離散,分裂格式CPR 格式仍能保證數(shù)值穩(wěn)定性。

    圖7 總自由度約為 643時(shí),采用SplitFe-LG-g2-Roe 方案的動(dòng)能、動(dòng)能耗散率、擬渦能計(jì)算結(jié)果Fig.7 Computational results of the SplitFe-LG-g2-Roe scheme with the total degree of freedom of 643: (a) kinetic energy;(b) kinetic energy dissipation rate;(c) enstrophy

    圖8 和圖9 是基于LGL 點(diǎn)的散度形式與分裂形式的對(duì)比,與LG 點(diǎn)所得結(jié)論一致,分裂形式極大地提高了格式的穩(wěn)定性。

    圖8 總自由度約為 643時(shí),采用Div-LGL-g2-Roe 方案的動(dòng)能、動(dòng)能耗散率、擬渦能計(jì)算結(jié)果Fig.8 Computational results of the Div-LGL-g2-Roe scheme with the total degree of freedom of 643: (a) kinetic energy;(b) kinetic energy dissipation rate;(c) enstrophy

    圖9 總自由度約為 643時(shí),采用SplitFe-LGL-g2-Roe 方案的動(dòng)能、動(dòng)能耗散率、擬渦能計(jì)算結(jié)果Fig.9 Computational results of the SplitFe-LGL-g2-Roe scheme with a total degree of freedom of 643: (a) kinetic energy;(b) kinetic energy dissipation rate;(c) enstrophy

    2.2.2.3 解點(diǎn)類型的影響

    如圖6 和圖8 所示,首先對(duì)比散度形式下LG 點(diǎn)和LGL 點(diǎn)的計(jì)算結(jié)果。LG 點(diǎn)在p1~p4 方案下穩(wěn)定,而LGL 點(diǎn)僅能在p1~p2 方案下穩(wěn)定。因此,這一結(jié)果表明,在多項(xiàng)式階不是很高的條件下(p5 以下),基于LG 點(diǎn)比基于LGL 點(diǎn)的散度形式CPR 格式穩(wěn)定性更優(yōu)。

    圖7 和圖9 是分裂形式下LG 點(diǎn)和LGL 點(diǎn)的計(jì)算結(jié)果。兩種解點(diǎn)類型均能在所要求的多項(xiàng)式階下保證穩(wěn)定。從圖7(c)和圖9(c)可以看出,LG 點(diǎn)相較于LGL 點(diǎn)的優(yōu)勢(shì)在于:其計(jì)算精度更高,尤其是多項(xiàng)式階為p7 時(shí),LG 點(diǎn)的擬渦能與參考解吻合得更好。

    2.3 激波捕捉能力測(cè)試

    本小節(jié)考慮二維無(wú)黏Kelvin-Helmholtz 不穩(wěn)定性問題[45]。該算例對(duì)于CPR 方法來說很具有挑戰(zhàn)性,因?yàn)樗瑖?yán)重欠解析的渦結(jié)構(gòu),馬赫數(shù)約等于0.6。本測(cè)試的目的是將結(jié)合子單元限制技術(shù)的分裂形式CPR 格式應(yīng)用于求解欠解析流動(dòng)。初始條件由下式給出:

    其中B=tanh(15y+7.5)-tanh(15y-7.5) 。計(jì)算域[-1,1]2且均為周期邊界。我們使用 64×64個(gè)四邊形單元對(duì)計(jì)算域進(jìn)行均勻劃分,多項(xiàng)式階數(shù)為p7,計(jì)算時(shí)間T=10。

    在圖10 中,我們分別繪制了T=3.7 和T=10時(shí)分裂形式結(jié)合子單元限制求解Kelvin-Helmholtz 不穩(wěn)定性問題的密度等值線。在使用相同激波偵測(cè)器和相同自由度的情況下,相比于Ramírez 等[20]提出的基于LGL 點(diǎn)、結(jié)合子單元限制技術(shù)的DGSEM 格式的計(jì)算結(jié)果,這種新策略能分辨更多的小尺度特征且振蕩更少。圖11 分別給出了兩個(gè)時(shí)刻下的問題單元分布。

    圖10 不同時(shí)刻下Kelvin Helmholtz 不穩(wěn)定性模擬的密度云圖Fig.10 Density contours for the Kelvin Helmholtz instability simulation at different time instances

    圖11 不同時(shí)刻下Kelvin Helmholtz 不穩(wěn)定性模擬的問題單元分布Fig.11 Distribution of problematic cells for the Kelvin Helmholtz instability simulation at different time instances

    3 結(jié)論

    本文在CPR 格式的基礎(chǔ)上,比較了采用不同解點(diǎn)類型的分裂形式對(duì)求解欠解析流動(dòng)時(shí)的穩(wěn)定性和精度的影響,主要有以下結(jié)論:

    1)Sod 激波管和等熵渦算例結(jié)果表明,以LG 點(diǎn)作為解點(diǎn)時(shí),通過邊界通量修正的Fe 分裂形式滿足離散守恒律。

    2)精度測(cè)試結(jié)果表明,使用LG 點(diǎn)的散度形式和分裂形式的L1、L∞誤差基本相同;在相同散度形式或分裂形式下,使用LG 點(diǎn)的誤差明顯小于LGL 點(diǎn)(約一個(gè)量級(jí))。

    3)二維欠解析旋渦流動(dòng)算例用來考察數(shù)值通量的特性和格式的穩(wěn)定性。LLF 和HLL 通量在出口邊界處會(huì)產(chǎn)生偽反射和小尺度噪聲,導(dǎo)致散度形式CPR 格式在計(jì)算時(shí)崩潰,而Fe 分裂形式能夠保證穩(wěn)定計(jì)算。Roe 通量則能完全抑制這種偽反射。此外,使用LGL 點(diǎn)時(shí)放大了出口處產(chǎn)生的偽反射,而LG 點(diǎn)對(duì)流場(chǎng)中的微小結(jié)構(gòu)捕捉得更清楚。

    4)三維黏性Taylor-Green 渦算例結(jié)果表明,無(wú)論采用LG 點(diǎn)或LGL 點(diǎn),即使是很高階(p7)的離散,分裂形式CPR 格式對(duì)于保證計(jì)算穩(wěn)定依然是有效的,且LG 點(diǎn)精度更高,在p7 下與參考解匹配得更好。如果統(tǒng)一使用散度形式,在多項(xiàng)式階不是很高的條件下(p5 以下),采用LG 點(diǎn)的格式穩(wěn)定性優(yōu)于LGL 點(diǎn)。

    5)二維無(wú)黏Kelvin-Helmholtz 不穩(wěn)定性算例考察了分裂形式結(jié)合子單元限制技術(shù)的激波捕捉能力,相比于基于LGL 點(diǎn)的間斷譜元法子單元限制策略,前者在分辨率和穩(wěn)定性方面均有優(yōu)勢(shì)。

    總的來說,本文將基于LG 點(diǎn)的分裂形式CPR 格式用于求解欠解析流動(dòng)問題,不僅能提高格式的穩(wěn)定性,而且相比于LGL 點(diǎn),該方法對(duì)流場(chǎng)的分辨率更高;在結(jié)合子單元限制策略后,其對(duì)含激波的欠解析流動(dòng)也能得到很好的解。本文基于CPR 方法、采用“分裂+LG”策略時(shí)需要采用邊界通量修正保證守恒性,如果想將這一策略應(yīng)用于其他高階譜元方法,在如何滿足格式的守恒性方面有待進(jìn)一步探索。

    猜你喜歡
    散度激波高階
    帶勢(shì)加權(quán)散度形式的Grushin型退化橢圓算子的Dirichlet特征值的上下界
    有限圖上高階Yamabe型方程的非平凡解
    高階各向異性Cahn-Hilliard-Navier-Stokes系統(tǒng)的弱解
    滾動(dòng)軸承壽命高階計(jì)算與應(yīng)用
    哈爾濱軸承(2020年1期)2020-11-03 09:16:02
    一種基于聚類分析的二維激波模式識(shí)別算法
    基于HIFiRE-2超燃發(fā)動(dòng)機(jī)內(nèi)流道的激波邊界層干擾分析
    具有部分BMO系數(shù)的非散度型拋物方程的Lorentz估計(jì)
    斜激波入射V形鈍前緣溢流口激波干擾研究
    H型群上一類散度形算子的特征值估計(jì)
    適于可壓縮多尺度流動(dòng)的緊致型激波捕捉格式
    亚洲国产精品一区二区三区在线| 女的被弄到高潮叫床怎么办| 国产成人精品婷婷| 午夜福利网站1000一区二区三区| 青春草视频在线免费观看| 美女国产高潮福利片在线看| 色哟哟·www| 18在线观看网站| 九草在线视频观看| 亚洲av成人精品一区久久| 51国产日韩欧美| 美女xxoo啪啪120秒动态图| 五月伊人婷婷丁香| 你懂的网址亚洲精品在线观看| 秋霞伦理黄片| 国产精品一二三区在线看| 男男h啪啪无遮挡| 国产精品国产三级国产av玫瑰| 日韩人妻高清精品专区| 精品少妇黑人巨大在线播放| 少妇精品久久久久久久| 日韩 亚洲 欧美在线| 人人妻人人添人人爽欧美一区卜| 日韩伦理黄色片| 久久国产精品大桥未久av| 男人添女人高潮全过程视频| 日韩一区二区三区影片| 丰满少妇做爰视频| 桃花免费在线播放| 男女边摸边吃奶| 国产白丝娇喘喷水9色精品| 国产日韩欧美在线精品| 一级黄片播放器| av卡一久久| 精品久久久久久久久av| 久久99蜜桃精品久久| 久久女婷五月综合色啪小说| 一区二区日韩欧美中文字幕 | 国产精品国产三级国产av玫瑰| 久久午夜综合久久蜜桃| 国产免费一区二区三区四区乱码| 性色av一级| 一个人免费看片子| 精品少妇黑人巨大在线播放| 亚洲丝袜综合中文字幕| 18禁观看日本| 男人爽女人下面视频在线观看| 18禁在线无遮挡免费观看视频| 国产国拍精品亚洲av在线观看| 人人妻人人添人人爽欧美一区卜| 免费不卡的大黄色大毛片视频在线观看| 亚洲av国产av综合av卡| 亚洲av二区三区四区| 久久精品熟女亚洲av麻豆精品| 在线亚洲精品国产二区图片欧美 | 成年人免费黄色播放视频| tube8黄色片| 国产精品蜜桃在线观看| 老司机亚洲免费影院| 综合色丁香网| xxx大片免费视频| 一区在线观看完整版| 中文字幕人妻丝袜制服| 欧美激情国产日韩精品一区| 久久久国产欧美日韩av| 久久精品久久久久久久性| 美女cb高潮喷水在线观看| 中国国产av一级| 久久久久精品久久久久真实原创| 大话2 男鬼变身卡| 久久久久国产网址| 日韩精品免费视频一区二区三区 | 精品一区在线观看国产| 草草在线视频免费看| 视频在线观看一区二区三区| 亚洲人成77777在线视频| 热99久久久久精品小说推荐| 国产日韩欧美亚洲二区| 美女视频免费永久观看网站| 在现免费观看毛片| 国产免费视频播放在线视频| 91久久精品国产一区二区三区| 五月玫瑰六月丁香| 国产国语露脸激情在线看| freevideosex欧美| 久久人妻熟女aⅴ| 日韩一本色道免费dvd| 永久免费av网站大全| 人妻一区二区av| 久久精品熟女亚洲av麻豆精品| 亚洲国产精品国产精品| 飞空精品影院首页| 亚洲精品国产av蜜桃| 欧美精品亚洲一区二区| 美女内射精品一级片tv| 高清不卡的av网站| 国产日韩一区二区三区精品不卡 | 亚洲精品av麻豆狂野| 男女国产视频网站| 日本vs欧美在线观看视频| 一区在线观看完整版| 国产永久视频网站| 久久久久久久久久久免费av| av视频免费观看在线观看| 新久久久久国产一级毛片| 另类精品久久| 一本大道久久a久久精品| 久久国产亚洲av麻豆专区| 18禁裸乳无遮挡动漫免费视频| 最新的欧美精品一区二区| 欧美精品人与动牲交sv欧美| 国产成人精品在线电影| 免费大片黄手机在线观看| 久热这里只有精品99| 亚洲国产成人一精品久久久| 国产精品一二三区在线看| 久久精品久久久久久久性| 熟女人妻精品中文字幕| 另类精品久久| 美女中出高潮动态图| 国产精品.久久久| 久久亚洲国产成人精品v| 中文字幕人妻熟人妻熟丝袜美| 在线精品无人区一区二区三| 免费大片黄手机在线观看| 五月开心婷婷网| 日韩电影二区| 欧美一级a爱片免费观看看| 99久国产av精品国产电影| 男女高潮啪啪啪动态图| 五月伊人婷婷丁香| 天堂中文最新版在线下载| 国产免费又黄又爽又色| 亚洲精品国产av蜜桃| 丝袜喷水一区| 在线观看免费视频网站a站| 人妻一区二区av| 欧美+日韩+精品| 亚洲经典国产精华液单| 亚洲av欧美aⅴ国产| 精品99又大又爽又粗少妇毛片| 黑人巨大精品欧美一区二区蜜桃 | 人妻少妇偷人精品九色| 国内精品宾馆在线| 下体分泌物呈黄色| 麻豆精品久久久久久蜜桃| 下体分泌物呈黄色| 热re99久久精品国产66热6| 免费黄网站久久成人精品| 亚洲国产精品国产精品| 亚洲第一av免费看| 亚洲精品一二三| 久久久久久久久久久免费av| 亚洲第一区二区三区不卡| 亚洲国产精品一区二区三区在线| 日韩av不卡免费在线播放| 91aial.com中文字幕在线观看| 日本wwww免费看| 自拍欧美九色日韩亚洲蝌蚪91| 国产黄频视频在线观看| 国产成人精品久久久久久| 熟女电影av网| 国产成人午夜福利电影在线观看| 久久精品国产a三级三级三级| 熟妇人妻不卡中文字幕| 丝袜脚勾引网站| 丝袜脚勾引网站| 欧美日韩av久久| 热re99久久精品国产66热6| av线在线观看网站| 男人添女人高潮全过程视频| 在线精品无人区一区二区三| 人体艺术视频欧美日本| 久久热精品热| 国产一区二区在线观看日韩| 欧美另类一区| 国产成人a∨麻豆精品| 大片免费播放器 马上看| 高清欧美精品videossex| 久久人人爽人人片av| 久久久久久久久久久丰满| 日韩视频在线欧美| 午夜免费观看性视频| 国产白丝娇喘喷水9色精品| 久久99蜜桃精品久久| 777米奇影视久久| 国产男女内射视频| 欧美少妇被猛烈插入视频| 亚洲内射少妇av| 国产精品一区二区在线不卡| 最近中文字幕高清免费大全6| av卡一久久| 九色亚洲精品在线播放| 欧美丝袜亚洲另类| 成年av动漫网址| 午夜影院在线不卡| 女性生殖器流出的白浆| 亚洲精品aⅴ在线观看| 日韩不卡一区二区三区视频在线| 欧美国产精品一级二级三级| 国产高清有码在线观看视频| 国产在视频线精品| 亚洲精品成人av观看孕妇| 久久精品国产亚洲av天美| 夜夜爽夜夜爽视频| 亚洲精品美女久久av网站| 永久网站在线| 99热网站在线观看| 国国产精品蜜臀av免费| 中国国产av一级| 夜夜骑夜夜射夜夜干| 精品一区二区三区视频在线| 丝袜脚勾引网站| 99热这里只有是精品在线观看| 极品少妇高潮喷水抽搐| 国产成人精品一,二区| 秋霞伦理黄片| 夫妻性生交免费视频一级片| 婷婷色麻豆天堂久久| 亚洲伊人久久精品综合| 亚洲欧洲日产国产| 国产免费视频播放在线视频| 简卡轻食公司| 午夜精品国产一区二区电影| 亚洲精品一二三| 飞空精品影院首页| 国产成人aa在线观看| 一级二级三级毛片免费看| 国国产精品蜜臀av免费| 亚洲国产欧美在线一区| 色婷婷av一区二区三区视频| 精品久久国产蜜桃| 大码成人一级视频| 中文字幕av电影在线播放| 人妻少妇偷人精品九色| 国产一级毛片在线| av专区在线播放| 全区人妻精品视频| 色婷婷av一区二区三区视频| 亚洲欧美日韩卡通动漫| 18禁观看日本| 女人久久www免费人成看片| 免费高清在线观看视频在线观看| 日本av手机在线免费观看| 九九爱精品视频在线观看| 九九在线视频观看精品| 午夜激情久久久久久久| 免费大片18禁| 老司机亚洲免费影院| 欧美精品高潮呻吟av久久| 三级国产精品欧美在线观看| 国产精品三级大全| 最后的刺客免费高清国语| 中文字幕亚洲精品专区| 免费看不卡的av| 在线看a的网站| av免费在线看不卡| 国产综合精华液| 国产 精品1| 18禁动态无遮挡网站| 久久精品国产亚洲av天美| 国产不卡av网站在线观看| 人成视频在线观看免费观看| 18禁裸乳无遮挡动漫免费视频| 国产片内射在线| 国产精品99久久99久久久不卡 | 国产 精品1| 热99久久久久精品小说推荐| 母亲3免费完整高清在线观看 | 午夜老司机福利剧场| 亚洲精品成人av观看孕妇| 日日啪夜夜爽| 久久av网站| av.在线天堂| 亚洲美女视频黄频| 丰满迷人的少妇在线观看| 春色校园在线视频观看| 啦啦啦啦在线视频资源| 另类亚洲欧美激情| 国产免费福利视频在线观看| 亚洲综合色惰| 国产精品久久久久久av不卡| 狂野欧美白嫩少妇大欣赏| 91成人精品电影| 91精品三级在线观看| 久久久午夜欧美精品| 久久久精品94久久精品| 曰老女人黄片| 一级毛片aaaaaa免费看小| 一级毛片黄色毛片免费观看视频| 中国三级夫妇交换| av不卡在线播放| 欧美亚洲日本最大视频资源| 丝袜在线中文字幕| 国产又色又爽无遮挡免| 男女边摸边吃奶| 激情五月婷婷亚洲| 精品99又大又爽又粗少妇毛片| 美女脱内裤让男人舔精品视频| 日韩伦理黄色片| 亚洲精品久久午夜乱码| 成人二区视频| 精品人妻一区二区三区麻豆| 欧美激情 高清一区二区三区| 伊人亚洲综合成人网| 亚洲综合色网址| 中文乱码字字幕精品一区二区三区| 亚洲av成人精品一区久久| 国产av码专区亚洲av| 26uuu在线亚洲综合色| 久久99热这里只频精品6学生| 天堂俺去俺来也www色官网| kizo精华| 久久99蜜桃精品久久| 日本wwww免费看| 各种免费的搞黄视频| 激情五月婷婷亚洲| 在现免费观看毛片| 99热这里只有是精品在线观看| 久久狼人影院| 亚洲在久久综合| 免费观看在线日韩| 一区二区三区精品91| 欧美性感艳星| 51国产日韩欧美| 精品人妻在线不人妻| 亚洲精品亚洲一区二区| 亚洲精品国产色婷婷电影| 午夜激情av网站| 亚洲精品久久成人aⅴ小说 | 91aial.com中文字幕在线观看| 丰满乱子伦码专区| 国产精品99久久99久久久不卡 | 亚洲五月色婷婷综合| 欧美一级a爱片免费观看看| 寂寞人妻少妇视频99o| 美女国产视频在线观看| 久久97久久精品| 妹子高潮喷水视频| 国产精品国产av在线观看| 久久99一区二区三区| 国产成人免费观看mmmm| 国产成人a∨麻豆精品| 欧美 日韩 精品 国产| 午夜福利,免费看| 亚洲一区二区三区欧美精品| 2018国产大陆天天弄谢| 欧美 亚洲 国产 日韩一| 亚洲不卡免费看| 看十八女毛片水多多多| 国产成人91sexporn| 午夜老司机福利剧场| av黄色大香蕉| 免费高清在线观看日韩| 精品久久蜜臀av无| 亚洲色图综合在线观看| 九九爱精品视频在线观看| 成人亚洲欧美一区二区av| 桃花免费在线播放| 一本久久精品| 大又大粗又爽又黄少妇毛片口| 你懂的网址亚洲精品在线观看| 免费不卡的大黄色大毛片视频在线观看| 一本大道久久a久久精品| 大又大粗又爽又黄少妇毛片口| 久久精品久久久久久噜噜老黄| 最近中文字幕高清免费大全6| 久热久热在线精品观看| 亚洲久久久国产精品| 黄片无遮挡物在线观看| 男女国产视频网站| 国产精品久久久久久精品古装| 人人妻人人爽人人添夜夜欢视频| 国产精品偷伦视频观看了| 亚洲情色 制服丝袜| 亚洲精品成人av观看孕妇| 国产精品不卡视频一区二区| 国产一区二区三区av在线| 亚洲精品av麻豆狂野| 熟妇人妻不卡中文字幕| 色网站视频免费| 99久久中文字幕三级久久日本| 五月天丁香电影| 18+在线观看网站| 国产一区二区三区av在线| 男男h啪啪无遮挡| 高清在线视频一区二区三区| 欧美三级亚洲精品| 哪个播放器可以免费观看大片| 日日摸夜夜添夜夜爱| 久久精品夜色国产| 99re6热这里在线精品视频| 亚洲av欧美aⅴ国产| 国产片内射在线| 在线 av 中文字幕| 久热这里只有精品99| 亚洲欧美精品自产自拍| 大码成人一级视频| 欧美xxxx性猛交bbbb| 久久狼人影院| 大香蕉久久成人网| 久久 成人 亚洲| 涩涩av久久男人的天堂| 人成视频在线观看免费观看| 亚洲av电影在线观看一区二区三区| 午夜精品国产一区二区电影| 国产成人精品久久久久久| 精品一区在线观看国产| 午夜精品国产一区二区电影| 久久国产精品大桥未久av| 精品午夜福利在线看| 欧美日韩国产mv在线观看视频| 晚上一个人看的免费电影| 五月开心婷婷网| 97超视频在线观看视频| 免费黄色在线免费观看| 考比视频在线观看| 欧美亚洲 丝袜 人妻 在线| 9色porny在线观看| 婷婷成人精品国产| 亚洲成色77777| 爱豆传媒免费全集在线观看| 欧美日韩av久久| 亚洲无线观看免费| 极品人妻少妇av视频| 韩国高清视频一区二区三区| 高清午夜精品一区二区三区| 天美传媒精品一区二区| 久久久久久久大尺度免费视频| 一级毛片黄色毛片免费观看视频| 国产男女内射视频| 亚洲成人一二三区av| 欧美日韩一区二区视频在线观看视频在线| 日韩三级伦理在线观看| 人妻系列 视频| 免费高清在线观看日韩| 欧美精品国产亚洲| 国产精品一二三区在线看| 久久精品人人爽人人爽视色| √禁漫天堂资源中文www| 国产精品蜜桃在线观看| 伊人亚洲综合成人网| 少妇猛男粗大的猛烈进出视频| 一级毛片aaaaaa免费看小| 免费观看a级毛片全部| 嫩草影院入口| 成人综合一区亚洲| 街头女战士在线观看网站| 99九九在线精品视频| 精品人妻在线不人妻| 又大又黄又爽视频免费| 建设人人有责人人尽责人人享有的| 国产精品久久久久久久电影| 搡老乐熟女国产| 精品亚洲成a人片在线观看| 一区二区三区精品91| 国产成人aa在线观看| 在线看a的网站| 99九九在线精品视频| 国语对白做爰xxxⅹ性视频网站| 中文字幕免费在线视频6| 大话2 男鬼变身卡| 色视频在线一区二区三区| 亚洲美女搞黄在线观看| 亚洲精品aⅴ在线观看| 国产成人免费无遮挡视频| 男人添女人高潮全过程视频| 午夜av观看不卡| 大香蕉97超碰在线| a级片在线免费高清观看视频| 视频区图区小说| a 毛片基地| 亚洲久久久国产精品| 亚洲一级一片aⅴ在线观看| 国产乱来视频区| 国产熟女欧美一区二区| 黑人欧美特级aaaaaa片| 国产色爽女视频免费观看| 欧美日韩视频精品一区| 2021少妇久久久久久久久久久| 麻豆成人av视频| 日韩三级伦理在线观看| 免费看av在线观看网站| 国产免费福利视频在线观看| 久久99热这里只频精品6学生| 精品一品国产午夜福利视频| av福利片在线| 欧美日韩av久久| 性高湖久久久久久久久免费观看| 有码 亚洲区| 欧美成人午夜免费资源| 亚洲精品视频女| 99久国产av精品国产电影| 国产成人freesex在线| 国产男女内射视频| 亚洲不卡免费看| 在线观看www视频免费| 99精国产麻豆久久婷婷| 久久人人爽av亚洲精品天堂| 久久婷婷青草| 日韩制服骚丝袜av| 久久99热6这里只有精品| 99视频精品全部免费 在线| 免费大片18禁| 男人添女人高潮全过程视频| 免费播放大片免费观看视频在线观看| 狂野欧美白嫩少妇大欣赏| 国产av精品麻豆| 精品少妇久久久久久888优播| 成人综合一区亚洲| 人人妻人人澡人人看| 国产乱来视频区| 人妻夜夜爽99麻豆av| 精品少妇内射三级| 国产毛片在线视频| 777米奇影视久久| 亚洲图色成人| 国产精品一国产av| 97精品久久久久久久久久精品| 亚洲人成77777在线视频| 久久人人爽av亚洲精品天堂| 午夜视频国产福利| 国产精品国产三级国产av玫瑰| 飞空精品影院首页| 亚洲精品视频女| 日韩亚洲欧美综合| 看十八女毛片水多多多| 18禁在线播放成人免费| 久久婷婷青草| 爱豆传媒免费全集在线观看| 久久久久久人妻| av播播在线观看一区| 下体分泌物呈黄色| 国产有黄有色有爽视频| 久久韩国三级中文字幕| 人妻一区二区av| 精品一区在线观看国产| 国产成人精品一,二区| 最新中文字幕久久久久| 久久99热6这里只有精品| 天天躁夜夜躁狠狠久久av| 久久国内精品自在自线图片| 亚洲经典国产精华液单| 国产黄片视频在线免费观看| 中国美白少妇内射xxxbb| 国国产精品蜜臀av免费| 黄色欧美视频在线观看| 日韩大片免费观看网站| 久久精品国产亚洲av天美| 日本午夜av视频| 久久精品久久久久久久性| 亚洲精品乱码久久久久久按摩| 我要看黄色一级片免费的| 亚洲av二区三区四区| 国精品久久久久久国模美| 在线播放无遮挡| 97在线视频观看| 欧美一级a爱片免费观看看| 国产成人91sexporn| 熟妇人妻不卡中文字幕| 日韩不卡一区二区三区视频在线| 亚洲精品456在线播放app| 2022亚洲国产成人精品| 国产黄色视频一区二区在线观看| 人妻夜夜爽99麻豆av| 美女cb高潮喷水在线观看| 最近中文字幕2019免费版| 在线观看免费高清a一片| 成人免费观看视频高清| av国产精品久久久久影院| 在线精品无人区一区二区三| 美女xxoo啪啪120秒动态图| 欧美老熟妇乱子伦牲交| 99九九线精品视频在线观看视频| 国产一区二区三区综合在线观看 | 一区二区av电影网| 午夜精品国产一区二区电影| 精品国产乱码久久久久久小说| 亚洲国产欧美日韩在线播放| 欧美日韩在线观看h| 另类精品久久| 黄片无遮挡物在线观看| 国产视频内射| 91国产中文字幕| 妹子高潮喷水视频| 制服诱惑二区| 婷婷成人精品国产| 少妇被粗大猛烈的视频| 亚洲成人一二三区av| 午夜av观看不卡| 国产精品久久久久久久久免| av免费观看日本| 美女视频免费永久观看网站| 免费观看a级毛片全部| 欧美 日韩 精品 国产| 99热国产这里只有精品6| 九九爱精品视频在线观看| 亚洲美女视频黄频| 18禁观看日本| 国产一区二区在线观看av| 新久久久久国产一级毛片| 一本色道久久久久久精品综合| 欧美精品高潮呻吟av久久| 国产男人的电影天堂91| 免费人成在线观看视频色| 男女国产视频网站| 日韩成人伦理影院| 人人妻人人添人人爽欧美一区卜| 一个人看视频在线观看www免费| 国产黄片视频在线免费观看| 国产在线免费精品| 嫩草影院入口| 肉色欧美久久久久久久蜜桃| 热99国产精品久久久久久7| 中文天堂在线官网|