亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        基于傳熱學(xué)的爐溫曲線研究

        2020-05-18 14:22:20方灝航
        現(xiàn)代信息科技 2020年20期

        摘? 要:文章針對設(shè)計(jì)回流焊爐焊接系統(tǒng)各溫區(qū)對應(yīng)下的最優(yōu)溫度和傳送帶過爐速度,使用牛頓冷卻定律與傅里葉定律,構(gòu)建了沿傳送帶方向試件經(jīng)過各溫區(qū)溫度的差分模型;進(jìn)一步設(shè)計(jì)一種搜索算法,在允許的范圍內(nèi),尋求出試件最大傳送帶過爐速度;最終設(shè)立合理的指標(biāo),使用蒙特卡洛法和模擬退火算法全局尋優(yōu),得到了各溫區(qū)對應(yīng)下的最優(yōu)溫度和傳送帶過爐速度。

        關(guān)鍵詞:回焊爐;爐溫曲線;傳熱學(xué);差分方程;搜索算法;模擬退火算法

        中圖分類號:TP273;TN405? ? ? 文獻(xiàn)標(biāo)識碼:A 文章編號:2096-4706(2020)20-0018-05

        Furnace Temperature Curve Research Based on Heat Transfer

        FANG Haohang

        (School of Mechanical Engineering Yangtze University,Jingzhou? 434023,China)

        Abstract:In order to design the optimal temperature and the conveyor belt passing speed in each temperature area of the reflow furnace welding system,Newtons law of cooling and Fouriers law were used to construct the difference model of the temperature of specimens passing through each temperature area along the conveyor belt direction. A search algorithm is further designed to find the maximum conveyor belt passing speed within the allowable range. In the end,a reasonable index was set up,Monte Carlo method and simulated annealing algorithm were used to search for global optimization,and the optimal temperature and conveyor belt passing speed under each temperature region were obtained.

        Keywords:reflow furnace;furnace temperature curve;heat transfer;difference equation;search algorithm;simulated annealing algorithm

        0? 引? 言

        回流焊是目前應(yīng)用于集成電路板最廣的焊接工藝,本文對2020年全國大學(xué)生數(shù)學(xué)建模競賽題進(jìn)一步研究,為提高回流焊接工藝的效益,提出一種設(shè)計(jì)回焊爐內(nèi)各溫區(qū)最佳溫度的方法。通過測量的電路板在爐內(nèi)任一位置溫度(爐溫曲線),調(diào)整模型中各參數(shù)。提出一種搜索算法,尋求出電路板最大過爐速度。由于焊接過程需要滿足一定的要求,建立了合理的評估指標(biāo),最終使用蒙特卡洛法與模擬退火算法進(jìn)行全局尋優(yōu),尋求出回焊爐各個溫區(qū)的最佳溫度以及電路板的最佳過爐速度。作者基于參加2020全國大學(xué)生數(shù)學(xué)建模競賽的經(jīng)驗(yàn),對回焊爐有一定了解,希望通過建模的方法探尋最佳爐溫系統(tǒng),在整個過程中主要負(fù)責(zé)數(shù)學(xué)模型的推導(dǎo)與構(gòu)建以及代碼的實(shí)現(xiàn)。

        1? 回焊爐介紹

        1.1? 回焊爐組成

        回焊爐分成4個大溫區(qū):預(yù)熱區(qū)、恒溫區(qū)、回流區(qū)、冷卻區(qū),如圖1所示。本文中根據(jù)賽題中已知的回焊爐尺寸以及一次回流焊數(shù)據(jù)進(jìn)行模型的構(gòu)建,其中含有11個小溫區(qū)及爐前區(qū)域和爐后區(qū)域,小溫區(qū)長度為30.5 cm,間隙5 cm,爐前、爐后區(qū)域長25 cm。

        1.2? 回焊爐內(nèi)環(huán)境與生產(chǎn)要求

        回焊爐啟動后,爐內(nèi)溫度短時間穩(wěn)定。間隙不做特殊控制,其溫度與相鄰溫區(qū)的溫度有關(guān)。實(shí)際生產(chǎn)時各小溫區(qū)設(shè)定溫度可以進(jìn)行±10 ℃調(diào)整。調(diào)整時要求小溫區(qū)1~5中的溫度保持一致,小溫區(qū)8~9中的溫度保持一致,小溫區(qū)10~11中的溫度保持25 ℃。傳送帶的過爐速度調(diào)節(jié)范圍為65~100 cm/min。在回焊爐電路板焊接生產(chǎn)中,爐溫曲線應(yīng)滿足一定的要求,稱為制程界限,如表1所示。

        1.3? 回焊爐內(nèi)溫度分析

        小溫區(qū)空氣溫度可視為均勻分布。沿著傳送帶方向焊接區(qū)域溫度沿著傳送帶方向不盡相同。因此本文僅研究傳送帶中心環(huán)境的一維傳熱過程,將回焊爐分為如圖1所示的三類區(qū)域:小溫區(qū)對應(yīng)傳送帶上的區(qū)域Di(i=1,2,…,11),間隙對應(yīng)傳送帶上的區(qū)域Bj(j=1,2,…,10),單熱源影響區(qū)域?qū)?yīng)傳送帶上的區(qū)域Ck(k=1,2)。

        2? 模型的建立

        2.1? 模型構(gòu)建思路

        由于電路板溫度與軌道焊接區(qū)域溫度直接相關(guān),而焊接區(qū)域溫度又與兩邊小溫區(qū)溫度相關(guān),則本文建模主要分兩部分:

        第一部分構(gòu)建模型旨在設(shè)計(jì)回焊爐焊接區(qū)域不同位置的溫度分布,考慮在熱傳導(dǎo)與熱對流的兩種熱傳遞方式下,運(yùn)用牛頓冷卻定律以及傅里葉定律,構(gòu)建爐內(nèi)氣體傳熱模型以及間隙環(huán)境傳熱模型。分別構(gòu)建出小溫區(qū)環(huán)境溫度、間隙環(huán)境溫度、爐前爐后環(huán)境溫度模型并整合為環(huán)境溫度分布模型。

        第二部分構(gòu)建模型旨在探究電路板溫度與焊接區(qū)域溫度之間的關(guān)系,由于回焊爐主要加熱方式為小溫區(qū)持續(xù)噴射出的熱風(fēng)與電路板直接的熱傳遞,因此本文主要考慮對流傳熱,運(yùn)用牛頓冷卻定律,構(gòu)建電路板溫度變化差分模型。

        綜上所述繪制出本文模型構(gòu)建的思維導(dǎo)圖,如圖2所示。

        2.2? 環(huán)境溫度分布模型

        2.2.1? 小溫區(qū)Di環(huán)境溫度模型

        小溫區(qū)產(chǎn)生的熱量持續(xù)向空氣輸入最終趨于平衡。則傳送帶中心位置氣體溫度Uat=Ui。Ui為第i個小溫區(qū)溫度。

        2.2.2? 間隙Bj的環(huán)境溫度模型

        由熱力學(xué)定律可知,熱量可以自發(fā)地從高溫物體傳給低溫物體。工程上任何復(fù)雜的傳熱現(xiàn)象都是由傳導(dǎo)導(dǎo)熱、對流傳熱、熱輻射三種基本傳熱方式組成[1]。間隙Bj的環(huán)境溫度變化主要通過熱傳導(dǎo)和熱對流產(chǎn)生。為探求間隙的環(huán)境溫度變化,進(jìn)行如下推導(dǎo)。

        2.2.2.1? 爐內(nèi)氣體傳熱

        傳送帶不同區(qū)域上的溫度變化是以氣體作為介質(zhì),易熱傳導(dǎo)的方式熱量交換。設(shè)兩個小溫區(qū)i、i+1溫度分別為Ui、Ui+1(Ui+1>Ui),間隙的溫度還未平衡時,由熱力學(xué)第二定律,熱量能夠自發(fā)地從高溫物體傳遞給低溫物體。根據(jù)傅里葉定律:單位時間內(nèi)流過單位橫截面積的熱量q與溫度的下降成正比[2]。

        (1)

        其中,負(fù)號為熱量流動方向與溫度升高方向相反,λ為導(dǎo)熱系數(shù), 為溫度沿x方向變化率。

        如圖3所示,設(shè)任意位置x0處t時刻的溫度為U(x0,t),如圖做兩個垂直熱流方向的橫截面積ΔS。dt時間內(nèi)通過x0+dx處ΔS橫截面流入小柱體內(nèi)的熱量為:

        流出小柱體內(nèi)的熱量為:

        由上述兩式可得dt時間內(nèi)小柱體凈流入量為:

        (2)

        由于,故可化簡為:

        (3)

        記dV=ΔSdx,則dt時間內(nèi)小柱體單位體積內(nèi)的凈熱量為:

        (4)

        已知,熱傳導(dǎo)介質(zhì)的密度為ρ,比熱容為c,dt時間內(nèi)引起小柱體單位體積溫度升高所需的熱量為dq=cρ(Ut+ dt-Ut),而dt時間內(nèi)流入小柱體的熱量即為使小柱體溫度升高所需的熱量由于,整理可得一維的熱傳導(dǎo)方程為:

        (5)

        當(dāng)溫區(qū)間的溫度穩(wěn)定后,溫度將不隨時間變化,因此爐內(nèi)空氣溫度是關(guān)于位置的函數(shù)U(x),溫度梯度? 也不隨時間變化。

        2.2.2.2? 間隙氣體對流傳熱模型

        通過氣體對流進(jìn)行熱量交換的方式為熱對流,其能量交換滿足牛頓冷卻定律[2]:

        Q=Sα(Uat-U)? ? ? ? ? ? ? ? ? ? ? ? ? ? ?(6)

        其中,Uat為小溫區(qū)焊接區(qū)域氣體溫度,U為間隙焊接區(qū)域氣體溫度,S為垂直高溫氣流方向的截面積,α為換熱系數(shù),反映對流換熱程度的強(qiáng)弱。

        由式(6)可知單位時間小溫區(qū)焊接區(qū)域氣體吸收的熱量dQ=Sα(Uat-U)dt,與氣體溫度升高dU需要吸收的熱量dQ=ρcVdU相等,即:

        Sα(Uat-U)dt=ρcVdU? ? ? ? ? ? ? ? ? ? ? ? ?(7)

        式中,ρ為小溫區(qū)焊接區(qū)域氣體密度,c為小溫區(qū)焊接區(qū)域氣體比熱容,V為電路板體積。設(shè)時間從t到t+dt時,溫度從U1變化到U2,對式(7)兩邊同時積分得 ,兩邊同時求原函數(shù)并取對數(shù)得Uat-U1=(Uat-U2) 記 (τ>0),則:

        U2=U1+(Uat-U1)(1-)? ? ? ? ? ? ? ? ? (8)

        根據(jù)以上推導(dǎo),最終選擇直線U(x)為間隙焊接區(qū)域氣體溫度分布:

        U(x)=? ? ? ? ? ? ? ? ?(9)

        理由如下:式(9)滿足傳導(dǎo)傳熱規(guī)律,當(dāng)溫度不隨時間變化,溫度是關(guān)于位置的函數(shù)U(x),溫度梯度? 也不隨時間變化;式中U(x)滿足對流傳熱過程中的牛頓冷卻定律;對于回焊爐內(nèi)間隙長度為5 cm較小,用直線擬合產(chǎn)生的誤差也較小;氣體溫度分布滿足流體的連續(xù)介質(zhì)假說[3]。

        2.2.3? 單熱源影響區(qū)域Ck溫度

        單熱源影響區(qū)域是指受單個熱源影響的區(qū)域。由于單熱源影響區(qū)域一邊溫度較高,熱傳導(dǎo)的能量會從另外一邊散失的很快,該區(qū)域氣體溫度升高所需的能量的主要是通過對流換熱,其溫度變化應(yīng)滿足牛頓冷卻定律。U0為生產(chǎn)車間的溫度,Us為單熱源的溫度,Ux為距離單熱源x位置處的氣體溫度,由式(8)可得溫度在單熱源影響區(qū)域的分布為:

        Ux=Us+(U0-Us)(1-e-xψi)? ? ? ? ? ? ? ? ?(10)

        其中,ψi為第i個單熱源對鄰近區(qū)域溫度的影響程度,其大小反映了介質(zhì)導(dǎo)熱能力的強(qiáng)弱,不同單熱源影響的區(qū)域,參數(shù)ψi不同。

        2.2.4? 環(huán)境溫度分布模型

        根據(jù)上述模型的推導(dǎo),聯(lián)立式(9)(12)可得爐內(nèi)環(huán)境溫度在傳送帶上的一維分布:

        (11)

        2.3? 爐溫曲線差分模型

        由于電路板附近熱氣流從小溫區(qū)噴射出,最終到達(dá)電路板表面與其發(fā)生熱傳遞。傳送帶周圍的氣體與電路板之間熱傳遞的方式為熱對流,則電路板表面溫度應(yīng)滿足對流換熱式(8)。

        將式(8)寫成差分形式:

        (12)

        其中,U(t)為t時刻電路板的表面溫度,Uat(vt)為電路板在t時刻的位置所對應(yīng)的環(huán)境溫度,β為電路板與環(huán)境的換熱系數(shù)。

        實(shí)際電路板表面的溫度隨時間的變化有一定的滯后性[4],即:

        U′(t)=U(t-γΔt)=U(t-Δt-γΔt)+[Uat(vt-vγΔt)

        -U(t-γΔt-Δt)](1-)? ? ? ? ? ? ? ? ? ? ? (13)

        將此模型記為模型Ⅰ,式中,γ為電路板表面溫度的延時常數(shù)。

        3? 參數(shù)確定

        本文根據(jù)賽題中第一問數(shù)據(jù):各溫區(qū)設(shè)定的溫度分別為175 ℃(小溫區(qū)1~5)、195 ℃(小溫區(qū)6)、235 ℃(小溫區(qū)7)、255 ℃(小溫區(qū)8~9)及25 ℃(小溫區(qū)10~11),傳送帶的過爐速度為70 cm/min的實(shí)際與理論爐溫曲線的偏離程度作為評判,利用實(shí)驗(yàn)數(shù)據(jù)調(diào)整爐溫曲線差分模型中各參數(shù),最終得到的ψ1=0.160,ψ2=0.014,ψ3=0.019,β=50,γ=2,最終理論與實(shí)際爐溫曲線對比圖如圖4所示。

        4? 尋求最大過爐速度

        實(shí)際生產(chǎn)過程中,電路板的過爐速度越快,則電路板焊接數(shù)量越多,效益越好。然而過爐速度并不能無限制增大,還需滿足制程界限如表1所示。因此為了使效益最大化,本文選取各溫區(qū)溫度的設(shè)定值分別為182 ℃(小溫區(qū)1~5)、203 ℃(小溫區(qū)6)、237 ℃(小溫區(qū)7)、254 ℃(小溫區(qū)8~9)爐溫曲線滿足制程界限的情況下,尋求出符合要求的最大傳送帶過爐速度。

        4.1? 最大過爐速度模型建立

        約束條件可由制程界限推導(dǎo)出。目標(biāo)函數(shù):需尋求出符合要求的最大傳送帶的過爐速度。則目標(biāo)函數(shù)為:max v,v為傳送帶過爐速度。

        最大過爐速度模型(模型Ⅱ):

        max v

        (14)

        其中,t1為第一次達(dá)到150 ℃對應(yīng)的時間,t2為第一次達(dá)到190 ℃對應(yīng)的時間, 為降溫至150 ℃時間, 為降溫至190 ℃時間。

        4.2? 最大過爐速度模型求解

        4.2.1? 搜索算法的構(gòu)建

        模型Ⅱ本質(zhì)上是一個單變量最優(yōu)化問題,提出一種搜索算法其過程為:

        (1)確定初始搜索區(qū)間[v1,v2],搜索步長為Δv;

        (2)在搜索區(qū)間[v1,v2]內(nèi)按步長Δv從大到小搜索,記錄第一個滿足約束條件的速度v;

        (3)更新搜索區(qū)間[v-Δv,v+Δv],并縮小搜索步長Δv,使其變?yōu)? ;

        (4)重復(fù)步驟(2)、(3),直到Δv

        上述算法流程圖如圖5所示。

        以上算法初值與參數(shù)的設(shè)定:

        傳送帶過爐速度調(diào)節(jié)范圍為v∈[65,100],初始搜索步長設(shè)定為1。參數(shù)n=10,即每次更新后步長變?yōu)樵瓉淼? ;參數(shù)p=10-4,即最終輸出速度的精度為10-4。

        4.2.2? 模型Ⅱ運(yùn)算結(jié)果

        模型Ⅱ通過以上搜索算法利用Matlab軟件編程求解,得到允許的最大傳送帶過爐速度為:

        vmax=75.76 cm/min

        5? 最優(yōu)爐溫系統(tǒng)設(shè)計(jì)

        理想的爐溫曲線還應(yīng)使超過217 ℃到峰值溫度所覆蓋的面積盡可能小,進(jìn)一步還希望以峰值溫度為中心線的兩側(cè)超過217 ℃的爐溫曲線應(yīng)盡量對稱。因此合理調(diào)節(jié)回焊爐各個溫區(qū)溫度,既能保證電路板的焊接質(zhì)量,又能使得效益最大化。

        5.1? 最優(yōu)爐溫系統(tǒng)模型建立

        首先要求在制程界限求出最優(yōu)爐溫曲線以及其對應(yīng)各溫區(qū)設(shè)定溫度、傳送帶過爐速度即可得對應(yīng)約束條件。各溫區(qū)最合適溫度表示為Ub=[Ub1,Ub2,Ub3,Ub4],其中Ub1∈[165,185]表示小溫區(qū)1~5的設(shè)定溫度,Ub2∈[185,205]表示小溫區(qū)6的設(shè)定溫度,Ub3∈[225,245]表示小溫區(qū)7的設(shè)定溫度,Ub4∈[245,265]表示小溫區(qū)8~9的設(shè)定溫度。最合適的傳送帶速度為vb。

        目標(biāo)函數(shù)的確定:

        設(shè)A為超過217 ℃到峰值溫度所覆蓋的面積,ta為第一次到達(dá)217 ℃的時間,tb為溫度到達(dá)峰值的時間,則覆蓋面積為:

        令Δt=,則原式表達(dá)為:

        (15)

        為描述兩邊曲線的對稱程度,令(Ti,Ui)與(ti,ui)分別是差分方程U′(t)上距離對稱軸(峰值所在的x軸垂線)最近的第i個點(diǎn)。當(dāng)Δx較小時,Ti與ti重合,要使得峰值溫度兩側(cè)超過217 ℃爐溫曲線盡可能對稱,即是使得|Ui-ui|足夠小,則:

        Ui-Ui+1≈ui-ui+1? ? ? ? ? ? ? ? ? ? ? ? ? ?(16)

        由式(16)可得 ,εi為第i個方程的偏差。方程偏差越小,說明中心線兩邊曲線的對稱性越好。根據(jù)方差的定義,該方程每一個偏差越小說明偏差的方差也越小。即研究對稱性目標(biāo)函數(shù)為:

        (17)

        為建立更合理的評價指標(biāo),本文將綜合考察曲線對稱性和面積大小兩個因素。

        將爐溫曲線超過217 ℃到峰值所覆蓋的面積映射到[0,1],消除量綱的影響,避免面積的數(shù)值遠(yuǎn)大于斜率偏差方差所造成的“大數(shù)吃小數(shù)”現(xiàn)象,所以面積大小比率為:

        (18)

        此外,由于焊接材料熔點(diǎn)為217 ℃左右,溫度在217 ℃到峰值溫度覆蓋的面積意義在于:若電路板溫度高于217 ℃時間過短或者溫度峰值太低會導(dǎo)致焊接材料融化不完全,焊接效果較差;若電路板溫度高于217 ℃時間過長或溫度峰值太高會導(dǎo)致焊接材料最終融化為液體四處流動,焊接效果較差。然而,使得超過217 ℃的爐溫曲線盡可能對稱的意義在于焊接材料不會發(fā)生驟冷或者驟熱,從而避免焊接材料由于降溫太慢或太快導(dǎo)致其最終形成的焊接材料韌性太差或者脆性太強(qiáng),實(shí)際材料的驟冷驟熱對材料的性能雖有影響,但影響并不是很大。權(quán)衡上述兩方面因素,材料溫度條件的重要性高于超過熔點(diǎn)部分爐溫曲線的對稱條件。所以衡量爐溫曲線優(yōu)劣程度的指標(biāo)Assess,取權(quán)重為ω1,ω2作為衡量這兩個因素的重要性,本文取ω1=0.7,ω2=0.3。最終,得到綜合評價指標(biāo):

        Assess=0.7η+0.3σ? ? ? ? ? ? ? ? ? ? ? ? ? (19)

        其中,面積大小比率η由式(18)求得,斜率偏差方差σ由式(17)求得。

        則計(jì)算超過217 ℃到峰值溫度的最小覆蓋面積的模型Ⅲ為:

        minAssess=0.7η+0.3σ

        (20)

        5.2? 最優(yōu)爐溫系統(tǒng)模型求解

        本文使用退火算法求解,其優(yōu)點(diǎn)在于其會按照一定概率跳出局部最優(yōu),以達(dá)到尋找全局最優(yōu)解的目的,并且算法收斂速度很快,用于求解多變量組合優(yōu)化問題十分有效。通過退火算法,運(yùn)用計(jì)算機(jī)求解,可得各溫區(qū)對應(yīng)下的最優(yōu)溫度和傳送帶過爐速度如表2所示。

        計(jì)算出斜率偏差方差σ=0.533 6,面積比率η=0.583 8,模型的評價指標(biāo)為Assess=0.568 7。對應(yīng)爐溫曲線如圖6所示,其中虛線是電路板已經(jīng)離開回焊爐或已經(jīng)通過回焊爐。

        6? 結(jié)? 論

        本文針對一種型號的回焊爐的爐溫曲線進(jìn)行研究,首先根據(jù)機(jī)理構(gòu)建了一維的爐溫曲線差分模型,使用一次實(shí)驗(yàn)測量的數(shù)據(jù)進(jìn)行模型各參數(shù)的調(diào)整后,模型的理論爐溫曲線與實(shí)際爐溫曲線基本重合,因此本文模型較為精確。在模型精確的前提下,進(jìn)一步使用搜索算法,蒙特卡洛法以及模擬退火算法全局尋優(yōu),尋求出了使得經(jīng)濟(jì)效益最大時對應(yīng)的試件的最大過爐速度以及各溫區(qū)應(yīng)被控制穩(wěn)定的溫度。

        因此,若要探究其他型號回焊爐的最優(yōu)爐溫系統(tǒng),只需測量一次該種型號回焊爐爐溫曲線,然后根據(jù)實(shí)驗(yàn)數(shù)據(jù)調(diào)整本文模型內(nèi)各參數(shù),本文模型即可適用于該種回焊爐。若制程界限有所不同,對應(yīng)調(diào)節(jié)本文模型Ⅱ與模型Ⅲ的約束條件,采用相同的研究方法即可設(shè)計(jì)出不同型號回焊爐對應(yīng)的最優(yōu)爐溫系統(tǒng)。

        參考文獻(xiàn):

        [1] 汪學(xué)軍.多溫區(qū)自動測控系統(tǒng)控制模型的建立與研究 [D].長沙:中南大學(xué),2007.

        [2] 張玉民.熱學(xué):第2版 [M].北京:科學(xué)出版社,2006.

        [3] 360百科.連續(xù)介質(zhì) [EB/OL].[2020-09-13]. https://baike.so.com/doc/897475-948678.html.

        [4] 高金剛.表面貼裝工藝生產(chǎn)線上回流焊曲線的優(yōu)化與控制 [D].上海:上海交通大學(xué),2007.

        作者簡介:方灝航(1999.05—),男,漢族,湖北黃岡人,本科在讀,研究方向:機(jī)械工程。

        人妻忍着娇喘被中进中出视频| 开心久久婷婷综合中文字幕| 免费观看成人欧美www色| 日本黄页网站免费观看| 亚洲AV一二三四区四色婷婷| 日本红怡院东京热加勒比| 97中文字幕精品一区二区三区| 最新中文字幕av无码不卡| 亚欧AV无码乱码在线观看性色| 中文字幕一区二区人妻痴汉电车| 女主播啪啪大秀免费观看| 五月综合激情婷婷六月| 伊人久久无码中文字幕| 国产亚洲日本人在线观看| 久久精品国产在热亚洲不卡| 免费无码不卡视频在线观看| 亚洲av日韩av永久无码色欲| 久久99精品久久久久久国产人妖| 蜜桃传媒免费在线观看| 九九久久99综合一区二区| 香色肉欲色综合| 亚洲综合精品在线观看中文字幕 | 亚洲精品综合久久国产二区| 色欲av永久无码精品无码蜜桃| 中文字幕人妻av一区二区| 韩国日本亚洲精品视频| 亚洲国产精品自拍成人| 成人精品视频一区二区| 国产精品美女| 一本色道亚州综合久久精品| 欧美老妇牲交videos| 亚洲老妇色熟女老太| 天堂AV无码AV毛片毛| 日本午夜精品一区二区三区| 两个人看的www免费视频中文 | 男女视频网站在线观看| 四虎影视成人永久免费观看视频 | 人人爽人人爽人人爽人人片av| 亚洲一区二区欧美色妞影院| 久久99精品久久只有精品| 色综合久久精品亚洲国产|