A method of simulating a bullet impact

Abstract

The application relates to the field of impact dynamics and protective engineering, and discloses a projectile-target impact simulation method. The method comprises the following steps: generating initialization finite element model data of a projectile body and mixed discrete model data of a target plate respectively, wherein the target plate model is composed of a core SPH particle region and a peripheral finite element grid region through a fixed connection relationship; subsequently, taking the two as initial system state data, starting a time step iteration solving cycle to perform dynamic simulation; when a simulation termination condition is met, outputting a final target plate damage form, a projectile penetration result and a ballistic limit speed. Through self-adaptive mixed discrete and coupling technology, high-fidelity and high-efficiency numerical simulation of the whole impact process is realized. The technical problem that the existing single numerical method is difficult to balance the calculation accuracy and efficiency when simulating low-strength projectile impact on high-strength target plates is solved.

Description

Technical Field

[0001] This invention relates to the field of impact dynamics and protective engineering, and in particular to a method for simulating projectile impact. Background Technology

[0002] Numerical simulation studies of the physical processes of bullets and other penetrating agents impacting metal targets are of great significance for evaluating the performance of protective structures and optimizing their design. Among these, the impact of "low-strength projectiles on high-strength targets" is a typical engineering problem, often encountered when evaluating the protective capability of solid rocket motor casings against small arms fragments. In such problems, the projectile material is relatively soft, undergoing significant plastic upsetting, mushroom-shaped deformation, and even fragmentation during impact; while the target material is high-strength, potentially exhibiting various complex failure modes including bulging, petal-shaped cracking, and slugging. This mismatch between material strength and deformation mechanisms presents a significant challenge to achieving high-fidelity numerical simulations of this dynamic process.

[0003] Currently, numerical simulation research in this field mainly relies on single methods such as mesh-based finite element method (FEM) or particle-based smoothed particle hydrodynamics (SPD). While the FEM has advantages in computational efficiency, severe mesh distortion occurs when simulating extreme large deformations and fractures such as projectile bulging, target material tearing, and fragment separation, leading to decreased computational accuracy, premature termination, or complete failure. In contrast, meshless methods such as SPD can naturally handle large deformations and material fractures, but their computational cost is high. When simulating large-sized targets, the number of particles required to ensure accuracy is enormous, making it difficult to complete the calculation within an engineering-acceptable timeframe. Although some existing research attempts to couple different methods, designing an efficient and stable coupling mechanism to achieve high-precision simulation of the entire process from impact, penetration, to failure remains a pressing technical challenge when dealing with complex impact problems such as "low-strength projectile-high-strength target" scenarios involving large deformations in both materials and requiring precise capture of multiple failure modes. Summary of the Invention

[0004] Based on this, it is necessary to propose a projectile impact simulation method to address the technical problem that the single numerical method used in the existing technology cannot simultaneously take into account simulation accuracy and computational efficiency.

[0005] Firstly, a method for simulating the impact of a projectile target is provided, the method comprising: Read the input projectile geometry model data, the material parameter data of the projectile and the target plate, and the impact initial condition data including the initial velocity of the projectile; Based on the projectile geometric model data, generate three-dimensional finite element mesh model data of the projectile; Based on the projectile material parameter data, the corresponding material property data are assigned to the elements in the three-dimensional finite element mesh model data; Based on the initial impact condition data, the mesh nodes of the three-dimensional finite element mesh model data are assigned initial motion state data, thereby generating initialized projectile finite element model data. Based on the input target plate geometric model data and material parameter data, a hybrid discrete model data of the target plate is generated; wherein, the hybrid discrete model data consists of SPH particle model data of the core region and finite element mesh model data of the outer region, and a consolidation connection relationship data is established at the interface between the two. Using the initialized finite element model data of the projectile and the hybrid discrete model data of the target plate as the initial system state data, the time-step iterative solution loop is started; When the simulation termination condition is met, the time-step iterative solution loop is exited, and target plate damage morphology data, projectile penetration result data, and ballistic limit velocity are generated and output based on the final system state data.

[0006] The beneficial effects of this application are: The target plate is discretized into zones. SPH particles are used for discretization in the core region where large deformation and fracture are expected, while finite element meshes are used for discretization in the peripheral region where deformation is smaller. Consolidation is established, which significantly reduces the computational scale while ensuring high-fidelity capture of complex physical phenomena such as penetration and failure, thus overcoming the problem of excessive computational cost of the pure SPH method. The projectile is discretized using finite element methods and is allowed to be dynamically converted into SPH particles according to the degree of distortion during the calculation. This adaptively combines the high efficiency of the finite element method in the initial stage of calculation and small deformation with the robustness of the SPH method in handling extreme large deformation, overcoming the problem of computational interruption or distortion caused by mesh distortion in the pure finite element method. Finally, this integrated scheme achieves continuous, efficient and high-precision numerical reproduction of the entire stage of impact, penetration and failure / penetration in the impact process of low-strength projectile to high-strength target. Attached Figure Description

[0007] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0008] in: Figure 1 This is a diagram illustrating the application environment of a computer device in one embodiment; Figure 2 This is a flowchart of a projectile-target impact simulation method in one embodiment; Figure 3 This is a structural block diagram of a projectile impact simulation device in one embodiment; Figure 4 This is a structural block diagram of a computer device in one embodiment. Detailed Implementation

[0009] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0010] Figure 1 An exemplary network architecture is shown that can be applied to the logistics waybill management method or logistics waybill management device of this application.

[0011] like Figure 1 As shown, the network architecture may include computer device 101 and server 102. Computer device 101 and server 102 can communicate with each other via the network, which serves as the medium for providing communication links between the various units. The network may include various types of wired or wireless communication links, such as wired communication links including fiber optic cables, twisted-pair cables, or coaxial cables, and wireless communication links including Bluetooth communication links, Wi-Fi communication links, or microwave communication links.

[0012] It should be noted that computer device 101 and server 102 can be either hardware or software. When computer device 101 and server 102 are hardware, they can be implemented as a distributed server cluster consisting of multiple servers, or as a single server. When computer device 101 and server 102 are software, they can be implemented as multiple software programs or software modules (for example, used to provide distributed services), or as a single software program or software module; no specific limitations are made here.

[0013] The computer device described in this application can be equipped with various communication client applications, such as video recording applications, video playback applications, voice interaction applications, search applications, instant messaging tools, email clients, social platform software, etc.

[0014] Computer devices can be either hardware or software. When a computer device is hardware, it can be any computer device with a display screen, including but not limited to smartphones, tablets, laptops, and desktop computers. When a computer device is software, it can be installed on the computer devices listed above. It can be implemented as multiple software programs or software modules (e.g., used to provide distributed services) or as a single software program or software module; no specific limitation is made here.

[0015] When a computer device is used as hardware, it can also be equipped with a display device and a camera. The display device can be any device capable of displaying information, and the camera is used to capture video streams. Examples of display devices include cathode ray tube displays (CR), light-emitting diode displays (LED), e-ink screens, liquid crystal displays (LCD), and plasma display panels (PDP). Users can utilize the display devices on their computer devices to view displayed text, images, videos, and other information.

[0016] It should be understood that Figure 1 The number of computer devices, networks, and servers shown is for illustrative purposes only. The number of computer devices, networks, and servers can be any number, depending on the implementation requirements.

[0017] Please see Figure 2 As shown, Figure 2 A flowchart illustrating the projectile impact simulation method provided in this embodiment of the invention includes the following steps: S1. Read the input projectile geometry model data, the material parameter data of the projectile and the target plate, and the impact initial condition data including the initial velocity of the projectile.

[0018] It receives three sets of structured data provided by users or external systems through its input interface.

[0019] Projectile geometric model data refers to the digital information used to uniquely determine and reconstruct the three-dimensional solid shape and internal structure of a projectile, including all dimensions, topology, and assembly relationships. This data defines the spatial morphology of the projectile and forms the geometric basis for finite element mesh discretization.

[0020] For example, when simulating a 7.62mm rifle cartridge, the projectile geometry model data should include: a projectile diameter of 7.62 mm, a projectile length of 10.75 mm, and a projectile radius of curvature of 40.00 mm. It can also include layered dimensions such as a steel core diameter of 5.80 mm, a lead sheath thickness (which can be calculated), a jacket thickness of 0.5 mm, and a steel core semi-cone angle of 11.9 degrees. This data can generate a correctly assembled three-dimensional solid comprising the jacket, lead sheath, and steel core.

[0021] Material parameter data refers to the set of all physical constants and model parameters used in numerical simulations to quantitatively describe the mechanical behavior and constitutive response of projectile and target materials. These parameters determine the stress-strain relationship, damage accumulation, and failure process of the material under impact loads.

[0022] For example, for a 30CrMnSiA steel target plate, material parameters should include: density 7850 kg / m³, elastic modulus 200 GPa, yield stress A 525 MPa, hardening coefficient B 101 MPa, strain rate C 0.1739, damage constants D1=0.0705, D2=1.732, etc. For a lead sleeve for a projectile, elastic-plastic parameters such as density 11340 kg / m³ and yield strength 20 MPa should be provided.

[0023] Impact initial condition data refers to the set of all parameters that define the initial kinematic state and relative spatial position of the projectile impact event. This data sets the starting point for the entire dynamic simulation process.

[0024] The initial velocity of the projectile refers to the instantaneous velocity of the projectile just before it makes contact with the target plate. It is a vector, and its magnitude (scalar value, such as meters per second) and direction (such as the components of the velocity vector in the global coordinate system or defined by the impact angle) usually need to be specified.

[0025] For example, the initial impact conditions for a simulation can be defined as follows: the initial velocity of the projectile is 616 meters per second, the impact angle is 0 degrees (indicating that the velocity direction is perpendicular to the target surface, i.e., a normal impact), the projectile axis is parallel to the target normal, and the initial spacing is 0.8 millimeters. Based on this, an initial velocity vector of 616 meters per second, perpendicular to the target, will be assigned to all projectile nodes, and the projectile model will be placed 0.8 millimeters from the target surface.

[0026] S2. Based on the projectile's geometric model data, generate the projectile's three-dimensional finite element mesh model data.

[0027] After analyzing the projectile's geometric model data, the geometric modeling and meshing process is initiated to generate a complete set of data describing the projectile's discrete spatial morphology, namely, three-dimensional finite element mesh model data. This process transforms a continuous geometric shape into a set composed of a finite number of discrete elements and their nodes, laying the topological foundation for subsequent physical property assignment and mechanical calculations. Its core lies in accurately reconstructing the projectile's three-dimensional solid model in computational memory based on the specific dimensional parameters obtained from the analysis, using a built-in geometric engine. This reconstruction process is not simply drawing the outline, but follows the projectile's composite structure defined by the input data, distinguishing and constructing the geometric regions corresponding to different material layers such as the jacket, lead sheath, and steel core, and assigning these regions independent logical identifiers. After the solid model is established, its mesh generator is invoked to perform finite element meshing of the projectile solid, primarily using hexahedrons. During the meshing process, nodes are arranged in the solid space according to preset mesh accuracy requirements (such as global element size or curvature adaptability), and adjacent nodes are connected to form elements, while automatically recording the material layer region to which each element belongs. Finally, a structured dataset is constructed and output. The core content of this dataset includes the three-dimensional coordinates of all nodes, the connection relationships of all elements consisting of which nodes, and the correspondence between each element and a specific material layer identifier. These data together constitute the three-dimensional finite element mesh model data of the projectile, which is a "skeleton" model that only includes geometric and topological relationships and has not yet been assigned physical properties.

[0028] For example, based on the analyzed projectile geometry model data, which specifies that the projectile consists of three layers: a jacket, a lead sheath, and a steel core, and provides the precise dimensions of each layer, three coaxial geometric entities are first created in memory, representing a 0.5 mm thick copper-clad steel jacket cylinder, the inner lead sheath layer, and the innermost low-carbon steel core. Next, a meshing algorithm discretizes these three combined entities, potentially generating a total of 5872 hexahedral elements. During meshing, elements numbered 3050 to 3120 are located in the steel core region, elements numbered 1000 to 1500 are located in the lead sheath region, and the remaining elements are located in the jacket region. These element sets are labeled "Core," "Lead," and "Jacket," respectively. The final output 3D finite element mesh model data is a digital model including thousands of node coordinates, tens of thousands of element connections, and material layer labels for each element. For example, the coordinates of data record node 1234 are (1.2, 0.5, 10.1) mm, unit 5678 is composed of nodes [1234, 1235, 1240, 1239, ...], and the material layer associated with this unit is identified as "Jacket".

[0029] S3. Based on the projectile material parameter data, assign corresponding material property data to the elements in the three-dimensional finite element mesh model data.

[0030] After obtaining the three-dimensional finite element mesh model data of the projectile, the next step is to assign material properties. The goal of this step is to accurately associate the mechanical properties of different material layers described in the projectile material parameter data with each discrete element in the mesh model, thereby transforming the purely geometric mesh into a computational model with physical properties.

