CN111859645A - Improved MUSL format material dot method for shock wave solving - Google Patents

Improved MUSL format material dot method for shock wave solving Download PDF

Info

Publication number
CN111859645A
CN111859645A CN202010657277.7A CN202010657277A CN111859645A CN 111859645 A CN111859645 A CN 111859645A CN 202010657277 A CN202010657277 A CN 202010657277A CN 111859645 A CN111859645 A CN 111859645A
Authority
CN
China
Prior art keywords
flow field
momentum
shock wave
grid
material point
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
CN202010657277.7A
Other languages
Chinese (zh)
Other versions
CN111859645B (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.)
Nanjing University of Science and Technology
Original Assignee
Nanjing University of Science and Technology
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 Nanjing University of Science and Technology filed Critical Nanjing University of Science and Technology
Priority to CN202010657277.7A priority Critical patent/CN111859645B/en
Publication of CN111859645A publication Critical patent/CN111859645A/en
Application granted granted Critical
Publication of CN111859645B publication Critical patent/CN111859645B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/08—Fluids
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

The invention provides an improved MUSL format object particle method for solving shock waves, which comprises the steps of firstly improving the momentum correction step in the traditional MUSL format; then establishing a shock wave flow field model, and carrying out discretization treatment on the shock wave flow field; then giving initial conditions and boundary conditions of a shock wave flow field, and setting a total calculation length; then calculating the time step of the current calculation step based on the stability condition; solving a shock wave flow field by an improved MUSL format material dot method based on shock wave solving; and finally, performing visual processing on the impact wave flow field. Compared with the traditional MUSL format, the method can effectively reduce the noise when the object particles pass through the grid, has lower requirement on the precision of a shape function, stronger robustness and stability, higher calculation efficiency and less calculation energy consumption, and can provide a new method for solving shock waves in engineering.

Description

