CN103136450B - Aircraft surface erosion amount balancing method under supersonic speed - Google Patents

Abstract

The invention discloses aircraft surface erosion amount balancing method under a kind of supersonic speed, can based on the flow field situation on Euler method calculating aircraft surface and mass particle distribution situation, calculate aircraft surface erosion amount situation, the three-dimensional aircraft erosion amount that under the supersonic speed designed by the present invention, aircraft surface erosion amount balancing method calculates has good accuracy.

Description

Method for measuring surface erosion amount of aircraft at supersonic speed

Technical Field

The invention relates to a method for calculating the surface erosion amount of an aircraft at supersonic speed.

Background

When an aircraft flies in the atmosphere at supersonic speed, particles such as dust, raindrops, ice crystals, snowflakes and the like can collide with each other, on one hand, the particles influence the boundary layer characteristics of the collision surface to increase convective heat transfer, on the other hand, the particles directly collide at high speed to cause the quality loss of the collision surface, the phenomenon is called as erosion phenomenon, the safety of the aircraft flying at the hypersonic speed is seriously damaged by the erosion of the surface, and therefore the erosion phenomenon is an important subject of the hypersonic speed protection design of the aircraft.

According to classical erosion theory, the erosion receding rate is believed to be a function of particle velocity, particle density, angle of incidence, and the properties of the material. The erosion calculation model that can be seen at present can be written in the form:

S_e＝E_p/(ρ_tC_N)

$E_p=5\×{10}^{-4}\[{\mathrm{\ρ}}_{P_{\mathrm{\∞}}}{V}_{P_{\mathrm{\∞}}}\left(\frac{m_P}{m_{P_{\mathrm{\∞}}}}\right)V_P^2{\sin }^2\mathrm{\θ}\]$

wherein,which represents the density of the particles of the incoming flow,representing the incoming particle velocity, V_PThe speed of the particles at the object plane is shown,denotes the initial mass of the particle, m_PRepresenting the object plane particle mass, θ is the object plane angle.

The above formula can simulate the solid particle erosion condition more effectively, but also has limitations: (1) the formula is suitable for calculating the particle orbit by the Lagrange method, the Lagrange method has advantages when the trace is in two dimensions, but in the simulation of a three-dimensional entity, the Eulerian method is a frequently adopted method, and the formula is difficult to realize in a three-dimensional flow field; (2) sin theta in the formula has better adaptability to specific two-dimensional conditions, and needs to be changed under three-dimensional conditions; (3) the formula does not take into account the impact of the impingement surface on the gas flow and there is an error in calculating the erosion.

Disclosure of Invention

The technical problem to be solved by the invention is to provide a method for measuring the surface erosion amount of the aircraft at supersonic speed, which can calculate the surface erosion amount of the aircraft based on the flow field condition and the particle mass distribution condition of the aircraft surface calculated by the Euler method.

The invention adopts the following technical scheme to solve the technical problems: the invention designs a method for measuring the surface erosion amount of an aircraft at supersonic speed, which comprises the following specific steps:

step (1): acquiring the speed V, the density rho and the ambient pressure p of the air flowing on the surface of the aircraft in a supersonic speed state;

step (2): according to the mass and momentum conservation principle, the velocity V, the density rho and the ambient pressure p of the incoming air collected in the step (1) are used for listing the following three-dimensional Reynolds time-average equation, and the air velocity u at any position in the space is obtained by solving the equation:

$\frac{\∂\mathrm{\ρ}}{\∂t}+\frac{\∂}{\∂x_i}\left({\mathrm{\ρu}}_i\right)=0$

$\frac{\∂}{\∂t}\left({\mathrm{\ρu}}_i\right)+\frac{\∂}{\∂x_j}\left({\mathrm{\ρu}}_i{u}_j\right)=\frac{\∂\mathrm{\ρ}}{\∂x_i}+\frac{\∂}{\∂x_j}(\mathrm{\μ}\frac{\∂u_i}{\∂x_j}-\mathrm{\ρ}\stackrel{\‾}{u_i^{\′}{u}_j^{\′}})$

