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

超級質子-質子對撞機中束流熱屏的熱-結構耦合模擬分析*

2021-01-14 02:47:58范佳錕王潔高勇游志明王盛張靜胡耀程許章煉王斌
物理學報 2021年1期
關鍵詞:模型

范佳錕 王潔? 高勇 游志明 王盛?張靜 胡耀程 許章煉 王斌

1) (西安交通大學,陜西省先進核能技術重點實驗室&動力工程多相流國家重點實驗室&陜西省先進核能工程研究中心,西安 710049)

2) (西安交通大學機械學院,西安 710049)

束流熱屏(beam screen)是新一代高能粒子對撞機中的重要部件,用于將束流在管道中運行時產生的熱量轉移到冷卻系統中,同時通過束流熱屏上的排氣孔將殘余氣體輸送至冷管壁上,維持良好的真空度.然而,在轉移熱負載的過程中,溫度變化產生的形變會影響束流熱屏的結構穩定性.如何在保證束流熱屏良好傳熱性能的情況下,盡量減小形變是優化束流熱屏結構設計的關鍵問題之一.本文采用ANSYS軟件對束流熱屏模型的傳熱性能和力學性能進行了模擬,并優化了束流熱屏結構設計,增強其傳熱性能和結構穩定性.對于束流熱屏外屏的內表面,采用減小銅涂層厚度的方式來降低運行過程中產生的洛倫茲力.相關理論模型計算結果表明: 與厚度為 100 μm 的銅涂層工況相比,當銅涂層的厚度在 0到 100 μm 之間變化時,厚度為 75 μm的銅涂層可以使束流熱屏外屏的最大形變降低70.9%,同時使束流熱屏的最高溫度升高1.1%.對于束流熱屏內屏,采用間隔布置支撐肋片的設計方案對束流熱屏的結構進行加固處理,提高束流熱屏整體的結構穩定性.計算結果表明: 與未加支撐肋片的工況相比,當相鄰兩個支撐肋片之間的間隔為1個排氣孔時,束流熱屏內屏的最大形變可降低86.8%,同時使束流熱屏的最高溫度降低7.69%.研究成果為新一代高能粒子加速器真空系統中關鍵部件束流熱屏的設計提供重要的理論參考.

1 背景介紹

2012年,標準模型中希格斯粒子(Higgs boson)在歐洲核子中心被發現[1],但仍然存在著夸克-輕子對稱性、電荷量化、基本粒子數量以及暗物質等標準模型所不能解釋的現象[2-4].因此,探究超越標準模型的“新物理”成為粒子物理學的首要任務,也是超級質子對撞機(super proton-proton collider,SPPC)的主要目標之一[5].

高能粒子束在運行過程中產生的同步輻射等會在管道內壁產生能量沉積,造成管道溫度升高,繼而影響束流的穩定運行[6-9].束流熱屏作為高能粒子加速器真空系統中的關鍵部件之一,通過將束流熱屏內的熱載荷轉移到冷卻管道內,維持超導磁鐵運行所需的工作溫度[10,11].同步輻射轟擊在束流熱屏內壁會釋放出大量氣體,因此,束流熱屏上必須存在貫通的排氣孔,使同步輻射致解吸的氣體分子可以從中逸出,并且凝結在冷管道內壁上[12,13].

束流運行過程中產生的溫度場變化會在束流熱屏結構中產生位移和應變.同時,考慮到超導失超的情況[14,15],在進行結構分析時需要考慮洛倫茲力的因素,進行熱-結構耦合分析.因此,本文采用ANSYS軟件對束流熱屏設計方案進行改進和優化,以增強其結構性能和傳熱性能,從而保證束流運行的穩定性.

2 模型與方法

2.1 模型介紹

