999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

空腔流動的動量分解及能量輸運特性1)

2022-03-20 15:52:36韓帥斌武從海張樹海
力學學報 2022年2期
關鍵詞:模態

韓帥斌 羅 勇 李 虎 武從海 張樹海

(中國空氣動力研究與發展中心,空氣動力學國家重點實驗室,四川綿陽 621000)

(中國空氣動力研究與發展中心,計算空氣動力研究所,四川綿陽 621000)

引言

空腔是一種非常具有代表性的結構,廣泛應用于航空航天飛行器部件及地面交通工具中,如飛機的起落架艙、航空發動機內部凹槽以及高速列車受電弓腔、汽車天窗等[1].空腔流動結構復雜,存在剪切層和旋渦,對超聲速流動還會產生激波,是一種重要的噪聲源,也是相關工程設計中必須考慮的關鍵問題.

空腔流場中各種結構之間的相互作用會產生豐富的非線性物理過程,包括多模態共振、渦聲的產生和傳播、多波相互作用等.在亞聲速空腔流動中,剪切層、旋渦和聲波等流動結構之間的相互作用占據主導地位;超聲速空腔流動中則存在復雜的激波波系,激波之間、激波與其他流動結構之間相互作用會產生激波噪聲[2-3].空腔噪聲的產生涉及以下動力學過程:起始于空腔前緣的渦擾動在剪切層內不斷增長并向下游對流,對流至空腔后緣點處的渦與后緣點相互作用從而產生向外輻射和向上游傳播的聲波,聲波擾動到達前緣點并進一步誘發新的渦擾動.在以上動力學過程中,腔體內形成擾動反饋增長回路,主流中的動力學能量不斷轉化為聲能,產生的強噪聲持續向外輻射.

在空腔流動中,流體動力學模態(渦模態和熵模態)與聲模態之間的相互作用是空腔自持振蕩的重要過程,二者在近場的耦合作用機制對流動的發展和演化起著關鍵作用[4].在低速流動中,熵模態通常較弱,渦聲模態之間的相互作用構成了流動的主要機制.Kerschen 和Tumin[5]與Alvarez 等[6]基于空腔邊緣散射機制,建立了超聲速和亞聲速空腔聲模態與剪切層內渦層在空腔前后緣相互作用并轉化的理論模型,預測了渦波和聲波的幅值和相位變化;Tang和Rockwell[7]研究了渦與空腔后緣相互作用形式對瞬態壓力場的影響,發現壓力幅值及相位與渦-緣相互作用形式密切相關;萬振華[8]基于渦量和擬渦能分析了空腔內渦結構撞擊空腔后緣點產生脈動壓力的動力學過程;韓帥斌等[9]采用拉格朗日擬序結構分析了空腔流動中渦的生成、對流、撞擊和破裂等動力學過程.隨流動馬赫數增加,流動可壓縮性增強,熵模態擾動也隨之增強.Arya 和De[10]基于聲學擾動方程(acoustic perturbation equation,APE)分析了可壓縮空腔流動中的Lamb 矢量所代表的渦聲源和熵的時空擾動所代表的熵聲源對噪聲產生的貢獻,發現熵聲源對于遠場和近場聲壓均具有一定影響;Liu 和Gaitonde[11]基于預解分析(resolvent analysis)對亞聲速和超聲速空腔流動中渦熵動力學模態和聲模態的激勵響應特性進行了研究.

已有的研究基本上都采用不同的變量分別表征流聲模態,并未對同一流動變量進行流聲模態分解,無法深入理解流聲模態之間的相互作用及不同模態之間的能量轉化機制.準確識別并解耦空腔內的流體動力學模態和聲模態,是深入理解空腔流聲相互作用和能量轉化機制的關鍵.對于任意一個矢量場,亥姆霍茲分解可將其分解為無旋有散的標量梯度場和一個有旋無散的向量旋度場[12-13].針對速度場,前者表征渦模態速度,對應于旋渦結構;后者則是流體可壓縮性的體現,在等熵流中可表征聲模態速度,對應于聲波、激波、膨脹波等可壓縮流動結構.針對速度場的傳統亥姆霍茲分解在矢量氣動聲學[14]、可壓縮湍流[15]中得到了廣泛應用.然而傳統亥姆霍茲分解無法解耦速度場中的熵模態,僅能應用于等熵流動.Doak[16]提出了針對動量的亥姆霍茲分解,并建立了動量勢理論(momentum potential theory,MPT).MPT 將動量分解為渦(有旋無散量)、聲(無旋等熵量)和熵(無旋等壓量)3 類組分,系統分析了各類組分相互作用的能流特性.目前MPT 已被廣泛應用于波包動力學[17]、射流[18-19]、鈍體繞流[20-21]以及邊界層轉捩[22-23]等,揭示了相應流動中不同模態的流動特性.