wherein, subscript i, j takes on 1, 2, 3, respectively representing x, y, z direction, i.e. u₁Denotes the component velocity of the air in the x-direction, u₂Representing the component velocity of the air in the y-direction, u₃Represents the component velocity of the air in the z direction; x is the number of₁Representing a parameter on the x-axis, x₂Representing a parameter, x, on the y-axis₃Representing a parameter in the z-axis; μ represents the kinetic viscosity of air;

the derivative with respect to time is represented as,the derivative of x is represented as a function of,the derivative to y is represented as a function of,represents the derivative to z;

and (3) based on the air velocity calculation result in the step (2), listing mass and momentum equations for the particles and solving to obtain the volume fraction α of the particles and the velocity u of the particles_s：

$\frac{\∂\left({\mathrm{\ρ}}_s\mathrm{\α}\right)}{\∂t}+\▿\·\left({\mathrm{\ρ}}_s{\mathrm{\αu}}_s\right)=0$

$\frac{\∂\left({\mathrm{\ρ}}_s{\mathrm{\αu}}_s\right)}{\∂t}+\▿\·\left({\mathrm{\ρ}}_s{\mathrm{\αu}}_s{u}_s\right)=K\mathrm{\α}\mathrm{\ρ}(u-u_s)+\mathrm{\α}\mathrm{\ρ}F$

Where ρ is_sDefined as the density of the solid particles, u_sDefined as the particle velocity, F as the gravitational acceleration, K as the inertia factor, α as the particle volume fraction at any position in space, u_aDefined as the air velocity;

and (4): calculating the erosion receding rate by using the following formula according to the particle volume fraction alpha and the particle speed result obtained by calculation in the step (3):

Se＝kE_knV/(ρ_tC_N)

$E_{kn}={\mathrm{\α}}_n{V}_n^2$

wherein S is_eDefined as the erosion receding Rate, E_knDefined as the normal momentum, p, of the particle striking the surface_tDefined as the density of the impinging material, C_NDefined as the erosion coefficient of the material, and k is defined as a parameter related to the profile of the eroded surface, α_nDefined as the volume fraction of particles striking the surface, V, in the calculation of the trajectory of the solid particles_nDefined as the normal velocity of the particle when it strikes the surface.

Compared with the prior art, the invention has the following advantages:

the method for measuring the erosion amount of the surface of the aircraft at the supersonic speed can calculate the mass distribution condition of the particles on the surface of the aircraft by using the Euler method, is convenient to popularize to the three-dimensional condition, and has better accuracy.

Drawings

FIG. 1 is a comparison graph of erosion amount simulation calculation results.

Detailed Description

The invention is further described in detail below with reference to the accompanying drawings:

when calculating the three-dimensional erosion situation, the lagrangian method is inconvenient to calculate, and therefore, the main calculation method is the eulerian method, but a formula for calculating the erosion amount of the three-dimensional eulerian method is not seen yet, so that a new erosion amount calculation method is needed.

The invention designs a method for measuring the surface erosion amount of an aircraft at supersonic speed, which comprises the following specific steps:

step (1): acquiring the speed V, the density rho and the ambient pressure p of the air flowing on the surface of the aircraft in a supersonic speed state;

step (2): according to the mass and momentum conservation principle, the velocity V, the density rho and the ambient pressure p of the incoming air collected in the step (1) are used for listing the following three-dimensional Reynolds time-average equation, and the air velocity u at any position in the space is obtained by solving the equation:

$\frac{\∂\mathrm{\ρ}}{\∂t}+\frac{\∂}{\∂x_i}\left({\mathrm{\ρu}}_i\right)=0$

$\frac{\∂}{\∂t}\left({\mathrm{\ρu}}_i\right)+\frac{\∂}{\∂x_j}\left({\mathrm{\ρu}}_i{u}_j\right)=\frac{\∂\mathrm{\ρ}}{\∂x_i}+\frac{\∂}{\∂x_j}(\mathrm{\μ}\frac{\∂u_i}{\∂x_j}-\mathrm{\ρ}\stackrel{\‾}{u_i^{\′}{u}_j^{\′}})$