束流熱屏模型[16-19]的結構和尺寸如圖1所示.在本文的模擬中,SPPC的束流熱屏主體材料選擇316LN不銹鋼[20-23].其中,不銹鋼表面鍍了一層銅涂層,銅材料選擇剩余電導率RRR(residual conductivity ratio)為100的高電導率無氧銅[17].當溫度為20 K時,不銹鋼和無氧銅的物性參數如表1所列[20-26].

圖1 束流熱屏模型的結構和尺寸Fig.1.Structure and size of beam screen model.

2.2 熱負載與洛倫茲力

2.2.1 熱負載

束流熱屏上的熱負載主要有三個來源,分別是同步輻射、鏡像電流和電子云產生的熱負載[6,12].在SPPC中,輻射沿水平面與軌道相切地發射,是束流熱屏上熱負載的主要來源[10,11].管壁鏡像電流產生的熱負載與束流熱屏材料的電導率相關[6].本文采用在束流熱屏內壁上鍍一層高電導率無氧銅的方法來降低鏡像電流產生的熱負載.束流與殘余氣體相互作用引起的相關效應最終導致電子云的形成,在束流熱屏內壁產生熱負載[27,28].綜合考慮三種熱負載,SPPC的束流在運行過程中產生的總熱負載為17.08 W/m[16].

2.2.2 洛倫茲力

磁場引起的洛倫茲力作用在束流熱屏上引起形變,使其在水平方向上張緊,而在垂直方向上收縮.超導磁鐵失超時的磁場會在束流熱屏內部感應出渦流[14,15],渦流的強度由感生電場的強度以及材料的電阻率決定.渦流在磁場結構內部沿著閉合回路運動,在磁場的作用下產生洛倫茲力,作用在束流熱屏結構上.對SPPC的束流熱屏而言,超導磁鐵失超過程中的最大洛倫茲力為

其中r為束流熱屏結構的半徑;Δ1為316 LN不銹鋼的厚度;Δ2為無氧銅的厚度;ρss為 316 LN 不銹鋼的電阻率;ρCu為無氧銅的電阻率;B為磁感應強度,12 T;為磁感應強度的變化率,63 T /s .對于SPPC而言,束流熱屏內屏邊緣處和外屏中心處的洛倫茲力最大值分別為23.6 N/mm和15.2 N/mm.

2.3 網格獨立性分析

本文研究的束流熱屏模型采用ANSYS meshing 進行網格劃分.在冷卻劑管道溫度為 80 K,束流運行過程中產生的總的熱負載為17.08 W/m時,對束流熱屏模型進行網格獨立性分析.根據劃分網格幾何尺寸的不同,得到如圖2所示的五種網格數目的最高溫度計算結果.由圖2可知,在網格數大于248956時,束流熱屏的最高溫度變化小于0.01 K.考慮到計算精度和效率,束流熱屏所用的網格數為248956.

表1 溫度為 80 K 時,無氧銅和 316LN 不銹鋼的材料屬性Table 1.Material properties of oxygen-free copper and 316LN stainless steel at 80 K.

圖2 五種網格數目的束流熱屏最高溫度計算結果Fig.2.The calculation result of the maximum temperature of beam screen with five kinds of the grids.

2.4 可靠性分析

為了驗證所采用的束流熱屏模型以及網格劃分方式的正確性,使用COMSOL仿真軟件在相同的工況下對同一模型進行驗證.在冷卻劑管道溫度為 80 K,排氣孔開口面積比為 6% 時,COMSOL的模擬結果為83.756 K,ANSYS在相同條件下的模擬結果為 83.613 K,誤差為 0.17%,在可接受的范圍內,因此可以驗證本文所采用的束流熱屏模型以及網格劃分方式的正確性.

3 熱-結構耦合模擬

在80 K的工況下對束流熱屏模型進行溫度場模擬,所涉及關鍵參數如表1 所列,316 LN 不銹鋼和無氧銅的導熱系數分別為 8.12 W m—1K—1和 540 W m—1K—1,將得到的溫度場模擬結果作為載荷施加到靜力學結構分析模型上,所涉及關鍵參數如表1所列.根據實際情況,將所有接觸面的接觸關系定義為Bonded,對冷管壁外表面施加固定約束,限制其位移,對束流熱屏兩端橫截面施加無摩擦約束,限制其法向位移,允許進行徑向位移,其他表面不施加約束.

