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

ITER 裝置中等離子體旋轉和反饋控制對電阻壁模影響的數(shù)值研究*

2021-03-04 05:55:00曹琦琦劉悅王碩
物理學報 2021年4期

曹琦琦 劉悅? 王碩

1) (大連理工大學物理學院, 三束材料改性教育部重點實驗室, 大連 116024)

2) (核工業(yè)西南物理研究院, 成都 610041)

在托卡馬克等離子體中, 電阻壁模是非常重要的磁流體不穩(wěn)定性, 特征時間在毫秒量級.對長時間穩(wěn)態(tài)運行下的先進托卡馬克, 電阻壁模限制著聚變裝置的運行參數(shù)空間(放電時間和比壓), 影響經濟效益, 所以研究電阻壁模穩(wěn)定性至關重要.本文使用MARS 程序, 針對ITER 裝置上9 MA 先進運行平衡位形, 研究了等離子體旋轉和反饋控制對電阻壁模的影響.結果表明, 在沒有反饋控制時, 當比壓參數(shù) C β 取0.7, 等離子體環(huán)向旋轉頻率達到1.1%的阿爾芬頻率時, 可以完全穩(wěn)定電阻壁模; 在等離子體環(huán)向旋轉和反饋控制共同作用時, 比壓參數(shù) C β 取0.7, 反饋增益 | G| 取0.6 時, 穩(wěn)定電阻壁模所需要的等離子體旋轉頻率為0.2%的阿爾芬頻率.可見, 單獨靠等離子體環(huán)向旋轉穩(wěn)定電阻壁模所需的旋轉頻率較大; 而等離子體環(huán)向旋轉和反饋控制共同作用可以降低穩(wěn)定電阻壁模的旋轉頻率臨界值, 符合先進托卡馬克的運行.本文的研究結果對中國聚變工程試驗堆CFETR 的工程設計和運行具有一定指導意義.

1 引 言

托卡馬克實驗中可以通過提高裝置的比壓來實現(xiàn)聚變功率和自舉電流份額的提升.而隨著裝置比壓的增加, 由壓強或壓強梯度驅動的磁流體不穩(wěn)定性可能會出現(xiàn), 外扭曲模是其中最重要的磁流體不穩(wěn)定性之一.外扭曲模對托卡馬克裝置放電產生重大影響, 嚴重時甚至直接導致放電破裂或終止.早期的研究表明, 當托卡馬克等離子體外部有足夠靠近的理想導體壁時, 外扭曲模可以完全被穩(wěn)定住.所謂理想導體壁是指電阻為零的導體壁.在實際實驗中, 真空室或第一壁材料等都有一定的電阻, 有電阻的導體壁稱為電阻壁.這類電阻壁使外扭曲模演化為一種增長率緩慢的不穩(wěn)定模式, 由于這種不穩(wěn)定性與電阻壁有關, 被稱為電阻壁模(resistive wall mode, RWM)[1?3].電阻壁模的特征時間是磁場在電阻壁滲透的特征時間, 為毫秒量級[4].電阻壁模一旦產生, 仍然會造成等離子體的扭曲, 嚴重時也會導致放電破裂.因此, 對于需要長時間穩(wěn)態(tài)運行的先進托卡馬克, 電阻壁模的研究是非常重要的.

2004 年Liu 等[5]對ITER 裝置的反饋控制和等離子體環(huán)向旋轉穩(wěn)定電阻壁模進行了研究, 結果表明, 安裝在真空室內的線圈和極向傳感器組成的反饋控制系統(tǒng)能夠有效地抑制電阻壁模, 并指出流體效應下等離子體旋轉對電阻壁模穩(wěn)定的有效性.Xia 等[6]針對JET-like 的平衡位形研究了等離子體旋轉效應和反饋控制的協(xié)同作用, 結果表明等離子體旋轉與反饋控制的協(xié)同作用下能促使電阻壁模產生兩個穩(wěn)定窗, 并且在反饋控制下, 等離子體電阻對電阻壁模有很大影響.Xia 等[7]還研究了在反饋控制下動理學效應對電阻壁模的影響.Liu 等[8]研究了ITER 裝置中等離子體剪切流對電阻壁模的影響.Hao 等[9?11]研究了捕獲熱粒子對電阻壁模的影響, 捕獲能量粒子穩(wěn)定電阻壁模, 以及由黏滯引起的湍流和等離子體旋轉對電阻壁模不穩(wěn)定性的協(xié)同作用.

