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

永磁電動懸浮系統三維解析建模與電磁力優化分析

2021-03-16 08:36:14李冠醇
電工技術學報 2021年5期
關鍵詞:優化模型

巫 川 李冠醇 王 東

(海軍工程大學艦船綜合電力技術國防科技重點實驗室 武漢 430033)

0 引言

永磁電動懸浮系統(Permanent Magnet-Electrodynamic Suspension, PM-EDS)通過永磁體與導體板相對運動在導體板內感應出渦流,與源磁場相互作用產生電磁力使永磁體懸浮于空氣中,具有無摩擦、能耗低、懸浮氣隙大、結構簡單等優點,在磁懸浮交通、電磁發射和航天發射助推領域具有廣泛的應用前景[1-2]。以美國為主的發達國家正抓緊對高速永磁電動懸浮系統關鍵技術進行攻關,美國通用原子能公司(GA)建造了一條基于Inductrack技術的 200m長的全尺寸試驗系統,勞倫斯利弗莫爾國家實驗室(LLNL)在 NASA的資助下建立了應用于航天電磁發射的懸浮系統,并取得一系列的研究成果[3],國內的西南交通大學[2,4]、國防科技大學[3,5]、中科院電工所[6-7]等研究機構也對永磁電動懸浮技術開展了一系列研究,并建立了不同的試驗平臺。電磁力的準確計算是對 Halbach陣列永磁電動懸浮系統進行研究和設計的基礎,結構優化是提高系統性能和減少工程化建設成本的關鍵,但相關的文獻調研表明,在懸浮系統橫向端部效應量化分析和全尺寸優化方面還需要進一步研究。

由于 Halbach永磁陣列的長和寬均有限,直線型永磁電動懸浮系統的磁場在縱向和橫向都會產生畸變。這會引起在直線電機分析中廣泛考慮的端部效應[8],因此電磁場分析在系統研究設計過程中尤其重要。目前常用的電磁場分析方法為有限元法和解析法。其中,有限元法計算較為準確,尤其在一些必須考慮3-D電磁場分布的場合應用較多,但有限元法計算時間較長且對計算設備性能要求較高,不適合在系統分析設計階段使用。針對這個問題,許多學者嘗試采用解析法計算永磁電動懸浮系統的電磁場分布。目前國際上的研究者主要利用磁矢位和傅里葉空間分解法對系統的電磁場進行2-D建模分析[9-10],但基于周期無限長的 2-D模型忽略了系統橫向端部效應和縱向端部效應對電磁力的影響,這使得電磁力計算結果與實際相比有較大誤差;有學者采用導體板電導率修正系數等效橫向端部效應的方法求解懸浮系統所受的電磁力[11],但修正系數只對氣隙磁場的垂直分量進行了矯正,縱向分量未得到矯正,因此電磁力計算的最后結果仍存在較大誤差。陳殷和 Z. B. Jonathan博士等采用二階矢量勢法對系統的三維電磁場進行建模分析[4,12-14],但該模型基于導體板無限寬假設,導體板寬度無法得到優化。王厚生博士對有限寬導體板的三維電磁場進行了分析[6,15],但建模過程中源磁場計算存在較大誤差并且忽略了導體板橫向分界面磁場透射對導體板渦流的影響。P. Subhra博士研究了橫向邊界面源磁場對導體板的影響[16-17],較為準確地計算出了電磁力,但磁輪與直線形永磁陣列相比磁場分布差別較大,同時只研究了永磁體相對于導體板特定位置的電磁解析模型。目前關于系統結構參數對橫向端部效應影響的研究仍存在空白,尤其是導體板寬度和永磁體寬度對系統電磁力的影響尚不明確。由于電磁力解析模型計算的限制,以往對永磁電動懸浮系統優化方面的研究多基于 2-D模型和導體板無限寬假設的3-D模型,無法對全系統尤其是導體板寬度進行優化設計[18-20],而對于永磁電動懸浮系統,導體板鋪設在軌道上,其寬度顯著影響系統成本,因此明確永磁體陣列寬度和導體板寬度對系統橫向端部效應的影響,將導體板寬度納入優化過程,對系統基于工程應用背景下的輕量化和小型化研究具有重要意義。

