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

熱聲載荷作用下薄壁結構非線性響應分析和試驗驗證

2019-04-22 10:45:22沙云東張墨涵趙奉同朱付磊
航空學報 2019年4期
關鍵詞:模態振動結構

沙云東,張墨涵,趙奉同,2,*,朱付磊

1. 沈陽航空航天大學 遼寧省航空推進系統先進測試技術重點實驗室,沈陽 110136 2. 北京航空航天大學 能源與動力工程學院,北京 100083

熱聲載荷作用下結構動力學響應特性及疲勞預測是現代高速飛行器結構設計和安全性的重要因素?;鹧嫱脖诎濉l動機機匣、渦輪葉片和飛行器蒙皮等都會受到高溫及聲壓的影響。飛行器熱保護系統(TPS)表面溫度可達1 600 ℃以上,聲壓級(SPL)可達180 dB[1]。TPS在高聲壓級下會產生較大位移,并可能在高溫下出現屈曲。為提高航空飛行器的飛行性能即解決嚴重的氣動、聲學和熱環境下的疲勞問題,應統計大量薄壁結構呈現出的非線性響應[2-6],并進行疲勞壽命的預測[7-9]。因此,為了滿足航空薄壁結構的結構設計要求、降低研究成本,有效的熱聲動態響應分析成為當前的重要任務。

薄壁結構仿真對提高試驗的有效性和可靠性起重要作用。結構非線性響應模擬分析方法主要有:攝動法、 Fokker-Plank-Kolmogorov(FPK)方程、von Karman-Herrmann大撓度板方程、蒙特卡羅法、等價線性化(EL)方法、降階模型法(ROM)和有限元數值積分方法(FEM)。已證明攝動法[10]僅限于較弱的幾何非線性,而FPK方程僅對單自由度(DOF)系統有精確解。蒙特卡羅模擬[11-12]是最常用的方法,但偏微分方程和伽遼金法(Galerkin)僅適用于較簡單的結構[13-14]。等價線性化方法已被廣泛應用,因為該法能夠準確獲取大范圍的響應統計,同時保持較輕的計算負擔[15]。其缺點是假設的響應必須是高斯分布。對于有限元數值積分方法[16],計算成本是主要關注的因素。有限元模型通常包含數百個物理自由度,且非線性項在每個時間步都會更新及重新組合。Dhainaut等[17]提出了一個用于預測受熱聲載荷作用下薄板的非線性隨機響應的有限元公式。為了驗證公式,模擬的數值結果包括最大位移和應變的均方根值、位移和應變響應的時間歷程、概率分布函數及功率譜密度等。楊智春等[18]通過分析壁板跳變運動產生的機理定義了“穿零頻次”,將跳變運動定量進行分類,得出了在顫振臨界動壓前,動壓的增大會導致薄壁板勢能阱深度呈現先增大后減小趨勢的結論。楊雄偉等[19]使用混合有限元-統計能量方法分析了高溫環境、165 dB聲壓作用下的結構聲振特性,總結了熱應力對結構模態頻率的影響規律,以及材料物性對模態頻率與振動響應的作用。桂業偉等[20]從單向及雙向耦合兩方面研究了流-熱-固多場耦合問題,對國內研發的熱響應耦合分析平臺進行總結并討論了未來的發展趨勢。文獻[17-20]均以理論公式及仿真計算為主體,并未有試驗驗證作為支撐。Liu等[21]通過對比試驗測得的后屈曲偏轉與有限元模擬后屈曲分析得到的結果,由一階模態方程推導出了動態響應參數的變化,從而描述了從線性振動到完全非線性動態響應的轉變,得出了聲激勵、熱應力和結構剛度的組合對動態響應的影響。劉磊等[22]結合熱-流-固耦合試驗需求開展了風洞試驗,通過計算對比分析驗證了多場耦合計算平臺方法的有效性,并預估了高溫合金模型的試驗結果。侯薇等[23]建立了基于截面平均的準一維非線性熱聲模型,并進行數值求解,采用寬頻時域聲阻抗邊界條件,預測了臨界起振溫度及熱聲系統的瞬態壓力波。以上文獻[21-23]均為在同一模型上進行了仿真與試驗的對比及預測,而本文是在仿真與試驗對比的基礎上做了適當的拓展仿真,即在薄壁板的基礎上做了加筋薄板的仿真。

