A global discontinuity detection method for high-order finite difference schemes for both interior points and boundary points

Abstract

This invention proposes a global discontinuity detection method applicable to both interior points and boundary high-order finite difference schemes. The method includes: setting a computational domain; constructing a semi-uniformly distributed grid based on the computational domain; using the grid cells divided by the solution points as the cells to be detected, further dividing the cells into interior cells and boundary cells, and performing discontinuity detection on each type of cell; constructing polynomial functions based on left-facing and right-facing discontinuity detection templates; obtaining the simplified form of the first derivative of the polynomial functions; calculating the overdetermined extrema on the cells to be detected based on the simplified first derivative; and determining whether the cells to be detected are discontinuous based on the sign of the overdetermined extrema. This invention overcomes the shortcomings of existing discontinuity detection methods in effectively handling non-uniformly distributed grid points at the boundaries and missing template points at the boundaries, achieving accurate and robust partitioning of discontinuous and smooth regions throughout the entire computational domain.

Description

Technical Field

[0001] This invention relates to the field of engineering fluid simulation, and in particular to a global discontinuity detection method applicable to both interior points and boundary high-order finite difference schemes. Background Technology

[0002] In the design and performance evaluation of aviation, aerospace, and marine equipment, high-precision engineering fluid simulation technology has become a core tool supporting design decisions, shortening R&D cycles, and reducing testing risks. This technology provides reliable data for key engineering problems such as aircraft aerodynamic shape optimization, propulsion system performance analysis, and structural thermal load prediction through high-fidelity numerical simulation of complex flow phenomena (such as supersonic / hypersonic flow, shock wave interference, and separated flow). With the continuous improvement of performance requirements for aircraft and other equipment, and the increasing complexity of the geometry of design objects, the demand for accuracy, resolution, and efficiency in engineering fluid simulation is also growing.

[0003] The finite difference method (FDM) is a widely used numerical discretization technique in engineering fluid simulation, offering advantages such as high accuracy, high computational efficiency, and ease of implementation. However, when simulating complex flow phenomena such as supersonic flows over aircraft, the flow field often contains various discontinuous structures, including shock waves, contact discontinuities, and shear layers. This poses a significant challenge to the discontinuity capture capabilities of the numerical schemes used in engineering fluid simulation. Traditional high-precision FDM shock wave capture methods use nonlinear schemes, which, while capable of stably solving strongly discontinuous problems, introduce additional dissipation in smooth regions of the flow field, reducing the resolution of flow details. Furthermore, their complex nonlinear mechanisms significantly increase computational costs in smooth regions.

[0004] To balance high accuracy, high resolution, and high computational efficiency, an effective strategy is to construct a hybrid scheme. This involves using a nonlinear scheme with shock wave capture capabilities in discontinuous regions such as shock waves, while employing a linear scheme with high resolution and high computational efficiency in smooth regions. A prerequisite for constructing an efficient and reliable hybrid scheme is the development of robust and effective discontinuity detection methods. The role of discontinuity detection is to accurately identify and distinguish discontinuous regions such as shock waves and contact discontinuities in the flow field, numerically differentiating them from the smooth parts of the flow field to facilitate subsequent calculations.

[0005] Existing discontinuity probing methods are primarily designed for the computational domain. When addressing the "missing template points" problem near the boundary caused by template extension beyond the computational domain, virtual mesh point technology is typically employed, allowing the same discontinuity probing algorithm to be used across the entire computational domain. However, for globally conserved high-precision finite difference schemes, the computational domain boundary not only suffers from missing template points, but the solution points at the boundary are also no longer uniformly distributed as within the computational domain. Therefore, interior discontinuity probing methods cannot be directly applied to such boundary regions. In practical applications, the complex shapes of engineering simulations (e.g., large passenger aircraft, spacecraft, automobiles, etc.) commonly rely on multi-mesh technology, resulting in a very high demand for simulating mesh boundaries. Failure to guarantee the simulation accuracy and reliability at the boundaries will severely impact the credibility of the overall engineering design. Summary of the Invention

[0006] Based on this, the present invention provides a global discontinuity detection method applicable to both interior points and boundary high-order finite difference schemes. This discontinuity detection method is based on the solution point grid distribution and performs discontinuity detection cell-by-cell. Its core lies in overcoming the shortcomings of existing discontinuity detection methods in effectively handling non-uniform boundary points and missing template points, thereby achieving accurate and robust partitioning of discontinuous and smooth regions throughout the entire computational domain, including the computational domain boundary. This provides an indispensable key technological foundation for overcoming existing bottlenecks and successfully constructing high-performance hybrid schemes.

[0007] This invention proposes a global discontinuity detection method applicable to both interior point and boundary high-order finite difference schemes, comprising:

[0008] S1. Based on the task to be calculated, construct the hyperbolic conservation law equations of the simulation object and set the computational domain;

[0009] S2. Based on the computational domain, construct a semi-uniform grid. The semi-uniform grid refers to dividing the computational domain into intersecting solution points and flux points in each coordinate direction, where flux points are evenly spaced and solution points are arranged at the center of adjacent flux points. Then, the boundary solution points are replaced with non-uniformly conserved points. The boundary solution points refer to the solution points closest to the boundary. The specific number is determined according to the accuracy of the format. For a global (2p-1) order format, where 2≤p≤6, (p-1) non-uniformly conserved points are arranged near the boundary.

[0010] S3. The grid cells divided according to the solution points are used as the cells to be detected. The cells to be detected are divided into cells inside the computational domain and cells at the boundary of the computational domain. Discontinuous detection is performed on each of the two types of cells. For each cell to be detected, a left-facing discontinuous detection template and a right-facing discontinuous detection template are selected. The discontinuous detection template refers to the engineering control unit used to determine whether the detection cell contains discontinuities. It is a set of several continuous solution points around the cell to be detected. The specific number is determined according to the precision of the format. For the global (2p-1) order format, the left-facing discontinuous detection template and the right-facing discontinuous detection template are both composed of (2p-1) solution points. The presence or absence of discontinuities on the detection cell is determined according to the known physical quantities on these (2p-1) solution points and the positional relationship between these (2p-1) solution points and the detection cell.

[0011] S4. Based on the left windward intermittent detection template and the right windward detection template, construct polynomial functions Q1(X) and Q2(X) respectively;

[0012] S5. Find the first derivative of the polynomial function to obtain the first derivative functions Q′1(X) and Q′2(X) of the polynomial function;

[0013] S6. Simplify the first derivative of the polynomial function to obtain the simplified form of the first derivative of the polynomial function. and The simplification process refers to converting the fractional coefficients in Q′1(X) and Q′2(X) into integer coefficients through multiplication;

[0014] S7. Based on the simplified form of the first derivative function... and Calculate the overdetermined maximum and minimum values ​​on the unit to be detected respectively;

