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

螺旋型電磁超聲輻射聲場(chǎng)對(duì)缺陷檢測(cè)影響的計(jì)算機(jī)模擬?

2023-09-15 12:36:04劉婉婷羅云躍杜好陽(yáng)朱大銘
應(yīng)用聲學(xué) 2023年4期
關(guān)鍵詞:磁場(chǎng)有限元

劉婉婷 羅云躍 杜好陽(yáng) 朱大銘 吳 剛 王 凱

(1 國(guó)網(wǎng)吉林省電力有限公司電力科學(xué)研究院 長(zhǎng)春 130021)

(2 國(guó)網(wǎng)吉林省電力有限公司 長(zhǎng)春 130022)

(3 中國(guó)科學(xué)院金屬研究所 沈陽(yáng) 110016)

0 引言

早在1970 年代,人們對(duì)電磁超聲(Electromagnetic acoustic testing,EMAT)的研究就有著明確的目標(biāo),即開(kāi)發(fā)非接觸缺陷檢測(cè)的工具[1]。雖然對(duì)EMAT 的理論研究推進(jìn)了EMAT 探頭設(shè)計(jì)和儀器的發(fā)展,但目前還不能認(rèn)為EMAT是一種通用且有效的缺陷檢測(cè)技術(shù)。一般認(rèn)為,這歸咎于EMAT 的換能效率比較低。但隨著現(xiàn)代電子技術(shù)的發(fā)展,具有高功率發(fā)射電路和高信噪比、高增益的接收電路的EMAT檢測(cè)儀,可以在一定程度上彌補(bǔ)EMAT技術(shù)的這一不足。因此,從輻射聲場(chǎng)角度深入探究EMAT 未能有效應(yīng)用于缺陷檢測(cè)的原因也是十分必要的。

一般來(lái)講,EMAT 是由永磁體和線圈組成,分別產(chǎn)生偏置靜磁場(chǎng)和動(dòng)磁場(chǎng)。不同的偏置磁場(chǎng)方向和線圈形式就能組合出不同的EMAT形式,也就對(duì)應(yīng)著不同的輻射聲場(chǎng)指向性。通過(guò)研究EMAT輻射聲場(chǎng)能夠得到橫波和縱波的指向性,指導(dǎo)應(yīng)用EMAT 進(jìn)行缺陷檢測(cè),對(duì)不同的檢測(cè)場(chǎng)合選用不同的檢測(cè)參數(shù)和探頭形式。

前人已經(jīng)開(kāi)展了一些關(guān)于EMAT 輻射聲場(chǎng)的研究工作。Pardee等[2]利用格林函數(shù)方法(Green’s function method,GFM)研究了曲折型EMAT的聲場(chǎng)輻射,并在數(shù)值上得到了曲折型EMAT 的輻射聲場(chǎng)指向性。而后,Kawashima等[3]就螺旋型EMAT的輻射聲場(chǎng)展開(kāi)了討論,他們借用了渦流檢測(cè)中著名的Dodd &Deeds 模型得到了在金屬工件表面產(chǎn)生的渦流大小,但是他們將EMAT 線圈當(dāng)成單匝環(huán)形線圈考慮,與實(shí)際的線圈不太符合。另外,Jiang等[4]也給出了螺旋型EMAT 的聲場(chǎng)分布,指出在軸線上并沒(méi)有橫波波幅產(chǎn)生,并給出了臨界角的計(jì)算公式。Zhai 等[5]利用截?cái)鄥^(qū)域特征函數(shù)展開(kāi)(Truncatedregion eigenfunction expansion,TREE)法來(lái)處理具有多匝的螺旋線圈在工件表面產(chǎn)生的渦流,簡(jiǎn)化數(shù)值計(jì)算過(guò)程。通過(guò)前人的研究,螺旋線圈的渦流場(chǎng)解已經(jīng)達(dá)到了很高的精度。但是沒(méi)給出螺旋型EMAT 輻射聲場(chǎng)指向性以及在缺陷檢測(cè)中比較關(guān)心的一些特征參數(shù),如擴(kuò)散角、偏離角等參數(shù)。

