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

一種考慮微凸體接觸角分布的結(jié)合面?zhèn)冉佑|分形模型

2023-12-18 03:21:38王世軍劉鑫范凌松崔圣奇李鵬陽(yáng)
關(guān)鍵詞:有限元模型

王世軍, 劉鑫, 范凌松, 崔圣奇, 李鵬陽(yáng)

(西安理工大學(xué)機(jī)械與精密儀器工程學(xué)院, 710048, 西安)

機(jī)械結(jié)構(gòu)中兩個(gè)相互接觸的表面稱為結(jié)合面。宏觀上是光滑平整的零件表面,微觀上總是粗糙不平的。兩個(gè)零件表面的接觸,微觀上都是兩個(gè)粗糙表面的接觸,接觸性質(zhì)并不是兩個(gè)光滑表面的接觸性質(zhì)。結(jié)合面的形貌特征和材料特性決定了結(jié)合面的接觸特性。如何從接觸面的微觀形貌入手構(gòu)建準(zhǔn)確的接觸模型,對(duì)于研究結(jié)合面的接觸特性[1]乃至分析機(jī)械整機(jī)的靜態(tài)和動(dòng)態(tài)性能,都有重要意義。

在結(jié)合面的接觸問(wèn)題研究中,學(xué)者大多將結(jié)合面簡(jiǎn)化成粗糙表面和剛性平面的接觸,相當(dāng)于峰-峰接觸,此時(shí)兩微凸體的接觸角度為0°。雙粗糙表面微觀上微凸體之間除了正接觸,大多為肩-肩接觸,接觸角度范圍在0°~90°之間[2]。

針對(duì)側(cè)接觸問(wèn)題,Sepehri等[3]將微凸體定義為橢球體,研究了傾斜接觸力在法線和切線方向上的分量,建立了SF統(tǒng)計(jì)模型。Gorbatikh等[4]提出了一個(gè)接觸角度分布函數(shù),并用統(tǒng)計(jì)方法研究了微凸體之間的側(cè)向接觸性質(zhì),但是沒(méi)有給出這個(gè)角度分布函數(shù)的理論依據(jù)。Misra等[5]通過(guò)引入一個(gè)方位函數(shù)來(lái)模擬接觸角和各向異性表面,利用微力學(xué)模型來(lái)研究?jī)A斜接觸微凸體的滑動(dòng)行為,并將剪切方向與法向的接觸剛度耦合起來(lái),建立了考慮接觸方向分布和微滑移特性的側(cè)接觸統(tǒng)計(jì)模型。高志強(qiáng)等[6]引入文獻(xiàn) [4]的接觸角度函數(shù),建立了考慮微凸體斜接觸和相互作用的結(jié)合面接觸剛度及阻尼模型。Jamshidi等[7]根據(jù)微凸體接觸時(shí)的穿透量來(lái)界定接觸狀態(tài),同時(shí)對(duì)Misra的接觸角度函數(shù)進(jìn)行修正,由此建立了接觸界面法向載荷與剪切力耦合的HJ模型。文獻(xiàn) [4-7]在研究斜接觸問(wèn)題時(shí)雖然給出了接觸角的分布函數(shù),但對(duì)其分布規(guī)律的來(lái)源和依據(jù)沒(méi)有明確說(shuō)明。范凌松等[8-9]研究了粗糙表面上相鄰微凸體的水平距離,發(fā)現(xiàn)水平距離遵循高斯分布規(guī)律,建立了包含法向和切向接觸特性的FLS模型,研究了相關(guān)輪廓統(tǒng)計(jì)參數(shù)及相互作用對(duì)結(jié)合面接觸剛度的影響。

利用分形模型研究側(cè)接觸問(wèn)題時(shí),孫見(jiàn)君等[10]根據(jù)微凸體斜接觸時(shí)上下兩輪廓的空間坐標(biāo)與接觸角的聯(lián)系,建立了雙粗糙面的側(cè)接觸分形模型,分析了表面形貌參數(shù)、接觸壓力對(duì)孔隙度和實(shí)際接觸面積的影響,該模型采用的接觸角函數(shù)是由空間幾何關(guān)系獲得,與水平距離無(wú)關(guān)。Zhang等[11]基于MB模型和側(cè)接觸理論建立了考慮接觸角分布的ZK改進(jìn)分形模型,該模型假定接觸角的范圍在0°~45°之間,通過(guò)給定接觸面積與接觸角的函數(shù)關(guān)系來(lái)研究接觸角對(duì)接觸剛度的影響,該模型忽略了水平距離的存在,且假定的角度范圍與實(shí)際不符。以上模型都未對(duì)接觸角的理論來(lái)源加以研究。

