A method for simulating the mechanical behavior of a spine based on a multi-layer material structure

Abstract

The application discloses a kind of simulation methods of spinal mechanics behavior based on multilayer material structure, it is related to spinal biomechanics simulation technical field. Including: obtaining the medical image data of target spinal segment and constructing initial three-dimensional geometric model;The initial three-dimensional geometric model is carried out multilayer material region segmentation, and the spinal bone tissue layer, annulus fibrosus layer, intervertebral disc nucleus layer, ligament attachment layer and nerve root sheath layer are obtained;Different viscoelastic constitutive parameters are respectively given to form multilayer material composite structure;Pre-set physiological load boundary condition is applied, and nonlinear finite element iterative calculation is executed, and stress-strain response field is obtained;According to stress-strain response field, the interface stress concentration area of each material layer is reversely identified, and the interface stress concentration area is marked as mechanical behavior simulation output point.The method realizes the fine simulation of the mechanical behavior of spinal segment by multilayer material layering and interface stress concentration reverse identification.

Description

Technical Field

[0001] This invention relates to the field of spinal biomechanics simulation technology, specifically a method for simulating spinal mechanical behavior based on a multilayer material structure. Background Technology

[0002] Simulation of spinal biomechanics has significant applications in implant design, pathological research of degenerative diseases, and surgical plan evaluation. Current techniques mostly simplify spinal segments into finite element modeling using two or three homogeneous materials: vertebrae and intervertebral discs. Ligamentary structures are often replaced by line or spring elements, while delicate soft tissues such as nerve root sheaths are typically ignored. This simplified modeling approach fails to reflect the influence of non-uniform density distribution within the vertebrae, the orientation dependence of collagen fibers in the annulus fibrosus of the intervertebral disc, the nonlinear characteristics of ligament pre-stretching, and the transverse anisotropy of the nerve root sheath on overall biomechanical behavior. The stress transfer mechanism at the interface between multi-layered materials is closely related to the initiation location of internal damage. Current techniques only output displacement and stress contour maps after simulation, requiring manual experience to locate stress concentration points region by region within a large dataset. There is a lack of automatic reverse identification methods for stress concentration areas at multi-layered material interfaces, resulting in insufficient clinical diagnostic guidance from simulation results and a tendency to miss deep, minute stress anomalies. Therefore, it is necessary to solve the problems of how to simultaneously construct a multi-layered composite material structure containing vertebrae, intervertebral disc annulus fibrosus, intervertebral disc nucleus pulposus, ligament attachment layer and nerve root sheath layer in simulation and assign differentiated constitutive parameters that conform to the actual mechanical characteristics of the tissue, and how to automatically and accurately identify and calibrate the stress concentration area at the interface of the multi-layered material from the stress-strain response field. Summary of the Invention

[0003] This invention provides a method for simulating the mechanical behavior of the spine based on a multilayer material structure. The method aims to segment the vertebral bone tissue layer, intervertebral disc annulus fibrosus layer, intervertebral disc nucleus pulposus layer, ligament attachment layer, and nerve root sheath layer from medical imaging data, and assign tissue-specific viscoelastic constitutive parameters to each layer. This enables precise modeling of the multilayer material composite structure of the spinal segments. Furthermore, after nonlinear finite element load response calculation, the method automatically identifies stress concentration regions at the interfaces between layers based on the stress-strain response field, and marks these stress concentration regions as output points for the mechanical behavior simulation. This improves the biofidelity of the spinal mechanical simulation and the efficiency of result analysis.

[0004] To achieve the above objectives, the present invention provides the following technical solution: The present invention provides a method for simulating the mechanical behavior of the spine based on a multilayer material structure, comprising: Acquire medical imaging data of the target spinal segment, and construct an initial three-dimensional geometric model based on the medical imaging data; The initial three-dimensional geometric model is segmented into multiple material regions to obtain the vertebral bone tissue layer, the intervertebral disc annulus fibrosus layer, the intervertebral disc nucleus pulposus layer, the ligament attachment layer, and the nerve root sheath layer; Different viscoelastic constitutive parameters are assigned to the vertebral bone tissue layer, intervertebral disc annulus fibrosus layer, intervertebral disc nucleus pulposus layer, ligament attachment layer and nerve root sheath layer respectively to form a multi-layer material composite structure; A preset physiological load boundary condition is applied to the multilayer material composite structure, and a nonlinear finite element iterative calculation is performed under the physiological load boundary condition to obtain the stress-strain response field of the multilayer material composite structure under load. Based on the stress-strain response field, the stress concentration regions at the interfaces of each material layer are identified in reverse, and these stress concentration regions are marked as the output points for mechanical behavior simulation.

[0005] Preferably, the multi-layer material region segmentation specifically includes: acquiring the grayscale gradient distribution map of the initial three-dimensional geometric model, and determining the initial interface between the vertebral bone tissue region and the intervertebral disc region based on the grayscale gradient; extracting the mean grayscale value of the central region and the grayscale gradient change rate of the edge region within the intervertebral disc region, thereby dividing the intervertebral disc region into the annulus fibrosus layer and the nucleus pulposus layer; acquiring the attachment coordinates of the ligament origin and termination points, and peeling the ligament attachment layer from the surface of the vertebral bone tissue layer; and expanding outward by a predetermined radial distance along the internal contour line of the intervertebral foramen to obtain the nerve root sheath layer. Through this segmentation process, a multi-layer geometric model of the spine that is highly consistent with the actual anatomical layers can be obtained without manual intervention.

[0006] As a preferred embodiment of the present invention, the process of assigning viscoelastic constitutive parameters to each material layer includes: obtaining an apparent density distribution cloud map of the vertebral bone tissue layer and assigning anisotropic elastic modulus to each local region; obtaining a collagen fiber orientation distribution map of the intervertebral disc annulus fibrosus layer and assigning position-dependent tensile stiffening parameters to different radial positions; obtaining the fluid volume fraction of the intervertebral disc nucleus pulposus layer and assigning incompressible fluid parameters; obtaining the initial fiber pre-stretch rate of the ligament attachment layer and assigning nonlinear stress-strain curve parameters; and obtaining the axonal direction vector of the nerve root sheath layer and assigning transversely isotropic viscoelastic parameters. The resulting multilayer composite structure can realistically reflect the differentiated mechanical responses of the bony components, annulus fibrosus, nucleus pulposus, ligaments, and nerve tissues in the spinal column.

[0007] Preferably, the process of applying physiological load boundary conditions and performing nonlinear finite element iterative calculation includes: determining the node set of the superior vertebral endplate by screening the normal vector direction and axial position of the patch on the outer surface of the vertebral bone tissue layer, and setting the determined node set as the load application surface; gradually increasing the compressive load value on the load application surface, while applying a tension boundary condition dynamically updated with the ligament elongation at the attachment point of the ligament attachment layer. This tension is calculated based on the nonlinear stress-strain curve of the ligament attachment layer and the real-time distance change of the attachment point; establishing a fiber-reinforced hyperelastic constitutive update algorithm based on the local coordinate system of the intervertebral disc annulus fibrosus layer, updating the fiber direction stress and the hydrostatic pressure of the intervertebral disc nucleus pulposus layer in each increment step; monitoring the normalized error of the displacement increment of all nodes between two adjacent increment steps, terminating the iteration when the error value is lower than the convergence tolerance, and outputting the stress-strain response field. This scheme considers the load transfer of the superior vertebral body, the ligament dynamic tension, and the coupling effect of fiber-reinforced hyperelasticity and incompressibility, which can improve the prediction accuracy of spinal biomechanical behavior simulation under physiological load while ensuring convergence stability.

[0008] A further preferred method for identifying stress concentration regions at the interface includes: acquiring the equivalent stress values ​​of the nodes at the interface between the vertebral bone tissue layer and the intervertebral disc annulus fibrosus layer in the stress-strain response field; clustering nodes with equivalent stress values ​​exceeding a preset stress threshold into connected domains to form multiple connected clusters of stress concentration regions; for each connected cluster, using the product of the average equivalent stress value of the nodes within the cluster and the number of nodes as the stress concentration weight coefficient; using the geometric center coordinates of the connected cluster with the largest stress concentration weight coefficient as the output point of the mechanical behavior simulation, and outputting a list of node numbers occupied by the connected cluster. This method can automatically extract the most dangerous stress concentration sites from complex stress fields, providing an objective and quantitative reference for assessing the risk of intervertebral degeneration, implant failure, and developing individualized spinal surgery plans.

[0009] The technical effects and advantages provided by the present invention in the above technical solution are as follows: By performing grayscale gradient analysis, identifying grayscale differences between the center and edge of the intervertebral disc, extracting the attachment coordinates of the ligament origin and termination points, and radially expanding the internal contour of the intervertebral foramen on the initial three-dimensional geometric model, five material regions were segmented: the vertebral bone tissue layer, the intervertebral disc annulus fibrosus layer, the intervertebral disc nucleus pulposus layer, the ligament attachment layer, and the nerve root sheath layer. Anisotropic elastic modulus was assigned to the vertebral bone tissue layer based on the apparent density distribution cloud map, reflecting the changes in the density of the cancellous bone within the vertebra. Tensile stiffening parameters varying with position were assigned to the intervertebral disc annulus fibrosus layer based on the collagen fiber orientation distribution map, reproducing the gradual stiffness variation characteristics of different radial regions of the annulus fibrosus. Incompressible fluid parameters were assigned to the intervertebral disc nucleus pulposus layer based on the fluid volume fraction, maintaining the volume conservation behavior of the nucleus pulposus under load. Nonlinear stress-strain curve parameters were assigned to the ligament attachment layer based on the initial fiber pre-stretch rate, simulating the pre-tension state of the ligament before physiological stretching. Transversely isotropic viscoelastic parameters were assigned to the nerve root sheath layer based on the axonal direction vector, reflecting the mechanically dominant direction of the nerve tissue along the axonal direction. The above-mentioned material partitioning and parameter assignment methods construct a multi-layered composite material structure that is closer to physiological reality. The stress-strain response field obtained by nonlinear finite element calculation can simultaneously present the mechanical coupling effect between vertebrae, intervertebral discs, ligaments, and nerve root sheaths, avoiding the distortion of load transmission channels caused by simplified modeling, and making the simulation results more refined in depicting the mechanical behaviors such as intervertebral disc degeneration, ligament overload, and nerve root compression.

