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

基于電容層析測量的油氣兩相流動態信號時頻分析方法

2024-03-18 08:59:48徐一李軼馬志揚王海剛
化工進展 2024年2期
關鍵詞:測量信號

徐一,李軼,馬志揚,王海剛

(1 中國科學院工程熱物理研究所,煤炭高效低碳利用全國重點實驗室,北京 100190;2 中國科學院大學,北京 100190;3 清華大學深圳國際研究生院,廣東 深圳 518055)

多相流動廣泛存在于石油開采及油氣混相輸送管線中,其流動復雜性體現在各相非均勻混合且存在相互作用、流型變化復雜、特征參數多等特點。對這一非線性、高度隨機性的瞬態變化過程實現實時在線的準確測量具有重要的理論意義和工程應用價值。近年來,研究者們致力于發展非侵入性、流型適應性強兼顧測量準確性、體積及成本的在線多相流量計,以期彌補依靠氣液分離器和單相流量計的傳統測量系統在實際工業應用中的時間延遲、維護困難等不足,從而適應油氣開發的不同場景[1]。

電容層析成像(electrical capacitance tomography,ECT)通過重構感測截面的介電常數分布實現對非導電介質的可視化測量,且具有不產生壓降、安全無輻射、響應速度快的優點,近年來被嘗試應用于油氣兩相流動過程的實時監測[2-3]。然而,由于存在反演誤差以及難以準確獲得流體的真實數據,現有基于ECT 的成像方法難以有效提取多相流動態流動關鍵信息。層析成像的相關研究多以靜態仿真和實驗測量為主,缺乏動態數據信息,而基于數字孿生(digital twin, DT)技術輔助ECT成像的方法[4]有效拓展了現有研究手段,能夠高效獲取更加接近真實流場的動態測量數據。在數據分析方面,現有研究多側重于優化ECT 的圖像重建質量[5-6],而對于原始電容信號內耦合的流場動態信息關注較少。在可視化基礎上進一步提取層析測量信號中與多相流型和流動特性相關的特征信息,將有助于提升測量精度和對復雜多相流動機理的理解。

多相流測量信號是由表征不同流體動力學現象的多個分量共同作用的結果。常見的氣液兩相流壓力信號概括反映了流速、相含率、湍流度以及氣泡、液塞等不同尺度的局部和全局流動行為。應用不同的信號分析方法,主要包括時域、頻域、時-頻域和非線性分析[7],可以從不同視角提取信號中蘊含的特征信息。小波變換是一種具有代表性的時-頻域分析方法,能夠在高頻范圍使用高時間分辨率、在低頻范圍使用高頻率分辨率對信號展開靈活解析,適用于流態信息集中在低頻段、以高度隨機性的瞬態變化為特點的多相流動。Wang等[8]通過壓力信號的連續小波分析有效提取了氣固流態化過程中的流場局部特征。Nguyen等[9]對多通道阻抗空隙率計的測量信號應用連續小波分析,通過提取有效局部小波能量和有效尺度,建立水氣垂直兩相流的流型識別標準。李楓[10]基于電阻層析成像技術應用小波分析和人工神經網絡模型,開發了針對氣水兩相水平管流的在線流型識別軟件。陳飛和孫斌[11]通過多孔孔板測量,基于自適應線性調頻小波對氣液兩相流的差壓信號進行時頻聯合分析,并提取與流型變化相關的特征值。

針對動態研究中多相流動狀態在時間、空間維度上均發生顯著變化的復雜特性,以及由于感測信號中多信息融合難以通過單一變量或靜態測試中有效提取流型特征的兩方面難點,本文通過建立流電耦合仿真計算方法,獲取更為接近真實多相流動過程的動態模擬測量數據,并對4種典型流型下的原始壓力和電容模擬信號開展了多分辨率時頻分析,從而有效解析流型變化與復雜測量信號特征分量的對應關系。連續小波變換實現了一維信號在時-頻二維平面上的顯示,有助于揭示流量波動和特征頻率在時間上的變化規律。基于離散小波變換的多分辨率解析則將原始信號在相同的時間尺度上分解為多個不同頻率尺度下的分量,獲得各分量的能量占比,進而提取不同流型下的主導頻率和能量分布變化。在流型的演化過程中,壓力和電容測量信號的主導頻率呈現一致的遷移規律,為流型轉變的提供了有力判據。本文提出的信號分析方法在段塞流的實驗測量中得到了檢驗,有效識別了段塞特征的細節變化。本文將ECT 測量的可視化圖像與時頻信號分析相結合,為精準多相流過程監測和計量提供了可靠依據和技術手段。

