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

φ300 mm直拉硅單晶生長過程中的變晶現象及其影響因素

2022-08-12 02:10:14晶,劉
人工晶體學報 2022年7期
關鍵詞:界面

張 晶,劉 丁

(西安理工大學晶體生長設備及系統集成國家地方聯合工程研究中心,西安 710048)

0 引 言

作為一種可持續發展的綠色行業,光伏產業正受到越來越多的關注。近年來,隨著電池生產技術的提升,光伏產業不斷朝著降低成本、減少硅材料的損耗和提高生產效率的方向發展。直拉(Czochralski, Cz)硅單晶擁有生產成本低,晶體質量和性能好的優點,目前已廣泛應用于光伏產業領域[1-3]。在直拉法制備太陽能級硅單晶中,為滿足行業需求,晶體生長熱系統的設計朝著提升晶體生長速度和降低加熱器功率的方向發展。其中,提高晶體生長速度主要是通過提高硅單晶的提拉速度來實現的。據統計,國際上生長200 mm硅單晶的提拉速度已超過1.8 mm/min,直徑為300 mm硅單晶最高提拉速度已經達到1.5 mm/min[4],硅單晶的生產效率有了明顯提高。然而,過高的提拉速度會給晶體生長過程帶來問題,主要體現在晶體等徑生長階段[5]。

等徑生長階段硅單晶以均勻直徑生長,其生長時間占晶體生長全周期的80%以上,該階段所生長的晶體直徑近似相等、成分均勻、結構完整、缺陷較少,是晶體生長的主要階段。等徑精度越高,晶體的質量越好。然而,在拉晶實驗中觀察到,當提拉速度較快,等徑生長的晶體沒有按照圓柱形軸對稱的晶格結構生長,出現扭曲生長或者螺旋生長的現象[6],這種現象稱為晶體的“扭晶”。扭晶的出現使晶體的部分晶格在特定方向產生塑性變形,并且變形區原子與未變形區原子在交界處仍維持緊密接觸,即生長的晶體由單晶體變為多晶體,發生變晶[6],極大地降低了晶體品質。由晶體生長實驗獲得扭晶后的硅單晶棒如圖1所示,可以看到晶體直徑呈現劇烈變化,光滑的晶體側面呈現波浪狀,晶體生長直徑不均勻。在合理范圍內調控提拉速度是解決扭晶現象的有效方法,然而降低提拉速度不僅會降低生產效率,增加生產成本,也使晶體內部缺陷分布發生變化,不利于晶體的生長。因此,深入研究扭晶現象的成因并避免這一現象發生,對直拉硅單晶工藝中尋求最大提拉速度并降低成本具有一定指導意義。

與上述扭晶現象類似的情況還發生于化合物晶體生長研究中。文獻[7]對不同加熱條件下的BGO晶體形狀進行研究,結果表明降低熔體自由表面的熱輻射有利于改善晶體形狀。Okan等[8]解釋了DyAlO3晶體在生長過程中比HoAlO3更容易扭晶。文獻[9]表明在直拉法生長鈣鎂鋯摻雜釓鎵石榴石(SGGG)晶體過程中,提升坩堝位置有助于提高固液界面溫度梯度。Kamada等[10]探討了由于熔體對流不穩定所引起的Cz法生長Pr∶LuAG晶體時非對稱溫度分布現象。熔體對流的不穩定是影響晶體扭晶的一個原因,而建立熔體穩定性與晶體直徑的變化關系十分困難。主要原因是熔體對流的不穩定直接影響熔體內溫度分布,進而影響固液界面處的溫度梯度,從而影響晶體生長直徑的變化。

為解決此問題,本文采用有限元法[11]對硅單晶生長過程進行數值模擬,通過數值模擬和理論相結合的方法,分析提拉速度過大時晶體生長出現扭晶的成因[12],建立晶體直徑與熔體溫度分布的關系。通過對影響因素的分析,提出改進工藝參數措施。通過仿真和實驗驗證該方法的有效性,為在晶體生長熱場中設計最大提拉速度以及尋求最佳工藝條件提供參考。

1 數學模型

