Numerical simulation method of lander low-speed impact on small celestial body's gravel surface

Abstract

The application provides a kind of dynamics numerical simulation method for lander low-speed impact small celestial body gravel surface, which comprises: determining the state information of lander and irregular particles, constructing the hierarchical bounding box of lander and irregular particles;With the help of linear hierarchical bounding box root node, coarse-grained contact detection is carried out, and with the help of linear hierarchical bounding box, fine-grained contact detection is carried out, and the relative intrusion depth of lander and irregular particles is calculated;Based on the contact force and contact torque between lander and irregular particles, combined with the translational dynamics equation and rotational dynamics equation of lander and irregular particles, the state of lander and irregular particles is updated by means of second-order leapfrog integration.The present application can solve the problem of multi-point contact between irregular particles and lander by constructing Delaunay triangular mesh model in view of the shortcomings of spherical particle discrete element hypothesis, and is more in line with the physical process of mutual contact between irregular particles, and is more efficient and reliable.

Description

Technical Field

[0001] This invention belongs to the field of numerical simulation analysis technology of particle dynamics in interstellar space, and in particular, a numerical simulation method for the dynamics of a lander impacting a small celestial body's gravel surface at low speed. Background Technology

[0002] Deep space exploration has revealed that most small celestial bodies are composed of fragments ranging from millimeters to tens of meters in size, accreted by their own gravity. This results in irregular shapes, weak gravitational fields, and surfaces covered by layers of mixed gravel particles spanning more than three orders of magnitude. Humanity is no longer content with mere glimpses through telescopes; it yearns to collect samples of material from the surface and interior of these small bodies to provide more data for studies on the evolution of the solar system, the origin of life on Earth, planetary defense, and the development and utilization of small body resources. Therefore, the era of surface contact exploration of small celestial bodies has arrived. NASA, the European Space Agency (ESA), and Japan's JAXA have conducted impact, sampling, rover exploration, and drilling experiments on various types of small celestial bodies, and the scientific results obtained have significantly updated humanity's original understanding of these bodies.

[0003] However, due to insufficient understanding of the gravel layer on small celestial bodies, the exploration mission experienced unexpected events such as lander bounce and sample leakage from the sampling device. Particle systems are inherently complex, with numerous particles in contact, dissipating energy through collisions and friction. Due to their significant internal inhomogeneity, nonlinearity, and discreteness, particle systems exhibit self-organization and high sensitivity, and are highly susceptible to solid-liquid phase transitions under external perturbations. Furthermore, the coupling of the weak, irregular gravitational field of the small celestial body and various complex spatial interactions leads to even more complex dynamic behaviors in the gravel layer. Existing numerical simulation methods mostly employ discrete element method (DEM) simulations of spherical particles to model the interaction between the lander or sampler and the gravel layer, offering advantages such as simple contact detection, ease of parallelization, and high computational efficiency. However, the particles on the surface of small celestial bodies are not perfectly spherical, and most are irregularly shaped. While spherical particles exhibit interlocking and interlocking phenomena during contact due to their non-spherical shape, it is difficult to simulate irregular particles. Therefore, conclusions about spherical particle systems cannot be easily extrapolated to irregular particle systems. Therefore, in order to better serve my country's surface exploration missions, it is urgent and necessary to seek a dynamic numerical simulation method for the low-speed impact of a lander on the gravel surface of a small celestial body, so as to deepen the understanding of the dynamics of the gravel particle layer of the small celestial body. Summary of the Invention

[0004] This invention addresses the shortcomings of existing technologies by proposing a dynamic numerical simulation method for low-speed impact of a lander onto a gravelly surface of a small celestial body. The method includes determining the state information of the lander and irregular particles, constructing a hierarchical bounding box for the lander and irregular particles; performing coarse-grained contact detection using the root node of the linear hierarchical bounding box, performing fine-grained contact detection using the same box, and calculating the relative intrusion depth between the lander and the irregular particles; updating the state of the lander and irregular particles based on the contact force and contact torque between the lander and the irregular particles, combined with the translational and rotational dynamic equations of the lander and the irregular particles, using a second-order leapfrog integral. This invention addresses the shortcomings of the discrete element method assumption for spherical particles by constructing a Delaunay triangular mesh model, which solves the problem of multi-point contact between irregular particles and the lander, better reflects the physical process of contact between irregular particles, and is more efficient and reliable.

[0005] This invention provides a dynamic numerical simulation method for a lander impacting a small celestial body's gravel surface at low speed, comprising the following steps:

[0006] S1. Initialize the state of the lander and irregular particles: Determine the state information of the lander and irregular particles, and construct the hierarchical bounding box of the lander and irregular particles.

[0007] S11. Determine the status information of the lander and irregular particles: Read in the Delaunay triangular mesh model of the lander and irregular particles, and send it from the CPU memory to the GPU memory.

[0008] S12. Based on the Deloitte triangular mesh model, construct a hierarchical bounding box between the lander and the irregular particles;

