CN114169184A - High-precision self-adaptive finite volume-finite difference coupling numerical simulation method - Google Patents

High-precision self-adaptive finite volume-finite difference coupling numerical simulation method Download PDF

Info

Publication number
CN114169184A
CN114169184A CN202111217254.5A CN202111217254A CN114169184A CN 114169184 A CN114169184 A CN 114169184A CN 202111217254 A CN202111217254 A CN 202111217254A CN 114169184 A CN114169184 A CN 114169184A
Authority
CN
China
Prior art keywords
calculation
substance
interface
equation
numerical simulation
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
Application number
CN202111217254.5A
Other languages
Chinese (zh)
Other versions
CN114169184B (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.)
Beijing Institute of Technology BIT
Original Assignee
Beijing Institute of Technology BIT
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 Beijing Institute of Technology BIT filed Critical Beijing Institute of Technology BIT
Priority to CN202111217254.5A priority Critical patent/CN114169184B/en
Priority claimed from CN202111217254.5A external-priority patent/CN114169184B/en
Publication of CN114169184A publication Critical patent/CN114169184A/en
Application granted granted Critical
Publication of CN114169184B publication Critical patent/CN114169184B/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
    • G06F30/28Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00Details relating to CAD techniques
    • G06F2111/10Numerical modelling
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00Details relating to the application field
    • G06F2113/08Fluids
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14Force analysis or force optimisation, e.g. static or dynamic forces

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • General Engineering & Computer Science (AREA)
  • Geometry (AREA)
  • Evolutionary Computation (AREA)
  • Fluid Mechanics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computing Systems (AREA)
  • Algebra (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

The invention discloses a high-precision self-adaptive finite volume-finite difference coupling numerical simulation method, and belongs to the field of high-precision numerical simulation of multi-substance coupling. The method can realize the self-adaptive coupling of two numerical methods of finite volume and finite difference, select different numerical formats to calculate according to the self-adaptive algorithm judgment rule at different positions of a calculation domain, and couple different high-precision algorithm calculation regions by adopting the self-adaptive algorithm coupling method, thereby not only obtaining the high-resolution shock wave discontinuity capturing effect of the high-order finite difference method, but also taking the inherent good conservation characteristic of the finite volume method into account, remarkably improving the accuracy of the numerical simulation prediction result and further effectively predicting the engineering technical problems related to the high-precision numerical simulation field. The high-precision numerical simulation field comprises the fields of high-speed/ultrahigh-speed warhead penetration and protection, underwater explosion, bubble dynamics, aerospace and mechanical engineering.

Description

High-precision self-adaptive finite volume-finite difference coupling numerical simulation method
Technical Field
The invention relates to a self-adaptive finite volume-finite difference coupling high-precision numerical simulation method, belonging to the field of multi-substance coupling high-precision numerical simulation.
Background
The bubble dynamics has important application in various industries, such as the fields of ship damage, offshore ice breaking, marine exploration, biomedical treatment and the like, and therefore, the bubble dynamics research has important significance for the research of the bubble dynamics.
The research on underwater explosion bubble dynamics mainly comprises theoretical derivation, numerical simulation and experimental research. The explosion has the characteristics of high temperature, high pressure, high speed and the like, is a typical nonlinear problem, and the bubble jet flow relates to the processes of collision, tearing, fusion and the like of bubbles, has more complex motion characteristics and causes great difficulty in theoretical derivation. The experimental research is very intuitive, but the experimental method has the defects of long research period, high expenditure and the like, the explosion test has high danger, the difficulty in developing a large-equivalent prototype test is high, the measurement and observation of the explosion test also have high difficulty, only the surface phenomenon can be analyzed, and the research on the internal mechanism cannot be carried out. Compared with the former two methods, the numerical simulation method has the advantages of short research period, small expenditure burden, high safety and the like, and internal mechanism analysis and research can be carried out on the representation, so that the research is more deep, but the numerical simulation has higher requirements on numerical format, calculation scale, calculation precision and the like, otherwise, the numerical simulation is only a high-class and cannot be applied.
The existing numerical simulation method mainly comprises the following steps: finite element method, finite difference method and finite volume method. The finite element method has higher requirements on the grid quality, and the grid can move along with the material, so that the grid distortion problem can be caused by the large deformation problem such as explosion impact, the calculation error is increased, and even the calculation termination can be caused; the finite difference method has the main idea that the difference replaces the derivative, and the high-order format only needs to use more node structures and is easy to realize, so the finite difference method has better effect on shock wave discontinuity capture, but the finite difference method has poor conservation due to the adoption of a differential form equation set; the finite volume method is a new method combining finite elements and finite difference methods, adopts an integral form equation, has inherent good conservation characteristic, and is suitable for non-structural grids, but due to the limitation of the method, the realization of a high-order form is difficult, only the second-order precision exists, and the capture effect on shock wave discontinuity is poor. Different numerical simulation methods have defects, so that numerical simulation errors of different formats are different, the reliability is reduced, and the defects need to be compensated urgently, so that the method provided by the invention has important engineering value and application significance.
Disclosure of Invention
Aiming at the problems that numerical simulation errors are increased and the precision of a calculation result is influenced by adopting a single method in the existing numerical simulation method, the characteristics that the conservation of a high-precision finite difference method is poor but the shock wave discontinuity capturing effect is good, a finite volume method has inherent good conservation characteristic but a high-order method is complex to realize and the shock wave discontinuity capturing effect is poor are considered, and the technical problems to be solved by the self-adaptive finite volume-finite difference coupling high-precision numerical simulation method disclosed by the invention are that: the method has the advantages that the self-adaptive coupling of two numerical methods of finite volume-finite difference is realized, different numerical formats are selected for calculation at different positions of a calculation domain according to self-adaptive algorithm judgment rules, and the self-adaptive algorithm coupling method is adopted for coupling calculation regions of different high-precision algorithms, so that the high-resolution shock wave discontinuity capturing effect of the high-order finite difference method can be obtained, the inherent good conservation characteristic of the finite volume method can be considered, the accuracy of the numerical simulation prediction result can be remarkably improved, and further, the engineering technical problems related to the high-precision numerical simulation field can be effectively predicted.
The high-precision numerical simulation field comprises the fields of high-speed/ultrahigh-speed warhead penetration and protection, underwater explosion, bubble dynamics, aerospace and mechanical engineering.
The purpose of the invention is realized by the following technical scheme.
The invention discloses a high-precision self-Adaptive Finite Volume-Finite Difference coupling numerical simulation Method, which is a self-Adaptive coupling Method of a Finite Difference Method and a Finite Volume Method in a computational domain, combines a Weighted essential Non-oscillation (WENO) Finite Difference format with a Conservation-law monotonic center windward (MUSCL) Finite Volume format, and simultaneously couples a Level-Set interface tracking Method and a multi-material Riemann problem solving Method, namely a self-Adaptive high-precision Finite Volume-Finite Difference computing Method (VDM). The high-order finite difference method has a high-resolution shock wave discontinuity capturing effect, and the WENO finite difference format is adopted, so that numerical value oscillation can be reduced while the precision is ensured; the finite volume method has inherent good conservation characteristic, adopts MUSCL finite volume format, has high-resolution contact interruption capturing effect, and simultaneously ensures the conservation of the physical state of a multi-substance interface; the multi-substance Riemann solution decouples the multi-substance interaction problem into a single-substance problem of each substance, and is convenient to solve by adopting a uniform format; the Level-Set method can effectively track the complex multi-material interface evolution and topology change. The method can effectively realize high-precision numerical simulation on the problem of multi-substance interaction, remarkably improve the accuracy of the numerical simulation prediction result, and further effectively predict the engineering technical problems related to the high-precision numerical simulation field.
The invention discloses a high-precision self-adaptive finite volume-finite difference coupling numerical simulation method, which comprises the following steps of:
step 1: determining a calculation region according to actual physical problem requirements, establishing a Cartesian coordinate system, and dividing the calculation region in x, y and z directions into m, n and l grids respectively, wherein the total number of the m, n and l grids is m.
Step 2: defining a level-set function in the calculation area for distinguishing the initial state interface position of each substance, defining the initial physical quantity, the material state equation and the conservation type Euler control equation set of each substance according to the area divided by the level-set, and simultaneously setting the boundary condition according to the actual physical problem.
Step 2.1: and determining the position of the multi-material interface in the initial state according to the actual physical problem, and establishing a level-set function.
Step 2.2: and (3) dividing each substance area according to the step 2.1, and defining initial state physical quantities of substances in different areas.
The initial state physical quantity includes density ρ of each substance in the initial state, and velocities u, v, w in three directions of pressure p, x, y, z.
Step 2.3: and constructing a conservation-type Euler control equation set.
Step 2.4: and selecting a state equation function of the material aiming at the actual physical problem, and further determining the physical state of the material under deformation and motion.
The conservation type Euler control equation set comprises 5 equations and 6 unknowns, and in order to seal the equation set, a proper material state equation needs to be selected to seal the equation set.
Preferably, the JWL equation of state is chosen for the explosive without taking into account the chemical reaction, expressed as:
Figure BDA0003311229960000031
where ρ is0A, B, R1, R2, omega are JWL equation of state parameters for the initial density of the explosive. An ideal gas state equation is selected for air and expressed as
p=(γ-1)ρe (2)
Where γ is the ideal gas adiabatic index.
Tammann equation of state, expressed as
p=(γ-1)ρe-γ(B-A) (3)
Wherein, gamma, A, B are the parameters of Tammann equation of state.
In a low-pressure area formed by sparse waves, because the pressure is less than saturated vapor pressure, cavitation phenomenon can occur in water, in order to carry out more accurate numerical simulation on the cavitation phenomenon, a cavitation model is introduced, and preferably, an improved Schmidt model in a single-fluid cavitation model is selected and expressed as
Figure BDA0003311229960000041
Wherein p issatIs the saturated vapor pressure, pminAt a very small positive pressure, pgAnd ρlDensity of gas and water, respectively, cgAnd clAcoustic velocity, p, of gas and water, respectivelyglIs obtained by the following formula:
Figure BDA0003311229960000042
step 2.5: and setting boundary conditions for the calculation regions divided in the step 2.1.
According to the actual physical problem, a plurality of layers of virtual points are required to be arranged outside a calculation domain, so that different types of boundary conditions are realized.
For non-reflection boundary conditions, such as deep water area and air area, the influence of gravity is considered, so that the pressure gradient boundary condition, p, needs to be set-i=p0+ρg(h-i-h0) N, i is 1,2, i is that the virtual point pressure of the first layer, the second layer and the third layer outside the boundary of the calculation region is equal to the pressure p of the first layer in the boundary of the calculation region0Plus the pressure ρ g (h) due to gravity-i-h0) Boundary setting of other physical quantities adopts U-i=U0And i is 1,2 … N, namely virtual points of the first layer, the second layer and the third layer outside the boundary of the calculation region are equal to physical quantities of the first layer inside the boundary of the calculation region.
For reflecting boundary conditions, e.g. walls, rigid floors, using un,-i=-un,iI is 1,2 … N, namely the normal speeds of the virtual points of the first layer, the second layer and the third layer outside the boundary of the calculation region are respectively equal to the reverse directions of the normal speeds of the points of the first layer, the second layer and the third layer inside the boundary of the calculation region, and the boundary setting of other physical quantities adopts U-i=UiAnd i is 1,2 … N, namely, the virtual points of the first layer, the second layer and the third layer outside the boundary of the calculation region are respectively equal to the physical quantities of the points of the first layer, the second layer and the third layer inside the boundary of the calculation region.
And step 3: and (4) according to the stability condition solved by the hyperbolic equation, taking a calculation parameter CFL and calculating the time step length.
And 4, step 4: determining a self-adaptive algorithm judgment rule, and judging the numerical simulation algorithm types adopted in different areas in the calculation domain.
Step 4.1: the high-precision finite volume calculation method has inherent good conservation characteristic and better capture effect on contact discontinuity, and the high-precision finite difference calculation method has better capture effect on shock waves. According to the advantages of the two high-precision numerical calculation methods, in the numerical calculation process, a high-precision finite volume calculation method is adopted in a calculation region near an interface, and a high-precision finite difference calculation method is adopted in a calculation region of shock wave propagation. Therefore, according to the self-adaptive coupling requirements of two high-precision algorithms and the characteristics of the algorithms, an algorithm type judgment rule is specified:
Figure BDA0003311229960000051
where Δ d represents the distance from the point to the interface, dx, dy, and dz represent the lengths of the grid in the x, y, and z directions, respectively, and α is a constant greater than 3 according to the requirements of the number of nodes in the reconstruction grid in the WENO format and the mud format, and is preferably selected to be 10.
Step 4.2: and (4) judging and marking the algorithm types of different areas in the calculation domain according to the self-adaptive algorithm judgment rule obtained in the step (4.1).
Because the multi-substance interfaces are all located in the calculation area of the finite volume method, the finite volume method needs to process the multi-substance independently, the decoupling is carried out on the single substance for calculation in a unified format, the multi-substance interfaces are not involved in the finite difference calculation area, and the problem of the multi-substance is not needed to be considered.
And 5: in order to avoid non-physical oscillation in multi-substance interface calculation, the multi-substance interface is processed, and meanwhile, to avoid the defects that the GFM method is poor in conservation and the like, a boundary condition multi-substance processing method is adopted, and by applying different boundary conditions to the multi-substance interface, the interaction of multiple substances is decoupled into a single substance problem.
Step 5.1: aiming at fluid-fluid interaction, the processing of a multi-material interface is realized by setting balance boundary conditions on the interface, and a one-dimensional control equation is realized for a finite volume method
Figure BDA0003311229960000052
Wherein U represents a conservation variable, F represents a flux, and when a Level Set symbol distance function
Figure BDA0003311229960000053
Namely, the multi-material interface is positioned at the right side of a grid node i, and the initial condition of establishing two sides of the interface is Ui、Ui+1The multi-substance Riemann problem is solved by a multi-substance Riemann problem solving method, preferably by an HLLC method
Figure BDA0003311229960000054
Setting an equilibrium boundary condition at the interface, equation (7) may become:
Figure BDA0003311229960000055
similarly, when the Level Set symbol distance function
Figure BDA0003311229960000056
Equation (7) may be changed to:
Figure BDA0003311229960000057
step 5.2: the treatment of the multi-material interface is realized by setting immersion boundary conditions on the interface aiming at the fluid-rigid body interaction, and a one-dimensional problem is also taken as an example for explanation. For a one-dimensional control equation of a finite volume method, when a Level Set symbolic distance function
Figure BDA0003311229960000058
Namely, a multi-material interface is positioned at the right side of a grid node i, the grid node i is fluid, the grid node i +1 is a rigid body, the Riemann problem of wall surfaces is established at two sides of the interface, and the initial condition is
Figure BDA0003311229960000061
Wherein u iswFor rigid motion speed, the multi-substance Riemann problem solving method, preferably HLLC method, is adopted to obtain
Figure BDA0003311229960000062
Setting the immersion boundary condition on the interface, the equation (7) can be converted into the equation (8), and similarly, when the Level Set symbolic distance function
Figure BDA0003311229960000063
Establishing a Riemann problem on the wall surfaces at two sides of the interface, wherein the initial conditions are as follows:
Figure BDA0003311229960000064
solving the problem by multi-substance Riemann method, preferably HLLC method
Figure BDA0003311229960000065
The formula (7) is converted to the formula (9) by setting an immersion boundary condition at the interface.
Step 6: and in the finite volume calculation region, reconstructing by adopting a 3-order MUSCL format, solving a multi-material Riemann problem by adopting a preferred HLLC method, obtaining flux, and performing space dispersion on the Euler equation set.
And 7: in a finite difference calculation area, a numerical flux is obtained by adopting a high-precision WENO format, and the Euler equation set is subjected to space dispersion.
And 8: and (4) coupling the euler equation sets obtained in the steps 6 and 7 in a space discrete mode by adopting a self-adaptive algorithm coupling method, and carrying out integral time advancing.
Step 8.1: integrating the spatial discrete results of the Euler equation sets obtained in the steps 6 and 7 to obtain an ordinary differential equation set of the conservation variable with respect to the time derivative in a uniform form in the whole calculation region:
Figure BDA0003311229960000066
where F, G, H are the x, y, z directional fluxes, respectively, and the subscripts indicate the grid positions.
Step 8.2: in order to prevent flux inconsistency of coupling areas of a finite volume method and a finite difference method and further cause non-physical oscillation, virtual grids are arranged at the position of an interface of the two methods for calculating the areas, and a multi-material Riemann problem solving method, preferably an HLLC method, is adopted in the directions x, y and z to obtain the fluxes in the directions x, y and z.
Step 8.3: the time derivative term of the formula (17) is dispersed, preferably in a third-order TVD Runge-Kutta format, and the equal-sign right-end term of the formula (17) is marked as S (U)i,j,k) The propulsion time step is Δ t, so that a fully discrete format of a conservation-oriented euler control equation set is obtained:
Figure BDA0003311229960000071
and step 9: a level-set function is advanced to obtain tn+1The Level-Set function of each substance area at the moment is obtained as tn+1The position of the interface of each substance at the moment.
Step 10: t obtained according to step 9n+1At the moment, the position of each material interface is determined, and t is obtained in the step 8n+1Updating and integrating physical states of all substances at the moment to obtain tn+1The physical state of the entire computing area at the moment.
Since no virtual fluid is used during the multi-material interface process at step 5, the physical state near the interface needs to be updated by solving the multi-material riemann problem. When the interface is advanced from the time t to the time t + delta t, grid nodes with changed material types exist, the physical state of the grid nodes is updated by establishing and solving a local multi-material Riemann problem, and preferably by adopting an HLLC method.
Step 11: judging the current calculation time tn+1Whether the set end time t is exceededendIf t isn+1>tendIf so, ending the numerical calculation and outputting the physical state of the whole calculation area at the moment; if tn+1<tendAnd returning to the step 3, and continuing to perform the multi-substance interaction high-precision numerical simulation.
Further comprising the step 12: and (3) carrying out self-adaptive finite volume-finite difference coupling high-precision numerical simulation by utilizing the steps 1 to 11, capturing shock wave discontinuity and contact discontinuity at high resolution, improving the precision of numerical simulation prediction in the multi-substance interaction process, and further solving the engineering technical problem related to the multi-substance interaction field.
The fields of multi-substance interaction comprise high-speed/ultrahigh-speed penetration and protection, aerospace navigation and mechanical engineering fields.
Predicting a multi-substance interaction process in the engineering field, wherein the multi-substance interaction process comprises explosive air explosion, near-ground explosion, deep water explosion, near-water surface explosion, bubble jet, energy-gathering jet, ultra-high speed collision problem, supernova explosion and Rayleigh Taylor instability problem.
Has the advantages that:
1. the high-precision self-adaptive finite volume-finite difference coupling numerical simulation method disclosed by the invention can capture shock wave discontinuity and contact discontinuity with high resolution, keep the constancy of the calculation process, improve the numerical simulation precision, and reduce the error with the experimental result, so as to accurately predict the problems of explosive air explosion, near-ground explosion, deep water explosion, near-water surface explosion, bubble jet, energy-gathering jet, ultra-high speed collision, supernova explosion, Rayleigh Taylor instability and the like in the engineering field.
2. The invention discloses a high-precision self-adaptive finite volume-finite difference coupling numerical simulation method, which adopts a self-adaptive high-precision finite volume-finite difference coupling calculation method, adaptively selects proper calculation methods in different calculation regions according to self-adaptive algorithm judgment rules, adopts the self-adaptive algorithm coupling method to couple the calculation regions of different high-precision algorithms, can realize high-resolution capture of shock waves, also considers good conservation characteristics, improves the precision of numerical simulation calculation, improves the goodness of numerical simulation and test, and has obvious advantages in processing the problem of multi-substance interaction in the engineering field.
3. The invention discloses a high-precision self-adaptive finite volume-finite difference coupling numerical simulation method, which aims to better perform numerical simulation on complex cavitation generation, development and collapse phenomena, simultaneously avoid non-physical phenomena of negative pressure and density in a numerical simulation process and improve program robustness.
4. The invention discloses a high-precision self-adaptive finite volume-finite difference coupling numerical simulation method, which is characterized in that different boundary conditions such as a balance boundary condition and an immersion boundary condition are applied to an interface, then the physical state near the interface is updated after the interface is pushed by solving a multi-substance Riemann problem, so that the processing of a fluid-fluid multi-substance interface, a fluid-rigid body multi-substance interface and the like is realized, the defect of local non-conservation of the multi-substance interface method is overcome, the calculation efficiency is improved, and the use of virtual fluid methods such as RGFM, MGFM, WGFM and the like is avoided.
Drawings
FIG. 1 is a schematic diagram of pressure gradient boundary condition settings;
FIG. 2 is a schematic diagram of a finite volume method and finite difference method coupling;
FIG. 3 is a schematic diagram of a coupling method of the FVM and the FDM virtual grid nodes;
FIG. 4 is a schematic diagram of the update of the physical state of the interface propulsion of multiple substances;
FIG. 5 is a schematic view of an initial calculation model (pressure cloud) in example 1;
FIG. 6 is a graph comparing numerical simulation and experiment of bubbles formed by underwater explosion at 20ms in example 1;
FIG. 7 is a graph comparing numerical simulation and experiment of bubble morphology at various times in example 1;
FIG. 8 is a comparison graph of pressure curve numerical simulation and experiment at pressure monitoring points at different times in example 1;
fig. 9 discloses a flow of an adaptive finite volume-finite difference coupling numerical simulation method in embodiment 1.
Detailed Description
Example 1:
for better illustrating the objects and advantages of the present invention, the present invention will be further described with reference to the accompanying drawings and examples.
The embodiment is applied to underwater explosion bubble jet calculation.
The initial calculation model in the calculation example disclosed in this embodiment is presented in fig. 5, and the comparison between the numerical simulation of bubbles formed by underwater explosion at 18ms and the experiment is presented in fig. 6, as shown in fig. 9, the embodiment discloses a specific implementation method of an adaptive finite volume-finite difference coupling numerical simulation method, which is:
step 1: according to the problem requirement, a two-dimensional axisymmetric model is adopted for calculation, and the calculation area is 1.6 x 1.6m2The number of grids is 1600 × 1600 ═ 2560000.
Step 2: the calculation area is divided into two parts of explosive and water, the type of the explosive is PETN, the mass is 2.5g, the coordinates of the circle center are (0m,1m), the pressure monitoring point is located in the horizontal direction of the initial position of the explosive and is 1.4m away from the explosive, and accordingly a level-set function is established.
A conservation-type Euler control equation set is adopted as a control equation,
Figure BDA0003311229960000091
wherein the content of the first and second substances,
Figure BDA0003311229960000092
e is the total energy of unit mass of the material, E is the specific internal energy of the material,
Figure BDA0003311229960000093
g is the acceleration of gravity.
The water adopts a Tammann equation of state, an improved Schmidt cavitation model is adopted, the explosive adopts an JWL equation of state, and the initial states of the water and the explosive and the parameters of the equation of state take the values as shown in the following table:
table 1:
Figure BDA0003311229960000101
and applying a non-reflection pressure gradient boundary condition at the boundary position according to the actual problem requirement.
And step 3: the CFL parameter takes 0.3 and calculates the time step.
Figure BDA0003311229960000102
Wherein Δ x is the grid step length in the x direction in the calculation region; Δ y is the grid step length in the y direction in the calculation region; Δ z is the grid step length in the z direction in the calculation region; u, v and w are velocity components of the current grid in the x, y and z directions in the calculation region; c is the speed of sound of the current grid in the calculation region.
And 4, step 4: according to the type of the interface position judgment algorithm judged by the level-set function, taking alpha as 10, namely, in the distance interface
Figure BDA0003311229960000103
And a finite volume algorithm is adopted in the distance, and a finite difference algorithm is adopted in the rest positions of the calculation domain.
And 5: the explosive-water multi-substance interface is processed in the x direction and the y direction respectively, and the multi-substance interaction is decoupled into a single-substance problem by applying a balance boundary condition on the explosive-water interface.
Step 6: calculating the finite volume calculation region by adopting a 3-order MUSCL format combined with a Van Albada limiter to obtain an ordinary differential equation system of the conservation variable of the calculation region with respect to the time derivative:
Figure BDA0003311229960000104
and 7: calculating the finite difference calculation region by adopting a 5-order WENO format to obtain an ordinary differential equation set of the conservation variable of the calculation region with respect to the time derivative:
Figure BDA0003311229960000105
and 8: and (4) integrating the ordinary differential equation sets of the conservation variables obtained in the steps (6) and (7) with respect to the time derivative to obtain the ordinary differential equation sets of the time derivative in a unified form in the whole calculation region, and dispersing the time derivative items by adopting a third-order TVD Runge-Kutta format to obtain a full-discrete format of the conservation-type Euler control equation set.
And step 9: a Level-Set function is advanced to obtain tn+1The Level-Set function of each substance area at the moment is obtained as tn+1The position of the interface of each substance at the moment.
Adopting HJ-WENO format and TVD Runge-Kutta format to respectively disperse the spatial derivative term and the time derivative term of the Level-Set motion equation to obtain the fully-discrete format of the Level-Set motion equation,
Figure BDA0003311229960000111
and then obtain tn+1The Level-Set function of each material area at the moment.
Step 10: t obtained according to step 9n+1At the moment, the position of each material interface is determined, and t is obtained in the step 8n+1Updating and integrating physical states of all substances at any moment, solving the multi-substance Riemann problem through an HLLC method, and updating the physical state near the interface to obtain tn+1The physical state of the entire computing area at the moment.
And then, carrying out steps 11 to 12 to obtain a numerical simulation result of the underwater explosion bubble jet.
And (4) analyzing a calculation result:
fig. 6 is a numerical simulation and experiment comparison diagram of bubbles formed by underwater explosion at 18ms, and it can be seen from the diagram that after bubbles are formed by explosion of explosives, as the bubbles gradually expand, the interface error gradually increases due to poor conservation of identity in the finite difference method, and the comparison error between the adaptive finite volume-finite difference coupling numerical simulation method and the experiment disclosed by the patent is smaller. The figure (7) shows the comparison of the numerical simulation of the bubble shape at each moment with the experiment, and it can be seen from the figure that the bubble gradually shrinks just after the bubble gradually increases, and the bubble forms a bubble jet stream at the later stage of shrinkage due to the gravity. The attached figure (8) is a comparison graph of numerical simulation and experiment of a pressure curve of a pressure monitoring point, and it can be seen from the graph that the explosive is detonated to form shock waves, pressure jump occurs at the pressure monitoring point, the numerical simulation result of the pressure peak value of the initial shock wave is 7.1MPa, the experiment result is 7.3MPa, and the error is less than 3%; at the end of the first period of the bubbles, bubble jet flow is formed and generates a second shock wave, the numerical simulation result of the peak pressure of the shock wave is 1.7MPa, the test result is 1.9MPa, and the error is less than 10%; and at the end of the second period of the bubbles, bubble jet flow is formed again and a third shock wave is generated, the numerical simulation result of the peak pressure of the shock wave is 1.05MPa, the test result is 1.15MPa, and the error is less than 10%. The precision and the effectiveness of the invention are verified by comparing numerical simulation with experimental results.
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 (9)

