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

正交異性薄膜非線性振動分析

2018-06-28 13:29:50何澤青張冬輝栗穎思
振動與沖擊 2018年12期
關鍵詞:方向振動

何澤青, 張冬輝, 宋 林, 栗穎思, 王 生

(1. 中國科學院大學 材料科學與光電技術學院, 北京 100190; 2. 中國科學院 光電研究院, 北京 100094)

隨著浮空器研究的迅速發展和新興建筑領域的興起,現代復合膜材料由于其高強、輕質、柔軟、化學穩定等特性得到了廣泛的應用[1]。膜結構由于質量輕、張力大,普遍存在大變形、低頻率的振動特點,當薄膜結構的固有頻率與外界激勵頻率非常接近時,結構極容易發生共振甚至導致破壞,因此膜結構的動力學問題引起了人們的廣泛關注,對于膜結構自振特性的研究也就顯得尤為必要。

當前對于膜結構振動研究主要分為線性領域和非線性領域。在線性領域,文獻[2]對于各向同性張拉薄膜的自由振動進行了理論分析,并對具有規則外形的矩形和圓形薄膜的振動頻率和振型進行了理論求解;文獻[3]根據哈密頓原理建立了薄膜法向振動的動力學方程,求得周邊固定圓形薄膜、扇形薄膜自由振動的理論解;文獻[4]求解了各向同性的雙向受力不等的矩形、圓形、橢圓形平面薄膜的自振頻率與振型,并給出了任意外形邊界的平面薄膜的近似解,但分析中并未考慮薄膜正交異性情況,同時動力學方程中也未體現出非線性對結構振動的影響;文獻[5]構造了6節點三角形單元,根據哈密頓原理建立薄膜自由振動方程,并推導其單元剛度矩陣和單元質量矩陣。文獻[6]用Bessel函數作為薄膜橫向振動偏微分方程的解,建立復雜邊界固有頻率求解模型,并對圓形、L-形、分段圓弧邊界平面張 不拉薄膜固有頻率進行求解,分析了邊界幾何參數對薄膜固有頻率的影響。文獻[7]建立了薄膜二維振動的偏微分方程,并對無限大薄膜的二維受迫振動一般方程進行了求解。以上模態分析均是建立在初始幾何形態上,沒有考慮結構法向位移對其幾何形態的影響,這在小變形假定下是可以接受的。但是,膜結構本身由于具有較強的幾何非線性特點,其大變形效應對曲面本身有一定的修正作用,同時其形態的不穩定必然會對其剛度矩陣產生較大的影響。

在非線性領域,目前主要有能量法、范式理論方法、同倫攝動法和迭代攝動法等對非線性振動展開近似求解。文獻[8]基于能量原理求出非線性系統的一次近似解析解,然后利用牛頓迭代思想得到周期系數微分方程,最后根據諧波平衡原理及最小二乘法求其高次近似解。文獻[9]拓展了范式方法,改進了Nayfeh關于響應頻率的提取,使其可以研究強非線性振動問題。文獻[10]利用微分求積法計算了變密度薄膜的橫向振動,獲得了系統無量綱復頻率與密度系數、薄膜張力比的關系曲線。文獻[11]應用了一種推廣的平均法研究強非線性振子的共振響應,并與數值積分法進行了對比。文獻[12]通過兩個典型張拉膜結構工程實例采用Lanczos法研究了膜面初張力、膜面矢高和膜材質量對自振特性的影響,分析結果表明膜面張力與結構質量是影響膜結構自振特性的主要因素,大變形產生的位移修正效應對膜結構自振特性也有較大的影響。文獻僅是通過仿真分析對上述現象進行闡述,并未從理論上對其進行分析和說明。

本文將參考Qian的研究對于各向同性的雙向受力不等平面薄膜的振動分析,對正交異性薄膜振動進行數理分析,建立薄膜振動的平衡微分方程,同時更進一步的根據薄膜振動的大撓度理論[13-14]引入其動力學方程的非線性項,考慮薄膜的正交異性和振動中的非線性特點,對動力學方程進行求解,得出非線性振動頻率的解析解,目的是明確非線性因素對于振動頻率的影響。同時為了提高理論研究在工程應用中的便利性,通過對方程解析解進行合理簡化,得到一種簡單精確,適用于工程需要的非線性頻率計算公式。

