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

計及遮雨帽和隔音裝置的空心電抗器散熱性能分析及結構參數優化設計

2022-11-01 10:39:20陳煒胡勝盧鈴唐奇曹浩
南方電網技術 2022年9期

陳煒,胡勝,盧鈴,唐奇,曹浩

(1. 國網湖南省電力有限公司電力科學研究院,長沙 410007;2. 國網電力設施噪聲與振動實驗室,長沙 410007)

0 引言

干式空心電抗器因具備結構簡單、線性度高和重量輕等優點,在電力系統中得到了廣泛應用[1 - 2]。然而,空心電抗器在實際運行過程中易發生過熱、燒傷甚至起火燒毀等故障,研究表明包封線圈局部溫度過高是主要原因[3]。目前,為降低環境因素(特別是雨水)對電抗器的影響,常在電抗器端部加裝遮雨帽[4 - 5];為降低電抗器周圍噪音,常在電抗器周圍安裝隔音裝置。然而,電抗器線圈產生的熱量使周圍流體向上流動,遮雨帽和隔音裝置阻礙了流體的流動,導致電抗器包封線圈的散熱能力降低,溫升顯著增加。因此,為了提高電抗器的散熱能力,在準確計算空心電抗器溫升的基礎上開展遮雨帽和隔音裝置結構參數的優化研究至關重要。

在電抗器溫升計算方面,文獻[6 - 7]介紹了變壓器繞組的平均溫升計算方法,該方法精度較低。文獻[8]將電抗器包封線圈的散熱過程等效為豎直管道,給出了溫升解析計算方法,該方法僅適用于氣道兩側包封壁面的熱流密度。文獻[9]由傳熱學準則式推導出電抗器包封壁面沿軸向的對流傳熱系數表達式,進而計算出電抗器的繞組溫升分布,然而該方法無法獲得電抗器詳細的溫度場和流場分布。在此基礎上,文獻[10]推導出在加遮雨帽工況下的溫升解析式,但該方法受制于包封線圈壁面對流傳熱系數,計算精度較低。相比于平均溫升和解析計算方法,采用有限元法能夠獲得詳細的溫度場和流場分布,被廣泛應用于電抗器/變壓器等電工裝備的溫度場計算中[11 - 12]。文獻[12 - 15]建立了空心電抗器溫度場仿真模型,采用包封線圈溫度與周圍流體耦合的方法獲得了電抗器詳細的溫度分布。然而,上述方法未考慮遮雨帽和隔音裝置對電抗器溫度場分布和散熱特性的影響。

在電抗器熱優化方面:文獻[17 - 18]提出了通過調整包封線圈高度和氣道寬度等方法以提高線圈-氣道單元散熱能力,該方法適用于未加遮雨帽和隔音裝置的工況。在此基礎上,相關文獻開展了關于空心電抗器遮雨帽和隔音裝置的優化研究,考慮到遮雨帽對電抗器溫升的影響,文獻[19]將遮雨帽等效為傾斜擋板,通過調整包封線圈參數以提高線圈散熱能力。文獻[3]采用有限元法分析了在強制風冷條件下遮雨帽結構對氣道內流體流速的影響規律,通過調整遮雨帽結構參數能夠實現包封線圈間氣道內流體流速基本相同。考慮到隔音裝置對電抗器溫升的影響,文獻[14]建立了考慮隔音裝置下空心電抗器溫度場仿真模型,分析了隔音裝置結構參數對電抗器溫升的影響規律,獲得了最佳的設計參數,能夠在一定程度上降低電抗器的溫升。然而,遮雨帽和隔音裝置結構參數對電抗器散熱能力的影響復雜,上述方法僅從局部開展遮雨帽或隔音裝置的優化設計,無法實現電抗器整體散熱性能的優化,限制了其實際應用。

本文建立空心電抗器流場-溫度場三維仿真模型,分析了加遮雨帽和隔音裝置前后電抗器溫度場和流場分布特點,給出了加遮雨帽和隔音裝置工況下電抗器溫升顯著增加的原因。在此基礎上,采用拉丁方試驗設計和有限元仿真計算相結合的方法,構建了電抗器包封線圈最高溫度與遮雨帽、隔音裝置結構參數間的Kriging模型,獲得了各參數對電抗器最高溫度的影響規律。采用粒子群優化算法獲得最佳的遮雨帽和隔音裝置結構參數,結果表明優化方法能夠顯著降低空心電抗器的溫升。

