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

基于智能優化算法和有限元法的多線圈均勻磁場優化設計

2019-05-29 07:32:34呂志峰張金生王仕成趙欣李婷
北京航空航天大學學報 2019年5期
關鍵詞:有限元法磁場智能

呂志峰, 張金生, 王仕成, 趙欣, 李婷

(火箭軍工程大學精確制導與仿真實驗室, 西安 710025)

磁場模擬裝置在地磁導航等航空航天工程中具有重要的應用價值[1-3]。理想的磁場模擬裝置應能夠在足夠大的空間內產生高均勻度的磁場,因此,磁場均勻性是衡量磁場模擬裝置性能的關鍵技術指標之一。近年來,隨著磁場模擬裝置的應用領域越來越廣泛,其磁場分布的均勻性得到了越來越多科研工作者的重視[4-6]。文獻[4]基于亥姆霍茲線圈加工設計了一套磁場模擬裝置,并對磁場分布的均勻性進行了測試,結果表明該裝置基本能夠滿足地磁導航半實物仿真系統中的磁場模擬需求。文獻[5]針對磁傳感器標定環境中均勻磁場的要求,對四線圈磁場分布的均勻性進行了仿真計算。文獻[6]針對某型號光纖陀螺性能測試的要求,設計了一套滿足磁場分布均勻性要求的大型磁場模擬裝置。

目前,均勻磁場的產生主要依靠亥姆霍茲線圈[7-8]。隨著工程應用對磁場均勻性的要求不斷提高,傳統的亥姆霍茲線圈已經很難滿足實際需求,因此,有學者提出采用共軸多組線圈組合的方式來提高磁場的均勻性[9-11]。但是,由于線圈組數增多,如果各組線圈的位置、尺寸及匝數等參數設計不合理,就會導致磁場的均勻區域大大減小,這就需要對線圈的參數進行優化使其均勻區達到最大化。對于線圈參數的優化,目前通常采用數學求導的方式求取最優參數,但是線圈組數的增多會導致待優化變量也相應增加,這就需要進行高階求導,如文獻[11]對9線圈系統的參數進行優化,進行了16階求導;文獻[12]對四環圓形線圈的參數進行優化,進行了8階求導,其計算較為復雜。因此,如何避免對函數高階求導是多線圈參數優化需要解決的一個問題。對此,文獻[13-14]提出采用差分進化算法對線圈的均勻磁場進行優化設計,取得了較好的效果,但是智能優化算法普遍存在容易陷入局部最優的缺點,而上述研究未討論尋優得到的結果是否為全局最優。另外,線圈參數的優化都基于解析式,很多學者直接以解析優化的結果作為最終參數,并根據這些參數加工制作磁場模擬裝置,但是裝置的最終性能與理論分析存在不同程度的偏差[15-16]。分析原因,解析式是在理想化條件下推導出來的,基于解析式得到的優化結果的可信度如果未進行有效評估,一旦參數優化有誤就會導致磁場的均勻性不理想,使已制作好的磁場模擬裝置的均勻性無法滿足要求,造成人力物力財力的浪費。因此,如何對優化結果的可信度進行有效評估也是線圈參數優化的重要一環。

鑒于此,本文提出一種基于智能優化算法的共軸多線圈參數優化方法,并通過Ansoft Maxwell軟件,采用有限元法對優化結果的可信度進行檢驗與評估,從而為產生高均勻度磁場的磁場模擬裝置的設計提供一種高可信度優化方法。

1 方形線圈均勻磁場的基本理論

目前,關于圓形亥姆霍茲線圈構建均勻磁場的研究較多,但是圓形線圈加工、固定等不方便,故在實際工程中多采用方形線圈用于模擬磁場。本文就以方形線圈為研究對象,對其磁場分布的均勻性進行分析研究。

1.1 方形線圈軸線上的磁場計算

首先分析單個方形線圈軸線的磁感應強度。假設方形線圈的邊長為2l,其中心為坐標原點O,水平向右為x軸的正方向,垂直于紙面向內為y軸的正方向,z軸垂直于線圈平面且向上為正方向,電流方向為逆時針,與z軸正方向滿足右手螺旋關系,線圈匝數為N,如圖1所示。P(0,0,z)為z軸上任意一點,那么根據畢奧-薩法爾定律,可以

圖1 單個方形載流線圈Fig.1 Single square current-carrying coil