基于文獻[24],本文首先采用了一種新型航空發動機材料GH188,對熱環境下薄壁方板進行聲激振試驗,得到了加速度響應和應變響應的結果。然后,根據馮卡門大撓度理論,采用耦合的FEM/ROM方法對薄壁板的動力學響應進行計算。再將熱聲試驗測得的響應結果與數值仿真結果進行對比,驗證了熱聲仿真計算分析方法的精度與有效性。最后,利用該計算模型計算了在不同溫度場及聲載荷作用下四邊固支交叉加筋金屬板的動力學響應,并進行了詳細的對比與分析。

1 非線性響應理論

1.1 平板結構有限元控制方程

根據有限元法對各單元振動方程進行推導和求解運算,得到整體結構控制方程:

(1)

式中:M為質量矩陣;K為剛度矩陣;W為整體位移;FT為合外力。

設時間步長為eiωt,i為步數,ω為頻率,t為時間,則結構振動方程可記為

(2)

式中:Φn為結構法向振型;ωn為結構基頻,分析阻尼對結構響應特性的影響,將有限元與邊界元方法進行耦合分析,得到結構在模態坐標下的運動方程為

(3)

其中:dn為模態位移;ζn與Mn分別為阻尼系數和模態質量。對結構模態坐標運動方程式(3)進行推導,可以獲得響應函數表達式為

(4)

將式(3)與式(4)結合,則式(4)轉化為

HSd=ΦΤFT

(5)

式中:HS為響應函數;d為模態位移;Φ為模態矩陣。

將模態載荷FT分解為外力載荷Fe和聲壓載荷FP兩部分,控制方程變為

HSd=ΦΤFP+ΦΤFe

(6)

式中:ΦΤFP為聲壓向量,從邊界元系統傳遞到到有限元中,傳遞過程如圖1所示。

通過式(6),得到平板結構控制方程為

HSd=(LT)ΤFP+ΦΤFe

(7)

式中:LT為傳遞矩陣。

圖1 有限元與邊界元傳遞過程Fig.1 Transfer process of finite element and boundary element

1.2 薄壁構件大撓度方程

薄壁構件在熱聲載荷作用下大撓度方程的推導是在克?;舴蚍e分方程和馮卡門大撓度方程2個 理論假設的基礎上進行的。以薄壁結構中截面位置為參考,推導結構上任意位置處各方向應變的表達式為

(8)

針對結構應變和位移響應進行矢量微分計算,得到帶有撓度項的變形協調方程:

(9)

將艾里應力函數F引入,式(8)可表示為式(10) 的形式

(10)

將式(9)與式(10)聯立,可以轉化成式(11)中帶有應力函數項的變形協調方程:

(11)

(12)

式中:ρ為密度;ξ為阻尼系數;υ為泊松比;p(x,y,t)為模擬聲載荷的隨機應力;D為彎曲剛度;θ為板厚的溫度梯度。

1.3 薄壁構件熱屈曲理論

隨著溫度的逐漸升高,薄壁結構因受熱而產生的壓應力會引起熱屈曲,導致結構變形并喪失承載能力,此時的溫度定義為臨界屈曲溫度,用T*表示。四邊固支板的一階臨界屈曲溫度表示為

(13)

式中:β=b/a為薄壁構件寬與長之比。將溫度無量綱化,表示為:T0=T0/T*,T0為高出常溫的那部分溫度,故

T=T0+δvT0fv(x,y)+δgT0fg(x,y)

(14)

式中:T為任意一點溫度;δvT0和δgT0分別為整個板面上溫度變化值的幅值大小和沿板厚度變化的溫度的幅值;fv(x,y)和fg(x,y)分別為薄板面上溫度變化情況和任一點(x,y)處溫度變化情況。

2 薄壁板熱聲載荷試驗驗證

2.1 薄壁板熱聲試驗

該試驗是針對GH188板材開展高溫環境下聲激振的動力學響應試驗。試驗件的幾何尺寸與應變片的貼片位置如圖2所示,其中板厚1.5 mm。1# 與3#應變片在薄壁板的短邊中點,以測量試驗件x方向的應變;2#與4#應變片在薄壁板的長邊中點,以測量試驗件y方向的應變。不同溫度下的材料參數如表1所示,K為熱導率。通過口框夾具將試驗件壓緊,并采用雙排螺栓擰緊的方式對試驗件四周固定,以實現四邊固支約束。采用雙面加溫的方式對試驗件進行加熱如圖3所示,溫度載荷為50~250 ℃,間隔為50 ℃;采用聲場控制方法施加有限帶寬高斯白噪聲,頻率范圍為100~1 250 Hz,加載的總聲壓級范圍為142~154 dB,間隔為3 dB。