[0010] After obtaining the stress-strain response field, the stress tensor field at the interface of each material layer is extracted. The Mises equivalent stress value is calculated node by node, and threshold filtering and connected component clustering are performed to form multiple stress concentration clusters. The product of the average stress value and the number of nodes in each cluster is calculated as the stress concentration weight coefficient. The geometric center coordinates of the cluster with the largest weight coefficient are used as the mechanical behavior simulation output point, along with a list of node numbers for that cluster. This process eliminates the need for manual traversal and examination of stress contour maps, automatically capturing the most dangerous interface stress concentration locations from massive amounts of simulation data and directly converting the simulation results into evaluation output points with clear spatial coordinates. This automatic reverse identification method reduces the risk of omissions during manual interpretation, allowing the simulation output to directly point to potential structural failure initiation areas, providing objective and reproducible reference location information for surgical planning and implant placement assessment. Attached Figure Description

[0011] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments recorded in this invention. For those skilled in the art, other drawings can be obtained based on these drawings.

[0012] Figure 1This is a flowchart of a simulation method for spinal mechanical behavior based on multilayer material structures. Figure 2 This is a flowchart of multi-layer material region segmentation for a three-dimensional geometric model of spinal segments; Figure 3 This is a schematic diagram illustrating the process of assigning viscoelastic parameters to multilayer composite structures. Figure 4 This is a schematic diagram of the radial gray-level gradient distribution and preset gradient threshold in the intervertebral disc region; Figure 5 It is the nonlinear finite element iterative displacement convergence curve of a multilayer material composite structure; Figure 6 It is the stress-strain curve of ligament fibers; Figure 7 It is the Mises equivalent stress distribution curve at the interface between the vertebral bone tissue layer and the intervertebral disc annulus fibrosus layer. Detailed Implementation

[0013] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, 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, 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.

[0014] See Figure 1 This invention provides a method for simulating the mechanical behavior of the spine based on a multilayer material structure. The method includes acquiring medical imaging data of the target spinal segment and constructing an initial three-dimensional geometric model based on the medical imaging data; segmenting the initial three-dimensional geometric model into multilayer material regions to obtain a vertebral bone tissue layer, annulus fibrosus layer of the intervertebral disc, nucleus pulposus layer of the intervertebral disc, ligament attachment layer, and nerve root sheath layer; assigning different viscoelastic constitutive parameters to the vertebral bone tissue layer, annulus fibrosus layer of the intervertebral disc, nucleus pulposus layer of the intervertebral disc, ligament attachment layer, and nerve root sheath layer respectively to form a multilayer material composite structure; applying preset physiological load boundary conditions to the multilayer material composite structure and performing nonlinear finite element iterative calculations under these physiological load boundary conditions to obtain the stress-strain response field of the multilayer material composite structure under load; and identifying the stress concentration regions at the interfaces of each material layer in the multilayer material composite structure based on the stress-strain response field, marking these stress concentration regions as mechanical behavior simulation output points.

[0015] Example 1: In specific implementation, please refer to Figure 2After acquiring medical imaging data of the target spinal segment, an initial three-dimensional geometric model is constructed based on the medical imaging data. The process of multi-layer material region segmentation of the initial three-dimensional geometric model includes: acquiring the gray-level gradient distribution map of the initial three-dimensional geometric model, and determining the initial interface between the vertebral bone tissue region and the intervertebral disc region based on the gray-level gradient distribution map; extracting the mean gray-level value of the central region and the gray-level gradient change rate of the edge region within the intervertebral disc region, and dividing the intervertebral disc region into the annulus fibrosus layer and the nucleus pulposus layer based on the mean gray-level value of the central region and the gray-level gradient change rate of the edge region; acquiring the ligament attachment coordinates of the ligament origin and termination points in the initial three-dimensional geometric model, and peeling the ligament attachment layer from the surface of the vertebral bone tissue layer based on the ligament attachment coordinates of the ligament origin and termination points; and expanding outward by a predetermined radial distance along the internal contour line of the intervertebral foramen in the initial three-dimensional geometric model to obtain the nerve root sheath layer.

[0016] The specific method for obtaining the grayscale gradient distribution map of the initial 3D geometric model is as follows: The grayscale gradient vector is calculated point-by-point for the voxel data corresponding to the initial 3D geometric model. The grayscale gradient vector is calculated using the central difference method, that is, the grayscale difference between adjacent voxels is calculated along three orthogonal directions in the 3D image, and the differences in the three directions are combined into a grayscale gradient vector. The magnitude of the grayscale gradient vector is calculated for each voxel to obtain the grayscale gradient amplitude. The grayscale gradient amplitudes of all voxels constitute the grayscale gradient distribution map. When determining the initial interface between the vertebral bone tissue region and the intervertebral disc region based on the grayscale gradient distribution map, voxels with grayscale gradient amplitudes exceeding a preset gradient threshold are marked as boundary candidate voxels. The preset gradient threshold is set by analyzing the grayscale difference characteristics of vertebral bone tissue and intervertebral disc tissue in medical images. Specifically, known vertebral bone region and intervertebral disc region samples are selected, and the range of grayscale gradient amplitudes at the interface between the two types of tissues is statistically analyzed. Ninety percent of the lower limit of the range is taken as the preset gradient threshold. Connectivity analysis was performed on candidate boundary voxels, and the largest connected surface was retained as the initial interface between the vertebral bone tissue region and the intervertebral disc region. The extraction of the largest connected surface was achieved through a connected component labeling algorithm, that is, the candidate boundary voxels were labeled with six-neighbor or twenty-six-neighbor connected components, and the connected component containing the most voxels was selected as the initial interface.

[0017] The method for extracting the mean gray level of the central region and the gray level gradient change rate of the edge region within the intervertebral disc region is as follows: Within the intervertebral disc region, a central sampling sphere is constructed with the geometric center of the intervertebral disc region as the center and a radius of 15% of the equivalent radius of the intervertebral disc region. The average gray level of all voxels within the central sampling sphere is calculated, and the resulting value is recorded as the mean gray level of the central region. Simultaneously, within the intervertebral disc region, an annular sampling band is constructed with the geometric center of the intervertebral disc region as the axis, and the distance between the band and the outer contour of the intervertebral disc region is two voxel dimensions. The rate of change of the gray level gradient amplitude along the radial direction is calculated within the annular sampling band. The average gray level gradient amplitude of all voxels within the annular sampling band is then divided by the radial width of the annular sampling band to obtain the gray level gradient change rate of the edge region. When dividing the intervertebral disc region into the annulus fibrosus layer and the nucleus pulposus layer based on the mean gray value of the central region and the gray gradient change rate of the peripheral region, if the gray value of a certain local area is lower than the product of the mean gray value of the central region and a preset coefficient, and the gray gradient amplitude of the corresponding local area is greater than the product of the gray gradient change rate of the peripheral region and a preset multiple, then this local area is marked as belonging to the annulus fibrosus layer; the remaining area within the intervertebral disc region is marked as belonging to the nucleus pulposus layer. The preset coefficient ranges from 0.7 to 0.9, and the preset multiple ranges from 1.2 to 1.5. The specific values ​​of the preset coefficient and preset multiple are adjusted within the range according to the degree of intervertebral disc degeneration: for intervertebral discs without significant degeneration, the preset coefficient is 0.9 and the preset multiple is 1.2; for intervertebral discs with moderate degeneration, the preset coefficient is 0.8 and the preset multiple is 1.35; for intervertebral discs with severe degeneration, the preset coefficient is 0.7 and the preset multiple is 1.5.

[0018] To obtain the ligament origin and termination attachment coordinates in the initial 3D geometric model, the curvature abrupt change at the interface between the bone surface and soft tissue is identified. Curvature analysis is performed on the surface of the vertebral bone tissue layer, calculating the Gaussian curvature and mean curvature of each surface vertex. Surface vertices with Gaussian or mean curvature exceeding a preset curvature threshold are marked as candidate ligament attachment points. The preset curvature threshold is set according to the anatomical attachment characteristics of different ligaments, such as the anterior longitudinal ligament, posterior longitudinal ligament, and ligamentum flavum: for the anterior longitudinal ligament, the preset curvature threshold is set to 1.8 times the mean curvature of the anterior vertebral surface; for the posterior longitudinal ligament, the preset curvature threshold is set to 1.6 times the mean curvature of the posterior vertebral surface; and for the ligamentum flavum, the preset curvature threshold is set to twice the mean curvature of the lamina region. The candidate ligament attachment points are clustered according to spatial continuity and anatomical direction, with the two endpoints of each category serving as the ligament origin and termination attachment coordinates, respectively. When the ligament attachment layer is peeled from the surface of the vertebral tissue layer based on the coordinates of the ligament origin and termination points, the surface area of ​​the ligament attachment layer is generated by expanding the line connecting the coordinates of the ligament origin and termination points outwards by a predetermined ligament width. Then, this surface area is offset inwards along the normal direction of the vertebral tissue layer surface by a ligament thickness value, resulting in the three-dimensional solid structure of the ligament attachment layer. The predetermined ligament width and thickness values ​​are set based on the anatomical statistical average values ​​of each ligament.

