Methods for establishing a finite element model of craniocerebral blast injury

Abstract

This invention relates to the field of twin simulation technology, specifically a method for establishing a finite element model of craniocerebral blast injury, comprising the following steps: acquiring tomographic images, extracting two-dimensional boundary node coordinates to generate a three-dimensional geometric solid model; performing tetrahedral mesh generation and adjusting deformation node positions to generate a head solid mesh model; dividing the outer region into Eulerian elements and performing node pairing operations to establish a fluid-structure interaction mesh; applying temporal pressure to the outer boundary of the Eulerian elements and setting non-reflective physical constraints to construct an initial model. In this invention, by extracting mesh deformation parameters to compensate for vertex centripetal displacement to improve topological connectivity, pairing external flow field control points with internal solid field boundaries to establish a fluid-structure interaction network transmission link to restore the structural collision energy dissipation mechanism, establishing a non-reflective environment to accurately project pressure and prevent false interference, and constructing a collaborative physical response mechanism to approximate the true in-situ deformation mechanical properties under extreme impact conditions.

Description

Technical Field

[0001] This invention relates to the field of twin simulation technology, and in particular to a method for establishing a finite element model of cranial blast injury. Background Technology

[0002] The field of twin simulation technology mainly involves constructing a mapping relationship between physical entities and virtual models to achieve dynamic mapping and behavior reproduction of real systems in digital space. Its core aspects include multi-source data acquisition and fusion, geometric structure modeling, parameterized description of physical properties, boundary condition setting, and synchronous updates of the numerical solution process. Relying on computational mechanics, computational fluid dynamics, and multiphysics coupling analysis methods, it performs refined time-domain and spatial-domain simulations of the response process of complex structures under external loads, thereby reproducing processes such as impact, vibration, thermal effects, and nonlinear material deformation in a virtual environment. It is widely used in aerospace, biomedical engineering, and safety protection fields.

[0003] The traditional finite element model for craniocerebral blast injury refers to the numerical modeling process for the injury mechanism of the skull and brain tissue under explosive impact loads. The technical issues it addresses are: constructing a discretized computational model that can reflect the interaction between the skull, brain parenchyma, and cerebrospinal fluid; obtaining the geometric structure of the head using a three-dimensional reconstruction method based on medical imaging data; establishing a computational mesh by dividing the skull into shell elements and the brain tissue and cerebrospinal fluid into solid elements; assigning material properties based on the density, elastic modulus, and Poisson's ratio of different tissues; applying an impact pressure load described in the form of a time history outside the model; setting constraints at the skull base or neck position; and gradually calculating nodal displacements and element stresses using an explicit time integration method to complete the model construction process.

[0004] Traditional methods construct geometric structures based on medical images and perform isolated meshing of phasic and volumetric regions. They assign material properties based on tissue differences and apply time-series loads to the external boundaries in conjunction with end constraints. The discreteness of a single structure leads to fuzzy representation of the geometric features of the interface between soft and hard tissues. Stress loading in isolated fields ignores the dynamic interaction and collision energy decay process between the fluid medium and the solid structure. Rigid boundary constraints cause abnormal concentrated reflection of ripples under transient impact conditions, resulting in distortion of internal strain distribution and deviation in damage assessment. Summary of the Invention

[0005] To address the technical problems of traditional methods that construct geometric structures based on medical images and perform isolated meshing of spherical and volumetric regions, assign material properties based on tissue differences and apply time-series loads to external boundaries in conjunction with end constraints, the discreteness of a single structure leads to fuzzy representation of the geometric features at the interface between soft and hard tissues, stress loading in isolated fields ignores the dynamic interaction and collision energy attenuation process between the fluid medium and the solid structure, and rigid boundary constraints cause abnormal concentrated reflection of ripples under transient impact conditions, resulting in distortion of internal strain distribution and deviation in damage assessment, this invention provides a method for establishing a finite element model of craniocerebral blast injury.

[0006] To achieve the above objectives, this invention employs a method for establishing a finite element model of craniocerebral blast injury, comprising the following steps: S1: Acquire tomographic images of healthy cranial soft tissue and hard tissue from MRI and CT scanners, input them into the image registration and fusion algorithm to perform grayscale feature calculation and cross-sectional contour extraction, extract the spatial coordinates of the nodes of the two-dimensional contour boundary, and perform inter-layer normal interpolation on the spatial coordinates to generate a three-dimensional geometric solid model. S2: Perform tetrahedral meshing on the three-dimensional geometric solid model, extract the coordinates of the mesh unit nodes and calculate the unit volume deformation factor, and adjust the spatial displacement coordinates of nodes below the preset deformation threshold towards the geometric center of the unit to generate the head solid mesh model. S3: Extract the coordinates of the nodes on the outer surface of the head solid mesh model, divide the external space region into a hexahedral air Eulerian element structure, calculate the spatial distance between the coordinates of the internal element nodes and the outer surface nodes, and perform bidirectional contact node pairing operation to establish a fluid-structure interaction mesh structure. S4: Apply the shock wave pressure time series data collected by the pressure sensor in the experimental environment to the Euler outer boundary node of the fluid-structure interaction mesh structure, set the non-reflective physical boundary constraint, obtain the tissue material property data and assign it to the internal solid element node of the fluid-structure interaction mesh structure to construct the initial finite element model of craniocerebral blast shock. S5: Based on the non-reflective physical boundary constraints, the initial finite element model of craniocerebral blast injury is input into the explicit dynamic solver, and step difference operation is performed to extract the stress values ​​of the model nodes; the stress values ​​of the model nodes are compared with the experimental stress test values ​​to perform error comparison operation, and the state is solidified when the difference is lower than the preset convergence threshold to generate the target craniocerebral blast injury finite element model.

[0007] As a further aspect of the present invention, the three-dimensional geometric solid model includes the three-dimensional morphology of the skull, the anatomical configuration of the ventricles, and the surface contour of the scalp; the head solid mesh model includes a tetrahedral element set, mesh node coordinates, and element connection topology matrix; the fluid-structure interaction mesh structure includes an Eulerian computational domain, a Lagrangian solid domain, and a fluid-structure interaction interface; the initial finite element model of craniocerebral blast injury includes a hyperelastic constitutive model, bulk modulus, and state equation parameters; the target craniocerebral blast injury finite element model includes an experimentally calibrated mesh topology, corrected material constitutive parameters, and contact penalty function constraints.

[0008] As a further aspect of the present invention, the specific steps of S1 are as follows: S101: Acquire soft tissue tomographic images and hard tissue tomographic images using MRI and tomography equipment, perform three-dimensional translation and alignment, extract the single-channel grayscale components corresponding to the overlapping areas, multiply the components with the same coordinates and take the absolute value as the target grayscale feature, and map them one by one into matrix elements according to the spatial voxel coordinates corresponding to the target grayscale feature to establish a soft and hard tissue feature matrix. S102: Call the soft and hard tissue feature matrix, calculate the second derivative of the pixel gray level gradient and filter edge points that exceed the preset mutation threshold, perform spatial connected component labeling and aggregate pixels with the same label, construct a set of closed curves, extract discrete node coordinate sequences from the set of closed curves, and obtain a two-dimensional contour node coordinate set. S103: Based on the coordinate set of the two-dimensional contour nodes, calculate the normal vector of the corresponding node of the adjacent fault plane, determine whether the cosine value of the included angle is within the preset parallel judgment interval, if it is, use the corresponding node as the matching node for search, construct the cross-layer three-dimensional related line segment, perform linear geometric interpolation along the related line segment to supplement the support point, and generate a three-dimensional geometric solid model.