1 空心電抗器基本結構及等效模型

1.1 基本結構及參數

本文以額定電壓為220 kV,電流為1 800 A,電感為21 mH的干式空心電抗器作為研究對象。空心電抗器的本體結構是由多個同軸的包封線圈組成,各包封線圈在電氣上并聯連接,相鄰包封間采用氣道隔開,氣道內由引拔條撐起,起到絕緣和散熱的作用。線圈由多根并聯圓導線或扁導線組成,每根導線上包有聚酯薄膜作為匝絕緣,包封由內向外依次為導體、絕緣材料和環氧玻璃纖維,包封上下端由星形架連接,起到加固和分配電流的作用。基于上述的技術參數和結構型式,設計的電抗器主要參數如下:線圈高度為1.5 m,氣道寬度為0.025 m,線圈內外半徑分別為0.3 m和1.0 m,包封數量為12。電抗器線圈頂端和底端安裝星形架,其長度、寬度和厚度分別為2.0 m、0.1 m和0.1 m,沿圓周呈45 °排列,結構如圖1所示。

圖1 空心電抗器基本結構Fig.1 Basic structure of air core reactor

1.2 等效模型

根據空心電抗器結構參數,利用COMSOL仿真軟件建立了三維模型,考慮到計算的準確性及計算時長,仿真模型做了如下簡化和等效:1) 模型中僅考慮包封線圈穩態散熱過程,忽略包封間撐條和絕緣支柱對電抗器溫升的影響;2) 為了兼顧計算精度和時長,整個計算域為正方體,邊長為6 m,空心電抗器三維等效模型如圖2所示。

圖2 空心電抗器三維等效模型 Fig.2 Three dimensional equivalent model of air core reactor

2 電抗器流場-溫度場仿真計算

2.1 流場-溫度場仿真計算

熱源計算、控制方程、邊界條件設置和網格剖分是準確獲得電抗器溫度場的關鍵步驟。

1)熱源計算

空心電抗器的損耗主要由包封線圈損耗和星形架損耗構成,線圈損耗由電阻損耗和渦流損耗組成。根據電抗器包封線圈電流和周圍磁場分布,可以得到電抗器各包封線圈的總損耗。基于星形架的結構尺寸和周圍磁場分布可以得到星形架的損耗。將上述計算得到的損耗施加在仿真模型中,作為電抗器溫度場仿真模型的熱源。

2)控制方程

空心電抗器熱源主要來自包封線圈和星形架的損耗,線圈和星形架產生的熱量主要通過熱傳導、熱對流和熱輻射的方式向外傳遞。在電抗器線圈和星形架固體區域,熱量主要以熱傳導的方式從高溫區域向低溫區域傳遞;考慮到電抗器內部相鄰包封線圈的溫升差別不大,在電抗器內部包封線圈表面,熱量主要以熱對流的方式與周圍流體進行換熱;在最內包封線圈內表面、最外包封線圈外表面和星形架表面,主要通過熱對流和熱輻射進行散熱。其控制方程可參考文獻[20 - 21]。

3)邊界條件設置

根據圖2的空心電抗器三維模型,溫度場仿真的邊界條件按如下方式進行設置:溫度場采用層流和流體傳熱2個模塊,電抗器包封線圈和星形架為設置為靜止壁面,各方向速度均為0;整個模型下底面設置為入口,前后側面、左右側面和上表面設置為出口,法向速度為0,表面溫度設置為環境溫度;考慮到最內最外包封線圈和星形架表面的熱輻射過程,設置其表面發射率為0.9。其中環境溫度設置為20 ℃。

4)網格剖分

網格的疏密程度直接影響到溫度場仿真計算的準確性,為了兼顧計算精度和計算速度,采用自由剖分網格,在靠近包封線圈和星形架位置處網格剖分密集,在遠離包封線圈位置網格剖分稀疏。網格剖分結果如圖3所示。

圖3 網格剖分結果 Fig.3 Results of mesh generation

2.2 溫度場仿真結果

根據上述的計算方法,仿真總時長設置為10 h,步長設置為0.1 h,容差設置為0.05。按照上述方法,計算得到的電抗器溫度場仿真結果如圖4所示。