1 油氣兩相流動態測量過程仿真計算

1.1 流電耦合模型及計算平臺構建

數字孿生通過構建多尺度、高保真的虛擬孿生體以及實際物理系統與虛擬空間之間的信息交互,從而準確獲取單個部件乃至整個系統的動態信息,實現狀態監測、性能診斷、行為預測和高效的數據分析[12]。本文基于前期研究中構建的流電耦合數值模擬方法[13],通過耦合多相流動的計算流體力學(CFD)模型和ECT 的軟測量模型,將多物理場同步計算下的流體力學參數和電學參數進行融合,從而實現傳感器測量過程從物理空間向數字空間的精準映射。

多相流可以通過連續性方程和N-S方程來描述其質量和動量守恒,如式(1)、式(2)。

質量控制方程(連續性方程)

動量控制方程(N-S方程)

式中,Sm為源項;ρ為流體密度;u為速度矢量;p表示壓力;μ為流體黏度;f為體積力矢量。

本文采用VOF模型[14]對多相流體進行求解,對于水平管氣液兩相流的典型流態可以獲得較高精度的計算結果。

不同相被視作是共存且互相摻混的連續流體,通過引入相體積分數αq的概念來描述上述的守恒方程。采用固定網格的體積追蹤方法以模擬氣液界面,并重構相應的流場參數,例如n個組分構成的多相流體的密度可表示為式(3)。

ECT 測量對應靜電場的控制方程為簡化的Maxwell方程,如式(4)。

式中,ε為介電常數;φ為電勢;Q?為電荷量。

在現有的ECT 成像算法中,線性反投影(linear back-projection, LBP)算法原理最為簡單,且計算速度快[15]。考慮到成像的實時性和對不同流型的適應性,經過初步算法對比測試確保成像質量在可接受范圍內,本文統一采用LBP算法進行圖像重構,根據測量得到的歸一化電容向量λ和靈敏度矩陣S,反演出灰度值分布矩陣g,如式(5)。

依據上述氣液多相流和ECT 的物理方程,兩者耦合的關鍵在于重構流場中的介電常數分布。與式(3)類似地,在VOF 模型下建立多相流體的介電常數ε與相體積分數αq之間的數學關系,如式(6)。

由此將多相流場計算結果轉換為Maxwell方程求解的條件,從而實現多相流場與靜電場的同步計算。

基于上述流電耦合模型,本文借助Fluent 19.2、COMSOL Multiphysics 5.4 和MATLAB R2018b 多 平臺協同,實現了動態測量仿真計算。Fluent軟件中基于VOF 模型求解多相流控制方程,并借助用戶自定義函數(user defined function, UDF)重構感測區域內的介電常數分布。COMSOL軟件中根據介電常數分布求解靜電場方程以模擬ECT 的測量過程,從而獲得每一次測量的獨立電容值。MATLAB軟件中借助LBP算法根據獨立電容值進行圖像重建,同時實現多相流場與靜電場之間的數據交互。通過對比流場計算結果與重構圖像可以驗證計算模型和架構的有效性,并對仿真測量進行質量評估。

1.2 幾何模型、網格及工況設置

