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

圓柱形推拉式電磁鐵的電磁力數(shù)值算法

2023-11-23 07:20:16魯夢(mèng)昆張俊洪易祥烈袁志方
科學(xué)技術(shù)與工程 2023年30期
關(guān)鍵詞:有限元

魯夢(mèng)昆,張俊洪,易祥烈*,袁志方

(1.海軍工程大學(xué)電氣工程學(xué)院,武漢 430033; 2.湖北工程學(xué)院機(jī)械工程學(xué)院,孝感 432000)

推拉式電磁鐵是一種典型的電磁能與機(jī)械能能量轉(zhuǎn)換裝置,當(dāng)給線圈中施加勵(lì)磁電流時(shí),形成的電磁場(chǎng)將使線圈內(nèi)部軟磁性材質(zhì)的鐵芯磁化,在電磁力的作用下,鐵芯和推桿會(huì)對(duì)外輸出推力和直線運(yùn)動(dòng)。在一些需求短距離直線運(yùn)動(dòng)控制的領(lǐng)域,相比于旋轉(zhuǎn)電機(jī)及配套運(yùn)動(dòng)轉(zhuǎn)換裝置結(jié)合的方案,螺線管電磁鐵裝置在制造成本、體積、質(zhì)量、響應(yīng)速度等性能方面具有明顯優(yōu)勢(shì)。鑒于這些特點(diǎn),推拉式電磁鐵被廣泛應(yīng)用于工業(yè)自動(dòng)化領(lǐng)域,包括電磁閥、繼電器、電磁振動(dòng)器等多種形式,在各種應(yīng)用種準(zhǔn)確的電磁力控制都是至關(guān)重要的[1-2]。

推拉式電磁鐵作為一項(xiàng)基礎(chǔ)應(yīng)用技術(shù)一直以來都是中外科研工作者的重點(diǎn)研究方向[3]。近些年隨著工業(yè)科技的發(fā)展和對(duì)電能源應(yīng)用的愈發(fā)重視,螺線管電磁鐵相繼出現(xiàn)了一些新的應(yīng)用方式。其中包括文獻(xiàn)[4]中介紹的一種電磁鉚接機(jī),主要功能是為大飛機(jī)的蒙皮外殼板材進(jìn)行鉚釘鉚接;文獻(xiàn)[5]中介紹的一種電磁卷邊機(jī);文獻(xiàn)[6]中分析了一種直線柱塞泵,液壓動(dòng)力直接通過施加在活塞上的電磁力得到,而不是由旋轉(zhuǎn)動(dòng)力轉(zhuǎn)換而來,更適用于小型、結(jié)構(gòu)緊湊的工作環(huán)境;相比于化學(xué)能,電能在太空中具有可再生的優(yōu)勢(shì),是衛(wèi)星運(yùn)行的主要?jiǎng)恿υ?推拉式電磁鐵體積小的特點(diǎn)非常適合在太空中應(yīng)用[7],如文獻(xiàn)[8]中介紹了一種小型衛(wèi)星的太空對(duì)接方案,文獻(xiàn)[9]中提出了一種基于單線圈鐵芯螺線管的立方體納衛(wèi)星微型電磁對(duì)接硬件重構(gòu)設(shè)計(jì)方案。雖然推拉式電磁鐵應(yīng)用廣泛,但因?yàn)殍F芯是活動(dòng)的,電磁力、磁芯位置、電流等電磁數(shù)據(jù)之間存在強(qiáng)非線性關(guān)系,解析計(jì)算非常困難[10-11],在工程應(yīng)用中的動(dòng)態(tài)性能難以準(zhǔn)確控制。

