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

注塑成型有限元分析管理和成型工藝優(yōu)化系統(tǒng)

2018-01-10 04:17:24凌成智胡廣洪
精密成形工程 2018年1期
關(guān)鍵詞:有限元工藝優(yōu)化

凌成智,胡廣洪

(上海交通大學(xué) 塑性成形技術(shù)與裝備研究院,上海 200030)

注塑成型有限元分析管理和成型工藝優(yōu)化系統(tǒng)

凌成智,胡廣洪

(上海交通大學(xué) 塑性成形技術(shù)與裝備研究院,上海 200030)

目的 針對注塑成型工藝過程,開發(fā)有限元分析項目流程管理和工藝參數(shù)優(yōu)化的集成系統(tǒng)。方法 參考企業(yè)項目流程,開發(fā)項目流程管理模塊,使用多元回歸擬合配合改進(jìn)的模擬退火算法進(jìn)行工藝參數(shù)的優(yōu)化,通過塑料托盤翹曲變形量的優(yōu)化來驗證優(yōu)化方法的準(zhǔn)確性,進(jìn)行塑料托盤的CAE分析正交試驗,使用多元回歸擬合試驗數(shù)據(jù)并進(jìn)行擬合優(yōu)度檢驗,分別使用Isight和改進(jìn)的模擬退火算法進(jìn)行最優(yōu)解的求解并對比計算結(jié)果,最后對得到的最優(yōu)參數(shù)組合進(jìn)行CAE分析得到變形量,將變形量與正交試驗變形量對比。結(jié)果 在回歸方法中,使用懲罰系數(shù)為2的多項式嶺回歸效果較好,其優(yōu)度指標(biāo)R2為0.975,改進(jìn)的模擬退火算法計算結(jié)果與Isight軟件自適應(yīng)模擬退火算法的計算結(jié)果基本一致,最優(yōu)參數(shù)組合在Moldflow分析中的變形量為0.6805 mm,與之前的最小變形量0.7436 mm相比更小。結(jié)論 使用回歸擬合配合改進(jìn)的模擬退火算法開發(fā)的工藝參數(shù)優(yōu)化程序能起到較好的工藝參數(shù)優(yōu)化效果。

注塑成型;有限元分析;項目管理;工藝優(yōu)化

目前在注塑工藝設(shè)計過程中,使用CAE軟件進(jìn)行仿真分析已經(jīng)成為其中至關(guān)重要的一環(huán),采用有限元計算方法,模擬整個注塑過程,并分析成型工藝對注塑成型產(chǎn)品質(zhì)量的影響,以幫助設(shè)計人員及早發(fā)現(xiàn)并解決問題[1],因而,CAE軟件在提高研發(fā)效率,降低研發(fā)成本等方面發(fā)揮了重要作用。在處理大量的重復(fù)性的CAE計算以及成型工藝參數(shù)的優(yōu)化問題上,使用單一 CAE軟件就顯得效率低下了,因此,CAE集成系統(tǒng)應(yīng)運而生,成為目前 CAE技術(shù)發(fā)展的一個重要方向。國內(nèi)外的許多軟件公司、高校和研究機構(gòu)著力于研發(fā)圍繞CAE軟件的集成系統(tǒng)或平臺[2—3]。

文中描述的是針對注塑成型進(jìn)行有限元分析,開發(fā)了 CAE項目流程管理和成型工藝參數(shù)優(yōu)化的集成系統(tǒng),實現(xiàn)對注塑成型有限元分析項目流程進(jìn)行科學(xué)、系統(tǒng)的管理,將人力從具有繁重的重復(fù)性的有限元分析任務(wù)中解放出來,同時,為注塑成型工藝參數(shù)的優(yōu)化提供了高效、精確的研發(fā)工具。

1 軟件架構(gòu)設(shè)計