[0009] S2. Contact detection between the lander and irregular particles: coarse-grained contact detection is performed using the root node of the linear hierarchical bounding box, fine-grained contact detection is performed using the linear hierarchical bounding box, and the relative intrusion depth between the lander and the irregular particles is calculated.

[0010] S21, coarse-grained contact detection, quickly identifies all irregular particles that may have come into contact with the lander;

[0011] S22. Fine-grained contact detection: Using the logic of the thread hierarchy organized by the GPU, determine the leaf nodes of all irregular particle-level bounding boxes that are in contact with the internal nodes of the lander-level bounding box.

[0012] S23. Calculate the contact force and contact torque between the lander and the irregular particles: Based on the triangular facet pairs formed by the leaf nodes of the lander's hierarchical bounding box and the leaf nodes of the irregular particles' hierarchical bounding box within the linear hierarchical bounding box, calculate the contact force and contact torque between the lander and the irregular particles.

[0013] S3. Lander and Irregular Particle State Update: Based on the contact force and contact torque between the lander and the irregular particles, combined with the translational and rotational dynamic equations of the lander and the irregular particles, and using second-order leapfrog integrals, the state of the lander and the irregular particles is updated; the state includes translational motion state and rotational motion state, and the hierarchical bounding box of the lander and the irregular particles is updated according to the position information; the update of the translational motion state is as follows:

[0014]

[0015] Among them, v m v m+12 and v m+1 These represent the lander's velocity at the moment before the update, the first half of the update, and the moment after the update, respectively; a m and a m+1 These represent the lander's acceleration before and after the update, respectively; r m and r m+1 These represent the lander's positions before and after the update, respectively; Δt represents the time interval.

[0016] The update of the rotational motion state is as follows:

[0017]

[0018] Where, ω m and ω m+1 These represent the angular velocities of the lander before and after the update, respectively; q represents the updated angular acceleration of the lander; m and q m+1 Let f(q) represent the quaternions before and after the lander update, respectively; d represents the differential; t represents time; m () represents a differential equation involving quaternions;

[0019] The hierarchical bounding box of the lander or irregular particles is updated as follows:

[0020]

[0021] Among them, V B Represents the original vertex coordinates of the lander or the bounding box of the irregular particle hierarchy; Q represents the quaternion; V N This represents the updated vertex coordinates of the lander or irregular particle hierarchy bounding box.

[0022] Furthermore, step S12 specifically includes the following steps:

[0023] S121. Based on the Delaunay triangular mesh model of the lander and the irregular particles, construct bounding boxes of all levels respectively, and execute steps S122 and S123.

[0024] S122. Obtain the Morton code corresponding to the centroid of the triangular face: convert the decimal number in the coordinates of the centroid into a binary number, and then interleave the binary numbers to obtain the Morton code corresponding to the centroid.

[0025] S123. Based on the obtained Morton coding, construct a hierarchical bounding box:

[0026] S1231. Using the binary search method, determine the first distinct binary bit of the Morton code;

[0027] S1232. Divide the left and right subtrees: Nodes with the highest bit of 0 are assigned to the left subtree and become left nodes, and nodes with the highest bit of 1 are assigned to the right subtree and become right nodes.

[0028] S1233. Based on the left and right subtrees respectively, repeat steps S1311 and S1312 to determine the next different binary bit of the Morton code and divide the left and right subtrees of the next level until the last binary bit of the Morton code is processed, and the hierarchical bounding box is completed.

[0029] Preferably, step S21 specifically includes the following steps:

[0030] S211. Execute step S212 sequentially according to the numbering of the irregular particles.

[0031] S212. Determine whether the maximum bounding box of the corresponding irregular particle overlaps with the maximum bounding box of the lander. If so, store the number of the corresponding irregular particle in the array and execute step S213; otherwise, execute step S213.

[0032] S213. Repeat step S212 until the last irregular particle corresponding to the last number is executed;

[0033] Step S22 specifically includes the following steps:

[0034] S221. If the left and right subtrees of the lander's hierarchical bounding box are both leaf nodes and both intersect with the leaf nodes of the irregular particle's hierarchical bounding box, then determine whether the leaf nodes of the lander's hierarchical bounding box and the leaf nodes of the irregular particle's hierarchical bounding box intersect in the triangular faces within the linear hierarchical bounding box.

[0035] S222. If the left and right subtrees of the lander's hierarchical bounding box are both internal nodes, and only one of the left or right subtrees intersects with the leaf node of the irregular particle's hierarchical bounding box, or if one of the left and right subtrees of the lander's hierarchical bounding box is a leaf node and the other is an internal node, and the internal node intersects with the leaf node of the irregular particle's hierarchical bounding box, then recursively determine whether all internal child nodes of the next level of the internal node of the lander's hierarchical bounding box intersect with the leaf node of the irregular particle's hierarchical bounding box.