3.1 無洛倫茲力情況下的熱-結構耦合模擬結果

在無洛倫茲力情況下,當冷卻管道內溫度為80 K,排氣孔開口面積比為6%時,對束流熱屏的形變量進行分析.通過如圖3(a)所示的四個路徑,對束流熱屏上的形變量進行分析.在圖3(b)中,束流熱屏距離為0處是指束流熱屏內屏與上冷卻管道相切處,其他距離是指以相切處為起點,沿著路徑1上某一點與相切處的距離.束流熱屏內屏處形變量沿著路徑1的變化如圖3(b)所示.

由計算結果可以得到: 沿著路徑1,束流熱屏內屏形變量沿著路徑1基本呈線性變化,在接近底部的位置形變僅為0.006 mm,在頂端則達到了最大形變量 0.031 mm,總體形變量較小,相對最大形變量為 2.3%.在圖3(c)中,束流熱屏距離為0處是指束流熱屏外屏最左端,其他距離是指以最左端為起點,路徑3上某一點與起點間的距離.束流熱屏外屏的形變量隨距離的變化如圖3(c)所示,形變呈波狀分布,變化趨勢與排氣孔位置的分布相關.束流熱屏的向外形變位于排氣孔之間,最大形變量為0.074 mm; 束流熱屏的向內形變位于相鄰排氣孔的間隙處,最大形變量為0.065 mm.

3.2 有洛倫茲力情況下的熱-結構耦合模擬結果

在束流熱屏內屏邊緣處以及束流熱屏外屏中心處分別施加大小為23.6和15.2 N/mm的洛倫茲力,水平指向束流熱屏外側.圖4為冷卻管道內溫度為80 K,排氣孔開口面積比為6%的工況下,施加了洛倫茲力的熱-結構耦合模擬結果.可以看到形變的最大值出現在束流熱屏外屏受同步輻射處,最大形變量為0.428 mm.總形變可分解為沿X,Y,Z三個方向的定向形變,分別如圖4(b),圖4(c)和圖4(d)所示.

圖4(b)為X方向上的形變分布云圖,其中紅色代表和X方向同向的形變,藍色代表和X方向反向的形變.可以看到X方向上的形變集中在束流熱屏外屏中心處和內屏邊緣處,分別向兩側發生形變,最大形變量為 0.416 mm.束流熱屏內屏邊緣處的X方向上的形變方向和束流熱屏外屏形變方向一致,指向束流熱屏外側.這是由于洛倫茲力產生的形變大于溫度引起的形變引起的,從而使束流熱屏內屏的形變方向變為相反方向.在束流熱屏的其他位置X方向上的形變量很小,束流熱屏X方向上總體的形變呈現出向兩側擴張的趨勢.

圖3 不同路徑和束流熱屏形變隨不同路徑的變化 (a) 路徑1、路徑2、路徑3和路徑4的示意圖; (b) 束流熱屏內屏形變量沿著路徑1的變化; (c) 束流熱屏外屏形變量沿著路徑3的變化.Fig.3.Different paths and variations with distance along paths in beam screen.(a) Schematic of paths 1,2,3 and 4; (b) variations of the inner screen with distance along path 1 in beam screen; (c) variations of the outer screen with distance along path 3 in beam screen.

圖4(c)為Y方向上的形變分布云圖,與X方向上的形變分布云圖相似,紅色和藍色分別代表和Y方向同向和反向的形變.可以看到束流熱屏在Y方向上的形變集中在冷卻劑管道周圍,束流熱屏上部的形變量略大于下部的形變量,最大形變量為 0.273 mm.因此,重力加劇了束流熱屏上部的形變,束流熱屏下部的形變正好相反,束流熱屏兩側對下部進行拉扯引起的形變方向與重力作用方向相反.因此,重力作用在一定程度上緩解了束流熱屏下部的形變.在束流熱屏其他位置Y方向上的形變很少,束流熱屏Y方向上總體的形變呈現向內部壓縮的趨勢.