[0015] S8. Determine whether the unit to be detected is discontinuous based on the sign of the over-determined maximum and minimum values.

[0016] In one embodiment, in step S2, a semi-uniform grid is constructed based on the computational domain. The semi-uniform grid refers to dividing the computational domain into staggered solution points and flux points along each coordinate direction, where flux points are evenly spaced, solution points are arranged at the centers of adjacent flux points, and then boundary solution points are replaced with non-uniformly conserved points. The boundary solution points refer to the several solution points closest to the boundary, the specific number of which is determined based on the precision of the format. For a global (2p-1) order format, where 2≤p≤6, (p-1) non-uniformly conserved points are arranged near the boundary. This includes:

[0017] For the 5th-order precision format, the same meshing method is used in the x, y, and z directions;

[0018] For a semi-uniform grid along the x-axis: for the computational domain [a, b], flux points are evenly distributed across the computational domain, and the coordinates of the (j+1)th flux point are x. j+1 / 2 = a + jΔx, j = 0, 1, 2, ..., N, where N is an integer greater than 1, and the number of flux points is N+1; Δx = (ba) / N represents the spatial step size along the x-axis; N represents the total number of solution points, where the 3rd to N-2th solution points are located exactly between two adjacent flux points, with coordinates x... j = a + (j - 1 / 2)Δx, j = 3, ..., N-2. The first two solution points and the last two solution points are located near the left and right boundaries of the computational domain, respectively. By arranging these four solution points in specific positions, the engineering simulation can conform to the conservation characteristics of fluid physical quantities. The solution point positions in the x-axis direction are:

[0019]

[0020] Where j represents the index of the solution point along the x-axis, x j Let be the coordinates of the j-th solution point along the x-axis. The semi-uniform grid formed by the solution point and the flux point ensures that the engineering simulation conforms to the conservation characteristics of fluid physical quantities.

[0021] In one embodiment, in step S3, the mesh cells defined by the solution points are used as the cells to be detected, and the cells to be detected are divided into cells inside the computational domain and cells at the boundary of the computational domain. The specific method is as follows:

[0022] For the 5th-order precision scheme, the same discontinuous probing method is used in the x, y, and z directions. For the x-axis direction: the N solution points in the x-axis direction divide the computational domain into N+1 grid cells. These N+1 grid cells are treated as N+1 probed cells, and discontinuous probing is performed one by one using the global discontinuous probing method. The first probed cell is the interval [a, x1], and the last probed cell is the interval [x...]. N The remaining N-1 detection unit intervals are both ends of the solution points, and the i-th detection unit is represented as the interval [x,b]. i-1 ,x i ], i = 2, 3, 4, ..., N; among the N+1 discontinuous detection units, the 6th to N-4th detection units are internal units of the computational domain, and the remaining 10 units are boundary units of the computational domain.

[0023] In one embodiment, S3

[0024] The method for detecting discontinuities within the computational domain includes: for the 6th to N-4th detection units, the same discontinuity detection method is used for each detection unit. For the i-th detection unit I... i =[x i-1 ,x i ], where i is an integer and 6≤i≤N-4: Select the left-facing intermittent detection template S1={u i-3 ,u i-2 ,u i-1 ,u i ,u i+1}, and the right-facing intermittent detection template S2={u i-2 ,u i-1 ,u i ,u i+1 ,u i+2}; where u i-3 ,u i-2 ,u i-1 ,u i ,u i+1 ,u i+2 These represent the conserved physical quantities at the (i-3)th to (i+2)th solution points in the computational domain, respectively.

[0025] The method for detecting discontinuities in the boundary elements of the computational domain includes: for the left boundary elements I1 = [a, x1], I2 = [x1, x2], I3 = [x2, x3], I4 = [x3, x4], and I5 = [x4, x5], elements I3, I4, and I5 each use pairwise different discontinuity detection templates; the left-facing discontinuity detection templates for elements I1 and I2 are the same as the right-facing discontinuity detection templates; for the five elements on the right boundary, discontinuity detection templates corresponding to the left boundary elements are obtained through symmetry.

[0026] In one embodiment, in step S4, polynomial functions Q1(X) and Q2(X) are constructed based on the left windward intermittent detection template and the right windward detection template, respectively; the specific method is as follows:

[0027] Discontinuity probing of cells within the computational domain:

[0028] The i-th unit to be detected i =[x i-1 ,x i ], where i is an integer and 6≤i≤N-4, and the left-facing intermittent detection template S1={u i-3 ,u i-2 ,u i-1 ,u i ,u i+1} and right-facing intermittent detection template S2={u i-2 ,u i-1 ,u i ,ui+1 ,u i+2 The conserved physical quantities and coordinates at the solution points on the} are determined, and polynomial functions Q1(X) and Q2(X) are constructed using the Lagrange interpolation method;

[0029] On the left-facing windward intermittent detection template S1={u i-3 ,u i-2 ,u i-1 ,u i ,u i+1 On the above, taking the (i-1)th solution point as the reference, a polynomial function Q1(X) is constructed using the Lagrange interpolation method, with the specific form as follows:

[0030]

[0031] in, These are the basis functions of the Lagrange interpolation method. These are the coefficients of each basis function; the superscript L indicates that the polynomial was constructed from the left-facing discontinuous detection template.

[0032] The polynomial satisfies: Q1(x) j )=u j j = i-3,…,i+1

[0033] The specific forms of the basis functions are as follows:

[0034]

[0035] The specific form of the correlation coefficient is:

[0036]

[0037] Where h = (ba) / N is the spacing between solution points that are uniformly distributed internally;

[0038] On the right-facing windward intermittent detection template S2={u i-2 ,u i-1 ,u i ,u i+1 ,u i+2 On the}, taking the i-th solution point as the reference, the polynomial function Q2(x) is constructed using the Lagrange interpolation method, with the specific form as follows:

[0039]

[0040] The superscript R indicates that the polynomial was constructed from a right-facing, intermittent detection template.

[0041] The polynomial satisfies: Q2(x) j )=u j j = i-2,…,i+2

[0042] The specific forms of the basis functions are as follows:

[0043]

[0044] The specific form of the correlation coefficient is:

[0045]

[0046] Discontinuity detection for computational domain boundary cells:

[0047] For the unit under test, a polynomial function is constructed based on the left-facing intermittent detection template. Construct a polynomial function based on the right-facing windward intermittent detection template. For each unit under test, polynomial functions Q1(X) and Q2(X) need to be constructed to establish a unified expression: based on the template node set {(x... j,1 ,u j,1 ),(x j,2 ,u j,2 ),(x j,3 ,u j,3 ),(x j,4 ,u j,4 ),(x j,5 ,u j,5 The polynomial function Q(x) = θ0g0(x) + θ1g1(x) + θ2g2(x) + θ3g3(x) + θ4g4(x) is constructed using Lagrange interpolation.

[0048] The polynomial function Q(x) satisfies the condition Q(x) j,k )=u j,k Given k = 1, 2, 3, 4, 5, we obtain the following system of linear equations:

[0049]

[0050] Among them, y k =f(x) j,k The coefficients are determined by the coordinate mapping relationship; by solving this system of equations, the explicit expressions for the coefficients are obtained:

[0051]

[0052] In one embodiment, S5, the first derivative of the polynomial function is obtained to obtain the first derivative functions Q′1(X) and Q′2(X); the specific method is as follows:

[0053] Discontinuity probing for cells within the computational domain and discontinuity probing for cells at the boundary of the computational domain:

[0054] Taking the first derivatives of the polynomial functions Q1(x) and Q2(x) with respect to x, we obtain the specific forms of the first derivative functions Q'1(x) and Q'2(x) as follows:

[0055]

[0056] In one embodiment, S6, the first derivative of the polynomial function is simplified to obtain a simplified form of the first derivative of the polynomial function. and The specific method is as follows:

[0057] Discontinuity detection for cells within the computational domain and discontinuity detection for cells at the computational domain boundary:

[0058] Multiplying both Q′1(x) and Q′2(x) by a factor of h transforms the coefficients of the polynomial function from fractional to integer, resulting in the simplified form of the first derivative. and

[0059]

[0060] and The signs of the function values ​​at the corresponding positions will be the same as those of Q′1(X) and Q′2(X), respectively.

[0061] In one embodiment, S7, according to the simplified form of the first derivative function. and Calculate the overdetermined maximum and minimum values ​​on the unit to be detected respectively; the specific method is as follows:

[0062] Discontinuity probing of cells within the computational domain:

[0063] polynomial function and Each is obtained by linear combination of four monomial functions, and the maximum and minimum values ​​of each monomial function are obtained at the endpoints of the unit to be tested.

[0064] polynomial function The overdetermined maximum value is:

[0065]

[0066] polynomial function The overdetermined maximum value is:

[0067]

[0068] Discontinuity detection for computational domain boundary cells:

[0069] In calculation and When the overdetermined maximum / minimum value is reached, for the 4th and 5th units, and The formula for calculating the overdetermined maximum and minimum values ​​is consistent with the formula for calculating the internal elements mentioned above; the third element The overdetermined extrema are as follows:

[0070]

[0071] polynomial function The overdetermined maximum value is:

[0072]

[0073] Based on the symmetrical distribution of the computational domain grid, the method for calculating the overdetermined maximum and minimum values ​​of the 3rd cell is also applicable to the (N-1)th cell. During the calculation, the coordinates of the endpoints of the 3rd cell are replaced with the coordinates of the corresponding endpoints of the (N-1)th cell.

[0074] For the 1st and 2nd units on the left boundary, and the Nth and (N+1th)th units on the right boundary, the polynomial function for these units... An improved method for calculating overdetermined extrema is proposed: These four units are used in the calculation When finding the maximum and minimum values ​​of each monomial, they are coupled into... Maximum and minimum value information of the corresponding item: The local maxima of each item are superimposed to In the calculation of the overdetermined maximum value, at the same time The local minimum values ​​of each item are superimposed to The calculation of the overdetermined minimum value; the calculation method for the overdetermined maximum and minimum values ​​of these 4 elements is as follows: For element I k =[a k ,b k ], k = 1, 2, N, N+1, a k With b k Representing unit I respectively k The coordinates of the two endpoints of , and a k k , The overdetermined extrema are as follows:

[0075]

[0076] In polynomials In the overdetermined extrema, the polynomial is coupled into the polynomial. Information:

[0077]

[0078] ​In one embodiment, S8, the detection unit is determined to be discontinuous based on the sign of the over-determined maximum / minimum value; the specific method is as follows:

[0079] Discontinuity detection for cells within the computational domain and discontinuity detection for cells at the boundary of the computational domain;

[0080] Introduce a modified discriminant based on the uniform mesh size h = (ba) / N; if (MAX1+h)·(MIN1+h)<0 or (MAX2+h)·(MIN2+h)<0, the element to be measured is a discontinuous element; if (MAX1+h)·(MIN1+h)≥0 and (MAX2+h)·(MIN2+h)≥0, the element to be measured is a smooth element.

[0081] The advantages of this invention compared to existing technologies are as follows: This invention designs a global discontinuity detection method applicable to both interior points and boundary high-order finite difference schemes. It can effectively detect discontinuities in both uniform internal mesh cells and non-uniform mesh cells near the boundary, achieving efficient and accurate detection of discontinuities near the computational domain boundary. Furthermore, this discontinuity detection method does not rely on any empirical parameters specific to the problem, thus it can be easily applied to various computational examples. This provides a technical foundation for constructing globally efficient and high-precision hybrid schemes. Attached Figure Description

[0082] Figure 1 This is a technical roadmap for intermittent detection in one embodiment of the present invention;

[0083] Figure 2 This is a template used for detecting discontinuities in internal units in one embodiment of the present invention;

[0084] Figure 3 This is a template used for left boundary unit discontinuity detection in one embodiment of the present invention;

[0085] Figure 4 This is a schematic diagram of the discontinuity identification result of a one-dimensional Lax shock tube problem in another embodiment of the present invention;

[0086] Figure 5 This is a schematic diagram of the discontinuity identification result of a one-dimensional Sod shock tube problem in another embodiment of the present invention;

[0087] Figure 6 This is a schematic diagram of the discontinuity identification result of the two-dimensional Riemann problem density calculation in another embodiment of the present invention;

[0088] Figure 7 This is a schematic diagram of the discontinuity identification results in the density calculation of the dual Mach reflection problem in another embodiment of the present invention. Detailed Implementation

[0089] To make the objectives, technical solutions, and advantages of this application clearer, the following detailed description is provided in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the scope of this application.

[0090] Parameter description: Let the numerical solution u j For the corresponding exact solution u(x) j The numerical approximation of ,t).

[0091] In one embodiment, such as Figure 1 As shown, a global discontinuity detection method applicable to both interior points and boundary high-order finite difference schemes is provided, comprising the following steps:

[0092] S1. Based on the task to be calculated, construct the hyperbolic conservation law equations of the simulation object and set the computational domain.

[0093] Specifically, the tasks to be calculated can include, but are not limited to, the flow fields of fighter jets flying at supersonic speeds and spacecraft re-entering the atmosphere.

[0094] In a Cartesian coordinate system, a three-dimensional problem can be decoupled into a one-dimensional problem along the x-axis, y-axis, and z-axis directions respectively. The same method is used to process the problems along the x-axis, y-axis, and z-axis directions. Therefore, this embodiment will take the x-axis direction as an example.

[0095] The one-dimensional hyperbolic conservation law equation along the x-axis is expressed as follows:

[0096] u t +f(u) x =0

[0097] Where u = u(x,t) represents the conserved physical quantity at coordinate x at time t, including fluid density, velocity, and pressure; f(u) represents the flux corresponding to the conserved physical quantity at coordinate x at time t, where flux characterizes the transport rate of the conserved physical quantity through a unit area per unit time; u t Let f(u) represent the derivative of the conserved physical quantity u with respect to time t. x This represents the derivative of flux f(u) with respect to x. The computational domain refers to the mapping region of the physical space in the engineering simulation onto the discrete model. Preferably, the computational domain along the x-axis can be represented as the interval [a, b], where a and b are real numbers and a... <b。

[0098] S2. Construct a semi-uniform grid based on the computational domain. A semi-uniform grid means dividing the computational domain into staggered solution points and flux points in each coordinate direction. The flux points are evenly spaced, and the solution points are arranged at the center of adjacent flux points. Then, the boundary solution points are replaced with non-uniform conservation points. The boundary solution points refer to the number of solution points closest to the boundary. The specific number is determined according to the precision of the format. For a global (2p-1) order format (where p is an integer and 2≤p≤6), (p-1) non-uniform points are needed near the boundary.

[0099] Specifically, taking a 5th-order precision format as an example, the same meshing method is used in the x, y, and z directions. Here, we introduce a semi-uniform mesh distribution using the x-axis as an example: for the computational domain [a, b], flux points are evenly distributed across the computational domain, and the coordinates of the (j+1)th flux point are x... j+1 / 2 = a + jΔx, j = 0, 1, 2, ..., N, where N is an integer greater than 1, and the number of flux points is N+1; Δx = (ba) / N represents the spatial step size in the x-axis direction; correspondingly, there are N solution points in total, of which the 3rd to N-2th solution points are located exactly in the middle of two adjacent flux points, with coordinates x... j = a + (j - 1 / 2)Δx, j = 3, ..., N-2. The first two solution points and the last two solution points are located near the left and right boundaries of the computational domain, respectively. These four solution points are no longer located exactly in the middle of two adjacent flux points, but are located in specific positions. By arranging these four solution points in specific positions, the engineering simulation can conform to the conservation characteristics of fluid physical quantities. In general, the solution point positions in the x-axis direction are:

[0100]

[0101] Where j represents the index of the solution point along the x-axis, x j Let be the coordinates of the j-th solution point along the x-axis. The semi-uniform grid formed by the solution point and the flux point ensures that the engineering simulation conforms to the conservation characteristics of fluid physical quantities.

[0102] S3. The grid cells divided according to the solution points are used as the detection cells. The detection cells are divided into cells inside the computational domain and cells at the boundary of the computational domain. Discontinuous detection is performed on each of the two types of detection cells one by one. For each detection cell, a left-facing discontinuous detection template and a right-facing discontinuous detection template are selected. The discontinuous detection template refers to the engineering control unit used to determine whether the detection cell contains discontinuities. It is a set of several continuous solution points near the detection cell. The specific number is determined according to the precision of the format. For the global (2p-1) order format, the left-facing discontinuous detection template and the right-facing discontinuous detection template are both composed of (2p-1) solution points. The presence or absence of discontinuities on the detection cell is determined according to the known physical quantities (usually density, velocity, and pressure, etc.) on these (2p-1) solution points and the positional relationship between these (2p-1) solution points and the detection cell.

[0103] Specifically, taking a 5th-order precision format as an example, the same discontinuous probing method is used in the x, y, and z directions. The x-axis direction will be used as an example. In step S2, the N solution points divide the computational domain into N+1 grid cells. These N+1 grid cells will serve as N+1 probed cells, and the global discontinuous probing method will be used to perform discontinuous probing one by one. Specifically, the first probed cell is the interval [a, x1], and the last probed cell is the interval [x...]. N The remaining N-1 detection unit intervals are both ends of the solution points. The i-th detection unit can be represented as the interval [x,b]. i-1 ,x i Let i = 2, 3, 4, ..., N. In these N+1 discontinuous detection units, the 6th to N-4th detection units are internal units of the computational domain, and the remaining 10 units are boundary units of the computational domain. Due to the different distribution patterns of the solution points, the selection methods for the discontinuous detection templates and the discontinuous detection methods for these two types of units will also be different. The specific implementation of steps S3 to S8 will be described below, starting from the internal units and the boundary units respectively.

[0104] 1) Detection of discontinuities within the computational domain

