Rapid calculation method for propagation process of shock waves and detonation waves in inhomogeneous medium

Abstract

The invention belongs to the field of computational explosion mechanics, and particularly relates to a rapid calculation method for a shock wave and detonation wave propagation process in a non-uniform medium, and the rapid calculation method comprises the following steps: 1, initializing a calculation program; step 2, estimating and calculating step length in each time step; step 3, carrying out integration on the ordinary differential equation by adopting an explicit three-order TVD Runge-Kutta format; 4, solving a Mach number through an A-M relation; 5, data correction is carried out, wherein the data correction comprises regularization and smoothing processing of grids; and step 6, checking whether termination time is reached or not, if so, ending all processes to complete operation, and if not, returning to the step 2 to continue to circularly solve. According to the rapid calculation method, rapid numerical calculation of propagation of shock waves and detonation waves in a non-uniform medium is achieved, wherein the phenomena comprise refraction, transmission, Mach reflection and the like.

Description

Technical Field

[0001] The invention belongs to the technical field of computational explosion mechanics, and in particular relates to a fast calculation method for the propagation process of shock wave detonation wave in a non-uniform medium. Background Art

[0002] Detonation refers to the process of violent combustion of combustible media. This complex process is widely present in nature and industry, and has extremely high destructive power and important application value. Chemical energy is rapidly released in the form of heat and mechanical energy, forming a supersonic shock wave. The detonation phenomenon was first studied by scientists in the 19th century, and important breakthroughs were made in experiments and theories. Scientists gradually established detonation and deflagration. The key feature of detonation is the self-sustaining combustion reaction guided by the high-pressure shock wave. Even in the absence of an ignition source, the pressure wave can provide the necessary energy to keep the reaction going.

[0003] The study of shock waves and detonation phenomena in inhomogeneous media has important scientific and engineering significance. In basic science, this type of research can deeply reveal the complex mechanisms in fluid mechanics, thermodynamics, and chemical reaction kinetics, and enrich our understanding of violent physical processes in nature. In the application field, understanding and controlling shock waves and detonation phenomena plays a key role in aerospace, energy engineering, military defense and other fields. For example, in the aerospace industry, shock wave management directly affects the performance and safety of aircraft; in national defense science and technology, studying detonation phenomena helps to improve the effectiveness and safety of weapon systems. Therefore, the systematic study of shock waves and detonation phenomena in inhomogeneous media can not only advance the development of scientific theories, but also promote technological advancement and innovation in multiple high-tech industries.

[0004] However, the prior art lacks a calculation method for the propagation of shock waves and detonation waves in inhomogeneous media. Summary of the invention

[0005] In order to improve the efficiency of studying shock waves and detonation waves in inhomogeneous media and enrich research methods, the present invention provides a fast calculation method for the propagation process of shock waves and detonation waves in inhomogeneous media. The calculation method can accurately calculate the propagation deformation of the wave surface, so as to efficiently and accurately calculate the shock wave and detonation wave problems in different inhomogeneous media. The method is suitable for the field of explosion and impact, and can handle the problems of shock waves and detonation waves in inhomogeneous media with arbitrary boundary shapes, such as: the propagation of shock waves and detonation waves, the refraction of shock waves and detonation waves, etc.

[0006] In order to achieve the above object, the present invention adopts the following specific technical solutions:

[0007] A fast calculation method for the propagation process of shock wave detonation wave in a non-uniform medium, the fast calculation method comprising the following steps:

[0008] Step 1, initializing the calculation program;

[0009] Step 2, estimate the calculation step size in each time step;

[0010] Step 3, integrating the ordinary differential equation using an explicit third-order TVD Runge-Kutta format;

[0011] Step 4, solving the Mach number through the AM relationship;

[0012] Step 5, data correction, including grid regularization and smoothing;

[0013] Step 6, check whether the end time has been reached. If it has been reached, end all processes and complete the calculation. If it has not been reached, return to step 2 and continue the loop to solve.

[0014] Furthermore, step 1 specifically includes:

[0015] Step 1.1, discretize the initial discontinuous wave peak. At the initial moment, the wave surface is expressed in Cartesian coordinates as:

[0016] x(α0,β)=(x(α0,β),y(α0,β)) T (1);

