CN103136450A - Method for measuring erosion amount of aircraft surface under supersonic speed - Google Patents

Abstract

The invention discloses a method for measuring erosion amount of an aircraft surface under supersonic speed. Conditions of the erosion amount of the aircraft surface can be calculated on the basis of calculation of flow field conditions and grain mass distribution conditions of the aircraft surface in an Eulerian method. According to the method for measuring the erosion amount of the aircraft surface under the supersonic speed, the calculated erosion amount of a three-dimensional aircraft is good in accuracy.

Description

Aircraft surface erosion amount balancing method under supersonic speed

Technical field

The present invention relates to the computing method of aircraft surface erosion amount under a kind of supersonic speed.

Background technology

When aircraft flies in atmospheric envelope with supersonic speed, can run into dust, raindrop, ice crystal, the particle hits such as snowflake, these particles affect the impact surface behavior of boundary layer on the one hand, increase convection heat transfer, the direct high-speed impact of particle on the other hand, cause the mass loss of impact surface, we call erosion to this phenomenon, the surface the erosion serious harm security of aircraft in hypersonic lower flight, therefore, erosion is the important topic of the hypersonic protection Design of aircraft.

Theoretical according to the erosion of classics, think that the erosional retreat rate is particle rapidity, particle density, the function of incident angle and properties of materials.The erosion computation model that can see at present can be write as following form:

S _e＝E _p/(ρ _tC _N)

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

Wherein,

Expression incoming flow particle density,

Expression incoming flow particle rapidity, V _PExpression object plane particle rapidity,

Expression particle initial mass, m _PExpression object plane mass particle, θ is the object plane angle.

Above-mentioned formula can comparatively effectively be simulated the solids erosion condition, but limitation is arranged also: (1) formula is applicable to Lagrangian method and calculates the particle rail, mark is when two dimension, Lagrangian method has its advantage, but in the simulation of 3D solid, Euler method is the method that often adopts, and this formula is more difficult realization in three-dimensional flow field; (2) the sin θ in formula has adaptability preferably to specific two-dimensional case, needs change under three-dimensional situation; (3) this formula is not considered impact surface to the impact of air-flow, and calculating erosion condition is to have error.

Summary of the invention

Technical matters to be solved by this invention be to provide a kind of can be take the flow field situation on Euler method calculating aircraft surface and mass particle distribution situation as the basis, calculate aircraft surface erosion amount balancing method under the supersonic speed of aircraft surface erosion amount situation.

The present invention adopts following technical proposal to solve the problems of the technologies described above: the present invention has designed aircraft surface erosion amount balancing method under a kind of supersonic speed, comprises following concrete steps:

Step (1): gather under the supersonic speed state speed V of aircraft surface incoming flow air, density p and environmental pressure p;

Step (2): according to quality and principle of conservation of momentum, utilize the speed V of the incoming flow air that collects in step (1), density p, equal equation when environmental pressure p lists following three-dimensional Reynolds, and find the solution the air velocity u that obtains any position in the space:

$\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{\∂p}{{\∂x}_i}+\frac{\∂}{{\∂x}_j}(\mathrm{\μ}\frac{\∂u_i}{\∂x_j}-\mathrm{\ρ}\stackrel{\‾}{u_i^{\′}{u}_j^{\′}})$

Wherein, subscript i, the j value is 1,2,3, represents respectively x, y, z direction, i.e. u ₁Expression air component velocity in the x-direction, u ₂Expression air component velocity in the y-direction, u ₃Expression air component velocity in the z-direction; x ₁Parameter on expression x axle, x ₂Parameter on expression y axle, x ₃Parameter on expression z axle; μ represents the kinetic viscosity of air;

The derivative of expression to the time, The derivative of expression to x,

The derivative of expression to y, The derivative of expression to z;

Step (3): take step (2) air velocity result of calculation as the basis, particle is listed following quality and the equation of momentum, and find the solution and obtain particle fraction by volume α and particle rapidity u _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)=\mathrm{K\α\ρ}(u_a-u_s)+\mathrm{\α\ρF}$

Wherein, ρ _sBe defined as solids density, u _sBe defined as particle rapidity, F is defined as acceleration of gravity, and K is defined as inertial factor, and α is defined as particle fraction by volume, u _αBe defined as air velocity;

Step (4): according to the particle fraction by volume α that calculates in step (3) and particle rapidity result, utilize following formula to calculate the erosional retreat rate:

Se=kE _knV/(ρ _tC _N)

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

Wherein, S _eBe defined as the erosional retreat rate, E _knNormal direction momentum when being defined as the particle hits surface, ρ _tBe defined as the density of clashing into material, C _NBe defined as the material corrosion coefficient, k is defined as with erosion surface profile relevant parameters, α _nBe defined as in the calculating of solids track the particle fraction by volume of impact surface, V _nNormal velocity when being defined as the particle hits surface.

The present invention compared with prior art has following advantage:

Under the designed supersonic speed of this method, aircraft surface erosion amount balancing method can utilize Euler method calculating aircraft surface particle mass distribution situation, conveniently is generalized to three-dimensional situation, and the three-dimensional aircraft erosion amount that calculates has accuracy preferably.

Description of drawings

Fig. 1 is erosion amount analog computation comparison diagram as a result.

Embodiment

Below in conjunction with accompanying drawing, the present invention is done further and illustrates:

When the Calculation of Three Dimensional erosion condition, utilizing Lagrangian method to calculate has more inconvenience, and therefore, main computing method are because Euler method, and the formula that calculates for three-dimensional Euler method erosion amount is at present also never seen, therefore is necessary to propose a kind of new erosion amount computing method.

The present invention has designed aircraft surface erosion amount balancing method under a kind of supersonic speed, comprises following concrete steps:

Step (1): gather under the supersonic speed state speed V of aircraft surface incoming flow air, density p and environmental pressure p;

Step (2): according to quality and principle of conservation of momentum, utilize the speed V of the incoming flow air that collects in step (1), density p, equal equation when environmental pressure p lists following three-dimensional Reynolds, and find the solution the air velocity u that obtains any position in the space:

$\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{\∂p}{{\∂x}_i}+\frac{\∂}{{\∂x}_j}(\mathrm{\μ}\frac{\∂u_i}{\∂x_j}-\mathrm{\ρ}\stackrel{\‾}{u_i^{\′}{u}_j^{\′}})$

Wherein, subscript i, the j value is 1,2,3, represents respectively x, y, z direction, i.e. u ₁Expression air component velocity in the x-direction, u ₂Expression air component velocity in the y-direction, u ₃Expression air component velocity in the z-direction; x ₁Parameter on expression x axle, x ₂Parameter on expression y axle, x ₃Parameter on expression z axle; μ represents the kinetic viscosity of air;

The derivative of expression to the time,

The derivative of expression to x,

The derivative of expression to y,

The derivative of expression to z;

Set up three-dimensional non-steady Navier-Stokes equation according to quality, momentum, the large Conservation Relationship of energy three, write as following integrated form for a certain control volume Ω:

Wherein Q is the conservation variable, and FC is that Fv is the viscosity flux without sticking flux, and n controls the in vitro method vector, and S is the control volume border:

$W=\left[\begin{array}{c}\mathrm{\ρ}\\ \mathrm{\ρu}\\ \mathrm{\ρv}\\ \mathrm{\ρw}\\ \mathrm{\ρE}\end{array}\right],$ $F_c=\left[\begin{array}{c}\mathrm{\ρV}\\ \mathrm{\ρuV}+n_{x}P\\ \mathrm{\ρvV}+n_{y}P\\ \mathrm{\ρwV}+n_{z}P\\ \mathrm{\ρHV}\end{array}\right],$ $F_v=\left[\begin{array}{c}0\\ n_x{\mathrm{\τ}}_{\mathrm{xx}}+n_y{\mathrm{\τ}}_{\mathrm{xy}}+n_z{\mathrm{\τ}}_{\mathrm{xz}}\\ n_x{\mathrm{\τ}}_{\mathrm{yx}}+n_y{\mathrm{\τ}}_{\mathrm{yy}}+n_z{\mathrm{\τ}}_{\mathrm{yz}}\\ n_x{\mathrm{\τ}}_{\mathrm{zx}}+n_y{\mathrm{\τ}}_{\mathrm{zy}}+n_z{\mathrm{\τ}}_{\mathrm{zz}}\\ n_x{\mathrm{\Θ}}_x+n_y{\mathrm{\Θ}}_y+n_z{\mathrm{\Θ}}_z\end{array}\right],$ $Q=\left[\begin{array}{c}0\\ \mathrm{\ρ}{f}_{e,x}\\ \mathrm{\ρ}{f}_{e,y}\\ \mathrm{\ρ}{f}_{e,z}\\ \mathrm{\ρ}{\stackrel{\→}{f}}_e\·\stackrel{\→}{v}+{\stackrel{\·}{q}}_h\end{array}\right]$