本文針對來流馬赫數Ma=0.8 的二維亞聲速空腔流動進行數值模擬,獲得高精度流場數據.基于MPT 對空腔流動的動量進行分解,得到動量的渦熵動力學組分和聲組分,研究各組分動量的時空特性,試圖分析并揭示不同組分相關的能量輸運特性.

1 數值及理論方法

1.1 空腔流動高精度數值模擬

在文獻[9,24]中,已經對空腔的流場結構進行了高精度數值模擬,并進行了網格收斂性驗證和計算結果與文獻對比驗證.在數值模擬中,控制方程為二維Navier-Stokes 方程,其具體形式為

式中 ρ,u,p,E,σ,T,μ,Re,Pr,M 分別為密度、速度、壓力、能量、黏性應力張量、溫度、動力學黏性系數、雷諾數、普朗特數和馬赫數.

無量綱過程中的參考長度為空腔深度,速度的參考值為無窮遠處聲速,密度、壓力、溫度等流場變量的參考值采用無窮遠值.計算中,Navier-Stokes方程的對流項采用五階WENO 格式[25-27]離散,黏性項采用6 階中心差分格式離散,時間推進采用3 階TVD Runge-Kutta 方法.

空腔流動數值模擬的計算區域如圖1 所示,其中空腔長深比為2:1,各長度參數為:L=2D,Ll=7D,Lr=20D,Ly=20D,Li=3D,Lo=10D,Lt=10D.空腔來流馬赫數為Ma=0.8,普朗特數為Pr=0.7,雷諾數為Re=2500.在數值模擬氣動聲學問題時,由于聲學擾動是小量,遠場邊界處微弱的聲波反射即可引起較大的數值計算誤差,因此本文采用Bodony[28]的方法,在計算域外緣設置海綿層從而降低邊界上聲波的反射,避免聲場污染.數值計算網格量為:腔外網格810×300,腔內網格240×120,其中固壁邊界處的網格進行了加密處理以捕捉邊界層,遠場網格進行拉伸以降低反射[29].數值模擬所得的流場紋影如圖2 所示,清晰顯示了空腔流動的遠場聲場和近場剪切層及渦結構.對空腔(1.0D,1.0D)處的速度和壓力脈動采樣分析,可得流動的無量綱時間周期為T=3.75[9],相應的無量綱頻率St=fL/U=0.67,與Rossiter 半經驗公式給出的第二模態頻率St2=0.686 和Krishnamurty[30]的實驗結果St=0.656 吻合.網格收斂性驗證和計算結果的準確性驗證可參考文獻[9,24].

圖1 空腔流動數值模擬計算區域Fig.1 The computational domain of open cavity flow for numerical simulation

圖2 空腔流動的數值紋影Fig.2 Numerical schlieren of open cavity flow

1.2 動量勢理論

MPT 由Doak[16]提出,通過對動量應用亥姆霍茲分解,將動量分解為有旋無散的渦動力學量、有散無旋的聲學量和熵動力學量,具體為

其中B0為時間平均的動力學量,B′為脈動渦動力學量,ψ 為脈動標量勢,表征流動的可壓縮特性,包含了聲組分和熵組分.結合連續性方程

可得標量勢滿足泊松方程

基于以上分解,動量可視為時間平均動力學量、渦組分、聲組分和熵組分4 個部分的線性疊加.通過求解式(5)和式(6)中任意兩個泊松方程即可獲得流場中動量的不同組分.由于壓力和密度為流場的原始變量,較易通過求解獲得其時間導數,因此本文通過求解總標量勢泊松方程(5)和聲標量勢泊松方程(6b)進行動量分解,其具體實現步驟如下:

(1)對周期性或準周期性流動,依據流場變量的時間演化獲得其流動周期T,并對一個流動周期[t0,t0+T]內的動量 ρu 進行時間平均,由此獲得B0;

