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

應力作用下NiTi 形狀記憶合金微結構演化的相場模擬及其本征應變率敏感性*

2022-10-10 02:01:02席尚賓
爆炸與沖擊 2022年9期
關鍵詞:記憶合金

席尚賓,蘇 煜

(北京理工大學宇航學院,北京 100081)

NiTi 形狀記憶合金因具有優良的形狀記憶效應和超彈性變形能力而被廣泛應用于航空航天、微機電系統和生物醫療領域,其優良的熱力學性能主要源于不同熱-機環境下高溫奧氏體相與低溫馬氏體相之間的可逆相變。當環境溫度高于奧氏體完成溫度()且低于應力誘導馬氏體相變的最高溫度()(即<<)時,外載荷作用下材料表現出回復大變形的功能以及宏觀應力-應變曲線形成大的遲滯回環。NiTi 合金的超彈性力學行為具有明顯的應變率敏感性,并常被應用于高應變率環境(如NiTi 合金軸承),因此對其在動態加載(應變率為0.1~100 s)下的超彈性變形能力以及變形微觀物理機制的研究顯得十分必要。

NiTi 形狀記憶合金的應變率敏感性極大地影響著其力學性能,促使研究者開始關注這一問題。Chen 等采用分離式霍普金森壓桿脈沖成形技術,測試了NiTi 合金動態壓縮過程中的應變率效應,發現在動態載荷作用下材料在卸載完成后會存在殘余變形,應力-應變曲線形成的滯回環是開放的,而且動態加載過程中剛性到柔性轉變的開始應力遠高于準靜態加載過程的。Dayananda 等研究了應變率對NiTi 線超彈性行為的影響,發現在動態加載過程中,隨著應變率的升高,馬氏體相變的開始應力和完成應力以及應力-應變曲線的斜率均增大,尤其是相變階段的應力-應變響應變化更明顯。Nemat-Nasser 等研究了準靜態、動態以及高應變率下NiTi 合金的宏觀力學響應,發現準靜態和動態載荷作用下,材料主要表現為初始奧氏體變形機制和應力誘導的馬氏體形成機制,而在高應變率(大于1 000 s)加載過程中,馬氏體界面的運動則以界面阻尼效應為主,從而加工硬化速率隨著應變速率的升高而迅速升高。最近Jiang 等提出了一種熱力耦合的本構模型,數值計算了NiTi 合金超彈性變形的應變率敏感性,認為隨著應變率的升高材料相變響應由水平型向硬化型轉變,且隨著晶粒尺寸的減小,材料超彈性變形的應變率敏感性減弱。

以上實驗研究較全面地探索了NiTi 形狀記憶合金在不同加載環境下超彈性變形的應變率敏感性,并主要集中于高應變率下的宏觀應力-應變響應。在高應變率加載條件下,由于材料內部絕熱升溫造成的影響不容忽視,因而難以探究材料力學行為的本征應變率敏感性。而在中低應變率下,由于材料內部的溫升可以忽略,因而適于開展對材料力學行為本征應變率敏感性的研究。目前,對于NiTi 合金力學行為應變率敏感性微觀機理的研究較少,對于超彈性變形過程中奧氏體B2 相到馬氏體B19'相微結構動力學演化過程的實驗觀測方法也較有限。

為此,本文中,基于Ginzburg-Landau 動力學控制方程,建立NiTi 形狀記憶合金非等溫相場模型,實現對NiTi 合金內應力誘導馬氏體相變的計算模擬;同時將晶界能密度引入系統局部自由能密度,從而考慮多晶系統中晶界的重要作用;數值計算單晶和多晶NiTi 形狀記憶合金在單軸機械載荷作用下微結構的動態演化過程和宏觀力學行為,揭示該類材料力學行為的本征應變率敏感性的微觀機制。

1 相場模型及本構關系