本次試驗有2個目的,一是對變形高溫合金材料隨機疲勞性能測試試驗系統開展聯合調試,驗證系統的有效性;二是在給定的不同溫度及聲載荷條件下獲取變形高溫合金材料結構表面指定位置處高溫聲激振測試的動態應變響應。

圖2 試驗件尺寸及應變片位置Fig.2 Sizes of test piece and locations of strain gage

表1 GH188材料參數Table 1 Material parameters of GH188

ParameterT=50 ℃T=400 ℃T=600 ℃E/GPa208175156υ0.3010.3180.326α/(10-6 ℃-1)11.713.414.4ρ/(103kg·m-3)9.099.099.09K/(W·℃-1)15.719.627.1

圖3 試驗件加熱現場照片Fig.3 Photo of test piece heated on site

試驗結果表明,試驗系統能滿足試驗要求,且能實現噪聲載荷和熱載荷的聯合加載,并能給出有效測試結果。在給定的頻帶范圍內聲壓級可達157 dB,試驗件表面溫度可達500 ℃。同時也分析總結出了結構動態響應特征以及隨溫度和聲載荷變化規律。得出了結構屈曲前和屈曲后過程響應數據,為進行仿真對比驗證提供了有效數據。

2.2 仿真對比驗證

為了確保仿真的準確性,仿真對象的特性、邊界條件和熱聲載荷等均與試驗相同。同時采用FEM/ROM方法計算薄壁板的熱聲動力學響應,并提取與試驗相同位置的應變響應,進行仿真與試驗的對比驗證。

在試驗過程中,利用激光測振儀獲得了不同熱聲載荷組合作用下薄壁板中心位置的加速度響應,如圖4所示。通過對不同溫度下薄壁板的一階模態頻率分析可以得出加速度響應和應變響應。本項試驗對象為GH188板材樣件,選用BAB350-5AA250(11)-150型中溫應變計,該應變計為自補償應變計,可自動修正熱輸出,無需在結果上再進行修正。表2為不同溫度下結構基頻試驗結果,由于高階響應幅值遠低于基頻響應,所以只給出了試驗基頻。表3列出了仿真的前六階模態頻率,可以看出,100 ℃為薄壁板的臨界屈曲溫度。對比結果表明兩者的一階模態頻率具有一致性,均呈現出先降低后升高的趨勢。對于薄壁板在前屈曲區域保持軟化,并在后屈曲區域轉化為硬化得到了驗證。圖5表示試驗測試的聲壓級為151 dB時,在1#測點的代表性加速度響應譜,在前屈曲區域T=50 ℃時,一階模態頻率為347 Hz,而在后屈曲區域T=150 ℃時,一階模態頻率為306 Hz。圖6為薄壁板在SPL=151 dB時T=50,100,150,200 ℃的應力響應云圖,可看出各溫度下薄壁板的危險點位置,反映了結構整體響應特性。

圖4 試驗中的加速度測量Fig.4 Acceleration measurement in test