在推拉式電磁鐵電磁場(chǎng)計(jì)算模型的研究方面,無鐵芯的螺線管磁場(chǎng)容易根據(jù)Biot-Savart定律計(jì)算,而含鐵芯的螺線管線圈則計(jì)算過程較為復(fù)雜,成熟的理論成果較少,目前可用的方法包括半解析計(jì)算法和有限元法(finite element method,FEM)、其中半解析法包括集總參數(shù)磁路模型法和磁場(chǎng)子域方程法[12]。有限元法計(jì)算存在計(jì)算速度慢,變量之間對(duì)應(yīng)關(guān)系不清晰的問題[13];集總參數(shù)磁路模型法主要用以計(jì)算磁場(chǎng)的平均值,而不是一個(gè)區(qū)域的詳細(xì)分布,因而計(jì)算準(zhǔn)確度不高[14];文獻(xiàn)[15-16]中從理論層面解釋了子域磁矢量方程計(jì)算鐵芯螺線管的磁感應(yīng)強(qiáng)度、電感、電磁力的過程。

基于子域磁矢量方程的數(shù)值算法相比于有限元法具有計(jì)算速度快、可程序化等優(yōu)點(diǎn),能完成一些有限元不能直接處理的計(jì)算任務(wù)。現(xiàn)利用數(shù)值算法分析一種不帶鐵磁性外殼的推拉式電磁鐵,含特殊函數(shù)的磁矢量勢(shì)方程通過MATLAB編程計(jì)算數(shù)值解,在計(jì)算得到電感和電磁力后對(duì)比有限元仿真結(jié)果和實(shí)驗(yàn)結(jié)果說明該方法的計(jì)算準(zhǔn)確度。

1 電磁力數(shù)值算法建模

推拉式電磁鐵主要由勵(lì)磁線圈、軟磁性材質(zhì)的鐵芯和彈簧組成。其中彈簧的作用主要是使鐵芯具有復(fù)位功能,彈簧力的方向與電磁力方向相反。當(dāng)線圈中通以任意方向的勵(lì)磁電流時(shí),鐵芯受到電磁力的作用并產(chǎn)生向線圈中心移動(dòng)的趨勢(shì)。推拉式電磁鐵的結(jié)構(gòu)圖如圖1所示。

圖1 推拉式電磁鐵的結(jié)構(gòu)圖

由于鐵芯頂端的圓環(huán)結(jié)構(gòu)不會(huì)進(jìn)入線圈,距離線圈較遠(yuǎn),下方的推桿為滲碳無磁不銹鋼材質(zhì),因此可將推拉式電磁鐵的鐵芯和線圈視為標(biāo)準(zhǔn)圓柱體進(jìn)行分析。

1.1 磁矢量勢(shì)方程

磁場(chǎng)子域磁矢量勢(shì)方程的數(shù)學(xué)模型中假設(shè)了鐵芯的磁導(dǎo)率無窮大,鐵芯上的磁力線均垂直于鐵芯表面[17],即鐵芯邊界處的磁感應(yīng)強(qiáng)度切向分量為零。以鐵芯的中軸線為中心線建立柱狀旋轉(zhuǎn)坐標(biāo)系,以求解域圓柱體的3個(gè)面為外部邊界條件,推拉式電磁鐵磁場(chǎng)子域劃分及各邊界條件如圖2所示。

z為螺線管的軸向方向;r為螺線管的徑向方向;R1為鐵芯半徑;R2為線圈內(nèi)半徑;R3為線圈外半徑;R4為求解域半徑;Z1為線圈的左邊界;Z2為線圈的右邊界;Z3為鐵芯的左邊界;Z4為鐵芯的右邊界;Z5為求解域長(zhǎng)度;AⅠ、AⅡ、AⅢ、AⅣ、AⅤ分別為5個(gè)區(qū)域的磁矢量勢(shì)

為了使求解結(jié)果更準(zhǔn)確,各邊界應(yīng)遠(yuǎn)離線圈和鐵芯。由麥克斯韋方程可知

-?B=?2A=-μ0J

(1)