以圖1所示基于五模塊Halbach陣列的永磁電動懸浮系統為研究對象,建立包含系統橫向端部效應和縱向端部效應的三維電磁場解析模型,通過對比 2-D、3-D解析模型計算結果與 3-D有限元仿真結果,驗證了所建三維解析模型計算結果的準確可靠性。利用三維解析模型對永磁電動懸浮系統進行參數化分析,研究結構參數對浮阻比和浮重比兩個優化指標的影響。最后利用多目標粒子群優化算法(Multi-Objective Particle Swarm Optimization, MOPSO)對系統結構尺寸進行了迭代尋優,通過對比MOPSO設計解和基于參數分析的系統優化設計解,證明本文提出的優化設計方法快速可靠。

圖1 永磁電動懸浮系統3-D結構示意圖Fig.1 3-D structural schematic view of PM-EDS system

1 系統建模及驗證

1.1 模型分析

圖2為五模塊Halbach陣列永磁電動懸浮系統xOy和yOz截面展開圖。系統將所處的空間劃分為Ⅰ~Ⅴ五個區域,各結構的尺寸參數如圖2所示,其中,τ為永磁體極距,h為永磁體厚度,2w0為永磁體寬度,g為氣隙,d導體板厚度,wz為導體板橫向伸出寬度,l0為Halbach陣列周期延拓長度。為方便分析,對系統三維模型做如下假設[1,11,16]:①導體板在x方向無限延伸,在z方向寬度有限;②永磁體陣列位于導體板橫向寬度中心,導體板厚度有限;③Halbach永磁體陣列沿x軸以速度vx運動且遠小于光速,系統電磁場為似穩電磁場;④各區域的物理參數是均勻的、各向同性的。

圖2 永磁電動懸浮系統xOy和yOz截面圖Fig.2 The xOy and yOz section views of PM-EDS

基于Maxwell電磁場方程組和磁矢位[1-2]

經分析推導可得導體板區域的磁矢位微分方程為

式中,μ為磁導率;σ為電導率。對于非導體板域,由于電導率為0,所以其磁矢位方程為

對于導體板域,根據文獻[16-17],將導體板內的感應電流等效為源磁場分別通過邊界T、B和邊界e、r透射入導體板激發的感應電流的疊加結果,則導體板內磁矢位A可表示為

式中,下標 1、2表示源磁場分別通過T、B邊界、e、r邊界與導體板反應產生的感應場的磁矢位。

研究表明,垂直于邊界面的渦流分量并不會因為源磁場的變化有明顯改變[17],且不考慮漏磁情況下垂直于源場透射邊界面的渦流分量不存在[6],因此在理想情況下解析計算可以忽略垂直于源場透射邊界面的渦流分量,可以得出

1.2 邊界條件

假定導體板區域和空氣域的磁導率相同,在T邊界面,根據安培環路定理,磁場強度的切向分量相等,即有

根據磁通連續性定理有

式中,上標r表示導體板渦流在區域Ⅱ激發的感應磁場;上標s表示永磁體激發的源磁場。B、e、r邊界面的磁場銜接條件同邊界面T相同。

為了實現磁矢位A各分量解耦,引入庫侖規范,則有

根據文獻[16]在導體板各邊界面上有

其中,nⅡ為邊界面的方向向量,在x方向上有

其他各區域的外層邊界條件為

1.3 各區域磁矢位方程通解

將式(5)代入式(2),運用分離變量法[21]先求出源磁場通過上、下邊界T、B與導體板反應產生的渦旋磁場磁矢位,考慮到z方向的渦流分布情況及式(9)和式(10),則區域Ⅱ中各分量的雙重傅里葉級數解[16]為

式中,2.5τ+2l0為x方向上考慮系統縱向端部效應的修正周期,認為縱向端部磁場在l0長度內衰減為0;2wd=2(w0+wz)為導體板的橫向寬度。同理,區域Ⅱ中各分量的通解為

對于非導體Ⅰ區域,考慮導體板中渦流在此區域產生的磁場要與區域Ⅱ中的場一一對應[16-17],由式(11)可解得

區域Ⅲ、Ⅳ、Ⅴ內磁矢位解的形式同區域Ⅰ類似。

1.4 源場計算