聲壓級為145 dB和151 dB時,應變片測量了與試驗件長邊和短邊中點的基頻對應的應變響應值,表4所列的是試驗與仿真的長邊中點和短邊中點應變響應值。結果表明,仿真應變響應值與試驗結果非常吻合,驗證了薄壁板在熱聲載荷作用下的非線性動力學響應模型和計算方法的有效性。圖7給出了試驗測試的短邊中點(1#)聲壓級為151 dB時,溫度載荷為50 ℃與150 ℃這2組 代表性應變響應譜,屈曲前T=50 ℃時,應變響應值為10.7×10-6,屈曲后T=150 ℃時,應變響應值為16.6×10-6。

表2 不同溫度下結構基頻試驗結果

表3 不同溫度下結構前六階熱模態頻率計算結果Table 3 Calculation results of the first six thermal modal frequencies of structures at different temperatures

圖5 薄壁板中點的加速度響應Fig.5 Acceleration response on thin-walled panel midpoints

圖6 薄壁板應力響應云圖Fig.6 Stress response contours of thin-walled panel

綜上所述,通過本次對高溫合金薄壁板進行高溫環境聲激振試驗和數值仿真,分別對比了一階模態頻率和測點的應變響應值,充分驗證了熱聲載荷作用下薄壁板非線性動力學響應模型和計算方法的有效性。

表4 測點位置單元在基頻處的單向應變Table 4 One-dimensional strain at fundamental frequency of measuring point unit

圖7 試驗件應變頻域響應Fig.7 Frequency-domain response of test piece

3 熱聲載荷條件下非線性響應分析

數值仿真計算對象選擇圖8所示的四邊固支交叉加筋板,其性能與試驗的薄壁板保持一致,表5 列出了特定的幾何參數,其中h1為肋寬,h2為肋高。假設熱載荷均勻分布在板的表面,聲載荷是有限帶寬高斯白噪聲,帶寬頻率范圍為0~2 000 Hz。本文主要提取模型中間節點處的熱聲響應結果以分析其振動特性,如圖8所示。圖9為結構的有限元網格劃分。

圖8 交叉加筋板示意圖Fig.8 Cross stiffen-reinforced plate sketch

表5 幾何參數Table 5 Geometric parameters

Parametera/mmb/mmh/mmh1/ mmh2/ mmValue2101451.51.51.5

圖9 加筋板的有限元網格模型Fig.9 Finite element mesh model of cross reinforced plate

采用有限元法計算出加筋薄壁板的第1階臨界屈曲溫度Tcr=27.5 ℃。為便于表達,本文均以熱屈曲系數S來表示溫度,其中S=T/Tcr,T為實際溫度。以熱屈曲系數S的大小為標準,有3種典型狀態,S<1、S=1和S>1分別代表屈曲前、臨界屈曲和屈曲后狀態。常溫狀態下(20 ℃)加筋板的前六階模態頻率如表6所示。本文在不同S與SPL組合作用下獲得了加筋板的動力學響應,其中S為0~20,間隔為0.1,SPL的范圍為145~172 dB,間隔為3 dB。

表6常溫環境下加筋板前六階模態頻率

Table6Firstsixordermodalfrequenciesofstiffen-reinforcedplateinnormaltemperatureenvironment

Order123456Frequency/Hz490.36778.591273.61310.51509.81886.3

3.1 橫向位移響應分析

室溫條件下(S=0.2),加筋板圍繞初始平衡位置做隨機振動,隨著聲壓級由160 dB升高到172 dB, 振動響應加劇,位移響應幅值由0.5 mm增大到2.1 mm,如圖10(a)、(d)所示。當溫度升高到S=0.6時,位移響應幅值隨之增大,由0.5 mm增大到1.0 mm和由2.1 mm增大到2.3 mm,如圖10(b)、 (e)所示。當溫度達到臨界熱屈曲溫度S=1.0 時,位移響應幅值明顯增大,由1.0 mm增大到1.3 mm和由2.3 mm增大到2.6 mm,如圖10(c)、 (f)所示。圖10中清晰地顯示了屈曲前,隨著溫度的升高,位移時間歷程越來越稀疏;且固定條件下,位移響應與振動響應基本相同,不同振動方向的幅值基本一致。

圖11為加筋板在聲壓級為172 dB條件下屈曲后的橫向時間位移響應,其位移響應的運動趨勢取決于熱載荷與聲載荷的相對強度。當加載溫度稍高于臨界屈曲溫度時,聲載荷強度高于熱載荷強度,聲載荷足以克服熱載荷下剛度區域的驅離力,此時聲載荷大小起主導作用,加筋板位移時間響應在屈曲后的2個平衡位置間做連續跳變運動,如圖11(a)所示。隨著溫度逐漸升高,當聲載荷強度與熱載荷強度相當時,聲載荷剛好可以克服熱載荷下剛度區域的驅離力,此時聲載荷與熱載荷作用大小相當,位移時間響應在屈曲后的2個 平衡位置間做間歇跳變運動,如圖11(b)、(c)所示。當溫度繼續升高,聲載荷強度低于熱載荷強度時,聲載荷無法克服熱載荷下剛度區域的驅離力,此時熱載荷大小起主導作用,位移時間響應圍繞屈曲后的一個凸面平衡位置做非線性的隨機振動,如圖11(d)所示。圖11顯示了屈曲后,位移時間歷程由稀疏變得緊密;且固定條件下,位移響應與振動響應有較大差別,不同振動方向振幅相差也很大。

圖10 屈曲前到臨界屈曲不同聲壓級下橫向位移時間歷程Fig.10 Pre-buckling to critical buckling time history of lateral displacements under different SPLs

圖11 屈曲后橫向位移時間歷程(SPL=172 dB)Fig.11 Time history of post-buckling lateral displacements (SPL=172 dB)

如圖12所示,屈曲前,結構橫向位移的概率密度函數(PDF)基本服從正態分布,如圖12(a)、 (b)所示。隨熱載荷的增大,結構橫向位移PDF不再服從正態分布,非線性響應逐漸增加。在臨界屈曲溫度附近時,非線性響應明顯增強,加筋板開始在圖12(c)所示的2個振動平衡位置之間進行往復振動。屈曲后,橫向位移PDF明顯呈現雙峰現象,如圖12(d)所示,可以得出,加筋板在凸面平衡位置周圍的振動概率比在凹面平衡位置附近的振動概率更大,表明了聲載荷的增加會在凸面平衡位置較好地驅動勢能阱深度的增加。因此,加筋板更可能在凸面平衡位置周圍進行隨機振動。當熱載荷與聲載荷相當時,由于熱載荷的加劇,凹面平衡位置中的勢能阱深度增大,使加筋板在凹面平衡位置附近振動的可能性更大,如圖12(e) 所示。隨著溫度持續升高,由于熱載荷的增加,勢能阱深度更大,從而引起振動響應平均值增加而振幅減小,并在一個屈曲后平衡位置附近產生非線性振動,如圖12(f)所示,橫向位移的PDF逐漸趨于正態分布,而非線性響應減弱。

圖13給出了屈曲前和屈曲后結構橫向位移的功率譜密度(PSD)。分析表明,在前屈曲區域,隨著溫度的升高加筋板的PSD幅值增大,而基頻降低,即呈現峰值左移現象;PSD出現高頻響應峰值,峰值隨著熱載荷的增加而增加,在臨界屈曲附近達到最大值,如圖13(a)所示。顯然,在屈曲前,由于熱應力的增加,加筋板進入軟化區域,剛度降低,導致振動響應增大。然而,與屈曲前的響應結果相比,屈曲后PSD出現了低頻響應峰值,而且隨著熱載荷的增加,響應峰值降低,加筋板的基頻升高,如圖13(b)所示。這表明在屈曲后,由于熱應力的增加,加筋板進入硬化區域且剛度增強,因此振動響應也隨之減小。

圖12 橫向位移的概率密度函數(PDF)(SPL=172 dB)Fig.12 Probability Density Function (PDF) of lateral displacement(SPL=172 dB)

圖13 橫向位移功率譜密度(PSD)(SPL=151 dB)Fig.13 Power Spectral Density (PSD) of lateral displacement (SPL=151 dB)

3.2 應力分量響應分析

圖14顯示了在加筋板中心位置提取的x方向應力分量的時間歷程。它具有以下幾點特征:

1)x方向應力分量時間歷程的振動平衡位置逐漸偏離橫軸。從圖14(a)可見,在臨界屈曲溫度條件下,平衡位置偏離值為y=-20 MPa,且溫度越高,偏移量越大。拉應力和壓應力同時存在并交替作用在加筋板上,而此時壓應力為主導因素。