[0019] The method for obtaining the nerve root sheath layer by extending the internal contour line of the intervertebral foramen from the initial 3D geometric model outward by a preset radial distance is as follows: First, the internal contour line of the intervertebral foramen is extracted, specifically through 2D cross-sectional image segmentation. In multiple planes perpendicular to the direction of the intervertebral foramen, a region growing algorithm is used to segment the inner boundary of the intervertebral foramen, and the boundary points of each plane are connected to form the internal contour line of the intervertebral foramen. Then, the local normal vector of each point on the internal contour line of the intervertebral foramen is calculated. The local normal vector is set as the direction from the center of the intervertebral foramen to the boundary point within the plane containing the internal contour line of the intervertebral foramen. The points on the internal contour line of the intervertebral foramen are moved outward by a preset radial distance along the local normal vector to obtain the outer contour line of the nerve root sheath layer. The preset radial distance is set as the average value of the anatomical thickness of the nerve root sheath, ranging from 0.8 mm to 1.2 mm. The nerve root sheath layer consists of the region between the internal contour line of the intervertebral foramen and the outer contour line of the nerve root sheath layer. This region is stretched along the direction of the intervertebral foramen to form the 3D geometric structure of the nerve root sheath layer.

[0020] See Figure 4The horizontal axis in the figure represents the radial position of the intervertebral disc region, in millimeters, ranging from 0 to 10 mm. The left side of the vertical axis indicates the grayscale value, and the right side indicates the grayscale gradient amplitude. The solid curve represents the spatial distribution of the grayscale value, showing a rapid decrease from approximately 250 at the radial position of 0 mm to approximately 30 at 10 mm, exhibiting a monotonically decreasing trend. The decrease is particularly steep in the initial range (0–1.5 mm), followed by a slow stabilization. The dashed curve represents the trend of the grayscale gradient amplitude, starting from a low value at approximately 10 mm, increasing rapidly along the radial direction, reaching a steady level of approximately 50 after 2 mm, and fluctuating but remaining generally stable in subsequent ranges. A horizontal dashed line is also marked in the figure, representing the preset gradient threshold, with a value of approximately 50.

[0021] This figure, combined with the description of multi-layer material region segmentation of the initial three-dimensional geometric model in Example 1, reflects the determination of the initial interface between vertebral bone tissue and intervertebral disc region, as well as the boundary between the annulus fibrosus layer and the nucleus pulposus layer of the intervertebral disc, based on the spatial distribution of grayscale values ​​and their gradient amplitudes. Specifically, a significant decrease in grayscale value corresponds to the transition from vertebral bone tissue to the intervertebral disc region, while the location where the grayscale gradient amplitude approaches or exceeds the preset gradient threshold is the candidate region for the tissue boundary. At a radial position of approximately 1.5 mm, the grayscale gradient amplitude crosses the preset gradient threshold, indicating that this is the initial interface between the vertebral bone and the intervertebral disc. Subsequently, the grayscale gradient amplitude increases and tends to stabilize, reflecting the grayscale gradient characteristics of the annulus fibrosus layer and the nucleus pulposus layer within the intervertebral disc. Combined with the determination method of the mean grayscale value of the central region and the grayscale gradient change rate of the edge region in Example 1, this provides a quantitative basis for further subdivision of the intervertebral disc region. The overall curve trend clearly reveals the technical path and the basis for setting key parameters for multi-layer material region segmentation through grayscale values ​​and their gradient distribution.

[0022] Example 2: In specific implementation, please refer to Figure 3 The process of assigning different viscoelastic constitutive parameters to the vertebral bone tissue layer, intervertebral disc annulus fibrosus layer, intervertebral disc nucleus pulposus layer, ligament attachment layer, and nerve root sheath layer to form a multilayer material composite structure is achieved through the following steps: An apparent density distribution cloud map of the vertebral bone tissue layer was obtained, and anisotropic elastic moduli were assigned to each local region of the vertebral bone tissue layer based on the apparent density distribution cloud map. The apparent density distribution cloud map was obtained by extracting the grayscale value of each voxel in the initial 3D geometric model from the medical image, converting the grayscale value to an apparent density value using a pre-calibrated correspondence between grayscale and bone mineral density, and generating a 3D density distribution field within the vertebral bone tissue layer. The calibration relationship was established by scanning a calcified hydroxyapatite phantom of known density, and linearly regressing the grayscale value of the phantom under different scanning parameters to its known density value to obtain the conversion coefficient. At each finite element integration point of the vertebral bone tissue layer, the average apparent density of the voxels covered by the current element was taken as the apparent density value of the integration point, denoted as ρ. The anisotropic elastic modulus was assigned based on the empirical relationship between the apparent density and elastic modulus of bone tissue, and the elastic modulus of the principal direction was calculated using the following formula: ; In the formula, Let ρ represent the elastic modulus along the principal axis i, where i = 1, 2, and 3 correspond to the radial, tangential, and axial directions of the vertebrae, respectively; ρ is the apparent density at the integration point, in grams per cubic centimeter; k is the proportionality coefficient, the value of which is obtained through reverse calibration. Ex vivo bone samples of the same species and anatomical location as the target spinal segment are selected, and their density distribution is obtained through micro-CT scanning. Then, point-by-point microindentation tests are performed to obtain the elastic modulus. The power-law relationship between density and elastic modulus is fitted using least-squares to obtain the value of k, which ranges from 0.5 to 3.0, in megapascals per gram per cubic centimeter raised to the power of α; α is the power-law exponent, determined through fitting of the same batch of experimental data, with a value of 2.0. Meanwhile, the anisotropic direction of the vertebral bone tissue layer is determined through structural tensor analysis: the gradient of the apparent density distribution cloud map is calculated to obtain the local density gradient vector field. At each integration point, a spherical neighborhood with a specific radius centered on the integration point is constructed. The structural tensor of the density gradient vector in this neighborhood is calculated. The structural tensor is decomposed into eigenvalues. The eigenvector corresponding to the largest eigenvalue is taken as the dominant orientation direction of the trabecular bone. The dominant orientation direction is taken as the axial material principal axis. The radial and tangential directions are set in orthogonal directions, thereby forming a local anisotropic elastic matrix.

[0023] A collagen fiber orientation distribution map of the annulus fibrosus of the intervertebral disc was obtained. Based on this map, position-dependent tensile stiffening parameters were assigned to different radial locations within the annulus fibrosus. The collagen fiber orientation distribution map was obtained by scanning the target spinal segment using diffusion tensor magnetic resonance imaging (DTI) to obtain the diffusion tensor matrix for each voxel. Eigenvalue decomposition was performed on the diffusion tensor matrix, and the eigenvector corresponding to the largest eigenvalue was used as the principal orientation direction of the collagen fibers for that voxel. The diffusion tensor magnetic field image was spatially registered with the initial 3D geometric model, and trilinear interpolation was used to map the fiber orientation vectors in the voxels to the finite element integration points of the annulus fibrosus. Inside the annulus fibrosus, a normalized radial coordinate variable *r* was established from the nucleus pulposus boundary to the outer boundary of the annulus fibrosus. *r* was set to 0 at the nucleus pulposus boundary and 1 at the outer boundary of the annulus fibrosus. For each integration point, the *r* coordinate value was obtained, and the angle deflection between the collagen fiber orientation direction and the circumferential direction of the intervertebral disc at the integration point was calculated. This deflection was used as the fiber angle parameter for that integration point. Tensile hardening parameters are assigned by defining a fiber strain energy function, which is in exponential form, where scalar parameters... and Controlling the tensile hardening response characteristics of the fiber. Parameters vary depending on the radial position r. and Assign values ​​according to the linear change rule: , ,in The value is 0.5 MPa. The value is taken as 1.5 MPa. The value is 10.0. The value is 30.0, and these values ​​are derived from the fitted average of uniaxial tensile test data of the human lumbar intervertebral disc annulus fibrosus. The tensile hardening parameter corresponding to each integration point is calculated based on its r-coordinate and written into the material property card of that integration point.

[0024] The fluid volume fraction of the nucleus pulposus layer of the intervertebral disc was obtained, and incompressible fluid parameters were assigned to the nucleus pulposus layer based on the fluid volume fraction. The fluid volume fraction was obtained by quantitatively measuring the T2 relaxation time of the nucleus pulposus region on T2-weighted magnetic resonance imaging (MRI) images. Specifically, all voxels of the nucleus pulposus layer were delineated on the T2 mapping sequence images of the MRI workstation, and the T2 value of each voxel was read. Using a pre-established conversion relationship between T2 value and water content, the T2 mapping was converted into a water content distribution, which is the fluid volume fraction. The conversion relationship between T2 value and water content was obtained by establishing a linear calibration curve between Karl Fischer titration of water content in ex vivo intervertebral disc samples and the corresponding T2 value. The arithmetic mean of the fluid volume fraction of all voxels within the nucleus pulposus layer was then calculated to obtain the overall average fluid volume fraction of the nucleus pulposus layer. When assigning incompressible fluid parameters to the nucleus pulposus layer of the intervertebral disc, the nucleus pulposus material model is set as a hybrid model of hyperelastic material and hydrostatic pressure, where the hyperelastic constitutive model is the Neo-Hookean model, and the initial shear modulus G is given by the formula... Sure, The reference shear modulus under normal moisture content is taken as 0.2 MPa. The reference fluid volume fraction was set to 0.8, and β was the water content sensitivity coefficient, obtained through fitting experiments on the dynamic shearing of the nucleus pulposus under different water content conditions, with a value of 15.0. This represents the actual measured overall average fluid volume fraction. The bulk modulus of the nucleus pulposus material was set to 500 times the initial shear modulus to ensure approximately incompressible properties.