求得P點處的磁場強度大小為

(1)

式中:μ0=4π×10-7N·A-2為真空磁導率;I為電流大小。

1.2 方形線圈均勻磁場的傳統優化設計方法

目前,磁場的模擬主要依靠亥姆霍茲線圈。亥姆霍茲線圈是由一對相同線圈匝數、相同纏繞方式、同軸平行放置的線圈組成。首先以方形亥姆霍茲線圈軸線上的均勻磁場區域為研究對象,采用傳統的優化方法對線圈參數進行優化。設有2個方形載流線圈1和線圈2,其邊長均為2l,匝數均為N,所通電流I同向,線圈之間相距為2a,以兩線圈中心軸線的中點作為坐標原點O,如圖2所示。

根據磁場的疊加原理及式(1)可知,兩線圈在中心軸線上任意一點Q(0,0,z)處產生的磁感應強度為

(2)

對B(z)在中心點O附近進行泰勒展開,得

(3)

從式(3)中可以看出,應盡可能使z=0處的各階導數均為0,從而使中心點附近的磁感應強度B(z)→B(0),進而使中心點附近的均勻磁場范圍盡可能大。下面對式(2)進行一階、二階求導,分別為

圖2 方形亥姆霍茲線圈示意圖Fig.2 Schematic diagram of square Helmholtz coil

(4)

(5)

顯然,由于B(z)為偶函數,那么其在z=0處的奇數階導數均為0,故dB/dz=0自動滿足。如果要使z=0處的二階導數d2B/dz2=0,則根據式(5)需要滿足:

F(a)=6a6+18a4l2+11a2-5l6=0

(6)

經計算,可得a=0.544 5l,即當方形亥姆霍茲線圈參數之間的關系滿足a=0.544 5l時,線圈軸線中心點附近的磁場均勻性最好。由式(2)可知,對于方形亥姆霍茲線圈,B(z)在z=0處的各階導數是否為0僅由間距參數a和邊長參數l之間的關系決定,在滿足d2B/dz2=0的情況下(即a=0.544 5l時),d2nB/dz2n≠0(n≥2,n∈Z),這就從根本上限制了方形亥姆霍茲線圈磁場均勻區域的大小。

鑒于此,可以采用多組線圈,通過增加參數的個數,使各參數之間滿足一定的關系,使z=0處的各階導數盡可能多為0,從而提高線圈磁場的均勻性。以兩組方形亥姆霍茲線圈為例,4個線圈的邊長均為2l,通入電流均為I,且電流方向相同。中間一組線圈之間的距離為2a1,線圈匝數均為N1,兩側一組線圈之間的距離為2a2,線圈匝數均為N2,且線圈1、2之間的距離與線圈3、4之間的距離相等,如圖3所示。

圖3 兩組方形亥姆霍茲線圈示意圖Fig.3 Schematic diagram of two sets of square Helmholtz coils

圖3中線圈在中心軸線上任意一點Q(0,0,z)處產生的磁感應強度為

(7)

對于式(7),B(z)在z=0處的各階導數是否為0由間距參數a1、a2與邊長參數l之間的關系及匝數參數N1與N2之間的關系共同決定,故需要確定3組參數關系,即k1=a1/l、k2=a2/l、k3=N1/N2,這3組參數關系的確定至少需要3個方程。由于B(z)為偶函數,那么其在z=0處的奇數階導數均為0,即dB/dz=0、d3B/dz3=0、d5B/dz5=0自動滿足。因此,需要通過d2B/dz2=0、d4B/dz4=0、d6B/dz6=0這3個方程共同確定這3組參數關系,從而使線圈在中心點附近的均勻磁場區域達到最大。式(5)是單組亥姆霍茲線圈的二階導數,其求導已經較為復雜,那么對2組亥姆霍茲線圈(即式(7))求6階導數必然會更加復雜,如果需要更大的磁場均勻區域就要采用更多組線圈,參數增多就需要求更高階導數,其計算復雜性可想而知。不難看出,在多組線圈均勻磁場優化設計過程中,如果仍采用傳統的求導優化方法,勢必會導致問題求解復雜化,因此,需要尋求一種新的解決方法以避免高階求導。

