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

基于幾何非線性有限單元法的溫控器用片彈簧優化設計

2015-10-28 10:59:30陳文華李小輝賀青川陳曉英
中國機械工程 2015年5期
關鍵詞:變形優化設計

陳文華 李小輝 潘 駿 賀青川 陳曉英,2 王 孟

1.浙江理工大學浙江省機電產品可靠性技術研究重點實驗室,杭州,3100182.衢州學院,衢州,324000

基于幾何非線性有限單元法的溫控器用片彈簧優化設計

陳文華1李小輝1潘駿1賀青川1陳曉英1,2王孟1

1.浙江理工大學浙江省機電產品可靠性技術研究重點實驗室,杭州,3100182.衢州學院,衢州,324000

利用幾何非線性有限單元法對溫控器用片彈簧參數進行了數值計算,并通過與實際測試結果的對比,驗證了基于幾何非線性有限單元法求解溫控器用片彈簧參數的可行性和正確性。在此基礎上,運用MATLAB優化工具,對片彈簧的結構尺寸進行了優化設計。然后建立觸頭彈跳的ADAMS仿真模型,并利用激光位移測試方法對模型進行了驗證。最后利用ADAMS軟件對溫控器進行觸頭彈跳仿真分析。結果表明:基于幾何非線性有限單元法優化設計的片彈簧能有效抑制動觸頭的彈跳。

幾何非線性有限單元法;溫控器;片彈簧;觸頭彈跳;優化設計

0 引言

溫控器是實現溫度控制的電開關元件,被廣泛應用于溫度控制電路中,因其觸頭材料有一定剛性,在通電過程中,難以避免因動觸頭與靜觸頭間的碰撞而引起動觸頭的彈跳。動觸頭的彈跳是引起電弧,造成觸頭燒蝕的主要原因,會影響溫控器的工作可靠性,并降低壽命[1-2]。文獻[3-7]的研究結果表明,片彈簧的預壓力和剛度等參數對動觸頭的彈跳有較大影響。在片彈簧的參數計算方面,羅依金[8]建立了片彈簧的彎曲微分方程,采用二次積分求解的解析方法對可簡化為懸臂梁模型的直線形片彈簧進行了計算。該方法雖能精確得到片彈簧在力作用點處的位移,但僅適用于對形狀簡單、規則的片彈簧進行計算[9]。Fedder[10]和Peroulis等[11]通過對片彈簧的受力進行分析,計算出了簧片截面承受的彎曲力矩,建立了片彈簧沿長度方向的彈性位能計算模型。盧錦鳳等[12]、翟國富等[13]基于變形能法,提出了離散復雜形體片彈簧、分段求解片彈簧柔度再線性疊加的思路。

上述求解方法均建立于小變形假設的基礎之上,即在片彈簧的整個變形過程中,假設片彈簧的彈性變形量相對初始結構尺寸比較小,忽略因變形而引起的結構尺寸變化。但片彈簧屬于薄壁元件,在一定載荷作用下會產生較大的位移,出現幾何非線性問題,若繼續基于小變形理論對片彈簧進行研究,將會產生比較大的誤差[10,14]。目前研究片彈簧的幾何非線性問題時通常借助于ABAQUS、ANSYS等有限元軟件,但在將其應用在片彈簧的優化設計時,這些軟件則存在分析計算不便、計算過程文件難以查看、求解耗時較長等缺點。

本文以溫控器用片彈簧為研究對象,提出基于幾何非線性有限單元法對其參數進行求解的數值計算方法,并與實測值進行對比驗證。然后建立以片彈簧的預壓力為目標函數,以片彈簧的剛度和最大彎曲應力為約束條件的優化模型,利用MATLAB軟件中的fmincon優化函數進行優化設計。最后利用多體動力學分析軟件ADAMS對溫控器的觸簧系統進行動態仿真分析,所得結果可為片彈簧的設計研究提供技術支撐。

1 幾何非線性有限單元法

1.1物理模型