采用有限元數值模擬的方法研究硅熔體的溫度分布和對流情況,是目前分析硅單晶典型生長過程的有效方法[13-16]。在仿真過程中建立單晶爐熱場準穩態數學模型,為了方便數值計算,對其進行軸對稱簡化。對整個單晶爐進行二維(2D)建模,包括熔體、晶體、熱屏、氣氛、石英坩堝、石墨坩堝、保溫套、加熱器、爐體等,如圖2所示。

(1)

式中:ρl為硅熔體的密度;v表示速度矢量;μ為熔體的粘滯系數;β為溫度引起的體膨脹系數;T為溫度;下標s表示固體;g為重力加速度;P為壓力;式(2)為連續性方程:

(2)

熔體的熱傳輸方程如式(3)所示:

(3)

式中:c為熱容量;kl為硅熔體的熱傳導系數。晶體與熔體固液界面處的熱平衡滿足能量守恒方程,如公式(4),表示流出固液界面的熱量Q3等于流入界面的熱量Q2與界面處相變所釋放的潛熱Q1之和:

(4)

式中:ks為硅單晶熱傳導系數;ρs為硅單晶的密度;A為硅單晶橫截面積;n為固液界面的法線向量;(?T/?n)l和(?T/?n)s分別表示相變界面附近熔體的溫度梯度和相變界面附近晶體的溫度梯度;ΔH為晶體潛熱。

所建立的局部三維(3D)模型包括晶體、熔體、石英坩堝和石墨坩堝,對于三維局部模型來說,邊界條件的設定影響仿真的精度,其中需要設定的邊界條件主要包括:熔體自由表面、晶體外表面、石英坩堝內表面、石英坩堝頂表面、石墨坩堝頂表面。由于晶體的外表面與熔體自由表面溫度分布受諸多因素影響,邊界條件不確定。原因是晶體外表面和熔體自由表面的溫度分布一方面受熱輻射影響,另一方面受氬氣流動影響[17]。在熱系統熱屏的作用下,熱輻射及氬氣流向改變,因此在晶體外表面和熔體自由表面處采用第一類邊界條件進行溫度約束是不恰當的。為了解決此問題,建立熱流密度方程來描述晶體外表面和熔體自由表面由于受到氬氣吹拂和輻射產生的熱損耗。

對晶體外表面處建立的熱流密度邊界條件如式(5)所示:

(5)

式中q′s為:

q′s=qout,k-qin,k=σεT4-εqin,k,qin,k=sumj=1~N(Fkjqout,j)

(6)

熔體自由表面處的熱損耗如式(7)所示:

(7)

式中q′l為:

q′l=qout,k-qin,k=σεT4-εqin,k
qin,k=sumj=1~N(Fkjqout,j)

(8)

式(5)和式(7)左側表示熱流密度邊界條件,δ[T(x)-T0(x)]1.25表示因氣體對流而產生的熱損耗,?T/?r,?T/?z為徑向溫度梯度和軸向溫度梯度,q′s表示晶體表面由輻射而產生的熱損耗,q′l表示熔體液面由輻射而產生的熱損耗,ra和l分別為晶體直徑和熔體深度,qout,k是第k個表面上由于熱輻射向外散發的熱量,qin,k是第k個表面上由于外部熱輻射而吸收的熱量。式(6)和式(8)與輻射模型相關,在本文中考慮晶體表面與熔體自由表面的自身輻射熱流、反射熱流和入射輻射熱流。δ為對流的熱損耗系數,σ為Stefan-Boltzmann常數,ε為輻射系數,Fkj為k、j兩個表面之間的角系數[18]。

石英坩堝內表面、石英坩堝頂表面和石墨坩堝頂表面的邊界條件同晶體,均采用熱流密度邊界條件。此外,在仿真過程中,密度采用Boussinesq模型近似,該簡化下,控制方程中的密度項只在重力項中考慮密度變化,而其他項中的密度取常數,即平均值ρ0。可近似表示為如公式(9)所示:

Δρ=ρ0(1-βΔT)

(9)

式中:Δρ為密度;β為溫度引起的體膨脹系數;ΔT為溫差。

在熔體自由表面采用自由滑移壁面,其他壁面均采用非滑移。