[0009] As a further aspect of the present invention, the specific steps of S2 are as follows: S201: Obtain the three-dimensional geometric entity model, implant discrete seed nodes along the surface and interior of the entity space boundary, connect adjacent nodes to construct a tetrahedral structure according to the spatial topology rules, extract the absolute position vectors of the grid unit vertices and splice them into a numerical arrangement, and establish a three-dimensional coordinate matrix of the grid unit. S202: Based on the three-dimensional coordinate matrix of the mesh unit, read the absolute position vectors of the four vertices, calculate the absolute value of the mixed product of adjacent common vertex edge vectors to obtain the actual unit volume, extract the standard reference volume of the same circumscribed sphere tetrahedron, calculate the ratio of the two types of volumes and record it as the distortion state parameter, and obtain the set of unit volume deformation factors. S203: Call the set of unit volume deformation factors, extract the target deformation nodes whose distortion state parameters are lower than the preset deformation threshold, perform an arithmetic mean operation on the absolute position vectors of the four vertices of the unit to extract the geometric center coordinates, construct the displacement line segment of the node pointing to the geometric center coordinates and adjust the coordinates to generate the head solid mesh model.

[0010] As a further aspect of the present invention, the preset deformation threshold is determined based on the statistical distribution of the unit volume deformation factor set. The preset deformation threshold is obtained by sorting the unit volume deformation factor set and selecting the distortion state parameter value corresponding to the preset quantile position.

[0011] As a further aspect of the present invention, the specific steps of S3 are as follows: S301: Obtain the head entity mesh model, extract the discrete coordinate node set along the outer boundary of the head entity mesh model, construct an orthogonal hexahedral air Euler element array in the external space according to a preset scale, perform topological relative position alignment determination on the discrete coordinate node set and the element array, and generate the boundary interaction initial coordinate matrix. S302: Call the initial coordinate matrix of the boundary interaction, traverse the coordinate nodes of the unit nodes and the outer surface coordinate nodes of the head entity mesh model, quantify the node span according to the Euclidean distance calculation rule to extract the distance value, compare the distance value with the preset contact tolerance limit, filter the candidate associated node set and remove isolated nodes, and generate a near-end potential contact topology map. S303: Based on the near-end potential contact topology map, perform bidirectional data interaction mapping pairing operation on the candidate associated node set, between the flow field Euler control nodes and the coordinate nodes on the outer surface of the head entity mesh model, force the mapping pairing attributes into the global node sequence space, aggregate the internal and external field associated data elements, and establish a fluid-structure interaction mesh structure.

[0012] As a further aspect of the present invention, the specific steps of S4 are as follows: S401: Extract the spatial coordinates of the Euler nodes on the outer side of the fluid-structure interaction mesh structure, collect the time series data of shock wave pressure from the pressure sensor, map the pressure time series data into a unidirectional scalar field of Euler nodes based on the spatial position correspondence, aggregate the dynamic amplitude of the load and the spatial coordinates, and establish a spatiotemporal distribution array of the boundary load. S402: Based on the boundary load spatiotemporal distribution array, scan the extreme coordinate topological boundary domain of the fluid-structure interaction mesh structure, extract the outermost free state node set, apply unidirectional energy dissipation transmission property to the outermost free state node set and lock the wavefront reflection scalar flux, map the node index and constraint state association key value, and generate non-reflective physical boundary constraints. S403: Invoke the non-reflective physical boundary constraints, read the set of physiological tissue material property parameters, lock the solid entity unit nodes inside the fluid-structure interaction mesh structure and perform material parameter assignment operations, bind the density mechanical modulus and node sequence index, integrate the global node position, external boundary constraints and internal field physical property variables, and construct the initial finite element model of craniocerebral blast shock.

[0013] As a further aspect of the present invention, the specific steps of S5 are as follows: S501: Call the non-reflective physical boundary constraints and the initial finite element model of craniocerebral blast shock, perform time domain discretization, perform algebraic equation calculations of displacement strain matrix for each time step, and converge the three-dimensional spatial component vectors of the vertices of the internal solid elements of the time node to obtain the stress time series tensor of the model node. S502: Perform time point registration and alignment operation on the model node stress time series tensor and the experimental standard stress test sequence parameters, calculate the absolute value of the stress scalar difference at the corresponding time point, add the absolute values ​​of the difference and perform time domain cross-node integration operation to generate a set of node stress error scalars. S503: Based on the set of nodal stress error scalars, compare the error scalars with the preset convergence threshold, remove the corresponding dynamic incremental step sequences that exceed the preset convergence threshold, extract the associated geometric coordinates and deformation attributes that do not exceed the preset convergence threshold, perform mesh freezing, and establish a finite element model of the target craniocerebral blast injury.

[0014] As a further aspect of the present invention, the preset convergence threshold is based on the upper limit of a fixed numerical range determined by the statistical distribution of the nodal stress error scalar set. The removal of the corresponding dynamic incremental step sequence that exceeds the preset convergence threshold means that the node stress error scalar set is arranged in chronological order and the difference change rate is calculated for adjacent time nodes. The difference change rate is compared with the preset convergence threshold one by one. The dynamic incremental step sequence corresponding to the consecutive time nodes that satisfy the difference change rate is less than or equal to the preset convergence threshold is retained. The associated geometric coordinates and deformation attributes that do not exceed the preset convergence threshold refer to the combination data of the three-dimensional coordinates and displacement vectors of the corresponding entity unit vertices at the same time node. Consistency verification processing is performed on the combination data to filter the set of entity units that meet the spatial continuity.

[0015] Compared with the prior art, the advantages and positive effects of the present invention are as follows: In this invention, medical image registration and fusion are performed to extract the spatial coordinates of the contour boundary and reconstruct the refined solid configuration by normal interpolation to eliminate the expression deviation defects of the interface between soft and hard tissues. The deformation state parameters of discrete grid units are extracted and centripetal displacement compensation is performed on the distorted vertices to improve the quality of topological connection. The distance between the gas phase control point of the external flow field and the physical boundary of the internal solid field is calculated to complete the bidirectional mapping pairing interaction to open up the fluid-structure coupling transmission link and restore the real physical mechanism of shock wave energy dissipation. A reflection-free constraint environment is established and the temporal pressure data is accurately projected onto the outer correlation domain to prevent false wavefront interference. A multi-physics field collaborative response network is constructed to effectively approximate the in-situ deformation mechanical mechanism under high-risk injury conditions. Attached Figure Description

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

[0017] Figure 1 This is a schematic diagram of the steps of the present invention; Figure 2 This is a detailed schematic diagram of S1 of the present invention; Figure 3 This is a detailed schematic diagram of S2 of the present invention; Figure 4 This is a detailed schematic diagram of S3 of the present invention; Figure 5 This is a detailed schematic diagram of S4 of the present invention; Figure 6 This is a detailed schematic diagram of S5 of the present invention. Detailed Implementation

[0018] The technical solution of the present invention will now be described with reference to the accompanying drawings.

[0019] To make the technical problems, technical solutions and advantages of the present invention clearer, a detailed description will be given below in conjunction with the accompanying drawings and specific embodiments.