由式(6)和式(7)可以看出,在邊界條件中需要源場值,采用面電流法[22]對五模塊永磁體構成的 Halbach陣列的外部磁場進行建模求解,最后對源場進行雙重傅里葉級數變換。分析邊界條件發現,源場磁通密度的數學表達形式與磁矢位A的各分量有關。為使源場與感應場逐項關聯便于求解,源場磁通密度各分量的傅里葉級數形式必須和感應場有相同本征值,因此在T邊界面的源磁場經過傅里葉變換后的數學表達形式為

1.5 電磁力計算驗證

將計算出的源場和各域磁矢位的通解代入邊界條件,可求解出通解中的未知數,進而根據式(1)求出磁通密度分布,采用麥克斯韋應力張量法求解懸浮系統受到的電磁力,面T、面B受到的電磁力為F1,面e、r受到的電磁力為F2。

對于面T、B,有

對于面e、r,有

所以,導體板所受的總電磁力各方向分量為

為了驗證解析模型準確性,利用Ansys Maxwell軟件建立3-D有限元仿真模型,如圖3所示,永磁電動懸浮系統的結構參數見表1。

表1 永磁電動懸浮系統基本尺寸參數Tab.1 The basical parameters of PM-EDS

圖3 3-D有限元模型Fig.3 3-D finite element model

圖 4為二維解析模型(2D_model)[1,2,8,11,23]與本文建立的三維解析模型(3D_model)電磁力計算結果的對比,以三維有限元仿真結果(3D_FEA)作為參考標準。分析圖中的計算結果可以看出,本文建立的三維解析模型相較于傳統的二維解析模型能大幅度提高電磁力的計算精度,且計算結果與有限元結果基本一致,這驗證了三維解析模型的計算準確性和可靠性。圖5為導體板橫向中心處磁感應強度沿x軸的分布情況,三維解析模型對磁場的計算結果同有限元仿真結果能較好地吻合,進一步證明了模型的準確可靠性。

圖4 不同模型電磁力計算結果對比Fig.4 The comparison of electromagnetic force calculations using different models

圖6為本文建立的三維解析模型與三維有限元模型在不同導體板寬度下的電磁力計算結果。從圖中可以看出,在不同導體板寬度下解析模型的電磁力計算結果與有限元仿真結果相差不大,驗證了模型能夠計算導體板有限寬時系統所受的電磁力,這為下一步量化分析永磁體寬度和導體板寬度對系統橫向端部效應的影響奠定了基礎。

圖5 磁感應強度在導體板橫向中心沿x軸的分布情況Fig.5 The distribution of magnetic flux density along xaxis at center of the plate width.

圖6 不同導體板寬度下三維解析模型與有限元模型電磁力計算結果對比Fig.6 The comparison of electrimagnetic force calculations using 3D_ model and 3D_FEA in case of different conductive plate widths

2 參數分析

基于三維解析模型對系統進行參數分析,研究結構參數變化對永磁電動懸浮系統性能的影響。令永磁體陣列始終位于導體板橫向寬度中心,以式(19)定義的系統浮阻比和浮重比為優化指標,改變特定參數分析優化指標受到的影響,給出最優參數取值標準。

式中,Fy和Fx分別為懸浮系統所受的懸浮力和電磁阻力,已在式(18)中計算得出;MG為Halbach永磁體陣列所受的重力,可通過永磁體陣列的體積和密度直接計算。

已有的研究表明,Halbach永磁陣列單塊永磁體的最優長度與氣隙有關,即 Halbach的最優極距長度與氣隙的選取高度有關[23],為避免工作的重復性,本文主要研究除極距外的其他結構參數變化對系統優化指標的影響。單獨研究某個參數的取值是毫無意義的,這里主要分析各參數與極距之比的最優取值標準,以保證結果的通用性。在對永磁體參數進行分析時,先將導體板寬度設定為無限寬以消除其對系統性能的影響,然后參數化分析永磁體寬度與極距之比對優化指標的影響,接著研究永磁體寬度與厚度之間是否存在耦合關系。在上述分析的基礎上研究永磁體厚度與極距之比對優化指標的影響,然后分析導體板寬度與永磁體寬度之比變化對優化指標的影響,最后研究導體板厚度對優化指標的影響。最終,給出各參數的最優取值標準。

2.1 永磁體參數分析

圖7和圖8分別為不同極距下和不同永磁體厚度下優化指標受永磁體寬度與極距之比變化的影響。