本文以型號為KSD 301的常閉式溫控器為研究對象,其動作溫度為90 ℃,復位溫度為75 ℃,結構原理如圖1所示。利用鉚釘將片彈簧與底座固定,動觸頭鉚接在片彈簧上。在電路處于接通狀態時,片彈簧彎曲變形產生預壓力P維持動靜觸頭緊密接觸。當外部溫度升高到90 ℃左右時,溫控器的碟形雙金屬片突跳,通過陶瓷桿對片彈簧產生推力,使動靜觸頭分離,從而切斷電路,推力所產生的機械能儲存在片彈簧中。當溫度降低到75 ℃附近時,碟形雙金屬片反向突跳,通過陶瓷桿作用在片彈簧上的推力消失,動觸頭閉合,觸頭之間的碰撞會使動觸頭發生彈跳。

1.碟形雙金屬片 2.靜觸頭 3.動觸頭 4.陶瓷桿 5.片彈簧 6.鉚釘圖1 溫控器的結構原理圖

KSD301型溫控器所用片彈簧包括三條直線段s1、s2、s3,其結構簡圖見圖2,結構尺寸如表1所示。根據片彈簧的實際受力和約束情況,可認為鉚釘鉚接孔附近的平面未發生變形,所以可將其簡化為平面梁系統。

圖2 片彈簧的結構簡圖

l1(mm)l2(mm)l3(mm)寬度b(mm)厚度h(mm)折彎角β(°)折彎半徑r(mm)10.01.65.03.70.18130.2

采用有限元方法對片彈簧建立如圖3所示的有限元模型,圖中數字為有限元節點編號。鉚釘鉚接孔中心在xy平面的投影位置為節點1,陶瓷桿與片彈簧的接觸點為節點13,動觸頭鉚接孔中心位置為節點15。節點1固定,節點15沿y軸負方向作用有預壓力P,該點y方向的位移為δ,則節點15處的剛度:

(1)

圖3 片彈簧的有限元模型

碟形雙金屬片突跳,通過陶瓷桿對片彈簧產生的推力為P1,作用于節點13。在P1作用下,片彈簧的彎曲變形最大,片彈簧的s2段產生最大彎曲應力:

(2)

式中,Mmax為片彈簧所受的最大彎矩。

1.2平面梁幾何非線性有限單元法

對于節點i和j、長為l、截面積為A的梁單元,在局部坐標系下,其節點位移列陣可表示如下:

(3)

其中,u、v、θ對應著軸向位移、撓度和轉角。

選擇多項式作為位移的插值函數,插值函數的矩陣表達形式為[15]

(4)

式中,δe為單元節點位移向量。

略去梁單元的剪切變形,則其應變包括拉壓應變εa和彎曲應變εb[15],表示為

(5)

將式(5)代入式(4),得

(6)

在大位移情況下,梁單元的應變和位移的關系式是非線性的[15],可表達為

ε=(B0+BLδe)×δe

(7)

式中,B0為一般線性分析時的應變矩陣,與單元節點位移無關;BLδe為單元節點位移向量δe的函數。

由式(6)、式(7)可知:

(8)

(9)

根據虛位移原理:梁單元內應力由于虛應變dε所做的功等于外力在虛位移dδe上所做的功,若用φ表示內力與外力的矢量總和,在局部坐標系下,梁單元的平衡方程式為[15]

d(δe)Tφ=∫dεTσdV-d(δe)TR=0

(10)

式中,σ為單元應力;R為載荷列陣。

根據一般的線彈性關系,在不考慮初應力和初應變的情況下,可將式(10)簡化為

(11)

式(11)中右邊第一項可改寫為

(12)

其中,k0為線性剛度矩陣,與單元的節點位移無關;kL為大位移矩陣,反映了單元位移的變化對平衡方程的影響情況[15]。

式(11)中右邊第二項可表示為

(13)

式中,kσ為初應力矩陣。

kT=k0+kL+kσ

(14)

式中,kT為切線剛度矩陣。

根據式(12)~式(14),將式(11)簡化為

dφ=kTd δe

(15)

式(15)是梁單元在局部坐標系下的平衡方程,對于梁系統,應將式(15)中φ、kT和δe轉換到整體坐標系下,再組裝形成系統的失衡力向量ψ、切線剛度矩陣kT和節點位移向量δ,則對系統建立平衡方程式為

dψ=kTdδ

(16)

1.3方程的求解

牛頓-拉斐遜迭代法是求解非線性方程的常用方法[15],其迭代格式如下:

kTnΔδn+1=R-F δn

(17)

δn+1=Δδn+1+δn

式中,F為系統的節點力向量。

