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

受限亞音速氣流中倒置懸臂壁板靜氣彈穩(wěn)定性的理論及實驗研究1)

2020-03-26 02:51:04張德春楊翊仁
力學(xué)學(xué)報 2020年2期
關(guān)鍵詞:理論實驗

張德春 李 鵬 梁 森 楊翊仁

(西南交通大學(xué)力學(xué)與工程學(xué)院,成都 610031)

引言

壁板是常見的結(jié)構(gòu)形式,其廣泛地應(yīng)用在諸多工程領(lǐng)域.壁板結(jié)構(gòu)在氣流中的失穩(wěn)問題是得到廣泛關(guān)注的一個研究重點.航空航天中的超音速壁板氣動彈性問題已有了豐富的研究[1-4].最近針對亞音速壁板氣動彈性問題的研究也隨著高速列車的發(fā)展越來越受到重視[5-7],研究成果也越來越豐富.針對壁板失穩(wěn)問題,相關(guān)的理論及實驗研究主要是基于線性模型,以預(yù)測失穩(wěn)臨界參數(shù)為目標(biāo),考察不同結(jié)構(gòu)形式(如初始構(gòu)型[6-7],邊界條件[8-9],質(zhì)量比及長寬比[10-16]等)及不同流體力建模理論(如準(zhǔn)定常點源[5-6,8]及點渦理論[7,12,17],非定常頻域Theodorsen 理論[13]及時域Peters 有限入流模型[7],數(shù)值渦格法[14-15]及NS-結(jié)構(gòu)方程耦合求解法[16]等)對分析結(jié)果的影響.雖然線性建模方式簡單且缺少對非線性因素的描述,但對深入了解這類動力系統(tǒng)失穩(wěn)問題的本質(zhì)有著不可替代的優(yōu)勢[1].已有的研究結(jié)果表明[5,7-9,13-16],壁板的失穩(wěn)形式(發(fā)散或顫振失穩(wěn))與流體速度及壁板邊界約束形式密切相關(guān).兩端固定的壁板(簡支及固支)在超音速氣流中會呈現(xiàn)顫振失穩(wěn)而在亞音速氣流中卻僅會出現(xiàn)發(fā)散失穩(wěn).然而當(dāng)懸臂壁板置于亞/超音速流體中時卻會分別經(jīng)歷顫振/發(fā)散失穩(wěn).

亞音速氣流中的懸臂壁板相較于兩端固定的壁板而言,一般具有較低的失穩(wěn)流速,這使得其被廣泛應(yīng)用于俘能器的設(shè)計中[18-23].如何優(yōu)化壁板的結(jié)構(gòu)形式[19-20]或設(shè)計流動特征[21]使得其具有較低的失穩(wěn)速度、較強的運動幅度及較多的平衡狀態(tài)是這類能量采集器的研究重點[18-21].已有研究表明,一方面在壁板周圍流場中設(shè)計某些流動約束,如前置圓柱體[21]、設(shè)置壁面限制[22]等,都會有效地降低壁板系統(tǒng)的失穩(wěn)速度.事實上,將懸臂壁板倒置于氣流中可獲得一種簡單卻具有優(yōu)良動力學(xué)特征的結(jié)構(gòu)形式,這是由于倒置懸臂壁板在結(jié)構(gòu)[24]及氣動力非線性[25]的作用下存在多穩(wěn)態(tài)之間的跳躍遷移現(xiàn)象.從目前的文獻調(diào)研結(jié)果來看,針對倒置懸臂壁板氣動彈性問題的研究還比較有限.已有的研究主要是通過數(shù)值模擬及風(fēng)洞實驗手段對結(jié)構(gòu)失穩(wěn)后的多穩(wěn)態(tài)非線性特征進行定量分析[24-26].然而若要對多穩(wěn)態(tài)特征進行更靈活的設(shè)計則需要對其出現(xiàn)的誘因,即靜態(tài)失穩(wěn)問題,從理論及實驗方面進行有針對性的研究.相關(guān)的研究主要是將未受限氣流中的壁板視為升力線(面),采用已有的氣動力理論[26]進行建模分析,然而已有氣動力理論均只適用于非受限氣流,鮮有理論上考察壁面約束限制作用的研究報道.