[0020] Please see Figure 1 This invention provides a method for establishing a finite element model of craniocerebral blast injury, comprising the following steps: S1: Acquire tomographic images of healthy cranial soft tissue and hard tissue from MRI and CT scanners, input them into the image registration and fusion algorithm to perform grayscale feature calculation and cross-sectional contour extraction, extract the spatial coordinates of the nodes of the two-dimensional contour boundary, and perform inter-layer normal interpolation on the spatial coordinates to generate a three-dimensional geometric solid model. S2: Perform tetrahedral mesh generation on the three-dimensional geometric solid model, extract the coordinates of the mesh unit nodes and calculate the unit volume deformation factor, and adjust the spatial displacement coordinates of nodes below the preset deformation threshold towards the geometric center of the unit to generate the head solid mesh model. S3: Extract the coordinates of the nodes on the outer surface of the head solid mesh model, divide the external space region into a hexahedral air Eulerian element structure, calculate the spatial distance between the coordinates of the internal element nodes and the nodes on the outer surface, and perform bidirectional contact node pairing operation to establish a fluid-structure interaction mesh structure. S4: Apply the shock wave pressure time series data collected by the pressure sensor in the experimental environment to the Euler outer boundary node of the fluid-structure interaction mesh structure, set the non-reflective physical boundary constraint, obtain the tissue material property data and assign it to the internal solid element node of the fluid-structure interaction mesh structure to construct the initial finite element model of craniocerebral blast shock. S5: Based on the non-reflective physical boundary constraints, the initial finite element model of craniocerebral blast injury is input into the explicit dynamic solver, and the step difference operation is performed to extract the stress values ​​of the model nodes. The stress values ​​of the model nodes are compared with the experimental stress test values. When the difference is lower than the preset convergence threshold, the state is solidified to generate the target craniocerebral blast injury finite element model. The three-dimensional geometric solid model includes the three-dimensional morphology of the skull, the anatomical configuration of the ventricles, and the surface contour of the scalp. The head solid mesh model includes a set of tetrahedral elements, mesh node coordinates, and element connection topology matrix. The fluid-structure interaction mesh structure includes the Eulerian computational domain, the Lagrangian solid domain, and the fluid-structure interaction interface. The initial finite element model of craniocerebral blast injury includes a hyperelastic constitutive model, bulk modulus, and state equation parameters. The finite element model of the target craniocerebral blast injury includes an experimentally calibrated mesh topology, corrected material constitutive parameters, and contact penalty function constraints.

[0021] Please see Figure 2 The specific steps of S1 are as follows: S101: Acquire soft tissue tomographic images and hard tissue tomographic images using MRI and tomography equipment, perform three-dimensional translation and alignment, extract the single-channel grayscale components corresponding to the overlapping areas, multiply the components with the same coordinates and take the absolute value as the target grayscale feature, and map them one by one into matrix elements according to the spatial voxel coordinates corresponding to the target grayscale feature to establish a soft and hard tissue feature matrix. A sequence of soft tissue tomographic images of the head with a slice thickness of 1.0 mm, output from a high-field-strength MRI scanner, was acquired. Simultaneously, a sequence of hard tissue tomographic images with a slice thickness of 1.0 mm, output from a computed tomography (CT) scanner, was read. These two types of low-level image data were imported into the memory space of the image registration and fusion module. Bilinear interpolation resampling was performed on the original scanned images to unify the spatial resolution of the two sets of images and eliminate differences in the acquisition field of view. The three-dimensional coordinate points of the scalp and brain parenchyma periphery of the soft tissue tomographic images were extracted as the primary registration reference point set, and the coordinate points of the inner and outer plates of the skull of the hard tissue tomographic images were extracted as the secondary registration reference point set. A three-dimensional translation and alignment operation was performed on the primary and secondary registration reference point sets. The spatial translation vector that minimizes the sum of the squares of the Euclidean distance differences between the two sets of points in the same three-dimensional Cartesian coordinate system was calculated. This spatial translation vector was applied to perform a global three-dimensional coordinate shift on the hard tissue tomographic images, achieving a strict topological correspondence between the soft and hard tissue tomographic images in physical space. The tomographic image voxel grid after 3D translation and alignment is traversed row by row. Voxel values ​​of the two types of images at the same coordinate position are compared, and the 3D continuous space containing effective anatomical structure signals is identified and defined as the overlapping region. For each voxel coordinate within the overlapping region, the single-channel soft tissue grayscale component responding to that coordinate in the soft tissue tomographic image is extracted, and the single-channel hard tissue grayscale component responding to that coordinate in the hard tissue tomographic image is also extracted. The single-channel soft tissue grayscale component and the single-channel hard tissue grayscale component are multiplied to obtain a mixed image signal product. Then, the absolute value of this mixed image signal product is taken to eliminate instrument inversion artifacts and negative values. The final absolute value is used as the target grayscale feature of the voxel in that space. For example, at a voxel coordinate at the junction of the parietal bone and dura mater, the single-channel soft tissue grayscale component value acquired by an MRI scanner is 150, and the single-channel hard tissue grayscale component value acquired by a computed tomography (CT) scanner at the same coordinate is 210. Multiplying the single-channel soft tissue grayscale component 150 and the single-channel hard tissue grayscale component 210 yields a mixed image signal product of 31500. After absolute value extraction, the target grayscale feature is 31500. A three-dimensional empty matrix is ​​initialized in memory. Based on the voxel coordinate sequence, the target grayscale feature values ​​of all voxels within the overlapping region are mapped one by one to the corresponding row and column depth positions of the associative array as matrix elements, thereby establishing a global soft and hard tissue feature matrix.

[0022] S102: Call the soft and hard tissue feature matrix, calculate the second derivative of the pixel gray-level gradient and filter edge points that exceed the preset mutation threshold, perform spatial connected component labeling and aggregate pixels with the same label, construct a set of closed curves, extract discrete node coordinate sequences from the set of closed curves, and obtain a two-dimensional contour node coordinate set. The system calls upon the soft and hard tissue feature matrix and reads the two-dimensional slice data layer by layer within the matrix according to the Z-axis slicing order. The two-dimensional slice data represents the feature distribution on a single physical cross-section. The Laplacian second-order differential operator is applied to the read two-dimensional slice data to calculate the target gray-level feature difference between the current center pixel and its eight neighboring pixels. The feature values ​​of the eight neighboring pixels are summed, and then eight times the feature value of the center pixel is subtracted to calculate the second derivative of the pixel gray-level gradient of the center pixel. After performing the second derivative calculation on all pixels within the slice, the arithmetic mean and standard deviation of the second derivative values ​​for the entire slice are calculated. The sum of the arithmetic mean and standard deviation (3 times the standard deviation) is used as the preset mutation threshold. The second derivative of the pixel gray-level gradient of each pixel within the slice is compared one by one with this preset mutation threshold, and spatial coordinates points whose absolute second derivative value is greater than or equal to the preset mutation threshold are selected as edge points. For the selected unordered edge points, a spatial connectivity labeling operation is performed in the two-dimensional grid space. The 8-connected neighborhoods of the edge points are scanned, and adjacent edge points are assigned the same integer topological label. Pixels with the same label are aggregated to construct connected sets, and tiny, isolated noise connected sets with fewer than 50 pixels are removed. Connected sets forming closed loops are retained to construct a set of closed curves. For each closed curve in the set, the physical horizontal and vertical coordinates of the pixels on the curve are extracted sequentially in a clockwise direction to form a discrete node coordinate sequence. The discrete node coordinate sequences from all tomographic slices are summarized to obtain a complete two-dimensional contour node coordinate set describing the shape of various tissues in the head. For example, in a local slice analysis, the sum of the feature values ​​of the eight neighboring pixels around the center pixel is 160,000, and the feature value of the center pixel is 15,000. Multiplying the feature value of the center pixel by 8 gives 120,000. Subtracting 120,000 from the sum of the surrounding feature values ​​of 160,000 gives the second derivative of the pixel grayscale gradient of that point as 40,000. The preset mutation threshold obtained from the current slice statistics is 35,000. Since the second derivative of 40,000 of this point is greater than the preset mutation threshold of 35,000, this pixel is filtered and marked as an edge point.