具體迭代時,將以線性有限元法求解得到的節點位移向量δ0作為迭代求解的初始值,當節點15處y方向的位移變化率小于給定值時,迭代過程即可結束。

1.4片彈簧預壓力的求解

以上求解過程是在已知外載荷的情況下實現對節點位移的求解的。根據溫控器用片彈簧的結構特點和實際約束情況,當動靜觸頭處于閉合狀態時,片彈簧節點15處的y方向位移δ可通過計算近似得到,即δ≈l1sinβ,此時在節點15處產生的彈性變形力即為片彈簧的預壓力P。

選取較小的預壓力P0作為初始值,利用幾何非線性有限單元法得到節點15處的y方向位移δ0,按照以下迭代公式對預壓力P進行求解:

(18)

其中,kn為第n次求解的剛度。當|δn+1-δ|≤ξ時,迭代停止,ξ為給定的誤差,此時得到的預壓力Pn+1即是預壓力P。

1.5最大彎曲應力的求解

片彈簧在推力P1作用下產生最大彎曲應力。對于KSD301型溫控器,動靜觸頭分開0.24 mm即可有效地切斷電路。在P1作用下,片彈簧彎曲變形到穩定狀態時,片彈簧節點15處的y方向位移δ可通過計算近似得到,即δ≈l1sinβ+0.24。同理按照式(18)的迭代格式可求解出推力P1的大小,進而按照幾何非線性有限單元法確定出片彈簧所受的最大彎矩Mmax,由式(2)可求出最大彎曲應力σmax。

1.6理論計算與試驗結果的對比

為了驗證利用幾何非線性有限單元法對溫控器用片彈簧求解計算的可行性和準確性,現利用LDW-1型微機控制電子萬能材料試驗機對片彈簧進行壓縮試驗如圖4所示。試驗選取3個片彈簧試樣,采用3次試驗的平均值作為最終的試驗結果。

圖4 試驗測試

利用變形能法[12]、線性有限單元法、幾何非線性有限元方法,按表1所列的片彈簧尺寸參數,對節點15處的y方向位移δ進行求解計算。通過圖5所示的結果可以看出,柔度法和線性有限單元法的求解結果比較接近,利用幾何非線性有限單元法得到的求解結果最小。通過試驗值的變化曲線可以看出:隨著外載荷的不斷增大,片彈簧的節點位移呈現出非線性增大的特點。利用變形能法和線性有限單元法得到的節點位移隨著外載荷的增大保持線性增大的趨勢,在外載荷大于0.9N時,求解結果與試驗值的誤差開始變得明顯,并隨著力的增大誤差逐漸變大。利用幾何非線性有限單元法得到的結果表明,變形位移與外載荷保持非線性增大的特點,與試驗結果能夠較好地吻合。

圖5 求解結果與試驗結果的對比

2 片彈簧的優化設計

2.1設計變量

選取表1所示的結構參數作為設計變量X:

X=(x1,x2,x3,x4,x5,x6,x7)=(l1,l2,l3,b,h,β,r)

設置變量的初始值和上下限如表2所示。

表2 設計變量的初始值和上下限

2.2目標函數

片彈簧的剛度和預壓力等參數對動觸頭的彈跳有較大影響,在合適的范圍內增大片彈簧的剛度和預壓力能有效抑制動觸頭的彈跳,其中提高預壓力值更為明顯[7]。本文選取片彈簧的預壓力作為目標函數。對可分斷的電接觸系統,預壓力P與系統的電流i有關[1],確定出預壓力的范圍為0.6~1.5N。最終確定出目標函數為

minf(P)=|P-1.5|

2.3約束條件

(1)根據溫控器廠家的經驗數據,片彈簧合適的剛度k范圍為0.5~0.8 N/mm。為了縮短觸頭彈跳時間,應使片彈簧具有較大的剛度[7]。現選擇片彈簧剛度k范圍為0.7~0.8 N/mm。

(2)該型溫控器所用片彈簧材料為鈹青銅(QBe2)。片彈簧的最大彎曲應力應小于等于片彈簧材料的許用應力,即σmax≤[σ]=600 MPa。

(3)對于一般片彈簧,其寬度b和厚度h應滿足2≤b/h≤10;為了避免片彈簧在彎曲處存在過大的應力集中,規定h/(2r)≤0.5;為了保證在預壓力作用下,片彈簧節點15處y方向變形位移遠大于s2段和s3段的y方向位移,規定l1在xz平面上的投影長度大于l3,即:l3≤l1cosβ。

