Bionic thin-wall structure crashworthiness analysis method based on finite element simulation

Abstract

The invention discloses a bionic thin-walled structure crashworthiness analysis method based on finite element simulation, particularly relates to the technical field of computer-aided engineering simulation, and is used for solving the problem that the crashworthiness of a complex bionic thin-walled structure is difficult to systematically evaluate in the early stage of design in the prior art. The method comprises the following steps: establishing a three-dimensional geometric model, carrying out finite element grid division, constructing a finite element model for collision simulation, then applying a dynamic impact load to execute collision calculation, extracting node coordinate data, and converting a structure dynamic deformation process into topological characteristic evolution through persistent coherence analysis; a distance-based filtering process is adopted to generate a filtering complex sequence, a persistence graph is obtained by calculating homology group evolution, quantitative indexes, including an average life cycle and persistence entropy, representing topological feature persistence are extracted, then a performance balance graph is constructed, Pareto leading edge distribution is recognized through a non-dominated sorting algorithm in multi-objective optimization, and the performance balance graph is obtained. And the scheme screening of the collision resistance of the bionic thin-wall structure is realized.

Description

Technical Field

[0001] This invention relates to the field of computer-aided engineering simulation technology, and more specifically, to a method for analyzing the impact resistance of biomimetic thin-walled structures based on finite element simulation. Background Technology

[0002] In the field of structural engineering design, especially in the development of thin-walled components such as automotive crash beams and energy-absorbing boxes, biomimetic designs are often inspired by the superior structural features of organisms in nature, such as bamboo and bone, to balance lightweighting and safety. Currently, the traditional methods for assessing the crashworthiness of such biomimetic thin-walled structures mainly rely on the designer's experience and comparative analysis of a few simplified cross-sectional configurations, and require subsequent verification through physical bench tests.

[0003] The aforementioned existing technical methods are insufficient to effectively predict and compare the systematic crashworthiness of structural components with complex biomimetic morphology under collision loads in the early stages of design. This lack of early assessment capability leads to considerable blindness in the design process, prolongs the trial production and verification cycle, and ultimately restricts the selection and application of efficient and reliable biomimetic structural solutions. Summary of the Invention

[0004] In order to overcome the above-mentioned defects of the prior art, the present invention provides a biomimetic thin-walled structure impact resistance analysis method based on finite element simulation to solve the problems mentioned in the background art.

[0005] To achieve the above objectives, the present invention provides the following technical solution: A method for analyzing the impact resistance of biomimetic thin-walled structures based on finite element simulation includes: S1. Establish a three-dimensional geometric model of the biomimetic thin-walled structure to be evaluated; S2. Perform finite element mesh generation on the three-dimensional geometric model and assign material properties and boundary conditions to construct a finite element model for collision simulation. S3. Apply dynamic impact loads to the finite element model and perform collision simulation calculations; S4. Extract the deformation mode sequence of the biomimetic thin-walled structure from the results of the collision simulation calculation. The deformation mode sequence contains the node coordinate data of multiple time steps arranged in chronological order. S5. Based on node coordinate data, perform persistent cohomology analysis on the geometric configuration of the biomimetic thin-walled structure during the collision process, and extract at least one quantitative index to characterize the persistence of topological features. S6. Construct a performance equilibrium diagram based on at least one quantitative index to characterize the constraint relationship between the quantitative indices, and identify the Pareto front distribution in the performance equilibrium diagram to screen out the biomimetic thin-walled structure scheme located at the Pareto front.

[0006] Furthermore, a three-dimensional geometric model of the biomimetic thin-walled structure to be evaluated is established, including: Based on predefined biomimetic structural feature parameters, a three-dimensional geometric model of a biomimetic thin-walled structure with at least one of the following features—bamboo-like cavity partitioning, skeletal gradient porous structure, or honeycomb rib structure—is generated using a parametric modeling method. The parametric modeling method controls the size and spatial distribution of the geometric configuration by adjusting the feature parameters.

[0007] Furthermore, the three-dimensional geometric model is meshed using finite element methods and assigned material properties and boundary conditions to construct a finite element model for collision simulation, including: The three-dimensional geometric model is divided using a hybrid mesh of tetrahedrons and hexahedrons, and local mesh refinement is implemented in characteristic regions with bamboo-like cavity partitions, skeletal gradient porous or honeycomb rib structures. The material distribution properties based on the three-dimensional geometric model are assigned corresponding elastoplastic material parameters, and fixed constraint boundary conditions are applied at the end of the finite element model to simulate the actual installation state.

[0008] Furthermore, dynamic impact loads are applied to the finite element model to perform collision simulation calculations, including: A dynamic impact load along the axial direction is applied to one end of the finite element model. The dynamic impact load is achieved by giving the impact mass block an initial velocity or by using a rigid wall displacement loading method. An explicit dynamic algorithm is used to solve the transient response of the finite element model during the collision process, where the impact velocity is set to a constant value that conforms to the actual working conditions.

[0009] Furthermore, the explicit dynamic algorithm is used to solve the transient response of the finite element model during the collision process, which includes: solving the equation of motion using the explicit time integration method, calculating the nodal displacement, velocity and acceleration at each time step, and considering the nonlinear constitutive relation of the material and the contact interaction.

[0010] Furthermore, from the results of the collision simulation calculations, a deformation mode sequence of the biomimetic thin-walled structure is extracted. This deformation mode sequence contains nodal coordinate data for multiple time steps arranged in chronological order, including: Multiple time steps are selected from the collision simulation results at equal time intervals. The spatial coordinates of all nodes of the biomimetic thin-walled structure at each time step are read using the post-processing function of the finite element analysis software. The node coordinate data arranged in chronological order are integrated into a deformation mode sequence, where the node coordinate data is stored in the form of a three-dimensional coordinate matrix.

[0011] Furthermore, based on nodal coordinate data, a persistent cohomology analysis is performed on the geometric configuration of the biomimetic thin-walled structure during the collision process, extracting at least one quantitative index to characterize the persistence of topological features, including: Point cloud data representing geometric configurations is constructed based on node coordinate data; A distance-based filtering process is applied to point cloud data to generate a filtered complex sequence; Persistent graphs are obtained by calculating the homology group evolution of filtered complex sequences; Extract the lifetime distribution of topological features from the persistent graph, and calculate the average lifetime of the lifetime distribution as the first quantification index, while calculating the persistence entropy of the lifetime distribution as the second quantification index.

[0012] Furthermore, the process of generating a filtered complex sequence by applying a distance-based filtering process to point cloud data includes: defining a distance filtering function, progressively increasing the distance threshold, and constructing a Vietoris-Rips complex sequence.

