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

    Experimental study on flow excursion instability of supercritical hydrocarbon fuel in scramjet regenerative cooling parallel channels

    2023-02-09 08:59:16XioyongLIUZhuqingYANGRuipuMIAOFengLIUShujunZHAOQinchengBI
    CHINESE JOURNAL OF AERONAUTICS 2023年1期

    Xioyong LIU, Zhuqing YANG, Ruipu MIAO, Feng LIU, Shujun ZHAO,Qincheng BI

    a Beijing Power Machinery Research Institute, Beijing 100074, China

    b Liaoning Key Laboratory of Complex Energy Conversion and Utilization, School of Energy and Power Engineering, Dalian University of Technology, Dalian 116024, China

    c State Key Laboratory of Multiphase Flow in Power Engineering, Xi’an Jiaotong University, Xi’an 710049, China

    KEYWORDS Flow distribution;Ledinegg instability;Onset of flow excursion;Parallel channels;Regenerative cooling;Scramjet;Supercritical flow

    Abstract Flow instability of supercritical hydrocarbon fuel is a crucial issue in scramjet regenerative cooling structure.In this study,flow excursion instability and flow distribution in parallel tubes were experimentally studied for supercritical fluids. Two types of flow excursion occur in a single tube.Type I and Type II excursions,and they are corresponding to decreasing and increasing flow rate respectively. They can trigger flow maldistribution between parallel tubes and the hysteresis phenomenon of flow distribution. The effects of system parameters, including inlet temperature,system pressure, and heat flux, on flow distribution were analyzed. In addition, the relationship between flow excursion and the pseudo-critical interval proposed in the literature was established according to the heated tube outlet temperature at the onset of flow instability. Finally, the flow excursion instability boundary was obtained using two dimensionless parameters.These experimental results can provide helpful insight on the mechanism of Scramjet regenerative cooling.

    1. Introduction

    Hypersonic aircrafts have attracted global attentions because they are developed by utilizing the advanced propulsion technology. The high speed of the aircraft subjects the engine surface to intense aerodynamic heating as well as combustion thermal load.1Researchers found that the scramjet combustion chamber temperature is up to 3000 K when the Mach number is 6.2While the latest thermal protection composite material can bear a maximum temperature of 2200 K, which means only using the heat resistant material is not feasible to ensure the operation safety of the scramjet.3A method of effectively cooling scramjet engines has become a key technology for flight safety.4-6Regenerative active cooling technology can not only alleviate the scramjet chamber thermal load, but also preheat the fuel and enhance the combustion efficiency.7,8Scholars regard it as one of the most promising thermal protection technologies because of its practicability and superiority.9-11Hydrocarbon fuel was used as coolant of scramjet regenerative cooling structure in the X-51A flight test.12The results indicated that the regenerative cooling technology is feasible to scramjet combustion chamber thermal protection.As a trans-critical process exist in regenerative cooling channels,the flow and heat transfer process of the coolant is extremely complicated.13,14Supercritical hydrocarbon fuel flow instability and heat deterioration have a profound influence on the scramjet regenerative cooling technology.A multichannel structure design has been used in the heated surface of regenerative active cooling systems. The size and design of the structure are closely related to heat distribution and the hydrocarbon fuel cooling process.15Hydrocarbon fuel flows in channels under adverse conditions,such as high temperature and pressure,high heat flux,and high superheat.Hydrocarbon fuel cracks and becomes small molecules also brings a large heat sink. The fluid gas-liquid interface disappeared at a supercritical state. Thus, the fluid thermophysical properties have a drastic variation near the pseudo-critical point. The scramjet cooling system’s supercritical pressure is less prone to uneven flow distribution and heat transfer deterioration near the pseudo-critical temperature.16,17These characteristics of hydrocarbon fuel induce complex hydrodynamic characteristics and flow instability in parallel tubes. The pipe structure,thermal load distribution, and fuel status parameters are crucial factors that influence the flow process.

    There are many serious safety problems caused by the flow instability and uniform flow distribution of hydrocarbon fuel,such as a negative influence on the thermal capacity of the heat exchange structure and thermal stress imbalance.In our previous research, the reliability of large plate’s active regenerative cooling structure at rather high heat flux condition was studied. Superheated ablation phenomenon was observed because mass flow rate deviation was found in parallel channels, as showed in Figs. 1(a) and 1(b), where m1and m2are mass flow rate in Tube 1 and Tube 2, respectively, and mallrespresents mass flow rate in summational tube. In addition, the overheated area moved from the end to the middle of the channels suggested a seriously uneven flow maldistribution in the test.Insufficient cooling in scramjet engine will pose a desperate threat to the safety. Fundamental and mechanism research on supercritical fluid flow instability characteristics in parallel tubes is significantly needed to solve this engineering problem.

    Flow instability can be classified into two representative types:static instability and dynamic instability;it occurs when the fluid encounters a small disturbance, and the fluid flow shows a state of flow excursion or flow oscillation.Static instability was first proposed by Ledinegg,18which is defined by the cubic curve relation between the pressure drop and mass flow rate relation of a two-phase flow. It is also called the ‘‘hydrodynamic multivalued characteristic curve.” For the multivalued characteristics, scholars had chosen the minimum point as the feature point to study.19-21Many researchers have proposed experimental fitting correlation formulas and theoretical prediction formulas.20,22-25They considered that the minimum point in the hydrodynamic multivalued characteristic curve was approximately equal to the Onset of Flow Instability(OFI). Swapnalee et al.26and Mahmoudi27believed that a drastic change in density is the main factor that causes hydrodynamic multivalued characteristics.According to the effect of system parameters on the hydrodynamic multivalued characteristic curve,various measures have been proposed to restrain multivalued characteristics.28-31Through these measures,increasing the system pressure, increasing the inlet temperature, and decreasing the heat flux were found to cause an increase in the fluid stability. Most researchers have studied supercritical fluid flow instability based on ultra-supercritical boiler and supercritical water-cooled reactor. Supercritical hydrocarbon flow in the scramjet regenerative cooling channels process is similar to the boiler and water-cooled reactor.Many scholars have investigated flow excursion,flow distribution, and flow oscillation characteristics to study the flow instability of supercritical hydrocarbon fuel in regenerative cooling parallel channels both experimentally and numerically.Liu et al.32has defined Static Flow Instability(SFI)to describe hydrodynamic characteristics for supercritical fluid.SFI is similar to the Ledinegg instability caused by some disturbance and trigger to flow excursion.Yang et al.25,33experimentally investigated the flow characteristics of supercritical cyclohexane and concluded that OFI is the minimum point of the hydrodynamic multivalued characteristic curve. They also found that the fluid outlet temperature reached its pseudo-critical temperature at this point. A one-dimensional transient model has been developed in Jin et al.’s investigation.34Their calculation results showed that the negative slope area of the multi-value curve is unstable.Many researchers believe that flow excursion occurs as the system enters the negative slope area of hydrodynamic multi-value characteristic curve29,30,35.

    Fig. 1 Abnormal phenomena found in authors’ former regenerative cooling experiments.

    Flow instability investigations have been carried out for parallel tubes based on the features of a single tube. In previous flow instability analyses,some scholars used the simulation software to study flow instability. They established a flow instability relationship between single and parallel tubes and obtained a criterion.36-38Wang et al.39conducted an experimental study on boiling fluids in parallel heated tubes. It was concluded that the hydrodynamic characteristic curve could be used to describe the Ledinegg instability in parallel tubes.Pressure drop oscillation occurs when the system is in the negative slope area of the curve, and density drop oscillation occurs in the other area. From their perspective, this phenomenon is attributed to the coupling of static and dynamic instabilities. Chatoorgoon40first used the state equation to predict the flow characteristics. Experimental and simulation studies were conducted by Xiong et al.41,42.Their results show that Ledinegg instability is difficult to observe in parallel tubes.Chen et al.43used a simulation method to analyze supercritical water flow in parallel tubes,and their results show that the factors that impact parallel tubes are similar to those of a single tube. As the flow direction decreases, flow excursion occurs.The results have not been verified experimentally. The threshold of heat flux has also attracted attention.44-47When the heat flux threshold is increased, the flow system becomes more stable.The instability boundary is obtained using the heat flux threshold. By controlling the system parameters, the relationship between the system parameters and flow instability can be established.48-50Approaches to verify if tubes are uniformly heated have also been widely investigated51-54.

    The static (Ledinegg) flow instability is described that the flow excursion occurrs as a result of some disturbance, and the flow condition will not go back to the original state. Flow excursion is commonly considered to exist in two-phase flow,and it can be related to the hydrodynamic multivalued characteristic curve. As summarized in Table 1,32,41,42,44-46,48-50,54-58recently study of supercritical flow in parallel channels mostly use SCWR as background. Supercritical water was usually used as the working fluid in parallel and vertical channels.Previous researchers focused on the flow distribution caused by pressure and density oscillations. For supercritical fluid flow in parallel tubes, seldom researchers combine the hydrodynamic multi-value curve with flow distribution to investigate Ledinegg instability. There is no definitive conclusion about flow in parallel tubes—specifically, whether Ledinegg instability exists under supercritical conditions.In addition,compared with the flow in a single tube,there are more factors that affect the flow characteristics in parallel tubes. Further research on the deterioration mechanisms of heat transfer and flow instability is required.6Thus,the exploration of flow characteristics is becoming one of the major directions toward the feasibility and safety design of supercritical hydrocarbon fuel systems.

    In this study, the flow excursion was considered to be the phenomenon that Ledinegg instability in parallel channels.Based on the classification of the flow excursion phenomenon by Liu et al.,59-61the hydrodynamic and heat transfer characteristics of supercritical cyclohexane flow in two heated parallel channels were investigated experimentally. Type I and Type II excursions correspond to a decrease and increase in the mass flow rate as showed in Fig. 259respectively, where m is mass flow rate and Δp represent the pressure drop through the pipe.Using the hydrodynamic multivalued characteristic curve, the flow process of two tubes under decreasing and increasing mass flow rates is described.Two types of instability boundary are proposed based on the experimental results.The predictive correlation of the unstable boundary is also presented. These conclusions can provide theoretical guidance for the regenerative active cooling technology of scramjet engines.

    2. Experimental system and data reduction

    2.1. Experimental system

    The experimental system is shown in Fig. 3 and illustrated in Fig. 4. The working fluid through the piston pump (Elite P500+, 0-500 mL/min) flows into the preheated section and is heated to the desired temperature. The flow is divided into two parallel heat tubes through a distribution header. After the test section, the fluid flows back to the main road through the collection header and through the water cooler.Finally,the cooled fluid flows back into the fuel tank for the next circulation.

    The temperatures of the inlet and outlet fluids were measured using an armored K-type thermocouple (Omega, 0-590 ℃).Volumetric flowmeters (CX-M5.1-SS, 0-800 mL/min) were used to measure the mass flow rate in parallel tubes. The test section was made of a horizontal stainless-steel tube with an inner diameter of 2.0 mm. Eight K-type thermocouples were equally spaced on the test section of each channel to measure the outside wall temperature.The pressure drop in the test section was measured using a differential pressure sensor (Rosemount 3051CD4, 0-1000 kPa). The data acquisition system consisted of NI 9209,NI 9205,NI 9212,and Keithley 2700.

    Cyclohexane with a purity of 99.5wt% was used for the working fluid in this study.The flow excursion and flow distribution phenomena observed under multiple conditions and the operating conditions for all experiments are listed in Table 2.The outlet temperature of the bulk fluid was controlled within 420 °C, and the wall temperature was lower than 540 °C to avoid thermal cracking of the fuel.

    2.2. Data reduction

    Taking the incompressible potential flow round a cylinder for example, the stream function is to obtain an accurate expression of the heat flux, a preliminary experiment was conducted to calculate the heat loss. The heat flux transformed through the outside surface of the heated test section is obtained using Eq. (1).

    where IU is the total input heat,I is the heating electricity current, U represent the heating voltage, q is heat flux of tube, D and L represent the outer diameter and heating length of the tube in the test section, respectively, and Qlossis the heat loss of the test section.

    Table 1 Reported supercritical fluid flow in parallel tubes.

    Fig.2 Type I and Type II flow excursion defined by Liu et al.59.

    The heat loss from the test section was less than 10%of the total power input and was calculated using Eq.(2),and the calculate step is as follows. Firstly, the test section was heated without fluid and the wall temperature and the power applied(the heat loss at this time) was measured and recorded. Secondly, the wall temperature was corresponded to its corresponding heat loss one by one and the Qloss-Twallcurve can be obtained. Finally, Eq. (2) could be obtained from the formula fitting of this heat loss curve:

    To verify the reliability of the experiment, a heat balance test was conducted. When the difference in inlet and outlet enthalpy with the heat flux was examined to calculate the heat efficiency, the results indicated that the efficiency of the test section was ~90 %. Then, the experimental data were used to obtain Nuexp(Nusselt number calculated from experimental data), and the experimental results were compared with the classical empirical heat transfer formula.

    Fig. 3 Photos of experimental system and parallel test sections.

    Fig. 4 Schematic diagram of experimental system.

    Table 2 Experimental conditions.

    The average temperature of the inner wall of the tube can be calculated using the measured heat,tube outer wall temperature,and physical properties of the fluid.The arithmetic average temperature can be considered as the temperature of the fluid,62the heat transfer coefficient of the fluid can be calculated, and then Nuexpcan be obtained.

    where Q(x)is internal heat source;d calculated the inner diameter of the tube in the test section;Twi,Tworepresent the inner and outer surface temperatures,respectively;and λ is heat conductivity coefficient.

    The heat transfer coefficient can be calculated using the diameter of the inner wall of the tube and the fluid temperature, and then Nuexpcan be calculated using Eq. (5).

    where Re is Reynolds number; Prbrepresent the fluid bulk of Prandtl number;and Prwicalculated the Prandtl number of the fluid near the inner wall.

    Five groups of data were used for checking; the detailed working condition parameters and calculation results are listed in Table 3, where p is pressure, Qeffrepresent effective heat flux.The error between the experiment and the Gnielinski classical empirical heat transfer formula calculation was within 5 %.

    An experimental data deviation analysis was performed.The uncertainty of the directly measured parameters is the equipment precision. The indirect parameter R is calculated from the many directly measured data,Ri.The calculation process produces a transfer error.In this case,the R uncertainty is calculated using Eq. (9):

    The heat flux is calculated by the measured power voltage(U) and current (I):

    Table 3 Error between experimental data and Gnielinski classical formula.

    The uncertainties of all test parameters are shown in Table 4.

    2.3. Experimental procedure

    Two extremely similar tubes were chosen for the preliminary experiment, and the heat flux of the parallel tubes was almost the same during the experiment.First,the total mass flow rate was decreased (Steps 1-6) to observe Type I excursion, and then the total mass flow rate was increased(Step 7)to observe Type II excursion.Because the total mass flow rate is constant when flow excursion occurs,Type I excursion occurring in one tube results in a sudden mass increase in the other tube. Similarly,Type II excursion occurring in one tube results in a sudden mass decrease in the other tube; thus, the onset of flow excursion can be captured. The detailed experimental process is as follows:

    Step 1.Start the pump and ensure that the volume flow rate of the pump is steady.

    Step 2. Set the back-pressure valve to the target value.

    Step 3.Adjust the preheater AC power to set inlet fuel temperature at the target value.

    Table 4 Uncertainty of experimental parameters.

    Step 4. Adjust the test AC powers to make heat flux to the target value.

    Step 5. When the outlet fuel bulk temperature becomes steady,collect the experimental data,including mass flow rate,heating power, inlet and outlet bulk temperatures, working pressure, pressure drop, and wall temperature. The entire collection period is maintained at ~40 s.

    Step 6.Decrease the total mass flow rate and repeat Step 5.As the outlet bulk temperature becomes close to the pseudocritical temperature, minimize the change step of total mass flow rate and observe Type I excursion occurring in every tube.

    Step 7. Increase the total mass flow rate and repeat Step 5.As the outlet bulk temperature becomes close to the pseudocritical temperature,minimize the change step of the total mass flow rate and observe Type II excursion occurring in every tube.

    Step 8. Change the heat flux, inlet temperature, or system pressure, and repeat Steps 1-7.

    3. Results and discussion

    3.1. Flow process analysis

    To investigate the flow characteristics of flow excursion and uneven flow distribution in parallel tubes, the hydrodynamic multivalue characteristic curve was used to analyze the experimental phenomena. A 4.5 MPa system pressure, 30 °C inlet temperature, and 245 kW/m2heat flux were chosen. Dimensionless time t*is plotted on the X-axis on account of different operating times under every working condition,and t*=tnow/tall(tnowis the real-time corresponding to the experiment and tallis the total experiment time).

    There was a small difference between the parallel tubes in the experiment; thus, the internal flow resistances of the two tubes were not exactly the same.In this study,the small difference between the two tubes was ignored to analyze the flow process based on the hydrodynamic characteristic curves. As the total mass flow rate increases or decreases, the fluid in the tube undergoes the following process. Type I excursions are shown in Fig.5,where the number without‘‘'”represents the data points in Tube 1,and the number with‘‘'”represents the data points in Tube 2. Under the same pressure drop constraint between parallel tubes, for the process from 1(1') to 3(3'), the fluid flow was stable in two tubes, and the mass flow rate in both tubes decreased simultaneously within the right positive slope area of the hydrodynamic curves. As the flow reached the minimum point 3(3'),the mass flow rate in parallel started to exhibit a different trend.The mass flow rate in Tube 1 decreased slightly, whereas it increased slightly in Tube 2.When the flow reached 4(4'), Type I excursion occurred in Tube 1,and a sudden mass flow rate increase occurred in Tube 2. From 4(4') to 5(5'), the mass flow rate in Tube 1 continued to increase, and that in Tube 2 continued to decrease, indicating that the flow distribution became more uneven as the total mass flow rate decreased because Tube 1 was located on the negative slope area and Tube 2 on the right positive slope area.For the process from 5(5') to 6(6'), both mass flow rates decreased as the total flow decreased. When the flow in Tube 2 reached the minimum point 6', the mass flow rate in Tube 2 decreased slightly, while that in Tube 1 increased slightly.Type I excursion occurred in Tube 2 when it came to 7',which led to a sudden flow increase in Tube 1. After this excursion,from 8(8')to 9(9'),the mass flow rate in both tubes maintained the same decreasing trend within the left positive slope area.

    Following the previous step, the total mass flow rate was increased gradually, as shown in Fig. 6. For the process from 1'(1)to 3'(3),the mass flow rate in both tubes increased simultaneously. As the flow reached the maximum point 3'(3), the mass flow rate in both tubes started to exhibit the opposite trend. The mass flow in Tube 1 decreased slightly, while that in Tube 2 increased slightly, which made the flow distribution slightly uneven. When the flow reached 4'(4), Type II flow excursion occurred in Tube 2, and the flow rate decrease occurred in Tube 1. The uneven flow distribution started to increase after point 4'(4).From 4'(4)to 5'(5),the mass flow rate in Tube 2 continued to increase,while that in Tube 1 showed a slight decreasing trend.For the process from 5'(5)to 6'(6),the mass flow rate in both tubes increased simultaneously because of the positive slope area.When the flow in Tube 1 reached the maximum point 6,Type II excursion occurring in Tube 1 led to a sudden flow rate decrease in Tube 2.From 6(6')to 7(7'), the mass flow rate in Tube 1 kept increasing and that in Tube 2 kept decreasing. The degree of uneven flow distribution decreased distinctly. From 7(7') to 8(8'), the mass flow rate in both tubes increased simultaneously, and the flow distribution was quite uniform because of the right positive slope.

    The above-mentioned results of the flow process analysis reveal that the flow path between the decrease and increase in the total mass flow rate is discrepant, which results from the hysteresis phenomenon. The hysteresis phenomenon of flow distribution in this study can be interpreted on the basis of Liu et al.’s research.59Because the hysteresis phenomenon exists in Type I and Type II flow excursions,the flow distribution caused by flow excursion also induces the same phenomenon.Thus,different flow evolutionary processes have different flow distribution results in parallel tubes.

    Fig. 5 Flow analysis with hydrodynamic multivalued characteristic curve in decreasing flow rate process.

    Fig. 6 Flow analysis with hydrodynamic multivalued characteristic curve in increasing flow rate process.

    3.2. Flow characteristic in parallel tubes

    3.2.1. Onset of flow excursion

    According to the classification by Liu et al.,59OFIIcorresponds to the onset of Type I excursion, and OFIIIcorresponds to the onset of Type II excursion. In this section, the characteristics of Type I and Type II flow excursions are analyzed.

    Fig.7 shows the experimental results of different inlet fluid temperatures(Tin)at 4.5 MPa system pressure and 245 kW/m2heat flux. As the flow rate decreased, the fluid outlet temperature was Tout1=281°C at OFII.With the increasing flow rate,the fluid outlet temperature was Tout2=296°C at OFIII.The data imply that at supercritical pressure, flow instability exists when the outlet fluid temperature is approaching the pseudocritical temperature. The drastic property change near the pseudo-critical temperature is one of the elements that causes these phenomena. Density is considered the major parameter of flow instability. The considerable inlet and outlet density difference resulted in the two types of flow excursion and eventually led to flow instability in parallel tubes.

    Fig. 8 shows the experimental results of different system pressures at 30°C inlet fluid temperature and 245 kW/m2heat flux. With the increase in the system pressure, both OFIIand OFIIIincreased. As the pressures increased to supercritical pressures, the variation of the thermal property decreased and the pseudo-critical temperature increased. Thus, OFIIand OFIIImoved to a larger temperature according to a higher system pressure. Increasing the system pressure can improve the flow stability of the fluid at supercritical pressure.

    There is a drastic property change in the two-phase boiling process. Some scholars have defined the pseudo-critical interval to describe the process from a gaseous to a liquid-like state.Drastic property changes occur in the pseudo-critical interval.64-66Banuti64proposed a method for calculating the pseudo-critical interval. It uses the relationship between enthalpy(h) and temperature. Based on this method and the principle shown in Fig. 9, the following parameters were calculated.

    Fig. 7 Effect of inlet fluid temperature on fluid outlet temperature.

    In Fig.9,hiG(T )and hL(T )are the tangent lines of the fluid property curve, which are extremely far from the pseudocritical temperature. Therefore, kiGand h0,iGin Eq. (12) were calculated based on the two points chosen extremely far from the pseudo-critical temperature on the fluid property curve obtained from REFPROP 9.0. Furthermore, kLand h0,Lin Eq. (13) was calculated in the same way. Similarly, kpband h0,pbin Eq. (14) were calculated by choosing two extremely close points near the pseudo-critical temperature on the fluid property curve. Tpbin Eq. (14) is the fluid temperature in the quasi-boiling range.

    Combining Eqs. (12) - (14), the formula for T+and T-is obtained.

    The pseudo-critical interval of cyclohexane under different working conditions was calculated in this study based on the method described above.T-was considered as the outlet fluid temperature of OFIIand T+as the outlet fluid temperature of OFIII.Fig.10 shows the experimental T-and T+.The results of the experimental and computational comparisons are summarized in Table 5. Thus, it can be verified that, as the outlet fluid temperature reaches the onset of the pseudo-critical interval in the total mass flow rate decrease process, Type I excursion occurs.As the outlet fluid temperature reaches the end of the pseudo-critical interval during the process of total mass flow rate increase, Type II excursion occurs.

    3.2.2. Flow distribution ratio

    The flow distribution in the parallel channels had a sharp jump after flow excursion. The flow distribution ratio (Mi) was defined to understand the flow distribution phenomenon more clearly.

    In Eq. (17),the numerator part is the mass flow rate of the two tubes,and the denominator part is the value if the fluid in the two tubes is evenly distributed.Moreover,i=1 and 2 corresponding to the flow deviations of the two tubes.When Miis closer to 1, the system becomes more stable.

    Fig. 8 Effect of system pressure on fluid outlet temperature.

    Fig. 9 Pseudo-critical interval calculation for cyclohexane at 4.2 MPa based on Banuti’s method.

    Fig.11 displays the flow distribution ratio at different inlet fluid temperatures. The system pressure was 4.5 MPa,and the heat flux was 245 kW/m2. The results of the decreased total mass flow rate are shown in Fig. 11(a). The system was stable before the flow excursion, and the flow distribution ratio was almost the same, which suggested that the inlet fluid temperature had negligible influence at this stage. When Type I excursion occurred in Tube 1, there was a sharp change in the flow distribution ratio in the two channels, and the maximum flow distribution difference was reached at the end of the flow excursion. As the Type I excursion appeared in Tube 2, the flow distribution difference between the two channels was reduced. Fig. 11(b) illustrates the flow distribution ratio with increased mass flow rate.In the initial stage,the flow distribution ratio of the two tubes was in a certain range.Similar to the previous process, the flow distribution difference of the two tubes was larger when the Type II excursion occurred in Tube 2. The difference in flow distribution reached its maximum as the flow excursion faded away.The maximum value was different under every inlet fluid condition. When the Type II excursion occurred in Tube 1,there was another sharp change in the flow distribution ratio; in this case, the flow distribution in both tubes was close to 1.

    Fig.10 Outlet fluid temperature of OFII and OFIII corresponding to density under all conditions.

    Table 5 Deviations of pseudo-critical interval between experiment and theoretical calculations.

    Fig. 11 Effect of inlet fluid temperature on flow distribution ratio.

    Overall, with increasing inlet temperature, the range of the flow excursion narrowed, and the maximum flow distribution difference decreased.The temperature gradient along the cross section of the tube became flatter at a higher inlet temperature,and the density variation weakened.The increase in the incompressible high-density fluid decreased the flow excursion intensity and increased the system stability.

    The outcomes for heat fluxes of 205,245,and 299 kW/m2at a supercritical pressure of 4.5 MPa with an inlet temperature of 30 °C are displayed in Fig. 12. The amplitude of the flow distribution ratio increased as the heat flux increased.The flow distribution ratio changed non-monotonically with heat flux.Under the higher and lower heat flux conditions,the difference between M1and M2was ~0.4. In contrast, under the middle heat flux condition, the difference reached 0.8. The nonmonotonic relationship between flow distribution and heat flux suggests that heat flux not only contributes to changing fluid thermal properties but also influences the flow process, which requires further investigation.

    Experiments were conducted at 4.2, 4.8, and 5.5 MPa,30 °C inlet fluid, and 245 kW/m2, and the results are shown in Fig. 13. The values at the different pressures were almost the same. Compared with the fluid inlet temperature and heat flux, the system pressure had a weaker relationship with the flow distribution ratio at supercritical pressure.

    3.3. Boundaries of two types of flow excursion

    For the two-phase fluid, two dimensionless parameters, subcooling number Nsuband phase change number Npch, were defined to obtain the flow instability boundary. The dimensionless sub-pseudo-critical number, NSPC, and trans-pseudocritical number, NTPC, have similarly been defined for supercritical fluids.

    where βpc,cp,pc,hpcrepresent the volumetric expansion coefficient, specific heat capacity at constant pressure, enthalpy at critical pressure, respectively; hinand minrepresent the enthalpy and volume flow rate at the entrance of tube,respectively; NTPCand NSPCare used to obtain the boundaries of OFIIand OFIII. Fig. 14 shows the flow excursion boundary based on the experimental data. The OFIIboundary in the two tubes almost coincided, as for OFIII. This indicates that the instability boundary in a single tube can be used to predict flow excursion occurring in parallel tubes.Based on the experimental data, flow excursion instability boundaries were developed:

    The application ranges of the two formulas are as follows:4.2 MPa<p <5.5 MPa,30°C <Tin<70°C,and 205 kW/m2<q <299 kW/m2.

    Fig. 12 Effect of heat flux on flow distribution ratio.

    Fig. 13 Effect of system pressure on flow distribution ratio.

    In addition, the section between the boundary lines constitutes the entire flow instability area,which means a strong uniform flow distribution may appear in this section.When NSPCis a constant,the variation of NTPCon the two boundary lines can qualitatively determine the flow excursion interval. The flow excursion interval decreases with an increase in the system pressure and inlet temperature, which is consistent with the analysis results in Section 3.2.2. Set the scramjet regenerative cooling system’s design point out of the unstable area can enhanced the system flow stability and safety.

    The discussion shows that supercritical hydrocarbon fuel flow in regenerative cooling parallel channels comes to a severe uneven flow distribution result. Flow distribution uneven degree enhanced as flow excursion occurs inside the instability boundary.The hydrodynamic characteristics should be considered seriously as well as the fuel heat transfer pattern in the scramjet regenerative cooling system design. It’s a crucial factor to system stability and safety.Channels structure optimization can recede the negative slope area of the multi-value curve. Flow uneven distribution triggered by flow excursion restrained with enhanced system pressure and reduced inlet fluid subcooled degree.The improvement measures mentioned above were got relatively optimized results in our follow-up large active regenerative plate cooling design project. As shown in Fig. 15, the channels’ wall temperature was uniform in the experimental research,which indicated an even flow distribution in parallel channels.In addition, mass flow rate coupling with thermal load between parallel tubes is also a promising way to improve system stability.

    Fig. 14 Flow instability boundary based on OFII and OFIII.

    Fig.15 Experimental results of optimized cooling structure with our theoretical analysis.

    4. Conclusions

    Supercritical fluid flow instability is a crucial issue for regenerative active cooling technology. Based on our previous large plate structure design of regenerative cooling system, fundamental investigations on supercritical flow instability in parallel heated tubes were conducted experimentally. New characteristics concerning flow excursion, flow distribution,and OFI were observed. The findings are summarized as follows,and these conclusions have provided theoretical guidance and technical support for the regenerative active cooling design.

    (1)Flow excursion occurring in a single tube can trigger the flow maldistribution between parallel tubes. The system has different distribution results during the process of decreasing and increasing the total mass flow rate, indicating that a hysteresis phenomenon occurs during the flow distribution process because of the different paths of Type I and Type II excursions.

    (2) A higher inlet fluid temperature is favorable for a uniform flow distribution. The effect of system pressure on the flow distribution is limited.The results above mean that higher inlet fluid temperature and system pressure is favorable to a stable flow in the regenerative cooling system of scramjet.The relationship between the heat flux and flow distribution is nonmonotonic in the current heat flux range.

    (3) In view of the pseudo-critical interval proposed in previous research, the outlet fluid temperature at OFIIis close to the left boundary of the interval,and the outlet fluid temperature at OFIIIis close to the right boundary of the interval.Thus, the pseudo-boiling theory is proved to be applicable for interpreting the flow excursion instability of supercritical fluids. Through the calculation of pseudo-boiling temperature interval,the two types excursion boundary can be roughly estimated.The non-uniform flow distribution in parallel tubes can be prevented in the regenerative cooling technology of scramjet.

    (4) The flow excursion instability boundary was predicted based on OFIIand OFIIIusing two dimensionless parameters.The same instability boundary can be approximately used for two parallel tubes, which indicates that the instability boundary obtained using a single tube can be used as a prediction criterion for flow excursion in a parallel-tube system. The flow excursion boundaries via two dimensionless number can be used in reliable working conditions design of scramjet.

    Declaration of Competing Interest

    The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

    Acknowledgements

    This study was co-supported by the Open Fund of Key Laboratory of Power Research of China (No. 2017-III-0005-0029)and the National Natural Science Foundation of China (No.51776167).

    九九久久精品国产亚洲av麻豆 | 国产极品精品免费视频能看的| av天堂中文字幕网| www.自偷自拍.com| 午夜免费观看网址| 51午夜福利影视在线观看| 一级作爱视频免费观看| 99久久成人亚洲精品观看| 99国产精品99久久久久| 免费在线观看影片大全网站| 精品国产美女av久久久久小说| 亚洲精品国产精品久久久不卡| 国产精品 国内视频| 成人性生交大片免费视频hd| 亚洲欧美日韩高清在线视频| 免费无遮挡裸体视频| 国产毛片a区久久久久| 色综合婷婷激情| 国产又黄又爽又无遮挡在线| 一进一出抽搐动态| 热99re8久久精品国产| 久久久久精品国产欧美久久久| 亚洲一区二区三区不卡视频| 午夜福利在线在线| 99精品欧美一区二区三区四区| 日日干狠狠操夜夜爽| 99精品久久久久人妻精品| 成人鲁丝片一二三区免费| 啦啦啦韩国在线观看视频| 亚洲专区国产一区二区| 国产不卡一卡二| or卡值多少钱| 麻豆国产av国片精品| 日本在线视频免费播放| 欧美日韩精品网址| 男人舔女人下体高潮全视频| 两个人视频免费观看高清| 在线免费观看不下载黄p国产 | 神马国产精品三级电影在线观看| 亚洲熟妇熟女久久| 真人一进一出gif抽搐免费| 日本三级黄在线观看| 欧美乱色亚洲激情| 日本与韩国留学比较| 长腿黑丝高跟| 欧美色视频一区免费| 黄色丝袜av网址大全| 国产精品乱码一区二三区的特点| 成人三级黄色视频| 在线观看美女被高潮喷水网站 | 淫秽高清视频在线观看| 亚洲成人久久爱视频| 91老司机精品| 亚洲中文字幕一区二区三区有码在线看 | 狠狠狠狠99中文字幕| 国产成人av教育| 91老司机精品| 国产极品精品免费视频能看的| 久久精品国产亚洲av香蕉五月| 在线观看美女被高潮喷水网站 | 亚洲狠狠婷婷综合久久图片| 黄色片一级片一级黄色片| 男女那种视频在线观看| 成年女人看的毛片在线观看| 亚洲av免费在线观看| 亚洲五月天丁香| 99久久综合精品五月天人人| 亚洲一区二区三区不卡视频| 午夜免费成人在线视频| 亚洲九九香蕉| 欧美丝袜亚洲另类 | 亚洲av第一区精品v没综合| 免费在线观看成人毛片| 99久久99久久久精品蜜桃| 午夜视频精品福利| 欧美午夜高清在线| 亚洲五月婷婷丁香| 亚洲av成人精品一区久久| 变态另类成人亚洲欧美熟女| 国产精品久久久久久精品电影| 免费观看人在逋| 亚洲人与动物交配视频| 少妇的逼水好多| 三级国产精品欧美在线观看 | 久久精品国产综合久久久| 悠悠久久av| 国产精品,欧美在线| 亚洲五月天丁香| 午夜福利成人在线免费观看| 婷婷精品国产亚洲av| 最新中文字幕久久久久 | 欧美日本视频| 欧美大码av| 亚洲自偷自拍图片 自拍| 人人妻,人人澡人人爽秒播| 女生性感内裤真人,穿戴方法视频| 一本一本综合久久| 美女被艹到高潮喷水动态| 亚洲成人久久性| 美女高潮喷水抽搐中文字幕| 久99久视频精品免费| 精品国产乱子伦一区二区三区| 日韩中文字幕欧美一区二区| 国产精品一区二区三区四区久久| 99精品久久久久人妻精品| 亚洲在线自拍视频| 麻豆久久精品国产亚洲av| 国产91精品成人一区二区三区| 免费人成视频x8x8入口观看| 人妻久久中文字幕网| 免费看光身美女| 欧美av亚洲av综合av国产av| cao死你这个sao货| 欧美另类亚洲清纯唯美| 国产精品一区二区三区四区免费观看 | 欧美中文综合在线视频| 成人高潮视频无遮挡免费网站| 在线视频色国产色| 国产精品 国内视频| 国产精品永久免费网站| 身体一侧抽搐| 久久久精品欧美日韩精品| 99精品久久久久人妻精品| 国产探花在线观看一区二区| 日本与韩国留学比较| 免费观看人在逋| 国产单亲对白刺激| 变态另类成人亚洲欧美熟女| 18禁国产床啪视频网站| 午夜免费激情av| 观看免费一级毛片| 91字幕亚洲| 一个人免费在线观看的高清视频| 狂野欧美激情性xxxx| 国内精品一区二区在线观看| 久久中文看片网| 男插女下体视频免费在线播放| 国产麻豆成人av免费视频| 一区二区三区国产精品乱码| 天天躁狠狠躁夜夜躁狠狠躁| 欧美激情在线99| 午夜成年电影在线免费观看| 国内精品美女久久久久久| 麻豆国产97在线/欧美| 一进一出抽搐动态| 亚洲自偷自拍图片 自拍| 啦啦啦免费观看视频1| 午夜亚洲福利在线播放| 级片在线观看| 波多野结衣高清无吗| 国产精品爽爽va在线观看网站| 久久这里只有精品19| 9191精品国产免费久久| 国产主播在线观看一区二区| 两性夫妻黄色片| 亚洲人成网站高清观看| 少妇熟女aⅴ在线视频| 波多野结衣高清无吗| 国产 一区 欧美 日韩| 久久热在线av| 久久香蕉精品热| 国产蜜桃级精品一区二区三区| 欧美黑人巨大hd| 国产精品久久电影中文字幕| 欧美绝顶高潮抽搐喷水| 99久久久亚洲精品蜜臀av| 国产欧美日韩一区二区三| 亚洲性夜色夜夜综合| 亚洲人成电影免费在线| 精品国内亚洲2022精品成人| 成人特级黄色片久久久久久久| 久久这里只有精品中国| 国产一区二区三区在线臀色熟女| 18禁黄网站禁片午夜丰满| 欧美成狂野欧美在线观看| 色av中文字幕| 少妇的丰满在线观看| www.精华液| 啦啦啦韩国在线观看视频| 999精品在线视频| 国产伦精品一区二区三区视频9 | 99久久精品一区二区三区| 亚洲午夜理论影院| 国产精品久久久久久久电影 | 国产主播在线观看一区二区| 露出奶头的视频| 日韩高清综合在线| 三级毛片av免费| 亚洲欧洲精品一区二区精品久久久| 欧美高清成人免费视频www| 亚洲中文字幕日韩| 精品国产三级普通话版| 国产精品自产拍在线观看55亚洲| 欧美日韩亚洲国产一区二区在线观看| 色哟哟哟哟哟哟| 色综合站精品国产| 亚洲男人的天堂狠狠| 夜夜看夜夜爽夜夜摸| 国产一区二区在线av高清观看| 色视频www国产| 亚洲午夜精品一区,二区,三区| 国产私拍福利视频在线观看| 99久久国产精品久久久| 一个人看的www免费观看视频| 黄片小视频在线播放| 国产精品美女特级片免费视频播放器 | 男人舔女人下体高潮全视频| 在线观看免费午夜福利视频| 无限看片的www在线观看| 99久久精品热视频| 老司机午夜十八禁免费视频| 国产黄片美女视频| 视频区欧美日本亚洲| 99久国产av精品| 三级男女做爰猛烈吃奶摸视频| 国产探花在线观看一区二区| www.自偷自拍.com| 国模一区二区三区四区视频 | 搡老妇女老女人老熟妇| 一级毛片女人18水好多| 精品久久久久久,| 18美女黄网站色大片免费观看| 国产探花在线观看一区二区| 18禁黄网站禁片免费观看直播| 又大又爽又粗| 午夜日韩欧美国产| 99热只有精品国产| 在线免费观看不下载黄p国产 | 男人和女人高潮做爰伦理| 丁香六月欧美| 啦啦啦免费观看视频1| 午夜福利视频1000在线观看| 亚洲国产欧美网| avwww免费| 亚洲欧美日韩高清专用| 老司机深夜福利视频在线观看| 啦啦啦免费观看视频1| 给我免费播放毛片高清在线观看| 又黄又粗又硬又大视频| 99精品欧美一区二区三区四区| e午夜精品久久久久久久| 国产精品1区2区在线观看.| 19禁男女啪啪无遮挡网站| 亚洲精品456在线播放app | bbb黄色大片| 久久久久国产精品人妻aⅴ院| 久久久久久久久久黄片| x7x7x7水蜜桃| 亚洲人成网站在线播放欧美日韩| 日本免费a在线| 最新美女视频免费是黄的| 婷婷精品国产亚洲av在线| 欧美黑人巨大hd| 看片在线看免费视频| 国产真人三级小视频在线观看| 欧美成狂野欧美在线观看| 国产伦一二天堂av在线观看| 九色成人免费人妻av| 狂野欧美激情性xxxx| 精品一区二区三区视频在线观看免费| 亚洲欧美一区二区三区黑人| 啦啦啦韩国在线观看视频| 亚洲成人中文字幕在线播放| 精品一区二区三区视频在线 | 色噜噜av男人的天堂激情| 精品欧美国产一区二区三| 国产熟女xx| 18禁国产床啪视频网站| 成人三级做爰电影| 午夜福利在线观看免费完整高清在 | 亚洲国产精品久久男人天堂| 高潮久久久久久久久久久不卡| 99久久久亚洲精品蜜臀av| 国产一区二区激情短视频| 黄色片一级片一级黄色片| 精品国内亚洲2022精品成人| 国产精品影院久久| 99国产精品一区二区蜜桃av| 午夜成年电影在线免费观看| 国产亚洲av高清不卡| 免费av毛片视频| 十八禁人妻一区二区| 国产成人啪精品午夜网站| 免费看光身美女| 国产一区二区在线观看日韩 | 我的老师免费观看完整版| 亚洲美女黄片视频| 亚洲av免费在线观看| 日韩高清综合在线| 成人鲁丝片一二三区免费| 丰满人妻熟妇乱又伦精品不卡| 久久久久久久久免费视频了| 熟妇人妻久久中文字幕3abv| 999久久久国产精品视频| av欧美777| 中文字幕人妻丝袜一区二区| 19禁男女啪啪无遮挡网站| 国产精品日韩av在线免费观看| 国产精品自产拍在线观看55亚洲| 亚洲av熟女| 欧美xxxx黑人xx丫x性爽| 亚洲成av人片免费观看| 99久久精品国产亚洲精品| 欧美一级毛片孕妇| 国产激情久久老熟女| 精品国产乱子伦一区二区三区| 亚洲国产日韩欧美精品在线观看 | 欧美丝袜亚洲另类 | 国产真实乱freesex| 性欧美人与动物交配| 制服人妻中文乱码| 特级一级黄色大片| 1000部很黄的大片| 免费看a级黄色片| 国产精品久久久久久精品电影| 国产成人福利小说| 99国产精品一区二区蜜桃av| 黄片大片在线免费观看| 成人鲁丝片一二三区免费| 精品久久蜜臀av无| 亚洲在线自拍视频| 91在线观看av| 老熟妇乱子伦视频在线观看| 每晚都被弄得嗷嗷叫到高潮| 在线视频色国产色| 露出奶头的视频| 久久久久九九精品影院| 久久久色成人| 久久草成人影院| 国产在线精品亚洲第一网站| 久久精品亚洲精品国产色婷小说| 在线观看免费视频日本深夜| 天堂av国产一区二区熟女人妻| 亚洲自偷自拍图片 自拍| 亚洲真实伦在线观看| 久久天堂一区二区三区四区| av在线天堂中文字幕| 亚洲五月婷婷丁香| 日本在线视频免费播放| 99国产精品99久久久久| av福利片在线观看| 日本五十路高清| 日韩欧美三级三区| 91av网站免费观看| 桃红色精品国产亚洲av| 无人区码免费观看不卡| 法律面前人人平等表现在哪些方面| 观看免费一级毛片| 成人高潮视频无遮挡免费网站| 亚洲成人久久性| avwww免费| 极品教师在线免费播放| 日本免费a在线| 天天添夜夜摸| 最近最新中文字幕大全免费视频| 99国产综合亚洲精品| 精品久久久久久久末码| 99热6这里只有精品| www.熟女人妻精品国产| 久久精品国产综合久久久| 亚洲 欧美 日韩 在线 免费| av国产免费在线观看| 国产av不卡久久| 老鸭窝网址在线观看| x7x7x7水蜜桃| 最新美女视频免费是黄的| 无限看片的www在线观看| 亚洲一区二区三区色噜噜| 观看美女的网站| 亚洲av熟女| 黄片大片在线免费观看| 欧美中文日本在线观看视频| 脱女人内裤的视频| 欧美日韩瑟瑟在线播放| 又粗又爽又猛毛片免费看| 国产精品久久久人人做人人爽| 亚洲九九香蕉| 久久久色成人| 国产精品 欧美亚洲| 精品久久久久久久久久免费视频| 亚洲午夜精品一区,二区,三区| av在线蜜桃| 国产免费av片在线观看野外av| 国产亚洲欧美在线一区二区| 又大又爽又粗| 熟女少妇亚洲综合色aaa.| 中文资源天堂在线| 两个人的视频大全免费| 男女做爰动态图高潮gif福利片| 嫩草影院精品99| 久久久色成人| svipshipincom国产片| 国产成+人综合+亚洲专区| 91字幕亚洲| 国产午夜精品久久久久久| 丝袜人妻中文字幕| 色视频www国产| 亚洲欧美日韩高清在线视频| 成年女人看的毛片在线观看| 日韩欧美一区二区三区在线观看| 国产一区二区在线av高清观看| 免费观看的影片在线观看| 午夜福利成人在线免费观看| 国产成+人综合+亚洲专区| 日韩免费av在线播放| 亚洲av免费在线观看| 好男人在线观看高清免费视频| 亚洲中文日韩欧美视频| 男女下面进入的视频免费午夜| 国产精品亚洲一级av第二区| 深夜精品福利| 成人鲁丝片一二三区免费| 亚洲人成网站在线播放欧美日韩| 国产欧美日韩精品亚洲av| 少妇的丰满在线观看| 国产精品综合久久久久久久免费| 两个人视频免费观看高清| 国产亚洲精品综合一区在线观看| 黄片小视频在线播放| 色综合欧美亚洲国产小说| 12—13女人毛片做爰片一| 在线观看午夜福利视频| 无限看片的www在线观看| 国产激情偷乱视频一区二区| 国产极品精品免费视频能看的| 亚洲国产精品成人综合色| 中文字幕高清在线视频| 免费电影在线观看免费观看| 国产成人欧美在线观看| 精品一区二区三区视频在线观看免费| 我要搜黄色片| 看黄色毛片网站| 国产欧美日韩一区二区三| av女优亚洲男人天堂 | 中文字幕人成人乱码亚洲影| 亚洲自偷自拍图片 自拍| 久久久久国产精品人妻aⅴ院| 成人性生交大片免费视频hd| 大型黄色视频在线免费观看| 精华霜和精华液先用哪个| 一个人免费在线观看电影 | 中文资源天堂在线| 91av网一区二区| 啪啪无遮挡十八禁网站| 欧美中文综合在线视频| av福利片在线观看| 欧美黄色片欧美黄色片| 国产欧美日韩一区二区精品| 国产av麻豆久久久久久久| 欧美一区二区精品小视频在线| 国产一区二区三区在线臀色熟女| 国产亚洲欧美98| 757午夜福利合集在线观看| 色综合亚洲欧美另类图片| 18禁裸乳无遮挡免费网站照片| 看免费av毛片| 久久久久亚洲av毛片大全| 免费大片18禁| www.精华液| 亚洲电影在线观看av| 九九热线精品视视频播放| 国产日本99.免费观看| 日韩精品中文字幕看吧| 午夜免费观看网址| av女优亚洲男人天堂 | 91麻豆精品激情在线观看国产| 亚洲国产色片| 欧美日韩国产亚洲二区| 草草在线视频免费看| 日本一二三区视频观看| 欧美乱码精品一区二区三区| 丁香欧美五月| 国产亚洲精品久久久com| 免费一级毛片在线播放高清视频| 国产野战对白在线观看| 免费搜索国产男女视频| 精品国产乱子伦一区二区三区| 国产精品影院久久| 中国美女看黄片| 9191精品国产免费久久| 国产精品一及| 色综合欧美亚洲国产小说| 黑人巨大精品欧美一区二区mp4| 香蕉丝袜av| 天堂av国产一区二区熟女人妻| 久久天堂一区二区三区四区| 日本撒尿小便嘘嘘汇集6| 成人特级黄色片久久久久久久| 国产成人精品无人区| 亚洲午夜理论影院| 天堂√8在线中文| 亚洲精品乱码久久久v下载方式 | 久久人人精品亚洲av| 亚洲国产色片| 999久久久精品免费观看国产| 日日干狠狠操夜夜爽| 亚洲国产欧美网| 亚洲熟妇熟女久久| 亚洲第一电影网av| 亚洲国产精品sss在线观看| 欧美日韩福利视频一区二区| 黄色成人免费大全| 国产一区二区三区在线臀色熟女| 国产人伦9x9x在线观看| 午夜福利欧美成人| 免费在线观看成人毛片| 国产高清视频在线播放一区| 中亚洲国语对白在线视频| 天天一区二区日本电影三级| 叶爱在线成人免费视频播放| 少妇的逼水好多| 天堂动漫精品| 欧美大码av| 在线永久观看黄色视频| 成熟少妇高潮喷水视频| 人人妻,人人澡人人爽秒播| 在线观看日韩欧美| 在线观看66精品国产| 男女之事视频高清在线观看| 免费观看精品视频网站| 国产主播在线观看一区二区| 亚洲一区高清亚洲精品| av中文乱码字幕在线| 久久久国产成人免费| 全区人妻精品视频| 日日干狠狠操夜夜爽| 看片在线看免费视频| 亚洲av免费在线观看| 国产在线精品亚洲第一网站| 久久久成人免费电影| 午夜福利成人在线免费观看| 神马国产精品三级电影在线观看| 欧美乱妇无乱码| 免费av不卡在线播放| 男女之事视频高清在线观看| 久久精品国产综合久久久| 久久久久久国产a免费观看| 99久久国产精品久久久| 熟女人妻精品中文字幕| 久久国产精品影院| 亚洲精品中文字幕一二三四区| 成人永久免费在线观看视频| 久久久精品欧美日韩精品| 熟女少妇亚洲综合色aaa.| 婷婷精品国产亚洲av在线| 一本久久中文字幕| 国产精品98久久久久久宅男小说| 夜夜爽天天搞| 老司机深夜福利视频在线观看| 99国产综合亚洲精品| 欧美成人一区二区免费高清观看 | 亚洲欧美精品综合一区二区三区| 亚洲欧美日韩卡通动漫| 黑人操中国人逼视频| 欧美一区二区国产精品久久精品| 最近最新中文字幕大全电影3| 国产精品久久久久久精品电影| svipshipincom国产片| 亚洲美女视频黄频| 日本黄色片子视频| 又粗又爽又猛毛片免费看| av女优亚洲男人天堂 | 两个人视频免费观看高清| 国产亚洲精品综合一区在线观看| 久久久久久国产a免费观看| 成人欧美大片| 美女高潮喷水抽搐中文字幕| 丰满人妻熟妇乱又伦精品不卡| www日本黄色视频网| 久久精品国产99精品国产亚洲性色| 午夜福利在线在线| 别揉我奶头~嗯~啊~动态视频| 两性午夜刺激爽爽歪歪视频在线观看| 久久久国产精品麻豆| 国产激情偷乱视频一区二区| 99在线人妻在线中文字幕| 中文字幕最新亚洲高清| 免费av毛片视频| 精品午夜福利视频在线观看一区| 欧美av亚洲av综合av国产av| 美女 人体艺术 gogo| 国产精品综合久久久久久久免费| tocl精华| 日本黄大片高清| 后天国语完整版免费观看| 亚洲电影在线观看av| 18禁裸乳无遮挡免费网站照片| 变态另类成人亚洲欧美熟女| 久久精品91无色码中文字幕| 婷婷精品国产亚洲av| 69av精品久久久久久| 曰老女人黄片| av视频在线观看入口| 欧美乱色亚洲激情| 亚洲国产欧美人成| 99精品在免费线老司机午夜| 宅男免费午夜| 桃色一区二区三区在线观看| 男人舔女人的私密视频| 欧美在线一区亚洲| 国产激情久久老熟女| 最新中文字幕久久久久 | 中文资源天堂在线| 琪琪午夜伦伦电影理论片6080| av黄色大香蕉| 美女免费视频网站| 欧美色欧美亚洲另类二区| 久久国产精品影院| 亚洲国产精品成人综合色| 国产成+人综合+亚洲专区| a级毛片a级免费在线| 久久久精品大字幕| 久久精品国产清高在天天线| www.精华液| 午夜福利欧美成人|