建立NiTi 形狀記憶合金外載荷作用下微結構演化行為的相場理論模型,通過引入一組時空連續的標量場 φ (,) 來表征材料晶體幾何結構特征和物理性質。通常 φ 被稱為序參量或場變量,是關于空間位置和時間的函數。為了描述材料在熱-機環境中微結構的動態演化和力學響應,用序參量 φ(0≤φ≤1) 來表示NiTi 形狀記憶合金微觀組織結構的組成物(奧氏體、馬氏體等),構造恰當的自由能來描述馬氏體相變過程中微結構的演化。在納米晶NiTi 形狀記憶合金數值計算系統中,總的自由能可以表示為系統自由能密度 ψ 的體積積分。系統的自由能密度主要由3 部分組成,即局部自由能密度 ψ、梯度能密度 ψ和彈性能 密 度 ψ,具體表 達 式 為:

1.1 局部自由能

NiTi 形狀記憶合金中發生馬氏體相變的熱力學驅動力來自系統自由能的降低,即奧氏體相與馬氏體相之間的自由能之差為馬氏體相變提供了驅動力。多晶材料中存在相變惰性區域,因此設計NiTi 多晶形狀記憶合金的局部自由能由相變活躍區(晶粒內部)能量密度 ψ和相變惰性區(晶界區域)能量密度 ψ組成,即:

系統中相變活性區的局部自由能通常被相關的熱力學或物理變量來表示,如過冷度 Δ和相變潛熱。參考Landau 型局部自由能密度表達形式,采用二-三-四-四階多項式來表征NiTi 形狀記憶合金相變活性區局部自由能:

式中:、、和為溫度相關的Landau 系數。當序參量的值為0,即 φ=0 (=1, 2, ···, 6) 時,材料的微結構為奧氏體相;當序參量的值為1,即 φ=1(=1, 2, ···, 6) 時,材料的微結構為馬氏體相,φ表示系統中出現了第個馬氏體變體。

為了使得相場模型在不同的條件(溫度或應力)下形成穩定或者亞穩定的奧氏體相或馬氏體相,Landau 系數、、和的取值并非任意的,需要滿足以下條件:

此方程中,能量勢壘的構造主要考慮3 個方面:(1)在等效溫度以下(即<),當 φ =1 時,必須使局部自由能函數取得極小值,即在低溫下必須確保馬氏體相能夠穩定存在;(2)等效溫度以下,為了使僅溫度作用下的馬氏體相變能夠發生,奧氏體相必須處于亞穩態;(3)為了描述NiTi 形狀記憶合金在擬實熱-機環境中能量的連續演化,低溫下的能量勢壘依然是溫度相關的函數。由此,相變活性區局部自由能函數被建立。

多晶材料中晶界對材料力學性能的影響明顯,因此為了更準確地計算材料熱力學環境下的微觀相變行為,需要對晶界的屬性加以區分。通過對晶界屬性的深入了解,方程(2)中的晶界能密度項ψ可以表示為:

式中:為額外晶界能量勢壘系數,在晶粒內部被設置為零;為典型彈性應變能,與材料的剪切模量 μ 和典型無應力轉變應變 ε相 關,=με。晶界處的能量勢壘更高,因而不易發生相變現象。為此,采用與文獻[30]類似的關于序參量的二次函數作為晶界能密度的數學表達式。在選取恰當的能量系數后,=5,這個二次函數可以確保晶界能密度始終高于晶粒內部的自由能密度,如圖1 所示。

圖1 NiTi 形狀記憶合金局部自由能密度與序參量的關系Fig. 1 Relationship between local free energy density and order parameter of NiTi shape memory alloy

1.2 梯度能密度

為了描述材料微結構演化過程中相或變體之間的非尖銳界面,場變量相關的梯度能密度 ψ被引入到相場模型中,通常采用如下表達形式:

1.3 應變能密度及應力-應變本構關系

應變能是材料相變過程中儲存的變形能,是微觀結構形貌和熱力學行為的重要驅動力,其產生的物理機制是納米材料晶體點陣在相變過程中引起的晶格失配。應變能密度可以表示為:

式(13)中由不同馬氏體變體的相變應變構成了應力-應變關系中的非線性部分,在機械載荷作用下,由應力誘導的馬氏體變體逐漸形成并進行空間演化,這導致了NiTi 合金在宏觀上的超彈性應力-應變回線。馬氏體變體在空間上的演化速度與外載速率的競爭關系將導致材料力學行為的本征應變率依賴性,這一應變率相關性是由材料自身的晶體微結構相變特性決定的,并非是由絕熱升溫導致的,因此稱為本征應變率相關性。這一相關性在中低應變率下可更好地加以體現。

此外,在多晶材料模擬過程中為了計算全局坐標系下的轉變應變,旋轉張量分量(θ) 被用來表示局部坐標和全局坐標之間的關系。根據坐標轉換規則,彈性應變張量和彈性張量可表示為多晶材料模擬過程中為了計算全局坐標系下的轉變應變。根據坐標轉換規則,彈性應變張量分量和彈性張量分量可分別表示為:

式中:上標G 和L 分別代表全局和局部坐標系下的量。旋轉張量分量(θ) 通常被表示為:

式中: θ 為多晶材料中每個晶粒的取向角。

1.4 相場動力學方程

相場數值分析描述了模擬系統中自由能的動態最小化過程,場變量的連續演化由時間相關的Ginzburg-Landau 方程控制,該方程是基于場變量變化率與熱力學驅動力成正比假設的隨機相場動力學方程,即:

采用有限元方法對式(18)進行數值求解,將Ginzburg-Landau 動力學方程和力學平衡方程作為基本方程,數值求解獲得材料在力學邊界條件下的位移和相場信息。

1.5 數值計算參數

本研究中NiTi 形狀記憶合金的4 個臨界溫度分別設置為:馬氏體開始溫度=251.3 K,馬氏體完成溫度=213 K,奧氏體開始溫度=260.3 K,奧氏體完成溫度=268.5 K;相平衡溫度=(+)/2=255.8 K。

2 數值計算及結果討論

采用建立的相場動力學模型來研究NiTi 形狀記憶合金單晶、納米多晶動態加載過程中微結構的演化和宏觀應力-應變響應。如圖2(a)所示,對數值計算模型的左邊界約束水平方向的位移,即u=0,左上角節點施加固定,右邊界施加隨時間線性變化的位移載荷,如圖2(b) 所示。最大拉伸應變為6%,即ε=6%,應變率 ε˙ =15 s,系統的環境溫度為290 K。

圖2 數值計算邊界條件及外加載荷歷程Fig. 2 Boundary conditions and applied load history of the simulations

2.1 單晶NiTi 合金超彈性行為

NiTi 形狀記憶合金的超彈性行為源于單晶水平上的馬氏體相變,對其應力誘導的微觀相變行為進行數值研究。建立了幾何尺寸為60 nm×60 nm 的二維(2D)單晶NiTi 形狀記憶合金有限元模型,共包含28 800個平面三角形單元。為了盡可能全面地捕獲NiTi 合金在動態加載過程中的相變行為,對模型的上下和左右邊界施加序參量 φ的周期性邊界條件。

