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

基于Fluent網格變形的旋流器的形狀優化

2016-08-18 06:36:05蔣明虎譚放金淑芹龍桂蘭陳桂芹徐保蕊東北石油大學機械科學與工程學院黑龍江大慶68大慶采油三廠一礦黑龍江大慶山東華宸高壓容器有限公司山東濟南500
化工進展 2016年8期
關鍵詞:變形效率優化

蔣明虎,譚放,金淑芹,龍桂蘭,陳桂芹,徐保蕊(東北石油大學機械科學與工程學院,黑龍江 大慶 68;大慶采油三廠一礦,黑龍江 大慶 6;山東華宸高壓容器有限公司,山東 濟南 500)

研究開發

基于Fluent網格變形的旋流器的形狀優化

蔣明虎1,譚放1,金淑芹2,龍桂蘭2,陳桂芹3,徐保蕊1
(1東北石油大學機械科學與工程學院,黑龍江 大慶 163318;2大慶采油三廠一礦,黑龍江 大慶 163113;3山東華宸高壓容器有限公司,山東 濟南 250101)

目前,旋流器的流體結構優化主要局限于單參數的尺寸優化;大多采用多組網格分別模擬的方法;不僅費時,而且對于微小的尺寸變形常常難以操控。當旋流器的分離效率優化到一定峰值后,進一步提高則比較困難。為解決這一問題,提出了基于網格變形的形狀優化方法。以旋流器的分離效率為目標,采用ANSYS Fluent的網格變形方法,優化了旋流器的外形結構,將旋流器的分離效率由原來的95.69%提高到了99.18%,并對比分析了優化前后的壓降、速度和油相體積分布。研究結果表明:Fluent 網格變形器驅動網格節點的平滑變形,可以同時實現變形區域的結構參數的優化組合,以及優化迭代的自動化和智能化,具有良好的全局搜索能力和較強的魯棒性;Fluent網格變形優化算法能夠縮短CAD建模的時間,避免網格重構;可在單一尺寸優化的基礎上,在一定范圍內提高旋流器的分離效率。

形狀優化;網格變形;分離效率;軸流式旋流器;目標函數

旋流分離器是利用離心加速度進行輕相和重相分離的高效分離設備,具有結構緊湊、分離效率高、節能環保、系統配件少、維修費用低等特點,廣泛應用于石油、環保、化工、礦山、冶金等領域。其中,軸流式旋流器的分離效率更高且流場穩定[1],不易產生循環流和短路流,特別適合于含油率較低的油水混合液分離。

優化旋流器的結構參數、提高分離效率一直是研究人員面臨的重大課題。近年來,研究人員采用CFD方法,對軸流式分離器的結構參數進行了一系列優化研究工作,并且取得了一定的研究成果。聶濤[2]、俞接成[3]等初步探討了軸流式旋流器的速度、壓力和分離效率;GONG等[4]模擬了不同螺旋角和導流葉片的軸流式分離器的分離性能;王云峰等[5]對實體錐段和多孔介質錐段的軸流式旋流器進行了模擬研究;宋民航等[6]優化了軸流式脫水型旋流器的溢流口結構和出口形式;蔣明虎等[7-9]模擬了軸流式分離器的速度分布和壓降規律并進行了實驗研究,探討了含氣量對分離器性能的影響,優選了氣、液、砂三相分離器的處理量和分流比等操作參數;趙立新等[10-11]優選了軸流式油水分離器的入口結構以及處理量和分流比等操作參數;史仕熒等[12]優化了軸流式分離器的矩形切向開縫出油管的結構;黃龍等[13]對軸流式氣液旋流器的速度分布規律與壓降變化進行了數值模擬。以上文獻均采用的是基于控制變量法的單一尺寸優化;運用該方法開展旋流器的結構優化研究已進行得比較深入,對于進一步提高旋流器的分離效率已幾乎沒有提升空間;況且,單一的尺寸優化對旋流器的局部外形難以準確控制。為此,采用ANSYS Fluent的網格變形的方法[14],以旋流器的分離效率為單一目標函數,對旋流器的流體結構進行了形狀優化,進一步提高了旋流器的分離效率。