[0031] The execution of this process relies on the association data between the elements and material layer regions established in the previous step. First, the projectile material parameter data is read and parsed. This data is typically organized by material layer name, with each set of parameters defining an elastoplastic material model and its specific constants. Then, each element in the mesh model is traversed, and based on the material layer region identifier carried by that element, the corresponding set of material parameters is searched and matched in the material parameter database. Upon successful matching, a material property data object is instantiated for that element. This object encapsulates the specified material model type and the complete set of corresponding material constants. For example, for an element identified as a steel core, its material property data object will be filled with the elastic modulus, Poisson's ratio, density, and yield strength of the steel core. Material property data refers to the complete set of parameters used in numerical simulation to quantitatively define and describe the specific mechanical response laws followed by each material constituting the projectile and target plate under external loads (especially high-speed impact loads). These include: Elastic modulus: describing the linear proportional relationship between stress and strain during the elastic deformation stage of the material; Poisson's ratio: describing the ratio of lateral strain to axial strain when the material is under uniaxial tension or compression. Density: The mass per unit volume of the material, used to calculate inertial forces. Initial yield strength: The critical stress value at which the material begins to undergo irreversible plastic deformation. Finally, the geometric connection information of all elements and their corresponding material property data objects are integrated to form a more complete finite element mesh model of the projectile with enhanced physical properties, preparing for the next stage of dynamic state initialization.

[0032] For example, a pre-defined projectile mesh model already exists in memory, where cell 5678 is marked as belonging to the jacket layer and cell 3050 is marked as belonging to the steel core layer. The system simultaneously reads material parameter data, defining a material named "Jacket" with an elastic modulus of 200 GPa and a yield strength of 0.3 GPa, and a material named "Core" with an elastic modulus of 200 GPa and a yield strength of 0.3 GPa, but with potentially different densities. When traversing to cell 5678, the "Jacket" material parameter set is retrieved based on its "Jacket" identifier, and material property data is created and assigned to that cell. Similarly, when traversing to cell 3050, the material properties of the steel core are assigned based on its "Core" identifier. After traversing all cells, the original geometric mesh data is upgraded to a comprehensive model including geometric and material properties, where each cell explicitly knows which set of physical laws and constants its material behavior will follow. S4. Based on the initial impact condition data, the initial motion state data is assigned to the mesh nodes of the three-dimensional finite element mesh model data, thereby generating the initialized projectile finite element model data.

[0033] After completing the geometric construction and material property assignment of the projectile's finite element mesh, the motion state initialization step is executed. The purpose of this step is to set an initial state for the projectile mesh model, which already possesses physical properties, based on the initial impact conditions data, conforming to the actual physical scenario and allowing for the commencement of dynamic calculations. This process transforms the macroscopic impact conditions into specific kinematic values ​​for each node on the discrete model.

[0034] First, the initial impact condition data is analyzed to obtain key information about the projectile's initial motion vector, including velocity magnitude, direction, and initial positional relationship within the global computational space. Based on this information, an initial velocity field applicable to the entire projectile is calculated. For each finite element mesh node on the projectile, its position relative to the projectile's center of mass or reference point, combined with the impact angle, is used to determine and assign a corresponding initial velocity vector. This vector ensures that the projectile, as a rigid body, exhibits overall translational and rotational motion conforming to the input conditions at the initial moment. Simultaneously, the internal state variables of all finite elements are initialized, including but not limited to setting the stress tensor of each element to zero, the cumulative plastic strain to zero, the damage variable to zero, and the temperature to the ambient reference temperature. After assigning values ​​to all nodes and elements, the updated node motion data, element state data, and previously existing mesh geometry and element material property data are integrated and encapsulated to form a complete finite element model of the projectile, including all necessary initial information. This data model is ready and can be imported into a dynamics solver as the starting point for calculations.

[0035] For example, the initial impact conditions data specify an initial projectile velocity of 616 meters per second, an impact angle of 0 degrees (normal impact), and an initial position of the projectile tip 0.8 millimeters from the target surface. After analyzing this data, the initial velocity vector direction is determined to be perpendicular to the target surface. Subsequently, every node of the projectile's finite element mesh is traversed, assigning all nodes the same initial velocity vector with a magnitude of 616 meters per second and a direction perpendicular to the target. Simultaneously, the entire projectile mesh model is translated in space so that its tip is precisely positioned 0.8 millimeters from the target surface. For each element in the mesh, regardless of whether it belongs to the jacket, lead sleeve, or steel core, the program resets its internal equivalent stress, plastic strain, and damage values ​​to 0, and sets the temperature value to 293 Kelvin (room temperature). The final generated projectile finite element model data is a digital model where all nodes have a defined initial velocity, all elements are in a "natural" state without deformation or damage, and it is in the correct initial relative position to the target.

[0036] S5. Based on the input target plate geometric model data and material parameter data, generate hybrid discrete model data of the target plate; wherein, the hybrid discrete model data consists of SPH particle model data of the core region and finite element mesh model data of the outer region, and a consolidation connection relationship data is established at the interface between the two.

[0037] The process involves constructing a hybrid discrete model of the target plate based on the input geometric model data and material parameter data. The core innovation of this step lies in addressing the highly localized deformation of the target plate in impact simulations by decomposing a single geometric target plate into two coupled computational regions discretized using different numerical methods, thus balancing computational accuracy and efficiency.

[0038] First, the target plate's geometric model data is analyzed to obtain its overall size and shape. Combined with the expected impact point location, the geometric model is logically divided into two regions. Typically, a rectangular or circular region centered on the impact point projection is defined as the core region, and the portion of the target plate outside this region is defined as the outer region. After division, different discretization strategies are applied to these two regions. For the core region, the program generates a set of smoothed particle hydrodynamics (SPH) particles in a regular arrangement within three-dimensional space based on a preset initial particle spacing. This forms SPH particle model data, which includes the initial spatial coordinates of each particle and attributes such as mass calculated based on spacing and material density. For the outer region, the program uses the finite element method to mesh the target plate, generating finite element mesh model data. This process specifically requires that the size of the finite element mesh at the interface adjacent to the core region matches the spacing of the SPH particles to ensure coupling compatibility. Simultaneously, a mesh transition technique is employed to gradually increase the element size from the interface towards the outer edge of the target plate, effectively controlling the total mesh size.

[0039] Subsequently, a crucial coupling operation is performed: establishing a consolidation connection at the interface between the SPH particles and adjacent finite element meshes. This is achieved through a specific consolidation algorithm that creates a constraint relationship or mapping weight between the particles and neighboring finite element nodes, ensuring that forces, displacements, and energy can be continuously transferred across the interface in subsequent calculations. Then, the unified material constitutive model and its complete set of parameters defined in the target material parameter data are assigned to all SPH particles and all finite element elements in the surrounding region, and their internal state variables such as stress, strain, and damage are initialized. Finally, fixed constraints are applied to the nodes of the surrounding finite element mesh located at the physical boundary of the target plate, and the aforementioned SPH particle data, finite element mesh data, interface consolidation relationship data, material property data, and boundary condition data are integrated and encapsulated, outputting a complete and computable target plate hybrid discrete model data.

[0040] For example, reading the target plate's geometric model data reveals that the target plate is a square plate with a side length of 250 mm and a thickness of 4.1 mm. The program automatically defines a square region with a side length of 27 mm as the core region, using the center of the plate as the impact point. Within this core region, a layer of 27,225 SPH particles are generated at a particle spacing of 0.5 mm. The mass of each particle is calculated based on the density of 30CrMnSiA steel. Around the core region, the program generates a finite element mesh. The element size is also 0.5 mm immediately adjacent to the core region boundary, gradually increasing to several millimeters towards the outer edge of the target plate. At the contact surface between the core region particles and the outer mesh, a consolidation algorithm calculates the interaction weight between each boundary particle and the surrounding finite element nodes and stores it as connectivity data. Subsequently, all particles and elements were assigned the same set of modified Johnson-Cook model parameters for 30CrMnSiA steel. The Johnson-Cook model is a phenomenological constitutive model based on the Mises yield criterion and employing a multiplicative combination, describing the strain hardening, strain rate strengthening, and thermal softening effects of the material through three sets of separable multipliers. Finally, fixed constraints were applied to the finite element nodes on the four sides of the target plate, and all data were packaged into a final hybrid discrete model of the target plate. Material parameter data can be obtained through material experiments (such as Hopkinson bar tests) or by referencing publicly available material databases.

[0041] S6. Using the initialized finite element model data of the projectile and the mixed discrete model data of the target plate as the initial system state data, start the time-step iterative solution loop.

[0042] After initializing the projectile and target models, the integrated initial system state data is loaded into the explicit dynamics solver, formally initiating a time-stepped core calculation process, i.e., a time-stepped iterative solution loop. This loop is a strictly controlled, repeatedly executed logical sequence used to numerically reproduce the dynamic evolution of the impact event over time.

[0043] The loop begins with the current time step as zero or an initial value. At the start of each time step, it first checks whether the current system state meets the preset simulation termination condition. If not, it sequentially executes a series of physical calculation sub-steps defined within the loop. These sub-steps include contact interaction detection and force calculation, material constitutive response and state update, adaptive element type conversion based on deformation degree, failure handling, and integral progression of the system motion state. Each sub-step receives updated data from the previous sub-step and passes its calculation results to the next sub-step, forming a coherent data processing chain. After all sub-steps in a time step have been executed, the system time is advanced by a small time increment, and based on the updated state data of all nodal displacements, velocities, stresses, and damage, the next round of loop determination and calculation begins. This process repeats cyclically, dynamically simulating the entire process of projectile penetration of the target plate, target plate deformation, and failure, until the system state is detected to meet the simulation termination condition, at which point the loop exits.

[0044] Specifically, the time-step iterative solution loop of this application includes the following sub-steps: Sub-step S61: Based on the current spatial position data of the projectile and the target plate SPH particles, perform contact search and generate contact pair data, and calculate and output contact force data based on the contact pair data.

[0045] Within the time-step iterative solution loop, a contact search and contact force calculation sub-step is executed. The purpose of this step is to accurately detect whether physical contact has occurred between the projectile and the target plate at the current moment, and to quantitatively calculate the interaction force for all occurring contacts, providing boundary load input for subsequent material deformation and motion state updates.

[0046] First, the real-time spatial coordinates of all projectile surface finite element nodes, converted projectile SPH particles, and target plate SPH particles in the system are acquired at the current time step. Based on this position data, a contact search algorithm is invoked. This algorithm typically consists of two stages. The first stage is a global search, which quickly filters out all potential object pairs that may come into contact through spatial partitioning or bounding box techniques, such as a projectile node located in the vicinity of a target plate particle. The second stage is a local precise detection, which accurately determines the geometric relationship of each potential contact object pair. If the projectile object is the surface of a finite element element, the minimum distance and penetration depth from the target plate particle to the surface of that element are calculated, and the normal direction of the contact point is determined. If the projectile object is also an SPH particle, the distance between the centers of the two particles is directly calculated. Based on the results of the precise detection, a structured list of contact pair data is generated, where each record includes the identifiers of the two objects that are in contact, the contact type, and contact geometry information such as penetration depth and contact normal.

[0047] Subsequently, based on the contact pair data and the material property data of the associated objects, the contact force is calculated. For the contact between the surface of the projectile finite element element and the target plate SPH particles, the penalty function method is typically used for calculation. This method treats the penetration depth as a virtual linear spring compression, calculates the magnitude of the normal contact force based on the equivalent material stiffness of the contact point, and then calculates the tangential friction force based on the Coulomb friction model and tangential relative velocity. For the contact between the projectile SPH particles and the target plate SPH particles, the interaction force is directly calculated based on the particle approximation theory of the SPH method, solved using the constitutive stress, density, and kernel function gradient of the two particles. After the calculation is completed, the corresponding contact force data is generated, and according to the principle of action and reaction forces, the calculated contact force vector is distributed to both objects in the contact pair; that is, the force is applied to the target plate particles, and the equal and opposite reaction forces are applied to the projectile nodes or particles.

[0048] For example, at a certain time step, the projectile's head has been partially converted into SPH particles. After reading the position data, a global search reveals that a specific SPH particle in the projectile's head region is very close to dozens of SPH particles on the target plate surface, and is listed as a potential contact pair. Local precise detection confirms that the center-to-center distance between this projectile particle and three of the target plate particles is less than the sum of their smooth lengths, thus determining that contact has occurred. These three particle pairs and their distance information are recorded as SPH-SPH (Smoothed Particle Hydrodynamics - Smoothed Particle Hydrodynamics contact pair) type contact pair data. Simultaneously, the search also reveals that an unconverted finite element surface at the projectile's tail penetrates a target plate particle, recorded as an FEM-SPH (Finite Element Method - Smoothed Particle Hydrodynamics) type contact pair, and the penetration depth is calculated to be 0.02 mm.

[0049] For the contact between SPH particle pairs, the SPH interaction force model is invoked, and the repulsive force between them is calculated based on the current stress, density, and kernel function derivative of the two particles. For the contact between the finite element element and the particles, a normal contact force of 150 Newtons is calculated using the penalty function method based on the equivalent material stiffness at that point (combining the modulus of the projectile jacket and the target plate steel) and a penetration depth of 0.02 mm. Finally, the calculated repulsive forces are applied to the corresponding pairs of SPH particles; a 150-Newton force is applied to the target plate particles, while a 150-Newton force of the same magnitude but in the opposite direction is distributed to the corresponding nodes constituting the finite element element. These force data will be aggregated for the next step of updating the acceleration of all objects.