[0013] Furthermore, obtaining a persistent graph by calculating the homology group evolution of the filtered complex sequence includes: calculating the homology group of each filtered complex, tracking the birth and death of topological features, and generating a persistent graph.

[0014] Furthermore, a performance equilibrium diagram is constructed based on at least one quantitative index to characterize the constraints between the quantitative indices. The Pareto front distribution in the performance equilibrium diagram is identified to screen out biomimetic thin-walled structure schemes located at the Pareto front, including: A performance equilibrium diagram was constructed using the first and second quantitative indicators as two-dimensional coordinate axes, and the quantitative indicator values ​​corresponding to multiple biomimetic thin-walled structure schemes were plotted as a scatter distribution. The Pareto front distribution in the performance equilibrium graph is identified by the non-dominated sorting algorithm in multi-objective optimization. The Pareto front distribution is composed of Pareto optimal solutions that are not dominated by other solutions. The structural scheme corresponding to the Pareto optimal solution is selected from all biomimetic thin-walled structural schemes as the final preferred scheme.

[0015] Compared with the prior art, the present invention has the following beneficial effects: 1. By integrating finite element simulation and topology data analysis techniques, a systematic assessment of the impact resistance performance of biomimetic thin-walled structures is achieved in the early design stage. Based on parametric modeling, a three-dimensional geometric model with bio-inspired features is generated, and dynamic collision simulation is combined to accurately capture the transient response of the structure under impact loads. By extracting the node coordinate sequence during the deformation process and applying persistent cohomology analysis, the complex geometric evolution is transformed into quantifiable topological characteristic indicators, thereby overcoming the limitations of traditional experience-dependent methods. The computer simulation-based analysis method can reveal the inherent laws of energy absorption and deformation modes of different biomimetic configurations before physical prototyping, effectively reducing design blindness and providing reliable data support for structural optimization.

[0016] 2. By constructing a performance equilibrium graph representing the persistence of topological features and identifying the Pareto front, efficient scheme selection under multi-objective optimization is achieved. Abstract structural performance is transformed into an intuitive two-dimensional spatial distribution, enabling designers to quickly locate the most promising biomimetic thin-walled structure scheme in the crashworthiness trade-off. This significantly shortens the traditional development cycle that relies on iterative trial and error, reduces experimental costs, and ensures a balance between lightweighting and safety in the final design scheme. It provides a decision-making basis for engineering applications that combines theoretical rigor with practical feasibility. Attached Figure Description

[0017] Figure 1 This is a flowchart of the biomimetic thin-walled structure impact resistance analysis method based on finite element simulation according to the present invention. Detailed Implementation

[0018] 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 embodiments of the present invention, and not all embodiments. 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.

[0019] Example: Figure 1 This invention presents a method for analyzing the impact resistance of biomimetic thin-walled structures based on finite element simulation, including: S1. Establish a three-dimensional geometric model of the biomimetic thin-walled structure to be evaluated; S2. Perform finite element mesh generation on the three-dimensional geometric model and assign material properties and boundary conditions to construct a finite element model for collision simulation. S3. Apply dynamic impact loads to the finite element model and perform collision simulation calculations; S4. Extract the deformation mode sequence of the biomimetic thin-walled structure from the results of the collision simulation calculation. The deformation mode sequence contains the node coordinate data of multiple time steps arranged in chronological order. S5. Based on node coordinate data, perform persistent cohomology analysis on the geometric configuration of the biomimetic thin-walled structure during the collision process, and extract at least one quantitative index to characterize the persistence of topological features. S6. Construct a performance equilibrium diagram based on at least one quantitative index to characterize the constraint relationship between the quantitative indices, and identify the Pareto front distribution in the performance equilibrium diagram to screen out the biomimetic thin-walled structure scheme located at the Pareto front.

[0020] S1. Establish a three-dimensional geometric model of the biomimetic thin-walled structure to be evaluated, specifically as follows: In establishing the three-dimensional geometric model of the biomimetic thin-walled structure to be evaluated, the model is first generated based on predefined biomimetic structural feature parameters. These parameters include specific values ​​such as cavity diameter, spacing, porosity gradient change rate, and rib thickness. The parameters are defined through physical measurements and data analysis of the macroscopic structures of natural organisms such as bamboo, bone, or honeycomb. For example, a 3D scanner is used to acquire surface point cloud data of biological samples. The point cloud data is then denoised and meshed using point cloud processing software. A geometric feature extraction algorithm is then applied to calculate key dimensions. The cavity diameter is obtained by identifying the boundary points of the cavity regions in the point cloud and calculating the average distance, with a range, for example, between 5 mm and 20 mm. The spacing is determined by measuring the minimum distance between adjacent cavities. The range of parameters is, for example, between 3 mm and 15 mm. The porosity gradient change rate is determined by analyzing the trend of the change in pore size along the structural axis and fitting the change curve using a linear regression method. The change rate value is set to, for example, 0.1 to 0.5 mm per pore. The rib thickness is determined by extracting the width information of the rib structure and taking the average value, with a range of, for example, between 1 mm and 10 mm. These parameters are set according to actual engineering requirements such as lightweighting goals and impact resistance requirements. The lightweighting goal is determined by structural mass calculation, and the impact resistance requirement is determined by setting a safety factor by referring to historical collision test data. The specific values ​​of the parameters are adjusted by iterative optimization methods. The iterative process includes setting initial parameter values ​​and performing preliminary finite element analysis, and modifying the parameters according to the analysis results, such as stress distribution and deformation mode, until the performance indicators are met.

[0021] Next, a three-dimensional geometric model of a biomimetic thin-walled structure with at least one of the following characteristics—a bamboo-like cavity partition, a skeletal gradient porous structure, or a honeycomb ribbed structure—is generated using parametric modeling methods. The parametric modeling method is implemented using procedural tools in computer-aided design software, such as using a Python script to call the application programming interface (API) of SolidWorks or CATIA. The script writing includes defining input parameter variables and geometric generation logic. The input parameter variables are bound to predefined biomimetic structural feature parameters, and the geometric generation logic automatically creates geometric elements based on the parameter values. For the bamboo-like cavity partition, the generation process includes first creating a basic cylindrical or rectangular outline, with the outline dimensions set according to engineering requirements, such as a length of 200 mm and a diameter of 50 mm. Then, partition surfaces are added inside the outline based on the cavity diameter and partition spacing parameters. These partition surfaces are generated through extrusion or rotation operations. The arrangement of the partition surfaces includes equidistant or variable spacing. In the equidistant arrangement, the partition spacing is fixed at, for example, 8 mm. In the variable spacing arrangement, the partition spacing varies from 5 mm to a linear function. The thickness is gradually increased to 15 mm. The function expression is implemented through the mathematical module in the script. For skeletal gradient porous structures, the generation process uses the Voronoi algorithm to create a random hole distribution. The Voronoi algorithm is implemented by generating random seed points and calculating Voronoi cells. The number of seed points is set according to the model size, for example, 100 points. The hole size is controlled by the porous gradient change rate parameter. The change rate is calculated by the interpolation function to calculate the hole diameter at different locations, for example, decreasing from 10 mm at one end of the structure to 2 mm at the other end. The interpolation function uses linear or quadratic functions, and the function parameters are set according to the gradient requirements. For honeycomb ribbed structures, the generation process uses a honeycomb grid algorithm to create a hexagonal grid pattern. The grid cell size, for example, the hexagonal side length is between 5 mm and 20 mm. The rib thickness parameter controls the width of the grid lines, for example, set to 3 mm. After the grid is generated, a three-dimensional ribbed structure is formed by stretching. All geometric models are saved in standard formats such as STEP or IGES to ensure compatibility with subsequent finite element software.

