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

The Analysis of Thermal-Induced Phase Transformation and Microstructural Evolution in Ni-Ti Based Shape Memory Alloys By Molecular Dynamics

2019-08-13 06:16:40HsinYuChenandNienTiTsou

Hsin-Yu Chen and Nien-Ti Tsou

1Department of Materials Science and Engineering,National Chiao Tung University,Ta Hsueh Road,Hsinchu 30013,Taiwan.

Abstract:Shape memory alloys has been widely applied on actuators and medical devices.The transformation temperature and microstructural evolution play two crucial factors and dominate the behavior of shape memory alloys.In order to understand the influence of the composition of the Ni-Ti alloys on the two factors,molecular dynamics was adopted to simulate the temperature-induced phase transformation.The results were post-processed by the martensite variant identification method.The method allows to reveal the detailed microstructural evolution of variants/phases in each case of the composition of Ni-Ti.Many features were found and having good agreement with those reported in the literature,such as the well-known Rank-2 herringbone structures;the X-interface;Ni-rich alloys have lower transformation temperature than Ti-rich alloys.In addition,some new features were also discovered.For example,the Ti-rich alloys enabled an easier martensitic transformation;the nucleated martensite pattern determined the microstructural evolution path,which also changed the atomic volume and temperature curves.The results generated in the current study are expected to provide the design guidelines for the applications of shape memory alloys.

Keywords:Thermal-induced phase transformation,transformation temperature,molecular dynamics,crystal variants.

1 Introduction

Ni-Ti shape memory alloys (SMAs) have gained great attention recently because of its unique shape memory effect (SME) and superelasticity.They have been widely used in aerospace,sensors,actuators,and biomedical devices [Manfredi,Huan and Cuschieri (2017);Takaoka,Horikawa,Kobayashi et al.(2002)].The behavior of shape memory effect and the material properties are dominated by the transformation temperature and the evolution of microstructures in the crystal.There are four characteristic phase transformation temperatures:martensite transformation start temperatureMs,martensite transformation finish temperatureMf,austenite transformation start temperatureAsand austenite transformation finish temperatureAf.Ni-Ti SMAs adopt the austenite phase at high temperature and transform into the martensite phases when cooling belowMs[Shaw and Kyriakides (1995)].Crystals with the martensite phases typically adopt the twinned structures,which minimize the total energy and result in the anisotropic behavior of the crystals.Therefore,transformation temperature and the phase transition process influence the material properties significantly.

Mehrabi et al.[Mehrabi,Kadkhodaei and Elahinia (2014)] pointed out that the twinned and detwinned structures play an important role in the level of the shape recovery.Chen et al.[Chen,Liu,Li et al.(2018)] showed that the SMAs exhibit different superelastic and shape memory effect under different temperatures.The composition of materials also greatly influences the behavior of SMAs.Carl et al.[Carl,Smith,Van Doren et al.(2017)] found that with the higher concentration of nickel (Ni-rich) the transformation temperature decreases dramatically;however,with the lower concentration of nickel (Tirich),the transformation temperature is similar with it of the equiatomic Ni-Ti[Razorenov,Garkushin,Kanel et al.(2009)].Many researchers have made great contribution on the thermal-mechanical properties of SMAs,however,the corresponding microstructural evolution and the detailed phase transition are still subtle.Thus,in the current work,molecular dynamics (MD)simulations were performed to predict the material properties,such as transformation temperatures,and also to generate the detailed patterns of the austenite/martensite phases.

Potential used in the MD simulations allow to generate accurate atomistic unit cell distortion,and thus,the martensitic phase transformation of Ni-Ti alloys can be predicted.Two types of potential have been widely used for the simulation of SMAs.Lai et al.[Lai and Liu(2000)]proposed the embedded-atom method Finnis-Sinclair(EAM-FS)potential to describe the atomic interaction of Ni-Ti alloy.Mutter et al.[Mutter and Nielaba(2011)] improved the FS potential to examine the austenite-martensite interfaces in Ni-Ti nanoparticles.Then,Ko et al.[Ko,Grabowski and Neugebauer (2015)] optimized the Ni-Ti potential based on the second nearest neighbor modified embedded-atom method(2NN-MEAM),and provided the accurate temperature-induced phase transformation and the stress distribution in Ni-Ti crystal.Srinivasan et al.[Srinivasan,Nicola and Simone(2017)] compared the 2NN-MEAM potential with the EAM-FS potential in detail and found that the 2NN-MEAM potential could predict martensite phase transformation more accurately.As a result,2NN-MEAM potential has become a well-accepted potential form for describing the thermal-mechanical behaviors of Ni-Ti binary alloys.