方程組(1)~(4)采用有限元方法離散,ELEMENTS數量為8 237,NODES數量為17 993。熱場的數值計算以及對固液界面耦合邊界處理已發表在文獻[19-21]中,壓力-速度修正采用SIMPLEC算法。對于幾何形狀規則的區域采用結構化網格進行離散,對于氬氣流動域這樣不規則的復雜區域采用三角形非結構化網格離散。determinant方法是有限元分析中評價網格質量常用的方法之一,對其網格劃分的雅可比值進行比較,該比值基于雅可比矩陣的行列式。本文采用determinant方法檢測網格質量,在分析過程中,采用雅可比矩陣將單元矩陣從理論形狀轉換為實際形狀。理想各單元的雅可比值為1,該值離1越遠表示單元網格質量越差。對上述離散的網格進行評價,所得雅可比值最小處為0.885,最大處為0.907,表明網格劃分效果較好。

2 結果與討論

在實際拉晶過程中,提拉速度越大越有利于加快生產進程、提高生產效率。然而,提拉速度受到晶體生長環境的種種約束,由能量守恒定理可知[22],晶體生長速度與固液界面溫度梯度和結晶潛熱有關,由式(4)可知晶體生長速度為:

(10)

式中:Gl和Gs分別表示相變界面附近熔體的溫度梯度和相變界面附近晶體的溫度梯度;由式(10)可知,當Gl=0時,晶體的生長速度可以達到極限值[22]:

(11)

在受控的直拉晶體生長中,提拉速度決定著晶體的軸向生長速度,即提拉速度的快慢決定了晶體生長的快慢[23]。但是由于晶體生長環境熱場復雜,熔體處的溫度梯度不能滿足Gl=0的條件,晶體生長速度不可能取得極限值。并且過大的提拉速度會導致晶體扭晶,由式(10)可知,晶體能夠獲得的最大生長速度取決于晶體和熔體中的溫度梯度的大小,所以對不同提拉速度下熔體內溫度分布進行分析,以尋找扭晶出現的原因和最大提拉速度。為了分析這一問題,仿真過程以西安理工大學某型號單晶爐24英寸(1英寸=2.54 cm)熱場為模型。在等徑階段,測量得到晶體發生扭晶時的工藝參數l=180 mm,φ=300 mm,其中,l為晶體長度,φ為晶體直徑。針對不同提拉速度下的晶體生長過程中的熔體溫度、流形分布等參數進行分析,模擬過程所采用的晶體生長過程控制參數為:晶體直徑300 mm,晶體長度180 mm,提拉速度0.3~1.2 mm/min,晶體轉速4 r/min,坩堝轉速8 r/min,氬氣流量0.001 m3/s,爐內壓力約1 500 Pa。

由圖3可知,當提拉速度為0.45 mm/min時,熔體自由表面溫度均高于固液界面溫度(1 685 K)。隨著提拉速度逐漸增大,熔體自由表面的平均溫度逐漸降低。當提拉速度為0.5 mm/min時,熔體自由表面靠近三相點處的溫度達到晶體的熔點;當提拉速度為1.2 mm/min時,熔體自由表面靠近三相點處溫度為1 680 K,說明在熔體自由表面存在比硅熔點(1 685 K)低的溫度區域。由式(10)可知,當提拉速度加快時,為了維持晶體生長的相對穩定性,結晶速率也相對提高,導致晶體固液界面附近軸向溫度梯度增大。而熱系統無法提供足夠的熱量,于是在熔體的自由表面三相點附近的區域出現了過冷區。出現過冷區的區域開始結晶,影響晶體生長直徑變化,而持續的直徑變化就會使晶體產生扭晶現象。觀察圖3固液界面形狀,提拉速度由0.45 mm/min升高至1.2 mm/mm的過程中,固液界面形狀逐漸凸向晶體。這是由于熔體向晶棒側傳熱效率增強,晶體側壁生長速度快于中心位置,生長界面從凸向熔體變得凸向晶體。

提拉速度分別為0.45 mm/min和1.2 mm/min時,二維仿真和三維仿真下熔體自由表面的徑向溫度梯度如圖4所示。觀察圖4可得,在同一仿真方法下,不同提拉速度得到的熔體自由表面徑向溫度梯度分布趨勢相似。但是對比同一提拉速度下不同仿真方法得到的結果,當提拉速度為0.45 mm/min時,不論是二維仿真還是三維仿真,熔體自由表面處的溫度都大于晶體熔點,表示沒有過冷區的產生。但是當提拉速度為1.2 mm/min時,三維仿真結果出現了低于晶體熔點的過冷區,即提拉速度過大時熔體自由表面會產生過冷區。而二維仿真的結果中熔體自由表面沒有出現低于晶體熔點的溫度,即沒有過冷區的產生,與圖3中提拉速度為1.2 mm/min時產生過冷區的現象不符。該結果表明了三維仿真的必要性,通過三維仿真,在提拉速度過大時產生的過冷區才可以在熔體自由表面的徑向溫度梯度分布中被明顯地觀察到。