1 ANSYS Fluent的網格變形優化器

所謂優化就是在設計空間中尋找最佳的合理設計點。Fluent網格變形器(mesh morpher)恰是基于用戶自定義的控制點的網格運動進行光滑地變形的[15];在無需過多的人工干預的情況下,智能地尋找最佳設計點,可以很好地解決形狀優化問題[16-17]。ANSYS Fluent在用戶設定的變形區域中,設計點是由“控制體”限定的一系列控制點來實現,由這些控制點操控指定區域(控制體)網格的光滑變形;依據優化器的不同,其變形原理基于徑向基函數插值法[18]或張量積的伯恩斯坦多項式[19]。

Fluent網格變形器無需重新修改CAD模型和網格重構,節省了CAD建模和劃分網格的時間,對網格類型和物理模型沒有限制;不足之處是無法對重疊的區域同時進行變形,不支持任意形狀的變形,不能單獨用于多目標優化。

Fluent Mesh Morpher/Optimizer(MMO)一共內置6個優化器,分別是:The Compass Optimizer[20]、The NEWUOA Optimizer[21]、 The Simplex Optimizer[22]、The Torczon Optimizer[23]、The Powell Optimizer[24]和The Rosenbrock Optimizer[25]。其中,Compass優化器以目標函數的最小值為網格變形的調整方向,從一個給定的值開始,基于最小二乘法的原理[26],通過逐步改進控制點進行網格變形。本文主要采用Compass優化器進行優化。

2 模型尺寸及參數設置

2.1旋流分離原理

模擬采用的軸流式旋流分離器是經尺寸優化過的,其幾何結構模型如圖1所示;初始結構的主要參數:入口段長度L1=50mm,螺旋段長度L2=57mm,旋流腔長度L3=80mm,錐管段長度L4=336mm,尾管段長度L5=50mm,溢流管長度L7=16mm,旋流腔直徑D2=50mm,底流口直徑D4=25mm,溢流管直徑D5=8mm,螺旋流道數N=5,錐管段錐角α=6°。其工作原理是,油水混合介質由入口進入螺旋流道產生強旋流,沖入旋流腔后繼續作強螺旋流動,產生重力加速度幾百倍以上的離心加速度。由于密度的差異,在離心力的作用下,水相流向壁面并匯集到底流口流出,而油相則被擠向旋流器的中心區域,形成細長狀的油心,在背壓的作用下從溢流口流出,從而實現油水分離。

圖1 旋流器初始結構模型

2.2網格劃分

因六面體結構化網格數量較少,網格單元和流體流動方向對齊,既可以提高求解精度,又可以節約計算成本,特別適合于多相流的模擬計算;故采用ICEM CFD對旋流器進行全六面體網格劃分。為檢驗網格的無關性,劃分了3種不同數量的六面體網格單元;分別是390901個(分離效率為95.56%)、445188個(分離效率為95.69%)和852067個(分離效率為95.71%),網格質量均在0.3以上;因后兩個分離效率模擬數據比較接近,考慮到節省計算資源,最終確定為 445188個六面體單元,413203個節點(如圖2所示)。

圖2 旋流器的網格劃分

2.3參數設置

影響油水分離器分離性能的參數主要有:入口混合介質參數(油滴粒徑、油和水的黏度、油水密度差和含油濃度等)、操作參數(入口流量、分流比等)和結構參數(旋流腔尺寸、錐管段尺寸、尾管段尺寸、入口尺寸、溢流管尺寸和螺旋流道尺寸等)。本文只改變旋流腔、錐管段及尾管段的形狀和尺寸進行數值模擬,力求得到三者的最佳匹配參數。

油水混合介質的物性參數:水的密度為0.9982×103kg/m3,水的黏度為 1.003×103kg/(m·s),油的密度為889kg/m3,油的黏度為1.06kg/(m·s),油滴粒徑為0.4mm。操作參數:環境溫度為298K,環境壓力為101325N/m2,入口速度為0.8m/s,底流分流比為80%,溢流分流比為20%,水相的體積分數為98%,油相的體積分數為2%。