2.4溫控器用片彈簧的優化設計結果

現采用MATLAB軟件中的可求解單目標約束優化問題的fmincon函數對溫控器用片彈簧進行優化設計。應用幾何非線性有限單元法對片彈簧的預壓力、剛度等參數進行計算,優化后的結構參數如表3所示。

表3 優化設計前后結果對比

3 片彈簧優化設計結果的仿真驗證

3.1仿真分析與試驗驗證

本文利用多體動力學分析軟件ADAMS對突跳式溫控器的觸頭彈跳運動進行仿真分析。首先利用有限元分析軟件ANSYS生成片彈簧的模態中性文件,然后將其導入到ADAMS分析軟件中。在動靜觸頭之間添加接觸,并采用碰撞函數法對其進行求解[2],所得動觸頭運動參數求解結果如圖6所示。

圖6 動觸頭彈跳的仿真結果

由圖6可知,動觸頭在閉合過程中的最大速度為291.6 mm/s,最大接觸力為10.8 N,彈跳時間為4.3 ms。作用于片彈簧的外力消失后,儲存在片彈簧中的變形能釋放,動觸頭的動能不斷增大,速度也不斷增大。動靜觸頭首次接觸時,觸頭間的劇烈碰撞使動觸頭的速度急劇減小,觸頭間的接觸力從零變為最大,動觸頭反向彈起。隨著動觸頭動能的不斷減小,動觸頭的彈起高度逐漸降低,直到動能被完全消耗,動觸頭才達到靜止、穩定的閉合狀態。

為驗證仿真分析的準確性,對動觸頭的彈跳位移進行試驗測試,所用儀器設備如圖7所示。動觸頭閉合時的彈跳位移由LK-G80型激光位移傳感器進行采集,經LK-GD500型控制器進行設置和處理后,輸送到計算機中進行數據的記錄和存儲。

圖7 試驗測試示意圖

圖8 觸頭彈跳位移對比

3.2優化結果的仿真分析

為檢驗經過優化設計后的片彈簧能否改善動觸頭的彈跳,將優化后的片彈簧導入到仿真模型中,并使模型中其他部件和接觸力的參數與優化前的模型保持一致,得到的仿真結果如圖9所示。

觸頭在彈跳過程中,彈跳時間越長,彈跳引起的電弧能量越大,對觸頭的燒蝕越嚴重[7],由此引起的接觸失效越容易發生,因此彈跳時間是評價觸頭彈跳的重要參數[2-7]。如圖9所示,優化后的彈跳時間為2.3 ms,相對于優化前的片彈簧,彈跳時間縮短46.5%;對比圖6和圖9所示結果可知,優化前彈跳位移波動較大,位移變化劇烈,優化后的彈跳位移變化平緩,最大彈跳高度很小,說明動觸頭的彈跳得到抑制。

圖9 優化后的仿真結果

由圖9可知,對于優化后的觸簧系統,觸頭閉合過程中的最大速度為499.5 mm/s,最大接觸力為22.05 N,均大于優化前的觸簧系統。這主要是因為優化后的片彈簧的預壓力和剛度都變大,在動觸頭分開距離一定的情況下,優化后的片彈簧中儲存的變形能增大,動觸頭的動能更大,導致觸頭間的碰撞更加劇烈。大的預壓力和剛度會使觸頭的回彈變得更加困難,雖然觸頭間的瞬時沖擊力變大,但該沖擊力并不能大幅度地推動動觸頭,所以觸頭閉合瞬時的位移變化比較平緩,彈跳時間較短。

4 結論

(1)本文主要考慮片彈簧在外載荷作用下變形位移較大,引入幾何非線性有限單元法求解,能反映出變形位移對結構平衡的影響,以此為基礎對溫控器用片彈簧進行有限元建模和求解。結果和試驗值較為吻合,表明用幾何非線性有限單元法對片彈簧進行分析更符合片彈簧的實際工作情況。

(2)由于動觸頭的碰撞彈跳與片彈簧的預壓力有較大關系,為了改善動觸頭的碰撞彈跳性能,以片彈簧的預壓力為目標函數,其尺寸參數為設計變量,運用幾何非線性有限單元法,實現了對片彈簧結構的優化設計。

