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

移動加熱器法碲鋅鎘晶體生長系統熱場研究

2023-09-22 07:12:28徐哲人張繼軍曹祥智祁永武
人工晶體學報 2023年9期
關鍵詞:界面生長

徐哲人,張繼軍,曹祥智,盧 偉,劉 昊,祁永武

(上海大學材料科學與工程學院,上海 201900)

0 引 言

碲鋅鎘(Cd1-xZnxTe,簡稱CdZnTe或CZT)晶體由于其寬禁帶、高密度、低電離能等特性,成為一種典型的室溫核輻射探測器材料,在國土安全、醫學成像、空間探測等方面獲得了廣泛應用[1-3]。溶液生長法從富碲溶液中制備CZT單晶,改善了高溫熔體生長法導致的晶體缺陷,近年來在生長探測器級CZT單晶的研究上受到越來越多的關注[4-5]。目前廣泛使用的溶液生長方法主要有溫度梯度溶液生長法、移動加熱器法。其中移動加熱器法(traveling heater method, THM)由于具備區熔提純、低溫生長、籽晶生長等優點,成為生長高性能CZT晶體的熱門方法。

采用THM生長大尺寸、高性能CZT單晶的核心是通過調節晶體組分、生長爐溫場實現富Te熔區在生長界面均勻定向結晶,同時盡可能地減小晶體的熱應力。由于CZT晶體生長在高溫且封閉的設備中進行,缺乏直接觀測和參數監測手段,想要對晶體生長的溫度場進行嚴格的分析,并試圖找到在高溫下生長CZT單晶過程中最優的溫度梯度條件,計算機數值仿真顯得十分重要。

在過去的幾年里,一些數值研究致力于分析CdTe類半導體晶體生長過程的熱量和質量傳輸[6-7]。Chang等[8]使用了有限元穩態模型模擬分析了THM生長HgCdTe過程中的溫場輸運,通過沿坩堝壁面設置熱剖面分析了各種因素對界面形狀的影響。Kuppurao等[9]針對布里奇曼法生長CZT晶體建立了二維的準穩態過程和瞬態過程,并且利用準穩態過程研究了傳熱和對流現象,利用瞬態過程研究Zn的傳質過程[10]。Martinez-Tomas等[11]對垂直布里奇曼(vertical Bridgman, VB)法生長CdTe的整個生長系統進行了建模,采用三步法對坩堝和爐膛兩個系統進行研究,分析了生長系統中的熱量傳導。另外還有一些數值研究致力于分析在各種熱條件、坩堝旋轉或者是在靜磁場的影響下的生長界面形狀變化[12-14]。然而上述研究中,加熱器都是通過施加一個恒定速率上升的溫度分布曲線來近似模擬。實際的THM晶體生長過程中,需要考慮加熱器的位置和尺寸,以及坩堝和周圍空腔之間的熱交換對晶體生長穩態溫場的影響。因此,對THM法生長CZT晶體的溫場進行數值模擬時,需將包括爐膛結構的完整爐體建立模型,本文在這一前提下進行了數值模擬研究。

本文的研究內容主要分為兩部分:1)通過設立溫度控制點的方式逆模擬加熱源的功率;2)通過計算獲得的功率研究幾何參數等變量對THM生長CZT的影響規律。目標是通過數值模擬探究整個爐體結構(爐體幾何參數和加熱器分布),對生長系統由瞬態加熱過程轉變至偽穩態的生長過程后的爐內溫場分布,進一步對晶體生長實驗起到指導作用。

1 移動加熱器法生長過程的熱輸運原理

1.1 基本物理過程

在單晶生長爐的爐膛內存在不同的介質,如熔體、晶體、氣體、爐膛管(瑞典Canthal爐管,FeCrAl合金管)等。不同的介質具有不同的溫度,即使在同一介質中溫度的分布也不均勻,可以把在同一時刻下生長系統中的溫度的空間分布稱為溫場。單晶爐內的溫場總是存在溫度梯度,因而存在熱流。生長爐中作為熱源的加熱爐絲產生的熱流與耗散在環境中的熱流達到動態平衡時,才達到穩態溫場,否則爐溫總是向上或是向下的非穩態變化過程。本文討論了生長系統中非穩態溫場達到偽穩態狀態時爐體內的溫度分布。