根據上述仿真結果,可以根據熔體自由表面的溫度分布預測出晶體生長熱系統最大提拉速度,將熔體自由表面的溫度與提拉速度之間的關系作為判斷熱系統最大穩定提拉速度的依據。通過該方法可以判斷晶體生長到各個階段所能獲得的最大提拉速度,從而避免扭晶現象的發生,設計出更為優良的熱系統。

對不同提拉速度下熔體自由表面過冷區和固液界面形狀進行仿真,結果如表1所示。由表1可知,隨著晶體提拉速度增加,熔體自由表面附近的過冷區逐漸擴大,固液界面形狀從凸向熔體逐漸變為凹向熔體。當V>0.45 mm/min時過冷區出現,并隨著提拉速度的增加而逐漸擴大,固液界面形狀越凹向熔體;當V≤0.45 mm/min時,熔體自由表面不存在過冷區,隨著提拉速度的減小固液界面形狀越凸向熔體。于是,可以判定在晶體生長等徑長度為180 mm時可以獲得的最大提拉速度為0.45 mm/min。當提拉速度大于0.45 mm/min時,熔體自由表面產生過冷區,過冷區的持續存在使晶體直徑猛然增大,晶體內部溫度梯度分布劇烈變化,產生較大的熱應力。當熱應力超過硅的臨界應力后引起位錯,導致晶體變晶斷苞產生扭晶,等徑生長失敗。

表1 不同提拉速度下的熔體自由表面過冷區和固液界面形狀仿真結果Table 1 Simulation results of free melt surface supercooling zone and the solid-liquid interface shape at different values of pulling speed

根據上述分析,過大的提拉速度使得晶體產生扭晶現象,導致晶體變晶,拉晶失敗。為了得到結構完整的合格晶體,提拉速度不能超過各個生長階段的臨界值,該臨界值取決于熱系統中材料的性質和生長工藝參數。由表1可知,增大晶體提拉速度,固液界面的形狀逐漸凸向晶體,即凹向熔體。其原因是根據固液界面能量守恒定理,流入固液界面的熱量等于流出固液界面的熱量與晶體結晶時釋放的潛熱之和,如式(4)所示。當提拉速度增大,為了維持晶體的穩定生長,固液界面結晶速率也相應增大,導致單位時間內由于相變產生的結晶潛熱增多,進而影響固液界面形狀。并且晶體生長過程中固液界面的形狀與熔體熱傳輸均勻性密切相關,過大的提拉速度導致熔體溫度分布不均勻,也影響固液界面形狀,降低生長質量。因此,通過判定固液界面形狀是獲得晶體生長最大提拉速度的又一判據。合理優化熱系統的溫度梯度并且獲得晶體生長過程中的最大生長速度及其臨界值的判斷條件,為設計高效熱場提供依據,對實際拉晶實驗有重大指導意義。

3 改進措施

由上述分析可知,過大的提拉速度產生扭晶,是因為提拉速度過大時熔體自由表面產生了過冷區。當提拉速度較高時,固液界面前端的結晶速率也相應提高,由固液界面熱量平衡關系可得,晶體固液界面附近軸向溫度梯度增大,而熱系統無法提供足夠的熱量,于是在熔體自由表面三相點附近的區域出現了過冷區,晶體生長固液界面形狀由凸向熔體變得凹向熔體,使得晶體扭晶。所以扭晶現象與熔體自由表面三相點附近過冷區的出現和固液界面形狀有關,通過改善三相點附近溫度分布和固液界面形狀,可以避免扭晶的產生。在晶體生長工藝參數中,晶體旋轉會產生強迫對流,影響熔體自由表面溫度分布和固液界面的溫度分布,適當的提高晶體旋轉的速度,可以改善三相點附近溫度分布以及固液界面形狀[24]。