另外,實際的磁場模擬裝置都是由多匝線圈纏繞而成,故會形成一定的纏繞寬度和纏繞厚度,而傳統的解析理論分析將線圈抽象為一根沒有寬度和厚度的導線,即式(2)和式(7)是一種理想情況下的解析式,這就導致理論分析與工程實際存在一定的偏差。然而很多學者直接以解析式優化的結果作為最終參數,并根據這些參數加工制作磁場模擬裝置,那么可能會導致最終的實際系統性能與理論分析出現較大偏差,因此還需要對解析式的優化結果的正確性進行檢驗與評估。

2 多線圈均勻磁場優化設計方法

2.1 智能優化算法基本理論

經典優化理論的局限性在于過度依賴優化問題的數學特征,但是,隨著問題規模的擴大及問題種類的增多,實際優化問題的數學特征越來越難以得到。在數學知識尚未能提供理論支撐的情況下,針對復雜優化系統出現的大規模、強非線性、不可微等特征,經典優化理論的局限性越來越明顯。因此,學者們在觀察現實世界中存在的自發性優化現象的基礎上,進行數學抽象以模擬該現象中的智能行為,提出了一類無需依賴優化問題數學特征的優化算法,這類算法統稱為智能優化算法。目前,常用的智能優化算法有遺傳算法、粒子群優化算法、蟻群算法、人工免疫算法、模擬退火算法、禁忌搜索等[17],這些算法雖然模擬的個體或群體現象不同,但是其共同之處都是利用它們的智能特性和行為方式來解決實際優化問題。智能優化算法為解決復雜優化問題提供了一種新的解決途徑。后續仿真實驗中,本文將采用粒子群優化算法進行尋優。

2.2 磁場的有限元分析方法

麥克斯韋方程組是電磁場計算的基本出發點及核心。目前,磁場的計算主要有解析法和數值法。這2種方法的理論基礎都是麥克斯韋方程組,但由于數學理論的局限性,對于復雜情況的磁場計算,需要通過一些近似處理才能得到其解析式,故解析法的計算精度相對較低。隨著計算機硬件的迅猛發展,數值法被廣泛應用于磁場的計算,其中最具代表性的就是有限元法。

有限元法是以變分原理為基礎逐步發展起來的,其基本思想是把待求解場域剖分成若干個小單元,然后利用一定的線性算法來求解各單元,而求解這些小單元是很簡單的,最后將所有單元的求解結果進行綜合便可得到整個場域的數值解[18]。從有限元法的基本思想可以看出,如果剖分的單元足夠小,計算結果就會無限趨近于理論解,故有限元法的求解精度要高于解析法。從現有文獻也可以看出,與解析法相比,有限元法的計算結果與實際系統的吻合度也較高[19-20],正因如此,有限元法常常被用于校驗解析計算結果的準確性。目前,世界上有多款電磁場有限元分析軟件,如美國的Ansoft Maxwell、法國的Flux、加拿大的Infolytica和日本的JMAG等,其中,Ansoft Maxwell軟件具有功能強大、計算結果精確、易于使用等優點,本文將采用這款軟件用于磁場的有限元計算。

2.3 線圈均勻磁場優化設計流程

在實際情況中,對線圈的均勻磁場進行優化并不是要求磁場的均勻區域無限大,而是在滿足我們的工程技術指標的前提下,使磁場的均勻區域盡可能最大化。

由 1.2 節可知,對于多組線圈結構參數的傳統優化設計,其難點在于高階求導,且基于解析法得到的優化結果的可信度有待確認,為此,本文提出一種基于智能優化算法和有限元法相結合的高可信度參數優化方法。首先,確定線圈的組數和待優化參數,同時,確定磁場均勻度指標,并以此作為待優化的目標函數。然后,采用智能優化算法對目標函數進行尋優,由于智能優化算法在尋優過程中具有隨機性并容易陷入局部極值,可以采用蒙特卡羅方法進行多次智能尋優。若大部分優化結果較為集中,說明大部分尋優收斂到了全局最優附近,就以眾多優化結果中的最優值作為最優參數;若大部分優化結果分散,說明尋優可能陷入了局部最優,即使是眾多優化結果中的最優值,也無法保證其是全局最優,造成這種結果可能是因為種群規模太小,這時需要增加種群的數量,進而增加了潛在解的多樣性,進行重新尋優,直到得到最優參數。由于智能優化算法尋優得到的線圈匝數比很可能不是整數比,而線圈的實際加工是以整數匝纏繞的,將匝數比近似成整數比后,在該匝數比的情況下,原來的間距比很有可能不再是最優間距比,因此,需要在該匝數比固定的情況下,進行重新尋優以調整間距比。得到最優參數后,需要判斷該參數下磁場的均勻性指標是否滿足工程要求,如果不滿足要求,說明線圈組數不夠多,需要增加線圈組數進行重新優化。最后,在智能尋優結果滿足均勻性指標的情況下,采用有限元法對優化結果的可信度進行檢驗與評估。由于有限元仿真更加接近真實條件,故如果有限元計算結果不滿足均勻性指標,說明按照智能尋優得到的最優參數加工成型的線圈很大程度上難以滿足均勻性指標,因此仍然需要增加線圈組數重新優化,反之,即可確認智能優化算法尋優得到的參數為最優參數。整個優化流程如圖4 所示。