圖4(d)為Z方向上的形變分布云圖,形變集中在排氣孔兩側,呈擠壓排氣孔的趨勢.X,Y和Z方向上的形變量最大值分別為0.416,0.273和0.014 mm,形變主要出現在X和Y方向上.出現這種情況的原因是: 排氣孔周圍溫度分布平均且較低,并且上下兩側有束流熱屏其他部分結構上的支撐來保持排氣孔部分的結構穩定性.

有洛倫茲力的情況下,束流熱屏內屏形變量沿著路徑1的變化如圖5(a)所示,形變量的范圍是在0.2—0.387 mm之間.相比于無洛倫茲力作用下的最大形變量(0.031 mm),有洛倫茲力存在情況下的最大形變量提高了近12倍.除了數值上的變化,有無洛倫茲力兩種工況下,束流熱屏內屏處的形變另一個區別在于兩者的方向是相反的.因此,在洛倫茲力作用前后,束流熱屏內屏頂端形變變化最大,由指向內側的最大形變量0.031 mm變為指向外側的 0.387 mm,形變變化量為 0.418 mm,相對最大形變量為29.7%.

圖4 有洛倫茲力情況下的熱-結構耦合模擬結果 (a) 總體形變分布; (b) X 方向上的形變分布云圖; (c) Y 方向上的形變分布云圖; (d) Z 方向上的形變分布云圖Fig.4.Simulation results of thermal-structure coupling with Lorentz force: (a) Overall deformation distribution; (b) cloud map of deformation distributions in X direction; (c) cloud map of deformation distributions in Y direction; (d) cloud map of deformation distributions in Z direction.

圖5 束流熱屏形變隨不同路徑的變化 (a) 束流熱屏內屏處形變量沿著路徑1的變化; (b) 束流熱屏外屏邊緣處沿著路徑2的形變量; (c) 束流熱屏外屏形變量沿著路徑3的變化Fig.5.Variations with distance along paths in beam screen: (a) Variations of the inner screen along path 1 in beam screen; (b) variations of the edge of outer screen along path 2 in beam screen; (c) variations of the outer screen along path 3 in beam screen.

束流熱屏外屏邊緣處形變量沿著路徑2的變化如圖5(b)所示.路徑2的起點為束流熱屏外屏與上冷卻管道的相切點.在0—7 mm處即接近束流熱屏頂部冷卻劑管道處,Y方向上的形變占主要地位.隨著位置向下移動,在 7—18 mm 處即排氣孔處的形變主要由X和Y方向上的形變構成,兩者的形變數值較小.因此,在排氣孔處的形變小于排氣孔兩側的形變,在18—28 mm處即束流熱屏外屏中心上側,X方向上的形變占主要地位,并且在中心點處達到了最大值.束流熱屏外屏下側部分變化趨勢與上側總體呈對稱分布,其中數值上不對稱是由于重力作用的影響.

有洛倫茲力的情況下,束流熱屏外屏形變量沿著路徑3的變化如圖5(c)所示,其外屏形變量變化趨勢與無洛倫茲力的情況下的變化趨勢相似,但最大形變量由 0.074 mm 增加到了 0.379 mm,最小形變量則由0.065 mm增加到了0.370 mm.

4 束流熱屏的結構穩定性強化

在無洛倫茲力作用的情況下,束流熱屏的形變主要出現在束流熱屏內屏處,最大形變量為0.031 mm,方向指向束流熱屏內側; 在有洛倫茲力作用的情況下,束流熱屏的形變主要出現在束流熱屏內屏及束流熱屏外屏中心處,最大形變量為 0.387 mm,方向指向束流熱屏外側.束流熱屏內屏和束流熱屏外屏的相對最大形變量分別為29.7%和42.1%.因此,需要考慮對束流熱屏進行結構上的優化和改進,增強其機械性能,從而使束流穩定運行.