事實上,受限流體中壁板失穩(wěn)問題也廣泛地存在于實際工程中,如核反應(yīng)堆中的層疊板、隧道內(nèi)的列車蒙皮及地效飛行器等.目前相關(guān)的研究主要是采用理論方法(積分變換法及鏡像函數(shù)法[27]等)、數(shù)值計算(面元法[9]等)對受限流體中懸臂壁板的顫振穩(wěn)定性進行分析.迄今為止,還未見針對亞音速氣流中倒置懸臂壁板這一結(jié)構(gòu)形式開展的相關(guān)理論及實驗研究,本文旨在對這一問題進行研究以期更進一步豐富壁板氣動彈性失穩(wěn)問題的研究.

已有的研究表明,當(dāng)壁板長寬比小于1 時,可將其近似視為二維問題[10,22].本文針對這一類受限亞音速氣流中的二維倒置懸臂壁板的靜態(tài)失穩(wěn)問題進行理論及風(fēng)洞實驗研究,以期了解剛性壁面效應(yīng)對這類壁板結(jié)構(gòu)靜態(tài)失穩(wěn)特性的影響規(guī)律.首先,采用鏡像函數(shù)法描述壁面效應(yīng),基于算子理論對壁板氣動力進行研究;其次,將壁板失穩(wěn)方程轉(zhuǎn)化為定區(qū)間內(nèi)的函數(shù)逼近問題并進行求解;最后,依據(jù)壓桿穩(wěn)定原理設(shè)計壁板靜態(tài)失穩(wěn)的測試方法,完成風(fēng)洞實驗,對理論及實驗結(jié)果進行對比分析.

1 力學(xué)模型及數(shù)學(xué)方程

本文考慮如圖1(a)所示的軸向受限亞音速氣流中倒置的二維懸臂壁板模型,壁板長度為L,左端自由而右端受到固支約束.氣流沿板的軸向x方向流動并在壁板的一側(cè)受到剛性壁面的限制,剛性壁面與壁板平行且兩者相距H.本文旨在分析該壁板模型在氣流中的失穩(wěn)特性,因此僅考慮壁板線性梁式彎曲運動微分方程

圖1 Fig.1

及如下的邊界條件

其中,w,ρs,h分別為壁板的橫向變形,密度及厚度,D=Eh3/12(1-ν2)為壁板的抗彎剛度,ν 為壁板的泊松比,Δp=pH+-pH-為壁板上下表面擾動壓力差(向下為正).

考慮理想歐拉流體,通過弓入等熵條件、小擾動假設(shè)和小擾動速度勢可得壁板上下表面壓力的近似線化表達式[28]

式中,ρ∞為氣流密度,φ 為擾動速度勢函數(shù),其滿足如下的線化速度勢方程

式(5a)為剛性壁面約束條件;式(5b)為氣固相容條件,其中

定義為壁板表面的法洗速度;式(5c)和式(5d)分別為流動的連續(xù)及Kutta-Joukowski 條件.方程(4)的定解還需要給定遠場擾動量的條件和流場的初始條件

2 Possio 積分方程及其求解

注意到方程(4)為關(guān)于時間及空間的偏微分方程,直接求解有很大困難.將其左右兩邊關(guān)于時間作用L(Laplace)變換,并利用初始條件(6b),可得

整理式(8)并將其寫為

其中k=λ/U∞定義為系統(tǒng)的減縮頻率.

對本文研究的亞音速流體而言,A具有正實部.考慮邊界條件(6a),方程(9)的解可以寫作

由式(5b)及式(12)可得壁板表面氣流的橫向速度

由式(5b)及式(13)可得壁板上下表面勢函數(shù)之差

進一步考慮式(3)后可得

雖然式(15)給出了任意運動時壁板表面上的壓力,但其求解仍很困難.對于本文研究的壁板靜態(tài)失穩(wěn)問題,可令k=0 并由式(15)解得

其中,c=2Hβ.注意到函數(shù)f(t)的H(Hilbert)變換的F 變換

的形式后,相應(yīng)定義如下變換[29]

在式(16)兩邊首先作用L 和F 的逆變換,然后考慮式(5b),式(5c),應(yīng)用投影運算PR向[0,L]投影可得如下的積分方程

由式(5d),有PRΔp=Δp,考慮算子QH 滿足交換律(簡要證明見附錄),則式(18)變?yōu)?/p>

上式即為Possio 積分方程[29-31].

考慮式(5d)給出的壓力條件,參考薄翼理論中的S?ghen 解[28],弓入移位Tricomi 算子T[29-30]