圖3 為NiTi 單晶形狀記憶合金單軸加載-卸載過程的應力-應變響應,圖4 為對應于應力-應變曲線上點處的微結構形態。從圖3 可以看到,加載過程可近似分為3 個階段:初期近似線彈性階段()、中期應力平臺()階段和末期近似線彈性階段()。在階段,序參量逐漸偏移0,數值計算系統能量開始發生變化,加載到點時,奧氏體組織完全消失,奧氏體向馬氏體的轉變開始,點處微結構如圖4 中點所示。馬氏體相變主要發生在應力平臺()階段,在點時形成了馬氏體驅前體如圖4 中點所示。當加載到點時序參量 φ 接近1,形成穩定的馬氏體,如圖4 中點所示,相變完成。加載末期(段),應力-應變曲線呈線性,材料的微結構不發生變化。加載完成后以相同的應變率進行卸載,卸載過程應力-應變響應出現應力遲滯現象,加載-卸載過程應力-應變曲線形成一個環。卸載完成后材料內部保留了0.11 %的殘余應變,材料未能完全恢復變形,原因是應變率過大,材料微結構來不及完全演化。但卸載完成后材料恢復到奧氏體結構,僅有極少數材料點處的序參量不為0,但非常接近于0(小于0.05)。卸載過程中材料發生逆馬氏體相變,序參量 φ 逐漸從1 轉向0,馬氏體結構逐漸褪去,隨后奧氏體結構開始形核并逐漸趨于穩定,微結構的演化如圖4 中點~所示。從加載-卸載過程微結構的演化(見圖4)可以觀察到,單晶水平上的馬氏體相變是均勻的,加載過程馬氏體的形核以及卸載過程中奧氏體的形核是在整個單晶區域內同時進行的。此處需要說明,本算例中采用了理想單晶假設,即晶體中不存在材料缺陷或熱擾動源,因而無法觀察到計算域內出現的馬氏體變體形核以及相界擴展現象。一般來說,在這種理想單晶假設下的相場計算會得到均勻分布的序參量,并且序參量會隨著加載增大而從0 到1 連續演化,類似的相場計算結果此前也有報道。這一點不同于隨后所考察的多晶問題,在多晶問題中,晶界在一定程度上起到了類似材料缺陷的作用,可誘發馬氏體變體在其附近形核生長。

圖3 單晶NiTi 合金290 K 環境溫度下單軸拉伸-卸載過程的應力-應變曲線Fig. 3 Stress-strain curve of monocrystal NiTi shape memory alloy during tension-unloading at 290 K

圖4 對應于圖3 中應力-應變曲線上關鍵點處的微結構,圖中不同顏色代表不同的微結構形態Fig. 4 Microstructures corresponding to the representative points of the stress-strain curve in Fig. 3,and different colors represent different microstructural morphologies

單晶NiTi 形狀記憶合金超彈性動態拉伸加載-卸載過程,應力-應變響應與準靜態加載-卸載過程相似,都出現了線性和應力平臺變化階段,且卸載過程出現應力滯后。但與準靜態加載過程不同的是,動態卸載完成過程中逆馬氏體相變未轉化完全而使得卸載完成后材料內部殘留了少量的應變。此外,動態加載過程中材料微觀相變更加激烈,馬氏體相變開始應力和結束應力都更大,應力-應變曲線具有較大的滯回環,即超彈性變形能力強。

2.2 納米晶NiTi 合金超彈性行為

為了研究納米多晶NiTi 形狀記憶合金在超彈性動態加載過程中微結構的演化和宏觀應力-應變響應,建立了如圖5 所示晶粒尺寸為60 nm 的NiTi 多晶幾何模型,每個晶粒的取向角 θ 在0°~45°隨機設置,晶界厚度設置為4 nm。對納米多晶動態加載過程中加載方式和應變率與2.1 節中對NiTi 單晶的加載過程完全相同,環境溫度設置為290 K。為了消除邊界對相變過程的影響,在NiTi 多晶幾何模型的4 個邊界上分別施加周期性邊界條件,如圖5 所示。

圖5 晶粒尺寸為60 nm 的多晶NiTi 形狀記憶合金有限元模型Fig. 5 A finite element model of the polycrystalline NiTi shape memory alloy with the grain size of 60 nm