圖1中顯示了本文數值模擬中構建的幾何模型以及對應的網格結構,包括12m 長、內徑為50mm的水平圓管和位于距管道入口10m處的八電極ECT傳感器,ECT 傳感器的單個電極徑向角為37.5°。流場模擬中的三維網格共有817718 個單元,靜電場模擬中的二維截面包括1194 個網格,均經過事先驗證以確保網格的獨立性和仿真精度。依據實驗室測試結果設置了表1中的油、氣兩相物性參數設置,其中氣相為空氣,油相為實驗室常用的白油。

表1 仿真計算物性參數設置

圖1 幾何模型及網格結構

經驗流型圖是長期以來多相流動研究中的常用手段之一。鑒于表1中的物性參數的適用范圍,使用Mandhane 流型圖[16]參照選取了表2 中的模擬工況,初步計算后在水平管內形成了相應的4種典型流態,即波浪流、氣團流、段塞流和泡狀流。

表2 仿真計算進口條件及多相流體速度計算結果

1.3 計算結果及初步分析

參考壓力信號,本文以電容平均值μc作為特征參數,將ECT 原始測量數據轉化為一維時間序列信號以便分析。由于采用八電極傳感器和單電極激勵方式,每次測量可獲得C28= 28 個獨立電容值,故μc計算公式如式(7)。

4種模擬流型下的壓力、電容原始信號以及相應的ECT 成像如圖2 所示,由于壓力變化范圍較大,為便于觀察,壓力信號的縱坐標選用了對數坐標。ECT圖像包括基于LBP算法生成的瞬時橫截面圖像以及管道軸向的流態重建圖像,軸向圖像是通過多幀橫截面圖像按時間順序堆疊而成的,能夠更加清晰地將管道內部的多相流場進行可視化表達。所有測試流型下的ECT 成像精度均已經過事先測試,ECT圖像與Fluent流場計算結果得到的流態始終保持一致。基于ECT成像得到的4種流型下截面空隙率的時變規律也與流場計算結果一致,平均相對誤差依次為3.62%、1.93%、5.47%和18.60%。

圖2 不同流型下的仿真測量信號及ECT重構圖像

結合并對比ECT的反演圖像和測量信號,4種流型及其測量信號主要有以下特征規律。

(1)波浪流持液率較低且流動狀態相對穩定,壓力值和電容值明顯低于其他三種流型,但測量信號均能靈敏感知其液位變化。

(2)氣團流的液相流量顯著增加,持液率較高,氣流吹起液波幅度會間歇性達到管頂,阻塞氣體流通面積,形成液塞和氣團交替出現的流動狀態。電容和壓力信號均較波浪流有顯著提升,并呈現出類似于周期性的變化規律。ECT圖像能夠感知氣團的出現,并準確還原氣團的形態。

(3)段塞流與氣團流類似,但通常具有更高氣相速度,導致相對較低的段塞頻率,以及更加劇烈的兩相間相互作用因而形成的不規則油氣界面。壓力和電容信號均能反映出段塞頻率和液塞形態的相應變化,且壓力信號在更高流速下表現出了更加寬范圍的周期振蕩。ECT的重構圖像依然可以較為敏感地還原液塞和氣團的形態特征。

(4)泡狀流的液相流量進一步提升,液體的湍流脈動破碎管道頂部的氣團,使其變為小氣泡分散于液相上層。ECT圖像能夠很好反映管道頂部以氣相為主、下部以大量液相為主的兩相分布規律。由于其高流體速度和高液相含率,壓力和電容信號值明顯高于其他三種流型,但變化較為平穩。

此外,根據測量原理及計算結果可知,壓力信號值與相含率、流速以及湍流度等多因素均相關,而電容信號值僅與相含率和氣液橫截面分布相關。盡管測量信號的時間序列波動和ECT 重構圖像能夠反映基本的流態結構和變化規律,但無法解析出耦合在原始信號中的各項流動參數和流場信息,亦無法定量分析以建立更加明確的流態判據,需要進一步應用先進信號分析方法對其進行深入探究。

2 測量信號的多分辨率時頻分析

2.1 連續小波變換

