• <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)的緊致型激波捕捉格式
    五月伊人婷婷丁香| 欧美性猛交黑人性爽| 亚洲色图av天堂| 嫩草影院入口| 日本爱情动作片www.在线观看| 国产亚洲精品久久久久久毛片| 麻豆成人午夜福利视频| 高清日韩中文字幕在线| 国产一级毛片七仙女欲春2| 蜜桃亚洲精品一区二区三区| 嫩草影院精品99| 麻豆一二三区av精品| 日本免费a在线| av福利片在线观看| 1024手机看黄色片| 亚洲激情五月婷婷啪啪| 蜜桃久久精品国产亚洲av| 天美传媒精品一区二区| 亚洲五月天丁香| 国产黄片美女视频| 麻豆乱淫一区二区| 一区二区三区四区激情视频 | 一级毛片久久久久久久久女| 看非洲黑人一级黄片| 人人妻人人澡欧美一区二区| av在线观看视频网站免费| 久久99精品国语久久久| 久久精品人妻少妇| 全区人妻精品视频| 日韩成人伦理影院| 亚洲精品乱码久久久v下载方式| 国产真实乱freesex| 18禁裸乳无遮挡免费网站照片| 日日干狠狠操夜夜爽| 99九九线精品视频在线观看视频| 淫秽高清视频在线观看| 国产蜜桃级精品一区二区三区| 日本与韩国留学比较| 久久久午夜欧美精品| 蜜臀久久99精品久久宅男| 欧美三级亚洲精品| .国产精品久久| 国语自产精品视频在线第100页| 国产 一区 欧美 日韩| 久久精品国产99精品国产亚洲性色| 99久国产av精品| 国产日韩欧美在线精品| 人体艺术视频欧美日本| 亚洲精品乱码久久久久久按摩| 国产精品人妻久久久影院| av在线天堂中文字幕| 色播亚洲综合网| 午夜激情福利司机影院| 亚洲丝袜综合中文字幕| 日本黄色视频三级网站网址| 18禁黄网站禁片免费观看直播| 日韩三级伦理在线观看| 夜夜爽天天搞| 久久精品久久久久久久性| 亚洲精华国产精华液的使用体验 | 久久人人爽人人片av| 97超碰精品成人国产| 人人妻人人澡欧美一区二区| 久久韩国三级中文字幕| 国产亚洲av片在线观看秒播厂 | 成人美女网站在线观看视频| 日韩人妻高清精品专区| 国产男人的电影天堂91| 91久久精品国产一区二区成人| 丝袜美腿在线中文| 色吧在线观看| 久久久久久久久久成人| 国产老妇伦熟女老妇高清| 2022亚洲国产成人精品| 赤兔流量卡办理| 又爽又黄无遮挡网站| 日本av手机在线免费观看| 搡女人真爽免费视频火全软件| 亚洲国产精品合色在线| 色综合色国产| 欧美xxxx黑人xx丫x性爽| 国产老妇女一区| 99九九线精品视频在线观看视频| 精品一区二区三区视频在线| 亚洲内射少妇av| 黑人高潮一二区| 爱豆传媒免费全集在线观看| 国产国拍精品亚洲av在线观看| 国产精品久久久久久精品电影| 婷婷色综合大香蕉| 国产精品三级大全| 国产精品女同一区二区软件| 91久久精品国产一区二区成人| 男女视频在线观看网站免费| 性插视频无遮挡在线免费观看| 欧美性感艳星| 午夜老司机福利剧场| 久久午夜亚洲精品久久| 搡老妇女老女人老熟妇| 久久久久久久午夜电影| 色综合站精品国产| 久久精品国产亚洲av香蕉五月| 国产av在哪里看| 身体一侧抽搐| 亚洲图色成人| 成年女人永久免费观看视频| 少妇高潮的动态图| 国产色婷婷99| 亚洲第一电影网av| 久久久久网色| 51国产日韩欧美| 欧美xxxx性猛交bbbb| 99久久成人亚洲精品观看| 国内久久婷婷六月综合欲色啪| 国产亚洲精品久久久com| 免费电影在线观看免费观看| 中文欧美无线码| 人体艺术视频欧美日本| 久久久久久久久中文| 成人欧美大片| 中文资源天堂在线| 久久亚洲国产成人精品v| 赤兔流量卡办理| 久久久久免费精品人妻一区二区| 可以在线观看毛片的网站| av专区在线播放| 中文字幕av在线有码专区| 国产精品电影一区二区三区| 在线国产一区二区在线| 伊人久久精品亚洲午夜| 亚洲av二区三区四区| 黄色日韩在线| 精品久久久久久久久亚洲| 一级二级三级毛片免费看| 国产成人影院久久av| 久久99热6这里只有精品| 免费观看精品视频网站| 国产午夜福利久久久久久| 日韩欧美精品v在线| 欧美性猛交黑人性爽| 国产精华一区二区三区| 久久精品国产亚洲av天美| 色5月婷婷丁香| 中文亚洲av片在线观看爽| 美女大奶头视频| 哪里可以看免费的av片| 少妇猛男粗大的猛烈进出视频 | 最近最新中文字幕大全电影3| 久久亚洲国产成人精品v| 两个人的视频大全免费| a级毛片a级免费在线| 国产精品美女特级片免费视频播放器| 免费人成在线观看视频色| 国产亚洲精品av在线| 国产老妇女一区| 国产高清激情床上av| 我要看日韩黄色一级片| 久久99热6这里只有精品| 99热精品在线国产| 欧美区成人在线视频| 看非洲黑人一级黄片| 舔av片在线| 老司机影院成人| 国产v大片淫在线免费观看| 看片在线看免费视频| 大香蕉久久网| 国产一区二区三区av在线 | 午夜精品一区二区三区免费看| 亚洲国产精品合色在线| 亚洲中文字幕一区二区三区有码在线看| 日本黄大片高清| 搞女人的毛片| a级毛色黄片| 久久综合国产亚洲精品| 欧美又色又爽又黄视频| 亚洲在线观看片| 亚洲七黄色美女视频| 九九在线视频观看精品| 黄色配什么色好看| 91狼人影院| 两性午夜刺激爽爽歪歪视频在线观看| 中文字幕免费在线视频6| 免费大片18禁| 欧美一区二区亚洲| 日韩av在线大香蕉| 能在线免费看毛片的网站| 能在线免费看毛片的网站| 国产日本99.免费观看| 日本黄色片子视频| 天堂网av新在线| 免费人成在线观看视频色| 成人二区视频| 在线免费观看不下载黄p国产| 日本欧美国产在线视频| 99riav亚洲国产免费| 欧美丝袜亚洲另类| 国产精品久久电影中文字幕| 国产日韩欧美在线精品| 人妻久久中文字幕网| 国产一级毛片在线| 亚洲人与动物交配视频| 亚洲欧美日韩高清专用| 禁无遮挡网站| 欧美成人a在线观看| a级毛片a级免费在线| 国产亚洲91精品色在线| 免费电影在线观看免费观看| 菩萨蛮人人尽说江南好唐韦庄 | 欧美+亚洲+日韩+国产| 狂野欧美激情性xxxx在线观看| 国产成人aa在线观看| 欧美变态另类bdsm刘玥| 此物有八面人人有两片| 国产一区亚洲一区在线观看| 在线观看午夜福利视频| 男人狂女人下面高潮的视频| 日韩在线高清观看一区二区三区| h日本视频在线播放| 精品国内亚洲2022精品成人| 欧洲精品卡2卡3卡4卡5卡区| 欧洲精品卡2卡3卡4卡5卡区| 观看免费一级毛片| 欧美+日韩+精品| 国产午夜福利久久久久久| 成人特级黄色片久久久久久久| 在线观看66精品国产| 国语自产精品视频在线第100页| 九草在线视频观看| 人妻制服诱惑在线中文字幕| 久久久久免费精品人妻一区二区| 日韩欧美在线乱码| kizo精华| 国产精品99久久久久久久久| 精品不卡国产一区二区三区| 成人毛片a级毛片在线播放| 国产高清视频在线观看网站| 最近的中文字幕免费完整| 尾随美女入室| 少妇丰满av| 春色校园在线视频观看| 免费一级毛片在线播放高清视频| 一区二区三区高清视频在线| 久久这里只有精品中国| 国产精品久久久久久av不卡| 性色avwww在线观看| 精品无人区乱码1区二区| 国产精品人妻久久久久久| 国产欧美日韩精品一区二区| 又粗又爽又猛毛片免费看| av在线天堂中文字幕| 晚上一个人看的免费电影| 波野结衣二区三区在线| 亚洲国产欧美在线一区| 美女 人体艺术 gogo| 两性午夜刺激爽爽歪歪视频在线观看| 亚洲美女视频黄频| 国产亚洲av嫩草精品影院| 免费人成在线观看视频色| av在线观看视频网站免费| 日日摸夜夜添夜夜爱| 亚洲三级黄色毛片| 国产精品福利在线免费观看| 六月丁香七月| 内地一区二区视频在线| 成人高潮视频无遮挡免费网站| 91aial.com中文字幕在线观看| 免费大片18禁| 午夜福利成人在线免费观看| 在线观看66精品国产| 国产男人的电影天堂91| 狂野欧美白嫩少妇大欣赏| 亚洲欧美日韩东京热| 九九久久精品国产亚洲av麻豆| 午夜免费激情av| 国产色婷婷99| 欧美又色又爽又黄视频| а√天堂www在线а√下载| 日本三级黄在线观看| 亚洲欧美精品自产自拍| 少妇的逼好多水| 午夜激情欧美在线| 大香蕉久久网| 嫩草影院新地址| 国产精品野战在线观看| 国产黄a三级三级三级人| 亚洲国产精品国产精品| 婷婷六月久久综合丁香| 热99re8久久精品国产| 欧美日韩一区二区视频在线观看视频在线 | 小蜜桃在线观看免费完整版高清| 五月伊人婷婷丁香| 久久国内精品自在自线图片| 亚洲性久久影院| 国产一区亚洲一区在线观看| 少妇熟女欧美另类| 一个人看的www免费观看视频| 免费不卡的大黄色大毛片视频在线观看 | 91精品国产九色| 中文亚洲av片在线观看爽| 国产一区二区激情短视频| 免费黄网站久久成人精品| 99久久久亚洲精品蜜臀av| 亚洲av第一区精品v没综合| 黄色欧美视频在线观看| 五月伊人婷婷丁香| 亚洲最大成人手机在线| 国产成人精品一,二区 | 少妇人妻精品综合一区二区 | 99久久精品热视频| 久久久国产成人免费| 天天一区二区日本电影三级| 日韩欧美三级三区| 久久久久性生活片| 99久久九九国产精品国产免费| av又黄又爽大尺度在线免费看 | videossex国产| 国产欧美日韩精品一区二区| 免费观看a级毛片全部| 性插视频无遮挡在线免费观看| 免费人成在线观看视频色| 色噜噜av男人的天堂激情| 麻豆国产av国片精品| 国内精品一区二区在线观看| 亚洲精品色激情综合| 亚洲国产精品久久男人天堂| 男人舔奶头视频| videossex国产| 午夜精品国产一区二区电影 | 少妇熟女aⅴ在线视频| 久久国产乱子免费精品| 亚洲精品影视一区二区三区av| 国产亚洲精品久久久com| 韩国av在线不卡| 亚洲在线自拍视频| 欧美日韩国产亚洲二区| 亚洲成人精品中文字幕电影| 国产黄色小视频在线观看| 欧美色视频一区免费| 免费看光身美女| 深夜精品福利| 久久草成人影院| 亚洲中文字幕一区二区三区有码在线看| 在线观看美女被高潮喷水网站| 亚洲人成网站高清观看| 免费观看a级毛片全部| 久久久国产成人免费| 大香蕉久久网| 成人特级av手机在线观看| 亚洲国产精品国产精品| 自拍偷自拍亚洲精品老妇| 男插女下体视频免费在线播放| 内射极品少妇av片p| 我的女老师完整版在线观看| 精品欧美国产一区二区三| 午夜免费男女啪啪视频观看| 亚洲欧美日韩无卡精品| 久久久久网色| 老熟妇乱子伦视频在线观看| 精品久久国产蜜桃| 久久人人爽人人片av| 久久人人爽人人片av| 国产片特级美女逼逼视频| 日本熟妇午夜| 国产黄a三级三级三级人| 久久久久久久亚洲中文字幕| 国产高潮美女av| 久久综合国产亚洲精品| 国产色婷婷99| 久久久久久伊人网av| 国产精品久久久久久久久免| 伦理电影大哥的女人| 色噜噜av男人的天堂激情| 日韩一区二区三区影片| 精品少妇黑人巨大在线播放 | 久久国产乱子免费精品| 久久午夜福利片| 免费看光身美女| 成人特级av手机在线观看| 亚洲人成网站高清观看| 久久亚洲精品不卡| 久久精品久久久久久久性| 亚洲第一区二区三区不卡| 卡戴珊不雅视频在线播放| 大又大粗又爽又黄少妇毛片口| 黄片无遮挡物在线观看| 久久久久网色| 亚洲电影在线观看av| 亚洲高清免费不卡视频| 国产精品99久久久久久久久| 精品久久久久久久久av| 高清在线视频一区二区三区 | 一区二区三区高清视频在线| 97热精品久久久久久| 不卡视频在线观看欧美| 久久精品影院6| 高清在线视频一区二区三区 | 国语自产精品视频在线第100页| 国产成人福利小说| 亚洲精品乱码久久久v下载方式| 神马国产精品三级电影在线观看| 国产精品久久视频播放| 久久久久久久久久成人| 国国产精品蜜臀av免费| 天堂影院成人在线观看| 97人妻精品一区二区三区麻豆| 亚洲欧美中文字幕日韩二区| 白带黄色成豆腐渣| 中文字幕熟女人妻在线| 成人美女网站在线观看视频| 精品国产三级普通话版| 国产视频内射| 国产高潮美女av| 能在线免费看毛片的网站| 亚洲欧美日韩东京热| 亚洲无线在线观看| 精品久久久噜噜| 九九久久精品国产亚洲av麻豆| 青春草国产在线视频 | 婷婷色av中文字幕| 99久久九九国产精品国产免费| 变态另类丝袜制服| 级片在线观看| 精华霜和精华液先用哪个| 久久精品国产亚洲网站| 我的老师免费观看完整版| 插阴视频在线观看视频| 九草在线视频观看| 免费大片18禁| 国产精品人妻久久久久久| 九色成人免费人妻av| 最近最新中文字幕大全电影3| 国内揄拍国产精品人妻在线| 一级av片app| 亚洲精品乱码久久久久久按摩| 亚洲欧美日韩高清专用| 在线观看免费视频日本深夜| 我要搜黄色片| 青春草视频在线免费观看| 久久亚洲精品不卡| 啦啦啦啦在线视频资源| 人妻夜夜爽99麻豆av| 18禁裸乳无遮挡免费网站照片| 日韩,欧美,国产一区二区三区 | 免费看美女性在线毛片视频| 国产伦精品一区二区三区四那| 国产成年人精品一区二区| 国产91av在线免费观看| 黄片无遮挡物在线观看| 非洲黑人性xxxx精品又粗又长| 身体一侧抽搐| 久久久久九九精品影院| 少妇高潮的动态图| 观看美女的网站| 91精品一卡2卡3卡4卡| 99在线人妻在线中文字幕| 久久精品国产鲁丝片午夜精品| 国产一区二区在线观看日韩| 女人被狂操c到高潮| 蜜臀久久99精品久久宅男| 久久草成人影院| 嫩草影院新地址| 国产黄片美女视频| 久久久国产成人免费| 最好的美女福利视频网| 日本在线视频免费播放| 精品久久国产蜜桃| 亚洲最大成人中文| 18禁裸乳无遮挡免费网站照片| 一进一出抽搐gif免费好疼| 成年免费大片在线观看| 色综合亚洲欧美另类图片| 青春草视频在线免费观看| 深爱激情五月婷婷| 国产精品一区二区三区四区免费观看| 日韩成人av中文字幕在线观看| 成人永久免费在线观看视频| 国产成人精品婷婷| 国产av一区在线观看免费| 国产极品天堂在线| 精品久久久久久久末码| 97超碰精品成人国产| av免费在线看不卡| 最近最新中文字幕大全电影3| 亚洲精品亚洲一区二区| 丝袜美腿在线中文| 大又大粗又爽又黄少妇毛片口| 亚洲一级一片aⅴ在线观看| 国产91av在线免费观看| 欧美又色又爽又黄视频| 乱系列少妇在线播放| 九九在线视频观看精品| 亚洲七黄色美女视频| 亚洲性久久影院| 成人漫画全彩无遮挡| 大又大粗又爽又黄少妇毛片口| 一本久久精品| 欧美日韩一区二区视频在线观看视频在线 | 在线观看av片永久免费下载| 秋霞在线观看毛片| 日本免费一区二区三区高清不卡| 国产精华一区二区三区| 国产精品1区2区在线观看.| 少妇裸体淫交视频免费看高清| 久久这里有精品视频免费| 国产精品一区二区三区四区久久| 国产成人a区在线观看| 搡女人真爽免费视频火全软件| 欧美高清性xxxxhd video| 久久精品久久久久久噜噜老黄 | 高清在线视频一区二区三区 | 久久久午夜欧美精品| 三级毛片av免费| 国产精品久久久久久亚洲av鲁大| 国产真实乱freesex| 久久久久久九九精品二区国产| 成人av在线播放网站| 22中文网久久字幕| 国产黄色视频一区二区在线观看 | 色噜噜av男人的天堂激情| 蜜桃久久精品国产亚洲av| 99久久精品热视频| 国内久久婷婷六月综合欲色啪| 日韩欧美三级三区| 熟妇人妻久久中文字幕3abv| 国内精品宾馆在线| 精华霜和精华液先用哪个| 午夜精品一区二区三区免费看| 日本免费a在线| 成年女人看的毛片在线观看| 性欧美人与动物交配| 亚州av有码| 久久久久久久久久成人| 男人和女人高潮做爰伦理| 成人三级黄色视频| 亚洲欧美日韩高清在线视频| 狠狠狠狠99中文字幕| 成人毛片60女人毛片免费| 成熟少妇高潮喷水视频| 国产精品三级大全| 国产精品1区2区在线观看.| 中文字幕精品亚洲无线码一区| 久久亚洲精品不卡| 日韩欧美在线乱码| 亚洲va在线va天堂va国产| 自拍偷自拍亚洲精品老妇| 老女人水多毛片| 欧美日韩乱码在线| av在线蜜桃| 成人亚洲精品av一区二区| 国产一区二区激情短视频| 亚洲av二区三区四区| 久久人人爽人人片av| 99久久中文字幕三级久久日本| 国产午夜精品论理片| 日本欧美国产在线视频| 少妇被粗大猛烈的视频| 午夜视频国产福利| 亚洲精品久久国产高清桃花| 村上凉子中文字幕在线| av在线老鸭窝| 狠狠狠狠99中文字幕| 国内少妇人妻偷人精品xxx网站| 亚洲国产精品久久男人天堂| 激情 狠狠 欧美| 你懂的网址亚洲精品在线观看 | 免费看光身美女| 成人av在线播放网站| 男女做爰动态图高潮gif福利片| 尤物成人国产欧美一区二区三区| 亚洲最大成人av| 欧美日韩一区二区视频在线观看视频在线 | 日日摸夜夜添夜夜添av毛片| 午夜福利视频1000在线观看| 五月伊人婷婷丁香| 赤兔流量卡办理| 久久久久九九精品影院| av免费观看日本| 国产伦精品一区二区三区四那| 午夜福利在线观看吧| 亚洲欧美精品专区久久| 亚洲三级黄色毛片| av卡一久久| 99国产极品粉嫩在线观看| 亚洲av成人精品一区久久| 一区二区三区四区激情视频 | 免费黄网站久久成人精品| 激情 狠狠 欧美| 国产极品精品免费视频能看的| 午夜福利成人在线免费观看| 综合色av麻豆| 国产一级毛片七仙女欲春2| 精品99又大又爽又粗少妇毛片| 狂野欧美激情性xxxx在线观看| 女同久久另类99精品国产91| av女优亚洲男人天堂| 久久99蜜桃精品久久| 18禁裸乳无遮挡免费网站照片| 精品久久久久久成人av| 精品国产三级普通话版| 亚洲国产精品国产精品| 亚洲av成人av| 欧美变态另类bdsm刘玥| 免费观看人在逋| 少妇熟女aⅴ在线视频| 欧美在线一区亚洲| 高清午夜精品一区二区三区 | 免费av观看视频| 日日干狠狠操夜夜爽| 欧美日韩综合久久久久久| 高清毛片免费观看视频网站| 免费观看在线日韩| 亚洲人与动物交配视频| 91aial.com中文字幕在线观看| 一卡2卡三卡四卡精品乱码亚洲| 免费av毛片视频|