[0022] Parametric modeling methods control the size and spatial distribution of geometric configurations by adjusting feature parameters. The adjustment process involves setting parameter-driven relationships within the modeling software; for example, the cavity diameter parameter can be associated with the model's global dimensions. This association is achieved through mathematical expressions. When the cavity diameter is modified, the model automatically updates all relevant geometric elements. Size control is achieved through scaling functions, which calculate the actual dimensions based on input parameters. For instance, a scaling factor can be set between 0.5 and 2.0, determined by structural strength requirements calculated using material yield strength and expected loads. Spatial distribution control is achieved through distribution functions, including uniform, random, or gradient distributions. In uniform distributions, parameter values ​​are constant, while in random distributions, a random number generator is used. Parameter values ​​are defined, and a random seed is set according to the model identifier. Linear or nonlinear interpolation functions are used in the gradient distribution to calculate the parameter values ​​at different locations. For example, for porous structures, a linear interpolation function is used to calculate the hole size along the axis. The number of interpolation points is set to 50 points according to the model complexity to ensure continuous distribution. The basis for parameter adjustment includes preliminary finite element analysis results or experimental data. For example, stress cloud maps are obtained by performing static load simulation. After identifying high-stress areas, parameters are adjusted, such as increasing the thickness of ribs or modifying the cavity distribution. The optimization process is automated by scripts. The scripts read simulation results and modify parameter values, iterating until the stress value is lower than the material yield strength. The final generated three-dimensional geometric model is used for subsequent finite element analysis to ensure that the model conforms to both biomimetic morphology and engineering performance.

[0023] S2. Perform finite element mesh generation on the 3D geometric model and assign material properties and boundary conditions to construct a finite element model for collision simulation. The specific implementation is as follows: In the finite element mesh generation process of the 3D geometric model, a hybrid tetrahedral and hexahedral mesh is first used. The hybrid mesh is generated using mesh generation tools in finite element analysis software, such as the automatic mesh generation function of ANSYS or Abaqus. The meshing process includes importing the 3D geometric model obtained in step S1 and setting the mesh type and size parameters. Tetrahedral meshes are used to handle complex geometric regions such as curved surfaces or hole structures, while hexahedral meshes are used for regular regions to improve computational efficiency. The mesh size parameters are set according to the overall size and feature complexity of the model; for example, the overall mesh size is set to between 1 / 5 and 1 / 10 of the minimum feature size of the model. The minimum feature size of the model is determined by measuring the dimensions of the 3D geometric model. The minimum values ​​of edge and radius of curvature are obtained. For example, when the minimum feature size is calculated to be 2 mm using the measurement tools of CAD software, the mesh size is set between 0.4 mm and 0.2 mm. During the mesh generation process, mesh quality indicators such as element aspect ratio and twist are checked. The element aspect ratio threshold is set to less than 5:1, and the element twist threshold is set to less than 0.7. These thresholds are set by referring to the stability requirements in finite element analysis theory to ensure computational convergence and accuracy. After the mesh is generated, mesh smoothing and optimization are performed using the software's built-in tools. The number of smoothing iterations is set between 10 and 20. The optimization process includes adjusting the node positions to improve the element shape, for example, using the Laplace smoothing algorithm to move nodes to the center position of adjacent elements.

[0024] Local mesh refinement is performed in characteristic regions with bamboo-joint cavity partitions, skeletal gradient porous structures, or honeycomb rib structures. Local mesh refinement is achieved by identifying the geometric boundaries of the characteristic regions. For example, for bamboo-joint cavity partitions, the cavity interface and partition edge are identified; for skeletal gradient porous structures, the pore boundaries and gradient change regions are identified; for honeycomb rib structures, the rib intersections and mesh cell boundaries are identified. The identification method uses a geometric feature detection algorithm based on curvature changes or distance measurements. For example, a curvature threshold of 0.1 mm is set to detect high curvature regions. The curvature threshold is determined by analyzing the histogram of curvature distribution on the model surface. The top 10% of curvature values ​​are designated as high curvature regions. Local refinement parameters include refinement level and refinement size. The refinement level is set to 2 to 4, and the refinement size is set to 1 / 2 to 1 / 4 of the global mesh size. For example, when the global mesh size is 0.5 mm, the local refinement size is 0.25 mm to 0.125 mm. The refinement process is implemented through the local mesh control function in the software. For example, in Abaqus, seed points or bounding boxes are used to define the refinement region. After refinement, mesh continuity is checked to ensure that there are no overlapping or gap elements. The continuity check is completed through the mesh diagnostic tools in the software, such as checking element shared nodes and face connectivity.

[0025] Material distribution properties based on a 3D geometric model are assigned corresponding elastoplastic material parameters. These properties are obtained through region identifiers or layer information within the geometric model, such as assigning material labels to different structural regions in CAD or CAE software. Material parameters include elastic modulus, yield strength, hardening factor, and density. These parameters are set according to the actual material type, such as aluminum alloy or composite materials. For example, the elastic modulus of aluminum alloy is set to 70 GPa, the yield strength to 250 MPa, the hardening factor is defined using a power-law hardening model with a power-law exponent of 0.2, and the density to 2700 kg / m³. In finite element software, parameter assignment is achieved through a material library or user-defined input. For example, in ANSYS, a material model is created and parameter values ​​are input. The parameter values ​​are associated with the geometric region by selecting the corresponding region and applying material properties. For gradient material distributions, such as skeletal gradient porous structures, material parameters are controlled by interpolation functions based on positional variations. For example, the elastic modulus changes linearly from 100 GPa at one end of the structure to 50 GPa at the other end. The interpolation function uses linear equations, and the equation parameters are set according to the gradient distribution requirements. The gradient distribution requirements are determined by structural performance objectives such as stress homogenization.

