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

直管科氏質(zhì)量流量傳感器靈敏度分析軟件設(shè)計(jì)

2022-01-12 10:23:54韓明哲鄭德智
測控技術(shù) 2021年12期
關(guān)鍵詞:測量分析質(zhì)量

胡 純,韓明哲,鄭德智,3,彭 鵬

(1.北京航空航天大學(xué) 電子信息工程學(xué)院,北京 100191;2.北京航空航天大學(xué) 儀器科學(xué)與光電工程學(xué)院,北京 100191;3.北京航空航天大學(xué) 前沿科學(xué)技術(shù)創(chuàng)新研究院,北京 100191;4.重慶德新機(jī)器人檢測中心有限公司,重慶 400700)

流量測量、壓力測量與溫度測量被稱為工業(yè)測量領(lǐng)域的三大參數(shù)測量,其中以流量測量最為復(fù)雜、精度要求最高,其測量的經(jīng)濟(jì)性、安全性、可靠性對生產(chǎn)過程有著重要影響,在工業(yè)中具有重要地位[1]。流量傳感器分為體積流量傳感器與質(zhì)量流量傳感器。體積流量會(huì)受到壓力、溫度、黏度與密度等因素的影響,并且標(biāo)定過程煩瑣;質(zhì)量流量則不受上述因素的影響,并且適用于氣液兩相流,因此日益得到人們的青睞[2]。

目前針對靈敏度的理論較為成熟,對測量管振動(dòng)模型、集中質(zhì)量的影響、最佳檢測位置等的研究較為充分,針對U型管和直管都進(jìn)行了較多研究。文軍浩等[3]關(guān)注了科氏質(zhì)量流量傳感器中的非線性因素,并研究了其對U型管質(zhì)量流量傳感器的影響。樊剛等[4]就溫度對科氏質(zhì)量流量傳感器進(jìn)行了分析,并建立了流量傳感器溫度補(bǔ)償模型,對溫度補(bǔ)償進(jìn)行了仿真。曹勝強(qiáng)等[5]將污垢按集中質(zhì)量處理,通過微分方程和ANSYS仿真的方法分析了污垢對單直管科氏質(zhì)量流量傳感器靈敏度的影響。部分科氏質(zhì)量傳感器采用了在測量管兩端安裝波紋管的結(jié)構(gòu),波紋管剛度遠(yuǎn)小于鋼管,因此此類測量管狀態(tài)與一般測量管不同,徐偉國等[6]和紀(jì)彩虹等[7]分別就兩段波紋管和三段波紋管的直管科氏質(zhì)量流量傳感器進(jìn)行了靈敏度計(jì)算,但尚沒有用于科氏質(zhì)量流量傳感器靈敏度分析的軟件。為彌補(bǔ)此空缺,本文設(shè)計(jì)了一款具有靈敏度分析、頻率計(jì)算、最佳檢測位置分析等功能的軟件,為科氏質(zhì)量流量傳感器的設(shè)計(jì)與開發(fā)提供便利。

首先從科氏質(zhì)量流量傳感器的數(shù)學(xué)模型入手,對科氏質(zhì)量流量傳感器進(jìn)行了分析,建立了測量管的歐拉梁靜力模型,并分別計(jì)算了測量管在激振力與科氏力作用下的撓度曲線;建立了測量管的微分方程模型,并通過求解其振動(dòng)微分方程,計(jì)算出了測量管在激振力下的諧振頻率和科氏質(zhì)量流量傳感器的靈敏度;接下來分析了測量管的壓力損失,給出了在層流和湍流下的壓損計(jì)算公式;研究了溫度對科氏質(zhì)量流量傳感器的影響,對當(dāng)前測量管的常用材料進(jìn)行了調(diào)研,采用線性模型估計(jì)了溫度對彈性模量的影響,并用以對分析結(jié)果進(jìn)行補(bǔ)償。