圖1(a)展示了THM生長系統的主要結構,包括保溫層、加熱器、Canthal爐管、坩堝、支撐桿以及坩堝內的原料,生長系統呈圓柱軸對稱型。圖1(b)中包含了生長系統中主要的傳熱形式。按區域傳熱特點可以主要分兩部分討論。第一個是坩堝內部原料區,主要由上部的多晶區、中部的富Te溶劑區與下部的籽晶區組成。第二個是坩堝外部加熱區,由保溫層、加熱器、Canthal爐膛、散熱區等組成。實際生長過程中,整個生長系統密封于不銹鋼金屬外殼中,實際對流由金屬外殼內壁與爐膛外壁間溫差所引起,因此在模型中環境對流的影響可以忽略不計。

圖1 THM生長系統結構(a)和生長系統內熱輸運方式(b)示意圖Fig.1 Schematic diagram of THM growth system structure (a) and heat transport in the growth system (b)

本文研究THM生長系統中,加熱器的能量通過熱輻射、熱傳導以及熱對流的方式傳遞至坩堝中,CZT晶體溶解在富Te的CZT溶液中,并在溫度梯度和濃度梯度的驅動下形成流動和擴散,最終運輸至生長界面完成單晶生長過程。晶體生長包含吸熱及放熱過程,這些機制相互耦合最終將生長系統控制在一個穩定的溫度范圍內。

1.2 傳熱過程的物理和數學描述與分析

晶體生長系統內的溫度最高能達到將近1 200 K,且靠近加熱器的區域溫度最高。THM生長系統中的輻射傳熱主要存在于兩個區域:1)加熱器與Canthal爐管外壁以及保溫棉外壁之間;2)坩堝外壁以及支撐桿外壁與Canthal爐管內壁之間。加熱器與Canthal爐管外壁之間的輻射傳熱量(Qradiheat)可由公式(1)計算。

(1)

式中:Akan為Canthal爐管外壁受輻射面積,εkan為Canthal爐管壁面發射率,Theater為加熱器溫度,Tkan為Canthal爐管壁面溫度,σ為玻爾茲曼常數。同理可以估算坩堝外壁與Canthal爐管內壁之間的輻射熱量(Qradi_amp)為

(2)

式中:Aamp為坩堝外壁受輻射面積,εamp為坩堝壁面的表面發射率,Tamp為坩堝壁面溫度。

對于輻射計算,能量方程形式如公式(3)。

(3)

式中:ρcP為熱容,k為熱導率。在此模型中可將所有固體物質視為不透明物體,所有參與輻射的壁面均為漫灰表面,采用面對面輻射模型,輻射熱流以第二類邊界條件的形式參與計算,每若干個迭代步更新一次輻射在邊界處的熱流,因此需要計算每個參與輻射表面的凈輻射熱流密度。使用參與輻射表面的表面溫度計算凈輻射熱流密度q。

(4)

式中:δij為克羅內克函數(Kronecker delta);Fij為視角系數,表示從j面發射的能量被i面攔截的份額,用式(5)計算。

(5)

式中:Ai、Aj分別為微元面i、j的面積,Rl為i面與j面心連線長度,θi、θj為該連線與i面、j面的外法向量夾角。

晶體生長爐內氣相組分流動的熱輸運強度可用無量綱量瑞利數Ra(Rayleigh number)表征,表達式為

(6)

式中:cp為氣體恒壓比熱容,μ為動力黏度,g為重力加速度,ν為運動黏度,β≈ 1/T為理想氣體熱膨脹系數,D為Canthal爐管內徑,ΔT為Canthal爐管與坩堝之間的溫差。

THM生長CZT晶體主要包括兩個過程,即CZT多晶在熔解界面處的分解及CZT溶液在生長界面處凝固成CZT單晶。這兩個過程存在較為復雜的表面相變反應,因此可以用式(7)估算熔解界面與生長界面處的熱量(Qchem)變化。

Qchem≈Qlatent=ρcrystalvcrystalAinterfaceHVS

(7)

式中:Qlatent為晶體生長相變釋放的熱量,ρcrystal為晶體密度,vcrystal為CZT單晶生長速率,Ainterface為界面面積,HVS為相變潛熱。潛熱必須滿足公式(8)。

(8)

式中:Ks和Kl分別為固體及液體內的導熱系數,dT/dz為軸向溫度梯度,ρ為生長速率。

計算使用的參數如表1所示。