式(1)中:?為倒數(shù)算符;B為磁感應(yīng)強(qiáng)度;A為磁矢量勢(shì);μ0為空氣中的磁導(dǎo)率,取4π×10-7H/m;J為線圈中的電流密度。

注意磁場(chǎng)呈圓柱狀具有旋轉(zhuǎn)對(duì)稱性,磁矢量勢(shì)僅有周向θ方向的分量,其大小與r、z的坐標(biāo)相關(guān),可表示為

(2)

柱狀坐標(biāo)系下螺線管磁感應(yīng)強(qiáng)度的3個(gè)分量[17]分別為

(3)

聯(lián)立式(1)~式(3)得到磁矢量勢(shì)方程滿足

(4)

式(4)中:Aθ=A(r,z);電流密度J表示為與坐標(biāo)點(diǎn)相關(guān)的函數(shù)J(r,z),以便區(qū)分電流密度為零的空氣區(qū)域。

利用分離變量法求解齊次微分方程,可得磁矢量勢(shì)的通解為

(5)

式(5)中:c1~c4為待定常數(shù);I1為第一類一階虛宗量Bessel函數(shù);K1為第二類一階虛宗量Bessel函數(shù)[18](也被稱為第二類一階修正的Bessel函數(shù));n為級(jí)數(shù)解的諧波項(xiàng)數(shù);αn為與級(jí)數(shù)相關(guān)的系數(shù)。

根據(jù)圓柱坐標(biāo)系的磁矢量勢(shì)在外部邊界處I1和K1的斂散性可分別得到簡(jiǎn)化后的各子區(qū)域磁矢量勢(shì)通解形式為

(6)

式(6)中:

(7)

根據(jù)線圈的電感儲(chǔ)能和磁矢量勢(shì)的關(guān)系,即

(8)

式(8)中:Lm為帶鐵芯的電感;I為線圈中的電流;V為導(dǎo)體的體積。

可推導(dǎo)出電感的計(jì)算式,設(shè)線圈中的電流均勻分布于線圈的矩形截面上,即電流密度為電流除以矩形截面的面積,推導(dǎo)出含鐵芯的電感為

(9)

式(9)中:

(10)

式(10)中:L為線圈長(zhǎng)度;N為線圈匝數(shù);F為超幾何函數(shù)。

由式(9)和式(10)可知,當(dāng)線圈和鐵芯結(jié)構(gòu)尺寸已知時(shí),只要知曉鐵芯與線圈的相對(duì)位置和勵(lì)磁電流即可通過子域磁矢量勢(shì)方程求得對(duì)應(yīng)的電感值。

1.2 電磁力方程

在得到電感曲線后,可通過電感梯度變化計(jì)算鐵芯受到的電磁力。使用虛擬位移法求解運(yùn)動(dòng)學(xué)過程,即

(11)

式(11)中:W為磁感應(yīng)能;L(x)為與可動(dòng)鐵芯位置x(t)有關(guān)的電感值;i(t)為線圈的勵(lì)磁電流。

鐵芯在線圈軸向上受到的電磁力為

(12)

實(shí)際上得到電磁力曲線的同時(shí),也容易計(jì)算出鐵芯的加速度和速度,鐵芯的運(yùn)動(dòng)方程為

(13)

式(13)中:m為鐵芯的質(zhì)量;x(0)為彈丸與線圈的初始相對(duì)位置。

由電磁力方程可知,在得到電感的情況下可計(jì)算鐵芯的運(yùn)動(dòng)狀態(tài),因而可用于控制電磁閥閥芯[20]。不限制鐵芯的運(yùn)動(dòng)時(shí),可將電磁鐵視為磁阻型電磁加速器,可根據(jù)電流和電感梯度的變化計(jì)算出鐵芯的運(yùn)動(dòng)過程,也能根據(jù)電感的變化逆向計(jì)算電流,形成電磁加速過程的動(dòng)態(tài)閉環(huán)計(jì)算[21]。