[0023] S103: Based on the coordinate set of two-dimensional contour nodes, calculate the normal vector of the corresponding node of the adjacent fault plane, determine whether the cosine value of the included angle is within the preset parallel judgment interval, if it is, use the corresponding node as the matching node for search, construct the cross-layer three-dimensional related line segment, perform linear geometric interpolation along the related line segment to supplement the support point, and generate a three-dimensional geometric solid model. Based on the acquired 2D contour node coordinate set, discrete node sequences on two adjacent physical fault planes are extracted. For any reference node on the upper fault plane, the tangential vectors of its two adjacent nodes on the closed curve are calculated. Then, the direction vector perpendicular to the tangential vector and pointing towards the center of curvature is calculated as the normal vector of the reference node. Similarly, the corresponding node normal vectors of all nodes on the lower fault plane are calculated. The normal vectors of the upper reference node and the normal vectors of the lower node to be matched are extracted and subjected to a dot product operation. The dot product result is divided by the product of the magnitudes of the two normal vectors to obtain the cosine value of the angle between the two node normal vectors. It is determined whether the cosine value of the angle is within a preset parallel judgment interval, with a lower limit of 0.95 and an upper limit of 1.0. If the cosine value of the angle is within this parallel judgment interval, the normal vectors of the two nodes are determined to be parallel, and the lower node to be matched is selected as the search matching node. Connecting the upper-level reference node with the lower-level search matching node constructs a cross-layer 3D associated line segment. When the distance between adjacent fault planes is 1.0 mm, the length space of this cross-layer 3D associated line segment is equidistantly subdivided. The difference between the coordinates of the two endpoints along the associated line segment is calculated and divided by the number of subdivisions to obtain the interpolation step vector. Starting from the upper-level endpoint coordinates, the interpolation step vector is continuously accumulated to perform linear geometric interpolation, supplementing the space between the two fault planes with additional 3D support point coordinates. All original discrete nodes, search matching nodes, and newly generated support points are reconstructed using 3D spatial meshing to generate a 3D geometric solid model containing closed surfaces of the scalp, skull, and brain parenchyma. For example, after calculating the 3D vector and normal vector of the upper-level fault plane reference node, a dot product and modulus division operation is performed with the normal vector of a lower-level node to obtain an angle cosine value of 0.98. Since 0.98 is within the parallel judgment interval of 0.95 to 1.0, a corresponding matching relationship is established and a cross-layer 3D associated line segment is generated. The upper endpoint of the line segment has a height elevation of 10.0 mm, and the lower endpoint has a height elevation of 9.0 mm. If the number of subdivision segments is set to 2, the interpolation step size is obtained by dividing the difference in endpoint coordinates by 2. The height elevation of the support point supplemented by linear geometric interpolation is 9.5 mm.

[0024] Please see Figure 3 The specific steps of S2 are as follows: S201: Obtain a three-dimensional geometric solid model, implant discrete seed nodes along the surface and interior of the solid space boundary, connect adjacent nodes to construct a tetrahedral structure according to the spatial topology rules, extract the absolute position vectors of the grid unit vertices and splice them into a numerical arrangement, and establish a three-dimensional coordinate matrix of the grid unit. A 3D geometric solid model containing the boundaries of each physiological structure is obtained and imported into a 3D spatial array of the mesh generation module. Discrete seed nodes are implanted along the spatial boundary surfaces of the solid model and the internal volume space of each tissue, according to the set mesh size requirements of 1.0 mm to 3.0 mm. For the implanted set of discrete seed nodes, the Delaunay spatial topology rule is used to generate a spatial mesh. Any four non-coplanar discrete seed nodes are selected and connected to construct a tetrahedral structure. The circumsphere empty circle property is checked on the tetrahedral structure, that is, to ensure that the circumsphere formed by the four nodes does not contain any other discrete seed nodes. If the empty circle property is satisfied, the connection relationship constructed by these four nodes is retained as a valid mesh unit. After traversing all discrete seed nodes to complete the full coverage of the spatial tetrahedral structure, for each generated tetrahedral mesh unit, the absolute position vectors of its four vertices are extracted, that is, the horizontal, vertical and depth values ​​of the vertices in the global 3D coordinate system. These four sets of absolute position vectors are spliced ​​into a continuous one-dimensional numerical arrangement according to a fixed order. Each mesh cell is assigned a unique cell number index, and the concatenated numerical values ​​are mapped to the corresponding index entries to establish a three-dimensional coordinate matrix of mesh cells that records the geometric characteristics of the entire mesh domain, for subsequent deformation verification and quality screening. Table 1 lists the configuration of the tetrahedral mesh cell node coordinate data extracted at the boundaries of some solid models.

[0025] Table 1: Characteristics of Mesh Cell Node Coordinate Array

[0026] As shown in Table 1, by extracting the absolute spatial coordinates of each vertex of the mesh unit, the relative position and geometric dimensions of each tetrahedron in the solid space can be accurately defined.

[0027] S202: Based on the three-dimensional coordinate matrix of the mesh element, read the absolute position vectors of the four vertices, calculate the absolute value of the mixed product of adjacent common vertex edge vectors to obtain the actual element volume, extract the standard reference volume of the tetrahedron with the same circumscribed sphere, calculate the ratio of the two types of volumes as the distortion state parameter, and obtain the set of element volume deformation factors. Based on the established 3D coordinate matrix of the mesh elements, the absolute position vectors of the four vertices corresponding to each tetrahedral mesh element are read sequentially. One vertex is selected as the reference origin, and the coordinates of the reference origin are subtracted from the absolute position vectors of the other three vertices to obtain three adjacent edge vectors sharing the same vertex origin. A mixed product operation is performed on these three adjacent edge vectors: first, the cross product of two edge vectors is calculated to generate a normal vector; then, the dot product of this normal vector and the third edge vector is performed to obtain the mixed product value. The absolute value of the mixed product value is taken, and the result is divided by 6 to calculate the actual element volume of the tetrahedral mesh element. The geometric radius of the circumsphere of the current tetrahedral mesh element is calculated based on the four vertices. Then, based on the geometric radius of the circumsphere, the standard reference volume of a regular tetrahedron enclosed by an equal circumsphere is derived. The actual element volume and the standard reference volume of the regular tetrahedron are divided by a ratio to calculate the ratio of the two volumes, which is recorded as the distortion state parameter of the mesh element. After traversing hundreds of thousands of tetrahedral mesh elements in the global model to complete the above calculations, the distortion state parameters of all mesh elements are summarized to obtain the set of element volume deformation factors from a global perspective. For example, at the boundary of a complex curved surface at the bottom of the anterior cranial fossa of the head, three adjacent edge vectors sharing the same vertex are calculated based on the absolute position vectors of the four vertices. After performing the mixed product, outer product, and inner product operations, the absolute value is 12.0 cubic millimeters. Dividing this by 6 yields an actual element volume of 2.0 cubic millimeters. Based on the circumscribed sphere radius calculated from its four vertices, the standard reference volume of a regular tetrahedron at this radius is derived to be 2.5 cubic millimeters. Dividing the actual element volume of 2.0 cubic millimeters by the standard reference volume of 2.5 cubic millimeters yields a volume ratio of 0.8. This 0.8 is the distortion state parameter of the current element; the closer the value is to 1.0, the closer the mesh is to the ideal regular tetrahedron state.