[0026] Fixed constraint boundary conditions are applied to the ends of the finite element model to simulate the actual installation state. The fixed constraints are achieved by restricting all degrees of freedom at the ends of the model, including translational and rotational degrees of freedom. The constraint positions are determined according to the actual installation method. For example, fixed constraints are applied to the end face for cylindrical structures and to the end plane for rectangular structures. The constraint application process involves selecting the end node or face in the finite element software and applying fixed support conditions. For example, in Abaqus, the Encastre option is used to restrict all degrees of freedom. The constraint range is defined using the geometric selection tool. For example, a circular face with a radius of 5 mm or a rectangular face with a side length of 10 mm is selected for the end region. The constraint is verified by running a preliminary static analysis to check the displacement response, ensuring that the constraints are effective and there is no rigid body motion. After the boundary conditions are set, the finite element model is saved for subsequent collision simulations.

[0027] S3. Apply dynamic impact loads to the finite element model and perform collision simulation calculations. The specific implementation is as follows: During the application of dynamic impact loads to the finite element model, an axial dynamic impact load is first applied to one end of the finite element model. The application of the dynamic impact load is achieved in two ways. The first way is by imparting an initial velocity to the impact mass block. The mass of the impact mass block is set to be between 5 and 20 times the mass of the finite element model. The specific value is calculated based on the principle of energy conservation. The principle of energy conservation requires that the kinetic energy of the impact mass block matches the energy absorbed by the finite element model. For example, when the mass of the finite element model is 2 kg, the mass of the impact mass block is set to be between 10 kg and 40 kg. The initial velocity is set to between 5 m / s and 20 m / s based on the actual working conditions. The first method uses a constant value within a meter range. The velocity value is determined by referring to typical collision scenarios in actual engineering applications, such as the common speed range in car crash tests. The second method uses rigid wall displacement loading. The displacement function of the rigid wall is set as a function that increases linearly with time. The displacement rate is calculated based on the impact velocity equivalence principle. The impact velocity equivalence principle requires that the rigid wall displacement rate and the impact velocity be consistent at the contact point. For example, when the target impact velocity is 10 m / s, the displacement rate should be set to 10 m / s. The contact distance between the rigid wall and the finite element model is set between 1 mm and 10 mm. The contact distance is set by considering the model size and calculation stability.

[0028] An explicit dynamics algorithm is used to solve the transient response of the finite element model during the collision process. The explicit dynamics algorithm uses an explicit time integration method to solve the equations of motion, which include the equilibrium equations where the mass matrix multiplied by the acceleration equals the external force vector minus the internal force vector. The mass matrix is ​​obtained by assembling the element mass matrices, which are in lumped mass matrix form, achieved by diagonalizing the element mass matrix. The external force vectors include impact loads and contact forces. The impact loads are calculated based on the actual collision conditions, and the contact forces are calculated using a contact algorithm. The internal force vectors are calculated using element stress integration, which employs the Gaussian integration method, with 2 to 4 Gaussian integration points set within each element. The explicit time integration method employs the central difference method. The implementation of the central difference method includes calculating the nodal acceleration at each time step. The nodal acceleration is solved through the equation of motion, i.e., the acceleration is equal to the inverse of the mass matrix multiplied by the residual force vector. The residual force vector is equal to the external force vector minus the internal force vector. The inverse of the mass matrix is ​​obtained by inverting the diagonal mass matrix, and the elements of the diagonal mass matrix are the nodal masses. Then, the nodal velocity is updated based on the nodal acceleration. The nodal velocity is equal to the nodal velocity of the previous time step plus the product of the nodal acceleration of the current time step and the time step length. Finally, the nodal displacement is updated based on the nodal velocity. The nodal displacement is equal to the nodal displacement of the previous time step plus the product of the nodal velocity of the current time step and the time step length.

[0029] When calculating the nodal displacement, velocity, and acceleration at each time step, the time step size is automatically adjusted based on the minimum mesh size and the material wave velocity. The time step size satisfies the stability condition, which requires the time step size to be less than the minimum mesh size divided by the material wave velocity. The minimum mesh size is obtained by traversing all element sizes and taking the minimum value. The element size is determined by calculating the minimum value of the element side length. The material wave velocity is calculated based on the material's elastic modulus and density, which are obtained from the material properties. For example, aluminum alloy has an elastic modulus of 70 gigapascals, a density of 2700 kg / m³, and a material wave velocity of approximately 5000 m / s. When the minimum mesh size is 0.5 mm, the time step is set to 1e-7 seconds. The nonlinear constitutive relation of the material is considered in the calculation. The nonlinear constitutive relation of the material is described by an elastoplastic model, which includes an elastic stage and a plastic stage. The stress-strain relationship in the elastic stage follows Hooke's law, and the plastic stage adopts an isotropic hardening model. The yield criterion adopts the von Mises yield criterion. The hardening curve is obtained by fitting experimental data. The experimental data comes from the tensile test of the material. For example, the hardening curve of aluminum alloy is fitted by a power law function. The strength coefficient and hardening exponent in the power law function are determined by fitting the experimental data by the least squares method.

[0030] Considering contact interactions involves defining contact pairs and a contact algorithm. Contact pairs consist of an impact mass or rigid wall and the surface of the finite element model. The contact algorithm employs a penalty function method, which prevents penetration by applying a contact force. The magnitude of the contact force is proportional to the penetration depth, and the proportionality coefficient is called the penalty stiffness. The penalty stiffness is set between 10 and 100 times the element stiffness of the contact region. The element stiffness of the contact region is obtained by calculating the ratio of the element's Young's modulus to its element size. For example, if the element's Young's modulus is 70 gigapascals and the element size is 1 millimeter, the element stiffness is approximately 70 meganewtons per millimeter. The penalty stiffness is set between 700 meganewtons per millimeter and 7000 meganewtons per millimeter. Contact search... A combination of global and local search methods is employed. The global search is performed every 10 time steps, while the local search is performed once per time step. The global search threshold is set to 1.1 times the maximum model size, and the local search threshold is set to 1.5 times the average element size. The friction coefficient is set between 0.1 and 0.3. The friction model is the Coulomb friction model, and the friction coefficient is obtained through material surface property testing. The contact interaction also includes a failure criterion, which is based on equivalent plastic strain. The failure threshold is set between 0.2 and 0.5 and is determined through material failure experiments. When the equivalent plastic strain of an element exceeds the failure threshold, the element is deleted from the model.