wherein, subscript i, j takes on 1, 2, 3, respectively representing x, y, z direction, i.e. u₁Denotes the component velocity of the air in the x-direction, u₂Representing the component velocity of the air in the y-direction, u₃Represents the component velocity of the air in the z direction; x is the number of₁Representing a parameter on the x-axis, x₂Representing a parameter, x, on the y-axis₃Representing a parameter in the z-axis; μ represents the kinetic viscosity of air;

to representThe derivative with respect to time is that of,the derivative of x is represented as a function of,the derivative to y is represented as a function of,represents the derivative to z;

establishing a three-dimensional unsteady Navier-Stokes equation according to three conservation relations of mass, momentum and energy, and writing the equation into an integral form aiming at a certain control body omega as follows:

wherein Q is a conservative variable, FC is a non-viscous flux, Fv is a viscous flux, n is a control body in vitro vector, and S is a control body boundary:

in viscous flux, the expressions for stress and heat flux are:

${\mathrm{\Θ}}_x={\mathrm{u\τ}}_{xx}+{\mathrm{\ν\τ}}_{xy}+k\frac{\∂T}{\∂x}$

${\mathrm{\Θ}}_y={\mathrm{u\τ}}_{yx}+{\mathrm{\ν\τ}}_{yy}+k\frac{\∂T}{\∂y}$

${\mathrm{\τ}}_{xx}=\frac{2}{3}\mathrm{\μ}(2\frac{\∂u}{\∂x}-\frac{\∂v}{\∂y}-\frac{\∂w}{\∂z}){\mathrm{\τ}}_{xy}={\mathrm{\τ}}_{yx}=\mathrm{\μ}(\frac{\∂v}{\∂x}+\frac{\∂u}{\∂y})$

${\mathrm{\τ}}_{yy}=\frac{2}{3}\mathrm{\μ}(2\frac{\∂v}{\∂y}-\frac{\∂u}{\∂x}-\frac{\∂w}{\∂z}){\mathrm{\τ}}_{yz}={\mathrm{\τ}}_{zy}=\mathrm{\μ}(\frac{\∂w}{\∂y}+\frac{\∂v}{\∂z})$

${\mathrm{\τ}}_{zz}=\frac{2}{3}\mathrm{\μ}(2\frac{\∂w}{\∂z}-\frac{\∂u}{\∂x}-\frac{\∂v}{\∂y}){\mathrm{\τ}}_{zx}={\mathrm{\τ}}_{xz}=\mathrm{\μ}(\frac{\∂u}{\∂z}+\frac{\∂w}{\∂x})$

in order to make the N-S equation system closed, some physical relations are also needed.

For a complete gas, there is the equation of state:

P＝ρRT

$\mathrm{\ρ}=(\mathrm{\γ}-1)\mathrm{\ρ}\[e-\frac{1}{2}(u^2+v^2+w^2)\]$

let the spatial derivatives qx, qy, qz be solved by the following system of equations:

$\begin{array}{c}\left[\begin{array}{ccc}{x}_{c2}-x_{c1}& y_{c2}-y_{c1}& z_{c2}-z_{c1}\\ 0.5\left(x_{n2}+x_{n3}\right)-x_{n1}& 0.5\left(y_{n2}+y_{n3}\right)-y_{n1}& 0.5\left(z_{n2}+z_{n3}\right)-z_{n1}\\ 0.5\left(x_{n3}+x_{n1}\right)-x_{n2}& 0.5\left(y_{n3}+y_{n1}\right)-x_{n2}& 0.5\left(z_{n3}+z_{n1}\right)-z_{n2}\end{array}\right]\left[\begin{array}{c}{q}_x\\ q_y\\ q_z\end{array}\right]\\ =\left[\begin{array}{c}{q}_{c2}-q_{c1}\\ 0.5\left(q_{n2}+q_{n3}\right)-q_{n1}\\ 0.5\left(q_{n3}+q_{n1}\right)-q_{n2}\end{array}\right]\end{array}$