由于ITER 裝置的諸多物理參數(shù)與已有的其他聚變裝置不同, 如ITER 裝置幾何尺寸更大, 需在高比壓、大自舉電流、長脈沖下運行.ITER 裝置目標平衡對應的歸一化比壓 βN= 2.92, 對應的比壓參數(shù) Cβ= 0.44, 本文模擬中使用的平衡是在ITER 設計平衡的基礎上略調高了比壓值, Cβ= 0.7,βN= 3.18.ITER 裝置希望運行在更高的比壓參數(shù)下, 以獲得更高的聚變功率、自舉電流份額等.而電阻壁模會隨著 βN的升高, 變得更加不穩(wěn)定, 模式的增長率相對更高, 有利于提高模擬結果的精度(因為目標平衡在使用磁通-電壓控制模式且反饋增益為 | G| = 0 時, 模式的增長率已經進入準穩(wěn)定狀態(tài), 可見文獻[12]的圖3).對于ITER 這樣的超大型托卡馬克裝置, 放電破裂的破壞性更大, 對裝置中面向等離子體的材料將會造成毀滅性的傷害.因此, 在裝置運行中必須有效控制電阻壁模不穩(wěn)定性的增長.

對電阻壁模的控制一般可以分為被動控制和主動控制兩種手段.被動控制是利用模式與波發(fā)生共振阻尼[13,14], 把不穩(wěn)定模式的自由能耗散掉.主動控制是通過在真空室內或外部安裝的反饋線圈補償泄露的磁場[15], 從而控制電阻壁模.近期研究發(fā)現(xiàn)使用主動控制和被動控制的協(xié)同作用能更好地控制電阻壁模[6].

本文利用MARS-F 程序[16]對ITER 裝置電阻壁模不穩(wěn)定性進行了模擬, 分別研究了流體模型下等離子體旋轉、反饋控制及協(xié)同作用對電阻壁模不穩(wěn)定性增長率的影響.文章安排如下: 第二部分將簡要介紹使用的MARS-F 程序的模型, 第三部分將給出使用的ITER 平衡, 第四部分將對數(shù)值計算結果進行分析和討論, 最后進行總結.

2 計算模型

使用線性模型研究電阻壁模不穩(wěn)定性及其控制是目前所有模擬中公認的方法, 使用線性模型是因為電阻壁模的增長率尺度較小, 它的非線性項貢獻可以忽略不計.MARS[17]程序是一個線性代碼,在環(huán)幾何位形下自洽求解電阻磁流體方程組, 程序不僅考慮了環(huán)向等離子體旋轉, 而且把電阻導體壁、外加反饋線圈以及各種阻尼效應也引入到模型中.MARS 程序聯(lián)立等離子體區(qū)域的非理想磁流體線性化方程組、真空區(qū)域方程、薄壁近似以及外加反饋線圈方程求解得到不穩(wěn)定模式的本征值.下面簡單地介紹一下本文使用的MARS 程序模型.

等離子體的方程為[17,18]:

這里 γ 是復數(shù), 代表電阻壁模的増長率和實頻, 由多普勒位移 i n? 修正, n 是環(huán)向模數(shù), ? 是沿著環(huán)向 ? 的等離子體旋轉頻率.B, J 和P 分別為由平衡程序CHEASE[16]求解得到的平衡磁場、平衡電流密度和平衡壓強.v, b, j 和p 分別表示擾動速度、擾動磁場、擾動電流密度和擾動壓強.ρ 是擾動質量密度, R 是等離子體大半徑, Z? 是在垂直方向的單位矢量, Π 是黏性應力張量, 表示平行聲波黏滯阻尼[19],代表阻尼強度, k 表示聲波數(shù) ( n ?m/q)(1/R) ,表示離子熱速度.η 是等離子體電阻率, Γ 是比熱比系數(shù).

