CN115329646B - Shock wave simulation method, device, equipment and medium - Google Patents

Shock wave simulation method, device, equipment and medium Download PDF

Info

Publication number
CN115329646B
CN115329646B CN202211249477.4A CN202211249477A CN115329646B CN 115329646 B CN115329646 B CN 115329646B CN 202211249477 A CN202211249477 A CN 202211249477A CN 115329646 B CN115329646 B CN 115329646B
Authority
CN
China
Prior art keywords
boundary
grid
physical parameter
determining
shock wave
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202211249477.4A
Other languages
Chinese (zh)
Other versions
CN115329646A (en
Inventor
杨灿群
仲彦旭
王伟
黄颖杰
夏梓峻
李健增
郑伟龙
卢海林
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Haihe Laboratory Of Advanced Computing And Key Software Xinchuang
National Supercomputer Center In Tianjin
National University of Defense Technology
Original Assignee
Haihe Laboratory Of Advanced Computing And Key Software Xinchuang
National Supercomputer Center In Tianjin
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Haihe Laboratory Of Advanced Computing And Key Software Xinchuang, National Supercomputer Center In Tianjin filed Critical Haihe Laboratory Of Advanced Computing And Key Software Xinchuang
Priority to CN202211249477.4A priority Critical patent/CN115329646B/en
Publication of CN115329646A publication Critical patent/CN115329646A/en
Application granted granted Critical
Publication of CN115329646B publication Critical patent/CN115329646B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation
    • G06F30/23Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation
    • YGENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02TCLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Abstract

The embodiment of the disclosure relates to a shock wave simulation method, a shock wave simulation device, equipment and a medium, wherein the method comprises the following steps: determining a shock wave grid in a current time step calculation region; when the boundary area of the calculation area comprises the shock wave grids, constructing a boundary Riemann problem corresponding to the boundary of the calculation area according to the first physical parameter and the second physical parameter; determining a solving result of the boundary Riemann problem, and determining a first virtual physical parameter corresponding to the virtual grid outside the calculation region according to the solving result; determining a discrete solution result of a preset fluid control equation according to the first virtual physical parameter of the virtual grid and the intermediate physical parameter of the current time step of the grid in the calculation region, and determining a shock wave simulation result of the next time step according to the discrete solution result. The embodiment of the invention effectively processes the problem of boundary propagation of the shock wave, inhibits non-physical reflection of the shock wave which is inconsistent with the real situation and appears at the boundary of the calculation area, and improves the accuracy of the simulation result.

Description

Shock wave simulation method, device, equipment and medium
Technical Field
The present disclosure relates to the field of simulation computing technologies, and in particular, to a method, an apparatus, a device, and a medium for simulating a shock wave.
Background
Computer Aided Engineering (CAE) technology can simulate phenomena such as explosion and combustion. In the simulation, the calculation area is of a limited size, however, the real physical space is usually an infinite open area or a semi-infinite open area, in order to improve the fit degree between the simulation and the real experiment, a plurality of layers of virtual grids can be additionally arranged on the boundary of the calculation area, and corresponding boundary conditions are applied, so that the effective processing of the boundary is realized.
In the related art, the Boundary of the calculation region may be processed using a Far-Field Boundary (FFBC) method. However, the far-field boundary method may form resistance against inflow/outflow of substances in a calculation region that is inconsistent with a real situation, which affects accuracy of CAE simulation, and the far-field boundary method has a poor effect of processing strong shock wave boundary propagation, which may generate significant boundary numerical reflection, thereby causing an increase in CAE simulation error.
Disclosure of Invention
To solve the above technical problem or at least partially solve the above technical problem, the present disclosure provides a shock wave simulation method, apparatus, device, and medium.
The embodiment of the disclosure provides a shock wave simulation method, which comprises the following steps:
determining a shock wave grid in a current time step calculation region;
when the boundary area of the calculation area comprises the shock wave grid, constructing a boundary Riemann problem corresponding to the boundary of the calculation area according to the first physical parameter and the second physical parameter; the first physical parameter is a physical parameter of a grid at infinity outside the calculation area, and the second physical parameter is a physical parameter of a grid current time step adjacent to the boundary in the calculation area;
determining a solving result of the boundary Riemann problem, and determining a first virtual physical parameter corresponding to the virtual grid outside the calculation area according to the solving result;
and determining a discrete solution result of a preset fluid control equation according to the first virtual physical parameter of the virtual grid and the intermediate physical parameter of the current time step of the grid in the calculation region, and determining a shock wave simulation result of the next time step according to the discrete solution result.
The embodiment of the present disclosure further provides a shock wave simulation apparatus, the apparatus includes:
the first determining module is used for determining a shock wave grid in a current time step calculation region;
the construction module is used for constructing a boundary Riemannian problem corresponding to the boundary of the calculation region according to the first physical parameter and the second physical parameter when the boundary region of the calculation region comprises the shock wave grid; the first physical parameter is a physical parameter of a grid at infinity outside the calculation region, and the second physical parameter is a physical parameter of a grid current time step adjacent to the boundary in the calculation region;
the solving module is used for determining a solving result of the boundary Riemannian problem and determining a first virtual physical parameter corresponding to the virtual grid outside the calculation area according to the solving result;
and the second determination module is used for determining a discrete solution result of a preset fluid control equation according to the first virtual physical parameter of the virtual grid and the intermediate physical parameter of the grid in the calculation region at the current time step, and determining a shock wave simulation result of the next time step according to the discrete solution result.
An embodiment of the present disclosure further provides an electronic device, including: a processor; a memory for storing the processor-executable instructions; the processor is used for reading the executable instructions from the memory and executing the instructions to realize the shock wave simulation method provided by the embodiment of the disclosure.
The embodiment of the present disclosure also provides a computer-readable storage medium, which stores a computer program for executing the shock wave simulation method provided by the embodiment of the present disclosure.
Compared with the prior art, the technical scheme provided by the embodiment of the disclosure has the following advantages: the shock wave simulation scheme provided in the embodiment of the disclosure determines a shock wave grid in a current time step calculation region; when the boundary area of the calculation area comprises the shock wave grids, constructing a boundary Riemann problem corresponding to the boundary of the calculation area according to the first physical parameter and the second physical parameter; the first physical parameter is a physical parameter of a grid at infinity outside the calculation region, and the second physical parameter is a physical parameter of a grid adjacent to the boundary in the calculation region at the current time step; determining a solving result of the boundary Riemann problem, and determining a first virtual physical parameter corresponding to the virtual grid outside the calculation region according to the solving result; determining a discrete solution result of a preset fluid control equation according to the first virtual physical parameter of the virtual grid and the intermediate physical parameter of the current time step of the grid in the calculation region, and determining a shock wave simulation result of the next time step according to the discrete solution result. By adopting the technical scheme, the boundary Riemannian problem is constructed under the condition that the boundary area of the calculation area comprises the shock wave grids, the boundary Riemannian problem is avoided being constructed under the condition that the shock wave is far away from the boundary of the calculation area, so that the influence of the boundary Riemannian problem on the inflow/outflow of substances into/from the calculation area is reduced, the accuracy of the CAE simulation result is improved, the shock wave boundary propagation problem can be effectively processed by constructing the boundary Riemannian problem corresponding to the boundary, the non-physical reflection which is not in accordance with the real condition and appears on the boundary of the calculation area of the shock wave is inhibited, and the accuracy of the CAE simulation result is further improved.
Drawings
The accompanying drawings, which are incorporated in and constitute a part of this specification, illustrate embodiments consistent with the present disclosure and together with the description, serve to explain the principles of the disclosure.
In order to more clearly illustrate the embodiments or technical solutions in the prior art of the present disclosure, the drawings used in the description of the embodiments or prior art will be briefly described below, and it is obvious for those skilled in the art that other drawings can be obtained according to the drawings without inventive exercise.
Fig. 1 is a schematic flow chart diagram of a shock wave simulation method according to an embodiment of the present disclosure;
FIG. 2 is a schematic diagram of a computing area and a virtual grid provided by an embodiment of the present disclosure;
FIG. 3 is a schematic diagram of a boundary region provided by an embodiment of the present disclosure;
FIG. 4 is a schematic flow chart diagram illustrating another method for simulating a shock wave according to an embodiment of the present disclosure;
fig. 5 is a schematic diagram of a boundary riemann problem provided by an embodiment of the present disclosure;
FIG. 6 is a schematic view of multiple substances in contact with one another provided by embodiments of the present disclosure;
FIG. 7 is a schematic flow chart diagram illustrating yet another method for simulating a shock wave according to an embodiment of the present disclosure;
FIG. 8 is a schematic illustration of an inflow/outflow boundary condition provided by an embodiment of the present disclosure;
FIG. 9 is a schematic flow chart diagram illustrating yet another method for simulating a shock wave according to an embodiment of the present disclosure;
FIG. 10 is a schematic illustration of a simulated explosive underwater explosion provided by an embodiment of the disclosure;
FIG. 11 is a cloud of pressure for an explosive shock wave flow field provided by an embodiment of the disclosure;
FIG. 12 is a schematic view of a pressure versus time curve provided by an embodiment of the present disclosure;
fig. 13 is a schematic structural diagram of a shock wave simulation apparatus according to an embodiment of the present disclosure;
fig. 14 is a schematic structural diagram of an electronic device according to an embodiment of the present disclosure.
Detailed Description
In order that the above objects, features and advantages of the present disclosure may be more clearly understood, aspects of the present disclosure will be further described below. It should be noted that the embodiments and features of the embodiments of the present disclosure may be combined with each other without conflict.
In the following description, numerous specific details are set forth in order to provide a thorough understanding of the present disclosure, but the present disclosure may be practiced otherwise than as described herein; it is to be understood that the embodiments disclosed in the specification are only a few embodiments of the present disclosure, and not all embodiments.
Computer Aided Engineering (CAE) technology can simulate phenomena such as explosion and combustion. In the simulation, the calculation area is of a limited size, however, the real physical space is usually an infinite open area or a semi-infinite open area, in order to improve the fit degree between the simulation and the real experiment, a plurality of layers of virtual grids can be additionally arranged on the boundary of the calculation area, and corresponding boundary conditions are applied, so that the effective processing of the boundary is realized. The boundary conditions have direct influence on the CAE simulation result, the error of numerical simulation can be greatly reduced by adopting reasonable boundary conditions,
when the simulation problems (especially the high-energy material combustion/explosion problem, the high/ultrahigh sound velocity shock wave and other strong nonlinear problems) in an infinite open area, a semi-infinite open area and a water area are solved, the reasonable boundary processing method is set, and the non-reflection outgoing absorption of the inflow/outflow of materials and the intermittent disturbance of the shock wave and the like can be realized at the boundary of a calculation area. In the related art, the Non-Reflection Boundary Condition (NRBC) mainly includes a Perfect Matched Layer (PML) method, a Far-Field Boundary (FFBC) method, and the like. The perfect matching layer method can better realize effective transmission and absorption of shock waves, but in the PML method, a perfect matching layer calculation grid with a larger range needs to be arranged at the boundary of a calculation area, so that the calculation amount is increased. The far field boundary method has poor effect of processing boundary propagation of strong shock waves and can generate obvious boundary numerical reflection. In the far-field boundary method, the computational domain needs to be increased in order to achieve better reflection-free effect, and the increase of the computational domain causes the reduction of computational efficiency. In addition, the two non-reflection boundary conditions can form non-physical numerical resistance which is not consistent with a real experiment on the inflow/outflow of substances at the boundary of the calculation region, and great influence is caused on the CAE simulation result.
In order to solve the above problem, embodiments of the present disclosure provide a method for simulating a shock wave, and the method is described below with reference to specific embodiments.
Fig. 1 is a schematic flow chart of a shockwave simulation method provided by an embodiment of the present disclosure, which may be executed by a shockwave simulation apparatus, where the apparatus may be implemented by software and/or hardware, and may be generally integrated in an electronic device. As shown in fig. 1, the method includes:
step 101, determining a shock wave grid in a current time step calculation region.
The shock wave simulation method provided by the embodiment of the disclosure can be applied to a computing system such as a super computer, wherein the shock wave refers to propagation of a discontinuous peak in a medium, and the shock wave can also be understood as discontinuous disturbance. The shockwave simulation method can be used in different scenarios including, but not limited to: any one of water explosion simulation, air explosion simulation, high-altitude explosion simulation, near-ground explosion simulation, deep water explosion simulation, near-water surface explosion simulation, bubble dynamics analysis, offshore ice breaking analysis, high-sound-velocity collision simulation, hypersonic-velocity collision simulation, energy-gathered jet simulation and wing streaming simulation.
In the embodiment of the present disclosure, the calculation region may be a region for performing simulation on the shock wave in the simulation process, and the calculation region includes a plurality of grids. The grid of the computing area corresponds to an intermediate physical parameter, which is also called a conservation variable, and the intermediate physical parameter may be a parameter that characterizes the grid in the computing area from physical dimensions, and includes but is not limited to: at least one of density, velocity, pressure, internal energy. In order to enable a computing area with limited space to simulate an infinite open area or a semi-infinite open area, a multi-layer virtual grid may be additionally arranged outside the computing area, and the virtual grid may be a grid arranged outside the computing area and used for assisting the simulation computing in the computing area. The virtual grid corresponds to virtual physical parameters, which may be parameters characterizing the virtual grid from physical dimensions, including but not limited to: at least one of density, velocity, pressure, internal energy. The boundary of the calculation region may be a boundary between the mesh and the virtual mesh within the calculation region. Fig. 2 is a schematic diagram of a computing area and a virtual grid according to an embodiment of the present disclosure, as shown in fig. 2, where a grid outside a boundary of the computing area is the virtual grid, and a boundary between the grid and the virtual grid in the computing area is a boundary of the computing area. In the disclosed embodiment, first, the shockwave mesh within the current time step calculation region is determined.
In the process of performing the simulation, the shock wave simulation apparatus can perform the simulation calculation based on a time step, wherein the time step can be a step length advancing in a time dimension for each simulation, and the current time step can calculate a corresponding time step for the current simulation. The grid of shock waves may be a grid in which shock waves are present.
In the embodiment of the present disclosure, the shock wave simulation apparatus may determine whether shock waves exist in each grid according to the intermediate physical parameters of each grid in the current time step calculation region, and determine the grid in which the shock waves exist as the shock wave grid.
In some embodiments, determining a grid of shock waves within the current time step calculation region includes:
sequentially determining each grid in the calculation area as a grid to be processed; determining an approximate physical parameter of the grid to be processed according to the intermediate physical parameters of the grid to be processed and the current time step of the grid before the grid to be processed; if the absolute value of the difference value between the intermediate physical parameter of the current time step of the grid to be processed and the approximate physical parameter is greater than or equal to a preset threshold value, determining the grid to be processed as a shock wave grid; the preset threshold is the product of the multi-resolution parameter and the grid size parameter.
Wherein, the grid to be processed may be a grid for which whether the shockwave exists is to be determined. In the embodiment of the present disclosure, grids may be sorted according to a position relationship between grids in the calculation region, and thus there may be corresponding sequence numbers in the grids. The mesh before the mesh to be processed may be a mesh whose sequence number is one bit before the mesh to be processed, and if the meshes are sorted by positive integers starting from 1, the mesh before the mesh to be processed may be a mesh whose sequence number is 1 smaller than that of the mesh to be processed. The approximate physical parameters may be physical parameters determined by performing weighted summation on intermediate physical parameters of the current time step of the mesh to be processed and intermediate physical parameters of the current time step of the mesh before the mesh to be processed. The absolute value of the difference value between the intermediate physical parameter of the grid to be processed and the approximate physical parameter can represent the change value of the intermediate physical parameter of the grid to be processed relative to the intermediate physical parameter of the previous grid. The preset threshold may be a criterion for determining whether shock waves exist in the grid to be processed. The multi-resolution parameter may be a preset parameter for excluding an influence of the mesh size on the impulse wave mesh determination. The mesh size may be a parameter characterizing the size of the mesh size.
In this embodiment, the shock wave simulation apparatus may obtain the intermediate physical parameters of the current time step of the mesh to be processed, and the intermediate physical parameters of the current time step of the mesh that is previous to the mesh to be processed. The intermediate physical parameter may be set according to an application scenario, and the embodiment is not limited, for example, the intermediate physical parameter may be set as a pressure. Further, the intermediate physical parameters of the grid to be processed and the intermediate physical parameters of the previous grid are weighted and summed, and the weighted and summed result is determined as the approximate physical parameters of the grid to be processed. And further, performing difference calculation on the intermediate physical parameters of the grid to be processed and the approximate physical parameters obtained by calculation to obtain difference results, comparing the absolute values of the difference results with a preset threshold, and if the absolute values of the difference results are greater than or equal to the preset threshold, determining that the grid to be processed is a shock wave grid. And sequentially taking the grids in the calculation area as the grids to be processed, and judging whether each grid to be processed is a shock wave grid or not according to the method, so that all shock wave grids in the calculation area can be determined.
In an alternative embodiment, the approximation error is
Figure 906830DEST_PATH_IMAGE001
Comprises the following steps:
Figure 89550DEST_PATH_IMAGE002
wherein the content of the first and second substances,
Figure 400445DEST_PATH_IMAGE003
is an intermediate physical parameter of the current time step of the mesh to be processed,
Figure 61234DEST_PATH_IMAGE004
approximating physical parameters for the grid to be processed, the approximating physical parameters
Figure 875606DEST_PATH_IMAGE005
Comprises the following steps:
Figure 229227DEST_PATH_IMAGE006
wherein the content of the first and second substances,
Figure 761839DEST_PATH_IMAGE007
intermediate physical parameters of the current time step of the previous grid of the grid to be processed according to the approximation error
Figure 727784DEST_PATH_IMAGE008
Judging whether the grid to be processed is a shock wave grid or not, if so, judging whether the grid to be processed is a shock wave grid or not
Figure 396662DEST_PATH_IMAGE009
The grid to be processed is a shock wave grid, wherein
Figure 186764DEST_PATH_IMAGE010
For multi-resolution parameters, can be taken
Figure 206672DEST_PATH_IMAGE011
Figure 209264DEST_PATH_IMAGE012
Is the mesh size.
102, when the boundary area of the calculation area comprises a shock wave grid, constructing a boundary Riemann problem corresponding to the boundary of the calculation area according to the first physical parameter and the second physical parameter; the first physical parameter is a physical parameter of a grid at infinity outside the calculation region, and the second physical parameter is a physical parameter of a grid adjacent to the boundary of the calculation region at the current time step.
The boundary area may be an area composed of a plurality of grids determined based on the boundary of the calculation area. In some embodiments, the boundary region may include: and calculating the grids of the preset layer number adjacent to the boundary in the area. The predetermined number of layers may be related according to the order of the subsequent fluid control equation solution method. For example, the fluid control equation is solved by using a Weighted essential Non-oscillator (WENO) finite difference format, where if the order of the WENO is 5, the preset number of layers may be 4 or 3, and if the order of the WENO is lower, the preset number of layers may be 1 or 2. Exemplarily, fig. 3 is a schematic diagram of a boundary area provided by an embodiment of the present disclosure, as shown in fig. 3, the computation area includes t layers of grids, the boundary area includes N layers of grids adjacent to the boundary in the computation area, then the 1 st to nth layers of grids belong to the boundary area, and the (t + 1-N) th to t th layers of grids belong to the boundary area.
The first physical parameter is a physical parameter of the grid at infinity outside the calculation area, and understandably, the grid at infinity is not affected by the shock wave caused by the reasons such as explosion, and the like, so that the first physical parameter may be the same as the initial physical parameter of the grid when the explosion does not occur, which is preset by the user. The second physical parameter is a physical parameter of the grid within the calculation region and adjacent to the boundary of the calculation region at the current time step, and the second physical parameter may be understood as a physical parameter determined through iterative simulation of one or more time steps based on the initial physical parameter set by the user. The riemann problem, also called shock tube problem, by which the movement state of the contact surface of different substances can be described, in which different substances that are in contact with each other can be classified into a riemann problem left side and a riemann problem right side. The boundary Riemannian problem is a Riemannian problem constructed based on the boundary of a computation region.
In this embodiment, if the grids forming the boundary region include a shockwave grid, the boundary riemann problem may be constructed at the boundary of the calculation region by using the first physical parameter and the second physical parameter as physical parameters corresponding to both sides of the riemann problem.
Fig. 4 is a schematic flow chart of another method for simulating a shock wave according to an embodiment of the present disclosure. As shown in fig. 4, constructing a boundary riemann problem corresponding to the boundary of the calculation region according to the first physical parameter and the second physical parameter includes the following steps:
step 401, the first physical parameter is used as a first riemann parameter value corresponding to the first side in the boundary riemann problem.
Step 402, taking the second physical parameter as a second riemann parameter value corresponding to a second side in the boundary riemann problem; wherein, if the first side is the left side in the boundary Riemannian problem, the second side is the right side in the boundary Riemannian problem; if the first side is the right side in the boundary Riemannian problem, the second side is the left side in the boundary Riemannian problem.
In this embodiment, the first side and the second side correspond to a left side or a right side of the riemann problem, and correspondingly, the first riemann parameter value and the second riemann parameter value correspond to a left riemann parameter value or a right riemann parameter value of the riemann problem, respectively. Specifically, if the first physical parameter is set as a left riemann parameter value corresponding to the left side in the boundary riemann problem; the second physical parameter is set to the right riemann parameter value corresponding to the right side in the boundary riemann problem. If the first physical parameter is set as a right Riemann parameter value corresponding to the right side in the boundary Riemann problem; the second physical parameter is set to the left-side riemann parameter value corresponding to the left side in the boundary riemann problem.
In an alternative embodiment of the method according to the invention,
Figure 998228DEST_PATH_IMAGE013
representing the left-hand riemann parameter in the boundary riemann problem,
Figure 693652DEST_PATH_IMAGE014
representing the right-hand riemann parameter in the boundary riemann problem,
Figure 935277DEST_PATH_IMAGE015
which is indicative of a first physical parameter of the device,
Figure 741559DEST_PATH_IMAGE016
representing a second physical parameter. If in the calculation area, a row includes t grids, and the grid with sequence number 1 is located at the leftmost side and the grid with sequence number t is located at the rightmost side, then
Figure 650609DEST_PATH_IMAGE017
A physical parameter representing the current time step of the leftmost mesh adjacent to the boundary,
Figure 516934DEST_PATH_IMAGE018
a physical parameter representing the current time step of the rightmost mesh adjacent to the boundary.
Fig. 5 is a schematic diagram of a boundary riemann problem provided by the embodiment of the present disclosure, as shown in fig. 5, if a left side of the boundary riemann problem is a virtual grid, then:
Figure 511435DEST_PATH_IMAGE019
Figure 855829DEST_PATH_IMAGE020
if the right side of the boundary Riemannian problem is a virtual grid, then:
Figure 619385DEST_PATH_IMAGE021
Figure 892497DEST_PATH_IMAGE022
and 103, determining a solving result of the boundary Riemannian problem, and determining a first virtual physical parameter corresponding to the virtual grid outside the calculation region according to the solving result.
The solution result may include a left solution result determined according to a left riemann parameter on the left side of the boundary riemann problem and a right solution result determined according to a right riemann parameter on the right side of the riemann problem. The first virtual physical parameter may be a physical parameter corresponding to a virtual grid determined by solving a boundary riemann problem.
In the embodiment of the disclosure, after the boundary riemann problem corresponding to the boundary of the calculation region is constructed, the boundary riemann problem can be solved, so that the solving results of the left side and the right side of the boundary riemann problem are determined, and then the first virtual physical parameters corresponding to the virtual grid are determined according to the solving results of the left side and the right side.
In some embodiments, determining the solution to the boundary Riemannian problem comprises: solving the boundary Riemann problem by adopting a Riemann problem solver to obtain a solved result; wherein the Riemann problem solver comprises: the Rankine-Raugueh relationship, and the contact break speed versus pressure continuity relationship.
The Riemann problem solver is a functional module for solving the Riemann problem. In order to give consideration to the accuracy and robustness of solution, an HLLC (Harten-Lax-van Leer-contact) solver can be selected as the Riemann problem solver. The HLLC solver comprises a Rankine-Yugonniu relation formula and a continuous relation between contact break speed and pressure.
The Rankine-Yugonniu relation is also called R-H relation, which is:
Figure 374294DEST_PATH_IMAGE023
wherein F represents flux, s represents shock wave velocity, U represents intermediate physical parameters corresponding to the grid,
Figure 522378DEST_PATH_IMAGE024
which represents the flux after the shock wave,
Figure 140442DEST_PATH_IMAGE025
representing the intermediate physical parameter after the shock wave.
The continuous relation between the contact discontinuity speed and the pressure is as follows:
Figure 348569DEST_PATH_IMAGE026
wherein the content of the first and second substances,
Figure 317662DEST_PATH_IMAGE027
representing the velocity to the left of the boundary riemann problem after the shockwave,
Figure 3858DEST_PATH_IMAGE028
representing the velocity to the right of the boundary riemann problem after the shock wave,
Figure 742007DEST_PATH_IMAGE029
representing the pressure to the left of the boundary riemann problem after the shock wave,
Figure 855457DEST_PATH_IMAGE030
representing the pressure on the right side of the boundary riemann problem after the shock wave.
In this embodiment, solving the boundary riemann problem by the riemann problem solver can determine that the left physical parameter of the boundary riemann problem after the interaction of the substances on the two sides is as follows
Figure 46267DEST_PATH_IMAGE031
Right side physical parameters of
Figure 536154DEST_PATH_IMAGE032
In some embodiments, determining the first virtual physical parameter corresponding to the virtual grid outside the calculation region according to the solution result includes: and extracting the computational physical parameters corresponding to the first physical parameters in the solving result, and determining the computational physical parameters as the first virtual physical parameters corresponding to the virtual grid.
The calculating physical parameter may be a solution result corresponding to the first physical parameter after the first physical parameter is solved by the boundary riemann problem.
In the present embodiment, as shown in fig. 5, if the first physical parameter is a left-side parameter on the left side of the boundary riemann problem, the calculated physical parameter corresponding to the first physical parameter in the solution result
Figure 128809DEST_PATH_IMAGE033
(ii) a If the first physical parameter is the right-side parameter on the right side of the boundary Riemann problem, the calculated physical parameter corresponding to the first physical parameter in the solving result
Figure 678739DEST_PATH_IMAGE034
First virtual physical parameters corresponding to the virtual grid
Figure 622424DEST_PATH_IMAGE035
Comprises the following steps:
Figure 140169DEST_PATH_IMAGE036
. Namely, the first virtual physical parameter corresponding to the virtual grid outside the boundary of the calculation area is equal to the calculation physical parameter corresponding to the first physical parameter.
And 104, determining a discrete solution result of a preset fluid control equation according to the first virtual physical parameter of the virtual grid and the intermediate physical parameter of the current time step of the grid in the calculation region, and determining a shock wave simulation result of the next time step according to the discrete solution result.
The preset fluid control equation is an equation based on physical parameters of the grid and capable of performing simulation calculation according to time steps, the equation can be set according to an application scenario of the shock wave simulation, the embodiment is not limited, and for example, the preset fluid equation may be one of a navier-stokes equation, an euler equation, and an elasto-plastic fluid equation. The discrete solution result is obtained after the preset fluid equation is subjected to time and space discrete processing.
In this embodiment, after determining the first virtual physical parameter of the virtual grid, the discrete solution of the time and space dimensions may be performed on the preset fluid control equation according to the first virtual physical parameter and the intermediate physical parameter of the current time step to obtain a discrete solution result, the intermediate physical parameter of the grid in the next time step in the calculation region may be obtained by calculation according to the discrete solution result, and the intermediate physical parameter of the next time step may be used as a shock wave simulation result of the next time step
In some embodiments, determining a discrete solution result of the preset fluid control equation according to the first virtual physical parameter of the virtual grid and the intermediate physical parameter of the current time step of the grid in the calculation region, and determining a shock wave simulation result of the next time step according to the discrete solution result includes: decoupling multiple substances which are contacted with each other in the calculation area; solving a time dimension and a space dimension of a preset fluid control equation based on the decoupled multiple substances, the intermediate physical parameters and the first virtual physical parameters to obtain a discrete equation corresponding to the preset fluid control equation; and determining an intermediate physical parameter of the next time step of the grid in the calculation region based on the discrete equation.
It is understood that the substance causing the shock wave and the transmission medium of the shock wave are different substances in general, fig. 6 is a schematic diagram of multiple substances contacting with each other provided by the embodiment of the present disclosure, as shown in fig. 6, the substance causing the shock wave may be a substance B (e.g., explosive), the transmission medium of the shock wave is a substance a (e.g., water), a mixing unit is present in a contact area of the substance a and the substance B, and the mixing unit is a mixture of the substance a and the substance B, and a multi-substance interface can be determined by calculation in the mixing unit. The multiple substances in contact with each other in the calculation region can be understood as different substances in contact with each other. The discrete equation may be an equation obtained by performing discrete processing in a spatial dimension and a time dimension on a preset fluid control equation.
In this embodiment, a multi-material interface processing (texture coherence) method may be used to decouple multiple materials in contact with each other in the calculation region. Specifically, the physical parameters of the multi-substance interface can be processed by using the texture method, as shown in fig. 6, the multi-substance interface is decoupled into the single substance, so that non-physical oscillation of the multi-substance interface is avoided, and the calculation accuracy and stability of the shock wave simulation are improved. Further, a numerical value space discrete method is adopted to carry out space discrete solution on the preset fluid control equation based on the decoupled multiple substances, the intermediate physical parameters and the first virtual physical parameters. Specifically, in order to improve the precision of the shock wave simulation, a 5-order WENO finite difference format can be adopted for solving, and the numerical flux is obtained
Figure 852910DEST_PATH_IMAGE037
So as to complete the space discrete solution and obtain the ordinary differential equation system of the intermediate physical parameter (namely the conservation variable) with respect to the time derivative
Figure 573742DEST_PATH_IMAGE038
Comprises the following steps:
Figure 739144DEST_PATH_IMAGE039
wherein i is a mark corresponding to the grid one by one,
Figure 570834DEST_PATH_IMAGE040
the size of the grid is represented as,
Figure 872502DEST_PATH_IMAGE041
indicating correspondence of i node plus half node
Figure 764235DEST_PATH_IMAGE042
Figure 416933DEST_PATH_IMAGE043
Representing i nodes minus half a node
Figure 786734DEST_PATH_IMAGE044
And further, carrying out time dimension propulsion solution on the ordinary differential equation set related to the time derivative by adopting a numerical time discrete method based on the decoupled multiple substances, the intermediate physical parameters and the first virtual physical parameters. Specifically, a 3-order invariant to Total Variation (TVD) Runge-Kutta format may be adopted to discretize a time derivative term of a system of ordinary differential equations about time derivatives to obtain a complete discrete format of a preset fluid control equation, where the complete discrete format is a discrete equation, and then an intermediate physical parameter of each grid in a next time step calculation region is obtained based on the discrete equation, and a multi-material interface is simultaneously promoted.
In some embodiments, the shockwave simulation method further comprises: and if the current time length corresponding to the current time step is greater than the preset ending time length, determining the shock wave simulation result of the next time step as a target simulation result. And the current time length corresponding to the current time step is the accumulated sum of the current time step and the time step before the current time step. The preset ending time is a parameter value for judging whether the shock wave simulation is ended, and the preset ending time may be set according to user requirements and the like, which is not limited in this embodiment.
In this embodiment, if the current time step is
Figure 208488DEST_PATH_IMAGE045
The preset end time is
Figure 271122DEST_PATH_IMAGE046
Will be
Figure 411117DEST_PATH_IMAGE047
And
Figure 850188DEST_PATH_IMAGE048
and adding n +1 time steps in the previous time steps to obtain the current time length. Judgment ofWhether the current time length exceeds the preset ending time length
Figure 126449DEST_PATH_IMAGE049
If the current time length is longer than the preset ending time length, taking the shock wave simulation result of the next time step determined by calculation as a target simulation result; and if the current time length is less than the preset ending time length, continuing to perform shock wave simulation calculation based on the shock wave simulation result of the next time step.
The shock wave simulation scheme provided in the embodiment of the disclosure determines a shock wave grid in a current time step calculation region; when the boundary area of the calculation area comprises the shock wave grids, constructing a boundary Riemann problem corresponding to the boundary of the calculation area according to the first physical parameter and the second physical parameter; the first physical parameter is a physical parameter of a grid at infinity outside the calculation region, and the second physical parameter is a physical parameter of a grid adjacent to the boundary in the calculation region at the current time step; determining a solving result of the boundary Riemann problem, and determining a first virtual physical parameter corresponding to the virtual grid outside the calculation region according to the solving result; determining a discrete solution result of a preset fluid control equation according to the first virtual physical parameter of the virtual grid and the intermediate physical parameter of the current time step of the grid in the calculation region, and determining a shock wave simulation result of the next time step according to the discrete solution result. By adopting the technical scheme, the boundary Riemannian problem is constructed under the condition that the boundary area of the calculation area comprises the shock wave grids, the boundary Riemannian problem is avoided being constructed under the condition that the shock wave is far away from the boundary of the calculation area, so that the influence of the boundary Riemannian problem on the inflow/outflow of substances into/from the calculation area is reduced, the accuracy of the CAE simulation result is improved, the shock wave boundary propagation problem can be effectively processed by constructing the boundary Riemannian problem corresponding to the boundary, the non-physical reflection which is not in accordance with the real condition and occurs on the boundary of the calculation area is inhibited, and the accuracy of the CAE simulation result is further improved.
In addition, the increase of the calculation amount caused by setting the perfect matching layer in the perfect matching layer method is also avoided.
In some embodiments, before determining the shockwave mesh within the current time step calculation region, further comprising: initializing simulation parameters, wherein the simulation parameters comprise: calculating at least one of a region parameter, a grid parameter, a multi-material interface parameter, an initial condition parameter, a fluid control equation parameter, a parallel calculation parameter, and a time parameter; wherein the initial condition parameters comprise initial physical parameters of the grids in the calculation region.
The simulation parameters can be parameters which need to be configured by a user in advance in the process of shockwave simulation. In this embodiment, the simulation parameters may be set according to a specific application scenario of the shock wave simulation method. The concrete description is as follows:
the calculation region parameter may be a parameter characterizing a calculation region. The mesh parameter may be a parameter characterizing a mesh size. In this embodiment, an euler coordinate system may be established according to a specific application scenario, and the size of the calculation region and the size of the grid may be defined.
The multi-substance interface parameter may be a parameter characterizing the contact interface of different substances within the calculation region. In this embodiment, a multi-substance interface tracking method may be used to define the initial position of the multi-substance interface, and optionally, the initial position of the multi-substance interface may be determined according to the Volume fractions of different substances based on a Fluid Volume of Fluid (VOF) interface capture method and a Youngs interface reconstruction method.
The initial condition parameters may include initial physical parameters, initiation location parameters, and initiation type parameters of the grid within the calculation region. The initial physical parameters may be physical parameters of grids in a calculation region when no shock wave is generated, the detonation position parameters may be parameters representing explosive position characteristics, and the detonation type parameters may be parameters representing explosive type characteristics of explosives. In this embodiment, the initial physical parameters may be set according to an application scenario, which is not limited in this embodiment, for example, the initial physical parameters may include: material type, density, pressure, velocity, volume fraction. The detonation position parameters and the detonation type parameters can be set according to actual problem requirements.
In the present embodiment, the fluid control equation parameters include, but are not limited to: one of a Navier-Stokes equation, an Euler equation and an elasto-plastic fluid equation, wherein parameters of the fluid control equation can be set according to characteristics of an actual physical problem.
The parallel computing parameters are parameters characterizing the type of parallel computing. In this embodiment, the parallel computing parameters may be set according to the amount of computation.
The time parameters are parameters representing shock wave simulation time dimension characteristics, and the time parameters comprise time step length and preset ending time length. In this embodiment, the CFL (cournt Friedrichs Lewy) time advance step length may be obtained according to the hyperbolic conservation law stable condition, and the preset end duration may be set.
In the scheme, various simulation parameters used in the shock wave simulation process are set, so that the normal operation of shock wave simulation calculation is ensured.
In some embodiments, fig. 7 is a schematic flowchart of another shockwave simulation method provided by an embodiment of the present disclosure, and as shown in fig. 7, the shockwave simulation method further includes:
step 701, when the boundary area of the calculation area does not include the shockwave grid, determining the second physical parameter as a second virtual physical parameter corresponding to the virtual grid.
The second virtual physical parameter may be a physical parameter corresponding to the virtual grid determined according to the second physical parameter.
In this embodiment, when none of the grids in the calculation region is the shockwave grid, a boundary closest to the virtual grid may be determined as an adjacent boundary, and a second physical parameter corresponding to a grid in the calculation region that is located in the same row as the virtual grid and adjacent to the adjacent boundary may be determined as a second physical parameter of the virtual grid.
In some embodiments, an adaptive transition condition may be set, which may be:
Figure 861449DEST_PATH_IMAGE050
specifically, as shown in fig. 3, s represents the number of the shockwave meshes, t represents the number of meshes in the x, y, or z direction, N represents the boundary region as N layers of meshes adjacent to the boundary in the calculation region, and N may be a constant related to a numerical discrete format. When the shock wave grids are located in the Boundary region (namely s =1~N or s = t + 1-N-t), determining a Boundary Riemann problem based on the Boundary of the calculation region, and constructing a Riemann-Reflection Boundary Condition (NRBC) to solve the shock wave simulation result of the next time step; when the shock wave grid is located outside the Boundary region (i.e. s ≠ 1~N or s ≠ t + 1-N-t), determining the second physical parameter as a second virtual physical parameter corresponding to the virtual grid, and constructing an Inflow/Outflow Boundary Condition (IOBC).
When the shockwave mesh is outside the boundary region, the computational domain boundary is processed using the outflow/outflow boundary conditions, as shown in fig. 8. Specifically, the second physical parameter is determined as a second virtual physical parameter corresponding to the virtual grid, as shown below, the second virtual physical parameter of the virtual grid
Figure 223160DEST_PATH_IMAGE051
Comprises the following steps:
Figure 934764DEST_PATH_IMAGE052
wherein j represents the mesh number, t represents the mesh number, N represents the boundary area as the N layers of meshes adjacent to the boundary in the calculation area,
Figure 331110DEST_PATH_IMAGE053
a physical parameter representing the current time step of the leftmost mesh adjacent to the boundary,
Figure 735547DEST_PATH_IMAGE054
a physical parameter representing the current time step of the rightmost mesh adjacent to the boundary,
Figure 584554DEST_PATH_IMAGE055
representing the corresponding physical parameters of the left virtual grid,
Figure 99849DEST_PATH_IMAGE056
representing the corresponding physical parameters of the right virtual grid.
And 702, determining a discrete solution result of a preset fluid control equation according to the second virtual physical parameter and the intermediate physical parameter of the current time step of the grid in the calculation area, and determining a shock wave simulation result of the next time step according to the discrete solution result.
In this embodiment, after determining the second virtual physical parameter of the virtual grid, the discrete solution of the time and space dimensions may be performed on the preset fluid control equation according to the second virtual physical parameter and the intermediate physical parameter of the current time step to obtain a discrete solution result, the intermediate physical parameter of the grid in the next time step in the calculation region may be obtained by calculation according to the discrete solution result, and the intermediate physical parameter of the next time step may be used as a shock wave simulation result of the next time step
In some embodiments, determining a discrete solution result of the preset fluid control equation according to the second virtual physical parameter and an intermediate physical parameter of the current time step of the grid in the calculation region, and determining a shock wave simulation result of the next time step according to the discrete solution result includes: decoupling multiple substances which are contacted with each other in the calculation area; solving a time dimension and a space dimension of a preset fluid control equation based on the decoupled multiple substances, the intermediate physical parameters and the second virtual physical parameters to obtain a discrete equation corresponding to the preset fluid control equation; and determining an intermediate physical parameter of the next time step of the grid in the calculation region based on the discrete equation. The specific calculation method is similar to the above-mentioned shock wave simulation result based on the first virtual physical parameter and the intermediate physical parameter to determine the next time step, and is not described here again.
In the above scheme, when the boundary area of the calculation area does not include the shock wave grid, that is, the distance between the shock wave and the boundary of the calculation area is long, the second physical parameter is determined as the second virtual physical parameter corresponding to the virtual grid, and the inflow and outflow boundary condition is adopted, so that the effective inflow and outflow of the material at the boundary are ensured. When the boundary region of the calculation region comprises the shock wave grids, namely the shock wave is close to the boundary of the calculation region, the Riemann problem is constructed at the boundary of the calculation region, so that the non-physical numerical reflection of the shock wave which is inconsistent with the actual condition at the boundary of the calculation region is avoided, and the accuracy of shock wave simulation and the matching degree with the actual condition are improved.
Next, the method for simulating a shock wave in the embodiment of the present disclosure is further described with a specific example. Exemplarily, fig. 9 is a schematic flowchart of a further shockwave simulation method provided by an embodiment of the present disclosure, and as shown in fig. 9, the method includes:
and step 901, CAE numerical simulation initialization. The initialization process includes: establishing a calculation area, meshing, setting multiple substances, defining initial conditions, selecting a fluid control equation and solving an event step length.
Step 902, determining a shockwave grid according to a discontinuity identification method.
Step 903, judging whether the shock wave grids are in the boundary area. If so, go to step 904, otherwise, go to step 907.
At step 904, a boundary Riemann problem is constructed at the boundary of the computation region.
Step 905, solving the boundary riemann problem by using a riemann problem solver.
And step 906, determining the virtual physical parameters of the virtual grid according to the solving result.
Step 907, determine the virtual physical parameters of the virtual grid according to the second physical parameters.
At step 908, the multiple substances within the computing area that are in contact with each other are processed.
In step 909, the fluid control equations are solved for spatial dispersion and time marching.
Step 910, determining the intermediate physical parameter of the next time step of the grid in the calculation region.
And 911, judging whether the current time length is greater than a preset ending time length. If yes, go to step 912; otherwise, step 903 is executed.
And 912, outputting the intermediate physical parameters of the next time step as target shock wave simulation results.
Next, the shockwave simulation method in the embodiment of the present disclosure is further described by another specific example, in which an application scenario of the shockwave simulation method is to simulate explosive shockwaves of an explosive in water.
Firstly, initializing CAE numerical simulation according to the problem of propagation of underwater explosion shock waves of explosives. Establishing an Euler coordinate system, adopting a full-size three-dimensional model, positioning the explosive at the central position of a calculation domain, determining an initial interface by adopting a VOF interface capturing method and a Youngs interface reconstruction method, and determining initial physical parameters according to the characteristics of different substances, wherein the mass of the explosive is 137kg, and the balance is water. As shown in fig. 10, the length and width of the calculation region are 18.3 meters, the height is 30 meters, the number of grids is 915 × 915 × 1500, the size of the grids is 2 centimeters, the calculation region runs on a supercomputer system, an information transfer interface (MPI) parallel method is adopted, and the number of Central Processing Units (CPUs) is 1215. According to the underwater explosion characteristic, an Euler equation in a differential form is adopted as a preset fluid control equation to be solved, the time step length is obtained, and the calculation termination time is set to be 14 milliseconds.
Furthermore, the shock wave is identified by adopting a multi-resolution intermittent indicator, and the position of the shock wave in the calculation area is determined and marked based on the approximation error, namely the shock wave grid is determined.
Further, whether the boundary area comprises the shock wave grid or not is judged, and if yes, a Riemann non-reflection boundary condition is set according to the boundary Riemann problem; if not, an inflow/outflow boundary condition is set. Thereby achieving adaptive coupling of the reflection-free boundary conditions with the inflow/outflow boundary conditions. Specifically, when the shock wave is far away from the boundary, namely the distance between the shock wave and the boundary exceeds 3 grids, the outflow/outflow boundary condition is adopted, and the second physical parameter is set as the virtual physical parameter of the virtual grid outside the boundary, so that the effective inflow/outflow of the material is realized. When the shock wave is close to the boundary, namely the distance between the shock wave and the boundary is less than 3 grids, a Riemann non-reflection boundary condition is adopted, a boundary Riemann problem is constructed on the boundary, an HLLC approximate Riemann solving method is adopted, and the following solving formula is obtained according to the Rankine-Yugonniu relation and the continuous relation between the contact interruption speed and the pressure:
Figure 350702DEST_PATH_IMAGE057
where ρ denotes density, u denotes velocity, p denotes pressure, E denotes specific energy, subscript L denotes the left side of the riemann problem, R denotes the right side of the riemann problem,. Denotes a physical quantity after impact, for example,
Figure 926040DEST_PATH_IMAGE058
representing the velocity to the left of the riemann problem after the shock wave,
Figure 262343DEST_PATH_IMAGE059
Figure 581329DEST_PATH_IMAGE060
the values are as follows:
Figure 686688DEST_PATH_IMAGE061
wherein c represents the material sound velocity, the boundary Riemann problem is solved based on the solving formula, and the solving result corresponding to the physical parameters of the grid at the infinite distance outside the calculated area in the solving result is determined
Figure 432927DEST_PATH_IMAGE062
According to
Figure 990948DEST_PATH_IMAGE063
And determining virtual physical parameters of the virtual grid outside the boundary, and finishing effective processing on the non-reflection boundary, thereby effectively inhibiting the numerical problem of non-physical reflection of the shock wave at the boundary.
In the present embodiment, as shown in fig. 10, in order to better compare the effects of different boundary conditions, the front boundary of the calculation region may be set as a boundary based on the boundary riemann problem, the rear boundary may be set as a wall boundary, the left boundary may be set as a far-field boundary, the right boundary may be set as an inflow/outflow boundary, and pressure monitoring points may be set at corresponding positions to capture pressure changes at different times at the positions.
The method for setting the wall surface boundary is as follows: rigid floors, walls in practical physical problems can be set as wall boundaries:
Figure 113625DEST_PATH_IMAGE064
where u denotes velocity, W denotes a physical parameter other than velocity, subscript n denotes a normal direction, subscript g denotes a virtual grid, and subscript j denotes a grid number. Based on the above formula, for the velocity, the normal velocity of the nth layer virtual grid node outside the calculation region boundary is equal to the reverse direction of the normal velocity of the nth layer grid node inside the boundary, and for other physical quantities, the other physical quantities of the nth layer virtual grid node outside the calculation region boundary are equal to the other physical quantities of the nth layer grid node inside the boundary.
In the embodiment, the underwater explosion problem is a typical multi-substance strong interaction problem, the multi-substance problem is decoupled into a single-substance problem by using a Mixture theory method to solve, the multi-substance interface physical state problem is effectively processed, and non-physical oscillation of the multi-substance interface is avoided. According to the effective processing of the boundary and the multi-material interface in the steps, a 5-order WENO finite difference format and a 3-order TVD Runge-Kutta format are adopted to carry out space-time discrete solution on a control equation, so that the physical state of the next time step is obtained, and the multi-material interface is promoted.
And finally, obtaining a numerical simulation result of the interaction between the explosive blast wave and the protective structure.
Based on the analysis of the simulation results in fig. 11 and fig. 12, fig. 11 is a pressure cloud graph of the blast shock wave flow field provided by the embodiment of the present disclosure, where fig. 11 may be a pressure cloud graph 11.2 milliseconds after the explosion occurs. Fig. 12 is a schematic diagram of a pressure-time curve provided by an embodiment of the present disclosure, where fig. 12 may be a schematic diagram of a pressure-time curve at a distance of 18 meters from an explosion center. As shown in fig. 11, the results of numerical simulation of different boundary condition types have large differences after the shock wave reaches the boundary. The pressure is high because the wall boundary conditions directly reflect the shock waves. The inflow/outflow boundary has a poor effect on shock wave treatment, and generates obvious non-physical reflection in the area B of fig. 11, the reflection pressure is large, and also forms obvious reflection in other positions, and an obvious non-physical secondary pressure peak exists in the area D of fig. 12, which has a large influence on the simulation result, and the error is large. The far field boundary conditions are better for shock wave processing than the inflow/outflow boundary conditions, and only in region a of fig. 11, a more pronounced non-physical reflection occurs, and the reflected pressure is also small, no pronounced secondary pressure peak occurs, but in region D of fig. 12, a more pronounced non-physical phenomenon of pressure decay interruption also occurs. The boundary processing method based on the boundary Riemann problem does not generate obvious non-physical reflection of the shock wave in the same area, and the shock wave is better in outgoing absorption effect. As shown in fig. 12, the similarity between the curve based on the boundary riemann problem and the experimentally determined curve is high, and the overall shock wave processing effect is superior to that of the far-field boundary condition. Compared with a simulation result and a test result which adopt the boundary processing method, the peak value error of the pressure curve is smaller, the pressure attenuation trend of the shock wave is approximate, and the precision and the effectiveness of the boundary processing based on the boundary Riemann problem are proved.
The shock wave simulation method provided by the embodiment has the following advantages:
1. according to the embodiment, the self-adaptive coupling of the non-reflection boundary and the inflow/outflow boundary based on the boundary Riemannian problem is realized according to whether the boundary region comprises the shock wave grid or not, the advantages of two boundary conditions are combined, the discontinuous transmission absorption of the shock wave of the large-scale CAE simulation calculation domain boundary and the inflow/outflow of the material can be effectively processed, and the efficiency and the precision of boundary processing are improved.
2. The boundary Riemann problem is established at the boundary of the calculation region, the impact wave discontinuity and the boundary interaction can be efficiently and accurately predicted, the Riemann non-reflection boundary condition is realized, the boundary impact wave number value reflection is restrained, the non-reflection boundary processing precision is improved, meanwhile, the increase of calculation amount caused by the increase of the solution region is avoided, and the calculation efficiency is improved.
3. The embodiment is used as a boundary processing method of a calculation area, can be coupled and connected with a high-precision simulation calculation format of a preset fluid control equation, and can reduce CAE simulation errors and improve the goodness of fit of CAE simulation and test based on the efficient and accurate processing of the boundary.
4. The method can automatically select a boundary processing method based on the Riemann boundary problem or an inflow/outflow boundary processing method based on shock wave identification, and realizes high-precision simulation and prediction of large-scale airspace/water area shock wave propagation and interaction processes by combining with a high-precision calculation format of a preset fluid control equation, so that the engineering technical problem related to the infinite/semi-infinite complex scene high-precision CAE simulation field is solved.
Fig. 13 is a schematic structural diagram of a shock wave simulation apparatus provided in an embodiment of the present disclosure, where the apparatus 1300 may be implemented by software and/or hardware, and may be generally integrated in an electronic device. As shown in fig. 13, the apparatus includes:
a first determining module 1301, configured to determine a shockwave grid in the current time step calculation region;
a constructing module 1302, configured to, when the boundary area of the calculation area includes the shockwave mesh, construct a boundary riemann problem corresponding to the boundary of the calculation area according to the first physical parameter and the second physical parameter; the first physical parameter is a physical parameter of a grid at infinity outside the calculation region, and the second physical parameter is a physical parameter of a grid current time step adjacent to the boundary in the calculation region;
a solving module 1303, configured to determine a solving result of the boundary riemann problem, and determine, according to the solving result, a first virtual physical parameter corresponding to the virtual grid outside the calculation area;
a second determining module 1304, configured to determine a discrete solution result of a preset fluid control equation according to the first virtual physical parameter of the virtual grid and the intermediate physical parameter of the current time step of the grid in the calculation region, and determine a shock wave simulation result of the next time step according to the discrete solution result.
In an alternative embodiment, the apparatus further comprises:
an initialization module to initialize simulation parameters prior to determining a shockwave mesh within a current time step calculation region, the simulation parameters including: calculating at least one of a region parameter, a grid parameter, a multi-material interface parameter, an initial condition parameter, a fluid control equation parameter, a parallel calculation parameter, and a time parameter; wherein the initial condition parameters comprise initial physical parameters of the grids in the calculation region.
In an optional implementation manner, the first determining module 1301 is configured to:
sequentially determining each grid in the calculation area as a grid to be processed;
determining an approximate physical parameter of the grid to be processed according to the grid to be processed and an intermediate physical parameter of a current time step of a grid before the grid to be processed;
if the absolute value of the difference value between the intermediate physical parameter of the current time step of the grid to be processed and the approximate physical parameter is greater than or equal to a preset threshold value, determining the grid to be processed as a shock wave grid; wherein the preset threshold is the product of a multi-resolution parameter and a grid size parameter.
In an alternative embodiment, the boundary region comprises: and calculating a preset number of layers of grids adjacent to the boundary in the area.
In an alternative embodiment, the building module 1302 is configured to:
taking the first physical parameter as a first Riemann parameter value corresponding to a first side in the boundary Riemann problem;
taking the second physical parameter as a second Riemann parameter value corresponding to a second side in the boundary Riemann problem;
wherein if the first side is the left side in the boundary Riemann problem, the second side is the right side in the boundary Riemann problem; if the first side is the right side in the boundary Riemannian problem, the second side is the left side in the boundary Riemannian problem.
In an optional embodiment, the solving module 1303 is configured to:
solving the boundary Riemann problem by adopting a Riemann problem solver to obtain a solving result; wherein the Riemann problem solver comprises: the Rankine-Rainbow equation, and the contact break speed and pressure continuity relationship.
In an optional embodiment, the solving module 1303 is configured to:
and extracting the computational physical parameters corresponding to the first physical parameters in the solving result, and determining the computational physical parameters as the first virtual physical parameters corresponding to the virtual grid.
In an optional implementation, the second determining module 1304 is configured to:
decoupling multiple substances in contact with each other within the calculation region;
solving the time dimension and the space dimension of the preset fluid control equation based on the decoupled multiple substances, the intermediate physical parameters and the first virtual physical parameters to obtain a discrete equation corresponding to the preset fluid control equation;
and determining intermediate physical parameters of the next time step of the grids in the calculation region based on the discrete equation.
In an alternative embodiment, the apparatus further comprises:
a third determining module, configured to determine, when the boundary area of the calculation area does not include the shockwave grid, the second physical parameter as a second virtual physical parameter corresponding to the virtual grid;
and the fourth determination module is used for determining a discrete solution result of a preset fluid control equation according to the second virtual physical parameter and the intermediate physical parameter of the current time step of the grid in the calculation area, and determining a shock wave simulation result of the next time step according to the discrete solution result.
In an alternative embodiment, the apparatus further comprises:
and the fifth determining module is used for determining the shock wave simulation result of the next time step as a target simulation result if the current time corresponding to the current time step is greater than the preset ending time.
The shock wave simulation device provided by the embodiment of the disclosure can execute the shock wave simulation method provided by any embodiment of the disclosure, and has corresponding functional modules and beneficial effects of the execution method.
Fig. 14 is a schematic structural diagram of an electronic device according to an embodiment of the present disclosure. As shown in fig. 14, the electronic device 1400 includes one or more processors 1401 and memory 1402.
The processor 1401 may be a Central Processing Unit (CPU) or other form of processing unit having data processing capabilities and/or instruction execution capabilities, and may control other components in the electronic device 1400 to perform desired functions.
Memory 1402 may include one or more computer program products that may include various forms of computer-readable storage media, such as volatile memory and/or non-volatile memory. The volatile memory may include, for example, random Access Memory (RAM), cache memory (cache), and/or the like. The non-volatile memory may include, for example, read Only Memory (ROM), hard disk, flash memory, etc. One or more computer program instructions may be stored on the computer-readable storage medium and executed by processor 1401 to implement the shockwave simulation methods of the embodiments of the present disclosure described above and/or other desired functions. Various contents such as an input signal, a signal component, a noise component, etc. may also be stored in the computer-readable storage medium.
In one example, the electronic device 1400 may further include: an input device 1403 and an output device 1404, which are interconnected by a bus system and/or other form of connection mechanism (not shown).
The input device 1403 may also include, for example, a keyboard, a mouse, and the like.
The output device 1404 may output various information including the determined distance information, direction information, and the like to the outside. The output devices 1404 may include, for example, a display, speakers, a printer, and a communication network and remote output devices connected thereto, among others.
Of course, for simplicity, only some of the components of the electronic device 1400 relevant to the present disclosure are shown in fig. 14, omitting components such as buses, input/output interfaces, and the like. In addition, electronic device 1400 may include any other suitable components, depending on the particular application.
In addition to the above-described methods and apparatus, embodiments of the present disclosure may also be a computer program product comprising computer program instructions that, when executed by a processor, cause the processor to perform the shockwave simulation methods provided by embodiments of the present disclosure.
The computer program product may write program code for carrying out operations for embodiments of the present disclosure in any combination of one or more programming languages, including an object oriented programming language such as Java, C + + or the like and conventional procedural programming languages, such as the "C" programming language or similar programming languages. The program code may execute entirely on the user's computing device, partly on the user's device, as a stand-alone software package, partly on the user's computing device and partly on a remote computing device, or entirely on the remote computing device or server.
Furthermore, embodiments of the present disclosure may also be a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, cause the processor to perform the shockwave simulation method provided by embodiments of the present disclosure.
The computer-readable storage medium may take any combination of one or more readable media. The readable medium may be a readable signal medium or a readable storage medium. A readable storage medium may include, for example, but not limited to, an electronic, magnetic, optical, electromagnetic, infrared, or semiconductor system, apparatus, or device, or a combination of any of the foregoing. More specific examples (a non-exhaustive list) of the readable storage medium include: an electrical connection having one or more wires, a portable disk, a hard disk, a Random Access Memory (RAM), a read-only memory (ROM), an erasable programmable read-only memory (EPROM or flash memory), an optical fiber, a portable compact disc read-only memory (CD-ROM), an optical storage device, a magnetic storage device, or any suitable combination of the foregoing.
It is noted that, in this document, relational terms such as "first" and "second," and the like, may be used solely to distinguish one entity or action from another entity or action without necessarily requiring or implying any actual such relationship or order between such entities or actions. Also, the terms "comprises," "comprising," or any other variation thereof, are intended to cover a non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements does not include only those elements but may include other elements not expressly listed or inherent to such process, method, article, or apparatus. Without further limitation, an element defined by the phrase "comprising a … …" does not exclude the presence of another identical element in a process, method, article, or apparatus that comprises the element.
The foregoing are merely exemplary embodiments of the present disclosure, which enable those skilled in the art to understand or practice the present disclosure. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the generic principles defined herein may be applied to other embodiments without departing from the spirit or scope of the disclosure. Thus, the present disclosure is not intended to be limited to the embodiments shown herein but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.

Claims (9)

1. A method for simulating a shock wave, comprising:
determining a shock wave grid in a current time step calculation region;
when the boundary area of the calculation area comprises the shock wave grid, constructing a boundary Riemann problem corresponding to the boundary of the calculation area according to the first physical parameter and the second physical parameter; the first physical parameter is a physical parameter of a grid at infinity outside the calculation area, and the second physical parameter is a physical parameter of a grid current time step adjacent to the boundary in the calculation area;
determining a solving result of the boundary Riemann problem, and determining a first virtual physical parameter corresponding to the virtual grid outside the calculation region according to the solving result;
determining a discrete solution result of a preset fluid control equation according to the first virtual physical parameter of the virtual grid and an intermediate physical parameter of the current time step of the grid in the calculation region, and determining a shock wave simulation result of the next time step according to the discrete solution result;
when the boundary area of the calculation area does not comprise the shock wave grids, determining the second physical parameters as second virtual physical parameters corresponding to the virtual grids;
and determining a discrete solution result of a preset fluid control equation according to the second virtual physical parameter and the intermediate physical parameter of the current time step of the grid in the calculation region, and determining a shock wave simulation result of the next time step according to the discrete solution result.
2. The method of claim 1, wherein determining the shockwave mesh within the current time step calculation region comprises:
sequentially determining each grid in the calculation area as a grid to be processed;
determining an approximate physical parameter of the grid to be processed according to the grid to be processed and an intermediate physical parameter of a current time step of a grid before the grid to be processed;
if the absolute value of the difference value between the intermediate physical parameter of the current time step of the grid to be processed and the approximate physical parameter is greater than or equal to a preset threshold value, determining the grid to be processed as a shock wave grid; wherein the preset threshold is the product of a multi-resolution parameter and a grid size parameter.
3. The method according to claim 1, wherein constructing a boundary Riemann problem corresponding to the boundary of the calculation region according to the first physical parameter and the second physical parameter comprises:
taking the first physical parameter as a first Riemann parameter value corresponding to a first side in the boundary Riemann problem;
taking the second physical parameter as a second Riemann parameter value corresponding to a second side in the boundary Riemann problem;
wherein if the first side is the left side in the boundary Riemannian problem, the second side is the right side in the boundary Riemannian problem; if the first side is the right side in the boundary Riemann problem, the second side is the left side in the boundary Riemann problem.
4. The method of claim 1, wherein determining the solution to the boundary Riemannian problem comprises:
solving the boundary Riemann problem by adopting a Riemann problem solver to obtain a solving result; wherein the Riemann problem solver comprises: the Rankine-Rainbow equation, and the contact break speed and pressure continuity relationship.
5. The method according to claim 1, wherein the determining the first virtual physical parameter corresponding to the virtual grid outside the calculation region according to the solution result comprises:
and extracting the computational physical parameters corresponding to the first physical parameters in the solving result, and determining the computational physical parameters as the first virtual physical parameters corresponding to the virtual grid.
6. The method of claim 1, wherein the determining a discrete solution result of a preset fluid control equation according to the first virtual physical parameter of the virtual grid and an intermediate physical parameter of a current time step of the grid in the calculation region, and determining a shock wave simulation result of a next time step according to the discrete solution result comprises:
decoupling multiple substances in contact with each other in the calculation area;
solving the time dimension and the space dimension of the preset fluid control equation based on the decoupled multiple substances, the intermediate physical parameter and the first virtual physical parameter to obtain a discrete equation corresponding to the preset fluid control equation;
and determining an intermediate physical parameter of the next time step of the grid in the calculation region based on the discrete equation.
7. A shock wave simulation apparatus, comprising:
the first determination module is used for determining a shock wave grid in a current time step calculation region;
the construction module is used for constructing a boundary Riemannian problem corresponding to the boundary of the calculation region according to the first physical parameter and the second physical parameter when the boundary region of the calculation region comprises the shock wave grid; the first physical parameter is a physical parameter of a grid at infinity outside the calculation region, and the second physical parameter is a physical parameter of a grid current time step adjacent to the boundary in the calculation region;
the solving module is used for determining a solving result of the boundary Riemannian problem and determining a first virtual physical parameter corresponding to the virtual grid outside the calculation area according to the solving result;
the second determination module is used for determining a discrete solution result of a preset fluid control equation according to the first virtual physical parameter of the virtual grid and an intermediate physical parameter of the grid in the calculation area at the current time step, and determining a shock wave simulation result of the next time step according to the discrete solution result;
a third determining module, configured to determine, when the boundary area of the calculation area does not include the shockwave grid, the second physical parameter as a second virtual physical parameter corresponding to the virtual grid;
and the fourth determination module is used for determining a discrete solution result of a preset fluid control equation according to the second virtual physical parameter and the intermediate physical parameter of the current time step of the grid in the calculation area, and determining a shock wave simulation result of the next time step according to the discrete solution result.
8. An electronic device, characterized in that the electronic device comprises:
a processor;
a memory for storing the processor-executable instructions;
the processor is used for reading the executable instructions from the memory and executing the instructions to realize the shock wave simulation method of any one of the claims 1 to 6.
9. A computer-readable storage medium, characterized in that the storage medium stores a computer program for executing the shockwave simulation method of any of the above claims 1-6.
CN202211249477.4A 2022-10-12 2022-10-12 Shock wave simulation method, device, equipment and medium Active CN115329646B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211249477.4A CN115329646B (en) 2022-10-12 2022-10-12 Shock wave simulation method, device, equipment and medium

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211249477.4A CN115329646B (en) 2022-10-12 2022-10-12 Shock wave simulation method, device, equipment and medium

Publications (2)

Publication Number Publication Date
CN115329646A CN115329646A (en) 2022-11-11
CN115329646B true CN115329646B (en) 2023-03-10

Family

ID=83914109

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211249477.4A Active CN115329646B (en) 2022-10-12 2022-10-12 Shock wave simulation method, device, equipment and medium

Country Status (1)

Country Link
CN (1) CN115329646B (en)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115795282B (en) * 2023-01-30 2023-05-09 武汉工程大学 Shock tube dynamic pressure reconstruction method and device, electronic equipment and storage medium
CN116050196B (en) * 2023-04-03 2023-06-30 国家超级计算天津中心 Multi-dimensional simulation method, device, equipment and storage medium
CN116864041A (en) * 2023-06-05 2023-10-10 国家超级计算天津中心 Method and device for determining temperature in material collision, electronic equipment and storage medium

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108985000A (en) * 2018-10-15 2018-12-11 武汉海威船舶与海洋工程科技有限公司 Underwater far field explosion wave emulated computation method based on gradient grid
CN109902376A (en) * 2019-02-25 2019-06-18 北京理工大学 A kind of fluid structurecoupling high resolution numerical simulation method based on Continuum Mechanics
CN111859646A (en) * 2020-07-09 2020-10-30 南京理工大学 Shock wave variable step length solving method based on B spline mapping function material point method
CN112464335A (en) * 2020-11-10 2021-03-09 中国人民解放军63921部队 Visual simulation analysis method for dangerous goods explosion in tall and large space complex building structure
CN113724346A (en) * 2020-05-26 2021-11-30 北京安云世纪科技有限公司 Method and system for realizing two-dimensional shock wave effect, storage medium and computer equipment thereof
CN114169184A (en) * 2021-10-19 2022-03-11 北京理工大学 High-precision self-adaptive finite volume-finite difference coupling numerical simulation method
CN114611340A (en) * 2022-03-17 2022-06-10 北京理工大学 Coupling Euler-Lagrange method for accurately capturing shock wave propagation process

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7225173B2 (en) * 2002-09-09 2007-05-29 Carmel - Haifa University Economic Corporation Ltd. Apparatus and method for efficient adaptation of finite element meshes for numerical solutions of partial differential equations
CN104699892B (en) * 2015-01-22 2017-08-01 三峡大学 Study of Landslides swell propagation rule and its model predicted Dam life and method
CN110852005B (en) * 2019-10-21 2021-06-15 北京理工大学 Numerical simulation method for self-adaptive expansion of computational domain of large-scale parallel computation
CN113919059A (en) * 2021-09-16 2022-01-11 东风汽车集团股份有限公司 Vehicle body analysis method under action of blast-hole shock wave, terminal device and medium

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108985000A (en) * 2018-10-15 2018-12-11 武汉海威船舶与海洋工程科技有限公司 Underwater far field explosion wave emulated computation method based on gradient grid
CN109902376A (en) * 2019-02-25 2019-06-18 北京理工大学 A kind of fluid structurecoupling high resolution numerical simulation method based on Continuum Mechanics
CN113724346A (en) * 2020-05-26 2021-11-30 北京安云世纪科技有限公司 Method and system for realizing two-dimensional shock wave effect, storage medium and computer equipment thereof
CN111859646A (en) * 2020-07-09 2020-10-30 南京理工大学 Shock wave variable step length solving method based on B spline mapping function material point method
CN112464335A (en) * 2020-11-10 2021-03-09 中国人民解放军63921部队 Visual simulation analysis method for dangerous goods explosion in tall and large space complex building structure
CN114169184A (en) * 2021-10-19 2022-03-11 北京理工大学 High-precision self-adaptive finite volume-finite difference coupling numerical simulation method
CN114611340A (en) * 2022-03-17 2022-06-10 北京理工大学 Coupling Euler-Lagrange method for accurately capturing shock wave propagation process

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
A Framework of Runge–Kutta, Discontinuous Galerkin, Level Set and Direct Ghost Fluid Methods for the Multi-Dimensional Simulation of Underwater Explosions;Nan Si et al.;《fluids》;20211229;全文 *
基于欧拉有限元方法的柱形药包水下爆炸冲击波载荷特性研究;许流逸等;《第十六届全国水动力学学术会议暨第三十二届全国水动力学研讨会论文集》;20211031;第672-678页 *
黎曼边界条件在高阶精度非结构有限体积方法中的验证与应用;孔令发等;《空气动力学学报》;20210630;第39卷(第3期);第21-30页 *

Also Published As

Publication number Publication date
CN115329646A (en) 2022-11-11

Similar Documents

Publication Publication Date Title
CN115329646B (en) Shock wave simulation method, device, equipment and medium
Wang et al. A computational framework for the simulation of high‐speed multi‐material fluid–structure interaction problems with dynamic fracture
Liu et al. Modeling and simulation of ice–water interactions by coupling peridynamics with updated Lagrangian particle hydrodynamics
Wang et al. Computational algorithms for tracking dynamic fluid–structure interfaces in embedded boundary methods
JP5280863B2 (en) Piecewise meshing and neighborhood search for object interaction simulation
US8959005B2 (en) Building envelope determination
Bonfiglioli et al. Unsteady shock‐fitting for unstructured grids
Ge et al. Investigation of underwater explosion near composite structures using a combined RKDG-FEM approach
Lu et al. Runge–Kutta discontinuous Galerkin method with front tracking method for solving the compressible two-medium flow
Yu et al. Discussion of assumptions behind the external dynamic models in ship collisions and groundings
Alshammari et al. Mitigation of blast loading through blast–obstacle interaction
Puscas et al. A conservative embedded boundary method for an inviscid compressible flow coupled with a fragmenting structure
Talley et al. Coalescence prevention algorithm for level set method
CN113255245A (en) Method and device for predicting downhole perforation shock wave and readable storage medium
JP4877744B2 (en) Building sound simulation system
Prescott et al. Incorporating dynamic 3D simulation into PRA
Abdullahi et al. Simulation and detection of blockage in a pipe under mean fluid flow using acoustic wave propagation technique
Doebling The escape of high explosive products: an exact-solution problem for verification of hydrodynamics codes
Wang et al. Verification and validation of a detonation computational fluid dynamics model
AU2007208112B2 (en) Termination assessment of a computer simulation
Devilliers et al. A new adaptive mesh refinement to model water flow around fishing nets
Icardi et al. Free Vibration of flexible soft-core sandwiches according to layerwise theories differently accounting for the transverse normal deformability
Lu et al. Assessment of a new mesh generation and modeling approach for the surface ship far-field underwater explosion problem
Lu Computationally-effective Modeling of Far-field Underwater Explosion for Early-stage Surface Ship Design
McDonald Development of a high-order finite-volume method for unstructured meshes

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
TR01 Transfer of patent right
TR01 Transfer of patent right

Effective date of registration: 20230724

Address after: 300457 5102, building 5, No. 19, Xinhuan West Road, economic and Technological Development Zone, Binhai New Area, Tianjin

Patentee after: NATIONAL SUPERCOMPUTER CENTER IN TIANJIN

Patentee after: National University of Defense Technology

Patentee after: Haihe Laboratory of Advanced Computing and Key Software (Xinchuang)

Address before: 300457 5102, building 5, No. 19, Xinhuan West Road, economic and Technological Development Zone, Binhai New Area, Tianjin

Patentee before: NATIONAL SUPERCOMPUTER CENTER IN TIANJIN

Patentee before: Haihe Laboratory of Advanced Computing and Key Software (Xinchuang)