Improved MUSL format material dot method for shock wave solving
Technical Field
The invention belongs to the technical field of shock wave flow field solving methods, and particularly relates to an improved MUSL-format substance particle method for solving shock waves.
Background
Shock waves are discrete peaks that propagate in a fluid medium, resulting in a step-like change in physical properties of the flow field, such as pressure, temperature, density, etc. It is widely existed in engineering, such as train shock wave, mine blasting, explosion-proof vehicle, etc. when high-speed train runs, and may cause casualties and equipment damage. With the development of computer technology, the CFD technology is gradually applied to the simulation of shock waves, and the disturbance effect of the shock waves on the flow field is obtained by numerically calculating the propagation process of the shock waves, so that calculation parameters can be provided for the calculation of the aerodynamic load of an acting object of the shock waves. In addition, the shock wave problem is one of the classical examples of CFD, the parameters such as density, pressure and the like of a wave front flow field and a wave rear flow field are discontinuous, and the solution of the discontinuous problems is always a difficult point and a core problem of CFD development. Therefore, the novel shock wave solving method has wide engineering value.
The physical point method is a numerical calculation method mixing the lagrange method and the eulerian method, and has been widely applied to solving the high-speed impact problem and the explosion problem in the solid structure. The method uses material points to disperse objects, related physical parameters are carried by the material points, and a background grid is only used for integrating a momentum equation and solving a gradient, so that a structural grid covering a calculation domain is only needed to be simply generated. The material point method gets rid of the constraint of the grid to a certain extent, avoids the grid distortion problem caused by large deformation, can well treat the phenomena of material fracture failure, structural breakage and the like, does not cause mass loss, and is suitable for solving extreme working conditions. The MUSL format is a calculation format with momentum correction, in the traditional MUSL format in the material point method, the material point momentum of the next time step is mapped to the grid nodes based on the shape function of the current time step to obtain the corrected grid node momentum of the next time, and the grid node speed is calculated by using the corrected grid node momentum and the grid node quality before correction. Since the position of the material point is changed at this time, the shape function of the material point should be changed accordingly, but the change is not considered in the conventional MUSL format, when the material point passes through the grid, the correction method generates a large error, which affects the calculation accuracy and even causes the calculation to be not converged.
Disclosure of Invention
The invention aims to provide an improved MUSL format matter dot method for solving shock waves, and provides a method for predicting the influence of the shock waves in practical engineering.
The technical solution for realizing the purpose of the invention is as follows:
an improved MUSL-format material dot method for solving shock waves comprises the following steps:
step 1, improving the momentum correction step in the traditional MUSL format: mapping the material point momentum of the next time step to the grid nodes based on the shape function of the current time step to obtain the corrected grid node momentum of the next time, canceling the momentum correction step in the traditional MUSL format, calculating the shape function based on the material point position of the next time step, and mapping the material point mass and the momentum of the next time step to the grid nodes;
step 2, establishing a shock wave flow field model, and carrying out discretization treatment on the shock wave flow field: dividing a background grid and distributing material points;
step 3, setting initial conditions and boundary conditions of a shock wave flow field, and setting a total calculation length;
step 4, calculating the time step length of the current calculation step based on the stability condition;
step 5, solving a shock wave flow field by an improved MUSL format object particle method based on shock wave solving;
Step 6, performing visualization processing on the impact wave flow field: and outputting the density and pressure parameters of the flow field.
A shock wave solving system based on an improved MUSL format material dot method comprises a material dot momentum correction module, a shock wave flow field model establishing module, a total calculation length setting module, a time step length calculating module, a shock wave flow field calculating module and a visualization processing module;
the material point momentum correction module is used for mapping the material point momentum of the next time step to the grid nodes based on the shape function of the current time step to obtain the corrected grid node momentum of the next time, canceling the momentum correction step in the traditional MUSL format, calculating the shape function based on the material point position of the next time step, and simultaneously mapping the material point mass and the momentum of the next time step to the grid nodes; the shock wave flow field model establishing module is used for establishing a shock wave flow field model and carrying out discretization processing on a shock wave flow field; the total length calculating setting module is used for setting the total length of calculation under the initial condition and the boundary condition of a given shock wave flow field; the time step calculation module is used for calculating the time step of the current calculation step; the shock wave flow field calculation module is used for solving a shock wave flow field based on an improved MUSL format material particle method for solving shock waves; the visual processing module is used for carrying out visual processing on the impact wave flow field: and outputting the density and pressure parameters of the flow field.
Compared with the prior art, the invention has the following remarkable advantages:
(1) the shock wave problem is solved by using a material point method, and a regular structure grid covering a calculation domain is generated, so that the workload and time of grid division are greatly reduced; the mass of a single object point is always unchanged in calculation, the mass conservation law is automatically met, and a mass equation does not need to be solved; physical parameters such as pressure, density and the like are stored on object points, so that the numerical dissipation of the object points is reduced; the method adopts the Lagrange method to solve in a single calculation step, avoids the problem of difficult processing of the convection term in the Eulerian method, abandons a deformed background grid between two calculation steps, reestablishes the mapping relation between the object points and the grid, and avoids the grid distortion problem in the Lagrange method.
(2) The correction method improves the momentum correction step in the traditional MUSL format, avoids the defect that the shape function in the traditional MUSL format lags behind the node momentum by one time step, corrects the grid node momentum, corrects the grid node quality, reduces the error caused by the material point crossing the grid, can more accurately capture the discontinuous characteristics before and after the shock wave surface, can well inhibit oscillation, and can improve the calculation convergence speed.
Drawings
FIG. 1 is a flow chart of a method for solving shock wave material points according to the present invention.
Fig. 2 is a schematic of the Riemann problem.
Fig. 3 is a diagram of the results of solving the Riemann problem in accordance with the present invention and the conventional MUSL format.
Fig. 4 is a diagram comparing the analytic solutions of the Riemann problem of the present invention and the conventional MUSL format.
Detailed Description
The invention is further described with reference to the following figures and embodiments.
Step 1, improving a momentum correction step in a traditional MUSL format, mapping the material point momentum of the next time step to grid nodes based on a shape function of the current time step to obtain the corrected grid node momentum of the next time, canceling a momentum correction step in the traditional MUSL format, calculating the shape function based on the material point position of the next time step, and simultaneously mapping the material point mass and the momentum of the next time step to the grid nodes to obtain the corrected grid node mass and momentum as follows:
Figure BDA0002577215680000031
Figure BDA0002577215680000032
wherein
Figure BDA0002577215680000033
For the quality of the mesh node i in the next computation step,
Figure BDA0002577215680000034
for the momentum of grid node i in the next computation step, mpIs the mass of a point of matter,
Figure BDA0002577215680000035
is the velocity of the material point p in the next calculation step,
Figure BDA0002577215680000036
is the mapping function between the mesh node i and the material point p in the next calculation step.
Step 2, establishing a shock wave flow field model, carrying out discretization treatment on the shock wave flow field, carrying out geometric cleaning on the boundary of the flow field and establishing topology, covering the whole calculation domain with a uniform structural grid to generate a background grid, generating at least one layer of grid outside the boundary as a virtual grid, and uniformly distributing object points in the calculation range of the shock wave.
Step 3, giving initial conditions and boundary conditions of a shock wave flow field, setting and calculating total length, and giving density of each material point at initial time
Figure BDA0002577215680000037
Pressure of
Figure BDA0002577215680000038
Speed of rotation
Figure BDA0002577215680000039
Calculating the volume of each material point at the initial moment
Figure BDA00025772156800000310
Mass mpInternal energy of
Figure BDA0002577215680000041
Figure BDA0002577215680000042
Figure BDA0002577215680000043
Figure BDA0002577215680000044
Wherein VΩIs the volume of the entire shock wave flow field, npThe total number of material points, gamma is the specific heat ratio of gas, and the air is 1.4.
And applying boundary conditions on the background grids, wherein for a solid wall boundary, the speed, momentum and nodal force of the grids at the boundary and the virtual grids outside the boundary are always 0.
Step 4, calculating the time step of the current calculation step based on the stability condition, and setting the number of the database, wherein the time step of the current calculation step is calculated as follows:
Figure BDA0002577215680000045
Figure BDA0002577215680000046
where Δ ttFor the time step of the current calculation step, CCFLIs the number of the Kuran number,
Figure BDA0002577215680000047
is the sound velocity of the material point p at the present moment,
Figure BDA0002577215680000048
respectively the speed of the material point p in the x, y and z directions in the current calculation step,
Figure BDA0002577215680000049
Is the density of the material point p in the current calculation step,
Figure BDA00025772156800000410
l is the grid length, which is the pressure of the material point p in the current calculation step.
Step 5, solving the shock wave flow field based on an improved MUSL format substance particle method for solving shock waves, which comprises the following steps:
5.1: mapping the material point parameters to grid nodes to obtain the mass and momentum of the grid nodes as follows:
Figure BDA00025772156800000411
Figure BDA00025772156800000412
wherein
Figure BDA00025772156800000413
For the quality of mesh node i in the current computation step,
Figure BDA00025772156800000414
for the momentum of mesh node i in the current computation step,
Figure BDA0002577215680000051
is the velocity of the mass point p in the current calculation step,
Figure BDA0002577215680000052
is a mapping function between the grid node i and the material point p in the current calculation step.
5.2: the grid node forces are calculated as:
Figure BDA0002577215680000053
wherein f isi tFor the force value of the mesh node i in the current computation step,
Figure BDA0002577215680000054
the derivative of the mapping function between the mesh node i and the material point p in the current computation step.
5.3: the integral momentum equation is:
Figure BDA0002577215680000055
wherein
Figure BDA0002577215680000056
Is a mesh node i isThe momentum in the preceding calculation step is calculated,
Figure BDA0002577215680000057
the momentum in the next computation step for mesh node i.
5.4: calculating the position and the speed of the material point at the next moment as follows:
Figure BDA0002577215680000058
Figure BDA0002577215680000059
wherein
Figure BDA00025772156800000510
Is the velocity of the mass point p in the current calculation step,
Figure BDA00025772156800000511
for the velocity of the material point p in the next calculation step, niThe number of the nodes of the grid is,
Figure BDA00025772156800000512
For the position of the object point p in the current calculation step,
Figure BDA00025772156800000513
the position of the material point p in the next calculation step is calculated.
5.5: the step of correcting the momentum of the grid nodes in the traditional MUSL format is cancelled, the mass and the momentum of the object points are mapped to the grid nodes again, and the mass and the momentum of the grid nodes are corrected as follows:
Figure BDA00025772156800000514
Figure BDA00025772156800000515
5.6, the density of the updated material points is:
Figure BDA0002577215680000061
Figure BDA0002577215680000062
Figure BDA0002577215680000063
wherein
Figure BDA0002577215680000064
The strain increase of the material point p in the next calculation step is a square matrix with rows j and k,
Figure BDA0002577215680000065
is the density of the material point p in the current calculation step,
Figure BDA0002577215680000066
the density of the material point p in the next calculation step is calculated.
5.7: the pressure of the mass point is updated using the ideal gas equation of state as:
Figure BDA0002577215680000067
wherein
Figure BDA0002577215680000068
The internal energy of the material point p in the next calculation step is calculated.
5.8: and judging whether the time is calculated to the end point, if not, making t equal to t +1 and transferring to the step 4, calculating the time step again, and if so, ending the step and entering the next step.
And 6, visualizing the shock wave flow field, storing time-course data of parameters such as density on grid nodes and the like, and obtaining the propagation process of the shock wave in the flow field.
The invention also provides an improved MUSL-format object particle method and a computer program based on the shock wave solving, and the shock wave solving system also comprises an object particle momentum correction module, a shock wave flow field model establishing module, a total length calculating module, a time step calculating module, a shock wave flow field calculating module and a visualization processing module;
The material point momentum correction module is used for mapping the material point momentum of the next time step to the grid nodes based on the shape function of the current time step to obtain the corrected grid node momentum of the next time, canceling the momentum correction step in the traditional MUSL format, calculating the shape function based on the material point position of the next time step, and simultaneously mapping the material point mass and the momentum of the next time step to the grid nodes, wherein the specific process is shown in the step 1; the shock wave flow field model establishing module is used for establishing a shock wave flow field model and carrying out discretization treatment on the shock wave flow field, and the specific process is shown in the step 2; the total length calculating setting module is used for setting the total length of calculation according to the initial condition and the boundary condition of the given shock wave flow field, and the specific process is shown in the step 3; the time step calculation module is used for calculating the time step of the current calculation step, and the specific process is shown in the step 4; the shock wave flow field calculation module is used for solving the shock wave flow field based on an improved MUSL format material dot method for solving shock waves, and the specific process is shown in the step 5; the visual processing module is used for carrying out visual processing on the impact wave flow field: outputting the density and pressure parameters of the flow field, wherein the specific process is shown in the step 6; the program processing process is the same as the processing process of the method, so that the method is not described in detail.
Example 1
The process of the shock wave solving method does not specify any specific object, namely, the method is applicable to solving any shock wave propagation problem in engineering. In order to explain the operation flow and application effect of the method in detail, the application of the invention is explained in detail by taking the classic Riemann problem as an example. The Riemann problem, namely the decomposition problem of the initial discontinuity, is one of the classical arithmetic examples of CFD, which comprises various discontinuities in CFD, and the solution of the discontinuity problems is always a difficult and core problem in the development of CFD. The Riemann problem has an analytic solution, and can be used for checking the correctness and the precision of a numerical simulation method, and therefore almost all flow field solving methods are established on the basis of solving the Riemann problem and then are expanded towards multiple dimensions.
The Riemann problem describes a closed tube with different gases on the left and right sides, the gases on the two sides being separated by a thin film in the middle, and at the beginning of the calculation, the diaphragm in the middle is suddenly removed, the gases on the two sides will be mixed with each other under the action of pressure, and finally, an equilibrium state is reached, which is also called a shock tube problem. The total calculation time of the calculation example is t equal to 0.2, all parameters in the calculation example are dimensionless parameters, a virtual grid needs to be generated at a boundary when a B spline mapping function is used, the total number of background grids generated in the present example is 1002, 1 virtual grid is respectively arranged outside the left end and the right end of a shock tube, each grid body at the initial moment except the virtual grid contains 2 object points, so the number of the object points is 2000, the number of the library is set to 0.8, the gas specific heat ratio gamma is 1.4, a calculation domain is x e (0,1), and the intermittent distribution at the initial moment is as follows:
Figure BDA0002577215680000071
In order to verify the robustness of the method, a 2-time B-spline basis function (3-point interpolation) with lower precision and a 3-time B-spline basis function (4-point interpolation) with higher precision are respectively adopted as shape functions between the material points and the grid nodes, and the calculation result of the method is compared with the analytic solutions of the traditional MUSL format and Riemann problems.
Step 1, improving a momentum correction step in a traditional MUSL format, mapping the material point momentum of the next time step to grid nodes based on a shape function of the current time step to obtain the corrected grid node momentum of the next time, canceling a momentum correction step in the traditional MUSL format, calculating the shape function based on the material point position of the next time step, and simultaneously mapping the material point mass and the momentum of the next time step to the grid nodes to obtain the corrected grid node mass and momentum as follows:
Figure BDA0002577215680000072
Figure BDA0002577215680000081
wherein
Figure BDA0002577215680000082
For the quality of the mesh node i in the next computation step,
Figure BDA0002577215680000083
for the momentum of grid node i in the next computation step, mpIs the mass of a point of matter,
Figure BDA0002577215680000084
is the velocity of the material point p in the next calculation step,
Figure BDA0002577215680000085
is the mapping function between the mesh node i and the material point p in the next calculation step.
Step 2, establishing a shock tube model, regarding the shock tube model as a one-dimensional line model with the length of 1, carrying out discretization treatment on a flow field in the shock tube, uniformly generating 1000 grids in the tube, respectively arranging 1 virtual grid outside the left end point and the right end point of the tube, setting the total number of background grids to be 1002, and setting 2 material points in each grid except the virtual grid, wherein the total number of the material points is 2000.
Step 3, giving initial conditions and boundary conditions of the shock tube, and setting the density, the speed and the pressure of material points in the tube at the initial moment as follows:
Figure BDA0002577215680000086
the volume, mass and internal energy of each material point are calculated as follows:
Figure BDA0002577215680000087
Figure BDA0002577215680000088
Figure BDA0002577215680000089
here VΩ＝1，np＝2000。
And applying boundary conditions to the left end and the right end of the shock tube to enable the speed, momentum and force of the end points and the grid nodes outside the end points to be 0 all the time.
And 4, calculating the time step according to the stability condition as follows:
Figure BDA00025772156800000810
Figure BDA00025772156800000811
where Δ ttFor the time step of the current calculation step, CCFL＝0.8，
Figure BDA0002577215680000091
Is the sound velocity of the material point p at the present moment,
Figure BDA0002577215680000092
is the velocity of the mass point p in the x-direction in the current calculation step,
Figure BDA0002577215680000093
is the density of the material point p in the current calculation step,
Figure BDA0002577215680000094
is the pressure of the mass point p in the current calculation step.
Step 5, solving the shock wave flow field based on an improved MUSL format substance particle method for solving shock waves, which comprises the following steps:
5.1: mapping the material point parameters to grid nodes to obtain the mass and momentum of the grid nodes as follows:
Figure BDA0002577215680000095
Figure BDA0002577215680000096
wherein
Figure BDA0002577215680000097
For the quality of mesh node i in the current computation step,
Figure BDA0002577215680000098
for the momentum of grid node i in the current computation step, mpIs the mass of the point of matter p,
Figure BDA0002577215680000099
is the velocity of the mass point p in the current calculation step,
Figure BDA00025772156800000910
is a mapping function between the grid node i and the material point p in the current calculation step.
5.2: the grid node forces are calculated as:
Figure BDA00025772156800000911
wherein f isi tFor the force value of the mesh node i in the current computation step,
Figure BDA00025772156800000912
the derivative of the mapping function between the mesh node i and the material point p in the current computation step.
5.3: the integral momentum equation is:
Figure BDA00025772156800000913
wherein
Figure BDA00025772156800000914
For the momentum of mesh node i in the current computation step,
Figure BDA00025772156800000915
the momentum in the next computation step for mesh node i.
5.4: calculating the position and the speed of the material point at the next moment as follows:
Figure BDA00025772156800000916
Figure BDA00025772156800000917
wherein
Figure BDA0002577215680000101
Is the velocity of the mass point p in the current calculation step,
Figure BDA0002577215680000102
for the velocity of the material point p in the next calculation step, niThe number of the nodes of the grid is,
Figure BDA0002577215680000103
for the position of the object point p in the current calculation step,
Figure BDA0002577215680000104
the position of the material point p in the next calculation step is calculated.
5.5: the step of correcting the momentum of the grid nodes in the traditional MUSL format is cancelled, the mass and the momentum of the object points are mapped to the grid nodes again, and the mass and the momentum of the grid nodes are corrected as follows:
Figure BDA0002577215680000105
Figure BDA0002577215680000106
5.6, the density of the updated material points is:
Figure BDA0002577215680000107
Figure BDA0002577215680000108
Figure BDA0002577215680000109
wherein
Figure BDA00025772156800001010
The strain increase of the material point p in the next calculation step is a square matrix with rows j and k,
Figure BDA00025772156800001011
is the density of the material point p in the current calculation step,
Figure BDA00025772156800001012
the density of the material point p in the next calculation step is calculated.
5.7: the pressure of the mass point is updated using the ideal gas equation of state as:
Figure BDA00025772156800001013
Wherein
Figure BDA00025772156800001014
The internal energy of the material point p in the next calculation step is calculated.
5.8: and judging whether the calculation time reaches 0.2, if not, turning to the step 4, and if so, ending the next step.
And 6, performing visual processing on the flow field in the shock tube, and outputting the density of each grid node to obtain the density distribution characteristic of the flow field.
When a 2-time B-spline function is taken as a shape function, the calculation result is shown in FIG. 3, and it can be seen that the difference between the calculation result of the conventional MUSL format and the analytic solution is large, the dispersion and dissipation of the density curve are large, and the result is severely oscillated; the calculation result and the analytic solution of the method of the invention are very close, and the traditional MUSL format needs about 1220 steps of calculation to complete the calculation, and the method of the invention only needs 1040 steps. When a 3-order B spline function is taken as a shape function, the calculation result is shown in figure 4, along with the improvement of the precision of a mapping function, the calculation precision of the MUSL format and the traditional MUSL format is improved, but the dispersion and the dissipation of the MUSL format are smaller, particularly at the discontinuous part, the density curve of the MUSL format is closer to an analytical solution, and the convergence speed is higher.
In summary, under the same shape function, the calculation speed and the calculation accuracy of the method are higher than those of the traditional MUSL format; the method has better robustness and lower requirement on the precision of the shape function, can obtain better calculation results under a 3-point interpolation format, and can better capture shock wave characteristics under a 4-point interpolation format in the traditional MUSL format, thereby indicating that the method has lower calculation energy consumption.