2) 如圖14(b)所示,加筋板在x方向應力響應做間歇跳變運動過程中,當加筋板在凸面平衡位置周圍振動時,x方向應力分量的振幅較大而平均值較??;而在凹面平衡位置周圍振動時,x方向應力分量的振幅較小而平均值較大。如圖14(b)、 (c)所示,隨著溫度的升高,x方向應力分量在2個平衡位置處的振幅減小而平均值增大。從S=1.5到S=2.2,2 個平衡位置處的振幅最大值從約為130 MPa減小到約為90 MPa;平均值從約為90 MPa增大到160 MPa。

3) 圖14(d)表明,當熱載荷足夠大時,加筋板在凸面平衡位置附近振動,使x方向應力的幅值顯著降低而平均值增大,此時加筋板主要表現出拉應力。

圖14 屈曲后x方向應力時間歷程(SPL=172 dB)Fig.14 Time history of post-buckling x-direction stress(SPL=172 dB)

4 結 論

1) 選擇四邊固支薄壁金屬板進行熱聲載荷作用下動力學響應試驗,并比較仿真與試驗的熱模態頻率,證實了結果的一致性。然后比較了應變響應值,結果表明應變響應的仿真值與試驗值在相同的數量級下比較吻合,得以驗證了在熱聲載荷環境下薄壁板的非線性動力學響應模型和計算方法的有效性。