本系統(tǒng)整體上使用客戶端/服務(wù)器(Client/Server)結(jié)構(gòu)。客戶端為本軟件主體,服務(wù)器則為存放各類數(shù)據(jù)的數(shù)據(jù)庫。本系統(tǒng)主要有用戶管理、項目流程管理、集成研發(fā)工具3個功能模塊。3個模塊以及相應(yīng)的數(shù)據(jù)庫都相互獨立,方便根據(jù)需求進(jìn)行擴展和對數(shù)據(jù)的管理、維護(hù)。

本系統(tǒng)主要使用C#語言,在.NET平臺上進(jìn)行開發(fā)。前2個模塊通過在.NET上對關(guān)系型數(shù)據(jù)庫進(jìn)行操作的方式,實現(xiàn)相應(yīng)的功能。集成研發(fā)工具模塊中,使用讀取數(shù)據(jù)庫正交表的方式實現(xiàn)工藝參數(shù)的設(shè)置;使用VB腳本調(diào)用Moldflow SynergyAPI來驅(qū)動CAE軟件,然后在程序里重開進(jìn)程,運行VB腳本完成自動運行 Moldflow軟件的功能;調(diào)用數(shù)據(jù)計算的Python腳本完成工藝參數(shù)優(yōu)化功能。

2 用戶和項目流程管理

有限元分析是現(xiàn)代企業(yè)成型工藝研發(fā)的一個必不可少的步驟,在整個研發(fā)過程中與其他環(huán)節(jié)有緊密的聯(lián)系,同時,人員上也存在不同工程師之間的協(xié)調(diào),因此,在項目的開展過程中,對人員和任務(wù)進(jìn)行清晰、科學(xué)的管理很有必要。在用戶管理模塊中,對所有的用戶都進(jìn)行了權(quán)限管理,每個工程師分工明確,在登陸系統(tǒng)之后只能調(diào)用與權(quán)限對應(yīng)的功能模塊。本系統(tǒng)將人員分為4個類別:超級管理員、總工程師、高級工程師、工程師。超級管理員主要負(fù)責(zé)對用戶信息進(jìn)行管理,如用戶注冊審批、用戶信息修改、用戶賬號銷毀等,注冊管理界面見圖1;總工程主要負(fù)責(zé)項目的新建和撤銷,在新建項目后需要將項目細(xì)分為具體的任務(wù);高級工程師的工作是將項目的具體任務(wù)分配給具體工程師,審核工程師提交的任務(wù),項目完成時進(jìn)行提交等,任務(wù)分配界面見圖2;工程師執(zhí)行通過“我的任務(wù)”菜單欄查看自己任務(wù),執(zhí)行并提交任務(wù)。

圖1 用戶注冊管理界面Fig.1 Interface of user registration management

圖2 任務(wù)分配界面Fig.2 Interface of task assignment

3 注塑成型工藝集成研發(fā)工具

3.1 基于Moldflow API的CAE分析任務(wù)自動運行程序設(shè)計

在本系統(tǒng)中,首先通過正交試驗的方法設(shè)置成型工藝參數(shù),然后調(diào)用Synergy API對Moldflow分析軟件進(jìn)行驅(qū)動,實現(xiàn)復(fù)制模型,設(shè)置參數(shù),同步運行以及結(jié)果數(shù)據(jù)的自動讀取,從而極大提高了CAE分析的工作效率和后期數(shù)據(jù)處理的準(zhǔn)確性。整個流程見圖3。

圖3 CAE分析任務(wù)自動與運行程序流程Fig.3 Work flow of automatic operation analysis of CAE

所有的實驗及數(shù)據(jù)都是以塑料托盤為例,分析內(nèi)容包括冷卻、填充、保壓、翹曲,材料為PP(POLYFLAM RPP 374 ND CS1),具體模型見圖4。

圖4 塑料托盤模型Fig.4 Model of plastic tray