將算子T 作用在式(19)的兩邊,可得壓力和法洗速度之間的關(guān)系[29,31]

其中,I 為單位算子,與未受限流體相關(guān);而壁面效應(yīng)則表征為一復(fù)合算子P,定義為

由式(21)可得

3 失穩(wěn)方程及其最小二乘解

考慮壁板的發(fā)散失穩(wěn)方程

上式為關(guān)于w的四階偏微分方程,將其擴階寫作

采用常數(shù)變易法求解(24)可得

在式(25)兩端作用D 的逆算子可得

其中

將壁板的邊界條件(2)代入式(26)后整理得到

式(27)僅有零解對應(yīng)系統(tǒng)的穩(wěn)定狀態(tài);當(dāng)其系數(shù)行列式等于零時系統(tǒng)出現(xiàn)非零解,意味著系統(tǒng)處于失穩(wěn)狀態(tài),即

其中,函數(shù)y(x)滿足y(x)≡1.

因算子I-R2是緊的[31],對所求參數(shù)α,其逆都是存在的,那么

由于R2中含有需要確定的參數(shù)α,直接采用式(29)求解并不利于計算.但注意到

上式表明方程(30)的解等價于尋找目標(biāo)函數(shù)f(x)使

將R2代入式(31)有

因(I-P)為壓縮映射[29],(I-P)-1T [f(x)]可近似展開為

假設(shè)所求函數(shù)f(x)在[0,L]區(qū)間內(nèi)是連續(xù)的,并注意到f(L)=0,由Weierstrass 定理可將f(x)近似展開為

考慮式(32)在[0,L]區(qū)間內(nèi)任意M,(M?N)個離散點上均成立.那么求解方程式(32)等價于求解如下的矛盾方程組

采用最小二乘方法可獲得方程組(35)的最優(yōu)解為

為尋找臨界參數(shù)αc,定義關(guān)于α 的誤差函數(shù)

改變α 使得誤差|E| ≤1.0e-4可得臨界參數(shù)αc及其對應(yīng)的臨界氣流速度.

4 理論結(jié)果及對比驗證

圖2 本文理論分析和計算的流程Fig.2 The flowchart of the present theoretical analysis

在下面計算中選取無量綱參數(shù)ξ=x/L,rH=H/L.為避免間隙H過小而導(dǎo)致較強的粘性作用,本文僅考慮rH>0.05 的情況.基于圖2 給出分析及計算流程進行算例分析.首先考察rH→∞,即無壁面的情況.此時,‖P‖→0,(I-P)-1→I.圖3(a)給出了此工況下,M=200 時,誤差函數(shù)E 隨截斷階數(shù)N的變化關(guān)系.由該圖可知隨著N的增加誤差會最終會在αc=1.84 處滿足設(shè)置的條件(進入綠色陰影區(qū)域),此參數(shù)即為系統(tǒng)的臨界參數(shù).下面考察式(33)和式(34)中截斷參數(shù)P和N對本文解收斂性的影響.圖3(b)給出了不同rH時截斷階數(shù)P對結(jié)果收斂性的影響.由圖3 可知,選擇P=N=7,M=200 可以保證本文計算結(jié)果的收斂性.圖4 給出了不同rH對應(yīng)的臨界參數(shù)αc.與預(yù)期一致,αc隨rH的增大而增大并最終穩(wěn)定在αc?1.84.

圖3 Fig.3

圖5 給出了不同rH下壁板的失穩(wěn)模態(tài)ψ(ξ).從圖可知,盡管不同rH下的臨界動壓不同,但失穩(wěn)模式都以倒置懸臂梁一階模態(tài)φ1為主導(dǎo).考察失穩(wěn)模態(tài)

圖4 本文理論解與其他理論及數(shù)值結(jié)果的對比分析Fig.4 Comparison of the present theoretical results with other existing theory and the numerical method

圖5 系統(tǒng)的失穩(wěn)模態(tài)(上)及壓強系數(shù)(下)Fig.5 The instability model(top)and the pressure coefficient(bottom)of the plate

下壁板表面氣動壓力的分布,定義壓強系數(shù)

圖5 給出了不同rH對應(yīng)的壓力系數(shù)(由于前緣壓差為無窮大,圖中已做截斷處理).由圖可知,壁板前緣附近為正壓區(qū)而尾緣附近為負(fù)壓區(qū);隨rH的減小正壓區(qū)擴大且壓力值增加,而負(fù)壓區(qū)縮小且壓力值減小.壁板上合力增大且力矩中心前移,對應(yīng)的臨界失穩(wěn)動壓也就越低.