本文從EMAT輻射聲場(chǎng)角度出發(fā),利用有限元方法(Finite element method,FEM)與GFM,得到了螺旋型EMAT 在非鐵磁性材料中輻射聲場(chǎng)的基本特征,在FEM 和格林函數(shù)法聲場(chǎng)計(jì)算結(jié)果相一致的前提下,通過(guò)與壓電活塞聲場(chǎng)進(jìn)行實(shí)驗(yàn)對(duì)比,驗(yàn)證了螺旋型EMAT輻射聲場(chǎng)的特征參數(shù)。同時(shí)本文還分析探討了螺旋型EMAT 輻射聲場(chǎng)對(duì)非鐵磁性金屬材料中進(jìn)行缺陷檢測(cè)時(shí)的影響以及與壓電超聲的區(qū)別。

1 螺旋型EMAT輻射聲場(chǎng)的有限元分析

1.1 螺旋型EMAT模型

軸對(duì)稱結(jié)構(gòu)的螺旋型體波EMAT 的基本配置如圖1 所示。偏置靜磁場(chǎng)的永磁鐵為NdFeB(N52),其剩磁強(qiáng)度Br=1.4 T,相對(duì)磁導(dǎo)率μmag=1.05,高度hmag=20 mm,半徑rmag=15 mm。線圈采用銅導(dǎo)線繞制,導(dǎo)線半徑rcoil=0.1 mm,距試樣表面提離距離h=0.05 mm,線圈內(nèi)徑rin=2 mm,線圈外徑rout=11 mm,線圈匝數(shù)n=40,施加的電流I=50 A。試樣材料為鋁合金,相對(duì)磁導(dǎo)率μsample=1,電導(dǎo)率σsample=3.7×107S/m,橫波速度cS=3120 m/s,縱波速度cL=6120 m/s。為了得到輻射聲場(chǎng),采用有限元軟件進(jìn)行頻域分析,并且利用完美匹配層(Perfectly matched layer,PML)來(lái)模擬無(wú)限大半空間,在半空間界面施加應(yīng)力為0邊界條件。

圖1 螺旋型EMAT 結(jié)構(gòu)示意圖與線圈結(jié)構(gòu)示意圖Fig.1 Schematic diagram of spiral EMAT structure and coil structure

永磁鐵產(chǎn)生的沿r與z方向靜磁場(chǎng)表示為,在邊界上需要滿足靜磁場(chǎng)邊界條件。線圈通入脈沖交變電流后在工件中產(chǎn)生的動(dòng)磁場(chǎng)磁感應(yīng)強(qiáng)度沿著r與z方向表示為。在工件中產(chǎn)生的沿著旋轉(zhuǎn)方向的渦流密度為Je。由電磁學(xué)理論,由靜磁場(chǎng)和動(dòng)磁場(chǎng)的磁感應(yīng)強(qiáng)度產(chǎn)生的洛倫茲力大小分別表示為

其中,每個(gè)力分量都會(huì)在工件中產(chǎn)生縱波(L 波)和橫波(S 波),并且可以利用位移z分量和位移r分量表示,如圖2所示。如,由產(chǎn)生的橫波和縱波位移可以分別表示為,并且與位移z分量和r分量有如下轉(zhuǎn)換關(guān)系(見(jiàn)圖2):

圖2 位移r、z 分量與橫波(S)、縱波(L)位移轉(zhuǎn)換關(guān)系示意圖Fig.2 Schematic diagram of the transformation relationship between displacement r,z components and shear wave (S) and longitudinal wave(L) displacement

1.2 靜磁場(chǎng)與動(dòng)磁場(chǎng)的分析

利用有限元模型可以計(jì)算得到EMAT 模型中的靜磁場(chǎng)與動(dòng)磁場(chǎng)分布。圖3 是靜磁場(chǎng)的磁感應(yīng)強(qiáng)度沿r向的分布。動(dòng)磁場(chǎng)是施加交變電流時(shí)產(chǎn)生的,并在工件內(nèi)產(chǎn)生渦流。當(dāng)施加的激勵(lì)電流的頻率為1 MHz、幅度為50 A時(shí)得到的渦流密度幅度隨距離試樣表面深度的變化如圖4 所示,利用有限元計(jì)算得出的集膚深度與理論值δ一致。

圖3 永磁鐵產(chǎn)生的磁感應(yīng)強(qiáng)度z 分量與r 分量隨半徑r 的變化Fig.3 The z-component and r-component of the magnetic induction intensity generated by the permanent magnet as a function of the radius r