[0025] The initial fiber pretension ratio of the ligament attachment layer is obtained, and nonlinear stress-strain curve parameters are assigned to the ligament attachment layer based on the initial fiber pretension ratio. The initial fiber pretension ratio is obtained as follows: after constructing the three-dimensional geometry of the ligament attachment layer, the spatial straight-line distance between the attachment points at both ends of the ligament attachment layer is extracted and recorded as the current attachment distance. Simultaneously, the original length of the ligament attachment layer under unloaded and relaxed conditions was obtained. This original length was obtained by averaging multiple measurements taken with vernier calipers under unloaded conditions on the ex vivo specimen, and is denoted as . Initial fiber pre-stretch ratio according to Calculate. If If it is less than 1.0, then... A forced setting of 1.0 indicates that the ligament is in its natural, unstretched state. The parameters of the nonlinear stress-strain curve are assigned by defining the normalized one-dimensional fiber stress-tension relationship of the ligament. This relationship is set as a nonlinear elastic curve, and its functional form includes three segments: a low-stiffness stage in the toe region, a linear high-stiffness stage, and a post-yield stage. The transition between these three stages is determined by two critical elongation ratios. and control, The value is 1.03. The value is taken as 1.08. In the nonlinear finite element solution process, the actual tensile ratio of the fiber is calculated at each step. ,in The additional stretching ratio caused by the external load is substituted into the nonlinear stress-strain curve function to obtain the current stress of the fiber, thereby introducing a pre-stretching effect into the material model of the ligament attachment layer.

[0026] The axonal orientation vector of the nerve root sheath is obtained, and transversely isotropic viscoelastic parameters are assigned to the nerve root sheath based on this vector. The axonal orientation vector is obtained through high-angle-resolution diffusion imaging and deterministic fiber tracing techniques. Specifically, in the diffusion imaging data, each voxel obtains a fiber orientation distribution function, and the direction of the maximum peak of the fiber orientation distribution function is extracted as the principal axonal orientation vector. The fiber orientation distribution function data is registered with the nerve root sheath in the initial 3D geometric model, and the axonal orientation vector is assigned to each integration point of the nerve root sheath through nearest-neighbor interpolation. When assigning transversely isotropic viscoelastic parameters to the nerve root sheath, the principal axonal orientation vector is set as the principal axis direction of the material, and a transversely isotropic linear viscoelastic constitutive model is established, which includes five independent elastic parameters: transverse elastic modulus. axial elastic modulus Transverse shear modulus lateral Poisson's ratio and axial Poisson's ratio It also includes two viscosity parameters: volume relaxation time constant. and bias relaxation time constant These parameters were calibrated through a combination of dynamic nanoindentation and micro-stretching experiments on the sciatic nerve, including the transverse elastic modulus. The axial elastic modulus is calibrated to 1.2 MPa. The calibrated value is 12.0 MPa, and the transverse shear modulus is... The calibration is 0.4 MPa, and the transverse Poisson's ratio is... The axial Poisson's ratio is calibrated to 0.49. Calibrated to 0.3, volume relaxation time constant Calibrated to 0.5 seconds, bias relaxation time constant. The calibration time is 0.2 seconds. The material stiffness matrix at each integration point is rotated according to its axial direction vector to ensure that the axonal direction vector remains parallel to the material axis, thereby completing the assignment of viscoelastic parameters of the nerve root sheath.

[0027] Example 3: In practical implementation, the process of applying a preset physiological load boundary condition to the multilayer material composite structure and performing nonlinear finite element iterative calculation under the physiological load boundary condition to obtain the stress-strain response field of the multilayer material composite structure under load is achieved through the following steps: To obtain the coordinates of the nodes on the outer surface of the vertebral tissue layer in a multilayer composite structure, the set of nodes located at the upper endplate of the superior vertebral body is set as the load application surface. The specific method for obtaining all triangular facets on the outer surface of the vertebral tissue layer is as follows: traverse all triangular elements located on the outer surface of the vertebral tissue layer in the finite element mesh of the multilayer composite structure. The outer surface of the vertebral tissue layer is determined by the material number of the triangular element. Elements whose material number belongs to the vertebral tissue layer and whose triangular element has only one adjacent element are marked as outer surface elements. Extract all node coordinates of the outer surface elements and the vertex connection relationships of the triangular facets. Calculate the dot product of the normal vector of each triangular facet and the global vertical upward direction vector. The global vertical upward direction vector is defined as (0,0,1). The normal vector is obtained by the cross product of the two edge vectors of the triangular facet. The edge vector is obtained by subtracting the vertex coordinates of the triangular facet. The cross product result is normalized to a unit normal vector. The dot product of the unit normal vector and the global vertical upward direction vector is then performed to obtain the dot product value. The direction cosine threshold is set to 0.8. The rationale for this choice is that the angle between the surface normal vector of the upper endplate and the vertical upward direction vector is typically between 0 and 36.9 degrees, corresponding to a direction cosine value range of 0.8 to 1.0. Therefore, setting the direction cosine threshold to 0.8 ensures complete screening of the upper endplate surface even when there are physiological curvature changes. All triangular surfaces with a dot product value greater than the direction cosine threshold are selected as candidate endplate surfaces. The coordinate values ​​of all nodes in the candidate endplate surfaces are obtained, and the x, y, and z coordinates of each node are extracted. After removing duplicate nodes from all candidate endplate surfaces, a candidate endplate node set is formed. The arithmetic mean of the z coordinates of all nodes in the candidate endplate node set is calculated to obtain the axial average. Nodes with axial coordinates greater than the axial average are identified as upper endplate nodes; the axial coordinate is the z coordinate value of the node. The index number of the upper endplate node is stored in the load application node set. The index number is an integer value representing the sequential numbering of nodes in the finite element mesh. The load application node set is stored using a one-dimensional integer array. The triangular patch region occupied by the load application node set is marked as the load application surface. The marking method is as follows: traverse all candidate final plate patches. If the index numbers of the three vertices of the triangular patch all belong to the load application node set, then add the element number of the triangular patch to the load application surface element list.

[0028] The compressive load value is gradually increased on the load application surface, while a displacement-dependent tension boundary condition is applied at the attachment points of the ligament attachment layer in the multilayer composite structure. The current load step value of the load application surface is obtained. The current load step value is the external load increment cycle index of the nonlinear finite element solution process, starting from 1 and increasing, with the maximum load step value set to 100. The single-step load increment is obtained by dividing the preset total compressive load target value of 750 Newtons by the maximum load step value of 100, and the single-step load increment is 7.5 Newtons. The current load step value is multiplied by the single-step load increment to obtain the current total compressive load. The current total compressive load is evenly distributed to each node on the load application surface. The even distribution method is as follows: count the total number of nodes in the load application node set, divide the current total compressive load by the total number of nodes to obtain the compressive load value borne by each node, and apply the compressive load value in the negative z-axis direction of each node, perpendicular to the upper end plate and downward. Simultaneously, the initial spatial coordinates of the starting and ending attachment points of the ligament attachment layer are acquired. These initial spatial coordinates are extracted and stored after the ligament attachment layer is constructed without any applied load. In each incremental step of the nonlinear finite element iterative calculation, the current distance between the starting and ending attachment points is calculated. This current distance is obtained by reading the current displacement vectors of the starting and ending attachment points. The current position coordinates of the starting and ending attachment points are the initial spatial coordinates plus the current displacement vector, and the current position coordinates of the ending attachment point are also the initial spatial coordinates plus the current displacement vector. The Euclidean distance between these two current position coordinates is then calculated as the current distance. The ligament elongation is obtained by subtracting the initial distance between the starting and ending attachment points from the current distance between them. The initial distance is calculated from the initial spatial coordinates. The current tension value is obtained by substituting the ligament elongation into the parameters of the nonlinear stress-strain curve of the ligament attachment layer. The substitution method is as follows: divide the ligament elongation by the relaxation length of the ligament attachment layer to obtain the ligament strain; substitute the ligament strain into the nonlinear stress-strain curve function to calculate the current ligament stress; and multiply the current ligament stress by the cross-sectional area of ​​the ligament attachment layer to obtain the current tension value. The cross-sectional area of ​​the ligament attachment layer is set to 50 square millimeters based on anatomical measurements. The current tension value is then applied to both the initial and final attachment points along a direction from the initial attachment point to the final attachment point. The application method is as follows: calculate the unit direction vector from the initial attachment point to the final attachment point; apply a concentrated force vector at the initial attachment point with a magnitude equal to the current tension value and a direction opposite to the unit direction vector; and apply a concentrated force vector at the final attachment point with a magnitude equal to the current tension value and a direction in the same direction as the unit direction vector.