表1 幾何和物性等計算參數[15-17]Table 1 Material properties and parameters for computation[15-17]

本文算例的計算域均為軸對稱,半徑約500 mm,高約1 280 mm,采用商業軟件Meshing對計算域劃分結構化網格,網格尺度為2 mm。采用Fluent軟件對方程進行離散和求解,對稱軸邊界條件為溫度零梯度,外邊界為400 K定溫。

2 結果與討論

2.1 逆模擬模式

熱源功率的設置是生長爐溫度場全局模擬的難點問題。實際生長爐多是通過多段溫控單元,基于熱電偶實時測溫并在PID控制算法的作用下,對溫場分布進行控制,這表明熱源功率是與溫度有關的變量。德國Fraunhofer 集成系統與設備技術研究所晶體生長實驗室開發了CrysMAS代碼,使用一種逆模擬模式計算熱源的輸入功率,從而實現多溫區VB爐的傳熱模擬[18]。Stelian等[19]使用此方法研究了CdTe與Ge晶體生長實驗中具有代表性的兩種VB爐的溫度場。Derby等[20]使用此方法研究了電動梯度冷凝法生長CZT晶體。本文采用Fluent商用軟件通過逆模擬模式實現了多溫區THM生長爐的全局傳熱模擬,計算了不同加熱器數量對空爐溫場分布的影響,并與實際晶體生長爐測溫溫場分布進行對比。

圖2為THM晶體生長爐結構與溫控點位置示意圖。逆模擬模式主要通過設定溫度控制點的方式實現,T1、T2、T3、T4、T5分別設置為893、953、1 073、923以及913 K,通過Fluent User-Defined Function(udf) 功能對Fluent軟件進行二次開發,使其在每次迭代后將控制點的溫度實時反饋到主計算節點,再根據控制點溫度值對熱源的輸入功率進行調整。在后續的研究中使用距離中心z軸0、16、25 mm的線段d0、d1、d2軸向溫度(B點至A點溫度分布)近似表示爐膛內溫場分布。

圖2 THM晶體生長爐結構與溫控點位置示意圖,其中B(z=0.70)至A(z=0.30)為有效溫度區間,距離中心z軸0、16、25 mm的線段d0、d1、d2處的軸向(B→A)溫度用來表示爐膛內部溫場分布Fig.2 Schematic diagram of the structure of the THM crystal growth furnace and the location of temperature control points, where B(z=0.70) to A (z=0.30) is the effective temperature range, and the axes at d0, d1, and d2 of the line segments 0, 16 and 25 mm away from the z-axis toward (B→A) temperature is used to approximate the temperature field distribution inside the furnace

為了探究加熱器分布對爐膛內溫度分布的影響,對比設置了三種不同的算例,分別為:1)僅使用3號熱源;2)使用2、3、4號熱源;3)使用1、2、3、4、5號熱源。圖3為三種算例下爐膛軸心處從控制點1(z=0.96 m)至控制點5(z=0.32 m)距離的溫度分布圖。

圖3 在不同加熱器數量條件下爐膛軸心溫度分布Fig.3 Temperature distribution of furnace axis under the condition of different number of heaters

僅使用3號熱源的條件下,控制點3(z=0.64 m)處溫度達到設定溫度T3(1 073 K);使用2、3、4號熱源的條件下,控制點2(z=0.80 m)、3(z=0.64 m)、4(z=0.48 m)處溫度達到設定溫度T2(953 K)、T3(1 073 K)、T4(923 K);使用1、2、3、4、5號熱源的條件下,控制點1(z=0.96 m)、2(z=0.80 m)、3(z=0.64 m)、4(z=0.48 m)、5(z=0.32 m)處溫度達到設定溫度T1(893 K)、T2(953 K)、T3(1 073 K)、T4(923 K)、T5(913 K)。結果表明,在逆模擬模式下,通過增加加熱器數量、設立多個溫度控制點,可以達到理想的溫度分布曲線。