本文基于分形理論研究粗糙表面上相鄰微凸體之間的高度差和水平距離的分布規(guī)律,采用分形方法描述相鄰微凸體之間的距離分布,結(jié)合相鄰微凸體之間高度差的分形規(guī)律,遵循微凸體峰高與水平距離的統(tǒng)一性,得到粗糙表面上微凸體側(cè)向接觸時(shí)接觸角的分布函數(shù),同時(shí)考慮彈塑性邊界條件的連續(xù)性及域擴(kuò)展因子對(duì)真實(shí)接觸面積的影響,以此為根據(jù)建立考慮微凸體水平距離和接觸角分布的側(cè)接觸分形模型。通過(guò)仿真分析探究了相關(guān)接觸參數(shù)對(duì)結(jié)合面接觸剛度的影響。

1 微凸體中心水平距離分布規(guī)律

針對(duì)45鋼磨削加工表面,利用DCM3D萊卡顯微鏡分別沿著輪廓水平方向按1 μm分辨率和高度方向上按0.25 μm分辨率進(jìn)行采樣,根據(jù)采樣數(shù)據(jù)繪制的三維表面輪廓圖如圖1所示。

圖1 磨削表面三維輪廓

本文統(tǒng)一按“三點(diǎn)峰”的定義來(lái)獲取一組輪廓上n個(gè)峰值點(diǎn)的水平坐標(biāo)xi,i取1~n,以第1個(gè)峰為基準(zhǔn),之后每個(gè)峰相對(duì)上一個(gè)峰的坐標(biāo)差定義為水平距離si,同理把后續(xù)每個(gè)峰與前一個(gè)峰的高度差的絕對(duì)值定義為峰的高度差hi,如圖2所示。

圖2 微凸體水平和峰高距離定義

針對(duì)磨削試樣按照順紋理方向提取多組斷面的二維輪廓數(shù)據(jù),分別統(tǒng)計(jì)各組輪廓不同間隔下的水平距離數(shù)據(jù),擬合出的輪廓曲線如圖3所示。

圖3 不同間隔水平距離輪廓

已有研究表明,磨削加工得到的機(jī)械表面在微觀結(jié)構(gòu)上微凸體的接觸類型主要為側(cè)向接觸[10]。本文采用PSD法來(lái)獲取上述輪廓的分形參數(shù),首先利用MATLAB編寫程序,得到對(duì)數(shù)坐標(biāo)下的功率譜密度函數(shù),利用最小二乘法線性擬合求出其分形參數(shù),從而驗(yàn)證水平距離是否具有分形規(guī)律。

對(duì)文獻(xiàn) [12]中的PSD函數(shù)左右各取對(duì)數(shù)

lnP(ω)=(2D-5)lnω+(2D-2)lnG-ln(2lnγ)

(1)

式中:D為分形維數(shù);G為分形尺度參數(shù);ω為空間截止頻率;γ為尺度參數(shù),通常取1.5。PSD曲線斜率為2D-5,截距為(2D-2)lnG-ln(2lnγ)。按照該方法對(duì)圖3中的水平距離輪廓數(shù)據(jù)進(jìn)行計(jì)算。間隔1 μm下水平距離輪廓的PSD曲線和最小二乘法擬合的曲線如圖4所示。

圖4 1 μm間隔水平距離功率譜密度曲線

按照相同方法分析相鄰微凸體峰的高度差輪廓,峰的高度差滿足分形規(guī)律。根據(jù)三點(diǎn)峰兩組表面輪廓數(shù)據(jù)計(jì)算得到的分形參數(shù)如表1所示,分析表1中數(shù)據(jù)可以看出,不同采樣間隔下三點(diǎn)峰水平輪廓和峰的高度差輪廓的斜率k在 [-3,-1]之間[13],證明本文按三點(diǎn)峰定義的兩種輪廓滿足分形規(guī)律,因此可以利用二者的分形規(guī)律建立法向和切向剛度模型。