1 動力學方程

1.1 動力學微分方程

建立如圖1所示的坐標系,在振動的膜上選取一振動的微單元,l1和l2為微元表面的弧線,l1與XOZ平面平行,l2與YOZ平面平行。

圖1 振動的膜微元Fig.1 Vibrating micro-membrane unit

在微元dy邊界,單位長度的張力為Nx,沿微元切線方向,與X坐標軸夾角為α,因此作用在x端的單位長度的張力在Z軸方向的分量為Nxsinα,由于α較小,則有sinα≈tanα,用w表示膜面上一點在振動中的撓度,因此有

于是作用在整個dy邊上的Z軸方向的力為

而在x+dx端的Z軸方向的力應為

由此可得在該面元的x與x+dx邊上Z軸方向的合力為

同理,作用在dx邊界上的的張力在Z軸方向的合力可以表示為

所以作用在整個微元上Z軸方向的合力為

設膜的面密度為ρ,則微元的質量為ρdxdy,根據達朗貝爾原理可寫出微元的運動微分方程

整理后可得

(1)

式中:Nx為微元dy邊界單位長度的張力,沿微元切線方向;Ny為微元dx邊界單位長度的張力,沿微元切線方向;w為微元撓度;ρ為微元的面密度。

對于受預張力薄膜在法向載荷下的振動微分方程,在式(1)中加入預張力和法向載荷項,就可以得到在法向載荷作用下,正交異性雙向受力不等的預張力薄膜振動微分方程

(2)

式中:Nox為微元dy邊界單位長度的預張力,沿微元切線方向;Noy為微元dx邊界單位長度的預張力,沿微元切線方向;p(x,y,t)為微元法向載荷。

1.2 變形協調方程

根據彈性大撓度理論[15],薄膜變形后,其應變由線性和非線性兩部分組成。其中,線性應變由面內位移u,v引起,非線性應變由撓度w引起,因此總的應變為

(3)

式中:εx為微元沿l1方向的應變;εy為微元沿l2方向的應變;γxy為切應變;u為微元沿l1方向的面內位移;v為微元沿l2方向的面內位移。

將式(3)中的位移u,v消去后可得膜面應變與撓度必須滿足的變形連續條件,即相容方程

(4)

1.3 物理方程

根據膜材料的均勻性與正交異性假定,纖維方向為彈性主方向,令其與l1,l2方向一致。按正交異性理論,l1和l2兩個方向的彈性模量分別為E1和E2,剪切模量為G12,徑向泊松比為μ1,緯向泊松比為μ2,彈性模量與泊松比滿足下面的關系式

(5)

根據應力應變物理方程關系式有

(6)

即有

(7)

式中:σx為微元dy邊界的正應力,沿微元切線方向;σy為微元dx邊界的正應力,沿微元切線方向;τxy為微元切應力;h為膜材厚度。

將式(7)代入相容方程式(4)中,得到{N}用w表示的相容方程

(8)

1.4 控制方程組的建立

總結以上內容,可得正交異性薄膜非線性振動問題的控制方程組為

(9)

式中:Nox為微元dy邊界單位長度的預張力,沿微元切線方向;Noy為微元dx邊界單位長度的預張力,沿微元切線方向;Nx為微元dy邊界單位長度的張力,沿微元切線方向;Ny為微元dx邊界單位長度的張力,沿微元切線方向;Nxy為微元的面內剪切力;ρ為微元的面密度;h為微元厚度;E1為微元沿l1方向的彈性模量;E2為微元沿l2方向的彈性模量;G12為剪切模量;μ1為微元沿l1方向的泊松比;μ2為微元沿l2方向的泊松比;w=w(x,y,t)為微元撓度;p(x,y,t)為微元法向載荷。

在薄膜振動過程中,根據膜面無矩無剪力假定[16],同時引入應力函數φ(x,y,t)[17]

(10)

對于自由振動薄膜,取其法向力p(x,y,t)=0,將式(10)代入式(9),則控制方程組為

(11)

2 方程簡化

2.1 初始條件與邊界條件

以四周固支的矩形薄膜為例,如圖2所示。

圖2 四周固支的矩形薄膜Fig.2 Rectangular membrane with four edges clamped