驅(qū)動MOLDFLOW軟件進(jìn)行批量CAE運算的流程為:首先將系統(tǒng)計算生成的實驗數(shù)據(jù)格式化為參數(shù)字符串;然后在系統(tǒng)中重開進(jìn)程運行VB腳本,將前一步生成的字符串傳入腳本作為參數(shù);最后調(diào)用Synergy API提供的相應(yīng)接口依次循環(huán)執(zhí)行打開模型、設(shè)置參數(shù)、提交分析和復(fù)制模型等操作。VB腳本中程序運行邏輯具體情況見圖5。

圖5 批量運行CAE分析程序Fig.5 Block diagram for batch run of CAE analysis

CAE軟件分析結(jié)果數(shù)據(jù)的獲取,同樣也是通過VB腳本調(diào)用Synergy API進(jìn)行實現(xiàn)的。本系統(tǒng)可以獲取各目標(biāo)參數(shù)的最大、最小及平均值。圖6顯示了部分分析數(shù)據(jù)的自動獲取情況。

圖6 分析結(jié)果獲取界面Fig.6 Interface of retrieving analysis results

3.2 基于回歸分析和改良模擬退火算法(SA)的工藝優(yōu)化程序

成型工藝參數(shù)優(yōu)化程序是整個系統(tǒng)的核心。本系統(tǒng)的工藝優(yōu)化流程是:按照正交試驗設(shè)計程序,設(shè)計關(guān)于工藝參數(shù)的正交試驗,通過自動化驅(qū)動程序驅(qū)動CAE軟件(Moldflow)進(jìn)行計算,提取得到結(jié)果數(shù)據(jù),將正交試驗參數(shù)組和結(jié)果數(shù)據(jù)的某一特定分析項一起作為樣本數(shù)據(jù),對這些樣本數(shù)據(jù)進(jìn)行多元回歸擬合,最后對擬合的結(jié)果使用最優(yōu)解搜尋算法進(jìn)行最值求解。得到的結(jié)果即為最優(yōu)的工藝參數(shù)組合。

回歸過程使用的回歸方法為帶有二次項的多項式回歸模型。其本質(zhì)上是對多元線性回歸的一個擴展。首先介紹一下簡單的多元線性回歸,其回歸方程的一般形式為:

式中:β1,β2…βp為p個待求的回歸系數(shù);β0為常數(shù)項系數(shù);ε為隨機誤差;y為因變量,即目標(biāo)函數(shù);xi為第i個自變量。文中使用回歸方法為scikit-learn工具中的嶺回歸[4],這種回歸方法是在普通的“最小二乘法”的殘差平方和項中加入一個懲罰項(α|β|2)來解決共線性和病態(tài)數(shù)據(jù)的問題[4—6]。其中α為懲罰系數(shù),α值的確定是由以α為自變量,擬合效果指標(biāo)Score為因變量的函數(shù)的最大值處確定,見圖7。

圖7 懲罰系數(shù)與擬合效果之間的關(guān)系Fig.7 Relationship between penalty factor and fitting effect

由此,本系統(tǒng)使用擴展參數(shù)組的方式來將其擴展為多項式回歸[7—9]。當(dāng)只有想x1,x2兩個參數(shù)時,簡單線性擬合時所求的系數(shù)有常數(shù)項、x1的系數(shù)、x2的系數(shù),而多項式回歸時則需計算x12,x22,x1x2值,并將其加入?yún)?shù),所求系數(shù)增加x12的系數(shù)、x22的系數(shù)、x1x2的系數(shù)。這樣就可以求出帶有二次項的回歸結(jié)果,見式(2)。

文中使用改進(jìn)的模擬退火算法求解最優(yōu)參數(shù)組合,其算法流程見圖8。改進(jìn)之處包括。

1) 初始值的選取:與一般的隨機取初始值不同,文中使用的初始值為田口方法的極差分析得到的局部最優(yōu)組合。該初始值在正交試驗的水平間隔上已經(jīng)是最優(yōu)解,在該初始值基礎(chǔ)上開始搜尋,以提高效率。

2) 新解的產(chǎn)生方式:通常新解產(chǎn)生方式是將固定步長加在原始解的一個隨機參數(shù)上形成新解。文中算法則采用與溫度相關(guān)的步長,在降溫過程中步長也隨之變小,以增加算法的局部搜索能力和降低丟失最優(yōu)解的概率[10]。