表1 輪廓分形參數(shù)統(tǒng)計(jì)

2 單對(duì)微凸體側(cè)接觸模型

微凸體在法向力P和切向力T作用下發(fā)生側(cè)接觸時(shí),其接觸過(guò)程如圖5所示,S、R分別為微凸體的高度、曲率半徑,下標(biāo)1、2代表上、下兩微凸體,下標(biāo)n、τ代表法線和切線方向,結(jié)合面的平均間距用d表示,微凸體中心水平距離和接觸角度分別用r、θ表示。在Z′OX′坐標(biāo)系中,法向接觸力P沿上下微凸體接觸面的法線和切線方向可分解為Pn、Pτ,同理切向力T可分解為Tn、Tτ,上下微凸體在Z′軸方向的形變量為ωZ′,X′軸方向位移量為ζX′。在ZOX坐標(biāo)系中,Z軸方向變形為ωZ,X軸方向位移為ζ。

圖5 微凸體肩對(duì)肩接觸

單對(duì)微凸體在發(fā)生側(cè)接觸時(shí)存在以下幾何關(guān)系

(2)

(3)

δ=S1+S2-d-r2/2Rs

(4)

式中Rs=R1+R2[14]。

根據(jù)側(cè)接觸受力關(guān)系可得

P=PZ′cosθ-TX′sinθ

(5)

T=PZ′sinθ+TX′cosθ

(6)

式中:PZ′=TX′tan(α+θ),其中PZ′、TX′分別代表Z′、X′方向所受的合力。

2.1 彈性接觸變形階段

單對(duì)微凸體發(fā)生側(cè)接觸時(shí),實(shí)質(zhì)上為兩個(gè)半球體在側(cè)臂上的接觸,此時(shí)兩球體存在一個(gè)接觸角θ,定義單對(duì)微凸體接觸時(shí)Z′方向上的變形量為ωZ′,由Weierstrass-Mandelbrot分形函數(shù)可得單個(gè)微凸體的變形及等效曲率半徑[15]

ωZ′=2GD-1(lnγ)0.5(2r′)2-D

(7)

R=r′D/[24-DGD-1(lnγ)0.5]

(8)

式中:r′為微凸體的截?cái)喟霃健?/p>

根據(jù)圖中幾何關(guān)系可以得到微凸體在Z′方向上的變形

ωZ′=ωZcosθ+ζsinθ

(9)

當(dāng)微凸體達(dá)到臨界變形接觸狀態(tài)時(shí),在Z′方向上變形量為[16]

δZ′ec=(πKH/2E)2R

(10)

聯(lián)立式(7)、(8)和(10)可得微凸體在Z′方向上的臨界接觸面積

(11)

微接觸點(diǎn)的截面積a的面積分布函數(shù)[18]為

(12)

式中:al為最大微接觸點(diǎn)的截?cái)嗝娣e;ψ為與分形維數(shù)D有關(guān)的域擴(kuò)展因子,且近似滿足ψ=1.549D-1.253+1.069。

當(dāng)微凸體在Z′方向上變形量ωZ′<ωZ′ec時(shí),此時(shí)微凸體為完全彈性接觸狀態(tài),依據(jù)赫茲接觸理論及文獻(xiàn) [19],可得微凸體彈性接觸階段的法向合力,切向合力與變形量ωZ′關(guān)系[19]如下

pe=4ER1/2b3/2c1/3

(13)

te=4ER1/2b3/2c2/3

(14)

b=ωZcosθ+ζsinθ

(15)

c1=[cosθ-sinθ/tan(θ+α)]

(16)

c2=[sinθ+cosθ/tan(θ+α)]

(17)

彈性階段,兩微凸體之間的法向接觸剛度

kne=dpe/dω=2ER0.5b0.5c1cosθ

(18)

切向接觸剛度為

kτe=dte/dζ=2ER0.5b0.5c2sinθ

(19)

2.2 彈塑性變形階段

上下微凸體發(fā)生側(cè)接觸時(shí)在彈塑性階段的臨界變形范圍為ωZ′ec≤ωZ′≤110ωZ′ec。根據(jù)文獻(xiàn) [20],在ZOX坐標(biāo)系中,微凸體在彈塑性區(qū)法向合力、切向合力與變形量ωZ′之間的關(guān)系分別為