在理論分析的基礎(chǔ)上,通過Matlab軟件的APP Designer模塊制作了直管科氏質(zhì)量流量傳感器的分析軟件,實(shí)現(xiàn)了計(jì)算靈敏度、頻率、壓損等數(shù)據(jù)的功能;并使用ANSYS Workbench進(jìn)行有限元分析;對其進(jìn)行了靜態(tài)結(jié)構(gòu)分析,分別得到其在激振力和科氏力作用下的靜撓度曲線,由此計(jì)算出傳感器的最佳安裝位置;進(jìn)行了瞬態(tài)動(dòng)力學(xué)分析,計(jì)算測量管在激振力作用下的振動(dòng);進(jìn)行了模態(tài)分析得到其各階模態(tài)及諧振頻率,并分析了溫度對諧振頻率的影響;最后通過諧響應(yīng)分析得到測量管位移振幅曲線。有限元仿真的結(jié)果被用于驗(yàn)證模型與所設(shè)計(jì)軟件的計(jì)算結(jié)果,以確保所設(shè)計(jì)軟件的正確。

1 直管型科氏質(zhì)量流量傳感器的數(shù)學(xué)模型

為設(shè)計(jì)靈敏度分析軟件,首先要建立科氏質(zhì)量流量傳感器的數(shù)學(xué)模型。為此分別采用了靜力學(xué)和微分方程的方法進(jìn)行分析,得到了靈敏度等傳感器指標(biāo)的計(jì)算公式,為GUI軟件的設(shè)計(jì)提供支持。

1.1 直管型科氏質(zhì)量流量傳感器的靜力學(xué)模型

在科氏質(zhì)量流量傳感器工作時(shí),其激振力頻率應(yīng)當(dāng)接近測量管的固有頻率,此時(shí)為簡便分析過程,可以將測量管在恒定激振力下的靜撓度曲線近似為測量管撓度曲線進(jìn)行分析[8]。直管型測量管的雙端固定且長細(xì)比大,符合雙端固支梁模型,可視為雙端固支梁進(jìn)行計(jì)算。

雙端固支梁為超靜定結(jié)構(gòu),為對其進(jìn)行分析,需要解除多余的約束條件并采用等效的作用力和力矩建模。將激振力F=Ffsin(ωt+φ)視為恒定的力F,設(shè)測量管兩端約束力與約束力矩分別為Fa,Fb,Ma,Mb,由于測量管結(jié)構(gòu)對稱,因此有

(1)

Ma=Mb

(2)

為計(jì)算兩固支端的約束力矩,考慮其兩端和中心的轉(zhuǎn)角。由于兩端是固支端,其3個(gè)自由度都受到約束,其轉(zhuǎn)角和撓度均為0;又由于雙端固支梁為對稱結(jié)構(gòu),故在作用于中點(diǎn)的激振力作用下,梁的形變同樣具有對稱性。因此梁的中點(diǎn)轉(zhuǎn)角也為0。

通過材料力學(xué)的方法分析梁的彎曲變形,通常借助梁的撓曲線近似微分方程:

(3)

對于集中力P作用下的測量管,其彎矩為

(4)

使用積分法求解梁的中點(diǎn)轉(zhuǎn)角,將式(3)從梁的端點(diǎn)到中點(diǎn)積分,得

(5)

(6)

因此,結(jié)合式(3)~式(6),得到在激振力F作用下梁的撓度微分方程:

(7)

對其進(jìn)行積分得到梁的轉(zhuǎn)角方程和撓度方程:

(8)

(9)

激振力F作用下梁的靜撓度曲線如圖1所示。

圖1 激振力F作用下梁的靜撓度曲線

圖1中對橫坐標(biāo)x和撓度y都做了去量綱處理。

在工作過程中,測量管所受科氏力載荷為

dFk=2Qω(x)dx

(10)

式中,w(x)為角速度,即角度對時(shí)間的微分:

(11)

為對測量管上的科氏力載荷進(jìn)行分析,考慮在測量管上橫坐標(biāo)a處施加集中力dFk的情況。可以使用與分析激振力時(shí)相同的方法,將其視為恒定的力,解除其一邊約束,并加以等效的約束力F′(a)和約束力矩M′(a)進(jìn)行求解。再將求解出的撓度曲線關(guān)于dFk進(jìn)行積分,得到科里奧利力載荷作用下的撓度曲線。經(jīng)計(jì)算,前半段梁的撓度為

5.64×10-4x3L3+1.95×10-4x2L4)

(12)

后半段梁的撓度可由前半段關(guān)于中點(diǎn)對稱得到。其撓度曲線如圖2所示。

圖2 科里奧利力載荷作用下梁的靜撓度曲線

分析式(12),可以得到在x=0.287L處,梁的撓度取得最大值。因此,最佳檢測位置在0.287L處和0.713L處。

1.2 直管型科氏質(zhì)量流量傳感器的微分方程模型

對于等截面歐拉梁振動(dòng)力學(xué)中已有成熟的理論,根據(jù)達(dá)朗貝爾原理可以列出其動(dòng)力學(xué)方程:

(13)

式中,X,Y分別為測量管橫縱坐標(biāo);T為時(shí)間;V為待測流體流速;E為測量管的彈性模量;I為測量管的截面慣性矩;m1,m2,m3分別為科氏質(zhì)量流量傳感器的檢測器和激振器質(zhì)量,其中m1,m3為檢測器,m2為激振器,兩檢測器質(zhì)量相等,安裝在關(guān)于中點(diǎn)對稱的位置上;mp,mf分別為測量管線密度和待測流體的線密度。

為計(jì)算諧振頻率,需要對該微分方程進(jìn)行求解。首先為方便分析,將其變形為無量綱形式:

(14)

式中,x,y,t,σk,v,β分別為無量綱的橫坐標(biāo)、縱坐標(biāo)、時(shí)間、集中質(zhì)量、流速、流體密度。

記測量管震動(dòng)角頻率為Ω,測量管無量綱振動(dòng)角頻率為ω:

對于此類微分方程,已經(jīng)有成熟的求解理論,其解為[9]

y(x,t)=RΦ(x)eiωt

(15)

式中,R為常系數(shù);Φ(x)為振型函數(shù);i為單位虛數(shù)。對測量管的振型函數(shù),可以各階振型函數(shù)構(gòu)成的無窮級數(shù)來表述:

(16)

其振型函數(shù)為[10]

φr(x)=cosh(krx)-cos(krx)-λr(sinh(krx)-sin(krx))

(17)

本文關(guān)注的是梁在激振力作用下的振動(dòng)和在科里奧利力作用下的振動(dòng),分別對應(yīng)一階和二階振型,振型曲線如圖3與圖4所示,其參數(shù)k和λ分別為:k1=4.730,k2=7.853,λ1=0.9825,λ2=1.001。

圖3 一階振型函數(shù)曲線

圖4 二階振型函數(shù)曲線

注意到微分方程解的一、二階振型函數(shù),與前文所推導(dǎo)的測量管在激振力與科氏力作用下的撓度方程具有相同的物理意義,由于在分析過程中采用了不同的方法去量綱以簡化分析,因此二者應(yīng)當(dāng)具有線性關(guān)系。在此采用了相關(guān)系數(shù)的方法進(jìn)行驗(yàn)證。

首先對計(jì)算出的函數(shù)進(jìn)行采樣,將其離散化。采樣間隔設(shè)為0.001,取1000個(gè)采樣點(diǎn),然后分別計(jì)算一階振型與激振力撓度、二階振型與科氏力撓度的皮爾遜相關(guān)系數(shù),分別為r1=0.999052,r2=0.999849,這說明二者高度相關(guān),驗(yàn)證了上述分析結(jié)果。二階振型函數(shù)在x=0.290處取得最大值,與2.1節(jié)中分析得到的最大值位置0.287相接近。歸一化后的對比圖如圖5、圖6所示。

