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

    Characterization of size effect of natural convection in melting process of phase change material in square cavity?

    2021-10-28 07:12:22ShiHaoCao曹世豪andHuiWang王輝
    Chinese Physics B 2021年10期
    關(guān)鍵詞:王輝

    Shi-Hao Cao(曹世豪) and Hui Wang(王輝)

    1College of Civil Engineering,Henan University of Technology,Zhengzhou 450001,China

    2School of Civil Engineering and Architecture,Hainan University,Haikou 570228,China

    Keywords: phase change material,natural convection,size effect,equivalent thermal conductivity

    1. Introduction

    Because of high latent heat capacity, phase change materials (PCMs) exhibit great advantages over traditional materials for thermal energy storage. For example, 1-kg solid paraffin with latent heat 210.5 kJ/kg can absorb 210.5-kJ thermal energy in the melting process and the melting temperature almost keeps constant. This preferable property makes PCMs considerably useful to store or release thermal energy in a wide range of applications such as buildings,solar power plants,electronic components,spacecrafts,etc.[1–4]

    Over the past decades,many studies have been carried out to analyze the heat transfer characteristic of PCMs. For example,the enthalpy method was first developed by Eyreset al.for solving heat conduction in solids[5]and then it was extended to the heat conduction controlled phase change problems,which are also named Stefan problems.[6–8]Recently, Sultanaet al.compared several fixed grid techniques for solving the Stefan problem,thereby predicting moving solid–liquid interface in the phase change process.[9]However, the classic Stefan problem is not suitable for the analysis of melting process of PCM, because the shape of melting front and the total melting time of PCM are actually determined by the heat conductive and convective processes together.[10]Natural convection caused by temperature gradient or density variation in liquid phase significantly changes the thermofluid transport mode.To effectively predict the thermal and hydrodynamic process of a material experiencing phase change, the natural convection effect must be included.[11–14]For conveniently describing the whirlpool flow caused by natural convection, Boussinesq approximation,[15]and porous flow theory[13]are generally utilized to introduce the driving force into the momentum equation.

    All those researches reported that the actual solid–liquid phase change process is extremely complex, because it involves with heat transfer, solid–liquid phase transition, liquid phase flow, solid phase motion, volume expansion, and moving solid–liquid interface.[16]Moreover, the melting process is strongly dependent on the geometry shape and size of PCM domain,[12,17,18]the heating direction,[19,20]the enhancing strategy,[21–27]and the constrained mode of the solid phase.[13,17]Although it has been demonstrated that the natural convection in the liquid phase plays a positive role in accelerating the melting process of PCM, the description of natural convection behavior in PCM melting process may encounter some obstacles: (i) the natural convection process is too complex to be captured accurately, (ii) the simulation of natural convection is time-consuming, (iii)the effect of natural convection is size-scale dependent,(iv)the effect of natural convection is related to the applied boundary conditions, and(v)the contribution of natural convection is difficult to quantitively assess. Therefore,comprehensively understanding the melting process and the role of natural convection presents a considerable challenge in applying PCMs to thermal energy storage. The study of influence of natural convection on melting behavior is still insufficient, and the physical principle in it is not specifically explained.

    In the present study,we aim to extend the understanding of the performance of natural convection melting process of PCM to the cases of different scales and develop an equivalent thermal conductivity concept for including the contribution of natural convection and reducing the computational difficulty as well. To do so, the melting process of paraffin wax within a rectangular container is studied. The melting experiment of paraffin wax in the rectangular box is conducted to trace the moving of solid–liquid interface and validate the computational model. Then,the temperature field and velocity field in the melting process are simulated and the influences of natural convection on the liquid fraction in various time steps and the total melting time are discussed. Also, with this computational model, the size effects of natural convection are analyzed in different size cases and the critical value of domain size which corresponds to the case that the effect of natural convection can be ignored below this value is obtained.Simultaneously,the physical principle behind it is revealed. Finally,the effect of natural convection on the melting time of PCM is represented by the introduced equivalent thermal conductivity concept,and an empirical correlation of equivalent thermal conductivity is given for simplifying the complex calculation of natural convection in engineering applications.

    The rest of this paper is organized as follows. In Section 2, the melting experiment of paraffin wax within a rectangular container is described. In Section 3, the basic governing equations to simulate the melting behavior of the PCM are presented. The numerical model and boundary conditions are explained in detail in Section 4. The simulation results and discussion are given in Section 5. Finally, the important findings are discussed in Section 6.

    2. Melting experiment

    To further investigate the thermal performance and the energy storage capability of the PCM, the paraffin melting experiment in a rectangular container is carried out,and the schematic diagram of the experiment is shown in Fig.1. The rectangular container with interior dimensions of 100 mm×100 mm×300 mm is made of six transparent perspex plates,each of which has a thickness of 10 mm and good thermal insulation performance,and allows visual observation of the melting inside the container. To further reduce heat loss in the melting process and provide nearly adiabatic environment, the rectangular container is covered with the cotton layer and the polyfoam layer as indicated in Fig.1.The insulation layer can be uncovered when we take photos by a digital camera to record the instantaneous locations of solid–liquid melting front during the melting experiment. The camera is positioned appropriately in front of the container. The cavity is filled with solid PCM paraffin wax(melting temperature 300 K) purchased from Shanghai Joule wax Company: such a wax is adopted because of its stability in melting and freezing cycles. The physical properties of paraffin wax are listed in Table 1. To measure the temperature change in the melting process of paraffin wax, five thermocouples are arranged at different locations (K1 ath=2 cm, K2 ath=4 cm, K3 ath=6 cm, K4 ath=8 cm, and K5 directly exposed to the environment). The temperature data are collected by a 24-channel temperature measurer(Anbai,AT4524),and stored in a computer. Besides,prior to the melting experiment,the rectangular container filled with paraffin should be exposed to the laboratory environment for at least 24 h to ensure a uniform initial temperature in the paraffin wax. The initial temperature of paraffin wax is(16±0.1)°C.

    Fig.1. Schematic diagram of paraffin melting experiment. The numbers in the figure describe the fallowing units: 1-polyfoam,2-cotton,3-rectangular container, 4-paraffin, 5-hot plate, 6-thermocouples, 7-data acquisition unit,8-computer,9-camera.

    Table 1. Material properties of paraffin.

    Fig.2. Evolution of experimental melting front in entire melting process(front view)at f =0(a),0.02(b),0.06(c),0.14(d),0.24(e),0.42(f),0.58(g),0.74(h),and 0.96(i).

    During the test,the enclosed paraffin wax is heated from the left side by using an electric heating plate,which is made of highly conductive silicon paste. A digital temperature controller is used to control the heating temperature,which is set to be 70°C in the melting test. In order to ensure the accuracy of heating temperature, the sensitivity and reliability of digital temperature controller have been calibrated by Ktype thermocouples prior to the melting test. To achieve a uniform heating source, a 2-mm-thick copper plate is glued on the electric heater to contact the paraffin wax. In Fig. 2, the experimental results show the evolution of liquid fractionfin the entire melting process of paraffin. Typically,f=0 when the PCM is solid, andf=1 when the PCM is liquid. At the beginning of melting, the heat transfer is dominated by heat conduction due to the direct contact between the heating wall and the paraffin wax, and a very thin liquid layer is formed rapidly. In this region, because the layer thickness is small,the motion of molten PCM is extremely slight. As a result,the melting front of paraffin almost keeps parallel to the heating wall as shown in Fig. 2(b). Subsequently, the paraffin in the upper portion melts faster than in the lower part, so that the solid–liquid interface becomes inclined,which may be attributed to the accelerating effect of natural convection on the liquid paraffin. Moreover,as the melting process goes on,the effect of natural convection becomes more significant. This can be demonstrated by the more and more inclined melting interface.

    3. Computational model

    3.1. Mathematical formulations

    Compared with the experimental results, the numerical results can provide detailed field variations inside the cavity.Here,in order to further investigate the heat transfer behaviors in the paraffin in the melting process at different size levels of the cavity, it is necessary to establish an appropriate computational model to describe the complex melting process of paraffin in the rectangular cavity which may involve with heat conduction, phase change, liquid-phase flow caused by natural convection,solid-phase motion,and moving of solid–liquid interface. In this work,a two-dimensional(2D)computational model along the square cross-section of the cavity is developed based on the symmetrical configuration of the experimental model and the applied boundary conditions, and its details can be found in Fig. 3. The following basic assumptions on the paraffin are made in the mathematical model: i)the fluid phase is incompressible, ii) the physical properties for each phase keep unchanged,iii)mushy zone exists during the melting of PCM,iv)the volume expansion in the melting process is ignored,and v)there is no heat generation inside the PCM.

    Based on the above assumptions,the transient governing equations for the problem under consideration are written below.

    Continuity equation

    whereρis the density,?is the Hamiltonian operator,uis the velocity vector,pis the pressure,μis the dynamic viscosity,Fis the volume force vector,Tis the absolute temperature,Cpis the specific heat capacity,kis the thermal conductivity,andLmis the latent heat of fusion.

    In momentum equation(2), when the density fluctuation caused by temperature variations is small,we can assume that the flow of the liquid phase is driven by a buoyancy force,which can be given by the following Boussinesq approximation:

    whereTrefis the reference temperature and can be taken as the melting temperatureTmfor the phase change problem.αis the thermal expansion coefficient,andgis the acceleration of gravity.

    Besides,the solid phase can be regarded as a liquid with much higher viscosity than the liquid phase, so the dynamic viscosity in Eq.(2)practically is written in the form of

    whereμs=108μlis the dynamic viscosity of solid phase andμlis the dynamic viscosity of liquid phase. Hereafter,the subscript“s”represents the solid phase,and the subscript“l(fā)”denotes the liquid phase.

    In Eq.(3),the effective thermal conductivityk,the effective densityρ, and the specific heat capacityCpcan respectively be given respectively by

    whereβis the mass fraction defined as

    In practice, in order to ensure the solution stability , the phase change is assumed to occur in a small temperature range[Ts,Tl], whereTsandTlusually are set to beTm??T/2 andTm+?T/2,respectively,and ?T=Tl?Tsis a small temperature interval,[28]i.e., 1 K. For this case, the phase change interphase is modelled as a mushy zone in which the value of liquid fraction can be evaluated by the linear interpolation of temperature

    3.2. Boundary and initial conditions

    In order to solve the unsteady PCM melting problem,appropriate boundary conditions are applied to the boundary of cavity. Initially, the temperature of the whole PCM is set to be the room temperatureT0=16°C,the velocity filed in it is zero and the pressure is 1 atm(1 atm=1.01325×105Pa). Additionally,the heating temperature is maintained at 70°C and the remaining boundaries are insulated. Besides,the pressure point constraint is applied to point B as shown in Fig.3. Such an operation is helpful in improving the convergence of phase change analysis when the natural convection is considered.

    Fig.3. Finite element mesh for melting computation model of paraffin.

    3.3. Meshing operation

    In this study,the COMSOL Multiphysics 5.3a is utilized to simulate the melting process occurring in the paraffin enclosed in a rectangular cavity by considering the heat conduction in the solid phase,the heat convection in the liquid phase,and the motion of the solid phase. The unstructured quadrilateral mesh containing 90009 DOFs is used to model in the computational domain. In order to exactly describe the flow velocity nearby the boundary, a boundary layer mesh is created as shown in Fig.3. The number of boundary layers is 12,and the boundary layer stretching factor is 20%to smooth the transition in element size from the boundary layer mesh to the interior mesh. The setting of boundary layer properties can achieve a better mesh.

    4. Results and discussion

    In this section, a comparison between the experimental results and the numerical results is first presented for the purpose of validation. Then,a detailed parameter sensitivity analysis is discussed.

    4.1. Validation

    Figure 4 shows the simulated melting fronts in the paraffin PCM at different liquid fractions. The blue color and red color represent the solid and liquid phase, respectively, and the yellow arrow denotes the velocity direction. It is found from Fig.4 that the evolution of the melting front is basically consistent with the experimental results given in Fig. 2. One can see that from Fig.4 that the portion of the solid phase descends gradually as expected. Atf=0.02,the melting interface is nearly vertical due to the dominated conduction effect.Then, owing to the effect of natural convection, a clockwise whirlpool flow indicated by the yellow velocity arrow in the molten liquid region appears clearly and becomes significant in accelerating the melting of the solid phase near the upper part of the cavity. As a result,the melting front becomes more and more inclined as the melting process goes on. It can be concluded that the melting process is initially dominated by the heat conduction and then by the natural convection under one-dimensional(1D)heat transfer conditions.

    To accurately describe the melting process, the liquid fraction of paraffin is calculated, and the fesults are shown in Fig. 5. It is found that the liquid fraction increases with time going by, but the melting rate slows down. This can be attributed to the fact that as the distance between the melting front and the heating wall increases with time going by, the heating effect of the hot wall becomes weak and the intensity of the whirlpool flow decreases. Moreover,such a variation of the liquid fraction approximately follows a power law.

    To further investigate the time-dependent variation of temperature in the melting process, six monitoring points are selected as shown in Fig. 3. The temperature variation at each monitoring point is depicted in Fig. 6. From Fig. 6,it is observed that the temperatures at monitoring points increase monotonically before they are lower than the melting temperature of paraffin. Once the temperature exceeds 300 K,the temperature curves rise rapidly with time. After that, all the temperatures slightly decrease. Finally, each temperature curve gradually reaches the heating temperatureTw=70°C.

    Fig.4. Evolution of numerical melting fronts in entire melting process at f =0(a),0.02(b),0.06(c),0.14(d),0.24(e),0.42(f),0.58(g),0.74(h),and 0.96(i).

    Fig.5. Computed and fitted curves of variation of liquid fraction with time.

    The temperature change at each point can be further explained by Fig. 7. It is found from Fig. 7(a) that there is a high temperature gradient near the heating wall and the melting front. Moreover, as the melting process goes on, it is observed from Fig.7(b)that thermally stable and unstable fluid regions exist together, and the separating line between them is close to the upper segment of the melting front. In the stable layer, the thermal energy is dominated by conduction.Whereas, in the unstable layer, the heat transfer is dominated by convection, and the hot fluid flow is basically wrapped up in the boundary layer of the melting front. As illustrated in Fig. 7(b), a triangular low temperature region is formed due to the effect of coupling between conduction and convection.Hence,the temperature at point B located in the triangular region is slightly lower than the values at points A and C that are located outside of the low temperature region. Such a temperature discrepancy can be used to explain why all the curves in Fig.6 show a slightly descent trend after the rapid rise in the melting process. Furthermore, owing to the presence of thermal boundary layer and velocity boundary layer which closely follows the melting front,the higher grid density is needed to accurately capture the positions of melting front.

    Fig.6. Variations of temperature with melting time at six monitoring points in melting process.

    Fig.7. Temperature contours at different time moments.

    4.2. Effect of natural convection

    To further ascertain the features of the natural convection,the case without fluid flow is investigated for comparison.Figure 8 displays that the melting front of the paraffin gradually moves rightward and keeps almost vertical at each timemoment, when the heat transfer is controlled by the heat conduction. Hence, the temperature curves at the six monitoring points(see Fig.3)are completely identical as shown in Fig.9.It is found that the paraffin at each of these monitoring points begins to melt at aroundt=38300 s at which the temperature approximately reaches the melting temperature 300 K and the liquid fraction is 0.5. Subsequently, the temperature at each of these points continues to increase by the heat conduction in the liquid PCM. Untilt=133860 s, the cavity is fully filled with the liquid paraffin.

    By comparing the melting times with and without considering the natural convection, it is found that the contribution of natural convection can make the melting front inclined and, more importantly, can reduce the total melting time by 92.36%. Thus,the natural convection in the liquid phase plays an important role in accelerating heat transfer in the melting process and should be considered. However, how to quantitatively assess its accelerating effect is still an open issue. In the next section, the dependence of natural convection on the cavity size is studied and its influence on the total melting time is analyzed numerically.

    Fig.8. Snapshots of melting front at different time moments:t=3600 s(a),38300 s(b),133860 s(c)(fully melted).

    Fig.9. Variations of temperature with melting time at different monitoring points in melting process.

    4.3. Size-dependent natural convection

    To investigate the effect of the cavity size on the convection heat transfer, a comprehensive analysis is performed in this work, and the results are shown in Figs. 10 and 11,in which the computation stops and the total melting time is recorded when the melting of paraffin is finished. As can be seen from Fig.10(a)which corresponds to the case without the natural convection, the log–log curves of the liquid fractionversusmelting time are approximately straight lines for different cavity sizes, which means that the liquid fraction varying with time follows a power law. Comparatively, under the influence of natural convection,the log–log curves in Fig.10(b)shows different behaviors. These curves initially keep straight and then their slopes increase with time. This indicates that the heat transfer is initially dominated by heat conduction and then the natural convection becomes strong and accelerate the melting rate. However,the melting rate slows down when the liquid fraction is over 0.6. This can be attributed to the fact that the strength of whirlpool becomes weak as the distance between the heating wall and the melting front increases. This stage accounts for around 30%of the total melting time.

    Figure 11 depicts the total melting time against the cavity size in the cases without and with the natural convection. It is clearly indicated that a smaller computational domain will melt faster than a larger one for the cases with and without natural convection and the relationship between the total melting time and the cavity side exhibits a power law, especially for the case with natural convection. Also,it is seen that the total melting time is significantly reduced when the natural convection present,and such an accelerating effect becomes stronger as the cavity size becomes bigger. Simultaneously, one observes in Fig.11 that the total melting time is almost identical for the two cases when the cavity size is very smaller.So it can be concluded that the effect of natural convection is closely related to the cavity size.

    Fig.10. Variations of liquid fraction with melting time in the cases without(a)and with(b)natural convection.

    Fig.11. Curve of total melting times versus domain size for cases without and with natural convection.

    To clearly demonstrate the size-dependence natural convection, the ratio between the total melting times for the two cases is defined as

    wheretvandtddenote the total melting time corresponding to the cases with and without natural convection,respectively.

    Figure 12 shows the variation of the total melting ratio with the cavity size. From Fig. 12, it is found that the effect of natural convection is significant. For example, when the cavity size is 50 mm, the total melting ratio is 0.1205. This means that the total melting time with natural convection decreases 87.95% in comparison with that without natural convection. Moreover, the total melting ratio becomes more and more small as the domain size increases. This indicates that the contribution of natural convection becomes more and more significant with the increase of domain size. Also, it is seen that when the cavity size is 2 mm, the total melting time ratio is very close to 1. This means that the accelerating effect of natural convection can be neglected when the cavity size is 2 mm. This can be further explained by the shaping of the melting front when the liquid fraction is 0.5 as indicated in Fig. 13, from which it is observed that the melting front becomes more inclined with the increase of domain size. Partially,as the domain size is 2 mm,the melting front is almost vertical.

    Fig. 12. Ratios of melting times versus domain size for different domain sizes.

    Fig.13. Melting fronts for different domain sizes: L=2 mm(a),5 mm(b),10 mm(c),20 mm(d),50 mm(e),and 100 mm(f).

    5. Equivalent thermal conduction

    The results above clearly show that the natural convection plays a positive role in accelerating the melting of paraffin and reducing the total melting time, which is directly related to the absorbing rate of latent energy in the melting process of PCM. However, we must point out that when the natural convection of liquid phase is involved, the melting process of PCM becomes very complex and a large number of iterations are needed at each time step. As a result, completing each simulating job may take about several days. For example,with the same mesh for discretizing the 100 mm×100 mm cavity, the CPU computing time is about 20 h by using the desktop with I7-8700 CPU and 16 G memories, although the full melting of paraffin can practically finish within 3 h. Such a time-consuming computation severely limits the quick analysis of PCM melting and hinders the engineering applications of PCM.

    In order to simplify the calculation of PCM, the simplest treatment is to introduce an equivalent thermal conductivity concept which integrates the contribution of natural convection in the original PCM into the thermal conductivity of equivalent PCM and assumes that the heat conduction is a dominant energy transfer mode in the melting process of equivalent PCM. The other physical parameters of PCM like the density,heat capacity,melting temperature,and latent heat,keep unchanged. On the basis of the same total melting time,the equivalent heat conduction coefficientk?effof the equivalent PCM can be determined by the intersection between the actual total melting time curve and the changing curve plotted by changing the thermal conductivity of the equivalent PCM in a large range as indicated in Fig.14 for the case ofL=100 mm.In Table 2, the values of desired equivalent thermal conductivity are listed for different domain sizes. In Table 2,tcvandtcddenote the CPU computing time for the original PCM and the equivalent PCM, respectively. Obviously, for the equivalent PCM,the computing time dramatically decreases.In other words,the computing efficiency greatly increases. For example,forL=100 mm,the computation for the equivalent PCM is 165.7 times faster than that for the original PCM.Moreover,the results in Table 2 indicate that the equivalent thermal conductivity of the equivalent PCM is related to the domain size.Using the data fitting technique, an empirical correlation of equivalent heat conduction coefficient is written as

    Equation(11)is valid for the range of parameters used in simulation cases: 23.8

    Fig.14. Variation of total melting time with heat conduction coefficient of equivalent PCM for L=100 mm.

    Table 2. Equivalent heat conduction coefficient.

    Fig.15.Comparison between results calculated by FEM and Eq.(11),showing(a)effect of domain size and(b)effect of temperature difference.

    6. Conclusions

    To analyze the effect of natural convection on the heat transfer behavior and the energy storage capability in the paraffin and reveal the size-dependence of natural convection,the experimental and numerical investigations of the melting process in a rectangular cavity are carried out in this paper.The enthalpy-based numerical model is validated experimentally. The main conclusions are shown below.

    The heat conduction is dominant at the beginning of the melting process and the corresponding melting front keeps almost parallel to the heating wall. Then,as the melting region grows,the melting front gradually becomes inclined due to the appearance of clockwise whirlpool in the liquid phase caused by natural convection. However,the intensity of whirlpool becomes weak as the distance between the melting front and the heating wall increases.

    Compared with the case without natural convection,the natural convection can reduce the total melting time by 92.36% and increase the melting efficiency of paraffin forL=100 mm.

    The domain size plays a significant role in the convective heat transfer in the melting process, and the relation between total melting time and domain side exhibits a power law.When the domain size is less than 2 mm,the influence of natural convection can be ignored.

    In order to simplify the calculation,the equivalent thermal conductivity concept is proposed,and an empirical correlation of the equivalent thermal conductivity is obtained by data fitting. This simple solution can greatly reduce the computation time and ensure the absorbing rate of latent heat simultaneously. This solution provides an approximate method for the future analysis of complex engineering problems associated with phase change materials.

    In summary, the present study can enhance our understanding of the relation between the particle size and the latent energy storing function, and the proposed empirical correlation of equivalent heat conduction coefficient provides a great convenience for developing advanced cement-based composite incorporating PCM particles used for potential engineering applications in designing ecology-grain-storage structures.

    猜你喜歡
    王輝
    本期面孔
    遼河(2022年12期)2023-01-29 13:24:58
    Spin transport properties for B-doped zigzag silicene nanoribbons with different edge hydrogenations
    劇作家王輝
    火花(2021年10期)2021-11-04 09:23:52
    Chinese Acupuncture
    竹之韻
    Investigation of the performance of CF3I/c-C4F8/N2 and CF3I/c-C4F8/CO2 gas mixtures from electron transport parameters
    A novel male-associated marker for sex identifi cation in Ussuri catfishPseudobagrus ussuriensis*
    The Thought on PPP in China
    弟子規(guī)
    山西老西兒
    丝袜喷水一区| 久久精品国产综合久久久| 成人手机av| 久久精品国产亚洲av涩爱| 观看美女的网站| 国产国语露脸激情在线看| 色网站视频免费| 国产又色又爽无遮挡免| 丰满饥渴人妻一区二区三| 亚洲国产最新在线播放| 国产男人的电影天堂91| 街头女战士在线观看网站| 国产免费福利视频在线观看| 国产男女超爽视频在线观看| 丁香六月天网| 99久国产av精品国产电影| 亚洲精品国产色婷婷电影| 91aial.com中文字幕在线观看| 午夜福利在线免费观看网站| 亚洲国产毛片av蜜桃av| 水蜜桃什么品种好| 国产成人av激情在线播放| 日本vs欧美在线观看视频| 大话2 男鬼变身卡| 久久鲁丝午夜福利片| 国产精品国产av在线观看| 1024香蕉在线观看| 亚洲av综合色区一区| 一区二区三区四区激情视频| 一级毛片 在线播放| 日韩伦理黄色片| 最黄视频免费看| 国产成人精品久久二区二区91 | 久久ye,这里只有精品| 成人黄色视频免费在线看| 久久毛片免费看一区二区三区| 韩国精品一区二区三区| 97人妻天天添夜夜摸| 国产1区2区3区精品| 免费在线观看黄色视频的| 亚洲欧美精品自产自拍| 国产探花极品一区二区| 精品视频人人做人人爽| 成年av动漫网址| 丰满饥渴人妻一区二区三| 国产精品女同一区二区软件| 丰满乱子伦码专区| av又黄又爽大尺度在线免费看| 日韩不卡一区二区三区视频在线| 国产女主播在线喷水免费视频网站| 国产男女超爽视频在线观看| 18禁国产床啪视频网站| 巨乳人妻的诱惑在线观看| 国产日韩欧美亚洲二区| 亚洲精品一区蜜桃| 蜜桃国产av成人99| 国产精品一区二区在线不卡| 亚洲欧美成人综合另类久久久| 男女边摸边吃奶| 母亲3免费完整高清在线观看 | 成年美女黄网站色视频大全免费| 看十八女毛片水多多多| 国产不卡av网站在线观看| videos熟女内射| 日日啪夜夜爽| 国产亚洲一区二区精品| 丰满饥渴人妻一区二区三| 美女视频免费永久观看网站| 美女国产高潮福利片在线看| 精品少妇内射三级| 午夜福利视频精品| 丰满迷人的少妇在线观看| 一本—道久久a久久精品蜜桃钙片| freevideosex欧美| 免费观看av网站的网址| 国产一区二区三区综合在线观看| 久久国产精品男人的天堂亚洲| 中文字幕亚洲精品专区| 国产精品免费视频内射| 亚洲欧美清纯卡通| 免费观看无遮挡的男女| 丰满乱子伦码专区| 中文精品一卡2卡3卡4更新| 少妇被粗大猛烈的视频| 午夜日本视频在线| 天堂中文最新版在线下载| 大话2 男鬼变身卡| 亚洲精品美女久久av网站| 亚洲国产精品999| 人人妻人人爽人人添夜夜欢视频| 在线精品无人区一区二区三| 国产亚洲av片在线观看秒播厂| 国产成人精品久久二区二区91 | 丰满少妇做爰视频| 免费在线观看视频国产中文字幕亚洲 | 999久久久国产精品视频| 亚洲精品中文字幕在线视频| 国产1区2区3区精品| 亚洲一级一片aⅴ在线观看| 亚洲美女视频黄频| 国产成人免费观看mmmm| 制服诱惑二区| 欧美精品一区二区大全| 国产精品av久久久久免费| 黑丝袜美女国产一区| 90打野战视频偷拍视频| 1024香蕉在线观看| 考比视频在线观看| 午夜免费观看性视频| 九九爱精品视频在线观看| 国产欧美日韩一区二区三区在线| 午夜精品国产一区二区电影| 三级国产精品片| 欧美成人午夜精品| 国产成人aa在线观看| 亚洲国产色片| 一级黄片播放器| 韩国高清视频一区二区三区| 一级a爱视频在线免费观看| 国产男女内射视频| 成人国产麻豆网| 少妇被粗大猛烈的视频| 国产97色在线日韩免费| 久久久久久伊人网av| 亚洲美女视频黄频| 精品国产一区二区久久| 少妇的丰满在线观看| 日日撸夜夜添| 精品午夜福利在线看| 七月丁香在线播放| 一区二区av电影网| 女人精品久久久久毛片| 搡老乐熟女国产| 91精品伊人久久大香线蕉| 男女高潮啪啪啪动态图| 精品一区二区免费观看| 最新中文字幕久久久久| 男人操女人黄网站| 午夜av观看不卡| 如日韩欧美国产精品一区二区三区| 欧美日韩成人在线一区二区| 亚洲精品视频女| 午夜日本视频在线| 2022亚洲国产成人精品| 这个男人来自地球电影免费观看 | 性高湖久久久久久久久免费观看| 不卡视频在线观看欧美| 久久久久国产网址| 在线天堂最新版资源| 1024视频免费在线观看| 免费高清在线观看日韩| 久久久久久久久久久免费av| 五月开心婷婷网| 欧美最新免费一区二区三区| 满18在线观看网站| 天天躁夜夜躁狠狠久久av| 永久免费av网站大全| 成人手机av| 美女国产视频在线观看| av网站在线播放免费| 一本—道久久a久久精品蜜桃钙片| 熟女少妇亚洲综合色aaa.| 天堂8中文在线网| 成人影院久久| 精品少妇久久久久久888优播| 99九九在线精品视频| 狂野欧美激情性bbbbbb| 国产精品免费视频内射| 国产乱来视频区| 超碰97精品在线观看| 九九爱精品视频在线观看| 亚洲国产毛片av蜜桃av| 考比视频在线观看| 伦精品一区二区三区| 国产高清不卡午夜福利| 男女下面插进去视频免费观看| 欧美xxⅹ黑人| 看非洲黑人一级黄片| xxx大片免费视频| 欧美日韩视频高清一区二区三区二| 人人妻人人添人人爽欧美一区卜| 丝袜美足系列| 国产黄色视频一区二区在线观看| 少妇人妻精品综合一区二区| 午夜精品国产一区二区电影| 日韩电影二区| 亚洲成色77777| 在线看a的网站| 成年人免费黄色播放视频| 一级爰片在线观看| 久久狼人影院| 国产精品二区激情视频| 男女免费视频国产| 成人国产麻豆网| 中文精品一卡2卡3卡4更新| 一级毛片电影观看| 午夜91福利影院| 国产片特级美女逼逼视频| 韩国av在线不卡| 久久精品国产亚洲av天美| 精品亚洲成国产av| 亚洲人成网站在线观看播放| 亚洲av电影在线进入| a级毛片在线看网站| 男女边摸边吃奶| 2022亚洲国产成人精品| 在线观看免费日韩欧美大片| 久久精品亚洲av国产电影网| 日韩伦理黄色片| 18禁裸乳无遮挡动漫免费视频| 午夜久久久在线观看| 国产日韩一区二区三区精品不卡| av在线播放精品| 咕卡用的链子| 亚洲一区中文字幕在线| 人妻系列 视频| 18禁观看日本| 一级爰片在线观看| a级片在线免费高清观看视频| 亚洲国产日韩一区二区| 国产有黄有色有爽视频| 成人毛片a级毛片在线播放| 久久99蜜桃精品久久| 一级毛片 在线播放| 午夜av观看不卡| 中文字幕色久视频| 欧美日韩视频高清一区二区三区二| 男人爽女人下面视频在线观看| 国产精品免费视频内射| a 毛片基地| 国产黄色视频一区二区在线观看| 国产欧美日韩一区二区三区在线| 久久精品aⅴ一区二区三区四区 | 久久人妻熟女aⅴ| 欧美精品一区二区大全| 亚洲熟女精品中文字幕| 激情五月婷婷亚洲| 色婷婷久久久亚洲欧美| 99久久人妻综合| 日韩在线高清观看一区二区三区| 91精品三级在线观看| av.在线天堂| 日韩不卡一区二区三区视频在线| 久久久久人妻精品一区果冻| 水蜜桃什么品种好| 一区二区三区激情视频| 国产欧美日韩一区二区三区在线| 久久狼人影院| 国产成人精品在线电影| 亚洲精品视频女| 国产成人免费观看mmmm| 天堂8中文在线网| 我要看黄色一级片免费的| 少妇熟女欧美另类| 久久这里有精品视频免费| 蜜桃在线观看..| 大话2 男鬼变身卡| 日韩一卡2卡3卡4卡2021年| 黄色怎么调成土黄色| av电影中文网址| 免费看av在线观看网站| 欧美老熟妇乱子伦牲交| 大香蕉久久网| 2022亚洲国产成人精品| 久久久久视频综合| 日韩免费高清中文字幕av| 国产深夜福利视频在线观看| 在线观看免费视频网站a站| 久久av网站| 看十八女毛片水多多多| 男的添女的下面高潮视频| 中文字幕色久视频| 国产1区2区3区精品| 丰满饥渴人妻一区二区三| 国产av一区二区精品久久| 亚洲欧美一区二区三区黑人 | 国产日韩一区二区三区精品不卡| 91精品国产国语对白视频| 欧美日韩av久久| 欧美成人午夜免费资源| 人妻 亚洲 视频| 黑丝袜美女国产一区| 91aial.com中文字幕在线观看| av网站在线播放免费| 婷婷色综合www| 午夜影院在线不卡| 国产免费现黄频在线看| 欧美精品一区二区大全| 国产精品秋霞免费鲁丝片| 欧美另类一区| 久久韩国三级中文字幕| 免费黄色在线免费观看| 极品少妇高潮喷水抽搐| 高清视频免费观看一区二区| 亚洲欧美精品自产自拍| 观看av在线不卡| 国产亚洲最大av| 免费在线观看视频国产中文字幕亚洲 | 亚洲国产日韩一区二区| 中文字幕人妻丝袜制服| 日韩伦理黄色片| 中国国产av一级| 侵犯人妻中文字幕一二三四区| 国产亚洲精品第一综合不卡| 99国产综合亚洲精品| 国产熟女午夜一区二区三区| 黄色怎么调成土黄色| 久久久久久人人人人人| 午夜福利乱码中文字幕| 亚洲人成77777在线视频| 中文字幕人妻熟女乱码| 久久久久人妻精品一区果冻| 狂野欧美激情性bbbbbb| 久久久欧美国产精品| 婷婷色麻豆天堂久久| 亚洲成国产人片在线观看| 久久国产亚洲av麻豆专区| 中文字幕制服av| 岛国毛片在线播放| 成人国产麻豆网| 国产免费视频播放在线视频| 一区二区三区激情视频| 在线观看人妻少妇| 亚洲成色77777| 99热国产这里只有精品6| 天美传媒精品一区二区| 青春草国产在线视频| 丁香六月天网| 美女高潮到喷水免费观看| 午夜av观看不卡| 人人妻人人添人人爽欧美一区卜| 99热全是精品| 国产午夜精品一二区理论片| 肉色欧美久久久久久久蜜桃| 日韩伦理黄色片| 国产1区2区3区精品| 亚洲少妇的诱惑av| 免费观看在线日韩| 精品一区二区免费观看| 亚洲婷婷狠狠爱综合网| 涩涩av久久男人的天堂| 女性生殖器流出的白浆| 国产麻豆69| 国产精品久久久久久av不卡| 伊人久久国产一区二区| 欧美激情高清一区二区三区 | 久久久久精品人妻al黑| 久久久亚洲精品成人影院| 少妇的丰满在线观看| 我要看黄色一级片免费的| 成年女人在线观看亚洲视频| 麻豆av在线久日| 亚洲av电影在线观看一区二区三区| 欧美精品亚洲一区二区| 日韩电影二区| 日韩人妻精品一区2区三区| 免费不卡的大黄色大毛片视频在线观看| 国产免费一区二区三区四区乱码| 亚洲国产毛片av蜜桃av| 精品国产乱码久久久久久男人| 菩萨蛮人人尽说江南好唐韦庄| 美国免费a级毛片| 亚洲国产毛片av蜜桃av| 男女下面插进去视频免费观看| 一边亲一边摸免费视频| 热re99久久国产66热| 精品国产超薄肉色丝袜足j| 国产精品久久久久久久久免| 国产毛片在线视频| 国产在线视频一区二区| 国产精品国产三级专区第一集| 制服人妻中文乱码| 欧美成人午夜精品| 熟女av电影| 青草久久国产| 免费观看性生交大片5| 18禁动态无遮挡网站| 免费大片黄手机在线观看| 欧美bdsm另类| 嫩草影院入口| 亚洲少妇的诱惑av| 免费观看av网站的网址| av又黄又爽大尺度在线免费看| 老司机亚洲免费影院| 久久精品熟女亚洲av麻豆精品| a级毛片在线看网站| 亚洲欧洲精品一区二区精品久久久 | 人妻一区二区av| 午夜av观看不卡| 美女福利国产在线| tube8黄色片| 如日韩欧美国产精品一区二区三区| 人人妻人人添人人爽欧美一区卜| 久热久热在线精品观看| 精品少妇黑人巨大在线播放| 国产精品国产三级专区第一集| 亚洲精品国产色婷婷电影| 人妻系列 视频| 欧美最新免费一区二区三区| 欧美日韩亚洲国产一区二区在线观看 | 丝袜美腿诱惑在线| 久久久久视频综合| 欧美日韩精品网址| 一本一本久久a久久精品综合妖精 国产伦在线观看视频一区 | 日韩制服骚丝袜av| 9191精品国产免费久久| 亚洲国产看品久久| 街头女战士在线观看网站| 国产成人av激情在线播放| 18禁观看日本| 你懂的网址亚洲精品在线观看| 一边亲一边摸免费视频| 建设人人有责人人尽责人人享有的| 欧美激情极品国产一区二区三区| 美女国产视频在线观看| 久久久久精品性色| 亚洲国产精品999| 亚洲欧美成人综合另类久久久| 亚洲欧美一区二区三区久久| 免费黄频网站在线观看国产| 只有这里有精品99| 久久精品国产综合久久久| 91久久精品国产一区二区三区| 秋霞在线观看毛片| 深夜精品福利| 极品人妻少妇av视频| 在线免费观看不下载黄p国产| 国产精品二区激情视频| 日本wwww免费看| 欧美人与性动交α欧美软件| 国产又色又爽无遮挡免| 久久久精品国产亚洲av高清涩受| 国产精品99久久99久久久不卡 | 巨乳人妻的诱惑在线观看| 黄频高清免费视频| 欧美国产精品一级二级三级| 天堂俺去俺来也www色官网| 91精品三级在线观看| 国产成人免费观看mmmm| 亚洲欧美日韩另类电影网站| 欧美日韩综合久久久久久| 亚洲国产欧美在线一区| 国产精品一二三区在线看| av有码第一页| 激情五月婷婷亚洲| 18禁动态无遮挡网站| 亚洲欧美一区二区三区黑人 | 伊人久久国产一区二区| 久久人妻熟女aⅴ| 精品第一国产精品| 国产亚洲精品第一综合不卡| 亚洲精品,欧美精品| 春色校园在线视频观看| 亚洲,欧美精品.| 亚洲av成人精品一二三区| av免费在线看不卡| 日韩av免费高清视频| 热re99久久精品国产66热6| 美女脱内裤让男人舔精品视频| 两个人看的免费小视频| 一个人免费看片子| 国产 一区精品| 亚洲欧美中文字幕日韩二区| 亚洲国产日韩一区二区| 欧美另类一区| 亚洲伊人色综图| 老鸭窝网址在线观看| 精品少妇内射三级| 亚洲,欧美精品.| 国产成人精品婷婷| 丝瓜视频免费看黄片| 国产黄频视频在线观看| 又黄又粗又硬又大视频| 欧美人与善性xxx| 免费不卡的大黄色大毛片视频在线观看| 麻豆精品久久久久久蜜桃| 国产免费福利视频在线观看| 亚洲,一卡二卡三卡| 赤兔流量卡办理| 国产精品亚洲av一区麻豆 | 日本欧美国产在线视频| 久久鲁丝午夜福利片| 欧美中文综合在线视频| 久久人人97超碰香蕉20202| 视频在线观看一区二区三区| 王馨瑶露胸无遮挡在线观看| 久久精品aⅴ一区二区三区四区 | 欧美日韩视频高清一区二区三区二| 亚洲欧美色中文字幕在线| 69精品国产乱码久久久| 黄色视频在线播放观看不卡| 久久精品人人爽人人爽视色| 天堂8中文在线网| 久久久a久久爽久久v久久| 91国产中文字幕| 人妻少妇偷人精品九色| 精品国产一区二区三区久久久樱花| 人人妻人人添人人爽欧美一区卜| 午夜免费观看性视频| 久久久久视频综合| 亚洲图色成人| 超色免费av| 人妻少妇偷人精品九色| 亚洲中文av在线| 欧美另类一区| tube8黄色片| 日韩精品有码人妻一区| 免费黄色在线免费观看| 亚洲欧美精品自产自拍| 亚洲婷婷狠狠爱综合网| 看免费成人av毛片| 激情五月婷婷亚洲| 国产又爽黄色视频| 日本vs欧美在线观看视频| 精品国产超薄肉色丝袜足j| 人成视频在线观看免费观看| av不卡在线播放| 国产av码专区亚洲av| 久久久久人妻精品一区果冻| 国产精品久久久av美女十八| av福利片在线| 97精品久久久久久久久久精品| 国产精品av久久久久免费| 久久国产亚洲av麻豆专区| 亚洲精华国产精华液的使用体验| 色吧在线观看| www.自偷自拍.com| 国产一区二区三区av在线| 久久久久精品性色| 高清黄色对白视频在线免费看| av电影中文网址| 色94色欧美一区二区| 青青草视频在线视频观看| 久久精品熟女亚洲av麻豆精品| 色婷婷av一区二区三区视频| 欧美 日韩 精品 国产| 亚洲av国产av综合av卡| 老司机影院成人| 美女脱内裤让男人舔精品视频| 亚洲 欧美一区二区三区| 蜜桃在线观看..| 国产激情久久老熟女| av在线老鸭窝| 亚洲av成人精品一二三区| 婷婷成人精品国产| 日产精品乱码卡一卡2卡三| 久久久久精品性色| 日本-黄色视频高清免费观看| 精品卡一卡二卡四卡免费| 在线天堂中文资源库| 美女xxoo啪啪120秒动态图| 日韩欧美一区视频在线观看| 99久久中文字幕三级久久日本| 欧美成人午夜免费资源| 欧美人与性动交α欧美精品济南到 | 国产综合精华液| kizo精华| 男女高潮啪啪啪动态图| 男女无遮挡免费网站观看| 免费不卡的大黄色大毛片视频在线观看| 自线自在国产av| 日韩中字成人| 色视频在线一区二区三区| 亚洲精华国产精华液的使用体验| www.精华液| 久久精品熟女亚洲av麻豆精品| 国产一区亚洲一区在线观看| av在线app专区| 菩萨蛮人人尽说江南好唐韦庄| 国产深夜福利视频在线观看| 欧美bdsm另类| 欧美日韩亚洲高清精品| 国产精品一国产av| 国产无遮挡羞羞视频在线观看| 91aial.com中文字幕在线观看| 在线观看www视频免费| 亚洲人成电影观看| 欧美成人午夜精品| 成人手机av| 激情视频va一区二区三区| 黄片无遮挡物在线观看| 韩国高清视频一区二区三区| 青草久久国产| 老汉色∧v一级毛片| 亚洲精品,欧美精品| 侵犯人妻中文字幕一二三四区| 久久久亚洲精品成人影院| 天天躁日日躁夜夜躁夜夜| 99精国产麻豆久久婷婷| 人妻人人澡人人爽人人| 精品国产乱码久久久久久男人| 国产精品久久久久成人av| 亚洲欧美精品自产自拍| 国产片特级美女逼逼视频| 日韩av不卡免费在线播放| 97精品久久久久久久久久精品| 午夜日本视频在线| 在线免费观看不下载黄p国产| 最新的欧美精品一区二区| 人妻系列 视频| 制服诱惑二区| 纯流量卡能插随身wifi吗| 26uuu在线亚洲综合色| 亚洲美女搞黄在线观看| 亚洲国产av新网站| 丁香六月天网| 精品午夜福利在线看| 午夜影院在线不卡| 熟女少妇亚洲综合色aaa.| 午夜福利一区二区在线看| 亚洲国产精品一区三区| 中文字幕人妻丝袜一区二区 | 色哟哟·www| 一区二区日韩欧美中文字幕| 男女高潮啪啪啪动态图|