[0050] Sub-step S62: Input the contact force data and the current stress, strain, and damage variable data of the elements and particles into the material constitutive model for calculation, and output the updated stress data, strain data, damage variable data, and temperature data.

[0051] Within the time-step iterative solution loop, a sub-step for material constitutive response and state update is executed. The purpose of this step is to calculate the new mechanical and thermal state of the material at this moment based on the currently applied external contact force and the material's own historical state, according to physical laws, thereby describing its deformation, hardening, damage, and even heating processes.

[0052] First, the contact force data acting on each unit or particle is combined with the current internal stress data of the object to calculate the strain increment and equivalent plastic strain rate data for the current time step using constitutive relations. Then, a sequential physical quantity update process based on a specific material model is initiated. The cumulative plastic strain and temperature data from the previous time step of the unit or particle are read, and new cumulative plastic strain data are obtained based on the newly calculated plastic strain increment. Simultaneously, according to the adiabatic assumption, the plastic deformation work generated in this time step is converted into heat in a certain proportion, and the temperature rise is calculated and updated to new temperature data. Next, the damage evolution model is invoked. The damage evolution model takes the new cumulative plastic strain data, stress triaxiality data characterizing the stress state, equivalent plastic strain rate data, and new temperature data as input. Based on a predefined damage evolution law, it determines whether damage has begun to accumulate and calculates the damage increment, finally outputting the updated damage variable data.

[0053] Then, the updated damage variable data, new cumulative plastic strain data, equivalent plastic strain rate data, and new temperature data are input into the Johnson-Cook model for coupled damage. This model comprehensively calculates multiple factors such as strain hardening, strain rate effect, temperature softening, and damage weakening to solve for the instantaneous yield stress data of the material under the current deformation and damage state. Finally, using this yield stress data as the target, combined with the strain increment data, the stress tensor of the previous time step is updated using an elastoplastic stress integration algorithm (such as the radial back mapping algorithm) to obtain new stress tensor data that satisfies the yield condition and flow law, thus completing the update of all state variables at this time step.

[0054] For example, for a specific SPH particle on a target plate, its cumulative plastic strain in the previous step was read as 0.15, the damage variable as 0.1, and the temperature as 320 Kelvin. In the current time step, the calculated equivalent plastic strain increment is 0.005. First, its cumulative plastic strain is updated to 0.155. Assuming a plastic work-thermal conversion factor of 0.9, the calculated temperature rise is approximately 2 Kelvin, hence the new temperature is 322 Kelvin. Next, assuming its current stress triaxiality is 0.5 and the strain rate is 1000 per second, according to the damage evolution formula, since the cumulative plastic strain has exceeded the damage threshold, a damage increment of 0.002 is generated, hence the new damage variable is 0.102. Then, substituting the new plastic strain of 0.155, the new damage of 0.102, the strain rate of 1000 per second, and the new temperature of 322 Kelvin into the Johnson-Cook formula for coupled damage, the calculated yield stress after considering damage softening is 800 MPa. Finally, using this constraint, the stress tensor from the previous step is updated from 780 MPa to 790 MPa, which matches the current strain increment and satisfies the yield surface condition, through a stress integration algorithm. The output results are the updated stress, strain, damage, and temperature data for this particle.

[0055] Sub-step S63: Calculate the distortion measurement data of the elements in the finite element model of the projectile based on the updated strain data of the elements. When the distortion measurement data exceeds the preset threshold, generate element conversion instruction data.

[0056] Within the time-step iterative solution loop, a sub-step for evaluating and converting element distortion metrics is executed. The purpose of this step is to monitor the deformation state of each element in the projectile finite element model, determine through quantitative indicators whether it is no longer suitable for further calculation using the finite element method due to excessive distortion, and make an intelligent decision on whether to convert it into an SPH particle. This is a key control step in achieving adaptive numerical method coupling.

[0057] First, updated strain data for all finite element elements of the projectile are extracted from the output of the material constitutive response calculation, particularly the equivalent plastic strain data for each element. Simultaneously, based on the node coordinates of these elements at the current time step, geometric parameters reflecting the degree of element shape distortion can be calculated, such as the rate of change of the element's Jacobian determinant relative to its initial value or the element's aspect ratio. Using these physical strain and geometric distortion data as input, a comprehensive, scalarized distortion metric is calculated for each element through a predefined evaluation function. This metric quantifies the degree to which the element deviates from its initial regularity.

[0058] Subsequently, the distortion metric data of each element is compared one by one with a preset conversion threshold. This threshold is a critical value determined in advance through calibration based on numerical stability and accuracy requirements. If the distortion metric data of an element exceeds the preset threshold, the element is determined to have suffered severe distortion. For all elements determined to be severely distorted, a structured element conversion instruction is generated. This instruction data includes at least a list of unique identifiers of the elements to be converted, as well as the current time step information, and its purpose is to indicate which elements need to undergo type conversion operations from finite element to SPH particle in subsequent steps.

[0059] The determination of this preset threshold is a process of calibration and optimization based on the mechanical properties of the projectile material under impact, the type and size of the finite element elements used, and the stability requirements of numerical calculations. This is achieved through pre-calibrated numerical experiments combined with engineering experience. Specifically, firstly, based on the constitutive model of the projectile material, theoretical analysis or single-element numerical experiments are used to preliminarily estimate the theoretical reference value of the distortion metric when the element experiences severe distortion (such as the Jacobian determinant approaching zero or the equivalent plastic strain approaching the failure strain). Subsequently, a simplified two-dimensional axisymmetric or three-dimensional local model including typical impact velocities is established. In such verification simulations, the preset threshold is adjusted dynamically, and the projectile under different thresholds is compared and analyzed. The simulation results of the large deformation region of the projectile's head include the degree of element distortion, mass-momentum conservation error, and computational stability, to observe the impact of the threshold on simulation accuracy and efficiency. Finally, a threshold value is selected that can ensure timely triggering of the transformation before severe element distortion leads to deterioration of computational accuracy or iteration divergence, while avoiding premature or unnecessary transformations, as the preset threshold for formal full-scale 3D simulation. This threshold is usually expressed as a dimensionless geometric distortion ratio (e.g., 0.3 to 0.7) or a physically meaningful critical equivalent plastic strain value (e.g., 0.6 to 0.9), and its specific value ultimately depends on the actual discretization accuracy and material parameters of the projectile's finite element model.

[0060] For example, for a specific hexahedral finite element in the projectile's head, its updated equivalent plastic strain data is read as 0.85. Simultaneously, based on the coordinates of its current eight nodes, the ratio of its volume to the initial volume, i.e., the Jacobian ratio, is calculated to be 0.25. The preset comprehensive evaluation function takes the larger value between the equivalent plastic strain and the geometric distortion coefficient, where the geometric distortion coefficient is defined as 1 - the Jacobian ratio. The calculated equivalent plastic strain of this element is 0.85, and the geometric distortion coefficient is 1 - 0.25 = 0.75; therefore, its distortion metric data is determined to be 0.85. The system's preset conversion threshold is 0.80. Since 0.85 is greater than 0.80, the element is determined to require conversion. Subsequently, an element conversion instruction is generated, recording "At time step T=105, perform FEM to SPH conversion on element ID 2047". This instruction will be passed to the next sub-step, triggering the specific element conversion and data mapping operations.

[0061] It should be noted that when the distortion measurement data does not exceed the preset threshold, the distortion degree of the finite element will be determined to be within an acceptable range, and no type conversion is required. Therefore, no conversion instruction data will be generated for the element. The element will continue to be retained in the projectile finite element model data and will directly participate in the subsequent sub-step S65 (failure particle identification and marking) and subsequent calculations with its current state. The entire process will not generate branch operations because of it.

[0062] Sub-step S64: In response to the element conversion command data, the current mass, momentum, energy and stress state data of the finite element to be converted are mapped to the newly generated SPH particle set to form new projectile SPH particle data, and this part of the data is removed from the projectile finite element model data.

[0063] Within the time-step iterative solution loop, a dynamic element transformation and data mapping sub-step is executed. The purpose of this step is to seamlessly transform severely distorted finite element elements from their original mesh-based computational framework to a smooth particle-hydrodynamic computational framework, in response to the element transformation command data generated in the previous step, while ensuring the conservation and continuity of physical quantities during the transformation process. This data refers to the current mass, momentum, energy, and stress state data of the finite element to be transformed, as identified in the element transformation command data.

[0064] First, the element conversion instruction data is received and parsed to obtain a list of unique identifiers for all finite element elements to be converted. For each target element in the list, the following operations are performed: First, all physical state data of the element at the current time step are read, including the element's total mass, total momentum calculated based on nodal velocity and mass distribution, total energy characterized by internal energy, and stress tensor, strain, damage, and temperature data at the element integration points. Then, according to preset particle generation rules, such as at the element centroid or integration point location, a new set of SPH particles is created. The number of particles created is usually related to the number of element nodes or integration points to ensure a reasonable mass distribution.

[0065] Next, the core data mapping operation is performed. Based on the principles of mass, momentum, and energy conservation, the total mass, momentum, and energy of the finite element to be converted are allocated to the newly generated SPH particles according to a certain weight ratio. Simultaneously, the stress state, damage variables, and temperature of the element are also assigned to these new particles accordingly. This process is completed through a mapping algorithm, ensuring that key physical quantities remain strictly conserved before and after the conversion. After mapping, the model data structure is updated. The newly generated SPH particle data, carrying all physical properties, is formally added to the system's SPH particle set, enabling it to participate in the SPH method calculation in subsequent time steps. At the same time, the successfully converted finite element is removed from its finite element model data, including deleting its identifier from the element list and severing its connections with surrounding elements, thus completing the transition of the element from the finite element system to the SPH particle system.

[0066] For example, upon receiving an instruction to convert the hexahedral element with element ID 2047, the program reads the current data of that element and determines its mass to be 2.5e. -6 kilograms, momentum vector is (1.5e) -3 The energy level is 0.01 joules, the equivalent stress is 500 MPa, the damage is 0.1, and the temperature is 350 Kelvin. According to the rules, the program generates eight new SPH particles near the eight nodes of this unit. Then, the program will have a total mass of 2.5e... -6 The kilogram was evenly distributed among eight particles, with each particle receiving 3.125e. -7The mass is kilogram. According to the law of conservation of momentum, the total momentum is also distributed proportionally by mass, ensuring that the initial velocity vectors of all eight particles are aligned with the original momentum direction of the element. Internal energy, stress, damage, and temperature data are also interpolated or averaged onto these particles. Subsequently, the IDs and properties of these eight new particles are registered in the global projectile SPH particle management list. Simultaneously, element 2047 is removed from the finite element list, and its original node connection information is cleared. From this point onward, this portion of the material will exist and evolve as eight interacting SPH particles in subsequent calculations.

[0067] Sub-step S65: Based on the updated damage variable data and stress state data of the target plate SPH particles, identify and label the failure particle data; wherein, the failure particle data stops participating in the constitutive calculation but continues to retain its mass and momentum data.

[0068] Within the time-step iterative solution loop, a sub-step for identifying and processing SPH particle failures in the target plate is executed. The purpose of this step is to determine, based on the physical damage state of the material, which material micro-elements represented by SPH particles in the target plate have fractured and failed, and to perform special numerical processing on these failed particles, thereby reasonably reproducing the physical phenomena of crack initiation and propagation in the numerical simulation, while maintaining the overall stability and conservation of the computational system.

[0069] First, the damage variable data and stress state data of all target plate SPH particles are acquired after the update at the current time step. The damage variable data characterizes the degree of damage accumulation within the material element represented by the particle, with values ​​between zero and one; the closer to one, the more severe the damage. The stress state data is used to distinguish whether the particle is currently under compression or tension. Subsequently, a logical judgment process is performed on each target plate SPH particle. This process is based on two key criteria. The first criterion is the damage criterion, which checks whether the particle's damage variable data has reached or exceeded a pre-set damage threshold based on the material properties. The second criterion is the stress state criterion, which analyzes the particle's current stress tensor to determine whether it is under net tensile stress (usually determined by the first invariant stress or the maximum principal stress being greater than zero).

[0070] A particle is considered to have undergone macroscopic fracture failure if and only if it simultaneously satisfies both of the above criteria. For all particles determined to be failed, their state flags are modified in system memory to mark them as failed particle data, with a clear meaning for subsequent calculations. Particles marked as failed will no longer participate in any form of material constitutive response calculation in subsequent time steps. This means that their internal state variables, such as stress tensor, strain, damage, and temperature, will be frozen and no longer updated, and they are generally considered to no longer participate in interactions during contact searches. However, to ensure that the mass and momentum conservation of the entire simulation system are not affected by local material failure, the mass and momentum data of these failed particles will be retained and tracked separately. These failed particles, like discrete material points, continue to move under inertia and collide with other unfailed particles or projectiles, but they no longer generate new stresses themselves. In this way, the total mass and total momentum are strictly maintained during the process of material fracture, crack formation, or fragmentation.