圖4 渦流大小隨距離試樣表面深度的變化Fig.4 Variation of eddy current size with depth from specimen surface

圖5 給出了根據(jù)式(1)得到的試樣表面的洛倫茲力。可以看到,相比于靜磁場(chǎng)產(chǎn)生的洛倫茲力,由動(dòng)磁場(chǎng)產(chǎn)生的洛倫茲力很小,可以忽略不計(jì)。另外,靜磁場(chǎng)產(chǎn)生的r方向洛倫茲力是最主要的,在產(chǎn)生輻射聲場(chǎng)也起到了主導(dǎo)作用。

圖5 靜磁場(chǎng)與動(dòng)磁場(chǎng)產(chǎn)生的洛倫茲力大小沿徑向分布Fig.5 The magnitude of the Lorentz force generated by the static and dynamic magnetic fields is distributed in the radial direction

1.3 螺旋型EMAT產(chǎn)生的輻射聲場(chǎng)

根據(jù)靜磁場(chǎng)與渦流作用產(chǎn)生的r和z方向的洛倫茲力,可以得到螺旋型EMAT 的輻射聲場(chǎng)。沿著圖6(a)所示的從左端起始的圓弧線獲取聲場(chǎng)橫波位移uS和縱波位移uL,得到圖6(b)的結(jié)果。可見(jiàn),螺旋型EMAT 的輻射聲場(chǎng)幅度以橫波為主,縱波相對(duì)較小。

圖6 試樣圓弧線與沿著圓弧線上橫波與縱波位移變化曲線Fig.6 The sample arc line and the shear wave and longitudinal wave displacement curve along the arc line

聲束指向性是換能器的一個(gè)重要特征。螺旋型EMAT 的橫波(uS)聲束指向性與縱波(uL)聲束指向性如圖7 所示。與壓電超聲輻射聲場(chǎng)相似,螺旋型EMAT 產(chǎn)生的橫波聲場(chǎng)也可分為具有超聲波干涉作用的近場(chǎng)區(qū)和隨著距離增大位移幅度單調(diào)變化的遠(yuǎn)場(chǎng)區(qū)。圖7(b)和圖7(d)分別給出了橫波和縱波的輻射聲場(chǎng)指向性。可以看到,橫波遠(yuǎn)場(chǎng)聲束輻射指向性圖上中心軸線上位移幅度為0,聲場(chǎng)呈喇叭型,最大幅度處的指向角為6.3?。而縱波聲束輻射指向性在中心軸線處位移幅度不為0,并且在稍偏離軸線一定角度也出現(xiàn)極大的位移輻射。但是正如前面指出的那樣,縱波位移幅度相對(duì)于橫波位移來(lái)說(shuō)是較小的。螺旋型EMAT同時(shí)產(chǎn)生的這兩種聲場(chǎng)與壓電超聲的活塞聲場(chǎng)顯然不同。

圖7 螺旋型EMAT 橫波輻射聲場(chǎng)指向性與縱波輻射聲場(chǎng)指向性Fig.7 The directivity of the shear wave radiation field and the directivity of the longitudinal wave radiation field of the helical EMAT

圖8 分別施加r 和z 方向洛倫茲力的橫波位移幅度與縱波位移對(duì)比Fig.8 Comparison of shear wave displacement amplitude and longitudinal wave displacement with applied Lorentz force in r and z directions,respectively

2 螺旋型EMAT輻射聲場(chǎng)的格林函數(shù)法分析

2.1 遠(yuǎn)場(chǎng)半空間格林函數(shù)

格林函數(shù)法分析輻射聲場(chǎng)時(shí),首先要求出點(diǎn)力源產(chǎn)生的位移場(chǎng)響應(yīng),即格林函數(shù),然后任意一點(diǎn)的位移可以表示為力源分布函數(shù)與格林函數(shù)的卷積,也就是:

其中,G(R-R′,ω)表示R′處的點(diǎn)力源f(R′,ω)在點(diǎn)R的格林函數(shù)表達(dá)式。

在遠(yuǎn)場(chǎng)近似下,可以將半無(wú)限大空間遠(yuǎn)場(chǎng)位移場(chǎng)用遠(yuǎn)場(chǎng)半空間的格林函數(shù)Ghf(R-R′,ω)表示為