在真空區(qū)域中, 沒有擾動速度和擾動壓強, 且擾動磁場滿足:

導體電阻壁滿足薄壁近似:

其中 ηW和 hW分別表示導體壁的電阻率和厚度, br,bx分別代表是擾動磁場逆變量的分量, 在托卡馬克曲線坐標 ( r,χ,?) 中, 有

其中r 是徑向坐標, 定義為極向磁通的平方根, χ是極向坐標, 依賴于雅可比行列式的選擇.

反饋線圈方程為[12]

其中 ψf, If, Rf和 τf分別表示穿過反饋線圈的擾動磁通、線圈電流、電阻以及反饋線圈的響應時間,G 表示反饋系統(tǒng)增益, 反饋增益是工程參數(shù), 在控制過程中通過調整反饋控制系統(tǒng)的PID 參數(shù)來實現(xiàn)對系統(tǒng)的有效控制.本文中的G 相當于PID 系統(tǒng)的比例項系數(shù).這里G 是歸一化的參數(shù), 歸一化因子為無量綱量 Kp和 Kd表示反饋控制器中的比例和微分因子, bs為傳感器線圈內測量/響應的噪聲, 這里也是用磁場表示的.

3 平衡與比壓

采用ITER 裝置等離子體電流 IP= 9 MA 的先進運行平衡位形, 其能量增益Q = 5, 輸出功率P = 340 MW, 真空環(huán)向磁軸處磁場 B0= 5.3 T,等離子體大半徑 R0= 6.2 m, 小半徑 a = 2.0 m.ITER 裝置的極向截面如圖1 所示, 真空室分為內外兩層, 在內層真空室弱場側偏上、中間、偏下分別安裝了三組磁擾動線圈, 用于控制邊緣局域模和電阻壁模不穩(wěn)定性, 以及內側中心平面處的一組傳感器, 用來測量極向磁擾動.

圖1 ITER 裝置極向截面及反饋控制示意圖Fig.1.Geometry of poloidal cross section and feedback control in ITER.

圖2 給出了ITER 裝置9 MA 先進運行平衡的各平衡量剖面分布, 圖2(a)是安全因子, 在芯部有較弱的負磁剪切, 磁軸處的安全因子為 q0=2.44, 邊界處的安全因子達到 qa= 7.13, 安全因子最小值為 qmin= 1.60, 這種弱剪切的結構能產生較高的比壓, 是托卡馬克先進運行方案的一個重要特點.圖2(b)是等離子體壓強分布, 壓強使用磁壓歸一化, 圖2(c)是等離子體質量密度剖面, 圖2(d)是等離子體環(huán)向電流密度剖面, 使用B0/(μ0R0)進行了歸一化處理.

圖2 ITER 裝置9 MA 先進運行平衡 (a)安全因子剖面; (b)等離子體壓強剖面, 由 B 02/μ0 歸一化; (c)質量密度剖面, 磁軸處歸一化為1; (d)環(huán)向電流密度剖面, 由 B 0/(μ0R0) 歸一化, s = , ψ 是歸一化的極向通量Fig.2.Radial profiles of (a)safety factor; (b)equilibrium pressure normalized by factor B 02/μ0 ; (c)plasma density normalized to unity at the magnetic axis; (d)toroidal current density normalized by factor B 0/(μ0R0) and s = ψ is the normalized poloidal flux for ITER target plasma designed for 9 MA steady state scenario.