為驗證本文的理論結(jié)果,考慮無壁面情況,此時滿足Kutta-Joukowski 條件的壁板的壓強解(Theodorsen 解[28])為

其中

取倒置懸臂梁模態(tài)函數(shù)φk,(k=1,2,...,K),應(yīng)用伽遼金方法求解方程(1),通過特征值分析獲得系統(tǒng)臨界參數(shù).如圖3(b)中虛線所示,選取前三階模態(tài)便可得收斂結(jié)果αc?1.8464.本文理論結(jié)果與經(jīng)典Theodorsen 解得到的結(jié)果吻合很好.

為了進一步考察本文理論方法對rH的適用性,采用如圖1(b)所示的離散面元模型進行數(shù)值求解.壁板上的面元分布由線性分布的連續(xù)渦來表示以計算滿足邊界條件的流體勢函數(shù),在質(zhì)量點處采用有限差分方法求解壁板結(jié)構(gòu)運動方程,詳細(xì)的離散和計算過程請參見文獻[12].本文采用該離散模型計算獲得了不同rH時的臨界參數(shù),結(jié)果如圖4 所示.由圖可知本文理論與數(shù)值結(jié)果的變化趨勢一致且吻合較好(兩者相對誤差保持在5%之內(nèi),圖中綠色區(qū)域內(nèi)).上述對比結(jié)果均表明了本文氣動力理論及失穩(wěn)分析方法的有效性.

5 風(fēng)洞實驗驗證

采用風(fēng)洞實驗研究壁板失穩(wěn)的臨界速度并驗證本文理論的正確性.實驗中使用直徑1.2 m 的開口式風(fēng)洞,其具有5~40 m/s 的風(fēng)速調(diào)節(jié)范圍和低于0.3% 的湍流度.實驗中壁板模型采用長寬比為2/3(20 cm×30 cm)的鍍鋅鐵板(ρp=7.85 g/cm3,Ep=210 GPa).如圖6 所示,壁板豎直安裝,并用兩塊厚鋼板將其夾持固定在剛性支架上以保證固支約束并進行實驗,夾持段壁板的長度約為3 cm;壁板一側(cè)豎直放置厚木板以模擬壁面對氣流的約束限制.靠近壁板固定端的根部位置貼有應(yīng)變片,實驗開始前采用敲擊法測試了壁板在無風(fēng)狀態(tài)下的自振頻率(7.9 Hz),并與理論結(jié)果(8.2 Hz)進行了對比,檢驗了固支約束實現(xiàn)及測試設(shè)備的可靠性.實驗中緩慢增加風(fēng)速至目標(biāo)風(fēng)速,穩(wěn)定至少一分鐘后進行數(shù)據(jù)采集,以保證實驗數(shù)據(jù)的可靠性.

圖6 風(fēng)洞中模型安裝圖和不同風(fēng)速下細(xì)繩的狀態(tài)圖Fig.6 Pictures of the setup of the model in the wind tunnel,and of the string states with different flow velocities

靜態(tài)失穩(wěn)不能像顫振這類動態(tài)失穩(wěn)一樣可以依據(jù)信號等幅周期變化的特征而直接判定.本文借鑒材料力學(xué)壓桿失穩(wěn)原理設(shè)計了一種等效拉力測試方法,測試原理如圖7(a)所示.理論上給定壁板初始變形量為失穩(wěn)模態(tài)的任意小倍數(shù),當(dāng)氣流速度小于臨界值時,壁板需要外部作用(實驗中依靠張緊的細(xì)繩提供拉力)才能維持該狀態(tài);而當(dāng)流速達到臨界值時,則可不依靠任何外部作用.因此可由是否需要細(xì)繩提供拉力而維持給定的初始變形來等效判定系統(tǒng)是否達到了臨界狀態(tài).本文的理論分析表明壁板將呈現(xiàn)一階懸臂梁模態(tài)失穩(wěn)模式,而該模態(tài)形式的初始壁板靜變形在實驗中卻不易精確給出.依據(jù)文獻[26]給出的懸臂梁一階模態(tài)與其端部受集中載荷產(chǎn)生的靜變形形態(tài)相類似這一結(jié)論,本文實驗采用在壁板自由端施加集中力(張緊細(xì)繩使其具有一定預(yù)拉力)而實現(xiàn)壁板的初始變形.