[0017] Among them, α0 is the initial coordinate of the grid point on the αβ plane, β is the coordinate value of the grid point in the β direction, x is the horizontal coordinate of the discrete point on the xy physical plane, and y is the vertical coordinate of the discrete point on the xy physical plane;

[0018] Step 1.2, introduce the abscissa of the discrete curve along the wavefront:

[0019]

[0020] Among them, for have And Δs i =s i -s i-1 ,i=2,...,N,s i is the abscissa of the discrete curve of grid point i on the wave surface, N is the total number of grid points, and Δs is the difference in spacing between adjacent grid points;

[0021] Step 1.3, homogenize the discretization of the initial wavefront:

[0022]

[0023] Among them, s i (α0) is the s-coordinate of the i-th point on the initial wave surface.

[0024] Furthermore, step 2 specifically includes:

[0025] The time step is calculated based on the CFL condition and the Ray line non-intersection principle, where the CFL condition is:

[0026]

[0027] In the above formula, M i is the Mach number of the current time step, λ and k are the Mach number M i The function is in the form of:

[0028]

[0029] In the above formula, γ is the gas specific heat ratio, and μ is the post-shock Mach number:

[0030]

[0031] The ray line non-intersection condition is:

[0032]

[0033] in, is the wave surface normal vector at the i-th point on the wave surface, l i-1 and k i There are two abbreviation operators:

[0034]

[0035] Furthermore, step 3 integrates the ordinary differential equation using an explicit third-order TVD Runge-Kutta scheme:

[0036]

[0037] Furthermore, step 4 is: use the bisection method to solve M i (α+Δα), the goal is the upper limit of the integral of the following formula:

[0038]

[0039] Furthermore, step 5 specifically includes:

[0040] Step 5.1, use the cubic spline interpolation method to interpolate the wave surface and give a unit normal vector;

[0041]

[0042] In the above formula, Y'(s i ) is the ordinate of the i-th point, X'(s i ) is the horizontal coordinate of the i-th point;

[0043] Step 5.2, insert points in the expansion area and delete points in the compression area to control the grid spacing to meet:

[0044] d min Δs≤Δs i ≤d max Δs (9).

[0045] Furthermore, step 6 specifically includes:

[0046] Step 6.1, output the result of the current step when the output condition is met; if the current number of iterations and time meet the termination condition, the calculation will stop;

[0047] Step 6.2: If the end condition is not met, return to step 2 and start the calculation of the next time step.

[0048] Compared with the prior art, the technical solution of the present invention has the following beneficial effects:

[0049] 1. The existing technology cannot calculate the propagation of shock waves and detonation waves in inhomogeneous media. The present invention theoretically models inhomogeneous media and realizes the calculation of shock and detonation wave propagation in inhomogeneous media. The present invention uses monotone cubic spline interpolation to reconstruct the wavefront and uses high-precision third-order TVD Runge-Kutta format discrete differential equations to realize rapid numerical calculation of shock and detonation wave propagation in inhomogeneous media.

[0050] 2. The present invention adjusts the discrete form of the AM relation, introduces a non-uniform gradient source term into the traditional AM relation, realizes the numerical description of the non-uniform medium, enables it to be integrated step by step, but still maintains the conservation characteristics of the scheme. In the traditional geometric shock wave dynamics relation, a source term representing the non-uniformity is added, which realizes the rapid calculation of the propagation of shock waves and detonation waves in non-uniform media.

[0051] 3. The present invention utilizes the Lagrangian motion of waves and the decomposition characteristics in the ray tube, and manages the number of points on the wavefront through a dynamic process called regularization. It can effectively solve the grid distortion problem when the wavefront is greatly deformed, and can handle the propagation of shock waves and detonation waves in inhomogeneous media under arbitrary boundaries.

[0052] 4. The calculation method of the present invention is intended to provide a tool for simulation research of aerodynamic gas phase explosion problems, and is suitable for scientific researchers, industrial and civil gas explosion protection designers, as well as rescue personnel and software verifiers who handle and evaluate gas explosion accidents. BRIEF DESCRIPTION OF THE DRAWINGS

[0053] Figure 1 is a flow chart of the fast calculation method of the present invention;

[0054] Figure 2 It is a schematic diagram of the wave propagation along the Ray tube;

[0055] Figure 3 The displacement of the wave surface grid points in a single time step and the schematic diagram of the Ray tube;