[0028] S203: Call the unit volume deformation factor set, extract the target deformation node whose distortion state parameter is lower than the preset deformation threshold, perform an arithmetic mean operation on the absolute position vectors of the four vertices of the unit to extract the geometric center coordinates, construct the node pointing to the geometric center coordinate displacement line segment and adjust the coordinates to generate the head solid mesh model. The system retrieves the aggregated set of element volume deformation factors and sorts all distortion state parameter values ​​in ascending order. It locates the 5th quantile position in the sorted sequence and extracts the corresponding distortion state parameter value as a preset deformation threshold. The system scans the entire mesh, extracting mesh nodes whose distortion state parameters are strictly below the preset deformation threshold as target deformation nodes. For the inferior tetrahedral element containing the selected target deformation node, the system reads the absolute position vectors of the four vertices of the element. It then adds the lateral, longitudinal, and depth coordinates of these four vertices, divides by 4, and performs an arithmetic mean calculation to extract the geometric center coordinates within the element. For the distorted vertex identified as a target deformation node, it constructs a spatial displacement segment pointing from the vertex's current coordinates to the geometric center coordinates. The system calculates 20% of the total length of this displacement segment as an adjustment step size and translates the target deformation node along this displacement segment towards the geometric center coordinates by this adjustment step size. This updates and adjusts the node's 3D coordinates, completing the local mesh correction and generating the final head solid mesh model. For example, in a model containing 200,000 tetrahedral elements, after all distortion state parameters are sorted in ascending order, the distortion state parameter value at the 10,000th position (i.e., the 5th quantile) is 0.45, thus determining the preset deformation threshold as 0.45. A certain tetrahedral element has a distortion state parameter of 0.35, which is lower than 0.45, and is therefore locked as an element requiring processing. The arithmetic mean of the coordinates of the four vertices of this element yields the geometric center coordinates as 10.0 mm horizontally, 20.0 mm vertically, and 30.0 mm deep. The initial coordinates of the target deformation node are 8.0 mm horizontally, 18.0 mm vertically, and 28.0 mm deep. The spatial distance between this node and the geometric center is calculated to be 3.46 mm, and the adjustment step size is 0.69 mm (20% adjustment step). The node is then moved 0.69 mm along the vector direction pointing to the geometric center, overwriting the original coordinate data.

[0029] Please see Figure 4 The specific steps of S3 are as follows: S301: Obtain the head solid mesh model, extract the discrete coordinate node set along the outer boundary of the head solid mesh model, construct an orthogonal hexahedral air Euler element array in the external space according to the preset scale, perform topological relative position alignment judgment on the discrete coordinate node set and the element array, and generate the boundary interaction initial coordinate matrix. A head entity mesh model with internal corrections is obtained. The outermost envelope structure of the global model, i.e., the topological polygonal patches of the scalp's outer surface, is scanned. All discrete coordinate node sets exposed to the external physical space are extracted along the outer model patches. A computational domain is established in the infinitely extending space outside the head entity model based on the hydrodynamic characteristics of shock wave propagation. According to a preset spatial scale, i.e., the length, width, and height dimensions of the orthogonal hexahedral elements are all set to 2.0 mm, and the minimum distance between the outer boundary of the computational domain and the nearest head entity boundary is set to 100.0 mm, a uniformly distributed array of orthogonal hexahedral air Eulerian elements is constructed in the external space. The coordinates of each scalp boundary node in the discrete coordinate node set are extracted and topologically aligned with the spatial center coordinates of each orthogonal hexahedral element in the air Eulerian element array. The alignment logic is based on whether the node coordinates are enveloped within the three-dimensional orthogonal boundary box of the corresponding air element. The correspondence between the outer air elements and the inner head nodes that satisfy the envelope relationship is recorded. This correspondence is integrated to generate an initial coordinate matrix of boundary interactions with internal and external spatial correspondence identifiers. For example, a bounding box is generated by extending outward by 100.0 mm based on the head geometry. Within this box, millions of orthogonal hexahedral nodes are generated at 2.0 mm intervals. The discrete nodes on the scalp are traversed. For a given scalp node with coordinates of 50.0 mm horizontally, 80.0 mm vertically, and 120.0 mm deep, specific Eulerian elements covering a boundary range of 49.0 to 51.0 mm horizontally, 79.0 to 81.0 mm vertically, and 119.0 to 121.0 mm deep are searched in the orthogonal hexahedral air Eulerian element array. The binding relationship between these elements is recorded and filled into the initial matrix.

[0030] S302: Call the initial coordinate matrix of the boundary interaction, traverse the coordinate nodes of the unit nodes and the outer surface of the head entity mesh model, quantify the node span according to the Euclidean distance calculation rules to extract the distance value, compare the distance value with the preset contact tolerance limit, filter the candidate associated node set and remove isolated nodes, and generate the near-end potential contact topology map. The generated initial coordinate matrix for boundary interaction is invoked, and the sequence of air Eulerian element nodes mapped within the matrix and the sequence of coordinate nodes on the outer surface of the head model are traversed. For each pair of pre-matched nodes, the node span is quantified according to the Euclidean distance calculation rule, i.e., the coordinate difference between the two nodes in the horizontal, vertical, and depth dimensions is calculated, the coordinate difference in each dimension is squared, the three squared values ​​are summed, and finally the square root of the sum is taken to extract the absolute spatial distance value. The obtained distance value is compared with a preset contact tolerance limit, which is set to 1.5 mm based on a comprehensive consideration of the fluid-structure mesh scale. If the calculated distance value is less than or equal to 1.5 mm, the pair of nodes is included in the candidate associated node set; if the distance value is greater than 1.5 mm, it is determined that the distance between the two is too great to produce effective fluid-structure interaction. After completing the global screening, the topological connections in the candidate associated node set are examined. Isolated flow field nodes that fail to form a mechanical transmission link with any effective solid field mesh element are identified and eliminated. The remaining effective node pairs are then organized to generate the final near-end potential contact topology map characterizing the explosion impact interface. Table 2 lists the distance quantification test results between boundary contact nodes.

[0031] Table 2: Quantization Decision Table for Spatial Span of Contact Nodes

[0032] As shown in Table 2, through detailed Euclidean spatial distance calculations and limit comparisons, it is possible to accurately separate free node pairs that do not meet the physical contact conditions.

[0033] S303: Based on the near-end potential contact topology map, perform bidirectional data interaction mapping and pairing operations on the flow field Euler control nodes and the coordinate nodes on the outer surface of the head entity mesh model for the candidate associated node set. Force the mapping and pairing attributes into the global node sequence space, aggregate the internal and external field associated data elements, and establish a fluid-structure interaction mesh structure. Based on the generated near-end potential contact topology, complex coupling variable transfers are performed on the locked candidate associated node set within the graph. Within each tiny analysis step of the time evolution, a bidirectional data interaction mapping pairing operation is performed between the external flow field Euler control nodes and the internal solid field head boundary nodes. Specifically, this involves extracting the transient fluid pressure and tangential shear stress variables from the flow field Euler control nodes and applying these force variables equivalently to the corresponding head solid mesh model's outer surface coordinate nodes using an inverse distance weighting algorithm. Simultaneously, the spatial displacement and dynamic deformation velocity variables calculated at the head solid mesh model's outer surface coordinate nodes in the current time step are extracted and used as boundary motion constraint feedback, written into the geometric boundary conditions of the flow field Euler control nodes. The mapping pairing attributes of the pressure loading and displacement feedback are digitally encoded and forcibly written into a globally unified node sequence memory space for fixed association. Under this operation, the physical state evolution data of the external explosion impact flow field and the dynamic deformation mechanics data of the internal head solid structure achieve deep feature aggregation, establishing a complete fluid-structure interaction mesh structure that transmits detonation energy and structural response. For example, at the instant the shock wave contacts the front of the head, the instantaneous pressure calculated by a flow field Euler control node at a distance of 1.0 mm is 120.0 kPa. This 120.0 kPa pressure data is immediately mapped and distributed to the corresponding scalp solid field node to generate a force response. Subsequently, the scalp node undergoes an inward displacement of 0.5 mm and a moving velocity of 2.0 m / s under pressure. This displacement and velocity variables are then immediately fed back to the Euler space, forcing the flow field mesh to update its volume shape to comply with mass conservation.