2 數(shù)值計(jì)算模型的仿真驗(yàn)證

為了驗(yàn)證電感計(jì)算的準(zhǔn)確性,構(gòu)建如圖1所示的推拉式電磁鐵的ANSOFT/Maxwell有限元仿真模型,同時(shí)編寫MATLAB數(shù)值計(jì)算程序,電磁鐵的各結(jié)構(gòu)參數(shù)如表1所示。

表1 圓柱形推拉式電磁鐵的建模參數(shù)

2.1 磁矢量勢(shì)對(duì)比分析

從圖1可知,推拉式電磁鐵鐵芯下移8.2 mm時(shí),鐵芯和線圈的相對(duì)位置為-17.3 mm。根據(jù)式(6)使用MATLAB編程求解到各子域磁矢量勢(shì)方程的7個(gè)未知積分常數(shù)后,可繪制磁矢量勢(shì)分布圖。當(dāng)給定線圈中的勵(lì)磁電流為4 A時(shí),子域法計(jì)算的磁矢量勢(shì)分布圖與有限元結(jié)果對(duì)比如圖3所示,其中有限元中鐵芯的磁導(dǎo)率分別設(shè)置為電工軟鐵的B-H曲線和定值10 000進(jìn)行仿真。

圖3 FEM與子域法計(jì)算的磁矢量勢(shì)的對(duì)比

當(dāng)勵(lì)磁電流為4 A時(shí),3個(gè)磁矢量勢(shì)峰值數(shù)據(jù)分別為3.608 8×10-3、3.585 6×10-3和3.599 4×10-3Wb/m,相差較小。實(shí)際上勵(lì)磁電流的大小與磁飽和程度必然相關(guān),為了具體分析磁矢量勢(shì)的計(jì)算準(zhǔn)確度,給線圈施加不同大小的勵(lì)磁電流,對(duì)比有限元的分析結(jié)果如圖4所示。

圖4 磁矢量勢(shì)峰值對(duì)比

由圖4可知,隨著電流的增加,有限元中磁導(dǎo)率設(shè)置為10 000時(shí)和子域法計(jì)算的磁矢量勢(shì)峰值基本相等,這對(duì)應(yīng)了子域磁矢量勢(shì)方程建立時(shí)鐵芯磁導(dǎo)率為無窮大的假設(shè),說明數(shù)值算法的計(jì)算結(jié)果非常近似地等于有限元中磁芯磁導(dǎo)率設(shè)置很大的情況。而當(dāng)有限元中鐵芯材質(zhì)設(shè)置為電工軟鐵計(jì)算的磁矢量勢(shì)峰值不會(huì)線性增加,而是按照B-H曲線的趨勢(shì)增加,磁芯磁飽和程度逐漸加大。由此可以判斷電磁鐵鐵芯下移8.2 mm時(shí),想要子域磁矢量勢(shì)方程計(jì)算結(jié)果具有較準(zhǔn)確的計(jì)算結(jié)果,線圈中的電流不宜超過6 A。

子域磁矢量勢(shì)方程計(jì)算較準(zhǔn)確的電流閾值必然與鐵芯的位置相關(guān),可在推拉式電磁鐵的鐵芯工作行程內(nèi)確定,以便明確該數(shù)值算法能準(zhǔn)確計(jì)算電磁力的使用條件。

2.2 電感對(duì)比分析

研究中可通過電感曲線計(jì)算電磁力曲線,直接分析不同勵(lì)磁電流情況下子域磁矢量勢(shì)方程計(jì)算電磁力的準(zhǔn)確度。由于子域磁矢量勢(shì)方程中磁導(dǎo)率為無窮大時(shí),電感不會(huì)隨勵(lì)磁電流變化,計(jì)算電感曲線全部重合。因此重點(diǎn)對(duì)比鐵芯工作行程內(nèi)電工軟鐵材質(zhì)鐵芯的電感曲線即可知曉電感計(jì)算準(zhǔn)確度。鐵芯在線圈中心線上每移動(dòng)0.5 mm取1個(gè)數(shù)據(jù)點(diǎn),電感曲線對(duì)比如圖5所示。