就旋流器模擬方法而言,大渦模擬(large eddy simulation,LES)和雷諾應力(Reynolds stress model,RSM)模型最適合三維高旋流,但大渦模擬需要較高的網格質量且耗費的計算資源較多,而雷諾應力模型忽略了各項同性的渦黏性假設,較好地考慮了漩渦、旋轉流動的變化情況,具有較強的模擬螺旋復雜湍流的能力,計算結果比較可信;又因為Mixture多相流模型能夠考慮相間的速度滑移,故采用RSM模型和Mixture模型進行數值模擬。

模擬采用的是基于壓力基的計算穩健的SIMPLE算法、雙精度和二階迎風格式;因計算域具有較大的曲率且為高旋流,故壓力項設置為PRESTO??;殘差精度設置為10-7,完全可以滿足計算要求。為方便配置分流比,底流和溢流出口的邊界條件均設置為 outflow。數值模擬設置如表 1所示。

表1 數值模擬設置

3 優化過程

3.1變形原理[27]

用于本文的網格變形原理基于伯恩斯坦多項式。若在一個坐標系內建立一個控制體,控制體內的坐標控制點被映射到網格節點上,由坐標控制點帶動網格節點的運動。

假定q為網格內節點,其位移由局部三維控制點(uq,vq,wq)進行限定,在 l×m×n的控制體內,u、v、w分別與坐標 i、j、k方向相對應。第 q個網格節點位置與第i、j、k控制點位置之間的線性關系見式(1)。

式(1)中,Bi,l(u)為第l階伯恩斯坦多項式,見式(2)。)

控制體內的網格變形是由一個或多個控制點的移動映射到網格節點來進行平滑定位的,由網格節點定位的目標函數δτ的變化可按式(3)、式(4)估算。

為控制點的敏感區域,δζijk表示控制點位置的調整。通過簡單的梯度運算,令δζijk=λWi jk,即可求出光滑邊界的網格形變,從而改進設計(式中比例系數λ默認值為1);增加網格數量和單元階次都可以提高提高計算精度。

3.2變形優化過程

Fluent MMO的運行步驟如下:

(1)按照Fluent設置首先求得初始結構的模擬結果文件;

(2)加載MMO模塊,依次設置變形體區域、網格變形方向和大小、目標函數、優化器、收斂精度及離散步數等;

(3)優化求解,Fluent軟件會根據設定的精度自動智能開始優化計算,到達設計精度時運算終止。

因旋流腔、錐管段和尾管段(如圖1所示)是影響旋流器的分離效率的主要結構參數[28];因此,在保持旋流器整體長度不變的情況下,取8個優化設計點,對這三部分進行結構變形優化,變形方向均為徑向中心方向。網格變形控制體如圖3所示。

網格變形體共得到49個,從中選取4個對稱變形結構,如圖4所示,可看出大致變形趨勢。

圖3 網格的變形控制體

圖4 旋流器的網格變形過程

定義旋流器分離效率為式(5)。

式中,Ez為旋流器分離效率;m˙in為入口油相的質量流量,kg/s;m˙ov為溢流油相的質量流量,kg/s。以分離效率作為單一目標函數(如圖 5所示),收斂精度設置為 0.001,由于溢流口的質量流量為負值,故優化目標就是力求目標函數值為最小。

目標函數的優化設計階段如圖6所示。圖6中,目標函數曲線左端點為目標函數(分離效率)的初始值,右端點為最終優化值;在經過進行49次網格變形后達到所設定的計算精度,分離效率達到最大并趨于穩定。目標函數值越小分離效率越高,即曲線越靠近 x 軸分離效率越高。經計算,優化前旋流器的分離效率為95.69%,優化后的分離效率為99.18%,分離效率提高幅度較為明顯。

4 結果分析

4.1形狀變化

