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

    Numerical evaluation of acoustic characteristics and their damping of a thrust chamber using a constant-volume bomb model

    2018-04-19 08:29:04JianxiuQINHuiqiangZHANGBingWANG
    CHINESE JOURNAL OF AERONAUTICS 2018年3期

    Jianxiu QIN,Huiqiang ZHANG,Bing WANG

    School of Aerospace Engineering,Tsinghua University,Beijing 100084,China

    1.Introduction

    Combustion instabilities have been experienced in many development programs of liquid rocket engines,which are manifested as transient large-amplitude pressure oscillations in the thrust chamber exceeding around 10%of the averaged chamber pressure.Based on the characteristic frequency of pressure oscillations,combustion instabilities can be classified as low-frequency(chugging),intermediate-frequency(buzzing),and high-frequency(acoustical)modes.1High-frequency combustion instability generally results from thermalacoustic coupling,which is the most destructive for its extensive damage to the chamber and injector face.1,2

    Research has shown that high-frequency combustion instability can generally occur when oscillatory energy supplied by unsteady combustion heat release is sufficiently greater than the loss of oscillatory energy damped in the chamber.1Acoustic combustion instability can be retrained by increasing damping,as well as by decreasing or breaking down thermalacoustic coupling for a concerned acoustic mode.Thus,it is important to determine the acoustical modes and their damping characteristics for a thrust chamber.The popular method for predicting acoustic characteristics of a combustion chamber is to solve linear or nonlinear pressure wave equations,3–7in which acoustic modes can be obtained by the Helmholtz analysis.For a linear model,the pressure disturbance grows without a bound and is unphysical,while nonlinear models can obtain the limit cycle behavior in real situations.5–7However,both linear and nonlinear models require an accurate definition of the mean flow and response functions,which are limited to apply to cases with complex chamber geometries.The Finite Element Method(FEM)can be employed to analyze acoustic fields in complex geometries by solving the Helmholtz equation with a Neumann boundary condition.8–10However,it is difficult for the FEM to handle the flow and combustion process occurring in a chamber,which in fact has an important effect on acoustic characteristics.

    The detailed processes of flow,combustion,and propagation of pressure wave can be solved based on the Computational Fluid Dynamics(CFD),which is somehow superior to the Helmholtz analysis of the FEM.Both linear and nonlinear disturbances can be handled in this method,which can also predict the effects of geometry variations and mean flow on acoustic characteristics.Using the CFD method,pressure oscillations have to be excited in a chamber first,and then acoustic characteristics can be obtained by analyzing the recorded time series of the pressure oscillations.Therefore,artificial disturbance models are usually employed in CFD to excite pressure oscillations,which have been developed by different researchers known as models of oscillating velocity disturbance,11–13mass flow rate disturbance,14energy disturbance15,and pressure disturbances16–22as well as pressure pulse.20–22They have then been successfully applied to analyze chamber acoustic properties and acoustic damping of various passive devices by researchers such as Abdelkader16and Taro et al.,17and to stimulate combustion instability by Kim,18Grenda,19Habiballah,20Zhuang,21and Urbano et al.22However,it must be indicated that artificial disturbance is not indispensable for self-sustained combustion instabilities in a combustion chamber predicted by CFD.23–25

    In summary,there are two main approaches to obtain pressure oscillations in a combustion chamber using the CFD method:one is to excite decaying oscillations by an artificial disturbance model,and the other is to realize self-excited combustion instability without artificial perturbations.Based on such pressure oscillations,the acoustic characteristics and the phenomenon of combustion instability in a combustion chamber can be investigated numerically.However,because it is hard to achieve self-excited combustion instability,artificial disturbances are more feasible to evaluate the acoustic characteristics and their damping in a thrust chamber,and often used in industries.For the existing artificial disturbance models,such as the models of mass flow rate,energy,velocity,and pressure,a specific acoustic mode is needed,and a resonantmode oscillation is then trigged.The amplitudes of pressure pulse models are limited to a small value,which are capable to excite a specific acoustic mode,usually a single-frequency mode.In practice,an artificial disturbance method that can simultaneously excite the multi-acoustic modes of a chamber is a pressing need to determine which acoustic mode is the easiest to be excited and which one is the most difficult to be attenuated in the acoustic modes in a given chamber.Therefore,an artificial model must be applied with a large-amplitude pressure pulse and then drive the nonlinear acoustics process of pressure propagations.However,all the existing disturbance models do not possess all of the above traits,which will be further analyzed in Section 2 of the present paper.

    Thus,a numerical constant-volume bomb model for a highamplitude pressure pulse is proposed in this paper.The nonreactive turbulent flows in a chamber are numerically simulated by CFD.Multi-frequency pressure oscillations are excited by the present numerical constant-volume bomb model imposed on a limited region in the chamber.Acoustic modes and their corresponding damping characteristics are analyzed by both the Fast Fourier Transformation(FFT)analysis and the half-power bandwidth method.Furthermore,the effects of the forced geometrical regions on the excited pressure oscillations are further discussed.

    2.Numerical method

    Acoustic properties in a typical configuration of a thrust chamber in a liquid rocket engine are numerically investigated in this paper.The chamber is initially filled with quiescent air with a temperature of 298 K and a pressure of 0.1 MPa.A pressure pulse is then imposed on the mean pressure in a small region,and then fluid flows and pressure oscillations in the chamber are induced.Such turbulent flows with pressure oscillations are numerically simulated by using Unsteady Reynolds-Averaged Navier-Stokes(URANS)equations,and the acoustic properties of the chamber can be obtained by analyzing the excited pressure oscillations.Based on the k-ε twoequation turbulent model,the governing equations in a general form for the turbulent flows in the chamber are written as follows23:

    where φ,Γφ,Sφand ujrepresent the conservation variables,the convective flux coefficient,the source terms,and the velocity in the jth direction,respectively;ρ,t,xjdenote density,time and coordinate axis in the jth direction in the Cartesian coordinate system,respectively.These variables are shown in details in Table 1.

    In Table 1,p,T,k,ε,and Yirepresent pressure,temperature,turbulent kinetic energy,turbulent dissipation rate and mass fraction of the ith species,respectively;λ denotes heat conductivity coefficient,while cpis specific heat at constant pressure;μedenotes the effective viscous coefficient,which contains the laminar viscous coefficient μ and the turbulentviscous coefficient μt.The detailed expressions are shown as follows:

    Table 1 Variables,coefficients,and source terms in Eq.(1).

    where A1and A2are constants,taken as 1.457×10-5and 110,respectively.

    The expression of Gkis shown as follows:

    The other parameters are turbulence model coefficients,shown in Table 2.

    The governing equations are discretized by the finite volume method,in which the diffusion and convection terms are discretized by a second-order central difference scheme and a second-order upwind scheme,respectively.The numerical dispersion as well as numerical dissipations are minimized by these schemes.The Semi-Implicit Method for Pressure-Linked Equations(SIMPLE)is employed to solve the discretized equations.26Time marching is approximated by an explicit first-order Eulerian scheme.To guarantee a computation stability,the time step is specified as 5×10-7s,but it is set as 1×10-8s when a constant-volume bomb is imposed on the steady flow in the chamber.The grid resolution is specified as 1 mm in the simulation.Further refinement of the mesh does not improve the prediction results of eigenfrequencies of the acoustic mode and their damping characteristics significantly.

    Adiabatic and non-slip conditions are adopted to the injector face and sidewall of the thrust chamber.An open outletcondition is applied for the nozzle exit.For that,the outlet pressure is set as the ambient pressure,while the velocity components are set equal to those of the logical inside neighbor vertexes.A choked condition and supersonic flow are not achieved in the nozzle for the cases investigated in this paper,and an open outlet condition means that the pressure at the nozzle exit is set as the ambient pressure.The turbulent flow and pressure oscillations described by the above models are induced by the initial pressure pulse,so how to realize such a suitable pressure pulse is the key problem,which is described in the following section.

    Table 2 Turbulence model coefficients.

    3.Artificial constant-volume bomb model

    Firstly,the existing oscillating velocity disturbance,11–13mass flow rate disturbance,14energy disturbance,15and pressure disturbance16–18as well as the pressure pulse20–22are compared,which are shown in Table 3.The variables in Table 3 are explained in detail in the corresponding references.

    These models can be classified as two types.One is organized disturbance models,in which resonant acoustic frequencies are requested in prior such as velocity disturbance,11–13mass flow rate disturbance,14energy disturbance,15and organized pressure disturbance.16,18Only the resonant acoustic mode can be excited in most of the cases.Thus,the characteristics on which mode is the easiest to be excited and which mode is the most difficult to be damped cannot be determined.Moreover,it is hard to get the acoustics frequencies of the chamber in advance for a complex geometry.The other is pulse disturbance models,such as those proposed by Habiballah,20Zhuang,21and Urbano et al.22The overpressure coefficient(αpin the Habiballah’s pressure pulse model)is defined as the ratio of the pulse amplitude to the mean chamber pressure,usually taken as a small value in these models.

    Taking the pressure pulse model developed by Habiballah and Lourme20as an example,it is applied to a combustion chamber with a pressure of 1.15 MPa and a temperature of 3000 K.A pressure pulse is imposed on the region at the chamber axis and near the head of the combustor with a duration of 0.1 μs,and its overpressure coefficient is taken in a range from 0.1 to 1.5.As shown in Fig.1,a pressure peak appears while the density keeps constant,and a significant temperature peak is therefore induced.The larger the overpressure coefficient of the pressure pulse is,the higher the temperature peak is.It can be seen that the temperature peak is greater than 5000 K when the amplitude of the pressure pulse takes a value of 1.5 times of the mean chamber pressure.It is unreasonable for such a high local temperature,and dramatically leads to a divergent computation with an increasing overpressure coefficient.

    Such a pressure pulse is imposed on the flow field in a small region.It propagates spherically at the initial stage,and then transforms to propagate in specific directions through re flections on the chamber wall,and finally multi-mode acoustic pressure oscillations are formed.Two factors are very important for this process.One is that the pressure waves are still strong enough at the end of the initial propagation stage,and the other is that after the transformation to specific modes,the pressure oscillations can maintain long enough to observe their decay.It is very important to increase the initial peak value of the pressure pulse.However,it is difficult to excite multi-mode acoustic pressure oscillations using the existingpressure disturbance model due to the limited low value of the overpressure coefficient,especially when it is applied to a real case of the combustion chamber.Inconsistencies of pressure and temperature may happen in the imposed regions for the existing model,which can be explained as that such a pressure pulse is realized based on a constant-volume or constantdensity process.Energy in a form of pressure is introduced to the flow field without any carrier mass,so that only a small overpressure coefficient can be acceptable.

    Table 3 Artificial disturbance models in references.

    Fig.1 Time histories of pressure,temperature,and density in imposed regions with Habiballah’s pressure pulse model.

    Then,a model of the numerical constant-volume bomb is therefore proposed to achieve a large-amplitude disturbance.Differed from the existing models,the corresponding mass as a carrier of the pressure pulse is also introduced in the imposed regions.The model is applied as the following based on the CFD scheme:

    (1)An imposed region Ω(x,y,z)is determined at a time instant t0,usually according to the research purpose,and the pressure,temperature,and density are denoted as p0(x,y,z),T0(x,y,z),and ρ0(x,y,z),respectively;

    (2)A gas that has the same components as in the original gas in the region Ω(x,y,z)with pressure γp0(x,y,z),temperature αT0(x,y,z),and density βρ0(x,y,z)is determined in the same region Ω(x,y,z);

    (3)The imposed and original gases are immediately mixed in the constant-volume region Ω(x,y,z).The tuned temperatureT(x,y,z),pressurep(x,y,z),and density ρ(x,y,z)in the same region are obtained as below:

    (4)The temperature,pressure,and density of the imposed region Ω(x,y,z)are set as Eq.(6)to achieve a largeamplitude disturbance from the initial time t0to time t0+Δt in numerical simulations.

    Here α,β,and γ= αβ are model coefficients,which determine the peak value of the pressure pulse.t0is the initial imposed time,while Δt is the imposed time period.This model can realize a pressure pulse with a large overpressure coefficient,which is verified to excite pressure oscillations corresponding to multi-frequency acoustics modes.We call the present model as a numerical constant-volume bomb model,because the imposed and original gases in the specified region are supposed to be immediately mixed in the unchanged constant-volume region,and such a pressure pulse with a large amplitude is achieved.

    The pressure,temperature,and density in the imposed regions for α =1, γ=0.5–9,and β = γ are shown in Fig.2.Different from the case shown in Fig.1 where γ is taken as 1.5,the over-temperature phenomenon does not appear in the present model even for the overpressure coefficient γ=9.A sharp rise is observed on the density,which is characterized by the constant-volume mixing process,so the present model is much more reasonable with an advantage of achieving a pressure pulse with a very large amplitude.

    The model coefficients γ and β (or α),especially the overpressure coefficient γ,are now in posterior for this study.The values of the model coefficients are applicable if multimode acoustic pressure oscillations can be excited.The larger the overpressure coefficient is,the more easily the evolution of multi-mode acoustic pressure oscillations is triggered.The selected value of γ yields to the below limitation:no divergence of the flow field in numerical simulations and competition over a possible damping of the flow field.Therefore,the parameter can be estimated during the implementation of numerical simulations,which is several times or dozens of the mean pressure of the engine chamber.Therefore,γ can be determined inde-pendently based on how much the amplitude of the pressure pulse is required to excite multi-mode pressure oscillations,while β is used to adjust the temperature in the imposed region.

    Fig.2 Time histories of pressure,temperature,and density in imposed regions with constant-volume bomb model.

    4.Quantification of damping capacity of chamber

    Pressure oscillations can be excited by the above numerical constant-volume bomb,which can be expressed as27

    where An,maxis the initial maximum amplitude of each mode,αnis the corresponding damping rate,fnis the acoustic resonant frequencies,and φnis the initial phase of each mode.

    Based on the results of FFT analysis,the damping factor for a peak frequency fn,peakis defined by the half-power bandwidth method as follows:2,10

    where fn,peakis the resonance frequency of each mode,fn,1and fn,2are the frequencies(here,fn,2>fn,1),of which the corresponding amplitudes areis the amplitude of the frequency fn,peak).Once ηnis determined,the damping rate αncan be obtained by the following relation:28,29

    where Δfn=fn,2-fn,1.

    A broadened bandwidth indicates a higher damping factor,followed by a larger damping rate.2,10Thus,the damping factor and half-power bandwidth have the same physical interpretation of acoustic damping as the damping rate.2,30The damping factor ηnand half-power bandwidth Δfnof each mode correspond to the damping rate of that mode.

    Through FFT analyzing the recorded excited pressure oscillations,multi-acoustic modes can be obtained by identifying and matching the peak frequencies with the corresponding theoretical acoustic eigenfrequencies.The acoustic mode that is the easiest to be excited can be determined based on the amplitudes of the peak frequencies;the acoustic mode that is the most difficult to be damped can be determined through the damping factor ηnand half-power bandwidth Δfnfor that mode.Therefore,the acoustics characteristics and their damping of the chamber can be compared and evaluated after the numerical constant-volume bomb is applied.

    5.Results and discussion

    5.1.Validation and verification of numerical methods

    Firstly,the numerical dissipation is investigated for a cylindrical chamber with a diameter of 32 mm and a length of 100 mm.The chamber has two closed ends,and the wall condition is slip and adiabatic.The grid-resolution is taken as 1 mm here.The chamber is filled with non-viscous quiescent air with a pressure of 0.1 MPa and a temperature of 298 K.A forcing is provided with a frequency of the first longitudinal mode 1735 Hz and an amplitude of 400 Pa.After the forcing is switched off,the pressure oscillation in the chamber is recorded.The damping is theoretically zero for such a case.31The amplitude decays from 399.1 Pa to 305.8 Pa during 19 periods.Thus,the damping rate is estimated about 0.104 s-1.The corresponding half-power bandwidth is 0.033 Hz,which is sufficiently smaller than that in a thrust chamber in next section.It proves that the numerical dissipation is well controlled in the present numerical simulations through a second-order discrete scheme and fine grids.

    An acoustic test under an ambient condition for a smallthrust Liquid Rocket Engine(LRE)chamber is conducted as shown in Fig.3.The diameter of the cylindrical chamber is 2rc=32 mm,while that of the throat is 2r*=14.3 mm.The length of the cylindrical section Lchis 38.4 mm,while that of the contraction section Lcvis 27.9 mm.The injector face is a solid wall,while the chamber exit is open.A mimicked bomb is placed near the sidewall,5 mm away from the injector face,while the acoustic sensor is near the chamber wall on the opposite side to the bomb.The time trace of dynamical pressure is recorded by a sensor after the bomb is detonated.The bomb in the test is realized by an explosion of gunpowder.

    Fig.3 Schematic diagram of experimental chamber.

    Table 4 Comparison of theoretical,experimental and numerical results.

    Numerical simulation is also conducted for an experimental test case to verify the above numerical methods.The computational region is chosen from the injector face to the throat plane.The chamber is initially filled with quiescent air with a pressure of 0.1 MPa and a temperature of 298 K.The eigenfrequencies of acoustic modes for this chamber are approximately calculated theoretically,in which the sound speed is taken as 347 m/s.The theoretical,experimental and numerical results are compared and shown in Table 4,in which 1L1T represents the combinations of the First Tangential(1T)and First Longitudinal(1L)mode.

    An artificial constant-volume bomb with γ=19 and β = γ is imposed into the chamber for 0.1 μs in a cylindrical region with a diameter of 4 mm and a length of 3 mm.The imposed position is taken as the same as that in the experimental case.After the constant-volume bomb is switched off,pressure oscillations are recorded in the opposite direction of the imposed region and near the chamber sidewall,which is the same position as that of an acoustic sensor fixed in the experiments.The corresponding FFT analysis results of pressure oscillations are shown in Fig.4.Several peak frequencies are observed as 2950 Hz,6545 Hz,and 7376 Hz in the experiments as well as 3300 Hz,7000 Hz,and 8100 Hz in numerical simulations.Compared to the theoretical frequencies shown in Table 4,2950 Hz in the experimental case and 3300 Hz in the numerical case are identified as 1L mode,while 6545 Hz in the experimental case and 7000 Hz in the numerical case are identified as 1T mode.7376 Hz and 8100 Hz are 1L1T mode.The frequencies of acoustic modes obtained by the simulations and experiments are a little bit higher than those of the theoretical acoustic eigenfrequencies due to the effects of the converging section of the thruster on acoustic characteristics,which has also been verified by Farshichi et al.8The relative errors of eigenfrequencies between the experimental and simulation cases are less than 11%,which may be caused by the effects of holes for equipping with the acoustic sensors and the excitation device,the initial pressure peak,and the outlet boundary conditions.In general,it can be concluded that multi-mode acoustic pressure oscillations are stimulated by the present numerical constant-volume bomb model,and the resonant acoustic frequencies obtained by the present simulations agree well with the experimental and theoretical results.

    The half-power bandwidth method is used to evaluate the damping characteristics.For the 1T acoustic mode,the halfpower bandwidths obtained in the numerical simulations and experiments are 146 Hz and 159 Hz,respectively,while the damping factors are 2.23%and 2.27%,respectively.For the 1L acoustic mode,the half-power bandwidths obtained in the numerical simulations and experiments are 130 Hz and 50 Hz,respectively,while the damping factors are 1.69%and 3.94%,respectively.The outlet boundary condition really has an effect on the eigenfrequencies and their damping rates.The chamber exit is open in the experiments,while the pressure is set as the ambient pressure and velocity components are set equal to those of the logical inside neighbor vertexes at the nozzle exit in the simulations.There are differences on the outlet condition between the experimental and numerical cases,and such differences have effects on the longitudinal acoustic modes.That is the reason why there are obvious differences for the damping rate of the 1L acoustic mode between the experimental and numerical cases.In order to avoid such differences on the outlet boundary between the experimental and numerical cases,a large region downstream from the nozzle exit should be included in the simulations,and the pressure boundary condition can be set on the boundary of this additional region.However,it will induce much expensive computational consumption.The main purpose of this paper is to introduce the numerical bomb,which can excite multi-mode acoustic pressure oscillations.Then the acoustic characteristics of the combustion chamber is investigated.Therefore,we do not try to avoid such differences on the outlet boundary.Fortunately,such differences have a little effect on the 1T acoustic mode.Therefore,we ignore the difference on the acoustic characteristics of the 1L acoustic mode,and validate the present works only using the experimental data of the 1T acoustic mode.It can be found that the predicted half-power bandwidth and damping factor of the 1T acoustic mode agree well with those obtained by experiments,and the errors between them are less than 10%.Thus,the numerical method in this paper can capture acoustic and damping characteristics of the thrust chamber correctly.

    Fig.4 FFT analysis of pressure time series.

    5.2.Effects of imposed position on excited pressure oscillations

    In order to evaluate the acoustics characteristics and their damping of a considered thrust chamber,effects of the imposed position on excited pressure oscillations are discussed in this section.Here,the constant-volume bomb model is applied in a different geometry of an LRE combustor.The diameter and length of the cylindrical part are 35.25 mm and 7.1 mm,respectively,while the length of the contraction part is 25.63 mm.The diameter of the throat is 10.22 mm,and the total length of the thrust chamber is 89.73 mm.The theoretical acoustic eigenmodes are shown in Table 5,in which 2T and 1R represent the Second Tangential mode and the First Radial mode,respectively.

    For a numerical evaluation,the chamber is initially filled with quiescent air with a pressure of 0.1 MPa and a temperature of 298 K.Then,a constant-volume bomb with γ=19 and β = γ is imposed on a cylinder-shaped region with a radius of 3.6 mm and a length of 3 mm.The positions of the center of the imposed region for different cases are shown in Table 6.The imposed cylinder region is located on the axis of the chamber for Cases 1,2 and 3,but off-axis for Cases 4 and 5,respectively.The former and latter are called as the central and offcentered constant-volume bombs,respectively.

    Table 5 Acoustic eigenmodes by theoretical calculations.

    Table 6 Positions of center of imposed region.

    The pressure oscillations p′at an observation point in the imposed region are shown in Fig.5(a)for Cases 1,2,and 3,respectively,and their corresponding FFT analyses are shown in Fig.5(b).The amplitudes are normalized by the mean chamber pressure pmean.The peak frequencies and their amplitudes are shown in Table 7.Compared to theoretical eigenfrequencies,these peak frequencies are identified and the corresponding modes are marked in Fig.5(b).The pressure oscillations characterized by 1L and 1R acoustic modes are excited for all the three cases.The half-power bandwidths of the 1R acoustic mode are larger than those of the 1L acoustic mode,which means that the 1R pressure oscillations decay faster than the 1L pressure oscillations.The half-power bandwidths of the 1L and 1R acoustic modes for Case 3 are both larger than those for Cases 1 and 2.Because the constant-volume bomb for Case 3 is imposed on the nozzle convergent section which can attenuate pressure oscillations.The amplitudes of 1L and 1R are weaker with the imposed region away from the injector face.

    The effects of the radial position of the imposed region on the pressure oscillations are investigated through Cases 2,4 and 5.The pressure oscillations for these cases are shown in Fig.6(a)and their FFT analyses are shown in Fig.6(b).The shapes of pressure oscillations for Cases 4 and 5 are different from that for Case 2.Compared to Case 2,more peak frequen-cies are found for Cases 4 and 5,which are identified based on the theoretical eigenfrequencies,and the corresponding modes are marked in Fig.6(b).For the central constant-volume bomb,the major acoustic modes are 1L and 1R,while 1T for the off-centered constant-volume bomb.

    Table 7 Acoustic and damping characteristics of Cases 1–3.

    Fig.5 Pressure oscillations and their FFT analyses in Cases 1–3.

    Fig.6 Pressure oscillations and their FFT analyses in Cases 2,4 and 5.

    In order to reveal the propagating process of pressure waves,the temporal behaviors of the velocity vector in a transverse section of the chamber with a distance of 4 mm away from the injector face are given in Figs.7–9.The red color represents high pressure,while the blue and green colors represent low and middle pressures,respectively.The arrow displays the propagating direction of pressure waves.

    As shown in Fig.7(a)for Case 2,a high-pressure region is formed at the center of the transverse section,because the constant-volume bomb is imposed on the axis of the chamber.Pressure waves propagate to the chamber wall along the radial direction as shown in Fig.7(a)–(d).When the pressure waves reach the wall,they are re flected and propagate to the center as shown in Fig.7(e)and(f).Thus,the central high-pressure region is formed again.A typical radial acoustic mode is therefore established,which is marked as 1R in Fig.5(b).

    For Case 4,the center of the imposed region is located at the half of the chamber radius as shown in Fig.8(a).Fig.8(b)–(h)show the propagating process of the pressure waves due to a constant-volume bomb.They propagate in both the tangential and radial directions because of low pressure in the other regions of the chamber.The pressure waves propagate to the right side first and reach one side of the chamber wall.Then they propagate along the wall due to re flection.Tangential acoustic modes are triggered,shown as Fig.8(e)and(f).Thus,the 1T,2T,and 1R acoustic modes are observed for Case 4 as shown in Fig.6(b).

    For Case 5,a constant-volume bomb is imposed very close on the chamber wall as shown in Fig.9(a).Comparing Fig.9(b)to Fig.8(b),the pressure wave propagating along the wall in Case 5 is much stronger than that in Case 4.The tangential acoustic modes are therefore enhanced,and the radial acoustic mode is suppressed,which can also be observed in Fig.6(b).Fig.9(e)and(f)show the 1T and 2T modes occurring in this case.

    Fig.7 Temporal behavior of velocity vector for Case 2.

    Fig.8 Temporal behavior of velocity vector for Case 4.

    Fig.9 Temporal behavior of velocity vector for Case 5.

    A high-pressure region is formed where a numerical bomb is imposed.The pressure waves travel from the high-pressure region to the low-pressure region spherically.For the central constant-volume bomb,the high-pressure region is always axially symmetrical for the chamber.The pressure varies in the radial direction and keeps constant in the azimuthal direction.Thus,the pressure waves propagate only in the radial direction in the transverse section.Pure radial modes occur in the transverse section,and tangential modes do not appear.For an offcentered constant-volume bomb,the high-pressure region is not axially symmetrical in the transverse section.Thus,the pressure waves propagate in both the radial and tangential directions.The tangential modes are excited.The pressure waves also propagate in the longitudinal direction for thesetwo types of constant-volume bomb.Thus,the 1L mode occurs in all the cases.When the imposed constant-volume bomb moves from the axial center of the chamber to the sidewall,the tangential acoustic mode becomes stronger,while the radial and longitudinal acoustic modes become weaker.

    Table 8 Acoustic and damping characteristics of Cases 2,4 and 5.

    The detailed acoustic and damping characteristics are shown in Table 8.When the off-centered constant-volume bomb is imposed on the chamber,such as in Cases 4 and 5,the 1T and 2T modes can be stimulated.The 1T mode is the strongest acoustic mode compared to the other modes,while the 2T mode is very weak.Therefore,the most inspirable mode is the 1T mode with the off-centered constant-volume bomb.When the central constant-volume bomb is imposed,the dominant modes are 1L and 1R,and the amplitudes are higher than those in Cases 4 and 5.The half-power bandwidths of the transverse modes are larger than that of the 1L mode,so the 1L pressure oscillations keep for the longest time,and they are the most difficult to be attenuated.In the transverse modes,the half-power bandwidth of the 1R mode is larger than that of the 1T mode.

    The half-power bandwidth of the higher acoustic mode(2T)is larger than that of the low acoustic mode(1T).The damping capacities of tangential acoustic modes are almost the same in Cases 4 and 5.When the imposed region is closer to the chamber sidewall,pressure waves in the longitudinal direction are easier to be weakened by the nozzle,and thus the damping factor should be larger as shown in Table 8.

    From Case 1 to Case 5,three acoustic modes have been encountered in the LRE model:the longitudinal,tangential,and radial acoustic modes.The frequencies of the acoustic modes and the corresponding damping capacities numerically obtained by the present CFD method agree well with those by the experiments,but they are reasonably larger than those obtained by the theoretical method,because the effects of the converging section on the acoustic characteristics of the chamber are considered in the simulations.The 1L mode appears in all the cases,whose half-power bandwidth is the smallest.For these simulation cases,the 1L acoustic mode is the most difficult to be damped for the present chamber configuration,which is affected by its lowest eigenfrequency and the outlet boundary condition.The amplitudes of the 1R mode are larger for Cases 1–3,while they are smaller in the other cases.The tangential acoustic modes appear in Cases 4 and 5.The 1T mode has been proven to be the most violent and harmful mode,whose amplitude is found to be the largest in Cases 4 and 5.The half-power bandwidth of the 1R mode is wider than that of the 1T mode.In practice,pressure disturbances due to injection or violence combustion in the head region of the chamber usually are not located at the center of the chamber.This work explains why the 1T highfrequency combustion instability occurs frequently and is hard to be eliminated without passive devices or active control.

    6.Conclusions

    This paper proposes a new artificial constant-volume bomb model and numerically mathematizes it,considering the existing disturbance models with the shortcoming of a small overpressure coefficient,which cannot excite multi-frequency oscillations inside the chamber.Based on this model,the acoustic characteristics of a small-thrust LRE are investigated by the CFD method.A constant-volume bomb is imposed on the steady flow to activate pressure oscillations in the chamber.Peak frequencies are obtained through the FFT analyses of the temporally recorded pressure oscillations.Their corresponding acoustic modes are identified by comparing to the theoretical acoustic modes.The damping capacities of the identified acoustic modes are evaluated by the half-power bandwidth.The effects of the position of the constant-volume bomb on excited pressure oscillations are further investigated.

    It is found that the 1L mode can be triggered easily.The half-power bandwidth of the 1L mode is the smallest.The amplitude of the 1L mode is smaller but its half-power bandwidth is larger as the center of the imposed region is closer to the chamber sidewall.The 1R mode can be excited by a central constant-volume bomb,because the pressure disturbances keep constant in the azimuthal direction,such as in Cases 1,2,and 3.In the other cases,the 1R mode amplitudes are too small to be identified.The tangential modes can only be excited by an off-centered constant-volume bomb,in which the pressure disturbances are non-axisymmetric about the axis in the chamber.Moreover,the 1T mode is verified as the fiercest mode and therefore the most inspirable one.The damping factor of the 1T mode is independent of the radial position of the bombed region.The half-power bandwidth of the 2T mode is larger than that of the 1T mode,indicating a faster decay.

    This work establishes a CFD method coupled with a numerical constant-volume bomb model to evaluate acoustic characteristics,which can be applied to,but not limited to,small-thrust LREs.It somehow verifies and explains why the 1T mode in high-frequency combustion instabilities has been encountered frequently in engineering.

    Acknowledgements

    Financial support from the National Natural Science Foundation of China(Nos.51676111 and 11628206)is acknowledged.

    1.Harrje DT,Reardon FH.Liquid propellant rocket combustion instability.Washington,D.C.:NASA;1972.p.16–23.

    2.Yang V,Anderson WE.Liquid rocket engine combustion instability.Reston:AIAA;1995.p.3–12,380–5.

    3.Reardon FH,Crocco L,Harrje DT.Velocity effects in transverse mode liquid propellant rocket combustion instability.AIAA J 1964;2(9):1631–41.

    4.Crocco L,Sirignano WA.Behavior of supercritical nozzle under three dimensional oscillatory conditions.Paris:Advisory Group for Aerospace Research and Development;1969.Report No.:AGARDograph 117.

    5.Culick FEC.Some recent results for nonlinear acoustics in combustion chambers.AIAA J 1994;32(1):146–69.

    6.Zinn BT.A theoretical study of nonlinear combustion instability in liquid-propellant engines.AIAA J 1968;6(10):1966–72.

    7.Zinn BT,Lores ME.Application of Galerkin methods in solution of nonlinear axial combustion instability problems in liquid rockets.Combust Sci Technol 1972;4(1):269–78.

    8.Farshichi M,Mehrjou H,Salehi MM.Acoustic characteristics of a rocket combustor chamber:Radial baf fle effects.Appl Acoust 2009;70(8):1051–60.

    9.Hardi JS,Oschwald M,Dally B.Acoustic characterization of a rectangular rocket combustor with liquid oxygen and hydrogen propellants.Proc Inst Mech Eng Part G:J Aerosp Eng 2012;227(3):436–46.

    10.Sohn CH,Park IS,Kim SK,Kim HJ.Acoustic tuning of gasliquid scheme injectors for acoustic damping in a combustion chamber of a liquid rocket engine.J Sound Vib 2007;304(3–5):793–810.

    11.Ducriux S,Candel S.External flow modulation in computational fluid dynamics.AIAA J 2004;42(8):1550–8.

    12.Shimizu T,Shigeru T,Yoshida S,Hori D,Shingo M,Mizobuchi Y.Intense tangential pressure oscillations inside a cylindrical chamber.AIAA J 2011;49(10):2272–81.

    13.Shimizu T,Morii Y,Hori D.Numerical study on intense tangential oscillation in a simulated liquid rocket chamber.AIAA J 2014;52(8):1795–9.

    14.Shingo M,Junji S,Yasuhiro M.LES of high frequency combustion instability in a rocket combustor.Reston:AIAA;2013.Report No.:AIAA-2013-0564.

    15.Shingo M,Shinjo J,Ogawa S,Mizobuchi Y.LES of highfrequency combustion instability in a single element rocket combustor.Reston:AIAA;2012.Report No.:AIAA-2012-1271.

    16.Abdelkader F,Tom N,Francisco CL.Control of combustioninstabilities through various passive devices.Reston:AIAA;2005.Report No.:AIAA-2005-2832.

    17.Taro S,Dan H.Acoustic structure and damping estimation of a cylinder rocket chamber during oscillation.Reston:AIAA;2012.Report No.:AIAA-2012-4206.

    18.Kim YM,Chen ZJ,Chen CP,Ziebarth JP.Pressure-based method for combustion instability analysis.Int J Numer Meth Fluids 1994;19(11):981–95.

    19.Grenda JM,Venkateswaran S,Merkle CL.Multi-dimensional analysis of combustion instabilities in liquid rocket motors.Reston:AIAA;1992.Report No.:AIAA-1992-3764.

    20.Habiballah M,Lourme D.PHEDRE-Numerical model for combustion stability studies applied to the Ariane Viking engine.J Propul Power 1991;7(3):322–9.

    21.Zhuang F,Nie W,Zhao W,Liu W.Liquid rocket combustion instability analysismethodology–methodsand representative examples.Reston:AIAA;1998.Report No.:AIAA-1998-3690.

    22.Urbano A,Selle L,Staffelbach G,Cuenot B,Schmitt T,Ducruix S,et al.Exploration of combustion instability triggering using large eddy simulation of a multiple injector liquid rocket engine.Combust Flame 2016;169:129–40.

    23.Zhang HQ,Ga YJ,Wang B,W XL.Analysis of combustion instability via constant volume combustion in a LOX/RP-1 bipropellant liquid rocket engine.Sci China Ser E 2012;55(4):1066–77.

    24.Harvazinski ME,Anderson WE,Merkle CL.Analysis of selfexcited combustion instabilities using two-and three-dimensional simulations.J Propul Power 2013;29(2):396–409.

    25.Yuan L,Shen CB.Computational investigation on combustion instabilities in a rocket combustor. Acta Astronaut 2016;127:634–43.

    26.Anderson JD.Computational fluid dynamics.In:Wu SP,Liu ZM,translated.1st ed.Beijing:China Machine Press;2009.p.174–81[Chinese].

    27.Kim HJ,Lee KJ,Seo S,Han YM,Seol WS,Lee SY.Stability rating tests of KSR-III baf fled chamber using pulse gun.Reston:AIAA;2004.Report No.:AIAA-2004-3364.

    28.Buffum FG,Dehority GL,Slates RO,Price EW.Acoustic attenuation experiments on subscale,cold- flow rocket motors.AIAA J 1967;5(2):272–80.

    29.Oschwald M,Farago Z.Resonance frequencies and damping of a combustor acoustically coupled to an absorber.J Propul Power 2008;24(3):2524–33.

    30.Park IS,Sohn CH.Nonlinear acoustic damping induced by a halfwave resonator in an acoustic chamber.Aerosp Sci Technol 2010;14(6):442–50.

    31.Shimizu T,Hori D,Yoshida S.On acoustic damping of a cylindrical chamber in resonant modes.Fluid Dyn Res 2012;44(4):1–20.

    国产激情欧美一区二区| www国产在线视频色| 国产精品美女特级片免费视频播放器 | 亚洲成a人片在线一区二区| 亚洲欧美精品综合久久99| 免费在线观看日本一区| 午夜两性在线视频| 国产又色又爽无遮挡免费看| 日本a在线网址| 老司机午夜十八禁免费视频| 少妇的逼水好多| 精品国产超薄肉色丝袜足j| 午夜福利欧美成人| 国产黄色小视频在线观看| 无人区码免费观看不卡| 真人做人爱边吃奶动态| 不卡一级毛片| 长腿黑丝高跟| 俺也久久电影网| 欧洲精品卡2卡3卡4卡5卡区| 国产乱人伦免费视频| 18禁黄网站禁片免费观看直播| 国产黄a三级三级三级人| 99riav亚洲国产免费| 国产精品久久电影中文字幕| 亚洲av五月六月丁香网| 精品人妻1区二区| 亚洲无线在线观看| xxxwww97欧美| 成人永久免费在线观看视频| 2021天堂中文幕一二区在线观| 欧美日韩中文字幕国产精品一区二区三区| 热99re8久久精品国产| 99热精品在线国产| 欧美三级亚洲精品| 久久精品影院6| 午夜福利视频1000在线观看| 亚洲欧美日韩高清在线视频| ponron亚洲| 88av欧美| 精品一区二区三区视频在线 | 变态另类成人亚洲欧美熟女| 他把我摸到了高潮在线观看| 亚洲国产中文字幕在线视频| 老司机午夜福利在线观看视频| 香蕉久久夜色| 午夜福利在线观看免费完整高清在 | 中文资源天堂在线| 日韩欧美精品v在线| 成人特级av手机在线观看| 在线观看日韩欧美| 国产亚洲精品综合一区在线观看| 日韩免费av在线播放| 免费电影在线观看免费观看| 午夜免费激情av| 久久草成人影院| 日韩国内少妇激情av| 男女午夜视频在线观看| 国产高清有码在线观看视频| 老司机福利观看| 国产亚洲精品av在线| 99国产极品粉嫩在线观看| 综合色av麻豆| 精品一区二区三区四区五区乱码| 国产午夜精品论理片| 波多野结衣高清无吗| 九九热线精品视视频播放| 可以在线观看的亚洲视频| 看免费av毛片| 国产精品影院久久| 51午夜福利影视在线观看| www国产在线视频色| 国产精品乱码一区二三区的特点| 巨乳人妻的诱惑在线观看| 日日干狠狠操夜夜爽| 亚洲精品456在线播放app | 香蕉丝袜av| 夜夜夜夜夜久久久久| 久久精品国产综合久久久| 欧美又色又爽又黄视频| 国产黄色小视频在线观看| 精品午夜福利视频在线观看一区| 757午夜福利合集在线观看| 亚洲乱码一区二区免费版| 久久中文字幕人妻熟女| 婷婷丁香在线五月| 白带黄色成豆腐渣| 亚洲18禁久久av| 好男人在线观看高清免费视频| 国内精品久久久久精免费| 国产高清视频在线播放一区| 亚洲精品国产精品久久久不卡| 国产欧美日韩一区二区精品| 国产激情欧美一区二区| 在线视频色国产色| 国产精品久久视频播放| 可以在线观看毛片的网站| 国内少妇人妻偷人精品xxx网站 | 国产精品一区二区精品视频观看| 黄片小视频在线播放| 青草久久国产| 视频区欧美日本亚洲| 国内精品久久久久久久电影| 成人永久免费在线观看视频| 国产精品久久久久久精品电影| 亚洲精品久久国产高清桃花| 看片在线看免费视频| 日韩精品青青久久久久久| 中文字幕熟女人妻在线| 久久精品国产99精品国产亚洲性色| 亚洲自偷自拍图片 自拍| 美女午夜性视频免费| 欧美黑人巨大hd| 高清在线国产一区| 三级国产精品欧美在线观看 | 日本熟妇午夜| 午夜福利欧美成人| 国产亚洲精品久久久com| 夜夜躁狠狠躁天天躁| 日本免费一区二区三区高清不卡| 在线观看美女被高潮喷水网站 | 亚洲欧美一区二区三区黑人| 欧美高清成人免费视频www| 人人妻人人澡欧美一区二区| 国产成人精品久久二区二区91| 蜜桃久久精品国产亚洲av| 变态另类成人亚洲欧美熟女| 国产一区二区激情短视频| 免费在线观看成人毛片| 国产精品乱码一区二三区的特点| 88av欧美| 成人鲁丝片一二三区免费| 手机成人av网站| 91麻豆精品激情在线观看国产| 国产亚洲av嫩草精品影院| 国产1区2区3区精品| 午夜福利在线在线| 精品久久久久久久末码| 美女被艹到高潮喷水动态| 国产三级在线视频| 国产成人精品无人区| 丰满人妻熟妇乱又伦精品不卡| 国产极品精品免费视频能看的| 19禁男女啪啪无遮挡网站| x7x7x7水蜜桃| 免费av不卡在线播放| 深夜精品福利| 久久久国产成人精品二区| 老司机在亚洲福利影院| 国产一区二区激情短视频| 成人鲁丝片一二三区免费| 国产黄片美女视频| 亚洲av日韩精品久久久久久密| 国产日本99.免费观看| 婷婷亚洲欧美| 19禁男女啪啪无遮挡网站| 欧美黄色淫秽网站| 老司机在亚洲福利影院| 国产高清视频在线观看网站| 香蕉国产在线看| 中文字幕最新亚洲高清| 国产毛片a区久久久久| 人妻夜夜爽99麻豆av| 亚洲精品美女久久av网站| 性色avwww在线观看| 成人18禁在线播放| 国产毛片a区久久久久| 18禁黄网站禁片免费观看直播| 亚洲美女视频黄频| www.熟女人妻精品国产| 在线十欧美十亚洲十日本专区| 一a级毛片在线观看| 99久久无色码亚洲精品果冻| 国产精品 欧美亚洲| 一二三四在线观看免费中文在| 丰满人妻一区二区三区视频av | 久久久成人免费电影| 精品久久久久久久末码| 欧美色欧美亚洲另类二区| 欧美乱色亚洲激情| 亚洲美女黄片视频| 最近视频中文字幕2019在线8| 青草久久国产| 夜夜躁狠狠躁天天躁| 亚洲成人久久性| 国产黄a三级三级三级人| 亚洲欧美日韩高清专用| 成人三级做爰电影| 亚洲av中文字字幕乱码综合| av在线蜜桃| 亚洲av成人av| 桃红色精品国产亚洲av| 久久久久精品国产欧美久久久| 成人国产综合亚洲| 亚洲自拍偷在线| 国产一区二区三区在线臀色熟女| 黄色片一级片一级黄色片| 亚洲av中文字字幕乱码综合| 波多野结衣高清无吗| 99久久无色码亚洲精品果冻| 国产激情欧美一区二区| 国产蜜桃级精品一区二区三区| 国产成人aa在线观看| 日本熟妇午夜| 亚洲av电影在线进入| svipshipincom国产片| 国产乱人视频| 久久久久国内视频| 一个人看视频在线观看www免费 | 亚洲国产色片| 又紧又爽又黄一区二区| 亚洲狠狠婷婷综合久久图片| 日本一本二区三区精品| 美女午夜性视频免费| 一级作爱视频免费观看| 亚洲精品一卡2卡三卡4卡5卡| 国产精品日韩av在线免费观看| 成人亚洲精品av一区二区| 国产一级毛片七仙女欲春2| 久久中文字幕一级| 国产激情偷乱视频一区二区| 丰满的人妻完整版| xxxwww97欧美| 国产高清视频在线播放一区| 18禁裸乳无遮挡免费网站照片| 91在线观看av| 久久精品影院6| 99热这里只有精品一区 | 一边摸一边抽搐一进一小说| 欧美绝顶高潮抽搐喷水| 国产精品,欧美在线| 午夜免费激情av| 淫妇啪啪啪对白视频| 国产乱人伦免费视频| 搡老岳熟女国产| 99国产精品一区二区蜜桃av| 观看免费一级毛片| 色哟哟哟哟哟哟| 中文字幕熟女人妻在线| 精品日产1卡2卡| 国产精品久久久久久久电影 | 亚洲性夜色夜夜综合| 国产69精品久久久久777片 | 嫩草影院精品99| 老司机午夜十八禁免费视频| 亚洲av成人精品一区久久| 老汉色∧v一级毛片| 亚洲精品粉嫩美女一区| 夜夜夜夜夜久久久久| 免费无遮挡裸体视频| 黄片小视频在线播放| 国模一区二区三区四区视频 | 久久久成人免费电影| 麻豆av在线久日| www.精华液| 国语自产精品视频在线第100页| 国产精品一及| 黄片大片在线免费观看| 高潮久久久久久久久久久不卡| av视频在线观看入口| 欧美黄色淫秽网站| 又黄又爽又免费观看的视频| 黄色女人牲交| 这个男人来自地球电影免费观看| xxx96com| 此物有八面人人有两片| 精品久久久久久久毛片微露脸| 亚洲人成电影免费在线| 国产亚洲精品av在线| 亚洲18禁久久av| 国产精品久久电影中文字幕| 免费看a级黄色片| 后天国语完整版免费观看| 欧美激情久久久久久爽电影| 日本免费a在线| 丁香六月欧美| 亚洲熟妇熟女久久| 成人三级做爰电影| 三级男女做爰猛烈吃奶摸视频| 好男人电影高清在线观看| 成人一区二区视频在线观看| 搡老熟女国产l中国老女人| 亚洲 国产 在线| 一个人免费在线观看的高清视频| 精品午夜福利视频在线观看一区| 黄色视频,在线免费观看| 女警被强在线播放| 9191精品国产免费久久| 久久久国产欧美日韩av| 一个人免费在线观看电影 | 国产视频一区二区在线看| 九九久久精品国产亚洲av麻豆 | 一级毛片精品| 日韩欧美精品v在线| 国产精品久久久久久人妻精品电影| 观看免费一级毛片| 欧美+亚洲+日韩+国产| 18禁美女被吸乳视频| 亚洲国产精品久久男人天堂| 国产极品精品免费视频能看的| 欧美性猛交黑人性爽| 狂野欧美激情性xxxx| 亚洲午夜精品一区,二区,三区| 色吧在线观看| 亚洲欧美激情综合另类| 看片在线看免费视频| 少妇人妻一区二区三区视频| 国产成人精品久久二区二区91| 国产精品一区二区免费欧美| 99视频精品全部免费 在线 | 中文字幕高清在线视频| 日本与韩国留学比较| 国产欧美日韩一区二区精品| 亚洲成人免费电影在线观看| a在线观看视频网站| 国产乱人伦免费视频| 搡老妇女老女人老熟妇| 亚洲国产日韩欧美精品在线观看 | 国产成人aa在线观看| 国产午夜福利久久久久久| 成人无遮挡网站| 欧美激情久久久久久爽电影| 亚洲av熟女| 欧美乱色亚洲激情| 男女之事视频高清在线观看| 亚洲自偷自拍图片 自拍| or卡值多少钱| 一区二区三区国产精品乱码| 变态另类丝袜制服| 国产 一区 欧美 日韩| 国产免费av片在线观看野外av| 一卡2卡三卡四卡精品乱码亚洲| 免费在线观看影片大全网站| 午夜精品久久久久久毛片777| 全区人妻精品视频| 久久性视频一级片| 亚洲一区二区三区色噜噜| 午夜免费激情av| 小说图片视频综合网站| 国产黄色小视频在线观看| 欧美乱码精品一区二区三区| av在线天堂中文字幕| 人妻丰满熟妇av一区二区三区| 国产精品一及| 亚洲天堂国产精品一区在线| 国产精品免费一区二区三区在线| 丰满人妻熟妇乱又伦精品不卡| 一a级毛片在线观看| 91麻豆av在线| 午夜福利欧美成人| 色av中文字幕| 国产伦精品一区二区三区视频9 | 国产又黄又爽又无遮挡在线| 日本 av在线| 人妻久久中文字幕网| 久久人人精品亚洲av| 女警被强在线播放| 大型黄色视频在线免费观看| 亚洲欧美日韩高清专用| 日本免费一区二区三区高清不卡| 在线免费观看不下载黄p国产 | x7x7x7水蜜桃| 亚洲成人精品中文字幕电影| 国产97色在线日韩免费| 亚洲美女黄片视频| 狂野欧美激情性xxxx| 婷婷精品国产亚洲av在线| 亚洲精品国产精品久久久不卡| 91麻豆av在线| 97碰自拍视频| 国产亚洲精品久久久com| 日韩精品中文字幕看吧| 日本 av在线| 两人在一起打扑克的视频| 精品福利观看| 久久久久久大精品| 久久久久精品国产欧美久久久| 久久国产乱子伦精品免费另类| 给我免费播放毛片高清在线观看| 亚洲va日本ⅴa欧美va伊人久久| 国产亚洲欧美在线一区二区| 午夜福利成人在线免费观看| 人妻夜夜爽99麻豆av| 久久久久国内视频| 亚洲欧美日韩东京热| 女生性感内裤真人,穿戴方法视频| 成人18禁在线播放| 欧美av亚洲av综合av国产av| 色精品久久人妻99蜜桃| 欧美高清成人免费视频www| 1024香蕉在线观看| 欧美一级毛片孕妇| 嫩草影视91久久| 久久九九热精品免费| 婷婷亚洲欧美| 天天添夜夜摸| 真人做人爱边吃奶动态| 97超视频在线观看视频| 校园春色视频在线观看| 午夜免费成人在线视频| АⅤ资源中文在线天堂| 一级作爱视频免费观看| 小说图片视频综合网站| АⅤ资源中文在线天堂| 伦理电影免费视频| 一进一出抽搐gif免费好疼| 国产黄色小视频在线观看| 最新在线观看一区二区三区| 日韩国内少妇激情av| 日本 av在线| 在线观看免费午夜福利视频| 精品不卡国产一区二区三区| 色av中文字幕| 全区人妻精品视频| 舔av片在线| 精品一区二区三区四区五区乱码| 欧美激情久久久久久爽电影| 欧美黑人欧美精品刺激| 国产亚洲欧美98| 在线观看免费视频日本深夜| 日日干狠狠操夜夜爽| 少妇的逼水好多| 美女高潮喷水抽搐中文字幕| www.精华液| 亚洲欧美日韩无卡精品| 国产麻豆成人av免费视频| 亚洲国产高清在线一区二区三| 欧美日韩中文字幕国产精品一区二区三区| 制服人妻中文乱码| 久久精品夜夜夜夜夜久久蜜豆| 久久久国产成人精品二区| 日韩人妻高清精品专区| 亚洲色图 男人天堂 中文字幕| 日韩国内少妇激情av| 欧美性猛交╳xxx乱大交人| 夜夜躁狠狠躁天天躁| 久久天堂一区二区三区四区| 欧美中文日本在线观看视频| 日韩欧美国产一区二区入口| 成人av在线播放网站| 亚洲电影在线观看av| 日本 欧美在线| 久久久久国产精品人妻aⅴ院| 欧美性猛交╳xxx乱大交人| 一级毛片精品| 天堂动漫精品| 日韩欧美国产在线观看| 国产高清激情床上av| 中出人妻视频一区二区| av中文乱码字幕在线| 国产精品亚洲av一区麻豆| 日本 av在线| 国产精品一区二区免费欧美| 一个人免费在线观看的高清视频| 岛国在线观看网站| 又大又爽又粗| 18禁观看日本| 欧美3d第一页| 久久天堂一区二区三区四区| 一级毛片精品| 亚洲精品中文字幕一二三四区| 午夜激情欧美在线| 午夜视频精品福利| www国产在线视频色| 九九久久精品国产亚洲av麻豆 | 午夜免费成人在线视频| 法律面前人人平等表现在哪些方面| 在线免费观看的www视频| 国产av麻豆久久久久久久| 国产精品精品国产色婷婷| 欧美不卡视频在线免费观看| 精品欧美国产一区二区三| 人人妻,人人澡人人爽秒播| 91麻豆av在线| 99热只有精品国产| 午夜日韩欧美国产| 欧美zozozo另类| 亚洲中文字幕日韩| 狂野欧美白嫩少妇大欣赏| 丁香欧美五月| 国产精品一区二区免费欧美| 国产精品精品国产色婷婷| 免费无遮挡裸体视频| 午夜精品在线福利| 久久久国产欧美日韩av| 午夜久久久久精精品| 亚洲国产欧美人成| 欧美激情久久久久久爽电影| 久久久色成人| 亚洲成人精品中文字幕电影| 一本精品99久久精品77| 一个人免费在线观看的高清视频| 国内精品一区二区在线观看| 久久久色成人| 亚洲国产精品999在线| 一本精品99久久精品77| or卡值多少钱| 黄色日韩在线| 成人性生交大片免费视频hd| 欧美在线黄色| 一级毛片女人18水好多| 欧美日韩乱码在线| 久久精品91蜜桃| 免费在线观看成人毛片| 久久婷婷人人爽人人干人人爱| 国产精品电影一区二区三区| 欧美日韩乱码在线| 成年女人看的毛片在线观看| 欧美黄色淫秽网站| 在线免费观看的www视频| 9191精品国产免费久久| 一边摸一边抽搐一进一小说| 此物有八面人人有两片| 色尼玛亚洲综合影院| 亚洲精品粉嫩美女一区| 婷婷亚洲欧美| 亚洲人成电影免费在线| 色噜噜av男人的天堂激情| 中文字幕久久专区| 亚洲精品中文字幕一二三四区| 日本黄色视频三级网站网址| 国产69精品久久久久777片 | 国产亚洲精品av在线| 午夜福利在线观看吧| 好男人电影高清在线观看| 99热这里只有精品一区 | 久久久久精品国产欧美久久久| 成人特级av手机在线观看| 国产精品国产高清国产av| 99在线视频只有这里精品首页| 97人妻精品一区二区三区麻豆| 高清在线国产一区| 嫩草影院精品99| 麻豆成人av在线观看| 熟女人妻精品中文字幕| 亚洲美女视频黄频| 亚洲中文日韩欧美视频| 麻豆av在线久日| 国产精品久久久久久亚洲av鲁大| 高清在线国产一区| 国产精品自产拍在线观看55亚洲| 亚洲人成网站在线播放欧美日韩| 三级毛片av免费| 国产免费av片在线观看野外av| 中文字幕av在线有码专区| 91久久精品国产一区二区成人 | 美女大奶头视频| 高潮久久久久久久久久久不卡| 搡老妇女老女人老熟妇| 免费在线观看影片大全网站| 校园春色视频在线观看| 琪琪午夜伦伦电影理论片6080| 一个人看的www免费观看视频| 欧美日韩一级在线毛片| av欧美777| 亚洲av中文字字幕乱码综合| 久久天躁狠狠躁夜夜2o2o| 国产乱人视频| 亚洲精品乱码久久久v下载方式 | 在线免费观看不下载黄p国产 | 免费在线观看日本一区| 国产成人精品久久二区二区91| 老熟妇仑乱视频hdxx| 成人亚洲精品av一区二区| 国产一区二区激情短视频| 国产精品一区二区精品视频观看| 精品一区二区三区四区五区乱码| 国产精品一及| 亚洲自偷自拍图片 自拍| 欧美色视频一区免费| 国产真人三级小视频在线观看| 成人欧美大片| 日本 欧美在线| 精品人妻1区二区| 国产高清三级在线| 欧美日韩乱码在线| 99久久久亚洲精品蜜臀av| 亚洲男人的天堂狠狠| 成人特级av手机在线观看| 欧美日韩瑟瑟在线播放| 色精品久久人妻99蜜桃| 午夜福利在线观看免费完整高清在 | 国产欧美日韩一区二区三| 夜夜躁狠狠躁天天躁| 亚洲av电影在线进入| 亚洲国产精品合色在线| 欧美成狂野欧美在线观看| 亚洲真实伦在线观看| 欧美国产日韩亚洲一区| 麻豆国产97在线/欧美| 91在线观看av| 亚洲av成人av| 国产亚洲精品av在线| 国产成人精品无人区| 国产亚洲av高清不卡| 香蕉国产在线看| 日韩欧美国产一区二区入口| 午夜福利成人在线免费观看| 美女黄网站色视频| 99国产精品99久久久久| 精品国产乱码久久久久久男人| h日本视频在线播放| 国产精品亚洲av一区麻豆| 久久精品夜夜夜夜夜久久蜜豆| 精品一区二区三区av网在线观看| 综合色av麻豆| 成人高潮视频无遮挡免费网站| 日韩高清综合在线| 亚洲专区中文字幕在线| 国产69精品久久久久777片 | 免费在线观看亚洲国产| 男人舔女人下体高潮全视频| 国产高清videossex| 国产黄a三级三级三级人| 两性午夜刺激爽爽歪歪视频在线观看|