圖5 不同勵(lì)磁電流時(shí)的等效電感曲線對(duì)比

由圖5可知,當(dāng)電流為1~4 A時(shí),無論鐵芯在任何位置,各電感曲線高度重合,可以推斷此階段鐵芯未發(fā)生磁飽和。當(dāng)電流達(dá)到5 A時(shí),僅當(dāng)鐵芯中心在線圈中心的-10~10 mm區(qū)域內(nèi)有所偏差,且鐵芯與線圈相對(duì)位置越近誤差越大。鐵芯在-17.3 mm 位置時(shí),勵(lì)磁電流為6 A時(shí)磁芯發(fā)生了輕微的磁飽,但誤差只有2.085%,可以認(rèn)為該電流值即為子域磁矢量勢(shì)方程計(jì)算推拉式電磁鐵電感較準(zhǔn)確的閾值。

2.3 電磁力對(duì)比分析

得到電感曲線后可根據(jù)式(12)計(jì)算電磁力,與有限元中磁芯為電工軟鐵的仿真結(jié)果對(duì)比如圖6所示,由于有限元中默認(rèn)z軸向上方向?yàn)檎?所以鐵芯向下移動(dòng)過程中電磁力為先負(fù)后正。

圖6 不同勵(lì)磁電流時(shí)的電磁力曲線對(duì)比

從圖6(b)中可以看出,數(shù)值算法計(jì)算的電磁力曲線呈波浪狀,這符合磁矢量勢(shì)方程諧波次數(shù)和函數(shù)的特征。推拉式電磁鐵鐵芯相對(duì)位置-17.3 mm時(shí),有限元和子域法的電磁力對(duì)比及和相對(duì)誤差如圖7所示,為了方便與實(shí)驗(yàn)中力傳感器受到的壓力作對(duì)比,電磁力以向下方向?yàn)檎?/p>

圖7 電磁力對(duì)比

由圖7可知,當(dāng)電流為1~5 A時(shí),數(shù)值算法計(jì)算的電磁力相比于考慮鐵芯材質(zhì)B-H曲線的有限元結(jié)果相差不超過3%,當(dāng)電流為6 A時(shí)的計(jì)算偏差約為5%,繼續(xù)增大電流將導(dǎo)致誤差迅速增加。

3 數(shù)值計(jì)算模型的實(shí)驗(yàn)驗(yàn)證

3.1 實(shí)驗(yàn)平臺(tái)搭建

為了驗(yàn)證磁場(chǎng)子域磁矢量勢(shì)方程數(shù)值算法計(jì)算推拉式電磁鐵電磁力的準(zhǔn)確程度,搭建了如圖8所示的電磁力測(cè)試實(shí)驗(yàn)平臺(tái),在實(shí)驗(yàn)過程中拆卸彈簧以剔除彈簧力的影響。為了獲得準(zhǔn)確的測(cè)試效果,測(cè)試過程中應(yīng)當(dāng)確保鐵芯相對(duì)于線圈不會(huì)在電磁力作用下發(fā)生位移,力傳感器應(yīng)當(dāng)有足夠的安裝預(yù)緊力。在鋼質(zhì)安裝架與螺線管之間裝夾木墊可盡量避免安裝架的鐵磁性對(duì)螺線管磁場(chǎng)產(chǎn)生影響。此研究中的實(shí)驗(yàn)平臺(tái)布置為縱向,可有效減少摩擦力和磁芯偏心等問題給計(jì)算帶來的誤差,使分析更準(zhǔn)確。為了避免力傳感器自身的重力和安裝預(yù)緊力對(duì)測(cè)試結(jié)果帶來誤差,測(cè)試前先對(duì)變送器調(diào)零,該操作也同時(shí)剔除了鐵芯的重力。