Most of the MD simulation results are analyzed by the conventional post-processing methods,such as the common neighbor analysis(CNA)[Stukowski(2012)]or polyhedral template matching(PTM)[Larsen,Schmidt and Schi?tz(2016)].However,those methods identify the lattice structures,such as HCP and BCC,which provide the information of the interface between the austenite and martensite phases only.In order to understand more detailed information about the variants and twin microstructure of SMAs,Yang et al.[Yang and Tsou(2017)]developed a post-processing method for MD simulations,referred as martensitic variant identification method (MVIM) in the current work,that identify multiple crystal variants and phases by examining the transformation of each lattice.Lu et al.[Lu,Chen and Tsou(2019)]adopted this method to study the superelasticity of Ni-Ti SMAs and corresponding microstructural evolution.The detailed microstructure,such as transitional orthorhombic phases between the austenite and monoclinic phases,twin patterns,and volume fraction can then be revealed.

In the current work,the standard MD package (LAMMPS) were used to stimulate the thermal-induced process and the corresponding microstructural evolution in the Ni-Ti alloys with the composition of Ti-rich,Ni-rich and the equiatomic Ni-Ti alloys were revealed by the post-processing method proposed by Yang et al.[Yang and Tsou(2017)].Many features were predicted,having good agreement with those reported in the literature.The relationship between the microstructural evolution and the atomic volume and temperature curves of the three groups is discussed in the following sections.

2 Methodology

The simulation was performed by a well-known MD package,LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) [Plimpton (1995)].The 2NN MEAM potential [Lee,Baskes,Kim et al.(2001)] was adopted to describe the thermal-induced phase transition behavior of Ni-Ti alloys.The lattice constanta0of Ni-Ti was set as 2.999 ?[Ko,Grabowski and Neugebauer(2015)],the cutoff radiusrcwas set as 5.0 ?.Periodic boundary conditions(PBC)were applied in the three directions,x,y,andzin all the cases.Atoms Ni and Ti were placed at a relative position at(0 0 0)and(0.5 0.5 0.5),respectively in the Cartesian coordinate([100],[010],[001]).The size effect of the simulation box was studied.Five different size of boxes with the length of 10 nm,15 nm,20 nm,25 nm,and 30 nm were generated,in order to obtain the converged results with efficiency.Then,different compositions of Ni and Ti were considered here.They were 49%and 49.5%of nickel for modeling Ti-rich alloys,50.5% and 51% of nickel for modeling Ni-rich alloys,and 50%for modeling the equiatomic Ni-Ti alloys.The simulation box was a 20.1×20.1×20.1×nm3cube as shown in Fig.1 which was filled with 67 unit cells along each side,giving 601,526 atoms in total in the box.The atomistic model with different composition of Ni-Ti was created by randomly replacing the atoms to match the given concentration.

2.1 Model setting

The model was initially relaxed at 0 K by the conjugate gradient method in order to minimize the total energy and achieve the dynamic balance.The time step was set as 2 fs,and the temperature was controlled by the isothermal-isobaric(NPT)canonical ensemble along with a pressure of 1.013 bar.Then the model underwent the thermal equilibrium at 500 K for 20 ps (10,000 time steps),where the velocities were assigned to all the atoms by the Gaussian distribution based on the given temperature.Noted that the temperature 500 K,which is above the austenite finish temperature (Af),was chosen here to obtain a stable initial austenite state.After the thermal equilibrium,the cooling and heating processes were performed at the rate 1 K/ps.Where the cooling process started from 500 K to 100 K and the heating process started from 100 K to 500 K.The thermal-induced microstructural evolution between austenite to martensite phases can be observed.

Figure 1:The molecular-dynamics model containing 601,526 atoms

2.2 Method of microstructure analysis

