• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    剪切流場(chǎng)中聚乙烯結(jié)晶過(guò)程的建模與模擬

    2016-08-22 02:46:12阮春蕾劉春太
    化工學(xué)報(bào) 2016年5期
    關(guān)鍵詞:球晶結(jié)晶剪切

    阮春蕾,劉春太

    (1河南科技大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院,河南 洛陽(yáng) 471023;2橡塑模具國(guó)家工程研究中心,河南 鄭州 450001)

    ?

    剪切流場(chǎng)中聚乙烯結(jié)晶過(guò)程的建模與模擬

    阮春蕾1,2,劉春太2

    (1河南科技大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院,河南 洛陽(yáng) 471023;2橡塑模具國(guó)家工程研究中心,河南 鄭州 450001)

    基于Eder模型推導(dǎo)了剪切流場(chǎng)中球晶、串晶形態(tài)演化的數(shù)學(xué)模型,將第一法向應(yīng)力差作為串晶成核的驅(qū)動(dòng),并引入兩相懸浮模型描述體系,認(rèn)為其由無(wú)定形相和半結(jié)晶相組成,分別用FENE-P模型和剛性啞鈴模型描述?;谏鲜鰯?shù)學(xué)模型,分別構(gòu)造了捕捉球晶、串晶生長(zhǎng)的Monte Carlo法與體系控制方程求解的有限差分法,成功模擬了二維剪切流場(chǎng)中聚乙烯的結(jié)晶過(guò)程,給出了球晶、串晶的形態(tài)演化,分析了剪切時(shí)間、剪切速率對(duì)串晶成核密度、結(jié)晶速率、流體黏度等的影響。數(shù)值結(jié)果表明:所構(gòu)造的Monte Carlo法合理有效,不僅成功捕捉了晶體的生長(zhǎng)與碰撞,而且較為準(zhǔn)確地預(yù)測(cè)了結(jié)晶速率。此外,提高剪切時(shí)間或剪切速率,將增加串晶成核密度、提高結(jié)晶速率、使流體黏度突增的時(shí)間點(diǎn)提前。

    流動(dòng)誘導(dǎo)結(jié)晶;形態(tài)學(xué);動(dòng)力學(xué);Monte Carlo模擬;結(jié)晶

    DOI:10.11949/j.issn.0438-1157.20151472

    引 言

    聚合物的結(jié)晶過(guò)程及其形態(tài)演化是決定制品力學(xué)性能的關(guān)鍵因素[1]。廣義而言,聚合物必須經(jīng)過(guò)成型加工才能成為有用的制品。成型過(guò)程中,聚合物受到流場(chǎng)的作用,高分子鏈發(fā)生運(yùn)動(dòng),在低溫條件下形成凝聚態(tài)結(jié)構(gòu),并最終形成球晶、串晶等復(fù)雜結(jié)晶形態(tài)[1]。由于成型中以剪切流場(chǎng)居多,因此,本文將研究限定在簡(jiǎn)單剪切流場(chǎng)下,并以聚乙烯為代表,研究結(jié)晶形態(tài)與動(dòng)力學(xué),流體性質(zhì)等的改變,揭示外場(chǎng)因素與結(jié)晶過(guò)程間的相關(guān)規(guī)律,并為復(fù)雜成型條件下的結(jié)晶過(guò)程提供指導(dǎo)。

    流場(chǎng)中的聚合物結(jié)晶又被稱為流動(dòng)誘導(dǎo)結(jié)晶(flow induced crystallization,F(xiàn)IC)[1-2]。流動(dòng)誘導(dǎo)結(jié)晶的實(shí)驗(yàn)研究表明[1]:與靜態(tài)結(jié)晶相比,流場(chǎng)中的結(jié)晶發(fā)生了實(shí)質(zhì)性變化,不僅提高了結(jié)晶速率,而且改變了結(jié)晶形態(tài)(不同于靜態(tài)條件下的球晶,流場(chǎng)作用下可獲得高度異相的串晶結(jié)構(gòu))?;谏鲜鰧?shí)驗(yàn)結(jié)果,許多學(xué)者提出了不少關(guān)于流動(dòng)誘導(dǎo)結(jié)晶的動(dòng)力學(xué)模型,其建立大多基于Nakamura方程[3]和Kolmogorov方程[4]。如,Doufas等[5]、Tanner等[6-7]通過(guò)修正Nakamura方程得到了流動(dòng)誘導(dǎo)結(jié)晶動(dòng)力學(xué)模型;Koscher等[8]、Zheng等[9]認(rèn)為流動(dòng)提高了晶體的成核數(shù),并將其反映在Kolmogorov方程中,得到了流動(dòng)誘導(dǎo)結(jié)晶動(dòng)力學(xué)模型。該類方法的主要缺點(diǎn)在于不能提供較為詳細(xì)的微觀結(jié)構(gòu)信息。Eder等[10]首次提出了基于結(jié)晶形態(tài)的數(shù)學(xué)模型。他們把球晶看成是不斷長(zhǎng)大的球,串晶看成是不斷長(zhǎng)大的圓柱體,得到一系列微分方程組。Zuidema等[11]對(duì)其進(jìn)行了改進(jìn),用可恢復(fù)形變替代Eder模型中流場(chǎng)應(yīng)變率作為串晶成核的驅(qū)動(dòng)。他們的模型較為成功地反映了串晶的微觀信息,但相關(guān)研究并沒(méi)有顯式地采用模擬方法給出其成核、生長(zhǎng)與碰撞的演化細(xì)節(jié)。因此,仍然擺脫不了使用結(jié)晶動(dòng)力學(xué)模型。Guo等[12]引入了大分子鏈形變因子,他們認(rèn)為當(dāng)形變因子大于某臨界值時(shí),結(jié)晶按串晶形式生長(zhǎng);反之,按球晶形式生長(zhǎng)。他們利用該方法成功預(yù)測(cè)了iPP注塑制件的皮層厚度,但未得到結(jié)晶形態(tài)細(xì)節(jié)。

    國(guó)內(nèi)關(guān)于流動(dòng)誘導(dǎo)結(jié)晶數(shù)值模擬方面的研究也不多見(jiàn)。王錦燕等[13]、榮彥等[14]對(duì)剪切流場(chǎng)中的聚合物結(jié)晶進(jìn)行了模擬,定性分析了剪切對(duì)結(jié)晶過(guò)程及流體性質(zhì)等的影響。他們的研究并未涉及到結(jié)晶形態(tài)。Zhou等[15]采用Guo等[12]的方法對(duì)成型中的結(jié)晶形態(tài)進(jìn)行了顯示,區(qū)分了球晶和串晶,但相關(guān)研究并沒(méi)有將其與結(jié)晶動(dòng)力學(xué)建立聯(lián)系,而只是給出了結(jié)晶形態(tài)的可能表示。

    本文基于Eder模型推導(dǎo)了剪切流場(chǎng)中球晶、串晶的形態(tài)演化模型,并采用Monte Carlo法[16]成功模擬了二維剪切流場(chǎng)中聚乙烯的結(jié)晶過(guò)程,給出了球晶、串晶的形態(tài)演化,分析了剪切時(shí)間、剪切速率對(duì)串晶成核密度、結(jié)晶速率、流體黏度、體系應(yīng)力的影響。

    1 模型與算法

    1.1結(jié)晶動(dòng)力學(xué)模型

    Avrami方程[17]是聚合物結(jié)晶動(dòng)力學(xué)模型中應(yīng)用最廣的模型,它也可用于描述流場(chǎng)中的聚合物結(jié)晶,其具體表達(dá)式為

    其中,0φ為球晶的“虛幻體積”,0ψ為串晶的“虛幻體積”。Eder等[10]給出了0φ和0ψ的表達(dá)式,其中,描述三維球晶生長(zhǎng)的方程組為

    這里,Ns、Rtot、Stot、Vtot分別為球晶的晶核數(shù)、總體直徑長(zhǎng)度、總體表面積、總體體積;而a為成核速率相關(guān)的參數(shù),Gs為球晶的生長(zhǎng)速率。三維串晶生長(zhǎng)的方程組為

    這里,Ns-k、Ltot、S~tot、V~tot分別為串晶的晶核數(shù)、總體長(zhǎng)度、總體表面積、總體體積;而τn為與溫度相關(guān)的串晶成核率參數(shù),為串晶成核的驅(qū)動(dòng)參數(shù),γ.為應(yīng)變率,為擬合參數(shù);τl為與溫度相關(guān)的串晶軸向生長(zhǎng)率參數(shù),為串晶軸向生長(zhǎng)率的驅(qū)動(dòng)參數(shù),為擬合參數(shù),Gs-k,r為串晶徑向生長(zhǎng)速率。

    Eder模型中將球晶視為不斷長(zhǎng)大的球,而串晶視為不斷長(zhǎng)大的圓柱體。本文研究二維情況下的晶體生長(zhǎng),此時(shí),認(rèn)為球晶按圓形生長(zhǎng),而串晶按長(zhǎng)方形生長(zhǎng)。

    由式(2)可導(dǎo)出三維球晶的等價(jià)微分方程組,將其降為二維,則有

    從中可知,球晶的生長(zhǎng)由成核數(shù)sN及生長(zhǎng)速率sG所決定。采用如下的成核密度公式[2]

    其中,sN是過(guò)冷度TΔ的函數(shù),過(guò)冷度TΔ定義為為平衡熔點(diǎn),0N和?是經(jīng)驗(yàn)參數(shù)。球晶的生長(zhǎng)速率用Hoffman-Lauritzen方程表示[18]

    其中,G0、Kg為參考因子,U?為聚合物的分子活化能,Rg為氣體常數(shù),T∞=Tg-30,Tg為玻璃化轉(zhuǎn)變溫度,

    由式(3)可導(dǎo)出三維串晶的等價(jià)微分方程組,同樣,將其降為二維,有

    這里,認(rèn)為參數(shù)nτ=∞,lτ=∞[11],則有

    由式(8)可知,串晶的生長(zhǎng)由串晶成核密度Ns-k、串晶軸向生長(zhǎng)速率Gs-k,l以及串晶徑向生長(zhǎng)速率Gs-k,r決定。串晶軸向生長(zhǎng)速率Gs-k,l,可記為[10]

    一般情況下,認(rèn)為串晶徑向生長(zhǎng)速率Gs-k,r與球晶生長(zhǎng)速率Gs一致[11],即

    問(wèn)題的關(guān)鍵在于求解串晶成核密度。采用如下公式[8]

    其中,C為參數(shù),N1為體系的第一法向應(yīng)力差。

    1.2無(wú)定形相與半結(jié)晶相的模型

    式(11)中出現(xiàn)了第一法向應(yīng)力差,因此,必須對(duì)晶體所在體系做一下闡述。這里,采用Zheng等[9]的思想,認(rèn)為半結(jié)晶相懸浮于無(wú)定形相之中,無(wú)定形相采用FENE-P彈性啞鈴模型描述,而半結(jié)晶相采用棒狀啞鈴模型來(lái)描述。

    無(wú)定形相的FENE-P彈性啞鈴模型為[19]

    其中,C為大分子鏈構(gòu)型張量,λa(T)為分子松弛

    其中,,0aλ為溫度為0T時(shí)的分子松弛時(shí)間,ag/ER為擬合常數(shù)。

    無(wú)定形相對(duì)應(yīng)力的貢獻(xiàn)為[9,19]

    其中,aτ為無(wú)定形相應(yīng)力,n為分子密度,k為Botzmann常數(shù)。

    半結(jié)晶相中的分子不能拉伸,只能取向,故采用硬棒模型來(lái)描述[9],其動(dòng)力學(xué)方程為

    其中,為分子構(gòu)型張量,γ.為形變率張量,λsc(x,T)為棒狀啞鈴的松弛時(shí)間,采用如下公式[9]α為相對(duì)結(jié)晶度,A為經(jīng)驗(yàn)值。式(15)中的四階張量<>RRRR,采用如下的二次封閉近似[9]半結(jié)晶相對(duì)體系應(yīng)力的貢獻(xiàn)為[9]

    這里,ηsc(x,T)為半結(jié)晶相的黏度,與無(wú)定形相的黏度間滿足如下關(guān)系[9]

    至此,整個(gè)體系就完整描述出來(lái)。由式(14)、式(18)可得到體系的整體應(yīng)力表達(dá)式而式(11)中的第一法向應(yīng)力差N1=τ11-τ22由式(20)計(jì)算得到。

    1.3算法分析

    1.3.1Monte Carlo法本文采用Monte Carlo法[16]來(lái)捕捉結(jié)晶形態(tài)的演化。這里,研究對(duì)象設(shè)為空間中某塊特定區(qū)域上[0,1]mm×[0,1]mm的聚合物。給定溫度、剪切速率及剪切時(shí)間,并按式(5)、式(6)計(jì)算球晶成核密度Ns、球晶生長(zhǎng)速率Gs;按式(9)~式(11)計(jì)算串晶軸向生長(zhǎng)速率Gs-k,l、串晶徑向生長(zhǎng)速率Gs-k,r以及隨時(shí)間變化的串晶成核密度Ns-k。

    實(shí)施時(shí),先將區(qū)域剖分為Ntot個(gè)小單元,這里取Ntot=106,并對(duì)每個(gè)單元賦底色。在t=0時(shí)生成密度為Ns的球晶晶核和密度為Ns-k的串晶晶核,并對(duì)晶核賦顏色。在t時(shí)刻,計(jì)算球晶的半徑R=Gst,對(duì)于每一個(gè)球晶,生成足夠多的隨機(jī)點(diǎn),若該點(diǎn)的顏色為底色,則判斷其坐標(biāo)是否落入該球晶范圍內(nèi),若是則認(rèn)為該點(diǎn)被此球晶覆蓋,賦上該球晶的顏色;對(duì)于串晶,由于其成核與時(shí)間相關(guān),先由式(11)計(jì)算串晶的成核密度,隨機(jī)產(chǎn)生新生核的坐標(biāo)并賦值,同時(shí),計(jì)算串晶的長(zhǎng)度l=2Gs-k,l及半徑R=Gs-k,r,對(duì)于每一個(gè)串晶,生成足夠多的隨機(jī)點(diǎn),判斷其是否落入某串晶長(zhǎng)度及半徑覆蓋的長(zhǎng)方形中,若是則賦上該串晶的顏色。如此循環(huán),直到所有單元都被晶體覆蓋,則結(jié)束。需要指出的是,上述t~的計(jì)算與串晶的晶核生成時(shí)間相關(guān),從晶核生成后開(kāi)始計(jì)時(shí)。此外,相對(duì)結(jié)晶度的計(jì)算,這里采用式(21)來(lái)表示。

    其中,晶體單元數(shù)可根據(jù)單元顏色來(lái)統(tǒng)計(jì),而總體單元數(shù)為totN。Monte Carlo法的優(yōu)勢(shì)在于可以避開(kāi)使用結(jié)晶動(dòng)力學(xué)模型。常見(jiàn)的結(jié)晶動(dòng)力學(xué)模型,如Avrami方程和Kolmogorov方程,在處理球晶、串晶共混體系時(shí)遇到困難,很難確定形狀因子mC和n。而Monte Carlo法直接從形態(tài)演化中統(tǒng)計(jì)獲得可信的相對(duì)結(jié)晶度,避免了由結(jié)晶動(dòng)力學(xué)模型中參數(shù)的不確定而導(dǎo)致的誤差。

    1.3.2有限差分法無(wú)定形相與半結(jié)晶相方程的計(jì)算采用有限差分法。演化方程式(12)和式(15)采用時(shí)間一階向前進(jìn)行離散,即

    2 結(jié)果與討論

    2.1材料參數(shù)

    本文以聚乙烯為代表,考察結(jié)晶過(guò)程及流體性質(zhì)。模擬所用參數(shù)為:N0=17.4×106,?=0.155, G0=2.83×102m/s,,Kg=5.5×105K2,Tg=269 K,=2.69×10-7, C=105Pa-1·s-1·m-1,λa,0=4.00×10-2s,T0=476.15 K,Ea/Rg=5.602×103K,b=5,n= 1.26×1026/m3,k=1.38×10-23,β=9.2,β1=0.05,A=0.44。其中,結(jié)晶形態(tài)參數(shù)見(jiàn)文獻(xiàn)[8,11],無(wú)定形相及半結(jié)晶相參數(shù)見(jiàn)文獻(xiàn)[9]。

    2.2算法有效性驗(yàn)證

    為了驗(yàn)證算法的有效性,將模擬結(jié)果與Eder模型預(yù)測(cè)結(jié)果進(jìn)行比較。這里,假定球晶和串晶晶核瞬時(shí)產(chǎn)生,且有成核密度Ns=107m-2,Ns-k=107m-2,球晶按生長(zhǎng)速率Gs=10-6m·s-1進(jìn)行生長(zhǎng),串晶按徑向生長(zhǎng)速率Gs=10-6m·s-1及軸向生長(zhǎng)速率Gs-k,l=10-5m·s-1進(jìn)行生長(zhǎng)。圖1給出了模擬結(jié)果與理論模型預(yù)測(cè)所得相對(duì)結(jié)晶度的比較。由圖可知,模擬結(jié)果與理論預(yù)測(cè)值吻合得較好。模擬結(jié)果主要在結(jié)晶后期發(fā)生偏離,這可能是由于Eder模型的建立是基于Avrami方程的,而該方程在建模時(shí)并未考慮到晶體間的碰撞。這也是Avrami方程的主要局限性。

    圖1 數(shù)值解與Eder模型的比較Fig.1 Comparison of simulation result with Eder model

    2.3剪切時(shí)間的影響本節(jié)討論剪切時(shí)間的影響。為了便于分析,這里固定剪切速率γ.=10 s-1,溫度T=137℃。

    2.3.1結(jié)晶速率及形態(tài)的變化圖2為剪切時(shí)間t=0、1、5、10 s時(shí)所生成的串晶成核密度。其中,t=0代表不對(duì)流體進(jìn)行剪切,即靜態(tài)條件。施加剪切后,串晶的成核密度將有所增加;且剪切時(shí)間越大,串晶成核密度越大。當(dāng)剪切停止后,成核密度將不再隨時(shí)間的變化而變化。

    圖2 不同剪切時(shí)間下的串晶成核密度Fig.2 Nucleation density of shish-kebabs with different shear time

    圖3 不同剪切時(shí)間下的相對(duì)結(jié)晶度比較Fig.3 Relative crystallinity comparison with different shear time

    圖3為不同剪切時(shí)間下的相對(duì)結(jié)晶度演化。與靜態(tài)條件(t=0)相比,剪切時(shí)間越大,結(jié)晶速率越快。結(jié)晶的加速作用主要是由剪切誘導(dǎo)的串晶所貢獻(xiàn)的。由于流場(chǎng)的剪切作用,促使串晶晶核的產(chǎn)生,并提供軸向的生長(zhǎng)速率,使串晶成長(zhǎng),進(jìn)而加速結(jié)晶過(guò)程。在剪切時(shí)間較小的情況下,所提供的串晶晶核相對(duì)較少,其加速作用不甚明顯。文獻(xiàn)[9]也報(bào)道了大剪切時(shí)間對(duì)結(jié)晶速率的加速作用。一些學(xué)者[13-14]認(rèn)為這種加速不是無(wú)限制的,當(dāng)剪切時(shí)間達(dá)到某臨界值時(shí),加速作用會(huì)變得緩慢。

    圖4給出了剪切時(shí)間為0、5、10 s時(shí)的結(jié)晶形態(tài)比較。由圖可知,隨著剪切時(shí)間的增加,結(jié)晶形態(tài)發(fā)生較大變化:不僅提高了串晶成核密度,而且提高了串晶的各向異性。這種變化趨勢(shì)與Zhou等[15]的實(shí)驗(yàn)結(jié)果也是一致的。

    2.3.2流體黏度及應(yīng)力的變化圖5給出了剪切時(shí)間t=0、1、5、10 s時(shí)流體黏度的變化。黏度隨著時(shí)間的增加在某臨界值時(shí)會(huì)發(fā)生突增現(xiàn)象。這是由相對(duì)結(jié)晶度引起的。由式(19)知,ηsc(x,T)=,當(dāng)α接近常數(shù)A時(shí),ηsc將趨向于無(wú)窮;故體系黏度急劇變化。此外,剪切時(shí)間越大,流體黏度越早發(fā)生突增。這與Zheng等[9]的結(jié)果一致。

    圖6為剪切時(shí)間t=10 s時(shí)的剪切應(yīng)力和第一法向應(yīng)力差。應(yīng)力值在前期變化較為緩慢,而到臨界值附近變化劇烈,該臨界值對(duì)應(yīng)于相對(duì)結(jié)晶度。當(dāng)結(jié)晶完成時(shí),這些應(yīng)力將鎖在材料內(nèi)部,形成殘余應(yīng)力[20]。

    2.4剪切速率的影響

    本節(jié)討論剪切速率的影響。為了便于分析,這里固定剪切時(shí)間t=10 s,溫度T=137℃。

    2.4.1結(jié)晶速率及形態(tài)的變化圖7給出了剪切速率0γ=.、1、5、10 s-1時(shí)串晶成核密度的比較。其中,γ=.0代表不對(duì)流體進(jìn)行剪切,即靜態(tài)條件。由圖可知,剪切速率越大,串晶成核密度越大。當(dāng)剪切停止后,成核密度將不再隨時(shí)間的變化而變化。

    圖8為剪切速率γ=.0、1、5、10 s-1時(shí)的相對(duì)結(jié)晶度發(fā)展。由圖可知,由于剪切作用的存在,結(jié)晶速率將比靜態(tài)條件下(0γ=.)的結(jié)晶速率有所增加,但這種增加在小剪切速率下(γ=.1 s-1)較微弱。同樣,結(jié)晶的這種加速作用也是由剪切誘導(dǎo)的串晶所貢獻(xiàn)的。

    結(jié)晶形態(tài)的演化受剪切速率的影響與受剪切時(shí)間的影響類似。隨著剪切速率的增加,串晶成核密度將有所增加,并且串晶的各向異性也更為明顯。這里將不再作圖顯示,圖形與圖4類似。

    2.4.2流體黏度的變化圖9為不同剪切速率下流體黏度的變化。流體隨時(shí)間的變化首先較為緩慢,但到某臨界值時(shí)發(fā)生突增,這也是由結(jié)晶度所引起的。此外,剪切速率越大,流體黏度突增的時(shí)間點(diǎn)越早。這與Zheng等[9]的結(jié)果也是一致的。

    圖4 不同剪切時(shí)間下結(jié)晶形態(tài)的比較Fig.4 Morphology comparison with different shear time

    3 結(jié) 論

    本文給出了一個(gè)二維簡(jiǎn)單剪切流場(chǎng)下聚乙烯結(jié)晶過(guò)程的建模與模擬?;贓der模型推導(dǎo)了球晶、串晶的形態(tài)演化模型;并用Zheng的兩相懸浮模型描述聚合物流體,將體系的第一法向應(yīng)力差作為串晶成核的驅(qū)動(dòng);與此同時(shí),構(gòu)建了模擬球晶、串晶形態(tài)演化的Monte Carlo法,分析了剪切時(shí)間、剪切速率對(duì)結(jié)晶速率、結(jié)晶形態(tài)、體系黏度等的影響。得出如下結(jié)論。

    (1)本文所建立的Monte Carlo法合理有效,成功區(qū)分了晶體生長(zhǎng)模式,捕捉了晶體生長(zhǎng)前沿,較為直觀地給出了晶體的演化發(fā)展;同時(shí),統(tǒng)計(jì)獲得了可信的相對(duì)結(jié)晶度。該方法的優(yōu)勢(shì)在于成功避開(kāi)了使用結(jié)晶動(dòng)力學(xué)模型,避免由結(jié)晶動(dòng)力學(xué)模型中參數(shù)的不確定而導(dǎo)致預(yù)測(cè)的失真。

    (2)剪切時(shí)間對(duì)結(jié)晶速率、結(jié)晶形態(tài)影響重大;剪切時(shí)間越長(zhǎng)結(jié)晶速率越快,串晶對(duì)結(jié)晶形態(tài)的影響越明顯,各向異性更為顯著。剪切時(shí)間對(duì)流體黏度也有顯著影響,這種影響主要體現(xiàn)在對(duì)結(jié)晶度的依賴上,當(dāng)結(jié)晶度接近于某臨界值時(shí),體系黏度將發(fā)生突增現(xiàn)象。

    (3)剪切速率對(duì)結(jié)晶速率、結(jié)晶形態(tài)、體系黏度等的影響與剪切時(shí)間的影響類似;增加剪切速率能提高結(jié)晶速率、提高串晶的各向異性以及使體系黏度發(fā)生突增的時(shí)間點(diǎn)提前。

    圖5 不同剪切時(shí)間下流體黏度的變化Fig.5 Evolution of melt viscosity with different shear time

    圖6 剪切應(yīng)力和第一法向應(yīng)力差(剪切時(shí)間t=10 s)Fig.6 Shear stress and first normal stress (shear time t=10 s)

    圖7 不同剪切速率下的串晶成核密度Fig.7 Nucleation density of shish-kebabs with different shear rate

    圖8 不同剪切速率下的相對(duì)結(jié)晶度比較Fig.8 Relative crystallinity comparison with different shear rate

    圖9 不同剪切速率下流體黏度的變化Fig.9 Evolution of melt viscosity with different shear rate

    符號(hào)說(shuō)明

    b——啞鈴的有限拉伸參數(shù)

    C——彈性啞鈴構(gòu)型張量,m

    Gs,Gs-k,l,Gs-k,r——分別為球晶、串晶軸向、串晶徑

    向生長(zhǎng)速率,m·s-1

    Ns,Ns-k——分別為球晶、串晶成核密度,m-2

    N1——第一法向應(yīng)力差,N·m-2

    R——?jiǎng)傂詥♀徣∠蛳蛄浚琺

    T——熱力學(xué)溫度,K

    t——時(shí)間,s

    Δt——時(shí)間步長(zhǎng),s

    α——相對(duì)結(jié)晶度

    ηa,ηsc——分別為無(wú)定形相、半結(jié)晶相黏度,Pa·sa

    λ,scλ——分別為無(wú)定形相、半結(jié)晶相分子

    松弛因子時(shí)間,s

    τ,aτ,scτ——分別為總體、無(wú)定形相、半結(jié)晶

    相應(yīng)力,N·m-2

    ?——梯度算子

    〈〉——系綜平均

    下角標(biāo)

    a ——無(wú)定形相

    s ——球晶

    s-k ——串晶

    sc ——半結(jié)晶相

    References

    [1]KENNEDY P, ZHENG R. Crystallization//Flow Analysis of Injection Molds [M]. Munich: MunichCarl Hanser Verlag, 2013: 141-153.

    [2]ZHE M, LUIGI B, GIUSEPPE P, et al. Flow induced crystallization in isotactic polypropylene during and after flow [J]. Polym., 2014, 55:6140-6151.

    [3]NAKAMURA K, WATANABE T, KATAYAMA K, et al. Some aspects of nonisothermal crystallization of polymers (Ⅰ):Relationship between crystallization temperature, crystallinity, and cooling conditions [J]. J. Appl. Polym. Sci., 1972, 16: 1077-1091.

    [4]KOLMOGOROV A N. On the statistic of crystallization development in metals [J]. Bull. Akad. Sci. USSR, Class Sci., Math. Nat., 1937, 1:355-359.

    [5]DOUFAS A K, MCHUGH A J, MILLER C. Simulation of melt spinning including flow-induced crystallization (1): Model development and predictions [J]. J. Non-Newtonian Fluid Mech., 2000, 92: 27-66.

    [6]TANNER R I, QI F Z. A comparison of some models for describing polymer crystallization at low deformation rates [J]. J. Non-Newtonian Fluid Mech., 2005, 127: 131-141.

    [7]WO D L, TANNER R I, FLETCHER D F. A numerical treatment of crystallization in tube flow [J]. Prog. Polym. Sci., 2012, 56: 1356-1366.

    [8]KOSCHER E, FULCHIRON R. Influence of shear on polypropylene crystallization: morphology development and kinetics [J]. Polym.,2002, 43: 6931-6942.

    [9]ZHENG R, KENNEDY P. A model for post-flow induced crystallization: general equations and predictions [J]. J. Rheol., 2004, 48: 823-842.

    [10]EDER G, JANESCHITZ-KRIEGL H, LIEDAUER S. Crystallization processer in quiescent and moving polymer melts under heat transfer conditions [J]. Prog. Polym. Sci., 1990, 15: 629-714.

    [11]ZUIDEMA H, PETERS G W M, MEIJER H E H. Development and validation of a recoverable strain-based model for flow induced crystallization of polymers [J]. Macromol. Theory Simul., 2001, 10:447-460.

    [12]GUO X, ISAYEV A I, GUO L. Crystallinity and microstructure in injection modelings of isotactic polypropylenes (Ⅰ): a new approach to modeling and model parameters [J]. Polym. Eng. Sci., 1999, 39(10): 2096-2114.

    [13]王錦燕, 陳靜波, 劉春太, 等. 聚合物流動(dòng)誘導(dǎo)結(jié)晶數(shù)值模擬 [J].化工學(xué)報(bào), 2011, 62 (4): 1150-1156. WANG J Y, CHEN J B, LIU C T, et al. Numerical simulation of shear-induced crystallization of polymer [J]. CIESC Journal, 2011, 62(4): 1150-1156.

    [14]榮彥, 賀惠萍, 曹偉, 等. 基于兩相模型的聚合物流動(dòng)誘導(dǎo)結(jié)晶數(shù)值模擬 [J]. 化工學(xué)報(bào), 2012, 63 (7): 2252-2257. RONG Y, HE H P, CAO W, et al. Numercial simulation of flow-induced crystallization of polymer based on two-phase model [J]. CIESC Journal, 2012, 63 (7): 2252-2257.

    [15]ZHOU Y G, TURNG L, SHEN C Y. Modeling and prediction of morphology and crystallinity for cylindrical-shaped crystals during polymer processsing [J]. Polym. Eng. Sci., 2010: 1-10.

    [16]楊玉良, 張紅東. 高分子科學(xué)中的Monte Carlo方法 [M]. 上海:復(fù)旦大學(xué)出版社, 1993. YANG Y L, ZHANG H D. Monte Carlo Method in Polymer Science[M]. Shanghai: Fudan University Press, 1993.

    [17]AVRAMI M. Kinetics of phase change (Ⅰ): General theory [J]. J. Chem. Phys., 1939, 7: 1103-1112.

    [18]HOFFMAN J D, MILLER R L. Kinetics of crystallization from the melt and chain folding in polyethylene fractions revisited: theory and experiment [J]. Polym., 1997, 38: 3151-3212.

    [19]NASIR A, MUHAMMAD A J, MUHAMMAD S. Theoretical analysis of the exiting thickness of sheets in the calendaring of FENE-P fluid [J]. J. Non-Newtonian Fluid Mech., 2015, 225: 28-36.

    [20]ABOUHAMZEH M, SINKE J, JANSEN K, et al. Closed form expression for residual stresses and warpage during cure of composite laminates [J]. Compos. Struct., 2015, 133: 902-910.

    Numerical simulation of morphology and kinetics of polyethylene in shear flow

    RUAN Chunlei1,2, LIU Chuntai2
    (1School of Mathematics and Statistics, Henan University of Science and Technology, Luoyang 471023, Henan, China;2National Engineering Research Center for Advanced Polymer Processing Technology, Zhengzhou 450001, Henan, China)

    The mathematical model of morphology evolution of spherulites and shish-kebabs in the shear flow is deduced based on the Eder model. The model considers that the nucleation and growth of spherulites are determined by the static temperature, while the nucleation and growth of shish-kebabs are determined by the flow,which is depended on the first normal stress difference of the system and the shear rate respectively. In order to calculate the nucleation density of the shish-kebabs, the two-phase suspension model of Zheng is introduced. The model treats the stress as the combination of the amorphous phase and semi-crystalline phase. The amorphous phase is described by FENE-P model while the semi-crystalline phase is depicted by a rigid dumbbell model. Based on the mathematical model, the Monte Carlo method and the finite difference method are constructed,respectively. The former is to capture the crystal growth while the latter is to calculate the equation of the system. By using these methods, the 2D crystallization in the shear flow is simulated. The evolution of spherulites and shish-kebabs is given. Also, the effects of shear time and shear rate on the nucleation density of shish-kebabs,crystallization rate, the viscosity of the fluid and the system stress are discussed. Numerical results show that the Monte Carlo method is valid which not only captures the morphology evolution of crystals successfully, but alsopredicts the crystallization rate well. In addition, the increase of the shear time or the shear rate will increase the nucleation density of shish-kebabs and the crystallization rate.

    date: 2015-09-21.

    RUAN Chunlei, ruanchunlei622@mail.nwpu. edu.cn

    supported by the National Natural Science Foundation in China (11402078, 51375148).

    flow induced crystallization; morphology; kinetics; Monte Carlo simulation; crystallization

    O631.+3

    A

    0438—1157(2016)05—2144—08

    2015-09-21收到初稿,2016-01-11收到修改稿。

    聯(lián)系人及第一作者:阮春蕾(1983—),女,博士,副教授。

    國(guó)家自然科學(xué)基金項(xiàng)目(11402078,51375148)。

    猜你喜歡
    球晶結(jié)晶剪切
    球晶的結(jié)構(gòu)和形成機(jī)理
    “長(zhǎng)大”的結(jié)晶
    Concurrence of Merkel Cell Carcinoma and Squamous Cell Carcinoma in A Patient with Generalized Actinic Keratosis: A Case Report
    淀粉球晶的制備及其理化性質(zhì)
    寬厚板剪切線控制系統(tǒng)改進(jìn)
    山東冶金(2018年5期)2018-11-22 05:12:46
    混凝土短梁斜向開(kāi)裂后的有效剪切剛度與變形
    共聚甲醛的自成核結(jié)晶行為
    土-混凝土接觸面剪切破壞模式分析
    sPS/PBA-aPS共混物的結(jié)晶與熔融行為
    Xanomeline新型衍生物SBG-PK-014促進(jìn)APPsw的α-剪切
    404 Not Found

    404 Not Found


    nginx
    天堂中文最新版在线下载| 国产黄色免费在线视频| 久久精品国产a三级三级三级| 狠狠婷婷综合久久久久久88av| 自线自在国产av| 午夜老司机福利剧场| 久久午夜综合久久蜜桃| 我的女老师完整版在线观看| 亚洲精品,欧美精品| 九草在线视频观看| 国产高清三级在线| 国产老妇伦熟女老妇高清| 一级a做视频免费观看| 亚洲精品乱久久久久久| 久久久久久久久久人人人人人人| 亚洲成人一二三区av| 午夜福利影视在线免费观看| 亚洲情色 制服丝袜| 26uuu在线亚洲综合色| 观看av在线不卡| 久久人人爽人人爽人人片va| 亚洲欧美日韩卡通动漫| 国产乱来视频区| 97在线人人人人妻| 人人妻人人澡人人爽人人夜夜| 久久精品久久久久久久性| 欧美亚洲日本最大视频资源| 18+在线观看网站| 免费女性裸体啪啪无遮挡网站| 又大又黄又爽视频免费| 久久99热这里只频精品6学生| 97人妻天天添夜夜摸| 少妇精品久久久久久久| 九九在线视频观看精品| 99热网站在线观看| 最后的刺客免费高清国语| 性色av一级| 日韩欧美一区视频在线观看| 久久免费观看电影| 欧美精品av麻豆av| 最后的刺客免费高清国语| 99久久人妻综合| 久久久久国产网址| 人人妻人人澡人人看| 久久热在线av| 国产黄色视频一区二区在线观看| 两个人免费观看高清视频| 最近最新中文字幕大全免费视频 | 人妻系列 视频| 久久这里有精品视频免费| 精品熟女少妇av免费看| 亚洲综合色网址| 欧美人与性动交α欧美软件 | 久久久国产一区二区| 咕卡用的链子| 在线观看一区二区三区激情| 大话2 男鬼变身卡| 免费av中文字幕在线| www日本在线高清视频| av在线app专区| 女性生殖器流出的白浆| 日韩大片免费观看网站| 春色校园在线视频观看| 亚洲伊人色综图| 国产在线免费精品| 国产成人91sexporn| 久久韩国三级中文字幕| 女人被躁到高潮嗷嗷叫费观| 看免费成人av毛片| 精品一区二区免费观看| 国产亚洲精品久久久com| 国产色爽女视频免费观看| 大陆偷拍与自拍| 极品人妻少妇av视频| 亚洲 欧美一区二区三区| 啦啦啦啦在线视频资源| 欧美人与善性xxx| 超色免费av| www.熟女人妻精品国产 | 国产一区二区三区av在线| 黑人欧美特级aaaaaa片| 国产一级毛片在线| 国产免费现黄频在线看| 成人国产av品久久久| 亚洲国产av影院在线观看| 美女福利国产在线| 精品亚洲乱码少妇综合久久| 国产免费福利视频在线观看| 国产激情久久老熟女| 老司机影院成人| 国产一区有黄有色的免费视频| 看十八女毛片水多多多| 久久久久精品久久久久真实原创| 热re99久久国产66热| 日韩制服丝袜自拍偷拍| 91精品三级在线观看| 成年美女黄网站色视频大全免费| 欧美xxⅹ黑人| 国产av国产精品国产| 曰老女人黄片| 婷婷成人精品国产| tube8黄色片| 日本黄色日本黄色录像| 伊人亚洲综合成人网| 男人操女人黄网站| 色吧在线观看| 各种免费的搞黄视频| 91aial.com中文字幕在线观看| 美女主播在线视频| 一级毛片 在线播放| 成人国产av品久久久| 制服诱惑二区| av福利片在线| 91国产中文字幕| 少妇的逼好多水| 精品久久久久久电影网| 精品一区二区三区四区五区乱码 | 欧美精品av麻豆av| 中文乱码字字幕精品一区二区三区| 亚洲精品第二区| a级毛色黄片| 亚洲天堂av无毛| 亚洲精品乱码久久久久久按摩| 亚洲精品乱码久久久久久按摩| 一区在线观看完整版| 免费在线观看黄色视频的| 国产精品成人在线| 亚洲国产av影院在线观看| 亚洲国产av影院在线观看| 久久久久视频综合| 久久久精品94久久精品| 视频中文字幕在线观看| 国产精品久久久久久av不卡| 午夜免费鲁丝| av卡一久久| 午夜激情av网站| 涩涩av久久男人的天堂| 伊人亚洲综合成人网| 成人18禁高潮啪啪吃奶动态图| 欧美亚洲日本最大视频资源| 免费观看在线日韩| av天堂久久9| 国产69精品久久久久777片| 亚洲久久久国产精品| 一区在线观看完整版| 精品少妇久久久久久888优播| 免费黄网站久久成人精品| 成人亚洲精品一区在线观看| 九色成人免费人妻av| 在线免费观看不下载黄p国产| 亚洲在久久综合| 欧美日韩精品成人综合77777| 精品少妇内射三级| 午夜免费观看性视频| av在线观看视频网站免费| 免费观看性生交大片5| 乱人伦中国视频| 午夜福利在线观看免费完整高清在| 精品国产一区二区久久| 91精品国产国语对白视频| 国产成人a∨麻豆精品| 国产国语露脸激情在线看| 少妇的丰满在线观看| 久久青草综合色| 国产一级毛片在线| 青春草视频在线免费观看| 久久久久国产网址| 中国美白少妇内射xxxbb| 99re6热这里在线精品视频| 免费人妻精品一区二区三区视频| 一级毛片电影观看| 在线免费观看不下载黄p国产| 亚洲第一av免费看| 天堂中文最新版在线下载| 视频区图区小说| 18禁裸乳无遮挡动漫免费视频| 国产成人a∨麻豆精品| 国产一级毛片在线| 蜜桃在线观看..| 国产av码专区亚洲av| 草草在线视频免费看| 在线观看美女被高潮喷水网站| 午夜91福利影院| 欧美成人午夜精品| 久久久久久人妻| 看免费成人av毛片| 日本vs欧美在线观看视频| 99久国产av精品国产电影| 涩涩av久久男人的天堂| 免费观看av网站的网址| 搡女人真爽免费视频火全软件| 亚洲 欧美一区二区三区| 啦啦啦啦在线视频资源| 一本一本久久a久久精品综合妖精 国产伦在线观看视频一区 | 建设人人有责人人尽责人人享有的| 成年女人在线观看亚洲视频| 成人毛片60女人毛片免费| 精品国产乱码久久久久久小说| 中文字幕最新亚洲高清| 日本av手机在线免费观看| 日韩 亚洲 欧美在线| 亚洲伊人色综图| 国语对白做爰xxxⅹ性视频网站| 久久人人爽人人片av| 精品人妻在线不人妻| 成人亚洲欧美一区二区av| 亚洲av日韩在线播放| 99热网站在线观看| 国产免费视频播放在线视频| 美国免费a级毛片| 丝瓜视频免费看黄片| 老熟女久久久| 人妻 亚洲 视频| 99国产综合亚洲精品| 香蕉丝袜av| 精品熟女少妇av免费看| 日本色播在线视频| 国产极品粉嫩免费观看在线| 亚洲国产欧美在线一区| √禁漫天堂资源中文www| 黄片播放在线免费| 亚洲婷婷狠狠爱综合网| 欧美日韩av久久| 精品亚洲成a人片在线观看| 久久久久久久亚洲中文字幕| 日韩视频在线欧美| 欧美激情 高清一区二区三区| 免费观看在线日韩| 91精品三级在线观看| 亚洲精品美女久久av网站| 黄色一级大片看看| 中文字幕另类日韩欧美亚洲嫩草| 日本与韩国留学比较| 欧美日韩视频精品一区| 中文字幕人妻熟女乱码| 少妇被粗大猛烈的视频| 男女高潮啪啪啪动态图| 日本av手机在线免费观看| 久久午夜福利片| 亚洲精品中文字幕在线视频| 18禁国产床啪视频网站| 99视频精品全部免费 在线| 欧美人与性动交α欧美精品济南到 | 伊人亚洲综合成人网| 亚洲伊人久久精品综合| 黄色怎么调成土黄色| 9191精品国产免费久久| 97在线人人人人妻| 亚洲在久久综合| 久久精品人人爽人人爽视色| 久久久久久久久久久久大奶| 中文欧美无线码| 精品午夜福利在线看| 日本av手机在线免费观看| 国语对白做爰xxxⅹ性视频网站| 韩国av在线不卡| 久久久欧美国产精品| 天堂中文最新版在线下载| 日韩中字成人| 国产精品一区二区在线不卡| 黄色怎么调成土黄色| 久久精品夜色国产| 美女中出高潮动态图| 18禁国产床啪视频网站| 色5月婷婷丁香| 亚洲图色成人| www.av在线官网国产| 国产男女超爽视频在线观看| 亚洲成av片中文字幕在线观看 | 日韩三级伦理在线观看| av天堂久久9| 韩国高清视频一区二区三区| 午夜免费观看性视频| 国产精品久久久久久久电影| 久久久国产一区二区| www日本在线高清视频| 亚洲三级黄色毛片| 九九在线视频观看精品| 七月丁香在线播放| 水蜜桃什么品种好| 久久久久久人人人人人| h视频一区二区三区| 一级,二级,三级黄色视频| 亚洲伊人久久精品综合| xxx大片免费视频| 一级毛片电影观看| 日日摸夜夜添夜夜爱| 久久久久精品性色| 国产乱来视频区| 高清毛片免费看| 精品视频人人做人人爽| 97超碰精品成人国产| 在线精品无人区一区二区三| 久久久久久久久久人人人人人人| 黑丝袜美女国产一区| 欧美bdsm另类| 日韩三级伦理在线观看| 9色porny在线观看| 在现免费观看毛片| 亚洲国产毛片av蜜桃av| 国产精品久久久久久久电影| 欧美精品一区二区免费开放| 免费av中文字幕在线| 97人妻天天添夜夜摸| 伦精品一区二区三区| 日本色播在线视频| √禁漫天堂资源中文www| 日韩 亚洲 欧美在线| 高清在线视频一区二区三区| 久久这里有精品视频免费| 色网站视频免费| 婷婷色麻豆天堂久久| 青春草亚洲视频在线观看| 丝袜脚勾引网站| 国产不卡av网站在线观看| 99国产精品免费福利视频| 亚洲第一区二区三区不卡| 国产一区二区在线观看日韩| 精品少妇久久久久久888优播| 亚洲欧美中文字幕日韩二区| 三级国产精品片| 日韩制服丝袜自拍偷拍| 久久av网站| 亚洲第一区二区三区不卡| 亚洲国产av影院在线观看| 中文字幕av电影在线播放| 亚洲精品一区蜜桃| 女人久久www免费人成看片| 亚洲国产看品久久| 国产男女超爽视频在线观看| 亚洲美女黄色视频免费看| 日日撸夜夜添| 欧美少妇被猛烈插入视频| 好男人视频免费观看在线| 另类精品久久| 夜夜爽夜夜爽视频| kizo精华| 亚洲精品乱码久久久久久按摩| 99久久中文字幕三级久久日本| 国产日韩一区二区三区精品不卡| 在线 av 中文字幕| 精品人妻一区二区三区麻豆| 人人妻人人澡人人爽人人夜夜| 国产欧美日韩综合在线一区二区| 一级毛片我不卡| 伊人久久国产一区二区| 女性生殖器流出的白浆| www.色视频.com| 人人澡人人妻人| 久久久久久久大尺度免费视频| 美女中出高潮动态图| 日韩,欧美,国产一区二区三区| 精品福利永久在线观看| 亚洲av福利一区| 亚洲婷婷狠狠爱综合网| 啦啦啦啦在线视频资源| 秋霞伦理黄片| 满18在线观看网站| 日韩中字成人| 午夜免费鲁丝| 纯流量卡能插随身wifi吗| 国产日韩欧美在线精品| 精品少妇黑人巨大在线播放| 男女边摸边吃奶| 成人国产av品久久久| 中文字幕亚洲精品专区| 亚洲三级黄色毛片| 精品熟女少妇av免费看| 一本—道久久a久久精品蜜桃钙片| a级毛片黄视频| 久久国产精品男人的天堂亚洲 | 欧美精品一区二区大全| 自拍欧美九色日韩亚洲蝌蚪91| 视频中文字幕在线观看| 91精品国产国语对白视频| 如何舔出高潮| 亚洲精品一二三| 国产女主播在线喷水免费视频网站| 久久人人爽人人爽人人片va| 国产免费视频播放在线视频| 丝袜在线中文字幕| 午夜福利乱码中文字幕| 男女边摸边吃奶| 国产熟女午夜一区二区三区| 亚洲综合色网址| 美国免费a级毛片| 精品人妻在线不人妻| 日本-黄色视频高清免费观看| 七月丁香在线播放| 色婷婷久久久亚洲欧美| 少妇精品久久久久久久| 免费观看在线日韩| 男女边吃奶边做爰视频| 欧美人与性动交α欧美软件 | 高清视频免费观看一区二区| 男人操女人黄网站| 国产精品一区二区在线不卡| 国产 精品1| 看免费成人av毛片| 成人黄色视频免费在线看| 中国国产av一级| 啦啦啦啦在线视频资源| 亚洲成人av在线免费| 久久久久久伊人网av| 久久久久国产精品人妻一区二区| 午夜福利视频在线观看免费| 香蕉丝袜av| 伦理电影免费视频| 国产精品久久久久久精品电影小说| 日韩电影二区| 亚洲精品,欧美精品| 综合色丁香网| 女的被弄到高潮叫床怎么办| 90打野战视频偷拍视频| 如日韩欧美国产精品一区二区三区| 久久精品国产综合久久久 | 免费人妻精品一区二区三区视频| 男女高潮啪啪啪动态图| 建设人人有责人人尽责人人享有的| 午夜免费男女啪啪视频观看| av.在线天堂| 国产精品欧美亚洲77777| 久久久久网色| 精品国产一区二区三区四区第35| 欧美精品一区二区免费开放| 久久99热这里只频精品6学生| 内地一区二区视频在线| 91精品三级在线观看| 男人操女人黄网站| 黄片播放在线免费| 97在线人人人人妻| 亚洲第一区二区三区不卡| 亚洲av电影在线观看一区二区三区| 青春草国产在线视频| 国产在线视频一区二区| 亚洲综合色网址| 最近的中文字幕免费完整| 久久av网站| 亚洲精品色激情综合| 最近2019中文字幕mv第一页| 少妇高潮的动态图| 亚洲国产精品专区欧美| 亚洲国产色片| 国产成人精品在线电影| 国产一区亚洲一区在线观看| 午夜福利影视在线免费观看| 人人妻人人爽人人添夜夜欢视频| 国产精品一国产av| 免费观看无遮挡的男女| 亚洲三级黄色毛片| 永久免费av网站大全| 最近中文字幕高清免费大全6| 9热在线视频观看99| 十八禁高潮呻吟视频| 久久久久久久久久成人| 又大又黄又爽视频免费| 毛片一级片免费看久久久久| 伊人久久国产一区二区| 午夜av观看不卡| 在线 av 中文字幕| 在线天堂最新版资源| xxxhd国产人妻xxx| 18禁国产床啪视频网站| 亚洲av男天堂| 国产深夜福利视频在线观看| 91午夜精品亚洲一区二区三区| 有码 亚洲区| 久久人人97超碰香蕉20202| 国产欧美日韩综合在线一区二区| 久久久久网色| 久久这里有精品视频免费| 美女国产视频在线观看| 青春草国产在线视频| 久久精品国产自在天天线| 欧美3d第一页| 女人精品久久久久毛片| 91成人精品电影| 亚洲成国产人片在线观看| 三级国产精品片| 精品福利永久在线观看| 日日摸夜夜添夜夜爱| 在线精品无人区一区二区三| 国产激情久久老熟女| 国产日韩欧美亚洲二区| 91午夜精品亚洲一区二区三区| 午夜91福利影院| 久久婷婷青草| 亚洲成国产人片在线观看| av线在线观看网站| 国产69精品久久久久777片| 亚洲欧洲国产日韩| 丝袜喷水一区| 777米奇影视久久| 欧美老熟妇乱子伦牲交| 亚洲情色 制服丝袜| 波野结衣二区三区在线| 亚洲国产最新在线播放| 99精国产麻豆久久婷婷| 亚洲婷婷狠狠爱综合网| 最近中文字幕2019免费版| 一二三四在线观看免费中文在 | 在线天堂最新版资源| 五月天丁香电影| 久久久国产精品麻豆| 高清欧美精品videossex| 日产精品乱码卡一卡2卡三| 日韩免费高清中文字幕av| 人人妻人人澡人人爽人人夜夜| 午夜91福利影院| 1024视频免费在线观看| 黄网站色视频无遮挡免费观看| 日韩制服骚丝袜av| 精品人妻一区二区三区麻豆| av天堂久久9| 18+在线观看网站| av在线播放精品| 亚洲精品美女久久av网站| 久久国产亚洲av麻豆专区| 男女免费视频国产| 国产成人aa在线观看| 国产极品粉嫩免费观看在线| 国产在线一区二区三区精| 一级毛片电影观看| 岛国毛片在线播放| 99久久综合免费| 男女啪啪激烈高潮av片| 亚洲四区av| 国产 一区精品| 久久99精品国语久久久| 女人被躁到高潮嗷嗷叫费观| 午夜激情久久久久久久| 免费大片黄手机在线观看| 国产亚洲精品久久久com| 欧美精品一区二区免费开放| 亚洲国产毛片av蜜桃av| 午夜免费鲁丝| 亚洲一码二码三码区别大吗| 亚洲国产成人一精品久久久| av又黄又爽大尺度在线免费看| 国产毛片在线视频| 国产亚洲精品久久久com| 两性夫妻黄色片 | 午夜日本视频在线| 亚洲少妇的诱惑av| 久久久久人妻精品一区果冻| 国产av码专区亚洲av| 韩国高清视频一区二区三区| 国产精品国产三级国产av玫瑰| 一本一本久久a久久精品综合妖精 国产伦在线观看视频一区 | 亚洲婷婷狠狠爱综合网| 1024视频免费在线观看| 97在线人人人人妻| 9色porny在线观看| 亚洲欧美一区二区三区黑人 | 中国美白少妇内射xxxbb| 午夜福利影视在线免费观看| 亚洲美女黄色视频免费看| 天天影视国产精品| 人人妻人人添人人爽欧美一区卜| 亚洲国产精品一区二区三区在线| 国产精品蜜桃在线观看| 国产精品三级大全| 国产男女超爽视频在线观看| 久久99热6这里只有精品| 爱豆传媒免费全集在线观看| 日本黄色日本黄色录像| 下体分泌物呈黄色| 日本wwww免费看| 国产精品不卡视频一区二区| 又黄又粗又硬又大视频| 男人操女人黄网站| h视频一区二区三区| 亚洲高清免费不卡视频| 亚洲精品日韩在线中文字幕| 亚洲少妇的诱惑av| 黄片播放在线免费| 成人毛片60女人毛片免费| 国产精品蜜桃在线观看| 免费av中文字幕在线| 欧美日韩综合久久久久久| 欧美3d第一页| 九色成人免费人妻av| 亚洲成色77777| 天堂俺去俺来也www色官网| 久久久久视频综合| 精品酒店卫生间| 久久人人爽人人片av| 成人国产av品久久久| 韩国高清视频一区二区三区| xxx大片免费视频| 性高湖久久久久久久久免费观看| 亚洲精品成人av观看孕妇| 一本色道久久久久久精品综合| 亚洲精品美女久久av网站| 丝袜脚勾引网站| 精品熟女少妇av免费看| 婷婷色av中文字幕| 一级片'在线观看视频| 亚洲经典国产精华液单| 午夜老司机福利剧场| 日本av免费视频播放| 在线观看免费视频网站a站| 国产精品女同一区二区软件| 九九在线视频观看精品| 国产黄色视频一区二区在线观看| 亚洲熟女精品中文字幕| 久久久精品免费免费高清| 午夜福利影视在线免费观看| 九九在线视频观看精品| 女人被躁到高潮嗷嗷叫费观| 91久久精品国产一区二区三区| 捣出白浆h1v1| 免费黄色在线免费观看| 国产精品人妻久久久影院| 少妇的逼水好多|