(2)對流場進行時間微分獲得 ?ρ′/?t,求解泊松方程(5)獲得標量勢 ψ ;

(3)對流場進行時間微分獲得 ?p′/?t,并計算當地流體聲速 c2=γp/ρ,求解泊松方程(6b)獲得聲標量勢 ψA;ψT=ψ-ψA

(4)計算熵標量勢 ;B′=ρu-B0+?ψ

(5)計算無散的渦組分 .

在空腔流動的動量分解過程中,泊松方程的求解采用Gauss-Seidel 迭代方法,每次迭代包括交替方向的兩次掃描.初始條件和邊界條件設置可參考文獻[18,20-21],具體為:標量勢 ψ 的初始條件設為ψ0=0;邊界條件設置為:左邊界、上邊界和右邊界采用遠場邊界(由圖3 中不同高度處動量脈動的旋度 ?×(ρu)′的分布可知,?×(ρu)′的幅值從近場到邊界處由 O(1) 衰減至 O(10-4),說明在邊界處動量脈動幾乎不存在有旋組分),可令有旋無散脈動量B′=0[18,21],于是由式(3a)有(ρu)′=-?ψ,因此沿邊界積分即可實現對遠場邊界的賦值,即遠場邊界施加狄利克雷邊界條件;對于物面邊界,設置法向零梯度邊界條件,即諾依曼邊界條件.求解 ψA時,遠場邊界賦值為 ψA=ψ,物面邊界則賦法向零梯度邊界條件.

圖3 不同高度處動量脈動的旋度分布Fig.3 The distribution of the curl of momentum perturbation

在以上動量分解的基礎上,能量守恒方程的時間平均可寫為如下形式

其中 H=cpT+u·u/2 為總焓;α′為加速度脈動,由流動中的渦量擾動、熵空間不均勻性以及黏性應力所貢獻.上式左端項表示動量的渦聲熵各組分所攜帶的脈動總焓的時均流量,右端則是輸運過程中各組分動量與加速度脈動相互作用所形成的源,表示對流體微團做功的功率,影響左端平均能流的產生(正源)或耗散(負源).當流場中不存在某一類動量組分時,其對應的流量及源相應為零,因此基于式(7)可分析各動量組分相關的能量的產生與耗散機制以及能量的輸運與轉化特性.

2 方法驗證

首先驗證泊松方程求解方法.圖4 給出了空腔內不同站位處泊松方程(5)左右端數值結果,兩者吻合很好,且方程兩端殘差的數量級為 O(10-7),數值求解取得了較好的收斂效果.

圖4 不同位置處泊松方程(5)左右端項分布Fig.4 The distribution of LHS and RHS of Eq.(5) at different locations

其次驗證動量分解所得各組分的物理屬性.圖5 給出了空腔上方y=1.2 處的時間平均動量、渦組分動量以及可壓縮組分動量的散度分布.由散度分布可見,時均動量和渦組分動量的散度均幾乎為零,其中在空腔附近[0,5]的流動旺盛區域,由于流動的空間不均勻性,時均動量和渦組分動量的散度均存在非零擾動,但二者最大幅值均為 O(10-3),遠遠小于 ?2ψ 的量級O (10-1),因此可認為采用MPT成功實現了動量的無散動力學量和有散可壓縮量的分離.

圖5 y=1.2 處不同組分的散度分布Fig.5 The divergence of different components of momentum at y=1.2

采用快速傅里葉變換(FFT)對密度、速度、原始動量以及MPT 分解所得各組分動量的頻譜特性進行對比驗證.對位于空腔(1.0D,1.0D)處的各變量進行FFT 分析,結果如圖6 所示,各流場變量以及原始動量的頻譜為離散頻譜,主頻均為St=0.67,其余峰值為主頻的諧頻,說明流動具有強周期性.MPT分解所得渦聲熵組分動量的主頻與原始變量一致,也為St=0.67,進一步證明了采用MPT 分解所得各組分的準確性.

圖6 頻譜分布Fig.6 Frequency spectrum

3 結果與分析

3.1 動量的流聲組分特性