[0036] S223. If the left and right subtrees of the lander's hierarchical bounding box are both internal nodes and both intersect with the leaf nodes of the irregular particle's hierarchical bounding box, then the right subtree is temporarily stored, and the two child nodes of the next level of the left subtree are recursively judged to see if they intersect with the leaf nodes of the irregular particle's hierarchical bounding box. This process continues until all child nodes of all levels of the left subtree have been judged. Then, the temporarily stored right subtree is read, and the recursive judgment continues to see if the two child nodes of the next level of the right subtree intersect with the leaf nodes of the irregular particle's hierarchical bounding box.

[0037] Step S23 specifically includes the following steps:

[0038] S231. The lander and the irregular particles are separated by a number of triangular facet pairs, wherein the mutual intrusion depth u of the k-th triangular facet pair is... nk for:

[0039]

[0040] Where, n i This represents the normal vector of the triangle corresponding to the lander; n j T represents the normal vector of the triangular face corresponding to the irregular particle; i and T j These represent the centroids of the triangular facets corresponding to the lander and the triangular facets corresponding to the irregular particles, respectively. and These represent the locations of the centroids of the lander and the irregular particle triangular facet, respectively.

[0041] S232. Calculate the normal contact force F of the k-th pair of triangular facets. cnk for:

[0042] F cnk =k nk u nk n k +C nk v nk (5)

[0043] Where, n k v represents the normal vector of the triangle corresponding to the k-th lander; nkThe relative velocity of the k-th pair of triangular facets is represented by k; nk C represents the stiffness coefficient. nk Indicates the damping coefficient;

[0044] S233. Calculate the tangential contact force F of the k-th pair of triangular facets. csk for:

[0045]

[0046] Among them, u s This represents the relative displacement of the k-th triangular element with respect to the tangential direction; v s The relative velocity of the k-th triangular element with respect to the tangential direction is represented by μ; μ represents the coefficient of friction; k sk C represents the tangential elastic modulus; sk Indicates the tangential damping coefficient;

[0047] S234, F based on each thread cnk and F csk To obtain the contact force and contact torque between the lander and irregular particles.

[0048] Preferably, in step S3, f(q) m ) is represented as:

[0049]

[0050] Where, (ω m ×) represents ω m The antisymmetric matrix;

[0051] V B Represented as:

[0052]

[0053] in, and These represent the x, y, and z components of the original vertex coordinates of the bounding box.

[0054] Q is represented as:

[0055] Q = [q0 q1 q2 q3] T (9)

[0056] Where q0, q1, q2, and q3 represent the first to fourth components of the quaternion, respectively;

[0057] V N Represented as:

[0058]

[0059] in, and These represent the x, y, and z components of the vertex coordinates after the hierarchical bounding box is updated.

[0060] Preferably, in step S12, the linear hierarchical bounding box of the first layer of the hierarchical bounding box is the largest; the second layer of the hierarchical bounding box divides the Delaunay triangular mesh in two, evenly distributing it into the linear hierarchical bounding boxes of the left and right subtrees, and so on downwards until the linear hierarchical bounding box of the last layer of the hierarchical bounding box contains only one triangular face; the linear hierarchical bounding box of the first layer of the hierarchical bounding box is called the root node, the linear hierarchical bounding boxes of the second to penultimate layers of the hierarchical bounding box are called internal nodes, and the linear hierarchical bounding box of the last layer of the hierarchical bounding box is called the leaf node; all Delaunay triangular meshes are contained in the root node, and all nodes except the leaf nodes are collectively called internal nodes.

[0061] Preferably, the logic for organizing the GPU thread hierarchy in step S22 is that the number of thread blocks of the GPU is equal to the number of irregular particles that may be in contact with the lander, and the number of threads allocated in each thread block is the same as the number of leaf nodes in the hierarchical bounding box of the irregular particle; that is, one thread block is responsible for handling the contact detection between one irregular particle and the lander, and one thread is allocated a leaf node of the hierarchical bounding box of the irregular particle.

[0062] Preferably, in step S11, after pre-allocating video memory to the GPU using the memory allocation function cudaMalloc, the status information of the lander and irregular particles on the CPU side is sent to the GPU side using the memory allocation function cudaMalloc.

[0063] Preferably, the state information in step S11 includes position, velocity, angular velocity, and quaternion.

[0064] Compared with the prior art, the technical effects of the present invention are as follows:

[0065] 1. This invention presents a dynamic numerical simulation method for the low-speed impact of a lander on a gravelly surface of a small celestial body. Addressing the shortcomings of the discrete element method assumption for spherical particles, it constructs a Deloitte triangular mesh model, capable of simulating the contact between irregularly shaped particles and the lander, closely resembling the actual surface environment of a gravelly particle layer on a small celestial body. Based on the Deloitte triangular mesh, the contact areas between irregular particles are determined, and the contact force and contact torque between irregular particles are calculated using a distributed continuous contact force model. This method solves the problem of multi-point contact between irregular particles and the lander, better reflects the physical process of contact between irregular particles, and offers higher reliability.