The MD results were post-processed by the scientific visualization software OVITO,which provides powerful visualization algorithms [Stukowski (2009)].For example,polyhedral template matching (PTM) [Larsen,Schmidt and Schi?tz (2016)] and dislocation analysis(DXA),were applied to reveal the martensitic transformation and dislocation evolution in the SMAs.Although these conventional algorithms provide the information about the phase transition and lattice structure,the information of the coexistence of multi-phases and multi-variants in SMAs was absent.Therefore,the martensite variant identification method(MVIM)proposed by Yang et al.[Yang and Tsou(2017)]was used to study more detailed information about microstructure in SMAs.Typically there are three crystal systems can be observed in Ni-Ti alloys,which are austenite(A),orthorhombic(O),and monoclinic(M).To be more specific,there are 6 crystal variants in the orthorhombic crystal system and 12 variants in the monoclinic system.According to the literature[Bhattacharya(2003)]each crystal variant has its theoretical transformation matrices UX,whereXis the labeled of the name of each crystal variant.For example,the theoretical UM5are illustrated as

whereα= 1.2043,γ= 0.9563,δ= 0.058,and?=-0.0427 in the case of Ni-Ti alloys [Knowles and Smith (1981)].Note that the crystal numbering used in the current work was based on it in the literature[Bhattacharya (2003)].There are 15 conditions for describing the relationship between components in the matrices [Yang and Tsou (2017);Lu,Chen and Tsou (2019)].For example,in the case of a theoretical variantM5,its transformation matrix is required to satisfy the conditions thatU11=U33,U12=U23,and so on so forth based on the arrangement of the transformation matrix shown in Eq.(1).Now consider one deformed lattice generated by the current MD model as an example.Its transformation matrices can be determined,say

It is obvious that the determined matrix matched with it ofM5 shown in Eq.(1).The similarity of determined matrix can be checked by comparing with the components in the theoretical transformation matrices.The lattice will be classified as the crystal variant with the theoretical transformation matrix with the highest similarity.The detailed settings,such as tolerance and threshold,are mentioned in the literature[Yang and Tsou(2017);Lu,Chen and Tsou(2019)].In the cases of the similarity lower than 50%or more than one crystal variant is identified,the lattice will be marked as"Other".By using the martensite variants identification method,detailed arrangement of the phases/crystal variants in the crystal can be identified.

3 Result and discussion

3.1 Thermal-induced phase transformation in the equiatomic Ni-Ti

Fig.2 shows that the atomic volume and temperature curves of the equiatomic Ni-Ti predicted by the MD model with the length of the simulation box of 10 nm,15 nm,20 nm,25 nm,and 30 nm.The temperature was cooling from 500 K to 100 K and then heating back to 500 K.The results show that the atomic volume changed dramatically during the phase transition,indicating the transformation temperatures in the figure.Many studies [Ko,Grabowski and Neugebauer (2015);Chen,Liu,Li et al.(2018);Prokoshkin,Khmelevskaya,Korotitskij et al.(2003)]pointed out that the sharp increased of the atomic volume predicted the austenite transforming into martensite.Where the atomic volume was calculated by the average volume occupied by a single atom.It can be observed that the curves converged,especially the value ofAf,when the length of the simulation box was greater than 20 nm.Thus,the box length of 20 nm is adopted in the current work for better efficiency.

Next consider the microstructural evolution corresponding to cooling and heating process.Fig.3 shows that the cooling process of the equiatomic Ni-Ti from 500 K to 100 K.The crystal adopted the austenite phase (A) initially (Fig.3(a)) and transformed into the martensite phases at around 220 K (Ms),leading to a dramatic increase of the atomic volume.There were three stages of transformation process during the cooling test in the equiatomic Ni-Ti.At the first stage as shown in Figs.3(a)-3(b),some regions in the crystal underwent the phase transformation fromAtoO3 andO4.Then,the monoclinic twin grew by expending theOvariants.There were two types of monoclinic twin in the crystal:M5-M6 twin andM7-M8 twin,as shown in Figs.3(b)-3(c).It can be observed that the atomic volume increased slower than the other two stages.At the third stage (Figs.3(c)-3(d)),the atomic volume increased sharply again.It is worth to note thatO3 transformed intoM5-M6 twin,whileO4 transformed intoM7-M8 twin,respectively.As most of the austenite transformed into the martensite phases,i.e.,temperatureMf ~200 K,a Rank-2 herringbone structure,consisting of the four monoclinic variants,was adopted by the crystal,as shown in Fig.3(e).