(20)

(21)

在ZOX坐標(biāo)系中,微凸體的法向接觸剛度、切向接觸剛度可表示為

(22)

(23)

式中:m=ln(330/K)/ln(110)。

2.3 塑性變形階段

當(dāng)微凸體在Z′方向上變形量ωZ′≥110ωZ′ec時(shí),微凸體進(jìn)入完全塑性接觸狀態(tài),此時(shí)微凸體的平均接觸載荷與兩接觸材料中較軟一方的硬度相等。完全塑性變形階段Z方向上的法向合力與變形量ωZ′的關(guān)系[21]如下

pp=2πHRbc1

(24)

此時(shí)法向接觸剛度

(25)

3 結(jié)合面的接觸剛度分形模型

由于水平距離和峰的高度差都具有分形規(guī)律,因此可用W-M函數(shù)來(lái)表征相鄰微凸體的輪廓

(26)

(27)

式中:Dh,s、Gh,s分別為分形維數(shù)、分形尺度參數(shù);γn為表征粗糙表面的頻譜;n為頻率指數(shù);xs為水平距離;xh為峰-峰相對(duì)距離。

在坐標(biāo)系中,接觸角θ與兩個(gè)輪廓的正切值有關(guān),因此接觸角的分布函數(shù)為

θ(xs,xh)=arctan[H(xh)/S(xs)]

(28)

根據(jù)接觸角幾何關(guān)系及式(28)可得微凸體接觸角的概率密度函數(shù)

f(θ)=

(29)

式中:θ為肩對(duì)肩接觸角,范圍為0°~90°;L為接觸表面加載次數(shù)。

根據(jù)式(29)擬合出分布函數(shù)曲線,分布規(guī)律如圖6所示,可以發(fā)現(xiàn):接觸角度函數(shù)值隨著接觸角θ的增加而減小;隨著加載次數(shù)增加,接觸面逐漸壓平,接觸角度函數(shù)值變化逐漸緩慢。

圖6 側(cè)接觸角度分布規(guī)律

微凸體發(fā)生彈性、彈塑性和全塑性變形的3個(gè)區(qū)域面積之和構(gòu)成了結(jié)合面的真實(shí)接觸面積[20]

(30)

將式(12)代入式(30),可得

amax=Ar(2-D)Ar/[Dψ(1-0.5D)]

(31)

根據(jù)文獻(xiàn) [21]及式(7)、(8)、(10)、(12)和(31),結(jié)合面無(wú)量綱法向接觸載荷P*可表示為

((2-D)/D)0.5Df(θ)dadθ

(32)

無(wú)量綱切向接觸載荷T*為

a[(0.5D-1)+(1-D)m]c2f(θ)dadθ

(33)

(34)

(35)

4 結(jié)果與分析

(a)D對(duì)的影響

對(duì)上述仿真曲線進(jìn)行分析可得如下結(jié)果。

圖8 材料參數(shù)對(duì)切向接觸剛度的影響

(3)接觸參數(shù)對(duì)真實(shí)接觸面積的影響。真實(shí)接觸面積隨法向載荷的變化曲線如圖9所示,可知D、塑性指數(shù)φ的增大會(huì)使真實(shí)接觸面積不斷變大,D越大,表面輪廓越趨近于光滑,而φ越大微凸體越不易發(fā)生塑性變形,從而使臨界的接觸面積不斷變小,而真實(shí)接觸面積的占比不斷變大。因此,改善分形維數(shù)和塑性指數(shù)有助于提高真實(shí)接觸面積。

圖9 接觸參數(shù)對(duì)真實(shí)接觸面積的影響

(4)法向載荷系數(shù)的影響。圖10為結(jié)合面接觸剛度與法向載荷系數(shù)的變化關(guān)系。h分別取2~50時(shí),無(wú)量綱法向和切向接觸剛度均隨著法向載荷系數(shù)h的增大而增大,隨著分形尺度參數(shù)G*的增大而減小。固定分形尺度參數(shù)G*時(shí),保持結(jié)合面法向載荷P不變,h的增大意味著結(jié)合面切向接觸載荷T減小,則接觸面抵抗切向變形的能力增強(qiáng),從而結(jié)合面法向和切向接觸剛度都會(huì)增大。