4.1 束流熱屏外屏的結構穩定性強化

保持束流熱屏的結構穩定性有三個方面: 一是改變材料,采用強度更大的材料作為束流熱屏的主體材料; 二是對形變較大處進行加固處理,如增加厚度等; 三是降低束流熱屏上的洛倫茲力,從而減小形變.本文對束流熱屏外屏的結構穩定性強化主要基于后面兩個方面.對于束流熱屏外屏,其與外部冷管壁之間的距離僅為1 mm,采用加固處理的方式有局限性且效果不明顯.因此,通過降低運行過程中產生的洛倫茲力來對束流熱屏外屏進行結構穩定性強化.

同等條件下,無氧銅內洛倫茲力遠大于不銹鋼內的洛倫茲力.因此,降低銅涂層的厚度是降低洛倫茲力的有效方式,但銅涂層厚度的增加有利于增強傳熱.在80 K的溫度下,銅和不銹鋼的導熱系數分別為 540 W/mK 和 8.12 W/mK,兩者相差近67倍.因此,需要分析不同厚度的銅涂層對傳熱性能的影響.對于束流熱屏內屏,在沒有洛倫茲力作用的情況下,也存在著一定的形變.單純地降低洛倫茲力不能很好地保持其結構穩定性,考慮到束流熱屏內屏與束流熱屏外屏之間空間較大,可以采用加固處理的方式增加其機械強度.

4.1.1 不同銅涂層厚度對束流熱屏傳熱性能和結構穩定性的影響

在冷卻管道內溫度為80 K,排氣孔開口面積比為6%的情況下,分別對銅涂層厚度為0,25,50,75和100 μm的束流熱屏模型進行溫度場模擬分析.由圖6 可知,與 100 μm 厚銅涂層工況相比,銅涂層厚度為 0,25,50 和 75 μm 的束流熱屏模型最高溫度分別增加了 31.6%,7.9%,3.2% 和 1.1%.銅涂層的存在對束流熱屏的傳熱能力有著很大的影響,但溫度分布趨勢則基本相同.

圖6 不同銅涂層厚度下的束流熱屏模型最高溫度Fig.6.Maximum temperature of beam screen model under different copper plating thickness.

由于銅涂層厚度發生改變,作用在束流熱屏外屏上的洛倫茲力也發生了改變.根據(1)式,銅涂層厚度為 0,25,50,75 和 100 μm 時,洛倫茲力分別為 0.58,4.23,7.88,11.53 和 15.18 N/mm.在冷卻管道內溫度為80 K,排氣孔開口面積比為6%的情況下,分別對上述五種束流熱屏模型進行熱-結構耦合分析.由圖7 可知,相比于 100 μm 厚銅涂層的束流熱屏外屏最大形變量,銅涂層厚度為0,25,50和75 μm的束流熱屏外屏最大形變量分別為 0.059,0.062,0.087,0.11 和 0.379 mm,即分別降低了 84.4%,83.6%,77.0% 和 70.9%.當銅涂層厚度為100 μm時,形變量最大且遠高于其余四種厚度的束流熱屏外屏形變量,其余四種厚度的束流熱屏外屏形變量在0.059—0.11 mm之間.綜合考慮傳熱性能和結構穩定性的因素,選擇使用厚度為75 μm的銅涂層比較合適.

圖7 不同銅涂層厚度下的束流熱屏外屏最大形變量Fig.7.Maximum deformation of beam screen outside screen under different copper coating thickness.

4.2 束流熱屏內屏的結構穩定性強化

束流熱屏內外屏之間的空隙在束流運行過程中起到了氣體分子逸散通道的作用.因此,采用在束流熱屏內屏外側增加支撐肋片的方式來對束流熱屏內屏進行結構穩定性強化,如圖8(a)所示.