Claims (10)

1. An improved MUSL-format material dot method for solving shock waves is characterized by comprising the following steps of:
step 1, improving the momentum correction step in the traditional MUSL format: mapping the material point momentum of the next time step to the grid nodes based on the shape function of the current time step to obtain the corrected grid node momentum of the next time, canceling the momentum correction step in the traditional MUSL format, calculating the shape function based on the material point position of the next time step, and mapping the material point mass and the momentum of the next time step to the grid nodes;
step 2, establishing a shock wave flow field model, and carrying out discretization treatment on the shock wave flow field: dividing a background grid and distributing material points;
step 3, setting initial conditions and boundary conditions of a shock wave flow field, and setting a total calculation length;
step 4, calculating the time step length of the current calculation step based on the stability condition;
step 5, solving a shock wave flow field by an improved MUSL format object particle method based on shock wave solving;
step 6, performing visualization processing on the impact wave flow field: and outputting the density and pressure parameters of the flow field.
2. The improved MUSL format object point method for shockwave resolution as recited in claim 1, wherein step 1 improves the momentum correction step in the conventional MUSL format, and the obtained grid node quality and momentum after correction are:
Figure FDA0002577215670000011
Figure FDA0002577215670000012
Wherein
Figure FDA0002577215670000013
For the quality of the mesh node i in the next computation step,
Figure FDA0002577215670000014
for the momentum of grid node i in the next computation step, mpIs the mass of a point of matter,
Figure FDA0002577215670000015
is the velocity of the material point p in the next calculation step,
Figure FDA0002577215670000016
is the mapping function between the mesh node i and the material point p in the next calculation step.
3. The improved MUSL-format object point method for shock wave solving as recited in claim 1, wherein step 2 is to create a shock wave flow field model, discretize the shock wave flow field, geometrically clean the boundaries of the flow field and create topology, cover the whole computational domain with a uniform structural grid to create a background grid, create at least one layer of grid outside the boundaries as a virtual grid, and uniformly lay object points within the computational range of the shock wave.
4. The method of claim 1, wherein step 3 provides initial conditions and boundary conditions of the shockwave flow field, sets the total length of the calculation, and provides the density of each material point at the initial time
Figure FDA0002577215670000017
Pressure of
Figure FDA0002577215670000018
Speed of rotation
Figure FDA0002577215670000019
Calculating the volume of each material point at the initial moment
Figure FDA00025772156700000110
Mass m of matter pointpInternal energy of
Figure FDA0002577215670000021
Figure FDA0002577215670000022
Figure FDA0002577215670000023
Figure FDA0002577215670000024
Wherein VΩIs the volume of the entire shock wave flow field, npIs the total number of material points, and gamma is the specific heat ratio of gas;
And applying boundary conditions on the background grids, wherein for the solid wall boundary, the speed, momentum and nodal force of the grids at the boundary and the virtual grids outside the boundary are always 0.
5. The improved MUSL-format material point method for shockwave resolution as recited in claim 1, wherein step 4 is performed by calculating the time step of the current calculation step based on the stability condition, and the set Curian number, wherein the time step of the current calculation step is calculated as follows:
Figure FDA0002577215670000025
Figure FDA0002577215670000026
where Δ ttFor the time step of the current calculation step, CCFLIs the number of the Kuran number,
Figure FDA0002577215670000027
is the sound velocity of the material point p at the present moment,
Figure FDA0002577215670000028
respectively the speed of the material point p in the x, y and z directions in the current calculation step,
Figure FDA0002577215670000029
is the density of the material point p in the current calculation step,
Figure FDA00025772156700000210
is the pressure of the material point p in the current calculation step, gamma is the gas specific heat ratio, L is the grid length, npIs the total number of material dots.
6. The improved MUSL format material particle method for shockwave resolution as recited in claim 1, wherein step 5 is for solving the shockwave flow field based on the improved MUSL format material particle method for shockwave resolution, comprising:
step 5.1, mapping the material point parameters to grid nodes to obtain the quality and momentum of the grid nodes as follows:
Figure FDA00025772156700000211
Figure FDA0002577215670000031
wherein
Figure FDA0002577215670000032
For the quality of mesh node i in the current computation step,
Figure FDA0002577215670000033
For the momentum of mesh node i in the current computation step,
Figure FDA0002577215670000034
is the velocity of the mass point p in the current calculation step,
Figure FDA0002577215670000035
for a mapping function between grid node i and material point p in the current computation step, mpIs the mass of a material point;
step 5.2, calculating the grid node force as follows:
Figure FDA0002577215670000036
wherein f isi tFor the force value of the mesh node i in the current computation step,
Figure FDA0002577215670000037
the derivative of the mapping function between the grid node i and the material point p in the current calculation step is obtained;
and 5.3, an integral momentum equation is as follows:
Figure FDA0002577215670000038
wherein
Figure FDA0002577215670000039
For the momentum of mesh node i in the current computation step,
Figure FDA00025772156700000310
the momentum of the grid node i in the next calculation step;
and 5.4, calculating the position and the speed of the material point at the next moment:
Figure FDA00025772156700000311
Figure FDA00025772156700000312
wherein
Figure FDA00025772156700000313
Is the velocity of the mass point p in the current calculation step,
Figure FDA00025772156700000314
for the velocity of the material point p in the next calculation step, niThe number of the nodes of the grid is,
Figure FDA00025772156700000315
for the position of the object point p in the current calculation step,
Figure FDA00025772156700000316
the position of the material point p in the next calculation step;
and 5.5, canceling the step of correcting the momentum of the grid nodes in the traditional MUSL format, and mapping the mass and the momentum of the object points to the grid nodes again, wherein the step of correcting the mass and the momentum of the grid nodes is as follows:
Figure FDA00025772156700000317
Figure FDA0002577215670000041
and 5.6, updating the density of the object points as follows:
Figure FDA0002577215670000042
Figure FDA0002577215670000043
Figure FDA0002577215670000044
Wherein
Figure FDA0002577215670000045
The strain increase of the material point p in the next calculation step is a square matrix with rows j and k,
Figure FDA0002577215670000046
is the density of the material point p in the current calculation step,
Figure FDA0002577215670000047
the density of the material point p in the next calculation step;
and 5.7, updating the pressure of the material point by using an ideal gas state equation to be as follows:
Figure FDA0002577215670000048
wherein
Figure FDA0002577215670000049
The internal energy of the material point p in the next calculation step is shown, and gamma is the specific heat ratio of the gas;
and 5.8, judging whether the terminal time is calculated or not, if not, making t equal to t +1 and turning to the step 4, and if so, ending the next step.
7. The improved MUSL-format material dot method for solving shock waves according to claim 1, wherein step 6 is used for visualizing the shock wave flow field, storing time-course data of parameters such as density on grid nodes and the like, and obtaining the propagation process of the shock waves in the flow field.
8. A shock wave solving system based on an improved MUSL format material dot method is characterized by comprising a material dot momentum correction module, a shock wave flow field model establishing module, a total calculation length setting module, a time step length calculating module, a shock wave flow field calculating module and a visualization processing module;
the material point momentum correction module is used for mapping the material point momentum of the next time step to the grid nodes based on the shape function of the current time step to obtain the corrected grid node momentum of the next time, canceling the momentum correction step in the traditional MUSL format, calculating the shape function based on the material point position of the next time step, and simultaneously mapping the material point mass and the momentum of the next time step to the grid nodes; the shock wave flow field model establishing module is used for establishing a shock wave flow field model and carrying out discretization processing on a shock wave flow field; the total length calculating setting module is used for setting the total length of calculation under the initial condition and the boundary condition of a given shock wave flow field; the time step calculation module is used for calculating the time step of the current calculation step; the shock wave flow field calculation module is used for solving a shock wave flow field based on an improved MUSL format material particle method for solving shock waves; the visual processing module is used for carrying out visual processing on the impact wave flow field: and outputting the density and pressure parameters of the flow field.
9. The shock wave solving system of claim 8, wherein the modified grid node masses and momentums processed by the object particle momentum modification module are:
Figure FDA0002577215670000051
Figure FDA0002577215670000052
wherein
Figure FDA0002577215670000053
For the quality of the mesh node i in the next computation step,
Figure FDA0002577215670000054
for the momentum of grid node i in the next computation step, mpIs the mass of a point of matter,is the velocity of the material point p in the next calculation step,
Figure FDA0002577215670000056
is the mapping function between the mesh node i and the material point p in the next calculation step.
10. The shockwave solution system of claim 8, wherein the step-size computation module processes the current computation step to obtain a step size of time:
Figure FDA0002577215670000057
Figure FDA0002577215670000058
where Δ ttFor the current calculation stepStep of time of, CCFLIs the number of the Kuran number,
Figure FDA0002577215670000059
is the sound velocity of the material point p at the present moment,
Figure FDA00025772156700000510
respectively the speed of the material point p in the x, y and z directions in the current calculation step,
Figure FDA00025772156700000511
is the density of the material point p in the current calculation step,
Figure FDA00025772156700000512
is the pressure of the material point p in the current calculation step, gamma is the gas specific heat ratio, L is the grid length, npIs the total number of material dots.
CN202010657277.7A 2020-07-09 2020-07-09 Improved MUSL format material dot method for shock wave solving Active CN111859645B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010657277.7A CN111859645B (en) 2020-07-09 2020-07-09 Improved MUSL format material dot method for shock wave solving

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010657277.7A CN111859645B (en) 2020-07-09 2020-07-09 Improved MUSL format material dot method for shock wave solving