圖5 一階振型與激振力撓度對比圖

圖6 二階振型與科氏力撓度對比圖

由圖5和圖6可以看出,兩種計(jì)算方法所得到的振型與撓度曲線高度相似,此結(jié)果對上述分析的正確性提供了支持。

使用伽遼金法求解微分方程(14),并取其前2階振型作為該微分方程的近似解,解得:

(18)

(19)

在科氏質(zhì)量流量傳感器工作時(shí),測量管同時(shí)存在激振力引起的振動(dòng)和科氏力引起的振動(dòng),這使得測量點(diǎn)存在振動(dòng)的相位差。因?yàn)榭剖腺|(zhì)量流量傳感器通過測量此相位差以得到待測流體的質(zhì)量流量,所以將科氏質(zhì)量流量傳感器的靈敏度定義為單位流量的相位差,即:

(20)

由式(19),可以得到測量管的振動(dòng)函數(shù):

(21)

其相位差Δφ為

(22)

進(jìn)一步計(jì)算出科氏質(zhì)量流量傳感器的靈敏度K:

(23)

其中,

1.3 壓損計(jì)算

壓損按照產(chǎn)生的機(jī)理可以分為沿程壓力損失和局部壓力損失。對于直管型測量管,不存在閥口、彎管等通流截面變化或流動(dòng)方向改變的結(jié)構(gòu),因此局部壓力損失可以忽略不計(jì)。

圓管層流是流體運(yùn)動(dòng)中較為簡單的一種情況,也是能夠得出流速分布及壓力損失解析解的為數(shù)不多的情況之一。經(jīng)計(jì)算,在層流中沿程壓力損失系數(shù)為[11]

(24)

湍流運(yùn)動(dòng)的情況較為復(fù)雜,沿程壓力損失系數(shù)的計(jì)算往往依賴經(jīng)驗(yàn)公式或半經(jīng)驗(yàn)公式。1913年德國水力學(xué)家布拉修斯(Blasius)提出了計(jì)算湍流沿程壓力損失系數(shù)的經(jīng)驗(yàn)公式:

(25)

布拉修斯公式適用于Re<105的紊流光滑區(qū),計(jì)算簡單方便,且具有較高精度。

沿程壓力損失的理論計(jì)算公式為

(26)

式中,l為測量管長度;ρ為測量管密度;v為流體流速,且滿足

(27)

在層流和湍流下測量管的壓損計(jì)算公式為

(28)

式中,Q為流體質(zhì)量流量(kg/s);η為流體動(dòng)力黏度(Pa·s);ρ為流體密度(kg/m3)。

1.4 溫度對科氏質(zhì)量流量傳感器的影響

科氏質(zhì)量流量傳感器在工業(yè)領(lǐng)域中的應(yīng)用逐步廣泛,其所面臨的工作環(huán)境也多種多樣。因此,在各種環(huán)境下都能保持較高的精度已成為科氏質(zhì)量流量傳感器的迫切需要。為此,科氏質(zhì)量流量傳感器通常會(huì)設(shè)置溫度傳感器,測量被測介質(zhì)的溫度,通過溫度補(bǔ)償?shù)姆椒ǖ窒麥囟葞淼挠绊懀员WC測量數(shù)據(jù)的精度。

科氏質(zhì)量流量傳感器最常用的材料是316L不銹鋼,除此之外哈氏合金c-22也被用于科氏質(zhì)量流量傳感器,在直管科氏質(zhì)量流量傳感器上,有時(shí)以硅合金作為測量管材料。筆者選取了最常見的這3種金屬作為備選材料,對其屬性進(jìn)行了調(diào)研,結(jié)果如表1~表4所示。