1. A self-adaptive finite volume-finite difference coupling high-precision numerical simulation method is characterized by comprising the following steps: comprises the following steps of (a) carrying out,
step 1: determining a calculation region according to actual physical problem requirements, establishing a Cartesian coordinate system, and dividing the calculation region into m, n and l grids in the directions of x, y and z respectively, wherein the total number of the m, n and l grids is m;
step 2: defining a level-set function in a calculation area, distinguishing the interface position of the initial state of each substance, defining the initial physical quantity, a material state equation and a conservation type Euler control equation set of each substance according to the area divided by the level-set, and simultaneously setting a boundary condition according to an actual physical problem;
and step 3: according to the stability condition solved by the hyperbolic equation, taking a calculation parameter CFL and calculating the time step length;
and 4, step 4: determining a self-adaptive algorithm judgment rule, and judging the numerical simulation algorithm types adopted in different areas in a calculation domain;
and 5: in order to avoid non-physical oscillation in the calculation of a multi-substance interface, the multi-substance interface is processed, and meanwhile, to avoid the defects of poor conservation and the like of GFM methods, a boundary condition multi-substance processing method is adopted, and by applying different boundary conditions on the multi-substance interface, the interaction of multiple substances is decoupled into a single substance problem;
step 6: reconstructing in a finite volume calculation region by adopting a 3-order MUSCL format, solving a multi-substance Riemann problem, preferably adopting an HLLC method to obtain flux, and performing space dispersion on an Euler equation set;
and 7: in a finite difference calculation area, a numerical flux is obtained by adopting a high-precision WENO format, and the Euler equation set is subjected to space dispersion;
and 8: coupling the euler equation sets obtained in the step 6 and the step 7 in a space discrete mode by adopting a self-adaptive algorithm coupling method, and carrying out integral time advancing;
and step 9: a level-set function is advanced to obtain tn+1The Level-Set function of each substance area at the moment is obtained as tn+1The position of each material interface at the moment;
step 10: t obtained according to step 9n+1At the moment, the position of each material interface is determined, and t is obtained in the step 8n+1Updating and integrating physical states of all substances at the moment to obtain tn+1The physical state of the whole calculation area is obtained at the moment;
step 11: judging the current calculation time tn+1Whether the set end time t is exceededendIf t isn+1>tendIf so, ending the numerical calculation and outputting the physical state of the whole calculation area at the moment; if tn+1<tendAnd returning to the step 3, and continuing to perform the multi-substance interaction high-precision numerical simulation.
2. The adaptive finite volume-finite difference coupling high-precision numerical simulation method of claim 1, wherein: further comprising the step 12: and (3) carrying out self-adaptive finite volume-finite difference coupling high-precision numerical simulation by utilizing the steps 1 to 11, capturing shock wave discontinuity and contact discontinuity at high resolution, improving the precision of numerical simulation prediction in the multi-substance interaction process, and further solving the engineering technical problem related to the multi-substance interaction field.
3. The adaptive finite volume-finite difference coupling high-precision numerical simulation method of claim 2, wherein: the multi-substance interaction field comprises the fields of high-speed/ultrahigh-speed penetration and protection, aerospace navigation and mechanical engineering;
predicting a multi-substance interaction process in the engineering field, wherein the multi-substance interaction process comprises explosive air explosion, near-ground explosion, deep water explosion, near-water surface explosion, bubble jet, energy-gathering jet, ultra-high speed collision problem, supernova explosion and Rayleigh Taylor instability problem.
4. The adaptive finite volume-finite difference coupling high-precision numerical simulation method of claim 1, wherein: the step 2 is realized by the method that,
step 2.1: determining the position of a multi-material interface in an initial state according to an actual physical problem, and establishing a level-set function;
step 2.2: dividing each substance area according to the step 2.1, and defining initial state physical quantities for substances in different areas;
the initial state physical quantity comprises density rho of each substance in the initial state, and velocities u, v and w in three directions of pressure p, x, y and z;
step 2.3: constructing a conservation-oriented Euler control equation set;
step 2.4: selecting a state equation function of the material aiming at the actual physical problem, and further determining the physical state of the material under deformation and motion;
the conservation type Euler control equation set comprises 5 equations and 6 unknowns, and in order to seal the equation set, a proper material state equation needs to be selected to seal the equation set;
step 2.5: setting boundary conditions for the calculation areas divided in the step 2.1;
according to the actual physical problem, a plurality of layers of virtual points are required to be arranged outside a calculation domain to realize boundary conditions of different types;
for non-reflection boundary conditions, such as deep water area and air area, the influence of gravity is considered, so that the pressure gradient boundary condition, p, needs to be set-i=p0+ρg(h-i-h0) N, i is 1,2, i is that the virtual point pressure of the first layer, the second layer and the third layer outside the boundary of the calculation region is equal to the pressure p of the first layer in the boundary of the calculation region0Plus the pressure ρ g (h) due to gravity-i-h0) Boundary setting of other physical quantities adopts U-i=U0I is 1,2 … N, that is, virtual points of the first layer, the second layer and the third layer outside the boundary of the calculation region are equal to the physical quantity of the first layer inside the boundary of the calculation region;
for reflecting boundary conditions, e.g. walls, rigid floors, using un,-i=-un,iI is 1,2 … N, namely the normal speeds of the virtual points of the first layer, the second layer and the third layer outside the boundary of the calculation region are respectively equal to the reverse directions of the normal speeds of the points of the first layer, the second layer and the third layer inside the boundary of the calculation region, and the boundary setting of other physical quantities adopts U-i=UiI 1,2 … N, i.e. the first layer, the second layer outside the boundary of the calculation regionThe virtual points of the second layer and the third layer are respectively equal to the physical quantities of the points of the first layer, the second layer and the third layer in the boundary of the calculation region.
5. The adaptive finite volume-finite difference coupling high-precision numerical simulation method of claim 4, wherein: without considering the chemical reaction, JWL equation of state was chosen for the explosive and expressed as:
Figure RE-FDA0003466903800000031
where ρ is0A, B, R1, R2, omega are JWL equation of state parameters for the initial density of the explosive;
an ideal gas state equation is selected for air and expressed as
p=(γ-1)ρe (2)
Wherein γ is the ideal gas adiabatic index;
tammann equation of state, expressed as
p=(γ-1)ρe-γ(B-A) (3)
Wherein, gamma and A, B are Tammann equation of state parameters;
in a low-pressure area formed by sparse waves, because the pressure is less than saturated vapor pressure, cavitation phenomenon can occur in water, in order to carry out more accurate numerical simulation on the cavitation phenomenon, a cavitation model is introduced, and preferably, an improved Schmidt model in a single-fluid cavitation model is selected and expressed as
Figure RE-FDA0003466903800000032
Wherein p issatIs the saturated vapor pressure, pminAt a very small positive pressure, pgAnd ρlDensity of gas and water, respectively, cgAnd clAcoustic velocity, p, of gas and water, respectivelyglIs obtained by the following formula:
Figure RE-FDA0003466903800000041
6. the adaptive finite volume-finite difference coupling high-precision numerical simulation method of claim 5, wherein: step 4, the method is realized by the following steps,
step 4.1: the high-precision finite volume calculation method has inherent good conservation characteristic and better capture effect on contact discontinuity, and the high-precision finite difference calculation method has better capture effect on shock waves; according to the advantages of the two high-precision numerical calculation methods, in the numerical calculation process, a high-precision finite volume calculation method is adopted in a calculation region near an interface, and a high-precision finite difference calculation method is adopted in a calculation region of shock wave propagation; therefore, according to the self-adaptive coupling requirements of two high-precision algorithms and the characteristics of the algorithms, an algorithm type judgment rule is specified:
Figure RE-FDA0003466903800000042
wherein, Δ d represents the distance from the point to the interface, dx, dy, dz represent the lengths of the grid in the x, y, z directions respectively, and α is a constant greater than 3 according to the requirements of reconstructing the number of grid nodes in WENO format and MUSCL format;
step 4.2: judging and marking the algorithm types of different areas in the calculation domain according to the self-adaptive algorithm judgment rule obtained in the step 4.1;
because the multi-substance interfaces are all located in the calculation area of the finite volume method, the finite volume method needs to process the multi-substance independently, the decoupling is carried out on the single substance for calculation in a unified format, the multi-substance interfaces are not involved in the finite difference calculation area, and the problem of the multi-substance is not needed to be considered.
7. The adaptive finite volume-finite difference coupling high-precision numerical simulation method of claim 6, wherein: step 5 the method is realized by the following steps,
step 5.1: aiming at fluid-fluid interaction, the processing of a multi-material interface is realized by setting balance boundary conditions on the interface, and a one-dimensional control equation is realized for a finite volume method
Figure RE-FDA0003466903800000043
Wherein U represents a conservation variable, F represents a flux, and when a Level Set symbol distance function
Figure RE-FDA0003466903800000044
Namely, the multi-material interface is positioned at the right side of a grid node i, and the initial condition of establishing two sides of the interface is Ui、Ui+1The multi-substance Riemann problem is solved by adopting a multi-substance Riemann problem solving method
Figure RE-FDA0003466903800000051
Setting an equilibrium boundary condition at the interface, equation (7) may become:
Figure RE-FDA0003466903800000052
similarly, when the Level Set symbol distance function
Figure RE-FDA0003466903800000053
Equation (7) may be changed to:
Figure RE-FDA0003466903800000054
step 5.2: aiming at the interaction of fluid and rigid body, the processing of the multi-material interface is realized by setting immersion boundary conditions on the interface, and the one-dimensional problem is also taken as an example for explanation; for a one-dimensional control equation of a finite volume method, when a Level Set symbolic distance function
Figure RE-FDA0003466903800000055
Namely, a multi-material interface is positioned at the right side of a grid node i, the grid node i is fluid, the grid node i +1 is a rigid body, the Riemann problem of wall surfaces is established at two sides of the interface, and the initial condition is
Figure RE-FDA0003466903800000056
Wherein u iswFor rigid motion speed, the multi-substance Riemann problem solving method, preferably HLLC method, is adopted to obtain
Figure RE-FDA0003466903800000057
Setting the immersion boundary condition on the interface, the equation (7) can be converted into the equation (8), and similarly, when the Level Set symbolic distance function
Figure RE-FDA0003466903800000058
Establishing a Riemann problem on the wall surfaces at two sides of the interface, wherein the initial conditions are as follows:
Figure RE-FDA0003466903800000059
solving the problem by means of multi-substance Riemann
Figure RE-FDA00034669038000000510
The formula (7) is converted to the formula (9) by setting an immersion boundary condition at the interface.
8. The adaptive finite volume-finite difference coupling high-precision numerical simulation method of claim 7, wherein: step 8 is realized by a method comprising the following steps of,
step 8.1: integrating the spatial discrete results of the Euler equation sets obtained in the steps 6 and 7 to obtain an ordinary differential equation set of the conservation variable with respect to the time derivative in a uniform form in the whole calculation region:
Figure RE-FDA0003466903800000061
wherein F, G, H are the x, y, z flux respectively, and the subscripts denote the grid positions;
step 8.2: in order to prevent the flux inconsistency of the coupling areas of the finite volume method and the finite difference method and further cause the non-physical oscillation, virtual grids are arranged at the interface positions of the calculation areas of the two methods, a multi-material Riemann problem solving method is adopted in the x direction, the y direction and the z direction respectively, and the HLLC method is preferably adopted to obtain the fluxes in the x direction, the y direction and the z direction;
step 8.3: the time derivative term of the formula (17) is dispersed, and the equal sign right end term of the formula (17) is marked as S (U)i,j,k) The propulsion time step is Δ t, so that a fully discrete format of a conservation-oriented euler control equation set is obtained:
Figure RE-FDA0003466903800000062
9. the adaptive finite volume-finite difference coupling high-precision numerical simulation method of claim 8, wherein: the step 10 is realized by that since the virtual fluid is not used in the multi-substance interface processing process in the step 5, the physical state near the interface needs to be updated by solving the multi-substance riemann problem; when the interface is advanced from the time t to the time t + delta t, grid nodes with changed material types exist, the physical state of the grid nodes is updated by establishing and solving a local multi-material Riemann problem, and the physical state of the grid nodes is updated by establishing the local multi-material Riemann problem.
CN202111217254.5A 2021-10-19 High-precision self-adaptive finite volume-finite difference coupling numerical simulation method Active CN114169184B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111217254.5A CN114169184B (en) 2021-10-19 High-precision self-adaptive finite volume-finite difference coupling numerical simulation method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111217254.5A CN114169184B (en) 2021-10-19 High-precision self-adaptive finite volume-finite difference coupling numerical simulation method