由圖4可知,電抗器最高溫度為73.6 ℃,環境溫度為20 ℃,最高溫升為53.6 ℃。根據上述的仿真結果,提取了第5、6和7包封線圈的溫度分布結果,如圖5所示。

圖5 部分包封線圈溫度分布Fig.5 Temperature distributions of partially encapsulated coils

由圖5可知,電抗器包封線圈溫升分布規律基本相同,沿軸向方向隨高度增加呈逐漸上升的趨勢,最高溫度位于包封線圈頂端位置。其中,電抗器各包封線圈的最高溫度結果如表1所示。由表可知,電抗器內部各包封線圈最高溫度基本相同。

表1 各包封線圈溫度仿真結果Tab.1 Simulation results of temperature of each encapsulated coil

3 加裝遮雨帽和隔音裝置的電抗器溫度場仿真結果與散熱特性分析

3.1 溫度場仿真結果

考慮到遮雨帽和隔音裝置會阻礙氣道內的流體流動,包封線圈散熱條件惡劣,因此需分析加遮雨帽和隔音裝置下電抗器溫升分布特點,獲得其散熱特性。其中:電抗器頂端安裝遮雨帽,通過絕緣支柱支撐;電抗器最外包封線圈周圍安裝隔音裝置,呈圓筒結構,上部和下部中心開孔。遮雨帽結構參數:半徑為1.2 m,高為0.3 m,遮雨帽底端與隔音罩頂端的距離為0.2 m;隔音裝置參數:整體高度為1.8 m,上下中心孔半徑分別為0.5 m和0.5 m,厚度為0.1 m,材料為熱塑性聚酯薄膜。根據上述參數,建立了涵蓋遮雨帽和隔音裝置的電抗器三維模型,如圖6所示。

圖6 加遮雨帽和隔音裝置的電抗器等效模型Fig.6 Equivalent model of reactor with rain cover and sound arrester

基于電抗器流場-溫度場仿真計算方法給出了加遮雨帽和隔音裝置下電抗器溫度場仿真結果,如圖7所示。

圖7 加遮雨帽和隔音裝置下溫度場仿真結果Fig.7 Simulation result of temperature fields with rain cover and sound arrester

由圖7可知,加遮雨帽和隔音裝置后電抗器的最高溫度為93.6 ℃,最高溫升為73.6 ℃。對比圖4的仿真結果可知,加遮雨帽和隔音裝置后電抗器溫升增加20.0 ℃,溫升顯著增加。

3.2 包封線圈散熱特性分析

為了分析遮雨帽和隔音裝置對電抗器包封線圈散熱特性的影響規律,提取了沿不同路徑下的溫度場和流場分布,其中提取的路徑如圖8所示。

圖8 路徑的選取Fig.8 Selection of the paths

1)頂端位置溫度分布

根據加/未加遮雨帽和隔音裝置兩種工況下電抗器的溫度場仿真結果,繪制得到在頂端位置處各包封溫度分布曲線,如圖9所示。

圖9 加遮雨帽和隔音裝置前后電抗器溫度分布Fig.9 Results of temperature fields with and without the rain cover and sound arrester

由圖9可知,未加遮雨帽和隔音裝置時,電抗器內部各包封線圈最高溫升基本相同,最內最外包封線圈溫度明顯低于內部包封線圈,主要由于最內最外包封線圈一側向氣道對流換熱,另一側與周圍大空間進行對流換熱和熱輻射,相比于內部線圈具有更好的散熱條件。加裝遮雨帽和隔音裝置后,電抗器各包封線圈溫升均有所增加,增加幅度為8~23 ℃,熱點溫升位于第10個包封線圈。

2)包封中軸線溫度分布

為了分析電抗器包封線圈內部溫度分布,提取包封線圈5、6和7中軸線位置處的溫度,如圖10所示。

由圖10可知,加和未加遮雨帽和隔音裝置時各包封線圈的溫升分布趨勢基本相同,隨著高度的增加呈現逐步增加的趨勢,在包封線圈頂端位置處略微下降,主要由于在線圈端部位置對流散熱充分。加裝遮雨帽和隔音裝置后,包封線圈溫升顯著增加,位于線圈中上部位置。

圖10 包封線圈中軸線溫度分布Fig.10 Temperature distributions along the central axis of the encapsulation coil

3)氣道中軸線溫度分布