圖7表明隨著永磁體寬度的增加,系統的浮重比先增大后保持不變,浮阻比基本不受影響,極距變化對最優永磁體寬度與極距之比取值無影響,2w0/τ的最優值取為2.5。通過改變永磁體厚度掃描分析永磁體寬度與厚度對系統性能的耦合影響,圖8表明永磁體厚度變化對最優 2w0/τ取值基本無影響,認為永磁體寬度與厚度對系統性能的影響是彼此解耦的,因此對永磁體厚度優化的研究可采用二維解析模型以加快分析進程。

圖7 不同極距下優化指標受永磁體寬度與極距之比變化的影響Fig.7 The impact of the ratio of permanent magnet width to pole pitch on optimization objectives in case of different pole pitches

圖8 不同永磁體厚度下優化指標受永磁體寬度與極距之比變化的影響Fig.8 The impact of the ratio of permanent magnet width to pole pitch on optimization objectives in case of different magnet thicknesses

不同極距下浮重比受永磁體厚度與極距之比變化的影響如圖9所示。圖9表明隨著永磁體厚度的增加,系統的浮重比先增大后減小,永磁體厚度與極距之比(h0/τ)存在最優取值。將圖中不同極距下浮重比極值和該點對應的h0/τ提取出來,繪制圖10所示的最優h0/τ和極距之間的關系圖,分析可以看出,隨著極距的增長,h0/τ的最優取值逐漸趨近于0.4。當極距較小時,最優h0/τ較大,這是因為:在極距較小時,永磁體厚度與極距的比值要進一步增大才能使氣隙磁場的縱向分量達到最大值,以獲得最大浮重比。當處于浮重比極值附近時,把h0/τ的最優取值修正為0.4,畫出圖10中的修正線,可以看出修正值與實際值之間的浮重比最大差值不超過5%,完全可將修正值代替實際值,因此綜合分析最優h0/τ取為 0.4。

圖9 不同極距下浮重比受永磁體厚度與極距之比變化的影響Fig.9 The impact of the ratio of permanent magnet thickness to pole pitch on lift to weight in case of different pole pitches

圖10 極距對永磁體厚度與極距最優比值的影響Fig.10 The impact of pole pitch on the optimal ratio of permanent magnet thickness to pole pitch

2.2 導體板參數分析

永磁電動懸浮系統本質為動態感應渦流場模型,導體板的結構尺寸對渦流的分布有很大的影響,進而影響系統的電磁力。有關研究表明當磁極的橫向寬度與極距之比大于一定值時,電機的端部效應可以忽略不計[8]。圖11所示為不同永磁體寬度下浮重比受導體板寬度與永磁體寬度之比變化的影響。圖11表明隨著永磁體寬度的增加,導體板寬度對系統橫向端部效應的影響逐漸變小。但考慮 2w0/τ的最優取值為2.5,在此前提下必須考慮系統的橫向端部效應,分析可得最優導體板寬度與永磁體寬度之比wd/w0可取為1.5。

在渦流場分析中趨膚效應的存在會影響導體板厚度的選取,圖12所示為導體板厚度對電磁阻力隨速度變化的影響。圖12的二維解析模型計算結果表明在低速時隨著導體板厚度的增加,電磁阻力峰值會逐漸降為一個穩定值,高速時系統阻力受導體板厚度變化的影響不大,為減小推進損耗,在速度小于 50m/s的低速區間可將導體板厚度取為15mm,在中高速運動區間段可根據實際運行速度進一步減小厚度,降低建設成本。

圖12 導體板厚度對電磁阻力隨速度變化的影響Fig.12 The impact of the conductive plate thickness on drag force changing along vx

2.3 參數分析系統優化設計法

2.1節和2.2節參數分析首次探究了導體板寬度對系統橫向端部效應的影響,研究了結構參數對系統橫向端部效應的影響,給出的結構參數最優取值標準可有效加快系統優化的速度,減少優化時間成本。基于參數分析系統優化設計法的永磁電動懸浮系統優化步驟總結如下:根據氣隙和實際應用場景選定系統 Halbach陣列的極距,根據參數分析將永磁體的寬度取為極距的2.5倍,厚度取為極距的0.4倍。確定永磁體結構后,可將導體板寬度取為永磁體寬度的 1.5倍,對于導體板厚度,在速度小于50m/s的低速區間可將導體板的厚度取為15mm,中高速運動區間結合二維計算方法選取合適的厚度。