Publications (2)

Publication Number Publication Date
CN114169184A true CN114169184A (en) 2022-03-11
CN114169184B CN114169184B (en) 2024-06-21

Family

ID=

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114840974A (en) * 2022-03-25 2022-08-02 中国气象局地球系统数值预报中心 Advection mode system suitable for complex terrain and operation method thereof
CN115329646A (en) * 2022-10-12 2022-11-11 国家超级计算天津中心 Shock wave simulation method, device, equipment and medium
CN116205152A (en) * 2022-12-12 2023-06-02 中广核风电有限公司 Numerical simulation method and device for offshore floating fan
CN116663373A (en) * 2023-07-25 2023-08-29 中国科学技术大学 High-precision numerical simulation method suitable for gas phase combustion and explosion

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20130231907A1 (en) * 2010-11-23 2013-09-05 Yahan Yang Variable Discretization Method For Flow Simulation On Complex Geological Models
CN109902376A (en) * 2019-02-25 2019-06-18 北京理工大学 A kind of fluid structurecoupling high resolution numerical simulation method based on Continuum Mechanics
CN112861263A (en) * 2021-02-22 2021-05-28 西北工业大学 Calculation simulation method suitable for compressible two-phase flow
CN112861445A (en) * 2020-12-29 2021-05-28 中国船舶重工集团公司第七研究院 Simulation method of grid-free numerical model
CN113408168A (en) * 2021-06-18 2021-09-17 北京理工大学 High-precision numerical simulation method based on Riemann problem accurate solution

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20130231907A1 (en) * 2010-11-23 2013-09-05 Yahan Yang Variable Discretization Method For Flow Simulation On Complex Geological Models
CN109902376A (en) * 2019-02-25 2019-06-18 北京理工大学 A kind of fluid structurecoupling high resolution numerical simulation method based on Continuum Mechanics
CN112861445A (en) * 2020-12-29 2021-05-28 中国船舶重工集团公司第七研究院 Simulation method of grid-free numerical model
CN112861263A (en) * 2021-02-22 2021-05-28 西北工业大学 Calculation simulation method suitable for compressible two-phase flow
CN113408168A (en) * 2021-06-18 2021-09-17 北京理工大学 High-precision numerical simulation method based on Riemann problem accurate solution

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
王成;SHU CHI-WANG;: "爆炸力学高精度数值模拟研究进展", 科学通报, no. 10, 10 April 2015 (2015-04-10) *
赵海涛;王成;宁建国;: "可压缩多介质问题的高精度数值模拟", 高压物理学报, no. 02, 15 April 2013 (2013-04-15) *

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114840974A (en) * 2022-03-25 2022-08-02 中国气象局地球系统数值预报中心 Advection mode system suitable for complex terrain and operation method thereof
CN115329646A (en) * 2022-10-12 2022-11-11 国家超级计算天津中心 Shock wave simulation method, device, equipment and medium
CN115329646B (en) * 2022-10-12 2023-03-10 国家超级计算天津中心 Shock wave simulation method, device, equipment and medium
CN116205152A (en) * 2022-12-12 2023-06-02 中广核风电有限公司 Numerical simulation method and device for offshore floating fan
CN116205152B (en) * 2022-12-12 2024-06-07 中广核风电有限公司 Numerical simulation method and device for offshore floating fan
CN116663373A (en) * 2023-07-25 2023-08-29 中国科学技术大学 High-precision numerical simulation method suitable for gas phase combustion and explosion
CN116663373B (en) * 2023-07-25 2023-11-17 中国科学技术大学 High-precision numerical simulation method suitable for gas phase combustion and explosion