[0031] The entire collision simulation calculation process is implemented in finite element analysis software, such as LS-DYNA or Abaqus / Explicit solver. The total calculation time is set according to the model size and impact velocity, for example, between 0.01 seconds and 0.1 seconds. The total calculation time is determined by estimating the collision duration, which is equal to the model length divided by the impact velocity. For example, when the model length is 0.5 meters and the impact velocity is 10 meters per second, the collision duration is approximately 0.05 seconds. During the calculation, the nodal displacement, velocity, and acceleration data for each time step are saved. The data saving frequency is set to once every 0.0001 seconds, which is set according to the analysis accuracy requirements. The saved data is used for subsequent deformation mode sequence extraction.

[0032] S4. From the results of the collision simulation calculation, extract the deformation mode sequence of the biomimetic thin-walled structure. The deformation mode sequence contains nodal coordinate data of multiple time steps arranged in chronological order. Specifically, the implementation is as follows: When extracting the deformation mode sequence of the biomimetic thin-walled structure from the collision simulation results, the first step is to extract data from the collision simulation results completed in step S3. The extraction process is based on the result file generated by the finite element analysis software. The result file contains simulation data for multiple time steps, and each time step stores the physical state information of all nodes of the biomimetic thin-walled structure. The first step in extracting the deformation mode sequence is to determine the time step selection strategy, using an equal time interval selection method. The interval value of the equal time interval is determined based on the total collision simulation time and the required analysis accuracy. The total collision simulation time is obtained from the simulation settings in step S3, for example, 0.1 seconds. The required analysis accuracy is determined by the rate of change of the deformation features, which is determined through preliminary analysis of adjacent... The node displacement difference is calculated by measuring the coordinate change of the same node at different time steps. When the node displacement difference exceeds a set threshold, the time interval needs to be reduced. The set threshold is usually one percent of the model feature size. For example, if the model feature size is 100 mm, the threshold is set to 1 mm. The specific value of the time interval is determined through iterative trial calculation. The iterative trial calculation process includes setting the initial time interval to 0.001 seconds, extracting the node coordinates, and calculating the maximum node displacement difference. If the maximum node displacement difference exceeds the set threshold, the time interval is halved and the extraction is repeated until the maximum node displacement difference is less than the set threshold. The final time interval value is generally between 0.0001 seconds and 0.001 seconds, corresponding to 100 to 1000 time steps.

[0033] The spatial coordinates of all nodes of the biomimetic thin-walled structure at each time step are read using the post-processing function of the finite element analysis software. The reading operation is completed in the post-processing module of the finite element analysis software, such as the Result module of ANSYS software or the Visualization module of Abaqus software. The specific operation process includes opening the result file, selecting the time step sequence, arranging the time step sequence at equal time intervals, selecting all nodes, and exporting the coordinate data. The coordinate data includes the X, Y, and Z coordinate values ​​of each node. The unit of the coordinate values ​​is consistent with the unit system of the finite element model, which is determined when the finite element model is established in step S2. For example, when using the millimeter unit system, the coordinate value is in millimeters. When reading the node coordinates, it is necessary to ensure data integrity. The total number of nodes is consistent with the total number of nodes in the finite element mesh in step S2. For example, if the total number of nodes is 10,000, 10,000 sets of coordinate values ​​need to be read at each time step. The coordinate data is stored in floating-point format. The floating-point precision is selected as single precision or double precision according to the analysis requirements. For example, when double precision is selected, the coordinate values ​​are retained to 6 decimal places.

[0034] The node coordinate data arranged in chronological order is integrated into a deformation pattern sequence. The integration process includes two stages: data sorting and data formatting. Data sorting is performed according to the chronological order of time steps, from the initial time step to the end of the final time step. The node coordinate data of each time step is treated as an independent data block. Data formatting organizes the sorted data blocks into a three-dimensional coordinate matrix. The first dimension of the three-dimensional coordinate matrix represents the time step number, the second dimension represents the node number, and the third dimension represents the coordinate components. For example, a deformation pattern sequence containing 500 time steps and 10,000 nodes can be represented as a 500×100 matrix. A 0×3 three-dimensional matrix is ​​generated, where the element a[i,j,k] represents the k-th coordinate component value of the j-th node at the i-th time step. The value of i ranges from 1 to 500, the value of j ranges from 1 to 10000, and the value of k ranges from 1 to 3, corresponding to the X, Y, and Z coordinates. After the matrix is ​​generated, data verification is required. Data verification includes checking whether the matrix dimensions are consistent with expectations and randomly selecting several nodes to check whether the coordinate values ​​are within a reasonable range. The reasonable range is determined according to the geometric dimensions of the model. For example, when the model size is 200 mm × 200 mm × 200 mm, the coordinate values ​​should be between -100 mm and 100 mm.

[0035] Node coordinate data is stored in the form of a three-dimensional coordinate matrix. The storage format is selected as binary or text format according to subsequent processing requirements. The binary format uses the IEEE floating-point standard for storage, which saves storage space and has fast read and write speeds. The text format uses plain text files separated by commas or spaces, which is convenient for manual viewing and cross-platform use. The storage process includes data writing and index building. Data writing saves the three-dimensional coordinate matrix to the file in row-major order. Row-major order means storing the coordinate values ​​of each node line by line according to the time step order. To improve data access efficiency, the index building records the matrix dimensions and the starting position of each time step data in the file header. The storage file name is associated with the model identifier, for example, using the model number and timestamp as the file name. The stored data is used for persistent cohomology analysis in step S5. The storage format of the three-dimensional coordinate matrix matches the input requirements of the persistent cohomology analysis algorithm. For example, if the persistent cohomology analysis algorithm requires the input data to be double-precision floating-point numbers, the storage format is set to double precision accordingly.

[0036] S5. Based on nodal coordinate data, perform persistent cohomology analysis on the geometric configuration of the biomimetic thin-walled structure during the collision process, and extract at least one quantitative index to characterize the persistence of topological features. Specifically, the implementation is as follows: When constructing point cloud data representing geometric configurations based on node coordinate data, node coordinate data is first extracted from the deformation pattern sequence obtained in step S4. The node coordinate data is stored in the form of a three-dimensional coordinate matrix, which contains the spatial position information of nodes at multiple time steps. The point cloud data is constructed by treating the node coordinates at each time step as a set of points in three-dimensional space. The number of point sets is consistent with the number of time steps. Each point set contains the three-dimensional coordinates of all nodes at that time step. The point cloud data is organized by reconstructing the three-dimensional coordinate matrix into a point list. Each element in the point list corresponds to the coordinates of a node, and the coordinate values ​​are directly derived from the elements in the three-dimensional coordinate matrix. The coordinate system is consistent with the coordinate system of the finite element model. For example, when using a rectangular coordinate system, the point cloud data includes three components: X, Y, and Z coordinates. The unit of the point cloud data is consistent with the unit of the node coordinate data. For example, when using millimeters, the point cloud coordinate values ​​are in millimeters. After the point cloud data is constructed, data verification is required. Data verification includes checking whether the number of points in the point cloud is consistent with the total number of nodes, and checking whether the coordinate values ​​are within a reasonable spatial range. The reasonable spatial range is determined according to the geometric dimensions of the biomimetic thin-walled structure. For example, when the outer contour dimensions of the biomimetic thin-walled structure are 200 mm × 200 mm × 200 mm, the coordinate values ​​should be within the range of -100 mm to 100 mm.