對空腔流動的流場應用MPT 進行動量分解,分析動量的聲熵渦各組分空間分布特性、時間演化特性.圖7 給出了標量勢以及聲熵渦各組分動量在一個流動周期內任一時刻的空間分布.其中聲熵組分為動量的可壓縮部分,其空間分布特性可由標量勢表征.由圖7(a)~ 圖7(c)可知,總標量勢 ψ 所表征的脈動主要分為3 部分:空腔剪切層附近同時包含聲組分和熵組分的周期性脹縮脈動,空腔左上方強烈的聲脈動和空腔右側相對較弱的聲脈動.聲標量勢ψA與總標量勢 ψ 在遠場幾乎一致,說明動量的可壓縮部分在遠場僅包含聲組分,而在近場由于熵組分的存在,聲標量勢與總標量勢在空腔剪切層附近有所差異.熵標量勢 ψT集中分布在空腔剪切層附近及空腔后方尾跡內,并呈現出周期性的壓縮膨脹特征.

圖7 動量不同組分的空間分布Fig.7 Spatial distribution of different components of momentum

動量的聲熵渦組分在x 和y 方向的分量如圖7(d)~ 圖7(i)所示.聲組分動量分布于空腔內以及遠場.空腔內的聲組分動量揭示了空腔自持振蕩的聲學脈動,其波陣面以近似垂直于空腔下壁面的平面波形式在腔體內傳播,強度與腔外聲波為同一量級,但高于腔外聲波,與文獻[8]中空腔內聲壓級高于腔外相一致;空腔外的遠場波陣面分布則揭示了聲波由空腔后緣點和前緣點向外輻射并傳播的過程,其中空腔左側聲波波系由于入口均勻來流產生了明顯的多普勒效應;空腔上方的波系較強,是空腔非定常流動產生的主要噪聲,空腔右側存在較弱波系,呈現為球面波形狀.熵組分動量的空間分布模式與渦組分相似,兩者均主要分布于空腔剪切層及空腔后緣點后的尾跡中,與剪切層內卷起的渦層及空腔右側的主渦密切相關,但是熵組分動量強度比渦組分低1~ 2 個數量級,與聲組分數量級接近.因此從強度看,在空腔流動的近場,渦組分動量為主要部分,主導著剪切層內渦層的發展演化;在遠場則只存在聲組分,揭示遠場噪聲的傳播過程.

由動量的聲熵渦組分在近場的空間分布可知盡管三者均呈現沿流向的周期性正負交替分布,但分布特性仍有差異.其中聲組分以近似垂直于空腔下壁面的平面波形式分布,并在空腔內振蕩傳播,而熵渦組分則集中于剪切層,以波包形式存在,隨主流向下游對流.為了揭示三者的傳播速度,對空腔唇口線上(x ∈(0,2D),y=1D) 的各組分動量在一個流動周期內的時空演化特性進行分析.圖8(a)~ 圖8(c)分別為聲熵渦流向動量的的x-t 分布圖,其中只有聲組分動量呈現向上游傳播的擾動,這是由空腔后緣點產生并向上游傳播的聲波,由于來流速度影響,其傳播速度為為唇口線上流向速度的空間平均值,c∞=1 為單位無量綱聲速;熵組分和渦組分動量則沿來流方向隨著剪切層向下游傳播,其傳播速度均為 |dx/dt|=0.5 ≈Uc,這里 Uc≈L/T=0.53 為剪切層的平均對流速度.

圖8 各組分流向動量的時空分布Fig.8 The spatial-temporal distribution of different components of momentum

3.2 動量流聲組分相關的能量輸運特性

渦聲熵組分動量相關的能量輸運特性可由式(7a)描述,其中左端項依次代表渦聲熵組分動量所攜帶的總焓流量,右端項則是相應的源.對能量方程(7a)的左端項和右端項分別計算并對比,結果如圖9所示,兩者完全吻合,且計算的殘差最大為 O(10-5),遠小于各組分動量所攜帶的總焓流量的大小,證明了計算的準確性.同時由圖9 可知,時均總焓流量分布可劃分為5 個區域:I 區為渦卷起并充分發展的剪切層核心區,總焓流量為正;位于剪切層外側的II 區和位于空腔后緣點左上附近的III 區總焓流量均為負,表明在剪切層內,動量所攜帶的總焓不斷流出,隨著剪切層對空腔后緣點的撞擊,總焓不斷向空腔后緣點附近匯聚,并向剪切層外側傳遞能量.IV 區為空腔內右半側的主渦區,在主渦內,總焓流量呈正負交替分布,主渦內整體呈現能量平衡狀態,因而可保持穩定.V 區為空腔后緣點后的尾跡區,總焓流量沿垂直于流向方向正負交替分布,邊界層內的能量不斷向邊界層外傳遞.