[0034] Please see Figure 5 The specific steps of S4 are as follows: S401: Extract the spatial coordinates of the Euler nodes on the outside of the fluid-structure interaction mesh structure, collect the time series data of shock wave pressure from the pressure sensor, map the pressure time series data into a unidirectional scalar field of Euler nodes based on the spatial position correspondence, aggregate the dynamic amplitude of the load and spatial coordinates, and establish a spatiotemporal distribution array of boundary loads. The spatial coordinate set of the outermost Euler nodes on the blast-facing surface of the constructed fluid-structure interaction mesh structure is extracted. Shock wave pressure time-series data collected during a live-fire free-field explosion experiment is used. This data is captured and recorded in real-time by pencil-shaped high-frequency pressure sensors distributed at a predetermined blast distance. Its characteristic feature is a numerical curve sequence where the pressure rises sharply to a peak within an extremely short microsecond time and then decays exponentially. Based on the relative orientation between the physical three-dimensional location of the sensors and the spatial coordinate system of the global model, the collected pressure time-series data is mapped and assigned values ​​to the unidirectional scalar field values ​​on the outer blast-facing Euler nodes according to time progression logic. The dynamic pressure amplitude of the shock wave load changing over time is aggregated with the spatial coordinate data of each receiving node. The pressure values ​​are bound to their timestamps and spatial horizontal and vertical depth positions in a high-dimensional array to establish a spatiotemporal distribution array of boundary loads describing the global external impact boundary. For example, time-series data of shock wave pressure at a measured free-field burst distance of 5.0 meters shows that after arriving at zero time, it reaches a peak of 150.0 kPa at 0.2 milliseconds, and then decays to the atmospheric pressure noise level at 5.0 milliseconds. The Euler node coordinates at the corresponding shock wavefront normal direct point are extracted as 0.0 mm horizontally, 0.0 mm vertically, and 100.0 mm deep. This position parameter is then bound to the 150.0 kPa pressure amplitude at 0.2 milliseconds, forming a structured array record written into the space. This process is repeated to aggregate all records for the entire wavefront scan duration.

[0035] S402: Based on the spatiotemporal distribution array of boundary loads, scan the topological boundary domain of the extreme coordinates of the fluid-structure interaction mesh structure, extract the outermost free state node set, apply unidirectional energy dissipation transmission properties to the outermost free state node set and lock the wavefront reflection scalar flux, map the node index and the constraint state association key value, and generate non-reflective physical boundary constraints. Based on the spatiotemporal distribution array of boundary loads recording the propagation state of shock waves, the extreme coordinate topological boundary domain of the fluid-structure interaction mesh structure is scanned, i.e., the boundary coordinates of the six outermost planes in the horizontal, vertical, and depth directions of the search space are extracted, and the set of outermost free-state nodes on these six planes is extracted. A unidirectional energy dissipation and transmission property is applied to this set of outermost free-state nodes, the core action of which is to calculate the acoustic impedance parameter of the outermost air element. The density parameter of air under standard atmospheric conditions is extracted as 1.225 kg / m³, and the sound velocity parameter in air is extracted as 340.0 m / s. The density parameter and the sound velocity parameter are multiplied to obtain the acoustic impedance parameter of 416.5 Reichs. This impedance parameter is used as a penalty force coefficient and assigned to the outermost node. When a pressure shock wave propagating from the interior collides with the boundary node, an equal wavefront energy transfer is calculated based on this impedance coefficient, and this energy is directly removed from the computational domain memory, thereby forcibly locking and eliminating the pressure scalar flux reflected into the internal computational domain. The external node index number, which is assigned a dissipation attribute, is recorded and bound one-to-one with a specific state-related key representing acoustic impedance constraints, generating a complete non-reflective physical boundary constraint in the configuration. The advantage of this computational logic is that by generating acoustic impedance through the product of external air density and sound speed and assigning node energy removal attributes, it successfully simulates the energy dissipation characteristics of an infinite open space on a finite three-dimensional computational domain truncated surface.

[0036] S403: Call the non-reflective physical boundary constraints, read the set of physiological tissue material property parameters, lock the solid entity unit nodes inside the fluid-structure interaction mesh structure to perform material parameter assignment operations, bind the density mechanical modulus and node sequence index, integrate the global node position, external boundary constraints and internal field physical property variables, and construct the initial finite element model of craniocerebral blast shock. The pre-defined non-reflective physical boundary constraint configuration is invoked, and the set of physiological tissue material property parameters describing the biomechanical behavior of each layer of the head is read. Within the fluid-structure interaction mesh structure, based on the pre-segmented anatomical region labels, the corresponding internal solid entity unit nodes are locked, and specific material parameter assignment operations are performed. For unit nodes labeled as hard tissue skull, the bone density value of 1900.0 kg / m³ and the mechanical modulus value describing elastic deformation resistance of 15.0 GPa are bound to the node sequence index of that region, and its Poisson's ratio value is set to 0.22 to characterize the lateral strain characteristics of hard materials; for unit nodes labeled as soft tissue brain parenchyma, its tissue density value of 1040.0 kg / m³ and the hyperelastic mechanical modulus value of 10.0 kPa are bound to the corresponding node sequence index, and its Poisson's ratio value is set to be close to the incompressible limit of 0.49. The global node position coordinates, externally set non-reflective Eulerian boundary constraints, and internal field physical property variables of each anatomical region are integrated and written into a general explicit finite element input instruction set to construct an initial finite element model of craniocerebral blast for final dynamic numerical calculation. Table 3 shows the specific settings for binding material properties of different physical tissues.

[0037] Table 3: Configuration of Mechanical Properties of Physiological Tissue Materials

[0038] As shown in Table 3, the density, modulus and Poisson's ratio parameters of different tissue regions vary greatly, reflecting the strong non-uniform biomechanical characteristics of the head structure.

[0039] Please see Figure 6 The specific steps of S5 are as follows: S501: Call the non-reflective physical boundary constraints and the initial finite element model of craniocerebral blast shock, perform time domain discretization, perform algebraic equation calculations of displacement strain matrix for each time step, and converge the three-dimensional spatial component vectors of the vertices of the solid elements inside the time node to obtain the stress time series tensor of the model node. A well-configured, non-reflective physical boundary constraint and an initial finite element model of craniocerebral blast shock are invoked and submitted to the explicit dynamics calculation kernel for time-domain discretization. Based on the explicit computational stability criteria set by mesh size and material sound velocity, the total 10.0 milliseconds of explosion impact time is divided into tens of thousands of tiny time steps, with a single time step size of 1.0 microsecond to ensure high-frequency impact response capture. For each discrete time step, the algebraic equations of the central differential displacement-strain matrix are performed on all solid element nodes. First, the resultant force vector of the external impact load and the elastic restoring force of the internal element at the current microsecond time step is read. The resultant force vector is divided by the corresponding mass value bound to the node to calculate the node acceleration vector at that moment. Then, this acceleration vector is multiplied by the 1.0 microsecond time step and added to the velocity vector of the previous moment to update the velocity vector of the current time step. Finally, the updated velocity vector is multiplied by the 1.0 microsecond time step and accumulated onto the historical displacement coordinates to obtain the current node displacement spatial components. Based on the displacement of the node, the three-dimensional normal strain and shear strain caused by the volume change of adjacent elements are obtained, and these are substituted into the material constitutive deviatoric stress equation to calculate a new stress tensor. The three-dimensional spatial stress component vectors of all vertices in the transverse, longitudinal, and depth directions of each solid element within the same time node are aggregated to obtain the model node stress time-series tensor reflecting the local stress state. For example, at the 5000th time step (i.e., 5.0 milliseconds), a brain parenchyma node in the frontal lobe experiences a net force divided by its infinitesimal mass, resulting in an acceleration of 1500.0 m / s². Given that the node's velocity at the previous time step was 1.1985 m / s, multiplying the current acceleration of 1500.0 m / s² by a step size of 1.0 microseconds yields a velocity increment of 0.0015 m / s². Adding this to the previous velocity of 1.1985 m / s², the current node velocity is calculated to be 1.2 m / s². Subsequently, the current velocity of 1.2 m / s is multiplied again by a step size of 1.0 microseconds, and the newly generated displacement within this tiny time step is calculated to be 0.0012 mm. Based on this local deformation increment of 0.0012 mm and its spatial topological span, a shear strain of 0.015 is derived for this element. Substituting this into the set modulus of 10.0 kPa, the shear stress at this moment is finally obtained as 0.15 kPa.