[0056] Figure 4 Schematic diagram of the process of adding and deleting grid points;

[0057] Figure 5 Schematic diagram of vertical wave surface condition at the boundary;

[0058] Figure 6 This is the Mach number cloud diagram of a detonation wave with a Mach number of 4.5 passing through a non-uniform interface with a sound speed ratio of 1:2 and an inclination angle of 30°;

[0059] Figure 7 It is the wave surface grid point distribution and grid point normal vector result diagram;

[0060] Figure 8 This is the Mach number cloud diagram of a shock wave with a Mach number of 3 passing through a non-uniform interface with a sound speed ratio of 1:2 and an inclination angle of 30°;

[0061] Fig. 9 A grid diagram with new grid points added after wavefront expansion. DETAILED DESCRIPTION

[0062] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0063] The fast calculation method of the present invention theoretically models the inhomogeneous medium and realizes the calculation of the propagation problem of shock wave and detonation wave in the inhomogeneous medium. The fast calculation method of the present invention is based on Whitham's Geometrical ShockDynamics theory. In order to describe the thermodynamic state changes before and after the shock wave and the detonation wave, the software uses monotone cubic spline interpolation to reconstruct the wavefront, and uses a high-precision third-order TVD Runge-Kutta format to discretize the differential equation. In order to construct the inhomogeneous conditions in the medium, the idea of ​​the present invention is to add an inhomogeneous gradient source term to the basic AM relationship, and use the bisection method to solve the AM relationship containing the source term to achieve the purpose of calculating the propagation of shock waves and detonation waves in the inhomogeneous medium. The present invention realizes the rapid numerical calculation of the propagation of shock waves and detonation waves in inhomogeneous media, including phenomena such as refraction, transmission, and Mach reflection.

[0064] The main ideas and contents of the present invention are as follows:

[0065] 1. Consider the shock wave as a constantly changing curve or surface, and use a few parameters such as the shock wave's Mach number and velocity vector to describe its geometry and evolution. The schematic diagram is as follows: Figure 2 shown.

[0066] 2. Based on the laws of conservation of mass, momentum and energy, a set of partial differential equations describing the evolution of shock waves is derived, called the Whitham equations. This set of equations describes the change of the local curvature and normal velocity of the shock wave with time.

[0067] 3. The Whitham equations are hyperbolic and can be solved using the characteristic line method. The characteristic lines represent the information propagation path on the shock wave surface.

[0068] 4. By numerically solving the AM relationship, the evolution process of the shock wave surface can be simulated and the position and shape of the shock wave surface can be predicted.

[0069] 5. Introduce the concept of shock ray, which is the trajectory of the normal line on the shock wave surface.

[0070] The present invention has high calculation efficiency and can quickly simulate the evolution of the shock wave surface without the need to calculate the flow state behind the shock wave in detail.

[0071] The core formula of the present invention is the area-Mach number relationship AM relationship:

[0072]

[0073] in:

[0074]

[0075] In the above formula, M is the wavefront Mach number, and the reaction heat release To set the initial incident wave to be a shock wave or a detonation wave.

[0076] After each iterative calculation, the normal direction and Mach number of each point on the wave surface will be calculated, and the position x of each point on the shock wave surface and detonation wave surface at the next moment is obtained by the following ordinary differential equation:

[0077]

[0078] In the above formula, t is the arrival time of the impact, n is the unit normal vector of the wave surface, c0 is the medium sound velocity, and U is the propulsion velocity of the discrete point. The schematic diagram is as follows: Figure 3 shown.

[0079] The idea of ​​setting the inhomogeneous zone in the present invention is innovative. The implementation method of the inhomogeneous zone is discussed as follows. It is necessary to consider the propagation of shock waves in inhomogeneous media. In this case, the shock wave bends and twists as it adapts to the inhomogeneity of the fluid and is affected by the curvature of the shock wave and the existence of the wall. The only additional feature of the inhomogeneous medium is the numerical treatment of the differential AM relationship. The traditional AM relationship is written as follows by adding the source term:

[0080]

[0081] Among them, M is the wavefront Mach number, A is the cross-sectional area of ​​the flow tube, and S is a source term that includes the gradient of the medium inhomogeneity. In the expression of the source term, a0 represents the medium sound velocity, p0 represents the medium initial pressure, and γ represents the specific heat ratio of the gas. Specifically:

[0082]

[0083] Among them, the coefficients of the sound velocity gradient term g(M,x), the pressure gradient term h(M,x) and the specific heat gradient term k(M,x) are as follows:

[0084]

[0085] In each iteration, the AM relationship is solved and the Mach number at the next moment is solved by integration:

[0086]

[0087] The key to solving the above integral is the integral of the non-uniform gradient. This integral is calculated using the following approximate polynomial summation formula:

[0088]

[0089] Embodiment 1

[0090] This embodiment provides a fast calculation method for the shock wave and detonation wave propagation process in a non-uniform medium. Figure 1 As shown, the following steps are included:

[0091] Step 1, initialize the calculation program; the initialization calculation program is the process of building a simulation calculation model, and a simplified simulation calculation model is built for the specific problem to be analyzed; the objects of initialization include: the shape of the boundary of the calculation domain, the type of initial discontinuity (shock wave or detonation wave), the position of the initial discontinuity, the intensity (speed) of the initial discontinuity, the shape of the inhomogeneous region, the medium properties of the inhomogeneous region, and the calculation control parameters such as the number of initial wave surface grid points, output interval, end time, and smoothing frequency; the specific process is:

[0092] Step 1.1, discretize the initial discontinuous wave peak. At the initial moment, the wave surface is expressed in Cartesian coordinates as:

[0093] x(α0,β)=(x(α0,β),y(α0,β)) T (1);

[0094] Among them, α0 is the initial coordinate of the grid point on the αβ plane, β is the coordinate value of the grid point in the β direction, x is the horizontal coordinate of the discrete point on the xy physical plane, and y is the vertical coordinate of the discrete point on the xy physical plane;

[0095] Step 1.2, introduce the abscissa of the discrete curve along the wavefront:

[0096]

[0097] Among them, for have And Δs i =s i -s i-1 ,i=2,...,N,s i is the abscissa of the discrete curve of grid point i on the wave surface, N is the total number of grid points, and Δs is the difference in spacing between adjacent grid points;

[0098] Step 1.3, homogenize the discretization of the initial wavefront:

[0099]

[0100] Among them, s i (α0) is the s-coordinate of the i-th point on the initial wave surface.

[0101] Step 2: Estimate the calculation step length in each time step, that is, at the beginning of each time step, recalculate the time step length in this time step according to the CFL condition of computational mechanics and the principle that Ray lines do not intersect. The specific process is:

[0102] Step 2.1, calculate the threshold of the time step Δα according to the CFL condition:

[0103]

[0104] Where M i is the Mach number of the current time step, λ and k are the Mach number M i The function is in the form of:

[0105]

[0106] In the above formula, γ is the gas specific heat ratio, and μ is the post-shock Mach number:

[0107]

[0108] Step 2.2, calculate the time step threshold according to the Ray line non-intersection principle:

[0109]

[0110] Among them, l i-1 and k i They are two abbreviation operators respectively.

[0111] Step 2.3, determine the final time step based on the above two conditions:

[0112] Δα=δmin(Δα s ,Δα c );

[0113] Among them, δ is the control factor, and its value is 0.5.

[0114] Step 3: Integrate the ordinary differential equation using the explicit third-order TVD Runge-Kutta format and iterate the solution of the ordinary differential equation. The main purpose of this step is to solve the integral of the AM relation and the Mach number M at the next time step. i (α+Δα) is obtained by integrating the following formula along the flow tube:

[0115]

[0116] Step 3.1, determine whether the flow tube is deleted or a new flow tube is added in this time step, as shown in the schematic diagram Figure 4 shown.

[0117] Step 3.2, if the grid points are deleted, the third-order Runge-Kutta method is used:

[0118]

[0119] Among them, n i The normal vector direction of the discrete point.

[0120] For the integral starting from α, the integrand formula is as follows:

[0121]

[0122] The letters with wavy lines on top represent the variable values ​​before adding and deleting.

[0123] Step 3.3, if the grid points are increased, the third-order Runge-Kutta method is used to integrate the following integrand:

[0124]

[0125] Step 3.4, use the monotone cubic spline interpolation method to interpolate the shock wave surface to avoid generating new extreme values. Calculate the wavefront unit normal vector accordingly:

[0126]