針對圖2 給出的平衡位形, 固定電阻壁的位置在 d /a=1.29 和環(huán)向模數(shù) n =1 , 運用MARS 程序計算得到了有理想壁和無壁的時候外扭曲模的增長率和比壓的關系, 如圖3 所示.所得到的計算結果表明, 在無壁時, 比壓極限為; 在有理想壁存在時, 比壓極限為=3.48 (其中βN=β(%)a(m)B0(T)/Ip(MA) 是 歸 一 化 比 壓, β 是等離子體熱壓與磁壓的比值, a 是等離子體小半徑,B0是環(huán)向等離子體磁軸處的磁場, Ip是等離子體總電流).這里定義一個等離子體比壓參量表示無壁時比壓極限, Cβ=1 表示有理想壁時比壓極限, 我們在研究電阻壁模的時候, Cβ會在0 和1 之間取值.

圖3 在無壁和有理想壁時, 不同比壓下的外扭曲模增長率Fig.3.Growth rate of external kink versus β N with and without an ideal wall.

4 數(shù)值結果與討論

4.1 流體效應對電阻壁模的穩(wěn)定作用

對于上面給出的ITER 平衡, 取等離子體比壓參量為 Cβ=0.7 , 研究等離子體旋轉頻率對電阻壁模的影響.對于平行黏滯 κ//分別取0.5, 1 和1.5,用MARS 程序計算得到了不同等離子體旋轉頻率下電阻壁模的增長率, 如圖4 所示, 其中等離子體旋轉頻率使用阿爾芬頻率 ?A=υA/R0歸一化的.所得到的計算結果表明, 在平行黏滯 κ//不變時, 電阻壁模增長率會隨著等離子體環(huán)向旋轉頻率增大而減小, 當旋轉頻率達到臨界頻率時, 電阻壁模的增長率為零.這是因為模式(波)與離子聲波在有理面附近發(fā)生共振阻尼, 耗散掉驅動不穩(wěn)定模式的自由能, 從而達到穩(wěn)定.所得到的計算結果還表明,平行黏滯越大, 穩(wěn)定電阻壁模所需的旋轉頻率越小, 穩(wěn)定效果越好.這是因為黏滯越大, 對壓強驅動的不穩(wěn)定模式的自由能耗散越大, 模式的增長率越小.

圖4 等離子體比壓參量為 C β =0.7 、平行黏滯 κ // 分別取0.5, 1 和1.5 時, 不同等離子體旋轉頻率下的電阻壁模增長率Fig.4.Growth rate of resistive wall mode versus plasma flow with equilibrium pressure scaling factor C β =0.7 ,parallel viscous coefficient κ // = 0.5, 1, 1.5.

當?shù)入x子體比壓參量 Cβ在0.1—0.9 之間變化,等離子體旋轉頻率在0— 1.2×10?2范圍內變化,用MARS 程序計算得到了電阻壁模的增長率, 計算中設置平行黏滯 κ//=1.5.計算結果如圖5 所示.在 Cβ比較小且等離子體具有較高旋轉頻率的情況下, 可以得到電阻壁模的穩(wěn)定區(qū)間.隨著 Cβ的增加, 電阻壁模表現(xiàn)得越不穩(wěn)定, 為了擴大穩(wěn)定電阻壁模的穩(wěn)定窗口, 在等離子體流體模型下, 僅利用等離子體旋轉實現(xiàn)ITER 裝置電阻壁模的控制是比較困難的, 需要反饋控制的加入.

4.2 反饋控制對電阻壁模的穩(wěn)定作用

用等離子體旋轉結合平行黏滯來穩(wěn)定電阻壁模需要很大的臨界等離子體旋轉閾值, 像ITER 等這種大型托卡馬克難以達到, 需要主動控制(反饋控制)參與.如圖1 所示, ITER 裝置的反饋系統(tǒng)是由三組反饋線圈和一組傳感器組成的, Li 等[20]分析了電壓-電流、電壓-電壓、磁通-電流以及磁通-電壓四種反饋方式, 本工作采用的反饋控制是磁通-電流反饋控制.