圖7 Fig.7

如圖6 及圖7(a)所示,懸臂板緊靠前緣正中心位置連接有細(xì)棉繩,細(xì)繩的另一端連接在拉力傳感器上(量程5 N,靈敏度0.01 N).調(diào)整細(xì)繩的長度可使壁板前緣產(chǎn)生不同的初始撓度Δ,其在傳感器上表現(xiàn)為不同的預(yù)拉力值.當(dāng)流速小于臨界值時,細(xì)繩會處于張緊狀態(tài)并為壁板提供拉力而維持其初始變形,如圖6(a)所示;而當(dāng)流速接近于臨界值時,細(xì)繩則會處于松弛狀態(tài),如圖6(b)所示.注意到細(xì)繩僅受氣流作用也會導(dǎo)致拉力(Fw)的存在,傳感器的實測拉力(F)則包含(Fw)及壁板對繩子的拉力(Fp).圖7(b)給出了風(fēng)速12 m/s 時時間1 s 內(nèi)兩種拉力的實測值.測試?yán)w時,將細(xì)繩保持自然狀態(tài),兩端分別固定在剛性支架及夾持端并單獨置于風(fēng)洞中進行測試.從圖7(b)可知,雖然兩種拉力均表現(xiàn)出明顯的波動性,但Fw的變化幅值較小(±0.005 N 之內(nèi)).若考慮兩種拉力之間的弱關(guān)聯(lián)性則可由下式

作為臨界狀態(tài)的近似判定條件.

理論上,對于任意的自由端部撓度Δ,實驗中均應(yīng)測試得到相同的臨界風(fēng)速.然而實驗中發(fā)現(xiàn),較小的值會因較差的抗干擾性而導(dǎo)致測量精度不易滿足,而較大的Δ 則會導(dǎo)致壁板產(chǎn)生幾何大變形非線性[26].因此,實驗中首先以無壁面的情形進行多組重復(fù)測試,通過對比已有理論及計算結(jié)果來確定最佳的Δ 給定范圍.圖8(a)給出了不同Δ 值下F及Fw在不同風(fēng)速下的測試值.依據(jù)式(40)可判定當(dāng)拉力Fm曲線處于圖8(a)的綠色陰影區(qū)域時,系統(tǒng)處于失穩(wěn)狀態(tài).依據(jù)圖8(a)得到臨界速度Ucr與Δ 之間的變化關(guān)系如圖8(b)所示.通過與理論值(39)及數(shù)值結(jié)果(圖4)對比可知,當(dāng)Δ 在0.3~0.5 cm 之間取值時,實驗結(jié)果與理論及數(shù)值結(jié)果吻合較好(實驗值與理論值(39)相差2%之內(nèi),圖中綠色區(qū)域內(nèi)),因此本文實驗選取該區(qū)間的值作為拉出撓度值完成了不同rH下的實驗,結(jié)果如圖9 所示.圖9(a)中無風(fēng)時細(xì)繩中的初始拉力并不完全相等,這是由于為了保證數(shù)據(jù)的多樣性而在0.3~0.5 cm 內(nèi)給定不同的Δ 值而導(dǎo)致的.

圖8 Fig.8

圖9 Fig.9

圖10 本文理論結(jié)果與試驗結(jié)果及已有理論結(jié)果的對比Fig.10 Comparison of the present theoretical results with the experiment and the other existing theory

圖10 給出了本文理論與試驗結(jié)果的對比,由圖可知兩者吻合很好,各風(fēng)速下的實驗結(jié)果與理論結(jié)果之間的相對誤差均保持在2%之內(nèi)(綠色區(qū)域內(nèi)),這充分表明了本文的理論計算及風(fēng)洞試驗的有效性和準(zhǔn)確性.

6 結(jié)語