[0066] 2. This invention presents a dynamic numerical simulation method for low-speed impact of a lander onto a gravelly surface of a small celestial body. Utilizing a heterogeneous computing framework combining CPU and GPU, the method performs contact detection between gravel particles and the lander, calculation of contact forces, and state updates of irregular particles and the lander within the GPU. In the most time-consuming contact detection stage, the method divides contact detection into coarse-grained contact detection, fine-grained contact detection, and relative intrusion depth calculation. By rationally organizing the GPU's thread hierarchy, the computational burden on each thread can be reduced, allowing computational tasks to be allocated to a finer parallel granularity, thereby improving program computational efficiency and supporting the simulation of the contact and collision process between the lander and a large-scale irregular particle system. Attached Figure Description

[0067] Other features, objects, and advantages of this application will become more apparent from the following detailed description of non-limiting embodiments with reference to the accompanying drawings.

[0068] Figure 1 This is a flowchart of the dynamic numerical simulation method for low-speed impact of a lander on a gravelly surface of a small celestial body according to the present invention.

[0069] Figure 2 This is a schematic diagram of the key algorithm flow of the present invention;

[0070] Figure 3 This is a demonstration diagram of Morton coding in a specific embodiment of the present invention;

[0071] Figure 4 This is a schematic diagram of constructing a hierarchical bounding box based on Morton coding in a specific embodiment of the present invention;

[0072] Figure 5 This is a schematic diagram of constructing a hierarchical bounding box for irregular particles in a specific embodiment of the present invention;

[0073] Figure 6 This is a schematic diagram of a pair of triangular facets in contact with each other according to the present invention;

[0074] Figure 7 This is a schematic diagram simulating the motion of a cube-shaped lander in a small celestial body gravel layer in a specific embodiment of the present invention. Detailed Implementation

[0075] The present application will now be described in further detail with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and not intended to limit it. Furthermore, it should be noted that, for ease of description, only the parts relevant to the invention are shown in the accompanying drawings.

[0076] Figure 1 and Figure 2This paper presents a numerical simulation method for the dynamics of a lander impacting a small celestial body's gravel surface at low speed, as described in this invention. The gravel surface is assumed to consist of several irregular particles, and low-speed impact typically refers to an impact at a velocity of 0.001-1 m / s. The method includes the following steps:

[0077] S1. Initialize the state of the lander and irregular particles: Determine the state information of the lander and irregular particles, and construct the hierarchical bounding box of the lander and irregular particles.

[0078] S11. Determine the status information of the lander and irregular particles: Read in the Delaunay triangular mesh model of the lander and irregular particles, and send it from the CPU memory to the GPU memory; the status information includes position, velocity, angular velocity and quaternion.

[0079] Data stored on the CPU and GPU cannot be accessed directly; data must be transferred via the API function cudaMemcpy for access to be possible. Therefore, after the CPU reads the position, velocity, angular velocity, quaternions, and other relevant status information of the lander and irregular particles from an external file, as well as the Delaunay triangular mesh model of the lander and irregular particles, it needs to pre-allocate video memory for the GPU using the cudaMalloc function. Then, the cudaMalloc function is used to send the status information of the lander and irregular particles from the CPU to the GPU.

[0080] S12. Based on the Deloitte triangular mesh model, construct a hierarchical bounding box between the lander and the irregular particles.

[0081] S121. Based on the Delaunay triangular mesh models of the lander and irregular particles respectively, construct bounding boxes at all levels, and execute steps S122 and S123.

[0082] S122. Obtain the Morton code corresponding to the centroid of the triangle: Convert the decimal number in the coordinates of the centroid to a binary number, and then interleave the binary numbers to obtain the Morton code corresponding to the centroid.

[0083] In one specific embodiment, such as Figure 3 As shown, when the centroid coordinates of the triangular face are (4,9,3), the binary representation is (0100,1001,0011), and the corresponding Morton code is 010100001011.

[0084] S123. Construct a hierarchical bounding box based on the obtained Morton code.

[0085] S1231. Using the binary search method, determine the first distinct binary bit of the Morton code.

[0086] S1232. Divide the left and right subtrees: Nodes with the highest bit of 0 are assigned to the left subtree and become left nodes, while nodes with the highest bit of 1 are assigned to the right subtree and become right nodes.

[0087] S1233. Based on the left and right subtrees respectively, repeat steps S1311 and S1312 to determine the next different binary bit of the Merton code and divide the left and right subtrees of the next level until the last binary bit of the Merton code is processed, and the hierarchical bounding box is completed. In a specific embodiment, a schematic diagram of constructing a hierarchical bounding box based on the Merton code is shown below. Figure 4 As shown.