[0071] For example, consider a target plate SPH particle, designated P12345. The particle's current updated damage variable is 0.31, while the preset damage threshold is 0.30. Simultaneously, its stress tensor is analyzed, and its mean stress (the first invariant of stress) is calculated to be +50 MPa (a positive value indicates net tension). Since the particle's damage variable of 0.31 exceeds the threshold of 0.30, and it is under tensile stress, particle P12345 is determined to have failed. The program then changes its status from "active" to "failed." In the next time step, when performing material constitutive calculations, particle P12345 will be skipped; its stress value will remain at the pre-failure level of 500 MPa and will not be recalculated, and its damage value will also remain at 0.31 and will not increase further. However, the particle still has a mass of 1.2e-7 kg, and its momentum data, determined by its pre-failure velocity vector (e.g., 200 m / s along the X direction), will continue to be retained. The failed particle will continue to move based on this momentum, and it can still transfer momentum if it comes into contact with other objects, but no new stress response will be generated inside it.

[0072] Sub-step S66: Based on the current total contact force data and internal stress data, solve the system motion equations by integration, output the displacement and velocity data of all calculated nodes in the next time step, and update the system state data.

[0073] Within the time-step iterative solution loop, the system motion equations integration and state advancement sub-steps are executed. The purpose of this step is to calculate the acceleration of all discrete objects based on Newton's second law, according to all external and internal forces acting on the system at the current time step, and advance their motion state to the next time step through numerical integration, thereby driving the dynamic evolution of the entire projectile-target impact system.

[0074] First, the force data is summarized. The total contact force data obtained from the contact calculation stage is combined with the nodal internal force data composed of internal stress gradients obtained from the constitutive response stage. For each calculation node in the system, whether it belongs to a finite element node of the projectile or target plate, or an SPH particle, the net resultant force vector acting on it is calculated. This resultant force vector is equal to the vector sum of all external contact forces and the nodal forces contributed by the internal stress of the elements or particles associated with that node.

[0075] Subsequently, the discrete form of the equations of motion for each node is solved. Specifically, the calculated net resultant force vector of the node is divided by the mass represented or allocated by that node to obtain the instantaneous acceleration vector of that node at the current time step. Next, an explicit time integration algorithm, such as the central difference method, is used to integrate the acceleration over time. This algorithm uses the velocity data of the current time step, the newly calculated acceleration data, and the known time step size to first update the velocity vector of the node, obtaining the intermediate or final velocity data for the next time step. Then, the updated velocity data is used to integrate over time to calculate and update the spatial position coordinates of the node, thereby obtaining the displacement data for the next time step. The system equations of motion refer to the discretized numerical equations based on Newton's second law, describing the quantitative relationship between the resultant force and the acceleration of all finite element nodes and SPH particles in the projectile-target coupling system. It is the core mathematical expression that transforms the mechanical principles into an executable time integration algorithm to advance the dynamic evolution of the entire system.

[0076] After updating the velocity and displacement of all nodes, a crucial system state update operation is performed. The newly calculated node displacement and velocity data, along with the previously updated element and particle stress, strain, damage variable, and temperature data, are written to the system state data storage area, overwriting the state of the old time step. Simultaneously, the current physical time of the simulation is incremented by one time step. This completes all state transitions from one physical moment to the next, and the updated system state data will become the initial conditions for all calculations at the start of the next iteration.

[0077] For example, for a specific SPH particle on a target plate, force composition is first performed. Assume the particle experiences a normal force of 150 Newtons (along the negative Z-axis) from the contact with the projectile, while the stress field of its neighboring particles generates an internal nodal force of 20 Newtons (along the positive Z-axis) on this particle. Therefore, the net resultant force acting on this particle is 130 Newtons (along the negative Z-axis). The mass of this particle is 1.0e. -7 Kilograms. According to Newton's second law, its acceleration is calculated to be 1.3 e. 9Meters per square second (along the negative Z-axis). Assume the current time step is 1.0e. -9 The particle's previous velocity was -200 meters per second (along the Z-axis). Using the central difference method, its new velocity is updated to: [-200 + 1.3e] 9 * 1.0e -9 The velocity is -199.87 meters per second. Its new position is updated to: original Z-coordinate + (-199.87 * 1.0e) -9 (meters). This new velocity and position, along with the particle's existing stress, damage, and other data, are recorded as its system state at the next time step. The system's global time is then updated from 10.0 microseconds to 10.001 microseconds. After all nodes and particles have completed this type of calculation, the entire simulation advances one time step.

[0078] S7. When the simulation termination condition is met, exit the time-step iterative solution loop, and generate and output the target plate damage morphology data, projectile penetration result data and ballistic limit velocity based on the final system state data.

[0079] In this step, after the iterative solution loop exits due to the preset simulation termination condition, the simulation result generation and output step is executed. The purpose of this step is to post-process and analyze the system state data of the final time step, extract and generate key result data that can intuitively and quantitatively reflect the physical outcome of this impact simulation, and complete the transformation from numerical calculation to engineering conclusion.

[0080] First, the final system state data at the end of the loop is read and parsed. This data fully includes all physical field information such as the spatial position, velocity, stress, strain, and damage of all computational nodes, finite element elements, and SPH particles at the end of the simulation. Based on this data, multiple analysis tasks are executed in parallel to generate three types of core results.

[0081] When generating target plate failure morphology data, the main processing steps involve the SPH particle data and peripheral finite element node data of the target plate. Based on the spatial coordinates of the particles, their damage state, and indicators of failure, the three-dimensional geometry of the target plate after impact is reconstructed. By analyzing the distribution of the particle swarm, missing areas (regions where failed particles have been removed), and the deformation of the element mesh, typical failure modes such as disc-shaped indentations, petal cracks, or plug perforations can be identified, and the morphology data is output in the form of three-dimensional point clouds, surface meshes, or specific identifier files.

[0082] When generating projectile penetration data, the focus is on the state of all finite element elements and SPH particles within the projectile. The final spatial position of the projectile's center of mass, average velocity vector, and remaining kinetic energy are calculated. By analyzing the overall shape of the projectile material (especially severely bulged or fragmented sections), the final deformation morphology data of the projectile can be output. These data collectively describe key results such as the projectile's penetration depth, whether penetration occurred, and the remaining velocity after penetration.

[0083] Generating the ballistic limit velocity is an inference process based on a series of simulations or a single simulation result. It may be possible to infer the critical velocity value that gives the projectile a 50% probability of penetrating the target plate, i.e., the ballistic limit velocity, by querying the initial velocity and final penetration result of the current simulation (such as "not penetrated" or "fully penetrated"), and combining it with a preset velocity threshold search algorithm (such as the binary search method), or by fitting a series of existing simulation results with different initial velocities. This predicted value is then used as the data output.

[0084] For example, the simulation terminates when the projectile completely penetrates the target plate. Reading the final state data reveals a through-hole formed by numerous failed SPH particles in the central region of the target plate, with surrounding particles exhibiting an outward-curving petal-like shape. The outer finite element mesh undergoes overall bulging deformation, generating a "plug-type penetration" failure morphology data file. For the projectile, the average velocity of all its components is calculated, yielding a remaining velocity of 150 meters per second, in the same direction as the initial velocity. Simultaneously, the projectile's mushroom-shaped morphology is reconstructed based on the particle and element coordinates, and this information is output as the projectile penetration result data. If this simulation is part of a series, with an initial velocity of 737 meters per second and successful penetration, while the previous simulation with an initial velocity of 665 meters per second failed to penetrate, the ballistic limit velocity can be preliminarily determined to be between the two. The predicted value is calculated to be approximately 670 meters per second, and this value is output as the ballistic limit velocity.

[0085] In one embodiment, S2, based on the projectile's geometric model data, generates three-dimensional finite element mesh model data of the projectile, including: S21. Analyze the projectile geometric model data to obtain specific dimensional parameters characterizing the layered structure of the projectile jacket, lead sleeve, and steel core. S22. Based on specific size parameters, a three-dimensional projectile solid model with a layered structure is generated through geometric reconstruction, and different material layer regions are identified. S23. Perform hexahedral finite element mesh generation on the three-dimensional projectile solid model to generate node coordinate data and element connection relationship data; S24. Based on the spatial position of the element in the three-dimensional projectile solid model, establish the association data between the element and the material layer region; S25. Based on node coordinate data, element connection relationship data, and the association data between elements and material layer regions, generate and output the three-dimensional finite element mesh model data of the projectile.

[0086] In step S21, the projectile geometry model data parsing sub-step is executed. The purpose of this step is to accurately identify, separate, and extract all independent parameters required to define the projectile's composite geometry from the structured information source of the input projectile geometry model data, particularly the key dimensional parameters describing its multi-layered material composition. The input data is scanned and interpreted according to predefined data format specifications. Fields, keywords, or section divisions in the data are identified, and dimensional parameters describing the overall outline and those describing the internal structure of each layer are categorized and read separately. For example, the program will identify and extract values ​​that identify the overall length of the projectile, and also identify and extract values ​​that specifically describe the diameter of the steel core, the thickness of the lead sheath, etc., belonging to specific material layers. All extracted parameters are converted into internal floating-point or integer variables and logically grouped and stored according to the geometric features they describe, providing accurate numerical input for the next step of geometric reconstruction.

[0087] For example, the projectile geometry model data is read from a text file. One line in the file states "overall_length = 10.75". The parsing program identifies "overall_length" as a keyword, extracts its value of 10.75, and stores it in a variable. Another line states "core_diameter = 5.80". "core_diameter" is identified as a keyword for the steel core layer, and its value of 5.80 is extracted and stored in another variable. Continuing the parsing, all other parameters are extracted, including the projectile surface radius of 40.00, the steel core semi-cone angle of 11.9, and the jacket thickness of 0.5. These parameters are then categorized into different logical sets such as "overall dimensions," "steel core dimensions," and "jacket dimensions," forming a complete and usable set of specific dimensional parameters.

[0088] In step S22, the geometric reconstruction and identification sub-step of the projectile's three-dimensional solid model is performed. The purpose of this step is to utilize the specific dimensional parameters obtained from the analysis, through assisted design of the geometric kernel, to accurately construct a three-dimensional, solid projectile model in digital space that reflects the real physical structure, and to logically distinguish regions composed of different materials. The geometric contours of each structural layer (such as the steel core, lead sheath, and jacket) are calculated based on the parameters, for example, the radius and height of a cylinder, the dimensions of a cone, etc. Three-dimensional solids representing the steel core, lead sheath, and jacket are created sequentially from the inside out. These solids are combined into a complete projectile solid model through Boolean operations (such as addition and fusion). When creating each solid or after combination, a unique, user-defined attribute identifier is assigned to the geometric surfaces or volumes belonging to the same material layer, such as "MaterialLayer_Core" or "MaterialLayer_Jacket". This identifier is not a geometric feature, but rather metadata attached to the geometry, used to trace the physical material represented by that part of the geometry in subsequent steps.

[0089] For example, using dimensional parameters, first create a solid with a bottom diameter of 5.80 mm and a conical section as the steel core, and immediately assign the attribute label "Core" to all its surfaces. Next, using the outer surface of the steel core as a reference, offset outwards to determine the thickness of the lead sleeve, generating a lead sleeve solid and assigning it the label "Lead". Finally, using the outer surface of the lead sleeve as a reference, offset outwards by 0.5 mm to generate the jacket solid and assign it the label "Jacket". Ultimately, these three solids are combined into a seamless, single 3D solid model of the projectile, but different regions within the model carry the labels "Core", "Lead", and "Jacket" respectively.

[0090] In step S23, the hexahedral finite element mesh generation sub-step is executed. The purpose of this step is to spatially discretize the reconstructed 3D projectile solid model, transforming it into a mesh composed of numerous small hexahedral elements connected by nodes, providing a topological basis for finite element calculations. A mesh generation algorithm is invoked, which automatically arranges a series of nodes inside and on the surface of the solid model based on the complexity of the model's geometry and user-defined mesh accuracy requirements (such as global element size or curvature control parameters). Then, the algorithm connects these nodes according to certain rules (such as sweeping, mapping, or multi-domain meshing) to form hexahedral elements. This process is entirely automated, ultimately generating two sets of core data: one is a list of the 3D spatial coordinates of all nodes, i.e., node coordinate data; the other is a list defining which nodes are connected to each hexahedral element in what order, i.e., element connection relationship data. These data collectively describe the discrete geometry of the model.

[0091] For example, a sweep mesh generation process is initiated on the assembled projectile solid model. The program divides the model into several layers along the projectile's axis, and then further divides the cross-section of each layer into a ring-shaped mesh. Ultimately, the entire projectile is discretized into 5872 hexahedral elements, which are composed of more nodes (e.g., approximately 7000). The program records the node number and its corresponding X, Y, Z coordinates for each element, as well as the element number and which 8 nodes (for hexahedral elements) constitute it.

[0092] In step S24, a sub-step is performed to establish the association data between the execution element and the material layer region. The purpose of this step is to link the purely geometric mesh elements with the physical material they represent, bridging the gap for subsequent assignment of material properties. Each hexahedral element generated after meshing is traversed. For each element, the program determines its material layer region by calculating the spatial location of its geometric center point in the 3D projectile solid model, or by detecting which identified material layer solid contains the majority of the element's volume. This determination is typically performed through a geometric spatial inclusion test or by utilizing the original geometric association information retained during meshing. After determining the region, an association entry is created for each element, mapping its unique identifier to the identifier of the material layer region it belongs to (e.g., "Core", "Jacket"), forming the association data between the element and the material layer region.

