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ⅴ无码专区亚洲av综合网| 日本精品视频一区二区| 国产福利在线免费| 2048国产精品原创综合在线| 色首页AV在线| 国产精品久久久久婷婷五月| 老司机精品99在线播放| 亚洲AⅤ综合在线欧美一区| 成年看免费观看视频拍拍| 日韩人妻无码制服丝袜视频| 97青青青国产在线播放| 精品午夜国产福利观看| 久久精品国产亚洲麻豆| 在线精品亚洲国产| 噜噜噜久久| 人妻免费无码不卡视频| 麻豆AV网站免费进入| 4虎影视国产在线观看精品| 亚洲综合第一区| 亚洲欧美在线综合一区二区三区| 亚洲国产精品国自产拍A| 欧美高清视频一区二区三区| 丁香五月婷婷激情基地| 日韩无码视频播放| 精品天海翼一区二区| 激情亚洲天堂| 五月婷婷中文字幕| 国产亚洲欧美在线专区| 亚洲精品视频免费看| 国产欧美日韩另类精彩视频| 伊人91视频| 超清人妻系列无码专区| 日韩在线观看网站| 亚洲欧美成人网| 国产精品密蕾丝视频| 欧美激情视频二区三区| 亚洲欧美成人综合| 亚洲国产日韩在线成人蜜芽| 日本高清有码人妻| 国产尤物jk自慰制服喷水| 国内毛片视频| 亚洲清纯自偷自拍另类专区| 久久成人国产精品免费软件 | 国产白丝av| av色爱 天堂网| 欧美日韩在线第一页| 国产一级视频久久| 久久狠狠色噜噜狠狠狠狠97视色| 精品视频一区二区观看| 免费国产黄线在线观看| 亚洲a免费| 青草国产在线视频| 国产亚洲欧美日韩在线一区二区三区| 国产成人a在线观看视频| 欧美成人午夜视频免看| 日韩午夜片| 少妇人妻无码首页| 亚洲AⅤ波多系列中文字幕| 亚洲女同一区二区| 啊嗯不日本网站| 国产香蕉在线| 三上悠亚一区二区| 亚洲天堂网视频| 国产欧美日韩免费| 丁香六月综合网| a亚洲视频| 五月婷婷综合在线视频| 国产无码制服丝袜| 久久精品一品道久久精品| 在线观看无码a∨| 波多野一区| 在线播放国产一区| 国产亚洲美日韩AV中文字幕无码成人 | 久久午夜夜伦鲁鲁片不卡| 伊人久热这里只有精品视频99| 老司机久久99久久精品播放| 国产99视频在线|