[0088] The first layer of the hierarchical bounding box has the largest linear hierarchical bounding box. The second layer of the hierarchical bounding box divides the Delaunay triangular mesh in two, evenly distributing it into the linear hierarchical bounding boxes of the left and right subtrees. This continues downwards until the last layer of the hierarchical bounding box contains only one triangular face. The linear hierarchical bounding box of the first layer is called the root node, the linear hierarchical bounding boxes of the second to penultimate layers are called internal nodes, and the linear hierarchical bounding box of the last layer is called a leaf node. All Delaunay triangular meshes are contained within the root node, and all nodes except the leaf nodes are collectively called internal nodes. In a specific embodiment, taking the construction of a hierarchical bounding box with irregular particles as an example, such as... Figure 5 As shown, the first layer of the hierarchical bounding box is the root node, the second layer is the internal node, and the seventh layer (the last layer) is the leaf node.

[0089] S2. Contact detection between the lander and irregular particles: coarse-grained contact detection is performed using the root node of the linear hierarchical bounding box, fine-grained contact detection is performed using the linear hierarchical bounding box, and the relative intrusion depth between the lander and the irregular particles is calculated.

[0090] S21, coarse-grained contact detection, quickly identifies all irregular particles that may have come into contact with the lander.

[0091] S211. Execute step S212 sequentially according to the numbering of the irregular particles.

[0092] S212. Determine whether the maximum bounding box of the corresponding irregular particle overlaps with the maximum bounding box of the lander. If so, store the number of the corresponding irregular particle in an array and execute step S213; otherwise, execute step S213.

[0093] S213. Repeat step S212 until the last irregular particle corresponding to the last number is reached.

[0094] S22. Fine-grained contact detection: Using the logic of organizing the thread hierarchy on the GPU side, the leaf nodes of all irregular particle bounding boxes that are in contact with the internal nodes of the lander's bounding box are determined. The logic of organizing the thread hierarchy on the GPU is that the number of thread blocks of the GPU is equal to the number of irregular particles that may be in contact with the lander, and the number of threads allocated in each thread block is the same as the number of leaf nodes of the irregular particle's bounding box. That is, one thread block is responsible for handling the contact detection between one irregular particle and the lander, and one thread is allocated a leaf node of the irregular particle's bounding box.

[0095] S221. If the left and right subtrees of the lander's hierarchical bounding box are both leaf nodes, and both intersect with the leaf nodes of the irregular particle's hierarchical bounding box, then determine whether the leaf nodes of the lander's hierarchical bounding box and the leaf nodes of the irregular particle's hierarchical bounding box intersect in the triangular faces within the linear hierarchical bounding box.

[0096] S222. If the left and right subtrees of the lander's hierarchical bounding box are both internal nodes, and only one of the left or right subtrees intersects with the leaf node of the irregular particle's hierarchical bounding box, or if one of the left and right subtrees of the lander's hierarchical bounding box is a leaf node and the other is an internal node, and the internal node intersects with the leaf node of the irregular particle's hierarchical bounding box, then recursively determine whether all internal child nodes of the next level of the internal node of the lander's hierarchical bounding box intersect with the leaf node of the irregular particle's hierarchical bounding box.

[0097] S223. If the left and right subtrees of the lander's hierarchical bounding box are both internal nodes and both intersect with the leaf nodes of the irregular particle's hierarchical bounding box, then temporarily store the right subtree and recursively determine whether the two child nodes of the next level of the left subtree intersect with the leaf nodes of the irregular particle's hierarchical bounding box. After all child nodes of all levels under the left subtree have been determined, read the temporarily stored right subtree and continue to recursively determine whether the two child nodes of the next level of the right subtree intersect with the leaf nodes of the irregular particle's hierarchical bounding box.

[0098] S23. Calculate the contact force and contact torque between the lander and the irregular particles: Based on the triangular facet pairs formed by the leaf nodes of the lander's hierarchical bounding box and the leaf nodes of the irregular particles' hierarchical bounding box within the linear hierarchical bounding box, calculate the contact force and contact torque between the lander and the irregular particles.

[0099] S231, The lander and the irregular particles include several pairs of triangular facets, such as Figure 6 The diagram shows a pair of contacting triangular facets, where the mutual intrusion depth u of the k-th triangular facet pair is... nk for:

[0100]

[0101] Where, n i This represents the normal vector of the triangle corresponding to the lander; n j This represents the normal vector of the triangular face corresponding to the irregular particle; and These represent the locations of the lander's and the centroids of the irregular particle's triangular facet, respectively; T i and T j B represents the centroid of the triangular facet corresponding to the lander and the triangular facet corresponding to the irregular particle, respectively; i and B j These represent the centroids of the triangular facets corresponding to the lander and the triangular facets corresponding to the irregular particles, respectively.

[0102] S232. Calculate the normal contact force F of the k-th triangular element pair. cnk for:

[0103] F cnk =k nk u nk n k +C nk v nk (5)

[0104] Where, n k v represents the normal vector of the triangle corresponding to the k-th lander; nk This represents the relative velocity of the k-th triangular element pair; k nk C represents the stiffness coefficient. nk This represents the damping coefficient.