其中,cS,L和ZS,L分別為橫、縱波聲速和聲阻抗,MS(θ)和ML(θ)為半空間橫、縱波Miller-Pussy 因子[6],若力源f(R′,ω)沿著r方向,有

其中,κ為縱波與剪切波波數(shù)比,κ=kP/kS。

可以看到,公式(5)相當(dāng)于對(duì)力源f(R′,ω)進(jìn)行空間傅里葉變換。因此,求解洛倫茲力產(chǎn)生的聲束指向性之前,需要對(duì)力源進(jìn)行空間傅里葉變換。

2.2 表面等效洛倫茲力的空間傅里葉變換

洛倫茲力的計(jì)算,需要得到永磁鐵的磁感應(yīng)強(qiáng)度和渦流的強(qiáng)度。在理論上可以借助Dodd &Deeds 模型得到解析的結(jié)果[7],或者利用Tree方法得到級(jí)數(shù)解[5]。但是這兩種方法在具體操作過(guò)程中也不能得到絕對(duì)的解析結(jié)果,對(duì)于Dood &Deeds模型,需要用到無(wú)窮積分,這個(gè)積分不能夠解析得到,因此只能進(jìn)行數(shù)值積分。而采用Tree 方法,涉及到無(wú)窮項(xiàng)級(jí)數(shù),因此也需要進(jìn)行數(shù)值求解才能夠得到相應(yīng)的結(jié)果。因此在這里直接采用有限元計(jì)算得出的洛倫茲力分布結(jié)果,然后將該洛倫茲力進(jìn)行數(shù)值傅里葉變換。另外考慮到洛倫茲力的分布特點(diǎn),采用門(mén)函數(shù)替代有限元計(jì)算得到的洛倫茲力結(jié)果,同樣也進(jìn)行傅里葉變換,并給出在這種簡(jiǎn)化處理下得到的EMAT聲場(chǎng)指向性的遠(yuǎn)場(chǎng)解析表達(dá)式,并比較兩者的差別。

有限元和理論分析都表明,當(dāng)頻率為1 MHz,在鋁板試樣中產(chǎn)生的渦流的深度(<0.1 mm)相對(duì)于波長(zhǎng)是較小的。因此,可以將本來(lái)作為體力作用的洛倫茲力,沿著板厚方向進(jìn)行積分,得到作為表面力作用的洛倫茲力進(jìn)行等效簡(jiǎn)化。洛倫茲力作為體力施加和作為表面力施加得到的位移幅度的有限元計(jì)算結(jié)果表明兩者相差無(wú)幾。通過(guò)將體力近似替代為表面力,在對(duì)洛倫茲力進(jìn)行空間傅里葉變換的時(shí)候,就從三維變換簡(jiǎn)化為了二維變換,大大簡(jiǎn)化了計(jì)算量。

J1為一階第一類(lèi)貝塞爾函數(shù)。

圖9 分布及其極坐標(biāo)傅里葉變換Fig.9 distribution and its Polar Fourier Transform

2.3 格林函數(shù)簡(jiǎn)化解析結(jié)果

在利用格林函數(shù)法分析橫波輻射聲場(chǎng)過(guò)程中,可以看到,實(shí)際上洛倫茲力分布與門(mén)函數(shù)十分相似,因此可以考慮用門(mén)函數(shù)表達(dá)式替代有限元結(jié)果得到的洛倫茲力分布,從而得到簡(jiǎn)化的格林函數(shù)解析模型。圖11(a)給出了簡(jiǎn)化的門(mén)函數(shù)分布的洛倫茲力,圖11(b)給出了有限元洛倫茲力傅里葉變換與簡(jiǎn)化洛倫茲力傅里葉變換的結(jié)果,可以看到簡(jiǎn)化的洛倫茲力與FEM 計(jì)算得到的洛倫茲力傅里葉變換基本是一致的。

圖11 簡(jiǎn)化洛倫茲力和空間傅里葉變換結(jié)果Fig.11 Simplified Lorentz force and spatial Fourier transform result

將洛倫茲力分布看成是幅度為F0,門(mén)位分別為a、b的門(mén)函數(shù),那么可以通過(guò)上述變換得到