為了驗證不同晶體旋轉速度對熔體自由表面溫度分布和固液界面形狀的影響,通過數值模擬得到了不同轉速下晶體和坩堝內的溫度分布曲線。圖5為坩堝旋轉速度為0,晶體旋轉速度(ωs)分別為0 r/min、5 r/min和10 r/min時晶體和坩堝內的溫度分布。觀察仿真結果可得,晶體旋轉使熔體溫度分布逐漸均勻。隨著晶體旋轉速度的增加,熔體內部等溫線向晶體下方移動,熔體內部強迫對流增強,使得晶體下方高溫熔體運動增強,熔體內部溫度分布更加均勻。并且隨著晶體旋轉速度的增加,晶體固液界面的形狀得到改善,固液界面變得更加平坦,更有利于晶體生長。

在坩堝轉速為0,晶體旋轉速度分別為0 r/min、5 r/min和10 r/min時固液界面附近硅熔體徑向溫度分布曲線如圖6所示。由圖6可以觀察到,不論晶體旋轉速度多大,固液界面附近熔體溫度分布呈現兩側溫度高,中心軸處溫度低的現象,符合晶體生長工藝。在晶體旋轉速度ωs=0 r/min時,相變界面附近熔體溫度低于ωs=5 r/min和ωs=10 r/min。在晶體旋轉速度ωs=5 r/min時,強迫對流使得晶體下方熔體運動,從而實現熱量傳遞,使得相變界面附近熔體溫度升高。在晶體旋轉速度ωs=10 r/min時,熔體溫度升高更為明顯,固液界面相變溫度場溫度分布均勻性得到改善。

基于上述仿真結果,增加晶體旋轉速度可以增強晶體下方熔體流動,改善熔體內溫度分布和固液界面形狀,使熔體自由表面溫度上升,改善扭晶產生的過冷區現象。在提拉速度為0.8 mm/min時,分析不同晶體旋轉速度下熔體中的對流形態及溫度分布,如圖7所示。由圖3分析,當提拉速度為1.2 mm/min時,熔體自由表面出現了低于晶體熔點1 685 K的過冷區,該過冷區的溫度為1 680 K,與晶體熔點溫差為5 K,所以圖7顯示了等溫線間隔為5 K左右的溫度分布情況。觀察可知,隨著晶體旋轉速度的增加,固液界面形狀更加凸向晶體,熔體自由表面的過冷區有所改善。

由圖7可知,晶體旋轉能夠控制固液界面的形狀,當晶體旋轉速度為8 r/min時,如圖7(a)所示,與圖3(c)中晶體旋轉速度為4 r/min進行對比,固液界面形狀平坦,是因為晶體旋轉速度較低時,加熱器從坩堝外部加熱,固液界面上的溫度梯度呈現出中間低,四周高的溫度分布。當晶體旋轉速度增加,晶體下方傳熱加快,使得固液界面下方的溫度梯度增加,彌補了中間區域較低的溫度梯度。當晶體旋轉速度為12 r/min時,固液界面形狀過于凸向晶體,容易產生較大的熱應力,不利于晶體生長,如圖7(b)所示。因此,在提拉速度過大時,適當增加晶體旋轉速度可以使得固液界面穩定生長,避免熔體自由表面產生過冷區和固液界面形狀發生變化導致晶體扭晶,有利于硅單晶生長。

按照上述理論分析和有限元數值模擬結果,在晶體生長國家及地方聯合工程中心進行了12英寸硅單晶生長實驗,圖8(a)所示為12英寸直拉硅單晶生長爐實驗平臺。為了驗證提拉速度、晶體旋轉速度對直拉硅單晶生長的影響,設置其余工藝參數相同,在不同的晶體長度、提拉速度和晶體旋轉速度下,比較了6組晶體生長。其中,L是晶體長度,V是提拉速度,ωs是晶體旋轉速度。實驗中設置的工藝參數如表2所示。

表2 拉晶實驗參數Table 2 Experimental parameters of pulling crystal

通過上述方法,獲得拉晶實驗中直徑均勻的晶體,如圖8(b)所示。結果表明,較大的提拉速度下適當的提高晶體旋轉速度可維持晶體穩定生長,獲得直徑均勻的硅單晶。

4 結 論