[0105] S233. Calculate the tangential contact force F of the k-th triangular element pair. csk for:

[0106]

[0107] Among them, u s This represents the relative displacement of the k-th triangular element with respect to the tangential direction; v s The relative velocity of the k-th triangular element with respect to the tangential direction is represented by μ; the coefficient of friction is represented by k. sk C represents the tangential elastic modulus; sk This represents the tangential damping coefficient.

[0108] S234, F based on each thread cnk and F csk To obtain the contact force and contact torque between the lander and irregular particles.

[0109] S3. Lander and Irregular Particle State Update: Based on the contact force and contact torque between the lander and the irregular particles, and combining the translational and rotational dynamic equations of the lander and the irregular particles, the state of the lander and the irregular particles is updated using second-order leapfrog integrals. The state includes translational motion state and rotational motion state, and the hierarchical bounding box of the lander and the irregular particles is updated according to the position information. The update of the translational motion state is as follows:

[0110]

[0111] Among them, v m v m+12 and v m+1 These represent the lander's velocity at the moment before the update, the first half of the update, and the moment after the update, respectively; a m and a m+1 These represent the lander's acceleration before and after the update, respectively; r m and r m+1 These represent the positions of the lander before and after the update, respectively; Δt represents the time interval.

[0112] The rotational motion state is updated as follows:

[0113]

[0114] Where, ω m and ω m+1 These represent the angular velocities of the lander before and after the update, respectively; q represents the updated angular acceleration of the lander; m and q m+1 Let f(q) represent the quaternions before and after the lander update, respectively; d represents the differential; t represents time; m Let ) denote a differential equation involving quaternions, and express it as:

[0115]

[0116] Where, (ω m ×) represents ω m An antisymmetric matrix.

[0117] The bounding box for the lander or irregular particles is updated as follows:

[0118]

[0119] Among them, V B Represents the original vertex coordinates of the lander or the bounding box of the irregular particle hierarchy; Q represents the quaternion; V N This represents the updated vertex coordinates of the lander or irregular particle hierarchy bounding box.

[0120] VB Represented as:

[0121]

[0122] in, and These represent the x, y, and z components of the original vertex coordinates of the bounding box.

[0123] Q is represented as:

[0124] Q = [q0 q1 q2 q3] T (9)

[0125] Where q0, q1, q2, and q3 represent the first to fourth components of the quaternion, respectively.

[0126] V N Represented as:

[0127]

[0128] in, and These represent the x, y, and z components of the vertex coordinates after the hierarchical bounding box is updated.

[0129] In one specific embodiment, a cubic lander is simulated moving through a layer of gravel particles on a small celestial body, as illustrated in the example. Figure 7 As shown in the figure, there are 19,023 irregular particles. Figure 7 The color of irregularly shaped particles indicates the velocity of the gravel particles; a darker color indicates a slower migration velocity, while a lighter color indicates a faster migration velocity. From... Figure 7 As can be seen from the initial moment of the simulation up to 0.25s, as the lander gradually penetrates downwards into the gravel layer, its energy is gradually transferred to the gravel layer, causing the velocity of the gravel layer to increase. After the irregular particles dissipate energy through mutual contact and friction, the velocity of the gravel layer gradually decreases. After the lander's velocity decays to 0, it begins to bounce upwards and leave the gravel layer.

[0130] This invention presents a dynamic numerical simulation method for the low-velocity impact of a lander onto a gravelly surface of a small celestial body. Addressing the shortcomings of the discrete element method assumption for spherical particles, it constructs a Delloy triangular mesh model, capable of simulating the contact between irregularly shaped particles and the lander, closely resembling the real surface environment of a gravelly particle layer on a small celestial body. Based on the Delloy triangular mesh, the contact areas between irregular particles are determined, and a distributed continuous contact force model is used to calculate the contact forces and moments between irregular particles. This method solves the problem of multi-point contact between irregular particles and the lander, and better reflects the physical characteristics of contact between irregular particles. The process is more reliable; by leveraging the heterogeneous computing framework of CPU and GPU, the contact detection, contact force calculation, and state update of irregular particles and the lander are completed in the GPU; in the most time-consuming contact detection stage, the contact detection is divided into coarse-grained contact detection, fine-grained contact detection, and relative intrusion depth calculation. The GPU thread hierarchy is reasonably organized, which can reduce the computational burden in each thread, classify the computational tasks into finer parallel granularity, improve the program's computational efficiency, and support the simulation of the contact and collision process between the lander and a large-scale irregular particle system.

[0131] Finally, it should be noted that the above embodiments are for illustration only and not for limiting the technical solutions of the present invention. Although the present invention has been described in detail with reference to the above embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the present invention without departing from the spirit and scope of the present invention. Any modifications or partial substitutions should be covered within the scope of the claims of the present invention.

Claims