2) 本文完成了不同熱聲載荷大小組合情況下交叉加筋板的動力學響應計算。通過分析橫向位移響應,可以得出以下結論:由于熱應力的增大,在屈曲前,基頻在軟化區域逐漸降低,由于響應的高斯分布,SPL=172 dB時位移峰值由2.1 mm逐漸增大到2.6 mm,使加筋板在初始平衡位置周圍隨機振動;在屈曲后,基頻在硬化區域逐漸增加,加筋板產生了由熱載荷與聲載荷的相對強度決定的跳變運動,同時考慮到響應的非高斯分布,加筋板顯示出具有2個勢能阱的雙峰現象。

3) 當加筋板在2個平衡位置間振動時,拉應力和壓應力交替作用,S=2.2, SPL=172 dB,結構x方向應力分量由拉應力(45~255 MPa)跳變為壓應力(150~300 MPa),x方向應力分量顯示加筋板跳變形態。當受到較強的熱載荷時,加筋板將圍繞2個平衡位置中的一個振動,S=2.8, SPL=172 dB,x方向應力分量為50~250 MPa,加筋板完全進入受拉區域。

航空航天薄壁結構熱聲載荷作用下非線性響應分析十分復雜,本文只是初步探究,需要大量理論、模擬,特別需要進行針對更多典型結構的復雜熱聲載荷試驗,獲取大量試驗數據,為進一步工程應用提供參考依據。

猜你喜歡
模態振動結構
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
振動與頻率
天天愛科學(2020年6期)2020-09-10 07:22:44
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
中立型Emden-Fowler微分方程的振動性
論《日出》的結構
國內多模態教學研究回顧與展望
創新治理結構促進中小企業持續成長
現代企業(2015年9期)2015-02-28 18:56:50
基于HHT和Prony算法的電力系統低頻振蕩模態識別
UF6振動激發態分子的振動-振動馳豫
計算物理(2014年2期)2014-03-11 17:01:44
主站蜘蛛池模板: 久综合日韩| Jizz国产色系免费| 福利在线免费视频| 欧美精品H在线播放| 无码免费视频| 人妖无码第一页| 重口调教一区二区视频| 精品精品国产高清A毛片| 欧美日韩精品一区二区在线线| 亚洲精品制服丝袜二区| 国产99热| 色婷婷久久| 亚洲精品无码抽插日韩| 亚洲无码免费黄色网址| 精品剧情v国产在线观看| 久久久久久久蜜桃| 日本精品视频| 国产女人综合久久精品视| 日本成人不卡视频| 一级毛片视频免费| 国产日产欧美精品| 欧美精品影院| 狂欢视频在线观看不卡| 毛片手机在线看| 影音先锋亚洲无码| 日韩最新中文字幕| 国产福利2021最新在线观看| 激情在线网| 日韩福利在线观看| 午夜精品区| 国产精品美女免费视频大全| 婷婷亚洲天堂| 日韩无码视频专区| 5555国产在线观看| 欧美中文字幕在线二区| 久久久久久久久久国产精品| 国产91全国探花系列在线播放| 一区二区在线视频免费观看| 99精品在线视频观看| 日韩欧美国产精品| 色婷婷色丁香| 欧美日韩亚洲国产主播第一区| 日本欧美在线观看| 91在线精品免费免费播放| 亚洲国产在一区二区三区| 欧美三级视频在线播放| 91精品国产自产在线老师啪l| 国产一级视频在线观看网站| 97av视频在线观看| jizz国产在线| a色毛片免费视频| 97视频免费在线观看| 三上悠亚一区二区| 色男人的天堂久久综合| 国产高清国内精品福利| 精品国产美女福到在线不卡f| 国产午夜无码专区喷水| 粗大猛烈进出高潮视频无码| 国产精品高清国产三级囯产AV| 亚洲欧美日韩成人在线| av无码一区二区三区在线| 国内精品一区二区在线观看| 国产成人高清精品免费5388| 5555国产在线观看| 亚洲区欧美区| 久久综合丝袜长腿丝袜| 亚洲国产欧洲精品路线久久| 无码电影在线观看| 中文字幕啪啪| 国产精品流白浆在线观看| 欧美一级在线播放| 女人18一级毛片免费观看| 在线看片中文字幕| 亚洲自偷自拍另类小说| 刘亦菲一区二区在线观看| 亚洲欧美日本国产综合在线| 国产综合亚洲欧洲区精品无码| 国产午夜福利亚洲第一| 欧美一区国产| 狠狠干欧美| 高h视频在线| 91www在线观看|