Publications (2)

Publication Number Publication Date
CN111859645A true CN111859645A (en) 2020-10-30
CN111859645B CN111859645B (en) 2023-03-31

Family

ID=73151958

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010657277.7A Active CN111859645B (en) 2020-07-09 2020-07-09 Improved MUSL format material dot method for shock wave solving

Country Status (1)

Country Link
CN (1) CN111859645B (en)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112364362A (en) * 2020-11-16 2021-02-12 宁波九寰适创科技有限公司 Parallel multilayer self-adaptive local encryption method facing fluid simulation direction
CN113673185A (en) * 2021-08-25 2021-11-19 北京理工大学 Method for accurately capturing shock wave discontinuous surface through weighted bidirectional mapping
CN113701981A (en) * 2021-09-14 2021-11-26 佛山奇正电气有限公司 Near-wall motion shock wave identification method

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20100250205A1 (en) * 2009-03-31 2010-09-30 Airbus Espana, S.L. Method and system for a quick calculation of aerodynamic forces on an aircraft in transonic conditions
CN107766287A (en) * 2017-10-26 2018-03-06 哈尔滨工程大学 A kind of Stochastic Dynamics analysis method based on thing particle method being applied in blast impulse engineering

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20100250205A1 (en) * 2009-03-31 2010-09-30 Airbus Espana, S.L. Method and system for a quick calculation of aerodynamic forces on an aircraft in transonic conditions
CN107766287A (en) * 2017-10-26 2018-03-06 哈尔滨工程大学 A kind of Stochastic Dynamics analysis method based on thing particle method being applied in blast impulse engineering

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112364362A (en) * 2020-11-16 2021-02-12 宁波九寰适创科技有限公司 Parallel multilayer self-adaptive local encryption method facing fluid simulation direction
CN112364362B (en) * 2020-11-16 2023-12-29 宁波九寰适创科技有限公司 Parallel multi-layer self-adaptive local encryption method oriented to fluid simulation direction
CN113673185A (en) * 2021-08-25 2021-11-19 北京理工大学 Method for accurately capturing shock wave discontinuous surface through weighted bidirectional mapping
CN113673185B (en) * 2021-08-25 2024-01-30 北京理工大学 Accurate shock wave discontinuities capturing method based on weighted bidirectional mapping
CN113701981A (en) * 2021-09-14 2021-11-26 佛山奇正电气有限公司 Near-wall motion shock wave identification method