圖5 平行黏滯 κ // 取1.5 時, 不同等離子體比壓參數(shù)和不同等離子體旋轉頻率下的電阻壁模增長率Fig.5.With parallel viscous coefficient κ // =1.5 , growth rate of resistive wall mode versus plasma flow for different equilibrium pressure scaling factor C β.

首先給出主動反饋線圈產生的磁場分布, 如圖6 所示, 以 Cβ= 0.7 為例, 可以看出不加等離子體旋轉, 只考慮主動反饋線圈時, 在等離子體邊界有比較大的磁場分布, 芯部基本沒有線圈產生的磁場.

圖6 沒有等離子體旋轉頻率, C β =0.7 時, 計算得到的主動線圈產生的磁場分布Fig.6.Without plasma flow and equilibrium pressure scaling factor C β =0.7 , the calculated magnetic field distribution of active coil.

然后, 在研究三組線圈增益前, 為達到上下兩組線圈相位的最佳組合, 利用MARS-F 程序對上下兩組線圈的相位進行了掃描.首先固定上下兩組線圈增益為 | G| = 0.5, 固定一組線圈相位, 去掃描另一組線圈相位, 以找到最佳的相位組合, 如圖7所示, 電阻壁模的增長率隨上下兩組線圈的相位組合不同而不同, 在上下兩組線圈的相位為 ?u=150°, ?L= –150°時, 反饋控制電阻壁模的效果最好, 把這一相位組合稱為最佳相位組合, 和Wang等[12]采用磁通-電壓的反饋控制方式相位掃描的結果一致, 以下所有上下線圈的研究中, 上下線圈的相位均取 ?u= 150°, ?L= –150.

圖7 線圈增益幅值為 | G| = 0.5 時, 上下兩組線圈不同相位下電阻壁模的增長率Fig.7.Growth rate of resistive wall mode with varying phase of feedback gains for upper and lower sets of active coils, feedback gain amplitude | G| = 0.5.

接著, 在沒有等離子體旋轉頻率, 平行黏滯κ∥=1.5時, 分別研究了中間反饋線圈、上下兩組反饋線圈及三組反饋線圈的增益對電阻壁模的影響, 如圖8—圖10 所示.對比發(fā)現(xiàn), 等離子體比壓越大, 所需要穩(wěn)定電阻壁模的臨界增益越大; 在不同等離子體比壓下, 當反饋增益逐漸增大時, 增長率慢慢變小, 當增益幅值足夠大時, 能使得電阻壁模的增長率為零, 穩(wěn)定電阻壁模, 這是因為通過導體壁泄露的磁場已由中間反饋線圈得到補償, 此時導體壁等效為理想壁.定義完全穩(wěn)定電阻壁模的增益為臨界反饋增益 Gcri, 當G > Gcri時, 電阻壁模是穩(wěn)定的, 當G < Gcri時, 電阻壁模仍然會增長,最終影響高比壓等離子體放電.例如在 Cβ= 0.6時, 中間一組線圈穩(wěn)定電阻壁模所需的臨界增益為|G|= 1.2, 上下兩組線圈穩(wěn)定電阻壁模所需的臨界增益為 | G| = 1.6, 上中下三組線圈穩(wěn)定電阻壁模所需的臨界增益為 | G| = 0.7, 可見, 在只有反饋控制中, 上中下三組反饋線圈的控制效果更好.

4.3 電阻壁模的協(xié)同控制

圖8 在沒有等離子體旋轉頻率、平行黏滯 κ // =1.5 時,不同的等離子體比壓參量, 不同中間線圈的增益下電阻壁模的增長率變化Fig.8.Without plasma flow and with parallel viscous coefficient κ // =1.5 , growth rate of resistive wall mode with varying equilibrium pressure scaling factor versus feedback gains for middle sets of active coils.

圖9 在沒有等離子體旋轉頻率、平行黏滯 κ // =1.5 時,不同的等離子體比壓參量, 不同上下兩組線圈的增益下電阻壁模的增長率變化Fig.9.Without plasma flow and with parallel viscous coefficient κ // =1.5 , growth rate of resistive wall mode with varying equilibrium pressure scaling factor versus feedback gains for upper and lower sets of active coils.