Figure 2:The hysteresis curves of the thermal-induced phase transformation in the equiatomic Ni-Ti,where the length of the simulation boxes were 10,15,20,25,and 30 nm

Figure 3:The evolution of the microstructure in the Ni-Ti equiatomic under the cooling process

The crystal was then heated back to 500 K.The atomic volume-temperature curve and the corresponding microstructural evolution are shown in Fig.4.The crystal remained the Rank-2 herringbone structure,until theAs ~410K was reached (Fig.4(a)).Then,the atomic volume curve decreased sharply.It is worth noting that the phase transformation path between the austenite and martensite phases was different with it during the cooling process.Here,the monoclinic twins phases firstly transformed into theOphase,and then,the region of theOphase further transformed into the austenite phase.More specifically,M5-M6 twins transformed intoO3 variants,whileM7-M8 twins transformed intoO4 variants (Fig.4(b)).There was no direct interfaces between the austenite andMphases observed during the heating process,as shown in Fig.4(c).Finally,the microstructure returned to the pure austenite phase when the temperature was aboveAf ~440 K.

Figure 4:The evolution of the microstructure in the Ni-Ti equiatomic under the heating process

In order to reveal the detail of the microstructure at the beginning and the end of the phase transformation between the martensite and austenite phases,two microstructures at 420 K and 430 K during the heating process were chosen here to be analyzed by PTM,DXA,and MVIM as shown in Fig.5.It is also worth to verify the resulting patterns and to demonstrate the power of the MVIM by comparing with those generated by the two commercial approaches.At 420 K,PTM revealed a HCP band structure separated by the interfaces,which were mainly composed with the BCC structure.These interfaces between the bands of HCP accumulated dislocations as shown in Fig.5(b).Fig.5(c)shows the result generated by MVIM,revealing a Rank-2 herringbone structure,i.e.,twins within twins,consisting of the four monoclinic variants,which is consistent with the finding in the experiment [Petersmann,Pranger,Waitz et al.(2017)].Based on the literature [Bhattacharya (2003)],both the orthorhombic (Ophases) and monoclinic (Mphases) were with the HCP lattice structure,while the austenite phase adopted the BCC structure.By comparing with the microstructure shown in Fig.5(a),it can be observed that the region of theOandMvariants identified by MVIM matched it of the HCP identified by PTM;the austenite phase at the interfaces were also accurately identified by MVIM.Although PTM can capture the martensitic transformation in SMAs,detailed variant information is absent.MVIM demonstrates its capability to reveal the coexistence of different phases and variants,providing an insight of the microstructural evolution in SMAs.

Figure 5:The microstructures at 420 K and 430 K post-processed by (a)(d) PTM,(b)(e)DXA (dislocation analysis) methods and (c)(f) MVIM (martensite variant identification method),respectively

The microstructure at 420 K contained interfaces with dislocations,and also a segment of the interface identified as the FCC structure,which referred as stacking faults in the crystal [Chen,Liu,Li et al.(2018)].These interfaces led to the high energy area where initialized the phase transformation.This can be seen by comparing Figs.5(a) and 5(d),where all the thin interfaces at 420 K expanded toward the directions which decrease the volume fraction of the HCP regions (monoclinic twins) at 430 K.This phenomenon was also reported in the literature[Ko,Maisel,Grabowski et al.(2017)].The monoclinic twin boundaries adopted two orientations,vertical and horizontal,as marked with red lines in Fig.5(f).They are also known as X-interfaces reported in the literature [Basinski and Christian (1954);Seiner (2015);Lu,Chen and Tsou (2019)].The X-interfaces allow the switch between the austenite and martensite phase along the low energy evolution path.Two types of monoclinic twins presented here,M5-M6 andM7-M8,which were surrounded by theO3 andO4 variants,respectively.It can be observed that the twoOphases served as the transition region between the monoclinic and austenite phases.Such transitional regions had no dislocation accumulated as shown in Fig.5(e).It is worth noting that this feature can be captured by the MVIM only,as theOandMphases are indistinguishable for the PTM method.

3.2 Thermal-induced phase transformation in ni-ti with different composition