(3)為驗證優化設計后的片彈簧是否能抑制動觸頭的彈跳,利用多體動力學軟件ADAMS對動靜觸頭閉合過程中的彈跳進行了仿真分析。結果表明,基于幾何非線性有限單元法設計的片彈簧能有效抑制動觸頭的彈跳。

[1][1]程禮椿.電接觸理論及應用[M].北京:機械工業出版社,2004.

[2]熊軍.航天繼電器觸簧系統接觸彈跳及其影響因素研究[D] .武漢:華中科技大學,2008.

[3]Mcbride J W,Sharkh S M.Electrical Contact Phenomena during Impact[J].IEEE Transactions on Components, Hybrids and Manufacturing Technology,1992, 15(2):184-192.

[4]Li Y, Meguid S A, Fu Y, et al. Unified Nonlinear Quasistatic and Dynamic Analysis of RF-MEMS Switches[J]. Acta Mechanica, 2013,224(8): 1-15.

[5]LaRose III R P, Murphy K D. Impact dynamics of MEMS Switches[J]. Nonlinear Dynamics, 2010, 60(3): 327-339.

[6]McCarthy B, Adams G G, McGruer N E, et al. A Dynamic Model, Including Contact Bounce of an Electrostatically Actuated Micro-switch[J]. Microelectromechanical Systems, Journal of, 2002, 11(3): 276-283.

[7]紐春萍,陳德桂,李興文 等.交流接觸器觸頭彈跳的仿真及影響因素[J].電工技術學報,2007,22(10):85-90.

Niu Chunping,Chen Degui,Li Xingwen. et al. Simulation of Contact Bounce of AC Contactor and Study of Its Influence Factors[J].Transactions of China Electrotechnical Society,2007,22(10):85-90.

[8]羅依金.小型密封電磁繼電器[M].王蓉芳,譯.北京:人民郵電出版社,1979.

[9]葉雪榮,梁慧敏,翟國富,等.電磁繼電器簧片反力特性的通用計算方法[J].低壓電器,2007(3):10-12.

Ye Xuerong,Liang Huimin,Zhai Guofu,et al.Universal Calculation Method for Electromagnetic Relay Spring Force Characteristic[J].Low Voltage Apparatus,2007(3):10-12.

[10]Fedder G K,Simulation of Microelectromechanical Systems[D].California:University of California,1994.

[11]Peroulis D,Pacheco S P,Sarabandi K.Electromechanical Considerations in Developing Low-voltage RF MEMS Switches[J].IEEE Transactions on Microwave Theory and Techniques,2003,51(1):259-270.

[12]盧錦鳳,翟國富,劉茂凱.變形能法在計算繼電器簧片柔度中的應用[J].機電元件,1995(2/3):14-22.

Lu Jinfeng,Zhai Guofu,Liu Maokai.Using the Strain Energy Method to Calculate the Relay’s Flat Spring[J].Electromechanical Component,1995(2/3):14-22.

[13]翟國富,樊薇薇,梁慧敏.基于均勻試驗設計的航天電磁繼電器觸點滑移長度優化[J].電工技術學報,2009,24(10):59-64.

Zhai Guofu,Fan Weiwei,Liang Huimin.An Optimization Method for the Contact Slip Length of Space Electromagnetic Relay Based on Uniform Experimental Design[J].Transactions of China Electrotechnical Society,2009, 24(10):59-64.

[14]Fatola B O,Keogh P,Hicks B.Modelling Flat Spring Performance Using FEA[J].Journal of Physics: Conference Series.IOP Publishing,2009,181(1):1-9.

[15]王勖成. 有限單元法基[M]. 北京:清華大學出版社, 2003.

(編輯王艷麗)

Optimal Design of Flat Spring Used in Thermostat Based on Geometrically Nonlinear Finite Element Method

Chen Wenhua1Li Xiaohui1Pan Jun1He Qingchuan1Chen Xiaoying1,2Wang Meng1

1.Zhejiang Province’s Key Laboratory of Reliability Technology for Mechanical &Electrical Product,Zhejiang Sci-Tech University,Hangzhou,310018 2.Quzhou University,Quzhou,Zhejiang,324000