[0029] A local coordinate system for the annulus fibrosus of the intervertebral disc is obtained. Based on this local coordinate system, a fiber-reinforced hyperelastic constitutive update algorithm is established. Within each incremental step, the fiber-direction stress of the annulus fibrosus is updated first, followed by the hydrostatic pressure of the nucleus pulposus. The radial and circumferential vectors at each integration point in the annulus fibrosus are obtained. The radial vector is obtained by subtracting the geometric center coordinates of the nucleus pulposus from the global coordinates of the integration point and then normalizing the result. The geometric center coordinates of the nucleus pulposus are obtained by calculating the arithmetic mean of the coordinates of all nodes in the nucleus pulposus. The circumferential vector is obtained by normalizing the cross product of the radial vector and the global vertical vector (0,0,1). If the radial vector is parallel to the global vertical vector, the circumferential vector is obtained by cross-product of the radial vector and the global anterior-posterior vector (0,1,0). A local coordinate system is constructed based on the radial and circumferential direction vectors. The three orthogonal basis vectors of the local coordinate system are: the first basis vector is the circumferential direction vector, the second basis vector is the radial direction vector, and the third basis vector is obtained by the cross product of the first and second basis vectors, ensuring that the local coordinate system satisfies the right-hand rule and that the three axes are pairwise orthogonal. The current tensile ratio in the fiber direction of the intervertebral disc annulus fibrosus is calculated in the local coordinate system. The fiber direction is determined by the principal fiber orientation direction at the corresponding integration point in the collagen fiber orientation distribution map. The principal fiber orientation direction is projected onto the plane formed by the first and third basis vectors of the local coordinate system. The magnitude of the projected vector is the ratio of the current fiber length to the original length, denoted as the current tensile ratio. The corresponding fiber stress increment is found based on the current tensile ratio. The lookup method is as follows: a discrete data table of fiber tensile ratio and fiber stress increment is pre-stored. The data table is generated by piecewise cubic Hermite interpolation of fiber uniaxial tensile experimental data. The interval containing the current tensile ratio is found in the data table, and the fiber stress increment is obtained through interval linear interpolation. The fiber stress increment is then superimposed on the fiber stress of the previous increment step to obtain the current fiber stress. The current fiber stress is transformed to the global coordinate system, and the global stress tensor of the intervertebral disc annulus fibrosus is updated. The transformation method is as follows: the fiber stress tensor in the local coordinate system is multiplied by the coordinate transformation matrix. The coordinate transformation matrix is ​​a 3×3 matrix formed by arranging the components of the three orthogonal basis vectors of the local coordinate system in the global coordinate system in column order, and its transpose is used in the operation. The global stress tensor is calculated using the following formula: ; In the formula, The stress tensor in the global coordinate system is a 3×3 symmetric matrix, and its unit is megapascals (MPa). R represents the fiber stress tensor in the local coordinate system, which is a 3×3 symmetric matrix, containing non-zero normal stress components only in the fiber direction, with all other components being zero; R represents the coordinate transformation matrix, which is a 3×3 matrix, with the three columns corresponding to the components of the first, second, and third basis vectors of the local coordinate system in the global coordinate system. This represents the transpose of the coordinate transformation matrix R. The current volume change rate of the intervertebral disc nucleus pulposus layer is obtained. This rate is defined as the volume of the finite element at the current moment divided by its volume under initial no-load conditions; this ratio is calculated using the determinant of the deformation gradient tensor. The hydrostatic pressure value of the intervertebral disc nucleus pulposus layer is updated based on the current volume change rate. The update method is: subtract 1 from the volume change rate, multiply by the bulk modulus of the nucleus pulposus, and then multiply by -1 to obtain the hydrostatic pressure value. The bulk modulus of the nucleus pulposus was set to 500 times the initial shear modulus during the material parameter assignment stage. The updated hydrostatic pressure value is added to the diagonal terms of the stress tensor of the intervertebral disc nucleus pulposus layer, specifically to the elements in the first row and first column, the second row and second column, and the third row and third column of the stress tensor.

[0030] When the displacement increment of any node in the multilayer composite structure between two adjacent increment steps is less than the convergence tolerance, the iteration stops and the stress-strain response field corresponding to the current increment step is output. The first displacement vector of all nodes at the end of the current increment step is obtained. This is achieved by extracting the displacement components of all nodes from the output vector of the nonlinear finite element solver. Each node's first displacement vector contains its displacement values ​​in the x, y, and z directions, with a dimension of 3 times the total number of nodes. The second displacement vector of all nodes at the end of the previous increment step is obtained. This is stored in the same way as the first displacement vector, having been saved to memory after convergence in the previous increment step. The node displacement increment vector is obtained by subtracting the first and second displacement vectors, with the subtraction operation performed one component at a time. The Euclidean norm of the node displacement increment vector is calculated by squaring the value of each component in the vector, summing all the squares, and then taking the square root of the sum. The total nodal displacement is defined as the square root of the sum of the squares of all nodal displacement vector components in the current increment step. The normalized displacement error is obtained by dividing the Euclidean norm by the total nodal displacement. The preset convergence threshold is set to 1×10⁻⁶. -5This value originates from the standard convergence criterion of finite element analysis, indicating that the iteration is considered fully converged when the ratio of the displacement increment to the total displacement is less than one ten-thousandth. When the normalized displacement error value is less than the preset convergence threshold, the current increment step is determined to meet the convergence condition, and the stress-strain response field corresponding to the current increment step is written to the output file. The stress-strain response field includes six stress components and six strain components at all integration points, stored in binary VTK file format. When the normalized displacement error value is greater than or equal to the preset convergence threshold, the Jacobian matrix of the multilayer material composite structure is updated, and the next iteration step is entered again. The Jacobian matrix update is achieved through the global tangent stiffness matrix assembly algorithm built into the finite element solver. The updated Jacobian matrix is ​​substituted into the Newton-Raphson iterative scheme to solve for the displacement correction vector of the next iteration step.

[0031] See Figure 5 In the figure, the vertical axis uses a logarithmic scale to represent the magnitude of the normalized displacement error, and the horizontal axis represents the number of iterations in the nonlinear finite element iterative calculation, with a total of 500 iterations. The curve shows the trend of normalized displacement error as a function of the number of iterations, and the dashed line represents the preset convergence threshold of 1×10⁻⁶. -5 As can be seen from the figure, the initial value of the normalized displacement error is relatively large, exceeding 1×10⁻⁶. -2 The error generally decreases with increasing iteration count, but decreases more rapidly in the first 100 iterations, with the error range decreasing from approximately 1×10⁻⁶. -2 Decreased to approximately 1×10 -3 Subsequently, within the interval of 100 to 500 iterations, the normalized displacement error fluctuated within 1×10⁻⁶. -3 The amplitude fluctuated slightly, exhibiting multiple small oscillations, failing to significantly decrease below the convergence threshold. The convergence threshold corresponding to the dashed line is 1×10⁻⁶. -5 The error level is significantly lower than that shown by the curve, indicating that the simulation has not yet met the preset convergence condition within 500 iterations. This curve reflects the convergence behavior of the multilayer composite structure in Example 3 under the applied preset physiological load boundary conditions during nonlinear finite element iterative calculations. It shows that the iteration process can rapidly reduce errors in the initial stage, but the convergence speed slows down and oscillates in subsequent iterations, suggesting the possible existence of complex nonlinear material responses or numerical stability issues. The overall trend and numerical range verify the iterative characteristics of the nonlinear finite element solver and the rationality of the preset convergence criterion.

[0032] Example 4: In practice, the process of gradually increasing the compressive load value on the load application surface while applying a displacement-dependent tension boundary condition at the attachment point of the ligament attachment layer of the multi-layer composite structure is achieved in the following way.

[0033] The current load step value of the load application surface is obtained. This value is derived from the external load incrementing loop variable of the nonlinear finite element solution program. This loop variable is initialized to zero at the start of the solution and automatically increments by 1 after each complete external load step. The current load step value is multiplied by the single-step load increment to obtain the current total compressive load. The single-step load increment is determined by: setting a target total compressive load value, dividing the target total compressive load value by the maximum load step value. The target total compressive load value is set to 750 Newtons. This value is calculated based on the load on the upper torso mass borne by an adult male of standard weight in an upright position in the L4-L5 segment, and rounded down. The calculation method is: upper torso mass approximately 40 kg multiplied by gravitational acceleration approximately 9.8 m / s², then multiplied by the spinal load distribution factor of 0.5, resulting in approximately 196 Newtons. Considering the amplification effect of daily activities such as bending over and lifting objects, the load is increased to approximately 3.8 times, resulting in 750 Newtons. The maximum load step value is set to 100. This is based on the principle that dividing the total compressive load target value into 100 equidistant incremental steps ensures that the structural deformation amplitude between adjacent incremental steps is sufficiently small during quasi-static loading, avoiding convergence difficulties, while also considering computational efficiency. The current total compressive load is evenly distributed to each node on the load application surface. The even distribution method is as follows: count the total number of nodes in the load application node set, divide the current total compressive load by the total number of nodes to obtain the compressive load value borne by each node, and apply a concentrated force equal to the compressive load value in the negative z-axis direction of the global coordinate system of each node.

[0034] Simultaneously, the initial spatial coordinates of the initiation and termination points of the ligament attachment layer are obtained. The initial spatial coordinates are extracted from the finite element node coordinate array and stored as a static array after the 3D geometry of the ligament attachment layer is constructed and no external load is applied. In each increment step of the nonlinear finite element iterative calculation, the current distance between the initiation and termination points is calculated. The current position coordinates of the initiation point are obtained by adding the three components of the initial spatial coordinates of the initiation point to the three components of the displacement vector of the initiation point obtained in the current increment step. The current position coordinates of the termination point are obtained in the same way from the corresponding initial spatial coordinates and displacement vector. The Euclidean distance between the current position coordinates of the initiation and termination points is calculated and used as the current distance between them. The ligament elongation is obtained by subtracting the initial distance between the initiation and termination points from the current distance between them. The initial distance is obtained by calculating the Euclidean distance from the initial spatial coordinates of the initiation and termination points. The current tension value is obtained by substituting the ligament elongation into the parameters of the nonlinear stress-strain curve of the ligament attachment layer. Specifically, the substitution process is as follows: the ligament elongation is divided by the relaxation length of the ligament attachment layer to obtain the ligament engineered strain, with the relaxation length taken as the initial distance. The ligament engineered strain is then input into the nonlinear stress-strain curve function of the ligament attachment layer. This function uses a piecewise approach; when the ligament engineered strain is less than the first critical tensile ratio... When the output is reduced by 1, the low stiffness stress in the toe region is reached, and the ligament engineering strain is between and The output is linear high-stiffness stress between these values, and when the ligament engineering strain is greater than... The system outputs the post-yield stress, which is the ligament fiber stress. Multiplying this stress by the cross-sectional area of ​​the ligament attachment layer yields the current tension value. The cross-sectional area of ​​the ligament attachment layer is set to 50 square millimeters based on statistical averages from anatomical specimens. The current tension value is applied to both the initial and final attachment points along a direction from the initial attachment point to the final attachment point. The application method is as follows: a unit direction vector is calculated from the initial attachment point to the final attachment point; a concentrated force vector is applied to the initial attachment point, with a magnitude equal to the current tension value and a direction opposite to the unit direction vector; a concentrated force vector is also applied to the final attachment point, with a magnitude equal to the current tension value and a direction the same as the unit direction vector.

