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

    Effect of Radial Porosity Oscillation on the Thermal Performance of Packed Bed Latent Heat Storage

    2021-08-30 07:05:08LiuZhao
    Engineering 2021年4期

    H.B.Liu,C.Y.Zhao*

    Institute of Engineering Thermo-Physics,School of Mechanical Engineering,Shanghai Jiao Tong University,Shanghai 200240,China

    Keywords:

    Packed bed latent heat storage Radial porosity oscillation Fluid flow Heat transfer Diameter ratio

    ABSTRACT Owing to its high heat storage capacity and fast heat transfer rate,packed bed latent heat storage(LHS)is considered as a promising method to store thermal energy.In a packed bed,the wall effect can impact the packing arrangement of phase change material(PCM)capsules,inducing radial porosity oscillation.In this study,an actual-arrangement-based three-dimensional packed bed LHS model was built to consider the radial porosity oscillation.Its fluid flow and heat transfer were analyzed.With different cylindrical sub-surfaces intercepted along the radial direction in the packed bed,the corresponding relationships between the arrangement of capsules and porosity oscillation were identified.The oscillating distribution of radial porosity led to a non-uniform distribution of heat transfer fluid(HTF)velocity.As a result,radial temperature distributions and liquid fraction distributions of PCMs were further affected.The effects of different dimensionless parameters(e.g.,tube-to-capsule diameter ratio,Reynolds number,and Stefan number)on the radial characteristics of HTF and PCMs were discussed.The results showed that different diameter ratios correspond to different radial porosity distributions.Further,with an increase in diameter ratio,HTF velocity varies significantly in the near wall region while the non-uniformity of HTF velocity in the center region will decrease.The Reynolds and Stefan numbers slightly impact the relative velocity distribution of the HTF—while higher Reynolds numbers can lead to a proportional improvement of velocity,an increase in Stefan number can promote heat storage of the packed bed LHS system.

    1.Introduction

    Recently,solar energy has been extensively used for its advantages of cleanness,safety,and inexhaustibility[1,2].According to REN21,the capacity of a solar photovoltaic(PV)and concentrating solar power(CSP)was up-regulated from 409.9 GW in 2017 to 510.5 GW in 2018,occupying 60.97%of the renewable energy growth[3].However,solar energy is intermittent and fluctuating,and requires energy storage technology to address its defects and improve energy efficiency.Considering the thermal utilization of solar energy,latent heat storage(LHS)can exhibit a higher energy density and maintain an almost constant temperature during phase change,which brought widespread attention[4,5].The main constraint of LHS is the low thermal conductivity of phase change materials(PCMs)(0.2-0.8 W·m-1·K-1),and requires effective heat transfer enhancements[6].Considering the heat storage unit,a packed bed is considered an efficient approach when compared with a shell-and-tube unit.As reported by Li et al.[7],the charging and discharging rates of a packed bed LHS system is 1.8-3.2 times that of the shell-and-tube heat storage system.

    In the packed bed LHS system,as an heat transfer fluid(HTF)directly flows through the PCM capsules,the arrangement of PCM capsules considerably impacts the flow field[8,9].For a sphere packing,the densest arrangements of the spheres include hexagonal packing and face-centered cubic packing,exhibiting a highest packing density as 74.05%[10].However,when spheres are randomly packed in a container,they cannot be densely arranged because of the wall,which is called the wall effect[11].Mueller[12]reported that radial porosity would oscillate along the radial direction in the packed bed.Moreover,the oscillation distribution of radial porosity causes non-uniformity of the flow filed,which would further affect the heat transfer inside the packed bed[13].Considering a cylindrical container,the wall effect on the average packing density is prominent with a decrease in diameter ratio of the container and spheres(D/dp)[14].When a packed bed is used for latent heat storage,higher heat storage capacity leads to smaller tank size and diameter ratio.Accordingly,the radial porosity variation should be considered in a numerical study on packed bed LHS systems,particularly for small diameter ratio(D/dp<10)[15].

    Several numerical models were built to assess the thermal performance of the packed bed LHS system.Based on the requirement for the specific arrangement of PCM capsules,the models can be divided into representative elementary volume(REV)scale model and pore scale model.In REV scale model,governing equations are adopted to calculate the thermal performance of HTF and PCM.Benmansour et al.[16]adopted the two-phase model to express the energy equations of HTF and PCMs,while the porosity employed in the equations is a constant;thus,the radial porosity variation was not considered[17].Subsequently,a monotonic exponential equation,where the porosity refers to a function of radial distance from the wall,was proposed to replace the constant porosity[18,19].To obtain the thermal gradients inside the PCM capsules,the dispersion-concentric(DC)model was proposed.Karthikeyan and Velraj[20]compared the two-phase model with the DC model;they reported that the results of the DC model complied better with their experimental results as the internal heat transfer inside the PCM capsules was considered.However,a packed bed with a small diameter ratio(D/dp<10)cannot be considered as a continuous porous medium.Since representative control volume should be small enough to gain good calculation results,a control volume would not contain even more than one capsule,which fails to satisfy the requirement of the REV scale model[21,22].

    In the pore scale model,the concrete arrangement of PCM capsules should be described;subsequently,the flow of HTF and the heat transfer between HTF and PCM capsules can be calculated.Xia et al.[23]developed an effective packed model,where the three-dimensional(3D)packed bed is converted to a twodimensional packed bed.The corresponding criteria was to ensure that the two packed beds exhibit identical porosity.However,the radial porosity distribution is not considered during the transformation.Moreover,PCM capsules in the effective model are not in contact with each other,leading to the change of flow channels among PCM capsules.Therefore,the flow field is not accurate and heat conduction between capsules is ignored.

    Thus,complying with the actual packing process of spheres,a 3D packed bed LHS model was built to consider the radial porosity oscillation.The mechanism of the radial porosity oscillation and its impacts on the radial velocity of HTF,and radial temperature and liquid fraction of PCMs are analyzed.Further,the effects of different dimensionless parameters(e.g.,diameter ratio,Reynolds number,and Stefan number)on the fluid flow,heat transfer,and heat storage performance of the packed bed LHS are discussed.

    2.Model description

    2.1.Establishment of physical model

    In studying a 3D packed bed,a difficulty faced is building an effective sphere packing model;the falling,collision,and friction process of the spheres should be simulated.The existing methods of building a packed bed consist of the discrete equation method(DEM)method[8,24,25],Blender method[26,27],and Monte Carlo method[15,28].In this study,the open source software Blender was adopted to build a 3D packed bed model.The Blender software was equipped with a prominent physical simulator and bullet physics library;thus,the physical process of spheres when they fall under the influence of gravity can be effectively calculated.In the modelling process,a cylindrical container was first built;then,small spheres were constantly and randomly generated to fall above the container,as illustrated in Fig.1(a).The position of every sphere was calculated during the packing process.If the positions of all the spheres are stable,the packing process can be considered to be attaining a steady state.Then,the 3D packed bed model is established.

    Porosity(φv)is the ratio of void volume(Vvoid)to total volume(V),while surface porosity(φA)is the ratio of void area(Avoid)to total area(A)on surface,as expressed in Eqs.(1)-(3)[29].

    where R is the radius of container and r is an independent variable changing from 0 to R.Radial porosity,which is discussed in this study,is a phenomenon of surface porosity.To calculate the radial porosity,a series of cylindrical sub-surfaces along the radial direction were cut off.For the corresponding sub-surface,a ray-casting algorithm was adopted to determine the areas pertaining to the interior of the spheres.Then,the areas of all the spheres are added and total area of spheres is divided by the area of the sub-surface,as presented in Fig.1(b).

    Fig 1.Schematic diagrams of(a)packing process of sphere and(b)generation of sub-surfaces along the radial direction.

    Fig 2.Comparison of radial porosity distributions when diameter ratio is 7.99.

    In Fig.2,the radial porosity of the model built in the Blender was compared with the results of the Mueller,and the corresponding diameter ratio obtained is 7.99[12].The two curves oscillate and correspond well with each other,verifying the reliability of the model.Moreover,the oscillation distribution of radial porosity shows that the packed bed cannot be treated as a uniform porous media,further confirming the necessity of establishing a 3D model.

    In the packed bed,the contact type between the spheres is point contact,which is not conducive to the mesh generation.Thus,four methods were proposed to address contact points:overlaps,gaps,caps,and bridges.The bridge method can achieve the optimal results in the situation of heat transfer[30].Accordingly,small cylinders with a diameter of dp/10 were inserted between every two spheres in contact as heat bridges,thereby leading to changes in the average porosity of the packed bed by 0.05%-0.08%.Moreover,the gaps method was used to address the contact points between the spheres and container wall.After the diameter of the container was increased by 1 mm,a narrow gap was identified between the spheres and container wall.

    Further,a 3D packed bed LHS model was built,as presented in Fig.3.In the model,205 PCMcapsules were packed into a cylindrical tank with an inner diameter of 240 mm,and the diameter of PCM capsule is 48 mm;thus,the diameter ratio of the packed bed is 5.Ternary carbonate Li2CO3-K2CO3-Na2CO3(32 wt%-35 wt%-33 wt%)was selected as the PCM,and its thermal properties are listed in Table 1.The packed bed was placed in the middle section,while extended parts were added on the sides of the inlet and outlet to eliminate the impact of the sides.

    2.2.Mathematical model

    2.2.1.Assumptions

    (1)The effect of gravity on HTF flow is ignored;

    (2)The convection effect of the PCM in the capsules is ignored;

    (3)Radiation heat transfer of the HTF and PCM capsules is ignored.

    2.2.2.Governing equation of HTF

    In the packed bed LHS system,the pore Reynolds number based on mean velocity is expressed as[9]

    whereρis density,μis viscosity,and uinis the inlet velocity.

    Fig.3.Schematic diagram of 3D packed bed LHS system.

    For a packed bed with low diameter ratio,the flow is normally turbulent flow.Considering that the heat transfer surface of the PCM capsule is spherical,renormalization group(RNG)k-εturbulence model was selected to facilitate the wall treatment.Subsequently,the governing equations of the HTF in the packed bed LHS system can be written as follows[31]:

    Continuity equation:

    where Cμis the model constant;γis the RNG k-εturbulence model coefficient;γ0is the RNG k-εturbulence model coefficient in the initial state;S denotes the mean strain-rate of the flow;and Sijis the deformation tensor.

    Table 1 Thermal properties of PCM.

    The constants in the aforementioned equations areσT=1.0,σk=σ=0.72,Cμ=0.0845,C1=1.42,C2=1.68,andγ0=4.38.

    2.2.3.Governing equation of the PCMs in the capsules

    For the PCMs in the capsules,as the internal natural convection is ignored,the energy equation can be expressed as[23]

    where Hpcmis the enthalpy of PCM and defined as the sum of sensible heat and latent heat;?is the gradient operator.

    where H0,pcmis the enthalpy of the PCM at the initial temperature.ΔH represents the latent heat released during the charging process and defined asΔH=βL.L is the latent heat andβis the liquid fraction of PCM,which,at different temperatures,is written as follows:

    where Tpcm,T0,Ts,and Tldenote the temperature of the PCM,initial temperature,onset melting temperature,and ending melting temperature,respectively.Furthermore,the charging process of the PCM demonstrates a ratio of sensible heat to latent heat,which is defined as the Stefan number[33]:

    where Tinrefers to the temperature of the inlet HTF.

    2.2.4.Boundary conditions and initial conditions

    At the initial state,the packed bed LHS system was maintained at a constant temperature and heat transfer fluid temperature Tf=Tpcm=T0.During the charging process,the HTF at the inlet was maintained at a constant temperature and mass flow rate(q),while heat flux of the HTF in the direction normal to the outlet was assumed to be zero.Therefore,the boundary conditions of the HTF in the packed bed are expressed as follows:

    where z is the height.

    To understand the effect of radial porosity distribution on the radial characteristics of the HTF,the container wall is set to be adiabatic.

    2.3.Heat storage of packed bed LHS system

    The total heat storage in the packed bed LHS system consists of the heat stored in the PCM,PCM capsule shells,and heat stored in the tank wall.The heat stored in the PCMincludes the sensible heat in the solid and liquid state,and latent heat during phase change.The heat stored in the PCM capsule shells and tank refers to the sensible heat.The charging time of the system is defined as the time when the central temperature of the PCM capsule at the outlet reaches 1 K lower than the inlet temperature.Hence,the total heat storage(Qtotal)in the packed bed LHS system can be expressed as follows:

    where Qpcm,Qshell,and Qtankare the heat storage in the PCM,shell of PCMcapsules,and the tank,respectively;mpcm,mshell,and mtankare the quality of PCM,shell of PCM capsules,and the tank,respectively;cp,sand cp,lare the specific heat capacity of the PCM in solid state and liquid state,respectively;Tm,T0,and Tinare the melting temperature,the initial temperature,and the inlet temperature,respectively;cp,shelland cp,tankare the heat capacity of the shell and tank,respectively.

    During the charging process,heat storage of the packed bed LHS system (Qstored) can be obtained using the first law of thermodynamics:

    where mfis the mass flow rate of HTF;cpf,inand cpf,outare the specific heat capacity of HTF at the inlet and outlet,respectively;Tf,inand Tf,outare the temperature of HTF at the inlet and outlet,respectively.

    The average charging power Pavercan be defined as[7]

    whereτchargeis the charging time.

    2.4.Numerical procedure and validation

    As there are hundreds of PCM capsules randomly filled in the packed bed,unstructured grids were adopted to generate the meshes.The shrink-wrap method was employed to generate the surface mesh,and then,the volume was filled with the polyhedral meshes[34].Subsequently,the effects of grid size on the charging time and pressure drop in the packed bed were determined,and a grid size of dp/18 was selected for the mesh generation,as shown in Fig.4.During the calculation,pressure and velocity fields were computed using the SIMPLE algorithms.The spatial discretization and transient formulation were calculated using the second order upwind scheme and second order implicit scheme,respectively.

    Fig.4.Independence verification of mesh size.

    Li et al.[7]built a packed bed LHS system filled with 385 PCM capsules.The inner diameter of the container was 260 mm,while the capsule diameter was 34 mm;thus,a diameter ratio of 7.65 was achieved.The PCM in the capsules was ternary carbonate Li2CO3-K2CO3-Na2CO3(32 wt%-35 wt%-33 wt%);its thermal properties are listed in Table 1.The PCM temperature variations at different heights during the charging process were selected to verify the established 3D packed bed LHS model.The locations of the PCM capsules along the height were z/h=0.25,0.5,0.75,as shown in Fig.5(a).During the charging process,the mass flow rate qinwas 260 kg·h-1,while the initial temperature T0of the packed bed was 325°C and inlet temperature Tinof the HTF reached 465°C.

    As shown in Fig.5(b),the simulation data of the model complied well with the experimental results,and the results at z/H=0.5 were optimally matched.At the beginning of the charging process,the simulation was faster than the experiment,whereas,after the phase change process of the PCM capsule was achieved,the simulation data and experimental results were synchronized.On one hand,the heating process of the HTF to be heated from initial temperature(e.g.325°C)to the inlet temperature(e.g.465°C)was time-consuming,whereas,this transformation was attained immediately in the simulation.On the other hand,when the PCMs were in the liquid phase,the differences of the heat storage attributed to the different inlet temperature transition mode will be small when compared with the overall heat storage of the system.Accordingly,the simulation becomes more accurate over time.The 3D packed bed LHS system model exhibits a good accuracy.

    3.Results and discussion

    3.1.Flow and thermal performance analysis

    In this part,the initial temperature T0of the packed bed is 325°C,inlet temperature Tinof HTF is 465°C,and mass flow rate reaches 260 kg·h-1.As the average porosity is 0.437 in the packed bed exhibiting the diameter ratio of 5,the Reynolds number Repis calculated to be 4998.72,thus confirming that the flow in the packed bed is turbulent.

    3.1.1.Reason for radial porosity oscillation

    For the packed bed LHS system model with a diameter ratio of 5,the distribution of the radial porosity is illustrated in Fig.6.Several cylindrical sub-surfaces at different radial positions inside the packed bed were intercepted.The distance of the sub-surfaces from the tank wall reaches 0,0.25dp,0.5dp,0.75dp,dp,1.25dp,1.5dp,1.75dp,and 2dp,respectively.Considering the surface porosity,the porosity of the sub-surfaces(0,dp,and 2dpaway from the wall)is suggested to be higher than that of the other sub-surfaces.The oscillating distribution of the radial porosity is related to the process of packing capsules and spherical shape of the capsules.During the packing process,capsules tend to fill the space close to the wall first and subsequently fill the space in the center[35].At the wall surface,the PCM capsules exhibit point contacts with the tank wall.Hence,the porosity of this position is close to 1,indicating that the sub-surface is almost HTF.At the position 0.5dpaway from the wall,the sub-surface passes through the center of the capsules,revealing that the PCMarea will occupy most of the area of the subsurface,and the porosity will decrease to 0.15.Similarly,at the subsurface dpaway from the wall,the capsules near the wall are in point contacts with the inner capsules.Gaps are formed near the contact points,and the porosity increases with a value of 0.67.Thus,the porosity of the sub-surfaces in the radial direction complies with an oscillating distribution.

    3.1.2.Flow velocity distribution

    Similar to radial porosity,radial velocity is the average value of HTF velocity on the sub-surface,defined as Eq.(29).In this study,the relative velocity(U)is determined by the ratio of the radial velocity(u(r))to the inlet velocity(uin).

    whereΩfis the area of HTF on the sub-surface and Afis the area value of HTF.The relative velocity distribution of the HTF at different radial positions is presented in Fig.7,also indicating an oscillation trend.It is observed that although the porosity at the wall is close to 1,the flow velocity is 0 due to wall viscosity.However,in the radial direction away from the wall,velocity increases significantly and subsequently varies with the oscillation of radial porosity.When the HTF ascends into the packed bed,it flows through the gaps among the PCM capsules.Locations with more gaps allow more HTF to flow through;thus,the HTF velocity is faster where the porosity is larger.To facilitate the analysis and comparison of the velocity distribution,the relative velocity distribution was divided into two parts:near wall region and center region.The distance of the near wall region varies from the wall to 0.5dp,where velocity varies drastically,while the center region is from 0.5dpto the center,and the velocity varies more regularly.

    Fig.5.(a)Position illustrations of thermocouples;(b)comparison of calculation results with experimental data.

    In the near wall region,velocity varies drastically.Hence,the most important parameter here is the maximumrelative velocity near the wall.In this case,the maximum relative velocity at the near wall region is 3.49.According to the continuity equation,under the same flow rate,the variation of velocity is inversely proportional to the variation of flow channel area,as expressed in Eq.(31).Accordingly,the velocity of the incompressible HTF increases with a sudden contraction of flow channel.Fig.8 illustrates the velocity distributions on several cross-sections in the packed bed.The velocity increases in the area near the wall for the contraction of flow channel.

    The relative velocity distribution in the center region exhibits a more regular oscillation that is consistent with the radial porosity distribution.To assess the non-uniformity of the HTF,the standard deviation of relative velocity was adopted.A higher standard deviation complies with a less uniform velocity distribution.For the packed bed with a diameter ratio of 5,when the Reynolds number Repis 4998.72,the standard deviation of the relative velocity(sU)in the center region is 0.39.

    Fig.6.Radial porosity distribution and sub-surfaces along the radial direction.

    Fig.7.Relative velocity distribution along the radial direction.

    Fig.8.Flow velocity distribution inside the packed bed.

    where,N is the number of sub-surfaces,Uiis the relative velocity at the No.i sub-surface;ˉU is the average relative velocity of all the sub-surfaces.

    3.1.3.Radial temperature and liquid fraction distribution

    For the corresponding sub-surface,the radial temperature and liquid fraction of the PCMs are calculated as Eq.(33)and Eq.(34),whereΩpcmis the area of the PCMs on the sub-surface and Apcmis the area of the PCMs.

    Fig.9.Radial distributions of PCMs at different time.(a)Temperature;(b)liquid fraction.

    To illustrate the radial distribution variations of PCMs during the charging process,the temperature and liquid fraction of the PCMs on the sub-surfaces at different time are presented in Figs.9(a)and(b),fluctuating along the radial direction.When t=200 s,except for the near wall region,the PCMs’temperatures at most positions do not attain the melting point,and the liquid fractions remain 0.Thus,the charging process remains in the solidstate heating stage.When t=1200 s,the PCMs’temperatures are within the melting temperature range,and the liquid fraction differs significantly at each position.The melting temperature range of PCM is 395.1-413°C.Considering the enthalpy method,an increase temperature for every 1°C leads to an elevation of the PCM liquid fraction by 0.056.Accordingly,although the temperature differences of the PCMs are small,the liquid fractions variations are high.When t=2400 s,the PCMs’temperatures at most positions are greater than the final melting temperature,and the liquid fractions become 1 indicating that the charging process is in the liquid-state heating stage.

    In Fig.10,the liquid fractions of PCMs on the radial sub-surfaces are observed when t=600 s.The PCMs near the entrance melt faster,and the melting is initiated from the flow side of the HTF.Moreover,the PCMs melt faster on the sub-surfaces of dpand 2dpaway from the wall,complying with the liquid fraction distribution in Fig.9(b).For the low porosity on the sub-surfaces,that is 0.5dpand 1.5dpaway from the wall,the flow velocity is slow,while the convection heat transfer is weak.The proportions of the PCMs are higher on the sub-surfaces corresponding to 0.5dpand 1.5dpaway from wall;thus,the amount of heat required for melting increases leading to a longer charging process.

    3.2.Effect of diameter ratio on the flow and thermal performance

    In this section,the effects of various diameter ratios(4,5,6)on the flow and thermal performance of the packed bed LHS system are studied;at the inlet temperature Tinof 465°C,the mass flow rate is 260 kg·h-1.The inner diameter of the heat storage tank is 240 mm,and the height is 500 mm.To determine the number of PCM capsules in each packed bed,the heat storage of the three packed bed LHS systems was set to be nearly equal.Based on the aforementioned criterion,the parameters of the three packed beds are listed in Table 2.In Fig.11,the radial porosity distributions of the packed beds with different diameter ratios are compared,and it is observed that the radial porosity distributions oscillate but exhibit different trends.Over a distance of 1.5dpfrom the wall,the radial porosity distributions are almost constant,whereas,as the distance gets closer to the center,more significant differences are observed.For the packed bed with a diameter ratio of 4,the diameter can be filled with 4 capsules,and the capsules are in point contacts or no contact at the tank center,in which gaps are formed.Consequently,the porosity is large.For the packed bed with a diameter ratio of 5,the diameter can be filled with 5 capsules;thus,the tank center is filled with the capsules,and the central porosity is approximately zero.For the packed bed with a diameter ratio of 6,the central porosity ought to be as large as that of the packed bed with a diameter ratio of 4,theoretically.However,as the tank center is slightly far away from the tank wall,and the wall effect decreases,the packing of capsules at the tank center tends to be random,as shown in Fig.11.

    Fig.10.Liquid fractions of PCMs on the radial sub-surfaces when t=600 s.

    The flow relative velocity distributions inside the three packed bed LHS systems are illustrated in Fig.12(a).It is observed that the relative velocity distributions along the radial direction comply with the radial porosity distributions.At the near wall region,the maximum relative velocity increases with an increase in the diameter ratio,which is 2.87,3.49,and 3.99,respectively.With a decrease in the diameter of the PCMcapsule,gaps formed between the capsules and container wall are smaller.As shown in Eq.(31),when hot air flows through a narrower channel,the flow velocity will increase.For the center region,the standard deviation of relative velocity distribution in the three packed beds reaches 0.53,0.39,and 0.27,respectively,thus demonstrating that an increase of diameter ratio can decrease the velocity non-uniformity in the center region.

    The radial liquid fraction distributions of the packed beds when t=1200 s are shown in Fig.12(b).For the packed bed with a diameter ratio of 4,the PCMs near the wall and tank center melt faster than the other positions.However,for the packed bed with the diameter ratio of 5 and 6,the charging processes of PCMs near the wall are faster than those at the tank center,which is consistent with the radial porosity distributions.

    In Fig.13,the heat storage and average charging power of the packed bed LHS systems with different diameter ratios were calculated.With an increase in the diameter ratio,the charging times of the three systems are 121.03,89.73,and 79.80 min,respectively;thus,the average charging powers of the packed bed LHS systems increase.The capsule diameter in the system with a diameter ratio of 4 is 60 mm,while that with the diameter ratio of 5 or 6 achieves a diameter of 48 or 40 mm.A smaller diameter indicates a shorter heat transfer distance,whereas it will further result in a denser packing arrangement and promotes a decrease in pressure.The calculated results show that the decrease in pressure in the packed beds are 142.97,213.34,and 257.13 Pa,respectively.Hence,higher pump power is required to maintain the operation of the packed bed LHS system.Accordingly,when the packed bed LHS system is employed for practical use,the charging time,charging power,and pressure drop should be comprehensively considered.

    Fig.11.Radial porosity distributions of the packed bed models with different diameter ratios.

    Table 2 Parameters of the packed bed LHS systems with different diameter ratios.

    Fig.12.Radial distributions of the packed bed LHS systems with different diameter ratios.(a)Relative velocity;(b)liquid fraction when t=1200 s.

    Fig.13.Heat storage and the average charging power of the packed bed LHS systems with different diameter ratios.

    3.3.Effect of Reynolds number on the flow and thermal performance

    Three different Reynolds numbers(4229.69,4998.72,and 5767.75)are compared,which comply with inlet flow rates of 220,260,and 300 kg·h-1,respectively.The diameter ratio of packed bed is 5,and the Stefan number is 0.31.Fig.14(a)shows that the relative velocities of the HTF along the radial direction for different Reynolds numbers comply with a similar distribution.The maximum relative velocities of the HTF at the near wall region are 3.40,3.49,and 3.55,while the standard deviations of relative velocity at the center region are 0.40,0.39,and 0.39,demonstrating that an increase in the Reynolds number would result in a proportional improvement in velocity.Correspondingly,convection heat transfer in the packed bed is enhanced,and the charging process of the PCMs would accelerate(Fig.14(b)).The PCMs melt faster with an increase in the Reynolds number when t=1200 s.

    Fig.14.Flow and thermal properties in the packed bed LHS systems with different Reynolds numbers.(a)Radial relative velocity distribution;(b)heat storage and the average charging power;(c)radial liquid fraction distribution when t=1200 s.

    Fig.14(c)shows that an increase in the Reynolds number will not promote heat storage of the system,while the charging time can be shortened and charging power can increase accordingly.The charging time of the system with Rep=4229.69 is 97.53 min,while those of the systems with Rep=4998.72 and Rep=5767.75 are 89.73 and 83.63 min,and 8%and 14.25%faster,respectively.However,an increase in the Reynolds number would increase the dissipation of the HTF when flowing through capsules and improve the pressure drop at the inlet and outlet.The corresponding pressure drops in the three conditions reach 162.84,213.34,and 270.06 Pa,respectively.

    3.4.Effect of Stefan number on the flow and thermal performance

    Similar to the previous section,the diameter ratio of the packed bed is 5 and Reynolds number is 4998.72.Subsequently,the HTF inlet temperatures are set to 445,465,and 485°C,leading to Stefan numbers of 0.19,0.31,and 0.43,respectively.A comparison of the relative velocity distributions is drawn in Fig.15(a).It is observed that the increase in Stefan number slightly impacts the flow velocity and its radial distribution is almost constant.Hence,the convection effect between the HTF and PCM capsules is the same.However,an increase in the Stefan number can increase the temperature differences between the HTF and PCM capsules,which enhances the heat transfer to the PCMs.Thus,the charging processes of the PCMs is facilitated and PCMs melt faster,as presented in Fig.15(b).

    An increase in the Stefan number further indicates that the sensible heat stored in the PCMs increase,thus,the heat storage of the packed bed LHS system is improved(Fig.15(c)).When the Stefan number is 0.19,the charging time is 99.2 min,and the charging time becomes 89.73 and 84.7 min with an increase in the Stefan number to 0.31 and 0.43,respectively.Accordingly,an increase in the Stefan number can shorten the charging time while increasing the heat storage and average charging power.

    4.Conclusion

    In this study,a 3D packed bed LHS model was built,where the radial porosity oscillation caused by the wall effect was considered.The model was validated by comparing the temperature variations of PCM capsules with the experimental results.The impact of the radial porosity oscillation on the radial velocity of HTF,radial temperature of PCMs,and liquid fraction of PCMs were studied.Moreover,the effects of different dimensionless parameters(e.g.,diameter ratio,Reynolds number,and Stefan number)on the radial characteristics of the HTF and PCMs were discussed.The main conclusions are as follows:

    (1)The oscillating distribution of the radial porosity is analyzed by intercepting different cylindrical sub-surfaces along the radial direction.For the sub-surface where the capsules are in point contact with the container wall or inner capsules,the gaps formed result in high porosity.However,for the sub-surfaces that pass through the capsules,the PCMs occupy most of the area on the surface,and a lower porosity is observed.

    (2)The oscillating distribution of the radial porosity leads to a non-uniform distribution of the HTF velocity.The velocity would change sharply at the near wall region while an oscillating distribution at the center region is observed.The radial temperature distribution of the PCMs is consistent with the radial relative velocity distribution of the HTF.Hence,the PCMs melt faster at the radial positions with higher porosity.

    (3)With an increase in the diameter ratio,the maximum velocity at the near wall region increases.However,as the packing of the PCM capsules in the center of the tank becomes more random,the velocity non-uniformity in the center region decreases.Furthermore,the charging time of the packed bed LHS system decreases drastically and average charging power can be improved,whereas larger pressure drops are observed.

    (4)An increase in the Reynolds number would result in a proportional improvement of the velocity and accelerate the charging processes of the PCMs,which further shortens the charging time and increases the average charging power.This can cause larger pressure drops.

    Fig.15.Flow and thermal properties in the packed bed LHS systems with different Stefan numbers.(a)Radial relative velocity distribution;(b)radial liquid fraction distribution when t=1200 s;(c)heat storage and the average charging power.

    (5)An increase in the Stefan number slightly impacts the velocity distribution of the HTF,whereas it can enhance the temperature differences between the HTF and PCMs,and the sensible heat stored in PCMs.Thus,the charging time can be shortened and heat storage of the packed bed LHS system increases.

    Acknowledgements

    This work is supported by the Foundation for Innovative Research Groups of the National Natural Science Foundation of China(51521004)and the National Natural Science Foundation of China(51906150).

    Compliance with ethics guidelines

    H.B.Liu and C.Y.Zhao declare that they have no conflict of interest or financial conflicts to disclose.

    Nomenclature

    V volume(m3)

    A area(m2)

    dpdiameter of PCM capsule(m)

    R radius of container(m)

    D diameter of container(m)

    T temperature(°C)

    H specific enthalpy(kJ·kg-1)

    cpspecific heat capacity(J·kg-1·K-1)

    L latent heat(kJ·kg-1)

    E energy(J)

    Ste Stefan number

    q mass flow rate(kg·s-1)

    u velocity(m·s-1)

    u′fluctuating velocity(m·s-1)

    U relative velocity

    Reppore Reynolds number

    t time(s)

    x distance(m)

    z height(m)

    p pressure(Pa)

    k turbulent kinetic energy(m2·s-2)

    c1,c2model constants in turbulent kinetic energy equation

    c model constant in turbulent viscosity correlation

    Q heat storage(J)

    P charging power(W)

    s standard deviation of relative velocity

    φ porosity

    ρ density(kg·m-3)

    λ thermal conductivity(W·m-1·K-1)

    μ viscosity(kg·m-1·s-1)

    μtturbulent viscosity(kg·m-1·s-1)

    σTPrandtl number in energy equation

    σkPrandtl number in turbulent kinetic energy equation

    σεPrandtl number in turbulence dissipation rate equation

    ε turbulence dissipation rate(m2·s-3)

    β liquid fraction

    γ RNG k-εturbulence model coefficient

    Ω area of HTF and PCM

    Subscripts

    i,j coordinate direction

    f heat transfer fluid

    PCM phase change material

    total total amount

    shell shell of PCM capsule

    tank tank of container

    stored stored heat

    charge charging process

    aver average value

    in inlet

    0 the initial state

    l liquid state

    s solid state

    人人澡人人妻人| 91成人精品电影| 亚洲性久久影院| av网站免费在线观看视频| 少妇人妻久久综合中文| 啦啦啦在线观看免费高清www| 精品少妇内射三级| 成年av动漫网址| 99国产精品免费福利视频| 男的添女的下面高潮视频| av免费在线看不卡| 日韩免费高清中文字幕av| 亚洲av综合色区一区| 国产淫语在线视频| 国产成人精品无人区| 日韩大片免费观看网站| 嫩草影院入口| 精品亚洲成a人片在线观看| 久久鲁丝午夜福利片| 男女啪啪激烈高潮av片| 国产精品三级大全| 最近2019中文字幕mv第一页| 久久久精品免费免费高清| 久久久国产精品麻豆| 成人黄色视频免费在线看| 成人漫画全彩无遮挡| 狂野欧美白嫩少妇大欣赏| 国产免费又黄又爽又色| 久久精品国产亚洲av涩爱| 永久免费av网站大全| 18禁裸乳无遮挡动漫免费视频| 大片免费播放器 马上看| 中文字幕亚洲精品专区| 久久国产精品大桥未久av | 高清欧美精品videossex| 九九爱精品视频在线观看| 精品人妻熟女毛片av久久网站| 蜜臀久久99精品久久宅男| a级毛片免费高清观看在线播放| 国产精品久久久久久精品电影小说| 又黄又爽又刺激的免费视频.| 欧美最新免费一区二区三区| 久久久久精品久久久久真实原创| 男男h啪啪无遮挡| 极品教师在线视频| 国产男女超爽视频在线观看| 精品久久久久久电影网| 日韩免费高清中文字幕av| 丁香六月天网| 大话2 男鬼变身卡| 又粗又硬又长又爽又黄的视频| 美女cb高潮喷水在线观看| 亚洲精品一区蜜桃| 国产白丝娇喘喷水9色精品| 五月伊人婷婷丁香| 国产高清有码在线观看视频| 国产在线一区二区三区精| 国产白丝娇喘喷水9色精品| 又大又黄又爽视频免费| 最近手机中文字幕大全| 97在线人人人人妻| 精品国产一区二区三区久久久樱花| 成人国产麻豆网| 男男h啪啪无遮挡| 亚洲不卡免费看| 国产精品.久久久| 精品久久久久久久久亚洲| 久久99热这里只频精品6学生| 亚洲av综合色区一区| 不卡视频在线观看欧美| 日本wwww免费看| 国产精品久久久久久av不卡| 免费黄频网站在线观看国产| 男的添女的下面高潮视频| 色哟哟·www| 国产精品一区www在线观看| 80岁老熟妇乱子伦牲交| 久久久亚洲精品成人影院| 多毛熟女@视频| 亚洲欧美一区二区三区国产| 国产精品99久久久久久久久| 国产 一区精品| 日韩成人av中文字幕在线观看| 免费看av在线观看网站| 日韩,欧美,国产一区二区三区| 国产真实伦视频高清在线观看| 天美传媒精品一区二区| 在线亚洲精品国产二区图片欧美 | av免费观看日本| av国产精品久久久久影院| 亚洲国产成人一精品久久久| 色94色欧美一区二区| 亚洲国产精品专区欧美| 亚洲国产色片| 亚洲人成网站在线观看播放| 少妇丰满av| 嫩草影院入口| 国产欧美日韩精品一区二区| 黄片无遮挡物在线观看| 好男人视频免费观看在线| 又爽又黄a免费视频| 日本免费在线观看一区| 激情五月婷婷亚洲| 久久久久久久大尺度免费视频| 少妇高潮的动态图| 亚洲av成人精品一二三区| 男男h啪啪无遮挡| 国产亚洲av片在线观看秒播厂| 一级毛片黄色毛片免费观看视频| 2018国产大陆天天弄谢| 亚洲欧美日韩另类电影网站| 看免费成人av毛片| 七月丁香在线播放| 少妇人妻精品综合一区二区| 在线播放无遮挡| 国产精品福利在线免费观看| 亚洲成人手机| 丝袜喷水一区| 少妇高潮的动态图| 大码成人一级视频| 国产老妇伦熟女老妇高清| 99久久精品一区二区三区| 亚洲精品乱码久久久久久按摩| 亚洲精品亚洲一区二区| 亚洲色图综合在线观看| 伊人久久国产一区二区| 在线看a的网站| 一级毛片 在线播放| 亚洲综合精品二区| 亚洲国产毛片av蜜桃av| 夜夜看夜夜爽夜夜摸| 国产国拍精品亚洲av在线观看| 亚洲国产色片| 国产成人精品福利久久| 免费久久久久久久精品成人欧美视频 | 亚洲情色 制服丝袜| 性色av一级| 大香蕉久久网| 日本vs欧美在线观看视频 | 欧美精品一区二区大全| 成人影院久久| 亚洲图色成人| 欧美xxⅹ黑人| 在线观看美女被高潮喷水网站| 熟妇人妻不卡中文字幕| 精品国产一区二区三区久久久樱花| 自拍欧美九色日韩亚洲蝌蚪91 | 观看av在线不卡| 18禁在线无遮挡免费观看视频| 高清在线视频一区二区三区| 成人毛片60女人毛片免费| 国产亚洲最大av| 丝袜在线中文字幕| 插逼视频在线观看| 青春草亚洲视频在线观看| 超碰97精品在线观看| 久久精品国产亚洲av涩爱| 精品久久久久久电影网| 欧美xxxx性猛交bbbb| 内射极品少妇av片p| 人人妻人人澡人人看| 大香蕉97超碰在线| 国产av一区二区精品久久| 亚洲欧美精品专区久久| 热99国产精品久久久久久7| 亚洲人成网站在线观看播放| 视频区图区小说| 久久国内精品自在自线图片| 99国产精品免费福利视频| 亚洲国产av新网站| 欧美xxxx性猛交bbbb| 成年人免费黄色播放视频 | 青春草国产在线视频| 伦理电影大哥的女人| 国产一区亚洲一区在线观看| 亚洲国产欧美日韩在线播放 | 国产片特级美女逼逼视频| 日日摸夜夜添夜夜添av毛片| 精品久久久久久电影网| 99热6这里只有精品| 岛国毛片在线播放| 久久久久久伊人网av| a 毛片基地| 日韩,欧美,国产一区二区三区| 99久久中文字幕三级久久日本| 亚洲欧美中文字幕日韩二区| 亚洲久久久国产精品| 成人国产av品久久久| 欧美xxxx性猛交bbbb| 内射极品少妇av片p| 女人精品久久久久毛片| 狂野欧美激情性bbbbbb| 在线播放无遮挡| 免费播放大片免费观看视频在线观看| 亚洲人与动物交配视频| 99精国产麻豆久久婷婷| 人人妻人人添人人爽欧美一区卜| 国产日韩欧美在线精品| 日本与韩国留学比较| 黄色毛片三级朝国网站 | 黑人猛操日本美女一级片| 国产精品一二三区在线看| 亚洲,欧美,日韩| 国产成人91sexporn| 精品人妻熟女av久视频| 国产成人精品福利久久| 亚洲欧洲精品一区二区精品久久久 | a级毛片免费高清观看在线播放| 极品教师在线视频| 亚洲一区二区三区欧美精品| 婷婷色综合www| 欧美老熟妇乱子伦牲交| 亚洲国产精品专区欧美| 国产一级毛片在线| 我的老师免费观看完整版| 成人18禁高潮啪啪吃奶动态图 | 男人添女人高潮全过程视频| 少妇人妻 视频| 亚洲欧美中文字幕日韩二区| 成人亚洲欧美一区二区av| 97在线视频观看| 纯流量卡能插随身wifi吗| 制服丝袜香蕉在线| 大片免费播放器 马上看| 亚洲人成网站在线观看播放| 18禁裸乳无遮挡动漫免费视频| 妹子高潮喷水视频| 亚洲精品乱码久久久久久按摩| 亚洲欧美日韩卡通动漫| 中文资源天堂在线| 各种免费的搞黄视频| 国产无遮挡羞羞视频在线观看| 寂寞人妻少妇视频99o| 少妇裸体淫交视频免费看高清| 丁香六月天网| 国产精品无大码| 国产精品久久久久成人av| 色5月婷婷丁香| 丰满人妻一区二区三区视频av| 国产69精品久久久久777片| av福利片在线观看| 亚洲,一卡二卡三卡| 亚洲丝袜综合中文字幕| 久久久久精品久久久久真实原创| a 毛片基地| 亚洲av二区三区四区| 一级毛片电影观看| 午夜影院在线不卡| 精品久久久久久电影网| 国产日韩欧美视频二区| 大陆偷拍与自拍| 人人澡人人妻人| 欧美97在线视频| 国产熟女午夜一区二区三区 | 精品久久久噜噜| 婷婷色综合大香蕉| 国产69精品久久久久777片| 亚洲成人av在线免费| 美女xxoo啪啪120秒动态图| 一级爰片在线观看| 色婷婷av一区二区三区视频| 亚洲国产精品一区二区三区在线| 国产一区二区在线观看av| 亚洲va在线va天堂va国产| 国产欧美日韩精品一区二区| 国产真实伦视频高清在线观看| 亚洲国产最新在线播放| 成年女人在线观看亚洲视频| 久久99蜜桃精品久久| 在线观看美女被高潮喷水网站| 国产日韩欧美亚洲二区| 伊人亚洲综合成人网| 国精品久久久久久国模美| 老熟女久久久| 男的添女的下面高潮视频| 日本爱情动作片www.在线观看| 精品99又大又爽又粗少妇毛片| 成年人免费黄色播放视频 | 亚洲高清免费不卡视频| 日韩精品免费视频一区二区三区 | 国产黄片美女视频| 青春草亚洲视频在线观看| 亚洲av成人精品一区久久| 又大又黄又爽视频免费| 特大巨黑吊av在线直播| 久久99精品国语久久久| av女优亚洲男人天堂| 精品酒店卫生间| 青春草亚洲视频在线观看| 中文字幕久久专区| 成人国产av品久久久| 欧美成人午夜免费资源| 国产真实伦视频高清在线观看| 熟女人妻精品中文字幕| 亚洲欧美精品自产自拍| 国产在视频线精品| 大陆偷拍与自拍| 啦啦啦中文免费视频观看日本| 人人妻人人看人人澡| 精品午夜福利在线看| 熟妇人妻不卡中文字幕| 日韩人妻高清精品专区| 丝袜脚勾引网站| 免费观看av网站的网址| 六月丁香七月| 久久久久国产网址| 日本av手机在线免费观看| 午夜免费男女啪啪视频观看| 黑人猛操日本美女一级片| 国产色爽女视频免费观看| 久久精品久久久久久噜噜老黄| 韩国高清视频一区二区三区| 91久久精品电影网| 国产淫片久久久久久久久| .国产精品久久| 校园人妻丝袜中文字幕| 99re6热这里在线精品视频| 男女边吃奶边做爰视频| 久久久精品免费免费高清| 五月玫瑰六月丁香| tube8黄色片| 搡女人真爽免费视频火全软件| 成人二区视频| 中文欧美无线码| 久久久久国产精品人妻一区二区| 91久久精品电影网| 国产片特级美女逼逼视频| 全区人妻精品视频| 中文字幕精品免费在线观看视频 | 人人妻人人看人人澡| 一边亲一边摸免费视频| 久久婷婷青草| 日韩av免费高清视频| 国产成人精品一,二区| 国产伦精品一区二区三区视频9| 欧美 日韩 精品 国产| 汤姆久久久久久久影院中文字幕| 麻豆乱淫一区二区| 男人舔奶头视频| 下体分泌物呈黄色| 久久免费观看电影| 欧美xxxx性猛交bbbb| 日韩欧美 国产精品| 春色校园在线视频观看| 天堂中文最新版在线下载| 久久 成人 亚洲| 尾随美女入室| 观看美女的网站| 色5月婷婷丁香| 成人影院久久| 人人妻人人爽人人添夜夜欢视频 | 国产成人精品一,二区| 十分钟在线观看高清视频www | 亚洲精品日韩av片在线观看| 亚洲av综合色区一区| 蜜臀久久99精品久久宅男| 亚洲欧美日韩另类电影网站| av在线观看视频网站免费| 少妇熟女欧美另类| 日日啪夜夜撸| 韩国av在线不卡| a级毛片在线看网站| 麻豆成人午夜福利视频| 22中文网久久字幕| av女优亚洲男人天堂| 国产成人精品无人区| 中文乱码字字幕精品一区二区三区| 99视频精品全部免费 在线| 国产成人a∨麻豆精品| 亚洲欧美日韩东京热| 久久青草综合色| 在线精品无人区一区二区三| 观看美女的网站| 在线精品无人区一区二区三| 日日爽夜夜爽网站| 国产男女超爽视频在线观看| 日韩免费高清中文字幕av| 高清午夜精品一区二区三区| 亚洲四区av| 九九在线视频观看精品| 亚洲精品国产av成人精品| 一级毛片黄色毛片免费观看视频| 在现免费观看毛片| 亚洲精品亚洲一区二区| 丰满饥渴人妻一区二区三| 久久99蜜桃精品久久| 另类精品久久| 又大又黄又爽视频免费| 亚洲精品一区蜜桃| 我的女老师完整版在线观看| 中文字幕久久专区| 黑人猛操日本美女一级片| 国产在线男女| 搡老乐熟女国产| 亚洲性久久影院| 日本av手机在线免费观看| 黄色配什么色好看| 王馨瑶露胸无遮挡在线观看| 制服丝袜香蕉在线| 亚洲欧美一区二区三区国产| 久久久午夜欧美精品| 建设人人有责人人尽责人人享有的| 这个男人来自地球电影免费观看 | 另类亚洲欧美激情| 天堂8中文在线网| 午夜视频国产福利| 少妇裸体淫交视频免费看高清| 伊人久久精品亚洲午夜| 99久久综合免费| 色网站视频免费| 中文字幕精品免费在线观看视频 | 国产白丝娇喘喷水9色精品| 热re99久久国产66热| 亚洲成色77777| 能在线免费看毛片的网站| 免费观看av网站的网址| 曰老女人黄片| 国产亚洲5aaaaa淫片| 免费看av在线观看网站| 亚洲欧洲国产日韩| 又黄又爽又刺激的免费视频.| 自拍欧美九色日韩亚洲蝌蚪91 | 国产综合精华液| 日韩,欧美,国产一区二区三区| 卡戴珊不雅视频在线播放| 国产白丝娇喘喷水9色精品| 丰满少妇做爰视频| 成年美女黄网站色视频大全免费 | 日韩一区二区视频免费看| 亚洲精品国产色婷婷电影| 亚洲性久久影院| 日韩不卡一区二区三区视频在线| 国产乱人偷精品视频| 各种免费的搞黄视频| 一级毛片aaaaaa免费看小| 狂野欧美激情性xxxx在线观看| 18+在线观看网站| 欧美日韩视频高清一区二区三区二| 免费大片黄手机在线观看| 久久午夜福利片| 午夜福利在线观看免费完整高清在| 男人和女人高潮做爰伦理| 欧美三级亚洲精品| 夜夜骑夜夜射夜夜干| 国产免费一级a男人的天堂| 少妇 在线观看| 亚洲综合色惰| 久久午夜福利片| 亚洲第一区二区三区不卡| 国产 精品1| 久久99蜜桃精品久久| 纯流量卡能插随身wifi吗| 日日撸夜夜添| 欧美+日韩+精品| 日韩成人av中文字幕在线观看| 久久女婷五月综合色啪小说| 在现免费观看毛片| 国产亚洲5aaaaa淫片| 特大巨黑吊av在线直播| 一本—道久久a久久精品蜜桃钙片| 日韩一区二区视频免费看| 亚洲情色 制服丝袜| 亚洲中文av在线| 婷婷色综合www| 日韩一本色道免费dvd| 免费看av在线观看网站| 亚洲综合色惰| 两个人免费观看高清视频 | 80岁老熟妇乱子伦牲交| 日日爽夜夜爽网站| 久久人妻熟女aⅴ| 精品人妻熟女av久视频| 一区在线观看完整版| 99久久人妻综合| 啦啦啦在线观看免费高清www| 日韩在线高清观看一区二区三区| 91精品伊人久久大香线蕉| 成年美女黄网站色视频大全免费 | 97超碰精品成人国产| 一区二区三区四区激情视频| 国产熟女欧美一区二区| 在线免费观看不下载黄p国产| 3wmmmm亚洲av在线观看| 亚洲第一av免费看| 观看美女的网站| 最近中文字幕2019免费版| 中文乱码字字幕精品一区二区三区| 国产亚洲5aaaaa淫片| 国产深夜福利视频在线观看| 国产亚洲欧美精品永久| 精品少妇黑人巨大在线播放| 美女cb高潮喷水在线观看| 久久久午夜欧美精品| 成年av动漫网址| 欧美精品国产亚洲| 久久av网站| videossex国产| 两个人的视频大全免费| 久久精品久久久久久久性| 亚洲欧美成人综合另类久久久| 精品久久久噜噜| 久久人人爽av亚洲精品天堂| 99视频精品全部免费 在线| 欧美xxxx性猛交bbbb| 最近中文字幕高清免费大全6| 亚洲欧美成人综合另类久久久| 国产视频首页在线观看| 亚洲中文av在线| 精品一区二区三区视频在线| 哪个播放器可以免费观看大片| 久久久a久久爽久久v久久| 在线观看国产h片| 国产精品一区二区在线观看99| 乱系列少妇在线播放| 国产一区二区在线观看日韩| 日韩电影二区| 日日摸夜夜添夜夜爱| 边亲边吃奶的免费视频| 三上悠亚av全集在线观看 | 亚洲精华国产精华液的使用体验| 亚洲婷婷狠狠爱综合网| 亚洲精品自拍成人| 伊人久久国产一区二区| av天堂久久9| 大陆偷拍与自拍| 免费在线观看成人毛片| 十分钟在线观看高清视频www | 九色成人免费人妻av| 黄色配什么色好看| 人妻少妇偷人精品九色| 精品国产乱码久久久久久小说| 国产淫语在线视频| 久久 成人 亚洲| 22中文网久久字幕| av线在线观看网站| av天堂久久9| 国产av码专区亚洲av| 国产精品一二三区在线看| 色吧在线观看| 99re6热这里在线精品视频| av不卡在线播放| 一区在线观看完整版| 精品一区二区三卡| 日日摸夜夜添夜夜爱| 男的添女的下面高潮视频| 免费黄网站久久成人精品| 国产男女内射视频| 国产黄片美女视频| 久久精品国产亚洲网站| 国产精品国产三级国产av玫瑰| 高清视频免费观看一区二区| 久久国内精品自在自线图片| 日本色播在线视频| 成年女人在线观看亚洲视频| av卡一久久| 蜜臀久久99精品久久宅男| 日本91视频免费播放| 色视频在线一区二区三区| 久久99一区二区三区| 国产91av在线免费观看| 久久久精品94久久精品| 人体艺术视频欧美日本| 超碰97精品在线观看| 又粗又硬又长又爽又黄的视频| 国产精品一区二区三区四区免费观看| 黄色视频在线播放观看不卡| 菩萨蛮人人尽说江南好唐韦庄| 日韩成人伦理影院| 在现免费观看毛片| 视频区图区小说| 亚洲精品一二三| 2021少妇久久久久久久久久久| 美女主播在线视频| av视频免费观看在线观看| 最后的刺客免费高清国语| 国产精品无大码| 国产淫语在线视频| 精品卡一卡二卡四卡免费| 如日韩欧美国产精品一区二区三区 | 免费久久久久久久精品成人欧美视频 | 最近最新中文字幕免费大全7| 国产精品不卡视频一区二区| 春色校园在线视频观看| 婷婷色综合大香蕉| 欧美变态另类bdsm刘玥| 春色校园在线视频观看| 国产一区二区在线观看日韩| 少妇被粗大的猛进出69影院 | 七月丁香在线播放| 十八禁网站网址无遮挡 | 美女主播在线视频| 少妇被粗大的猛进出69影院 | 最近最新中文字幕免费大全7| 免费看光身美女| 欧美xxxx性猛交bbbb| 我的女老师完整版在线观看| 99re6热这里在线精品视频| a 毛片基地| 国产伦精品一区二区三区视频9| 日日摸夜夜添夜夜添av毛片| 国内精品宾馆在线| 搡女人真爽免费视频火全软件| 中国美白少妇内射xxxbb| 女性被躁到高潮视频| 亚洲av成人精品一二三区| 国产成人一区二区在线| 麻豆成人午夜福利视频| 色婷婷av一区二区三区视频| 亚洲在久久综合| 51国产日韩欧美| 亚洲av男天堂| 成年人免费黄色播放视频 | 亚洲成人av在线免费| 嘟嘟电影网在线观看| av一本久久久久|