3) 新的降溫規(guī)則:配合新的初始值選取方式(見式(3)),文中采用較低的初始溫度以及緩慢的降溫方式。(θ為降溫系數(shù),文中的算法中取0.9)。

4) 記錄歷史最優(yōu)解,在不增加時間、空間復(fù)雜度的情況下保證算法結(jié)果優(yōu)于或等于算法結(jié)束時的當(dāng)前解。

圖8 改進(jìn)的模擬退火算法Fig.8 Improved simulated annealing algorithm

4 塑件變形量的優(yōu)化分析

本節(jié)以制件翹曲變形量為目標(biāo)函數(shù),利用本系統(tǒng)進(jìn)行優(yōu)化分析,從而降低制件的最終變形量[11]。首先在本系統(tǒng)下進(jìn)行2組五因素五水平的正交試驗[12],兩組試驗一共50次CAE分析,將第一組的數(shù)據(jù)用于擬合各工藝參數(shù)與變形量最大值之間的關(guān)系,第二組數(shù)據(jù)用于驗證該擬合結(jié)果的效果。最后對本系統(tǒng)得出的最優(yōu)參數(shù)組合進(jìn)行CAE分析實驗驗證,其中設(shè)置的參數(shù)及其水平區(qū)間見表1。

表1 注塑成型CAE分析參數(shù)表Tab.1 CAE analysis parameters of injection molding

通過本系統(tǒng)進(jìn)行正交試驗的設(shè)計,得到的參數(shù)表界面見圖 9。調(diào)用 Moldflow計算后提取的結(jié)果界面見圖6,Moldflow分析結(jié)果第一組的最大變形量分別為4.52, 3.25, 1.98, 1.01, 0.74, 2.78, 1.56, 0.77, 2.17,3.77, 1, 2.96, 1.69, 3.26, 1.12, 2.5, 0.9, 1.85, 1.18, 3.08,0.95, 1.55, 4.01, 1.71, 0.95 mm;第二組變形量分別為4.78, 3.56, 2.32, 1.15, 0.77, 3.04, 1.9, 0.81, 2.67, 4.1,1.06, 3.37, 2.08, 3.69, 1.22, 2.87, 0.98, 2.07, 1.31, 3.49,1.03, 1.69, 4.33, 2.25, 1.09 mm。

圖9 正交試驗參數(shù)Fig.9 Orthogonal experimental parameters

根據(jù)分析結(jié)果,本系統(tǒng)對第一組實驗數(shù)據(jù)分別使用普通最小二乘法線性回歸、普通最小二乘法多項式回歸和懲罰系數(shù)為2的多項式嶺回歸擬合,將第二組的正交試驗參數(shù)代入第一組數(shù)據(jù)的各個擬合結(jié)果進(jìn)行計算預(yù)測值(predict value),得到的結(jié)果與第二組數(shù)據(jù)在 Moldflow中的分析得到的值(true value)進(jìn)行對比,見圖10—12,其中橫坐標(biāo)為試驗組數(shù)。

圖10 普通最小二乘法線性回歸Fig.10 Ordinary least squares linear regression

圖11 普通最小二乘法多項式回歸Fig.11 Ordinary least squares polynomial linear regression

圖12 α值為2時的多項式嶺回歸Fig.12 Polynomial ridge regression when α=2

從圖10—12可知,使用懲罰系數(shù)為2的多項式嶺回歸效果最好,擬合效果指標(biāo)達(dá)到 0.975。說明這個回歸模型是可靠的。且回歸結(jié)果模型見式(4)。

式中:x1,x2,x3,x4,x5分別為注射時間、V/P切換、保壓壓力、熔體溫度、“注射+保壓+冷卻”時間。式(4)各項系數(shù)見表2。

表2 回歸多項式系數(shù)表Tab.2 Coefficients of polynomial regression

