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

    Thermal state calculation of chamber in small thrust liquid rocket engine for steady state pulsed mode

    2019-02-27 08:59:06AlexeyGenndievichVOROBYEVSvtlnSergeevnVOROBYEVALihuiZHANGEvgeniyNikolevichBELIAEV
    CHINESE JOURNAL OF AERONAUTICS 2019年2期

    Alexey Genndievich VOROBYEV,Svtln Sergeevn VOROBYEVA,Lihui ZHANG,Evgeniy Nikolevich BELIAEV

    aDepartment of Rocket Engines,Moscow Aviation Institute,Moscow 125993,Russia

    bSchool of Astronautics,Beihang University,Beijing 100191,China

    Abstract This paper presents a method of thermal state calculation of combustion chamber in small thrust liquid rocket engine.The goal is to predict the thermal state of chamber wall by using basic parameters of engine:thrust level,propellants,chamber pressure,injection pattern, film cooling parameters,material of wall and their coating,etc.The difficulties in modeling the startup and shutdown processes of thrusters lie in the fact that there are the conjugated physical processes occurring at various parameters for non-design conditions.A mathematical model to predict the thermal state of the combustion chamber for different engine operation modes is developed.To simulate the startup and shutdown processes,a quasi-steady approach is applied by replacing the transient process with time-variant operating parameters of steady-state processes.The mathematical model is based on several principles and data commonly used for heat transfer modeling:geometry of flow part,gas dynamics of flow,thermodynamics of propellants and combustion spices,convective and radiation heat flows,conjugated heat transfer between hot gas and wall,and transient approach for calculation of thermal state of construction.Calculations of the thermal state of the combustion chamber in single-turn-on mode show good convergence with the experimental results.The results of pulsed modes indicate a large temperature gradient on the internal wall surface of the chamber between pulses and the thermal state of the wall strongly depends on the pulse duration and the interval.

    KEYWORDS Combustion chamber;Film cooling;Mathematical model;Nonstationary thermal mode;Small thrust liquid rocket engine;Steady pulse mode;Thermal state

    1.Introduction

    Thrusters are used as the executive devices of the control system for spacecraft's angular position.There is development mainly with the radiation cooling system,so the question of assessing the thermal state of such an engine becomes very important.Pulsed mode is the main operation mode of the thrusters when it is used as an executive element.In Ref.1,the requirements to the executive thrusters of control system,and the classification of operating modes,taking into account the thermal effect of engine wall heating,are considered.In many respects,the efficiency of the engine in the pulsed mode determines the quality of the operation and the lifetime of the entire spacecraft.2An big issue for the developer of thrusters is to ensure the permissible thermal state of the wall in each of its operation modes.

    An analysis of the thermal state of the combustion chamber,operating on a pulsed operating mode,is a rather difficult task to calculate.The greatest difficulty is the simulation of unsteady processes during engine startup and shutdown.3,4This is due to the complex thermochemical processes of conversion of fuel into combustion products,and the hydrodynamics of the flow under the non-calculated operating conditions of the injectors when the pressure in the combustion chamber is absent.Models for calculating processes in chamber at unsteady regimes are based on empirical relationships obtained during experimental engine processing.Mathematical models of the engine operation in the pulsed modes are compiled,and the thermal characteristics of the considered engines are shown.5,6The availability of reliable experimental data on intra chamber processes in the modes of filling and emptying the mixing head and chamber is the key to adequate modeling for unstable operation modes of engine.

    The difficulty in calculating the unsteady regime makes it necessary to use a quasi-stationary method when unsteady operating modes of the engine(startup and shutdown)are replaced by a set of stationary steady-state regimes with intermediate constant parameters of the working process in the combustion chamber.7,8

    Theoretical modeling of steady pulsed operating modes of the engine,when the influence of unsteady processes is small,is possible without involving a large number of experimental data,but with subsequent verification of the results obtained.The complexity of such calculations is as follows:

    (1)Use a large number of design relationships describing the physical processes occurring in the engine;

    (2)Use a large amount of information of a thermodynamic database for various engine operating modes,taking into account the distribution of the component ratio in the cross section of the chamber,including film cooling and near wall layers;

    (3)Develop a software program for calculating nonstationary thermal processes in the pulsed mode of operation of the engine.

    Despite the large number of works devoted to the prediction of the thermal state of the combustion chamber wall,the engineer sometimes needs a simple fast 2D tool which in the early stages of design would allow calculating and comparing engine parameters without hardware producing,without performing CFD calculations,because the exact geometry has not yet been determined.

    The prediction of the thermal state of the combustion chamber wall in pulsed operating modes allows analysis and optimization of engine design parameters during the development and testing phase.9,10

    The purpose of the work is,based on the refined mathematical model,to calculate the thermal state of the thrusters at the steady pulse-operating mode.

    2.Physical processes in combustion chamber of thrusters working at a single pulse

    The processes in the engine with a single pulse are represented in the following order:after the electric command signal is applied,the process of opening the valves starts with a certain delay,the components enter the chamber,and their ignition and pressure increase in the combustion chamber.Simultaneously,the beginning of the pressure up,the process of the outflow of the vapor-gas phase,and then that of the combustion products from the engine nozzle occur.When the gas moves along the wall of the combustion chamber,the boundary layers on the wall of the chamber and the nozzle increase.Consequently,when the pressure in the chamber changes,the interaction of the gas flow with the wall leads to the formation of an unsteady boundary layer,through which the heat transfer process from the gas to the wall takes place.After the pressure becomes steady state,the boundary layer stabilizes according to its dynamic parameters,because the core of the flow becomes stable and the gas parameters in the core of the flow can be assumed to be stationary.

    In view of the large thermal capacity and the thermal inertia of the combustion chamber wall,its temperature will vary.Because of the interaction of the boundary layer with the wall,the parameters of the boundary layer will change,as well as the heat flux to the wall.Thus,despite the stabilization of the flow in the core of the flow,stabilization of the boundary layer does not occur.The nonstationary process of the heat transfer process will take place until heat flux from products of combustion to wall will equal heat flux removed from the walls by heat radiation processes into space and insignificantly heat dissipation to the mixing head,valve assembly and the engine attachment points.Because heat radiation process will play an appreciable role in the heat balance only at high temperatures,the time to enter the stationary thermal state is comparatively large.Therefore,if the pulse duration is small,the heat transfer process will be nonstationary,even in the absence of stationary parameters of the working process in the combustion chamber.

    After switching off the fuel supply and the pressure decrease in the chamber,the process of engine cooling due to radiation begins.If the interval time between the pulses is large,the temperature of the chamber and the pressure become equal to the temperature and pressure of the external conditions.In this case,the next pulse starts with the parameters corresponding to the external parameters.In this case,the effect of the previous pulse on the subsequent pulse does not occur-the engine operates in the mode of single or disconnected pulses.If the interval time between pulses is small,the next pulse starts at a chamber temperature different from the ambient temperature,in which case the effect of the engine operation in the previous cycle on the initial conditions of the next pulse will be affected.This will be the mode of operation with temperature-related pulses.If during a single pulse,the wall temperature has time to stabilize for a sufficiently long time-this will be a continuous stationary mode of operation.

    In the mode of single pulse and pulsed operation,the boundary layer on the walls of the combustion chamber and the nozzle is nonstationary,because the temperature dependent parameters(density,viscosity,thermal conductivity of the gas)vary by wall temperature.This issue will in turn change the heat transfer to the wall.Thus,there is a relationship between the parameters of the boundary layer and the wall temperature-the conjugate problem,which considerably complicates the calculation of heat transfer.

    The design parameters of the combustion chamber also affect the heat 1transfer process.For example,the axial extension of the combustion chamber leads to a significant turbulence of the boundary layer,which in turn leads to an increase in the heat flux to the wall.The material of wall also determines the wall temperature state in its main sections.

    3.Description of mathematical model

    The mathematical model is a set of functional calculation modules and an interface for setting the initial data and outputting the calculation results.The mathematical model11,12includes the following calculation modules:

    (1)Module for calculating the main parameters of the thruster;

    (2)Module for constructing the internal and external profile of the wall of the engine;

    (3)Module for calculating the thermodynamic and thermophysical properties of combustion products;

    (4)Module for calculating internal and external heat fluxes acting on the wall;

    (5)Module for calculating dynamic boundary layer parameters;

    (6)Module for calculating the evaporation and mixing of film cooling;

    (7)Module for calculating the stationary thermal state of the thruster;

    (8)Module for calculating the nonstationary thermal state of the thruster for pulse operation modes.

    Convective heat transfer is determined by the parameters of the combustion products in the wall layer.In addition,the dissociation-recombination processes,the chemical reactions of combustion,evaporation and decomposition of liquid components in the boundary layer influence the convective heat transfer in the conditions of liquid rocket engines.It is difficult to take into account all these factors;therefore,when choosing the general calculated equations of convective heat transfer in engine conditions, only the effect of dissociation recombination is taken into account.

    In practice,the approximate formula of Ievlev method13is used to calculate the convective heat flux,taking into account the change in the parameters with respect to time,which has the form:

    wherexis axial coordinate,tis operation time,Bis coefficient taking into account the transition from viscosity μ0Gto viscosity μ1000Gat a temperature of 1000 K, τ(λ)is gas-dynamics function,Pcc(t) is combustion chamber pressure,S(Kmw(x),Tw(x,t))is function of thermophysical parameters of gas,which depends on the type of fuel,the ratio of components in the near wall layerKmw(x)and temperature of wallTw(x,t),is relative diameter,dthis throat diameter of combustion chamber,andPris Prandtl number.

    The model cannot predict whole temperature field of combustion chamber because they do not operate with Navier-Stokes equations,and there is no need to calculate the whole field.Instead of this,the model predicts the near wall parameters that directly affect heat transfer.The model uses solutions of the boundary layer equations including continuity equation,equations of momentum,energy,friction,chemical equilibrium,enthalpy and viscosity and others.The differential equations of the turbulent boundary layer are used to compile integral equations and then transform it to question(1).

    The convective flow is determined by the parameters of the combustion products near the wall,taking into account the mixing of the wall layer with the film cooling layer.When there is internal film cooling,it is necessary to use the parameters of the combustion products,determined by the ratio of the components near the wallKmw.

    The radiant heat flux from the combustion products to the inner wall surface and from the outer wall surface to the environment is calculated according to the Stefan-Boltzmann law the radiation power is directly proportional to the fourth power of the body temperature:

    where φfcis the coefficient the absorption of the radiant flux from the core of the combustion products through the near wall layer,εw.ef.is the effective coefficient of blackness of the wall material,εgis the coefficient of emissivity of gas,C0is the Stefan-Boltzmann constant,TkandTware the static core and wall temperatures(changing by time)on the gas side,respectively.

    If there is internal cooling of the wall by organizing the film cooling with a fuel or an oxidizer,it is necessary to calculate the length at which the film remains in the liquid state.The length of the evaporation section of the film can be determined from the balance:the heat transferred from the combustion products to the liquid film goes to its heating and evaporation14:

    where η is the coefficient that takes into account the process of partial tearing of the film from its wall surface(0.5-0.9),˙mfcis the mass flow of liquid film cooling,πDis the circumference of the combustion chamber,Tinitis the initial temperature of film liquid,TSis the boiling temperature of the liquid at a given pressure in the combustion chamber,clis the specific heat capacity of liquid at medium temperature,Tave=(Tinit+Ts)/2,Qs,is the heat of evaporation of component,and αkis the convective heat transfer coefficient,calculated without allowance for curtain cooling:

    whereTG0is the temperature of the gas in the wall layer andTwGis the wall temperature,which depend on the ratio of the components in current calculation section,if the film cooling is absent.The coefficient η depends on the hydrodynamics of the film flow and its interaction with the main flow η=f(Refc),Refc=ρ1W1δ1/μ1=m˙fc/(πDμ1),where ρlis the film liquid density,Wlis the average velocity of film cooling liquid,δlis the thickness of film cooling,and μlis the liquid viscosity.

    When modeling the thermal state of the combustion chamber wall,it is necessary to correctly set the conditions for external convection.When modeling the thermal state with vacuum conditions,external convection is excluded.When modeling under atmospheric conditions,it is necessary to analyze the external conditions acting on the heat exchange process.The external convective heat flux is calculated according to the following equation:qk.ext= αk.ext(T-T∞)where αk.extdetermines the intensity of heat exchange between the outer wall surface and the environment,Tis the temperature of the outer wall surface,T∞is the temperature of the environment.For natural convection,αk.ext=2÷25( J/s·m2·K),and for forced convection, αk.ext=25÷250( J·s·m2·K).The exact definition of αk.extis very difficult because in the test conditions and the working place there are a large number of factors affecting the heat transfer:ambient temperature,engine location,air velocity on the surface of the combustion chamber,etc.In this calculation, αk.ext=2÷ 25( J/s·m2·K),which is typical for turbulent natural convection.

    Except formula(1),heat convection in thrust chamber can be calculated by formula Bartz(5).The formula describes the relationship between heat convection intensity and fluid's physical features as well as flow behavior.Gas-side heat transfer coefficienthgis expressed by

    WhereDtis throat diameter,μ is the viscosity of gas,ns means nozzle stagnation condition,Cpis the specific heat capacity at constant pressure,P0is chamber pressure,C*is characteristic velocity,R=(Rtin+Rtout)/2 is the radius of curvature of nozzle contour at throat,AtandAare area of throat and calculated section respectively,and σ is the correction factor for property variations across the boundary layer and can be evaluated in terms of nozzle stagnation temperatureT0,local gasside chamber wall temperatureTw,specific heat γ,and local Mach numberMa:

    Quantity of heat convection is indicated by effective temperature:

    whereTgthe temperature of the gas in the wall layer.

    The calculation shows that Ievlev and Bartz equations for the same conditions have most likely results;the results of convective heat flux are shown in Fig.1.The difference is less than 10%for throat section and less than 1%for other section.In the future,we will use Ievlev approach for heat flux calculation.

    When simulating engine cooling after the tests are completed,it is necessary to take into account the possible purging procedure and,consequently,the presence of forced convection between the internal wall of the combustion chamber and the purge gas,in this case a convective heat transfer coefficient value is required.

    To solve the nonstationary heat transfer problem,the differential Fourier-Kirchhoff thermal conductivity equation is used in the case of a stationary flow and the absence of internal heat sources in cylindrical coordinates15:

    where λw(T),ρw(T),Cw(T)are thermal conductivity,density and heat capacity of the wall material of the combustion chamber as a function of temperature respectively,ris the geometrical size of the dimension(radius),zis the geometric size of the longitudinal dimension(step along the radius).

    To obtain a difference analogue,we will write the partial derivatives with allowance wherekis the time step,iis the number of the grid node along the axis of the shell,jis the grid node number perpendicular to the axis,Δzis the step along the axis,and Δriis the step along the radius in the cross section:

    The difference analogue of the equation will have the form:

    Fig.1 Ievlev and Bartz approaches in result of convective heat flux.

    For the solution,initial conditions are set up in the form of a temperature distribution on the surfaces and in the wall of the combustion chamber,and the boundary conditions in the form of the dependences of the change in the heat flux on the inner and outer surfaces of the combustion chamber and also at the ends on the side of the mixing head and the nozzle section.The methodology and calculation formulas used in the mathematical model are considered in more details in previous works,11,12 and the structural scheme is shown in Fig.2.

    As any kind of mathematical models,this model has limitation and assumptions.Foremost,as mentioned above,only near wall layer is calculated exactly,and the kernel layer seems to be stationary and uniform.

    4.Modeling pulsed mode of operation

    A typical experimental cyclogram of a single pulse of a thruster with a command signal duration of 0.1 s is shown in Fig.3.

    The operation of the engine can be divided into the following characteristic areas and points:

    (1)I-supply of voltage to the valves,the beginning of their opening,the flow of liquid components into the mixing head,the full opening of the valves and the final filling of the mixing head by the components;pointa-receipt of the first portions of fuel components into the combustion chamber volume.

    (2)II-entering of the first portions of the components in the combustion chamber,pre- flame processes and the increase of pressure to the level of 0.1 Pkdue to the fuel injection without the combustion processofthe propellants.pointb-ignition of propellants,engine thrust uprising;

    Fig.3 Cyclogram of operation of thruster with a command signal duration 0.1 s.

    (3)III-combustion of fuel and oxidizer with the formation of the main products of the reaction, filling the combustion products with a volume of combustion chamber and increasing the pressure up to a level of 0.9 from the nominal pressure Pk;pointc-engine output for full thrust pulse.

    (4)IV-engine operation in steady state;pointd-engine output level of nominal thrust.pointe-command to close the valves;pointf-actual closing of valves.

    (5)V-burn-out of components flowing from the trapped cavities of the mixing head-the end of the combustion process;pointg-pressure in the combustion chamber drops down to 0.9 Pk.

    Fig.2 Structural diagram of mathematical model.

    (6)VI-emptying the volume of the combustion chamber with the remaining propellant components and gases of the combustion products.

    Characteristic parameters of the pulsed operating mode are16:

    τST-electric signal time on thruster valves;

    τP-pause time between starts of single thruster;

    τPB-post burning time to 0.9Pk;

    τPS-post electric signal time to 0.1Pk;

    T-cycle time between starts of single thruster T= τST+ τP.

    The time parameters related to the time of pulse interlacing form the dimensionless parameters of the pulsed mode of the thruster:

    k=τST/T- filling factor;

    S=τP/T-duty cycle;

    f=1/T-frequency of working of thruster17,18

    According to the pulsed mode of the thruster,the problem of modeling the conjugate heat transfer is divided into two problems5,6,19:modeling heat transfer in the steady-state mode of the pulse and modeling heat transfer in unsteady mode in a quasi-stationary formulation.For the second case,equations for the steady-state mode are used with correction of changing thermodynamic parameters of combustion products' flow:pressure and temperatures.

    To solve the task,the process of growth and pressure drop is replaced by a set of dependencies:

    (1)The process between the points of supply of voltage to valve 0 and pointbis not modeled,the process of calculating the thermal state will start according to the calculated time delay;

    (2)The process of pressure growth is modeled by a linear relationship with known points (τ0.1;0,1Pk)-(τ0.9;0,9Pk),which are user-defined;

    (3)The steady-state process is modeled according to the stationary approach;

    (4)The process between the pointeand pointgis the conditional time of post-burning τPB,modeled according to the stationary approach,and τPBis set by the user;

    (5)The process of pressure drop to the end of the pulse is modeled by quadratic or exponential interpolation from known points( τST+ τPB;0,9Pk)- (τST+ τPS0,1Pk);

    (6)The process after τPSis not modeled.

    This approach allows avoiding the connection between the process of electromechanical opening-closing of the solenoid valve and the processes occurring in the mixing head and combustion chamber at the time of engine startup and shutdown.Values of the parameters τ0.1,τ0.9,τPBand τPSare set directly by user,using theoretical or previously obtained statistical characteristics of the operation of a particular thruster or analogs.

    The pressure drop in startup period when combustion chamber comes to the steady-state operating mode and pressure fluctuations during the engine operation can be neglected,assuming that short-term pressure pulsations do not exert a large influence on the thermal state of the chamber wallthe calculation is based on the steady-state average pressure in chamber.

    Such an assumption restricts the applicability of the mathematical model,which is used only for steady-state pulse engine operation.For unsteady pulse modes of operation,modes with cut pulse,where the proportion of nonstationary processes is large enough,regression models based on statistical testing of the results of fire tests of thrusters should be used to determine the thermal state of the chamber wall and the engine's energy parameters.

    5.Object of investigation and initial data for calculation

    The initial data for the calculation were the parameters of the low-thrust rocket engine developed in the Moscow Aviation Institute,nitrogen tetroxide and asymmetric dimethyl hydrazine,200 N of thrust.The engines of this dimension,on these fuel components,are most often used as executive controls for the reactive spacecraft control system.20-22

    The experimental engine DMT MAI-200-1(Fig.4)has a removable mixing head,which is fixed to the body of the combustion chamber by a flange connection.The sealing of the flange-mixing head was carried out using a copper ring.The head of the engine DMT MAI 200-1 has one two-component centrifugal jet and a slotted film cooling along the circumference,which provides protection of the walls of the combustion chamber(Fig.5).

    The main characteristics of the engine MAI-200-1 are as follows:

    (1)Thrust-200 N;

    (2)Working time-50 s;

    (3)Components: nitrogen tetroxide+unsymmetrical dimethyl hydrazine(AT+UDMH);

    (4)Specific pulse without film cooling-292 s;

    (5)Coefficient of completeness of combustion-0.92;

    (6)Working pressure in the combustion chamber-0.94 MPa;

    (7)Degree of expansion of the nozzle by pressure-1060;

    (8)Mixing ratio of components-2.6;

    (9)Total mass flow rate-0.066 kg/s;

    (10)Film cooling component-fuel;

    Fig.4 Combustion chamber with mixing head and geometry profile of LREST MAI-200-1.

    Fig.5 Design of single-jet mixing head of engine MAI-200-1.

    (11)Relative mass flow for film cooling(by sum mass flow)-20%.

    6.Results of experiments in single-turn-on mode

    For the simulation of the pulsed mode,the initial data for the calculation and the time parameters are set(τ0.1,τ0.9,τPBandτPS)which are taken based on the results of the fire tests of the MAI-200-1 in the single-turn-on mode of different duration(from 2 to 50 s).

    During the tests,thruster with a shortened nozzle is used.The test scheme without a nozzle part allows checking the efficiency of the combustion chamber itself for long-term work without the involvement of expensive vacuum equipment.As shown by the calculation of the thermal state of the chamber with the nozzle,the presence of the nozzle part does not significantly affect the thermal characteristics of the combustion chamber.

    The stand(Fig.6)is equipped with a special working area,which includes a frame 2 for securing the engine in a vertical position,a protective casing 7 and a gas-dynamic cooled pipe 1 for flue gas discharge to the stand neutralization system.The combustion chamber of the engine 4 with the installed head and thermocouples is fixed to the frame of the working area with the aid of three pins passing through the head flange.

    The main sensors used to analyze the working process in the engine are:

    (1)Two sensors for measuring the pressure in the chamber;

    (2)Two turbine sensors for measuring mass flow of fuel;

    (3)Two turbine sensors for measuring mass flow of oxidizer;

    (4)Thermocouples mounted on the outer surface of the chamber(4 section with 2 reverse thermocouples for each section);

    (5)Two pressure sensors at the engine inlet(in front of the valves);

    (6)Two thermocouples for measuring the temperature of the components at the entrance to the engine;

    (7)Two current measurement of electromagnet valve.

    Comparison between the results of experiments and calculations of the thermal state of the combustion chamber wall in the single-start-up mode is shown in Figs.7 and 8.

    Calculations of the thermal state of the combustion chamber show good convergence with the experimental results of LREST MAI-200-1 at the level of 5-10%.The discrepancies are explained by the following reasons:

    (1)The real boundary conditions differ from the functions adopted in the calculation of the boundary conditions qk(x,t)and qrad(x,t),and the real boundary conditions are more complex and then adopted in mathematical model;

    (2)At the initial moment of time,the exact definition of the boundary conditions is problematic because of the unsteady combustion and near wall processes in the chamber.

    Fig.6 Schematic diagram of fire stand.

    Fig.7 Comparison between experimental(interpolation points)and calculated(lines)temperature of external surface of combustion chamber by time(current time in seconds displayed in gray square).

    Fig.8 Comparison between experimental(interpolation points)and calculated(lines)temperature of external surface of combustion chamber by time in different sections.

    (3)Thex=0.12 m section was near exit section of nozzle and near gas-dynamic cooled pipe(see Fig.6,position 1).When engine works,the nitrogen passes through protective casing(see Fig.6,position 7)into cooled pipe near exit section of nozzle right near thermocouplex=0.12 m that can be the reason of decrease of real temperature wall.

    The calculation results and test data show that the most heat-stressed place is the wall near critical section.The highest rate of growth of temperature is in the same place.It agrees with the theoretical issue.

    One can draw a conclusion about the adequacy of the numerical solution of the nonstationary heat conduction problem for a single launch.This enables us to use the mathematical model for further calculations of the pulsed mode in the quasi-stationary formulation of the problem.

    7.Results of calculation of a thermal state in a pulsed mode of operation

    Figs.9 and 10 show the results of mathematical modeling of the thermal state of the combustion chamber wall of the engine in a pulsed steady-state mode,for the case of 20%of the fuel film cooling and α=0.85 for different pulse duration and pause time at a duty ratioS=0.5.In the calculation,the operating conditions of the engine with a cut nozzle are simulated under the conditions of an atmospheric fire stand.It is assumed that a vertical installation of the engine,and mutual radiation between the combustion chamber and the elements of the stand,are excluded.The ambient temperature is 300 K.The values of the characteristic time are taken as follows: τ0.1=0.01 s; τ0.9=0.02 s; τPB=0.01 s; τPS=0.03 s.Fig.9 also shows the change of pressure with calculation time.

    In the case of pulsed operation,a cyclic thermal action is observed on the inner surface of the chamber.When the engine is switched off,the heat is redistributed from the more heated layers(inner surface of the wall)to the less heated layers(outer surface)and from the more heated zones(subsonic part of the nozzle)to the less heated(cylindrical part of the combustion chamber).The layers that are closer to the inner surface of the chamber work with a large temperature amplitude,which decreases with subsequent heating of the wall.Thus,coatings that are used in the production of thruster must have sufficient strength and adhesion properties,allowing them to work in the conditions of cyclic exposure to temperatures and combustion products.

    Figs.11 and 12 show the simulation results of the operation of the engine on a firing atmospheric stand by a composite cyclogram:a single pulse of 10 s,5 pulses of 2 s with 2 s interval,5 pulses of 1 s with 1 s interval,and purging the engine by nitrogen.Fig.12 shows the change of the temperature of the outer surface of the combustion chamber wall in different sections for this operating mode.

    The results of calculation(Figs.11 and 12)clearly show the influence of the cycle time between single inclusions on the thermal state of the combustion chamber wall for different sections.The critical sectional area during the transition to the pulse mode stops being heated,and as the interval time increases,the wall temperature decreases.The zone of the cylindrical part of the combustion chamber continues to be heated during the transition to the pulse mode of the engine,There is no experimental verification of this result,but we will plan in future.

    Fig.9 Variation of wall temperature of combustion chamber in critical section with time.(Pulse duration 0.05 s,interval time 0.05 s,20%of fuel film cooling,α=0.85).

    Fig.10 Variation of wall temperature of combustion chamber in critical section with time.(Pulse duration 1 s,interval time 1 s,20%of fuel film cooling,α=0.85).

    Fig.11 Variation of wall temperature of combustion chamber at critical section with time when modeling composite engine operation cycle.

    Fig.12 Variation of temperature of external surface of combustion chamber wall in different sections with time for modeling composite operation cycle of engine.(x=0.04 m-combustion chamber,x=0.09 m-critical section,x=0.12 m-out of short nozzle).

    In fact,it is possible to use model to predict thermal state of any thruster.For doing this,the engineer should set up(or calculate)geometry and thickness of wall,injector pattern, film cooling parameters,materials,propellants,chamber pressure,mixture ratio and others.If any single-turn-on mode exists,it allows the adjusted model to have more accuracy for thermal prediction.The new results using this model show that it can be possible to get dependencies of thermal state of combustion chamber from any general parameter of engine and working mode.There are more than 50 parameters,excluding dataset of materials,chemical elements,propellants and others.The model is a hopeful tool for engineer of space propulsion system and can be extended.For example,based on this model,the calculation of regenerative cooling of liquid rocket engine(including transient heat transfer during filling up the cooling jacket of coolant)is done.For verifying the model,the experimental data are now collected and systemized.

    8.Conclusions

    (1)A mathematical model that allows predicting the thermal field of a combustion chamber wall in various engine operating modes is developed.The physical processes taking place in the engine for startup and shutdown mode are considered.

    (2)A theoretical study of the thermal state of a model liquid rocket engine operating in a pulse steady-state mode has been investigated.The cases of engine operation in pulsed modes with different cycle time between single inclusions are considered,as well as a mode simulating fire tests at an atmospheric stand with varying time of single pulses.

    (3)A large temperature amplitude is obtained on the internal surface of the combustion chamber when calculating the pulsed mode.On the first pulse,the amplitude of the temperature is approximately 400 K,and then the amplitude decreases noticeably and reaches values of up to 100 K,which indicates that the thermal state of the wall in a pulsed steady-state mode strongly depends on the pulse duration and the interval between the operations.

    (4)The mathematical model can be used for modeling thermal state of thruster with various propellants,mixing scheme,geometry of combustion chamber and nozzle,and material of wall,including wall coating.

    欧美av亚洲av综合av国产av| 成人永久免费在线观看视频| 高清黄色对白视频在线免费看| 国产91精品成人一区二区三区| 国产乱人伦免费视频| 狂野欧美激情性xxxx| 在线观看免费视频网站a站| 亚洲七黄色美女视频| 国产单亲对白刺激| 日韩有码中文字幕| 巨乳人妻的诱惑在线观看| 亚洲精品一卡2卡三卡4卡5卡| 亚洲中文日韩欧美视频| 色尼玛亚洲综合影院| 日本一区二区免费在线视频| 一边摸一边抽搐一进一小说| 长腿黑丝高跟| 久久国产精品男人的天堂亚洲| 日韩欧美国产一区二区入口| av天堂在线播放| 欧美日本亚洲视频在线播放| 久久国产精品影院| 亚洲一区中文字幕在线| 国产成人精品久久二区二区免费| 亚洲欧洲精品一区二区精品久久久| 99久久综合精品五月天人人| 日本一区二区免费在线视频| 免费搜索国产男女视频| 欧美激情极品国产一区二区三区| 国产精品,欧美在线| 波多野结衣av一区二区av| 嫩草影视91久久| 黄色 视频免费看| 亚洲熟女毛片儿| 欧美精品啪啪一区二区三区| 啦啦啦韩国在线观看视频| 最新在线观看一区二区三区| 精品熟女少妇八av免费久了| 丝袜人妻中文字幕| 日本免费a在线| 露出奶头的视频| 国产欧美日韩精品亚洲av| 一进一出好大好爽视频| 亚洲国产看品久久| 亚洲精品一区av在线观看| 91成人精品电影| 天天躁夜夜躁狠狠躁躁| 久久欧美精品欧美久久欧美| 亚洲av熟女| 午夜a级毛片| 久久国产精品男人的天堂亚洲| 国产xxxxx性猛交| 在线观看免费视频网站a站| 一进一出好大好爽视频| 久久中文看片网| 国内毛片毛片毛片毛片毛片| 欧美日本亚洲视频在线播放| 丝袜在线中文字幕| 国产视频一区二区在线看| 麻豆成人av在线观看| 成年版毛片免费区| 色老头精品视频在线观看| 成人国产综合亚洲| 欧美不卡视频在线免费观看 | 中文字幕人妻丝袜一区二区| 99在线人妻在线中文字幕| АⅤ资源中文在线天堂| 90打野战视频偷拍视频| 日韩欧美免费精品| 久久久久久久久久久久大奶| 老鸭窝网址在线观看| 久久久久国内视频| 麻豆av在线久日| 国产精华一区二区三区| 国产成人av教育| 精品午夜福利视频在线观看一区| 亚洲成国产人片在线观看| 亚洲国产中文字幕在线视频| 成人18禁高潮啪啪吃奶动态图| 久久国产亚洲av麻豆专区| 在线观看舔阴道视频| 国产熟女xx| 欧美不卡视频在线免费观看 | 亚洲精品中文字幕一二三四区| 俄罗斯特黄特色一大片| 91av网站免费观看| 精品免费久久久久久久清纯| 免费人成视频x8x8入口观看| 夜夜爽天天搞| 中亚洲国语对白在线视频| 久久性视频一级片| 一级作爱视频免费观看| 很黄的视频免费| 久久草成人影院| www.精华液| 中文字幕最新亚洲高清| 国产精品亚洲av一区麻豆| 亚洲avbb在线观看| 级片在线观看| 国产成人av激情在线播放| 一级毛片精品| 亚洲中文日韩欧美视频| 精品久久久精品久久久| 亚洲国产日韩欧美精品在线观看 | 日本免费a在线| 欧美成人性av电影在线观看| 一级,二级,三级黄色视频| 欧美成狂野欧美在线观看| 99香蕉大伊视频| 国产午夜福利久久久久久| 色尼玛亚洲综合影院| 欧洲精品卡2卡3卡4卡5卡区| 最近最新中文字幕大全电影3 | 一边摸一边抽搐一进一小说| 女人高潮潮喷娇喘18禁视频| 久久人人精品亚洲av| 久久国产精品男人的天堂亚洲| 久久精品国产亚洲av香蕉五月| 18禁国产床啪视频网站| 国产精品免费一区二区三区在线| 黄色视频不卡| 国内毛片毛片毛片毛片毛片| 亚洲一码二码三码区别大吗| 午夜免费激情av| 欧美成狂野欧美在线观看| 在线观看日韩欧美| 最新在线观看一区二区三区| www.精华液| 欧美一级a爱片免费观看看 | 99在线人妻在线中文字幕| 亚洲性夜色夜夜综合| 亚洲aⅴ乱码一区二区在线播放 | 99国产极品粉嫩在线观看| 美女免费视频网站| 午夜a级毛片| 99在线人妻在线中文字幕| 亚洲九九香蕉| 中亚洲国语对白在线视频| 欧美色视频一区免费| 国产精品亚洲av一区麻豆| 丰满人妻熟妇乱又伦精品不卡| 国产成人免费无遮挡视频| 亚洲va日本ⅴa欧美va伊人久久| 国产亚洲av高清不卡| 国产精品久久久久久亚洲av鲁大| 国产一区在线观看成人免费| 欧美午夜高清在线| 女生性感内裤真人,穿戴方法视频| 制服丝袜大香蕉在线| 最近最新免费中文字幕在线| 精品熟女少妇八av免费久了| 精品久久蜜臀av无| 黑人欧美特级aaaaaa片| 男男h啪啪无遮挡| 亚洲最大成人中文| 成人特级黄色片久久久久久久| 亚洲国产高清在线一区二区三 | 国产激情欧美一区二区| 久久精品亚洲精品国产色婷小说| 精品不卡国产一区二区三区| 午夜两性在线视频| 亚洲成人精品中文字幕电影| 日韩欧美三级三区| 两个人视频免费观看高清| 久久久久九九精品影院| x7x7x7水蜜桃| 正在播放国产对白刺激| 好看av亚洲va欧美ⅴa在| 欧美不卡视频在线免费观看 | 亚洲av电影不卡..在线观看| 91老司机精品| 一本大道久久a久久精品| 精品不卡国产一区二区三区| 18禁观看日本| 亚洲成人精品中文字幕电影| 亚洲第一青青草原| 嫩草影院精品99| 久久天躁狠狠躁夜夜2o2o| 18禁黄网站禁片午夜丰满| 午夜成年电影在线免费观看| 91九色精品人成在线观看| 夜夜看夜夜爽夜夜摸| 久久人人精品亚洲av| 亚洲男人天堂网一区| 51午夜福利影视在线观看| 在线观看免费日韩欧美大片| 免费在线观看亚洲国产| 精品久久久久久成人av| 脱女人内裤的视频| 欧美国产精品va在线观看不卡| 亚洲中文字幕日韩| 午夜福利18| 男女床上黄色一级片免费看| 欧美成人性av电影在线观看| 久久久国产成人免费| 成人国语在线视频| 久久影院123| 国产成+人综合+亚洲专区| 亚洲精品在线美女| 国产精品久久久av美女十八| 亚洲av成人av| 国产成人精品无人区| 女警被强在线播放| av天堂久久9| 女人爽到高潮嗷嗷叫在线视频| 成人特级黄色片久久久久久久| 久久国产乱子伦精品免费另类| 一二三四社区在线视频社区8| 亚洲aⅴ乱码一区二区在线播放 | 黄色 视频免费看| 亚洲国产毛片av蜜桃av| 国产成年人精品一区二区| 亚洲一码二码三码区别大吗| 人人妻人人澡人人看| 国内久久婷婷六月综合欲色啪| 中文字幕久久专区| 人人妻人人澡欧美一区二区 | 免费在线观看黄色视频的| 欧美一级毛片孕妇| 看免费av毛片| 亚洲精品在线美女| 欧美激情 高清一区二区三区| 性色av乱码一区二区三区2| 好男人在线观看高清免费视频 | 亚洲国产欧美一区二区综合| 国产麻豆成人av免费视频| 美女扒开内裤让男人捅视频| 性少妇av在线| 日本在线视频免费播放| 国产亚洲精品av在线| 好看av亚洲va欧美ⅴa在| 欧美乱码精品一区二区三区| 欧美日韩瑟瑟在线播放| 国产熟女xx| 大型黄色视频在线免费观看| 一区福利在线观看| 99国产综合亚洲精品| 69av精品久久久久久| av片东京热男人的天堂| 制服人妻中文乱码| 亚洲专区字幕在线| 色哟哟哟哟哟哟| 丰满的人妻完整版| 亚洲av熟女| 国产精品秋霞免费鲁丝片| 国内久久婷婷六月综合欲色啪| 亚洲自偷自拍图片 自拍| 又紧又爽又黄一区二区| 黄色 视频免费看| 一a级毛片在线观看| 女性生殖器流出的白浆| 免费看十八禁软件| 自拍欧美九色日韩亚洲蝌蚪91| 在线国产一区二区在线| 亚洲 欧美 日韩 在线 免费| 丁香欧美五月| 美女高潮到喷水免费观看| 这个男人来自地球电影免费观看| 欧美丝袜亚洲另类 | 亚洲视频免费观看视频| 视频区欧美日本亚洲| 男人操女人黄网站| 黄色丝袜av网址大全| 99国产精品一区二区三区| 久久人人97超碰香蕉20202| 国产欧美日韩精品亚洲av| 一区二区三区精品91| 成人18禁高潮啪啪吃奶动态图| 亚洲精品一卡2卡三卡4卡5卡| 久久狼人影院| 免费搜索国产男女视频| 亚洲第一欧美日韩一区二区三区| 中文字幕最新亚洲高清| 亚洲午夜精品一区,二区,三区| 午夜影院日韩av| 中文字幕人妻丝袜一区二区| 香蕉国产在线看| 免费在线观看视频国产中文字幕亚洲| 免费无遮挡裸体视频| 国产野战对白在线观看| 亚洲av日韩精品久久久久久密| 欧美色视频一区免费| 久久久久久久久免费视频了| 中文亚洲av片在线观看爽| 国产一卡二卡三卡精品| av超薄肉色丝袜交足视频| 亚洲天堂国产精品一区在线| 久久国产精品影院| 午夜免费观看网址| 禁无遮挡网站| 亚洲情色 制服丝袜| 国产精品免费视频内射| 国产精品久久久久久精品电影 | 日韩av在线大香蕉| 欧美另类亚洲清纯唯美| 亚洲色图av天堂| 国产又色又爽无遮挡免费看| 狠狠狠狠99中文字幕| 夜夜夜夜夜久久久久| 国产成+人综合+亚洲专区| av在线播放免费不卡| www.www免费av| 精品国内亚洲2022精品成人| 99久久久亚洲精品蜜臀av| 在线天堂中文资源库| bbb黄色大片| 国产熟女xx| 天天躁狠狠躁夜夜躁狠狠躁| 亚洲九九香蕉| 精品午夜福利视频在线观看一区| 欧美国产精品va在线观看不卡| 美女免费视频网站| 免费人成视频x8x8入口观看| 亚洲国产精品久久男人天堂| 欧美精品啪啪一区二区三区| 国产又色又爽无遮挡免费看| 亚洲精品久久国产高清桃花| 亚洲情色 制服丝袜| www日本在线高清视频| 999久久久精品免费观看国产| 国产成+人综合+亚洲专区| 亚洲国产精品久久男人天堂| 久久这里只有精品19| 国产成人精品无人区| 一二三四社区在线视频社区8| 美女大奶头视频| 露出奶头的视频| 十分钟在线观看高清视频www| 色播在线永久视频| 日韩精品免费视频一区二区三区| 狠狠狠狠99中文字幕| 中出人妻视频一区二区| 老司机在亚洲福利影院| av视频免费观看在线观看| 午夜福利在线观看吧| 欧美国产日韩亚洲一区| 国产亚洲av高清不卡| av视频免费观看在线观看| 男女之事视频高清在线观看| 美女扒开内裤让男人捅视频| 欧美绝顶高潮抽搐喷水| 十分钟在线观看高清视频www| 国产精品一区二区在线不卡| 不卡一级毛片| 欧美日本中文国产一区发布| 亚洲黑人精品在线| 久久久久国产一级毛片高清牌| 久久青草综合色| 亚洲专区中文字幕在线| 精品一区二区三区av网在线观看| 91在线观看av| 午夜精品国产一区二区电影| 一a级毛片在线观看| 久久中文字幕一级| 亚洲av电影不卡..在线观看| 日日爽夜夜爽网站| 精品卡一卡二卡四卡免费| 午夜影院日韩av| 国产黄a三级三级三级人| 国产精品爽爽va在线观看网站 | 国产成人欧美在线观看| xxx96com| 久久久久久人人人人人| 午夜福利,免费看| 黄频高清免费视频| 美女扒开内裤让男人捅视频| 国产亚洲精品久久久久久毛片| 午夜福利18| 长腿黑丝高跟| 国产欧美日韩综合在线一区二区| 正在播放国产对白刺激| 欧美日韩黄片免| 亚洲,欧美精品.| 久久午夜亚洲精品久久| 国产熟女午夜一区二区三区| 狂野欧美激情性xxxx| 日日摸夜夜添夜夜添小说| xxx96com| 91字幕亚洲| 色综合亚洲欧美另类图片| 18禁黄网站禁片午夜丰满| 韩国av一区二区三区四区| ponron亚洲| 亚洲第一青青草原| 欧美色视频一区免费| 9色porny在线观看| 黑人巨大精品欧美一区二区蜜桃| 日本精品一区二区三区蜜桃| bbb黄色大片| 久久热在线av| 国产精品98久久久久久宅男小说| 亚洲,欧美精品.| 麻豆av在线久日| 久久 成人 亚洲| 国产精品亚洲一级av第二区| 亚洲欧美日韩无卡精品| 国产成人精品在线电影| 97超级碰碰碰精品色视频在线观看| 伦理电影免费视频| 老司机在亚洲福利影院| 亚洲,欧美精品.| 午夜福利在线观看吧| 50天的宝宝边吃奶边哭怎么回事| 少妇 在线观看| 涩涩av久久男人的天堂| 国产精品久久视频播放| 级片在线观看| 国产99白浆流出| 久久天堂一区二区三区四区| 欧美乱妇无乱码| 国产av一区在线观看免费| 韩国av一区二区三区四区| 啦啦啦免费观看视频1| 日韩欧美免费精品| 亚洲av电影不卡..在线观看| 国产麻豆69| 亚洲第一电影网av| 午夜精品国产一区二区电影| 在线免费观看的www视频| 亚洲国产精品成人综合色| 多毛熟女@视频| 麻豆成人av在线观看| 亚洲成人精品中文字幕电影| 日本 欧美在线| 国产亚洲欧美精品永久| 91九色精品人成在线观看| 久久婷婷人人爽人人干人人爱 | 国产av一区在线观看免费| 国产精品一区二区精品视频观看| 美女高潮喷水抽搐中文字幕| 波多野结衣一区麻豆| 91在线观看av| 亚洲自拍偷在线| 香蕉丝袜av| 久久午夜综合久久蜜桃| 亚洲电影在线观看av| 大陆偷拍与自拍| 亚洲成国产人片在线观看| 亚洲精品在线美女| 在线天堂中文资源库| 99在线视频只有这里精品首页| 国产精品亚洲美女久久久| 午夜免费成人在线视频| 满18在线观看网站| 国产精品久久电影中文字幕| 两性午夜刺激爽爽歪歪视频在线观看 | 国产一区二区三区视频了| 国产又爽黄色视频| 亚洲一码二码三码区别大吗| 波多野结衣一区麻豆| 国产成人免费无遮挡视频| 一级a爱片免费观看的视频| 大型av网站在线播放| 午夜免费鲁丝| 少妇粗大呻吟视频| 国产精品,欧美在线| x7x7x7水蜜桃| 国产亚洲精品av在线| 禁无遮挡网站| 自拍欧美九色日韩亚洲蝌蚪91| 亚洲av电影在线进入| 国产区一区二久久| 成人永久免费在线观看视频| 91老司机精品| 中文字幕精品免费在线观看视频| 又紧又爽又黄一区二区| 国产av精品麻豆| 在线观看日韩欧美| 国语自产精品视频在线第100页| 丰满的人妻完整版| 国产区一区二久久| 少妇熟女aⅴ在线视频| 女生性感内裤真人,穿戴方法视频| 久久亚洲真实| 色尼玛亚洲综合影院| 免费看a级黄色片| 国产精品av久久久久免费| 高潮久久久久久久久久久不卡| 12—13女人毛片做爰片一| 一区二区日韩欧美中文字幕| 日韩有码中文字幕| 欧美乱码精品一区二区三区| 亚洲精品美女久久久久99蜜臀| 制服丝袜大香蕉在线| 精品国产一区二区久久| 丝袜在线中文字幕| 18美女黄网站色大片免费观看| 人妻丰满熟妇av一区二区三区| 精品电影一区二区在线| 天天添夜夜摸| 人妻丰满熟妇av一区二区三区| 国产精品影院久久| 一区二区三区精品91| 人妻丰满熟妇av一区二区三区| 国产91精品成人一区二区三区| 女人被狂操c到高潮| 一二三四在线观看免费中文在| 嫩草影院精品99| 级片在线观看| 美女午夜性视频免费| 亚洲五月天丁香| 亚洲va日本ⅴa欧美va伊人久久| 国产亚洲欧美在线一区二区| 狂野欧美激情性xxxx| av天堂在线播放| a级毛片在线看网站| 中文字幕精品免费在线观看视频| 久久天堂一区二区三区四区| 亚洲欧美精品综合久久99| 又黄又爽又免费观看的视频| 久久婷婷人人爽人人干人人爱 | 午夜福利视频1000在线观看 | 757午夜福利合集在线观看| av视频免费观看在线观看| 99久久精品国产亚洲精品| 国语自产精品视频在线第100页| 免费看美女性在线毛片视频| 中文字幕最新亚洲高清| 美女免费视频网站| 两性夫妻黄色片| 免费在线观看视频国产中文字幕亚洲| 日韩欧美一区二区三区在线观看| 首页视频小说图片口味搜索| 夜夜看夜夜爽夜夜摸| 黄频高清免费视频| 人妻丰满熟妇av一区二区三区| 一进一出好大好爽视频| 一级作爱视频免费观看| 日本黄色视频三级网站网址| 天天添夜夜摸| 日日夜夜操网爽| 99国产精品免费福利视频| 亚洲欧美激情综合另类| 亚洲 国产 在线| 亚洲av熟女| 女性被躁到高潮视频| 日本欧美视频一区| 成人国语在线视频| 欧美成人午夜精品| 一边摸一边抽搐一进一小说| 窝窝影院91人妻| 99在线人妻在线中文字幕| 亚洲欧美激情综合另类| 不卡一级毛片| 欧美成人性av电影在线观看| 国产激情欧美一区二区| 一a级毛片在线观看| 成年女人毛片免费观看观看9| 国产视频一区二区在线看| 国产午夜精品久久久久久| av福利片在线| 亚洲va日本ⅴa欧美va伊人久久| 免费高清在线观看日韩| www.熟女人妻精品国产| 久久人人97超碰香蕉20202| 久久久水蜜桃国产精品网| 欧美丝袜亚洲另类 | 免费观看精品视频网站| 国产高清视频在线播放一区| tocl精华| 夜夜看夜夜爽夜夜摸| 免费看美女性在线毛片视频| 女人爽到高潮嗷嗷叫在线视频| 久久中文看片网| 国产aⅴ精品一区二区三区波| 人人妻人人澡欧美一区二区 | 91大片在线观看| 久久久久亚洲av毛片大全| 99国产精品一区二区蜜桃av| 人妻丰满熟妇av一区二区三区| 久久精品亚洲熟妇少妇任你| 免费看美女性在线毛片视频| 搞女人的毛片| 久久久久久国产a免费观看| 国产99久久九九免费精品| 最好的美女福利视频网| 三级毛片av免费| 欧美不卡视频在线免费观看 | 最近最新中文字幕大全免费视频| 亚洲专区字幕在线| 99精品欧美一区二区三区四区| 高潮久久久久久久久久久不卡| 男女做爰动态图高潮gif福利片 | 不卡av一区二区三区| 欧美日本中文国产一区发布| 国产成人av激情在线播放| 亚洲片人在线观看| 多毛熟女@视频| 999久久久精品免费观看国产| 欧美日韩一级在线毛片| 天天躁夜夜躁狠狠躁躁| 99riav亚洲国产免费| 动漫黄色视频在线观看| 免费看美女性在线毛片视频| 国产精品亚洲美女久久久| 国产三级黄色录像| 一级毛片精品| 免费av毛片视频| 女人精品久久久久毛片| 免费在线观看影片大全网站| 操出白浆在线播放| 国产乱人伦免费视频| 琪琪午夜伦伦电影理论片6080| 亚洲男人的天堂狠狠| 91麻豆精品激情在线观看国产| 狂野欧美激情性xxxx| 啦啦啦免费观看视频1| 69精品国产乱码久久久| 大型黄色视频在线免费观看| 后天国语完整版免费观看| 久久久久久人人人人人| 久久精品成人免费网站| 精品国产超薄肉色丝袜足j| 90打野战视频偷拍视频| 曰老女人黄片|