[0037] When applying distance-based filtering to point cloud data to generate a filtered complex sequence, a distance filtering function is first defined. This function uses the Euclidean distance function, which calculates the straight-line distance between any two points in the point cloud. The input to the distance filtering function is the coordinates of point pairs in the point cloud data, and the output is the distance value between point pairs. The filtering process is achieved by progressively increasing a distance threshold, which linearly increases from zero to the point cloud diameter. The point cloud diameter is obtained by calculating the maximum distance between any two points in the point cloud. The increment step of the distance threshold is set according to the point cloud density, which is determined by calculating the average nearest neighbor distance of points in the point cloud. For example, if the point cloud contains 10,000 points... If the average nearest neighbor distance of a point is 0.5 mm, then the step size for increasing the distance threshold is set between 0.1 mm and 0.5 mm. The value of the distance threshold ranges from 0 mm to the diameter of the point cloud, for example, a point cloud diameter of 200 mm. A Vietoris-Rips complex sequence is constructed under each distance threshold. The Vietoris-Rips complex is constructed by adding an edge between two points when the distance between them is less than or equal to the current distance threshold, adding a triangle when the distance between any two of three points is less than or equal to the current distance threshold, and so on to construct higher-order simplexes. The filtered complex sequence consists of a series of complexes that gradually expand as the distance threshold increases.

[0038] When obtaining a persistent graph by calculating the homology group evolution of the filtered complex sequence, the homology group of each filtered complex is first calculated. The homology group calculation adopts a linear algebra method, which is implemented by constructing a boundary matrix and calculating the rank of the matrix. The rows of the boundary matrix correspond to the low-dimensional simplex, and the columns correspond to the high-dimensional simplex. The matrix elements represent the boundary relationship between simplexes. The calculation of the homology group includes calculating the quotient group of the kernel space and image space of the boundary matrix. The specific calculation process uses a matrix reduction algorithm, which transforms the boundary matrix into Smith canonical form. The number of non-zero diagonal elements in the Smith canonical form corresponds to the Betti number. The birth and death of topological features are tracked by comparing the differences in the homology groups of adjacent filtered complexes. When a new topological feature appears at a certain distance threshold, it is recorded as a birth. When a topological feature disappears, it is recorded as a death. The distance thresholds for birth and death constitute the life cycle of the topological feature. When generating the persistent graph, each topological feature is represented as a point on a two-dimensional plane. The x-coordinate of the point represents the birth distance threshold, and the y-coordinate of the point represents the death distance threshold. The coordinate range of the persistent graph is consistent with the range of the distance thresholds.

[0039] When extracting the lifecycle distribution of topological features from a persistent graph, the lifecycle distribution is obtained by collecting the differences between the horizontal and vertical coordinates of all points in the persistent graph. The difference between the horizontal and vertical coordinates is the death distance threshold minus the birth distance threshold. The lifecycle distribution is organized in a list structure, where each element corresponds to the lifecycle value of a topological feature. When calculating the average lifecycle of the lifecycle distribution as the first quantification indicator, the average lifecycle is obtained by summing all lifecycle values ​​and dividing by the number of topological features, which is the total number of points in the persistent graph. When calculating the persistence entropy of the lifecycle distribution as the second quantification indicator, the persistence entropy is calculated by multiplying the negative logarithm of each lifecycle value divided by the total lifecycle by the sum of the probabilities of each lifecycle value divided by the total lifecycle. The total lifecycle is the sum of all lifecycle values. The persistence entropy ranges from zero to positive infinity. The larger the persistence entropy value, the more dispersed the persistence of the topological features.

[0040] The specific implementation of generating a filtered complex sequence by applying distance-based filtering to point cloud data includes defining a distance filtering function, which adopts the Euclidean distance formula in the Minkowski distance function. The Euclidean distance formula calculates the distance between two points as the square root of the sum of the squares of the differences of each coordinate component. The distance threshold is gradually increased through iterative iteration. The initial distance threshold of the iteration is set to zero, and the final distance threshold is set to the point cloud diameter. The iteration step size is set according to the required analysis accuracy. For example, when the point cloud diameter is 200 mm, the iteration step size is set to 1 mm. When constructing the Vietoris-Rips complex sequence, a corresponding complex is generated at each iteration step size. The complex is stored using an abstract simplex complex data structure, which records the dimension of each simplex and the index of the points it contains.

[0041] The specific implementation of obtaining a persistent graph by calculating the homology group evolution of the filtered complex sequence includes calculating the homology group of each filtered complex. The homology group calculation uses standard calculation methods in persistent homology algorithm libraries, such as persistent homology calculation functions in JavaPlex or Dionysus libraries. Tracking the birth and death of topological features is achieved by maintaining a feature tracking table. The feature tracking table records the identifier of each topological feature and its birth and death distance thresholds. When generating the persistent graph, the records in the feature tracking table are converted into a graphical representation. The graphical representation is in the form of a scatter plot, with the horizontal axis of the scatter plot being the birth distance threshold and the vertical axis being the death distance threshold. The data structure of the persistent graph is a list of point coordinates.