在矩形薄膜邊界處,其撓度在振動過程中始終為零,并且運動的速度也為零,即薄膜應滿足如下邊界條件:

x=0或a時

(12)

y=0或b時

(13)

t=0時

(14)

以上三組邊界和初始條件即為薄膜在振動過程中的定解條件,在確定薄膜的模態參數時將會用到。

2.2 控制方程組簡化

根據函數的性質,對w(x,y,t)采用分離變量法進行求解,設滿足邊界條件的位移函數為[18]

式中:w(x,y,t)為膜面撓度;m為x方向的節線數,取正整數;n為y方向的節線數,取正整數;a,b為矩形膜面的邊長;X(t)為膜面撓度最大值的時間函數。

將上式只取一項進行求解,再將最終結果按上式進行求和計算,并不影響最終結果,則可以得到

(15)

將式(15)代入式(11)的第二式可以得到

(16)

(17)

將式(17)代入式(16)可得

將式(17)代入邊界條件式(12)和式(13)可得

將上面所得的系數ki代入式(17)可得

(18)

將式(18)代入式(10)的前兩式可得

(19)

(20)

將式(15)、式(19)、式(20)代入式(11)的第一式可得

(21)

(22)

3 方程求解

根據式(22)的推導得到薄膜振動頻率方程的非線性表達式,對于式(22)可以采用利茲-伽遼金法進行求解[19],取權函數為W,構造積分方程為

(23)

其中,

式(23)為一個典型的非線性振動Duffing方程[20],由于X3(t)的系數B并非為一個微小量,因此式(23)為一個強非線性振動方程。

3.1 方程的解析解

根據方程(23)可得

(24)

由式(15)可知,當t=0時有

令X(0)=A,表示初始時刻(t=0)膜面撓度的最大值,有

根據邊界條件式(14),當t=0時

將上式再次代入式(24),有

(25)

對上式進行變量分離,得

薄膜從平衡位置到最大位移A處需要經歷1/4周期,因此對上式展開積分

則可以得到薄膜振動的周期T為

(26)

將上式代入式(26)并進行積分可得

將p和q代入上式,則可以得出薄膜振動的頻率f為

(27)

式(27)即為正交異性薄膜非線性振動頻率的解析解的表達式。由上式可見,薄膜振動頻率的非線性解不僅包含其頻率的線性部分,同時還與薄膜經緯方向的彈性模量有關,并且受薄膜振動初始位移的影響。

3.2 方程的近似解

根據邊界條件式(14),可以設式(23)的解為

X(t)=Acos(ωt)=Acosφ=X

取權函數為cosφ,代入式(23),并經由伽遼金法在一個周期T內積分得

(28)

將X=Acosφ代入式(28),可得

(29)

由式(29)可得

得出

(30)

式(30)即為正交異性薄膜非線性振動頻率的近似解的表達式。

3.2 非線性項對頻率的影響分析

對于薄膜振動微分方程式(23),若令其非線性項系數B=0,則可以得到標準的線性振動微分方程

(31)

對式(31)進行求解則可以得到薄膜振動頻率的線性頻率f0的表達式為

(32)

由式(30)和式(32)可得

(33)

從式(23)可知,B,K>0,因此ε>0,則式(33)可轉化為

(34)

從式(34)可以看出,薄膜的非線性自由振動仍然為周期振動,但振動的頻率值隨薄膜初始位移A的變化而變化,而不同于線性系統的固有頻率,如圖3所示。

圖3 非線性自振頻率-初始位移曲線Fig.3 Nonlinear frequency according to different initial displacement

4 算例分析

以圖2所示的四周固支矩形薄膜為例,分別對其自振頻率的線性解、非線性解析解和非線性近似解進行計算,并展開比較分析。薄膜參數為:a=2 m,為薄膜長邊長度;b=1 m,為薄膜短邊長度;ρ=0.12 kg/m2,為薄膜面密度;Nox=5 000 N/m,為膜面x向預張力;Noy=2 500 N/m,為膜面y向預張力;h=0.16 mm,為薄膜厚度;E1=0.8 GPa,為薄膜x向彈性模量;E2=0.6 GPa,為薄膜y向彈性模量。

4.1 自振頻率的線性解

將K代入式(32),則可以得到薄膜自振頻率的線性解公式

(35)