圖8 電磁力測(cè)試實(shí)驗(yàn)平臺(tái)

3.2 數(shù)據(jù)對(duì)比分析

實(shí)驗(yàn)過程中利用可調(diào)直流電源給線圈施加1~8 A不同大小的直流電(螺線管的額定工作電流約7.2 A),電流數(shù)據(jù)通過直流電流傳感器檢測(cè),力傳感器檢測(cè)的電壓與力轉(zhuǎn)換比為10 N/V。示波器中采集到的電流和電磁力與有限元仿真和子域法計(jì)算結(jié)果對(duì)比如表2所示,其中有限元和MATLAB程序中的電流激勵(lì)按照實(shí)驗(yàn)中實(shí)測(cè)電流給定。另外,圖9中顯示了勵(lì)磁電流分別為約為2.982 A和6.997 A兩種情況的示波器波形。

表2 仿真、計(jì)算及實(shí)驗(yàn)獲得的電磁力對(duì)比

當(dāng)電流較小時(shí),傳感器等本身的誤差相對(duì)較大,因而不適合以該數(shù)據(jù)結(jié)果為標(biāo)準(zhǔn)分析基于子域磁矢量勢(shì)方程的半解析數(shù)值算法的誤差。有限元中分析模型則過于理想化,以實(shí)驗(yàn)結(jié)果為標(biāo)準(zhǔn)分析誤差更合理。實(shí)驗(yàn)中當(dāng)電流為2~5 A時(shí),數(shù)值算法計(jì)算誤差不超過3%,當(dāng)電流約為6 A時(shí)誤差約為7.25%,隨著電流繼續(xù)增大數(shù)值算法的誤差迅速增加。因此判斷數(shù)值算法的適用電流閾值應(yīng)不高于 5~6 A,在此條件下數(shù)值算法的計(jì)算結(jié)果具有可信度。

4 結(jié)論

推拉式電磁鐵具有結(jié)構(gòu)簡(jiǎn)單、制造成本低、體積小、響應(yīng)速度快等優(yōu)點(diǎn),雖然其已被廣泛應(yīng)用于各行各業(yè),但由于磁場(chǎng)的非線性變化,解析計(jì)算困難,現(xiàn)階段的計(jì)算基本上都是靠有限元來處理。與有限元法相比,基于數(shù)學(xué)方程的半解析法具有計(jì)算過程清晰、可程序化、算法編輯靈活等優(yōu)點(diǎn),可實(shí)現(xiàn)一些有限元無法處理的計(jì)算,是理論與應(yīng)用研究的發(fā)展趨勢(shì)。

研究了柱狀坐標(biāo)系下的磁場(chǎng)子域方程,通過MATLAB編程實(shí)現(xiàn)了一種推拉式電磁鐵電磁特性數(shù)據(jù)的計(jì)算,包括磁矢量勢(shì)、電感和電磁力等。數(shù)值算法計(jì)算的結(jié)果與有限元中設(shè)置鐵芯磁導(dǎo)率為10 000時(shí)的結(jié)果高度重合,但當(dāng)鐵芯線設(shè)置為電工軟鐵的B-H曲線時(shí),隨著電流的加大鐵芯會(huì)發(fā)生磁飽和。當(dāng)磁飽和問題不嚴(yán)重時(shí),數(shù)值算法的計(jì)算結(jié)果具有高準(zhǔn)確度,可等效于有限元法;隨著電流的增加,發(fā)生磁飽和后誤差迅速增加。在實(shí)際使用時(shí),可對(duì)比空載情況下有限元計(jì)算的磁矢量勢(shì)或電感,提前獲取電流大小與誤差的變化規(guī)律,明確數(shù)值算法計(jì)算較準(zhǔn)確的電流閾值,以此避免磁飽和導(dǎo)致計(jì)算誤差過大的問題。