與h的變化曲線

(5)接觸方式的影響。圖11為兩種模型的接觸剛度和真實(shí)面積隨法向載荷變化的曲線。其中LZT模型[16]的接觸類型為正接觸,本文模型則考慮了側(cè)接觸,給定相同的無(wú)量綱法向載荷,新模型的接觸剛度和真實(shí)接觸面積小于LZT模型,根據(jù)圖6及文獻(xiàn) [6]的研究可知,正接觸的接觸角度函數(shù)值相比側(cè)接觸較大,而LZT模型忽略了側(cè)接觸的影響,導(dǎo)致該模型的剛度偏高。

(a)真實(shí)接觸面積

5 有限元仿真與試驗(yàn)驗(yàn)證

采用文獻(xiàn) [9]的試驗(yàn)裝置來(lái)驗(yàn)證本文提出的側(cè)接觸分形模型,結(jié)合部由兩半環(huán)及螺栓構(gòu)成,半環(huán)材料均為45鋼,半環(huán)接觸面為磨削表面,兩材料的厚度和半徑分別為35、275 mm,表面經(jīng)丙酮清洗,其他材料參數(shù)如表2所示。

表2 整環(huán)材料參數(shù)

截取部分二維輪廓數(shù)據(jù),經(jīng)過(guò)PSD功率譜計(jì)算可得D=1.425,G=2.928×10-10m,半環(huán)聯(lián)結(jié)面上每個(gè)螺栓能承受的最大擰緊力矩為80 N·m。本文采用錘擊法獲得整環(huán)試驗(yàn)裝置的振型和固有頻率,同時(shí)建立整機(jī)的有限元接觸模型并進(jìn)行模態(tài)分析,綜合比較試驗(yàn)與有限元模型的結(jié)果。

由于包含結(jié)合部的機(jī)械結(jié)構(gòu)的接觸剛度具有非線性的特征,因此需要將包含接觸的非線性模型線性化,再通過(guò)有限元模態(tài)分析進(jìn)行線性化求解,具體過(guò)程如下[21]:基于虛擬材料法建立整環(huán)有限元接觸模型;通過(guò)靜力分析獲得接觸層在工作靜載荷下的剛度;如果整環(huán)工作時(shí)的載荷波動(dòng)遠(yuǎn)小于靜載荷,可認(rèn)為整環(huán)在工作時(shí)接觸層的剛度保持不變,即用靜力分析獲得的接觸層剛度建立整環(huán)的線性模型,將非線性模型線性化;通過(guò)整環(huán)線性模型的模態(tài)分析獲得結(jié)合部的振型及固有頻率。

依據(jù)式(34)、(35)及LZT模型計(jì)算的法向、切向接觸剛度數(shù)據(jù)可得虛擬材料層的非線性接觸參數(shù),試樣和整環(huán)有限元模型的材料參數(shù)一致。采用APDL語(yǔ)言將虛擬接觸層的材料參數(shù)[22]、非線性接觸參數(shù)寫入整環(huán)有限元接觸模型,整環(huán)模型共使用80 798個(gè)SOLID185單元,其中結(jié)合部虛擬接觸層單元有1 858個(gè),厚度為1 mm[23],如圖12所示。

圖12 整環(huán)有限元模型

對(duì)上述有限元接觸模型進(jìn)行靜力分析可以求得虛擬接觸層在靜載荷下的接觸剛度,再利用APDL語(yǔ)言將虛擬接觸層剛度代入模態(tài)分析,可得各階固有頻率及振型[22-27]。本文計(jì)算了法向面壓為0.6 MPa時(shí)本文模型與LZT模型[16]的有限元模態(tài),圖13為本文的試驗(yàn)與有限元前3階振型匯總,表3為3種模型不同面壓下的模態(tài)分析結(jié)果。

(a)第一階

由表3可知LZT模型的最大相對(duì)誤差絕對(duì)值為8.12%,而本文模型為4.32%,與LZT模型相比精度提高了46.80%。利用本文研究的接觸角分布函數(shù)及側(cè)接觸建模方法來(lái)計(jì)算有限元模型的固有頻率和振型,與試驗(yàn)結(jié)果相比誤差最小,本文模型的精度得到了驗(yàn)證。

6 結(jié) 論