圖9 能量方程左右端項Fig.9 Left hand side (LHS) and right hand side (RHS) of the averaged energy equation

圖10 給出了能量方程左端渦聲熵各組分動量所攜帶的總焓流量與相應源的分布.渦組分源分布主要集中于剪切層,呈現垂直于流向的正負交替的雙層分布,并在后緣點附近形成負源,因此渦組分動量所攜帶的總焓由空腔前緣點起始,跟隨剪切層的運動,從剪切層內不斷輸運至剪切層外側以及空腔后緣點處.空腔內的主渦處交替分布了相對較弱的正負渦組分源,其能量的平衡使得空腔右側主渦保持相對穩定狀態.聲組分動量所攜帶的總焓流量在剪切層內呈現出強度近似的周期性正負分布,空腔后緣點左側區域總焓流量為負值,表明聲組分動量攜帶的總焓流量由后緣點右側不斷流向該區域,并向上游和遠場傳播.聲組分源僅存在于近場,在剪切層內沿流向呈周期性正負分布,且正源顯著強于負源,因此聲組分動量攜帶的總焓從剪切層內不斷流出,以聲能形式向外輻射;同時空腔右側主渦附近也存在著強度相近的正負聲組分源,表明空腔內的聲自持振蕩可以保持能量的穩定狀態.熵組分動量與渦組分動量空間分布相似,均與剪切層密切相關;另一方面熵組分動量與聲組分動量均具備可壓縮特性,因此熵組分動量攜帶的總焓流量分布同時呈現出了類似于渦組分的垂直于流向的層狀分布特性和類似于聲組分的沿流向的周期性正負交替分布特性.熵組分動量攜帶的總焓流量由剪切層下方輸運至剪切層核心區域,剪切層上方則沿流向周期性流入流出,保持平衡狀態.熵組分源主要分布于剪切層附近,且主要為負源,因此熵組分動量攜帶的總焓不斷向剪切層內輸運,并在剪切層內耗散.

圖10 能量方程各項空間分布Fig.10 The spatial distribution of various terms in energy equation

對比3 類動量組分所攜帶的總焓流量的強度可知,渦組分高于聲組分和熵組分至少一個數量級,因而渦組分相關的能量輸運過程在空腔流動中占據主導地位;聲組分和熵組分數量級接近,盡管兩者強度較小,但在空腔內的聲能輸運以及能量耗散過程中均起著重要作用.值得注意的是,渦聲熵組分相關的能量方程左右端分布均有顯著差異,其差值resi=RHSi-LHSi 如圖11 所示,resi>0 表示相應組分的源大于該組分的總焓流量,反之則相反,因此動量的渦聲熵組分所攜帶的總焓流量與相應的源并非一一對應關系,在不同組分之間存在著能量的轉化.為了量化分析各組分之間的能量輸運特性,在圖10(a)中的流動旺盛區域A 內對能量方程(7a)中的各項進行空間積分,獲得該區域內各組分動量相關的能量的凈流出/流入和凈產生/耗散,結果如表1 所示.由左端項的凈流量可知,渦組分動量攜帶動力學能凈流出區域A,輸運至空腔尾跡的邊界層和空腔內的主渦內;聲組分動量攜帶聲能凈流出,向遠場持續輻射聲波;熵組分動量則攜帶能量凈流入至區域A 內.由右端項的凈源可知,渦組分和聲組分源均為正源,不斷產生相應的能量,熵組分源則為負源,耗散流動中的能量.對比渦聲熵各組分的凈源和凈流量可知,渦組分和聲組分能量的凈產生均高于其凈流出,而熵組分的凈耗散高于其凈流入,因此由能量平衡可知,在剪切層內存在著渦聲組分能量向熵組分能量的轉化.

圖11 各組分源與流量殘差Fig.11 The residue for different components of source and energy flux

表1 能量方程7(a)中各組分的空間積分Table 1 The spatial integration of each term in the energy equation 7(a)

4 結論

本文針對來流馬赫數Ma=0.8 的二維亞聲速空腔流動進行數值模擬,獲得其高精度流場數據,采用基于亥姆霍茲分解的動量勢理論,對空腔流動進行動量分解,并分析了聲熵渦各組分動量的物理特性以及能量輸運特性.