從式(35)中可以看出,由于并未考慮正交方向不同的彈性模量的影響以及撓度引起的薄膜面內應力的變化所帶來的影響,而是按均勻各向同性材料來考慮,因此其線性頻率值僅與初始條件以及邊界條件有關。下面將對其前八階線性頻率進行計算,如表1所示。

表1 自振頻率的線性解

4.2 自振頻率的非線性理論解

根據式(27)計算振動頻率的解析解。首先以一階(m=1,n=1)頻率為例,確定其收斂性。

編制程序,取初始位移A=0.02 m,計算f11與j的關系如表2所示。

表2 解析解的收斂性

從表2可以看出,當j≥2時,計算結果即收斂于穩定值,可以判定頻率的解析解是快速收斂并趨于穩定的。

根據式(27),取j=100,則正交異性薄膜非線性振動前八階頻率的解析解,如表3所示。

表3 自振頻率的非線性解析解

自振頻率的解析解隨初始位移變化情況如圖4所示。

從表3可以看出,當初始位移為0時,薄膜的自振頻率與其按小撓度計算得到的線性解相同;而從圖4可以看出,隨著初始位移的增加,薄膜的自振頻率也逐漸變大,這是由于隨著初始位移的增加,薄膜面內應力增大,導致薄膜剛度增大,從而引起其振動頻率變大。

圖4 前八階頻率-初始位移曲線Fig.4 The first eight modal frequencies according to different initial displacement

4.3 自振頻率的非線性近似解

式(30)為經由伽遼金算法得出的正交異性薄膜非線性振動頻率的近似解,其頻率公式為

其前八階計算頻率如表4所示。

表4 自振頻率的非線性近似解

薄膜自由振動前三階理論解與近似解的對比如圖5~圖7所示。

圖5 一階頻率-初始位移曲線Fig.5 The first order frequency according to different initial displacement圖6 二階頻率-位移曲線Fig.6 The second order frequency according to different initial displacement圖7 三階頻率-位移曲線Fig.7 The third order frequency according to different initial displacement

圖5~圖7顯示了薄膜一~三階頻率的理論解與近似解隨初始位移的變化情況。從圖中可以看出,近似解與理論解吻合情況良好,二者的誤差隨著初始位移的增加而逐漸增大,表4中的最大誤差僅為1.46%。更進一步的,當A→∞時,頻率的近似解與理論解的比值為即對于任意的初始位移,近似解與理論解的最大誤差不超過3.86%。

5 結 論

本文通過建立正交異性薄膜的非線性振動微分方程并對其進行求解,得到了其振動頻率的理論解,同時利用伽遼金法對方程進行了簡化求解,并得到了其近似解。利用此兩種解分別對四邊簡支矩形正交異性薄膜結構的線性頻率和非線性自由振動頻率分別進行了求解計算和對比,結果表明:

(1)薄膜結構線性振動分析僅適用于小撓度情況。

(2)隨著薄膜初位移的增加,導致薄膜面內剛度變大,進而引起其頻率變大。

(3)本文得到的薄膜自振頻率的近似解簡單實用,且具有較高的精度,對于任意的初始位移,按近似解表達式計算的結果與理論解的最大誤差不超過3.86%,完全能夠滿足工程精度的要求。

綜上所述,文中通過解析方法對正交異性薄膜的非線性振動展開分析,建立了結構的振動微分方程,方程在推導過程中并未對薄膜的結構形態進行特殊設定,因此其表達式是通用的,對于復雜形態的薄膜結構同樣適用,當系統的結構形式和邊界條件明確時,可以通過解析方法得出非線性系統的理論解或近似解。相較于數值方法,解析法由于不存在截斷誤差和舍入誤差,因此具有更高的精度,同時避免了數值方法在迭代計算過程中可能出現的不易收斂問題。為了提高分析結果的工程適用性,文中利用伽遼金法對非線性振動的Duffing方程進行求解并對結果展開誤差分析,得到的近似解不僅適用于線性和弱非線性系統,對強非線性系統同樣適用,因此具有較大的適用范圍和較高的精度,豐富和完善了非線性振動的研究,其研究結果為浮空器和建筑膜結構的工程設計提供了理論計算依據。

參 考 文 獻

[ 1 ] 陳務軍. 膜結構工程設計[M]. 北京: 中國建筑工業出版社, 2004: 134-139.