1. A dynamic numerical simulation method for a lander's low-speed impact on a small celestial body's gravelly surface, characterized in that, It includes the following steps: S1. Initialize the state of the lander and irregular particles: Determine the state information of the lander and irregular particles, and construct the hierarchical bounding box of the lander and irregular particles. S11. Determine the status information of the lander and irregular particles: Read in the Delaunay triangular mesh model of the lander and irregular particles, and send it from the CPU memory to the GPU memory. S12. Based on the Deloitte triangular mesh model, construct the hierarchical bounding box between the lander and the irregular particles; S2. Contact detection between the lander and irregular particles: coarse-grained contact detection is performed using the root node of the linear hierarchical bounding box, fine-grained contact detection is performed using the linear hierarchical bounding box, and the relative intrusion depth between the lander and the irregular particles is calculated. S21, coarse-grained contact detection, quickly identifies all irregular particles that may have come into contact with the lander; S22. Fine-grained contact detection: Using the logic of the thread hierarchy organized by the GPU, determine the leaf nodes of all irregular particle-level bounding boxes that are in contact with the internal nodes of the lander-level bounding box. S23. Calculate the contact force and contact torque between the lander and the irregular particles: Based on the triangular facet pairs formed by the leaf nodes of the lander's hierarchical bounding box and the leaf nodes of the irregular particles' hierarchical bounding box within the linear hierarchical bounding box, calculate the contact force and contact torque between the lander and the irregular particles. S3. Lander and Irregular Particle State Update: Based on the contact force and contact torque between the lander and the irregular particles, combined with the translational and rotational dynamic equations of the lander and the irregular particles, and using the second-order leapfrog integral, the translational and rotational motion states of the lander and the irregular particles are updated, and the hierarchical bounding boxes of the lander and the irregular particles are updated according to the position information; the update of the translational motion state is as follows: Among them, v m v m+12 and v m+1 These represent the lander's velocity at the moment before the update, the first half of the update, and the moment after the update, respectively; a m and a m+1 These represent the lander's acceleration before and after the update, respectively; r m and r m+1 These represent the lander's positions before and after the update, respectively; Δt represents the time interval. The update of the rotational motion state is as follows: Where, ω m and ω m+1 These represent the angular velocities of the lander before and after the update, respectively; q represents the updated angular acceleration of the lander; m and q m+1 Let f(q) represent the quaternions before and after the lander update, respectively; d represents the differential; t represents time; m () represents a differential equation involving quaternions; The hierarchical bounding box of the lander or irregular particles is updated as follows: Among them, V B Represents the original vertex coordinates of the lander or the bounding box of the irregular particle hierarchy; Q represents the quaternion; V N This represents the updated vertex coordinates of the lander or irregular particle hierarchy bounding box.

2. The dynamic numerical simulation method for low-speed impact of a lander onto a gravelly surface of a small celestial body according to claim 1, characterized in that, Step S12 specifically includes the following steps: S121. Based on the Delaunay triangular mesh model of the lander and the irregular particles, construct bounding boxes of all levels respectively, and execute steps S122 and S123. S122. Obtain the Morton code corresponding to the centroid of the triangular face: convert the decimal number in the coordinates of the centroid into a binary number, and then interleave the binary numbers to obtain the Morton code corresponding to the centroid. S123. Based on the obtained Morton coding, construct a hierarchical bounding box: S1231. Using the binary search method, determine the first distinct binary bit of the Morton code; S1232. Divide the left and right subtrees: Nodes with the highest bit of 0 are assigned to the left subtree and become left nodes, and nodes with the highest bit of 1 are assigned to the right subtree and become right nodes. S1233. Based on the left and right subtrees respectively, repeat steps S1311 and S1312 to determine the next different binary bit of the Morton code and divide the left and right subtrees of the next level until the last binary bit of the Morton code is processed, and the hierarchical bounding box is completed.