[0040] S502: Perform time-point registration and alignment operations on the model nodal stress time series tensor and the experimental standard stress test sequence parameters, calculate the absolute value of the stress scalar difference at the corresponding time point, add the absolute values ​​of the difference and perform time-domain cross-node integration to generate a set of nodal stress error scalars. The model node stress time series tensor, continuously output by the explicit computation kernel, is extracted, and the experimental standard stress test sequence parameter stored in the external validation database is retrieved. This parameter originates from the measured stress history of the pressure sensor implanted inside the standard head model within the physical detonation shock tube. Time point registration and alignment are performed on the two sets of time-dimensional data sequences, selecting a complete comparison period from 0.0 ms when the shock wave contacts the leading edge of the head model to 10.0 ms when the shock wave completely decays. At each discrete analysis time point of the alignment, the transient stress value of the model node calculated by the simulation model at the corresponding spatial coordinates of the specific intracranial sensor is extracted, and the calibration stress value corresponding to the experimental standard stress sequence at the same physical time is also extracted. The transient stress value of the model node is subtracted from the experimental calibration stress value, and the stress scalar difference between the two at the corresponding time points is calculated. Then, the absolute value operation is performed on this difference to eliminate the error cancellation effect caused by positive and negative fluctuations. The absolute values ​​of the scalar differences at all discrete time points within the period are continuously added together, and cross-node integration and summation are performed across the entire time domain to generate a scalar set of nodal stress errors that covers the cumulative deviations across different sensor locations and the entire process. For example, at the brainstem region comparison point, the model-calculated stress at 1.2 ms is 45.0 kPa, while the experimental standard test stress is 40.0 kPa. Subtracting the two gives a difference of 5.0 kPa, the absolute value of which is 5.0 kPa. At 1.3 ms, the model-calculated stress is 38.0 kPa, while the experimental value is 42.0 kPa. Subtracting the two gives a difference of -4.0 kPa, the absolute value of which is 4.0 kPa. These 5.0 and 4.0 are added together and integrated, and this process is repeated to calculate the total error at all times from 0.0 to 10.0 ms, ultimately obtaining a value of 185.0 kPa representing the overall degree of deviation, which is recorded in the error scalar set.

[0041] S503: Based on the nodal stress error scalar set, compare the error scalar with the preset convergence threshold, remove the corresponding dynamic incremental step sequence that exceeds the preset convergence threshold, extract the associated geometric coordinates and deformation attributes that do not exceed the preset convergence threshold, perform mesh freezing, and establish a finite element model of the target craniocerebral blast injury. Based on the generated set of nodal stress error scalars, the error accumulation in each local region is evaluated in depth. According to the tolerance statistical distribution established in large-scale basic cranial experiments, the upper limit of the acceptable overall fixed numerical range is set to a total integral of 150.0 kPa, and this 150.0 kPa is defined as the preset convergence threshold used for evaluation. Furthermore, the accumulated errors at each time point in the nodal stress error scalar set are arranged chronologically, and the absolute error values ​​at two adjacent time points are extracted, subtracted, and then divided by the time step to calculate the rate of change of the difference. The obtained rate of change of the difference for each time period is compared one by one with another preset convergence threshold of 5.0 kPa per millisecond. If the rate of change of the difference is greater than 5.0 kPa per millisecond or the total integral of the local error exceeds 150.0 kPa, the mechanical calculation in that region is considered unstable, and the dynamic incremental step sequence corresponding to the convergence threshold and its subsequent calculation data are immediately discarded. Only when consecutive time nodes with a difference rate of change less than or equal to a preset convergence threshold of 5.0 kPa per millisecond are observed will their corresponding dynamic incremental step sequences be retained as valid analysis results. For the retained time node data, the associated geometric coordinates and deformation attributes that do not exceed the preset convergence threshold are extracted; that is, the 3D coordinate values ​​of each vertex of the effective solid element at the same time node and their corresponding spatial displacement vector values ​​are extracted and combined. Subsequently, this combined data undergoes a consistency check of the displacement values ​​of adjacent elements sharing vertices, i.e., it is determined whether the absolute difference of the displacement vectors of adjacent elements at common nodes is zero. After screening, a set of solid elements that do not exhibit coordinate tearing at common nodes and satisfy the spatial continuity law is selected. Based on this selected set of high-quality, non-penetrating solid elements, structural solidification is implemented, and a mesh freeze command is executed to block all subsequent coordinate updates. Finally, a stable and usable finite element model of the target craniocerebral blast injury is established and output. For example, if the rate of change of the difference between adjacent time nodes within a certain period is 3.2 kPa per millisecond, it is compared with the preset convergence threshold of 5.0 kPa per millisecond and determined to be less than the threshold. Therefore, the incremental step sequence is retained. The original coordinates of the vertex at that moment and the displacement vector data combination are extracted. The displacement components of two adjacent brain parenchymal units at the shared node are both 0.5 mm laterally and 0.2 mm longitudinally. The difference is zero, which satisfies the consistency. The set of physical units is directly frozen and solidified.

[0042] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in the present invention should be included within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of protection of the described technical solutions.

Claims

1. A method for establishing a finite element model of craniocerebral blast injury, characterized in that, Includes the following steps: S1: Acquire tomographic images of healthy cranial soft tissue and hard tissue from MRI and CT scanners, input them into the image registration and fusion algorithm to perform grayscale feature calculation and cross-sectional contour extraction, extract the spatial coordinates of the nodes of the two-dimensional contour boundary, and perform inter-layer normal interpolation on the spatial coordinates to generate a three-dimensional geometric solid model. S2: Perform tetrahedral meshing on the three-dimensional geometric solid model, extract the coordinates of the mesh unit nodes and calculate the unit volume deformation factor, and adjust the spatial displacement coordinates of nodes below the preset deformation threshold towards the geometric center of the unit to generate the head solid mesh model. S3: Extract the coordinates of the nodes on the outer surface of the head solid mesh model, divide the external space region into a hexahedral air Eulerian element structure, calculate the spatial distance between the coordinates of the internal element nodes and the outer surface nodes, and perform bidirectional contact node pairing operation to establish a fluid-structure interaction mesh structure. S4: Apply the shock wave pressure time series data collected by the pressure sensor in the experimental environment to the Euler outer boundary node of the fluid-structure interaction mesh structure, set the non-reflective physical boundary constraint, obtain the tissue material property data and assign it to the internal solid element node of the fluid-structure interaction mesh structure to construct the initial finite element model of craniocerebral blast shock. S5: Based on the non-reflective physical boundary constraints, the initial finite element model of craniocerebral blast injury is input into the explicit dynamic solver, and step difference operation is performed to extract the stress values ​​of the model nodes; the stress values ​​of the model nodes are compared with the experimental stress test values ​​to perform error comparison operation, and the state is solidified when the difference is lower than the preset convergence threshold to generate the target craniocerebral blast injury finite element model.