3 多目標優化

3.1 優化函數

利用多目標粒子群優化算法對系統的結構進行優化,優化目標為系統浮重比和浮阻比最大化,其中優化目標函數定義為

優化時,系統的運行速度vx=200m/s,g=10mm,系統結構參數組成的參數集合矢量如式(20)所示,每一個結構參數的變化范圍由第 2節參數分析確定,見表2,其他參數固定不變(見表1)。

表2 優化設計參數取值范圍Tab.2 The ranges of optimization design

3.2 多目標粒子群優化算法

多目標粒子群優化算法(MOPSO)目前已發展為比較成熟的優化算法之一,其算法中的粒子具備記憶功能,粒子之間可以相互交流自己存儲的信息,在代表全局最優解的粒子和代表局部最優解的粒子的共同引導之下迭代尋優,最終使整個粒子群不斷收斂到最優解集上。多目標粒子群優化算法的主要關系函數為[24-25]

式中,n為目前為止的迭代次數;i為種群的第i個粒子;xi為粒子的位置,即式(21)中的優化參數矢量;vi為粒子的速度,即迭代尋優時參數集合矢量的增量;Pbesti和Gbest分別為迭代過程中每個粒子的局部最優解和整個種群的全局最優解;ω為慣性權重,代表該粒子當前速度對迭代得到的該粒子的下一代位置的影響;c1為自我學習因子;c2為社會學習因子;r1和r2為0~1之間的隨機數。整個種群通過迭代尋優找出合適的Pareto最優解集,畫出Pareto前沿面軌跡,最后找出符合設計指標的最優解。

3.3 優化結果

將粒子的初始種群數量設置為 100,通過隨機數在參數取值范圍內生成100個初始化粒子,如果導體板寬度小于永磁體寬度,則將導體板寬度直接設置為永磁體寬度的 1.5倍,Pareto最優解集的大小設置為100,最大迭代次數設置為300。為了防止粒子進入局部尋優,在迭代過程中設置與迭代次數相關的突變因子,當系統生成的隨機數小于突變因子時,將該粒子的本次迭代解變異為與突變因子有關的新解。

經過優化以后,Pareto前沿面軌跡如圖13所示。選擇其中的A、B、C三個點作為優化設計解(優化指標f1和f2),并與初始設計進行對比,其參數及性能對比見附表1中MOPSO設計解。優化點A的浮重比(55.17)和浮阻比(16.51)皆大于初始設計點(49.04, 10.04),C點的浮重比(74.12)大于A點,但浮阻比(11.14)比A點小且與初始設計點接近,B(62.25,14.47)處于整個軌跡的中點平衡處,對優化設計指標的權衡達到最佳,它及其周圍的解是最佳設計目標選取區域。

圖13 多目標優化Pareto前沿面Fig.13 Pareto front of the muti-objective optimization

從附表1中可以看出MOPSO優化解的導體板橫向寬度較大,由優化算法尋找出來的設計解是指標最優解。此時保持A、B、C三點的極距不變,其他參數取值采用 2.3節參數分析系統優化法進行計算,得出附表1中的基于參數分析的系統優化設計解Ac、Bc、Cc,并計算出設計點的優化指標,可以看出后者的優化指標與前者基本接近,誤差不超過10%,分析得到誤差主要來源于參數分析系統優化法給出的結構參數計算規則,規則基礎為相關指標變化的拐點值,拐點以后的增長被忽略,但優化算法計算了拐點后優化指標的增長值,同時相應的結構尺寸也大幅增大,所以最后會產生一定的誤差。但從附表1的結果對比中可以看出,基于參數分析的系統優化設計解在對系統性能影響較小的情況下能大幅減小導體板的寬度,降低建設成本,且優化設計時間相比于MOPSO可以忽略不計,這驗證了參數分析給出的最優參數取值標準的可靠性和實用性,減少了工程設計中的計算量,節約了設計成本,最優參數取值標準見表3。

表3 最優參數取值標準Tab.3 The standard of optimization parameters calculating

4 結論

準確計算電磁力是對永磁電動懸浮系統進行深入研究和分析的基礎,本文建立了直線型永磁電動懸浮系統三維電磁解析模型,考慮了系統的橫向端部效應和縱向端部效應計算問題,基于構建的三維解析模型分析了結構參數對系統性能指標的影響,主要結論如下:

1)建立的三維電磁解析模型考慮了直線型永磁電動懸浮系統的橫向端部效應和導體板橫向邊界面磁場透射對系統性能的影響,能夠有效分析導體板寬度對懸浮系統性能的影響。解析模型適用于直線型永磁電動懸浮系統永磁體與導體板任意位置關系,更具通用性。與二維計算模型相比,大幅提高了電磁力計算精度,與3-D有限元仿真相比,將單個設計解的電磁力計算時間縮短為20s左右,大大節約了計算時間。

2)基于參數分析對直線型永磁電動懸浮系統的橫向端部效應進行了量化研究,提出一種系統優化設計方法,得到了系統結構參數最優取值標準,見表 3。通過對比參數分析系統優化設計解和MOPSO的優化解,表明在對系統性能影響較小的情況下參數分析系統優化設計法能大幅減小導體板的寬度,降低系統建設成本,且優化設計時間相較于MOPSO計算時間可以忽略不計,減少了工程設計中的計算量,節約了設計成本。

目前作者所在的課題組正有序推進桶式旋轉試驗平臺的建設,下一步將利用試驗平臺驗證理論模型的準確性。

附 錄

附表1 MOPSO設計解與參數分析系統優化設計解對比App.Tab.1 The comparison of design given by optimization based on parameter analysis and that given by MOPSO

猜你喜歡
優化模型
一半模型
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
由“形”啟“數”優化運算——以2021年解析幾何高考題為例
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
主站蜘蛛池模板: 日韩视频精品在线| 一本大道在线一本久道| 亚洲中文字幕日产无码2021| 免费观看精品视频999| 亚洲中文字幕在线精品一区| 亚洲AV电影不卡在线观看| 精品福利国产| 亚洲综合色婷婷| 久久中文电影| 茄子视频毛片免费观看| 国产成人精品免费av| 91丝袜在线观看| 国产视频a| 免费国产在线精品一区| 精品国产成人av免费| 重口调教一区二区视频| 一区二区无码在线视频| 97亚洲色综久久精品| 久久精品娱乐亚洲领先| 日本欧美在线观看| 日韩小视频在线播放| 久久久久久国产精品mv| 一级毛片无毒不卡直接观看| 欧美日韩另类在线| 国产美女在线免费观看| 天天视频在线91频| 国产成人高清亚洲一区久久| 日韩天堂在线观看| 欧美在线三级| 欧美日韩国产成人在线观看| 欧美五月婷婷| 天堂网亚洲系列亚洲系列| 亚洲无码高清视频在线观看 | 国产日韩AV高潮在线| 日本免费a视频| 91麻豆久久久| 狠狠做深爱婷婷综合一区| 国产后式a一视频| 国产亚洲精久久久久久久91| 亚洲女人在线| 热re99久久精品国99热| 亚洲天堂精品在线| 亚洲swag精品自拍一区| 欧美成人a∨视频免费观看| 手机在线看片不卡中文字幕| 国产91线观看| 中文字幕人妻无码系列第三区| 亚洲av无码成人专区| 欧美成人精品高清在线下载| 精品91在线| 在线观看亚洲精品福利片| 国产美女主播一级成人毛片| 亚洲精品成人片在线观看| 国内丰满少妇猛烈精品播| 亚洲欧美成人综合| 精品久久综合1区2区3区激情| 久久不卡国产精品无码| 不卡网亚洲无码| 国产精品亚洲片在线va| 日韩小视频网站hq| 色香蕉影院| 日本不卡视频在线| 久久成人国产精品免费软件| 福利视频一区| 国产成人精品视频一区视频二区| 九九九国产| 老熟妇喷水一区二区三区| 91成人试看福利体验区| 国产精品私拍在线爆乳| 国产精品成人不卡在线观看 | 色噜噜综合网| 国产精品欧美激情| 久热精品免费| 久久精品中文无码资源站| 美女无遮挡拍拍拍免费视频| a在线亚洲男人的天堂试看| 五月婷婷激情四射| 一区二区自拍| 毛片一级在线| 久久精品国产电影| 日本少妇又色又爽又高潮| 亚洲经典在线中文字幕|