[0035] The local coordinate system of the annulus fibrosus layer of the intervertebral disc is obtained, and a fiber-reinforced hyperelastic constitutive update algorithm is established based on the local coordinate system. The process of updating the fiber orientation stress of the annulus fibrosus layer and then updating the hydrostatic pressure of the nucleus pulposus layer of the intervertebral disc in each incremental step is achieved in the following way.

[0036] Obtain the radial and circumferential direction vectors for each integration point in the annulus fibrosus of the intervertebral disc. The radial direction vector is obtained by normalizing the global coordinate vector of the integration point after subtracting the geometric center coordinate vector of the nucleus pulposus layer. The geometric center coordinates of the nucleus pulposus layer are obtained by calculating the arithmetic mean of the coordinates of all nodes within the nucleus pulposus layer. The circumferential direction vector is obtained by performing a cross product operation between the radial direction vector and the global vertical direction vector (0,0,1). If the magnitude of the cross product vector is greater than 10 to the power of -8, the cross product vector is normalized to obtain the circumferential direction vector. If the magnitude of the cross product vector is less than or equal to 10 to the power of -8, i.e., the radial direction vector is parallel to the global vertical direction vector, the circumferential direction vector is obtained by performing a cross product operation between the radial direction vector and the global anterior-posterior direction vector (0,1,0) and normalizing the result. A local coordinate system is constructed based on the radial and circumferential direction vectors. The three orthogonal basis vectors of the local coordinate system are set as follows: the first basis vector is the circumferential direction vector, the second basis vector is the radial direction vector, and the third basis vector is obtained by the cross product of the first and second basis vectors, ensuring that the local coordinate system satisfies the right-hand rule and that the three basis vectors are pairwise orthogonal.

[0037] The current tensile ratio in the fiber direction of the annulus fibrosus of the intervertebral disc is calculated in a local coordinate system. The fiber direction is given by the principal fiber orientation direction at the corresponding integration point in the collagen fiber orientation distribution map, which was obtained and mapped to the integration point during the material parameter assignment stage using diffusion tensor magnetic resonance imaging. The principal fiber orientation direction vector is projected onto the plane formed by the first and third basis vectors of the local coordinate system, and the magnitude of the projected vector is the current tensile ratio in the fiber direction. The corresponding fiber stress increment is found based on the current tensile ratio. The lookup method is as follows: a discrete data mapping table of fiber tensile ratio and fiber stress increment is pre-generated. This mapping table is constructed by segmented cubic Hermite interpolation of uniaxial tensile test data of isolated intervertebral disc annulus fibrosus. The tensile ratio range in the test data covers 0.8 to 2.0, and the stress increment range covers 0 MPa to 50 MPa. The interval in which the current tensile ratio is located in the mapping table is located, and the fiber stress increment is calculated by linear interpolation of the two endpoints of this interval. The fiber stress increment is superimposed on the fiber stress of the previous increment step to obtain the current fiber stress. The superposition method is to directly add the values.

[0038] Transform the current fiber stress to the global coordinate system and update the global stress tensor of the intervertebral disc annulus fibrosus. The transformation is achieved using the following formula: ; In the formula, The fiber stress tensor in the global coordinate system is a 3x3 symmetric matrix, and the unit of each matrix element is megapascal. Let represent the fiber stress tensor in the local coordinate system. It is a 3x3 symmetric matrix, where non-zero normal stress components are only present at the diagonal positions corresponding to the fiber direction. These components are equal to the current fiber stress, and all other matrix elements are zero. Q represents the coordinate transformation matrix, which is also a 3x3 matrix. The first column of the matrix consists of the three components of the first basis vector of the local coordinate system in the global coordinate system arranged in rows. The second column consists of the three components of the second basis vector of the local coordinate system in the global coordinate system arranged in rows. The third column consists of the three components of the third basis vector of the local coordinate system in the global coordinate system arranged in rows. The value of each element of the matrix ranges from -1 to 1. This represents the transpose of the coordinate transformation matrix Q. The global stress tensor is calculated. Then, the six independent stress components of the global stress tensor are updated in the global stress storage array at that integration point.

[0039] The current volume change rate of the intervertebral disc nucleus pulposus is obtained. This rate is defined as the volume of a finite element in the current deformed configuration divided by the volume of that element in the initial unloaded configuration. This ratio is obtained through the determinant of the deformation gradient tensor, which is calculated from the node displacements of the current increment step. The hydrostatic pressure value of the intervertebral disc nucleus pulposus is updated based on the current volume change rate. The update method is as follows: subtract 1 from the current volume change rate, multiply by the bulk modulus of the intervertebral disc nucleus pulposus, and then multiply by -1 to obtain the hydrostatic pressure value. The bulk modulus of the intervertebral disc nucleus pulposus is set to 500 times the initial shear modulus of the nucleus pulposus during the material parameter assignment stage. The initial shear modulus is 0.2 MPa, therefore the bulk modulus is 100 MPa. The updated hydrostatic pressure value is added to the diagonal terms of the stress tensor of the intervertebral disc nucleus pulposus, that is, the hydrostatic pressure value is added to the first row and first column, the second row and second column, and the third row and third column of the stress tensor, respectively, completing the stress state update of the nucleus pulposus in the current increment step.

[0040] See Figure 6 The figure shows the stress-engineered strain curve of the ligament fiber. The horizontal axis represents the ligament's engineered strain, ranging from 0 to 0.2, and the vertical axis represents the ligament fiber stress, in megapascals (MPa). The curves exhibit a typical nonlinear stress-strain relationship, reflecting the mechanical response characteristics of the ligament fiber during loading.

[0041] The curve can be divided into three stages from left to right: The first stage is from engineering strain 0 to about 0.03 (corresponding to the vertical dashed line at lambda1=1.03 in the figure). In this stage, the stress of the ligament fibers increases slowly, showing the low stiffness characteristics of the toe area, reflecting the relaxation and initial buffering capacity of the ligament fiber microstructure; The second stage is from engineering strain about 0.03 to 0.08 (corresponding to the vertical dashed line at lambda2=1.08 in the figure). In this stage, the stress of the ligament fibers increases significantly with strain, showing a linear high stiffness increase, representing that the fibers have entered the elastic tensile state and the ligament load-bearing capacity has increased rapidly; The third stage is the part after the engineering strain is greater than 0.08. The increase in stress on the curve slows down and is accompanied by fluctuations, reflecting that the ligament has entered the post-yield stage, which is manifested as micro-damage or nonlinear plastic deformation of the fibers.

[0042] The curve's twists and turns reflect, to some extent, the nonlinear iterative characteristics and material response complexity during the numerical calculation process. The peak stress on the curve reaches approximately 4.8 MPa, demonstrating the maximum load-bearing capacity of the ligament fibers under high strain conditions.

[0043] The data in this figure closely follows the description of the nonlinear stress-strain curve parameters of the ligament attachment layer in Example 4, clarifying the three stages of ligament engineering strain and the corresponding stress response characteristics. It verifies the simulation results of the mechanical behavior of the ligament attachment layer during load application and demonstrates the accurate simulation capability of the nonlinear finite element iterative calculation method of this invention for the mechanical response of ligament fibers.

[0044] Example 5: In practice, when the displacement increment of any node in the multilayer material composite structure between two adjacent increment steps is less than the convergence tolerance, the process of stopping the iteration and outputting the stress-strain response field corresponding to the current increment step is achieved in the following way.

[0045] Obtain the first displacement vector of all nodes at the end of the current increment step. The first displacement vector is obtained by reading the displacement solution corresponding to the latest converged increment step from the global displacement array of the nonlinear finite element solver. The global displacement array is a one-dimensional double-precision floating-point array, and the array length is equal to the total number of nodes in the finite element mesh multiplied by the spatial dimension 3. The 3×k-2th element in the array represents the displacement component of the k-th node in the x-direction, the 3×k-1th element represents the displacement component of the k-th node in the y-direction, and the 3×kth element represents the displacement component of the k-th node in the z-direction. k is the node index number from 1 to the total number of nodes. Obtain the second displacement vector of all nodes at the end of the previous increment step. The second displacement vector has been copied from the global displacement array to a separate memory buffer after the previous increment step converged. The data format in the buffer is the same as that of the first displacement vector. The first displacement vector is subtracted from the second displacement vector to obtain the node displacement increment vector. The subtraction operation is performed one by one according to the elements with the same index position. That is, the j-th element of the node displacement increment vector is equal to the j-th element of the first displacement vector minus the j-th element of the second displacement vector. The value of j ranges from 1 to 3 times the total number of nodes.

[0046] The Euclidean norm of the nodal displacement increment vector is calculated as follows: Square each element of the nodal displacement increment vector, sum all the squared values, and then take the square root of the sum to obtain the Euclidean norm. The total nodal displacement is then calculated, defined as the Euclidean norm of all nodal displacements in the current increment step, which is obtained by squaring, summing, and taking the square root of all elements of the first displacement vector. Dividing the Euclidean norm by the total nodal displacement yields the normalized displacement error value. The preset convergence threshold is set to 1 × 10⁻⁶. -5 The value is set based on the following: In the Newton-Raphson iteration convergence criterion for finite element calculation, when the relative error of the displacement norm is less than one part per hundred thousand, it indicates that the displacement correction of the current iteration step is much smaller than the current overall deformation of the structure, and continuing the iteration will not have an observable effect on the stress-strain distribution. Moreover, this threshold is widely used as a convergence criterion in structural nonlinear analysis.

