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

    Turbulent characteristics and rotation correction of wall function in rotating channel with high local rotation parameter

    2018-10-15 02:44:06ZhiTAOHuijieWURuquanYOUHaiwangLIKuanWEI
    CHINESE JOURNAL OF AERONAUTICS 2018年10期

    Zhi TAO,Huijie WU,Ruquan YOU,Haiwang LI,*,Kuan WEI

    aNational Key Laboratory of Science and Technology on Aero-Engine Aero-Thermodynamics,Beihang University,Beijing 100083,China

    bCollaborative Innovation Center for Advanced Aero-Engine,Beihang University,Beijing 100083,China

    KEYWORDS Coriolis force;Entry section;Relaminarization;Rotating channel;Turbulence

    Abstract The turbulent fluctuation and the rotation correction of wall function law are investigated in the entrance section of a rotating channel.The one-dimensional hot wire probe and the X-type probe are utilized to measure the boundary layer at four streamwise stations.Through the analysis on the boundary layer near the leading side and trailing side,it is found that the turbulent fluctuation is promoted in the trailing side whereas suppressed in the leading side.This difference is attributed to the Coriolis instability near the trailing side.In addition,considering the local rotation parameter Rc,whose maximum absolute value is 0.014,is larger than that in previous research,whose maximum value is 0.007,the whole process of the relaminarization is captured.To understand this phenomenon better,the effects of the generation term and the Coriolis term in the transport equation of the Reynolds stress are discussed.In addition,the rotation correction of the viscous-Coriolis region and the Coriolis region are discussed,a new revising method for the wall function is proposed.

    1.Introduction

    As we all know,the turbine blades are operated in the high temperature environment to improve the thermal efficiency of modern gas turbine engine.Although the melting point of the turbine blade material increases year by year,this extremely high temperature is far beyond the material’s tolerance.Accordingly,lots of cooling techniques are applied to protect the turbine blade from the hot gases.Internal cooling technique is one of the most prevalent methods used to protect the turbine blade.For instance,serpentine passage,as an effective method,was developed in the middle section of a turbine blade.1,2To improve the performance of serpentine passages,most of the previous researches focused on the heat transfer instead of the flow behavior in the serpentine passages.However,without knowledge about the flow behavior under rotating condition and the interaction between heat transfer and fluid flow,the principle behind the heat transfer phenomenon may not be understood comprehensively.Therefore,some researchers investigated the flow behavior in rotating ducts.These researches can be divided into five branches:the effects of rotation on the primary velocity profiles in a straight channel without ribs;3–5the effects of rotation on the secondary flow fields in a straight channel without ribs;6–8the effects of rotation on the turbulent fluctuation;9,10the effects of rotation on the flow fields in a channel with complex geometry,such as U turn and ribs;11–15the effects of rotation on the flow fields in a channel with heated walls.16–18Nevertheless,some important issues are not understood clearly,such as the rotation correction of the wall function and the relaminarization process discovered in the analysis of the turbulent fluctuation.Considering both turbulent fluctuation and the rotation correction of the wall function have relationship with the primary velocity profiles,the present work also investigates the primary flow in a smooth rotating duct.Therefore,the rotation correction of the wall function and the turbulent fluctuation are the main concern of the present paper.

    Koyama et al.19measured the primary velocity profiles on the mid-span of the duct with aspect ratio equaling 7.A universal rotation correction of the log law formula is obtained after the analysis of the primary flow of the turbulent boundary layer.There are two correction methods of the standard log law formula.The first one only changes the empirical coefficients of the log law formula without changing the basic form,the second one adds an additional term to the basic form of the standard log law formula.

    Watmuff et al.3measured the primary flow on the mid-span of the duct with aspect ratio equaling 4.After analyzing primary velocity profile at different regions of the turbulent boundary layer,the viscous region is not influenced by rotation,whereas the logarithmic region needs rotation correction.They only changed the empirical coefficients of the log law formula instead of changing the basic form of log law.

    Nakabayashi and Kitoh4measured the fully developed flow of a straight rotating duct with aspect ratio equaling 8.They mainly focused on the effects of rotation on the flow under low Reynolds number.Using y/L,L/δν,δν/δcas dimensionless parameters,they divided the turbulent boundary layer into several regions:viscous region,buffer region,logarithmic region,Coriolis region and core region.Here,y means the normal distance from the wall,δν(δν= ν/Uτ)means the viscous length scale,Uτmeans the friction velocity,ν is the kinematic viscosity coefficient,δc(δc=Uτ/|Ω|)means the Coriolis length scale and L refers to the channel half-width,Ω means the angular velocity of the rotating duct.In the process of the rotation correction of the log law formula,they changed the empirical coefficients of the log law formula instead of changing the basic log law form as rotation correction.However,the expression of the intercept C,namely the intercept of the log law formula,is so complicated,influenced not only by Ων/U2τ,but also by Reynolds number,that its correction formula is not given.

    Nickels and Joubert20obtained the basic form of log law formula with rotation correction through dimensional analysis.The corrected formula is valid for different streamwise pressure gradients and the intensities of secondary flow.

    Apart from the rotation correction of the wall function,the turbulent fluctuation is also important,because the rotation correction of the wall function law has much relationship with the turbulent fluctuation.For example,by contrasting the turbulent fluctuation in the leading side and the trailing side,we can see whether the so-called relaminarization occurs,which is adverse to the rotation correction of the wall function.Besides,the turbulent fluctuation profiles,rather than the mean primary profiles,are used to judge the penetration of the Coriolis force,10which is essential for the division of the turbulent boundary layer under rotating conditions,thus influencing the rotation correction of the wall function.Therefore,it is necessary to pay attention to the turbulent fluctuation.

    Koyama et al.19measured the fluctuation of the turbulent boundary layer at a fixed streamwise location in the rotating straight duct with aspect ratio equaling 7.The experiment results show that the transition of the boundary layer near the leading side is found to be suppressed,while promoted on the trailing side.

    Watmuff et al.3measured the distribution of all components of the Reynolds stress of the turbulent boundary layer in the mid-span of the rotating straight duct with aspect ratio equaling 4.The experimental results suggest that the component of Reynolds stress is affected by the rotation only outside the viscous region.〈uv〉(u means the streamwise fluctuating velocity,v means the normal fluctuating velocity)remains to be nearly constant in the logarithmic region,though affected by the rotation.They also analyzed the streamwise energy spectrum at the same dimensionless distance under different rotating conditions.The results show that only low wavenumber spectrum is affected by the rotation.

    Macfarlane and Joubert21investigated the distribution of all components of the Reynolds stress of the turbulent boundary layer in the mid-span of the rotating straight duct with aspect ratio equaling to 4,2,1,respectively.The results show that the Reynolds stress is affected by rotation when y+>10,which is different from the mean primary velocity profile that is not influenced by rotation until y+>40.They investigated the influence of the secondary flow on the Reynolds stress by changing aspect ratio.It turns out that the distribution of the Reynolds stress at downstream locations of the ducts with aspect ratio equaling 4 is almost the same with that at downstream locations of the duct with aspect ratio equaling 2,whereas obviously different from that at downstream locations of the duct with aspect ratio equaling 1.

    Macfarlane et al.8also investigated the interaction between the secondary flow and the turbulent fluctuation of the turbulent boundary layer in the rotating straight duct with aspect ratio equaling 4,2,1,respectively.They contrasted the streamwise energy spectra at a certain location of the boundary layer in the ducts with different aspect ratios.Spectral results for aspect ratio equaling 2 show almost no variation at the high wavenumber,whereas spectra results for aspect ratio equaling 1 show that the whole spectrum shift towards low wavenumber on the trailing side and high wavenumber on the leading side.

    Nakabayashi and Kitoh10measured the distribution of the parameters of turbulent flow in the mid-span of rotating channel with aspect ratio equaling 8.Using dimensional analysis,they divided the turbulent boundary layer near the wall into several regions such as viscous region,buffer region,Coriolis region,etc.And they find the ranges of these regions are functions of Re′(Re′=UL/ν)and|Ω|ν/.Here,U means the primary velocity.Through measuring 〈uv〉 and performing quadrant analysis,they found the strong ejection which occurs periodically on the trailing side even though the ejection frequency is low.In the energy spectrum of the turbulent fluctuation,a rise in the spectrum is seen in the corresponding ejection frequency.In this way,it is demonstrated that the increase of the ejection is responsible for the increase of 〈 uv〉near the trailing side.

    In a word,the researches mentioned above have concerned the fundamental issues such as the rotation correction of the wall function law and the turbulent fluctuation,but those researches present some limitations.For example,the local rotation parameteris low in the researches above,the maximum value of which is 0.007.Low local rotation parameter often narrows down the applicable range of the revised wall function law and limits the occurrence of new phenomena of turbulent fluctuation,such as the relaminarization.In addition,the researches mentioned above focuses on the rotation correction of logarithmic region,but the intercept of the revised log law cannot be obtained,thus the log law corrected by rotation is of no use.

    Thus,the authors investigate the turbulent flow in the entrance section of the straight duct without ribs under rotating conditions.The present work focuses on the turbulent fluctuation to reveal some new phenomena because of the high local rotation parameter,such as the relaminarization process,which is adverse for us to correct the wall function near the leading side.Besides,the previous work only focuses on rotation correction on the log law region,which would result in a problem that the rotation correction of the intercept C cannot be obtained.Thus,after figuring out that there is no relaminarization,a new method to correct the wall function of the whole inner layer of the turbulent boundary layer is developed.

    2.Apparatus,techniques and experimental conditions

    2.1.Rotating facility and test section

    The experimental apparatus and methods utilized are the same with the ones in Wei et al.22The channel geometry and coordinate system are presented in Fig.1.In Fig.1,there are tripping wires attached both on the trailing side and leading side,which are 70 mm from the inlet of the channel.The measuring stations are the same with those in Wei et al.,22but the flow is measured by the X-type probe,which can obtain the twodimensional flow.It is different from that in Wei et al.,22which only obtain one-dimensional flow.

    2.2.Hotwire technique

    2.2.1.X-type hot wire

    The two-dimensional flow field is measured with the X-type hot wire anemometer,the Dantec 55P61 probe.But this hot wire is unable to measure the two-dimensional flows near the wall,owing to its structure.Considering there is no twodimensional boundary layer probe,the authors have to change the support of the Dantec 55H25.To change the support of the two-dimensional probe,three requirements need to be guaranteed.First,the support of the probe is located behind the probe when measuring the velocity,so as to minimize the influence of the support on the flow;second,the change of the support has to be made to guarantee that the sensing element of the probe can be as close to the wall as possible,so as to measure the viscous region of the boundary layer;third,the direction of the incoming flow direction must be within the acceptance angle of the probe.To satisfy the three requirements,the authors add a long support with a turning angle equaling 93°to the short support of the Dantec 55H25,which is shown in Fig.2.The measuring locations of the twodimensional hot wire are the same with the one-dimensional hot wire.

    Fig.1 Channel geometry and coordinate system.

    Fig.2 Installation structure of changed X-type hot wire.

    In addition,the measurement of the wall shear stress,the distance from the hotwire to the wall for the one dimensional hot wire and the measurement of the wall shear stress for X-type probe are presented by Wei et al.22As for the distance from the X-type hot wire to the wall,the velocity measured by the X-type probe is interpolated to the velocity profile measured by the one-dimensional hot wire at each station,then the distance from the X-type probe to the wall and the wall shear stress τwis obtained at the same time.Finally,the friction velocity is obtained through the equationHere,ρ means the density of the gas.

    2.2.2.Calibration of X-type hot wire

    For the X-type hot wire,the calibrating process is more complicated than the one-dimensional process.The absolute value of the velocity and the flow angle should be obtained at the same time from the voltages measured by the X-type probe.To do this,the X-type should be installed on a device,through which the angle of the X-type probe can be changed,as shown in Fig.3.

    The calibrating curve of the X-type hot wire is presented in Fig.4,the abscissa and the ordinate represent the voltage obtained by the two probes of the X-type hot wire.Different calibrating curve corresponds to different flow angle θ,and different point in the same curve corresponds to different velocity.

    Look-up matrix method is utilized to process the voltage through the calibrating curve listed in Fig.4.Then the relationship between Ut–θ coordinate system(Utis the absolute value of total velocity)and E1–E2coordinate system is obtained.In this way,the velocity and the flow angle can be obtained through the voltage E1,E2acquired by the X-type probe respectively.The concrete process of the look-up process is presented by Lueptow et al.23

    Fig.3 Holder of X-type hot wire.

    Fig.4 Calibrating curves of X-type hot wire.

    2.3.Experimental conditions

    The governing parameters for the boundary layer in the rotating channel of this paper are Reynolds number Re(Re=UD/ν)(D is the hydraulic diameter of the rotating channel)and rotation number Ro(Ro=ΩD/U).The local Reynolds number Reδ(Reδ=Uτδ/ν), δ is the boundary layer thickness and the local rotation parameter Rc are useful in the boundary layer analysis.The running parameters are presented in Table 1.Since the high Reynolds number will make the boundary layer so thin that we cannot measure the velocity within the boundary layer;whereas there is no logarithmic law within the boundary layer when the Reynolds number are small.Thus,we take Reynolds number equal 19000.

    In the first column of Table 1,the number 1,2,3,4 are the streamwise stations in the channel,and the number 30,60,90,120,150 are the rotating speed Ω,whose unit is r/min.The letter ‘‘T”refers to the trailing side whereas the letter ‘‘L” means the leading side.The word ‘‘stationary” means the rotating channel does not rotate.

    2.4.Discussion of uncertainties

    For any experiment,errors are inevitable in all components of the measurement system.For the X-type probe,the primary velocity measured by the X-type probe is obtained not through transforming the measured voltage into the velocity directly,but through the look-up table,23which firstly transforms the voltage into Utand θ, then utilizes the formula U=Utcosθ,V=Utsinθ to obtain the primary velocity and normal velocity.Therefore,the measurement errors of the velocity using X-type probe come from not only the calibration process and the variation of the temperature,but also the interpolation error and fitting error of look-up table.After the temperature correction,the calibration error and the error caused by the temperature variation only influence the absolute value of Ut.The fitting error and the interpolation error of the look-up table influence the value of Utand the angle of the flow direction.According to the error analysis of theone-dimensional hot wire,the sum of the calibration error and the error caused by the temperature variation is approximately±2%±0.02 m/s,both of which only influence the value of Ut.Look-up table includes two fitting processes and one interpolation process,and the error of E1and E2in the fitting process and interpolation process is approximate±0.00365 V.After the calculation of θ ( E1,E2)and Ut(E1,E2),it is found that the absolute error of θ is ±0.3°and the relative error of Utis±2%.In a word,the maximum relative error of Utis±4% ± 0.02 m/s,whereas,the maximum relative error of θ is±0.3°.In the present work,the range of Utis 1.3-3.6 m/s,whereas,the range of θ is-7.5°to 2.3°.Considering θ may equal 0,which may lead the relative error of normal velocity approaches to in finite,the authors only calculate the absolute error of the parameters related to V.Through calculation,the relative error of U measured by X-type probe is 5.8%,whereas,the absolute value of V is 0.006-0.019 m/s.

    Table 1 Experimental conditions.

    The third important parameter is the friction velocity Uτ,the uncertainty of which is also discussed by Wei et al.22

    The fourth parameter is U+(U+=U/Uτ).According to the error propagation formula,

    According to the Wei et al.22

    Thus,the relative uncertainty of the friction velocity is obtained:

    2.5.Verification of techniques

    It is necessary to conduct the verification of the techniques before the measurement in the present work.

    For the X-type probe,the measurements are conducted at Station 3,with Re=19000 and Ro=0.The mean velocity profile measured by X-type probe and the standard log law

    are presented in Fig.5.

    According to the uncertainty analysis of U+,the minimum relative uncertainty is 11%.It can be seen in Fig.5 that the velocity profile measured in the present work almost overlaps with the standard logarithmic law,except for some slight difference.And the small difference can be explained as follows:first,there is some inevitable errors in the measurement of the wall shear stress;second,the standard log law is obtained in the duct with zero pressure gradient,whereas,the present work is performed in the duct which has a slightly positive pressure gradient.

    3.Results and Discussion

    3.1.Inlet conditions

    The turbulent intensity is discussed in the entrance section of a rotating channel flow;thus,the inlet conditions is very crucial,as presented in our previous paper.24

    3.2.Distribution of turbulent intensities

    Fig.5 Distribution of U+at Station 3.

    In our previous paper we investigated the distribution of themeasured by the hot wire anemometer,through the contrast of those between trailing side and leading side.It was found that the Coriolis force can suppress the turbulent fluctuation in the leading side,which is named the relaminarization process by Koyama et al.19When relaminarization process occurs,the value of the peak of the turbulent intensity decreases and the location of that moves away from the wall,which is in good agreement with the experimental result of Kim et al.25In addition,with the help of X-type probe,more details of turbulent intensities are listed in the present work.Although it has been found by Koyama19and Nakabayashi10et al.,the local rotation parameter Rc equals-0.007(i.e.the maximum absolute value of the Rc number in previous work),which is smaller than that in the present work,the maximum value of which is 0.0136.Thus the process of relaminarization in the leading side is captured,as is conveyed in Section 3.3.

    Because the X-type probe can measure U and V,〈uu〉/measured in the present work is the real 〈uu〉/Whereas,themeasured by Wei et al.22refers toIn addition,some differences can be seen in contrast with 〈 uu〉/in the leading side of present work and that of Wei et al.When y+equals about 100,the 〈 uu〉/in the leading side measured by Wei et al.is larger than that measured by the present work.This is becausemeasured by Wei et al. refers towhereasin Fig.6 refers to the true 〈 uu〉/Thus,it can be concluded that〈vv〉/is responsible for the difference of 〈 uu〉/measured by Wei et al.and the present work.The variation of 〈 vv〉/with streamwise location and rotation number is presented in Fig.7.It can be seen in Fig.7 that 〈 vv〉/also has the attenuation trend in the leading side,like〈 uu 〉/Because the relative error of 〈vv〉/is larger than that of 〈uu〉/the distribution of 〈 vv〉/is more scattered than that of 〈 uu〉/

    Fig.6 Variation of 〈uu〉/with streamwise location and rotation number.

    Fig.7 Variation of 〈vv〉/ with streamwise location and rotation number.

    In Eq.(5),Ujmeans one of the three components of the velocity whereas ujor uimeans one of the three components of the fluctuating velocity,means the production term of turbulent kinetic energy,Sijmeans the strain rate tensor and p means pressure.When 〈 uv〉/=0,the production of turbulent kinetic energy stops,which is also the essential characteristic of relaminarization process.The attenuation of 〈 uu〉/and〈vv 〉/is the result of stop of production of turbulent kinetic energy,but the turbulent kinetic energy equation is not only produced by the production term,but also transported by convection termand transportation termThus 〈uu〉/and〈vv〉/do not equal 0.

    Fig.8 Variation of 〈uv〉/ with streamwise location and rotation number.

    3.3.Quadrant analysis of velocity fluctuation

    In the investigation of Wei et al.,22the relaminarization phenomenon is found through the contrast ofIn order to have a better understanding of the characteristics of the turbulent fluctuation of the relaminarization and the effects of the Coriolis instability on the turbulent fluctuation,the authors analyzed the quadrant distribution of fluctuating velocity u and v.The coordinate system and quadrant used here are shown in Fig.9.

    Fig.9 Distribution of fluctuation in four quadrants under different rotation conditions at Station 3(y+=20).

    As we all know,there are a series of streaks in the inner layer of the turbulent boundary layer,and burst is one of its characteristics.Burst,the streak flows downstream with the fluid,will also move away from the wall at the same time.When y+is larger than 10,the movement of streaks away from the wall becomes faster;subsequently,the streaks vibrate violently,break down,with the vortex shedding from the wall.Burst is considered to be the primary form of the production of the turbulent kinetic energy.Besides,the movement of streaks away from the wall is named ejection.Ejection causes the relatively high speed fluid to fill the vicinity of the wall,and the process in which fluid flow to the wall is named sweep.In a word,ejection and sweep are considered to be the characteristics of the production of the turbulent kinetic energy.In the description above,we can see that ejection is the process in which the relatively low speed fluid moves away from the wall;therefore,in the four quadrant analysis,ejection lies in the second quadrant(i.e.u<0,v>0).Besides,sweep is the process in which the relatively high speed fluid comes to the wall;therefore,in the four quadrant analysis,sweep lies in the fourth quadrant(i.e.u>0,v<0).

    The distribution of the fluctuating velocity in the quadrant analysis measured at Station 3,y+=20,in a sampling period(i.e.10 s)under three rotating conditions is presented in Fig.9.It can be seen that there is a greater possibility in Fig.9(a)and Fig.9(b)for the fluctuating velocity lying in the second quadrant,and this is because the ejection happens at the measuring location.Comparing the distribution of the fluctuating velocity near the trailing side and that near the stationary wall,we can see the possibility that fluctuating velocity lies in the second quadrant near the trailing side is much higher than that near the stationary wall.This is because the ejection happens more frequently near the trailing side,resulting in a frequent production of turbulent kinetic energy,which is in good agreement with the conclusion that the Coriolis instability can promote the turbulent intensity.

    Comparing the distribution of the fluctuating velocity in the quadrant analysis near the leading side and that near the stationary wall,we can see that the distribution of fluctuation takes on the following characteristics: first,the fluctuation decreases sharply,which means that the kinetic energy of turbulent fluctuations dissipates very quickly after relaminarization.Second,the shape of the distribution of the fluctuation near the leading side is more flat than stationary cases;that is to say,the absolute value of v is much smaller than the absolute value of u,which means the momentum transfer is weakened in the normal direction.Third,the distribution of fluctuation is more uniform in four quadrants,which means ejection does not happen as frequently as it should do.

    3.4.Effect of Coriolis force on turbulent fluctuation

    Nakabayashi and Kitoh10simplified the Reynolds transport equation in the boundary layer and normalized the equation,and finally obtained the transport equation as follows:

    Through the analysis of Eqs.(6)–(8),the origin of the turbulent fluctuation near the leading side is the first term on the right-hand side of the Eq.(7),which is also the Coriolis term.Coriolis term transports the turbulent kinetic energy from the term 〈vv〉/to 〈uu〉/resulting in the decrease of the〈 vv 〉/and the increase of the term〈 uu 〉/Turbulent kinetic energy transferred from 〈 vv〉/to 〈 uu〉/results in the decrease of〈-uv〉/in two different ways.The first way is that the first term on the right-hand of the Eq.(8),the production term,is responsible for the decrease of-〈uv〉/;the second way is that the second term on the right-hand,the Coriolis term,leads to the decrease of-〈uv〉/By contrast,we can see from Eq.(6)that the first term on the right-hand of the equation(i.e.the production term)results in the decrease of〈uu 〉/while the second term on the right-hand of the equation(i.e.the Coriolis term)results in the increase of the〈uu〉/.In order to make clear the mechanism of the relaminarization,it is necessary to make it clear which term,Coriolis term or production term,has a predominant role.

    The distribution of the production term.and the

    3.5.Rotation correction of wall function

    It has been presented in the research of Wei et al.22that the turbulent boundary layer in the present case can be divided into several layers:viscous region,viscous-Coriolis region,Coriolis region and core region.Considering there are so many influencing factors in the core region,it is difficult to discuss the issues of rotation correction of the velocity formula at this region.Besides,the viscous region is free from the impact of the rotation,there is no need to discuss the rotation correction of the velocity formula at viscous region.Therefore,we mainly discuss the rotation correction of velocity formula at viscous-Coriolis region and Coriolis region.

    Fig.10 Coriolis term and production term.

    For the rotation correction of the velocity formula at the Coriolis region,considering the Coriolis region is derived from the logarithmic region,the authors try to correct the velocity formula from the very beginning of the deduction of the log law;whereas for the viscous-Coriolis region,considering this region is derived only from the viscous region,the authors also try to correct the velocity formula from the very beginning of the deduction of the wall function for the viscous region.By the way,the most common expression of the wall function in the viscous region is Van Driest formula,thus the correction of the velocity formula at the viscous-Coriolis region starts at the beginning of the deduction of the Van Driest formula.What’s more,log law and Van Driest formula has the same parameter,the von Karman coefficient κ.Therefore,the authors first correct the log law and get the von Karman coefficient κ with rotation correction;based on this κ,the Van Driest formula with rotation correction is finally obtained.

    The mean primary velocity profiles under different rotation numbers at Stations 1-4 are presented by Wei et al.,22and the difference of mean primary velocity profiles between the rotating cases and stationary cases,namely the δU+,at Station 2 and Station 3 is presented in Fig.11.As is shown in Fig.11,the difference near the trailing side has a logarithmic distribution in the log region,and the difference increases monotonously with the rise of the rotation number.The difference neartheleading side,apart from the casein which Ro=0.072,does not have the logarithmic distribution,and the difference near the leading side does not increase monotonously with the rise of rotation number.According to the analysis above,the authors speculate that relaminarization in the leading side under rotation conditions is responsible for the contrast of the differences between the trailing side and leading side.Therefore,the authors confirm the speculation by demonstrating that relaminarization is responsible for the difference of the velocity profiles near the leading side not having the logarithmic distribution.

    Fig.11 Difference of mean primary velocity profile measured at Station 2 and Station 3 between stationary cases and rotating cases.

    In the deduction of log law,there is an important assumption:local equilibrium hypothesis is valid in the log law region,which means that production term of turbulent fluctuation equals dissipation term and the equationis valid.can be neglected,compared to,and then there is,and then there is-〈uv〉=in log law region.

    The distribution of the dimensionless shear stress at different locations under different rotating conditions when 30< y+< 100 is presented in Fig.12.The formula of the dimensionless shear stress is

    Fig.12 Distribution of dimensionless shear stress under different rotating cases when 30< y+< 100.

    Fig.12 shows that the equationis not valid for the inner layer of the boundary layer in stationary cases.This is mainly because the small Reτ(Reτ=Uτδ/ν)in the present case leads to the attenuation of the-〈uv〉when y+<100.Therefore,the mean primary velocity profile also has the log law whensuggesting thatis not the prerequisite for the existence of log law.The authors consider that the log law of the mean primary velocity profile can be obtained when the distribution ofat the range of 30<y+<100 under rotating conditions is consistent with that in stationary cases(i.e.decreases monotonously with the rise of y+).Apart from the cases of Ro=0.290,Ro=0.362 at the leading side of Station 2 and the cases of Ro=0.145,0.217,0.290,0.362 at the leading side of Station 3,the distribution ofat the range of 30<y+<100 under other rotating conditions is consistent with that under stationary conditions.increases with y+in the rotating cases in which the distribution ofis not consistent with that stationary cases.Through the contrast of Fig.8 and Fig.12,it is not difficult to see thatincreasing with y+is the result of relaminarization.Meanwhile,the differences of the velocity profiles having the logarithmic distribution in Fig.11 are the cases in Fig.12,in which the distribution ofat the range of 30<y+<100 is consistent with that under stationary cases.And the differences not having the logarithmic distribution are the cases in Fig.12,in which the distribution ofat the range of 30<y+<100 is not consistent with that under stationary cases.Thus,it is demonstrated the relaminarization results in the difference of the primary velocity profile not having the logarithmic distribution.

    The primary velocity does not have a logarithmic distribution when relaminarization influences the velocity profile in the log region,not to mention rotation correction on the wall function of the Coriolis region which has a logarithmic distribution.Therefore,it is important to make sure whether the relaminarization influences the wall function within the log law region in the inner layer of the boundary layer.And thus,the critical conditions for the relaminarization must be obtained before the rotation correction for the wall function of the turbulent boundary layer.It is not difficult to find that the critical rotation number,at which the logarithmic distribution of the velocity profile within the log region is affected by relaminarization,decreases with the rise of the radius of rotation.This can be explained as follows:relaminarization is the process in which the inner layer of the boundary layer deviates from the equilibrium turbulent flow;in other words,the dissipation of the turbulent kinetic energy is larger than the generation.Therefore,the kinetic energy of turbulent flow decreases along the stream-wise location and the relaminarization becomes severer and severer.

    Fig.13 gives the critical condition for the relaminarization when Re=19000.In Fig.13,X/D equals 4.0625,5.3125,6.5625 respectively for the four streamwise locations shown in Fig.1.The point in Fig.13 means the maximum rotation number at which the velocity differences of the primary flow between the rotating conditions and the stationary conditions are not influenced by the relaminarization,which can be named the limit of the relaminarization.Fig.13 shows that the critical rotation number for the limit of relaminarization has a linear relationship with X/D.

    Fig. 13 Critical condition for relaminarization when Re=19000.

    With the limit of the relaminarization,we can perform the rotation correction of the wall function within the Coriolis region in the cases without the influence of relaminarization.It is known that all the forces have little influence on the log law region.Therefore,there should be no direct relationship between ?U/?y and the viscous force;in other words,?U/?y should not include viscous coefficient ν.However,considering the revised von Karman κ is the function of Rc,which is one of the governing parameters in the turbulent boundary layer under rotating conditions and includes the viscous coefficient,it seems that revising von Karman coefficient will lead?U/?y to includeν.But the von Karman coefficient is the constant related to the turbulent mixing length scale lm,considering there is lm=κy,which is valid in the region not dominated by viscous force.Therefore,in any given case,von Karman coefficient is still the constant in the physical essence,and revising von Karman coefficient can guarantee that ?U/?y does not directly include viscous coefficient ν.

    The variation of 1/κ with Rc at different streamwise locations is presented in Fig.14,and only the cases in which the flow near the leading side is not affected by relaminarization are presented.1/κ has the linear relationship with Rc on the whole.Therefore,the present paper takes the following fitting formula to describe the von Karman coefficient in the Coriolis region.

    Here,M and N are all fitting parameters of the curve in Fig.14.Through fitting the experimental data into Eq.(10),then it becomes

    Fig.14 Variation of 1/κ with Rc at different streamwise locations.

    It should be pointed out that the Eq.(11)is not valid for the flows in all kinds of channel geometries.First,we have studied the entry section of the turbulent flow.As we all know,the entry section of the turbulent flow is greatly affected by the boundary conditions.For example,if the crossing-section changes from square into trapezoid,the secondary flow in the cross-section may change,and Macfarlane et al.8have discussed the interaction between the secondary flow and the primary flow.

    Frankly speaking,there is some limitation in the generality of the conclusion;besides,the formula present in our paper is applied in the square cross-section.However,this paper is in the preliminary stage of our research.In the subsequent research,we will discuss the rotation correction of the boundary layer in different channel structures.However,the basic formula form=M+NRc may not change,only the value of the undetermined coefficient M and N may vary with the channel geometry.

    After revising von Karman coefficient κ,we can revise the intercept of the standard log law.However,the distribution of intercept is so sophisticated,as shown in Fig.15.Watmuff et al.3and Nakabayashi and Kitoh4do not give the revising formula of intercept C either.The reason why the influence of rotation on κ and C is so different is that κ is related to dU+/dy+,whereas,C is related to the integral of dU+/dy+.Therefore,the viscous and viscous-Coriolis region can also influence the intercept C.So,the present paper does not intend to revise the intercept C with rotation.

    Failing to find the rotation correction of the intercept C,the authors cannot revise the formula of the log law as the formula of Coriolis region.However,the authors find the revised formula of von Karman coefficient κ,and we can use the Van Driest formula to revise the wall function of the inner layer of the boundary layer based on the revised von Karman coefficient.The Van Driest formula is presented as follows:

    Fig.15 Variation of intercept C with Rc at different streamwise locations.

    where A+=UτA/ν and the A means the damping length scale in the stationary cases,which is an empirical parameter.Through fitting the mean velocity profiles,which is presented by Wei et al.,26into Eq.(12),the revised A+can be obtained as shown in Fig.16.It is not difficult to see that the variation of A+with Rc at different streamwise locations is similar to the intercept C,which means that there may be some implicit relationship between A+and intercept C.Therefore,it is also difficult to get a universal revised formula with rotation.In addition,seen from Fig.16,A+can be very large when Rc>0,which is physically contradictory.The physical meaning of A+is the dimensionless length scale dominated by viscousforce,and the larger A+means the larger are a dominated by viscous force.However,when Rc>0,the turbulent fluctuation near the trailing side increases,the area dominated by viscous force should shrink.

    But why is the A+obtained near the trailing side physically contradictory?This can be explained as follows:

    The deduction of Van Driest formula is based on an assumption:

    where κy[ 1-exp(-y/A)]means the turbulent mixing length scale.When the flow is far away from the wall,the scale equalsκy.Considering the damping effect near the wall,a damping function[1-exp(-y/A)]is added.However,for the boundary layer in the rotating channel,it is inappropriate to use κy[ 1-exp(-y/A)]to represent the turbulent mixing length scale.Instead,we should use the κ of the rotating cases in the Coriolis region,use the κ of the stationary cases in the viscous region.Only in this way the velocity profile of rotating cases and stationary cases can overlap near the wall without A+contradictory to its physical essence.Otherwise,the κ near the trailing side of rotating cases is larger than that in stationary cases,and the term A should be enlarged so that the velocity profile near the wall of rotating cases overlaps with that in stationary cases,resulting in the A+contradictory to its physicalessence near the trailing side.In a word,only κy[ 1-exp(-y/A)]cannot describe the distribution of the turbulent mixing length of the boundary layer of rotating cases,which also leads to the result why it is difficult to obtain a universal A+through Eq.(12).

    Fig.16 Variation of A+with Rc at different streamwise locations.

    Fig.17 Distribution of dimensionless mixing length for rotating cases.

    Based on the analysis above,the authors consider it is necessary to define the turbulent mixing length scale again to obtain the universal wall function in the rotating channels.The newly-defined turbulent mixing length scale must have the following characteristics.First,the new mixing length is κ1y[ 1-exp(-y/A)]in viscous region,here,κ1means the von Karman coefficient in the boundary layer in stationary cases.Second,from the inferior limit of the region influenced by the Coriolis force,namely theis a coefficient that determines the inferior limit of the influencing region of Coriolis force),the turbulent mixing length should vary from κ1y[ 1-exp(-y/A)]to κ2y smoothly,ending up with κ2y(κ2means the von Karman coefficient in the rotating channel)in the Coriolis region,as shown in Fig.17.To realize the smooth variation,it is necessary to add an additional damping length scale B.In a word,the newly defined turbulent mixing length scale is sophisticated:First of all,it is a segment function;second,an additional damping length scale B that should include penetration depth of Coriolis force is needed;third,the smooth variation from κ1y[ 1-exp(-y/A)]to κ2y is pretty sophisticated.Considering the complexity of the newly defined turbulent mixing length scale,the authors discuss some important characteristics it should have,rather than give the specific form of the newly defined turbulent length scale.

    4.Conclusions

    From the distribution of the component of Reynolds stress,it can be seen that relaminarization occurred near the leading side.Because the absolute value of local rotation parameter Rc is almost twice as large as that in previous researches,the whole process of relaminarization is presented.

    Through the analysis of the transport equation of Reynolds stress,it can be seen that though the Coriolis term indeed contributes to the attenuation of the fluctuation in the leading side,the influence of Coriolis term on attenuation is slight,compared to the production term.However,the Coriolis term is very important for the attenuation of turbulent fluctuation,considering the Coriolis term is the origin of the attenuation of turbulent fluctuation.The influence of the Coriolis term is amplified by the production term.

    Last,the present paper revises the wall function of the turbulent boundary layer not affected by the relaminarizaiton.The revised formula of von Karman coefficient is also presented.In the rotation correction of the wall function of the Coriolis region and the Van Driest formula in the viscous region,the authors fail to obtain the intercept of the log law or the universal A+in the Van Driest formula,thus a new revising method is proposed.

    Acknowledgments

    The present work is supported by the National Natural Science Foundation of China(No.51541605).

    看片在线看免费视频| 91麻豆av在线| 少妇丰满av| 国产午夜福利久久久久久| 免费av观看视频| 日本在线视频免费播放| 欧美xxxx黑人xx丫x性爽| 我的老师免费观看完整版| 亚洲一区二区三区不卡视频| xxx96com| 久久午夜亚洲精品久久| 有码 亚洲区| 亚洲一区二区三区色噜噜| 免费观看精品视频网站| 听说在线观看完整版免费高清| 男女午夜视频在线观看| 欧美丝袜亚洲另类 | 亚洲一区二区三区色噜噜| 人人妻人人看人人澡| 51国产日韩欧美| 51国产日韩欧美| 日韩av在线大香蕉| 久久精品国产清高在天天线| 国内精品美女久久久久久| 国产探花在线观看一区二区| 最近在线观看免费完整版| 九九热线精品视视频播放| 国产私拍福利视频在线观看| 全区人妻精品视频| 国产探花极品一区二区| 成人一区二区视频在线观看| 欧美大码av| 亚洲男人的天堂狠狠| 亚洲av日韩精品久久久久久密| 1000部很黄的大片| 欧美+日韩+精品| 亚洲中文日韩欧美视频| 女人被狂操c到高潮| 免费人成视频x8x8入口观看| 国产欧美日韩一区二区三| 亚洲一区二区三区色噜噜| 日韩国内少妇激情av| 亚洲av不卡在线观看| 久久九九热精品免费| 日日干狠狠操夜夜爽| 久久久国产精品麻豆| 天堂动漫精品| 国产高潮美女av| 国产伦在线观看视频一区| 色尼玛亚洲综合影院| 久99久视频精品免费| 亚洲五月天丁香| 午夜亚洲福利在线播放| 国产精品女同一区二区软件 | 日本黄色片子视频| 欧美日韩乱码在线| 国产不卡一卡二| 国产精品乱码一区二三区的特点| 精品人妻一区二区三区麻豆 | 免费看光身美女| 蜜桃亚洲精品一区二区三区| 操出白浆在线播放| 中亚洲国语对白在线视频| 精品乱码久久久久久99久播| 日本免费一区二区三区高清不卡| 国产精品1区2区在线观看.| 亚洲av中文字字幕乱码综合| 99久久成人亚洲精品观看| 桃色一区二区三区在线观看| 十八禁人妻一区二区| 我要搜黄色片| 99国产综合亚洲精品| 国产三级黄色录像| 亚洲片人在线观看| 神马国产精品三级电影在线观看| 日韩精品中文字幕看吧| 免费av观看视频| 看免费av毛片| 国产精品久久电影中文字幕| 99热这里只有精品一区| 中文在线观看免费www的网站| 99久久精品一区二区三区| 国产色爽女视频免费观看| 精品99又大又爽又粗少妇毛片 | 十八禁人妻一区二区| 综合色av麻豆| 亚洲国产欧美人成| АⅤ资源中文在线天堂| 他把我摸到了高潮在线观看| 3wmmmm亚洲av在线观看| 亚洲成人精品中文字幕电影| 午夜福利18| 日本与韩国留学比较| 日本成人三级电影网站| 欧美不卡视频在线免费观看| 国产伦在线观看视频一区| 免费观看精品视频网站| 亚洲在线观看片| 午夜精品久久久久久毛片777| 免费av毛片视频| 欧美绝顶高潮抽搐喷水| 99热6这里只有精品| 91麻豆精品激情在线观看国产| 国产亚洲精品久久久久久毛片| 一本精品99久久精品77| 日韩高清综合在线| 亚洲精品成人久久久久久| 婷婷亚洲欧美| 国产色爽女视频免费观看| 亚洲国产精品成人综合色| 国产不卡一卡二| 欧美日韩精品网址| 亚洲avbb在线观看| 免费看美女性在线毛片视频| 国产色婷婷99| 淫秽高清视频在线观看| www.www免费av| 成人亚洲精品av一区二区| 免费人成在线观看视频色| 激情在线观看视频在线高清| 欧美三级亚洲精品| 尤物成人国产欧美一区二区三区| 可以在线观看的亚洲视频| 少妇的丰满在线观看| 久久久久精品国产欧美久久久| 真人做人爱边吃奶动态| 久久精品亚洲精品国产色婷小说| 中文亚洲av片在线观看爽| 久久精品夜夜夜夜夜久久蜜豆| 欧美一区二区亚洲| 午夜精品在线福利| 一本精品99久久精品77| 日本熟妇午夜| 小说图片视频综合网站| 亚洲第一电影网av| 欧美激情久久久久久爽电影| 国产不卡一卡二| 99久久无色码亚洲精品果冻| 亚洲国产精品久久男人天堂| 老熟妇乱子伦视频在线观看| 久久亚洲精品不卡| 国产黄色小视频在线观看| 欧美中文日本在线观看视频| 午夜免费男女啪啪视频观看 | 国产成人啪精品午夜网站| 欧美+日韩+精品| 中文字幕av成人在线电影| 美女黄网站色视频| 欧美成人性av电影在线观看| 无人区码免费观看不卡| 天堂影院成人在线观看| 成年女人看的毛片在线观看| 丰满人妻一区二区三区视频av | 制服丝袜大香蕉在线| 床上黄色一级片| 舔av片在线| 国产黄色小视频在线观看| 国产成人影院久久av| 久久性视频一级片| 国产中年淑女户外野战色| 中出人妻视频一区二区| 女警被强在线播放| 美女 人体艺术 gogo| 久久草成人影院| 日本在线视频免费播放| 啦啦啦观看免费观看视频高清| 国产激情偷乱视频一区二区| 国产精品99久久久久久久久| 成人高潮视频无遮挡免费网站| 欧美乱色亚洲激情| 久久久久久久亚洲中文字幕 | 狂野欧美白嫩少妇大欣赏| 国产欧美日韩一区二区精品| 欧美黄色淫秽网站| 观看免费一级毛片| 99久久精品一区二区三区| 757午夜福利合集在线观看| 91av网一区二区| 免费搜索国产男女视频| 国产成人啪精品午夜网站| 在线a可以看的网站| 精品国产三级普通话版| 观看美女的网站| 女人被狂操c到高潮| 99久久精品一区二区三区| 日韩欧美国产一区二区入口| 嫩草影院精品99| 一级作爱视频免费观看| 麻豆一二三区av精品| 国产爱豆传媒在线观看| 欧美性猛交黑人性爽| 欧美色欧美亚洲另类二区| 熟妇人妻久久中文字幕3abv| 大型黄色视频在线免费观看| 色噜噜av男人的天堂激情| 在线观看午夜福利视频| 高清在线国产一区| 亚洲电影在线观看av| 免费高清视频大片| 国产中年淑女户外野战色| 丰满人妻熟妇乱又伦精品不卡| 国产精品嫩草影院av在线观看 | 亚洲av中文字字幕乱码综合| 在线观看舔阴道视频| 成人av一区二区三区在线看| 久久久成人免费电影| 90打野战视频偷拍视频| 国产中年淑女户外野战色| 成人精品一区二区免费| 亚洲精品亚洲一区二区| 国产精品日韩av在线免费观看| 亚洲av一区综合| 俄罗斯特黄特色一大片| tocl精华| 国产精品永久免费网站| 老司机深夜福利视频在线观看| 亚洲精品在线观看二区| 十八禁网站免费在线| 久久久久久久午夜电影| 欧美在线一区亚洲| 级片在线观看| 老司机在亚洲福利影院| 欧美又色又爽又黄视频| www日本黄色视频网| 色综合站精品国产| 日本a在线网址| 美女高潮的动态| 一区二区三区国产精品乱码| 亚洲人成网站在线播放欧美日韩| 国产精品影院久久| 亚洲,欧美精品.| 岛国在线免费视频观看| 一本一本综合久久| 色av中文字幕| 亚洲国产日韩欧美精品在线观看 | 色综合站精品国产| 观看美女的网站| 成年免费大片在线观看| 免费av观看视频| 好男人电影高清在线观看| 亚洲国产欧洲综合997久久,| 麻豆国产av国片精品| 九色国产91popny在线| 亚洲精品一卡2卡三卡4卡5卡| 99久久久亚洲精品蜜臀av| 国产精品av视频在线免费观看| 黄色视频,在线免费观看| 男插女下体视频免费在线播放| 亚洲av不卡在线观看| 男女那种视频在线观看| 免费在线观看成人毛片| 夜夜夜夜夜久久久久| 久久精品91蜜桃| 内射极品少妇av片p| 久久国产精品影院| 国产高潮美女av| 搡女人真爽免费视频火全软件 | 最好的美女福利视频网| 精品乱码久久久久久99久播| www.色视频.com| 搡老妇女老女人老熟妇| 啪啪无遮挡十八禁网站| 免费观看的影片在线观看| 三级毛片av免费| 免费一级毛片在线播放高清视频| 女同久久另类99精品国产91| 久久精品国产亚洲av香蕉五月| 淫秽高清视频在线观看| 午夜免费激情av| 国产精品 欧美亚洲| 91九色精品人成在线观看| 国产三级中文精品| 老汉色av国产亚洲站长工具| 在线观看免费午夜福利视频| 99久久精品一区二区三区| 成年女人看的毛片在线观看| 免费高清视频大片| 久久午夜亚洲精品久久| 国产久久久一区二区三区| 麻豆一二三区av精品| 最好的美女福利视频网| 天美传媒精品一区二区| 成年免费大片在线观看| 免费人成在线观看视频色| 亚洲中文字幕一区二区三区有码在线看| 一级毛片女人18水好多| 国产在视频线在精品| 色综合婷婷激情| 女同久久另类99精品国产91| 国产精品av视频在线免费观看| 深爱激情五月婷婷| 女警被强在线播放| 国产成人av激情在线播放| 免费无遮挡裸体视频| 日本与韩国留学比较| 国产欧美日韩精品亚洲av| 欧美色视频一区免费| 午夜久久久久精精品| 波野结衣二区三区在线 | 欧美+亚洲+日韩+国产| 黄色视频,在线免费观看| 美女高潮的动态| 18+在线观看网站| av天堂中文字幕网| 男女下面进入的视频免费午夜| 偷拍熟女少妇极品色| e午夜精品久久久久久久| 中文字幕高清在线视频| 午夜福利高清视频| 国产男靠女视频免费网站| 免费看a级黄色片| 神马国产精品三级电影在线观看| 1000部很黄的大片| 亚洲电影在线观看av| 全区人妻精品视频| 成年版毛片免费区| 国产探花极品一区二区| 真人做人爱边吃奶动态| 中文字幕久久专区| 国产 一区 欧美 日韩| 真实男女啪啪啪动态图| 国产精品一区二区免费欧美| 亚洲av电影不卡..在线观看| 在线免费观看不下载黄p国产 | 午夜福利18| 国产av不卡久久| 亚洲欧美激情综合另类| 在线天堂最新版资源| 国产一区二区在线av高清观看| 成人亚洲精品av一区二区| 日韩欧美 国产精品| 国产三级黄色录像| 欧美丝袜亚洲另类 | 一个人看视频在线观看www免费 | 日本a在线网址| 在线观看日韩欧美| 一本精品99久久精品77| 国产免费男女视频| 久久久久九九精品影院| 一边摸一边抽搐一进一小说| 观看免费一级毛片| 国产精品亚洲一级av第二区| 久久久成人免费电影| 亚洲成人久久爱视频| 免费一级毛片在线播放高清视频| 99国产极品粉嫩在线观看| 久久久久久久精品吃奶| 国产久久久一区二区三区| 欧美乱码精品一区二区三区| 午夜福利18| 叶爱在线成人免费视频播放| 午夜福利成人在线免费观看| av福利片在线观看| 一a级毛片在线观看| 国产伦人伦偷精品视频| 亚洲成a人片在线一区二区| 精品久久久久久久久久免费视频| 少妇熟女aⅴ在线视频| 欧美日本亚洲视频在线播放| 亚洲av成人av| 亚洲最大成人手机在线| 在线播放国产精品三级| 在线十欧美十亚洲十日本专区| 神马国产精品三级电影在线观看| 国产欧美日韩一区二区三| 精品欧美国产一区二区三| 特级一级黄色大片| 国产高清激情床上av| 真人一进一出gif抽搐免费| 亚洲国产中文字幕在线视频| 国产精品国产高清国产av| 男女做爰动态图高潮gif福利片| 久久天躁狠狠躁夜夜2o2o| 日韩欧美在线二视频| 精品不卡国产一区二区三区| 88av欧美| 国产中年淑女户外野战色| 特级一级黄色大片| 免费av毛片视频| 国产成人影院久久av| 免费一级毛片在线播放高清视频| 国产男靠女视频免费网站| 香蕉丝袜av| 美女 人体艺术 gogo| 内射极品少妇av片p| 他把我摸到了高潮在线观看| 99精品久久久久人妻精品| 精品一区二区三区av网在线观看| 淫妇啪啪啪对白视频| 一级a爱片免费观看的视频| 无人区码免费观看不卡| 变态另类成人亚洲欧美熟女| 天天添夜夜摸| 国产成人啪精品午夜网站| 久久人人精品亚洲av| 精品熟女少妇八av免费久了| 变态另类成人亚洲欧美熟女| 中文字幕人妻丝袜一区二区| 亚洲精品456在线播放app | 少妇裸体淫交视频免费看高清| 精品人妻一区二区三区麻豆 | 久久久久九九精品影院| 午夜免费男女啪啪视频观看 | 国产中年淑女户外野战色| 精品国产超薄肉色丝袜足j| 69av精品久久久久久| 国产成+人综合+亚洲专区| 亚洲精品色激情综合| 亚洲第一电影网av| 97人妻精品一区二区三区麻豆| 国产日本99.免费观看| 久久久久久人人人人人| 男人和女人高潮做爰伦理| 女同久久另类99精品国产91| 夜夜爽天天搞| 久久精品国产亚洲av涩爱 | 夜夜看夜夜爽夜夜摸| 国产91精品成人一区二区三区| 午夜福利在线观看吧| bbb黄色大片| 一区二区三区激情视频| 国产精品久久久久久人妻精品电影| 91九色精品人成在线观看| 中文字幕人妻丝袜一区二区| 在线播放国产精品三级| 美女大奶头视频| 桃色一区二区三区在线观看| 中文亚洲av片在线观看爽| 69人妻影院| 亚洲中文日韩欧美视频| 国产精品,欧美在线| 少妇的丰满在线观看| 一区二区三区高清视频在线| 一个人观看的视频www高清免费观看| 国产伦一二天堂av在线观看| 国产一区二区三区在线臀色熟女| 中文字幕久久专区| 91麻豆精品激情在线观看国产| 偷拍熟女少妇极品色| 琪琪午夜伦伦电影理论片6080| 男女那种视频在线观看| 在线免费观看的www视频| 国产高清videossex| 国产高清激情床上av| 99热这里只有精品一区| 搡老岳熟女国产| 中文字幕精品亚洲无线码一区| 身体一侧抽搐| 亚洲成av人片免费观看| 亚洲成人免费电影在线观看| 午夜福利在线观看吧| 老熟妇乱子伦视频在线观看| 久久6这里有精品| 欧美日韩瑟瑟在线播放| 制服丝袜大香蕉在线| 亚洲国产日韩欧美精品在线观看 | 亚洲精品456在线播放app | av专区在线播放| 欧美午夜高清在线| 欧美另类亚洲清纯唯美| 亚洲真实伦在线观看| 久99久视频精品免费| 国产午夜福利久久久久久| 国产精品日韩av在线免费观看| 国产精品久久视频播放| 亚洲va日本ⅴa欧美va伊人久久| 亚洲精品国产精品久久久不卡| 丰满人妻熟妇乱又伦精品不卡| 午夜福利18| 欧美日韩一级在线毛片| 中文亚洲av片在线观看爽| 深爱激情五月婷婷| 国产在视频线在精品| 男人舔奶头视频| 欧美成人一区二区免费高清观看| 欧美乱码精品一区二区三区| 波多野结衣高清作品| www国产在线视频色| 我的老师免费观看完整版| 两个人视频免费观看高清| 夜夜躁狠狠躁天天躁| 久久精品国产自在天天线| 麻豆国产av国片精品| 日本熟妇午夜| 久久久久久久久久黄片| 一卡2卡三卡四卡精品乱码亚洲| 国产亚洲欧美98| 欧美区成人在线视频| 一夜夜www| 丁香六月欧美| 亚洲熟妇中文字幕五十中出| 天堂影院成人在线观看| 美女大奶头视频| 一本久久中文字幕| 少妇熟女aⅴ在线视频| 琪琪午夜伦伦电影理论片6080| 最近最新中文字幕大全电影3| 亚洲精品一区av在线观看| 国产淫片久久久久久久久 | 狠狠狠狠99中文字幕| 国产高清videossex| 一本综合久久免费| 亚洲国产欧洲综合997久久,| 蜜桃久久精品国产亚洲av| 国产一区二区三区视频了| 欧美成人性av电影在线观看| 亚洲成av人片在线播放无| 国产精品,欧美在线| 成年免费大片在线观看| 老熟妇仑乱视频hdxx| 成人无遮挡网站| 亚洲国产精品合色在线| 女生性感内裤真人,穿戴方法视频| 国产在视频线在精品| 男人舔女人下体高潮全视频| 无人区码免费观看不卡| 99久久久亚洲精品蜜臀av| 亚洲五月婷婷丁香| 欧美日韩综合久久久久久 | 亚洲成a人片在线一区二区| 欧美三级亚洲精品| 精品久久久久久久末码| 成人永久免费在线观看视频| 国产av不卡久久| 精品人妻一区二区三区麻豆 | 午夜两性在线视频| 免费人成视频x8x8入口观看| 国产精品99久久久久久久久| 国产一区二区三区在线臀色熟女| 中文字幕高清在线视频| АⅤ资源中文在线天堂| 国产高清三级在线| 美女黄网站色视频| 成人18禁在线播放| 最新中文字幕久久久久| 蜜桃亚洲精品一区二区三区| 免费看a级黄色片| 日韩有码中文字幕| 老汉色∧v一级毛片| 午夜精品在线福利| 亚洲国产高清在线一区二区三| 免费看十八禁软件| 少妇高潮的动态图| 69人妻影院| 又黄又爽又免费观看的视频| 18禁美女被吸乳视频| 99国产精品一区二区三区| 我要搜黄色片| 欧美+日韩+精品| 天美传媒精品一区二区| 特级一级黄色大片| 亚洲av一区综合| 99国产精品一区二区三区| 亚洲七黄色美女视频| 特级一级黄色大片| 搡老妇女老女人老熟妇| 老司机深夜福利视频在线观看| 叶爱在线成人免费视频播放| 91在线观看av| 淫秽高清视频在线观看| 91在线观看av| 俄罗斯特黄特色一大片| 香蕉久久夜色| 亚洲av一区综合| 亚洲无线观看免费| 欧美性猛交╳xxx乱大交人| 在线观看舔阴道视频| 国内毛片毛片毛片毛片毛片| 久久精品国产亚洲av涩爱 | 波野结衣二区三区在线 | 成熟少妇高潮喷水视频| 五月伊人婷婷丁香| 成人亚洲精品av一区二区| 国产一区二区在线观看日韩 | 人人妻,人人澡人人爽秒播| 热99re8久久精品国产| 欧美+亚洲+日韩+国产| 无人区码免费观看不卡| 伊人久久大香线蕉亚洲五| 亚洲人与动物交配视频| 欧美一区二区精品小视频在线| 午夜福利视频1000在线观看| 亚洲电影在线观看av| 亚洲五月婷婷丁香| 精品国产超薄肉色丝袜足j| 免费观看的影片在线观看| 欧美三级亚洲精品| 给我免费播放毛片高清在线观看| 99视频精品全部免费 在线| 国产精品久久视频播放| 日本成人三级电影网站| 免费看十八禁软件| 在线观看免费视频日本深夜| 一区二区三区免费毛片| 中文亚洲av片在线观看爽| 99精品久久久久人妻精品| 亚洲精品一区av在线观看| 成人鲁丝片一二三区免费| 少妇人妻一区二区三区视频| 国产97色在线日韩免费| 日本五十路高清| 欧美3d第一页| 国产真实伦视频高清在线观看 | 日韩精品青青久久久久久| 中文字幕精品亚洲无线码一区| 又粗又爽又猛毛片免费看| x7x7x7水蜜桃| 99国产极品粉嫩在线观看| 国产亚洲欧美98| 一区福利在线观看| 99久久成人亚洲精品观看| 一进一出好大好爽视频| 国产伦精品一区二区三区四那| 欧美性感艳星| 白带黄色成豆腐渣|