• <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.

    国产伦在线观看视频一区| 亚洲人成77777在线视频| 日本 欧美在线| aaaaa片日本免费| 黄色a级毛片大全视频| 精品福利观看| 久久人人精品亚洲av| 国产av不卡久久| 一个人观看的视频www高清免费观看 | 亚洲avbb在线观看| 香蕉久久夜色| 好男人电影高清在线观看| 亚洲成av片中文字幕在线观看| 中亚洲国语对白在线视频| 中出人妻视频一区二区| 热re99久久国产66热| 黄色 视频免费看| 国产私拍福利视频在线观看| 大型av网站在线播放| 国产成人一区二区三区免费视频网站| 成人亚洲精品av一区二区| 欧美在线一区亚洲| 亚洲精品美女久久av网站| 免费av毛片视频| 国产精品一区二区三区四区久久 | 制服丝袜大香蕉在线| 色精品久久人妻99蜜桃| 欧美日韩福利视频一区二区| 国产视频内射| 亚洲精品国产一区二区精华液| 黑人欧美特级aaaaaa片| 怎么达到女性高潮| 丁香六月欧美| 亚洲avbb在线观看| 欧美三级亚洲精品| 亚洲国产欧美日韩在线播放| 国产一卡二卡三卡精品| 午夜福利在线观看吧| 国产精品久久电影中文字幕| 亚洲最大成人中文| 人人妻人人澡欧美一区二区| 亚洲av日韩精品久久久久久密| 麻豆成人午夜福利视频| 91老司机精品| 亚洲电影在线观看av| 亚洲成人国产一区在线观看| 最近最新免费中文字幕在线| 国产激情久久老熟女| 无限看片的www在线观看| 神马国产精品三级电影在线观看 | 最近在线观看免费完整版| 91老司机精品| 色哟哟哟哟哟哟| 免费在线观看完整版高清| 99riav亚洲国产免费| 一本精品99久久精品77| av福利片在线观看| 国产一区二区在线观看日韩| 日日啪夜夜撸| av在线老鸭窝| 天堂网av新在线| 99久久成人亚洲精品观看| 真人做人爱边吃奶动态| 国产成人一区二区在线| 高清午夜精品一区二区三区 | 极品教师在线视频| 卡戴珊不雅视频在线播放| 日韩欧美免费精品| 欧美一区二区亚洲| 观看美女的网站| 性插视频无遮挡在线免费观看| 禁无遮挡网站| 麻豆乱淫一区二区| 成年免费大片在线观看| 国产乱人偷精品视频| 看黄色毛片网站| 色综合色国产| av在线蜜桃| 国产一区二区三区av在线 | 国产一区二区亚洲精品在线观看| av在线播放精品| 精品一区二区三区视频在线| 99精品在免费线老司机午夜| 在线观看av片永久免费下载| 日韩强制内射视频| 日本撒尿小便嘘嘘汇集6| 小蜜桃在线观看免费完整版高清| 天天躁夜夜躁狠狠久久av| 久久久久久久午夜电影| 国产成人aa在线观看| 赤兔流量卡办理| 日本五十路高清| 国产精品女同一区二区软件| 国产精品久久久久久久久免| 日韩精品中文字幕看吧| 国产伦精品一区二区三区视频9| 亚洲高清免费不卡视频| 国产色爽女视频免费观看| 午夜福利在线观看免费完整高清在 | 深夜精品福利| 欧美色视频一区免费| 亚洲av中文字字幕乱码综合| 国产av在哪里看| 我要搜黄色片| 校园春色视频在线观看| 成人亚洲欧美一区二区av| 99热网站在线观看| 国产亚洲精品久久久com| 国产成人a区在线观看| 老熟妇仑乱视频hdxx| a级毛片免费高清观看在线播放| 中文字幕精品亚洲无线码一区| 91在线观看av| 亚洲欧美日韩无卡精品| 精品免费久久久久久久清纯| 一级av片app| 黄色一级大片看看| www.色视频.com| 国产亚洲av嫩草精品影院| 日韩大尺度精品在线看网址| 久久久久久九九精品二区国产| 久久久久久大精品| 成年免费大片在线观看| 国产精品国产高清国产av| 99热网站在线观看| 国产淫片久久久久久久久| 在线国产一区二区在线| 亚洲欧美日韩无卡精品| 性插视频无遮挡在线免费观看| 男女那种视频在线观看| 麻豆久久精品国产亚洲av| 成年女人永久免费观看视频| 中文字幕精品亚洲无线码一区| 国产探花极品一区二区| 久久亚洲精品不卡| 午夜免费激情av| 免费人成在线观看视频色| 蜜臀久久99精品久久宅男| 成人午夜高清在线视频| 乱码一卡2卡4卡精品| 免费在线观看影片大全网站| 嫩草影院精品99| 色av中文字幕| 99视频精品全部免费 在线| 天堂√8在线中文| 桃色一区二区三区在线观看| 伦理电影大哥的女人| 我的女老师完整版在线观看| 成人av在线播放网站| 国产精品不卡视频一区二区| 亚洲精品成人久久久久久| 国产精品久久视频播放| 午夜精品一区二区三区免费看| 国产爱豆传媒在线观看| 一夜夜www| 国产乱人偷精品视频| 国产高清视频在线播放一区| av在线天堂中文字幕| 亚洲第一区二区三区不卡| 久久精品国产自在天天线| 午夜福利视频1000在线观看| 国产激情偷乱视频一区二区| 午夜福利18| 日韩欧美精品免费久久| 一级毛片aaaaaa免费看小| 此物有八面人人有两片| 69人妻影院| 又爽又黄a免费视频| 成人av一区二区三区在线看| 黄色视频,在线免费观看| 午夜精品国产一区二区电影 | 12—13女人毛片做爰片一| 亚洲性夜色夜夜综合| 欧美极品一区二区三区四区| 全区人妻精品视频| 国产精品日韩av在线免费观看| 成人二区视频| 精品不卡国产一区二区三区| 变态另类丝袜制服| 91精品国产九色| 日本-黄色视频高清免费观看| 亚洲精品粉嫩美女一区| 久久精品国产99精品国产亚洲性色| 亚洲精品一卡2卡三卡4卡5卡| 亚洲av成人av| 男人和女人高潮做爰伦理| 久久精品国产鲁丝片午夜精品| 人人妻人人看人人澡| 亚洲激情五月婷婷啪啪| 欧美bdsm另类| 国产精品爽爽va在线观看网站| 欧美xxxx黑人xx丫x性爽| 国产精品久久久久久久久免| 欧美激情国产日韩精品一区| 日韩精品青青久久久久久| 99在线人妻在线中文字幕| 午夜免费激情av| 尤物成人国产欧美一区二区三区| 成年版毛片免费区| 美女免费视频网站| 久久综合国产亚洲精品| 国产69精品久久久久777片| 久久久久免费精品人妻一区二区| 亚洲精品亚洲一区二区| 露出奶头的视频| 精品99又大又爽又粗少妇毛片| 麻豆久久精品国产亚洲av| 一夜夜www| 成人美女网站在线观看视频| 97人妻精品一区二区三区麻豆| 一本久久中文字幕| 亚洲成av人片在线播放无| 一进一出好大好爽视频| 99热这里只有是精品在线观看| 在线观看免费视频日本深夜| a级毛片免费高清观看在线播放| 久久精品国产亚洲av天美| 国产高清视频在线播放一区| 99国产精品一区二区蜜桃av| 精品国内亚洲2022精品成人| 97人妻精品一区二区三区麻豆| 日韩亚洲欧美综合| 一个人看视频在线观看www免费| 丰满人妻一区二区三区视频av| 韩国av在线不卡| 在线天堂最新版资源| 国产精品精品国产色婷婷| 99久久成人亚洲精品观看| 99国产极品粉嫩在线观看| 国产精品久久视频播放| 亚洲成av人片在线播放无| 国产精品,欧美在线| 桃色一区二区三区在线观看| 成年女人看的毛片在线观看| 观看美女的网站| 免费观看精品视频网站| 国产一区二区在线观看日韩| 精品人妻熟女av久视频| 免费无遮挡裸体视频| 国产一区二区亚洲精品在线观看| 在线播放国产精品三级| 人人妻,人人澡人人爽秒播| 男插女下体视频免费在线播放| 熟妇人妻久久中文字幕3abv| 十八禁网站免费在线| 亚洲欧美精品综合久久99| 国产免费一级a男人的天堂| 久久6这里有精品| 日韩大尺度精品在线看网址| 亚洲欧美成人综合另类久久久 | 插阴视频在线观看视频| 国产精品99久久久久久久久| 少妇的逼好多水| 欧洲精品卡2卡3卡4卡5卡区| 久久久午夜欧美精品| 最后的刺客免费高清国语| 久久久久久久久久成人| 国产一区二区在线av高清观看| 直男gayav资源| av天堂在线播放| 午夜福利视频1000在线观看| 黄色一级大片看看| 中文字幕熟女人妻在线| 真实男女啪啪啪动态图| 久久人人精品亚洲av| 别揉我奶头~嗯~啊~动态视频| 国产精品国产高清国产av| or卡值多少钱| 国产在线精品亚洲第一网站| 国内揄拍国产精品人妻在线| 国产中年淑女户外野战色| 久久精品91蜜桃| 两个人视频免费观看高清| a级一级毛片免费在线观看| 男女做爰动态图高潮gif福利片| 最新在线观看一区二区三区| 成人av一区二区三区在线看| 午夜老司机福利剧场| 国产精品不卡视频一区二区| 十八禁网站免费在线| 亚洲图色成人| 亚洲美女黄片视频| 麻豆国产av国片精品| 国产精品1区2区在线观看.| 一级毛片久久久久久久久女| 真人做人爱边吃奶动态| 亚洲精品色激情综合| 成人亚洲精品av一区二区| 亚洲真实伦在线观看| 国产片特级美女逼逼视频| 精品不卡国产一区二区三区| 日日干狠狠操夜夜爽| 成人精品一区二区免费| 久久久久久久久久黄片| 成年av动漫网址| 无遮挡黄片免费观看| 成人性生交大片免费视频hd| 夜夜看夜夜爽夜夜摸| 免费人成视频x8x8入口观看| 欧美激情国产日韩精品一区| a级毛片a级免费在线| 久久久久久久亚洲中文字幕| 热99在线观看视频| 啦啦啦啦在线视频资源| 久久久久久九九精品二区国产| 亚洲精品国产av成人精品 | 久久午夜亚洲精品久久| 国产高清不卡午夜福利| 最近中文字幕高清免费大全6| 欧美高清成人免费视频www| 亚洲av成人av| av在线老鸭窝| 午夜日韩欧美国产| 亚洲熟妇熟女久久| 国产精品一区二区免费欧美| 色综合亚洲欧美另类图片| 亚洲av不卡在线观看| 日产精品乱码卡一卡2卡三| 亚洲国产精品成人久久小说 | 久久精品国产亚洲av香蕉五月| 91精品国产九色| 国产三级在线视频| 美女 人体艺术 gogo| 寂寞人妻少妇视频99o| 天天躁日日操中文字幕| 久久久久国产精品人妻aⅴ院| 欧美一区二区国产精品久久精品| 青春草视频在线免费观看| 丰满人妻一区二区三区视频av| 别揉我奶头 嗯啊视频| 美女大奶头视频| 久久精品国产亚洲网站| 在线观看av片永久免费下载| av女优亚洲男人天堂| 亚洲中文字幕日韩| 最后的刺客免费高清国语| h日本视频在线播放| 一夜夜www| 成人永久免费在线观看视频| 国产精品精品国产色婷婷| 波多野结衣巨乳人妻| 亚洲熟妇熟女久久| 人妻丰满熟妇av一区二区三区| 亚洲av.av天堂| 欧美国产日韩亚洲一区| 草草在线视频免费看| 亚洲人成网站在线播| 欧美日韩精品成人综合77777| 亚洲人与动物交配视频| 亚洲专区国产一区二区| 99热网站在线观看| 国产淫片久久久久久久久| 少妇熟女aⅴ在线视频| 3wmmmm亚洲av在线观看| 午夜精品在线福利| 亚洲欧美成人精品一区二区| 婷婷色综合大香蕉| 亚洲最大成人手机在线| 国产av麻豆久久久久久久| 精华霜和精华液先用哪个| 久久午夜亚洲精品久久| 国产伦在线观看视频一区| 国产精品免费一区二区三区在线| 亚洲成人久久爱视频| 不卡视频在线观看欧美| 亚洲欧美成人综合另类久久久 | 国产蜜桃级精品一区二区三区| 午夜激情福利司机影院| 我的老师免费观看完整版| 禁无遮挡网站| 亚洲一区二区三区色噜噜| 日本一二三区视频观看| 日本-黄色视频高清免费观看| 久久国内精品自在自线图片| 高清午夜精品一区二区三区 | 亚洲欧美日韩高清在线视频| 国产欧美日韩精品亚洲av| 免费观看精品视频网站| 国产一区二区在线av高清观看| av专区在线播放| 成年女人永久免费观看视频| 看黄色毛片网站| 国产一区二区三区在线臀色熟女| 国产精品不卡视频一区二区| 国产精品女同一区二区软件| 天天一区二区日本电影三级| 日韩av在线大香蕉| 99久久九九国产精品国产免费| 毛片一级片免费看久久久久| 小说图片视频综合网站| 免费观看的影片在线观看| 不卡一级毛片| 麻豆国产av国片精品| 国产高潮美女av| 狂野欧美白嫩少妇大欣赏| 91午夜精品亚洲一区二区三区| 狂野欧美白嫩少妇大欣赏| 午夜福利在线在线| 欧美中文日本在线观看视频| 亚洲精品日韩av片在线观看| 精品日产1卡2卡| 国产高清激情床上av| 高清毛片免费看| 亚洲av.av天堂| 99热这里只有是精品50| 亚洲内射少妇av| 久久久欧美国产精品| 中文字幕熟女人妻在线| 成人永久免费在线观看视频| 女人十人毛片免费观看3o分钟| 亚洲欧美日韩东京热| 色播亚洲综合网| 亚洲av二区三区四区| 性色avwww在线观看| 你懂的网址亚洲精品在线观看 | 欧美xxxx黑人xx丫x性爽| 国产视频一区二区在线看| 国内精品宾馆在线| 国产高潮美女av| 国产中年淑女户外野战色| 村上凉子中文字幕在线| 极品教师在线视频| 久久99热这里只有精品18| 一本精品99久久精品77| 免费高清视频大片| 欧美高清成人免费视频www| 亚洲av一区综合| 伦理电影大哥的女人| 日韩欧美精品免费久久| 成年女人看的毛片在线观看| 国产精品美女特级片免费视频播放器| 亚洲av中文av极速乱| 性欧美人与动物交配| 男女下面进入的视频免费午夜| 99久久成人亚洲精品观看| 国产精品人妻久久久影院| 国产精品野战在线观看| 黄色视频,在线免费观看| 国产精品乱码一区二三区的特点| 一边摸一边抽搐一进一小说| 人妻夜夜爽99麻豆av| 国产精品美女特级片免费视频播放器| 亚洲欧美日韩无卡精品| 99久久成人亚洲精品观看| 国产高清视频在线观看网站| 欧美日本视频| 天天躁日日操中文字幕| av在线蜜桃| 麻豆乱淫一区二区| 精品国产三级普通话版| 日韩大尺度精品在线看网址| 偷拍熟女少妇极品色| 老熟妇乱子伦视频在线观看| 在现免费观看毛片| 亚洲乱码一区二区免费版| 精品午夜福利视频在线观看一区| 美女黄网站色视频| 精品久久久久久久人妻蜜臀av| 嫩草影院精品99| 国产精品国产高清国产av| 一本久久中文字幕| 亚洲精品亚洲一区二区| 成人亚洲精品av一区二区| 丰满人妻一区二区三区视频av| 看片在线看免费视频| 搡老妇女老女人老熟妇| 全区人妻精品视频| 精品无人区乱码1区二区| 精品久久久久久久久亚洲| 精品99又大又爽又粗少妇毛片| 在线免费观看的www视频| 美女被艹到高潮喷水动态| 超碰av人人做人人爽久久| 亚洲av二区三区四区| 最近在线观看免费完整版| 亚洲精品一区av在线观看| 免费黄网站久久成人精品| 国产男人的电影天堂91| 欧美人与善性xxx| 欧美成人一区二区免费高清观看| 日韩人妻高清精品专区| 亚洲av电影不卡..在线观看| 精品午夜福利在线看| 午夜免费男女啪啪视频观看 | 亚洲国产高清在线一区二区三| 国产麻豆成人av免费视频| 深夜a级毛片| 中文字幕人妻熟人妻熟丝袜美| 久久久久久国产a免费观看| 伦精品一区二区三区| 国产成人a∨麻豆精品| 国产精品一区二区性色av| 日本黄色片子视频| av免费在线看不卡| 国产亚洲精品av在线| 性插视频无遮挡在线免费观看| 国产免费一级a男人的天堂| 干丝袜人妻中文字幕| 国产精品亚洲美女久久久| 嫩草影院入口| 日韩一本色道免费dvd| 成熟少妇高潮喷水视频| 长腿黑丝高跟| 特级一级黄色大片| 禁无遮挡网站| 午夜福利18| 一本精品99久久精品77| 嫩草影院精品99| 欧美最新免费一区二区三区| 人妻少妇偷人精品九色| 五月伊人婷婷丁香| 亚洲三级黄色毛片| 一夜夜www| 亚洲中文字幕一区二区三区有码在线看| 99热这里只有是精品在线观看| 色综合亚洲欧美另类图片| 午夜视频国产福利| 欧美在线一区亚洲| 你懂的网址亚洲精品在线观看 | 亚洲av免费高清在线观看| 日产精品乱码卡一卡2卡三| 国内久久婷婷六月综合欲色啪| 国国产精品蜜臀av免费| h日本视频在线播放| 国产精品美女特级片免费视频播放器| 最新中文字幕久久久久| 免费人成在线观看视频色| 两个人视频免费观看高清| 亚洲一区二区三区色噜噜| 一级av片app| 日本精品一区二区三区蜜桃| 日韩国内少妇激情av| 男人和女人高潮做爰伦理| 久久国产乱子免费精品| 天堂网av新在线| 最近中文字幕高清免费大全6| 寂寞人妻少妇视频99o| 成年女人毛片免费观看观看9| 国产精品一区二区三区四区久久| 免费看光身美女| 国内精品宾馆在线| 婷婷色综合大香蕉| 国产av麻豆久久久久久久| 卡戴珊不雅视频在线播放| 天天躁夜夜躁狠狠久久av| 亚洲不卡免费看| 成年女人永久免费观看视频| 波野结衣二区三区在线| 好男人在线观看高清免费视频| 悠悠久久av| 午夜激情欧美在线| 91午夜精品亚洲一区二区三区| 日韩 亚洲 欧美在线| 午夜福利在线在线| 韩国av在线不卡| 精品人妻偷拍中文字幕| 久久人人爽人人片av| 18禁黄网站禁片免费观看直播| 啦啦啦观看免费观看视频高清| 亚洲欧美日韩卡通动漫| 在线观看av片永久免费下载| 久久久久久久久久久丰满| 91久久精品国产一区二区三区| 欧美日韩综合久久久久久| 麻豆av噜噜一区二区三区| 五月玫瑰六月丁香| 黄色配什么色好看| 51国产日韩欧美| 99久久精品热视频| 色5月婷婷丁香| 在现免费观看毛片| 日韩精品有码人妻一区| 看片在线看免费视频| 欧美成人一区二区免费高清观看| 欧美日本亚洲视频在线播放| 欧美一级a爱片免费观看看| 欧美激情在线99| 欧美+亚洲+日韩+国产| 乱码一卡2卡4卡精品| 校园人妻丝袜中文字幕| 悠悠久久av| 国内久久婷婷六月综合欲色啪| 18禁在线播放成人免费| 三级国产精品欧美在线观看| 精品日产1卡2卡| 久久久久久大精品| 九九在线视频观看精品| 国产精品久久久久久久久免| 成年女人永久免费观看视频| 欧美区成人在线视频| 日本黄色片子视频| 国产国拍精品亚洲av在线观看| 亚洲av二区三区四区| 成年免费大片在线观看| 天堂√8在线中文| 波野结衣二区三区在线| 免费看av在线观看网站| 免费高清视频大片| 亚洲精品一卡2卡三卡4卡5卡| 欧美丝袜亚洲另类| 波多野结衣高清作品| 色综合色国产| a级毛色黄片| 长腿黑丝高跟| av天堂中文字幕网| 男人狂女人下面高潮的视频| 国产一区二区亚洲精品在线观看| 国产精品亚洲一级av第二区| 如何舔出高潮| 色综合站精品国产| 高清午夜精品一区二区三区 | 精品福利观看| 女人被狂操c到高潮| 哪里可以看免费的av片| 人人妻人人澡人人爽人人夜夜 |