[ 2 ] TIMOSHENKO S, YOUNG S H, WEAVER W. Vibration problems in engineering[M]. 4th ed. New York: John Wiley & Sons, 1974: 327-334.

[ 3 ] 林文靜, 陳樹輝, 李森. 圓形薄膜自由振動的理論解[J]. 振動與沖擊, 2009, 28(5): 84-86.

LIN Wenjing, CHEN Shuhui, LI Sen. Analytical solution of the free vibration of circular membrane[J].Journal of Vibration and Shock, 2009, 28(5): 84-86.

[ 4 ] QIAN Guozhen. Solution for free vibration problem of membrane with unequal tension in two directions[J]. Applied Mathematics and Mechanics, 1982, 3(6): 885-892.

[ 5 ] 林文靜, 陳樹輝. 平面薄膜自由振動的有限元分析[J]. 動力學與控制學報, 2010, 8(3): 202-206.

LIN Wenjing, CHEN Shuhui. Free vibration analysis of plane membranes by finite element method[J]. Journal of Dynamics and Control, 2010, 8(3): 202-206.

[ 6 ] 劉充, 李玉宇, 保宏, 等. 邊界幾何參數對空間平面張拉膜結構固有頻率影響研究[J]. 振動與沖擊, 2015, 34(20): 198-202.

LIU Chong, LI Yuyu, BAO Hong, et al. Natural frequencies of pre-tensioned membrane structure with different boundary geometrical parameters [J]. Journal of Vibration and Shock, 2015, 34(20): 198-202.

[ 7 ] 張俊生. 薄膜二維振動數理方程的推導與求解[J]. 榆林學院學報, 2006, 16(6): 29-31.

ZHANG Junsheng. Inferential reasoning solution of thin film 2D vibration M&P equation [J]. Journal of Yulin University, 2006, 16(6): 29-31.

[ 8 ] 周一峰. 強非線性系統周期解的能量法[J]. 力學季刊, 2002, 23(4): 514-520.

ZHOU Yifeng. Energy iteration method for analytic periodic solutions of full strongly nonlinear vibration systems[J]. Chinese Quarterly Mechanics, 2002, 23(4): 514-520.

[ 9 ] 張琪昌, 郝淑英, 陳予恕. 用范式理論研究強非線性振動問題[J].振動工程學報, 2000, 13(3): 481-486.

ZHANG Qichang, HAO Shuying, CHEN Yushu. Study on strongly non-linear vibration systems by normal form theory [J]. Journal of Vibration Engineering, 2000,13(3): 481-486.

[10] 武吉梅, 陳媛, 王硯, 等. 基于微分求積法的印刷運動薄膜動力穩定性分析[J]. 振動與沖擊, 2015, 34(20): 57-60.

WU Jimei, CHEN Yuan, WANG Yan, et al. Dynamic stability of printing moving membrane based on differential quadrature method [J]. Journal of Vibration and Shock, 2015, 34(20): 57-60.

[11] 徐兆, 詹杰民. 強非線性振子的受迫振動[J]. 中山大學學報(自然科學版), 1995, 34(2): 1-6.

XU Zhao, ZHAN Jiemin. Forced oscillations of strongly nonlinear oscillators [J]. Journal of Sun Yat-sen University(Natural Science), 1995, 34(2): 1-6.

[12] 余志祥, 趙雷. 張拉膜結構自振特性研究[J]. 西南交通大學學報, 2004, 39(6): 734-739.

YU Zhixiang, ZHAO Lei. Research on free vibration properties of membrane structure [J]. Journal of Southwest Jiaotong University, 2004, 39(6): 734-739.

[13] CHEN Shanlin, ZHENG Zhoulian. Large deformation of circular membrane under the concentrated force [J]. Applied Mathematics and Mechanics, 2003, 24(1): 28-31.

[14] WANG Jing. Nonlinear free vibration of the circular plate with large deflection [J]. Journal of South China University of Technology, 2001, 29(8): 4-6.

[15] LIU Changjiang, ZHENG Zhoulian, YANG Xiaoyan, et al. Nonlinear damped vibration of pre-stressed orthotropic membranestructure under impact loading [J]. International Journal of Structural Stability and Dynamics, 2014, 14(1): 1091-1115.

