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

    Mechanism analysis and evaluation of thermal effects on the operating point drift of servo valves

    2022-05-10 10:07:44JianKANGZhaohuiYUANJingchaoLI

    Jian KANG,Zhao-hui YUAN,Jing-chao LI

    School of Automation,Northwestern Polytechnical University,Xi'an 710072,China

    Abstract:Operating point drift over large temperature spans can significantly degrade the performance of servo valves.The direction and magnitude of the deviation of the operating point are uncertain.To analyze and evaluate the mechanism of this complex system with a multi-level structure and multi-variables,it is necessary to construct a theoretical model with a clear physical concept to describe it.However,since the physical processes contain complex variations of structural parameters and flow properties,there is a problem of simplifying approximations in deriving analytical mathematical relations.The advantages of multi-physics field numerical analysis can compensate for this shortcoming of analytical formulations.Based on this,we constructed a whole-valve transfer function model to realize the mechanism analysis and evaluate the operating point drift when a thermal effect acts on a servo valve.The results show that the asymmetric fit relationship between the armature-nozzle assemblies is an important reason for the drift of the operating point caused by the thermal effect.Differences in structural parameters and fluid medium characteristics at different temperatures lead to nonlinear changes in the operating point.When the deviation angle reaches ±1°,an increase in temperature will cause the absolute value of the tangent slope of the displacement deviation of the spool to decrease from 1.44×10-5 m/℃to 1.25×10-6 m/℃.The influence of the deviation angle is reflected in the change in the absolute value of the tangent slope of the pressure deviation from 1.14×103 Pa/℃to 110 Pa/℃.

    Key words:Servo valve;Operating point drift;Mathematical model;Numerical analysis;Thermal effect

    1 Introduction

    Pressure control valves are widely used in aircraft braking systems because of their high precision,high sensitivity,and fast response.However,the large temperature variations during aircraft takeoff and landing pose a significant challenge to the performance of pressure servo valves.Therefore,the mechanism of thermal effects on the operating point drift of pressure servo valves has become a current issue for analysis.

    Servo valves are composed of a multi-stage structure,and the thermal effect on the output is the result of the combined effect of several subsystem characteristics.The study of the characteristics of each subsystem forms the basis of the system analysis.The torque motor,as the first stage of the servo valve,was first modeled and analyzed at the theoretical level.To obtain the ideal mathematical model more conveniently,Merritt (1967) considered the magnetic circuit as a 2D plane.Based on this study,an improved torque motor model was derived by considering the reluctance and leakage of permanent magnets(Urata,2000;Zhang QF et al.,2020).Subsequent theoretical modeling attempted to make the model more consistent with the actual physical process by improving the model of the equivalent magnetic circuit (Yan et al.,2017).Equating a 3D problem to a 2D model inevitably ignores some details.Since then,researchers have tried to explain the working process using 3D numerical simulations.The results showed that this method is more accurate than the equivalent model,but the internal pressure characteristics and flow state are still poorly understood.To observe the flow state inside the pilot stage,Mchenya et al.(2011) observed the jet shape,cavitation,and vorticity using a high-speed camera and computer-aided image measurement techniques.However,the experimental method is limited by the precise mating dimensions,which makes it difficult to achieve reproduction of the actual state and accurate data acquisition.To gain insight into the velocity and vortex distribution inside the pilot stage under different flapper positions and inlet pressure conditions,the use of computational fluid dynamics (CFD) techniques is a more suitable approach(Jacob et al.,2011).The results of numerical simulations revealed that cavitation occurs where vortexes are generated (Li et al.,2013).A CFD method was used to try to verify the suppression of this phenomenon by using new flapper and nozzle shapes (Aung and Li,2014;Yang et al.,2019).The high-speed fluid entering the pilot stage impacts the flapper,which is subjected to the main flow force.The vortex that appears after the fluid interferes with the solid flapper makes it subject to a secondary flow force(Aung et al.,2014).To grasp the principles and characteristics of the flow dynamics and the source of forced vibrations of the armature assembly,Li et al.(2018) derived and validated mathematical models of the flow dynamics and forced vibrations.For complex flow states that are difficult to calculate theoretically,numerical simulations can be an effective alternative to experiments (Pan et al.,2011).The numerical simulation technique can effectively compensate for the shortcomings of traditional theoretical analysis and better analyze the system characteristics of actual 3D problems.

    For the study of system-level output characteristics,the construction of a mathematical model with clear physical concepts is a more intuitive and desirable approach.Merritt (1968) proposed a third-order mathematical model of a servo valve,which is widely used for its high accuracy.Kim and Tsao (2000) developed a fifth-order system model of a servo valve based on consideration of some nonlinear effects.Based on differential geometry theory,Mu and Li (2011)considered the higher-order nonlinear terms in the system and applied input and output linearization methods to the mathematical modeling of a servo valve.Liu and Jiang(2014)derived a seventh-order model of the dynamic response of an electro-hydraulic servo valve from the nonlinear equations considering two degrees of freedom of the armature-flapper assembly.For high-precision servo valves,it is difficult to accurately represent the dynamic performance of the system by relying only on theoretical modeling methods.Numerical analysis has been used to achieve more accurate predictions of complex flow states in the flow field,which can be used to refine the theoretical models (Pountney et al.,1989).Numerical analysis also enables finite element analysis of critical and precision components,such as feedback rod assemblies and spring tubes,to obtain the required simulation parameters for the actuator system (Somashekhar et al.,2007).Yu et al.(2014) used CFD techniques to evaluate the flow dynamics and valve frequency response characteristics of the new rotary direct drive servovalve spool.Conventional 2D simulations and theoretical calculations always estimate the receiving pressure lower than the experimental data.Yan et al.(2019) proposed two hypotheses about the pilot stage flow distribution based on which a more complete theoretical model can be constructed.It has been demonstrated that the flow field is too complex,and thus an accurate mathematical model cannot be established.Numerical simulations make complex hydrodynamic problems simple,and the results describe the actual physical processes more accurately than do the theoretical models(Li,2016).

    The various analytical methods described above are based on performance analysis and characterization of the system in a normal temperature environment.The evaluation of the performance of servo valves in environments with large temperature spans by introducing temperature variables into the system analysis model introduces new problems to be faced.The direct effects of high-and low-temperature environmental conditions on electro-hydraulic servo valves are reflected in structural clearances,fluid characteristics,spool leakage,and material hardness(Zhang et al.,2010).These changes affect the flow coefficient of the throttle orifice and the change of zero control pressure(Zhao et al.,2015).The effect of hydraulic fluid temperature variations on the control characteristics of servo valves can be analyzed based on AMESim software (Chen and Wang,2018).To find out the causes of servo valve current fluctuation due to temperature drift,Li et al.(2017)investigated the servo valve current fluctuation due to temperature drift by simulation.Temperature also affects the spool operating force(Zhang Y et al.,2020).

    Therefore,to analyze and evaluate the mechanisms of thermal effects on the operating point of a pressure servo valve,the following work was carried out.Firstly,under the basic assumptions,a fundamental analysis of the forces on the spool was carried out through a problem description,and the various components affecting the spool motion were identified.A theoretical model of the components affecting the spool state was then established based on the working principle of the servo valve,and an expression for the moment containing environmental variables and structural parameters was derived.Expressions with clear physical concepts can explain the mechanism of thermal effects on output.Based on this,the complex physical process was simulated numerically based on a 3D physical model of actual dimensions,taking into account the variation of structural parameters and fluid properties.Finally,a mathematical model of the entire valve was constructed using the mathematical model fitted by numerical analysis.The analysis and evaluation of the mechanism of the thermal effect at the operating point of the pressure servo valve were realized by analyzing the trends of the factors affecting the spool motion.

    2 Description of problem

    2.1 Composition and working principle of a servo valve

    An electro-hydraulic servo valve is usually composed of an electrical-mechanical converter,pilot stage,spring tube,and spool.The double-nozzle flapper pressure servo valve in this study consists of a torque motor,nozzle flapper pilot stage,and spool (Fig.1).The torque motor converts an input electrical signal into mechanical quantities,such as torque and rotation angle,which are then converted by the spring tube into a linear displacement or rotation angle to drive the motion of the pilot stage.

    Fig.1 Structure diagram of a nozzle flapper pressure servo valve

    The pilot stage of the hydraulic control element has an important influence on the steady-state control accuracy and dynamic response characteristics of a servo valve.The double-nozzle flapper pilot stage uses a frictionless variable throttle orifice amplifier to regulate the control cavity pressure on each side of the spool through the variable throttle orifice overflow area.The spool described here uses integrated pressure feedback to achieve pressure output through the combined effect of differential pressure in the feedback,control cavity,and other forces.

    2.2 Basic assumptions

    Our purpose was to make the analysis more general and more convenient for subsequent theoretical derivation and numerical analysis.The basic assumptions were as follows:

    (1) The shape of the spool and sleeve is ideal cylindrical.

    (2) The radial clearance between the spool and sleeve,and the lateral force on the spool,are not considered.

    (3) The inhomogeneity of the structure is not considered.

    (4)The effect of pressure pulsation generated by the fluid is not considered.

    2.3 Analysis of the force balance relationship of the spool and its influencing factors

    As the last stage of the servo valve,the movement of the spool directly affects the output pressure of the whole valve.The spool is subjected to the control cavity differential pressure force,feedback cavity differential pressure force,flow force,inertia force,and damping force.The differential pressure in the control cavity directly causes the balance of forces on the spool to change and affects the size of the valve opening.

    The kinetic equation of the spool can be expressed as

    where

    p

    and

    p

    are the pressures of the control cavities on the left and right sides,respectively,

    A

    is the area of action of the control pressure,

    p

    and

    p

    are thepressures of the feedback cavities on the left and right sides,respectively,and

    A

    is the area of action of the feedback pressure.

    F

    is the inertia force,

    F

    is the damping force,and

    F

    is the flow force.The force analysis of the valve spool is shown in Fig.2.

    Fig.2 Schematic diagram of the force analysis of valve spool. Fcl and Fcr are the driving forces of the control cavities on the left and right sides,respectively

    The drift of the operating point under the action of the thermal effect is due to the destruction of the spool force balance relationship.The factors that affect the control pressure at both ends of the spool include mainly the torque motors,spring tube,and pilot stage components.Therefore,the mechanism analysis and evaluation of the thermal effect on the drift of the servo valve operating point need to be analyzed comprehensively from the system level.

    Changes in temperature are most likely to affect the structural parameters of the torque motor and spring tube assembly (Fig.3).The electromagnetic field of the torque motor is used to obtain the electromagnetic torque coefficient,and the structural field of the spring tube assembly is used to obtain the stiffness coefficient.For the subsystems involving the flow field,due to the strong interference between the fluid and the structure,the flow coefficient and the flow force need to be analyzed.

    Fig.3 Schematic diagram of the effect of temperature loading on the pressure in the control cavity

    In this study,based on the working principle of the servo valve,intensive parameter modeling of each main subsystem was carried out.By constructing a theoretical model with a clear physical concept to describe the relationship between structural parameters,environmental variables,input parameters,and performance,we attempt to explain the mechanism of the thermal effect on the drift of the servo valve operating point.Multiphysics numerical analysis was used to numerically fit the input-output relationship of each main subsystem.On this basis,a full-valve transfer function model was constructed to realize the evaluation of the operating point drift during the thermal effect of the servo valve.

    3 Theoretical model of thermal effects acting on the components of the servo valve

    The main purpose of this section is to obtain clear mathematical expressions of physical concepts at each stage,based on consideration of the influence of thermal effects on structural parameters and medium characteristics.Analysis of the mechanism of the thermal effect on the drift of the servo valve operating point is based on the physical process.

    3.1 Electromagnetic torque of the torque motor

    The function of the torque motor is to convert the input control signal into a rotation angle or displacement proportional to the input signal.In this study,we adopted a permanent magnet moving iron type torque motor composed of a permanent magnet,magnetic conductor,energized coil,and armature (Fig.4).The armature is located at the top of the spring tube and can swing around its top in four working air gaps.When the control signal is input,the control magnetic flux generated by the control coil passes through the four air gaps in different directions,so that the sum of the magnetic flux of the air gaps is no longer equal.The force balance of the armature is broken,and the armature generates a torque proportional to the control signal.A schematic diagram and equivalent circuit diagram of the magnetic circuit with intrinsic flux are shown in Fig.5.

    Fig.4 Schematic diagram of the torque motor structure and magnetic circuit

    Fig.5 Schematic diagram (a) and equivalent circuit diagram(b)of the magnetic circuit with intrinsic flux(the parameters are explained in the text)

    Part of the intrinsic flux generated by the permanent magnet returns to the magnet after passing through the magnetic conductor and the working air gap,while the other part returns to the magnet through the air.The effective fluxes through the magnetic conductor and the working air gap are denoted by

    φ

    and

    φ

    ,respectively,while the leakage fluxes are indicated by

    φ

    and

    φ

    .To complete the flux calculation using the mesh current method,let the fluxes inside the armature and conductor circuit be

    φ

    and

    φ

    ,respectively.

    Based on the basic method of circuit calculation,the following expressions for the relationship between magnetic flux,magnetoresistance,and magnetic potential can be obtained:

    where

    R

    (

    i

    =1,2,3,4) is the reluctance of the four working air gaps,

    R

    is the internal resistance of the permanent magnet,

    R

    is the reluctance of the magnetic leakage,

    R

    is the reluctance of the magnet on the magnetic circuit,

    R

    is the reluctance of the armature,and

    M

    is the magnetomotive force of the permanent magnet.

    Ignoring the reluctance change of the permeable magnet,the air gap reluctance and magnetomotive force will directly affect the magnetic flux of the air gap.We assumed that the flux strength of the magnetic conductor and armature would not saturate,and neglected the hysteresis of the magnetic circuit.

    When the temperature rises,the thickness of the air gaps 1 and 2 at the neutral position can be expressed as

    where

    H

    is the height of the center of the armature,which is related to the thermal expansion of the spring tube.

    H

    is the thickness of the armature,and

    H

    is the height of the magnetic conductor.

    α

    ,

    α

    ,and

    α

    are the coefficients of thermal expansion of the spring tube,armature,and magnetic conductor,respectively.

    T

    is the temperature.

    Similarly,the thickness of air gaps 3 and 4 can be expressed as

    where

    H

    is the height of the permanent magnet,

    H

    is the thickness of the magnetic conductor,and

    α

    is the coefficient of thermal expansion of the permanent magnet.

    The magnetic flux through the air gap is influenced not only by the magnetic resistance,but also by the magnetic potential of the permanent magnets.The magnetic momentum provided by the permanent magnet to the external magnetic circuit after removing the internal reluctance to consume magnetic momentum can be expressed as

    where

    l

    is the height of the permanent magnet,

    α

    is the temperature coefficient of remanence,

    μ

    is the magnetic permeability,

    R

    is the irreversible loss rate,and

    B

    is the remanence at room temperature.The relationship between the magnetic field strength and flux density of a permanent magnet is represented by the

    B

    -

    H

    characteristic curve.When the

    B

    -

    H

    curve of the permanent magnet is approximated as linear,the magnetic permeability

    μ

    of the permanent magnet can be regarded as a constant.

    where

    F

    is the force acting on the armature’s cantilever,

    φ

    is the magnetic flux passing through the working air gap,

    μ

    is the permeability of air,and

    A

    is the effective area of the air gap.

    The expression for the electromagnetic moment can be expressed as

    where

    a

    is half the length of the armature,and

    F

    (

    i

    =1,2,3,4) are the electromagnetic forces generated by the four working air gaps,respectively.

    Therefore,the thermal effect will act on the structural parameters of the magnetic conductor and the performance of the permanent magnet.Variations in the properties of permanent magnets affect the intrinsic magnetic flux in the magnetic circuit.At the same time,under the premise of changing the structural parameters of the permeable magnet,the deflection movement of the armature under the action of the control signal will cause a difference in the magnetic flux of the four air gaps.Based on the combined effect of these two main factors,the angle of rotation of the armature and the electromagnetic torque output in different temperature environments will differ from those at room temperature.

    3.2 Bending moment of the spring tube

    The spring tube is a thin-walled round tubular structure made from a highly elastic beryllium copper alloy with a high elastic modulus.The spring tube is assembled on the valve body between the torque motor and the pilot stage (Fig.6).A flapper located in the center of the armature fits into the spring tube using interference assembly.The armature is activated by the electromagnetic driving force to counteract the torsion moment of the spring tube and the flow force of the pilot stage to deflect the flapper.

    The stiffness of the spring tube can be expressed as

    where

    D

    and

    d

    are the outer diameter and inner diameter of the spring tube,respectively,

    l

    is the length of the spring tube,and

    E

    is the elastic modulus.Thus,The elastic modulus

    E

    decreases with the increasing temperature,andwhere

    η

    is the temperature coefficient of the modulus of elasticity.Metal physics research shows that the ratio of the temperature coefficient of elastic modulus

    η

    to the linear expansion coefficient

    α

    is a constant value of 0.04.The expression for the change in stiffness with respect to the change in temperature can be written as

    where

    K

    is the stiffness of the spring tube at room temperature.

    Therefore,the thermal effect will act on the elastic modulus and structural parameters of the spring tube,thereby affecting its stiffness.Under different thermal loads,the same action will result in different bending moments and deflection angles.

    3.3 Hydrodynamic torque of fluid

    The thermal effect on the pilot stage and the valve body is directly manifested in the thermal expansion and fluid characteristics of the structure.Previous studies have found that the influence of temperature on fluids is reflected mainly in two properties:viscosity and density.The viscosity and density characteristics of fluids can be modeled by the following equations:

    where

    υ

    and

    ρ

    are the kinematic viscosity and fluid density at normal temperature,respectively.

    α

    and

    α

    are the temperature coefficients,and

    β

    and

    β

    are the pressure coefficients.

    p

    is the pressure of the oil.Temperature has a more significant effect on fluid properties than does pressure.Therefore,we considered only temperature in this study and ignored the influence of pressure.

    The flow coefficient of the throttle orifice is related to the Reynolds number.The Reynolds number characterizes the ratio of the inertial and viscous forces of a fluid,which determines whether the flow pattern of the fluid is laminar or turbulent.The Reynolds number of a fixed throttle orifice can be expressed as

    where

    υ

    is the viscosity,

    r

    is the radius of the fixed throttle orifice,and Δ

    p

    is the pressure difference between the two sides of the orifice.

    The Reynolds number of the nozzle can be expressed as

    where

    D

    is the nozzle diameter,

    x

    is the initial distance between the nozzle and the flapper,

    x

    is the displacement of flapper,and Δ

    p

    is the pressure difference across the nozzle.

    Thermal effects will directly affect the Reynolds number through density,viscosity,and structural parameters.Based on previous studies,we conclude that the flow coefficient is a function of the Reynolds number.

    The nozzle ejects fluid that acts on the flapper flow force in two main ways.On the one hand,flow force is generated by the stagnation pressure in the nozzle projection area.This is due to the high-speed impact of the fluid on the flapper surface,converting the kinetic energy of the fluid into static pressure energy.On the other hand,flow force is generated by the change in fluid momentum on the remaining flapper surface outside the nozzle projection area.

    where

    F

    and

    F

    are the flow forces on each side of the flapper,

    p

    and

    p

    are the pressures of the control cavity on each side of the spool,respectively,

    C

    is the flow coefficient of the nozzle,and

    p

    is the return pressure.

    Therefore,the thermal effect will cause changes in fluid properties and structural parameters,which will directly affect the flow force.The force balance relationship of the spool and the working point of the whole valve is determined mainly by the flapper-nozzle gap flow resistance.The position of the flapper depends on the combined effects of the electromagnetic moment,the bending moment of the spring tube,and the hydrodynamic moment.These critical moments to which the armature-flapper assembly is subjected are all functions of temperature.Therefore,such complex system problems need to be evaluated with the help of mathematical models of the whole valve to assess their performance.However,there is the problem of simplifying the approximation in deriving the analytical mathematical relationships,and a multi-physics field numerical analysis is needed to analyze the forces on the armature-flapper assembly.

    4 Simulation results and discussion

    Within complex systems containing electromagnetic,fluid,and solid interactions,it is difficult to study the thermal effects of servo valves by theoretical analysis alone,especially the interference between fluids and solids present in the pilot stage and valve body.To obtain a more accurate model of the local characteristics of each stage and to further evaluate the effect of thermal effects on the operating point drift of the servo valve,numerical simulations of the electromagnetic field of the torque motor,the structural field of the spring tube,and the flow field of the pilot stage are described in this section.

    4.1 Numerical simulation of electromagnetic torque on the armature

    In the theoretical analysis section,the output torque equation was derived by analyzing the magnetic circuit of a simplified 2D torque motor.Due to the non-axisymmetric structure of the torque motor and the complexity of the magnetic circuit,if the structure is considered as a 2D problem then many details related to the structure are ignored.In addition,there are approximate estimates for the calculation of the air gap reluctance,component reluctance,and leakage flux,which eventually lead to some deviations in the coefficients in the electromagnetic torque equations.

    To understand the change law of electromagnetic torque and the influence of thermal effects more clearly,numerical simulations of electromagnetic torque under different angles and temperatures were carried out.The model of the torque motor is shown in Fig.6.The permanent magnet material is LNGT60,which is characterized by a relatively small temperature coefficient and a small change in magnetic properties caused by temperature changes.At 40 ℃,the permanent magnet has a remanence of 0.9 T,a coercivity of-11 kA/m,and a relative permeability of 6.51.Both the magnet conductor and armature are made of iron-nickel soft magnetic alloy 1J50,which has relatively high permeability and low coercive force in low magnetic field environments.The structure of this alloy does not change significantly under the effect of temperature.There are 1100 turns of wire loops on each arm of the armature.

    Fig.6 Three-dimensional model of the torque motor used for numerical analysis

    As shown in Fig.7,the electromagnetic torque grows in an essentially linear manner during the deflection of the armature within the limits.At the same deflection angle,the electromagnetic torque decreases with increasing temperature.In the case of an armature deflection angle of 0.1°,the electromagnetic torque is 0.019835 N·m when the temperature is 40 ℃,decreasing to 0.018974 N·m when the temperature increases to 120 ℃.The amount of change in torque will be even smaller when the armature is within a deflection angle of 0.1°.According to the theoretical analysis,the armature deflection angle and the electromagnetic torque are linearly related,and the coefficient can be defined as the electromagnetic spring torque coefficient

    K

    .

    Fig.7 Simulation results of electromagnetic torque

    The expression for the electromagnetic moment derived from the theoretical analysis can be expressed as

    where

    K

    =-0.004837 N·m/(rad·℃) is the temperature coefficient of the electromagnetic moment,and coefficient

    K

    =11.61 N·m/rad.

    Accordingly,thermal effects can influence the electromagnetic torque by changing the structural parameters and the permanent magnet characteristics.The structural parameters directly affect the reluctance of the magnet conductor and the operating air gap,while the permanent magnet characteristics determine the magneto dynamic potential.

    4.2 Numerical simulation of the flow coefficient of the orifice

    The nozzle flapper servo valve relies mainly on a multi-stage throttling structure to achieve control of the output.In particular,the variable throttle orifice of the pilot stage realizes the adjustment of flow resistance by changing between the flapper and the nozzle.Analysis of the flow characteristics of the throttle structure is necessary.The first factor that should be considered for the numerical analysis of the flow field is the fluid properties.A nonlinear least-square method based on the trust-region algorithm is used to fit the available oil property data,and temperature coefficients of fluid viscosity and density are 0.0626 and 0.000575,respectively.

    For the analysis of the influence of the thermal effects of the flow field in the valve,it is necessary to consider not only the characteristics of the fluid,but also the variation of the key structural dimensions.Among the structural parameters are the diameter of the fixed throttle orifice,the diameter and length of the nozzle,and the width of the flapper.The variation patterns of these structural parameters are based on the equations in the theoretical calculation section.The 3D model created is shown in Fig.8,in which as many structural details as possible are preserved.

    Fig.8 Three-dimensional model of the fluid domain inside the valve

    Considering the symmetry of the structure,the calculation domain adopts a single nozzle structure.This can effectively reduce the computational burden caused by the excessive number of grids,and can collect data more conveniently.Considering irregular geometric shapes and large-scale changes,polyhedra elements are more suitable for meshing.By using the fixed type size function,a finer grid is constructed around the flapper-nozzle gap area and the fixed orifice.The results of grid independence analysis show that the adopted grid can obtain reasonable and acceptable results.The number of element grids was 1972034,and the skewness values were all less than 0.7(Fig.9).

    Fig.9 Meshing of the computational domain

    The corresponding boundary conditions were set according to the actual working conditions of the servo valve.There are three types of boundaries set in the calculation domain,namely the pressure inlet,wall boundary,and pressure outlet.Wall type boundaries confining the fluid or solid regions are defined at all other surfaces.For viscous flow,the no-slip stationary condition is applied by default.Simulated by FLUENT,the flow medium is aviation hydraulic oil,and the inlet pressure is set as 21 MPa during the calculation,while the outlet pressure is kept constant at 0.4 MPa.For turbulent behaviors,the set of turbulent kinetic energy and dissipation rate equations are solved using the standard

    k-ε

    model.In the pressure-velocity coupling,the SIMPLEC scheme is used since it allows the higher pressure-correction under relaxation factor that helps to speed up convergence.In spatial discretization,the body“PRESTO!”scheme is used for pressure discretization to work flexibly in a complex computational domain.The first-order upwind is applied for the momentum and the set of standard

    k-ε

    equations.

    Since the structure of the flow channel in the valve body is symmetrical,only half of each structure is used as the calculation domain to facilitate data collection,while also saving energy,computer memory,and calculation time.The inlet pressure and outlet pressure are set boundary conditions.The pressure loss between the control cavity pressure and the nozzle cavity can be ignored,and the control cavity pressure can be easily collected.Based on the flow continuity,the flow rate of the inlet is equal to that of the nozzle and fixed orifice,and its value can be obtained.According to Eqs.(22)and(23),the flow coefficients of the fixed orifice and the nozzle can be easily obtained.According to the flow equation at the throttle orifice,the curved surface of the flow coefficient of the fixed orifice and the nozzle under the combined action of different pressure differences and thermal effects can be obtained,as shown in Fig.10.

    Fig.10 Curved surface of the flow coefficient of the fixed orifice and the nozzle under the combined action of different pressure differences and thermal effects:(a) flow coefficient of the fixed throttle orifice (Cdf);(b) flow coefficient of the nozzle(Cdn)

    Based on the numerical analysis,expressions for the flow coefficients of fixed throttle orifices under the influence of thermal effects can be obtained:

    where the coefficients

    c

    =4.196×10,

    c

    =-1.226×10,

    c

    =1.202×10,and

    c

    =0.3875.

    With the change of the flapper displacement,the nozzle flow shows basically a linear trend.However,when the flapper is at the same offset,the thermal effects make the flow show nonlinear behavior.From the maximum flow rate at 80 ℃,temperature increases or decreases will reduce the flow rate through the nozzle.Based on the fitted curve,an expression for the nozzle flow coefficient under the influence of thermal effects can be obtained:

    where the coefficients

    k

    =-2.26×10,

    k

    =-5.11,

    k

    =1.47×10,

    k

    =-4×10,

    k

    =3.49×10,and

    k

    =0.5841.

    Therefore,the influence of environmental variables and structural parameters on flow characteristics can be obtained through numerical simulation.Thermal effects can affect the fluid flow state by changing fluid characteristics and structural parameters.The effect can be evaluated by numerical simulation of explicit expressions with physical concepts.

    4.3 Numerical simulation of the flow force on the flapper

    After the high-speed fluid in the nozzle enters the pilot stage volume cavity,the impact force on the flapper surface is the flow force.The components of the flow force are complex,and the theoretical analysis can be divided roughly into the static pressure on the nozzle projection area and the flow force generated by the vortex.

    The length of the pipeline is appropriately shortened without affecting the flow field solution.According to the results of the theoretical analysis,the stagnation pressure is generated by the flow force,nozzle diameter,and the pressure of the nozzle cavity,while the fluid momentum generated by the flow force is also related to the flow coefficient and the flapper-nozzle gap width.The nozzle flapper structure under the influence of thermal effects,and the deflection of the flapper have been parameterized in Solidworks.To obtain a clearer overview of the pressure distribution on the flapper surface,the pressure distribution on each side of the surface of the flapper was collected at different offset positions of the flapper(Fig.11).

    Fig.11 Pressure distribution on each side of the surface of the flapper at different temperatures when the displacement of the flapper xf=10 μm:(a)T=40 ℃,left surface;(b)T=40 ℃,right surface;(c)T=80 ℃,left surface;(d)T=80 ℃,right surface;(e)T=120 ℃,left surface;(f)T=120 ℃,right surface

    There is a distinct circular peak in the central region of the nozzle and a circular low-pressure area around the peak.When the flapper moves to the left,the pressure in the left nozzle chamber is higher and the flow rate at the outlet increases.Therefore,the surface stagnation pressure is also greater when the kinetic energy of the nozzle jet fluid increases.As the flapper offset increases,the peak on the left gradually increases and the peak on the right gradually decreases.Between the peak and the annular low pressure area is a pressure gradient area,in which the pressure drops smoothly and rapidly.The periphery of the annular low-pressure zone,that is,the remaining area of the flapper,can be called a flat area.The pressure in the flat area is slightly higher than the return pressure.

    The areas of the peak,ring,and flat regions remain essentially constant at different temperatures.However,the amplitudes of the peak and flat region gradually increase with rising temperature.In contrast,the pressure in the annular low pressure region reduces with increasing temperature.

    The numerical simulation results of the flow force for different displacements of the flapper and different temperatures are shown in Fig.12.At the same temperature,the variation of the flow force appears essentially linear.However,at the same offset,the effect of temperature on the flow dynamics is relatively large.In particular,at an offset of 10 μm,the flow force decreases rapidly from 0.575 to 0.152 N when changing the temperature from 40 to 60 ℃.The flow force is the lowest at 60 ℃and then shows a nonlinear increase with the gradient of the gradual decrease in growth.The exponential change between fluid viscosity and temperature characteristics greatly affects the fluid motion characteristics.This effect is more obvious when the gap between the nozzle and the flapper is large.

    Fig.12 Numerical simulation results of the flow force at different displacements of the flapper and different temperatures

    Therefore,the expression for the relationship between the flow force and temperature can be written as

    where the coefficients

    c

    =1.1×10,

    c

    =-4.09,

    c

    =

    547,

    c

    =-3.158×10,and

    c

    =6.823×10.

    The numerical analysis shows that the effect of temperature and flapper offset on the flow state of the fluid is significant.

    5 Analysis and evaluation

    In the previous sections,the numerical simulation results of complex links have realized the numerical fitting and parameter correction of the theoretical model.The purpose of this section is to synthesize the results of the above analyses and analyze the influence of thermal effects on the output of the servo valve by modeling the transfer function of the entire valve.A modular block diagram of the whole valve performance analysis is shown in Fig.13,which is divided into three main parts:the torque motor,pilot stage,and spool.The nomenclature of the simulation parameters of the servo valve is shown in Table 1.

    Table 1 Nomenclature table of simulation parameters

    Fig.13 A modular block diagram of the whole valve performance analysis

    The equilibrium equation for the torque on the armature-flapper assembly can be obtained by applying Newton’s second law and rigid body dynamics:

    where

    K

    is the coefficient of neutral electromagnetic torque,and

    F

    is the flow force on the flapper.

    Since the angle of rotation is relatively small,the relationship between the angle of rotation of the armature and the displacement can be approximated as

    The nozzle flapper pilot stage consists of the flapper,fixed orifice,nozzle,and oil return orifice.The nozzle and the movable flapper create two-variable gaps,forming two structures similar to variable orifices(Fig.14).According to the throttling equation,the flow rate through the fixed throttle orifice can be expressed as

    Fig.14 Schematic diagram of the structure and flow distribution of the valve body.pLv and pRv are the pressures in the control cavities on the left and right sides of the spool,respectively

    where

    A

    is the area of the throttle orifice.Assuming that the flapper is offset from the neutral position by

    x

    distance,the flow rate from each side of the nozzle can be expressed as

    Since the volume of the cavity has little effect on the pressure of the control cavity,the effect of temperature on the volume of the cavity is ignored.Considering the compressibility of the closed cavity at both ends of the valve spool,the pressure in the control cavity can be expressed as

    where

    v

    is the speed of the spool movement.

    Therefore,the control driving force of the valve spool can be expressed as

    The change of control force will affect the force balance state of the spool,which is also the factor with the most direct influence on the operating point drift of the servo valve under the thermal effect.

    Since the spool used in this study has an integrated feedback structure,a dynamic equilibrium equation of the spool can be established:

    Considering the oil compressibility,the relationship between the output flow rate and the displacement of flow rate of the servo valve output can be expressed as

    where

    C

    is the flow coefficient of the orifice of the spool,

    W

    is the area gradient,and

    x

    is the overlap of the slide valve.

    When there is no input control current signal,the static equilibrium equation of the armature-flapper assembly can be written as

    where

    K

    is the flow force stiffness coefficient obtained by numerical analysis.Assuming that the armatureflapper assembly maintains a vertical matching relationship,whether it is in the neutral position or the overall deflection,its position will be stable at the neutral position.

    We assumed that the armature,spring tube,and nozzle flapper had a deviation angle during processing or assembly.The deviation angle of these elements can be converted into an additional deflection angle of the armature.Therefore,the above formula can be written as

    When the pressure difference of the control cavity is enough to drive the spool to produce output pressure,the spool will close again under the action of the feedback force.To facilitate the analysis of the spool movement range,it can be assumed that the spool has a large overlap.When there is no input signal,the analysis of the displacement of the spool is realized by setting different deviation angles and temperatures.

    In the two cases of positive and negative deviation angles,the effect of temperature on the spool displacement is completely different.The

    z

    -axis in Fig.15 is the relative value of the steady-state position of the spool at the zero deviation angle.The direction of the deviation angle is the same as the movement direction of the spool.The greater the absolute value of the deviation angle,the greater the deviation of the spool deviation.This means that when the deviation angle of the armature assembly is positive,the dead zone of the servo valve will gradually increase as the temperature increases.Conversely,when the deviation angle is negative,the increase in temperature will cause the dead zone of the servo valve to decrease.When the deviation angle reaches 1°,in the range of 40-50 ℃,the displacement deviation decreases from 5.68×10m to 4.24×10m,and the tangent slope of the displacement deviation reaches 1.44×10m/℃.In the range of 110-120 ℃,the displacement deviation is reduced from 1.60×10m to 1.48×10m,and the tangent slope of the displacement deviation is only 1.25×10m/℃.

    Fig.15 Variation in the displacement deviation of the spool with the temperature and deviation angle

    When the torque motor has a 20-mA input signal,the deviation of the output pressure varies with the temperature and deviation angle as shown in Fig.16.The

    z

    -axis is the relative value of the output pressure at the zero deflection angle,which can be defined as the deviation pressure.A larger deviation angle will result in a greater deviation pressure.The deviation pressure is more obvious at lower temperatures.When the deviation angle is positive,the output pressure drop gradually decreases as the temperature increases.The reduction process is nonlinear,and the gradient is relatively large in the relatively low temperature stage.When the deviation angle is negative,the temperature increase will cause the output pressure to increase.This indicates that the variation of the operating point caused by the deviation angle in the positive and negative directions under the action of the thermal effect is opposite.When the deviation angle reaches 1°,in the range of 40-50 ℃,the pressure deviation decreases from 4.49×10Pa to 3.35×10Pa,and the tangent slope of the pressure deviation reaches 1.14×10Pa/℃.In the range of 110-120 ℃,the pressure deviation is reduced from 1.27×10Pa to 1.16×10Pa,and the tangent slope of the pressure deviation is only 110 Pa/℃.

    Fig.16 Variation in the deviation of the output pressure with the temperature and deviation angle

    6 Conclusions

    Based on the physical process of the servo valve,we first established a mathematical model of the key links under the action of the thermal effect to illustrate the mechanism of the thermal effect acting on the operating point of the servo valve.To evaluate the thermal effect,a mathematical model of the entire valve was established based on the results of numerical simulation.Based on the above analysis results,the main conclusions were as follows:

    (1)The thermal effect will affect the electromagnetic torque of the torque motor and the bending torque of the spring tube by changing the structural parameters.In addition,the matching relationship of the structure and the difference in the characteristics of the fluid will affect the flow force acting on the flapper.

    (2) Through the numerical simulation of each stage,the expression of the fitting relationship between input,output,and temperature was obtained.The influence of temperature on the electromagnetic torque output by the torque motor is about 5%,and its changing trend with temperature is relatively linear.The change trend of the discharge coefficient is more complicated,and the influence of temperature on the discharge coefficient of the nozzle is about 2%at most.Temperature will not only change the fluid characteristics,but also have a significant impact on the pressure difference and structural coordination.In addition,the difference of the above factors will also affect the flow force,which will be reduced to about 60%of the normal temperature.

    (3)According to the working principle of the servo valve,a mathematical model of the whole valve was established.On this basis,the influence of the thermal effect on the drift of the operating point was evaluated.When the deviation angle reaches±1°,within the range of 40-50 ℃,the tangent slope of the displacement deviation reaches 1.44×10m/℃,and the tangent slope of the pressure deviation within this range can reach 1.14×10Pa/℃.An increase in temperature will gradually reduce the slope of the tangent.In the range of 110-120 ℃,the above values will be reduced to 1.25×10m/℃and 110 Pa/℃,respectively.Therefore,we conclude that the deviation angle of the armatureflapper is an important factor affecting the drift of the operating point.

    In conclusion,the mechanism of the working point drift of servo valves under thermal effects is related to the structural and fluid characteristics.The deviation angle between components generated during the process of machining or assembly will directly lead to the working point drift.Moreover,the direction of the deviation angle will affect the drift direction of the operating point and the output amplitude.

    Acknowledgments

    This work is supported by the Civil Aircraft Research Project(No.MJ-2016-S-54),China.

    Author contributions

    Jian KANG designed the research and wrote the first draft of the manuscript.Jing-chao LI processed the corresponding data and helped to organize the manuscript.Zhao-hui YUAN revised and edited the final version.

    Conflict of interest

    Jian KANG,Zhao-hui YUAN,and Jing-chao LI declare that they have no conflict of interest.

    老鸭窝网址在线观看| 电影成人av| 免费黄色在线免费观看| 一级片'在线观看视频| 天天躁夜夜躁狠狠久久av| 成年女人在线观看亚洲视频| 精品久久久久久电影网| 欧美人与善性xxx| 乱人伦中国视频| 久久久久久久精品精品| 电影成人av| 黄片小视频在线播放| 欧美精品高潮呻吟av久久| 精品国产超薄肉色丝袜足j| 成人国语在线视频| 午夜福利在线观看免费完整高清在| 午夜福利乱码中文字幕| 亚洲国产欧美日韩在线播放| 免费看av在线观看网站| 欧美 日韩 精品 国产| 国产亚洲一区二区精品| 国产成人aa在线观看| 国产精品久久久久久精品古装| 久久久精品区二区三区| www.精华液| 最新中文字幕久久久久| 亚洲精品aⅴ在线观看| 亚洲经典国产精华液单| 亚洲欧美成人精品一区二区| 国产亚洲精品第一综合不卡| 亚洲成人一二三区av| 天美传媒精品一区二区| 中文字幕人妻熟女乱码| 国产av精品麻豆| 欧美av亚洲av综合av国产av | 亚洲国产色片| 性色avwww在线观看| 男女国产视频网站| 搡女人真爽免费视频火全软件| 久久av网站| 纯流量卡能插随身wifi吗| 天天躁日日躁夜夜躁夜夜| 久久久久人妻精品一区果冻| 精品一区二区三区四区五区乱码 | 亚洲视频免费观看视频| 日韩中文字幕欧美一区二区 | 欧美黄色片欧美黄色片| 国产亚洲精品第一综合不卡| 中文字幕亚洲精品专区| 中文字幕精品免费在线观看视频| 2022亚洲国产成人精品| 午夜福利乱码中文字幕| 999久久久国产精品视频| 国产成人aa在线观看| 不卡av一区二区三区| 一区二区三区乱码不卡18| 可以免费在线观看a视频的电影网站 | 在线观看免费视频网站a站| 美女福利国产在线| www.熟女人妻精品国产| 久久久国产一区二区| 欧美精品国产亚洲| 99国产综合亚洲精品| 国产1区2区3区精品| 精品国产乱码久久久久久小说| 极品少妇高潮喷水抽搐| 又黄又粗又硬又大视频| 日日摸夜夜添夜夜爱| av电影中文网址| 国产精品免费视频内射| 欧美人与性动交α欧美软件| 91久久精品国产一区二区三区| 亚洲精品国产一区二区精华液| 久久久久国产一级毛片高清牌| 国产精品一区二区在线观看99| 人人妻人人爽人人添夜夜欢视频| 在线观看免费视频网站a站| 一级爰片在线观看| 国产精品国产三级国产专区5o| 国产极品天堂在线| 热99久久久久精品小说推荐| 女人高潮潮喷娇喘18禁视频| 赤兔流量卡办理| 久久国产精品男人的天堂亚洲| 人人妻人人爽人人添夜夜欢视频| 精品国产乱码久久久久久男人| 成人亚洲精品一区在线观看| 女性生殖器流出的白浆| 尾随美女入室| 国产成人午夜福利电影在线观看| 久久久久久久亚洲中文字幕| 久久国产精品男人的天堂亚洲| 涩涩av久久男人的天堂| 亚洲国产看品久久| 国产在线一区二区三区精| 日韩中文字幕欧美一区二区 | 日本午夜av视频| 亚洲av成人精品一二三区| 大片电影免费在线观看免费| 亚洲内射少妇av| 美女国产视频在线观看| 在线精品无人区一区二区三| 国产精品久久久久久精品古装| 国产 精品1| 五月天丁香电影| 日本欧美国产在线视频| 在线 av 中文字幕| 91精品伊人久久大香线蕉| 你懂的网址亚洲精品在线观看| 成年美女黄网站色视频大全免费| 日韩人妻精品一区2区三区| 国产午夜精品一二区理论片| 精品国产一区二区三区四区第35| 欧美97在线视频| 在线观看www视频免费| 天堂8中文在线网| 亚洲av男天堂| 成人国产麻豆网| 最近最新中文字幕大全免费视频 | 亚洲精品视频女| av.在线天堂| 精品国产乱码久久久久久男人| 亚洲精品aⅴ在线观看| 最近最新中文字幕大全免费视频 | 免费久久久久久久精品成人欧美视频| 婷婷色综合大香蕉| 日本av免费视频播放| 一区二区三区精品91| 丝袜人妻中文字幕| 成人免费观看视频高清| 在线精品无人区一区二区三| 成年人免费黄色播放视频| 欧美激情极品国产一区二区三区| 亚洲,欧美精品.| 精品亚洲成国产av| 亚洲伊人色综图| 日韩一区二区三区影片| 亚洲一区中文字幕在线| 99久久精品国产国产毛片| 天堂8中文在线网| 色吧在线观看| 一区二区av电影网| 人成视频在线观看免费观看| 在线观看人妻少妇| 免费黄色在线免费观看| 亚洲欧美精品综合一区二区三区 | 免费黄色在线免费观看| 在线免费观看不下载黄p国产| 欧美精品亚洲一区二区| 精品久久久精品久久久| 国产免费一区二区三区四区乱码| 免费观看a级毛片全部| 久久午夜福利片| 男女高潮啪啪啪动态图| 成年女人毛片免费观看观看9 | 午夜日韩欧美国产| 中国国产av一级| 一级毛片电影观看| 大码成人一级视频| 中文字幕色久视频| 精品国产乱码久久久久久小说| 麻豆av在线久日| 国产精品国产三级专区第一集| 国产有黄有色有爽视频| 国产成人午夜福利电影在线观看| 老熟女久久久| 日韩中文字幕欧美一区二区 | 午夜福利乱码中文字幕| 国产在线一区二区三区精| 久久精品夜色国产| 国产乱来视频区| 黑人巨大精品欧美一区二区蜜桃| 建设人人有责人人尽责人人享有的| 免费看av在线观看网站| 精品少妇内射三级| 欧美日韩亚洲国产一区二区在线观看 | 国产男人的电影天堂91| 99国产精品免费福利视频| 人妻系列 视频| 亚洲婷婷狠狠爱综合网| 桃花免费在线播放| 亚洲美女视频黄频| av在线观看视频网站免费| videosex国产| 看免费成人av毛片| av视频免费观看在线观看| 精品人妻在线不人妻| 亚洲成人一二三区av| 国产乱来视频区| 婷婷成人精品国产| 精品99又大又爽又粗少妇毛片| 午夜影院在线不卡| 考比视频在线观看| 飞空精品影院首页| 黑人巨大精品欧美一区二区蜜桃| 飞空精品影院首页| 飞空精品影院首页| 久久精品国产a三级三级三级| 天天躁狠狠躁夜夜躁狠狠躁| 成人午夜精彩视频在线观看| 99国产精品免费福利视频| 久久久久久久久久人人人人人人| 制服丝袜香蕉在线| 亚洲激情五月婷婷啪啪| 视频在线观看一区二区三区| 母亲3免费完整高清在线观看 | 丝瓜视频免费看黄片| 天天躁日日躁夜夜躁夜夜| 国产综合精华液| 精品国产国语对白av| a级毛片黄视频| 免费黄色在线免费观看| 欧美精品亚洲一区二区| 亚洲国产成人一精品久久久| 亚洲国产看品久久| 黄片小视频在线播放| 久久国产精品大桥未久av| 伊人亚洲综合成人网| 中文字幕人妻熟女乱码| 国产欧美日韩综合在线一区二区| 日韩三级伦理在线观看| 国产精品一区二区在线不卡| 日韩成人av中文字幕在线观看| 久久影院123| 国产精品久久久久久av不卡| 亚洲av免费高清在线观看| 美女国产高潮福利片在线看| 免费看av在线观看网站| 久久久久久人人人人人| 亚洲av在线观看美女高潮| 亚洲欧美色中文字幕在线| 新久久久久国产一级毛片| 亚洲第一av免费看| 精品少妇黑人巨大在线播放| 国产av精品麻豆| 精品国产乱码久久久久久男人| a级毛片在线看网站| xxx大片免费视频| 天天躁日日躁夜夜躁夜夜| 亚洲av在线观看美女高潮| 七月丁香在线播放| 国产 一区精品| av福利片在线| 精品一区二区三区四区五区乱码 | 国产97色在线日韩免费| 少妇人妻久久综合中文| 国产精品人妻久久久影院| 免费久久久久久久精品成人欧美视频| 亚洲av电影在线观看一区二区三区| 丝袜美足系列| 国产精品久久久久久久久免| 女人高潮潮喷娇喘18禁视频| 亚洲精品一二三| 亚洲中文av在线| 亚洲精品久久久久久婷婷小说| 亚洲美女视频黄频| 国产精品人妻久久久影院| 国产综合精华液| 精品国产国语对白av| 亚洲成人一二三区av| 久久久久久人人人人人| 一本一本久久a久久精品综合妖精 国产伦在线观看视频一区 | 欧美av亚洲av综合av国产av | 亚洲精品一二三| 亚洲精品aⅴ在线观看| 男女边吃奶边做爰视频| 精品少妇一区二区三区视频日本电影 | 18在线观看网站| 成人国产av品久久久| 秋霞伦理黄片| 欧美日韩一区二区视频在线观看视频在线| 99热国产这里只有精品6| 黑人欧美特级aaaaaa片| 久久午夜综合久久蜜桃| 久久久久国产一级毛片高清牌| 一边亲一边摸免费视频| 91在线精品国自产拍蜜月| 水蜜桃什么品种好| 国产乱人偷精品视频| 国产欧美亚洲国产| freevideosex欧美| 国产男人的电影天堂91| 精品卡一卡二卡四卡免费| 欧美成人午夜免费资源| 蜜桃在线观看..| 国产精品香港三级国产av潘金莲 | av免费在线看不卡| 一本—道久久a久久精品蜜桃钙片| 久久午夜综合久久蜜桃| 日本欧美国产在线视频| 韩国精品一区二区三区| 伦精品一区二区三区| freevideosex欧美| 中文字幕色久视频| 国产视频首页在线观看| 久久国产亚洲av麻豆专区| 亚洲国产精品成人久久小说| av在线app专区| 精品酒店卫生间| 2022亚洲国产成人精品| 自拍欧美九色日韩亚洲蝌蚪91| 王馨瑶露胸无遮挡在线观看| 一区二区三区激情视频| 日本91视频免费播放| 国产一区二区三区综合在线观看| 久久免费观看电影| 中国三级夫妇交换| 免费观看性生交大片5| 免费黄色在线免费观看| 精品一区二区三卡| 色网站视频免费| 最近最新中文字幕免费大全7| 毛片一级片免费看久久久久| 伊人久久国产一区二区| 国产成人精品久久久久久| 一区二区三区精品91| 久久99蜜桃精品久久| 国产一区二区在线观看av| 搡老乐熟女国产| 国产一区二区 视频在线| 午夜激情久久久久久久| 国产精品久久久久成人av| 韩国精品一区二区三区| 国产一区二区三区综合在线观看| 最新的欧美精品一区二区| 免费久久久久久久精品成人欧美视频| 不卡av一区二区三区| 中国三级夫妇交换| 欧美另类一区| 又黄又粗又硬又大视频| 亚洲少妇的诱惑av| 夫妻性生交免费视频一级片| 日本av免费视频播放| 国产一区二区 视频在线| 久久人人爽人人片av| 人妻 亚洲 视频| 日本-黄色视频高清免费观看| 精品久久久久久电影网| 少妇被粗大猛烈的视频| www.熟女人妻精品国产| 国产精品免费视频内射| 久久这里有精品视频免费| 欧美精品一区二区大全| 亚洲成国产人片在线观看| 王馨瑶露胸无遮挡在线观看| 在线观看人妻少妇| 久久av网站| 1024视频免费在线观看| 亚洲国产av新网站| 精品第一国产精品| 一区二区三区激情视频| 热99久久久久精品小说推荐| 亚洲一区二区三区欧美精品| 激情五月婷婷亚洲| 在线免费观看不下载黄p国产| 亚洲人成网站在线观看播放| 婷婷色麻豆天堂久久| 多毛熟女@视频| 啦啦啦啦在线视频资源| 免费看不卡的av| 午夜av观看不卡| 色视频在线一区二区三区| 国产精品久久久久久精品古装| 国产av码专区亚洲av| 激情五月婷婷亚洲| 久久人人97超碰香蕉20202| 亚洲国产精品一区二区三区在线| 色视频在线一区二区三区| 久久久国产欧美日韩av| av有码第一页| 亚洲国产av影院在线观看| 免费观看性生交大片5| 999精品在线视频| 伊人久久国产一区二区| 色网站视频免费| 免费不卡的大黄色大毛片视频在线观看| 欧美人与性动交α欧美软件| 欧美97在线视频| 久久免费观看电影| 国产高清国产精品国产三级| 韩国高清视频一区二区三区| 亚洲第一青青草原| 日韩av不卡免费在线播放| 国产激情久久老熟女| 国产黄色视频一区二区在线观看| 国产午夜精品一二区理论片| 欧美日韩成人在线一区二区| 日韩成人av中文字幕在线观看| 亚洲美女搞黄在线观看| 精品国产露脸久久av麻豆| 久久久久视频综合| 在线亚洲精品国产二区图片欧美| 成人毛片60女人毛片免费| 中文字幕av电影在线播放| 亚洲第一av免费看| 天堂中文最新版在线下载| 国产男女内射视频| 水蜜桃什么品种好| 欧美+日韩+精品| 成人国产麻豆网| 免费在线观看视频国产中文字幕亚洲 | 日本-黄色视频高清免费观看| 搡女人真爽免费视频火全软件| 精品少妇久久久久久888优播| 日韩一区二区三区影片| 人人妻人人爽人人添夜夜欢视频| 久久精品国产自在天天线| 一级毛片我不卡| 午夜影院在线不卡| 美女大奶头黄色视频| 一级毛片电影观看| 蜜桃国产av成人99| 午夜福利在线免费观看网站| 亚洲av成人精品一二三区| 亚洲国产精品成人久久小说| 免费少妇av软件| 亚洲欧美精品自产自拍| 黑人欧美特级aaaaaa片| 亚洲视频免费观看视频| 欧美中文综合在线视频| 精品一品国产午夜福利视频| 精品久久久久久电影网| 久久久久久久国产电影| 日产精品乱码卡一卡2卡三| 久久韩国三级中文字幕| 久久国产精品男人的天堂亚洲| 亚洲人成电影观看| 最近中文字幕高清免费大全6| 久久久久久免费高清国产稀缺| 美女高潮到喷水免费观看| 美女福利国产在线| 卡戴珊不雅视频在线播放| 九九爱精品视频在线观看| 寂寞人妻少妇视频99o| 亚洲av电影在线观看一区二区三区| 男人舔女人的私密视频| 亚洲国产精品一区二区三区在线| 丰满迷人的少妇在线观看| 国产深夜福利视频在线观看| 99国产综合亚洲精品| 免费久久久久久久精品成人欧美视频| 青春草国产在线视频| 人人妻人人澡人人爽人人夜夜| 午夜激情久久久久久久| 在线亚洲精品国产二区图片欧美| 伊人亚洲综合成人网| 中文字幕av电影在线播放| 日本欧美视频一区| 国产精品成人在线| 咕卡用的链子| 日韩欧美一区视频在线观看| 制服诱惑二区| 男男h啪啪无遮挡| 亚洲国产看品久久| 久久精品国产自在天天线| 亚洲欧美清纯卡通| 亚洲视频免费观看视频| 亚洲美女视频黄频| 一级a爱视频在线免费观看| 十八禁网站网址无遮挡| 可以免费在线观看a视频的电影网站 | 亚洲成av片中文字幕在线观看 | 黄片播放在线免费| 午夜福利视频在线观看免费| 亚洲精品日韩在线中文字幕| 国产激情久久老熟女| 在线观看www视频免费| 精品国产乱码久久久久久男人| 中国国产av一级| 丝袜喷水一区| 免费播放大片免费观看视频在线观看| 国产精品嫩草影院av在线观看| 91精品三级在线观看| 777米奇影视久久| 晚上一个人看的免费电影| 亚洲欧美一区二区三区黑人 | 少妇被粗大的猛进出69影院| 人人妻人人澡人人看| 亚洲av免费高清在线观看| 女性生殖器流出的白浆| 日韩成人av中文字幕在线观看| 人妻一区二区av| 男女高潮啪啪啪动态图| h视频一区二区三区| 夫妻午夜视频| 久久久精品免费免费高清| 中文字幕色久视频| 男女午夜视频在线观看| 一本色道久久久久久精品综合| 少妇精品久久久久久久| 成人国语在线视频| 成年av动漫网址| 中文字幕人妻熟女乱码| 丝袜在线中文字幕| av在线观看视频网站免费| 少妇熟女欧美另类| 日本wwww免费看| 天天影视国产精品| 综合色丁香网| av.在线天堂| 青春草国产在线视频| 日韩制服骚丝袜av| 国产午夜精品一二区理论片| 精品少妇内射三级| 91午夜精品亚洲一区二区三区| 国产精品蜜桃在线观看| 亚洲色图 男人天堂 中文字幕| 婷婷成人精品国产| 看非洲黑人一级黄片| 妹子高潮喷水视频| 日本午夜av视频| 麻豆精品久久久久久蜜桃| 日韩一本色道免费dvd| 亚洲精品日本国产第一区| 亚洲国产精品一区二区三区在线| 在现免费观看毛片| 亚洲成人av在线免费| 99九九在线精品视频| 人体艺术视频欧美日本| videossex国产| 国产熟女欧美一区二区| 一本大道久久a久久精品| 一级毛片黄色毛片免费观看视频| 中文字幕精品免费在线观看视频| 香蕉丝袜av| 国产精品熟女久久久久浪| 国产福利在线免费观看视频| 曰老女人黄片| 亚洲一区二区三区欧美精品| 久久人人爽人人片av| 一区二区三区精品91| 999精品在线视频| 久久精品国产自在天天线| 男人舔女人的私密视频| 性少妇av在线| 天天躁狠狠躁夜夜躁狠狠躁| 中文精品一卡2卡3卡4更新| 午夜免费观看性视频| 亚洲欧洲精品一区二区精品久久久 | 亚洲伊人久久精品综合| 精品亚洲成a人片在线观看| 日韩中文字幕视频在线看片| 国产成人精品一,二区| 国产精品免费大片| 亚洲欧美成人精品一区二区| 国产一区二区 视频在线| 黄色一级大片看看| 咕卡用的链子| 美女大奶头黄色视频| 人妻系列 视频| 黄片无遮挡物在线观看| 一级黄片播放器| 日韩av不卡免费在线播放| 成年人午夜在线观看视频| 欧美日韩综合久久久久久| 久久97久久精品| 在线精品无人区一区二区三| 搡老乐熟女国产| av视频免费观看在线观看| www.精华液| 成年动漫av网址| 又粗又硬又长又爽又黄的视频| 午夜福利视频在线观看免费| 日韩三级伦理在线观看| 亚洲三级黄色毛片| 国产亚洲午夜精品一区二区久久| 午夜福利,免费看| 人人澡人人妻人| 免费在线观看完整版高清| 国产高清国产精品国产三级| 日韩欧美一区视频在线观看| 又粗又硬又长又爽又黄的视频| 一级爰片在线观看| 亚洲一级一片aⅴ在线观看| 国产野战对白在线观看| 精品久久久久久电影网| 久久人人97超碰香蕉20202| 最新中文字幕久久久久| 18禁动态无遮挡网站| 狂野欧美激情性bbbbbb| 国产av精品麻豆| 免费日韩欧美在线观看| 亚洲成人av在线免费| a级毛片黄视频| 亚洲精品视频女| 欧美成人精品欧美一级黄| 免费不卡的大黄色大毛片视频在线观看| 国产不卡av网站在线观看| 久久97久久精品| 久久人妻熟女aⅴ| 少妇 在线观看| 日韩中字成人| 国产亚洲精品第一综合不卡| 欧美中文综合在线视频| 午夜福利,免费看| 黄色怎么调成土黄色| 五月开心婷婷网| 亚洲精品国产av成人精品| 亚洲av中文av极速乱| 国产无遮挡羞羞视频在线观看| 亚洲av成人精品一二三区| 男人添女人高潮全过程视频| 两个人看的免费小视频| 国产成人精品在线电影| 男人爽女人下面视频在线观看| 丝袜人妻中文字幕| 久久精品国产a三级三级三级| 在线观看三级黄色| 日韩在线高清观看一区二区三区| 色视频在线一区二区三区| 最近中文字幕2019免费版| 成年人午夜在线观看视频| 国产麻豆69| 亚洲欧美成人精品一区二区| 亚洲成人手机|