in the above equation set, the flow parameters to the respective grid nodes are used, which solve the formula as follows:

$q_n=\left(\underset{i=1}{\overset{N}{\Σ}}\frac{q_{0,i}}{r_i}\right)/\left(\underset{i=1}{\overset{N}{\Σ}}\frac{1}{r_i}\right)$

wherein: r is_i＝(x_0,ix_n)²+(y_0,i-y_n)²+(z_0,i-z_n)²]^1/2

q_0,,iIs the flow parameter at the grid center where node n is located, ri is the distance from point n to the grid center point.

The viscosity coefficient is composed of a laminar viscosity coefficient and a turbulent viscosity coefficient, and the determination of the vortex viscosity coefficient is illustrated by taking an RNG two-equation turbulence model as an example. Assuming that the Boussinesq vortex viscosity relationship is satisfied between the turbulent shear stress and the average strain rate tensor, the Reynolds stress can be written as:

$-\mathrm{\ρ}\stackrel{\‾}{u_i^{\′}{u}_j^{\′}}={\mathrm{\μ}}_l\left(\frac{\∂u_i}{\∂x_j}+\frac{\∂u_j}{\∂x_i}\right)-\frac{2}{3}{\mathrm{\ρ\δ}}_{ij}k$

according to dimensional analysis, μ_tCan be represented by k and:

${\mathrm{\μ}}_t={\mathrm{\ρC}}_{\mathrm{\μ}}\frac{k^2}{\mathrm{\ϵ}}$

wherein C is_μAn empirical constant of 0.09; k is turbulent pulsation kinetic energy; is the dissipation ratio of turbulent pulsating kinetic energy. Transport method of turbulent flow pulsation kinetic energy kThe process is as follows:

$\frac{\∂\left(\mathrm{\ρ}k\right)}{\∂t}+\frac{\∂\left({\mathrm{\ρu}}_{i}k\right)}{\∂x_j}=\frac{\∂}{\∂x_j}\[{\mathrm{\α}}_k(\mathrm{\μ}+{\mathrm{\μ}}_l)\frac{\∂k}{\∂x}\]+G_k+\mathrm{\ρ}\mathrm{\ϵ}$

$\frac{\∂\left(\mathrm{\ρ}\mathrm{\ϵ}\right)}{\∂t}+\frac{\∂\left({\mathrm{\ρ\ϵu}}_i\right)}{\∂x_i}=\frac{\∂}{\∂x_i}\[{\mathrm{\α}}_{\mathrm{\ϵ}}(\mathrm{\μ}+{\mathrm{\μ}}_l)\frac{\∂\mathrm{\ϵ}}{\∂x_i}\]+C_{\mathrm{\ϵ}1}\frac{\mathrm{\ϵ}}{k}{P}_k-C_{\mathrm{\ϵ}2}\mathrm{\ρ}\frac{{\mathrm{\ϵ}}^2}{k}$

after introducing a transport equation, an N-S equation is closed, a central difference format is adopted in spatial dispersion, and a multi-step Runge-Kutta format is adopted in time dispersion. The solution adopts a famous Jameson center difference method, and has two characteristics: on spatial dispersion, the flux at the inner field edge is obtained by averaging the flux values at the center of the left and right cells; in order to better capture discontinuity points and improve the stability of a format, an artificial dissipation term is added in flux calculation, the artificial dissipation term consists of second-order and fourth-order difference terms of conservation variables, the calculation of a mass, momentum and energy three-large fluid control equation is carried out on the basis of a structural grid, a constant solution is obtained by multi-step Runge-Kutta iteration in time, in addition, an accelerated convergence and residual fairing technology is adopted, the CFL number can be increased to 6-8, the convergence speed is improved by 50%, the convergence precision is controlled to be 10-8 orders of magnitude, and the far-field undisturbed boundary condition is set to be 15-fold chord length, namely a non-reflection boundary condition.