Fig.6 shows the atomic volume-temperature curves in the cases of 49%,49.5%,50%,50.5%,and 51%-Ni in the Ni-Ti alloys.It can be observed that the overall atomic volume of the system increased with the higher concentration of Ni.The characteristic temperaturesMsandAfcan also be obtained by examining the onset and the end of thephase transformation,respectively,as shown in Tab.1.The results showed that both theMsandAfdecreased with the increase of the concentration of Ni.This feature has very good agreement with the experimental results in the literature [Razorenov,Garkushin,Kanel et al.(2009);Otsuka and Ren(2005);Frenzel,George,Dlouhy et al.(2010)].In order to observe the trend of transformation temperature,here we used the estimated valueT0which was proposed by Frenzel ea al,whereT0(XNi) = (Ms(XNi)+Af(XNi))/2 [Frenzel,George,Dlouhy et al.(2010)].Fig.7 shows the comparison result ofT0between the experimental data and our present work.The trend of the simulation results were fit with the experimental data.Furthermore,the current simulation results also show two more features described as follows.(1)The decreases ofMsdue to the increasing concentration of Ni were less significant compared with those ofAf.More specifically,the range ofMsin the five cases considered here was within 17 K only,while it ofAfwas 95 K.(2)The composition effect to alter theAfwas more significant in the cases of the Ni-rich alloys compared with the Ti-rich alloys.There was only a 2 K temperature increment as the composition of Ni decreased from 49.5%to 49.0%,while a 52 K temperature increment as the composition of Ni decreased from 51.0%to 50.5%.

Table 1:Transformation temperature of Ni-Ti shape memory alloys with different compositions

Figure 6:The atomic volume-temperature curves in the cases of 49.0%,49.5%,50.0%,50.5%,and 51.0%-Ni in the Ni-Ti alloys

The MD results of the five cases were post-processed by MVIM.It is of interest to note that the microstructures and their evolution paths were similar in all the cases described here.Only slight differences were observed,and yet,made such significant changes on the characteristic temperature as shown in Tab.1.Here we take the cases of 50.5%-Ni and 49.5%-Ni as examples for the illustration purposes.

Figure 7:The comparison result of T0 between experimental data[Frenzel,George,Dlouhy et al.(2010)]and our present work

Fig.8 shows the atomic volume-temperature curve and the corresponding microstructures in the case of 50.5%-Ni.Firstly,the crystal was cooled down from the austenite phase.After theMstemperature was reached,variantsO1 andO2 presented in the crystal randomly,as shown in Fig.8(a).They expanded further into a Rank-1,alternating band structure,and then,two types of twin,M1-M2 andM3-M4,occurred in the region ofO1 andO2,respectively,as shown in Fig.8(b).As the temperature kept decreasing,theOphases were replaced by the two types of twins.These twin were separated by oblique interfaces marked as red lines in Fig.8(c).It is of interest to note that the exchange of the volume fraction between monoclinic twins andOphases remained the overall atomic volume,giving a small plateau between stage(b)and(c).However,a dramatic change of atomic volume occurred when the twin interfaces rotated to a horizontal one(marked as the red line),forming a Rank-2 herringbone twin consisting the four monoclinic variants,as shown in Fig.8(d).This herringbone structure was stable during the remaining cooling and the heating process until the temperature reached around 390 K.The evolution path was very similar to it described in the case of 50.0%-Ni as shown in Figs.5(c) and 5(f).The phase transformation initiated from the interfaces and theOphase served as the transition region between the monoclinic twin and austenite phase,as shown in Figs.8(e) and 8(f).At stage (g),the crystal returned to the austenite phase,however,many unknown phases(identified asOther) remained.It can be observed that all the microstructures show in Fig.8 contained noise and much moreOtherregions compared with those in Figs.3 and 4.This is because there were many Ti atoms replaced by Ni atoms,resulting in local lattice distortion.

Figure 8:The evolution of the microstructure and the hysteresis curves of the thermalinduced phase transformation in the 50.5%-Ni

Figure 9:The evolution of the microstructure and the hysteresis curves of the thermalinduced phase transformation in the 49.5%-Ni