表1 316L不銹鋼(Acerinox公司)屬性(自20 ℃起)

表2 哈氏合金c-22(Hastelloy公司)屬性(自20 ℃起)

表3 鈦合金Ti 6Al-4V(NeoNickel公司)屬性(自20 ℃起)

表4 常見材料密度

上述各材料的平均線膨脹系數(shù)都在10-5量級,因此在本研究中,可以忽略溫度對密度和測量管結(jié)構(gòu)參數(shù)的影響,溫度的影響主要體現(xiàn)在彈性模量上。對溫度與彈性模量進(jìn)行線性回歸分析,其結(jié)果如圖7~圖9所示。

圖8 哈氏合金c-22彈性模量與溫度的線性回歸模型

圖9 鈦合金Ti 6Al-4V彈性模量與溫度的線性回歸模型

上述回歸分析相關(guān)系數(shù)R2均大于0.997,下面使用F檢驗(yàn)法進(jìn)行回歸方程顯著性檢驗(yàn)。

相關(guān)系數(shù)R2與統(tǒng)計(jì)量F有如下關(guān)系:

(29)

由式(29)計(jì)算得,上述回歸分析統(tǒng)計(jì)量F分別為

F1=5795.89,F(xiàn)2=4340.53,F(xiàn)3=2896.18

對于316L不銹鋼、哈氏合金c-22和鈦合金Ti 6Al-4V的數(shù)據(jù),其樣本容量分別為n=6,5,4。對F檢驗(yàn)進(jìn)行查表,得

F0.01(1,4)=21.20

F0.01(1,3)=34.12

F0.01(1,2)=98.49

因此上述回歸均在0.01的水平上高度顯著,故接受回歸分析的結(jié)果,并在靈敏度計(jì)算中采用上述模型進(jìn)行彈性模量的計(jì)算。在軟件設(shè)計(jì)中還加入了自定義材料的功能,支持直接輸入材料的彈性模量和密度,以適應(yīng)使用其他材料或用戶具有準(zhǔn)確的材料屬性的情況。

2 GUI軟件設(shè)計(jì)

本GUI軟件基于Matlab的APP designer進(jìn)行設(shè)計(jì),由參數(shù)設(shè)置、材料屬性設(shè)置、集中質(zhì)量設(shè)置、最佳檢測位置分析、指定位置靈敏度分析、壓損分析和曲線繪制7個(gè)面板構(gòu)成,根據(jù)以上各節(jié)所推導(dǎo)計(jì)算的公式,完成了研究目標(biāo)中所提出的靈敏度分析、諧振頻率分析、壓損分析等目標(biāo),實(shí)現(xiàn)了輸入測量管結(jié)構(gòu)參數(shù)、材料屬性參數(shù)、傳感器集中質(zhì)量等數(shù)據(jù),進(jìn)行分析計(jì)算最佳檢測位置、靈敏度、頻率、壓損等數(shù)據(jù)的功能,并繪制相關(guān)曲線圖。軟件還具有根據(jù)溫度對材料屬性進(jìn)行補(bǔ)償?shù)墓δ堋TO(shè)計(jì)的GUI軟件功能架構(gòu)圖如圖10所示,GUI軟件界面如圖11所示。

圖10 軟件功能架構(gòu)圖

圖11 GUI軟件界面

Matlab的APP Designer提供了豐富而強(qiáng)大的回調(diào)函數(shù),極大地便利了軟件的設(shè)計(jì)。本軟件充分利用了APP Designer的回調(diào)函數(shù),基本上全部功能都通過回調(diào)函數(shù)實(shí)現(xiàn)。在程序運(yùn)行的初期,先進(jìn)行GUI界面的初始化,生成GUI界面、初始化數(shù)據(jù)并繪制表格和坐標(biāo)軸;然后等待用戶對界面進(jìn)行操作,根據(jù)用戶的操作執(zhí)行相對應(yīng)的回調(diào)函數(shù)。程序流程如圖12所示。