In order to obtain reliable flow field calculation parameters, the developed flow field solver is mainly based on an Euler or Navier-Stokes equation, an S-A model, A k-model and the like can be selected when A viscous flow field is solved, an explicit format-Runge-KuttA (Runge-KuttA) multi-step format and A multi-step fourth-order Runge-KuttA method can be adopted for solving the constant differential equation of time, and the format is as follows:

$Q_{I,0}=Q_I^k$

Q_I,1＝Q_I,0+α₁ΔtR_I,0

Q_I,2＝Q_I,0+α₂ΔtR_I,1

Q_I,3＝Q_I,0+α₃ΔtR_I,2

Q_I,4＝Q_I,0+α₄ΔtR_I,3

$Q_I^{k+1}=Q_{I,4}$

wherein: ${\mathrm{\α}}_1=\frac{1}{4},{\mathrm{\α}}_2=\frac{1}{3},{\mathrm{\α}}_3=\frac{1}{2},{\mathrm{\α}}_4=1.$

for the explicit propulsion format of the Runge-Kutta method, the time step is limited for stability requirements. However, since the unsteady equation is used to solve the unsteady flow problem, and the time advance and the space dispersion are independent in the solution process, the unsteady solution is independent of the time step, so that the local time step can be adopted on each grid unit, and the flow field is advanced with the time step close to the stability limit everywhere, thereby accelerating the convergence process of the whole iterative computation.

The time step per grid cell can be calculated using the following formula:

${\mathrm{\Δt}}_i=\frac{1}{4}\underset{n=1}{\overset{N}{\Σ}}\frac{CFL\·{\mathrm{\Ω}}_i}{|{\mathrm{uS}}_x+{\mathrm{vS}}_y+{\mathrm{wS}}_z|+c_i\sqrt{{S_x}^2+{S_y}^2+{S_z}^2}}$

in the above formula, Ω is the volume of the grid, u, v, w are velocity components, SX, SY, SZ are projection components of the area along the coordinate axis, c is the local sound velocity, and the value of the CFL number is determined by the specific flow condition. The explicit format has the advantages of small calculation amount and storage amount and simple program every time step is advanced. The disadvantages are that the time step length is limited by stability conditions, the number of calculation is too small, the efficiency is low, the number of steps required for convergence is large, and the total calculation time is very long, so that the characteristics of the explicit format can be described by 'small step fast running'.

And (3) based on the air velocity calculation result in the step (2), listing mass and momentum equations for the particles and solving to obtain the volume fraction α of the particles and the velocity u of the particles_s：

$\frac{\∂\left({\mathrm{\ρ}}_s\mathrm{\α}\right)}{\∂t}+\▿\·\left({\mathrm{\ρ}}_s{\mathrm{\αu}}_s\right)=0$

$\frac{\∂\left({\mathrm{\ρ}}_s{\mathrm{\αu}}_s\right)}{\∂t}+\▿\·\left({\mathrm{\ρ}}_s{\mathrm{\αu}}_s{u}_s\right)=K\mathrm{\α}\mathrm{\ρ}(u-u_s)+\mathrm{\α}\mathrm{\ρ}F$

Where ρ is_sDefined as the density of the solid particles, u_sDefined as the particle velocity, F as the gravitational acceleration, K as the inertia factorα is defined as the volume fraction of particles, u, at any position in space_aDefined as the air velocity;

and (4): calculating the erosion receding rate by using the following formula according to the particle volume fraction alpha and the particle speed result obtained by calculation in the step (3):

S_e＝kE_knV/(ρ_tC_N)

$E_{kn}={\mathrm{\α}}_n{V}_n^2$

wherein S is_eDefined as the erosion receding Rate, E_knDefined as the normal momentum, p, of the particle striking the surface_tDefined as the density of the impinging material, C_NDefined as the erosion coefficient of the material, and k is defined as a parameter related to the profile of the eroded surface, α_nDefined as the volume fraction of particles striking the surface, V, in the calculation of the trajectory of the solid particles_nDefined as the normal velocity of the particle when it strikes the surface.