[0105] In step S3, the discontinuity detection unit within the computational domain includes the 6th to N-4th detection units. The same discontinuity detection method is used for each detection unit, with the i-th detection unit I... i =[x i-1 ,x i Let's take [i] as an example, where i is an integer and 6 ≤ i ≤ N-4. For example... Figure 2 As shown, for I i =[x i-1 ,x i Select the left-facing intermittent detection template S1 = {u i-3 ,ui-2 ,u i-1 ,u i ,u i+1}, and the right-facing intermittent detection template S2={u i-2 ,u i-1 ,u i ,u i+1 ,u i+2}, where u i-3 ,u i-2 ,u i-1 ,u i ,u i+1 ,u i+2 These represent the conserved physical quantities at the (i-3)th to (i+2)th solution points in the computational domain. It is easy to observe that templates S1 and S2 use conserved physical quantities at different solution points, which increases the accuracy of discontinuous detection.

[0106] S4. Based on the left-facing intermittent detection template and the right-facing intermittent detection template, construct polynomial functions Q1(X) and Q2(X) respectively.

[0107] Specifically, for the i-th probe unit I i =[x i-1 ,x i ], where i is an integer and 6≤i≤N-4, and the left-facing intermittent detection template S1={u i-3 ,u i-2 ,u i-1 ,u i ,u i+1} and right-facing intermittent detection template S2={u i-2 ,u i-1 ,u i ,u i+1 ,u i+2 The conserved physical quantities and coordinates at the solution points are determined using the Lagrange interpolation method, and polynomial functions Q1(X) and Q2(X) are constructed. The conserved physical quantities at the solution points include density, velocity, and pressure. Preferably, the density at the solution points on the template is chosen to construct the polynomial functions Q1(X) and Q2(X).

[0108] The idea behind discontinuity detection is to check whether there are extreme points of polynomial functions Q1(X) and Q2(X) in the unit to be tested. If either polynomial function Q1(X) or Q2(X) has an extreme point in the unit to be tested, then the unit is considered to be a discontinuous unit.

[0109] On the left-facing windward intermittent detection template S1={u i-3 ,u i-2 ,u i-1 ,u i ,u i+1On the above, taking the (i-1)th solution point as the reference, a polynomial function Q1(X) is constructed using the Lagrange interpolation method, with the specific form as follows:

[0110]

[0111] in, Let be the basis functions of the Lagrange interpolation method. The coefficients of each basis function are given; the superscript L indicates that the polynomial is constructed from a left-facing, intermittent detection template. This polynomial satisfies: Q1(x j )=u j j = i-3,…,i+1

[0112] The specific forms of the basis functions are as follows:

[0113]

[0114] The specific form of the correlation coefficient is:

[0115]

[0116] Where h = (ba) / N is the spacing between solution points that are uniformly distributed internally.

[0117] Similarly, in the right-facing intermittent detection template S2={u i-2 ,u i-1 ,u i ,u i+1 ,u i+2 On the}, taking the i-th solution point as the reference, the polynomial function Q2(x) is constructed using the Lagrange interpolation method, with the specific form as follows:

[0118]

[0119] The superscript R indicates that the polynomial was constructed using a right-facing, intermittent detection template. This polynomial satisfies: Q²(x) j )=u j j = i-2,…,i+2

[0120] The specific forms of the basis functions are as follows:

[0121]

[0122] The specific form of the correlation coefficient is:

[0123]

[0124] S5. Take the first derivatives of the polynomial functions Q1(x) and Q2(x) to obtain the first derivative functions Q′1(x) and Q′2(x) of the polynomial functions.

[0125] Based on the relationship between the extreme points of the original function and the zeros of the first derivative, it is possible to determine whether the original function has extreme points within the unit under test by detecting whether the first derivative has zeros.

[0126] Specifically, by taking the first derivative of the polynomial functions Q1(x) and Q2(x) with respect to x, we obtain the specific forms of the first derivative functions Q'1(x) and Q'2(x) as follows:

[0127]

[0128] S6. Simplify the first derivative of the polynomial function to obtain the simplified form of the first derivative of the polynomial function. and Simplification refers to converting the fractional coefficients of Q'1(x) and Q'2(x) into integer coefficients through multiplication without affecting the sign of the function value.

[0129] Both Q'1(x) and Q'2(x) are continuous functions. Taking Q'1(x) as an example, if the maximum and minimum values ​​of Q'1(x) on the test cell have opposite signs, then Q'1(x) must have a zero on the test cell; otherwise, it can be approximated that Q'1(x) has no zero on the test cell. Based on this criterion, we only need to consider the sign of the corresponding function value of the first derivative, without considering the magnitude of the specific function value. Multiplying a function by a positive number does not change the sign of the function value; therefore, multiplication can be used to simplify the forms of Q'1(x) and Q'2(x), making subsequent calculations more convenient.

[0130] Specifically, by multiplying both Q′1(x) and Q′2(x) by a factor of h, the coefficients of the polynomial function are transformed from fractional to integer coefficients, resulting in the simplified form of the first derivative function. and

[0131]

[0132] and The signs of the function values ​​at the corresponding positions will be exactly the same as those of Q′1(X) and Q′2(X), respectively, but and The function form is simpler and more intuitive.

[0133] S7. Based on the simplified form of the first derivative function and Calculate the overdetermined maximum and minimum values ​​on the unit under test respectively.

[0134] Specifically, polynomial functions and Each is obtained by linear combination of four monomials, and the maximum and minimum values ​​of each monomial can be obtained on the unit under test. For example, The maximum value on the unit under test must be less than or equal to the sum of the maximum values ​​of the four monomials that make it up. The minimum value on the unit under test must be greater than or equal to the sum of the minimum values ​​of the four monomials that make it up. The sum of the maximum values ​​of the four monomials is defined as... The overdetermined maximum value, the sum of the minimum values ​​of the four monomials is The minimum value of overdeterminedness. For The same overdetermined extrema definition method is adopted. Compared with the actual maximum and minimum values ​​of the function, the overdetermined extrema are larger, enhancing the sensitivity of the discontinuity detector to extreme points.

[0135] As can be easily seen from the specific forms of the aforementioned basis functions, each basis function is a strictly monotonic function on the unit under test. The maximum and minimum values ​​of the function are obtained at the two endpoints of the unit under test, and are 0 and 1 respectively. Therefore, the overdetermined extrema can be calculated using the following formula.

[0136] polynomial function The overdetermined maximum value is:

[0137]

[0138] polynomial function The overdetermined maximum value is:

[0139]

[0140] S8. Determine whether the detection unit is discontinuous based on the sign of the over-determined maximum / minimum value.

[0141] Specifically, as mentioned earlier, the presence of zeros in the first derivative function on the element under test is determined based on the sign of the overdetermined extrema. To avoid spurious oscillations in the extremum sign caused by floating-point rounding errors, a uniform mesh size h = (ba) / N is introduced to correct the discriminant. If (MAX1+h)·(MIN1+h)<0 or (MAX2+h)·(MIN2+h)<0, the element under test is a discontinuous element; if (MAX1+h)·(MIN1+h)≥0 and (MAX2+h)·(MIN2+h)≥0, the element under test is a smooth element.

[0142] 2) Boundary discontinuity detection

[0143] Near the boundary of the computational domain, due to the lack of template points, step S3 requires designing a discontinuous detection template specifically for the boundary detection unit. Furthermore, since there are non-uniformly distributed points at the boundary, step 4 requires re-establishing the mapping relationship specifically for the boundary detection unit. This demonstrates that the discontinuous detection method for boundary units differs from that for internal units. Boundary units include the first 5 units of the left boundary and the last 5 units of the right boundary. The following description uses the first 5 units of the left boundary as an example to illustrate the discontinuous detection method. The last 5 units of the right boundary can be discontinuously detected in the same way by symmetrically flipping the template.

[0144] Specifically, in step 3, the discontinuous detection templates of the left boundary units i1=[a,x1], i2=[x1,x2], I3=[x2,x3], I4=[x3,x4], and I5=[x4,x5] are as follows: Figure 3 As shown in the figure. Units I3, I4, and I5 each use two different templates; while for units I1 and I2, since the templates cannot be extended outside the computational domain, the templates are moved inside the computational domain to keep the total number of template points unchanged. This restriction also results in the left windward templates of units I1 and I2 being the same as the right windward templates.

[0145] In step 4, for the unit under test, a polynomial function is constructed based on the left-facing intermittent detection template. Construct a polynomial function based on the right-facing windward intermittent detection template. There are 5 boundary test units. For each test unit, polynomial functions Q1(X) and Q2(X) need to be constructed. This step can be summarized into a unified expression: from the template node set {(x... j,1 ,u j,1 ),(x j,2 ,u j,2 ),(x j,3 ,u j,3 ),(x j,4 ,u j,4 ),(x j,5 ,u j,5 The polynomial function Q(x) = θ0g0(x) + θ1g1(x) + θ2g2(x) + θ3g3(x) + θ4g4(x) is constructed using Lagrange interpolation.

[0146] like Figure 3 As shown in the figure, the arrows indicate the interpolation base points when constructing the polynomial function on the corresponding template. For the i-th probe unit, the specific form of the Lagrange interpolation basis function still satisfies the basis function form of the i-th probe unit located inside, but the corresponding coefficients need to be recalculated. The polynomial function Q(x) needs to satisfy the condition Q(x) j,k )=u j,kSince k = 1, 2, 3, 4, 5, we obtain the following system of linear equations:

[0147]

[0148] Among them, y k =f(x) j,k The coefficients can be determined by coordinate mapping relationships. By solving this system of equations, explicit expressions for the coefficients can be obtained:

[0149]

[0150] In the coefficient expression, the parameter y k See Table 1 for the specific format:

[0151] Table 1 shows the coefficients of the polynomial functions corresponding to the first 5 units on the left boundary, including the parameter y. k Specific value

[0152]

[0153] In step S5, for the cell to be detected, the first derivative of the constructed polynomial function is calculated to obtain the first derivative function. and

[0154] In step S6, consistent with the internal unit, the first-order derivatives Q′1 and Q′2 can be simplified. The simplification method is the same: multiply both Q′1 and Q′2 by h in the denominator of the coefficients, resulting in a simplified form of the first-order derivative with integer coefficients. and This simplification does not change the sign of the function value.

[0155] Calculate in step S7 and When the overdetermined maximum / minimum value is reached, for the 4th and 5th units, and The formula for calculating the overdetermined maximum and minimum values ​​is completely consistent with the formula for calculating the internal elements mentioned above; simply substitute the corresponding coefficients from Table 1. For the third element I3 = [x2, x3], and The calculation method for the overdetermined extrema remains the same as the internal element calculation method described above. However, since point x2 is a non-uniform point, the interval length of the third element is no longer h. Therefore, the values ​​of the basis functions at the endpoints of the interval are no longer 0 and 1, and need to be calculated separately. The specific values ​​of the basis functions at the endpoints of the interval are: Substituting into the formula, we can obtain the third unit. and The overdetermined maximum value.

[0156] Specifically, the third unit The overdetermined extrema are as follows:

[0157]

[0158] polynomial function The overdetermined maximum value is:

[0159]

[0160]

[0161] Based on the symmetrical distribution of the computational domain grid, the method for calculating the overdetermined maximum and minimum values ​​for the 3rd element is also applicable to the (N-1)th element.

[0162] The first and second elements on the left boundary, and the Nth and (N+1)th elements on the right boundary, are located on the outermost edge of the computational domain and are most constrained by boundary conditions. Therefore, the polynomial functions for these elements... An improved method for calculating overdetermined extrema is proposed. This method differs from the method used for internal element calculations. and Employing a decoupled computation strategy, these four units perform computation... When finding the maximum and minimum values ​​of each monomial, they are coupled into... Maximum and minimum value information of the corresponding item: The local maxima of each item are superimposed to In the calculation of the overdetermined maximum value, at the same time The local minimum values ​​of each item are superimposed to The calculation of the overdetermined minimum value. This coupling mechanism increases the... The overdetermined maximum value was simultaneously reduced. The overdetermined minimum value is obtained, thereby expanding the extreme response range of the boundary discontinuity detector.

[0163] Specifically, taking the second unit on the left boundary, I2 = [x1, x2], as an example, its polynomial function The overdetermined maximum / minimum value calculation scheme is the same as that described above:

[0164]

[0165] In polynomials In the overdetermined extrema, the polynomial is coupled into the polynomial. Information:

[0166]

[0167] For the 1st, Nth, and N+1th elements, the same coupling strategy as for the 2nd element above is used to calculate the overdetermined maximum and minimum values. Simply replace the coordinates of the element endpoints in the formula with the corresponding coordinates of the element endpoints.

[0168] For step S8, the sign of the overdetermined maximum / minimum value is used to determine whether the element to be detected is discontinuous, with the boundary discontinuity detector and the internal discontinuity detector remaining consistent. Specifically, the sign of the overdetermined maximum / minimum value is used to determine whether the first derivative function has a zero point on the element to be detected. To avoid false oscillations of extremum signs caused by floating-point rounding errors, a uniform grid size h = (ba) / N is introduced to correct the discriminant. If (MAX1+h)·(MIN1+h)<0 or (MAX2+h)·(MIN2+h)<0, the element to be detected is a discontinuous element; if (MAX1+h)·(MIN1+h)≥0 and (MAX2+h)·(MIN2+h)≥0, the element to be detected is a smooth element.

[0169] In a verification embodiment, the discontinuity identification results of the global discontinuity detection method are verified through different one-dimensional and two-dimensional examples. These examples are all calculated using the fifth-order multi-resolution WCNS format, and the global discontinuity detection method is used to detect discontinuities in the calculated waveforms, with discontinuity units labeled.

[0170] 1) Testing the problem of one-dimensional shock tube

[0171] Next, we will test the accuracy of the discontinuity detection method for discontinuity capture using two one-dimensional shock tube problems. The computational domain for both shock tube problems is set to x∈[0,1], and discretized into 200 grid points.

[0172] Solving the Lax shock tube problem up to t = 0.15, with the following initial conditions.

[0173]

[0174] Solving the Sod shock tube problem up to t = 0.2, with the following initial conditions.

[0175]

[0176] Where ρ represents density, u represents velocity, and p represents pressure. Figure 4 The results of density calculation for the Lax shock tube problem and the spatial distribution of discontinuous elements identified by the discontinuity detection method are presented. Figure 5 The results of density calculations for the Sod shock tube problem and the spatial distribution of discontinuous elements identified by the discontinuity detection method are presented. It can be seen that the discontinuity detection method can accurately capture key characteristic elements such as shock waves and contact discontinuities.

[0177] 2) Two-dimensional Riemann problem test