其中,J0、J1和H0、H1分別為Bessel 和Struve 函數(shù)。進(jìn)一步,可以得到聲場(chǎng)的遠(yuǎn)場(chǎng)指向性解析表達(dá)式:

利用簡(jiǎn)化格林函數(shù)方法(Simplified Green’s function method,SGFM)可以計(jì)算螺旋型EMAT產(chǎn)生的橫波聲束指向性。圖12對(duì)比了用3 種方法計(jì)算得到的聲束指向性,可以看到SGFM 與其他兩種方法求解結(jié)果一致性很好。

圖12 SGFM、GFM 和FEM 計(jì)算得到的橫波聲束指向性對(duì)比Fig.12 Comparison of shear wave beam directivity calculated bySGFM,GFM and FEM

SGFM 相較于基于有限元的GFM 和FEM 來(lái)說(shuō),在計(jì)算聲束指向性需要的計(jì)算時(shí)間更少,計(jì)算更方便。利用簡(jiǎn)化格林函數(shù)模型可以很方便地計(jì)算偏離角θ隨著一些線圈參數(shù)的變化情況,如圖13所示。可以看出,通過(guò)增加頻率和加大線圈的內(nèi)徑可以減小聲束的偏離角,進(jìn)而使得聲束相對(duì)集中一些。

圖13 偏離角θ 隨頻率的變化關(guān)系與隨線圈內(nèi)徑的變化關(guān)系Fig.13 Deviation angle θ versus frequency versus coil inner diameter

3 螺旋型EMAT輻射聲場(chǎng)實(shí)驗(yàn)驗(yàn)證

利用EMAT檢測(cè)厚度為50 mm、材料為鋁合金的長(zhǎng)條形試樣,在試樣的下方有一平底孔。在試樣上方橫向移動(dòng)EMAT,以平底孔中心正上方試樣表面為原點(diǎn),平底孔反射回波大小與橫向位置的關(guān)系如圖14 所示。實(shí)驗(yàn)中所用的頻率為1.6 MHz,其他參數(shù)與有限元仿真一致。

圖14 實(shí)驗(yàn)所用平底孔試樣Fig.14 Flat-bottom hole sample used in the experiment

圖15的實(shí)驗(yàn)結(jié)果顯示,在偏離平底孔兩側(cè)進(jìn)行檢測(cè)的時(shí)候,會(huì)分別出現(xiàn)回波極大值,兩極大值之間的距離為X=11.7 mm,而平底孔反射面距離上表面距離為h=45 mm,利用三角函數(shù)關(guān)系可以計(jì)算得到偏離角:

圖15 平底孔反射回波隨橫向位置變化Fig.15 Flat-bottom hole reflected echo as a function of lateral position

采用上述幾種分析方法得到的偏離角計(jì)算結(jié)果,見(jiàn)表1。

表1 實(shí)驗(yàn)與幾種計(jì)算方法得到的偏離角θ 對(duì)比Table 1 Comparison of deviation angles θ obtained from experiments and several calculation methods

從圖16的B掃圖上還可以看出,螺旋型EMAT產(chǎn)生了較強(qiáng)的橫波一次底波信號(hào),同時(shí)也有次強(qiáng)的縱波一次底波信號(hào),以及更弱的縱波在底面經(jīng)一次波型轉(zhuǎn)換成為橫波的回波信號(hào),這些固有回波都可能比缺陷回波信號(hào)強(qiáng),干擾缺陷的識(shí)別;而且由于EMAT的喇叭形聲場(chǎng)特征,一個(gè)缺陷在掃查時(shí)形成了兩個(gè)“缺陷”信號(hào)。實(shí)驗(yàn)結(jié)果與理論分析結(jié)果一致。因此,螺旋型EMAT 的這種復(fù)雜的聲場(chǎng)特征會(huì)使之在缺陷探傷方面的應(yīng)用受到局限。

圖16 平底孔反射回波及B 掃描結(jié)果Fig.16 Flat-bottomed hole reflected echo and B-scan results

4 結(jié)論

通過(guò)FEM 和GFM 研究了螺旋型EMAT 的輻射聲場(chǎng),并提出了將洛倫茲力簡(jiǎn)化為門(mén)函數(shù)的SGFM解析模型,提高了模擬計(jì)算的速度。