[0093] For example, consider element number 5678. The centroid coordinates of this element are calculated to be (1.0, 0.2, 8.5) mm. Based on these coordinates, the program performs a spatial lookup within the identified projectile solid model and finds that this coordinate point lies within the geometry labeled "Jacket". Therefore, the program creates a record in the association data table: Element 5678 -> Material layer "Jacket". After performing this operation on all elements, a complete mapping table is created, explicitly indicating which physical material layer each mesh element corresponds to.

[0094] In step S25, the sub-step of generating and outputting 3D finite element mesh model data is executed. The purpose of this step is to integrate all intermediate data generated in previous steps, assemble and output a structurally complete and information-rich 3D finite element mesh model data object for use in subsequent steps. The nodal coordinate data, element connection relationship data, and element-material layer region association data are encapsulated and integrated. This data is organized according to a predetermined data structure (which may be custom-defined or conform to a general simulation software interface specification), ensuring that the reference relationships between them are correct. For example, it is ensured that the element ID referenced in the association data exists in the element connection relationship data, and the node ID referenced in the element connection data exists in the nodal coordinate data. Finally, this composite data object integrating geometric topology information and material layer attribution information is output to a specific variable in memory or written to an intermediate file, marking that the 3D finite element mesh model data of the projectile is ready.

[0095] For example, create a data object named "Bullet_Mesh_Model". This object consists of three main data arrays: the first array is a list of nodes, storing their X, Y, and Z coordinates in order of node ID. The second array is a list of elements, storing the sequence of node IDs they contain in order of element ID. The third array is an association list, storing the corresponding material layer identifier string (such as "Jacket") in order of element ID. After the program performs integrity checks, it passes this data object to the next step (S3, material property assignment) as input.

[0096] In one possible embodiment, S5, based on the input target plate geometric model data and material parameter data, generates hybrid discrete model data of the target plate, including: S51. Based on the target plate geometric model data and the expected impact point position, define the core region data and the outer region data surrounding the core region; S52. Based on the core region data and the preset particle spacing data, generate SPH particle set basic data within the core region; wherein, the SPH particle set basic data includes the initial coordinates and mass data of the particles. S53. Based on the peripheral region data, the core region boundary data, and the interface unit size data that matches the particle spacing data, generate the basic data of the FEM mesh, wherein the unit size of the FEM mesh gradually increases from the interface with the core region to the periphery. S54. Based on the particle data of the adjacent interface in the SPH particle set basic data and the cell node data of the interface in the FEM mesh basic data, establish the SPH-FEM consolidation connection relationship data. S55. Assign the target plate material parameter data to all particles in the SPH particle set base data and all elements in the FEM mesh base data, and initialize their stress, strain and damage variable data to form particle data and element data with material properties. S56. Apply fixed constraints to the nodes located at the geometric boundary of the target plate in the FEM mesh base data, and output the hybrid discrete model data of the target plate based on the particle data, element data and consolidation connection relationship data with material properties.

[0097] In step S51, a region definition sub-step is performed based on the target plate geometric model data and the expected impact point location. This step aims to logically divide the target plate geometric space into two computational regions using different discretization strategies, according to the characteristics of the impact physics process, laying the foundation for subsequent hybrid modeling. First, the target plate geometric model data is parsed to obtain its overall shape, size, and spatial orientation. Simultaneously, the coordinates of the expected impact point, specified by the user or determined based on the impact scenario, are read. This point is typically located at the geometric center of the target plate or a specific location. Subsequently, a continuous sub-region is defined within the three-dimensional entity of the target plate, centered on this impact point, according to preset rules (such as fixed side length, radius, or a thickness-based ratio). This region is marked as the core region. The core region is typically designed to encompass all large deformations, damage, and even failures that may be caused by the projectile impact. The remaining parts within the target plate geometry, excluding the core region, are uniformly defined as the outer region. Data structures describing the spatial extent of these two regions are generated, for example, through bounding box coordinates or geometric Boolean operation results.

[0098] For example, if the target plate is read as a square plate with dimensions of 250 mm x 250 mm and a thickness of 4.1 mm, and the expected impact point is located at the center of the plate, the program, according to preset rules, defines a square prism region with a side length of 27 mm and a thickness that extends through the center point as the core region. The rest of the target plate is automatically defined as the outer region surrounding this square prism. Two sets of data will be generated: one set describes the boundary of the core region (e.g., the range in the X and Y directions is -13.5 mm to +13.5 mm, and in the Z direction is 0 to 4.1 mm), and the other set describes the outer region formed by the difference between the entire target plate and this core region.

[0099] In step S52, based on the defined core region data and preset particle spacing data, the SPH particle set generation sub-step is executed. The goal of this step is to geometrically discretize the continuous core region into a set of discrete points, i.e., SPH particles, that follow the calculation rules of the smooth particle hydrodynamics method. First, based on the spatial range and shape of the core region, combined with the preset particle spacing, the total number of particles required for a regularly arranged particle array within the region and their initial position coordinates are calculated. Typically, particles are initialized in three-dimensional space in a cubic or hexagonal close-packed manner. For each generated particle, the mass assigned to that particle is calculated based on its represented discrete volume (determined by the particle spacing) and the density of the target material. Finally, the basic SPH particle set data is generated. This data is a structured list, where each record corresponds to a particle, including its unique identifier, three-dimensional initial coordinates, and calculated mass value.

[0100] For example, the core region is 27x27x4.1 cubic millimeters. The preset particle spacing is 0.5 millimeters. First, it is calculated that 54 particles can be arranged in the X direction, 54 in the Y direction, and 9 in the thickness direction, generating a total of 26244 particles. The program assigns an initial coordinate to each particle; for example, the first particle's coordinates are (-13.25, -13.25, 0.25) millimeters. Assuming the target material density is 7850 kg / m³, and the volume represented by a single particle is (0.0005 m³), ​​its mass is approximately 9.81e-10 kg. The ID, coordinates, and mass information of all these particles are recorded as the SPH particle set basic data.

[0101] In step S53, based on the peripheral region data, core region boundary data, and interface element size data, a sub-step for generating the peripheral region FEM mesh is performed. This step aims to generate a finite element mesh for the peripheral region and ensure its compatibility and computational efficiency with the core SPH region at the interface. First, at the geometric interface between the core and peripheral regions, a fine-grained finite element surface mesh is generated based on the interface element size data matching the SPH particle spacing. Subsequently, starting from this interface, a mesh transition technique is used to expand the three-dimensional volume mesh towards the outer edge of the target plate. The mesh transition is achieved by gradually increasing the element size, for example, by increasing the element side length by a certain percentage gradient. Finally, the basic FEM mesh data is generated, which includes the coordinate data of all nodes in the peripheral region, the connection relationship data of all elements, and its mesh density exhibits a distribution characteristic from dense to sparse from the inside out.

[0102] For example, the interface between the core and outer regions is a 27*27 mm square ring. Quadrilateral surface elements with a size of 0.5 mm are generated on this ring. Starting from this surface mesh, three-dimensional hexahedral elements are generated by stretching or sweeping towards the outer side (four edges) and back of the target plate, while the element size gradually increases; for example, the element size increases to 2 mm at a distance of 10 mm from the interface and to 5 mm at the edge of the target plate. Finally, a finite element mesh for the outer region with a total of approximately 13056 elements is formed, and its node coordinates and element connection relationships are recorded as the basic FEM mesh data.

[0103] In step S54, a sub-step for establishing interface consolidation connections is performed based on the particle data of the adjacent interfaces in the SPH particle set base data and the element node data of the interface in the FEM mesh base data. This step is crucial to ensuring the continuity of the mechanical behavior of the hybrid discrete model at the interface, aiming to establish a physical quantity transfer mechanism between the two discretization methods. First, SPH particles located at the core region boundary, i.e., adjacent to the outer region, and element nodes in the outer region mesh that connect to the core region are identified. Then, a specific consolidation algorithm is invoked, such as a constraint-based or weight-based algorithm. This algorithm finds several (usually one or more) action nodes among the adjacent finite element nodes of each SPH particle near the interface and calculates a set of weight coefficients. These weights define how the motion and forces of the SPH particle are constrained or influenced by these finite element nodes, and vice versa. All calculated particle-node pairs and their corresponding weight coefficients are stored as SPH-FEM consolidation connection data. The specific consolidation algorithm is a numerical method for establishing multi-point constraint relationships between SPH particles and adjacent FEM element nodes. Its core lies in assigning a set of constraint weights to one or more adjacent finite element nodes for each SPH particle near the interface, based on its spatial position. This algorithm uses these weights to construct constraint equations between particle displacements and node displacements, ensuring that the displacements on both sides of the interface are consistent in subsequent solutions, and that action and reaction forces are transmitted between particles and nodes according to the same weight relationship. This achieves seamless mechanical coupling between two heterogeneous numerical regions at the interface, forming a unified and coherent computational system. The constraint weights are typically determined based on the distance between the SPH particle and the relevant finite element node or the shape function value of the element.

[0104] For example, consider an SPH particle located at the upper surface boundary of the core region, with coordinates (13.0, 0, 2.0) mm. In the adjacent peripheral finite element mesh, its four nearest surface nodes are found. The consolidation algorithm calculates the distances between the particle and these four nodes and determines a set of weights based on the inverse distance ratio or other functions, for example, weights of 0.4, 0.3, 0.2, and 0.1 for the four nodes. This relationship is recorded as a consolidation connection, indicating that in subsequent calculations, the particle's displacement will be constrained by the weighted average of the displacements of these four nodes, and the forces acting on the particle will also be distributed and reacted upon these four nodes according to this weight.

[0105] In step S55, the target plate material parameter data is assigned to all particles in the SPH particle set base data and all elements in the FEM mesh base data, and a state initialization sub-step is executed. This step allows the discrete geometric objects to acquire their material mechanical behavior definition and initial physical state. The target plate material parameter data is read, which includes a complete constitutive model (such as the modified Johnson-Cook model) and all its constants. Then, each particle in the SPH particle set base data is traversed, and this set of material parameters is completely associated with that particle. At the same time, each element in the FEM mesh base data is traversed, and the same set of material parameters is associated with that element. This ensures the consistency of the material between the core region and the outer region. After the association is completed, the initial mechanical state variables of all particles and elements are set. Typically, the components of the stress tensor are set to zero, the cumulative plastic strain is set to zero, the damage variable is set to zero, and the initial temperature value is set according to the reference temperature, thereby forming particle data and element data with material properties.

[0106] For example, the target material is 30CrMnSiA steel, and its material parameters include Johnson-Cook strength parameters A=525 MPa, B=101 MPa, n=0.081, C=0.1739, m=1.635, and damage parameter D1=0.0705. This complete parameter set is simultaneously assigned to all 26244 SPH particles and all 13056 FEM elements. Subsequently, the equivalent stress, plastic strain, and damage values ​​of all particles and elements are initialized to 0, and the initial temperature is 293 Kelvin.

[0107] In step S56, fixed constraints are applied to the nodes located on the geometric boundaries of the target plate in the FEM mesh base data, and the hybrid model data integration output sub-step is executed. This is the final step in generating the final computable model. Nodes located on the physical boundaries of the target plate in the finite element mesh of the outer region are identified, such as the four sides assumed to be clamped or fixed in the simulation, and the possible rear surface. Applying fixed constraints to these nodes typically means restricting their displacement degrees of freedom in the X, Y, and Z directions. Finally, the particle data with material properties, the element data with material properties and constraint information, and the SPH-FEM consolidation connection relationship data are encapsulated and integrated to form a unified, self-consistent data structure: the hybrid discrete model data of the target plate. This data object includes all the geometric, discretization, material, connection, and boundary information required for subsequent coupled dynamics simulations.

[0108] For example, all nodes in the finite element mesh of the outer region with X coordinates equal to -125 mm or +125 mm (target side) and Y coordinates equal to -125 mm or +125 mm are identified, and the UX, UY, and UZ degrees of freedom of these nodes are constrained to zero. Then, the program creates a data object named "Target_Hybrid_Model", which internally links to or includes an initialized SPH particle data set, a constrained FEM element data set, and an interface consolidation connection table. This data object, as the output of step S5, is used together with the projectile model data to initiate the impact simulation.

[0109] In one possible embodiment, S61, based on the current spatial position data of the projectile and the target particles, performs contact search and generates contact pair data, and calculates and outputs contact force data based on the contact pair data, including: S611. Based on the position data of the finite element nodes or SPH particles of the projectile and the position data of the SPH particles of the target plate, perform a global contact search to identify all potential contact object pairs. S612. Perform local precision detection on each potential contact object pair. If the projectile object is a finite element, calculate the penetration depth and contact normal of the target plate particle on the surface of the element. If the projectile object is an SPH particle, calculate the center distance between the two particles and generate contact pair data including contact type and contact geometry accordingly. S613. Based on contact pair data and the constitutive parameter data of the associated projectile and target plate materials, parallel calculations are performed for different contact types: For the contact between the projectile finite element and the target plate SPH particles, the contact force is calculated using the penalty function method based on the penetration depth and material stiffness; For the contact between the projectile SPH particles and the target plate SPH particles, the contact force is calculated using the SPH particle interaction model based on the particle spacing, kernel function derivative, and constitutive stress, thereby generating contact force metadata corresponding to each contact pair. S614. Distribute and integrate contact force metadata: For finite element-particle contact, distribute the contact force to the corresponding finite element nodes according to the shape function, and distribute the reaction force to the target plate SPH particles; for particle-particle contact, distribute the interaction force directly to the corresponding paired SPH particles, and finally summarize to generate contact force data acting on all relevant nodes and particles in the system.