A geometrically nonlinear finite element method was applied herein to the numerical calculation of the flat spring used in a thermostat. The calculation results were verified by comparison with experiments. Then the optimal design of the flat spring’s construction parameters were completed by using the MATLAB optimization toolbox. At last, the contact bounce of the thermostat was simulated by ADAMS. The results show that the flat spring which is designed optimally based on the geometrical nonlinear finite element method can prevent the contact bounce effectively.

geometrically nonlinear finite element method;thermostat;flat spring;contact bounce;optimal design

2014-02-27

長江學者和創新團隊發展計劃資助項目(IRT13097);浙江省重點科技創新團隊項目(2010R50005)

TM564.8< class="emphasis_italic">DOI

:10.3969/j.issn.1004-132X.2015.05.003

陳文華,男,1963年生。浙江理工大學機械與自動控制學院教授、博士研究生導師。主要研究方向為可靠性工程、機械傳動與機構。李小輝,男,1990年生。浙江理工大學機械與自動控制學院碩士研究生。潘駿,男,1974年生。浙江理工大學機械與自動控制學院教授。賀青川,男,1984年生。浙江理工大學機械與自動控制學院講師。陳曉英,女,1974年生。衢州學院機械工程學院副教授,浙江理工大學機械與自動控制學院博士。王孟,男,1989年生。浙江理工大學機械與自動控制學院碩士研究生。

猜你喜歡
變形優化設計
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
談詩的變形
中華詩詞(2020年1期)2020-09-21 09:24:52
瞞天過海——仿生設計萌到家
藝術啟蒙(2018年7期)2018-08-23 09:14:18
“我”的變形計
設計秀
海峽姐妹(2017年7期)2017-07-31 19:08:17
有種設計叫而專
Coco薇(2017年5期)2017-06-05 08:53:16
例談拼圖與整式變形
主站蜘蛛池模板: 国产欧美视频综合二区| 日韩AV无码免费一二三区| 亚瑟天堂久久一区二区影院| 日本a级免费| 一级毛片在线免费看| 国产成人精品免费av| 成人在线视频一区| 日韩高清一区 | 亚洲AV色香蕉一区二区| 国产视频资源在线观看| 伊人久热这里只有精品视频99| 欧美精品亚洲精品日韩专区va| 国产不卡一级毛片视频| 四虎在线观看视频高清无码| 国产美女91视频| 91福利免费视频| 欧美综合在线观看| 亚洲精品国产综合99| 综1合AV在线播放| 一级毛片基地| 国产日韩欧美黄色片免费观看| 成人免费网站久久久| 亚洲天堂在线视频| 中文字幕丝袜一区二区| 一边摸一边做爽的视频17国产| 国产在线日本| 日本三级黄在线观看| 成人午夜天| 亚洲国产日韩在线观看| 欧美日韩一区二区在线播放| 亚洲国产精品国自产拍A| 亚洲国产精品不卡在线| 亚洲毛片一级带毛片基地| 强乱中文字幕在线播放不卡| 亚洲 欧美 偷自乱 图片 | 日本人真淫视频一区二区三区| 中文字幕人妻av一区二区| 啪啪永久免费av| 国产成人一区二区| 丰满少妇αⅴ无码区| 国产综合精品一区二区| 大陆国产精品视频| 国产成人成人一区二区| 日日噜噜夜夜狠狠视频| 2018日日摸夜夜添狠狠躁| 九九热视频精品在线| 久久免费成人| 欧美激情第一欧美在线| 女人一级毛片| 色天堂无毒不卡| 久久亚洲精少妇毛片午夜无码| 欧美中文字幕第一页线路一| 五月激激激综合网色播免费| 欧美国产菊爆免费观看| 日韩免费毛片视频| 国产91小视频| 青青草久久伊人| 欧美在线三级| 熟女视频91| 亚洲 日韩 激情 无码 中出| 欧美一区精品| a级毛片免费网站| 国产久草视频| 成人a免费α片在线视频网站| 91香蕉国产亚洲一二三区| 色综合国产| 国产福利在线免费| 久青草免费在线视频| 日韩资源站| 亚洲第一成人在线| 99青青青精品视频在线| 国产一级无码不卡视频| 国产好痛疼轻点好爽的视频| 欧美第一页在线| 欧美综合中文字幕久久| 无码高潮喷水专区久久| 国产美女在线观看| 国产91在线免费视频| 久久亚洲国产视频| 成人av专区精品无码国产| 欧美午夜在线观看| 国产精品露脸视频|