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

    微重力雙向溫差作用下Czochralski法硅熔體中的熱毛細(xì)對(duì)流*

    2016-12-29 09:04:54魏程星云張全壯
    功能材料 2016年12期
    關(guān)鍵詞:毛細(xì)熔體對(duì)流

    魏程星云,彭 嵐,張全壯

    (重慶大學(xué) 動(dòng)力工程學(xué)院,低品位能源利用技術(shù)及系統(tǒng)教育部重點(diǎn)實(shí)驗(yàn)室,重慶 400044)

    微重力雙向溫差作用下Czochralski法硅熔體中的熱毛細(xì)對(duì)流*

    魏程星云,彭 嵐,張全壯

    (重慶大學(xué) 動(dòng)力工程學(xué)院,低品位能源利用技術(shù)及系統(tǒng)教育部重點(diǎn)實(shí)驗(yàn)室,重慶 400044)

    采用數(shù)值模擬的方法研究了微重力條件下Czochralski法生長硅晶體過程中熔體熱毛細(xì)對(duì)流的基本特征,探討了水平和垂直溫度梯度的耦合對(duì)熔體流動(dòng)的影響。熔體自由表面與外界輻射換熱,水平溫度梯度Marangoni(Ma)數(shù)選取(0~3 000),底部熱流Q選取(1.39×10-2~1.76×10-2)。結(jié)果表明,當(dāng)Q和Ma數(shù)均較小時(shí),流動(dòng)為穩(wěn)態(tài),液池內(nèi)產(chǎn)生3個(gè)流胞,熔體流動(dòng)由Q主導(dǎo),減小Q或增大Ma數(shù)可使流動(dòng)更穩(wěn)定。當(dāng)Ma數(shù)增大到一定值時(shí),流動(dòng)從穩(wěn)態(tài)轉(zhuǎn)變?yōu)榉欠€(wěn)態(tài),流動(dòng)的臨界Mac數(shù)隨Q的增大而顯著減小。流動(dòng)失穩(wěn)后,出現(xiàn)了新的流動(dòng)轉(zhuǎn)變方式,Ma數(shù)為影響表面波動(dòng)形式的關(guān)鍵因素,Q會(huì)改變熱流體波數(shù),是晶體附近的熱流體波產(chǎn)生的決定因素。隨著Ma數(shù)和Q的不斷增強(qiáng),自由表面最終形成彎曲條幅狀熱流體波。

    熱毛細(xì)對(duì)流;Cz結(jié)構(gòu);雙向溫度梯度;數(shù)值模擬;表面張力

    0 引 言

    硅晶體因其優(yōu)越的熱學(xué)、電學(xué)、光學(xué)等特性被廣泛應(yīng)用于光伏工程、電子技術(shù)、航空航天技術(shù)等領(lǐng)域。目前,Czochralski生長法(Cz法)是制備硅晶體的最主要的方式之一。在Cz法制備硅單晶的過程中,有兩種不同形式的表面張力驅(qū)動(dòng)硅熔體進(jìn)行流動(dòng),一種是水平溫度梯度形成的表面張力所驅(qū)動(dòng)的熱毛細(xì)對(duì)流;另一種是垂直溫度梯度所驅(qū)動(dòng)的Marangoni對(duì)流。這兩種流動(dòng)相互作用,對(duì)Cz法生成晶體的質(zhì)量有著很大的影響。近些年來,許多學(xué)者對(duì)液池內(nèi)的熱毛細(xì)對(duì)流進(jìn)行了大量的研究[1-5]。Kyung-Woo Yi等[6]采用數(shù)值模擬的方法,發(fā)現(xiàn)熱毛細(xì)對(duì)流會(huì)影響Czochralski結(jié)構(gòu)中的流體溫度和流動(dòng)結(jié)構(gòu)的不對(duì)稱性,并且認(rèn)為Marangoni-Bénard不穩(wěn)定性和Rayleigh-Bénard不穩(wěn)定性是在不同表面張力系數(shù)下引起流動(dòng)失穩(wěn)的主要原因。Nepomnyashchy[7]最早使用數(shù)值模擬方法及線性穩(wěn)定性理論對(duì)微重力條件下的矩形液池內(nèi)的Bénard-Marangoni對(duì)流進(jìn)行了研究,闡述了垂直和水平的溫度梯度的相對(duì)大小對(duì)液層內(nèi)的流胞運(yùn)動(dòng)方向的影響。Jing等[8]發(fā)現(xiàn)Marangoni效應(yīng)對(duì)大尺寸Cz結(jié)構(gòu)內(nèi)熔體自由表面的波動(dòng)形式有著決定性的作用。Yao等[9-10]用數(shù)值模擬的方法對(duì)液橋內(nèi)的熱毛細(xì)對(duì)流進(jìn)行研究,結(jié)果表明磁場可以很好的抑制熱毛細(xì)對(duì)流,且普朗特?cái)?shù)的大小對(duì)熱毛細(xì)對(duì)流起著重要的作用。Peng等[11]對(duì)方形池內(nèi)的硅油進(jìn)行了不同溫差下的實(shí)驗(yàn)研究,實(shí)驗(yàn)結(jié)果表明浮力對(duì)流可以抑制熱毛細(xì)對(duì)流。Y.Takagi等[12]發(fā)現(xiàn)坩堝旋轉(zhuǎn)和磁場相結(jié)合對(duì)熱毛細(xì)對(duì)流產(chǎn)生的的熱流體波有很好的抑制作用。T.Yamamoto等[13]發(fā)現(xiàn)通過改變液池的寬深比可以控制熔體內(nèi)熱毛細(xì)對(duì)流的流動(dòng)方向。Wang等[14]利用三維數(shù)值模擬研究了環(huán)形淺液池內(nèi)垂直溫度梯度對(duì)硅油熱毛細(xì)對(duì)流的影響,結(jié)果表明增加底部熱流可以改變自由表面振蕩波波形,流動(dòng)最不穩(wěn)定的區(qū)域在環(huán)形池內(nèi)壁附近。

    以往國內(nèi)外的大部分研究主要集中于Marangoni對(duì)流或者熱毛細(xì)對(duì)流單獨(dú)存在下的情況,很少有研究報(bào)道流體中兩種流動(dòng)同時(shí)存在的情況,然而在晶體的實(shí)際生長過程中由于存在雙向溫差使得熔體中同時(shí)存在著兩種流動(dòng),稱為Marangoni-熱毛細(xì)對(duì)流;并且大部分學(xué)者研究的液池均為淺液池,而Cz法制備晶體的坩堝深度一般較深。因此,對(duì)雙向溫差下深液池中的熱毛細(xì)對(duì)流的研究顯得尤為重要。為了了解實(shí)際過程中雙向溫差耦合作用對(duì)熔體熱毛細(xì)對(duì)流的影響,本文采用有限差分法對(duì)Cz法生長硅晶體過程中的熱毛細(xì)對(duì)流進(jìn)行了三維數(shù)值模擬,得出了不同溫差和底部熱流密度耦合的情況下的硅熔體流動(dòng)的基本特征,重點(diǎn)分析了水平和豎直方向的溫度梯度耦合作用對(duì)熔體內(nèi)部流動(dòng)的影響。

    1 數(shù)學(xué)物理模型

    物理模型如圖1所示,坩堝內(nèi)半徑rc=50 mm,晶體半徑rs=15 mm,液池深度h=8 mm,坩堝底部以恒定均勻熱流q加熱,自由邊界為不變形的平面且與環(huán)境進(jìn)行輻射換熱,坩堝壁維持恒溫Th,結(jié)晶界面溫度為Tc=Tm,水平溫差ΔT=Th-Tc。

    圖1 物理模型

    為簡化計(jì)算,作出如下假定:(1) 熔體為不可壓縮牛頓流體;(2) 流速較低,流動(dòng)為層流;(3) 所有固-液界面滿足無滑移條件;(4) 液體自由表面平整無變形;(5) 液體自由表面考慮熱毛細(xì)力的作用并且與外界進(jìn)行輻射換熱。

    在上述假定條件下,引入無量綱參數(shù)溫度

    和無量綱熱流密度

    ·V=0

    (1)

    (2)

    (3)

    邊界條件為:

    自由表面

    (4)

    (5)

    熔體-晶體界面

    Z=H,0≤R≤Rs,0≤θ≤2π

    Vθ=VR=VZ=0,Θ=1

    (6)

    坩堝底部

    Z=0,0≤R<1,0≤θ≤2π

    (7)

    坩堝側(cè)壁

    0

    Vθ=VR=VZ=0,Θ=Θh

    (8)

    初始條件(τ=0)

    VR=VZ=Vθ=0,

    (9)

    其中

    分別為普朗特?cái)?shù),毛細(xì)雷諾數(shù)和熱輻射數(shù)。定義Marangoni數(shù)為

    式中

    為表面張力溫度系數(shù),v為動(dòng)量擴(kuò)散率,α為熱擴(kuò)散率,σ為斯忒藩-玻耳茲曼常數(shù),ε為發(fā)射率,λ為導(dǎo)熱系數(shù),ρ為熔體密度。熔點(diǎn)溫度Tm=1 683 K下硅熔體的物性參數(shù)如表1所示。

    表1Tm=1 683 K時(shí)硅熔體物性參數(shù)

    Table 1 Physical properties of the silicon melt atTm=1 683 K

    參數(shù),符號(hào)單位數(shù)值表面張力溫度系數(shù),γTN/(m·K)-7.0×10-5密度,ρkg/m32530動(dòng)力粘度,μkg/(m·s)7.0×10-4動(dòng)量擴(kuò)散率,vm2/s2.767×10-7比熱,CpKJ/(kg·K)1000導(dǎo)熱系數(shù),λW/(m·K)64熱擴(kuò)散率,αm2/s2.53×10-5輻射率,ε0.318

    利用有限容積法對(duì)基本方程進(jìn)行離散,擴(kuò)散項(xiàng)采用中心差分,對(duì)流項(xiàng)采用QUICK格式,壓力-速度修正采用SIMPLE方法,計(jì)算采用80r×40z×60θ的非均勻網(wǎng)格,時(shí)間步長取為1.19×10-3s,在每個(gè)時(shí)間步長內(nèi),如果控制容積上的最大殘差<10-8,則認(rèn)為計(jì)算收斂。計(jì)算程序的正確性已被大量的計(jì)算結(jié)果所證實(shí)[15],網(wǎng)格無關(guān)性驗(yàn)證如表2所示,因此計(jì)算中選擇的網(wǎng)格為80r×40z×60θ。

    表2 網(wǎng)格收斂性的檢驗(yàn)(Ma=1 581,Q=1.39×10-2)

    Table 2 Validation of grid convergence (Ma=1 581,Q=1.39×10-2)

    網(wǎng)格(r×z×θ)振蕩頻率波數(shù)80×32×400.084480×40×600.0833480×60×800.08354

    2 結(jié)果與討論

    2.1 基本流動(dòng)

    當(dāng)?shù)撞繜崃髅芏群蜏夭钶^小時(shí),液池內(nèi)硅熔體的流動(dòng)較弱為穩(wěn)態(tài)軸對(duì)稱流動(dòng),豎直平面內(nèi)的流動(dòng)流型可用等流函數(shù)線表示(圖中正、負(fù)號(hào)分別代表逆時(shí)針和順時(shí)針方向)。無因次流函數(shù)定義如下

    (10)

    2.1.1Ma數(shù)對(duì)基本流的影響

    當(dāng)水平溫差較小時(shí),流動(dòng)為穩(wěn)態(tài)。自由表面的熔體在熱毛細(xì)力的驅(qū)動(dòng)下從坩堝壁向晶體流動(dòng),在熔體中產(chǎn)生逆時(shí)針流胞。圖2給出了Q=1.53×10-2時(shí),不同Ma數(shù)下液池內(nèi)R-Z平面的等流函數(shù)線和等溫線分布圖。當(dāng)Ma數(shù)較小時(shí),垂直溫差主導(dǎo)熔體內(nèi)部流動(dòng),產(chǎn)生了占據(jù)主流區(qū)的順時(shí)針大流胞。晶體附近的熔體在水平溫差作用下形成一個(gè)逆時(shí)針小流胞,熔體內(nèi)形成雙胞流。隨著Ma數(shù)增大,左側(cè)小流胞流動(dòng)增強(qiáng),右側(cè)大流胞流動(dòng)被抑制,液池右上角出現(xiàn)微小的逆時(shí)針小流胞。熔體內(nèi)流函數(shù)最大值隨Ma數(shù)增加而減小。由等溫線圖可知,隨Ma數(shù)的增加,熔體最高溫度和自由表面最高溫度增加。晶體下方等溫線保持水平,熔體幾乎無流動(dòng)。

    圖2 R-Z平面等流函數(shù)線和等溫線圖

    圖3給出了Q=1.53×10-2時(shí),自由表面徑向速度分布圖和R=0.4處的徑向速度在Z方向的分布圖。

    圖3 自由表面徑向速度分布和R=0.4處的徑向速度沿Z向分布

    Fig 3 Distribution of radial velocity along the free surface andZ

    由圖3可見,自由表面上的徑向速度出現(xiàn)了3個(gè)極值??拷w和坩堝壁的逆時(shí)針流胞使得自由表面形成了正極大值,而右側(cè)順時(shí)針流胞使得自由表面形成了負(fù)極大值。隨著Ma數(shù)的增加,自由表面正向流動(dòng)增強(qiáng)而負(fù)向流動(dòng)減弱,說明熔體中逆時(shí)針流胞流動(dòng)增強(qiáng)而順時(shí)針流胞流動(dòng)被抑制。從R=0.4處的徑向速度分布圖可見,隨著Ma數(shù)的增大徑向速度也增大。自由表面張力誘導(dǎo)的流動(dòng)主要發(fā)生在0.1

    2.1.2 熱流密度Q對(duì)基本流的影響

    底部熱流加熱和自由表面的輻射散熱形成了貫穿液池的垂直方向的溫度梯度,誘發(fā)了自下而上順時(shí)針流動(dòng)的右側(cè)大流胞。圖4給出了Ma=198時(shí),液池內(nèi)R-Z平面的等流函數(shù)線和等溫線分布圖。隨著Q增大,熔體流動(dòng)強(qiáng)度增強(qiáng),右上角小流胞消失,坩堝內(nèi)剩下兩個(gè)反向旋轉(zhuǎn)的流胞。底部和自由表面的最高溫度逐漸坩堝壁移動(dòng)。熔體中心等溫線發(fā)生劇烈變化而晶體下部等溫線保持水平。

    當(dāng)Ma=198時(shí),垂直和水平溫差在靠近晶體側(cè)的驅(qū)動(dòng)力方向一致,此處的徑向速度達(dá)到了最大值,如圖5所示。隨著Q增加,左側(cè)正向速度和右側(cè)負(fù)向速度增加,右側(cè)正向速度減小。這說明增大Q促進(jìn)了左側(cè)流胞和右側(cè)大流胞的流動(dòng),抑制了右側(cè)小流胞的流動(dòng)。R=0.4處熔體流動(dòng)為逆時(shí)針流動(dòng),自由表面溫度梯度驅(qū)動(dòng)的徑向內(nèi)流集中在0.1

    2.2 流動(dòng)轉(zhuǎn)變

    當(dāng)?shù)撞繜崃髅芏群妥杂杀砻鏈囟忍荻炔粩嘣龃笾钡匠^一定值時(shí),流場中任何一個(gè)微小的擾動(dòng)都會(huì)不斷地放大,使得流動(dòng)狀態(tài)從穩(wěn)態(tài)變?yōu)榉欠€(wěn)態(tài)。采用逼近法確定了流動(dòng)轉(zhuǎn)變的臨界Mac數(shù)。如圖6所示,增大Q可以使得流動(dòng)轉(zhuǎn)變的臨界值減小,且較小的Q的改變可以造成較大臨界Mac數(shù)變化。

    圖4 R-Z平面等流函數(shù)線和等溫線圖

    圖5 自由表面徑向速度分布和R=0.4處的徑向速度沿Z向分布

    Fig 5 Distribution of radial velocity along the free surface andZ

    圖6 臨界Mac數(shù)隨Q的變化

    2.3 三維振蕩流動(dòng)

    當(dāng)Ma數(shù)超過臨界值后,熔體的流動(dòng)轉(zhuǎn)變?yōu)槿S非穩(wěn)態(tài)流動(dòng)。此時(shí)硅熔體的流動(dòng)較為復(fù)雜,為了更直觀的了解熔體的周向速度的時(shí)相關(guān)特性,研究水平和垂直溫度梯度對(duì)流動(dòng)的影響,定義自由表面溫度波動(dòng)δT和自由表面速度波動(dòng)δV為

    (11)

    (12)

    2.3.1Ma數(shù)對(duì)振蕩流的影響

    由圖7可見,隨著Ma數(shù)的增大,自由表面的最大溫度波動(dòng)值不斷增加。圖8為Q=1.53×10-2,Ma=692時(shí),一個(gè)周期內(nèi)流動(dòng)演變的過程。由圖8可見,自由表面上速度波動(dòng)呈波數(shù)為5的單或雙層花瓣?duì)钫袷幉?,速度波?dòng)集中于晶體附近。振蕩由外向內(nèi)傳遞,為單、雙層交替產(chǎn)生的結(jié)構(gòu)。圖9為R=0.35處的自由表面速度波動(dòng)隨時(shí)間變化的STD圖,圖中可見一簇明暗相間的直線且每條線從下往上由亮變暗,說明該點(diǎn)振蕩在徑向上有周期性的強(qiáng)弱變化,周向上沒有運(yùn)動(dòng)。

    圖7 不同Q下的自由表面溫度波動(dòng)最大值隨Ma數(shù)的變化

    Fig 7 The max surface azimuthal temperature changes along the Ma at differentQ

    圖8 Ma=642, Q=1.53×10-2時(shí),一個(gè)周期內(nèi)表面速度波動(dòng)隨時(shí)間的變化

    圖9Ma=642時(shí)R=0.35處表面速度波動(dòng)的STD圖

    Fig 9 The space-time diagram (STD) of surface azimuthal velocity along the circumference atR=0.35 forMa=642

    圖10給出了Q=1.53×10-2時(shí)不同Ma數(shù)下自由表面的速度波動(dòng)圖和STD圖。當(dāng)Ma=791時(shí),表面波動(dòng)轉(zhuǎn)變?yōu)榫w附近振蕩的波數(shù)為5的花瓣?duì)顭崃黧w波。隨著Ma數(shù)的進(jìn)一步增大,振蕩從晶體附近擴(kuò)散到了整個(gè)熔體表面。當(dāng)Ma=988時(shí),晶體附近熱流體波的流動(dòng)區(qū)域被壓縮,在熱流體波外層出現(xiàn)了波數(shù)為7的振蕩波,由STD圖可見,外部振蕩為三維穩(wěn)態(tài)結(jié)構(gòu)。隨著Ma數(shù)增加到2 174,表面波動(dòng)轉(zhuǎn)變?yōu)椴〝?shù)為5的熱流體波。從STD圖可以看出,此時(shí)熱流體波振幅呈周期性的變化。由此可見,自由表面振蕩隨Ma數(shù)的增加從晶體附近擴(kuò)散到整個(gè)表面,隨著Ma數(shù)增大,流動(dòng)越來越不穩(wěn)定。

    2.3.2 熱流密度Q對(duì)振蕩流的影響

    在液池底部所加的熱流密度和表面散熱產(chǎn)生了垂直方向的溫度梯度,由圖7可見,自由表面的最大溫度波動(dòng)值隨Q的增加不斷增加。

    圖11給出了Q=1.62×10-2和Ma=692時(shí),自由表面速度波動(dòng)和STD圖??梢钥闯?,隨著Q的增大,熱流體波波數(shù)和圖10(a)相比減少到4個(gè),波動(dòng)集中在晶體附近,STD圖由四對(duì)明暗相間的傾斜條紋組成。

    圖12給出了Q=1.62×10-2,Ma=1 383時(shí)的表面速度波動(dòng)和STD圖,在晶體附近存在波數(shù)為4的熱流體波。其流動(dòng)區(qū)域較圖11明顯減小,在熱流體波外側(cè)出現(xiàn)了7條直幅狀振蕩波。

    圖11Q=1.62×10-2,Ma=988,R=0.35表面速度波動(dòng)圖和STD圖

    Fig 11 Snapshots of surface azimuthal velocity and STD

    由R=0.35處的STD圖可以看出,晶體附近波動(dòng)為熱流體波,從R=0.75處的STD圖可知,外圍振蕩波為三維穩(wěn)態(tài)結(jié)構(gòu)。與圖10(b)對(duì)比可以發(fā)現(xiàn),當(dāng)流動(dòng)轉(zhuǎn)變?yōu)橥獠咳S穩(wěn)態(tài)而內(nèi)部為熱流體波的雙層結(jié)構(gòu)時(shí),隨著Q的增大,晶體附近的熱流體波流動(dòng)區(qū)域增大,流動(dòng)強(qiáng)度增強(qiáng)。由圖13可見,當(dāng)Q較小而Ma數(shù)較大時(shí),自由表面波動(dòng)主要為三維穩(wěn)態(tài),而熱流體波則被壓縮到晶體附近區(qū)域幾乎消失,當(dāng)Ma數(shù)更大時(shí),熱流體波完全消失。

    圖12 Q=1.62×10-2, Ma=1 383表面速度波動(dòng)圖和STD圖

    圖13Q=1.46×10-2,Ma=691時(shí)的表面速度波動(dòng)圖和STD圖

    Fig 13 Snapshots of surface azimuthal velocity fluctuation and STD

    3 結(jié) 論

    采用數(shù)值模擬的方法對(duì)雙向溫度梯度下Czochralski法生長硅晶體中的熔體熱毛細(xì)對(duì)流進(jìn)行了研究。分析了垂直和水平溫差對(duì)熔體流動(dòng)的影響,確定了流動(dòng)轉(zhuǎn)變的臨界值。結(jié)果表明:

    (1) 當(dāng)Ma數(shù)和Q較小時(shí),流動(dòng)為穩(wěn)態(tài)。隨Ma數(shù)的增加和Q的減小,熔體右上角出現(xiàn)一個(gè)微小的流胞,晶體下部熔體幾乎沒有流動(dòng)。此時(shí)的熔體流動(dòng)Q起主導(dǎo)作用,右上角小流胞和右側(cè)大流胞相互抑制,增大Ma數(shù)和減小Q可使得熔體整體流動(dòng)減弱,流動(dòng)更穩(wěn)定;

    (2) 隨著Ma數(shù)的增大,流動(dòng)轉(zhuǎn)變?yōu)榉欠€(wěn)態(tài)。臨界Mac數(shù)隨Q的增大而減小,且Q微小的改變能造成臨界Mac數(shù)的顯著改變。硅晶體生長過程中使用較小的Q可以延緩流動(dòng)失穩(wěn)現(xiàn)象。

    (3) 當(dāng)流動(dòng)失穩(wěn)轉(zhuǎn)變?yōu)榉欠€(wěn)態(tài)后,Ma數(shù)為改變自由表面波形的關(guān)鍵因素。隨著Ma數(shù)的增加,自由表面波形由單、雙層花瓣互相轉(zhuǎn)換的振蕩波轉(zhuǎn)變?yōu)榛ò隊(duì)顭崃黧w波,接著轉(zhuǎn)變?yōu)橥鈱又狈ê蛢?nèi)層熱流體波同時(shí)存在的結(jié)構(gòu)并最終轉(zhuǎn)變?yōu)閺澢臈l幅狀熱流體波,隨著Ma數(shù)的增大熔體自由表面波動(dòng)不斷增強(qiáng);

    (4) 流動(dòng)失穩(wěn)后,隨著Q的增大,熱流體波波數(shù)減小,熔體自由表面波動(dòng)不斷增強(qiáng)。當(dāng)Q相對(duì)于Ma數(shù)較小時(shí),靠近晶體側(cè)熱流體波消失,晶體附近流動(dòng)較穩(wěn)定;而增大Q,晶體附近速度波動(dòng)增加,選用較大的Ma數(shù)和較小的Q有利于增加晶體附近流體流動(dòng)的穩(wěn)定性。

    綜上所述,底部熱流Q對(duì)熔體內(nèi)Marangoni-熱毛細(xì)對(duì)流的穩(wěn)定性有著極為重要的影響,Q過大容易導(dǎo)致流動(dòng)失穩(wěn)并且會(huì)造成結(jié)晶界面附近劇烈振蕩,影響晶體質(zhì)量。而選擇較大的Ma數(shù)可以抑制結(jié)晶界面附近的振蕩。因此在晶體生長過程中可以選用較小的Q搭配較大的Ma數(shù)來減弱結(jié)晶界面周圍的流體振蕩,提高晶體質(zhì)量。

    [1] Kumar V, Biswas G, Brenner G, et al. Effect of thermocapillary convection in an industrial Czochralski crucible: numerical simulation [J]. International Journal of Heat and Mass Transfer, 2003, 46(9):1641-1652.

    [2] Anilkumar A V, Grugel R N, Bhowmick J, et al. Suppression of thermocapillary oscillations in sodium nitrate floating half-zones by high-frequency end-wall vibrations[J]. Journal of Crystal Growth, 2004, 276(s1-2):194-203.

    [3] Teitel M, Schwabe D, Gelfgat A Y. Experimental and computational study of flow instabilities in a model of Czochralski growth [J]. Journal of Crystal Growth, 2008, 310(7-9):1343-1348.

    [4] Hintz P, Schwabe D, Wilke H. Convection in a Czochralski crucible-part1: non-rotating crystal [J]. Journal of Crystal Growth, 2001, 222(1-2):343-355.

    [5] Schwabe D. Buoyant-thermocapillary and pure thermocapillary convective instabilities in Czochralski systems [J]. Journal of Crystal Growth, 2002, 237-239:1849-1853.

    [6] Yi K W, Kakimoto K, Eguchi M, et al. Spoke patterns on molten silicon in Czochralski system[J]. Journal of Crystal Growth, 1994, 144(1): 20-28.

    [7] Nepomnyashchy A A, Simanovskii I B, Braverman L M. Stability of thermocapillary flows with inclined temperature gradient[J]. Journal of Fluid Mechanics, 2001, 442: 141-155.

    [8] Jing C J, Tsukada T, Hozawa M, et al. Numerical studies of wave pattern in an oxide melt in the Czochralski crystal growth[J]. Journal of Crystal Growth, 2004, 265(265):505-517.

    [9] Yao L, Zeng Z, Li X, et al. Effects of rotating magnetic fields on thermocapillary flow in a floating half-zone[J]. Journal of Crystal Growth, 2011, 316(1):177-184.

    [10] Yao L, Zeng Z, Chen J, et al. Investigation of convection control under the non-uniform RMF in a liquid bridge[J]. Procedia Engineering, 2012, 31(4):659-664.

    [11] Peng Z, Li D, Qi K. Transition to chaos in thermocapillary convection[J]. International Journal of Heat & Mass Transfer, 2013, 57(2):457-464.

    [12] Takagi Y, Okano Y, Minakuchi H, et al. Combined effect of crucible rotation and magnetic field on hydrothermal wave[J]. Journal of Crystal Growth, 2014, 385(1):72-76.

    [13] Yamamoto T, Takagi Y, Okano Y. Numerical investigation of the effect of free surface shape on the direction of thermocapillary flow in a thin circular pool[J]. Journal of Chemical Engineering of Japan, 2015, 48(6):407-417.

    [14] Wang F, Peng L, Zhang Q Z, et al. Effect of the vertical heat transfer on the thermocapillary convection in a shallow annular pool[J]. Microgravity Science & Technology, 2015, 27(2):107-114.

    [15] Li Y R, Imaishi N, Azami T, et al. Three-dimensional oscillatory flow in a thin annular pool of silicon melt[J]. Journal of Crystal Growth, 2004, 260(1):28-42.

    Thermocapillary convection in Czochralski silicon growth by bidirectional temperature gradients under microgravity

    WEI Chengxingyun,PENG Lan, ZHANG Quanzhuang

    (Key Laboratory of Low-Energy Utilization Technologies and Systems of Ministry of Education, College of Power Engineering, Chongqing University, Chongqing 400044,China)

    The bidirectional temperature gradients play an important part in the convection of crystal growth. However, most of the researches before only focused on unidirectional temperature. In this paper, under microgravity, a serious of numerical simulations was conducted to study the thermocapillary convection in Czochralski growth of silicon crystal. The coupling effects of vertical and horizontal temperature gradients were discussed. The radiation on the melt surface was considered, the range of horizontal and vertical temperature gradient wereMa=(0-3000) andQ=(1.39×10-2-1.76×10-2) respectively. Results show that flow is steady for smallQandMa. In this case, there are three flow cells in the melt, the melt flow is dominated byQ; the flow is more stable due to the decrease of theQand the increase of theMa. When the Ma exceeds a threshold value, an unsteady three dimensional flow is developed. With the increase of theQ, the critical Ma significantly decreased. When the flow is unsteady, a new way of the flow pattern transtion is found.Mais the key for the change of the surface fluctuation pattern. The number of the hydrothermal wave is changed byQ, and the generation of hydrothermal wave which is near the crystal is determined byQ. With the enhancement of theQandMa, the flow pattern enventually changes into curved spoke patterns.

    thermocapillary convection; Cz configuation; bidirectional temperature gradients; numerical simulation; surface tension

    1001-9731(2016)12-12090-07

    國家自然科學(xué)基金資助項(xiàng)目(51276203)

    2016-01-04

    2016-03-15 通訊作者:彭 嵐,E-mail:penglan@cqu.edu.cn

    魏程星云 (1991-),男,武漢人,在讀碩士,師承彭嵐教授,從事工程熱物理、熱毛細(xì)對(duì)流相關(guān)研究。

    TK124

    A

    10.3969/j.issn.1001-9731.2016.12.014

    猜你喜歡
    毛細(xì)熔體對(duì)流
    “毛細(xì)”政務(wù)號(hào)關(guān)停背后
    廉政瞭望(2024年5期)2024-05-26 13:21:07
    齊口裂腹魚集群行為對(duì)流態(tài)的響應(yīng)
    多孔建筑材料毛細(xì)吸水系數(shù)連續(xù)測量的方法
    能源工程(2020年6期)2021-01-26 00:55:16
    出現(xiàn)憋喘 可能是毛細(xì)支氣管炎!
    媽媽寶寶(2017年3期)2017-02-21 01:22:16
    聚合物熔體脈振傳遞過程的協(xié)同學(xué)研究
    中國塑料(2016年4期)2016-06-27 06:33:48
    注射保壓過程中O2/N2分子在PMMA熔體內(nèi)部的擴(kuò)散行為
    中國塑料(2016年3期)2016-06-15 20:30:01
    基于ANSYS的自然對(duì)流換熱系數(shù)計(jì)算方法研究
    高滲鹽水霧化吸入治療毛細(xì)支氣管炎的療效觀察
    含硅芳炔樹脂及其共混物熔體的流變性能
    二元驅(qū)油水界面Marangoni對(duì)流啟動(dòng)殘余油機(jī)理
    国产男人的电影天堂91| 国产白丝娇喘喷水9色精品| 婷婷色av中文字幕| 国产精品久久视频播放| 国产精品国产高清国产av| 国产亚洲欧美98| 亚洲aⅴ乱码一区二区在线播放| 久久亚洲精品不卡| 亚洲成人中文字幕在线播放| 国产精品无大码| 大香蕉久久网| 久久久久久久久中文| 国产精品电影一区二区三区| 国产在视频线在精品| 国产精品一二三区在线看| 一级黄片播放器| 狂野欧美白嫩少妇大欣赏| 99久久精品国产国产毛片| 国产精品久久视频播放| 偷拍熟女少妇极品色| 丰满乱子伦码专区| h日本视频在线播放| 精品久久久久久久久久免费视频| 一级av片app| 国产中年淑女户外野战色| 热99re8久久精品国产| 欧美色视频一区免费| 久久久久网色| 成人性生交大片免费视频hd| 18禁在线无遮挡免费观看视频| 最近中文字幕高清免费大全6| 91久久精品国产一区二区成人| 亚洲av第一区精品v没综合| 亚洲国产欧美人成| 亚洲成人中文字幕在线播放| 看非洲黑人一级黄片| 少妇被粗大猛烈的视频| av天堂在线播放| 亚洲一区二区三区色噜噜| 色哟哟哟哟哟哟| av在线蜜桃| 久久热精品热| 久久热精品热| 神马国产精品三级电影在线观看| 两个人的视频大全免费| 欧美精品一区二区大全| 国产精品综合久久久久久久免费| 成人毛片a级毛片在线播放| 午夜福利高清视频| 成人毛片a级毛片在线播放| 免费观看在线日韩| 国语自产精品视频在线第100页| 大香蕉久久网| 国产精品一区二区三区四区免费观看| 不卡视频在线观看欧美| 人人妻人人澡人人爽人人夜夜 | 男女下面进入的视频免费午夜| 日本-黄色视频高清免费观看| 国产探花极品一区二区| 欧美又色又爽又黄视频| 丝袜喷水一区| 大又大粗又爽又黄少妇毛片口| 成人一区二区视频在线观看| 中文字幕av在线有码专区| 国产伦精品一区二区三区四那| 日韩一区二区视频免费看| 日韩大尺度精品在线看网址| 在线国产一区二区在线| 在线播放无遮挡| 熟妇人妻久久中文字幕3abv| 欧美日本视频| 亚洲国产日韩欧美精品在线观看| 波多野结衣巨乳人妻| 赤兔流量卡办理| 97超碰精品成人国产| 亚洲性久久影院| 国产精品99久久久久久久久| 日本熟妇午夜| 国产高清不卡午夜福利| 欧美三级亚洲精品| 国产一级毛片七仙女欲春2| 久久精品91蜜桃| 国产精品综合久久久久久久免费| 波多野结衣高清无吗| 男人舔奶头视频| 村上凉子中文字幕在线| 欧美日本亚洲视频在线播放| 精品国产三级普通话版| 丰满人妻一区二区三区视频av| 国产精品,欧美在线| 亚洲av成人av| 国产亚洲91精品色在线| 嫩草影院精品99| 九九热线精品视视频播放| 日韩欧美 国产精品| 亚洲第一区二区三区不卡| 男的添女的下面高潮视频| 一区二区三区高清视频在线| 在线观看av片永久免费下载| 男人舔奶头视频| 高清日韩中文字幕在线| 神马国产精品三级电影在线观看| 亚洲成人久久性| 成人特级黄色片久久久久久久| 男女啪啪激烈高潮av片| 国产爱豆传媒在线观看| 亚洲中文字幕一区二区三区有码在线看| 国产高潮美女av| 欧美+亚洲+日韩+国产| 久久人人精品亚洲av| 亚洲av中文av极速乱| 亚洲av中文av极速乱| 亚洲自拍偷在线| 国产精品久久久久久亚洲av鲁大| 一级黄色大片毛片| 亚洲中文字幕一区二区三区有码在线看| 男女啪啪激烈高潮av片| 久久久午夜欧美精品| 久久欧美精品欧美久久欧美| 免费看a级黄色片| 久久99精品国语久久久| 国产 一区精品| av国产免费在线观看| 国产精品久久久久久精品电影| 婷婷六月久久综合丁香| 国内揄拍国产精品人妻在线| 精品人妻视频免费看| 日本与韩国留学比较| 欧美精品一区二区大全| 久久欧美精品欧美久久欧美| 亚洲欧美日韩无卡精品| 久久99热这里只有精品18| 爱豆传媒免费全集在线观看| 国产成人精品婷婷| 超碰av人人做人人爽久久| 性插视频无遮挡在线免费观看| 精品一区二区三区视频在线| 又爽又黄无遮挡网站| 尤物成人国产欧美一区二区三区| 啦啦啦韩国在线观看视频| 赤兔流量卡办理| 少妇高潮的动态图| 免费观看在线日韩| ponron亚洲| 国产 一区精品| 亚洲综合色惰| 亚洲国产欧美在线一区| a级毛片a级免费在线| 舔av片在线| 搡老妇女老女人老熟妇| 国产高清有码在线观看视频| 精品一区二区三区视频在线| 国产淫片久久久久久久久| 久久精品国产亚洲网站| 欧美成人免费av一区二区三区| 午夜久久久久精精品| 天堂av国产一区二区熟女人妻| 国产精品.久久久| 亚洲国产高清在线一区二区三| 精品熟女少妇av免费看| 国产极品精品免费视频能看的| 亚洲一区二区三区色噜噜| 此物有八面人人有两片| 日本五十路高清| 欧美日韩在线观看h| 亚洲精品自拍成人| 18禁黄网站禁片免费观看直播| 美女国产视频在线观看| 波多野结衣高清作品| 国产黄色小视频在线观看| 一边摸一边抽搐一进一小说| 国产美女午夜福利| 婷婷色av中文字幕| 尤物成人国产欧美一区二区三区| 中文字幕av在线有码专区| 久久久久久久久久久免费av| 寂寞人妻少妇视频99o| 日本-黄色视频高清免费观看| 亚洲精品日韩在线中文字幕 | 国国产精品蜜臀av免费| 欧美性猛交╳xxx乱大交人| 一级毛片电影观看 | 又粗又爽又猛毛片免费看| 久久婷婷人人爽人人干人人爱| 精品午夜福利在线看| 亚洲国产欧美在线一区| 一边亲一边摸免费视频| 狠狠狠狠99中文字幕| 亚洲欧美日韩高清专用| 久久人人爽人人片av| 亚洲人成网站在线播| 精品国内亚洲2022精品成人| 少妇熟女aⅴ在线视频| 国产精品一区二区性色av| 成人三级黄色视频| 我的女老师完整版在线观看| 欧美一区二区国产精品久久精品| 欧美激情在线99| 99在线视频只有这里精品首页| 亚洲精品456在线播放app| 亚洲精华国产精华液的使用体验 | 中文字幕熟女人妻在线| 嫩草影院入口| 在线免费观看不下载黄p国产| 国产色爽女视频免费观看| 国内精品久久久久精免费| 中文资源天堂在线| 边亲边吃奶的免费视频| 看非洲黑人一级黄片| 国国产精品蜜臀av免费| videossex国产| 免费看日本二区| 成人美女网站在线观看视频| 亚洲,欧美,日韩| 在线观看av片永久免费下载| 国产伦精品一区二区三区视频9| 久久久精品94久久精品| 一级毛片我不卡| 久久久久性生活片| 夜夜夜夜夜久久久久| 亚洲av成人精品一区久久| 网址你懂的国产日韩在线| 自拍偷自拍亚洲精品老妇| 亚洲内射少妇av| 18禁裸乳无遮挡免费网站照片| 国产精品蜜桃在线观看 | 三级男女做爰猛烈吃奶摸视频| 久久久午夜欧美精品| 亚洲一级一片aⅴ在线观看| 久久久久久大精品| 国模一区二区三区四区视频| 韩国av在线不卡| 99久久成人亚洲精品观看| 九草在线视频观看| 久久欧美精品欧美久久欧美| 九色成人免费人妻av| 韩国av在线不卡| 人妻制服诱惑在线中文字幕| 久久99热这里只有精品18| 国产精华一区二区三区| 校园春色视频在线观看| 亚洲国产日韩欧美精品在线观看| 天堂中文最新版在线下载 | 老熟妇乱子伦视频在线观看| 亚洲精品久久久久久婷婷小说 | 在线观看av片永久免费下载| 中文精品一卡2卡3卡4更新| 少妇猛男粗大的猛烈进出视频 | 一本精品99久久精品77| 亚洲色图av天堂| 日本与韩国留学比较| 丰满乱子伦码专区| 能在线免费观看的黄片| 天堂网av新在线| 在线观看免费视频日本深夜| 久久精品久久久久久噜噜老黄 | 国产精品女同一区二区软件| 国产黄色视频一区二区在线观看 | 久久久久久久久久久丰满| 黄色一级大片看看| 国产 一区精品| 超碰av人人做人人爽久久| 国产精品一二三区在线看| 人人妻人人看人人澡| 日本免费一区二区三区高清不卡| 丰满人妻一区二区三区视频av| 久久6这里有精品| 少妇熟女欧美另类| 中国美白少妇内射xxxbb| 2022亚洲国产成人精品| 国产一区二区在线av高清观看| 男人的好看免费观看在线视频| 亚洲图色成人| 亚洲av二区三区四区| 久久这里只有精品中国| 免费搜索国产男女视频| 免费看av在线观看网站| 在线播放无遮挡| 欧美一区二区精品小视频在线| 哪里可以看免费的av片| 99久久精品国产国产毛片| 高清在线视频一区二区三区 | 国产精品伦人一区二区| 国产午夜精品一二区理论片| 亚洲激情五月婷婷啪啪| 天堂√8在线中文| 国产精华一区二区三区| 国产精品麻豆人妻色哟哟久久 | 一进一出抽搐动态| 青春草视频在线免费观看| 国产熟女欧美一区二区| 男女做爰动态图高潮gif福利片| 欧美激情在线99| 欧洲精品卡2卡3卡4卡5卡区| 色综合色国产| 婷婷色综合大香蕉| 国产黄色小视频在线观看| 国产精品久久久久久久久免| 国产探花极品一区二区| 久久精品夜夜夜夜夜久久蜜豆| 97热精品久久久久久| 高清在线视频一区二区三区 | 国产91av在线免费观看| 老熟妇乱子伦视频在线观看| 国产蜜桃级精品一区二区三区| 99久久精品国产国产毛片| 亚洲三级黄色毛片| 日韩欧美 国产精品| 波多野结衣高清作品| av女优亚洲男人天堂| 中文在线观看免费www的网站| 欧美日韩综合久久久久久| 亚洲一级一片aⅴ在线观看| 国产精品久久久久久精品电影小说 | 最近最新中文字幕大全电影3| 成人欧美大片| 在线观看av片永久免费下载| a级毛片免费高清观看在线播放| 久99久视频精品免费| 午夜久久久久精精品| 寂寞人妻少妇视频99o| 欧美日韩精品成人综合77777| 亚洲国产精品成人综合色| 精品人妻一区二区三区麻豆| 婷婷精品国产亚洲av| 男女做爰动态图高潮gif福利片| 国产精品久久视频播放| 国产在线精品亚洲第一网站| 久久精品久久久久久久性| 久99久视频精品免费| 亚洲人成网站在线观看播放| 国产精品久久久久久亚洲av鲁大| 黄色视频,在线免费观看| 婷婷色av中文字幕| 色综合色国产| 亚洲综合色惰| 99在线视频只有这里精品首页| 精品久久久久久久人妻蜜臀av| 日产精品乱码卡一卡2卡三| 国产高清有码在线观看视频| 亚洲国产精品久久男人天堂| 欧美日韩一区二区视频在线观看视频在线 | 国产精品麻豆人妻色哟哟久久 | 久久九九热精品免费| 成人综合一区亚洲| 精品人妻偷拍中文字幕| 青春草亚洲视频在线观看| 午夜福利成人在线免费观看| 国产高清不卡午夜福利| 免费不卡的大黄色大毛片视频在线观看 | 直男gayav资源| 日本-黄色视频高清免费观看| 男女边吃奶边做爰视频| 国产中年淑女户外野战色| 国产成人freesex在线| 午夜a级毛片| 亚洲最大成人手机在线| 欧美高清成人免费视频www| av国产免费在线观看| 亚洲最大成人手机在线| 波多野结衣巨乳人妻| 少妇人妻一区二区三区视频| av卡一久久| 可以在线观看毛片的网站| 欧美日韩在线观看h| 99riav亚洲国产免费| 精品国内亚洲2022精品成人| 日本与韩国留学比较| 国产精品久久视频播放| 国产激情偷乱视频一区二区| 美女大奶头视频| av卡一久久| 国产色婷婷99| 国产爱豆传媒在线观看| av在线播放精品| 乱码一卡2卡4卡精品| 悠悠久久av| 日本av手机在线免费观看| 亚洲第一电影网av| 国产精品一区二区性色av| 欧美性感艳星| 亚洲av电影不卡..在线观看| 亚洲第一区二区三区不卡| 我要看日韩黄色一级片| 亚洲精品乱码久久久v下载方式| 91午夜精品亚洲一区二区三区| 日本爱情动作片www.在线观看| 亚洲美女搞黄在线观看| 国产极品精品免费视频能看的| 精品久久久久久久久av| 日本黄色片子视频| 在线a可以看的网站| 深夜精品福利| 少妇人妻精品综合一区二区 | 免费看a级黄色片| 日韩高清综合在线| 日韩欧美在线乱码| 婷婷色综合大香蕉| 99久久精品国产国产毛片| 在线播放无遮挡| 亚洲av男天堂| 欧美xxxx性猛交bbbb| 岛国毛片在线播放| 国产成人精品婷婷| 久久精品国产99精品国产亚洲性色| 美女脱内裤让男人舔精品视频 | 精品人妻熟女av久视频| www.av在线官网国产| 精品午夜福利在线看| 成人二区视频| 中文亚洲av片在线观看爽| 久久久精品94久久精品| 国产精品永久免费网站| 男女啪啪激烈高潮av片| 国产精品久久久久久av不卡| av专区在线播放| 男人舔女人下体高潮全视频| 99久久精品国产国产毛片| 免费av观看视频| 国产成人91sexporn| 熟女人妻精品中文字幕| 九草在线视频观看| 亚洲18禁久久av| 午夜福利成人在线免费观看| 插逼视频在线观看| or卡值多少钱| 久久精品国产鲁丝片午夜精品| 神马国产精品三级电影在线观看| av卡一久久| 免费一级毛片在线播放高清视频| 看非洲黑人一级黄片| 乱系列少妇在线播放| 白带黄色成豆腐渣| 五月玫瑰六月丁香| 国产精品蜜桃在线观看 | 99久久精品热视频| 最近2019中文字幕mv第一页| 黄色一级大片看看| 一边亲一边摸免费视频| 亚洲av不卡在线观看| 精品人妻偷拍中文字幕| 亚洲精品国产av成人精品| 在线播放无遮挡| www.色视频.com| 人妻系列 视频| 高清午夜精品一区二区三区 | 中文字幕久久专区| 在线观看一区二区三区| 亚洲国产色片| .国产精品久久| 舔av片在线| 精品久久久久久久久久免费视频| 中文欧美无线码| 久久国内精品自在自线图片| 亚洲在久久综合| 精品久久久久久久久av| 国产精品久久电影中文字幕| 九九在线视频观看精品| 免费av毛片视频| 亚洲内射少妇av| 亚州av有码| 国产高清视频在线观看网站| 可以在线观看的亚洲视频| or卡值多少钱| 久久久a久久爽久久v久久| 欧美又色又爽又黄视频| 午夜精品在线福利| 亚洲av二区三区四区| 免费观看人在逋| 久久久久久久久中文| 国产精品人妻久久久影院| 亚洲美女搞黄在线观看| 国产午夜精品久久久久久一区二区三区| 亚洲国产色片| 熟女电影av网| 欧美在线一区亚洲| 欧美日本视频| 久久精品综合一区二区三区| 午夜久久久久精精品| 精品久久久久久久久亚洲| 国产日本99.免费观看| 1000部很黄的大片| 日韩欧美精品免费久久| 国产毛片a区久久久久| 九九热线精品视视频播放| 久久久久久久亚洲中文字幕| 最新中文字幕久久久久| 免费看av在线观看网站| 午夜福利高清视频| 一本久久精品| 三级经典国产精品| 一个人免费在线观看电影| 99久久精品一区二区三区| 亚洲一级一片aⅴ在线观看| 婷婷色av中文字幕| 此物有八面人人有两片| 禁无遮挡网站| 欧美区成人在线视频| 欧美激情在线99| 熟妇人妻久久中文字幕3abv| 欧美性猛交黑人性爽| av黄色大香蕉| 99久久精品热视频| 晚上一个人看的免费电影| 免费看a级黄色片| 成人性生交大片免费视频hd| 熟女电影av网| 欧美xxxx黑人xx丫x性爽| 国产探花极品一区二区| 亚洲国产欧美人成| 啦啦啦观看免费观看视频高清| 免费电影在线观看免费观看| 成人亚洲精品av一区二区| 2021天堂中文幕一二区在线观| 亚洲国产色片| 国产午夜福利久久久久久| av视频在线观看入口| 国产精品不卡视频一区二区| 精品久久国产蜜桃| 亚洲在线自拍视频| 99久久中文字幕三级久久日本| 99热全是精品| 午夜福利高清视频| 1024手机看黄色片| 波野结衣二区三区在线| 日日摸夜夜添夜夜爱| 欧美激情国产日韩精品一区| 九九热线精品视视频播放| av专区在线播放| 中文精品一卡2卡3卡4更新| 欧美一级a爱片免费观看看| 亚洲国产欧美在线一区| 级片在线观看| 99久久精品热视频| 久久草成人影院| 欧美日韩在线观看h| 熟女人妻精品中文字幕| 悠悠久久av| 久久九九热精品免费| 亚洲欧美精品自产自拍| 最后的刺客免费高清国语| 欧美激情久久久久久爽电影| 在线观看美女被高潮喷水网站| 欧美激情国产日韩精品一区| 一级毛片久久久久久久久女| 久久久久久久久久成人| 我要搜黄色片| 日韩欧美在线乱码| 噜噜噜噜噜久久久久久91| kizo精华| 免费大片18禁| 日本熟妇午夜| 亚洲av第一区精品v没综合| 一个人看视频在线观看www免费| 女人十人毛片免费观看3o分钟| 国内精品一区二区在线观看| 校园人妻丝袜中文字幕| 日本熟妇午夜| 欧美在线一区亚洲| 亚洲av二区三区四区| 美女内射精品一级片tv| 看免费成人av毛片| 女人十人毛片免费观看3o分钟| 国产白丝娇喘喷水9色精品| 国产精品久久视频播放| 少妇高潮的动态图| 六月丁香七月| 少妇熟女aⅴ在线视频| 校园人妻丝袜中文字幕| 欧美激情久久久久久爽电影| 六月丁香七月| .国产精品久久| 亚洲va在线va天堂va国产| 亚洲电影在线观看av| 免费观看的影片在线观看| 在线观看免费视频日本深夜| av卡一久久| 免费av不卡在线播放| 一边摸一边抽搐一进一小说| 天堂中文最新版在线下载 | 人妻夜夜爽99麻豆av| 久久久国产成人免费| 天天躁夜夜躁狠狠久久av| 熟女人妻精品中文字幕| 日韩一本色道免费dvd| 日韩一区二区三区影片| 人人妻人人澡欧美一区二区| 亚洲欧美精品综合久久99| 深夜a级毛片| 成年av动漫网址| 国产精品av视频在线免费观看| 亚洲不卡免费看| 精品一区二区免费观看| 中文在线观看免费www的网站| 深夜精品福利| 亚洲欧美日韩高清在线视频| 亚洲精品粉嫩美女一区| 中文资源天堂在线| 亚洲七黄色美女视频| 少妇高潮的动态图| 亚洲五月天丁香| 99视频精品全部免费 在线| 欧美+日韩+精品| 亚洲国产精品成人综合色| 欧美区成人在线视频| 好男人在线观看高清免费视频| 中文字幕av成人在线电影| 国产单亲对白刺激| 久久精品夜夜夜夜夜久久蜜豆| 老女人水多毛片| 国产av不卡久久| 日韩成人伦理影院| 免费人成视频x8x8入口观看| 人体艺术视频欧美日本| 国产真实乱freesex|