前面兩節(jié)分別研究了被動控制和主動控制單獨對電阻壁模的穩(wěn)定作用, 當旋轉頻率超過一定閾值或反饋增益足夠大時, 電阻壁模能被完全穩(wěn)定住.然而在實際裝置運行中, 兩種控制手段同時存在, 例如裝置放電時給等離子體輔助加熱的中性束注入會產生等離子體旋轉, 控制其他不穩(wěn)定性模式的反饋線圈對電阻壁模也有穩(wěn)定作用等等.所以接下來研究ITER 裝置中等離子體旋轉與反饋對電阻壁模的協(xié)同作用.

圖10 在沒有等離子體旋轉頻率、平行黏滯 κ // =1.5 時,不同的等離子體比壓參量, 不同上中下三組線圈的增益下電阻壁模的增長率變化Fig.10.Without plasma flow and with parallel viscous coefficient κ // =1.5 , growth rate of resistive wall mode with varying equilibrium pressure scaling factor versus feedback gains for all three sets of active coils.

固定等離子體旋轉頻率 ?0/?A= 0.002,?0是等離子體中心處的旋轉頻率, 平行黏滯 κ//=1.5 ,黏滯提供了離子聲波阻尼穩(wěn)定機制, 考慮磁通-電流反饋控制, 發(fā)現(xiàn)反饋控制和旋轉共同作用比單獨反饋控制能夠更好地穩(wěn)定電阻壁模.例如 Cβ=0.6 時, 只有反饋控制存在(三組反饋線圈同時加入), 沒有等離子體旋轉頻率時, 完全穩(wěn)定電阻壁模所需要的反饋增益為 | G| = 0.7, 如圖10 所示; 當再加入等離子體旋轉頻率 ?0/?A= 0.002 時, 完全穩(wěn)定電阻壁模所需要的反饋增益為 | G| = 0.6, 如圖13 所示.可見, 反饋加入旋轉可以降低穩(wěn)定電阻壁模的臨界增益, 可以更好地控制電阻壁模.

圖11 在等離子體旋轉頻率 ? 0/?A = 0.002、平行黏滯κ// =1.5時, 不同的等離子體比壓參量, 加上旋轉后中間線圈的增益和增長率的變化Fig.11.With plasma flow ? 0/?A = 0.002 and parallel viscous coefficient κ // =1.5 , growth rate of resistive wall mode with varying equilibrium pressure scaling factor versus feedback gains for middle sets of active coils.

圖12 在等離子體旋轉頻率 ? 0/?A = 0.002、平行黏滯κ// =1.5時, 不同的等離子體比壓參量, 加上旋轉后上下兩組線圈的增益和增長率的變化Fig.12.With plasma flow ? 0/?A = 0.002 and parallel viscous coefficient κ // =1.5 , growth rate of resistive wall mode with varying equilibrium pressure scaling factor versus feedback gains for upper and lower sets of active coils.

圖13 在等離子體旋轉頻率 ? 0/?A = 0.002、平行黏滯κ// =1.5時, 不同的等離子體比壓參量, 加上旋轉后上中下三組線圈的增益和增長率的變化Fig.13.With plasma flow ? 0/?A = 0.002 and parallel viscous coefficient κ // =1.5 , growth rate of resistive wall mode with varying equilibrium pressure scaling factor versus feedback gains for all three sets of active coils.

接著分別研究在各個平衡壓強下, 中間線圈、上下兩組線圈和上中下三組線圈與等離子體旋轉的協(xié)同作用, 如圖11—圖13 所示.例如 Cβ= 0.6時, 在中間反饋線圈與等離子體旋轉的共同作用下, 穩(wěn)定電阻壁模需要的臨界增益為 | G| = 1.0; 在上下兩組線圈與等離子體旋轉的共同作用下, 穩(wěn)定電阻壁模需要的臨界增益為 | G| = 1.4; 在上中下三組線圈與等離子體旋轉的共同作用下, 穩(wěn)定電阻壁模需要的臨界增益為 | G| = 0.6.可見, 在三組線圈與等離子體旋轉的作用下穩(wěn)定效果最好.