[0110] In step S611, a global contact search sub-step is executed. The purpose of this step is to quickly and efficiently filter out all possible contact object pairs based on the instantaneous spatial positions of all objects involved in contact at the current time step, thereby narrowing the scope of subsequent precise detection and improving overall computational efficiency. First, the position data of all possible contact objects in the system is acquired, including the coordinates of the finite element nodes on the projectile surface, the coordinates of the transformed SPH particles of the projectile, and the coordinates of all target plate SPH particles. Then, a spatial partitioning algorithm is invoked, such as a bounding box-based algorithm or a spatial grid hashing algorithm. This algorithm divides the entire computational space into many small cells or regions and quickly classifies objects into corresponding cells based on their positions. Next, within each cell and between adjacent cells, the algorithm checks whether the distance between objects from different objects (e.g., a projectile node and a target plate particle) is less than a preset coarse search radius. All object pairs that meet this coarse distance condition are identified and recorded as potential contact object pairs. This process outputs a preliminary list including all candidate object pairs that require further precise geometric detection.

[0111] For example, at a certain time step, there are tens of thousands of target plate SPH particles and thousands of projectile surface nodes and particles. Using spatial grid hashing, the entire computational domain is divided into a cubic grid with sides of 1 mm. A projectile node located at coordinates (10.1, 0.5, 2.0) mm is placed in the corresponding grid (10,0,2), and a target plate particle located at (10.2, 0.5, 2.1) mm is placed in an adjacent grid. The algorithm checks the grid containing the projectile node and its 26 surrounding adjacent grids. It finds that the target plate particle is located within an adjacent grid, and the approximate distance between them is less than the preset search radius of 1.5 mm. Therefore, this pair (projectile node, target plate particle) is recorded as a potential contact pair. The system performs this operation on all objects, ultimately generating a list that may include thousands of candidate pairs.

[0112] In step S612, a sub-step for local precise detection and contact pair data generation is performed. The purpose of this step is to perform a fine-grained geometric relationship judgment on each potential contact object pair output by the global search, confirming whether actual contact has occurred and quantifying the geometric characteristics of the contact. Each potential contact object pair is traversed. For each pair, the program first determines the type of the projectile object. If the projectile object is a node of a finite element unit, the program needs to find the surface of the element to which the node belongs and calculate the shortest distance from the target particle to this element surface, i.e., the penetration depth, while simultaneously determining the unit normal vector at the contact point. If the projectile object itself is an SPH particle, the program directly calculates the Euclidean distance between the centers of the two particles. Based on the calculation results and the preset contact tolerance, the program determines whether the contact has actually occurred. For example, for finite element surface contact, if the penetration depth is positive; for inter-particle contact, if the center distance is less than the sum of the smooth lengths of the two particles. For all contacts determined to have actually occurred, a structured contact pair data record is generated. Each record includes unique identifiers of both parties in the contact, the contact type, and precise contact geometric quantities, such as penetration depth and contact normal vector, or inter-particle center distance.

[0113] For example, consider a potential contact pair where the projectile object is a finite element surface (a quadrilateral composed of four nodes), and the target object is an SPH particle. Calculate the projection point and distance of the particle onto the quadrilateral surface. If the particle passes through the surface with a penetration depth of 0.02 mm, and the contact point normal points towards the target particle, the program determines that contact has occurred and generates contact pair data, recorded as: Type=FEM-SPH, Projectile Element ID=1001, Target Particle ID=5500, Penetration Depth=0.02 mm, Normal Vector=(0,0,1). For another potential contact pair of SPH particles, the center-to-center distance is calculated to be 0.45 mm. Since the sum of the smooth lengths of the two particles is 1.0 mm, the center-to-center distance is less than this value, so contact is determined to have occurred, and data is generated as: Type=SPH-SPH, Projectile Particle ID=200, Target Particle ID=6000, Center-to-Center Distance=0.45 mm.

[0114] In step S613, a contact force calculation sub-step based on contact type is executed. The purpose of this step is to quantitatively calculate the interaction force at the contact point using a physical model that matches the contact type, based on the contact geometry described by the contact pair data and the material properties of both parties. Each contact pair data point and the associated constitutive parameter data of both materials are read. For the contact between the projectile finite element element and the target plate SPH particle, the program uses the penalty function method model. This model treats the penetration depth as a small, permissible geometric intrusion and assumes the existence of a virtual linear spring at the contact point. Based on the equivalent material stiffness of both parties and the penetration depth, the program calculates the normal contact force generated by this virtual spring. Simultaneously, the tangential friction force is calculated according to Coulomb's law of friction and the tangential relative velocity at the contact point. The vector sum of the normal and tangential forces is the total contact force of the contact pair.

[0115] For the contact between the projectile SPH particles and the target plate SPH particles, the SPH particle interaction model is invoked. This model is based on the kernel approximation principle of the SPH method. The program directly calculates the interaction force between the two particles according to the distance between them, the derivative of the kernel function at that distance, and the current constitutive stress tensor and density of the two particles, following the discretized form of the SPH momentum equation. This calculation naturally includes the material response. All calculated force vectors corresponding to specific contact pairs are stored as contact force metadata.

[0116] For example, for an identified FEM-SPH contact pair, the penetration depth is 0.02 mm. Assume the equivalent contact stiffness calculated based on the elastic modulus of the projectile jacket and target plate steel is 1e. 9 Newtons per meter. According to the penalty function method, the magnitude of the normal contact force is 1e. 9 N / m * 0.02e -3m = 20000 Newtons. Assuming a friction coefficient of 0.1, the tangential friction force is 2000 Newtons. Generate contact force metadata of 20100 Newtons, directed along the combined direction of the contact normal and tangential. For an SPH-SPH contact pair with a center-to-center distance of 0.45 mm, the program queries the kernel function table to obtain the kernel function gradient value at this distance. Combining the densities of the two particles and the current stress (e.g., particle A has a compressive stress of -300 MPa, and particle B has a tensile stress of 50 MPa), and substituting them into the SPH interaction formula, a repulsive force of 15000 Newtons, directed along the line connecting the centers of the two particles, is directly calculated as the contact force metadata for this pair.

[0117] In step S614, the contact force allocation and system integration steps are performed. The purpose of this step is to correctly allocate the calculated force metadata, which is still at the "contact pair" level, to each basic computational entity constituting the system, and integrate it into a complete list of node and particle forces, providing input for solving the equations of motion. All contact force metadata is traversed. For each data point, allocation is performed according to its contact type. For the contact between a finite element element and an SPH particle, the program needs to apply the calculated contact force to the target particle, while simultaneously applying equal and opposite reaction forces to the projectile finite element element. Since the force on the finite element element acts on its surface, the program needs to allocate the reaction force according to the weights of the nodes constituting the element, based on the local coordinates of the contact point on the element surface and using the element's shape function. For the contact between SPH particles, the calculated interaction forces directly correspond to the two particles, so the program directly assigns two equal and opposite forces to these two particles respectively.

[0118] After all contact force metadata has been assigned, integration is performed. A force array for all nodes and particles in the system is maintained. For each node or particle, the program may receive assigned force components from multiple different contact pairs. These force vectors from different contact pairs are summed to obtain the total contact force acting on that node or particle. Finally, a complete set of contact force data is output, clearly showing the total external contact force vector experienced by every node and every particle in the system at the current time step.

[0119] For example, for the 20,000 Newton normal force calculated for the aforementioned FEM-SPH contact pair (assuming the tangential force is zero), the 20,000 Newton force is first applied to target particle ID 5500. For projectile element ID 1001, assuming the contact point is located on the surface of this element near nodes A and B, the weights calculated by the shape function are: 0.7 for node A and 0.3 for node B. The program then distributes the -20,000 Newton reaction force, applying -14,000 Newtons to node A and -6,000 Newtons to node B. For the 15,000 Newton repulsive force calculated for the aforementioned SPH-SPH contact pair, the program directly applies a +15,000 Newton force to projectile particle ID 200 and a -15,000 Newton force to target particle ID 6000. Finally, assuming node A also receives a force of -5000 Newtons from other contacts, these are added together to obtain a total contact force of -19000 Newtons for node A. This accumulation operation is performed on all nodes and particles to generate the current total contact force data for the system.

[0120] In one possible embodiment, S62, the contact force data and the current stress, strain, and damage variable data of the elements and particles are input into the material constitutive model for calculation, and the updated stress data, strain data, damage variable data, and temperature data are output, including: S621. Based on the contact force data and the current stress data, the strain increment data and equivalent plastic strain rate data of the element and the particle are calculated. Based on the strain increment data, the current cumulative plastic strain data, and the current temperature data, new cumulative plastic strain data is obtained, and new temperature data is calculated based on the adiabatic model of plastic work-heat conversion. S622. Based on new cumulative plastic strain data, current stress triaxiality data, equivalent plastic strain rate data, and new temperature data, new damage variable data are calculated using a damage evolution model that includes damage threshold and fracture strain criterion. S623. Based on new damage variable data, new cumulative plastic strain data, equivalent plastic strain rate data, and new temperature data, the current yield stress data considering damage softening is calculated by the Johnson-Cook model coupled with damage variables. S624. Based on the current yield stress data and strain increment data, update the stress tensor through the elastoplastic stress integral algorithm, output the updated stress data, and simultaneously output the new cumulative plastic strain data, new damage variable data and new temperature data as the updated strain data, damage variable data and temperature data.

[0121] In step S621, the strain, temperature, and cumulative plastic strain update sub-steps are executed first. The purpose of this step is to calculate the basic increments of material deformation and temperature rise based on the mechanical action at the current time step, providing updated state input for subsequent damage and strength calculations. Based on the stress state jointly determined by contact force data and current stress data of unit particles, the strain increment data and equivalent plastic strain rate data for the current time step are calculated through constitutive relations. The strain increment data describes the amount of deformation change from the previous time step to the current time step, while the equivalent plastic strain rate data describes the rate of plastic deformation. Subsequently, the state quantities are updated. The calculated plastic strain increment is added to the cumulative plastic strain data from the previous time step, thereby updating the new cumulative plastic strain data, which records the total amount of irreversible plastic deformation of the material to date. Simultaneously, based on the adiabatic assumption, the plastic deformation work is largely converted into heat. Based on the specific value of the plastic work generated in the current time step, combined with the material density and specific heat capacity parameters, the temperature rise data is calculated using the plastic work-heat conversion formula. This temperature rise is then added to the temperature data of the previous time step to update and obtain new temperature data.

[0122] For example, for a target plate SPH particle, the equivalent plastic strain increment at the current time step is calculated to be 0.005 based on the force acting on it, and the equivalent plastic strain rate is 2000 per second. The program reads the cumulative plastic strain of the particle from the previous step as 0.15 and the temperature as 320 Kelvin. First, the cumulative plastic strain is updated to 0.155. Then, assuming the plastic work generated at this time step is 0.001 joules, the particle mass is 1.0e-7 kg, and the material specific heat capacity is 452 joules per kilogram of Kelvin, the temperature rise is calculated to be approximately 22 Kelvin according to the formula, thus updating the particle temperature to 342 Kelvin.

[0123] In step S622, the damage evolution and damage variable update sub-steps are then executed. The purpose of this step is to determine and calculate the cumulative amount of internal damage based on the material's current deformation and stress state. The new cumulative plastic strain data obtained in the previous step is read, and it is determined whether it exceeds a pre-set damage threshold. If it does not exceed the threshold, the damage does not evolve, and the damage variables remain unchanged; if it does exceed the threshold, damage calculation is initiated. The core of the damage evolution model used for the calculation is to determine the fracture strain of the material in its current state. This fracture strain is a function of the stress triaxiality data, the equivalent plastic strain rate data, and the new temperature data, and is typically described by a formula including constants D1 to D5. The stress triaxiality data characterizes the ratio of hydrostatic pressure to deviatoric stress in the stress state and has a decisive influence on damage evolution. Based on the current stress triaxiality, strain rate, and temperature, the fracture strain value in the current state is calculated using this formula. Then, according to the damage evolution law, the damage increment is calculated from the relationship between the cumulative plastic strain increment and the fracture strain. Finally, this damage increment is added to the damage variable data from the previous time step, and the updated damage variable data is output.

[0124] For example, the program determines that the cumulative plastic strain of the aforementioned particle, 0.155, has exceeded the damage threshold (e.g., 0.05). The program reads its current stress triaxiality data as 0.33, strain rate as 2000 per second, and temperature as 342 Kelvin. Substituting these values ​​into the fracture strain calculation formula, it assumes the fracture strain in the current state is 0.80. According to the damage evolution law, the damage increment is approximately the current plastic strain increment of 0.005 divided by the fracture strain of 0.80, multiplied by the critical damage value, assumed to be 0.00625. If the damage variable for this particle in the previous step was 0.10, then the updated damage variable data is 0.10625.