[0178] The problem was calculated up to t = 0.8, with the initial conditions being:

[0179]

[0180] Where ρ represents density, u represents velocity along the x-axis, v represents velocity along the y-axis, and p represents pressure. The grid resolution is 450×450. For two-dimensional problems, the discontinuity detection method operates independently along the x-axis and y-axis, and finally merges the discontinuous cells detected in the two directions to present the overall discontinuous distribution of the computational domain. This example involves a relatively complex flow structure. Figure 6 The diagram shows the discontinuous cells identified by the discontinuity detection method at the final state time. It demonstrates that the method accurately captures the discontinuous structures in the flow field. Notably, no false positives were detected at the boundaries.

[0181] 3) The problem of double Mach reflection

[0182] This problem describes a Mach 10 shock wave impacting an infinitely thick plate at a 60° incident angle, generating reflected waves on both sides of the plate and forming a double Mach reflected wave behind the plate. It allows for the examination of the scheme's ability to capture strong shock waves and its resolution on small-scale structures in the recirculation region. The computational domain is (x,y)∈[0,4]×[0,1]. To verify the effectiveness of the boundary discontinuity detection method, the flow field is divided into two subdomains along the line x=2, each with a resolution of 1000×500. The initial conditions are...

[0183]

[0184] Figure 7 The discontinuous elements captured by the discontinuity detection method at the final state time are shown. It can be seen that the discontinuity detection method accurately captures the discontinuous elements, especially in the complex flow field at the interface of the grid blocks (straight line x=2). The boundary discontinuity detection method accurately distinguishes discontinuous elements from smooth elements, proving the effectiveness of the invention.