本文選取以下五種不同的支撐肋片分布方式,即每兩個支撐肋片之間分別間隔1個、2個、3個、4個和5個排氣孔.圖8(b)和圖8(c)分別是相鄰支撐肋片間隔1個排氣孔時的支撐肋片分布以及束流熱屏模型溫度分布圖,其中藍色部分為支撐肋片.接下來將分析不同的支撐肋片分布對束流熱屏傳熱性能和結構穩定性的影響.

4.2.1 支撐肋片結構對束流熱屏傳熱性能的影響

在冷卻劑管道溫度為80 K的情況下,分別對五種束流熱屏模型進行熱-結構耦合分析,計算結果如圖9所示.在沒有支撐肋片存在的情況下,束流熱屏的最高溫度為83.613 K.支撐肋片間隔排氣孔數量分別為1個、2個、3個、4個和5個的束流熱屏模型的最高溫度分別為83.335,83.388,83.542,83.577 和 83.588 K,這五種工況下的最高溫度變化量相比于無支撐肋片存在的束流熱屏模型的最高溫度變化量分別降低了7.69%,6.22%,1.96%,0.99% 和 0.69%.

4.2.2 不同支撐肋片分布對束流熱屏內屏結構穩定性的影響

在冷卻管道內溫度為80 K,排氣孔開口面積比為 6%,洛倫茲力分別為 23.62和 15.18 N/mm的工況下,對五種不同支撐肋片分布的束流熱屏模型進行熱-結構耦合分析.圖10為支撐肋片間隔排氣孔數量分別為0個、1個、2個、3個、4個和5個的束流熱屏內屏邊緣處沿著路徑4的形變分布.

圖8 (a)增加支撐肋片的束流熱屏模型; (b)相鄰支撐肋片間隔排氣孔數量為1個時的支撐肋片分布; (c)相鄰支撐肋片間隔排氣孔數量為1個時束流熱屏模型的溫度分布圖Fig.8.(a) Beam screen model with added support ribs; (b) distributions of support fins when the number of exhaust holes between adjacent support fins is one; (c) temperature distribution diagram of the beam screen model when the number of exhaust holes between adjacent support fins is one.

圖9 不同支撐肋片分布下的束流熱屏最高溫度Fig.9.Maximum temperature of beam screen under different supporting ribs distribution.

圖10 6種不同支撐肋片束流熱屏內屏邊緣處沿著路徑4的形變分布Fig.10.Deformation distributions along path 4 at the edge of inner screen of six beam screens with different number of support ribs.

圖11 六種束流熱屏的內屏邊緣處最大形變量Fig.11.Maximum deformation at the edge of the inner screen of six beam screens.

由圖11可知,無支撐肋片時,束流熱屏的內屏邊緣處最大形變量為0.387 mm,支撐肋片間隔排氣孔數量分別為1個、2個、3個、4個和5個時的束流熱屏模型最大形變量分別為0.051,0.179,0.307,0.371 和 0.395 mm.相比無支撐肋片的束流熱屏模型的最大形變量,這五種工況的最大形變量分別降低了 86.8%,53.7%,20.6%,4.1% 和 2%.

由圖10和圖11可知,當束流熱屏模型的支撐肋片間隔排氣孔數量為1個時,束流熱屏內屏的形變量最低且遠低于其他五種模型的束流熱屏內屏形變量,而且分布更加均勻.因此,選擇支撐肋片間隔排氣孔數量為1個的支撐肋片分布模式來增強束流熱屏內屏的結構穩定性.

5 結 論