本文針對直拉硅單晶實驗中提拉速度過大時觀察到的扭晶現象,通過有限元數值模擬得到熔體的溫度分布和熱分布,并結合晶體生長理論分析了扭晶的成因,提出一種最大提拉速度估計的方法,并且通過實驗驗證了該方法的準確性。過高的提拉速度會影響熔體內溫度分布的均勻性,主要表現為熔體自由表面出現過冷區。隨著提拉速度的增加,熔體自由表面過冷區面積逐漸擴大,導致晶體生長直徑不均勻產生扭晶現象,極大地降低了晶體生長品質。基于理論分析以及仿真和多次拉晶實驗驗證,提出了提高晶體旋轉速度的工藝措施,在維持較高的提拉速度的同時能夠有效避免扭晶現象,從而獲得直徑均勻的硅單晶,提高晶體生長效率。這一方法的提出對實際晶體生長工藝參數設定具有一定的指導作用。

猜你喜歡
界面
聲波在海底界面反射系數仿真計算分析
微重力下兩相控溫型儲液器內氣液界面仿真分析
國企黨委前置研究的“四個界面”
當代陜西(2020年13期)2020-08-24 08:22:02
基于FANUC PICTURE的虛擬軸坐標顯示界面開發方法研究
西門子Easy Screen對倒棱機床界面二次開發
空間界面
金秋(2017年4期)2017-06-07 08:22:16
鐵電隧道結界面效應與界面調控
電子顯微打開材料界面世界之門
人機交互界面發展趨勢研究
手機界面中圖形符號的發展趨向
新聞傳播(2015年11期)2015-07-18 11:15:04
主站蜘蛛池模板: 免费一级毛片在线观看| 青青草国产精品久久久久| 国产成人免费高清AⅤ| 国产成人午夜福利免费无码r| 五月天香蕉视频国产亚| 国产不卡一级毛片视频| 亚洲色图狠狠干| 91在线无码精品秘九色APP| 精品少妇人妻一区二区| 精品一区二区三区自慰喷水| 亚州AV秘 一区二区三区| 中文字幕无码电影| 国产成人久久综合777777麻豆| 国产国语一级毛片在线视频| 国产精品3p视频| 久久五月天国产自| 亚洲欧洲自拍拍偷午夜色| 伊人五月丁香综合AⅤ| 亚洲欧美不卡| 2021精品国产自在现线看| 国产午夜无码片在线观看网站 | 热伊人99re久久精品最新地| 国产精品思思热在线| A级毛片高清免费视频就| 国产屁屁影院| 最新国产午夜精品视频成人| 欧美激情视频一区| 国产情侣一区| 亚洲中文字幕无码mv| 精品福利网| 国模私拍一区二区三区| 欧美一级高清片欧美国产欧美| 国语少妇高潮| 日韩成人在线网站| 亚洲国产精品一区二区第一页免| 成人福利免费在线观看| 国产丰满大乳无码免费播放| 亚洲欧美另类日本| 夜色爽爽影院18禁妓女影院| 亚洲成人精品在线| 粗大猛烈进出高潮视频无码| 久久精品一品道久久精品| 亚洲人成网站在线播放2019| 青青草原国产| 亚洲天堂自拍| 日韩午夜伦| 久久鸭综合久久国产| 欧美一区二区精品久久久| 亚洲swag精品自拍一区| 特级做a爰片毛片免费69| 亚洲成av人无码综合在线观看| 国产福利不卡视频| 亚洲区视频在线观看| 成人va亚洲va欧美天堂| 亚洲自拍另类| 国产成人精品2021欧美日韩| 婷婷六月在线| 精品国产自在在线在线观看| 刘亦菲一区二区在线观看| 欧美激情福利| 中文字幕亚洲乱码熟女1区2区| 4虎影视国产在线观看精品| 国产1区2区在线观看| 亚洲国产精品日韩av专区| 中日韩欧亚无码视频| 亚洲国产欧美自拍| 欧美天堂在线| 国产女人水多毛片18| 国产毛片片精品天天看视频| 看看一级毛片| 国产成人午夜福利免费无码r| 国产精品入口麻豆| 国产成人午夜福利免费无码r| 无码日韩精品91超碰| 激情综合激情| 美女无遮挡被啪啪到高潮免费| 国产香蕉国产精品偷在线观看 | 五月激情婷婷综合| 亚洲人成网线在线播放va| 红杏AV在线无码| 91精品啪在线观看国产| 久久综合伊人77777|