圖7中(a)和(b)為變形前后的幾何形狀變化圖,(c)為變形前后幾何形狀重疊圖,透明部分為變形前的形狀。從圖7中可以看出,變形后得到的幾何結構(旋流腔、錐管段和尾管段)較原來纖細,旋流腔和尾管段均由圓柱形變成圓臺形,錐管段的錐角變小了。表2為幾何變形前后各段錐角的變化情況。

圖5 目標函數的設置

圖6 目標函數的優化設計階段

圖7 變形前后的形狀比較(單位:m)

表2 變形前后的錐角    單位:(°)

從以上數據可知,旋流器的幾何形狀的變化不大,優化后的旋流腔和錐管段平滑過渡并融為一體。這種幾何形狀的微小變化是傳統的單一尺寸優化難以達到的。

4.2壓降變化

由Fluent求解入口面、出口面及溢流口面的平均靜壓力及壓降如表3所示。

表3 平均壓降變化    單位:Pa

由表3可知,變形后的壓降均比變形前有所增大,但變化幅度均較小。變形前溢流口的壓力為負值,對于溢流液的流出極為有利;由于流體壓強的相互傳遞,變形后入口的壓力增加了,溢流口的壓力為雖為正值,但變化較小,約為7196.45Pa,對溢流液的流出影響不大。變形后的出口壓降較變形前增加了,約為2720.85Pa,對于分離液的流出影響更小。圖8為錐管段Z= -100mm處的變形前后靜壓力分布圖,可知變形后混合液的壓力有一定的損失,大約在104Pa數量級,約為0.1個大氣壓。

4.3速度變化

圖8 變形前后的Z= -100mm處的壓力

圖9為錐管段截面Z= -100mm處的變形前后合速度對比圖,可以看出,變形前后的圖形基本形狀未變,均為M形。變形后,外旋流的速度變大,內旋流的速度反而變小,二者的速度差有利于油水混合液的分離。

圖9 變形前后的Z= -100mm處的合速度

切向速度決定旋流器中離心加速度和離心力的大小,是決定旋流器分離性能的重要因素[29]。圖10為變形前后Z= -100mm處切向速度的變化圖。由圖10可知,變形前后的切向速度形狀相似,且速度最大值重合;但變形后外旋流的速度變小且最值的差值較大,比變形前更靠近壁面;這種變化可促進油水分離,提高分離效率。

圖10 變形前后的Z= -100mm處的切向速度

綜合壓降和速度變化可知,由旋流器幾何形狀的纖細改變而增加的旋流速度,是由增大的壓降轉化而來的;優化后的壓力能比優化前降低較多,補充了油水分離所需要的動能,符合能量守恒定律,正是這種能量的轉化提高了旋流器的分離效率。

4.4油相體積的變化

圖11為油相體積的軸截面圖。從圖11中可以看出,在軸向中心和溢流口處,變形后的油心長度和寬度均比變形前的大,這說明油相的分布與變形前相比具有較大的體積分數,分離效率較變形前有所提高。

總之,經尺寸優化過的流體結構,再進行形狀優化,分離效率仍有一定提升空間。分別改變網格的最小正交質量或選擇不同的優化器(前文已提及,MMO共包含6個優化器),模擬結果是相同的;這說明MMO是穩健可靠的。

圖11 變形前后油相的體積分數(單位:m)

5 模擬與實驗對比分析

因入口流量為入口速度與入口截面積之積,為方便在實驗[30]中控制流量,把入口速度0.7m/s、0.8m/s、0.9m/s、1.0m/s、1.1m/s換算為對應的進口流量,分別為5.0m3/h、5.7m3/h、6.4m3/h、7.1m3/h、7.8m3/h。

5.1分離效率的對比