實際的THM生長CZT晶體生長過程中,石英坩堝的移動范圍在A(z=0.30)至B(z=0.70)這40 cm長的區間,因此重點模擬這部分區域的溫度分布。圖4(a)~(d)分別為逆模擬模式下控制點溫度隨迭代步變化曲線,加熱器功率隨迭代步變化曲線,d0、d1、d2處軸向溫度計算值與d0、d1、d2處軸向溫度實驗測量值。圖4(a)、(b)結果表明,在同一迭代步下控制點的溫度值同步影響著加熱器的功率值,隨著迭代步的增進,加熱器功率以及控制點的溫度都趨于收斂并最終穩定在了固定值。圖4(c)、(d)結果表明,當功率分別為225.6、343.7、1 045.9、92.5、199.6 W時,使用5個加熱器的溫度場模擬結果與實驗測得的爐膛內部溫場分布在趨勢上基本一致。由模擬結果可知,與d0處的溫度曲線相比,d1與d2處的溫度最高點的位置分別向爐膛頂端偏移了19、25 mm,并且最高溫度也分別上升了11.61、36.84 K,表明爐膛軸心處溫場發生偏移。

圖5(b)、(c)分別為THM晶體生長爐爐膛中部的溫度云圖和空氣流動速度矢量圖。圖5(b)顯示,爐膛壁面處靠近加熱器部分溫度高于靠近散熱區部分,在區域交接處溫度梯度達到最大。圖5(c)中空氣流動速度矢量圖則顯示在爐膛中部存在兩種流動方向不同的渦流。這是由于空氣密度受溫度影響出現差異而產生了對流。Canthal爐管靠近加熱器部分溫度高于靠近散熱區部分溫度,且爐膛軸心處溫度介于兩者之間,因此在爐膛上下部出現了流動方向相反的兩種渦流。可以得出,對流效應是使爐膛軸心處溫場發生偏移最主要的原因。爐膛中空氣的流動對于爐膛內溫場的分布有著重要的影響,在逆模擬模式下充分考慮了空氣流動對溫場的影響。

2.2 CZT晶體生長系統熱場研究

本節基于數值模擬方法對THM生長CZT晶體的熱場進行研究,提出熱場設計原則,以模擬生長直徑45 mm、長度150 mm的CZT晶體為目標,探究結構參數如加熱模塊與Canthal爐管間距離、散熱區長度等對熱場的影響規律。

在THM生長CZT晶體過程中,軸向溫度梯度與晶體生長速度有關。界面前沿處液體由于溶質的析出使凝固溫度偏離液相線溫度,能夠穩定生長的最大生長速度(v)可由式(9)計算。

(9)

式中:D為擴散系數,m為液相線的斜率,Cl為界面處液相的濃度,K為分凝系數,Gl為液相區域垂直于界面的溫度梯度。最大晶體生長速率取決于液相區域中溫度梯度的大小,提高液相區域中的溫度梯度可以加快晶體生長速率。然而,Wang等[21]研究發現,較大的溫度梯度對晶體質量產生不利的影響。因此,THM晶體生長模擬需要研究富Te溶液和CZT晶體中的溫度梯度與生長系統中關鍵部位的尺寸參數之間的關系。

CZT晶體生長系統模型如圖6所示。坩堝內部熔解界面和結晶界面處溫度梯度主要受加熱源和散熱區影響,因此將加熱源與Canthal爐管距離S、散熱區長度l作為研究變量。熱源與Canthal爐管間距離基礎值S0=10 mm,散熱區寬度基礎值l0=50 mm。采用單一控制變量法進行研究,設置3組算例,其中算例1為供對照的基礎算例,算例2~3分別研究不同設計參數對坩堝內溫度分布的影響,具體參數設置如表2所示,參數變化范圍見圖6中虛線框域。計算結果重點考察生長界面和熔解界面處晶體中的溫度梯度GS1、GS2,熔體中的溫度梯度Gl1、Gl2,坩堝壁面處溫度分布,及熔區內最大溫差變化。

表2 算例設計說明Table 2 Settings of the cases

圖6 生長系統結構和變量說明,散熱區長度l變化范圍為30 mm至80 mm,加熱器與Canthal爐管距離變化范圍為5 mm至15 mm。坩堝壁面溫度分布為A點至D點溫度曲線,Gs1、Gl1、Gs2、Gl2分別為B點C點兩側溫度梯度Fig.6 Growth system structure and variable description, the length D of the cooling zone varies from 30 mm to 80 mm, and the distance between the heater and the Canthal furnace tube varies from 5 mm to 15 mm. The temperature distribution of the crucible wall is a temperature curve from point A to point D, and Gs1, Gl1, Gs2, and Gl2 are the temperature gradients on both sides of point B and point C respectively