Fig.9 shows the atomic volume-temperature curve and the corresponding microstructures in the case of 49.5%-Ni.The Ti-rich crystal adopted a microstructural evolution path which was very similar with it of the Ni-rich crystal described in the previous paragraph.The major difference between the two cases can be found in Fig.9(a),where more volume fraction ofO1 andO2 presented,and most importantly,they seemed to be patterned rather than randomly occurred in the crystal.This indicates that the lattice distortion in the Ti-rich case enabled an easier phase transformation between the austenite toOphase.In addition,this patternedOstructure,which is different from the random structure shown in Fig.8(a),allow a direct phase transformation from variantO1 toM1-M2 twins and variantO2 toM3-M4 twins.This resulted in a gradual change between stage (b) to (c),and also less twin interface rotation was required,as shown in Figs.9(c) to 9(d).The microstructures corresponded to the remaining heating process were similar to those in the previous case.However,the noise andOtherregion were not that obvious compared with the case of 50.5%-Ni.

4 Conclusion

In the current work,the thermal-induced phase transformation and the corresponding microstructural evolution in a Ni-Ti single crystal was studied by MD simulation and various post-processing tools,PTM,DXA and MVIM to analyze the detailed microstructural evolution.The effect of the simulation box size was studied,and the length of 20 nm was chosen for efficiency and accuracy.The relationship between the characteristic temperature and the different composition of Ni-Ti alloys was also revealed.The feature that theMsandAftemperature decreased with the increase of the concentration of Ni was accurately predicted,having good agreement with the experimental results in the literature.Many features were also discovered,such as Rank-2 herringbone structure;the transitionalOphases region between the austenite andMphases;the dislocation accumulated at the twin boundaries where initiated the phase transformation;the X-interfaces etc.Moreover,the different nucleated martensite pattern lead to different microstructural evolution path.That also altered the atomic volume and temperature curves.The results presented here may provide fundamental understanding to the microstructure behavior in SMAs,and provide design guideline for their applications.

主站蜘蛛池模板: 欧洲高清无码在线| 国产91丝袜在线播放动漫 | 亚洲国产日韩欧美在线| 无码中文字幕乱码免费2| 囯产av无码片毛片一级| 成人国产精品一级毛片天堂| 91久久偷偷做嫩草影院| 99久久精品国产精品亚洲| 欧美日韩一区二区三区在线视频| 国产日产欧美精品| 久久毛片基地| 91久久偷偷做嫩草影院精品| 99热在线只有精品| 国产在线视频导航| 欧美色视频在线| 欧美色图第一页| 亚洲一区色| 国产亚洲现在一区二区中文| 老汉色老汉首页a亚洲| 亚洲成人一区二区三区| 青青久视频| 日本不卡在线播放| 福利在线免费视频| 国产精品久线在线观看| 国产91精品久久| 伊人大杳蕉中文无码| 国产性爱网站| 久久女人网| 又猛又黄又爽无遮挡的视频网站| 国产色婷婷| 国产欧美网站| 少妇极品熟妇人妻专区视频| 日本黄色a视频| 制服无码网站| 熟妇丰满人妻| 亚洲欧美日韩另类| 久久亚洲综合伊人| 精品免费在线视频| 在线欧美一区| 91久久性奴调教国产免费| 国产午夜小视频| 日韩在线成年视频人网站观看| 人妻91无码色偷偷色噜噜噜| julia中文字幕久久亚洲| 国产精品久久自在自2021| 亚洲第一成人在线| 国产人人射| 欧美色图第一页| 日本一区二区三区精品AⅤ| 亚洲V日韩V无码一区二区 | 色婷婷视频在线| 亚洲一级毛片免费观看| 亚洲第一香蕉视频| 91国内在线视频| 亚洲无码视频图片| 成人自拍视频在线观看| 色视频国产| 国产福利不卡视频| 国产一区二区三区在线精品专区| 男女男免费视频网站国产| 片在线无码观看| 欧美激情视频一区| 一区二区三区成人| 国产日韩欧美在线视频免费观看 | 亚洲天堂视频在线观看免费| 日本不卡视频在线| 国产乱子伦精品视频| 国产成人精品无码一区二| 亚洲日韩欧美在线观看| 日韩精品无码一级毛片免费| 青青青视频免费一区二区| 日韩欧美色综合| 亚洲天堂成人在线观看| 自慰网址在线观看| 中文字幕人成乱码熟女免费| 久久77777| 亚洲欧美自拍视频| 在线观看av永久| 亚洲日韩精品伊甸| 久久久久亚洲av成人网人人软件| 草逼视频国产| www欧美在线观看|