現實中的測量信號通常是非平穩信號。常見的傅里葉變換被廣泛應用于解析信號中出現的頻率分量[17],但其無法提供不同頻率發生的時間定位。故在此基礎上發展出了短時傅里葉變換(short-time Fourier transform, STFT)[18],將信號拆分成不同單元,從而獲取信號的局部特征。然而,STFT 中的窗口大小在整個計算過程中是固定的[圖3(a)],時間分辨率和頻率分辨率不能兩者兼顧。連續小波變換(continuous wavelet transform, CWT)[7]則克服了STFT 中的固有分辨率問題[圖3(b)],在高頻區域使用高時間分辨率,在低頻段則使用高頻率分辨率,這將能適應現實中很多信號在長時間內具有緩慢振蕩、在高頻則多以瞬態突變居多的物理特性。

圖3 短時傅里葉變換和連續小波變換的時頻拆分特性

信號f在尺度s下的連續小波變換公式由式(8)表示[9]。

式中,ψ為小波基函數;上角標*表示復共軛;通過調節尺度s和位移τ的大小來改變基本小波的形狀以獲得不同時頻尺度下的信息。

連續小波變換公式也可寫為式(9)的內積形式,表明小波變換系數C(s,τ)是原函數與小波基函數相似性的度量。圖4中顯示了取絕對值后不同流型下基于式(8)生成的CWT系數圖。

圖4 不同流型下的CWT系數圖

不同流型下的測量信號CWT 結果呈現出鮮明的特征差異。

(1)波浪流以中低頻段、長時間穩定的小幅震蕩為主,CWT 圖像可以清晰反映液位變化的劇烈程度,并定位其出現的時間。在0.5~1s 的時間段中,系數絕對值明顯增大,且伴有高頻響應,表明此時的液位變化較其他時段更為劇烈。

(2)對于氣團流和段塞流這兩種間歇流型,系數較大的部分仍分布在較低頻段,并在液塞通過的時間區域內有顯著響應,且段塞流的液塞頻率更低但能量更為集中。

(3)泡狀流的波動則主要集中在高頻段,反映了以氣泡行為為主導的流態特征。與波浪流類似,電容信號的小波變換系數整體偏低,表明氣液比隨時間變化很小,流態相對穩定。但壓力信號受高流體速度的影響,系數較波浪流顯著升高。

在模擬的4種流型下,小波變換系數的時頻分布在壓力和電容信號上都呈現出一致的規律性,且對于氣液結構的瞬態變化能夠提供準確的時間定位。能量集中的頻率范圍會隨流型變化而產生顯著差異,這意味著對原始信號進行頻段分解可能有助于解析不同流型特征在測量信號上的相關表征,從而建立更為量化的流型判據。

2.2 小波多分辨率分析

信號通常由多個具有物理意義的分量構成。多分辨率分析方法(multi-resolution analysis, MRA)將信號可逆地分解為不同頻段上的多個分量,從而能夠在與原始信號一致的時間尺度上單獨研究每一個分量的變化,常用方法包括小波多分辨率分析[19]、經驗模態分解(empirical mode decomposition, EMD)[20]、變分模態分解(variational mode decomposition,VMD)[21]等。小波MRA 通過小波的固定函數來分離信號,對于在較低頻率上分離多個頻率接近的振蕩分量以及捕捉數據中的瞬態變化通常有較好的效果。參考連續小波變換分析結果,表征多相流型變化的信號通常集中于較低頻率,壓力信號響應多集中在30Hz 以內,電容信號大多在10Hz 以內。同時,對于間歇流,液塞出現對應的信號突變是研究中關注的重點之一。因此,小波MRA 方法在原理及特性上適用于本文研究中的多相流動信號分析。

信號f對于小波函數Ψ和尺度函數φ的小波級數展開定義為式(10)。

式中,c和d分別為近似系數和細節系數。