2.2.1 加熱器與Canthal爐管距離變化對坩堝溫場變化影響

圖7(a)為加熱器和Canthal爐管距離變化對坩堝壁面溫度的影響規律圖,熔體區與晶體區呈現明顯的溫度梯度差異,熔解界面與結晶界面溫度保持在1 083 K,山峰型溫度分布符合THM生長CZT單晶的要求。相較于基礎算例,由于視角系數Fij的變化,加熱器與Canthal爐管距離增大時整體溫度下降,相反,當加熱器與Canthal爐管距離減小后整體溫度升高。當距離由5 mm增大至10、15 mm,中間加熱器相對Canthal爐管的視角系數Fij分別為0.865、0.778和0.725,1 083 K長晶溫度之間的Te熔區寬度分別減小20.5%、37.5%,最高溫度減小3.7%、5.6%。圖7(b)為加熱器和Canthal爐管距離變化對熔解界面晶體側和熔體側溫度梯度Gs1、Gl1,生長界面晶體側和熔體側溫度梯度Gs2、Gl2,熔區內最大溫差的影響規律圖。從圖7(b)可見,隨著Canthal爐管與加熱器間距離增大,生長界面與熔解界面的溫度梯度不斷減小。相較于熔體區,晶體區的溫度梯度變化更為明顯。并且Canthal爐管與加熱器間距離越大,熔體區最大溫差越低。當Canthal爐管與加熱器間距離增大10 mm,熔體區內溫度相差40.5 K。具體數據如表3所示。

表3 加熱器和Canthal爐管距離改變算例參量變化Table 3 Investigated parameter variation when change the distance between the heater and the Canthal furnace tube

圖7 加熱器和Canthal爐管距離變化對不同溫場參量的影響規律。(a)坩堝壁面溫度;(b)熔解界面晶體側和熔體側溫度梯度Gs1、Gl1,生長界面晶體側和熔體側溫度梯度Gs2、Gl2,熔區內最大溫差Fig.7 Effect of the change of the distance between the heater and the Canthal furnace tube on different temperature field parameters. (a) temperature of the crucible wall; (b) crystal side temperature gradient Gs1 at melting interface, melt side temperature gradient Gl1, growth interface crystal side temperature gradient Gs2, melt side temperature gradient Gl2 and maximum temperature difference in melting zone

2.2.2 散熱區長度對坩堝溫場變化影響

圖8(a)、(b)分別為散熱區長度變化對坩堝壁面溫度的影響規律與散熱區長度變化對熔解界面晶體側溫度梯度Gs1、熔體側溫度梯度Gl1,生長界面晶體側溫度梯度Gs2、熔體側溫度梯度Gl2以及熔區內最大溫差的變化規律。由圖8(a)可見,散熱區寬度減小后整體溫度同步上升,散熱區寬度增加后整體溫度同步下降。計算結果顯示,當散熱區寬度由30 mm增大至50、80 mm,1 083 K長晶溫度之間的Te熔區寬度分別減小32.7%、50.0%,最高溫度分別減小3.6%、4.6%。由圖8(b)可見,溫度梯度呈現與圖7(b)相反的變化規律,隨著散熱區寬度增大,生長界面與熔解界面兩側溫度梯度不斷增大。相較于更改加熱器與Canthal爐管之間的距離,改變散熱區寬度對溫度梯度的影響較小。同時,散熱區寬度越大,熔區內溫差越小。當散熱區寬度增大50 mm,熔體區內溫差減小7.2 K。具體數據如表4所示。

表4 散熱區長度改變算例參量變化Table 4 Investigated parameter variation when change the heat dissipation zone length

圖8 散熱區長度變化對不同溫場參量的影響規律。(a)坩堝壁面溫度;(b)熔解界面晶體側溫度梯度Gs1、熔體側溫度梯度Gl1,生長界面晶體側溫度梯度Gs2、熔體側溫度梯度Gl2,以及熔區內最大溫差Fig.8 Effect of the change of the heat dissipation zone length on different temperature field parameters. (a) Temperature of the crucible wall; (b) crystal side temperature gradient Gs1 and melt side temperature gradient Gl1 at melting interface, crystal side temperature gradient Gs2 and melt side temperature gradient Gl2 at growth interface, and maximum temperature difference in melting zone

3 結 論