圖4 參數優化流程圖Fig.4 Flowchart of parameter optimization

3 仿真分析

3.1 有效性驗證

為了說明智能優化算法尋優與求導尋優的結果具有一致性,本文以單組亥姆霍茲線圈(如圖2所示)的參數優化為例,采用2種不同的方法進行優化,將2種方法得到的結果進行對比。

首先,確定待優化參數及目標函數,由1.2節式(2)可知,待優化參數為k=a/l,定義磁場偏差率f=|B(z)-B(0)|/B(0),并以此作為目標函數,即以中心點附近的磁感應強度B(z)與中心點處的磁感應強度B(0)之間的相對偏差為目標函數,當二者之間的偏差最小時,線圈產生的磁場最均勻。然后,采用智能優化算法對目標函數進行尋優,這里,智能優化算法選用粒子群優化算法,該算法的基本思想是模擬鳥類的群體行為,具有群體智能、并行計算、快速收斂等優點[17],在MATLAB環境下編程運行,求得最優參數為k=a/l=0.546 4。

由1.2節可知,求導尋優得到的最優參數為k=a/l=0.544 5,經計算,智能優化算法得到的最優參數與求導尋優得到的最優參數的相對偏差僅為0.35%。這里需要說明的是,智能優化算法是一種近似算法,其在有限的時間開銷、資源消耗等成本約束下求得的是近似最優解,因此,2種結果出現一定的偏差是必然的,但是二者偏差僅為0.35%,可以看做2種方法的求解結果具有一致性。假設線圈的的半邊長l=0.5 m,那么在2種參數情況下,線圈中心軸線1 m長度范圍內,以磁場偏差率f=|B(z)-B(0)|/B(0)作為磁場均勻性評價指標進行定量計算,計算結果如表1所示。

從表1可以看出,本文方法和求導方法得到的參數下,相同范圍內的線圈的磁場偏差率保持在同一數量級且基本相近,這說明本文提出的方法是有效的,能夠代替傳統的求導尋優。

表1 單組線圈中心軸線磁場偏差率Table 1 Magnetic field deviation rate of single set of coils along central axis

3.2 優勢性驗證

本文以2組亥姆霍茲線圈(如圖3所示)的參數優化為例,通過與傳統的求導方法進行對比說明本文所提方法具有一定的優勢性。這里,假設均勻性指標為:沿著軸線方向,以軸線中點為中心的0.6倍線圈長度范圍內,磁場偏差率小于0.05%。

3.2.1 兩組亥姆霍茲線圈均勻磁場優化設計

首先采用本文提出的方法進行參數優化。以磁場偏差率f=|B(z)-B(0)|/B(0)為目標函數,采用粒子群優化算法對如圖3所示的2組亥姆霍茲線圈的結構參數進行優化,得到最優參數為k1=a1/l=0.293 4,k2=a2/l=1.101 4,k3=N1/N2=0.465 9。由于線圈的匝數只能為整數,為了盡可能接近k3=N1/N2=0.465 9,將匝數比設置為整數比k3=N1/N2=7/15≈0.466 7。固定該匝數比,按照圖4所示流程圖,對間距參數k1和k2進行再次尋優,得到最終的最優參數為k1=a1/l=0.294 2,k2=a2/l=1.103 0,k3=N1/N2=7/15。

最后,在結構參數固定的情況下,調整電流的大小,使2種情況下中心點處產生的磁感應強度大小一致,從而更直觀地體現2種情況下磁場均勻性的好壞。2種情況的參數設置如表2所示。

按照表2的結構參數,將其代入式(7),通過MATLAB軟件計算得到2種參數情況下線圈中心軸線上的磁場如圖5所示。