In the viscosity flux, the expression formula of stress and heat flux is:

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

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

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

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

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

In order to make the sealing of N-S system of equations, also need replenish some physical relation formulas.

For perfect gas, state equation is arranged:

P=ρRT

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

If space derivative qx, qy, qz are found the solution by following system of equations:

$\left[\begin{array}{ccc}{x}_{c2}-x_{c1}& y_{c2}-y_{c1}& z_{c2}-z_{c1}\\ 0.5(x_{n2}+x_{n3})-x_{n1}& 0.5(y_{n2}+y_{n3})-y_{n1}& 0.5(z_{n2}+z_{n3})-z_{n1}\\ 0.5(x_{n3}+x_{n1})-x_{n2}& 0.5(y_{n3}+y_{n1})-y_{n2}& 0.5(z_{n3}+z_{n1})-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(q_{n2}+q_{n3})-q_{n1}\\ 0.5(q_{n3}+q_{n1})-q_{n2}\end{array}\right]$

In above-mentioned system of equations, use the flow parameter at each grid node place, its solution formula is as follows:

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

Wherein: r _i=(x _{0, i}x _n) ²+ (y _{0, i}-y _n) ²+ (z _{0, i}-z _n) ²] ^1/2

q _{0,, i}Be the flow parameter at the grid lattice heart place at node n place, ri is that a n is to the distance of grid lattice heart point.

Coefficient of viscosity is comprised of laminar flow coefficient of viscosity and coefficient of eddy viscosity, and determining of eddy viscosity is described as an example of RNG two Equation Turbulence Models example here.Suppose all to satisfy Boussinesq eddy viscosity relational expression between the rate of strain tensor in turbulent flow shear stress peace, eddy stress can be written as:

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

According to dimensional analysis, μ _tCan represent with k and ε:

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

C wherein _μBe empirical constant 0.09; K is turbulence pulsation kinetic energy; ε is the dissipative shock wave of turbulence pulsation kinetic energy.The transport equation of turbulence pulsation kinetic energy k is:

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

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

The N-S equation of having introduced the transport equation rear enclosed adopts central difference schemes on spatial spreading, adopt multistep Runge-Kutta form on time discrete.adopt famous Jameson central difference method when finding the solution, it has two characteristics: on space discrete, the flux on limit, internal field is to be averaged acquisition by the amount of flux with unit center place, left and right, stability for better capturing discontinuous point and raising form, added artificial dissipation's item when calculating flux, and artificial dissipation's item is comprised of second order and the Four order difference item of conservation variable, quality, momentum, the large fluid governing equation calculating of energy three is based on structured grid and carries out, obtain Stationary Solutions by multistep Runge-Kutta iteration in time, accelerating convergence and residual smooth technology have also been adopted in addition, make the CFL number can be increased to 6～8, speed of convergence improves 50%, convergence precision is controlled at the 10-8 magnitude, far field undisturbed boundary condition is set to 15 times of wing chord lengths, it is nonreflecting boundary condition.

In order to obtain reliable Flow Field Calculation parameter, the flow field solver of developing is mainly based on Euler or Navier-Stokes equation, can select S-A model, k-ε model etc. when finding the solution Viscous Flow, explicit form-Runge Kutta (Runge-Kutta) multistep form, ordinary differential equation about the time, can adopt the quadravalence Runge one Kutta method of multistep to find the solution, its form 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},$ α ₄=1。

For the explicit propelling form of Runge one Kutta method, for the requirement of stability, its time step is conditional.But because being utilizes non-permanent equation to find the solution the Steady Flow problem, time propelling and spatial spreading are independent of each other in solution procedure simultaneously, so Stationary Solutions is irrelevant with the time step that adopts, therefore can adopt the local time step-length on each grid cell, the flow field is pushed ahead with the time step near limit of stability everywhere, thereby accelerate the convergence process of whole iterative computation.

Time step on each grid cell can adopt following formula to calculate:

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

In following formula, Ω is the volume of grid, and u, v, w are speed component, SX, SY, SZ be area along the projection components of coordinate axis, c is local velocity of sound, CFL counts value and is decided by concrete mobility status.The advantage of explicit form is every propelling one time step, and calculated amount and memory space are all smaller, and program is simple.Shortcoming is that time step is limited by stability condition, calculates number too small, and efficient is lower, restrains required step number very large, causes overall computing time very long, therefore, and the characteristics that can use " small step is hurried up " to describe explicit form.