2. The method for establishing a finite element model of craniocerebral blast injury according to claim 1, characterized in that, The three-dimensional geometric solid model includes the three-dimensional morphology of the skull, the anatomical configuration of the ventricles, and the surface contour of the scalp; the head solid mesh model includes a set of tetrahedral elements, mesh node coordinates, and element connection topology matrix; the fluid-structure interaction mesh structure includes an Eulerian computational domain, a Lagrangian solid domain, and a fluid-structure interaction interface; the initial finite element model of craniocerebral blast injury includes a hyperelastic constitutive model, bulk modulus, and state equation parameters; the target craniocerebral blast injury finite element model includes an experimentally calibrated mesh topology, corrected material constitutive parameters, and contact penalty function constraints.

3. The method for establishing a finite element model of craniocerebral blast injury according to claim 1, characterized in that, The specific steps of S1 are as follows: S101: Acquire soft tissue tomographic images and hard tissue tomographic images using MRI and tomography equipment, perform three-dimensional translation and alignment, extract the single-channel grayscale components corresponding to the overlapping areas, multiply the components with the same coordinates and take the absolute value as the target grayscale feature, and map them one by one into matrix elements according to the spatial voxel coordinates corresponding to the target grayscale feature to establish a soft and hard tissue feature matrix. S102: Call the soft and hard tissue feature matrix, calculate the second derivative of the pixel gray level gradient and filter edge points that exceed the preset mutation threshold, perform spatial connected component labeling and aggregate pixels with the same label, construct a set of closed curves, extract discrete node coordinate sequences from the set of closed curves, and obtain a two-dimensional contour node coordinate set. S103: Based on the coordinate set of the two-dimensional contour nodes, calculate the normal vector of the corresponding node of the adjacent fault plane, determine whether the cosine value of the included angle is within the preset parallel judgment interval, if it is, use the corresponding node as the matching node for search, construct the cross-layer three-dimensional related line segment, perform linear geometric interpolation along the related line segment to supplement the support point, and generate a three-dimensional geometric solid model.

4. The method for establishing a finite element model of craniocerebral blast injury according to claim 3, characterized in that, The specific steps of S2 are as follows: S201: Obtain the three-dimensional geometric entity model, implant discrete seed nodes along the surface and interior of the entity space boundary, connect adjacent nodes to construct a tetrahedral structure according to the spatial topology rules, extract the absolute position vectors of the grid unit vertices and splice them into a numerical arrangement, and establish a three-dimensional coordinate matrix of the grid unit. S202: Based on the three-dimensional coordinate matrix of the mesh unit, read the absolute position vectors of the four vertices, calculate the absolute value of the mixed product of adjacent common vertex edge vectors to obtain the actual unit volume, extract the standard reference volume of the same circumscribed sphere tetrahedron, calculate the ratio of the two types of volumes and record it as the distortion state parameter, and obtain the set of unit volume deformation factors. S203: Call the set of unit volume deformation factors, extract the target deformation nodes whose distortion state parameters are lower than the preset deformation threshold, perform an arithmetic mean operation on the absolute position vectors of the four vertices of the unit to extract the geometric center coordinates, construct the displacement line segment of the node pointing to the geometric center coordinates and adjust the coordinates to generate the head solid mesh model.

5. The method for establishing a finite element model of craniocerebral blast injury according to claim 4, characterized in that, The preset deformation threshold is determined based on the statistical distribution of the unit volume deformation factor set; the preset deformation threshold is obtained by sorting the unit volume deformation factor set and selecting the distortion state parameter value corresponding to the preset quantile position.

6. The method for establishing a finite element model of craniocerebral blast injury according to claim 4, characterized in that, The specific steps for S3 are as follows: S301: Obtain the head entity mesh model, extract the discrete coordinate node set along the outer boundary of the head entity mesh model, construct an orthogonal hexahedral air Euler element array in the external space according to a preset scale, perform topological relative position alignment determination on the discrete coordinate node set and the element array, and generate the boundary interaction initial coordinate matrix. S302: Call the initial coordinate matrix of the boundary interaction, traverse the coordinate nodes of the unit nodes and the outer surface coordinate nodes of the head entity mesh model, quantify the node span according to the Euclidean distance calculation rule to extract the distance value, compare the distance value with the preset contact tolerance limit, filter the candidate associated node set and remove isolated nodes, and generate a near-end potential contact topology map. S303: Based on the near-end potential contact topology map, perform bidirectional data interaction mapping pairing operation on the candidate associated node set, between the flow field Euler control nodes and the coordinate nodes on the outer surface of the head entity mesh model, force the mapping pairing attributes into the global node sequence space, aggregate the internal and external field associated data elements, and establish a fluid-structure interaction mesh structure.

7. The method for establishing a finite element model of craniocerebral blast injury according to claim 6, characterized in that, The specific steps of S4 are as follows: S401: Extract the spatial coordinates of the Euler nodes on the outer side of the fluid-structure interaction mesh structure, collect the time series data of shock wave pressure from the pressure sensor, map the pressure time series data into a unidirectional scalar field of Euler nodes based on the spatial position correspondence, aggregate the dynamic amplitude of the load and the spatial coordinates, and establish a spatiotemporal distribution array of the boundary load. S402: Based on the boundary load spatiotemporal distribution array, scan the extreme coordinate topological boundary domain of the fluid-structure interaction mesh structure, extract the outermost free state node set, apply unidirectional energy dissipation transmission property to the outermost free state node set and lock the wavefront reflection scalar flux, map the node index and constraint state association key value, and generate non-reflective physical boundary constraints. S403: Invoke the non-reflective physical boundary constraints, read the set of physiological tissue material property parameters, lock the solid entity unit nodes inside the fluid-structure interaction mesh structure and perform material parameter assignment operations, bind the density mechanical modulus and node sequence index, integrate the global node position, external boundary constraints and internal field physical property variables, and construct the initial finite element model of craniocerebral blast shock.

8. The method for establishing a finite element model of craniocerebral blast injury according to claim 7, characterized in that, The specific steps of S5 are as follows: S501: Call the non-reflective physical boundary constraints and the initial finite element model of craniocerebral blast shock, perform time domain discretization, perform algebraic equation calculations of displacement strain matrix for each time step, and converge the three-dimensional spatial component vectors of the vertices of the internal solid elements of the time node to obtain the stress time series tensor of the model node. S502: Perform time point registration and alignment operation on the model node stress time series tensor and the experimental standard stress test sequence parameters, calculate the absolute value of the stress scalar difference at the corresponding time point, add the absolute values ​​of the difference and perform time domain cross-node integration operation to generate a set of node stress error scalars. S503: Based on the set of nodal stress error scalars, compare the error scalars with the preset convergence threshold, remove the corresponding dynamic incremental step sequences that exceed the preset convergence threshold, extract the associated geometric coordinates and deformation attributes that do not exceed the preset convergence threshold, perform mesh freezing, and establish a finite element model of the target craniocerebral blast injury.

9. The method for establishing a finite element model of craniocerebral blast injury according to claim 8, characterized in that, The preset convergence threshold is based on the upper limit of a fixed numerical range determined by the statistical distribution of the nodal stress error scalar set; The removal of the corresponding dynamic incremental step sequence that exceeds the preset convergence threshold means that the node stress error scalar set is arranged in chronological order and the difference change rate is calculated for adjacent time nodes. The difference change rate is compared with the preset convergence threshold one by one. The dynamic incremental step sequence corresponding to the consecutive time nodes that satisfy the difference change rate is less than or equal to the preset convergence threshold is retained. The associated geometric coordinates and deformation attributes that do not exceed the preset convergence threshold refer to the combination data of the three-dimensional coordinates and displacement vectors of the corresponding entity unit vertices at the same time node. Consistency verification processing is performed on the combination data to filter the set of entity units that meet the spatial continuity.

Images