Also Published As

Publication number Publication date
CN111859645B (en) 2023-03-31

Similar Documents

Publication Publication Date Title
CN111859645B (en) Improved MUSL format material dot method for shock wave solving
CN108446445B (en) Composite material wing optimization design method based on aerodynamic reduced order model
Blom A monolithical fluid-structure interaction algorithm applied to the piston problem
CN108563843B (en) Method for updating disturbance area of steady compressible flow
Mavriplis Directional agglomeration multigrid techniques for high-Reynolds-number viscous flows
CN113673027B (en) Agent model-based hypersonic aircraft pneumatic load optimization design method
CN116151084B (en) Simulation method and device based on structural grid, terminal equipment and storage medium
CN116245049B (en) Node type non-structural grid boundary correction method, device, equipment and medium
CN110795869B (en) Numerical calculation method and device for flow field data
CN108897716B (en) Data processing device and method for reducing calculation amount through memory read-write operation
CN113609597B (en) Method for updating time-space hybrid propulsion disturbance domain of supersonic aircraft streaming
CN115618498B (en) Prediction method, device, equipment and medium for cross-basin flow field of aircraft
CN114492250A (en) Curved surface mesh generation method and system based on recursive decomposition and computer equipment
CN116738891A (en) LU-SGS improvement method for enhancing simulation stability of aircraft flow field
CN113361173B (en) Compressible modification method of transition model completely based on local flow field parameters
CN111027250A (en) Special-shaped curved surface reinforced shell modeling method based on grid deformation technology
CN105718619A (en) Method for determining fuel quality characteristics of aircraft based on finite element method
CN112307669B (en) Multi-medium level set method based on actual physical point iterative correction
CN113378440A (en) Fluid-solid coupling numerical simulation calculation method, device and equipment
CN116341421B (en) Hypersonic flow field numerical simulation method, hypersonic flow field numerical simulation system, electronic equipment and storage medium
Baumeister Time-dependent difference theory for noise propagation in a two-dimensional duct
CN111859646A (en) Shock wave variable step length solving method based on B spline mapping function material point method
JP2017224299A (en) Optimal pressure-projection method for incompressible transient and steady-state navier-stokes equations
Roy et al. A finite volume method for viscous incompressible flows using a consistent flux reconstruction scheme
CN110909511B (en) Non-viscous low-speed streaming numerical simulation method without curved surface volume division

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