本文基于ANSYS 軟件,針對束流熱屏冷卻管道內溫度為80 K,排氣孔開口面積比為6%的工況,對束流熱屏的傳熱性能和結構穩定性進行了綜合分析.在溫度場和洛倫茲力共同作用下,束流熱屏X方向上最大形變量為0.416 mm.Y方向上的最大形變量為0.273 mm.Z方向上的形變量最大值僅為0.014 mm.為了增強束流熱屏的的結構穩定性,采用減小銅涂層厚度的方式來降低運行過程中產生的洛倫茲力,同時采用加固處理的方式來增加束流熱屏內屏的結構穩定性.銅涂層厚度分別為 0,25,50,75 和 100 μm 時,最高溫度分別為110.09,90.302,86.295,84.569 和 83.613 K,銅涂層的厚度對束流熱屏傳熱性能有著很大的影響,相應的最大形變量分別為 0.059,0.062,0.087,0.11和0.379 mm.銅涂層厚度的變化帶來的洛倫茲力的變化造成了束流熱屏外屏最大形變量的變化.綜合考慮傳熱性能和結構穩定性的因素,選擇使用75 μm的銅涂層.在沒有支撐肋片存在的情況下,束流熱屏的最高溫度為83.613 K.支撐肋片間隔排氣孔數量分別為1個、2個、3個、4個和5個的束流熱屏模型的最高溫度變化量相比于無支撐肋片存在的束流熱屏模型的最高溫度變化量分別降低了 7.69%,6.22%,1.96%,0.99% 和 0.69%.當束流熱屏模型的支撐肋片間隔排氣孔數量為1個時,束流熱屏內屏的形變量最低且遠低于其他五種模型的束流熱屏內屏形變量.因此,合理設置支撐肋片能夠增強束流熱屏的傳熱能力和結構穩定性.

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 天天综合网色中文字幕| 国产久操视频| 99在线视频免费观看| 99免费视频观看| 欧美一区二区三区国产精品| 国产成人乱码一区二区三区在线| 在线亚洲天堂| 国产三级成人| AV在线天堂进入| 夜夜操国产| 欧美乱妇高清无乱码免费| 日韩无码视频播放| 白丝美女办公室高潮喷水视频| 97在线免费| 日韩毛片免费观看| 91精品国产一区自在线拍| 依依成人精品无v国产| 日韩视频免费| 深爱婷婷激情网| 久久久受www免费人成| 国产熟女一级毛片| 亚洲无码91视频| 国产91精品久久| 中文字幕永久在线看| 欧洲高清无码在线| 日韩人妻少妇一区二区| 亚洲成在人线av品善网好看| 国产精品午夜电影| 一边摸一边做爽的视频17国产| 成人福利免费在线观看| 国产高清又黄又嫩的免费视频网站| 免费看美女毛片| 免费可以看的无遮挡av无码| 国产人人射| 亚洲av日韩综合一区尤物| 亚洲精选高清无码| 91精品国产丝袜| 欧美激情综合| 日本免费a视频| 99人体免费视频| 亚洲免费人成影院| 国内精品视频| 国产精品第页| a级毛片免费播放| 欧美精品v| 日韩成人免费网站| 动漫精品啪啪一区二区三区| 久久久成年黄色视频| 欧美午夜性视频| 国产网站在线看| 不卡午夜视频| 91青青草视频在线观看的| 99国产在线视频| 激情综合网激情综合| 亚洲AV无码不卡无码 | 国产一区亚洲一区| 久久香蕉国产线| 另类专区亚洲| 免费网站成人亚洲| 亚洲人精品亚洲人成在线| 国产精品第5页| 日韩毛片免费观看| 日韩av高清无码一区二区三区| 亚洲无线一二三四区男男| 国产精品极品美女自在线看免费一区二区| 无码电影在线观看| 网友自拍视频精品区| 国产青青操| 99国产精品国产| 人妻无码中文字幕一区二区三区| 亚洲色图综合在线| 四虎影视永久在线精品| 久久永久免费人妻精品| 99免费在线观看视频| aⅴ免费在线观看| 一级一级特黄女人精品毛片| 午夜丁香婷婷| 国产精品一区在线麻豆| 99re热精品视频中文字幕不卡| 国产乱子伦一区二区=| 狼友视频国产精品首页| 九九久久精品免费观看|