CN110852005A - Numerical simulation method for self-adaptive expansion of computational domain of large-scale parallel computation - Google Patents
Numerical simulation method for self-adaptive expansion of computational domain of large-scale parallel computation Download PDFInfo
- Publication number
- CN110852005A CN110852005A CN201911001549.1A CN201911001549A CN110852005A CN 110852005 A CN110852005 A CN 110852005A CN 201911001549 A CN201911001549 A CN 201911001549A CN 110852005 A CN110852005 A CN 110852005A
- Authority
- CN
- China
- Prior art keywords
- calculation
- domain
- grid
- level
- pos
- 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.)
- Granted
Links
Images
Landscapes
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
The invention relates to a numerical simulation method of a large-scale parallel computing self-adaptive expanding computing domain, belonging to the field of numerical simulation. The method is an adaptive grid amplification method based on an algorithm of an essential non-oscillation finite difference method, a level set method and a real virtual fluid method, and can greatly improve the calculation efficiency of large-scale problems. The WENO finite difference format converts flux solution in the Euler equation into discrete solution by using the weight of the discontinuity as zero, and reduces numerical value oscillation under the condition of ensuring precision; the Level-Set method can implicitly track the position evolution of a complex multi-medium interface; the RGFM method utilizes a virtual flow field to process the interface state, and maintains the continuous characteristics of pressure and speed of the interface; the self-adaptive grid amplification method can perform high-precision calculation on a large-range flow field of a meter and kilometer scale under the condition of certain calculation resources, and remarkably improves the calculation efficiency. By the method, the explosion flow field propagation process can be simulated efficiently and accurately.
Description
Technical Field
The invention relates to a numerical simulation method of a large-scale parallel computing self-adaptive expanding computing domain, belonging to the field of numerical simulation.
Background
The explosion damage experiment test can quickly and accurately obtain the shock wave pressure value, but the shock wave speed is thousands of meters per second, the near-field pressure is generally higher than ten atmospheric pressures, and the requirement on the measurement technology of a sensor, an oscilloscope and the like is higher. Foreign scholars measure and calibrate pressure values at different dosages and different positions by using experimental means, and fit widely-used Henrych empirical formulas (J.Henrych, The Dynamics of expansion and its Use, Elsevier, New York,1979), Mills empirical formulas (C.Mills,1st int. Conf. Hazard procedure. Edinburgh, UK,1987), Brode empirical formulas (H.Brode, Phys. fluids,1959,2: 217-; the empirical formula has coefficients determined by experiments, and the coefficients obtained by different researchers are different and are not suitable for the actual engineering problem; for large explosive quantities or different boundary conditions, the actual conditions need to be considered in selecting the formula. In addition, the experimental tests have the disadvantages that: (1) it is difficult for the measuring instrument to obtain a correct response at high speed and high voltage. (2) The near-field pressure is extremely high, and a high-temperature fireball is generated, so that the sensor has high requirements on the performances of impact resistance and high temperature resistance and is easy to damage. (3) The error of the experiment carried out for many times is large, and various uncertain factors such as man-made factors, environment factors and the like exist.
With the development of numerical simulation technology and the advent of high-performance computers, numerical simulation has become a main technical means for explosion damage research. The technology has the advantages that: (1) compared with experimental research, the numerical simulation is very safe and low in cost. (2) The efficiency of researching the explosion damage mechanism is high, and the operation is simple. The explosion damage problem relates to high speed, high pressure, high temperature, phase change, mutual coupling and even mixing of various media such as gas, liquid and solid, and the material not only can be seriously deformed and broken, but also can be melted and even vaporized. Foreign scholars have conducted numerical simulation studies on explosions in pipelines, for example, Yanez et al (J.Yanez, International Journal of Hydrogen Energy,2011,36: 2613-; a Heidari (A.Heidari, International journal of Hydrogen Energy,2011,36: 2538-. Therefore, when large-scale high-precision parallel computation is performed, a numerical simulation result cannot be obtained efficiently and accurately. The invention provides a numerical simulation method of a large-scale parallel computing self-adaptive enlarged computing domain, which can quickly solve a large-scale explosion flow field under the condition of certain computing resources, and the three-dimensional solving speed is far higher than that of a common method. In summary, the method provided by the invention has important technical background.
Disclosure of Invention
The invention aims to provide a numerical simulation method of a self-adaptive expanded calculation domain of large-scale parallel calculation; the method meets the requirement of high precision under the condition of certain computing resources, can efficiently simulate the large-scale explosion flow field problem, obviously improves the efficiency, and further solves the engineering problem in the technical field of large-scale high-precision numerical simulation.
The large-scale high-precision numerical simulation technical field comprises the fields of explosion damage, aerospace and mechanical engineering.
The purpose of the invention is realized by the following technical scheme:
the invention provides a self-adaptive calculation domain expansion method for large-scale high-precision parallel calculation. Specifically, the method is an adaptive grid amplification method based on an algorithm of a finite difference method without shock in essence (WENO), a Level Set method (Level-Set) and a real virtual fluid method (RGFM), and can greatly improve the calculation efficiency of large-scale problems. The WENO finite difference format converts flux solution in the Euler equation into discrete solution by using the weight of the discontinuity as zero, and reduces numerical value oscillation under the condition of ensuring precision; the Level-Set method can implicitly track the position evolution of a complex multi-medium interface; the RGFM method utilizes a virtual flow field to process the interface state, and maintains the continuous characteristics of pressure and speed of the interface; the self-adaptive grid amplification method can perform high-precision calculation on a large-range flow field of a meter and kilometer scale under the condition of certain calculation resources, and remarkably improves the calculation efficiency. By the method, the explosion flow field propagation process can be simulated efficiently and accurately.
A numerical simulation method of a self-adaptive expansion calculation domain of large-scale parallel calculation comprises the following steps:
step 1: according to the explosion flow field scale to be simulated and the calculation resources (CPU core number and memory size), determining the final calculation domain size, the number of parallel partitions in the x direction, the y direction and the z direction, and giving the expansion times and the expansion ratio in the self-adaptive expansion calculation domain method and the grid number in each partition; giving serial numbers of CPUs in a calculation domain; and determining the grid widths in the x, y and z directions when the initial calculation domain l is 0, the sizes of the corresponding calculation domains and the coordinates of each calculation grid point.
Step 1.1: determining final calculated field size [ pi ]x1,πx2]*[πy1,πy2]*[πz1,πz2]The parallel partition numbers cutx, cuty and cutz in the three directions of x, y and z, and the grid number N in each partitionx×Ny×NzIn which N isxDenotes the number of x-direction grids, NyRepresenting the number of grids in the y-direction, NzIndicating the number of grids in the z direction. The total amount of CPU used in parallel is cutx × cuty × cutz, and the flow field is exploded to calculate the central point O (x)O,yO,zO) Expanding in the directions of x-, x +, y-, y +, z-and z +, expanding times in the directions of x, y and z, namely largex, largey and largez, and expanding the ratio r to obtain the total grid number of cutx cut Nx*Ny*Nz。
Step 1.2: the serial number of each CPU in the computational domain is given.
Defining myid as a CPU serial number (starting from 0); posx、posy、poszPos indicating that the CPU is in the x directionxPos in the direction of (y)yPos in the z directionzA CPU. posxIs in the range of 1 to cutx, posyIs in the range of 1 to cut, poszThe value of (A) is in the range of 1 to cutz.
CPU serial numbers myid and posx、posy、poszThe relationship of (1) is:
posx＝mod(myid,cutx)+1
posy＝myid/cutx+1 (1)
posz＝myid/(cutx×cuty)+1
where mod is the remainder symbol.
Step 1.3: when the initial calculation field l is 0, the grid widths of the calculation region in the x, y and z directions are respectively
Calculating the domain size omega when the initial calculation domain l is 00Comprises the following steps:
wherein r islargexIs the largex power of r, rlargeyIs the largey power of r, rlargezIs the power of r to the largez power, where pix1For the final calculation of the leftmost coordinate, π, in the x-direction of the domainx2For the final calculation of the rightmost coordinate, π, in the x-direction of the domainy1For final calculation of the leftmost coordinate, π, in the y-direction of the domainy2For final calculation of the rightmost coordinate, π, in the y-direction of the fieldz1For the final calculation of the leftmost coordinate, π, in the z-direction of the domainz2The domain z-direction rightmost coordinate is finally calculated.
Thus, the x, y, z coordinates (x) of each grid in the domain are initially computedi,yjzk) Comprises the following steps:
wherein i, j, k are the ith grid in the x direction, the jth grid in the y direction and the kth grid in the z direction of the grid in the CPU.
Step 2: defining a level-set function for distinguishing detonation products and air interface positions at the moment when the initial state t is 0 in the initial calculation region; assigning initial values to physical quantities of two regions of detonation products and air; establishing a conservation type Euler control equation set for describing detonation products and air; and establishing an equation of state function of the detonation product and the air.
Step 2.1: defining a level-set function for distinguishing detonation products and air interface positions at the moment when the initial state t is 0 in a minimum (initial) calculation region, giving the position of a multi-material interface at any moment by using a symbol distance function, and tracking the evolution process of the level-set function in real time, wherein the interface gamma (t) at the moment t is positioned in the level-set function phi (x)i,yj,zkT) position of zero:
Γ(t)＝{(xi,yj,zk)|φ(xi,yj,zk,t)＝0} (5)
at the moment when t is 0, the interface position of the detonation product and the air can be obtained by the distance between the position of each calculation grid point in the calculation domain and the interface position:
wherein x, y, z are physical space coordinates of the calculation grid points, d (Γ (0), (x), respectivelyi,yj,zk) Denotes a calculation grid point (x)i,yj,zk) Distance to initial interface Γ (0), SDetonation productsAnd SAir (a)Respectively, regions on both sides of the material interface.
Step 2.2: the initial state of the detonation product and air obtained in step 2.1 is given.
Initial density ρ of a given detonation product1Three directional velocities u1、v1、w1Pressure p1(ii) a Or detonationInitial density of product ρ1Three directional velocities u1、v1、w1Internal energy per unit mass e1；
And initial density ρ of air2Three directional velocities u2、v2、w2Pressure p2(ii) a Or initial density ρ of air2Three directional velocities u2、v2、w2Internal energy per unit mass e2。
Step 2.3: and constructing a system of Euler control equations describing the conservation of detonation products and air.
The conservation-oriented Euler control equation system is as follows:
wherein the content of the first and second substances,
U＝[ρ ρu ρv ρw ρE]T，
F＝[ρu ρu2+p ρuv ρuw u(ρE+p)]T，
G＝[ρv ρuv ρv2+p ρvw v(ρE+p)]T，
H＝[ρw ρuw ρvw ρw2+p w(ρE+p)]T，
and E is the total energy per unit mass,if the detonation product is calculated, the rho, u, v, w, p and e in the formula (7) are used as rho1、u1、v1、w1、p1、e1Substituting; if calculating air, p, u, v, w, p, e in the formula (7) is used as p2、u2、v2、w2、p2、e2And (6) substituting.
Step 2.4: and selecting a state equation function of the material to further determine the physical state of the material under deformation and motion.
Since there are only five equations in equation (6), but there are six unknowns, the state equation needs to be selected to form the whole control equationThe row is closed. Internal energy per unit mass e of detonation product1And pressure p1JWL equation of state relationship (sqj) expressed as:
where ρ is0Is the initial density of the explosive, A, B, R1、R2And ω are the equation of state parameters. The ideal gas equation of state relationship for air is expressed as:
p2＝(γ-1)ρ2e2(9)
where γ is the polytropic gas index.
And step 3: and setting boundary conditions of each stage of area.
Because the computation domains differ in size at different times, the computation domain boundary conditions are different at each level. For the difference of the spatial discrete formats, each computation grid point needs physical quantity computation flux of the surrounding N points, so N layers of virtual points are needed at the boundary of the computation domain.
(1) For an outgoing boundary, at each level of computational domain boundary, there is
① x direction and when posxWhen 1, there are
i＝1,2,...,N,j＝0,1,...,Ny,k＝0,1,...,Nz；
② x direction and when posxWhen n x, there is
i＝1,2,...,N,j＝0,1,...,Ny,k＝0,1,...,Nz；
③ y direction and when posyWhen 1, there are
i＝0,1,...,Nx,j＝1,2,...,N,k＝0,1,...,Nz④ y direction and when posyNy, there are
i＝0,1,...,Nx,j＝1,2,...,N,k＝0,1,...,Nz；
⑤ z direction and when poszWhen 1, there are
i＝0,1,...,Nx,j＝0,1,...,Ny,k＝1,2,...,N；
⑥ z direction and when poszWhen n is nz, there is
i＝0,1,...,Nx,j＝0,1,...,Ny,k＝1,2,...,N
(2) There are a slip boundary and a symmetric boundary for the fixed wall,
① x direction and when posxWhen 1, there are
φ(x-i,yj,zk,t)＝φ(xi,yj,zk,t)，
i＝1,2,...,N,j＝0,1,...,Ny,k＝0,1,...,Nz；
② x direction and when posxWhen n x, there is
φ(xnx+i,yj,zk,t)＝φ(xnx-i,yj,zk,t)，
i＝1,2,...,N,j＝0,1,...,Ny,k＝0,1,...,Nz；
③ y direction and when posyWhen 1, there are
φ(xi,y-j,zk,t)＝φ(xi,yj,zk,t)，
i＝0,1,...,Nx,j＝1,2,...,N,k＝0,1,...,Nz；
④ y direction and when posyNy, there are
φ(xi,yny+j,zk,t)＝φ(xi,yny-j,zk,t)，
i＝0,1,...,Nx,j＝1,2,...,N,k＝0,1,...,Nz；
⑤ z direction and when poszWhen 1, there are
φ(xi,yj,z-k,t)＝φ(xi,yj,zk,t)，
i＝0,1,...,Nx,j＝0,1,...,Ny,k＝1,2,...,N；
⑥ z direction and when poszWhen n is nz, there is
φ(xi,yj,znz+k,t)＝φ(xi,yj,znz-k,t)，
i＝0,1,...,Nx,j＝0,1,...,Ny,k＝1,2,...,N
The outflow boundary comprises, for example, an infinite volume of air;
the fixed wall with a sliding boundary comprises a rigid ground;
the symmetric boundary comprises an 1/2 model, a 1/4 model, or a 1/8 model;
and 4, step 4: and selecting the l-th-level calculation time step.
Selecting the l-th level to calculate the time step delta tl：
Wherein the parameter CFL is a constant between 0 and 1, Δ xlCalculating the grid width in the x direction in the area for the l level; Δ ylCalculating the grid width in the y direction in the area for the l level; Δ zlCalculating the grid width in the z direction in the area for the l level; u, v and w are the speeds of the three directions of the current calculation grid point, and if the current calculation grid point is in the detonation product region, u is taken as1、v1、w1If in the air region, it is taken as u2、v2、w2(ii) a c is the speed of sound of the current calculated grid point.
And 5: and solving a normal vector of a material interface in each Level of region by using a Level-set function, solving a normal Riemann problem at the interface of a detonation product and air, determining the physical state of N layers of grids near the interface by using a real virtual fluid method, and updating the physical state of the full flow field.
Using Level-set equation to advance the interface position in the flow field, using the following formula to advance for each calculation grid point in the calculation domain, and using phi (x) to simplify the formulai,yj,zkT) is abbreviated as φ;
if the calculated grid point is in the detonation product region, u, v and w in the formula (11) are used as u1、v1、w1Substituting; if the calculation grid point is in the air region, u, v, w in the formula (11) are used as u2、v2、w2And (6) substituting.
The solving formula of the normal vector n of each calculation grid point is as follows:
the method comprises the steps of firstly decomposing the speed of each calculation grid point according to a normal vector of the point to obtain a normal speed and a tangential speed, then utilizing physical quantity states of detonation products, air density, normal speed and pressure close to an interface, and utilizing a Riemann solver to give the normal speed and the pressure at the interface and the density of the detonation products and the air at two sides of the interface according to the normal speed continuity, the pressure continuity and the conservation law at two sides of the interface, wherein the normal speed and the tangential speed are physical states under a local coordinate system, the normal speed of the calculation grid points close to the interface needs to be kept unchanged, and the speed component is reversely transformed to obtain the speed under a global coordinate system. Respectively establishing a layer of banded regions with 20 grid widths in detonation products and air regions at two sides of the interface, wherein calculation grid points in the banded regions are virtual points; the density, three directional velocities and pressure at the virtual point are obtained by continuation. The continuation equation is:
where I represents physical quantities at a virtual point, including density, three directional velocity, and pressure. And finally, solving the detonation product and the air respectively by utilizing a high-precision WENO format.
Step 6: and judging whether the current calculation domain needs to be expanded or not.
After a period of time, detonation products in the initial calculation domain expand outwards, namely reach the boundary of the initial calculation domain, and are respectively analyzed by adopting shock waves (physical quantity mutation) and rarefaction waves (physical quantity slow rising or falling). Judging whether expansion is needed; if the expansion is judged not to be needed, continuing to operate the step 3, and if the expansion is judged to be about to be carried out, executing the steps 7-9;
when the shock wave is about to propagate to reach the boundary of the initial calculation domain, the shock wave is located at the position xMIs a critical point; if the critical point satisfies the following threshold condition, the calculation domain is expanded,
|IM-I0|＞Icr(14)
wherein, I can represent density, pressure; i is0Is the constant density, constant pressure, I, of the airMDensity, pressure at critical point; i iscrThe critical threshold value is generally a small number, so that when the shock wave reaches a critical point, the calculation domain is automatically expanded, and important information is prevented from being lost due to the fact that the shock wave flows out of the boundary of the calculation domain.
When the sparse wave is about to propagate to reach the boundary of the initial calculation domain, the position x is locatedNIs a critical point; when this point satisfies the following threshold condition, the calculation domain is about to expand,
wherein I may represent density, pressure; i isNDensity, pressure at critical point; i isc'rThe critical gradient threshold is typically taken to be a small number and is related to the grid width. If the shockwave threshold (as in equation 14) is still used, within a finite computation step, at the boundary point xNA "step" is created, thereby introducing an error. Therefore, by determining whether or not the size is increased using equation (15), it is possible to prevent the loss of important information due to the sparse wave flowing out of the boundary of the calculation domain.
And 7: in the calculation process, detonation products are spread outwards, pressure and density are spread outwards, the boundary of the first-level calculation domain meets the expansion threshold condition in the step 6, and at this time, l' ═ 1 is the expanded calculation domain, and the calculation domain of the expanded calculation domain, the corresponding grid width and the coordinates of each calculation grid point in each calculation domain need to be given.
The computation domain is ranked as l 0, 1., maxlage, where maxlage is the maximum of largex, largey, and largez, where l 0 represents the initial computation domain and l maxlage represents the final computation domain. Therefore, the grid widths in the x, y and z directions of the expanded l' th-level calculation region are respectively
The l-th order computing domain size omegal’Comprises the following steps:
wherein pix1To calculate the leftmost coordinate, π, in the x-direction of the fieldx2To calculate the rightmost coordinate, π, of the field x-direction y1 is the leftmost coordinate of the y direction of the calculation field, piy2To calculate the rightmost coordinate, π, of the y-direction of the fieldz1To calculate the leftmost coordinate, π, in the z-direction of the domainz2To calculate the rightmost coordinate in the z-direction of the field.
The x, y, z coordinates (x) of each calculation grid point in the l' th calculation domaini,yjzk) Comprises the following steps:
wherein i, j, k are the ith grid in the x direction, the jth grid in the y direction and the kth grid in the z direction of the grid in the CPU.
And 8: and assigning functions of density, three-direction speed, pressure and Level-set of all grid points of each CPU in the expanded calculation domain (the l' th Level). The grids within the expanded range are classified into three categories: newly calculating domain points, namely the grid points of the heart obtained after the expansion; the original calculation domain grid points are grid points which exist in the original calculation domain before and after the expansion and can be superposed; new points of the original calculation domain, namely grid points which exist in the original calculation domain before and after the expansion and are not coincident;
for the grid points of the original computing domain, a method of directly assigning values to the l-th computing domain can be adopted.
Aiming at new points of the original calculation domain, because the space discrete format is 2N-1 order precision, interpolation points with 2N-1 order precision also need to be provided at the points of the calculation grid, and the consistent precision of the whole calculation domain is ensured. Therefore, the points are assigned using 2N-1 order WENO interpolation.
And directly assigning the density, the speed and the pressure of the air at normal temperature and normal pressure aiming at the new calculation domain point.
And step 9: and (4) performing multi-CPU joint solution by adopting an MPI parallel computing module of the self-adaptive enlarged computing domain.
By adopting the steps, the serial program of the single CPU can be subjected to the fast and efficient solution of the self-adaptive expanded calculation domain. But in order to be able to compute large scale problems, parallel computations need to be introduced. In MPI (multi-processor interface) parallelism, each CPU occupies a unique memory space, the density, three-direction speed, pressure and Level-set function values of all grid points in all CPUs in a pre-expansion area (an l-Level calculation area) need to be shared in each CPU in an expanded area (an l '-Level calculation area), and a plurality of adjacent points in the l' -Level calculation area and the l-Level calculation area need to be searched.
By using the parallel block search method, the calculation time of the parallel self-adaptive expanded calculation domain module can be greatly reduced, and the parallel efficiency is improved. The parallel block search method comprises the following steps: due to different boundary conditions, for example, the outflow boundary and the fixed wall in step 3 have a sliding boundary, different parallel block search methods are adopted.
For the outflow boundary, dividing the outflow boundary into 4 blocks (two-dimensional) or 8 blocks (three-dimensional) according to a symmetry axis, and repeating the step 8 in each block region to obtain the density, three-direction speed, pressure and Level-set function values of all grid points of the CPU at the l' th Level;
and (3) assigning values to the partial area of the l 'Level aiming at the sliding boundary and the symmetrical boundary of the fixed wall to obtain the expanded density, three-direction speed, pressure and Level-set function values of all grid points of the CPU of the l' Level, so that the calculation amount of 3/4 (three-dimensional 7/8) can be reduced. The partial areas are: a portion where the l' th-order region coincides with the l-order region; at the moment, the calculation domain is expanded, new sizes in the x direction, the y direction and the z direction are generated, the boundary conditions of each level of area are reset, and the steps 3-9 are repeated continuously; ending the cycle until the specified time is reached;
further comprising the step 10: the steps 1 to 9 are utilized to carry out high-precision large-scale numerical simulation based on the self-adaptive expansion calculation domain, so that the shock wave dissipation is reduced, and the large-scale and ultra-large-scale practical engineering problems can be solved with high precision and high efficiency;
the large-scale high-precision numerical simulation technical field comprises the fields of high-speed killer weapons, aerospace and mechanical engineering.
And predicting a large-scale calculation process in the technical field of large-scale high-precision numerical simulation, wherein the large-scale calculation process comprises shock wave bubbling, explosive air explosion, near-ground explosion, underwater explosion, energy-gathered jet flow, interaction of explosion shock waves and buildings, multi-component mixed gas chemical explosion and supernova explosion.
Has the advantages that:
1. the self-adaptive calculation domain expansion method for large-scale parallel calculation disclosed by the invention can meet the requirements of high precision and can solve the large-scale explosion flow field problem efficiently under the condition of certain calculation resources, the solution efficiency is obviously improved, and further the engineering problems in the technical fields of large-scale high-precision numerical simulation, such as shock wave bubbling, explosive air explosion, near-ground explosion, underwater explosion, energy-gathering jet flow, interaction of explosion shock waves and buildings, multi-component mixed gas chemical explosion, supernova explosion and the like, can be accurately predicted.
2. The invention discloses a self-adaptive calculation domain expansion method for large-scale parallel calculation, which is characterized in that a WENO format is adopted to disperse a Euler equation spatial derivative term, and a TVD Runge-Kutta format is adopted to disperse a Euler equation time derivative term.
3. The invention discloses a self-adaptive calculation domain expansion method for large-scale parallel calculation, which adopts a Levelset and RGFM combined Riemann solution method to process the strong interruption problem of the interaction of multiple material interfaces and has the interruption characteristic of high resolution of the material interfaces; compared with the classical GFM method, the method can effectively treat the complex flowing process, effectively inhibit the non-physical oscillation phenomenon near the interface and ensure the stable operation of the numerical simulation process.
4. The invention discloses a self-adaptive calculation domain expansion method for large-scale parallel calculation, which solves the large-scale explosion problem with high precision and high efficiency by utilizing a shock wave/sparse wave expansion threshold condition, an expanded grid high-precision assignment method and a parallel self-adaptive expansion module.
Drawings
FIG. 1 is a diagram of an exemplary numerical model and observation points;
FIG. 2 is a schematic diagram of a three-dimensional model of each stage of the calculation example;
FIG. 3 is a schematic diagram of a shockwave threshold condition;
FIG. 4 is a schematic diagram of sparse wave threshold conditions;
FIG. 5 is a diagram illustrating an adaptive expanding computation domain assignment method;
FIG. 6 is a diagram illustrating a parallel search method under outflow boundary conditions;
FIG. 7 is a schematic diagram of a parallel search method under the condition of a symmetric/fixed wall sliding boundary;
FIG. 8 is a pre-dilation and post-dilation pressure cloud plot for an example t of 0.025 ms;
FIG. 9 is a three-dimensional pressure cloud of the example at times 0.05ms, 0.1ms, 0.2ms, 0.5ms, 1.0ms, 5.0ms, 15ms, and 25 ms;
FIG. 10 is a graph of numerical simulation and experimental comparison of pressure-time curves of examples at distances of 3m, 5m, 7m, 9m, and 11m from the center of the explosion.
Detailed Description
The invention is further described with reference to the following figures and examples.
Example 1:
the present embodiment discloses a method for adaptively expanding a computation domain for massively parallel computing, and for better illustrating the objects and advantages of the present invention, the present invention is further described below with reference to the accompanying drawings and embodiments.
This example is the application in the near-surface detonation problem where 2kgTNT explosive is placed 1.6 meters above the ground.
The initial geometric model in the calculation example disclosed in the embodiment is shown in fig. 1 and fig. 2, the propagation effects of the calculated explosion shock wave, the ground reflected wave and the mach wave formed by superposing the explosion shock wave and the ground reflected wave at different times are shown in fig. 8 and fig. 9, and the comparison result of the typical position pressure-time curve and the experiment is shown in fig. 10.
First, the information required in steps 1,2 is given: as shown in FIG. 1 and FIG. 2, the present embodiment adopts 1/4 model calculation, where x-and y-directions are symmetric boundaries, and x +, y +, and z + are outflow edgesThe boundary, z-is the outflow boundary when the computation domain is not in contact with the ground (i.e. l is 0,1,2 before the third-stage computation domain), and the computation domain is the fixed-wall sliding boundary after the approach ground z is-1.6 (i is 3,4, 5). Largex ═ largey ═ largez ═ 5, i.e., 5-order expansion, and the expansion ratio was taken as r ═ 2. The explosive is placed 1.6 meters from the ground, so the influence of the ground needs to be considered when setting the calculation domain and the boundary conditions thereof. Final calculated Domain size is set to [0,12.0 ]]×[0,12.0]×[-1.6,11.2]The number of parallel partitions in the x, y, and z directions cutx-cuty-6 and cutz-4, and the number of grids in each CPU is Nx×Ny×NzThe total number of the CPUs is 144, namely 50 × 50 × 80, and the total number of the grids is 2880 ten thousand. The calculation stop timing is set to 40 ms. The explosion center is selected as O (0,0,0), and when the calculation domain is not in contact with the ground (i.e., before the third-stage calculation domain, i.e., l ═ 0,1,2), the calculation domain is expanded in the directions of x +, y +, z-, and z +, and after the calculation domain is close to the ground, z ═ 1.6 (i ═ 3,4,5), the calculation domain is expanded only in the directions of x +, y +, and z +.
The CPU serial number of the calculation domain is myid, which is taken as 0-143, and the CPU serial number myid and posx、posy、poszThe relationship of (1) is:
posx＝mod(myid,6)+1
posy＝myid/6+1
posz＝myid/36+1
grid width of initial calculation domain is Deltax0＝Δy0＝Δz0＝0.00125m。
The size of the initial computational domain is Ω0＝[0,0.375]×[0,0.375]×[-0.2,0.2]
The x, y, z coordinates (x) of each grid in the domain are initially computedi,yjzk) Is composed of
xi＝i·Δxl+(posx-1)·50·Δxl
yj＝j·Δyl+(posy-1)·50·Δyl
zk＝k·Δzl+(posz-1)·80·Δzl-0.2
The TNT explosive is positioned at the center of (0,0,0), has the mass of 2kg and the radius of 0.0664m, and the other substances around the TNT explosive are air at normal temperature and normal pressure. The Level set signed distance function (explosive defined as negative and air defined as positive) for each point at the initial time (l ═ 0 computational domain) is
Initial state of explosive: density p1＝1630kg/m3Speed u1＝v1＝w 10, pressure p1＝8600MPa；
Initial state of air: density p2＝1.225kg/m3Speed u2＝v2＝w 20, pressure p2＝0.101MPa；
Establishing a conservation type Euler control equation set as a formula (7), and substituting the density, the speed and the pressure of the two substances into the formula (7) according to the region where the explosive and the detonation product are located.
The explosive adopts JWL equation of state, and the parameters are as follows: a-373800 MPa, B-3747 MPa and R1＝4.15、R2The ideal gas state equation is adopted for air with the parameters of gamma being 1.4, and the parameters of omega being 0.9 and 0.35.
Thereby, steps 1 and 2 are completed.
And then, circularly performing the step 3 to the step 5. Since the present embodiment adopts 1/4 model calculation and considers the existence of the ground, x-and y-directions are set as symmetric boundaries, x +, y +, and z + are set as outflow boundaries, z-is set as an outflow boundary when the calculation domain is not in contact with the ground (i.e., l is 0,1,2 before the third-stage calculation domain), and z is set as a solid wall sliding boundary after the calculation domain is close to the ground z is-1.6 (i is 3,4, 5). L-th stage calculation time step Δ tlIs given by the formula (10). The positions of interfaces in the flow field can be pushed by using the formulas (11) and (12), and the physical state of the full flow field can be updated by using a high-precision WENO format and a real virtual fluid method. Thus, step 3-5 is completed.
And 5, judging in step 6 after step 5, and automatically expanding in steps 7-9 if the pressure at the boundary of the calculation domain reaches the following conditions.
When the shock wave is about to propagate up to the initial calculation domain boundary, when the boundary is right and only one grid point (point x in FIG. 3) awayM-1) If the critical point satisfies the following threshold condition, the calculation domain is expanded,
wherein the overpressure Δ pn(＝pn-0.101) overpressure in the x +, y +, z-and z + directions, px+,py+,pz-And pz+To calculate the pressure on the right, back, bottom, top of the field.
When a rarefaction wave is about to propagate to the right of the boundary and only one grid point (point x in fig. 4) awayM-1) When this point satisfies the following threshold condition, the calculation domain is about to expand,
wherein the determination of the pressure on the ground p is due to the presence of the ground* z-Effective only in level 0,1,2 computational domains, the floor pressure can be written
In step 7, the expanded computational domain grid width is
Δx1＝Δy1＝Δz1＝0.0025m
Δx2＝Δy2＝Δz2＝0.005m
Δx3＝Δy3＝Δz3＝0.01m
Δx4＝Δy4＝Δz4＝0.02m
Δx5＝Δy5＝Δz5＝0.04m
The size of the expanded computational domain is
Ω1＝[0,0.75]×[0,0.75]×[-0.4,0.4]
Ω2＝[0,1.5]×[0,1.5]×[-0.8,0.8]
Ω3＝[0,3.0]×[0,3.0]×[-1.6,1.6](close to the ground, not expanding in the z-direction at this time and thereafter)
Ω4＝[0,6.0]×[0,6.0]×[-1.6,4.8]
Ω5＝[0,12.0]×[0,12.0]×[-1.6,11.2](Final calculation Domain)
X, y, z coordinates (x) of each grid in the 1,2, 3-th computation domaini,yjzk) Comprises the following steps:
xi＝i·Δxl+(posx-1)·50·Δxl
yj＝j·Δyl+(posy-1)·50·Δyl
x, y, z coordinates (x) of each grid in the 4,5 th-level computation domaini,yjzk) Comprises the following steps:
xi＝i·Δxl+(posx-1)·50·Δxl
yj＝j·Δyl+(posy-1)·50·Δyl
zk＝k·Δzl+(posz-1)·80·Δzl-1.6
for the assignment method in step 8, the following three types of points are analyzed, such as the original calculation domain grid point (blue circle), the original calculation domain new point (blue cross) and the new calculation domain point (red circle) in the l + 1-th calculation domain shown in fig. 5.
For the grid points of the original computing domain, a method of directly assigning values to the l-th computing domain can be adopted.
For the new points of the original calculation domain, because the spatial discrete format is 5-order precision, interpolation points with 5-order precision also need to be provided at the blue fork points, and the consistent precision of the whole calculation domain is ensured. Therefore, these points are assigned with a WENO interpolation value of 5 th order.
And directly assigning the density, the speed and the pressure of the air at normal temperature and normal pressure aiming at the new calculation domain point.
For the method in step 9, in this embodiment, when the computation domain is not in contact with the ground (i.e., l is 0,1,2 before the third-level computation domain), the outflow boundary parallel block search method is adopted, as shown in fig. 6, the computation domain may be divided into 8 blocks according to the symmetry axis, and data storage, search distance, and assignment are performed in each block region, which can improve the efficiency of the three-dimensional problem 7/8. After the calculation domain is close to the ground z to-1.6 (l is 3,4,5), a fixed wall sliding boundary parallel block search method is adopted, as shown in fig. 7, a partial region where the l ' th Level region and the l ' th Level region coincide is assigned, and the density, the three-direction speed, the pressure and the Level-set function value of all grid points of the CPU of the l ' th Level are obtained after expansion, so that the calculation amount of 3/4 (three-dimensional 7/8) can be reduced.
Therefore, the whole physical process of the near-surface explosion can be solved jointly by the self-adaptive expansion calculation method and the WENO format, Level set interface tracking and real virtual fluid interface processing method.
And (4) analyzing results:
fig. 8 is a graph of the result at the time t equal to 0.025ms (only x and z directions are shown), the left side is a calculation domain of l equal to 0, and the right side is a calculation domain of l equal to 1 after expansion, at this time, the pressure disturbance just reaches the top end and the bottom end of the calculation domain of l equal to 0, which leads to the first adaptive expansion, and as can be seen from comparison of the left graph and the right graph, the pressure results are consistent. FIG. 9 is a pressure cloud plot from t 0.05ms to t 25ms, lower right corner Ω1～Ω5Representing at which level the computational domain is now. It can be seen from the figure that the computational domain is slowly expanding and that at 0.05ms the pressure peak decays and the shock wave propagates outward, the maximum pressure of the second stage computational domain being significantly lower than the first stage. The dilation threshold condition is reached at 0.074ms, 0.217ms, 0.730, 6.531 ms. The shock wave reaches the ground at about 1.0ms and forms a ground reflection wave, and the ground reflection wave has a high speed and slowly catches up with the front shock wave to form a Mach wave by superposition.
The numerical simulation and experiment places observation points and sensors at positions 3m, 5m, 7m, 9m and 11m away from the center of explosion, the calculated pressure-time curve and the experiment curve with the same dosage are shown in the attached drawing 10, the numerical simulation result can embody the formation of shock waves, reflected waves and Mach waves, the pressure-time curve at each point is basically consistent with the experiment, and the effectiveness of the numerical simulation technology of the self-adaptive expanding calculation domain is proved.
The document (J.Yanez, International Journal of Hydrogen Energy,2011,36: 2613-. While the structured cubic grid is adopted in the embodiment, the grid width is 0.00125m at the minimum, and the calculation is performed by using the 144CPU +41472000 grid, which takes about 138 hours. In the embodiment, a fine grid is adopted for calculation, an 943718400000 grid is needed in a common method, and only a 41472000 grid is needed in the method of the present invention, so that if the problem of the embodiment is solved by adopting the above three types of software, under ideal conditions (without considering factors such as parallel data transfer), about 30782 hours, 102608 hours and 1659104 hours are needed, which are far longer than 138 hours in the embodiment. Therefore, compared with the existing simulation method, the method provided by the invention can be used for rapidly solving the explosion problems and calculating the large-scale explosion problems under the condition of certain calculation resources.
The shock wave, the reflected wave and the Mach wave which are obtained as the result of the embodiment can damage personnel and equipment, the equipment can be deformed until being damaged by the high-pressure shock wave, and the phenomena of eardrum rupture, death and the like can be generated when the personnel are subjected to the high-pressure shock wave and the low-pressure shock wave. It is therefore necessary that the present invention be used in this embodiment.
The above detailed description is intended to illustrate the objects, aspects and advantages of the present invention, and it should be understood that the above detailed description is only exemplary of the present invention and is not intended to limit the scope of the present invention, and any modifications, equivalents, improvements and the like made within the spirit and principle of the present invention should be included in the scope of the present invention.
Claims (1)
1. A numerical simulation method of a self-adaptive expansion calculation domain of large-scale parallel calculation is characterized by comprising the following steps of: the method comprises the following steps:
step 1: determining the size of a final calculation domain, the number of parallel partitions in the x direction, the y direction and the z direction, and the expansion times and the expansion ratio in the self-adaptive expansion calculation domain method and the number of grids in each partition according to the explosion flow field scale and calculation resources needing to be simulated; giving serial numbers of CPUs in a calculation domain; determining the grid widths in the x, y and z directions when the initial calculation domain l is 0, the sizes of the corresponding calculation domains and the coordinates of each calculation grid point; the computing resources comprise the number of CPU cores and the size of a memory;
step 1.1: determining final calculated field size [ pi ]x1,πx2]*[πy1,πy2]*[πz1,πz2]The parallel partition numbers cutx, cuty and cutz in the three directions of x, y and z, and the grid number N in each partitionx×Ny×NzIn which N isxDenotes the number of x-direction grids, NyRepresenting the number of grids in the y-direction, NzRepresents the number of grids in the z direction; the total amount of CPU used in parallel is cutx × cuty × cutz, and the flow field is exploded to calculate the central point O (x)O,yO,zO) Expanding in the directions of x-, x +, y-, y +, z-and z +, expanding times in the directions of x, y and z, namely largex, largey and largez, and expanding the ratio r to obtain the total grid number of cutx cut Nx*Ny*Nz；
Step 1.2: giving serial numbers of CPUs in a calculation domain;
defining myid as a CPU serial number from 0; posx、posy、poszPos indicating that the CPU is in the x directionxPos in the direction of (y)yPos in the z directionzA CPU; posxIs in the range of 1 to cutx, posyIs in the range of 1 to cut, poszThe value range of (1) to cutz;
CPU serial numbers myid and posx、posy、poszThe relationship of (1) is:
wherein mod is a remainder symbol;
step 1.3: when the initial calculation field l is 0, the grid widths of the calculation region in the x, y and z directions are respectively
Calculating the domain size omega when the initial calculation domain l is 00Comprises the following steps:
wherein r islargexIs the largex power of r, rlargeyIs the largey power of r, rlargezIs the power of r to the largez power, where pix1For the final calculation of the leftmost coordinate, π, in the x-direction of the domainx2For the final calculation of the rightmost coordinate, π, in the x-direction of the domainy1For final calculation of the leftmost coordinate, π, in the y-direction of the domainy2For final calculation of the rightmost coordinate, π, in the y-direction of the fieldz1For the final calculation of the leftmost coordinate, π, in the z-direction of the domainz2Calculating the coordinate of the rightmost side of the domain in the z direction finally;
thus, the x, y, z coordinates (x) of each grid in the domain are initially computedi,yjzk) Comprises the following steps:
wherein i, j, k is the ith grid in the x direction, the jth grid in the y direction and the kth grid in the z direction of the grid in the CPU;
step 2: defining a level-set function for distinguishing detonation products and air interface positions at the moment when the initial state t is 0 in the initial calculation region; assigning initial values to physical quantities of two regions of detonation products and air; establishing a conservation type Euler control equation set for describing detonation products and air; establishing a state equation function of a detonation product and air;
step 2.1: defining a level-set function for distinguishing detonation products and air interface positions at the moment when the initial state t is 0 in an initial calculation region, giving the position of a multi-material interface at any moment by using a symbol distance function, tracking the evolution process of the level-set function in real time, and locating the interface gamma (t) at the moment t in the level-set function phi (x)i,yj,zkT) position of zero:
Γ(t)＝{(xi,yj,zk)|φ(xi,yj,zk,t)＝0} (5)
and at the moment when t is 0, calculating the distance between the position of each point in the calculation domain and the interface position according to the interface position and the detonation product and the air interface position:
wherein x, y, z are physical space coordinates of the calculation grid points, d (Γ (0), (x), respectivelyi,yj,zk) Denotes a calculation grid point (x)i,yj,zk) Distance to initial interface Γ (0), SDetonation productsAnd SAir (a)The areas on both sides of the material interface are respectively;
step 2.2: giving the initial states of the detonation product and the air obtained in the step 2.1;
initial density ρ of a given detonation product1Three directional velocities u1、v1、w1Pressure p1(ii) a Or the initial density ρ of detonation products1Three directional velocities u1、v1、w1Internal energy per unit mass e1；
And initial density ρ of air2Three directional velocities u2、v2、w2Pressure p2(ii) a Or initial density ρ of air2Three directional velocities u2、v2、w2Internal energy per unit mass e2；
Step 2.3: constructing a conservation type Euler control equation set for describing detonation products and air;
the conservation-oriented Euler control equation system is as follows:
wherein the content of the first and second substances,
U＝[ρ ρu ρv ρw ρE]T，
F＝[ρu ρu2+p ρuv ρuw u(ρE+p)]T，
G＝[ρv ρuv ρv2+p ρvw v(ρE+p)]T，
H＝[ρw ρuw ρvw ρw2+p w(ρE+p)]T，
and E is the total energy per unit mass,
if the detonation product is calculated, the rho, u, v, w, p and e in the formula (7) are used as rho1、u1、v1、w1、p1、e1Substituting; if calculating air, p, u, v, w, p, e in the formula (7) is used as p2、u2、v2、w2、p2、e2Substituting;
step 2.4: selecting a state equation function of the material, and further determining the physical state of the material under deformation and motion;
therefore, a state equation needs to be selected to seal the whole control equation set; internal energy per unit mass e of detonation product1And pressure p1JWL equation of state relationship (sqj) expressed as:
where ρ is0Is the initial density of the explosive, A, B, R1、R2ω is the equation of state parameter; the ideal gas equation of state relationship for air is expressed as:
p2＝(γ-1)ρ2e2(9)
wherein γ is the polytropic gas index;
and step 3: setting boundary conditions of each level of area;
because the sizes of the calculation domains at different moments are different, the boundary conditions of the calculation domains at each level are different; aiming at different spatial discrete formats, each calculation grid point needs the physical quantity calculation flux of N surrounding points, so N layers of virtual points are needed at the boundary of a calculation domain; rho, u, v, w, p and E appearing below are taken as rho if the components are in a detonation product region1、u1、v1、w1、p1、E1In the air region, ρ is2、u2、v2、w2、p2、E2。
(1) For an outgoing boundary, at each level of computational domain boundary, there is
① x direction and when posxWhen 1, there are
i＝1,2,...,N,j＝0,1,...,Ny,k＝0,1,...,Nz；
② x direction and when posxWhen n x, there is
i＝1,2,...,N,j＝0,1,...,Ny,k＝0,1,...,Nz；
③ y direction and when posyWhen 1, there are
i＝0,1,...,Nx,j＝1,2,...,N,k＝0,1,...,Nz；
④ y direction and when posyNy, there are
i＝0,1,...,Nx,j＝1,2,...,N,k＝0,1,...,Nz；
⑤ z direction and when poszWhen 1, there are
i＝0,1,...,Nx,j＝0,1,...,Ny,k＝1,2,...,N；
⑥ z direction and when poszWhen n is nz, there is
i＝0,1,...,Nx,j＝0,1,...,Ny,k＝1,2,...,N
(2) There are a slip boundary and a symmetric boundary for the fixed wall,
① x direction and when posxWhen 1, there are
φ(x-i,yj,zk,t)＝φ(xi,yj,zk,t)，
i＝1,2,...,N,j＝0,1,...,Ny,k＝0,1,...,Nz；
② x direction and when posxWhen n x, there is
φ(xnx+i,yj,zk,t)＝φ(xnx-i,yj,zk,t)，
i＝1,2,...,N,j＝0,1,...,Ny,k＝0,1,...,Nz；
③ y direction and when posyWhen 1, there are
φ(xi,y-j,zk,t)＝φ(xi,yj,zk,t)，
i＝0,1,...,Nx,j＝1,2,...,N,k＝0,1,...,Nz；
④ y direction and when posyNy, there are
φ(xi,yny+j,zk,t)＝φ(xi,yny-j,zk,t)，
i＝0,1,...,Nx,j＝1,2,...,N,k＝0,1,...,Nz；
⑤ z direction and when poszWhen 1, there are
φ(xi,yj,z-k,t)＝φ(xi,yj,zk,t)，
i＝0,1,...,Nx,j＝0,1,...,Ny,k＝1,2,...,N；
⑥ z direction and when poszWhen n is nz, there is
φ(xi,yj,znz+k,t)＝φ(xi,yj,znz-k,t)，
i＝0,1,...,Nx,j＝0,1,...,Ny,k＝1,2,...,N
The outflow boundary comprises, for example, an infinite volume of air;
the fixed wall with a sliding boundary comprises a rigid ground;
step 3 the symmetric boundary comprises an 1/2 model, a 1/4 model, or a 1/8 model;
and 4, step 4: selecting the l-level calculation time step;
selecting the l-th level to calculate the time step delta tl：
Wherein the parameter CFL is a constant between 0 and 1, Δ xlCalculating the grid width in the x direction in the area for the l level; Δ ylCalculating the grid width in the y direction in the area for the l level; Δ zlCalculating the grid width in the z direction in the area for the l level; u, v and w are the speeds of the three directions of the current calculation grid point, and if the current calculation grid point is in the detonation product region, u is taken as1、v1、w1If in the air region, it is taken as u2、v2、w2(ii) a c is the sound velocity of the current calculated grid point;
and 5: solving a normal vector of a material interface in each Level of region by using a Level-set function, solving a normal Riemann problem at the interface of a detonation product and air, determining the physical state of N layers of grids near the interface by using a virtual fluid method, and updating the physical state of the full flow field;
using Level-set equation to advance the interface position in the flow field, using the following formula to advance for each calculation grid point in the calculation domain, and using phi (x) to simplify the formulai,yj,zkT) is abbreviated as φ;
if the calculated grid point is in the detonation product region, u, v and w in the formula (11) are used as u1、v1、w1Substituting; if the calculation grid point is in the air region, u, v, w in the formula (11) are used as u2、v2、w2Substituting;
the solving formula of the normal vector n of each calculation grid point is as follows:
decomposing the speed of each calculation grid point according to a normal vector of the point to obtain a normal speed and a tangential speed, then giving the normal speed and the pressure at the interface and the densities of the detonation products and the air at two sides of the interface by a Riemann solver according to the physical quantity states of the detonation products, the density of the air, the normal speed and the pressure near the interface, the continuity of the normal speed and the pressure at two sides of the interface and the conservation law, wherein the normal speed and the tangential speed are physical states under a local coordinate system, the normal speed of the calculation grid points near the interface needs to be kept unchanged, and the speed component is reversely transformed to obtain the speed under a global coordinate system; respectively establishing a layer of banded regions with 20 grid widths in detonation products and air regions at two sides of the interface, wherein calculation grid points in the banded regions are virtual points; the density, the three-direction speed and the pressure at the virtual point are obtained by continuation; the continuation equation is:
wherein I represents a physical quantity at a virtual point, the physical quantity at the virtual point including density, three directional velocity, and pressure; finally, solving the detonation product and the air respectively by using a high-precision WENO format;
step 6: judging whether the current calculation domain needs to be expanded or not;
after a period of time step is calculated, detonation products in the initial calculation domain expand outwards, namely reach the boundary of the initial calculation domain, and shock waves (physical quantity mutation) and rarefaction waves (physical quantity slow rising or falling) are respectively adopted for analysis; judging whether expansion is needed; if the expansion is judged not to be needed, continuing to operate the step 3, and if the expansion is judged to be about to be carried out, executing the steps 7-9;
when the shock wave is about to propagate to reach the boundary of the initial calculation domain, the shock wave is located at the position xMIs a critical point; if the critical point satisfies the following threshold condition, the calculation domain is expanded,
|IM-I0|＞Icr(14)
wherein, I can represent density, pressure; i is0Is the constant density, constant pressure, I, of the airMDensity, pressure at critical point; i iscrIs a critical threshold value, and is generally a small number, so that when the shock wave reaches a critical point, the calculation domain is automatically expanded to prevent the shock wave from flowing out of the boundary guide of the calculation domainLoss of important information;
when the sparse wave is about to propagate to reach the boundary of the initial calculation domain, the position x is locatedNIs a critical point; when this point satisfies the following threshold condition, the calculation domain is about to expand,
wherein I may represent density, pressure; i isNDensity, pressure at critical point; i'crA critical gradient threshold, typically taken as a small number and related to the grid width; if the shockwave threshold (as in equation 14) is still used, within a finite computation step, at the boundary point xNA "step" is created, thereby introducing an error; therefore, whether the calculation domain is enlarged is judged by adopting the formula (15), so that important information loss caused by sparse waves flowing out of the boundary of the calculation domain can be prevented;
and 7: in the calculation process, detonation products are spread outwards, pressure and density are spread outwards, the boundary of the first-level calculation domain meets the expansion threshold condition in the step 6, at the moment, l' ═ l +1 is the calculation domain after expansion, and the calculation domain of the calculation domain after expansion, the corresponding grid width and the coordinate of each calculation grid point in each calculation domain need to be given;
ranking the computation domain as l ═ 0, 1., maxlage, wherein maxlage is the maximum of largex, largey, and largez, wherein l ═ 0 represents the initial computation domain, and l ═ maxlage represents the final computation domain; therefore, the grid widths in the x, y and z directions of the expanded l' th-level calculation region are respectively
The l-th order computing domain size omegal'Comprises the following steps:
wherein pix1To calculate the leftmost coordinate, π, in the x-direction of the fieldx2To calculate the rightmost coordinate, π, of the field x-directiony1To calculate the leftmost coordinate, π, in the y-direction of the fieldy2To calculate the rightmost coordinate, π, of the y-direction of the fieldz1To calculate the leftmost coordinate, π, in the z-direction of the domainz2Calculating the coordinate of the rightmost side of the domain in the z direction;
the x, y, z coordinates (x) of each calculation grid point in the l' th calculation domaini,yjzk) Comprises the following steps:
wherein i, j, k is the ith grid in the x direction, the jth grid in the y direction and the kth grid in the z direction of the grid in the CPU;
and 8: assigning density, three-direction speed, pressure and Level-set functions of all grid points of each CPU in the expanded calculation domain (the l' th Level); the grids within the expanded range are classified into three categories: newly calculating domain points, namely the grid points of the heart obtained after the expansion; the original calculation domain grid points are grid points which exist in the original calculation domain before and after the expansion and can be superposed; new points of the original calculation domain, namely grid points which exist in the original calculation domain before and after the expansion and are not coincident;
aiming at the grid points of the original computing domain, a method of directly assigning values to the l-level computing domain can be adopted;
aiming at new points of an original calculation domain, because the space discrete format is 2N-1 order precision, interpolation points with 2N-1 order precision also need to be provided at the points of the calculation grid, and the consistent precision of the whole calculation domain is ensured; therefore, 2N-1 WENO interpolation is adopted to assign values to the points;
aiming at the new calculation domain point, the density, the speed and the pressure of the air at normal temperature and normal pressure can be directly assigned;
and step 9: adopting an MPI parallel computing module of a self-adaptive enlarged computing domain to carry out multi-CPU combined solution;
by adopting the steps, the serial program of the single CPU can be subjected to fast and efficient solution of the self-adaptive expanded calculation domain; however, in order to be able to compute large scale problems, parallel computation needs to be introduced; in MPI (message passing interface) parallelism, each CPU (Central processing Unit) occupies a unique memory space, the density, three-direction speed, pressure and Level-set function values of all grid points in all CPUs in a pre-expansion area (an l-Level calculation area) need to be shared into each CPU in an expanded area (an l '-Level calculation area), and a plurality of adjacent points in the l' -Level calculation area and the l-Level calculation area need to be searched;
by using the parallel block search method, the calculation time of the parallel self-adaptive expanded calculation domain module can be greatly reduced, and the parallel efficiency is improved; the parallel block search method comprises the following steps: due to different boundary conditions, if the outflow boundary and the fixed wall in the step 3 have a sliding boundary, different parallel block searching methods are adopted;
for the outflow boundary, dividing the outflow boundary into 4 blocks (two-dimensional) or 8 blocks (three-dimensional) according to a symmetry axis, and repeating the step 8 in each block region to obtain the density, three-direction speed, pressure and Level-set function values of all grid points of the CPU at the l' th Level;
aiming at the fact that a fixed wall has a sliding boundary and a symmetrical boundary, the I 'Level partial area is assigned, and the density, the three-direction speed, the pressure and the Level-set function value of all grid points of each CPU of the I' Level after expansion are obtained, so that the calculation amount of 3/4 (three-dimensional 7/8) can be reduced; the partial areas are: a portion where the l' th-order region coincides with the l-order region; at the moment, the calculation domain is expanded, new sizes in the x direction, the y direction and the z direction are generated, the boundary conditions of each level of area are reset, and the steps 3-9 are repeated continuously; ending the cycle until the specified time is reached; the method can reduce shock wave dissipation and can solve the actual engineering problems of large scale and ultra-large scale with high precision and high efficiency.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201911001549.1A CN110852005B (en) | 2019-10-21 | 2019-10-21 | Numerical simulation method for self-adaptive expansion of computational domain of large-scale parallel computation |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201911001549.1A CN110852005B (en) | 2019-10-21 | 2019-10-21 | Numerical simulation method for self-adaptive expansion of computational domain of large-scale parallel computation |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110852005A true CN110852005A (en) | 2020-02-28 |
| CN110852005B CN110852005B (en) | 2021-06-15 |
Family
ID=69597824
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201911001549.1A Active CN110852005B (en) | 2019-10-21 | 2019-10-21 | Numerical simulation method for self-adaptive expansion of computational domain of large-scale parallel computation |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110852005B (en) |
Cited By (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111783365A (en) * | 2020-06-04 | 2020-10-16 | 三多(杭州)科技有限公司 | Virtual medium method, device and equipment applied to low-voltage interface processing |
| CN111859530A (en) * | 2020-06-11 | 2020-10-30 | 北京航空航天大学 | Iterative propulsion disturbance domain updating method for aircraft dynamic aerodynamic characteristic simulation |
| CN111859819A (en) * | 2020-06-16 | 2020-10-30 | 空气动力学国家重点实验室 | Construction method of high-order WENO format |
| CN111859529A (en) * | 2020-06-11 | 2020-10-30 | 北京航空航天大学 | Multi-grid disturbance domain updating acceleration method for aircraft streaming numerical simulation |
| CN113408168A (en) * | 2021-06-18 | 2021-09-17 | 北京理工大学 | High-precision numerical simulation method based on Riemann problem accurate solution |
| CN113836659A (en) * | 2021-09-22 | 2021-12-24 | 西安石大派普特科技工程有限公司 | Method for characterizing performance characteristic color gradation chart of centrifugal compressor |
| CN115329646A (en) * | 2022-10-12 | 2022-11-11 | 国家超级计算天津中心 | Shock wave simulation method, device, equipment and medium |
| CN116050196A (en) * | 2023-04-03 | 2023-05-02 | 国家超级计算天津中心 | Multi-dimensional simulation method, device, equipment and storage medium |
| CN118070574A (en) * | 2024-04-24 | 2024-05-24 | 国家超级计算天津中心 | Parallel simulation method, device and storage medium |
| WO2024119956A1 (en) * | 2022-12-07 | 2024-06-13 | 腾讯科技（深圳）有限公司 | Control method and apparatus for virtual projectile, device, medium, and program product |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104375975A (en) * | 2014-12-01 | 2015-02-25 | 天津工业大学 | One-dimensional vacuum Crank-Nicolson complete matching layer implementation algorithm based on bilinear transformation |
| CN104461466A (en) * | 2013-09-25 | 2015-03-25 | 广州中国科学院软件应用技术研究所 | Method for increasing computing speed through parallel computing based on MPI and OpenMP hybrid programming model |
| US20150127311A1 (en) * | 2013-11-06 | 2015-05-07 | Weidlinger Associates, Inc. | Computer Implemented Apparatus and Method for Finite Element Modeling Using Hybrid Absorbing Element |
| CN105260249A (en) * | 2015-09-19 | 2016-01-20 | 中国地质大学（武汉） | Method for extracting calculation intensity features of spatial calculation domain |
| CN109902376A (en) * | 2019-02-25 | 2019-06-18 | 北京理工大学 | A kind of fluid structurecoupling high resolution numerical simulation method based on Continuum Mechanics |
-
2019
- 2019-10-21 CN CN201911001549.1A patent/CN110852005B/en active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104461466A (en) * | 2013-09-25 | 2015-03-25 | 广州中国科学院软件应用技术研究所 | Method for increasing computing speed through parallel computing based on MPI and OpenMP hybrid programming model |
| US20150127311A1 (en) * | 2013-11-06 | 2015-05-07 | Weidlinger Associates, Inc. | Computer Implemented Apparatus and Method for Finite Element Modeling Using Hybrid Absorbing Element |
| CN104375975A (en) * | 2014-12-01 | 2015-02-25 | 天津工业大学 | One-dimensional vacuum Crank-Nicolson complete matching layer implementation algorithm based on bilinear transformation |
| CN105260249A (en) * | 2015-09-19 | 2016-01-20 | 中国地质大学（武汉） | Method for extracting calculation intensity features of spatial calculation domain |
| CN109902376A (en) * | 2019-02-25 | 2019-06-18 | 北京理工大学 | A kind of fluid structurecoupling high resolution numerical simulation method based on Continuum Mechanics |
Non-Patent Citations (1)
| Title |
|---|
| 马天宝,等: "爆炸与冲击问题的大规模高精度计算", 《力学学报》 * |
Cited By (15)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111783365A (en) * | 2020-06-04 | 2020-10-16 | 三多(杭州)科技有限公司 | Virtual medium method, device and equipment applied to low-voltage interface processing |
| CN111783365B (en) * | 2020-06-04 | 2023-09-22 | 三多(杭州)科技有限公司 | Virtual medium method, device and equipment applied to low-voltage interface processing |
| CN111859529B (en) * | 2020-06-11 | 2022-04-08 | 北京航空航天大学 | Multi-grid disturbance domain updating acceleration method for aircraft streaming numerical simulation |
| CN111859529A (en) * | 2020-06-11 | 2020-10-30 | 北京航空航天大学 | Multi-grid disturbance domain updating acceleration method for aircraft streaming numerical simulation |
| CN111859530B (en) * | 2020-06-11 | 2022-04-22 | 北京航空航天大学 | Iterative propulsion disturbance domain updating method for aircraft dynamic aerodynamic characteristic simulation |
| CN111859530A (en) * | 2020-06-11 | 2020-10-30 | 北京航空航天大学 | Iterative propulsion disturbance domain updating method for aircraft dynamic aerodynamic characteristic simulation |
| CN111859819A (en) * | 2020-06-16 | 2020-10-30 | 空气动力学国家重点实验室 | Construction method of high-order WENO format |
| CN113408168A (en) * | 2021-06-18 | 2021-09-17 | 北京理工大学 | High-precision numerical simulation method based on Riemann problem accurate solution |
| CN113408168B (en) * | 2021-06-18 | 2022-11-08 | 北京理工大学 | High-precision numerical simulation method based on Riemann problem accurate solution |
| CN113836659A (en) * | 2021-09-22 | 2021-12-24 | 西安石大派普特科技工程有限公司 | Method for characterizing performance characteristic color gradation chart of centrifugal compressor |
| CN113836659B (en) * | 2021-09-22 | 2024-01-30 | 西安石大派普特科技工程有限公司 | Method for characterizing performance characteristic color level chart of centrifugal compressor |
| CN115329646A (en) * | 2022-10-12 | 2022-11-11 | 国家超级计算天津中心 | Shock wave simulation method, device, equipment and medium |
| WO2024119956A1 (en) * | 2022-12-07 | 2024-06-13 | 腾讯科技（深圳）有限公司 | Control method and apparatus for virtual projectile, device, medium, and program product |
| CN116050196A (en) * | 2023-04-03 | 2023-05-02 | 国家超级计算天津中心 | Multi-dimensional simulation method, device, equipment and storage medium |
| CN118070574A (en) * | 2024-04-24 | 2024-05-24 | 国家超级计算天津中心 | Parallel simulation method, device and storage medium |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110852005B (en) | 2021-06-15 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110852005B (en) | Numerical simulation method for self-adaptive expansion of computational domain of large-scale parallel computation | |
| US10296672B2 (en) | Generating inviscid and viscous fluid-flow simulations over an aircraft surface using a fluid-flow mesh | |
| Smiljanovski et al. | A capturing-tracking hybrid scheme for deflagration discontinuities | |
| Rizzi et al. | Computation of flow around wings based on the Euler equations | |
| Daxini et al. | A review on recent contribution of meshfree methods to structure and fracture mechanics applications | |
| Hannemann | High enthalpy flows in the HEG shock tunnel: experiment and numerical rebuilding | |
| Arai et al. | Ultra-large scale fracture mechanics analysis using a parallel finite element method with submodel technique | |
| Mohler | Wind-US flow calculations for the M2129 S-duct using structured and unstructured grids | |
| Sansica et al. | Global stability analysis of full-aircraft transonic buffet at flight Reynolds numbers | |
| Saito et al. | Development and application of high-resolution adaptive numerical techniques in Shock Wave Research Center | |
| Yin et al. | Lagrangian analysis of unsteady partial cavitating flow around a three-dimensional hydrofoil | |
| Zhan et al. | Three-dimensional high-order finite-volume method based on compact WENO reconstruction with hybrid unstructured grids | |
| Gillberg et al. | Accuracy and efficiency of stencils for the eikonal equation in earth modelling | |
| Deng et al. | A new fully coupled method for computing turbulent flows | |
| Landsberg et al. | An efficient parallel method for solving flows in complex three-dimensional geometries | |
| Margha et al. | Dynamic transition from regular to mach reflection over a moving wedge | |
| Rahmani et al. | Large eddy simulation of the Sandia axisymmetric transonic hump using a high-order method | |
| Manzke | Development of a scalable method for the efficient simulation of flows using dynamic goal-oriented local grid-adaptation | |
| Li et al. | Numerical computations of resonant sloshing using the modified isoAdvector method and the buoyancy-modified turbulence closure model | |
| Bernard et al. | Dispersion analysis of discontinuous Galerkin schemes applied to Poincaré, Kelvin and Rossby waves | |
| Le | Numerical simulation of shock (blast) wave interaction with bodies | |
| Fang et al. | Predictions of flow separation at the valve seat for steady-state port-flow simulation | |
| Tseng et al. | Numerical simulation of vorticity production in shock diffraction | |
| Di Giacinto et al. | Shock detection and discontinuity tracking for unsteady flows | |
| Shah et al. | ANSYS Fluent Scale-Resolving Simulations with SBES & Validation of a Re-Entry Capsule at Hypersonic Speed |
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 |