Similar Documents

Publication Publication Date Title
CN109948301B (en) Near-water surface sliding jump fluid-solid coupling numerical value prediction method based on grid control
Zhang et al. Nonlinear interaction between underwater explosion bubble and structure based on fully coupled model
Farhat et al. Robust and provably second‐order explicit–explicit and implicit–explicit staggered time‐integrators for highly non‐linear compressible fluid–structure interaction problems
Wang et al. A computational framework for the simulation of high‐speed multi‐material fluid–structure interaction problems with dynamic fracture
Liu et al. Investigation of hydrodynamics of water impact and tail slamming of high-speed water entry with a novel immersed boundary method
Saksono et al. An adaptive remeshing strategy for flows with moving boundaries and fluid–structure interaction
CN110795869B (en) Numerical calculation method and device for flow field data
CN110516342B (en) Propeller compressible cavitation flow numerical prediction method based on OpenFOAM platform
Wang et al. Simulations of the dynamics and interaction between a floating structure and a near-field explosion bubble
CN111709196B (en) Cavitation erosion resistance assessment method for underwater high-speed navigation body
Muralidharan et al. Simulation of moving boundaries interacting with compressible reacting flows using a second-order adaptive Cartesian cut-cell method
CN115758569A (en) Fluid-solid coupling simulation method for high-speed water-entering slamming and cavitation of aircraft
Li et al. Experimental study of underwater explosions below a free surface: Bubble dynamics and pressure wave emission
Cayzac et al. Computational fluid dynamics and experimental validations of the direct coupling between interior, intermediate and exterior ballistics using the Euler equations
CN111695282A (en) Liquid tank sloshing prediction and control analysis method based on fluid-solid coupling simulation
Liu et al. The modified ghost fluid method applied to fluid-elastic structure interaction
CN114169184A (en) High-precision self-adaptive finite volume-finite difference coupling numerical simulation method
CN114169184B (en) High-precision self-adaptive finite volume-finite difference coupling numerical simulation method
Zhou et al. Design of scaled model for dynamic characteristics of stiffened cylindrical shells based on equivalent similar method
CN113591345B (en) Explosion reaction flow high-precision prediction method based on generalized Riemann solver
Zheng et al. Comparative study of different SPH schemes on simulating violent water wave impact flows
KR101910742B1 (en) High Order Accurate Meshless Method for solving the compressible flow analysis
CN113408168B (en) High-precision numerical simulation method based on Riemann problem accurate solution
Zarchi et al. Computational Fluid Dynamics Validation and Post-test Analysis of Supersonic Retropropulsion in the Ames 9x7 Unitary Tunnel
Balabel RANS modeling of gas jet impinging onto a deformable liquid interface

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