在回歸模型基礎(chǔ)上,使用改進(jìn)的模擬退火算法進(jìn)行最優(yōu)參數(shù)組合的搜尋,搜尋過程當(dāng)前值隨著每次降溫的變化曲線見圖13,得到的最優(yōu)解組合見表3。

為驗證本算法的準(zhǔn)確性,在Isight軟件上也使用這一回歸模型(4)進(jìn)行了計算[13],具體系數(shù)見表 2,參數(shù)的取值范圍為表 1中的第一組中實驗參數(shù)的范圍,求出的最優(yōu)組合見表3。

圖13 改進(jìn)的SA算法當(dāng)前解變化曲線Fig.13 Current solution curve of improved SA algorithm

表3 模擬退火算法求出的最優(yōu)參數(shù)組合Tab.3 Optimal combination for parameters of simulated annealing algorithm

由表4可知,兩種模擬退火算法得出最優(yōu)解的參數(shù)組合幾乎一致,除V/P切換和熔體溫度有0.01%的誤差外,其余幾個參數(shù)都一致,說明了本系統(tǒng)改進(jìn)的模擬退火算法準(zhǔn)確有效。

將得到的最優(yōu)參數(shù)組合導(dǎo)入 Moldflow進(jìn)行分析,結(jié)果見圖14,得到的最大變形量為0.6805 mm,比正交試驗25組CAE分析的最優(yōu)結(jié)果0.7436更小,優(yōu)化力度為 9.3%,進(jìn)一步說明本算法準(zhǔn)確有效,且說明回歸與最優(yōu)解搜尋整個程序精確有效。

圖14 最優(yōu)參數(shù)組合的分析結(jié)果Fig.14 Analysis result of the best parameter combination

5 結(jié)論

本系統(tǒng)開發(fā)了3個功能模塊:用戶管理、項目流程管理和集成研發(fā)工具,其能有效地實現(xiàn)對企業(yè)用戶、項目流程的管理以及注塑成型工藝參數(shù)的優(yōu)化。用戶管理和項目流程管理有利于企業(yè)明確分工、加強協(xié)作、降低研發(fā)成本、提高工作效率。集成研發(fā)工具模塊一方面為大量重復(fù)的CAE分析工作提供了一個簡單、高效的解決方案,另一方面基于多元回歸和改進(jìn)的模擬退火算法實現(xiàn)了對成型工藝參數(shù)的優(yōu)化。通過與專業(yè)優(yōu)化軟件的結(jié)果對比,驗證了本系統(tǒng)計算結(jié)果的精確性和有效性。

[1] 王小明, 趙明娟, 陳炳輝, 等. Moldflow 在注塑成型翹曲優(yōu)化中的應(yīng)用[J]. 塑料工業(yè), 2011, 39(4): 49—51.WANG Xiao-ming, ZHAO Ming-juan, CHEN Bing-hui,et al. Application of Moldflow in Warp Optimization of the Injection Molding[J]. China Plastics Industry, 2011,39(4): 49—51.

[2] MATIN I, HADZISTEVIC M, HODOLIC J, et al. A CAD/CAE-integrated Injection Mold Design System for Plastic Products[J]. The International Journal of Advanced Manufacturing Technology, 2012, 63(5): 595—607.

[3] 陳華書, 胡廣洪. 金屬體積成形有限元分析管理與工藝優(yōu)化平臺[J]. 熱加工工藝, 2013, 42(23): 125—128.CHEN Hua-shu, HU Guang-hong. Finite Element Analysis Management and Process Optimization Platform for Metal Volume Forming[J]. Hot Working Technology,2013, 42(23): 125—128.

[4] PEDREGOSA F, VAROQUAUX G, GRAMFORT A, et al. Scikit-learn: Machine Learning in Python[J]. Journal of Machine Learning Research, 2011, 12(10): 2825—2830.

[5] LIU X Q, GAO F, XU J W. Linearized Restricted Ridge Regression Estimator in Linear Regression[J]. Communications in Statistics-Theory and Methods, 2012, 41(24):4503—4514.