圖12為溢流分流比為0.2時,分離效率隨入口流量變化的模擬與實驗值對比圖。由圖12可知,隨著流量的增大,分離效率有逐漸增加的趨勢,但流量越大,分離效率增加的越緩慢且實驗值均低于模擬值;當流量約為6.77m3/h時,變形后分離效率的實驗值最高,超過7.1m3/h時,分離效率的實驗值逐漸下降。這是因為實際分離過程中,隨著流量的增加,剪切應力也隨之增大,導致一部分油滴破碎,增大了油滴的乳化作用,致使分離效率降低。

圖12 分離效率對比圖

5.2壓降的對比

圖 13為變形后壓降隨入口流量變化的趨勢對比圖。

由圖13可知,隨著入口流量的增大,模擬和實驗結果的壓降也隨之增大,溢流壓降值均大于出口壓降值;實驗值均大于模擬值,這主要是旋流器實體結構制造質量的影響、模擬時對旋流器的壁面進行近似處理所致。

圖13 壓降對比圖

6 結 論

(1)Fluent網格變形技術,可以在保證流體網格質量的前提下實現幾何結構的平滑變形,避免受到幾何重新建模和網格重新劃分的制約。

(2)Fluent網格變形可以同時對控制體內的幾何參數進行智能匹配,從全局角度系統地對結構參數進行目標驅動優化。優化過程是自動智能的,無需過多人工干預,具有很好的全局搜索能力。

(3)提高分離效率所需要的旋流速度,是由進口、出口和溢流口之間的壓降來補償的;壓降的轉化量是提高分離效率的基本前提。

[1]趙立新,宋民航,蔣明虎,等.軸流式旋流分離器研究進展[J].化工機械,2014,41(1):20-25.

[2]聶濤. 軸流式液液旋流器內流場的數值模擬[D]. 東營:中國石油大學(華東),2008:47-51.

[3]俞接成,陳家慶,韓景. 軸向入口油水分離水力旋流器及其數值模擬[J]. 北京石油化工學院學報,2009,17(2):19-23.

[4]GONG Guangcai,YANG Zhouzhou,ZHU Shaolin. Numerical investigation of the effect of helix angle and leaf margin on the flow pattern and the performance of the axial flow cyclone separator[J]. Applied Mathematical Modelling,2012,36(8):3916-3930.

[5]王云峰,劉仁桓,王金花,等. 2種錐段結構的軸流式旋流器內的流場模擬[J]. 流體機械,2013,41(6):22-26.

[6]宋民航,韓佳軒. 軸流式脫水型旋流器的流場分析與結構優化[J].石油礦場機械,2012,41(12):56-60.

[7]蔣明虎,陳世琢,李楓,等. 緊湊型軸流式除油旋流器模擬分析與實驗研究[J]. 油氣田地面工程,2010,29(9):18-20.

[8]蔣明虎,宮磊磊,徐保蕊,等. 含氣條件對井下油水分離旋流器性能影響的數值模擬[J]. 化工機械,2014,41(5):629-632.

[9]蔣明虎,李永山,徐保蕊,等. 軸流式脫氣除砂三相旋流分離器操作參數優選[J]. 化工機械,2015,42(1):68-71.

[10]趙立新,代佳鑫,郭現臣. 葉片式水力旋流器操作參數優選[J]. 流體機械,2013,41(10):7-9.

[11]趙立新,宋民航,蔣明虎,等. 新型軸入式脫水型旋流器的入口結構模擬分析[J]. 石油機械,2013,41(1):68-71.

[12]史仕熒,吳應湘,許晶禹,等. 導流片型管道式分離器油水分離結構優化[C]//第十三屆全國水動力學學術會議暨第二十六屆全國水動力學研討會,青島,2014:1063-1067.

[13]黃龍,鄧松圣,王斌,等. 新型軸流導葉式氣液旋流器流場數值模擬[J]. 機床與液壓,2015,43(1):164-167.

[14]楊磊,李寶軍,胡平. 網格變體方法的工程應用與進展綜述[C]//中國計算力學大會2014暨第三屆錢令希計算力學獎頒獎大會,貴陽,2014.