[0042] S6. Construct a performance equilibrium diagram based on at least one quantitative index to characterize the constraints between the quantitative indexes, identify the Pareto front distribution in the performance equilibrium diagram to screen out biomimetic thin-walled structure schemes located at the Pareto front, specifically as follows: When constructing a performance equilibrium diagram based on the first and second quantification indicators obtained in step S5, the first quantification indicator is the average lifetime, and the second quantification indicator is the persistent entropy. The performance equilibrium diagram is constructed by using the first quantification indicator as the horizontal axis and the second quantification indicator as the vertical axis. The range of the axes is determined based on the quantification indicator values ​​of all biomimetic thin-walled structure schemes. Specifically, the minimum and maximum values ​​of the first and second quantification indicators are calculated for all schemes. For example, if the first quantification indicator value is distributed between 0.1 and 0.9, and the second quantification indicator value is distributed between 1.5 and 3.5, then the horizontal axis range is set to 0 to 1, and the vertical axis range is set to 1 to 4. The axis scale interval is set according to the data precision. The data precision is calculated by the distribution density of the quantification indicator values, and the distribution density is obtained by statistically analyzing the frequency of the quantification indicator values ​​in different intervals. For example, when the first quantification indicator value is evenly distributed in the range of 0.1 to 0.9, the scale interval is set to 0.1, and when the second quantification indicator value is evenly distributed in the range of 1.5 to 3.5, the scale interval is set to 0.5. Multiple biomimetic thin-walled structures are then considered. When the quantitative index values ​​corresponding to the schemes are plotted as a scatter distribution, each biomimetic thin-walled structure scheme corresponds to a scatter point. The x-axis of the scatter point is the first quantitative index value of the scheme, and the y-axis of the scatter point is the second quantitative index value of the scheme. The scatter distribution is implemented using data visualization tools, such as the Matplotlib library in the Python programming language or the plotting function of MATLAB software. The specific implementation process is to create a two-dimensional coordinate system, traverse all schemes, and plot the quantitative index value pairs of each scheme as coordinate points in the coordinate system. The color and shape of the scatter points are distinguished according to the scheme type. Scheme types include bamboo-joint cavity partition, skeletal gradient porous and honeycomb rib structure. For example, the bamboo-joint cavity partition scheme is represented by red circles, the skeletal gradient porous scheme is represented by blue squares, and the honeycomb rib structure scheme is represented by green triangles. The scatter distribution plot includes a legend and coordinate axis labels. The legend explains the scheme type corresponding to different colors and shapes, and the coordinate axis labels indicate the physical meaning of the first and second quantitative indices. The label for the first quantitative index is average lifetime, and the label for the second quantitative index is persistent entropy.

[0043] When identifying the Pareto front distribution in the performance equilibrium graph using the non-dominated ranking algorithm in multi-objective optimization, the input of the non-dominated ranking algorithm is a solution set consisting of the quantized index values ​​of all biomimetic thin-walled structure schemes. Each solution in the solution set corresponds to a pair of quantized index values ​​for a scheme. The core steps of the non-dominated ranking algorithm include calculating the dominance relationship of each solution. The dominance relationship is defined as follows: if a solution is no worse than another solution in all objective functions and is strictly better than another solution in at least one objective function, then the solution dominates the other solution. The objective function is the index to be optimized. In this scenario, the objective function is to maximize the first quantized index and minimize the second quantized index. The specific calculation process involves traversing all solution pairs in the solution set and checking whether each solution is dominated by other solutions. The algorithm records the number of times each solution is dominated and the set of dominated solutions. Then, it performs hierarchical sorting. The first layer is the set of non-dominated solutions, i.e., Pareto optimal solutions, which are dominated zero times. The second layer consists of non-dominated solutions remaining after removing solutions from the first layer, and so on until all solutions are hierarchically sorted. The Pareto front distribution is composed of the non-dominated solutions in the first layer. These solutions form a front curve in the performance equilibrium graph. The front curve is extracted by connecting adjacent solutions in the first layer. The connection method uses straight line segments or curve fitting, such as piecewise linear interpolation or polynomial fitting. Piecewise linear interpolation is achieved by sequentially connecting adjacent points in the first layer solutions after sorting by the first quantization index value. Polynomial fitting uses the least squares method to fit the data points of the first layer solutions to generate a smooth curve.

[0044] When selecting the Pareto optimal solution from all biomimetic thin-walled structure schemes as the final preferred scheme, the selection process is based on the Pareto optimal solution set obtained by the non-dominated sorting algorithm. The Pareto optimal solution set corresponds to scatter points on the Pareto front distribution in the performance equilibrium diagram. Each scatter point is associated with an identifier of a biomimetic thin-walled structure scheme. The scheme identifier is assigned when the three-dimensional geometric model is built in step S1, for example, using the scheme number or a combination of characteristic parameters as the identifier. During the selection process, schemes corresponding to all solutions in the Pareto optimal solution set are selected. These schemes constitute the final preferred scheme set. The preferred scheme set may contain multiple schemes, and the number of schemes depends on the number of Pareto optimal solutions. For example, when there are 5 Pareto optimal solutions, there are 5 preferred schemes. The output form of the final preferred scheme is a scheme list, which contains the identifier of each scheme and the corresponding quantification index value. The list sorting is determined according to the position of the scheme on the Pareto front, for example, sorting from the smallest to the largest first quantification index, or sorting according to the second quantification index value. The final preferred scheme is used to guide the design selection of biomimetic thin-walled structures to ensure that the scheme achieves the optimal balance in terms of topological feature persistence.

[0045] All calculations involved in the embodiments are dimensionless numerical calculations, and the preset parameters and thresholds in the calculations are set by those skilled in the art according to the actual situation.

[0046] It should be noted that this invention can be deployed on the device itself to realize embedded applications, or it can run on a PC or other terminal with a user interface, thereby meeting various hardware environments and usage requirements.

[0047] The above embodiments can be implemented, in whole or in part, by software, hardware, firmware, or any other combination thereof. When implemented using software, the above embodiments can be implemented, in whole or in part, as a computer program product. The computer program product includes one or more computer instructions or computer programs. When the computer instructions or computer programs are loaded or executed on a computer, all or part of the processes or functions described in the embodiments of this application are generated. The computer can be a general-purpose computer, a special-purpose computer, a computer network, or other programmable device. The computer instructions can be stored in a computer-readable storage medium or transmitted from one computer-readable storage medium to another. For example, the computer instructions can be transmitted from one website, computer, server, or data center to another website, computer, server, or data center via wireless or wired transmission; wired transmission methods include optical fiber, twisted pair, coaxial cable, etc.; wireless transmission includes infrared, microwave, etc. The computer-readable storage medium can be any available medium that a computer can access or a data storage device such as a server or data center containing one or more sets of available media. The available medium can be a magnetic medium (e.g., floppy disk, hard disk, magnetic tape), an optical medium (e.g., DVD), or a semiconductor medium. A semiconductor medium can be a solid-state drive.

[0048] Those skilled in the art will understand that, for the sake of convenience and brevity, the specific working processes of the systems, devices, and modules described above can be referred to the corresponding processes in the foregoing method embodiments, and will not be repeated here.

[0049] In the several embodiments provided in this application, it should be understood that the disclosed systems, apparatuses, and methods can be implemented in other ways. For example, the apparatus embodiments described above are merely illustrative; for instance, the division of modules is only a logical functional division, and in actual implementation, there may be other division methods. For example, multiple modules or components may be combined or integrated into another system, or some features may be ignored or not executed. Furthermore, the coupling or direct coupling or communication connection shown or discussed may be through some interfaces; the indirect coupling or communication connection between apparatuses or modules may be electrical, mechanical, or other forms.

[0050] The modules described as separate components may or may not be physically separate. The components shown as modules may or may not be physical modules; they may be located in one place or distributed across multiple network modules. Some or all of the modules can be selected to achieve the purpose of this embodiment according to actual needs.