[0185] Those skilled in the art should recognize that the embodiments described herein are intended to help the reader understand the principles of the invention, and should be understood that the scope of protection of the invention is not limited to such specific statements and embodiments. Those skilled in the art can make various other specific modifications and combinations based on the technical teachings disclosed in this invention without departing from the spirit of the invention, and these modifications and combinations are still within the scope of protection of this invention.

Claims

1. A global discontinuity detection method applicable to both interior point and boundary high-order finite difference schemes, characterized in that, It includes: S1. Based on the task to be calculated, construct the hyperbolic conservation law equations of the simulation object and set the computational domain; S2. Based on the computational domain, construct a semi-uniform grid. The semi-uniform grid refers to dividing the computational domain into intersecting solution points and flux points in each coordinate direction, where flux points are evenly spaced and solution points are arranged at the center of adjacent flux points. Then, the boundary solution points are replaced with non-uniformly conserved points. The boundary solution points refer to the solution points closest to the boundary. The specific number is determined according to the accuracy of the format. For a global (2p-1) order format, where 2≤p≤6, (p-1) non-uniformly conserved points are arranged near the boundary. S3. The grid cells divided according to the solution points are used as the cells to be detected. The cells to be detected are divided into cells inside the computational domain and cells at the boundary of the computational domain. Discontinuous detection is performed on each of the two types of cells. For each cell to be detected, a left-facing discontinuous detection template and a right-facing discontinuous detection template are selected. The discontinuous detection template refers to the engineering control unit used to determine whether the detection cell contains discontinuities. It is a set of several continuous solution points around the cell to be detected. The specific number is determined according to the precision of the format. For the global (2p-1) order format, the left-facing discontinuous detection template and the right-facing discontinuous detection template are both composed of (2p-1) solution points. The presence or absence of discontinuities on the detection cell is determined according to the known physical quantities on these (2p-1) solution points and the positional relationship between these (2p-1) solution points and the detection cell. S4. Based on the left windward intermittent detection template and the right windward detection template, construct polynomial functions Q1(X) and Q2(X) respectively; S5. Find the first derivative of the polynomial function to obtain the first derivative functions Q′1(X) and Q′2(X) of the polynomial function; S6. Simplify the first derivative of the polynomial function to obtain the simplified form of the first derivative of the polynomial function. and The simplification process refers to converting the fractional coefficients in Q′1(X) and Q′2(X) into integer coefficients through multiplication; S7. Based on the simplified form of the first derivative function... and Calculate the overdetermined maximum and minimum values ​​on the unit to be detected respectively; S8. Determine whether the unit to be detected is discontinuous based on the sign of the over-determined maximum and minimum values.

2. The method according to claim 1, characterized in that, In step S2, a semi-uniform grid is constructed based on the computational domain. The semi-uniform grid refers to dividing the computational domain into intersecting solution points and flux points in each coordinate direction, where flux points are evenly spaced and solution points are arranged at the center of adjacent flux points. Then, the boundary solution points are replaced with non-uniformly conserved points. The boundary solution points refer to the number of solution points closest to the boundary. The specific number is determined according to the accuracy of the format. For a global (2p-1) order format, where 2≤p≤6, (p-1) non-uniformly conserved points are arranged near the boundary. It includes: For the 5th-order precision format, the same meshing method is used in the x, y, and z directions; For a semi-uniform grid along the x-axis: for the computational domain [a, b], flux points are evenly distributed across the computational domain, and the coordinates of the (j+1)th flux point are x. j+1 / 2 = a + jΔx, j = 0, 1, 2, ..., N, where N is an integer greater than 1, and the number of flux points is N+1; Δx = (ba) / N represents the spatial step size along the x-axis; N represents the total number of solution points, where the 3rd to N-2th solution points are located exactly between two adjacent flux points, with coordinates x... j = a + (j - 1 / 2)Δx, j = 3, ..., N-2. The first two solution points and the last two solution points are located near the left and right boundaries of the computational domain, respectively. By arranging these four solution points in specific positions, the engineering simulation can conform to the conservation characteristics of fluid physical quantities. The solution point positions in the x-axis direction are: Where j represents the index of the solution point along the x-axis, x j Let be the coordinates of the j-th solution point along the x-axis. The semi-uniform grid formed by the solution point and the flux point ensures that the engineering simulation conforms to the conservation characteristics of fluid physical quantities.

3. The method according to claim 2, characterized in that, In step S3, the mesh cells defined based on the solution points are used as the cells to be probed. The cells to be probed are divided into cells within the computational domain and cells at the boundary of the computational domain. The specific method is as follows: For the 5th-order precision scheme, the same discontinuous probing method is used in the x, y, and z directions. For the x-axis direction: the N solution points in the x-axis direction divide the computational domain into N+1 grid cells. These N+1 grid cells are treated as N+1 probed cells, and discontinuous probing is performed one by one using the global discontinuous probing method. The first probed cell is the interval [a, x1], and the last probed cell is the interval [x...]. N The remaining N-1 detection unit intervals are both ends of the solution points, and the i-th detection unit is represented as the interval [x,b]. i-1 ,x i ], i = 2, 3, 4, ..., N; among the N+1 discontinuous detection units, the 6th to N-4th detection units are internal units of the computational domain, and the remaining 10 units are boundary units of the computational domain.

4. The method according to claim 3, characterized in that, In S3, the methods for detecting discontinuities between cells within the computational domain include: For the 6th to N-4th detection units, the same intermittent detection method is used for each detection unit. For the i-th detection unit I... i =[x i-1 ,x i ], where i is an integer and 6≤i≤N-4: Select the left-facing intermittent detection template S1={u i-3 ,u i-2 ,u i-1 ,u i ,u i+1 }, and the right-facing intermittent detection template S2={u i-2 ,u i-1 ,u i ,u i+1 ,u i+2 }; where u i-3 ,u i-2 ,u i-1 ,u i ,u i+1 ,u i+2 These represent the conserved physical quantities at the (i-3)th to (i+2)th solution points in the computational domain, respectively. The method for detecting discontinuities in the boundary elements of the computational domain includes: for the left boundary elements I1 = [a, x1], I2 = [x1, x2], I3 = [x2, x3], I4 = [x3, x4], and I5 = [x4, x5], elements I3, I4, and I5 each use pairwise different discontinuity detection templates; the left-facing discontinuity detection templates for elements I1 and I2 are the same as the right-facing discontinuity detection templates; for the five elements on the right boundary, discontinuity detection templates corresponding to the left boundary elements are obtained through symmetry.