為了分析電抗器包封線圈氣道內流體的溫度分布,提取了包封線圈5-6和6-7間的氣道內中軸線流體溫度分布,如圖11所示。

圖11 氣道中軸線溫度分布Fig.11 Temperature distributions on central axis of airway

由圖11可知,加和未加遮雨帽和隔音裝置時,氣道中軸線位置處溫度分布趨勢基本相同,沿軸向方向隨高度增加逐漸增加的趨勢;且沿著軸向方向的高度增加至1.2 m以后,溫度呈急劇上升的趨勢。

為分析加遮雨帽和隔音裝置后對電抗器包封線圈間流體流速的影響,提取了包封線圈5-6和6-7間氣道中軸線位置處軸向流速分布,如圖12所示。

圖12 氣道中軸線軸向流速分布Fig.12 Distributions of fluid velocity on central axis of airway

由圖12可知,加/未加遮雨帽和隔音裝置時,包封線圈間氣道內流體流速沿軸向方向流速分布規律基本相同,在高度的20%和80%區間范圍內流速分布基本相同,由于路徑選取為星形架下端,在兩端位置速度為0。同時,加遮雨帽和隔音裝置后,氣道流速分布明顯降低。

提取加/未加遮雨帽和隔音裝置下氣道內流體最高流速分布,如圖13所示。

圖13 中間位置氣道內流體最高流速分布Fig.13 Distributions of the highest fluid velocity of airway in the middle position

由圖13可知,未加遮雨帽和隔音裝置時,各包封線圈間的流速分布及平均流速基本相同,最高速度約為1.2 m/s;加遮雨帽和隔音裝置后,氣道內整體流速都大幅度下降且分布不均,下降的幅度在0.3~0.6 m/s之間,且最高流速為0.8 m/s,明顯低于未加遮雨帽和隔音裝置的工況下,主要是由于遮雨帽與隔音裝置阻礙了氣道內流體的流速。同時,由于遮雨帽和隔音裝置對靠近外側的影響較大,靠近外側包封線圈的溫度較高。

綜上所述,通過對比分析可知,在加遮雨帽和隔音裝置下電抗器內部包封線圈最高溫升不相同,其主要原因為由于隔音裝置對流速的影響程度迥異,內部包封線圈溫升顯著增加。

4 遮雨帽和隔音裝置結構參數優化設計

根據上述的分析可知,遮雨帽和隔音裝置對電抗器包封線圈間氣道內流體流動具有明顯的阻礙作用,造成各包封線圈溫升顯著增加。考慮到影響電抗器溫升的遮雨帽和隔音裝置結構參數眾多,因此,有必要開展遮雨帽和隔音裝置結構參數的優化研究,以提高電抗器線圈的散熱能力。

4.1 試驗設計

根據加遮雨帽和隔音裝置的電抗器等效模型,結合實際工程經驗,確定影響電抗器溫升的遮雨帽和隔音裝置主要結構參數,其中影響因素總共8個,分別為遮雨帽半徑x1、遮雨帽高度x2、隔音裝置上端到遮雨帽底端距離x3、隔音裝置頂端中心孔半徑x4、隔音裝置頂端與線圈頂端距離x5、最外包封線圈與隔音裝置外側的距離x6、隔音裝置底端到線圈底端距離x7,隔音裝置底端中心孔半徑x8,如圖14所示。

圖14 遮雨帽和隔音裝置優化參數Fig.14 Optimum parameters of the rain cover and sound arrester

考慮到電抗器實際的絕緣要求、散熱和隔音效果,各參數選取的范圍如下:x1的范圍為1.0~1.3 m,x2的范圍為0.1~0.3 m,x3的范圍為0.1~0.3 m,x4的范圍為0.3~0.7 m,x5的范圍為0.1~0.3 m,x6的范圍為0.2~0.4 m,x7的范圍為0.1~0.3 m,x8的范圍為0.3~0.7 m。

采用拉丁方試驗設計方法[22],可以獲得試驗設計表,總數為50組。結合電抗器溫度場仿真計算方法,可以得到不同遮雨帽和隔音裝置結構參數下的溫度場仿真結果,如表2所示。

表2 拉丁方試驗設計及仿真結果Tab.2 Latin square experimental design and the simulation results