[0051] In addition, the functional modules in the various embodiments of this application can be integrated into one processing module, or each module can exist physically separately, or two or more modules can be integrated into one module.

[0052] If the aforementioned functions are implemented as software functional modules and sold or used as independent products, they can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of this application, in essence, or the part that contributes to the prior art, or a portion of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the steps of the methods described in the various embodiments of this application. The aforementioned storage medium includes various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks.

[0053] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

[0054] In conclusion, the above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for analyzing the impact resistance of biomimetic thin-walled structures based on finite element simulation, characterized in that, include: S1. Establish a three-dimensional geometric model of the biomimetic thin-walled structure to be evaluated; S2. Perform finite element mesh generation on the three-dimensional geometric model and assign material properties and boundary conditions to construct a finite element model for collision simulation. S3. Apply dynamic impact loads to the finite element model and perform collision simulation calculations; S4. Extract the deformation mode sequence of the biomimetic thin-walled structure from the results of the collision simulation calculation. The deformation mode sequence contains the node coordinate data of multiple time steps arranged in chronological order. S5. Based on node coordinate data, perform persistent cohomology analysis on the geometric configuration of the biomimetic thin-walled structure during the collision process, and extract at least one quantitative index to characterize the persistence of topological features. S6. Construct a performance equilibrium diagram based on at least one quantitative index to characterize the constraint relationship between the quantitative indices, and identify the Pareto front distribution in the performance equilibrium diagram to screen out the biomimetic thin-walled structure scheme located at the Pareto front.

2. The method for analyzing the impact resistance of biomimetic thin-walled structures based on finite element simulation according to claim 1, characterized in that, A three-dimensional geometric model of the biomimetic thin-walled structure to be evaluated is established, including: Based on predefined biomimetic structural feature parameters, a three-dimensional geometric model of a biomimetic thin-walled structure with at least one of the following features—bamboo-like cavity partitioning, skeletal gradient porous structure, or honeycomb rib structure—is generated using a parametric modeling method. The parametric modeling method controls the size and spatial distribution of the geometric configuration by adjusting the feature parameters.

3. The method for analyzing the impact resistance of biomimetic thin-walled structures based on finite element simulation according to claim 1, characterized in that, The 3D geometric model is meshed using finite element methods, and material properties and boundary conditions are assigned to it to construct a finite element model for collision simulation, including: The three-dimensional geometric model is divided using a hybrid mesh of tetrahedrons and hexahedrons, and local mesh refinement is implemented in characteristic regions with bamboo-like cavity partitions, skeletal gradient porous or honeycomb rib structures. The material distribution properties based on the three-dimensional geometric model are assigned corresponding elastoplastic material parameters, and fixed constraint boundary conditions are applied at the end of the finite element model to simulate the actual installation state.

4. The method for analyzing the impact resistance of biomimetic thin-walled structures based on finite element simulation according to claim 1, characterized in that, Dynamic impact loads are applied to the finite element model, and collision simulation calculations are performed, including: A dynamic impact load along the axial direction is applied to one end of the finite element model. The dynamic impact load is achieved by giving the impact mass block an initial velocity or by using a rigid wall displacement loading method. An explicit dynamic algorithm is used to solve the transient response of the finite element model during the collision process, where the impact velocity is set to a constant value that conforms to the actual working conditions.

5. The method for analyzing the impact resistance of biomimetic thin-walled structures based on finite element simulation according to claim 4, characterized in that, The explicit dynamic algorithm is used to solve the transient response of the finite element model during the collision process. This includes solving the equations of motion using the explicit time integration method, calculating the nodal displacement, velocity and acceleration at each time step, and considering the nonlinear constitutive relations of the material and the contact interaction.

6. The method for analyzing the impact resistance of biomimetic thin-walled structures based on finite element simulation according to claim 1, characterized in that, From the results of the collision simulation, a sequence of deformation modes for the biomimetic thin-walled structure is extracted. This sequence contains nodal coordinate data for multiple time steps arranged in chronological order, including: Multiple time steps are selected from the collision simulation results at equal time intervals. The spatial coordinates of all nodes of the biomimetic thin-walled structure at each time step are read using the post-processing function of the finite element analysis software. The node coordinate data arranged in chronological order are integrated into a deformation mode sequence, where the node coordinate data is stored in the form of a three-dimensional coordinate matrix.

7. The method for analyzing the impact resistance of biomimetic thin-walled structures based on finite element simulation according to claim 1, characterized in that, Persistent cohomology analysis of the geometric configuration of the biomimetic thin-walled structure during the collision process is performed based on nodal coordinate data, and at least one quantitative index for characterizing the persistence of topological features is extracted, including: Point cloud data representing geometric configurations is constructed based on node coordinate data; A distance-based filtering process is applied to point cloud data to generate a filtered complex sequence; Persistent graphs are obtained by calculating the homology group evolution of filtered complex sequences; Extract the lifetime distribution of topological features from the persistent graph, and calculate the average lifetime of the lifetime distribution as the first quantification index, while calculating the persistence entropy of the lifetime distribution as the second quantification index.

8. The method for analyzing the impact resistance of biomimetic thin-walled structures based on finite element simulation according to claim 7, characterized in that, The process of generating a filtered complex sequence by applying distance-based filtering to point cloud data includes: defining a distance filtering function, progressively increasing the distance threshold, and constructing a Vietoris-Rips complex sequence.

9. The method for analyzing the impact resistance of biomimetic thin-walled structures based on finite element simulation according to claim 7, characterized in that, Obtaining a persistent graph by calculating the homology group evolution of the filtered complex sequence includes: calculating the homology group of each filtered complex, tracking the birth and death of topological features, and generating a persistent graph.

10. The method for analyzing the impact resistance of biomimetic thin-walled structures based on finite element simulation according to claim 1, characterized in that, A performance equilibrium diagram is constructed based on at least one quantitative index to characterize the constraints between the quantitative indices. The Pareto front distribution in the performance equilibrium diagram is identified to screen for biomimetic thin-walled structure schemes located at the Pareto front, including: A performance equilibrium diagram was constructed using the first and second quantitative indicators as two-dimensional coordinate axes, and the quantitative indicator values ​​corresponding to multiple biomimetic thin-walled structure schemes were plotted as a scatter distribution. The Pareto front distribution in the performance equilibrium graph is identified by the non-dominated sorting algorithm in multi-objective optimization. The Pareto front distribution is composed of Pareto optimal solutions that are not dominated by other solutions. The structural scheme corresponding to the Pareto optimal solution is selected from all biomimetic thin-walled structural schemes as the final preferred scheme.

Images