從圖5中可以定性看出,本文方法得到的磁場均勻性明顯優于傳統求導方法的磁場均勻性。為了定量說明本文方法具有優勢,在圖5所示的線圈中心軸線1 m長度范圍內,以磁場偏差率f=|B(z)-B(0)|/B(0)作為磁場均勻性評價指標進行定量計算,計算結果如表3所示。

圖5 兩種參數情況磁場分布均勻性對比Fig.5 Comparison of magnetic field distribution uniformity between two different parameters

軸線長度/m磁場偏差率本文方法傳統求導方法0.21.8471×10-51.2616×10-40.46.5529×10-56.2984×10-40.62.6851×10-40.00390.80.00640.020710.03560.0702

從表3可以定量看出,求導尋優得到的參數在0.6 m范圍內無法滿足磁場偏差率小于0.05%的均勻性指標,而智能優化算法尋優得到的參數能夠滿足。在軸線長度0.2、0.4和0.6 m的范圍內,傳統求導方法與本文方法的磁場偏差率不在同一數量級,二者的比值分別為6.83、9.61、14.52,采用本文方法求得的結構參數明顯優于傳統的求導尋優結果,且由中點沿軸線向兩側延伸,本文方法的優勢越發明顯。

3.2.2 最優參數的可信度檢驗與評估

下面采用有限元法對智能優化算法得到的最優參數的可信度進行檢驗與評估,同時,為了進一步說明智能尋優優于傳統的求導尋優,本文在這部分對求導尋優得到的參數也進行有限元仿真,比較二者在更真實的仿真情況下的磁場分布。

按照表2所示的結構參數,在Ansoft Maxwell軟件中建立如圖6所示的有限元數值仿真模型,其中,為了更符合工程實際情況,線圈的纏繞寬度設定為10 mm,纏繞厚度設定為4 mm,方形線圈4個直角處均以半徑為50 mm、R角為90°的圓弧平滑過渡。通過有限元計算,得到本文方法參數和求導方法參數下xOy平面和xOz平面的磁場分布情況,如圖7所示。

從圖7可以看出,在xOy平面和xOz平面中心區域的很大范圍內,無論是本文方法還是傳統的求導優化方法,磁場都呈現高度均勻分布,為了更直觀地比較2種方法磁場分布的均勻性,取z軸1 m長度范圍內的磁場進行考察,將Ansoft Maxwell軟件中的仿真結果導出,在MATLAB中作圖,結果如圖8所示。

對比圖8和圖5可以看出,有限元法求得的中心點處的磁感應強度大小與解析式求得的結果存在一定的偏差:在智能優化參數下,二者的相對偏差為5.51%;在求導優化參數下,二者的相對偏差為4.55%。這說明純解析計算與接近實際情況的有限元計算存在一定的偏差,通過優化算法得到的結果的可信度必須經過檢驗與評估。

圖7 不同平面磁場分布Fig.7 Magnetic field distribution in different planes

這里,主要關注的是磁場的均勻性指標。從圖8可以看出,單就磁場的均勻性而言,有限元法計算得到的軸線中心點附近的磁場依然具有良好的均勻性。根據圖8中的仿真數據,將磁場均勻性進行量化,計算磁場偏差率,將本文方法參數和傳統求導方法參數下的有限元計算結果進行比較,如表4所示。

從表4可以看出,在考慮了線圈纏繞寬度、纏繞厚度及R角的情況下,智能優化算法參數下的有限元仿真結果依然滿足磁場偏差率小于0.05%的均勻性指標,因此,該結構參數可以得到確認。另外,在軸線長度0.2 m和0.4 m的范圍內,傳統求導方法與智能優化算法的磁場偏差率相差很小,可以看做基本一致,但是在軸線長度0.6 m的范圍內,二者的比值達到了4.29,這與3.2.1節中仿真結論“由中點沿軸線向兩側延伸,本文方法的優勢越發明顯”相一致,說明智能優化算法優于傳統的求導尋優。

圖8 z軸磁場分布Fig.8 Magnetic field distribution in z axis

軸線長度/m磁場偏差率本文方法傳統求導方法0.28.3268×10-59.2408×10-50.43.2363×10-43.7401×10-40.63.9642×10-40.00170.80.00560.017410.03590.0678

4 結 論