[16] 喬磊, 譚峰, 楊慶山. 薄膜結構的動力反應分析[J]. 振動與沖擊, 2011, 30(6): 109-113.

QIAO Lei, TAN Feng, YANG Qingshan. Dynamic analysis of membrane structures [J]. Journal of Vibration and Shock, 2011, 30(6): 109-113.

[17] 徐芝綸. 彈性力學[M]. 4版. 北京: 高等教育出版社, 2006: 32-36.

[18] ZHENG Zhoulian, LIU Changjiang, HE Xiaoting, et al. Free vibration analysis of rectangular orthotropic membranes in large deflection [J]. Mathematical Problems in Engineering, 2009(4): 634362.

[19] 周紀卿, 朱因遠. 非線性振動[M]. 西安: 西安交通大學出版社, 2001: 140-144.

[20] 劉延柱, 陳立群. 非線性振動[M]. 北京: 高等教育出版社, 2001: 59-60.

猜你喜歡
方向振動
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
噴水推進高速艇尾部振動響應分析
2022年組稿方向
計算機應用(2022年2期)2022-03-01 12:33:42
2022年組稿方向
計算機應用(2022年1期)2022-02-26 06:57:42
2021年組稿方向
計算機應用(2021年4期)2021-04-20 14:06:36
2021年組稿方向
計算機應用(2021年3期)2021-03-18 13:44:48
2021年組稿方向
計算機應用(2021年1期)2021-01-21 03:22:38
This “Singing Highway”plays music
振動攪拌 震動創新
中國公路(2017年18期)2018-01-23 03:00:38
中立型Emden-Fowler微分方程的振動性
主站蜘蛛池模板: 五月婷婷综合网| 久久久久亚洲av成人网人人软件| 在线视频亚洲色图| 福利片91| 2048国产精品原创综合在线| 成人年鲁鲁在线观看视频| 欧美日本激情| 91欧美在线| 91福利免费| 狠狠色丁婷婷综合久久| 一本色道久久88综合日韩精品| 国产精品免费福利久久播放| 欧美综合成人| 伊人久久福利中文字幕| 久久99国产乱子伦精品免| 国产白浆在线| 伊人久久久久久久久久| 永久免费无码日韩视频| 欧美日本在线观看| 久久精品无码一区二区日韩免费| 中文字幕无码制服中字| 免费A级毛片无码免费视频| 亚洲自偷自拍另类小说| 免费日韩在线视频| 国产尤物视频在线| 日韩在线2020专区| 丁香婷婷激情网| 国产精品真实对白精彩久久| 国产国拍精品视频免费看 | 国产成人精品男人的天堂| 精品亚洲欧美中文字幕在线看| 亚洲不卡影院| 亚洲伦理一区二区| 亚洲乱码在线视频| 日韩成人午夜| 国产极品嫩模在线观看91| 色偷偷一区| 亚洲二区视频| 毛片网站免费在线观看| 亚洲中文字幕国产av| 伊人久久久久久久| 不卡无码网| 九九久久精品国产av片囯产区| 亚洲三级网站| 国产va在线观看免费| 国产精品护士| 久久 午夜福利 张柏芝| 自偷自拍三级全三级视频| 十八禁美女裸体网站| 国产在线小视频| 国产91小视频| 欧美人与牲动交a欧美精品 | 欧美h在线观看| 欧美另类图片视频无弹跳第一页| 少妇人妻无码首页| 国产精品一区在线观看你懂的| 亚洲资源站av无码网址| 国产在线观看一区精品| 在线免费无码视频| AV在线天堂进入| 欧美天堂在线| 99国产精品国产高清一区二区| 国产一区成人| 自拍偷拍欧美| 亚洲精品国产日韩无码AV永久免费网 | 精品国产欧美精品v| 日韩午夜伦| 亚洲无码91视频| 毛片网站观看| 色婷婷视频在线| 暴力调教一区二区三区| 欧美曰批视频免费播放免费| 欧美啪啪精品| 国产黄在线免费观看| 97亚洲色综久久精品| 亚洲一区免费看| a毛片免费在线观看| 熟妇人妻无乱码中文字幕真矢织江| 曰韩人妻一区二区三区| 香蕉视频在线观看www| 成人午夜天| 国产1区2区在线观看|