3. The dynamic numerical simulation method for low-speed impact of a lander onto a gravelly surface of a small celestial body according to claim 1, characterized in that, Step S21 specifically includes the following steps: S211. Execute step S212 sequentially according to the numbering of the irregular particles. S212. Determine whether the maximum bounding box of the corresponding irregular particle overlaps with the maximum bounding box of the lander. If so, store the number of the corresponding irregular particle in the array and execute step S213; otherwise, execute step S213. S213. Repeat step S212 until the last irregular particle corresponding to the last number is executed; Step S22 specifically includes the following steps: S221. If the left and right subtrees of the lander's hierarchical bounding box are both leaf nodes and both intersect with the leaf nodes of the irregular particle's hierarchical bounding box, then determine whether the leaf nodes of the lander's hierarchical bounding box and the leaf nodes of the irregular particle's hierarchical bounding box intersect in the triangular faces within the linear hierarchical bounding box. S222. If the left and right subtrees of the lander's hierarchical bounding box are both internal nodes, and only one of the left or right subtrees intersects with the leaf node of the irregular particle's hierarchical bounding box, or if one of the left and right subtrees of the lander's hierarchical bounding box is a leaf node and the other is an internal node, and the internal node intersects with the leaf node of the irregular particle's hierarchical bounding box, then recursively determine whether all internal child nodes of the next level of the internal node of the lander's hierarchical bounding box intersect with the leaf node of the irregular particle's hierarchical bounding box. S223. If the left and right subtrees of the lander's hierarchical bounding box are both internal nodes and both intersect with the leaf nodes of the irregular particle's hierarchical bounding box, then the right subtree is temporarily stored, and the two child nodes of the next level of the left subtree are recursively judged to see if they intersect with the leaf nodes of the irregular particle's hierarchical bounding box. This process continues until all child nodes of all levels of the left subtree have been judged. Then, the temporarily stored right subtree is read, and the recursive judgment continues to see if the two child nodes of the next level of the right subtree intersect with the leaf nodes of the irregular particle's hierarchical bounding box. Step S23 specifically includes the following steps: S231. The lander and the irregular particles are separated by a number of triangular facet pairs, wherein the mutual intrusion depth u of the k-th triangular facet pair is... nk for: Where, n i This represents the normal vector of the triangle corresponding to the lander; n j T represents the normal vector of the triangular face corresponding to the irregular particle; i and T j These represent the centroids of the triangular facets corresponding to the lander and the triangular facets corresponding to the irregular particles, respectively. and These represent the locations of the centroids of the lander and the irregular particle triangular facet, respectively. S232. Calculate the normal contact force F of the k-th pair of triangular facets. cnk for: F cnk =k nk u nk n k +C nk v nk (5) Where, n k v represents the normal vector of the triangle corresponding to the k-th lander; nk The relative velocity of the k-th pair of triangular facets is represented by k; nk C represents the stiffness coefficient. nk Indicates the damping coefficient; S233. Calculate the tangential contact force F of the k-th pair of triangular facets. csk for: Among them, u s This represents the relative displacement of the k-th triangular element with respect to the tangential direction; v s The relative velocity of the k-th triangular element with respect to the tangential direction is represented by μ; μ represents the coefficient of friction; k sk C represents the tangential elastic modulus; sk Indicates the tangential damping coefficient; S234, F based on each thread cnk and F csk To obtain the contact force and contact torque between the lander and irregular particles.

4. The dynamic numerical simulation method for low-speed impact of a lander onto a gravelly surface of a small celestial body according to claim 1, characterized in that, In step S3, f(q) m ) is represented as: Where, (ω m ×) represents ω m The antisymmetric matrix; V B Represented as: in, and These represent the x, y, and z components of the original vertex coordinates of the bounding box. Q is represented as: Q=[q0 q1 q2 q3] T (9) Where q0, q1, q2, and q3 represent the first to fourth components of the quaternion, respectively; V N Represented as: in, and These represent the x, y, and z components of the vertex coordinates after the hierarchical bounding box is updated.

5. The dynamic numerical simulation method for low-speed impact of a lander onto a gravelly surface of a small celestial body according to claim 1, characterized in that, In step S12, the first layer of the hierarchical bounding box has the largest linear hierarchical bounding box; the second layer of the hierarchical bounding box divides the Delaunay triangular mesh in two, evenly distributing it into the linear hierarchical bounding boxes of the left and right subtrees, and so on downwards until the last layer of the hierarchical bounding box contains only one triangular face; the linear hierarchical bounding box of the first layer of the hierarchical bounding box is called the root node, the linear hierarchical bounding boxes of the second to penultimate layers of the hierarchical bounding box are called internal nodes, and the linear hierarchical bounding box of the last layer of the hierarchical bounding box is called the leaf node; all Delaunay triangular meshes are contained in the root node, and all nodes except the leaf nodes are collectively called internal nodes.

6. The dynamic numerical simulation method for low-speed impact of a lander onto a gravelly surface of a small celestial body according to claim 1, characterized in that, The logic of organizing the GPU thread hierarchy in step S22 is that the number of thread blocks of the GPU is equal to the number of irregular particles that may have contact with the lander, and the number of threads allocated in each thread block is the same as the number of leaf nodes of the hierarchical bounding box of the irregular particle; that is, one thread block is responsible for handling the contact detection between one irregular particle and the lander, and one thread is allocated a leaf node of the hierarchical bounding box of the irregular particle.

7. The dynamic numerical simulation method for low-speed impact of a lander on a gravelly surface of a small celestial body according to claim 1, characterized in that, In step S11, after pre-allocating video memory to the GPU using the memory allocation function cudaMalloc, the status information of the lander and irregular particles on the CPU is sent to the GPU using the memory allocation function cudaMalloc.

8. The dynamic numerical simulation method for low-speed impact of a lander onto a gravelly surface of a small celestial body according to claim 1, characterized in that, The state information mentioned in step S11 includes position, velocity, angular velocity, and quaternion.

Images