圖6~7 為納米晶NiTi 形狀記憶合金動態加載過程的相場數值計算結果,圖6 為加載-卸載過程應力-應變響應,圖7 為對應于應力-應變曲線上關鍵點處的微結構形態。初始加載階段(圖6,段),材料應力-應變響應近似線性,數值計算模型中代表材料微結構的序參量 φ 開始偏離0,奧氏體結構逐漸消失。加載到點時,數值計算系統中部分材料點處序參量接近于1(圖7 中的點),應力誘導的馬氏體相變開始發生,馬氏體在容易產生局部應力集中的位置(晶界交匯處或施加載荷的邊界)形核,點對應的相變臨界應力 σ=796 MPa。這一臨界應力通常略高于實際形狀記憶合金中應力誘導馬氏體相變的臨界應力,因為實際材料內部往往存在有助于微結構形核的缺陷,如位錯、微觀孔洞等。段發生馬氏體相變,準靜態加載過程中點所對應的應力 σ(1.04 GPa)通常稱為相變完成應力,然而動態加載到點時,從應力-應變曲線上觀察到應力平臺階段結束,但從微結構的演化可以看到納米多晶系統中馬氏體相變還未完成(見圖7 中的點)。加載過程中應力平臺產生的原因是:當外加應力達到某一臨界值時,形狀記憶合金材料的微結構開始發生爆發式轉變,在幾十納秒內就可以從體心立方奧氏體結構(B2 相)轉變為單斜馬氏體結構(B19'相),在這極短時間內外加應力的變化很小,但相變過程產生了較大的局部內應變,從而導致宏觀應力-應變曲線表現出應力平臺。動態加載過程中加載速率與材料相變之間存在一定競爭關系,導致了應力平臺結束時材料中晶界附近還存在非馬氏體結構,如圖7 中點黑色虛線區域所示,直到加載完成(點)才完全轉變為馬氏體結構,如圖7 中點所示。

圖6 晶粒尺寸為60 nm 的NiTi 多晶形狀記憶合金290 K 環境溫度下單軸拉伸-卸載過程的應力-應變響應Fig. 6 Stress-strain response of the NiTi polycrystalline shape memory alloy with the grain size of 60 nm during tension-unloading at 290 K

卸載過程材料的應力-應變響應也表現出與加載過程相似的3 個階段。由圖7 中點~可以觀察到,當外加載荷卸載到達點時,材料晶界附近開始出現奧氏體結構,表明馬氏體到奧氏體的逆相變開始發生,馬氏體逐漸收縮成帶狀結構,奧氏體逐漸在晶界附近形核并向周圍擴展。因此,在納米晶NiTi 形狀記憶合金超彈性卸載過程中,逆馬氏體相變是通過晶粒內馬氏體結構的收縮-減少以及奧氏體結構的形核-擴展來實現的。然而,如圖7 中點所示,在卸載完成時刻,系統中有個別區域序參量 φ 的數值尚未完全達到0,即個別區域的奧氏體相變過程并未百分之百結束,導致此刻材料整體具有0.027%的殘余應變,這主要是由于卸載速率略高于微結構演化速度所致。此外,加載-卸載過程的應力-應變曲線形成的滯回環面積較小,相較于單晶材料來說,多晶NiTi 合金的吸能性能較弱。

圖7 對應于圖6 中應力-應變曲線上關鍵點處的微結構云圖Fig. 7 Microstructural morphologies corresponding to the representative points of the stress-strain curve in Fig.6

納米晶NiTi 形狀記憶合金在奧氏體完成溫度以上動態加載的相場計算表明,應力誘導馬氏體相變過程中材料微結構轉變速率與外部加載速率存在一定的相互作用。動態加載過程中NiTi 多晶系統內自由能迅速升高,驅動了立方奧氏體結構向單斜馬氏體結構轉變,但加載速率過高時,奧氏體-馬氏體界面前沿能量聚集,增大了馬氏體界面的擴展阻力,提高了相變開始應力。卸載過程的逆馬氏體相變也存在相同的微觀機制,因而卸載完成后系統中出現了殘余應變。

為了便于比較分析,將數值計算獲得的單晶和多晶NiTi 合金動態加載過程中應力-應變響應和Elibol 等的NiTi 形狀記憶合金超彈性實驗研究結果顯示在圖8 中。計算設置與實驗條件保持一致:平均晶粒尺寸為45 μm,加載應變為5%,應變率為20 s。從圖8 中可以觀察到:(1)相同應變率下,單晶系統的正、逆馬氏體相變開始應力和完成應力均高于多晶系統,表明晶界對晶粒內部的約束作用降低了馬氏體形核的臨界應力;(2)單晶系統應力-應變曲線滯回環面積更大,表明相同條件下單晶NiTi 合金具有更優的吸能性能;(3)多晶NiTi 合金的動態加載過程的數值模擬結果與實驗結果的趨勢吻合良好,這也驗證了相場模型的有效性。

圖8 NiTi 單晶、多晶系統在20 s-1 應變率下的應力-應變響應Fig. 8 Stress-strain curves of NiTi monocrystalline and polycrystalline at the strain rate of 20 s-1

為了研究NiTi 合金超彈性變形的應變率效應,模擬了晶粒尺寸為60 nm 的多晶系統在5×10~10s應變率下加載-卸載過程的應力-應變響應(見圖9(a)),并與Ahadi 等的實驗結果(見圖9(b))進行對比。數值計算結果表明,在納米晶NiTi 多晶系統中,應力-應變響應與應變率之間存在明顯的相關性??梢杂^察到,隨著應變速率的升高,奧氏體相到馬氏體相的轉變應力升高。較高的應變率會使正向和逆向馬氏體相變的開始應力升高,從而導致應力滯后回環面積的變化。由圖10 點處微結構圖可以看到,以5×10s的應變率加載到點時多晶系統中形成明顯的帶狀馬氏體結構,并隨著加載的進行迅速向相鄰的晶粒擴展。當外加應變達到4.7%(點)時,NiTi 合金中的馬氏體相變完成,多晶系統中奧氏體結構完全消失了。然而,隨著應變率的增加,加載到4.7%應變時NiTi 多晶系統中奧氏體的含量逐漸增加,如圖10 中點~所示。在較低應變率下多晶NiTi 形狀記憶合金奧氏體到馬氏體的轉變是通過馬氏體帶的空間擴展來實現的。在相變初期相鄰晶粒內部產生馬氏體變體,馬氏體疇可穿過晶界在空間內形成均勻分布的馬氏體條帶。隨著載荷的增加馬氏體條帶向周圍擴展,最終完成多晶系統從奧氏體相到馬氏體相的轉變。而在較高應變率作用下,材料內馬氏體轉變速度相對于加載速度產生滯后,因而在相同應力水平下材料內的大部分晶粒中未發生大面積馬氏體相變,難以形成廣泛分布的馬氏體條帶狀結構。只有當載荷進一步提升后,在晶界交匯處由于應力集中導致馬氏體變體在附近形核并擴展,直至馬氏體相變完成。因此隨著應變率的提高NiTi 多晶系統的臨界相變應力增大,且馬氏體形核位置較分散,難以形成明顯的馬氏體條帶疇結構。數值計算結果與實驗結果在趨勢上吻合良好。

圖9 納米晶NiTi 多晶系統在不同應變率下的應力-應變曲線Fig. 9 Stress-strain curves of NiTi polycrystalline system at different strain rates

圖10 對應于圖9(a)中應力-應變曲線上關鍵點處的微結構形態Fig. 10 Microstructure morphologies corresponding to the representative points on the stress-strain curves in Fig. 9(a)

3 結 論

以Ginzburg-Landau 動力學控制方程為基礎,通過在系統局部自由能密度中引入晶界能項,建立了非等溫的納米晶NiTi 形狀記憶合金相場模型,研究了單晶和多晶NiTi 形狀記憶合金在機械外載作用下微結構的演化特征和宏觀應力-應變響應,獲得的主要結論如下。

(1)單晶NiTi 形狀記憶合金系統內應力誘導的馬氏體相變是通過整體立方奧氏體結構向馬氏體結構的轉變而實現的,未形成奧氏體-馬氏體界面;而多晶系統中微結構演化通過局部馬氏體帶的形成-擴展來實現,通過馬氏體-奧氏體界面的向奧氏體結構區域的推移來完成相變。

(2)相同外載過程中,NiTi 單晶系統的宏觀應力-應變曲線具有更大的滯回環面積,擁有更優的超彈性變形能力;多晶系統卸載過程中奧氏體結構首先形成于晶界附近,并通過奧氏體-馬氏體界面向馬氏體疇結構區域的推移來完成逆馬氏體相變。

(3)在中低應變率下納米晶NiTi 形狀記憶合金應力-應變關系表現出較明顯的應變率相關性,隨著應變率的增高,正、逆馬氏體相變的開始和完成應力也隨之增大,卸載完成后材料內還保留一定的殘余應變。這一應變率相關性主要源于相場模型中外加載荷速率與馬氏體空間演化速度的相互競爭關系。

猜你喜歡
記憶合金
記憶合金在防護領域的應用
“記憶合金”科普儀器的改進與開發
化學教與學(2023年5期)2023-04-03 06:12:14
淺論SMA記憶合金
形狀記憶合金性能及結構加固應用綜述
廣東建材(2021年6期)2021-07-01 02:24:02
形狀記憶合金及其研究進展綜述
走進記憶合金
記憶合金內固定治療骨折臨床效果及存在問題研究
記憶合金骨卡環結合鋼板、髓內釘治療脛骨多段粉碎性骨折療效分析
西南軍醫(2015年3期)2015-04-23 07:28:21
基于形狀記憶合金的結構裂縫自修復研究
鎳鈦記憶合金環抱器內固定術后聯合中藥治療鎖骨骨折59例
中國藥業(2014年24期)2014-05-26 09:00:33
主站蜘蛛池模板: 国产高清免费午夜在线视频| 国产va在线| 成年女人a毛片免费视频| 无码国产伊人| 操国产美女| 欧美国产精品不卡在线观看| 国产欧美在线观看一区| 一级成人欧美一区在线观看| 精品无码国产一区二区三区AV| 欧美午夜在线观看| 无码啪啪精品天堂浪潮av| 国产精品视频第一专区| 国产高清在线丝袜精品一区| 久久久久久午夜精品| 老司机久久99久久精品播放| 日本免费高清一区| 91精品专区| 黄色不卡视频| 91精品亚洲| 国产精品亚洲专区一区| 1024国产在线| 国产精品人人做人人爽人人添| 九色在线观看视频| 亚洲一区二区三区麻豆| 成人午夜天| 五月婷婷欧美| av色爱 天堂网| 91成人在线观看视频| 国产制服丝袜91在线| 少妇露出福利视频| 91久久精品日日躁夜夜躁欧美| 亚洲精品视频网| 日本少妇又色又爽又高潮| 欧美国产在线看| 欧美影院久久| www中文字幕在线观看| 伊人久久久久久久久久| 99国产精品国产| 亚洲人成色在线观看| 国产主播一区二区三区| 国产91小视频| 日韩乱码免费一区二区三区| 成年av福利永久免费观看| 自拍亚洲欧美精品| 精品国产自| 国产哺乳奶水91在线播放| 国产黑丝视频在线观看| 国产在线观看人成激情视频| 国产精品网拍在线| 亚洲欧美成人在线视频| 精品视频在线观看你懂的一区| 国产一级毛片yw| 亚洲欧美另类专区| 欧美性精品| 国产精品视频白浆免费视频| 精品国产一区91在线| 免费国产黄线在线观看| 99偷拍视频精品一区二区| 在线毛片网站| 亚洲精品欧美日本中文字幕| 无遮挡一级毛片呦女视频| 日韩国产黄色网站| 亚洲 欧美 日韩综合一区| 国产精品网址在线观看你懂的| 大香网伊人久久综合网2020| 日韩欧美国产中文| 试看120秒男女啪啪免费| 色天天综合| 国产成人乱无码视频| 欧美综合在线观看| 精品色综合| 日韩av在线直播| 亚洲精品国产综合99| 日韩无码白| 日韩123欧美字幕| 亚洲人成在线精品| 日韩精品成人网页视频在线| 超碰aⅴ人人做人人爽欧美| 欧美啪啪网| 亚洲区第一页| 国产高清国内精品福利| 久久免费看片|