[15]BIANCOLINI Marco Evangelos. Mesh morphing accelerates design optimization[E/OL]. Rome,Italy,2010. http://www.ansys.com/ staticassets/ANSYS/staticassets/resourcelibrary/article/AA-V4-I1-Mesh-Morphing-Accelerates-Design-Optimization.pdf.

[16]JEMCOV Aleksandar,MARUSZEWSKI Joseph. Shape optimization based on downhill simplex optimizer and free-form deformation in general purpose CFD code[C]// 17th Annual Conference of CFD Society of Canada,Ottawa,Ontario,Canada,2009.

[17]MURUGAN S,WOODS B K S,FRISWELL M I. Hierarchical modeling and optimization of camber morphing airfoil[J]. Aerospace Science and Technology,2015,42(2):31-38.

[18]BIANCOLINI Marco Evangelos. Advanced mesh morphing for automotive applications using RBF Morph[C]//Automotive Simulation World Congress,Frankfurt am Main,Germany,2013.

[19]RAIDL Günther R,WUTM Christian. Approximation with evolutionary optimized tensor product Bernstein Polynomials[A]. Karlsplatz 13,1040 Vienna,Austria,2003.

[20]KOLDA Tamara G,LEWIS Robert Michael,TORCZON Virginia. Optimization by direct search:new perspectives on some classical and modern methods[J]. SIAM Review,2003,45(3):385-482.

[21]ZASLAVSKI Alexander,POWELL M J D. The NEWUOA software for unconstrained optimization without derivatives[M]. US:Springer,2004:255-297.

[22]MCKINNON K I M. Convergence of the Nelder-Mead simplex method to a nonstationary point[J]. SIAM J. Optimization,1998,9 (1):148-158.

[23]TORCZON Virginia. On the convergence of the multidirectional search algorithm [J]. SIAM J. Optimization,1991,1(1):123-145.

[24]PRESS William H,TEUKOLSKY Saul A,VETTERLING William T. Numerical recipes:the art of scientific computing[M]. 3rd Ed. Cambridge,England:Cambridge University Press,2007.

[25]ROSENBROCK H H. An automatic method for finding the greatest or least value of a function[J]. Computer Journal,1960,3(3):175-184.

[26]王仁芳,許秋兒,汪沁,等. 基于最小二乘網格的模型變形算法[J].計算機輔助設計與圖形學學報,2010,22(5):777-783.

[27]ELSAYED Khairy. Design of a novel gas cyclone vortex finder using the adjoint method[J]. Separation and Purification Technology,2015,142:274-284.

[28]OLSONT T J,VAN OMMEN R. Optimizing hydrocyclone design using advanced CFD model[J]. Minerals Engineering,2004,17(5):713-720.

[29]袁惠新,李雙雙,付雙成,等. 三相分離旋流器內流場及分離性能的研究[J]. 流體機械,2015,43(1):28-32.

[30]盛慶嬌. 新型螺旋入口水力旋流器模擬分析及實驗研究[D]. 大慶:東北石油大學,2012:39-45.

Sharp optimization of hydrocyclone based on Fluent mesh morphing

