• <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
    美女视频免费永久观看网站| 亚洲欧美成人综合另类久久久| 最近最新免费中文字幕在线| 精品久久久久久久毛片微露脸 | 欧美黑人欧美精品刺激| 欧美精品高潮呻吟av久久| 99久久人妻综合| 999精品在线视频| 欧美成人午夜精品| 国产精品久久久久久精品古装| 丰满饥渴人妻一区二区三| 久久精品国产综合久久久| 日日夜夜操网爽| 91精品国产国语对白视频| 久久久久久久精品精品| 咕卡用的链子| 亚洲专区中文字幕在线| 成年人免费黄色播放视频| 嫁个100分男人电影在线观看| 国产激情久久老熟女| 亚洲精品国产av成人精品| 嫩草影视91久久| 啦啦啦 在线观看视频| 超色免费av| av一本久久久久| 下体分泌物呈黄色| 久久精品国产亚洲av香蕉五月 | 岛国毛片在线播放| 午夜精品久久久久久毛片777| 久久热在线av| 亚洲欧美激情在线| 中文字幕人妻丝袜一区二区| 99国产精品一区二区三区| 丰满饥渴人妻一区二区三| 91麻豆精品激情在线观看国产 | 亚洲精品在线美女| 国产精品二区激情视频| 两性夫妻黄色片| av天堂久久9| 侵犯人妻中文字幕一二三四区| 日韩 欧美 亚洲 中文字幕| 叶爱在线成人免费视频播放| 一本—道久久a久久精品蜜桃钙片| 啦啦啦在线免费观看视频4| 亚洲精品中文字幕在线视频| 国产三级黄色录像| 日本vs欧美在线观看视频| 午夜福利视频精品| 新久久久久国产一级毛片| 免费高清在线观看日韩| 精品久久蜜臀av无| 欧美精品啪啪一区二区三区 | 男女边摸边吃奶| 成人国产av品久久久| 777久久人妻少妇嫩草av网站| 巨乳人妻的诱惑在线观看| 国产有黄有色有爽视频| 久久久国产一区二区| 91九色精品人成在线观看| av线在线观看网站| 亚洲第一青青草原| 日韩欧美免费精品| 成年人黄色毛片网站| 国产麻豆69| 18禁黄网站禁片午夜丰满| 国产精品久久久久久人妻精品电影 | 日本精品一区二区三区蜜桃| 一区福利在线观看| 美女视频免费永久观看网站| 蜜桃国产av成人99| 狂野欧美激情性bbbbbb| 人妻人人澡人人爽人人| 我的亚洲天堂| 天堂8中文在线网| 久久久久久久精品精品| 亚洲欧美日韩高清在线视频 | 狠狠婷婷综合久久久久久88av| 在线观看免费视频网站a站| 国产老妇伦熟女老妇高清| 国产精品一区二区免费欧美 | 国产三级黄色录像| a级毛片在线看网站| 精品国产乱码久久久久久小说| 亚洲avbb在线观看| 日日爽夜夜爽网站| 高清欧美精品videossex| 韩国高清视频一区二区三区| 丝袜喷水一区| 日本精品一区二区三区蜜桃| 男女下面插进去视频免费观看| 精品熟女少妇八av免费久了| 美女扒开内裤让男人捅视频| 亚洲av成人一区二区三| 亚洲熟女毛片儿| 波多野结衣av一区二区av| 国产一区二区三区综合在线观看| 老鸭窝网址在线观看| 国产成人免费观看mmmm| 窝窝影院91人妻| 久久人人爽人人片av| 国产主播在线观看一区二区| 久久人妻熟女aⅴ| 69精品国产乱码久久久| 亚洲国产欧美在线一区| 中文字幕av电影在线播放| 日韩制服骚丝袜av| 少妇被粗大的猛进出69影院| 十八禁高潮呻吟视频| www.自偷自拍.com| 成人18禁在线播放| 亚洲欧美一区二区三区黑人| 免费高清视频大片| 亚洲性夜色夜夜综合| 狠狠狠狠99中文字幕| 欧美日韩乱码在线| 搡老妇女老女人老熟妇| 国产欧美日韩一区二区精品| 亚洲免费av在线视频| 亚洲成av人片免费观看| 一本一本综合久久| 一级作爱视频免费观看| 国产亚洲精品第一综合不卡| 精品人妻1区二区| 美女黄网站色视频| 午夜免费激情av| 级片在线观看| 又紧又爽又黄一区二区| 美女大奶头视频| 午夜亚洲福利在线播放| 91麻豆av在线| 亚洲欧美日韩东京热| 日韩中文字幕欧美一区二区| 他把我摸到了高潮在线观看| 精品免费久久久久久久清纯| 日本一区二区免费在线视频| 制服人妻中文乱码| 狂野欧美白嫩少妇大欣赏| 成人亚洲精品av一区二区| 欧美高清成人免费视频www| 久久草成人影院| 精品久久久久久成人av| 十八禁人妻一区二区| 青草久久国产| 色综合欧美亚洲国产小说| 中文字幕av在线有码专区| 午夜两性在线视频| 搡老岳熟女国产| 亚洲成av人片在线播放无| 国产午夜精品论理片| 欧美日韩黄片免| 18禁观看日本| 18禁国产床啪视频网站| 午夜亚洲福利在线播放| 久久这里只有精品19| 免费一级毛片在线播放高清视频| 国产亚洲精品久久久久久毛片| 美女大奶头视频| 最近视频中文字幕2019在线8| 国产av又大| 91大片在线观看| 亚洲精品一卡2卡三卡4卡5卡| 天天躁狠狠躁夜夜躁狠狠躁| 90打野战视频偷拍视频| 在线观看日韩欧美| 国产成人av激情在线播放| 久热爱精品视频在线9| 嫩草影视91久久| 国产主播在线观看一区二区| 18禁黄网站禁片午夜丰满| 欧美绝顶高潮抽搐喷水| 国产高清视频在线观看网站| 99久久精品国产亚洲精品| 国产乱人伦免费视频| 男女床上黄色一级片免费看| 婷婷丁香在线五月| 国产一区二区三区视频了| 国产一区二区三区在线臀色熟女| 一本精品99久久精品77| 精品久久久久久,| 午夜福利在线观看吧| 国产精品电影一区二区三区| 一a级毛片在线观看| 国产亚洲av嫩草精品影院| netflix在线观看网站| 亚洲美女视频黄频| 在线看三级毛片| 国产精品免费视频内射| 琪琪午夜伦伦电影理论片6080| 国产成人欧美在线观看| 悠悠久久av| 亚洲av第一区精品v没综合| 亚洲第一欧美日韩一区二区三区| 中文资源天堂在线| 日韩 欧美 亚洲 中文字幕| 亚洲成av人片免费观看| 亚洲av五月六月丁香网| 变态另类丝袜制服| 亚洲欧洲精品一区二区精品久久久| 亚洲午夜精品一区,二区,三区| 精品少妇一区二区三区视频日本电影| 性色av乱码一区二区三区2| 免费在线观看完整版高清| 欧美在线一区亚洲| 日韩欧美国产在线观看| 久热爱精品视频在线9| 一进一出好大好爽视频| 给我免费播放毛片高清在线观看| 国产aⅴ精品一区二区三区波| 国产精品亚洲美女久久久| 啦啦啦免费观看视频1| 老鸭窝网址在线观看| 欧美成人免费av一区二区三区| 国产精品亚洲美女久久久| 亚洲精品久久成人aⅴ小说| 亚洲精品久久国产高清桃花| 夜夜躁狠狠躁天天躁| 欧美一区二区精品小视频在线| 国产一区二区激情短视频| 日本一区二区免费在线视频| 中文字幕高清在线视频| 国产av在哪里看| 香蕉av资源在线| 青草久久国产| 亚洲精品av麻豆狂野| 波多野结衣高清无吗| 日本免费a在线| 亚洲第一欧美日韩一区二区三区| 婷婷亚洲欧美| 国产精品,欧美在线| 亚洲国产欧美一区二区综合| www国产在线视频色| 麻豆一二三区av精品| 国产一区二区三区视频了| 久久久久久久久久黄片| 香蕉国产在线看| 日日干狠狠操夜夜爽| 亚洲色图 男人天堂 中文字幕| 欧美一级a爱片免费观看看 | 国产av麻豆久久久久久久| 老司机靠b影院| 日韩欧美三级三区| 国产成人影院久久av| 久久人妻av系列| 18禁观看日本| 国产高清有码在线观看视频 | 国产午夜福利久久久久久| 久久中文看片网| 久久久国产精品麻豆| 视频区欧美日本亚洲| www.999成人在线观看| 美女 人体艺术 gogo| 国产成人av激情在线播放| 亚洲一区二区三区色噜噜| 日本精品一区二区三区蜜桃| 国产成人一区二区三区免费视频网站| 日日摸夜夜添夜夜添小说| 国产视频一区二区在线看| 777久久人妻少妇嫩草av网站| 一a级毛片在线观看| 精品欧美国产一区二区三| 超碰成人久久| 国产亚洲精品久久久久久毛片| 亚洲全国av大片| 精品久久久久久久毛片微露脸| 成熟少妇高潮喷水视频| 亚洲av片天天在线观看| 午夜a级毛片| 欧美成狂野欧美在线观看| 色综合婷婷激情| 亚洲成人中文字幕在线播放| 这个男人来自地球电影免费观看| 窝窝影院91人妻| 国产精品影院久久| 亚洲精品粉嫩美女一区| 丝袜人妻中文字幕| 国产片内射在线| 精品少妇一区二区三区视频日本电影| 亚洲av美国av| 禁无遮挡网站| 可以在线观看毛片的网站| 亚洲专区中文字幕在线| 亚洲熟妇熟女久久| 免费在线观看成人毛片| 又黄又粗又硬又大视频| 久久亚洲精品不卡| 91国产中文字幕| 一级毛片精品| 亚洲美女黄片视频| 男女视频在线观看网站免费 | 一进一出抽搐gif免费好疼| av在线天堂中文字幕| 色综合欧美亚洲国产小说| 国产精品免费视频内射| 日韩欧美三级三区| 亚洲国产精品成人综合色| 成人亚洲精品av一区二区| 国内少妇人妻偷人精品xxx网站 | 亚洲欧美日韩高清在线视频| 女人高潮潮喷娇喘18禁视频| 日韩大码丰满熟妇| 国产精品1区2区在线观看.| 亚洲熟妇熟女久久| 精品午夜福利视频在线观看一区| 日本五十路高清| 亚洲色图 男人天堂 中文字幕| 国产精品免费视频内射| 亚洲国产欧美一区二区综合| 亚洲av五月六月丁香网| 成人午夜高清在线视频| 在线十欧美十亚洲十日本专区| 久久精品综合一区二区三区| 黄色毛片三级朝国网站| 亚洲精品美女久久av网站| 中文在线观看免费www的网站 | 日本一本二区三区精品| 黄片小视频在线播放| 欧美中文综合在线视频| 在线观看一区二区三区| 国产97色在线日韩免费| 久久人人精品亚洲av| 国产麻豆成人av免费视频| 亚洲中文字幕一区二区三区有码在线看 | 久久久国产精品麻豆| 国内精品久久久久精免费| 他把我摸到了高潮在线观看| 伦理电影免费视频| 亚洲激情在线av| 美女扒开内裤让男人捅视频| 一区福利在线观看| 欧洲精品卡2卡3卡4卡5卡区| 亚洲美女视频黄频| 色尼玛亚洲综合影院| 嫩草影视91久久| 夜夜爽天天搞| 国内精品久久久久精免费| 久久久久久免费高清国产稀缺| 五月伊人婷婷丁香| 亚洲精华国产精华精| 看片在线看免费视频| 最好的美女福利视频网| 日韩欧美免费精品| 国产aⅴ精品一区二区三区波| 亚洲专区中文字幕在线| 舔av片在线| 欧美日韩亚洲国产一区二区在线观看| 一级黄色大片毛片| 国产精品av久久久久免费| svipshipincom国产片| 国产三级中文精品| 亚洲黑人精品在线| 高潮久久久久久久久久久不卡| 国产av在哪里看| 亚洲成av人片免费观看| 亚洲av五月六月丁香网| 亚洲中文av在线| 色av中文字幕| 不卡一级毛片| 最近最新中文字幕大全免费视频| 真人做人爱边吃奶动态| 18禁黄网站禁片免费观看直播| 国产精华一区二区三区| 一级毛片高清免费大全| 午夜免费观看网址| 一进一出抽搐动态| 免费无遮挡裸体视频| 国产精品九九99| √禁漫天堂资源中文www| 国内少妇人妻偷人精品xxx网站 | 亚洲aⅴ乱码一区二区在线播放 | 欧美日韩中文字幕国产精品一区二区三区| www国产在线视频色| 麻豆久久精品国产亚洲av| 一边摸一边抽搐一进一小说| 少妇熟女aⅴ在线视频| 一个人免费在线观看的高清视频| 天堂影院成人在线观看| 国产精品98久久久久久宅男小说| 成年版毛片免费区| 琪琪午夜伦伦电影理论片6080| 九色国产91popny在线| 一进一出抽搐gif免费好疼| 久久性视频一级片| 激情在线观看视频在线高清| 亚洲精品国产精品久久久不卡| 婷婷丁香在线五月| 亚洲专区中文字幕在线| 岛国在线免费视频观看| xxxwww97欧美| 麻豆成人av在线观看| 舔av片在线| 91大片在线观看| 男女午夜视频在线观看| 真人一进一出gif抽搐免费| 中文字幕av在线有码专区| 亚洲精品久久成人aⅴ小说| 麻豆国产av国片精品| 国产亚洲精品一区二区www| 亚洲中文av在线| www.熟女人妻精品国产| 一级黄色大片毛片| 丝袜美腿诱惑在线| 日本一本二区三区精品| 欧美又色又爽又黄视频| 国产真实乱freesex| 亚洲电影在线观看av| 99国产综合亚洲精品| 亚洲第一电影网av| 国产亚洲精品第一综合不卡| 国产av又大| 成人特级黄色片久久久久久久| 午夜福利免费观看在线| 亚洲第一欧美日韩一区二区三区| 老司机在亚洲福利影院| 一个人免费在线观看的高清视频| 亚洲国产中文字幕在线视频| 男人舔奶头视频| 国产成人影院久久av| 精品久久久久久久人妻蜜臀av| 美女黄网站色视频| 夜夜夜夜夜久久久久| 嫩草影视91久久| 女警被强在线播放| 动漫黄色视频在线观看| 国产成人精品无人区| 久久精品亚洲精品国产色婷小说| 亚洲av美国av| www.精华液| 午夜福利在线在线| 很黄的视频免费| 伊人久久大香线蕉亚洲五| 禁无遮挡网站| 欧美日本视频| 黄色女人牲交| aaaaa片日本免费| 亚洲色图av天堂| 亚洲人成77777在线视频| 男男h啪啪无遮挡| 亚洲成人国产一区在线观看| 日韩欧美在线二视频| 亚洲成人中文字幕在线播放| bbb黄色大片| 不卡一级毛片| 国产精华一区二区三区| 成年版毛片免费区| 舔av片在线| 亚洲av成人一区二区三| 视频区欧美日本亚洲| 久久久久久久久中文| 久久久水蜜桃国产精品网| 日韩欧美在线二视频| 美女午夜性视频免费| 麻豆成人午夜福利视频| 亚洲中文字幕一区二区三区有码在线看 | 久久精品91蜜桃| tocl精华| 精品国内亚洲2022精品成人| www日本黄色视频网| 久久久久久亚洲精品国产蜜桃av| 午夜福利高清视频| 国产精品精品国产色婷婷| 丰满人妻熟妇乱又伦精品不卡| 国产又色又爽无遮挡免费看| 精品久久蜜臀av无| 国产精品综合久久久久久久免费| 国产伦人伦偷精品视频| 人人妻人人看人人澡| 深夜精品福利| 日韩大尺度精品在线看网址| 九色成人免费人妻av| 天堂av国产一区二区熟女人妻 | 午夜视频精品福利| 亚洲国产精品合色在线| 国产三级黄色录像| 亚洲第一电影网av| 亚洲成人中文字幕在线播放| 别揉我奶头~嗯~啊~动态视频| 每晚都被弄得嗷嗷叫到高潮| 国产av又大| 三级男女做爰猛烈吃奶摸视频| 成人午夜高清在线视频| 欧美久久黑人一区二区| 午夜福利欧美成人| 两人在一起打扑克的视频| av天堂在线播放| 欧美色欧美亚洲另类二区| 亚洲最大成人中文| 久久精品综合一区二区三区| 国产精品野战在线观看| 一区福利在线观看| 制服诱惑二区| 亚洲人成网站在线播放欧美日韩| 久久精品91无色码中文字幕| 在线观看舔阴道视频| 久久久国产精品麻豆| 精品久久久久久久人妻蜜臀av| 国产精品久久久久久亚洲av鲁大| 麻豆国产av国片精品| 欧美精品亚洲一区二区| 久久久国产精品麻豆| 在线国产一区二区在线| a级毛片a级免费在线| 91在线观看av| 国产精品自产拍在线观看55亚洲| 一级毛片精品| 宅男免费午夜| 国产高清videossex| 日韩精品免费视频一区二区三区| 亚洲成人久久爱视频| 久久婷婷成人综合色麻豆| or卡值多少钱| 熟妇人妻久久中文字幕3abv| 香蕉丝袜av| 久久精品国产综合久久久| 男女床上黄色一级片免费看| 国内精品久久久久精免费| 亚洲激情在线av| e午夜精品久久久久久久| 亚洲一区中文字幕在线| 国产成人精品无人区| 精品一区二区三区av网在线观看| 国产蜜桃级精品一区二区三区| xxxwww97欧美| av中文乱码字幕在线| 可以在线观看的亚洲视频| 岛国在线观看网站| 欧美成人性av电影在线观看| 欧美激情久久久久久爽电影| 一夜夜www| 国产熟女xx| 国产伦人伦偷精品视频| 午夜日韩欧美国产| 亚洲中文av在线| 成人欧美大片| 国产精品免费一区二区三区在线| 欧美日韩亚洲国产一区二区在线观看| 男人舔女人的私密视频| 久久久久性生活片| 免费在线观看影片大全网站| aaaaa片日本免费| 又紧又爽又黄一区二区| 亚洲自拍偷在线| 99热这里只有精品一区 | 怎么达到女性高潮| 中文字幕人成人乱码亚洲影| www国产在线视频色| 日本免费一区二区三区高清不卡| 1024视频免费在线观看| 国产午夜福利久久久久久| 搡老岳熟女国产| 国产精品久久久久久久电影 | 久久久久久久久中文| 丁香六月欧美| 国产男靠女视频免费网站| 久久中文看片网| 黄频高清免费视频| 亚洲午夜理论影院| 99精品欧美一区二区三区四区| 免费在线观看黄色视频的| 人人妻人人看人人澡| 不卡一级毛片| 久久中文看片网| 在线观看免费日韩欧美大片| 国产人伦9x9x在线观看| 人人妻,人人澡人人爽秒播| 亚洲一码二码三码区别大吗| 欧美日韩国产亚洲二区| x7x7x7水蜜桃| 在线观看免费日韩欧美大片| 欧美 亚洲 国产 日韩一| 欧美日韩中文字幕国产精品一区二区三区| 久久午夜综合久久蜜桃| 黄色女人牲交| 国产人伦9x9x在线观看| 精品乱码久久久久久99久播| 91老司机精品| 亚洲熟妇熟女久久| 免费观看精品视频网站| 国产精品,欧美在线| 欧美黄色淫秽网站| avwww免费| 女警被强在线播放| 欧美日韩精品网址| 叶爱在线成人免费视频播放| 久久人妻av系列| 午夜视频精品福利| 精品久久久久久久末码| 亚洲最大成人中文| 国产亚洲精品第一综合不卡| 亚洲人成电影免费在线| 精品国产亚洲在线| 啪啪无遮挡十八禁网站| 亚洲欧美日韩高清在线视频| 欧美日韩亚洲综合一区二区三区_| 欧美一级毛片孕妇| 成人三级黄色视频| 99久久国产精品久久久| 色哟哟哟哟哟哟| 99久久无色码亚洲精品果冻| 欧美日本视频| 日日干狠狠操夜夜爽| 国内精品一区二区在线观看| 999精品在线视频| 桃红色精品国产亚洲av| 麻豆久久精品国产亚洲av| 久久久久国内视频| 国产精品亚洲一级av第二区| 国产高清视频在线观看网站| 国产高清videossex| 天堂√8在线中文| 最新在线观看一区二区三区| 精品久久久久久,| 法律面前人人平等表现在哪些方面| 美女免费视频网站| 午夜福利成人在线免费观看| 午夜a级毛片|