(1)磨削表面微凸體大多為肩并肩的側(cè)接觸,本文對(duì)三點(diǎn)峰水平距離和峰的高度差進(jìn)行定義,經(jīng)過(guò)分析發(fā)現(xiàn),二者近似滿足分形規(guī)律,首次建立了考慮微凸體水平距離分布的結(jié)合面接觸剛度分形模型。

(2)本文構(gòu)造了微凸體水平距離和峰的高度差的分布函數(shù),根據(jù)側(cè)接觸相對(duì)關(guān)系確定了接觸角θ的分布函數(shù),與Gorbatikh模型的角度分布函數(shù)趨勢(shì)相同,接觸角θ和加載次數(shù)是影響接觸角度分布規(guī)律的主要原因。

(3)本文建立的有限元模型與試驗(yàn)結(jié)果最大相對(duì)誤差絕對(duì)值為4.32%,說(shuō)明考慮側(cè)接觸水平距離及接觸角影響的接觸剛度分形模型更加貼近磨削結(jié)合面的真實(shí)接觸特性,可為研究側(cè)接觸問(wèn)題提供新的思路。

猜你喜歡
有限元模型
一半模型
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
新型有機(jī)玻璃在站臺(tái)門的應(yīng)用及有限元分析
基于有限元的深孔鏜削仿真及分析
基于有限元模型對(duì)踝模擬扭傷機(jī)制的探討
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉(zhuǎn)換方法初步研究
磨削淬硬殘余應(yīng)力的有限元分析
基于SolidWorks的吸嘴支撐臂有限元分析
主站蜘蛛池模板: 色综合久久无码网| 国产精品对白刺激| 波多野结衣在线se| 在线色综合| 美女裸体18禁网站| 免费Aⅴ片在线观看蜜芽Tⅴ | 五月婷婷丁香色| 日本一本在线视频| 免费中文字幕在在线不卡| 免费一级α片在线观看| 中文字幕日韩久久综合影院| 国产美女主播一级成人毛片| 特级欧美视频aaaaaa| 欧美日韩一区二区在线免费观看| 国产亚洲精品无码专| 国产精品综合久久久| 福利姬国产精品一区在线| 中国毛片网| 强奷白丝美女在线观看| 中美日韩在线网免费毛片视频 | 综合人妻久久一区二区精品 | 四虎影院国产| 国产精品亚洲日韩AⅤ在线观看| 久青草免费在线视频| 99在线视频精品| 在线观看免费国产| 欧洲精品视频在线观看| 91外围女在线观看| 国产成人综合网| 亚洲A∨无码精品午夜在线观看| 97超碰精品成人国产| 亚洲精品麻豆| 亚洲区视频在线观看| 呦视频在线一区二区三区| 国产成人h在线观看网站站| AV无码无在线观看免费| 伊人狠狠丁香婷婷综合色| 国产精品真实对白精彩久久| 久久96热在精品国产高清| 伊人AV天堂| 91精品国产91久久久久久三级| 日韩毛片免费| 18禁影院亚洲专区| 久久亚洲高清国产| 久久伊人操| www.精品视频| 国产高清又黄又嫩的免费视频网站| 国产精品久久自在自线观看| 丝袜国产一区| 欧美日韩一区二区三区在线视频| 无码中文字幕精品推荐| 无码人中文字幕| 波多野结衣中文字幕一区二区| 欧美专区在线观看| 国产69囗曝护士吞精在线视频 | 97久久精品人人| 国产亚洲美日韩AV中文字幕无码成人 | 毛片在线看网站| 一级毛片在线播放免费| 国产美女在线观看| 四虎影视8848永久精品| 国产乱子伦视频在线播放| 福利一区在线| 亚洲综合国产一区二区三区| 国产成人在线无码免费视频| 男女性午夜福利网站| 亚洲精品制服丝袜二区| 欧美一区中文字幕| 99中文字幕亚洲一区二区| 一级成人a做片免费| 欧美人在线一区二区三区| 亚洲无码A视频在线| 亚洲精品国产首次亮相| 亚洲成人黄色在线观看| 国产精品部在线观看| 免费一极毛片| 亚洲成人网在线观看| 国产精品三级av及在线观看| 天天操天天噜| 91无码网站| 91精品国产综合久久不国产大片 | 亚洲精品中文字幕无乱码|