由表2可知,電抗器最高和最低溫度分別為109.9 ℃和81.5 ℃,因此,遮雨帽和隔音裝置結構參數對電抗器的溫升影響顯著。根據上述的仿真結果,給出最高和最低溫度的溫度場仿真結果,如圖15—16所示。

圖15 最高溫度仿真結果(樣本38)Fig.15 Results of temperature distribution (case 38)

圖16 最低溫度仿真結果(樣本7)Fig.16 Results of temperature distribution (case 7)

4.2 Kriging模型構建

Kriging方法即空間局部插值法,可根據已知數據點對未知數據點進行無偏最優估計[23]。未測點的估計值由相鄰已測點加權求和求得,如式(1)所示。

(1)

式中:x0為估值點;x1,x2,…,xm為已知數據點,所得樣本值相應為T(x1),T(x2),…,T(xm);未測點的估值為T(x0);μi為權重系數。

本文根據表2仿真結果建立了電抗器最高溫升與結構參數間的Kriging近似模型,以反映電抗器最高溫度與遮雨帽和隔音裝置結構參數的響應關系,如圖17所示。

圖17 Kriging近似模型Fig.17 Kriging approximation model

由圖17可知,隨著參數x7、x8的增加電抗器的最高溫度呈現逐漸下降的趨勢,且參數x7下降的趨勢相對于參數x8較為平緩;隨著參數x4、x6的增加,電抗器的最高溫度呈現先增加后下降的趨勢。

4.3 靈敏度分析

為了分析遮雨帽和隔音裝置結構參數對電抗器溫度的影響規律,本文采用靈敏度分析技術,分析各因素對電抗器溫度影響程度,其中,靈敏度指標定義為:

(2)

式中:Xi為設計變量;Y為各子系統狀態變量,Var(E(Y|Xi))為E(Y|Xi)的無條件方差;Var(Y)為Y的無條件方差。

根據式(2),可以獲得遮雨帽和隔音裝置結構參數對電抗器最高溫度的靈敏值,如圖18所示。

圖18 各優化參數對電抗器最高溫度的靈敏值Fig.18 Sensitivities of optimized parameters to the maximum temperature of reactor

根據靈敏度分析結果,將設計變量劃分為兩個層次,將敏感性指數Si>0.1劃分為敏感設計變量,并將Si≤0.1劃分為不敏感設計變量。因此參數x3、x4、x7和x8對電抗器最高溫度的影響較為顯著,其他因素影響較小。

根據近似模型數據,可以得到各參數對電抗器最高溫度的影響曲線,如圖19所示。

由圖19可知,參數x1、x3、x5、x6和x7對電抗器最高溫度的影響趨勢基本相同,隨著參數的增加最高溫度呈先增加后下降的趨勢。其中,參數x3對溫度的影響最大,參數x6對溫度的最小。此外,隨著參數x2、x4和x8的增加最高溫度大體呈現逐漸降低的趨勢,且參數x4的下降趨勢最快,影響最小的為參數x2。

圖19 各參數對電抗器最高溫度的影響曲線Fig.19 Impact curves of parameters on the maximum temperature of reactor

4.4 基于粒子群算法的優化結果

考慮到影響電抗器溫度的遮雨帽和隔音裝置結構參數眾多,為獲得最佳的結構參數,以提高電抗器散熱能力,采用粒子群優化算法(partical swarm optimization,PSO)。PSO算法來源于鳥類群體活動的規律性,它模仿鳥類捕食行為,將優化問題的搜索空間類比于鳥類的飛行空間,流程圖如圖20所示[24 - 25]。其中,設置最大迭代次數為100、粒子種群規模為10、慣性權重為1.0、學習因子c1=c2=1.5。

圖20 PSO-Kriging優化流程圖Fig.20 Optimization flow chart of PSO-Kriging

基于上述的Kriging近似模型,結合粒子群優化算法PSO,獲得了最佳的遮雨帽和隔音裝置結構參數,如表3所示。

根據表3中的最優結構參數,在保持電抗器包封線圈電氣和結構參數恒定的條件下,結合流場-溫度場仿真計算方法,得到在最佳遮雨帽和隔音裝置參數下電抗器溫度場仿真結果,如圖21所示。

表3 設計變量最優參數Tab.3 Optimal parameters of design variables

圖21 最優參數下電抗器溫度場仿真結果Fig.21 Temperature field simulation result of reactor with optimal parameters