本文考慮壁面效應(yīng)對亞音速氣流中倒置懸臂壁板的靜態(tài)失穩(wěn)問題進行了理論及實驗研究.文中采用鏡像函數(shù)法描述壁面約束條件,基于算子理論對壁板上的氣動力進行了理論分析,獲得了以Possio 積分方程為表征的氣動力表達式,壁面效應(yīng)則表征為一包含T,Q 及H 算子的復(fù)合算子;壁板的失穩(wěn)方程轉(zhuǎn)化為了定區(qū)間上的函數(shù)逼近問題,并利用Wererstrass 定理及最小二乘方法得到了最優(yōu)逼近函數(shù)確定臨界動壓,無需進行特征值計算;給出了壁面效應(yīng)對失穩(wěn)的定量影響并對壁面影響失穩(wěn)的原理進行了探討;設(shè)計了靜態(tài)失穩(wěn)測試方法,其相比于動態(tài)實驗對風(fēng)洞品質(zhì)要求更低且魯棒性較高;風(fēng)洞實驗分析結(jié)果與本文理論結(jié)果吻合較好,驗證了本文理論方法的正確性.

本文中以Possio 積分方程表征的氣動力物理含義清晰,依據(jù)微分算子理論將失穩(wěn)問題轉(zhuǎn)化為了定區(qū)間上的函數(shù)逼近問題,為該類問題的求解提供了新的思路.值得指出的是,本文方法在壁板顫振問題中的推廣及應(yīng)用正是下一步計劃開展的工作.

附錄:算子QH 的交換性

由本文算子的定義可知

其中

分別令t=τ+p及t=x+q,則上式變作

上式中的積分為瑕積分,注意到c>0,利用圍道積分進行計算后可得

猜你喜歡
理論實驗
記一次有趣的實驗
堅持理論創(chuàng)新
微型實驗里看“燃燒”
神秘的混沌理論
理論創(chuàng)新 引領(lǐng)百年
相關(guān)于撓理論的Baer模
做個怪怪長實驗
NO與NO2相互轉(zhuǎn)化實驗的改進
實踐十號上的19項實驗
太空探索(2016年5期)2016-07-12 15:17:55
理論宣講如何答疑解惑
主站蜘蛛池模板: 久久人人97超碰人人澡爱香蕉 | 国产极品美女在线播放| 国产精品久久久久久久久久久久| 国产欧美高清| 欧美日韩免费观看| 在线综合亚洲欧美网站| 91小视频版在线观看www| 88av在线| 亚洲精品成人片在线观看| 久久亚洲美女精品国产精品| 久久精品最新免费国产成人| 国产亚洲精品无码专| 精品国产电影久久九九| 国产成人精品在线1区| 亚洲精品国偷自产在线91正片| 国产日韩丝袜一二三区| 国产精品不卡永久免费| 九九热视频在线免费观看| 色爽网免费视频| 国产乱人乱偷精品视频a人人澡| 亚洲日韩AV无码精品| 欧美在线天堂| 香蕉久人久人青草青草| 狠狠五月天中文字幕| 91在线国内在线播放老师| 全部无卡免费的毛片在线看| 福利片91| 日本日韩欧美| 国产精品制服| 国产日本视频91| 尤物亚洲最大AV无码网站| 极品av一区二区| 国精品91人妻无码一区二区三区| 日韩黄色大片免费看| 在线观看免费国产| 91无码视频在线观看| 亚洲精品在线观看91| 国产成人综合久久精品下载| 国产成人资源| 天天综合色网| 亚洲首页在线观看| 亚洲欧洲日韩综合色天使| 伊人久综合| 爱做久久久久久| 国产乱人伦偷精品视频AAA| 亚洲一区免费看| 久久久久青草线综合超碰| 亚洲国产成熟视频在线多多| 国产在线拍偷自揄观看视频网站| 欧美日本在线一区二区三区| 四虎影视库国产精品一区| 一区二区三区国产精品视频| 91在线精品麻豆欧美在线| 精品一区二区三区波多野结衣 | 亚洲AV无码不卡无码 | 天天操天天噜| 经典三级久久| 性色生活片在线观看| 久久综合亚洲鲁鲁九月天| 欧洲一区二区三区无码| 在线欧美国产| 久久综合五月婷婷| 在线观看无码a∨| 国产精品一区二区在线播放| 99在线国产| 免费看a级毛片| 欧美精品亚洲精品日韩专区| 日韩精品资源| 亚洲中文字幕国产av| 91久久精品日日躁夜夜躁欧美| 亚洲成人一区二区三区| 亚洲无码熟妇人妻AV在线| 成人午夜视频免费看欧美| 99热这里只有精品免费国产| 亚洲无限乱码一二三四区| 亚洲视频一区在线| 国产精品思思热在线| 欧美精品在线视频观看| 国产第二十一页| 99er精品视频| 无码乱人伦一区二区亚洲一| 亚洲无码高清免费视频亚洲|