In order to verify the feasibility of the method, the erosion amount of a specific shape is calculated by the method and is compared with the test result under the same working condition and the equation simulation result in budget, as shown in fig. 1, the comparison graph of the erosion amount result calculated by the formula and the erosion amount result calculated by the original formula and the test result under the same working condition is shown in the graph, and the calculation result of the method is closer to the test value than the calculation result of the original method in the overall trend.

At the geometric leading edge point (Y is 0.00 m), the experimental value is 0.02544m, the calculated value of the original formula is 0.0299m, and the calculated value of the patent is 0.02722m, so that the method adopted by the patent is slightly smaller than the calculated value of the original formula, is closer to the experimental result, and is in good agreement.

Claims (1)

1. A method for measuring the surface erosion amount of an aircraft at supersonic speed is characterized by comprising the following specific steps:

step (1): acquiring the speed V, the density rho and the ambient pressure p of the air flowing on the surface of the aircraft in a supersonic speed state;

step (2): according to the mass and momentum conservation principle, the velocity V, the density rho and the ambient pressure p of the incoming air collected in the step (1) are used for listing the following three-dimensional Reynolds time-average equation, and the air velocity u at any position in the space is obtained by solving the equation:

$\frac{\∂\mathrm{\ρ}}{\∂t}+\frac{\∂}{\∂x_i}\left({\mathrm{\ρu}}_i\right)=0$

$\frac{\∂}{\∂t}\left({\mathrm{\ρu}}_i\right)=+\frac{\∂}{\∂x_j}\left({\mathrm{\ρu}}_i{u}_j\right)=\frac{\∂\mathrm{\ρ}}{\∂x_i}+\frac{\∂}{\∂x_j}(\mathrm{\μ}\frac{\∂u_i}{\∂x_j}-\mathrm{\ρ}\stackrel{\‾}{u_i^{\′}{u}_j^{\′}})$

wherein, subscript i, j takes on 1, 2, 3, respectively representing x, y, z direction, i.e. u₁Denotes the component velocity of the air in the x-direction, u₂Representing the component velocity of the air in the y-direction, u₃Represents the component velocity of the air in the z direction; x is the number of₁Representing a parameter on the x-axis, x₂Representing a parameter, x, on the y-axis₃Representing a parameter in the z-axis; μ represents the kinetic viscosity of air;

the derivative with respect to time is represented as,the derivative of x is represented as a function of,the derivative to y is represented as a function of,represents the derivative to z;

and (3) based on the air velocity calculation result in the step (2), listing mass and momentum equations for the particles and solving to obtain the volume fraction α of the particles and the velocity u of the particles_s：

$\frac{\∂\left({\mathrm{\ρ}}_s\mathrm{\α}\right)}{\∂t}+\▿\·\left({\mathrm{\ρ}}_s{\mathrm{\αu}}_s\right)=0$

$\frac{\∂\left({\mathrm{\ρ}}_s{\mathrm{\αu}}_s\right)}{\∂t}+\▿\·\left({\mathrm{\ρ}}_s{\mathrm{\αu}}_s{u}_s\right)=K\mathrm{\α}\mathrm{\ρ}(u-u_s)+\mathrm{\α}\mathrm{\ρ}F$

Where ρ is_sDefined as the density of the solid particles, u_sDefining as particle velocity, F as gravitational acceleration, K as an inertia factor, α as particle volume fraction at any position in space, u as air velocity;

and (4): calculating the erosion receding rate by using the following formula according to the particle volume fraction alpha and the particle speed result obtained by calculation in the step (3):

S_e＝kE_knV/(ρ_tC_N)

$E_{kn}={\mathrm{\α}}_n{V}_n^2$

wherein S is_eDefined as the erosion receding Rate, E_knDefined as the normal momentum, p, of the particle striking the surface_tDefined as the density of the impinging material, C_NDefined as the erosion coefficient of the material, and k is defined as a parameter related to the profile of the eroded surface, α_nDefined as the volume fraction of particles striking the surface, V, in the calculation of the trajectory of the solid particles_nDefined as the normal velocity of the particle when it strikes the surface.