3 科氏質(zhì)量流量傳感器的仿真分析

使用ANSYS Workbench對測量管進(jìn)行了有限元仿真。首先建立了直管型測量管的模型,然后分別對其進(jìn)行了靜態(tài)結(jié)構(gòu)分析、瞬態(tài)動(dòng)力學(xué)分析、模態(tài)分析和諧響應(yīng)分析,以驗(yàn)證理論計(jì)算結(jié)果。所使用的測量管結(jié)構(gòu)為單圓直管,并忽略了集中質(zhì)量,其結(jié)構(gòu)參數(shù)如表5所示。

表5 測量管結(jié)構(gòu)尺寸

3.1 靜態(tài)結(jié)構(gòu)分析

于測量管兩端施加固定約束,施加作用于測量管外表面上、作用位置為測量管中央的正上方、方向豎直向下的激振力F,以進(jìn)行激振力作用下的靜撓度仿真;為計(jì)算科氏力作用下的撓度,采用了分段施加載荷的方式,將測量管分為20段,在每段下施加相應(yīng)的科氏力載荷,以模擬實(shí)際的科氏力載荷分布。在網(wǎng)格劃分方面,采用了由ANSYS軟件智能劃分的自由網(wǎng)格。其仿真結(jié)果如圖13、圖14所示。

圖13 激振力作用下的靜態(tài)結(jié)構(gòu)仿真

圖14 科氏力作用下的靜態(tài)結(jié)構(gòu)仿真

通過將激振力作用下測量管的靜撓度曲線和以靜力學(xué)模型計(jì)算出撓度曲線、微分方程的振型函數(shù)進(jìn)行對比,驗(yàn)證了前文所進(jìn)行分析的正確性。

3.2 模態(tài)分析和諧響應(yīng)分析

同樣于兩端施加固定約束,無預(yù)應(yīng)力輸入,取20 ℃常溫下的材料屬性,同樣采用自由網(wǎng)格劃分。激振力作用下的振動(dòng)為一、二階模態(tài),科氏力作用下的振動(dòng)為三、四階模態(tài),因此選取前六階模態(tài)進(jìn)行模態(tài)分析。得到其前六階模態(tài)諧振頻率如表6所示。

表6 模態(tài)分析諧振頻率

直管型測量管結(jié)構(gòu)關(guān)于其中軸線對稱,因此其一二階、三四階、五六階的諧振頻率相同,其振型函數(shù)在空間上正交,且在繞測量管中軸線旋轉(zhuǎn)90°后重合,可以看作是一個(gè)二自由度振動(dòng)在空間上的分解。

上述仿真得到測量管的諧振頻率為131.53 Hz,運(yùn)行所設(shè)計(jì)的GUI軟件,輸入表1所示的材料屬性數(shù)據(jù)和表5所示的結(jié)構(gòu)尺寸參數(shù),選擇測量管內(nèi)流體為空,將溫度設(shè)置為20 ℃,計(jì)算得到的諧振頻率為132.2 Hz,其相對誤差為+0.5094%。

施加在測量管外表面中央正上方,方向向下的簡諧激振力作為激勵(lì),進(jìn)行諧響應(yīng)分析結(jié)果見圖15。

圖15 幅值響應(yīng)

通過修改模態(tài)分析中的溫度條件,對軟件的溫度補(bǔ)償功能進(jìn)行驗(yàn)證。計(jì)算與仿真結(jié)果如表7所示。

表7 諧振頻率仿真驗(yàn)證數(shù)據(jù)

上述兩組數(shù)據(jù)關(guān)于溫度的線性回歸分析都在0.01的水平上顯著,此結(jié)果支持了溫度對諧振頻率具有線性影響的結(jié)論,驗(yàn)證了前文中進(jìn)行的溫度對科氏質(zhì)量流量傳感器的影響的正確性。