5. The method according to claim 4, characterized in that, In step S4, polynomial functions Q1(X) and Q2(X) are constructed based on the left windward intermittent detection template and the right windward detection template, respectively; the specific method is as follows: Discontinuity probing of cells within the computational domain: The i-th unit to be detected i =[x i-1 ,x i ], where i is an integer and 6≤i≤N-4, and the left-facing intermittent detection template S1={u i-3 ,u i-2 ,u i-1 ,u i ,u i+1 } and right-facing intermittent detection template S2={u i-2 ,u i-1 ,u i ,u i+1 ,u i+2 The conserved physical quantities and coordinates at the solution points on the} are obtained, and polynomial functions Q1(X) and Q2(X) are constructed using the Lagrange interpolation method; On the left-facing windward intermittent detection template S1={u i-3 ,u i-2 ,u i-1 ,u i ,u i+1 On the above, taking the (i-1)th solution point as the reference, the polynomial function Q1(X) is constructed using the Lagrange interpolation method, with the specific form as follows: in, These are the basis functions of the Lagrange interpolation method. These are the coefficients of each basis function; the superscript L indicates that the polynomial was constructed from the left-facing discontinuous detection template. The polynomial satisfies: Q1(x) j )=u j j = i-3,…,i+1 The specific forms of the basis functions are as follows: The specific form of the correlation coefficient is: Where h = (ba) / N is the spacing between solution points that are uniformly distributed internally; On the right-facing windward intermittent detection template S2={u i-2 ,u i-1 ,u i ,u i+1 ,u i+2 On the}, taking the i-th solution point as the reference, the polynomial function Q2(x) is constructed using the Lagrange interpolation method, with the specific form as follows: The superscript R indicates that the polynomial was constructed from a right-facing, intermittent detection template. The polynomial satisfies: Q2(x) j )=u j j = i-2,…,i+2 The specific forms of the basis functions are as follows: The specific form of the correlation coefficient is: Discontinuity detection for computational domain boundary cells: For the unit under test, a polynomial function is constructed based on the left-facing intermittent detection template. Construct a polynomial function based on the right-facing windward intermittent detection template. For each unit under test, polynomial functions Q1(X) and Q2(X) need to be constructed to establish a unified expression: based on the template node set {(x... j,1 ,u j,1 ),(x j,2 ,u j,2 ),(x j,3 ,u j,3 ),(x j,4 ,u j,4 ),(x j,5 ,u j,5 The polynomial function Q(x) = θ0g0(x) + θ1g1(x) + θ2g2(x) + θ3g3(x) + θ4g4(x) is constructed using Lagrange interpolation. The polynomial function Q(x) satisfies the condition Q(x) j,k )=u j,k Given k = 1, 2, 3, 4, 5, we obtain the following system of linear equations: Among them, y k =f(x) j,k The coefficients are determined by the coordinate mapping relationship; by solving this system of equations, the explicit expressions for the coefficients are obtained:

6. The method according to claim 5, characterized in that, S5. Find the first derivative of the polynomial function to obtain the first derivative functions Q′1(X) and Q′2(X); the specific method is as follows: Discontinuity probing for cells within the computational domain and discontinuity probing for cells at the boundary of the computational domain: Taking the first derivatives of the polynomial functions Q1(x) and Q2(x) with respect to x, we obtain the specific forms of the first derivative functions Q'1(x) and Q'2(x) as follows:

7. The method according to claim 6, characterized in that, S6. Simplify the first derivative of the polynomial function to obtain the simplified form of the first derivative of the polynomial function. and The specific method is as follows: Discontinuity detection for cells within the computational domain and discontinuity detection for cells at the computational domain boundary: Multiplying both Q′1(x) and Q′2(x) by a factor of h transforms the coefficients of the polynomial function from fractional to integer, resulting in the simplified form of the first derivative. and and The signs of the function values ​​at the corresponding positions will be the same as those of Q′1(X) and Q′2(X), respectively.

8. The method according to claim 7, characterized in that, S7. Based on the simplified form of the first derivative function... and Calculate the overdetermined maximum and minimum values ​​on the unit to be detected respectively; the specific method is as follows: Discontinuity probing of cells within the computational domain: polynomial function and Each is obtained by linear combination of four monomial functions, and the maximum and minimum values ​​of each monomial function are obtained at the endpoints of the unit to be tested. polynomial function The overdetermined maximum value is: polynomial function The overdetermined maximum value is: Discontinuity detection for computational domain boundary cells: In calculation and When the overdetermined maximum / minimum value is reached, for the 4th and 5th elements, and The formula for calculating the overdetermined maximum and minimum values ​​is consistent with the formula for calculating the internal elements mentioned above; the third element The overdetermined extrema are as follows: polynomial function The overdetermined maximum value is: Based on the symmetrical distribution of the computational domain grid, the method for calculating the overdetermined maximum and minimum values ​​of the 3rd cell is also applicable to the (N-1)th cell. During the calculation, the coordinates of the endpoints of the 3rd cell are replaced with the coordinates of the corresponding endpoints of the (N-1)th cell. For the 1st and 2nd units on the left boundary, and the Nth and (N+1th)th units on the right boundary, the polynomial function for these units... An improved method for calculating overdetermined extrema is proposed: These four units are used in the calculation When finding the maximum and minimum values ​​of each monomial, they are coupled into... Maximum and minimum value information of the corresponding item: The local maxima of each item are superimposed to In the calculation of the overdetermined maximum value, at the same time The local minimum values ​​of each item are superimposed to The calculation of the overdetermined minimum value; the calculation method for the overdetermined maximum and minimum values ​​of these 4 elements is as follows: For element I k =[a k ,b k ], k = 1, 2, N, N+1, a k With b k Representing unit I respectively k The coordinates of the two endpoints of , and a k k , The overdetermined extrema are as follows:​ In polynomials In the overdetermined extrema, the polynomial is coupled into the polynomial. Information:

9. The method according to claim 1, characterized in that, S8. Based on the sign of the over-determined maximum / minimum value, determine whether the unit to be detected is discontinuous; the specific method is as follows: Discontinuity detection for cells within the computational domain and discontinuity detection for cells at the boundary of the computational domain; Introduce a modified discriminant based on the uniform mesh size h = (ba) / N; if (MAX1+h)·(MIN1+h)<0 or (MAX2+h)·(MIN2+h)<0, the element to be measured is a discontinuous element; if (MAX1+h)·(MIN1+h)≥0 and (MAX2+h)·(MIN2+h)≥0, the element to be measured is a smooth element.

Images