由圖21可知,電抗器最高溫度為81.2 ℃,最高溫升僅為61.2 ℃,相比于優化前電抗器的溫升降低了12.4 ℃,仿真結果驗證了優化方法的正確性。

5 結論

本文建立了空心電抗器流場-溫度場三維仿真模型,分析了加遮雨帽和隔音裝置前后電抗器的散熱特性,通過優化遮雨帽和隔音裝置結構參數以提高電抗器的散熱能力,可以得到如下結論。

1)獲得了加裝遮雨帽和隔音裝置前后電抗器包封線圈溫升分布特點,加裝遮雨帽和隔音裝置后氣道流體的最高流速從1.2 m/s下降到0.8 m/s,流體流速的降低造成電抗器溫升顯著增加。

2)將電抗器溫度場仿真計算方法和拉丁方試驗設計相結合,獲得了不同遮雨帽和隔音裝置結構參數下的溫度場仿真結果,構建了Kriging近似模型,采用靈敏度分析技術,獲得了不同結構參數對電抗器最高溫度的影響程度。其中,隔音裝置上端到遮雨帽底端距離x3、隔音裝置頂端中心孔半徑x4、隔音裝置底端到線圈底端距離x7和隔音裝置底端中心孔半徑x8對電抗器最高溫度影響顯著。

3)根據建立的Kriging近似模型,結合粒子群優化算法獲得了最佳的遮雨帽和隔音裝置結構參數,在最優的結構參數下電抗器最高溫升從73.6 ℃下降為61.2 ℃,優化方法能夠顯著降低電抗器溫升。

主站蜘蛛池模板: 欧美中文一区| 中文字幕波多野不卡一区| 国产精品视频第一专区| 国产黄在线免费观看| 在线免费观看a视频| 色天堂无毒不卡| 国产丝袜一区二区三区视频免下载| 国产福利免费在线观看| 欧美色伊人| 亚洲国产天堂久久综合| 亚洲中文字幕精品| 国产精品30p| 中文字幕丝袜一区二区| 91亚洲免费视频| 国产在线一区视频| 国产精品成| 91无码网站| 69综合网| 久热中文字幕在线| 99999久久久久久亚洲| 无码中文字幕精品推荐| 亚洲 欧美 偷自乱 图片 | 亚洲an第二区国产精品| 中文字幕精品一区二区三区视频 | 亚洲视频一区在线| 国产一级毛片在线| 麻豆精品视频在线原创| 国产香蕉一区二区在线网站| 国产白浆视频| 日韩欧美91| 亚洲a级毛片| 国产一级一级毛片永久| 国产亚洲精品va在线| 亚洲中文字幕无码mv| 欧美有码在线| 久久精品无码中文字幕| 午夜啪啪福利| 国产理论一区| 亚洲国产天堂在线观看| 亚洲人成网站观看在线观看| 欧美区一区二区三| 国产女人综合久久精品视| 亚洲Va中文字幕久久一区| 在线高清亚洲精品二区| 国产精品亚洲日韩AⅤ在线观看| 国产男女免费完整版视频| 欧美午夜一区| 国产成在线观看免费视频| 少妇露出福利视频| 亚洲码在线中文在线观看| 99久久精品免费看国产免费软件| 欧美国产日韩在线观看| 亚洲国产精品无码久久一线| 特黄日韩免费一区二区三区| 亚洲婷婷六月| 免费看一级毛片波多结衣| 四虎AV麻豆| 毛片免费在线视频| 国产激情国语对白普通话| 亚洲综合中文字幕国产精品欧美| 精品少妇人妻无码久久| 婷婷午夜天| 日韩在线2020专区| 97在线视频免费观看| 黄色在线不卡| 波多野结衣一区二区三视频 | 久久频这里精品99香蕉久网址| 欧美成人看片一区二区三区| 精品国产免费观看| 亚洲天堂日韩av电影| 99精品福利视频| 狠狠做深爱婷婷久久一区| 久久天天躁夜夜躁狠狠| 亚洲高清无码久久久| 亚洲国产欧美中日韩成人综合视频| 欧美日韩在线成人| 无码中文字幕精品推荐| 狠狠色婷婷丁香综合久久韩国| 一级毛片中文字幕| 日韩无码黄色网站| 亚洲综合色婷婷中文字幕| 国产区成人精品视频|