[0125] In step S623, the yield stress calculation sub-step for coupled damage is then executed. The purpose of this step is to comprehensively consider the material's strain hardening, strain rate strengthening, temperature softening, and recently evolved damage weakening effects to calculate the material's true yield stress at the current moment, which can resist deformation. The Johnson-Cook model for coupled damage variables is invoked. This model, based on the classic Johnson-Cook formula, incorporates damage effects through multiplication. New damage variable data, new cumulative plastic strain data, equivalent plastic strain rate data, and new temperature data, along with material constants A, B, n, C, and m, are substituted into the coupled model formula for calculation. The first part of the formula reflects the contribution of strain hardening, the second part reflects the contribution of strain rate strengthening, and the third part reflects the contribution of temperature softening. The coefficient multiplied by these parts reflects the reduction in the material's load-bearing capacity due to damage. This calculation outputs a current yield stress that considers all hardening and softening mechanisms.

[0126] For example, the state data of the aforementioned particles—damage variable 0.10625, cumulative plastic strain 0.155, strain rate 2000 kJ / s, temperature 342 Kelvin—and the Johnson-Cook constants of the target material (A = 525e6 Pa, B = 101e6 Pa, n = 0.081, C = 0.1739, m = 1.635, reference parameters omitted) are substituted into the Johnson-Cook formula for coupled damage. The formula first calculates the base stress before damage, then multiplies it by a factor (1 - 0.10625) = 0.89375 to reflect the strength degradation caused by damage. Finally, the current yield stress data considering damage softening is calculated, assumed to be 4.50 eJ / s. 8 Pa.

[0127] In step S624, the final step is to synchronize the stress integration and state data output. The purpose of this step is to determine the precise stress state of the material element or particle at the new time step based on the calculated current yield stress and strain increment, using a strict stress-strain constitutive relation, and to integrate all updated state variables for output. An elastoplastic stress integration algorithm is used, aiming to solve for a new stress tensor that satisfies the elastic stress-strain relationship without exceeding the yield surface defined by the current yield stress data. This algorithm typically employs a prediction correction strategy. First, based on the strain increment data, assuming the deformation is purely elastic, an elastic predicted stress is calculated. Next, it is checked whether the equivalent stress of this predicted stress exceeds the current yield stress. If it does not exceed it, the deformation is indeed elastic, and the predicted stress is the new stress. If it exceeds it, it indicates plastic flow, requiring plastic correction. The plastic correction algorithm, according to the flow law, pushes the stress exceeding the yield surface back onto the yield surface and simultaneously updates the plastic strain. After integration, the updated stress tensor is output as the new stress data. Simultaneously, the new cumulative plastic strain data, new damage variable data, and new temperature data finally determined at this time step are output as updated strain data, damage variable data, and temperature data. These data together constitute the complete initial state of the element or particle in the next time step.

[0128] For example, for the aforementioned particles, the strain increment results in an elastically predicted stress of 5.00e. 8 Pa. However, the current yield stress calculated in step S623 is 4.50e. 8 Pa. Since the predicted stress exceeds the yield stress, the program initiates a plastic correction. Based on the correlated flow rule, the algorithm calculates the plastic strain increment and corrects the stress tensor from 5.00e8 Pa back to the yield surface, ultimately obtaining a new stress value of 4.50e. 8Pa (equivalent stress value). At the same time, the cumulative plastic strain is finely adjusted from 0.155 to 0.156, the damage variable is kept at 0.10625, and the temperature is kept at 342 Kelvin. These values ​​are output together to complete all state updates for the particle in this time step.

[0129] In one possible embodiment, S63, calculating the distortion metric data based on the updated strain data of the elements in the projectile finite element model, including: S631. Extract the updated strain data of each element in the finite element model of the projectile. The strain data shall include at least the equivalent plastic strain data of the element. S632. Based on the deformed node coordinate data of each element, calculate the geometric parameter data characterizing the shape distortion of the element; S633. Using equivalent plastic strain data and / or geometric parameter data as input, calculate a distortion metric for each element that comprehensively characterizes its deformation and distortion degree.

[0130] In step S631, the computer device executes the sub-step of extracting strain data after element update. The purpose of this step is to accurately read and obtain the latest strain state of each element in the projectile finite element model at the current time step from the dataset where the material constitutive response calculation has been completed and the state has been updated, providing crucial physical quantity inputs for subsequent distortion assessment. First, a specific memory region or data structure storing the system state data is accessed. This region, indexed by element identifier, stores all updated state variables for each element. The program locates all element records contained in the projectile finite element model according to a predefined index relationship. Then, from the multiple state variables contained in each element record, the equivalent plastic strain data, one of the core judgment criteria, is identified and extracted. This data is a scalar representing the degree of irreversible deformation accumulated by plastic deformation in the element material. The computer temporarily stores or organizes these extracted equivalent plastic strain values, corresponding one-to-one with each element, into a new data list for direct retrieval in the next calculation.

[0131] For example, the finite element model of the projectile contains one thousand elements with element IDs ranging from 1001 to 2000. After completing the material constitutive calculation in S62, the updated state data (including stress, strain, damage, temperature, etc.) for each element has been saved. The computer program iterates through these one thousand elements, reading the "equivalent plastic strain" field from the record of element ID 1001, obtaining a value of 0.15; reading the same field from the record of element ID 1002, obtaining a value of 0.08; and so on, until element ID 2000. Finally, the program generates a list containing one thousand values, where the Nth item in the list is the equivalent plastic strain data corresponding to element ID (1000+N).

[0132] In step S632, the computer device executes a sub-step for calculating the geometric distortion parameters of the element. The purpose of this step is to quantify the degree of distortion of the element's shape relative to its initial regular state through geometric calculations, based on the nodal spatial coordinates updated by deformation at the current time step, providing another dimension for distortion measurement. The computer first acquires the coordinate data of all nodes of each element in the projectile finite element model at the current time step. This coordinate data has been recalculated in the system state update step. For each finite element to be evaluated, the program calls the corresponding geometric distortion evaluation function according to its type (e.g., hexahedral element). This function takes the current coordinates of the element nodes as input and obtains one or more scalar geometric parameters characterizing shape distortion through a series of vector operations and determinant calculations. Common calculations include: calculating the determinant value of the element's Jacobian matrix based on the nodal coordinates or its ratio to the initial determinant value to measure volume compression or expansion distortion; calculating the length ratio or angle change of each side of the element to measure shear distortion. The calculation results, i.e., the geometric parameter data characterizing the shape distortion of the element, are recorded and associated with the element.

[0133] For example, for a specific hexahedral element within a projectile, the computer reads the three-dimensional coordinates of its eight nodes at the current time step. The program first calculates the Jacobian matrix of the element in its current form, a matrix describing the coordinate mapping within the element. By calculating the determinant of this Jacobian matrix, a parameter measuring the change in the element's volume can be obtained. Assuming the element is initially a regular cube, the ideal value of the Jacobian determinant is 1. After deformation, the calculated current Jacobian determinant value is 0.3, so the program records the geometric parameter "Jacobi ratio" as 0.3. Furthermore, the program may also calculate the ratio of the longest to the shortest side of the element, obtaining an aspect ratio of 2.5. These values ​​(0.3 and 2.5) together constitute the geometric parameter data of the element.

[0134] In one possible embodiment, the simulation termination condition is that at least one of the following conditions is met: The spatial position data of the projectile or the target plug that moves with the projectile indicate that it has completely passed through the target plate and that its remaining velocity data is greater than zero. The total kinetic energy of the system is lower than the first preset threshold. The rate of change of the system's total energy over time is lower than the second preset threshold. The total simulated physical time data reached the preset maximum duration.

[0135] The simulation termination condition is a series of quantifiable and detectable physical or numerical criteria used during the runtime step-by-step iterative solution loop to determine whether to exit the loop and end the simulation. Its purpose is to ensure that the simulation automatically stops when the physical process has substantially completed or the computation has stabilized, avoiding unnecessary consumption of computational resources and outputting meaningful final results.

[0136] 1. Condition 1: The projectile or target plug completely penetrates the target plate and the remaining velocity is greater than zero.

[0137] This condition applies to scenarios where the projectile penetrates the target plate with sufficient kinetic energy. At each time step or every few time steps, the spatial position data of all parts of the projectile and the target plate impingement portion that adheres to and moves with the projectile are monitored in real time. A typical method for determining "complete penetration of the target plate" is to check whether the coordinates of all particles of the aforementioned moving body in the direction perpendicular to the thickness of the target plate have exceeded the coordinates of the target plate's back surface (i.e., the side opposite the impact surface). Simultaneously, the mass-weighted average velocity of the moving body, i.e., the residual velocity data, is calculated. When the spatial position data meets the geometric criterion for "complete penetration," and the calculated residual velocity data is greater than zero, it indicates that the penetrator has successfully penetrated the target plate and still possesses forward kinetic energy; the main physical process ends, and the simulation can be terminated.

[0138] For example, in a simulation, the target impact surface is located in the Z=0 mm plane, and the back plate surface is located in the Z=4.1 mm plane. When the Z coordinates of all units or particles representing the projectile and the adhering target plug are greater than 4.1 mm, it is determined as "complete penetration". At the same time, the average velocity of these units and particles is calculated, and the remaining velocity data is 150 meters per second (greater than zero). At this point, the condition is met, and the simulation terminates.

[0139] 2. Condition 2: The total kinetic energy of the system is lower than the first preset threshold.

[0140] This condition applies to scenarios where the projectile fails to penetrate the target plate, and its kinetic energy is completely absorbed by processes such as the plastic deformation of the target plate. At each time step, the total kinetic energy of the entire model (including the projectile and the target plate) is calculated, which is the sum of the kinetic energies of all finite element nodes and SPH particles. The first preset threshold is a very small positive value, set according to the initial total kinetic energy of the model; for example, it could be one millionth of the initial kinetic energy or an absolute value (such as 1.0e-6 joules). After the impact process ends, there is no more macroscopic motion, and the total kinetic energy approaches zero. When the calculated total kinetic energy data remains consistently below this threshold, it indicates that the dynamic process of the system has substantially stopped, and the simulation can be terminated.

[0141] The first preset threshold is determined through a calibration process combining pre-calibrated numerical simulations with engineering experience. Specifically, when the simulated system reaches a physically stationary state, its total kinetic energy will tend towards a theoretically near-zero minimum value. This threshold is set to capture this stable state, and its specific value is usually determined proportionally based on the initial total kinetic energy data of the projectile-target system. For example, it may be set between one ten-thousandth and one millionth of the initial kinetic energy, or an absolute value determined according to the model size (e.g., 1 × 10⁻⁶). -6 Joules up to 1×10 - (³ Joules); those skilled in the art can determine a specific threshold value by running a simplified model or using empirical formulas to estimate the magnitude when the system's kinetic energy is completely dissipated. When the total kinetic energy of the system calculated by simulation is consistently lower than this threshold, it can be reasonably determined that the macroscopic motion of the system has stopped, thus meeting the conditions for terminating the calculation.

[0142] For example, the initial total kinetic energy in the simulation is approximately 1500 joules. As the impact progresses, the kinetic energy is continuously converted into internal energy (plastic heat). Later in the simulation, the calculated total kinetic energy is 0.001 joules. The preset first threshold is 0.01 joules. Since 0.001 joules < 0.01 joules, the condition is met, and the simulation terminates.

[0143] 3. Condition 3: The rate of change of the system's total energy over time is lower than the second preset threshold.

[0144] This condition determines whether the simulation has reached a steady state from the perspective of energy conservation and stability. The total energy of the system includes all forms of energy, such as kinetic energy and internal energy (reflecting plastic work and elastic strain energy). The rate of change of total energy over time is monitored, i.e., the change in total energy per unit simulation time. The second preset threshold is a very small value used to characterize the "stagnation" level of energy change. In the early stages of the simulation, energy conversion is rapid, and the rate of change is large. When the simulation progresses to the point where the physical processes are essentially complete, the total energy tends to be constant, and its rate of change approaches zero. When this rate of change remains below the second preset threshold, it indicates that the system has reached a stable state of energy balance, and continuing the calculation will no longer change the physical results; the simulation can then be terminated.

[0145] The determination of the "second preset threshold" is accomplished by analyzing the energy balance state of the system at the end of the simulation and combining it with the convergence requirements of the numerical calculation. Specifically, this threshold is a minimal rate of change that characterizes the near-stagnant change in the total energy of the system (including kinetic energy and internal energy). Its setting is based on the following process: during the system relaxation phase after the completion of the main impact physical processes (such as projectile penetration cessation or penetration), by monitoring the change in total energy over a simulation period and calculating its average rate of change, and according to the required calculation accuracy (e.g., the allowable final energy fluctuation range), the magnitude of the rate of change reflecting that the energy curve has entered a stable plateau region is determined as the second preset threshold. This threshold can typically be set to 10 times the initial total energy of the system per second. -6 Up to 10 -4 The magnitude, or an absolute value based on the model size (e.g., 1 × 10⁻⁶ per second). - (3 joules to 1 joule) By running a simplified pre-simulation to observe the energy stabilization process, those skilled in the art can calibrate the proportional coefficient or absolute value applicable to a specific model, thereby ensuring that when the rate of change of the total energy of the system over time is consistently below this threshold, it can be determined that the system has reached a stable state of energy balance.

[0146] For example, the total energy of the system is monitored over multiple consecutive time steps, and its rate of change is calculated. In the initial stage of the impact, the rate of change may be as high as 1e9 joules per second. Towards the end of the simulation, the average energy change rate over a recent period is calculated to be 1.0e-3 joules per second. The preset second threshold is 1.0 joules per second. Since 1.0e-3 joules / second < 1.0 joules / second, the condition is met, and the simulation terminates.