[0047] When the normalized displacement error is less than the preset convergence threshold, the current increment step is deemed to meet the convergence condition, and the stress-strain response field corresponding to the current increment step is written to the output file. The stress-strain response field includes the stress and strain components of all integration points in the multilayer composite structure. The stress components contain six independent components: normal stress in the x-direction, normal stress in the y-direction, normal stress in the z-direction, shear stress in the xy plane, shear stress in the yz plane, and shear stress in the zx plane. The strain components also contain six independent components. The output file is stored in VTK format, and the coordinates, stress components, and strain components of each integration point are written to a binary data block in sequence. When the normalized displacement error is greater than or equal to the preset convergence threshold, the Jacobian matrix of the multilayer composite structure is updated, and the next iteration step is initiated. The Jacobian matrix of a multilayer composite structure is updated as follows: The finite element solver calls the global tangent stiffness matrix assembly function, iterates through all elements, calculates the tangent stiffness matrix at each integration point based on the current material constitutive relation and deformation state, and then assembles the element tangent stiffness matrices into the global Jacobian matrix according to the node connection relationship. The Jacobian matrix is ​​a sparse symmetric positive definite matrix, and its storage format adopts compressed row storage or compressed column storage format. The updated Jacobian matrix is ​​substituted into the Newton-Raphson scheme to solve the linear equation system to obtain the displacement correction vector, and the nodal displacements are updated in the next Newton iteration.

[0048] The process of identifying stress concentration regions at the interfaces of each material layer in a multilayer composite structure based on the stress-strain response field and marking these stress concentration regions as output points for mechanical behavior simulation is achieved in the following way.

[0049] Obtain the first stress tensor field at the interface between the vertebral bone tissue layer and the intervertebral disc annulus fibrosus layer in the stress-strain response field. The node set at the interface between the vertebral bone tissue layer and the intervertebral disc annulus fibrosus layer is determined as follows: traverse all elements in the multilayer composite structure, search for elements that simultaneously contain nodes belonging to both the vertebral bone tissue layer material and the intervertebral disc annulus fibrosus layer material. Extract nodes shared by both materials from these elements, and remove duplicate nodes to form the interface node set. Read the stress tensor corresponding to each node in the interface node set from the output stress-strain response field file. The stress tensor contains six independent components: , , , , , Each component is read from the file according to its node index number.

[0050] Calculate the Mises equivalent stress value at each node in the first stress tensor field. The formula for calculating the Mises equivalent stress value is as follows: ; In the formula, The value represents the Mises equivalent stress at the node, in megapascals (MPa). This represents the normal stress component of the node in the x-direction, in megapascals (MPa). This represents the normal stress component of the node in the y-direction, in megapascals (MPa). This represents the normal stress component of the node in the z-direction, in megapascals (MPa). This represents the shear stress components of the node in the xy plane, in megapascals (MPa). This represents the shear stress component of the node in the yz plane, in megapascals (MPa). This represents the shear stress components of the node in the zx plane, in megapascals (MPa). All stress components are directly extracted from the first stress tensor field.

[0051] Nodes with Mises equivalent stress values ​​greater than a preset stress threshold are marked as candidate stress concentration nodes. The preset stress threshold is set as the 90th percentile of the Mises equivalent stress values ​​of all nodes at the interface between the vertebral bone tissue layer and the intervertebral disc annulus fibrosus layer. This value is determined as follows: the Mises equivalent stress values ​​of all nodes at the interface are sorted in ascending order, and the value at the 90th percentile is taken as the preset stress threshold. The 90th percentile is chosen because it only includes the 10% of nodes with the highest stress in the candidate list, which can filter out non-concentrated background stress nodes while retaining the true stress concentration area.

[0052] Connectivity clustering analysis is performed on candidate stress concentration nodes to obtain multiple connected clusters of stress concentration regions. The connectivity clustering analysis is implemented using a region growing algorithm in three-dimensional space. Specifically, an access marker array of the same length as the total number of nodes is created, and all elements are initialized to an unvisited state. The candidate stress concentration node list is traversed, and for each unvisited candidate stress concentration node, region growing is performed: the node is marked as visited, added to the newly created connected cluster node list, and pushed into the growth queue. The head node is removed from the growth queue, and all its first-order adjacent nodes are searched in the topological connectivity relationships of the interface mesh. A first-order adjacent node is defined as a node directly connected to the current node through a triangle facet edge. For each first-order adjacent node, if it is both a candidate stress concentration node and unvisited, it is marked as visited, added to the current connected cluster node list, and pushed into the growth queue. This process is repeated until the growth queue is empty, thus completing the extraction of a connected cluster of stress concentration regions. The process continues to traverse the next unvisited candidate stress concentration node and repeat the region growing process until all candidate stress concentration nodes have been visited.

[0053] Obtain the average stress value and number of nodes in all nodes of each stress concentration region connected cluster. Calculate the stress concentration weight coefficient for each stress concentration region connected cluster based on the average stress value and number of nodes. Obtain the first Mises equivalent stress value list for all nodes in the first stress concentration region connected cluster. Sum all values ​​in the first Mises equivalent stress value list and divide the sum by the number of first nodes in the first stress concentration region connected cluster to obtain the first average stress value of the first stress concentration region connected cluster. Obtain the second Mises equivalent stress value list for all nodes in the second stress concentration region connected cluster. Sum all values ​​in the second Mises equivalent stress value list and divide the sum by the number of second nodes in the second stress concentration region connected cluster to obtain the second average stress value of the second stress concentration region connected cluster. Multiply the first average stress value by the number of first nodes to obtain the first weighted stress sum, and multiply the second average stress value by the number of second nodes to obtain the second weighted stress sum. For cases with more than two stress concentration region connected clusters, calculate the weighted stress sum for each connected cluster using the same method. Compare the sum of weighted stresses across all connected clusters and select the cluster with the largest sum. Calculate the arithmetic mean of the coordinates of all nodes in the connected cluster with the largest sum of weighted stresses. This arithmetic mean is obtained by summing the x-coordinates of all nodes in the connected cluster and dividing by the number of nodes, summing the y-coordinates and dividing by the number of nodes, and summing the z-coordinates and dividing by the number of nodes. The spatial point formed by these three coordinate averages is marked as the simulation output point for mechanical behavior. Simultaneously, output a list of node numbers occupied by the connected cluster with the largest sum of weighted stresses. The node numbers in the list are integer values ​​from the global node numbers of each node in the finite element mesh within that connected cluster. The node number list is written to the simulation output file in text format.

[0054] See Figure 7 In the figure, the horizontal axis represents the node number at the interface between the vertebral bone tissue layer and the intervertebral disc annulus fibrosus layer, ranging from 0 to 500, and the vertical axis represents the Mises equivalent stress value at the corresponding node, in megapascals (MPa). The solid curve depicts the distribution of Mises equivalent stress at each interface node calculated based on the multilayer material structure simulation method, while the dashed line represents the preset stress threshold, which is approximately 28 MPa.

[0055] Observing the trend of the curve, the Mises equivalent stress values ​​show significant fluctuations within the node number range, with two relatively significant stress peak regions: one is located between node numbers approximately 150 and 180, corresponding to a stress peak value reaching and exceeding the preset stress threshold of nearly 40 MPa; the other is located between node numbers approximately 360 and 390, where the stress peak value also exceeds the threshold, approaching 32 MPa. Apart from these, the stress values ​​at other nodes are mainly concentrated in the 10 to 25 MPa range, below the preset threshold.

[0056] This stress distribution reflects a significant stress concentration region at the interface between the vertebral bone tissue and the annulus fibrosus of the intervertebral disc in a multilayered composite structure under physiological load conditions. By filtering using a preset stress threshold, nodes located within the stress peak range in the figure are marked as candidate stress concentration nodes, conforming to the stress concentration region identification criteria based on Mises equivalent stress in Example 5, providing a basis for subsequent connected domain clustering analysis. This analysis helps to accurately locate potential mechanical weaknesses or high-risk areas in the interface, supporting the determination and marking of mechanical behavior simulation output points.

[0057] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any changes 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.

Claims

1. A simulation method for the mechanical behavior of the spine based on a multilayer material structure, characterized in that, include: Acquire medical imaging data of the target spinal segment, and construct an initial three-dimensional geometric model based on the medical imaging data; The initial three-dimensional geometric model is segmented into multiple material regions to obtain the vertebral bone tissue layer, the intervertebral disc annulus fibrosus layer, the intervertebral disc nucleus pulposus layer, the ligament attachment layer, and the nerve root sheath layer; Different viscoelastic constitutive parameters are assigned to the vertebral bone tissue layer, the intervertebral disc annulus fibrosus layer, the intervertebral disc nucleus pulposus layer, the ligament attachment layer, and the nerve root sheath layer, respectively, to form a multilayer material composite structure, specifically including: Obtain the apparent density distribution cloud map of the vertebral bone tissue layer, and assign anisotropic elastic modulus to each local region of the vertebral bone tissue layer according to the apparent density distribution cloud map; Obtain the collagen fiber orientation distribution map of the intervertebral disc annulus fibrosus layer, and assign position-dependent tensile hardening parameters to different radial positions of the intervertebral disc annulus fibrosus layer based on the collagen fiber orientation distribution map; Obtain the fluid volume fraction of the nucleus pulposus layer of the intervertebral disc, and assign incompressible fluid parameters to the nucleus pulposus layer of the intervertebral disc based on the fluid volume fraction; The initial fiber pre-stretch ratio of the ligament attachment layer is obtained, and nonlinear stress-strain curve parameters are assigned to the ligament attachment layer based on the initial fiber pre-stretch ratio. Obtain the axonal direction vector of the nerve root sheath, and assign transverse isotropic viscoelastic parameters to the nerve root sheath based on the axonal direction vector; A preset physiological load boundary condition is applied to the multilayer material composite structure, and a nonlinear finite element iterative calculation is performed under the physiological load boundary condition to obtain the stress-strain response field of the multilayer material composite structure under load. Based on the stress-strain response field, the stress concentration regions at the interfaces of each material layer in the multilayer composite structure are identified in reverse. These stress concentration regions at the interfaces are then marked as output points for mechanical behavior simulation, specifically including: Obtain the first stress tensor field at the interface between the vertebral bone tissue layer and the intervertebral disc annulus fibrosus layer in the stress-strain response field, and calculate the Mises equivalent stress value at each node in the first stress tensor field; Nodes whose Mises equivalent stress values ​​are greater than a preset stress threshold are marked as candidate stress concentration nodes. Connectivity clustering analysis is performed on the candidate stress concentration nodes to obtain multiple connected clusters of stress concentration regions. Obtain the average stress value and number of nodes in all nodes of each stress concentration region connected cluster. Calculate the stress concentration weight coefficient of each stress concentration region connected cluster based on the average stress value and the number of nodes. Use the geometric center coordinates of the connected cluster with the largest stress concentration weight coefficient as the output point of the mechanical behavior simulation, and simultaneously output the list of node numbers occupied by that connected cluster.