經(jīng)過研究,基于子域磁矢量勢(shì)方程的半解析數(shù)值算法的主要缺點(diǎn)包括:僅適用于標(biāo)準(zhǔn)圓柱狀的含鐵芯或空心螺線管;磁芯磁導(dǎo)率太低或發(fā)生磁飽和后計(jì)算準(zhǔn)確度降低。在后續(xù)的研究中有必要研究可動(dòng)鐵芯動(dòng)態(tài)電流情況下磁飽和修正算法。

猜你喜歡
有限元
基于擴(kuò)展有限元的疲勞裂紋擴(kuò)展分析
非線性感應(yīng)加熱問題的全離散有限元方法
TDDH型停車器制動(dòng)過程有限元分析
新型有機(jī)玻璃在站臺(tái)門的應(yīng)用及有限元分析
基于I-DEAS的履帶起重機(jī)主機(jī)有限元計(jì)算
基于有限元模型對(duì)踝模擬扭傷機(jī)制的探討
10MN快鍛液壓機(jī)有限元分析
磨削淬硬殘余應(yīng)力的有限元分析
基于SolidWorks的吸嘴支撐臂有限元分析
箱形孔軋制的有限元模擬
上海金屬(2013年4期)2013-12-20 07:57:18
主站蜘蛛池模板: 欧美人与牲动交a欧美精品 | 不卡视频国产| 亚洲综合在线网| 人妻丰满熟妇αv无码| 色婷婷亚洲十月十月色天| 国产极品粉嫩小泬免费看| 99久久国产精品无码| 亚洲床戏一区| 久久久亚洲色| lhav亚洲精品| 五月天天天色| www.亚洲一区二区三区| 欧美日韩亚洲综合在线观看| 99久久亚洲综合精品TS| 欧美a在线| 国产麻豆另类AV| 久久黄色免费电影| 日韩在线第三页| 日韩一区精品视频一区二区| 99re在线视频观看| 91在线播放国产| 日韩123欧美字幕| 日韩精品毛片| 国产日韩欧美在线播放| 国产99免费视频| 99视频国产精品| 久久男人资源站| 日本成人福利视频| 国产成人精品一区二区免费看京| 亚洲人成网18禁| 国产欧美性爱网| 精品久久人人爽人人玩人人妻| 精品日韩亚洲欧美高清a| 日韩精品亚洲一区中文字幕| 日本国产精品一区久久久| 欧美日韩理论| 麻豆精品在线| 国产欧美精品午夜在线播放| 国产在线精品网址你懂的| 亚洲高清无在码在线无弹窗| 亚洲国产精品美女| 久久精品人妻中文视频| 成人欧美日韩| 国产亚洲一区二区三区在线| 色综合久久88色综合天天提莫| 亚洲免费人成影院| 最新国产精品第1页| 毛片网站观看| 青青草原国产免费av观看| 国产视频a| 久久黄色一级视频| 动漫精品啪啪一区二区三区| 波多野吉衣一区二区三区av| 999精品色在线观看| 99久久人妻精品免费二区| 素人激情视频福利| 91精品专区| 国产男人天堂| 五月天久久婷婷| 真实国产乱子伦视频| 亚洲最大在线观看| 久久超级碰| 高清欧美性猛交XXXX黑人猛交| 欧美日韩91| 亚洲日韩图片专区第1页| 精品视频福利| 亚洲无码精品在线播放| 亚洲va欧美va国产综合下载| 亚洲激情区| 草逼视频国产| 国产视频欧美| 国产福利免费视频| 狼友视频国产精品首页| 国产丝袜第一页| 亚洲啪啪网| 在线视频亚洲色图| 又爽又大又黄a级毛片在线视频| 免费AV在线播放观看18禁强制| 天堂中文在线资源| 日本在线国产| 欧美国产在线看| 伊人久久大香线蕉成人综合网|