由式(10)的計算原理可知,若令初始尺度j0=0,原始信號f經過有限N層小波分解之后會獲得N個由小波函數產生的高頻“細節”組分(detail,D)和1個由尺度函數產生的低頻“近似”組分(approximation,A)。本文中對原始測量信號進行了6 層小波分解,Dk層的頻率范圍應為為采樣時間間隔。壓力信號的采集速度統一設置為100Hz,電容信號則為50Hz。表3中列出了分解后不同頻段的頻率范圍,并統一定義了本章后續分析中的高頻、中頻、低頻所指代的頻段。

表3 壓力和電容信號在多分辨率分析中的頻段劃分

不同流型下壓力和電容信號的多分辨率分析結果如圖5所示,每幅子圖中均包含6個“細節”組分和1個“近似”組分,隨著頻率從上至下逐漸降低,組分信號變得越來越平滑。圖中分別用藍、灰底色標記出了構成信號最主要和次主要的頻率分量,最主要頻率分量隨流型變化的遷移規律在壓力和電容信號的分解中展現出了令人滿意的一致性。當流型從波浪流到間歇流再到泡狀流的變化過程中,主導頻率先由中頻段轉向低頻段再轉至高頻段,這分別對應了波狀流中的中等頻率液面波動、間歇流中的低頻液塞以及泡狀流中的高頻氣泡行為。由于壓力和電容參數的影響因素和所表征的流動信息不同,信號分解結果在次主要頻段的分布上呈現出細節差異。例如波浪流下壓力信號的次主要分量位于比主分量更低的兩個中低頻段,而其電容信號的次主要分量則分別位于較主分量一高一低的兩個分量。較低頻率的分量反映了波浪流全局上的穩定流動,而中頻段則對應了氣流吹拂下的液位波動。在兩相物性不變的條件下,由于電容值僅與相含率和兩相分布相關,對于液位波動會更加敏感。

圖5 不同流型下的小波多分辨率信號分析

為了更加直觀量化地觀察各頻率分量,定義各組分的信號能量為該組分小波分解系數的平方和,并根據式(11)計算第L層細節組分在所有細節項中的能量占比。

式中,k為小波分解系數數量。

圖6中顯示了各頻段的小波能量比,進一步佐證了上述頻率遷移規律,并使得各頻段的能量分布差異更為直觀。而在后續實驗分析中,能量分布圖還將表現出在同一流型下展現不同流動細節的能力。

圖6 各頻段小波能量分布

3 動態實驗測量

在模擬測量信號分析的基礎上,進一步對段塞流的ECT 實驗測量信號展開多分辨率時頻分析,以驗證上述數據分析方法和模擬結論的可靠性,并著重分析在同種流型下小波分析方法對于流場細節變化的提取能力。

3.1 多相流實驗平臺

圖7中顯示了多相流實驗平臺的示意圖,不同的流體工質從各自的儲罐流出,混合后在測試管段形成氣液多相流,應用ECT 傳感器進行測量后分離回到各自儲罐,從而進行循環流動。鑒于本文的研究對象,僅使用油路和氣路進行實驗,且兩相物性參數與表1中仿真計算的設置一致,進口流量通過流量調節閥進行調節以形成不同的兩相流動狀態。實驗管道的橫截面結構也與仿真計算的幾何模型相同,均為內徑50mm的水平圓管。ECT傳感器為八電極結構,電極徑向角為31°。參考仿真計算結果,電容信號中表征流型變化的響應大多集中于10Hz 以內。綜合考慮信號分析的有效性、實驗室傳感器的采集能力和長時間測量的數據量,實驗測量中的數據采集頻率選取20Hz。

圖7 多相流實驗平臺示意圖

3.2 段塞流過程監測及信號分析

實驗在30m3/h的恒定氣相流量下進行,通過在0.5~2.5m3/h范圍內調節油相流量,在被測管段形成不同的段塞流形態。每一測試工況下取穩定運行后60s 內的測量信號進行分析,為便于觀察,表4 中對比了0~20s 內的ECT 圖像與小波變換結果。由ECT軸向重建圖像可知,隨著液相流量的增長,液塞呈現出先變長而后變短促的特征規律。這一特征變化在連續小波變換結果中表現為整體持液率提升導致的低頻段響應逐漸增強,以及液位波動頻率先變緩后加劇導致的較高頻段振蕩先減弱后增強。與2.1節中的仿真結果進行比對,CWT 的分析結果在流型基本特征的反映上具有一致性,在0~10Hz 的頻率范圍內能夠明顯捕捉段塞流的間歇振蕩并提供時間定位。但實驗信號整體上的時頻特性會更加復雜,尤其當觀測的時間尺度顯著增長時,僅從CWT 變換結果中直接辨別出流場細節的特征變化更為困難。