[6] 任維雅, 李國輝. 面向監(jiān)督學(xué)習(xí)的稀疏平滑嶺回歸方法[J]. 國防科技大學(xué)學(xué)報, 2015, 37(6): 121—128.REN Wei-ya, LI Guo-hui. Sparse Smooth Ridge Regerssion Method for Supervised Learning[J]. Journal of National University of Defense Technology, 2015, 37(6):121—128.

[7] KONG D, BONDELL H D, WU Y. Domain Selection for the Varying Coefficient Model Via Local Polynomial Regression[J]. Computational Statistics & Data Analysis,2015, 83: 236—250.

[8] 盧靜波, 吳藝能. 非線性回歸模型的線性變換和正交多項式回歸[J]. 統(tǒng)計與決策, 2009, 2009(23): 13—14.LU Jing-bo, WU Yi-neng. Linear Transformation and Orthogonal Polynomial Regression in Nonlinear Regression Model[J]. Statistics and Decision, 2009, 2009(23):13—14.

[9] 黎明, 陳穎, 楊楠. 應(yīng)用回歸分析[M]. 上海: 復(fù)旦大學(xué)出版社, 2008.LI Ming, CHEN Ying, YANG Nan. Application of Regression Analysis[M]. Shanghai: Fudan University Press,2008.

[10] 袁澎, 艾芊, 趙媛媛. 基于改進(jìn)的遺傳-模擬退火算法和誤差度分析原理的PMU多目標(biāo)優(yōu)化配置[J]. 中國電機工程學(xué)報, 2014, 34(13): 2178—2187.YUAN Peng, AI Qian, ZHAO Yuan-yuan. Research on Multi-objective Optimal PMU Placement Based on Error Analysis Theory and Improved GASA[J]. Proceedings of the CSEE, 2014, 34(13): 2178—2187.

[11] 孫寶壽, 陳哲, 吳真繁, 等. 薄壁注塑件翹曲影響因素分析及優(yōu)化研究進(jìn)展[J]. 機械制造, 2009, 47(12): 25—29.SUN Bao-shou, CHEN Zhe, WU Zhen-fan, et al. Analysis and Optimization of Warping Factors of Thin Wall Injection Molding Parts[J]. Machinery, 2009, 47(12): 25—29.

[12] 郎志正. 質(zhì)量管理及其技術(shù)和方法[M]. 北京: 中國標(biāo)準(zhǔn)出版社, 2003.LANG Zhi-zheng. Quality Management and Its Technical Methods[M]. Beijing: Standards Press of China, 2003.

[13] VAN D V A, KOCH P. Isight Design Optimization Methodologies[J]. ASM Handbook, 2010, 22: 79.

Finite Element Analysis Management and Process Optimization System of Injection Molding

LING Cheng-zhi,HU Guang-hong
(Institute of Plastic Forming Technology and Equipment, Shanghai Jiao Tong University, Shanghai 200030, China)

The paper aims to develop an integrated system of project work flow management and process parameter optimization for injection molding finite element analysis. The project process management module was developed in line with the company's project process. Multiple regression methods were used in combination with the improved simulated annealing algorithm to optimize process parameters. Cases of plastic tray warpage deformation optimization were used to verify the accuracy of the optimization method and carry out CAE analysis and orthogonal test for plastic tray. The multivariate regression method was used to fit the test data and calculate the score. The optimal solution was solved with Isight's ASA and the improved simulated annealing algorithm respectively. And their results were compared. CAE analysis was performed with parameters of the optimal solution; the max warpage of analysis result was obtained, and the warpage was compared with those of orthogonal test.Polynomial ridge regression with penalty coefficient of 2 has the best regression effect and itsR2was 0.975. The calculation result of the improved simulated annealing algorithm was consistent with that of Isight's ASA. The warpage of the optimal parameters in the Moldflow analysis was 0.6805 mm, which was smaller than the previous minimum warpage in orthogonal test of 0.7436 mm. The process parameter optimization program based on regression fitting and the improved simulated annealing algorithm can optimize parameters to some extent.