[0127] Step 3.5, calculate the grid information at the boundary according to the vertical condition of the wave surface at the boundary, as shown in the schematic diagram Figure 5 As shown. The spline interpolation formula at the boundary is:

[0128] X'(s bound )=N x ; Y'(s bound )N y ;

[0129] Among them, N x and N y are the tangent vector components of the boundary point curve respectively.

[0130] The normal vector of the wave surface grid point at the boundary is defined as:

[0131] n(s bound )=(N y' -N x ) T ;

[0132] The second derivative of the spline interpolant should be set to 0 at the boundaries:

[0133] X"(s bouond )=0;Y"(s bouond )=0;

[0134] Step 4: Use Newton's iterative root-finding method to solve the Mach number M at the next moment. i (α+Δα).

[0135] Step 4.1, choose an initial guess value M i (α+Δα);

[0136] Step 4.2, calculate the AM relationship to determine whether it is positive or negative;

[0137] Step 4.3, if the positive and negative conditions are not met, continue to iterate;

[0138] Step 4.4, repeat the above process until the guessed value is within the expected accuracy range.

[0139] Step 5, perform data correction, i.e. wave surface shaping correction, including grid regularization and smoothing;

[0140] Step 5.1, determine the spacing Δs between adjacent grid points on the wave surface iWhether it is too far or too close, the judgment criteria are as follows:

[0141] dminΔs≤Δs i ≤dmaxΔs;

[0142] Step 5.2, add or delete grid points as needed.

[0143] Step 5.3, use a smoothing procedure to suppress the high-frequency errors of the wave surface. After multiple iterations (between 10-50 in practical applications), let i = 2, ..., N-1:

[0144]

[0145] Among them, hr and hl are two abbreviation operators, as follows:

[0146]

[0147] Smooth the even points first, and then smooth the odd points using the above formula.

[0148] Step 6: Based on the set termination conditions, output the corresponding simulation results:

[0149] Step 6.1: If the current calculation step meets the output setting conditions, the calculation data will be output and the serial number of the output file will be recorded; the variable information contained in the grid point will be written into the file.

[0150] Step 6.2: If the end condition is not met, return to step 2 and start a new time step; if the end condition is met, end the calculation.

[0151] The above six steps realize fast and accurate numerical calculation of shock wave and detonation wave propagation problems in inhomogeneous media.

[0152] Embodiment 2

[0153] Taking the propagation of shock waves and detonation waves in inhomogeneous media with different sound speeds as an example, the propagation process of shock waves and detonation waves in inhomogeneous media is numerically simulated. The results show that shock waves and detonation waves have complex wavefront structures at inhomogeneous interfaces, including three wave points, oblique shock waves, sparse waves and other complex structures, indicating that the wavefront of the method of the present invention has high tracking accuracy and computing power; and can output the Mach number and wavefront normal vector at any point on the wavefront. In scientific research, the wavefront Mach number and normal vector can be used as the basis for identifying the wavefront structure to analyze the wavefront structure. The following takes the propagation process of shock waves and detonation waves refracted at an inclined inhomogeneous interface as an example. Figure 1 As shown, a fast calculation method for the shock wave and detonation wave propagation process in a non-uniform medium of the present invention is described in detail:

[0154] Step 1: construct a simplified simulation model for the specific problem to be analyzed. The process of constructing a simplified simulation model is called the model initialization process. The objects of initialization include: the shape of the boundary of the calculation domain, the type of initial discontinuity (shock wave or detonation wave), the position of the initial discontinuity, the intensity (speed) of the initial discontinuity, the shape of the inhomogeneous region, the medium properties of the inhomogeneous region, and the number of initial wave surface grid points, output interval, end time, smoothing frequency and other calculation control parameters; Step 1 specifically includes the following steps:

[0155] Step 1.1, set the calculation area to [0,3]×[0,1], and set the initial grid point spacing of the wave surface to 1.0 / 100=0.001;

[0156] Step 1.2, set the boundary to wall condition, the shape to a straight rectangular pipe, set the gas model to an ideal gas model, and the gas specific heat ratio to 1.4;

[0157] Step 1.3, the type of initial discontinuity is set to shock wave (Mach number 3) and detonation wave (Mach number 4.5) according to the heat release, and the discontinuity position is set to 1.0;

[0158] Step 1.4, the inhomogeneous region is set as a 30° inclined interface, and the sound velocity ratio of the front and rear media is set to 1:2;

[0159] Step 1.5, set the spacing threshold coefficient for adding and deleting grids, adding points when 1.5 times the initial spacing, and deleting points when 0.5 times the initial spacing;

[0160] Step 1.6, set the initial calculation control parameters, set the end time to 5, the step interval for saving the results to 30, and the smoothing frequency to perform interpolation smoothing every 50 iterations.

[0161] Step 2: A new time step starts. At the beginning of each time step, the time step length within the current time step is recalculated according to the CFL condition of computational mechanics and the Ray line non-intersection principle. Step 2 specifically includes the following steps:

[0162] Step 2.1, calculate the time step threshold according to the CFL condition.

[0163] Step 2.2: Calculate the time step threshold based on the Ray line non-intersection principle.

[0164] Step 2.3, determine the final time step based on the above two conditions.

[0165] Step 3, iteratively solve the ordinary differential equation. Solving the integral of the AM relation is a preparation for the next step of solving the Mach number; Step 3 specifically includes the following steps:

[0166] Step 3.1, determine whether the flow tube is deleted or a new flow tube is added in this time step.

[0167] Step 3.2, if the grid points are deleted, use the third-order Runge-Kutta method to integrate the following formula:

[0168]

[0169] The letters with wavy lines on top represent the variable values ​​before adding and deleting.

[0170] Step 3.3, if the grid points are increased, use the third-order Runge-Kutta method to integrate the following formula:

[0171]

[0172] Step 3.4, use the monotone cubic spline interpolation method to interpolate the shock wave surface to avoid generating new extreme values, and calculate the unit normal vector of the wavefront.

[0173] Step 3.5, calculate the grid information at the boundary according to the vertical condition of the wave surface at the boundary.

[0174] Step 4, traverse each grid point and use Newton's iterative root-finding method to solve the Mach number at the next moment. Step 4 specifically includes the following steps:

[0175] Step 4.1, choose an initial guess value.

[0176] Step 4.2, calculate the AM relationship to determine whether it is positive or negative.

[0177] Step 4.3, if the positive and negative conditions are not met, continue iterating.

[0178] Step 4.4, repeat the above process until the guessed value is within the expected accuracy range.

[0179] Step 5, wavefront shaping correction, specifically includes the following steps:

[0180] Step 5.1, determine whether the distance between adjacent grid points on the wave surface is too far or too close.

[0181] Step 5.2, add or delete grid points as needed.

[0182] In step 5.3, a smoothing procedure is used to suppress high-frequency errors in the wave surface.

[0183] Step 6, based on the set termination condition, output the corresponding simulation result, specifically including the following steps:

[0184] Step 6.1: If the current calculation step meets the output setting conditions, the calculation data will be output and the serial number of the output file will be recorded; the variable information contained in the grid point will be written into the file.

[0185] Step 6.2: If the end condition is not met, return to step 2 and start a new time step; if the end condition is met, end the calculation.

[0186] exist Figure 6 The flow field cloud diagram of a detonation wave with a Mach number of 4.5 passing through a non-uniform interface is shown in Figure 1. The global specific heat ratio of the flow field is set to a fixed value of 1.4. The inclination angle of the non-uniform interface is 30°, and the ratio of the sound speed of the incident and transmitted media is 1:2. When the detonation wave contacts the interface and transmits, the wave surface bends. In the transmitting medium, the detonation wave accelerates forward significantly, and an obvious Mach number discontinuity appears at the interface, and a three-wave point "shock-shock" appears. Figure 7 The distribution and movement of grid points on the wave surface, as well as the forward direction vector of each grid point, are shown. Figure 7 The Mach number distribution on the wave front is also shown. Figure 7 In the expansion region of the wave surface, the wave surface particles move away from each other, and near the triple wave point, the wave surface particles approach each other.

[0187] exist Figure 8 The flow field cloud map of a shock wave with a Mach number of 3 passing through an inhomogeneous interface is shown in the figure. The global specific heat ratio of the flow field is set to a fixed value of 1.4. The inclination angle of the inhomogeneous interface is 30°, and the ratio of the sound speed of the incident and transmitted media is 1:2. When the shock wave contacts the interface for transmission, the wavefront bends significantly. Compared with the detonation wave, the shock wave has a slower speed and weaker intensity, and the wavefront bending is more severe. The shock wave in the transmitting medium is convex, and an obvious Mach number discontinuity appears at the interface. Three wave points "shock-shock" appear above the interface. The oblique shock wave is significantly longer than the oblique shock wave in the detonation wave example, and the inclination angle of the three wave points is larger. More drastic wavefront deformation places higher requirements on the wavefront particle grid, and the grid spacing needs to be adjusted in time at the drastic changes. This method solves the problem of adaptively adjusting the grid spacing well, such as Fig. 9 shown.

[0188] It can be seen from the above embodiments that the fast calculation method of the present invention can well track the propagation process of the discontinuity surface in the inhomogeneous medium, and can record the Mach number and direction vector of each point on the wave surface, effectively realizing the numerical calculation of the shock wave and detonation wave propagation problems in the inhomogeneous medium, and providing a new and efficient calculation method for the research in the related gas dynamics field.

[0189] Obviously, those skilled in the art can make various changes and modifications to the embodiments of the present invention without departing from the spirit and scope of the present invention. Thus, if these modifications and variations of the present invention fall within the scope of the claims of the present invention and their equivalents, the present invention is also intended to include these modifications and variations.

Claims

1. A fast calculation method for the propagation process of shock wave and detonation wave in inhomogeneous medium, characterized in that: The following steps are involved: Step 1, initializing the calculation program; Step 2, estimate the calculation step size in each time step; Step 3, integrating the ordinary differential equation using the explicit third-order TVD Runge-Kutta format; Step 4, solving the Mach number through the AM relationship; Step 5, data correction, including grid regularization and smoothing; Step 6, check whether the end time has been reached. If it has been reached, end all processes and complete the calculation. If it has not been reached, return to step 2 and continue the loop to solve.

2. The fast calculation method according to claim 1, characterized in that: Step 1 specifically includes: Step 1.1, discretize the initial discontinuous wave peak. At the initial moment, the wave surface is expressed in Cartesian coordinates as: x(α0,β)=(x(α0,β),y(α0,β)) T (1); Among them, α0 is the initial coordinate of the grid point on the αβ plane, β is the coordinate value of the grid point in the β direction, x is the horizontal coordinate of the discrete point on the xy physical plane, and y is the vertical coordinate of the discrete point on the xy physical plane; Step 1.2, introduce the abscissa of the discrete curve along the wavefront: Among them, for have And Δs i =s i -s i-1 ,i=2,...,N,s i is the abscissa of the discrete curve of grid point i on the wave surface, N is the total number of grid points, and Δs is the difference in spacing between adjacent grid points; Step 1.3, homogenize the discretization of the initial wavefront: Among them, s i (α0) is the s-coordinate of the i-th point on the initial wave surface.

3. The fast calculation method according to claim 2, characterized in that: Step 2 specifically includes: The time step is calculated based on the CFL condition and the Ray line non-intersection principle, where the CFL condition is: In the above formula, M i is the Mach number of the current time step, λ and k are the Mach number M i The function is in the form of: In the above formula, γ is the gas specific heat ratio, and μ is the post-shock Mach number: The ray line non-intersection condition is: in, is the wave surface normal vector at the i-th point on the wave surface, l i-1 and k i There are two abbreviation operators:

4. The fast calculation method according to claim 3, characterized in that: Step 3 Integrate the ordinary differential equation using the explicit third-order TVD Runge-Kutta scheme:

5. The fast calculation method according to claim 4, characterized in that: Step 4: Use the binary search method to solve M i (α+Δα), the goal is the upper limit of the integral of the following formula:

6. The fast calculation method according to claim 5, characterized in that: Step 5 specifically includes: Step 5.1, use the cubic spline interpolation method to interpolate the wave surface and give a unit normal vector; In the above formula, Y'(s i ) is the ordinate of the i-th point, X'(s i ) is the horizontal coordinate of the i-th point; Step 5.2, insert points in the expansion area and delete points in the compression area to control the grid spacing to meet: d min Δs≤Δs i ≤d max Δs (9)。 7. The fast calculation method according to claim 6, characterized in that: Step 6 specifically includes: Step 6.1, output the result of the current step when the output condition is met; if the current number of iterations and time meet the termination condition, the calculation will stop; Step 6.2: If the end condition is not met, return to step 2 and start the calculation of the next time step.

Images