5 結 論

本文研究了在ITER 裝置上電阻壁模的主動控制和被動控制, 給出了ITER 裝置先進運行平衡位形; 用MARS 程序在流體模型下研究了等離子體旋轉對電阻壁模穩(wěn)定的物理機理; 使用磁通-電流反饋控制系統(tǒng), 控制系統(tǒng)包括三組反饋控制線圈和一組極向傳感器, 用MARS 程序掃描了上下兩組線圈穩(wěn)定電阻壁模的最佳組合相位; 研究了單獨反饋控制時, 穩(wěn)定電阻壁模的臨界增益及穩(wěn)定電阻壁模的物理機理; 研究了在等離子體旋轉和反饋控制共同作用下對電阻壁模的控制, 研究表明反饋控制的加入, 不僅可以降低穩(wěn)定電阻壁模的旋轉頻率, 而且能夠更快地達到穩(wěn)定.

本文的研究結果對中國聚變工程試驗堆CFETR 的工程設計和運行具有一定指導意義.

感謝房玉在MARS 程序運行過程中給予的幫助和支持, 感謝劉超博士對本工作的有益討論.

主站蜘蛛池模板: 第九色区aⅴ天堂久久香| 欧美综合区自拍亚洲综合天堂| 99热国产这里只有精品无卡顿"| 国产小视频在线高清播放| 激情网址在线观看| 亚洲另类色| 国产免费精彩视频| 国产探花在线视频| 久久精品人人做人人爽97| 日韩东京热无码人妻| 色噜噜狠狠色综合网图区| 天天色综合4| 成人精品区| 免费 国产 无码久久久| 免费人成视网站在线不卡| 国产一级裸网站| 美女无遮挡被啪啪到高潮免费| 成人一区在线| 中文字幕 欧美日韩| 91网址在线播放| 欧美亚洲国产一区| www.99精品视频在线播放| 国产亚洲男人的天堂在线观看| 无码免费试看| 精品一區二區久久久久久久網站| 欧美另类一区| 青青草国产一区二区三区| 精品久久高清| 久久99精品久久久久纯品| 自拍偷拍欧美日韩| 狠狠色噜噜狠狠狠狠色综合久| 99热国产这里只有精品无卡顿" | 国产精品成人不卡在线观看| 久久综合色视频| 性做久久久久久久免费看| 日韩无码黄色网站| 日本午夜三级| 亚洲欧美成人| 欧美色亚洲| 亚洲欧洲日产国产无码AV| 免费 国产 无码久久久| 特黄日韩免费一区二区三区| 青青网在线国产| 一级毛片免费观看不卡视频| 99re在线免费视频| 日韩无码精品人妻| 99re在线免费视频| 永久毛片在线播| 99爱在线| 精品无码日韩国产不卡av| 中文无码精品A∨在线观看不卡| 亚洲天堂免费| 性视频久久| 91色在线视频| 精品国产成人高清在线| 亚洲AV无码不卡无码| 波多野结衣亚洲一区| 国产精品一区在线麻豆| 91精品国产情侣高潮露脸| 国产乱子精品一区二区在线观看| 成人国内精品久久久久影院| 最新日韩AV网址在线观看| 午夜爽爽视频| 一级毛片中文字幕| 91在线激情在线观看| 欧美色亚洲| 亚洲天堂视频网站| 极品国产在线| 国产一级在线播放| 国产精品污污在线观看网站| 制服丝袜无码每日更新| 日本伊人色综合网| 婷婷综合亚洲| 99久久人妻精品免费二区| 91亚洲精品国产自在现线| 国产精品福利社| 在线国产毛片手机小视频| 国产精品无码影视久久久久久久 | 亚洲日韩欧美在线观看| 精品三级在线| 国产视频你懂得| 免费毛片a|