針對多線圈均勻磁場優化設計中的高階求導及優化結果可信度評估問題,本文提出一種基于智能優化算法和有限元法的多線圈均勻磁場優化設計方法,主要得到以下結論:

1) 較傳統的亥姆霍茲線圈,多線圈系統能夠擴大磁場的均勻區域,但是其參數優化設計如果仍采用傳統的求導優化方法,會增加問題求解的復雜度。

2) 在多線圈系統的參數優化問題中,本文方法要優于傳統的求導尋優,對于2組方形亥姆霍茲線圈,在其中心軸線長度為線圈邊長的0.2倍、0.4倍和0.6倍的范圍內,傳統求導方法與本文方法的磁場偏差率不在同一數量級,二者的比值分別為6.83、9.61、14.52,說明由中點沿軸線向兩側延伸,本文方法的優勢越發明顯。

3) 通過引入有限元法進行驗證,一方面能夠對本文方法求得的參數的可信度進行檢驗與評估,另一方面也驗證了本文方法優于傳統的求導尋優。

猜你喜歡
有限元法磁場智能
西安的“磁場”
當代陜西(2022年6期)2022-04-19 12:11:54
為什么地球有磁場呢
正交各向異性材料裂紋疲勞擴展的擴展有限元法研究
智能前沿
文苑(2018年23期)2018-12-14 01:06:06
智能前沿
文苑(2018年19期)2018-11-09 01:30:14
智能前沿
文苑(2018年17期)2018-11-09 01:29:26
智能前沿
文苑(2018年21期)2018-11-09 01:22:32
磁場的性質和描述檢測題
2016年春季性感磁場
Coco薇(2016年1期)2016-01-11 16:53:24
三維有限元法在口腔正畸生物力學研究中發揮的作用
主站蜘蛛池模板: 美女被操91视频| 国产高清无码第一十页在线观看| 内射人妻无码色AV天堂| 露脸真实国语乱在线观看| 九九热精品视频在线| 国产午夜一级毛片| 香蕉eeww99国产在线观看| 深爱婷婷激情网| 一本色道久久88| 国产91线观看| 高清不卡一区二区三区香蕉| 欧美亚洲国产精品久久蜜芽| 国产精品吹潮在线观看中文| 欧美激情视频一区| 99久久99视频| 中文字幕va| 二级毛片免费观看全程| 亚洲色欲色欲www网| 四虎永久免费地址| 激情综合五月网| 免费国产好深啊好涨好硬视频| 成人91在线| 在线播放真实国产乱子伦| 亚洲不卡影院| 亚洲成人网在线观看| 日本亚洲成高清一区二区三区| 国产一在线观看| 无码专区在线观看| 欧美日韩一区二区在线播放| 欧美翘臀一区二区三区| 亚洲视频影院| 亚洲精品无码久久毛片波多野吉| 亚洲侵犯无码网址在线观看| 欧美在线网| 亚洲综合经典在线一区二区| 精品视频一区二区观看| 97免费在线观看视频| a级毛片在线免费| 丝袜无码一区二区三区| 亚洲第一成年网| 99精品视频九九精品| 亚洲成肉网| 免费看a毛片| 91久久青青草原精品国产| 国内精品手机在线观看视频| 国产在线小视频| 亚洲欧美不卡视频| 久久亚洲精少妇毛片午夜无码 | 色综合色国产热无码一| 国语少妇高潮| 天天做天天爱夜夜爽毛片毛片| 色综合天天娱乐综合网| 国内精品久久九九国产精品 | 国产专区综合另类日韩一区| 国产精品一区二区无码免费看片| 澳门av无码| 极品av一区二区| 青青青国产精品国产精品美女| 99九九成人免费视频精品| 亚洲系列无码专区偷窥无码| 激情六月丁香婷婷| 91av成人日本不卡三区| 国产精品丝袜视频| 一级毛片免费观看不卡视频| 免费中文字幕在在线不卡| 美女被躁出白浆视频播放| 波多野结衣一二三| 最新精品国偷自产在线| 91在线高清视频| 国产欧美日韩综合在线第一| 性欧美精品xxxx| 欧美视频在线播放观看免费福利资源 | 高清欧美性猛交XXXX黑人猛交| 国产精品第一区在线观看| 超碰免费91| 婷婷成人综合| 国产97色在线| 国产色爱av资源综合区| 国产欧美日韩精品综合在线| 日本精品影院| 播五月综合| 久久夜夜视频|