Step (3): take step (2) air velocity result of calculation as the basis, particle is listed following quality and the equation of momentum, and find the solution and obtain particle fraction by volume α and particle rapidity u _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)=\mathrm{K\α\ρ}(u_a-u_s)+\mathrm{\α\ρF}$

Wherein, ρ _sBe defined as solids density, u _sBe defined as particle rapidity, F is defined as acceleration of gravity, and K is defined as inertial factor, and α is defined as particle fraction by volume, u _αBe defined as air velocity;

Step (4): according to the particle fraction by volume α that calculates in step (3) and particle rapidity result, utilize following formula to calculate the erosional retreat rate:

S _e=kE _knV/(ρ _tC _N)

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

Wherein, S _eBe defined as the erosional retreat rate, E _knNormal direction momentum when being defined as the particle hits surface, ρ _tBe defined as the density of clashing into material, C _NBe defined as the material corrosion coefficient, k is defined as with erosion surface profile relevant parameters, α _nBe defined as in the calculating of solids track the particle fraction by volume of impact surface, V _nNormal velocity when being defined as the particle hits surface.

In order to verify the feasibility of this method, we adopt this method to calculate the erosion amount of a given configuration, and with contrast with the test findings under operating mode and budgetary equation simulation result, as shown in Figure 1, provided under same operating in figure, erosion amount result and the erosion amount result of former formula calculating and the comparison diagram of experimental result that the present invention adopts formula to calculate, from overall trend, the more former method result of calculation of result of calculation of the present invention is more near experiment value.

At how much leading edge point places (Y=0.00m place), experiment value is 0.02544m, and former formula calculated value is 0.0299m, the this patent calculated value is 0.02722m, can find out that this patent employing method is more smaller than former formula calculated value, more near experimental result, and it is good to coincide.

Claims (1)

1. aircraft surface erosion amount balancing method under a supersonic speed, is characterized in that, comprises following concrete steps:

Step (1): gather under the supersonic speed state speed V of aircraft surface incoming flow air, density p and environmental pressure p;

Step (2): according to quality and principle of conservation of momentum, utilize the speed V of the incoming flow air that collects in step (1), density p, equal equation when environmental pressure p lists following three-dimensional Reynolds, and find the solution the air velocity u that obtains any position in the space:

$\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{\∂p}{{\∂x}_i}+\frac{\∂}{{\∂x}_j}(\mathrm{\μ}\frac{\∂u_i}{\∂x_j}-\mathrm{\ρ}\stackrel{\‾}{u_i^{\′}{u}_j^{\′}})$

Wherein, subscript i, the j value is 1,2,3, represents respectively x, y, z direction, i.e. u ₁Expression air component velocity in the x-direction, u ₂Expression air component velocity in the y-direction, u ₃Expression air component velocity in the z-direction; x ₁Parameter on expression x axle, x ₂Parameter on expression y axle, x ₃Parameter on expression z axle; μ represents the kinetic viscosity of air;

The derivative of expression to the time,

The derivative of expression to x,

The derivative of expression to y,

The derivative of expression to z;

Step (3): take step (2) air velocity result of calculation as the basis, particle is listed following quality and the equation of momentum, and find the solution and obtain particle fraction by volume α and particle rapidity u _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)=\mathrm{K\α\ρ}(u_a-u_s)+\mathrm{\α\ρF}$

Wherein, ρ _sBe defined as solids density, u _sBe defined as particle rapidity, F is defined as acceleration of gravity, and K is defined as inertial factor, and α is defined as particle fraction by volume, u _αBe defined as air velocity;

Step (4): according to the particle fraction by volume α that calculates in step (3) and particle rapidity result, utilize following formula to calculate the erosional retreat rate:

S _e=kE _knV/(ρ _tC _N)

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

Wherein, S _eBe defined as the erosional retreat rate, E _knNormal direction momentum when being defined as the particle hits surface, ρ _tBe defined as the density of clashing into material, C _NBe defined as the material corrosion coefficient, k is defined as with erosion surface profile relevant parameters, α _nBe defined as in the calculating of solids track the particle fraction by volume of impact surface, V _nNormal velocity when being defined as the particle hits surface.

Images