為了進一步清晰展現流型和流場細節變化對測量信號的影響,對電容信號進行小波多分辨率分析。由于均為段塞流狀態,不同液相流量下各信號分量的整體波動規律接近,故選擇圖8中更加簡潔和便于量化觀察的能量分布圖進行分析。由圖可知,低頻段D5始終為主導頻段,對應實際頻率范圍為0.31~1.25Hz,與第2 節中數值模擬的段塞流主頻范圍(0.39~1.56Hz)基本一致。而隨著液相流量增加,高頻段D1的能量比重明顯回落,而低頻段D4、D5的主導優勢變得更加明顯。進一步觀察主頻段D5,液相流量由0.5m3/h 提升至1.5m3/h 的過程中,能量占比明顯增加,這反映了穩定段塞頻率下液塞長度顯著增加。而液相流量繼續增加到2.5m3/h的過程中,D5的能量占比有所下降,而D3、D4的能量有所提升,這反映了液塞長度變短且段塞頻率增長的趨勢。

圖8 不同油相流量下的小波能量分布

通過動態實驗測量和上述相關數據分析,本文應用的信號時頻分析方法在實際測量中可以有效揭示流動狀態和氣液結構特征。作為流型的主要判據,實驗測量信號的主頻范圍與仿真計算結果一致,驗證了數值計算模型和流型判斷依據的有效性。連續小波變換結果能夠有效定位頻率變化出現的時間,這對于段塞流這類瞬時變化劇烈的間歇流的過程狀態監測能夠起到重要作用。而小波多分辨率分析結果一方面可以在低頻段提供較高的頻率分辨率,從而有效通過主導頻段確定流型;另一方面通過分析能量遷移規律,能夠簡潔快速地獲取氣液兩相結構的細節變化。將基于小波的多分辨率時頻分析與常規的ECT 可視化圖像相結合,可以僅在電容測量的單一模態下同時為管內多相流流動狀態提供定性和定量的特征描述,有效拓展了電容數據在多相流測量上的表征作用,從而為實時在線監測和過程控制提供更加豐富的多維信息。

4 結論及展望

基于流電耦合仿真和多相流實驗對油氣兩相流過程進行動態監測,并基于小波時頻分析方法將原始時間序列信號在時間域和頻率域上同步解析,在流動可視化的基礎上進一步挖掘了測量信號中蘊含的流場特征信息。主要研究結論如下。

(1)連續小波變換能夠有效提取油氣兩相流動中長時間的低頻波動和瞬態的高頻突變,并給予時間定位和強度表征。

(2)小波的多分辨率分析從頻率角度提供了不同流型的量化判據,當兩相流動在“波浪流—間歇流—泡狀流”的演化過程中,主頻帶相應呈現出“中頻段—低頻段—高頻段”的遷移規律。

(3)不同頻段上的小波能量分布變化能夠簡明反映氣液兩相結構的動態變化,在確定主導頻率的基礎上,有效識別段塞頻率和液塞形態的變化規律。

在本文研究的基礎上,可進一步開展以下工作。

(1)在仿真計算以及實驗、工業環境的實際測量中拓展數據集測試范圍,以適應不同流體的物性差異和多變的實測環境因素。

(2)借助連續小波變換的時間定位能力,修正現有的流量計算模型,以提高計算模型的適應性和準確性。

(3)將ECT 圖像和測量信號的時頻特征分量共同作為模式識別的有效輸入,提升流型判別模型的準確性和可解釋性,進而優化不同流型下圖像重建效果和關鍵參數的計量精度。