1)THM生長CZT晶體的生長系統中,對于考慮坩堝與坩堝外部環境進行熱交換的數值模擬,加熱器功率的設置可以通過逆模擬模式近似模擬,通過設定控溫點的方式對整體爐膛溫場分布有較好的數值模擬結果。本文將加熱器功率分別設置為225.6、343.7、1 045.9、92.5、199.6 W,得到的爐內溫場可以較好地模擬實際溫場分布。

2)THM生長CZT晶體的生長系統中,以模擬生長直徑45 mm、長度150 mm的CZT晶體為目標,改變Canthal爐管與加熱器間距離因改變了輻射傳熱中的視角因子,對于溫場的影響較大,在距離由5 mm增大至10、15 mm后,最高溫度減小3.7%、5.6%。而改變散熱區寬度對調控原料區熔區寬度則更具實際效果,散熱區寬度由30 mm增大至50、80 mm后,熔區寬度分別減小32.7%、50.0%。界面處溫度梯度與散熱區寬度成正相關,與Canthal爐管與加熱器間距離成負相關。

猜你喜歡
界面生長
碗蓮生長記
小讀者(2021年2期)2021-03-29 05:03:48
國企黨委前置研究的“四個界面”
當代陜西(2020年13期)2020-08-24 08:22:02
共享出行不再“野蠻生長”
生長在哪里的啟示
華人時刊(2019年13期)2019-11-17 14:59:54
野蠻生長
NBA特刊(2018年21期)2018-11-24 02:48:04
生長
文苑(2018年22期)2018-11-19 02:54:14
基于FANUC PICTURE的虛擬軸坐標顯示界面開發方法研究
空間界面
金秋(2017年4期)2017-06-07 08:22:16
電子顯微打開材料界面世界之門
人機交互界面發展趨勢研究
主站蜘蛛池模板: 91视频国产高清| 国产91小视频| 波多野吉衣一区二区三区av| 欧美日韩在线观看一区二区三区| 国产剧情国内精品原创| 国产女人爽到高潮的免费视频| 99人妻碰碰碰久久久久禁片| 全免费a级毛片免费看不卡| 在线观看国产黄色| 亚洲黄色网站视频| 日本影院一区| 亚洲制服中文字幕一区二区 | 中国丰满人妻无码束缚啪啪| 成人精品视频一区二区在线| 久久伊人色| 日韩一二三区视频精品| 日本不卡在线视频| 手机精品视频在线观看免费| 亚洲综合久久成人AV| 99热这里只有精品久久免费| 无码网站免费观看| 日韩精品少妇无码受不了| 日韩黄色大片免费看| 国产成人亚洲欧美激情| 日本一本在线视频| 片在线无码观看| 国产原创第一页在线观看| 永久毛片在线播| 亚洲国产成人精品一二区| 国产香蕉一区二区在线网站| 97一区二区在线播放| 高清色本在线www| 视频在线观看一区二区| 尤物精品视频一区二区三区| 成人国产精品2021| 亚洲精品va| 亚洲男人天堂网址| 国产精品自拍露脸视频| 亚洲欧美日韩久久精品| 性色一区| 国产91麻豆免费观看| 国产精品19p| 干中文字幕| 欧美怡红院视频一区二区三区| 久久综合亚洲鲁鲁九月天| 国产成年女人特黄特色毛片免 | 国产精品免费露脸视频| 欧美成a人片在线观看| 麻豆精品在线| 天堂在线亚洲| 国产成人麻豆精品| 亚洲熟妇AV日韩熟妇在线| 午夜成人在线视频| 又黄又爽视频好爽视频| 91区国产福利在线观看午夜| 国内精品91| 青草精品视频| 国国产a国产片免费麻豆| 在线观看热码亚洲av每日更新| 国产真实乱人视频| 91精品免费久久久| 成人国产免费| 中文字幕亚洲电影| 久久久精品久久久久三级| 成人看片欧美一区二区| 国产无码制服丝袜| 午夜日b视频| 久久国产精品77777| 九九九精品成人免费视频7| 六月婷婷精品视频在线观看| 婷婷六月激情综合一区| 六月婷婷精品视频在线观看| 大乳丰满人妻中文字幕日本| 日韩视频福利| 国产人成乱码视频免费观看| 毛片卡一卡二| 国产视频a| 精品亚洲国产成人AV| 久久性妇女精品免费| 亚洲综合九九| 欧美国产综合色视频| 亚洲天堂久久久|