injection molding; finite element analysis; project management; process optimization

2017-12-10

凌成智(1994—),男,碩士研究生,主要研究方向為計算機輔助工程技術(shù)。

胡廣洪(1973—),男,博士,副研究員,主要研究方向為高分子材料成型工藝、高分子材料合金化技術(shù)、材料成型CAD/CAE/CAM。

10.3969/j.issn.1674-6457.2018.01.021

TQ320

A

1674-6457(2018)01-0161-06

猜你喜歡
有限元工藝優(yōu)化
超限高層建筑結(jié)構(gòu)設(shè)計與優(yōu)化思考
民用建筑防煙排煙設(shè)計優(yōu)化探討
關(guān)于優(yōu)化消防安全告知承諾的一些思考
一道優(yōu)化題的幾何解法
轉(zhuǎn)爐高效復(fù)合吹煉工藝的開發(fā)與應(yīng)用
山東冶金(2019年6期)2020-01-06 07:45:54
5-氯-1-茚酮合成工藝改進(jìn)
一段鋅氧壓浸出與焙燒浸出工藝的比較
磨削淬硬殘余應(yīng)力的有限元分析
絡(luò)合鐵脫硫工藝在CK1井的應(yīng)用
基于SolidWorks的吸嘴支撐臂有限元分析
主站蜘蛛池模板: 国产精品夜夜嗨视频免费视频| 六月婷婷精品视频在线观看 | 国产在线视频二区| 天天摸夜夜操| 久久久久久久久久国产精品| a毛片免费在线观看| 欧美区在线播放| 国产高清在线丝袜精品一区| 中国毛片网| 欧美精品高清| 亚洲天堂网在线视频| 久久伊人操| 国产精品自在拍首页视频8| 国产区在线看| 亚洲天堂网在线观看视频| 国产精品成人免费综合| 色网在线视频| 国产一在线观看| 亚洲三级影院| 亚洲国产成人精品无码区性色| 任我操在线视频| 亚洲国产日韩在线观看| 99在线小视频| 免费中文字幕一级毛片| 亚洲中文字幕手机在线第一页| 亚洲午夜福利在线| av在线无码浏览| 国产成人免费高清AⅤ| 少妇露出福利视频| 久久99热66这里只有精品一| 亚洲一区二区约美女探花| 亚洲三级网站| 青青热久免费精品视频6| 日韩毛片免费视频| 国产美女精品一区二区| 丁香婷婷久久| 97超碰精品成人国产| 她的性爱视频| 国产成人高清精品免费软件| 精品无码国产自产野外拍在线| 99热这里只有精品在线播放| 日韩天堂网| 成人毛片免费在线观看| 国产精品深爱在线| 人妻免费无码不卡视频| 亚洲中文久久精品无玛| 欧美激情首页| 久久伊人操| 国内视频精品| 韩日无码在线不卡| 自慰网址在线观看| 亚洲伦理一区二区| 精品少妇三级亚洲| 久久久久亚洲av成人网人人软件| 人人艹人人爽| 精品国产成人国产在线| 99在线小视频| 97在线公开视频| 久久99久久无码毛片一区二区| 丁香婷婷激情综合激情| 伊人久久精品亚洲午夜| 亚洲色图欧美在线| 91年精品国产福利线观看久久| 看国产一级毛片| 国产 在线视频无码| 久久国产亚洲偷自| 国产精品专区第1页| 日本道综合一本久久久88| 亚洲精品无码日韩国产不卡| 美女被躁出白浆视频播放| 亚洲网综合| 99手机在线视频| 久久综合色播五月男人的天堂| 亚洲激情99| 制服丝袜在线视频香蕉| 成年人视频一区二区| 国产日本欧美亚洲精品视| 欧美自拍另类欧美综合图区| 91小视频在线观看免费版高清| 日本91视频| 国产一区二区三区在线观看免费| 露脸一二三区国语对白|