符號說明

A—— 小波多分辨率分解后的近似組分

C—— 小波變換系數

c—— 近似系數

D—— 小波多分辨率分解后的細節組分

d—— 細節系數

f—— 體積力矢量,N

g—— ECT反演圖像的灰度值分布矩陣

p—— 壓力,Pa

Q—— 電荷量,C

S—— ECT靈敏度矩陣

Sm—— 源項

s—— 小波函數尺度

u—— 速度矢量,m/s

αq—— 相體積分數

ε—— 介電常數,F/m

λ—— 歸一化電容向量,F

μ—— 黏度,Pa·s

μc—— 電容平均值,F

ρ—— 密度,kg/m3

τ—— 小波函數位移

φ—— 尺度函數

Ψ—— 小波函數

ψ—— 小波基函數

猜你喜歡
測量信號
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
完形填空二則
把握四個“三” 測量變簡單
滑動摩擦力的測量和計算
孩子停止長個的信號
滑動摩擦力的測量與計算
測量的樂趣
測量
基于LabVIEW的力加載信號采集與PID控制
一種基于極大似然估計的信號盲抽取算法
主站蜘蛛池模板: 色噜噜在线观看| 亚洲综合网在线观看| 久久婷婷五月综合色一区二区| 热思思久久免费视频| 亚洲乱码精品久久久久..| 国产成人艳妇AA视频在线| 亚洲一区二区三区在线视频| 国产特一级毛片| 91精品国产麻豆国产自产在线| 在线高清亚洲精品二区| 2021最新国产精品网站| 97人妻精品专区久久久久| 日本人真淫视频一区二区三区| 亚洲高清日韩heyzo| 露脸真实国语乱在线观看| 怡春院欧美一区二区三区免费| 亚洲人成成无码网WWW| 九九热这里只有国产精品| 无码啪啪精品天堂浪潮av| 国产精品污污在线观看网站| 国产精品太粉嫩高中在线观看| 久久性妇女精品免费| 欧美日韩精品在线播放| 精品视频一区二区三区在线播| 亚洲视频欧美不卡| 国产熟睡乱子伦视频网站| 免费无码AV片在线观看国产| 99视频在线免费| 四虎影视永久在线精品| 中文精品久久久久国产网址| 刘亦菲一区二区在线观看| 国产精品亚洲片在线va| 国产精品美女网站| 亚洲中文字幕在线一区播放| 国产成人一区| 国语少妇高潮| 伦精品一区二区三区视频| 日本欧美中文字幕精品亚洲| 日本在线欧美在线| 亚洲福利一区二区三区| 亚洲精品国产精品乱码不卞| 天天爽免费视频| 操操操综合网| 久久无码免费束人妻| 一级毛片在线播放| 亚洲人成网7777777国产| 亚洲国产91人成在线| 欧美专区在线观看| 久久国产精品国产自线拍| 全部无卡免费的毛片在线看| 亚洲黄网视频| 日韩AV无码免费一二三区| 97视频在线观看免费视频| 91精品国产丝袜| 狠狠色噜噜狠狠狠狠奇米777| 日韩亚洲高清一区二区| 人人91人人澡人人妻人人爽| 国产在线一二三区| 亚洲中文在线看视频一区| 成年网址网站在线观看| 久久国产乱子| 99热这里只有免费国产精品| 凹凸精品免费精品视频| 97超碰精品成人国产| 欧美一道本| 五月激情婷婷综合| 国产精品密蕾丝视频| 免费观看国产小粉嫩喷水| 亚洲永久免费网站| 欧美在线免费| 亚洲欧美日韩中文字幕一区二区三区| 国产在线拍偷自揄观看视频网站| a毛片免费在线观看| 97国产在线播放| 亚洲精品无码av中文字幕| 国产免费羞羞视频| 日韩精品毛片| 91色在线观看| 国产香蕉在线| 久久一色本道亚洲| 亚洲中文字幕无码爆乳| 日韩天堂在线观看|