JIANG Minghu1,TAN Fang1,JIN Shuqin2,LONG Guilan2,CHEN Guiqin3,XU Baorui1
(1School of Mechanical Science and Engineering,Northeast Petroleum University,Daqing 163318,Heilongjiang,China;2The First Oil Mine,Daqing No.3 Oil Production Company,Daqing 163113,Heilongjiang,China;
3Shandong Huachen High Pressure Vessel Co.,Ltd.,Ji'nan 250101,Shandong,China)

Currently,fluid-structure optimization of hydrocyclone has been mainly limited to single-parameter size optimization. Generally,the optimization process is dependent upon the method of separate simulation with multiple sets of grid. Nevertheless,the approach is time-consuming and not capable of controlling the micro size deformation. Besides,further improvement on the separation efficiency of cyclone is rather difficult once a certain maximum optimization is achieved. To solve this problem,mesh deformation-based shape optimization method was proposed. Aimed at improving the hydrocyclone separation efficiency with ANSYS Fluent mesh deformation method,the hydrocyclone sharp structure was optimized and the hydrocyclone separation efficiency was increased from 95.69% to 99.18%. A comparative analysis involved with pressure drop,velocity and oil phase volume distributions before and after optimization was then conducted. The results indicated that the smoothing mesh deformation driven grid nodes by Fluent mesh morpher could concurrently contribute to the optimal combination of structural parameters in the deformation zone,and the optimization of the automation and intelligentization of the iteration process. Besides,it had a considerable capability ofglobal searching and strong robustness. Fluent mesh morphing optimization algorithm can shorten CAD modeling time,avoid mesh reconstruction and further improve the separation efficiency of hydrocyclone within a certain range on the basis of single-parameter size optimization.

sharp optimization; mesh morphing;separation efficiency;axial cyclone;objective function

TQ 028.4

A

1000-6613(2016)08-2355-07

10.16085/j.issn.1000-6613.2016.08.08

2016-02-18;修改稿日期:2016-03-11。

國家高技術研究發展計劃(2012AA061303)。

及聯系人:蔣明虎(1962—),男,教授,博士生導師,主要研究方向為旋流分離技術。E-mail jmhdq@126.com。

猜你喜歡
變形效率優化
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
談詩的變形
中華詩詞(2020年1期)2020-09-21 09:24:52
提升朗讀教學效率的幾點思考
甘肅教育(2020年14期)2020-09-11 07:57:42
“我”的變形計
例談拼圖與整式變形
會變形的餅
跟蹤導練(一)2
主站蜘蛛池模板: 久久大香香蕉国产免费网站| 亚洲最新在线| 免费欧美一级| 亚洲精品日产AⅤ| 91视频99| 国产三级国产精品国产普男人| 亚洲免费福利视频| 久久99国产乱子伦精品免| 久久一日本道色综合久久| 亚洲免费人成影院| 538精品在线观看| 婷婷亚洲最大| 亚洲av日韩av制服丝袜| 999精品在线视频| 国产杨幂丝袜av在线播放| 久久人妻xunleige无码| 国产熟女一级毛片| 国产午夜不卡| 国产成人久视频免费| 国产视频大全| 免费无码AV片在线观看国产| 国产精品福利尤物youwu| 国产麻豆永久视频| 欧美成人一级| 欧美日韩午夜| 亚洲精品桃花岛av在线| 国产91视频免费观看| 亚洲 欧美 中文 AⅤ在线视频| 人妻丝袜无码视频| www.91在线播放| 五月婷婷综合网| 日韩天堂在线观看| 色综合久久久久8天国| 99精品欧美一区| 欧美激情视频二区| 波多野结衣在线se| 2021国产精品自产拍在线| 欧美另类视频一区二区三区| 亚洲欧美另类日本| 萌白酱国产一区二区| 欧美第一页在线| 国产国语一级毛片| 成AV人片一区二区三区久久| a毛片基地免费大全| 欧美黑人欧美精品刺激| 国产一区二区网站| 成年人福利视频| 欧美午夜网站| 97se亚洲| 国产精品自拍合集| 香蕉视频国产精品人| 一级毛片免费的| 91亚洲精品第一| 久久先锋资源| 26uuu国产精品视频| 伊人久综合| 国产亚洲日韩av在线| 四虎成人在线视频| 国内嫩模私拍精品视频| 久久性妇女精品免费| 国产性精品| 久久成人18免费| 老色鬼欧美精品| 亚洲无码不卡网| 欧美在线观看不卡| 亚洲国产精品日韩欧美一区| 亚洲人成影院在线观看| 亚洲午夜福利精品无码| 国产无遮挡猛进猛出免费软件| 亚洲综合狠狠| 97久久超碰极品视觉盛宴| 国产精品开放后亚洲| 一级毛片在线播放| 在线99视频| 九九九国产| 老熟妇喷水一区二区三区| 欧美日韩在线成人| 中文字幕在线观| 亚洲第一精品福利| 欧美一级夜夜爽| 成人在线不卡视频| 亚洲中文在线视频|