4 結(jié)束語

設(shè)計(jì)了直管型科氏質(zhì)量流量傳感器的靈敏度分析軟件,通過理論分析得到的公式進(jìn)行靈敏度、頻率、壓損的計(jì)算,并具有對溫度進(jìn)行補(bǔ)償?shù)墓δ堋2⑹褂肁NSYS進(jìn)行有限元仿真,對所設(shè)計(jì)軟件的計(jì)算結(jié)果進(jìn)行了驗(yàn)證,表明其具有較高的可靠性。

猜你喜歡
測量分析質(zhì)量
“質(zhì)量”知識鞏固
隱蔽失效適航要求符合性驗(yàn)證分析
質(zhì)量守恒定律考什么
把握四個(gè)“三” 測量變簡單
做夢導(dǎo)致睡眠質(zhì)量差嗎
滑動(dòng)摩擦力的測量和計(jì)算
電力系統(tǒng)不平衡分析
電子制作(2018年18期)2018-11-14 01:48:24
滑動(dòng)摩擦力的測量與計(jì)算
電力系統(tǒng)及其自動(dòng)化發(fā)展趨勢分析
測量
主站蜘蛛池模板: 又大又硬又爽免费视频| 亚洲免费成人网| 国产免费久久精品99re丫丫一| 波多野结衣一二三| 无码一区18禁| 国产精品网址在线观看你懂的| 538精品在线观看| 在线欧美日韩| 97超碰精品成人国产| 国产区免费| 欧美69视频在线| 亚洲日本一本dvd高清| 亚洲成人高清无码| 久久综合丝袜长腿丝袜| 日韩在线中文| 中文字幕永久视频| 狠狠v日韩v欧美v| 久久频这里精品99香蕉久网址| 日韩天堂网| 欧美a级在线| 亚洲国产清纯| 毛片在线区| 亚洲日韩精品伊甸| 91精品小视频| 午夜精品久久久久久久2023| 亚洲中文字幕无码爆乳| 国产一区二区精品高清在线观看| 亚洲欧美不卡| 啪啪永久免费av| 免费国产一级 片内射老| av在线无码浏览| 国产精品欧美亚洲韩国日本不卡| av在线无码浏览| 久久香蕉国产线看观看精品蕉| 2021天堂在线亚洲精品专区| 国产永久免费视频m3u8| 色视频国产| 91九色国产在线| 日本黄网在线观看| 欧美va亚洲va香蕉在线| 美女扒开下面流白浆在线试听| 手机在线免费毛片| 亚洲成人黄色网址| 3D动漫精品啪啪一区二区下载| 欧美在线精品一区二区三区| 亚洲欧美不卡视频| 国内丰满少妇猛烈精品播| 国产美女91呻吟求| yjizz视频最新网站在线| 国产靠逼视频| 精品99在线观看| 亚洲AⅤ综合在线欧美一区| 国产成人精品一区二区三区| 精品中文字幕一区在线| 国产在线八区| 色婷婷在线播放| 中文字幕日韩久久综合影院| 日韩精品亚洲精品第一页| 久久婷婷六月| 狠狠色狠狠综合久久| 2020最新国产精品视频| 欧美伦理一区| 精品一区二区三区中文字幕| 亚洲天堂视频在线播放| 亚洲性一区| 高清久久精品亚洲日韩Av| 四虎免费视频网站| 国产美女无遮挡免费视频| 日韩在线网址| 国产成年女人特黄特色毛片免| www.国产福利| a级毛片免费网站| 欧美午夜视频在线| 国产91无毒不卡在线观看| 国产高清无码麻豆精品| 一级爱做片免费观看久久 | 人人91人人澡人人妻人人爽| 日韩福利在线视频| 亚洲va视频| 亚洲福利视频网址| 精品中文字幕一区在线| 精品国产91爱|