[0147] 4. Condition 4: The total simulated physical time data reaches the preset maximum duration.

[0148] This condition is a safety precaution to prevent the simulation from looping indefinitely if the aforementioned physical conditions are not triggered. The total simulation physical time data is accumulated from the start of the simulation. Users preset a maximum duration sufficient to cover all possible physical processes (e.g., 500 microseconds after the impact) based on experience or estimation. Regardless of whether the impact penetrates or the energy stabilizes, the simulation is forcibly terminated once the accumulated simulation time reaches this preset maximum duration. This ensures the controllability of the computational task.

[0149] The determination of the "preset maximum duration" is based on an engineering estimate of the time required for the complete physical process of the projectile impact event, and on this basis, a sufficient safety boundary is set for the numerical simulation calculation. Specifically, this duration is determined by analyzing empirical formulas, experimental data under similar conditions (such as the total time of impact events recorded by high-speed photography), or preliminary simplified simulation results. The total physical time from the start of the impact to the system reaching stability or the projectile moving until it no longer has a significant impact on the target plate is first estimated, and then this time is multiplied by a safety factor (usually 2 to 10 times) to finally determine the duration. For example, for the case of a 7.62mm bullet impacting a 4.1mm steel plate, the main impact process observed in experiments is usually completed within hundreds of microseconds, so the preset maximum duration can be set to the order of several milliseconds. Those skilled in the art can determine a specific duration value (such as 5000 microseconds) for a specific simulation scenario using the above method to ensure that all possible dynamic processes are covered, while avoiding infinite simulation loops due to individual cases.

[0150] The user-preset maximum duration is 1000 microseconds. After the simulation starts, the internal clock continuously accumulates the simulation time. When the accumulated simulation physical time data reaches 1000 microseconds, the simulation will automatically terminate even if the projectile has not yet completely stopped or penetrated.

[0151] Please see Figure 3 As shown, in one embodiment, a projectile impact simulation device is provided, the device comprising: The reading module 301 is used to read the input projectile geometric model data, the material parameter data of the projectile and the target plate, and the impact initial condition data including the initial velocity of the projectile. The generation module 302 is used to generate three-dimensional finite element mesh model data of the projectile based on the projectile geometric model data; The assignment module 303 is used to assign corresponding material property data to the elements in the three-dimensional finite element mesh model data based on the projectile material parameter data. The generation module 302 is also used to assign initial motion state data to the mesh nodes of the three-dimensional finite element mesh model data based on the impact initial condition data, thereby generating initialized projectile finite element model data. The generation module 302 is also used to generate hybrid discrete model data of the target plate based on the input target plate geometric model data and material parameter data; wherein, the hybrid discrete model data consists of SPH particle model data of the core region and finite element mesh model data of the outer region, and a consolidation connection relationship data is established at the interface between the two. Iteration module 304 is used to start a time-step iterative solution loop using the initialized projectile finite element model data and the target plate hybrid discrete model data as initial system state data; The output module 305 is used to exit the time-step iterative solution loop when the simulation termination condition is met, and generate and output target plate damage morphology data, projectile penetration result data and ballistic limit velocity based on the final system state data.

[0152] For other details regarding the implementation of the above technical solutions by each module in the above-mentioned target impact simulation device, please refer to the description in the target impact simulation method provided in the above-mentioned invention embodiments, which will not be repeated here.

[0153] In one embodiment, a computer device is provided, the internal structure of which can be shown as follows: Figure 4 As shown, the computer device includes a processor, memory, network interface, display screen, and input devices connected via a system bus. The processor provides computational and control capabilities. The memory includes non-volatile storage media and internal memory. The non-volatile storage media stores the operating system and computer programs. The internal memory provides an environment for the operation of the operating system and computer programs in the non-volatile storage media. The network interface is used to communicate with an external server via a network connection. When the computer program is executed by the processor, it implements the functions or steps of a projectile impact simulation method.

[0154] In one embodiment, a computer device is provided, including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor, when executing the computer program, implements as follows: Figure 2 The steps are shown.

[0155] It should be noted that the functions or steps that can be implemented by the computer-readable storage medium or computer device described above can be referred to the relevant descriptions on the server side and client side in the foregoing method embodiments. To avoid repetition, they will not be described one by one here.

[0156] Those skilled in the art will understand that all or part of the processes in the methods of the above embodiments can be implemented by a computer program instructing related hardware. The computer program can be stored in a non-volatile computer-readable storage medium. When executed, the computer program can include the processes of the embodiments of the above methods. Any references to memory, storage, databases, or other media used in the embodiments provided in this application can include non-volatile and / or volatile memory. Non-volatile memory may include read-only memory (ROM), programmable ROM (PROM), electrically programmable ROM (EPROM), electrically erasable programmable ROM (EEPROM), or flash memory. Volatile memory may include random access memory (RAM) or external cache memory. By way of illustration and not limitation, RAM is available in a variety of forms, such as static RAM (SRAM), dynamic RAM (DRAM), synchronous DRAM (SDRAM), dual data rate SDRAM (DDRSDRAM), enhanced SDRAM (ESDRAM), synchronous link DRAM (SLDRAM), RAMbus direct RAM (RDRAM), direct memory bus dynamic RAM (DRDRAM), and memory bus dynamic RAM (RDRAM), etc.

[0157] Those skilled in the art will clearly understand that, for the sake of convenience and brevity, the above-described division of functional units and modules is used as an example. In practical applications, the above functions can be assigned to different functional units and modules as needed, that is, the internal structure of the device can be divided into different functional units or modules to complete all or part of the functions described above.

[0158] The above-described embodiments are only used to illustrate the technical solutions of the present invention, and are not intended to limit it. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention, and should all be included within the protection scope of the present invention.

Claims

1. A method for simulating projectile impact, characterized in that, The method includes the following steps: Read the input projectile geometry model data, the material parameter data of the projectile and the target plate, and the impact initial condition data including the initial velocity of the projectile; Based on the projectile geometric model data, generate three-dimensional finite element mesh model data of the projectile; Based on the projectile material parameter data, the corresponding material property data are assigned to the elements in the three-dimensional finite element mesh model data; Based on the initial impact condition data, the mesh nodes of the three-dimensional finite element mesh model data are assigned initial motion state data, thereby generating initialized projectile finite element model data. Based on the input target plate geometric model data and material parameter data, a hybrid discrete model data of the target plate is generated; wherein, the hybrid discrete model data consists of SPH particle model data of the core region and finite element mesh model data of the outer region, and a consolidation connection relationship data is established at the interface between the two. Using the initialized finite element model data of the projectile and the hybrid discrete model data of the target plate as the initial system state data, the time-step iterative solution loop is started; When the simulation termination condition is met, the time-step iterative solution loop is exited, and target plate damage morphology data, projectile penetration result data, and ballistic limit velocity are generated and output based on the final system state data.

2. The projectile impact simulation method according to claim 1, characterized in that, The time-step iterative solution loop includes the following sub-steps: Based on the current spatial position data of the projectile and the SPH particles of the target plate, a contact search is performed and contact pair data is generated, and contact force data is calculated and output based on the contact pair data. The contact force data, along with the current stress, strain, and damage variable data of the elements and particles, are input into the material constitutive model for calculation, and the updated stress data, strain data, damage variable data, and temperature data are output. The distortion measurement data is calculated based on the updated strain data of the elements in the finite element model of the projectile. When the distortion measurement data exceeds a preset threshold, element conversion instruction data is generated. In response to the unit conversion command data, the current mass, momentum, energy and stress state data of the finite element element to be converted are mapped to the newly generated SPH particle set to form new projectile SPH particle data, and this part of the data is removed from the projectile finite element model data. Based on the updated damage variable data and stress state data of SPH particles on the target plate, the failure particle data is identified and marked; wherein, the failure particle data stops participating in the constitutive calculation but continues to retain its mass and momentum data; Based on the current total contact force data and internal stress data, the system motion equations are solved by integration, and the displacement and velocity data of all calculated nodes in the next time step are output, and the system state data is updated.

3. The projectile impact simulation method according to claim 1, characterized in that, The process of generating a three-dimensional finite element mesh model of the projectile based on the projectile's geometric model data includes: Analyze the projectile geometric model data to obtain specific dimensional parameters characterizing the layered structure of the projectile jacket, lead sleeve, and steel core. Based on the specific size parameters, a three-dimensional projectile solid model with the layered structure is generated through geometric reconstruction, and different material layer regions are identified. The three-dimensional projectile solid model is meshed using a hexahedral finite element method to generate node coordinate data and element connection relationship data; Based on the spatial position of the unit in the three-dimensional projectile solid model, establish the association data between the unit and the material layer region; Based on the node coordinate data, the unit connection relationship data, and the association data between the unit and the material layer region, the three-dimensional finite element mesh model data of the projectile is generated and output.

4. The projectile impact simulation method according to claim 1, characterized in that, The hybrid discrete model data of the target plate, generated based on the input target plate geometric model data and material parameter data, includes: Based on the target plate geometric model data and the expected impact point position, define the core region data and the peripheral region data surrounding the core region; Based on the core region data and the preset particle spacing data, SPH particle set basic data is generated within the core region; wherein, the SPH particle set basic data includes the initial coordinates and mass data of the particles. Based on the peripheral region data, the core region boundary data, and the interface unit size data that matches the particle spacing data, FEM mesh basic data is generated, wherein the unit size of the FEM mesh gradually increases from the interface with the core region towards the periphery. Based on the particle data of the adjacent interface in the SPH particle set basic data and the cell node data of the interface in the FEM mesh basic data, SPH-FEM consolidation connection relationship data is established. The target plate material parameter data is assigned to all particles in the SPH particle set base data and all cells in the FEM mesh base data, and their stress, strain and damage variable data are initialized to form particle data and cell data with material properties. Fixed constraints are applied to the nodes located at the geometric boundary of the target plate in the FEM mesh base data, and the hybrid discrete model data of the target plate is output based on the particle data with material properties, the element data and the consolidation connection relationship data.

5. The projectile impact simulation method according to claim 2, characterized in that, The process of performing contact search and generating contact pair data based on the current spatial position data of the projectile and the target particles, and calculating and outputting contact force data based on the contact pair data, includes: Based on the position data of the finite element nodes or SPH particles of the projectile and the position data of the SPH particles of the target plate, a global contact search is performed to identify all potential contact object pairs. For each potential contact object pair, perform local precision detection. If the projectile object is a finite element, calculate the penetration depth and contact normal of the target plate particle to the surface of the element. If the projectile object is an SPH particle, calculate the center distance between the two particles and generate contact pair data including contact type and contact geometry accordingly. Based on the contact pair data and the constitutive parameter data of the associated projectile and target plate materials, parallel calculations are performed for different contact types: for the contact between the projectile finite element element and the target plate SPH particles, the contact force is calculated using the penalty function method based on the penetration depth and material stiffness; for the contact between the projectile SPH particles and the target plate SPH particles, the contact force is calculated using the SPH particle interaction model based on the particle spacing, kernel function derivative, and constitutive stress, thereby generating contact force metadata corresponding to each contact pair. The contact force metadata is allocated and integrated: for finite element-particle contact, the contact force is allocated to the corresponding finite element nodes according to the shape function, and the reaction force is allocated to the target plate SPH particles; for particle-particle contact, the interaction force is directly allocated to the corresponding paired SPH particles, and finally the contact force data acting on all relevant nodes and particles of the system is generated.

6. The projectile impact simulation method according to claim 2, characterized in that, The process of inputting the contact force data and the current stress, strain, and damage variable data of the elements and particles into the material constitutive model for calculation, and outputting updated stress data, strain data, damage variable data, and temperature data, includes: Based on the contact force data and the current stress data, the strain increment data and equivalent plastic strain rate data of the unit and the particle are calculated. Based on the strain increment data, the current cumulative plastic strain data and the current temperature data, new cumulative plastic strain data is obtained, and new temperature data is calculated based on the adiabatic model of plastic work-heat conversion. Based on the new cumulative plastic strain data, the current stress triaxiality data, the equivalent plastic strain rate data, and the new temperature data, new damage variable data are calculated using a damage evolution model that includes damage threshold and fracture strain criterion. Based on the new damage variable data, the new cumulative plastic strain data, the equivalent plastic strain rate data, and the new temperature data, the current yield stress data considering damage softening is calculated by the Johnson-Cook model coupled with the damage variables. Based on the current yield stress data and the strain increment data, the stress tensor is updated using the elastoplastic stress integral algorithm, and the updated stress data is output. Simultaneously, the new cumulative plastic strain data, new damage variable data, and new temperature data are output as the updated strain data, damage variable data, and temperature data.

7. The projectile impact simulation method according to claim 2, characterized in that, The calculation of distortion measurement data based on the updated strain data of the elements in the finite element model of the projectile includes: Extract the updated strain data of each element in the finite element model of the projectile, wherein the strain data includes at least the equivalent plastic strain data of the element; Based on the deformed node coordinate data of each element, calculate the geometric parameters that characterize the shape distortion of the element. Using the equivalent plastic strain data and / or the geometric parameter data as input, distortion measurement data that comprehensively characterizes the degree of deformation and distortion of each unit is calculated.

Images