對(duì)螺旋型EMAT的輻射聲場(chǎng)分析表明:該聲場(chǎng)中同時(shí)存在有橫波和縱波,并具有不同的指向性,在缺陷檢測(cè)時(shí)還會(huì)出現(xiàn)縱/橫波的模態(tài)轉(zhuǎn)換,這種復(fù)雜的聲場(chǎng)特性不利于對(duì)缺陷的檢測(cè);而且這種喇叭型聲場(chǎng)分布使得探頭正下方的缺陷不能被直接檢出,在掃描檢測(cè)時(shí),同一缺陷會(huì)反映出兩次缺陷信號(hào)。這些復(fù)雜的聲場(chǎng)特性局限了EMAT在缺陷檢測(cè)方面的應(yīng)用。

上述研究方法可進(jìn)一步用于EMAT 探頭的優(yōu)化設(shè)計(jì),從而得到所需要的聲場(chǎng)分布。

猜你喜歡
磁場(chǎng)有限元
西安的“磁場(chǎng)”
為什么地球有磁場(chǎng)呢
文脈清江浦 非遺“磁場(chǎng)圈”
新型有機(jī)玻璃在站臺(tái)門(mén)的應(yīng)用及有限元分析
《磁場(chǎng)》易錯(cuò)易混知識(shí)剖析
基于有限元的深孔鏜削仿真及分析
基于有限元模型對(duì)踝模擬扭傷機(jī)制的探討
磁場(chǎng)的性質(zhì)和描述檢測(cè)題
2016年春季性感磁場(chǎng)
Coco薇(2016年1期)2016-01-11 16:53:24
磨削淬硬殘余應(yīng)力的有限元分析
主站蜘蛛池模板: 亚亚洲乱码一二三四区| 国产99久久亚洲综合精品西瓜tv| 亚洲自拍另类| 久久 午夜福利 张柏芝| 制服丝袜一区| 日韩av无码精品专区| 国产日韩精品一区在线不卡| 9啪在线视频| 亚洲欧洲日产国码无码av喷潮| 国产成人无码播放| 亚洲—日韩aV在线| 国产清纯在线一区二区WWW| 欧美成人日韩| 亚洲Va中文字幕久久一区| 91久久偷偷做嫩草影院| 欧美日韩中文字幕在线| 日韩天堂视频| 亚洲av无码久久无遮挡| av色爱 天堂网| 日韩无码视频播放| 国产欧美视频在线| 青青青国产免费线在| 亚洲天堂日韩在线| 成人福利在线视频| 久热中文字幕在线观看| 欧美日本在线| 亚洲国产精品无码AV| 97超碰精品成人国产| 国产精品三级专区| 欧美一级在线播放| 国产午夜在线观看视频| 欧美精品aⅴ在线视频| 97视频在线精品国自产拍| 久久亚洲美女精品国产精品| 亚洲精品自在线拍| 在线观看欧美国产| 91成人在线免费视频| 中文字幕欧美日韩| 亚洲黄网在线| 亚洲日韩AV无码精品| 国产99视频精品免费视频7| 国产成人精品男人的天堂| 久久久久青草大香线综合精品| 无码精品国产dvd在线观看9久 | 欧美成人h精品网站| 114级毛片免费观看| 综合色在线| 77777亚洲午夜久久多人| 欧美午夜网站| 国产麻豆福利av在线播放| 人妻无码AⅤ中文字| 日韩精品无码免费一区二区三区| 亚洲成a人片77777在线播放| 久久五月视频| 国产视频a| 国产成人三级在线观看视频| 亚洲精品第五页| 国产成人啪视频一区二区三区| 亚洲免费播放| 亚洲国产成人综合精品2020| 黄色片中文字幕| 91小视频在线观看免费版高清| 亚洲一区免费看| 国模私拍一区二区三区| 日韩av手机在线| 国产精品视频3p| 日韩毛片基地| 国产精品女熟高潮视频| 精品三级在线| 麻豆国产精品视频| 国产探花在线视频| 伊人精品视频免费在线| 免费在线播放毛片| 欧美一区福利| a级毛片免费在线观看| 国产XXXX做受性欧美88| 欧美色亚洲| 亚洲欧美极品| 精品一区二区三区水蜜桃| 中文字幕在线免费看| 日韩黄色在线| 9cao视频精品|