2. The method for simulating the mechanical behavior of the spine based on a multilayer material structure according to claim 1, characterized in that, The initial three-dimensional geometric model is segmented into multiple material regions to obtain the vertebral bone tissue layer, intervertebral disc annulus fibrosus layer, intervertebral disc nucleus pulposus layer, ligament attachment layer, and nerve root sheath layer, specifically including: Obtain the grayscale gradient distribution map of the initial three-dimensional geometric model, and determine the initial interface between the spinal bone tissue region and the intervertebral disc region based on the grayscale gradient distribution map; The mean gray value of the central region and the gray gradient change rate of the edge region are extracted within the intervertebral disc region. Based on the mean gray value of the central region and the gray gradient change rate of the edge region, the intervertebral disc region is divided into the annulus fibrosus layer and the nucleus pulposus layer. Obtain the ligament origin attachment coordinates and ligament termination attachment coordinates in the initial three-dimensional geometric model, and peel the ligament attachment layer from the surface of the vertebral bone tissue layer according to the ligament origin attachment coordinates and ligament termination attachment coordinates. Based on the internal contour line of the intervertebral foramen in the initial three-dimensional geometric model, the nerve root sheath layer is obtained by extending outward by a predetermined radial distance along the internal contour line of the intervertebral foramen.

3. The spinal biomechanical behavior simulation method based on a multilayer material structure according to claim 2, characterized in that, A preset physiological load boundary condition is applied to the multilayer material composite structure, and a nonlinear finite element iterative calculation is performed under the physiological load boundary condition to obtain the stress-strain response field of the multilayer material composite structure under load. Specifically, this includes: Obtain the coordinates of the outer surface nodes of the vertebral bone tissue layer in the multilayer material composite structure, and set the set of nodes located at the upper endplate of the upper vertebral body in the outer surface node coordinates as the load application surface; The compressive load value is gradually increased on the load application surface, while a displacement-dependent tension boundary condition is applied at the attachment point of the ligament attachment layer of the multilayer material composite structure. The local coordinate system of the annulus fibrosus layer of the intervertebral disc is obtained, and a fiber-reinforced hyperelastic constitutive update algorithm is established based on the local coordinate system. In each incremental step, the fiber orientation stress of the annulus fibrosus layer of the intervertebral disc is updated first, and then the hydrostatic pressure of the nucleus pulposus layer of the intervertebral disc is updated. When the displacement increment of any node in the multilayer material composite structure between two adjacent increment steps is less than the convergence tolerance, the iteration stops and the stress-strain response field corresponding to the current increment step is output.

4. The spinal biomechanical behavior simulation method based on a multilayer material structure according to claim 3, characterized in that, Obtaining the coordinates of the outer surface nodes of the vertebral bone tissue layer in the multilayer material composite structure, and setting the set of nodes located at the upper endplate of the superior vertebral body in the outer surface node coordinates as the load application surface specifically includes: Obtain all triangular facets on the outer surface of the vertebral bone tissue layer, and calculate the dot product of the normal vector of each triangular facet and the global vertical upward direction vector; All triangular facets whose dot product value is greater than the direction cosine threshold are selected as candidate final plate facets; Obtain the coordinate values ​​of all nodes in the candidate endplate patch, calculate the axial average value of the coordinate values ​​of all nodes, and determine the nodes whose axial coordinates are greater than the axial average value as the upper endplate nodes. The index number of the upper end plate node is stored in the load application node set, and the triangular facet region occupied by the load application node set is marked as the load application surface.

5. The spinal biomechanical behavior simulation method based on a multilayer material structure according to claim 3, characterized in that, Gradually increasing the compressive load value on the load application surface, while simultaneously applying displacement-dependent tension boundary conditions at the attachment points of the ligament attachment layer in the multilayer composite structure, specifically includes: Obtain the current load step value of the load application surface, multiply the current load step value by the single-step load increment to obtain the current total compressive load, and evenly distribute the current total compressive load to each node of the load application surface; Obtain the initial spatial coordinates of the starting attachment point and the ending attachment point of the ligament attachment layer, and calculate the current distance between the starting attachment point and the ending attachment point in each incremental step of the nonlinear finite element iterative calculation. The ligament elongation is obtained by subtracting the initial distance between the starting attachment point and the ending attachment point from the current distance between the starting attachment point and the ending attachment point. The ligament elongation is then substituted into the nonlinear stress-strain curve parameters of the ligament attachment layer to obtain the current tension value. The current tension value is applied to the starting attachment point and the ending attachment point respectively in the direction from the starting attachment point to the ending attachment point.

6. The spinal biomechanical behavior simulation method based on a multilayer material structure according to claim 3, characterized in that, Obtain the local coordinate system of the annulus fibrosus layer of the intervertebral disc, and establish a fiber-reinforced hyperelastic constitutive update algorithm based on the local coordinate system. In each incremental step, first update the fiber orientation stress of the annulus fibrosus layer of the intervertebral disc, and then update the hydrostatic pressure of the nucleus pulposus layer of the intervertebral disc. Specifically, this includes: Obtain the radial direction vector and circumferential direction vector of each integration point in the annulus fibrosus layer of the intervertebral disc, and construct the local coordinate system based on the radial direction vector and the circumferential direction vector; Calculate the current tensile ratio in the fiber direction of the intervertebral disc annulus fibrosus in the local coordinate system, find the corresponding fiber stress increment based on the current tensile ratio, and add the fiber stress increment to the fiber stress of the previous increment step to obtain the current fiber stress. Transform the current fiber stress to the global coordinate system and update the global stress tensor of the intervertebral disc annulus fibrosus layer; Obtain the current volume change rate of the nucleus pulposus layer of the intervertebral disc, update the hydrostatic pressure value of the nucleus pulposus layer of the intervertebral disc according to the current volume change rate, and add the updated hydrostatic pressure value to the diagonal term of the stress tensor of the nucleus pulposus layer of the intervertebral disc.

7. The spinal biomechanical behavior simulation method based on a multilayer material structure according to claim 3, characterized in that, When the displacement increment of any node in the multilayer composite structure between two adjacent increment steps is less than the convergence tolerance, the iteration stops and the stress-strain response field corresponding to the current increment step is output, specifically including: Obtain the first displacement vector of all nodes at the end of the current increment step, obtain the second displacement vector of all nodes at the end of the previous increment step, and subtract the first displacement vector from the second displacement vector to obtain the node displacement increment vector. Calculate the Euclidean norm of the node displacement increment vector, and divide the Euclidean norm by the total displacement of all nodes to obtain the normalized displacement error value. When the normalized displacement error value is less than the preset convergence threshold, it is determined that the current increment step meets the convergence condition, and the stress-strain response field corresponding to the current increment step is written to the output file. When the normalized displacement error value is greater than or equal to the preset convergence threshold, the Jacobian matrix of the multilayer material composite structure is updated and the next iteration step is entered.

8. The spinal biomechanical behavior simulation method based on a multilayer material structure according to claim 7, characterized in that, Obtain the average stress value and number of nodes in all nodes of each stress concentration region connected cluster. Calculate the stress concentration weight coefficient of each stress concentration region connected cluster based on the average stress value and the number of nodes. Use the geometric center coordinates of the connected cluster with the largest stress concentration weight coefficient as the mechanical behavior simulation output point. Specifically, this includes: Obtain a list of the first Mises equivalent stress values ​​for all nodes in the first stress concentration region connected cluster. Sum the first Mises equivalent stress value list and divide it by the number of the first nodes in the first stress concentration region connected cluster to obtain the first average stress value of the first stress concentration region connected cluster. Obtain the list of second Mises equivalent stress values ​​for all nodes in the connected cluster of the second stress concentration region. Sum the list of second Mises equivalent stress values ​​and divide it by the number of second nodes in the connected cluster of the second stress concentration region to obtain the second average stress value of the connected cluster of the second stress concentration region. Multiply the first average stress value by the first number of nodes to obtain the first weighted stress sum, and multiply the second average stress value by the second number of nodes to obtain the second weighted stress sum; Compare the magnitudes of the first weighted stress sum and the second weighted stress sum, select the connected cluster with the larger weighted stress sum, calculate the arithmetic mean of the coordinates of all nodes in the connected cluster, and mark the spatial point corresponding to the arithmetic mean as the mechanical behavior simulation output point.

Images