(1)數值實現了動量勢理論,獲得并驗證了空腔流動的聲熵渦組分的動量分布,其中渦熵組分動量僅存在于近場,與剪切層的發展演化密切相關,隨主流向下游運動;而聲組分動量則同時存在于整個流場,并由近場傳播至遠場.

(2)渦熵組分動量攜帶的總焓從剪切層內不斷輸運至剪切層外側及空腔后緣點處;熵組分動量攜帶的總焓不斷向剪切層內輸運,并在剪切層內耗散;聲組分動量攜帶的總焓則以聲能形式由后緣點右側不斷向上游和遠場傳播.

(3)渦組分相關的能量輸運過程在空腔流動中占據主導地位,同時剪切層內存在著渦聲組分能量向熵組分能量的轉化和耗散.

猜你喜歡
模態
基于BERT-VGG16的多模態情感分析模型
跨模態通信理論及關鍵技術初探
一種新的基于模態信息的梁結構損傷識別方法
工程與建設(2019年1期)2019-09-03 01:12:12
多跨彈性支撐Timoshenko梁的模態分析
車輛CAE分析中自由模態和約束模態的應用與對比
國內多模態教學研究回顧與展望
基于HHT和Prony算法的電力系統低頻振蕩模態識別
由單個模態構造對稱簡支梁的抗彎剛度
計算物理(2014年2期)2014-03-11 17:01:39
利用源強聲輻射模態識別噪聲源
日版《午夜兇鈴》多模態隱喻的認知研究
電影新作(2014年1期)2014-02-27 09:07:36
主站蜘蛛池模板: 久久久久国产精品免费免费不卡| 最新加勒比隔壁人妻| 亚洲成a人片| 午夜天堂视频| 欧美一区精品| 91麻豆精品视频| 中文字幕久久波多野结衣| 日本AⅤ精品一区二区三区日| 99精品免费欧美成人小视频| 青青操国产视频| 国产一区亚洲一区| 天天躁日日躁狠狠躁中文字幕| 亚洲日本韩在线观看| 国产正在播放| 91视频99| 国产精品手机在线观看你懂的| 久久青草免费91观看| 久久香蕉国产线看精品| 欧美精品v日韩精品v国产精品| 欧美另类精品一区二区三区| 中文字幕不卡免费高清视频| 91免费片| 国产网友愉拍精品| 免费高清毛片| 精品无码国产一区二区三区AV| 亚洲中文字幕久久无码精品A| 美女毛片在线| 91高清在线视频| 美女视频黄又黄又免费高清| 国产成人精品2021欧美日韩| 色欲色欲久久综合网| 男女性午夜福利网站| 另类综合视频| 99re这里只有国产中文精品国产精品 | 综合久久五月天| 毛片免费视频| 日本爱爱精品一区二区| 国产黄网永久免费| 精品亚洲欧美中文字幕在线看| 亚洲精品午夜天堂网页| 好吊色国产欧美日韩免费观看| 国产成人永久免费视频| 欧美精品影院| 国产微拍精品| 国产免费人成视频网| 无码国内精品人妻少妇蜜桃视频| 伊人色综合久久天天| 香蕉在线视频网站| 亚洲午夜久久久精品电影院| 97国产精品视频人人做人人爱| 国产午夜一级毛片| 超级碰免费视频91| 国产后式a一视频| 欧美亚洲香蕉| 亚洲成人在线网| 超清人妻系列无码专区| 精品少妇人妻av无码久久| jijzzizz老师出水喷水喷出| 国产丝袜啪啪| 亚洲专区一区二区在线观看| 亚欧成人无码AV在线播放| 精品乱码久久久久久久| 97se亚洲综合在线天天| 国产玖玖视频| 亚洲国产成人在线| 日韩AV手机在线观看蜜芽| 激情国产精品一区| 国模私拍一区二区| 亚洲高清免费在线观看| 国产免费网址| 久久久久人妻精品一区三寸蜜桃| 国产精品天干天干在线观看 | 超清无码熟妇人妻AV在线绿巨人| 4虎影视国产在线观看精品| 国产偷国产偷在线高清| 天堂久久久久久中文字幕| 直接黄91麻豆网站| 人妻无码中文字幕第一区| 国产精品亚洲а∨天堂免下载| 日本人真淫视频一区二区三区| 波多野结衣中文字幕一区二区| 亚洲色图欧美在线|