An isogeometric analysis simulation method based on volumetric subdivision of a hyperelastic material model

Abstract

The application discloses an isogeometric analysis simulation method of a hyperelastic material model based on body subdivision. The isogeometric analysis adopts the division of a spline model parameter domain and the mapping of the parameter domain to a physical model, avoids the time consumption of mesh division in a finite element, and realizes an integrated process from modeling to simulation in combination with a CC subdivision method. The application comprises the following steps: converting a CC subdivision body into a spline body model, and establishing a mapping relationship of control points of the subdivision body and the spline body; establishing a hyperelastic material constitutive model according to geometric boundary conditions and material parameters; establishing a discrete balance equation of the model according to an isogeometric method; in combination with the balance equation, establishing a model motion equation and a nonlinear solving system; using a Newton iteration method to solve the nonlinear system, obtaining displacement solutions at each time step, and mapping the displacement solutions back to the subdivision body, so as to realize deformation simulation of the subdivision body.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of super-elastic material simulation, in particular to a super-elastic material model based on volume subdivision isometric analysis simulation method BACKGROUND

[0002] The present application relates to the field of super-elastic material model deformation simulation process, in particular to using isometric analysis method to solve Catmull-Clark (CC) subdivision model under stress deformation process.

[0003] Super-elastic model simulation has always been a research hotspot in computer graphics, and in industrial manufacturing, finite element method is usually used to perform numerical simulation calculation on the model, and the precision of simulation is determined by the number and quality of grid division. However, the divided grid is an approximation to the original accurate geometry, so the model has a certain degree of data distortion after grid division, and for complex models, to improve the precision, the grid is often refined to ensure the simulation quality, which also leads to a great increase in simulation time. Isometric analysis is based on the idea of finite element discretization, and uses non-uniform rational B-spline (NURBS) used in CAD model as the shape function of the analysis process, which directly simulates and analyzes the accurate geometry without going through the traditional grid division process, greatly improving the simulation precision. The NURBS function used in isometric analysis has high-order continuity, which can meet the needs of high-order differentiation in engineering analysis, and the accurate geometric description makes the coarsest discrete grid have better simulation analysis results than finite element method.

[0004] Volume spline modeling of complex geometric models is a key problem in isometric analysis, and the corresponding spline volume can be defined by a set of control grids and a set of refinement rules. These rules act on the given control grid to generate new refined grid. Catmull-Clark (CC) volume subdivision can directly generate an unstructured three-variable spline volume as an approximation of the subdivision limit volume from a given hexahedral grid through a simple two-step process. Compared with traditional approximation subdivision, this method is simple, easy to implement, and requires less resource consumption, and can quickly generate refined grids of any precision.

[0005] The research on isometric simulation analysis mainly focuses on statics, fluid mechanics, plate and shell problems, electromagnetism, etc. The deformation simulation of super-elastic material model belongs to nonlinear large deformation simulation analysis, and isometric analysis has less research in this field. SUMMARY

[0006] In order to provide a fast and effective super-elastic material deformation simulation method for CC subdivision, the present application uses isogeometric analysis instead of finite element analysis to simulate and analyze the spline model converted from the CC subdivision, and can calculate the displacement size of each time step for the complex three-dimensional geometric model, so as to realize the deformation simulation of the super-elastic material model.

[0007] The technical scheme adopted by the present application is as follows: converting the CC subdivision into a spline body model, establishing the mapping relationship between the control points of the subdivision and the spline body; establishing the constitutive model of the super-elastic material according to the geometric boundary condition and the material parameter; establishing the discrete balance equation of the isogeometric method; establishing the model motion equation and the nonlinear solution system; using the Newton iteration method to solve the nonlinear system, and mapping the displacement solution back to the subdivision, and the specific steps are as follows:

[0008] Step 1: converting the CC subdivision into a B-spline body model, and establishing the mapping relationship between the control points of the subdivision and the spline body, which is as follows:

[0009] Given a complex multi-piece hexahedral mesh model, denoted as H, the spline body of each hexahedral element is generated by the three-variable spline approximation method, and the spline body control points are H'. First, the coordinates of the inner points of the mesh v in are calculated, and the formula is as follows:

[0010]

[0011] Wherein υ is an arbitrary point in the hexahedral mesh, and the domain of the hexahedral mesh element where υ is located is L=v, e1, e2, e3, f1, f2, f3, c. Wherein e t , f t , c respectively represent the coordinates of the adjacent edge points, the adjacent face points and the adjacent block points of v. The adjacent edge points are the points on the same edge as υ, and there are 3 adjacent edge points in a hexahedral mesh element; the adjacent face points are the opposite vertices on the same face as υ, and there are 3 adjacent face points in a hexahedral; the adjacent block points are the opposite vertices in the same hexahedral mesh element as υ, and there is 1 adjacent block point in a hexahedral.

[0012] According to the calculated inner point coordinates, the new edge point, face point or corner point coordinates are obtained by averaging the adjacent inner point coordinates, that is,

[0013]

[0014] Wherein n is the number of hexahedral mesh elements containing the to-be-solved edge point, face point or corner point, is the coordinate of the inner point closest to the to-be-solved edge point, face point or corner point in the tth hexahedral mesh element containing the to-be-solved edge point, face point or corner point.

[0015] The above calculation formula is expressed by mapping matrix V, and the mapping relationship of two groups of control points is obtained:

[0016] H = VH' (3)

[0017] Step 2, according to the input material properties and boundary conditions, the stress-strain relationship is established, and the calculation form of the corresponding hyperelastic model is selected, which is as follows:

[0018] Input hyperelastic material properties, and establish stress-strain relationship according to material model. Unlike traditional linear elastic problems, the stress-strain of hyperelastic material is a nonlinear relationship, which has an elastic potential energy function related to strain. The model stress is obtained by deriving the strain from the potential energy function. Taking St. Venant-Kirchhoff hyperelastic material as an example, the corresponding elastic potential energy function is:

[0019] Ψ(E) = λtr(E) + μtr(E 2 ) (4)

[0020] In the formula, Ψ(E) is the size of elastic potential energy when the model is deformed, E represents the relative deformation size, i.e. strain, and λ and μ are input material coefficients. The stress component is obtained by taking the partial derivative of the elastic potential energy function with respect to the strain component, as follows:

[0021]

[0022] In the formula, S is the Kirchhoff stress tensor, and I is a unit matrix. According to the constitutive relationship of stress and strain, the derivative of stress with respect to strain is calculated to obtain the tangent material matrix D T :

[0023]

[0024] Strain refers to the relative deformation size of the object node, which can be represented as:

[0025]

[0026] Where represents the displacement gradient of the object.

[0027] The input boundary conditions are generally two kinds: displacement boundary condition (Dirichlet condition) and external force condition (Neumann condition). According to the two conditions, the node displacement values are constrained to ensure the uniqueness of the solution.

[0028] Step 3, establish discrete equilibrium equation by equal geometry method

[0029] The spline basis function input in step 1 is used as the basis function of the solution space by the isogeometric method, and the discrete points in the spline element are interpolated and solved to obtain the numerical solution under the accurate geometric modeling expression. First, the element balance equation is established. According to the expression of the spline basis function, the displacement of the discrete points in the element is linearly expressed by the displacement of each control point, as follows:

[0030] u(ξ)=∑ i N i u i =Nu e (8)

[0031] In the formula, u(ξ) represents the displacement size of any point of the model, u i represents the displacement size of the control point i, N i represents the basis function size of the control point i, and N is the matrix of the basis function combination. According to the displacement of the control point, the expression of the strain coefficient matrix B of the object can be derived from formula (7):

[0032] Bu e =E=LNu e (9)

[0033] L is a differential operator, and the above formula rewrites the strain into a tensor form. According to the Lagrange virtual work principle, the total virtual work of the system in the balanced state is zero, and the following element balance equation is established:

[0034] d e ∫ Ω B T SdΩ=d e ∫ Ω N T fdΩ+∫ Γ N T tdΓ (10)

[0035] In the formula, d e represents the virtual displacement, f represents the volume force, and t represents the area force. The left side of the equation represents the virtual work done by the internal force of the model on the virtual displacement, and the right side represents the virtual work done by the external force on the virtual displacement. The object remains in balance, so the two kinds of virtual work cancel each other out, and (10) is obtained. The element balance equations are assembled into the overall system balance equation, and after the virtual displacement is eliminated, the system balance equation is obtained. Since numerical calculation needs to be performed on the parameter domain, the above equation is subjected to Jacobian transformation, i.e., Ω→Ω′, to obtain the calculable system balance equation:

[0036] ∫ Ω′ B T S|J|dΩ′=∫ Ω′ N T f|J|dΩ′+∫ Γ′ N Tt|J|dΓ' (11)

[0037] where |J| is the Jacobian of the transformation from the physical domain to the parametric domain, and the above equation can also be written as:

[0038] Ku = f (12)

[0039] K is the stiffness matrix, f is the load vector, and u is the displacement solution to be solved. The displacement solution of the static analysis of the model can be obtained by solving the above equation.

[0040] Step 4, establish the model motion equation and use implicit time integration to establish the nonlinear solution system, as follows:

[0041] When dealing with the deformation process of the model, the changes of state variables and deformation over time must be considered. According to Newton's law, the general form of the object motion equation can be expressed as follows:

[0042]

[0043] where M represents the mass matrix, which is related to the volume density of the model; C is the damping matrix, which is related to the mass matrix and the stiffness of the model; R is the initial load vector, which is related to the initial deformation size of the model; and P is the size of the external force received by the model. and represent the second-order derivative and first-order derivative of displacement, i.e. acceleration and velocity, respectively.

[0044] After establishing the object motion equation, the Newmark implicit time integration is used for time discretization, and the velocity-displacement recursive relationship between time steps is established, which can be expressed as follows:

[0045]

[0046]

[0047] where Δt is the discretized time step, and the smaller the step size, the higher the solution accuracy; β and γ are self-defined parameters that determine the form of integral convergence. Equations (14) and (15) are the n-step displacement derived from the n+1-step velocity and acceleration. In the solution, the n-step displacement, velocity and acceleration are known, so the above equations need to be moved and combined to obtain a new derivation:

[0048]

[0049]

[0050] where Δu n+1 = u n+1 - u n represents the displacement difference between the previous and next two time steps, and α 1～6is the transformed parameter, and is a constant term. The solution of the next time step is determined by the solution of the previous time step and the displacement change, and the unknown quantity is only Δu n+1 By bringing (16) and (17) into the object motion equation (13), a nonlinear algebraic equation of unknown displacement is obtained:

[0051]

[0052] Equation (18) is an equation group composed of N equations of N time steps, and the displacement solution of each time step is obtained by solving the equation group. Then, each set of displacement solution is mapped back into the control points of the model, and the motion state and deformation degree of the spline body model of each time step are obtained.

[0053] Step 5, solving the nonlinear equation group and mapping the displacement solution of each time step back into the control points corresponding to the CC subdivision body:

[0054] For equation (18) of a certain time step, the Newton iteration method is used to solve the displacement change Δu n+1 The tangent stiffness matrix to be solved in each sub-iteration is as follows:

[0055]

[0056] The matrix can be obtained by deriving the balance equation (11) obtained in step 3. Then, the following scheme is executed for each Newmark time step iteration:

[0057]

[0058]

[0059] The initial value of the iteration is the convergence value of the last iteration: After the iteration converges, the displacement, velocity and acceleration of the time step are calculated and substituted into the next time step to continue solving, and the displacement solution of each time step is obtained. At this time, the displacement of each control point of the spline body is obtained, and then the subdivision body and the spline body control point mapping formula (3) in step 1 is used to finally obtain the motion state and deformation size of the subdivision body, and the motion simulation is completed.

[0060] The present application has the beneficial effects that the present application uses the isogeometric analysis method to realize the deformation simulation analysis of the hyperelastic model, and compared with the traditional finite element analysis, the present application has higher simulation precision and efficiency for the model with the same precision. The present application directly interfaces with the model after CC subdivision, realizes the deformation simulation and simulation of the complex CC subdivision body model, and perfects the integrated process from modeling to simulation. The present application constructs an isogeometric framework for numerical simulation of the CC subdivision body model, and the framework is also applicable to the structural optimization of the model.

[0061] The application converts the CC subdivision model into a spline body model with precise geometric representation, and then uses isogeometric analysis to simulate, directly connects the modeling and simulation processes, and realizes the integrated process of modeling and simulation. Moreover, for the simulation of large deformation of super-elastic materials, the finite element analysis often causes the precision to decrease due to the element distortion, or the solving time to be too long due to the dense mesh. The high precision and high continuity of isogeometry can well capture the deformation area, and only a small number of degrees of freedom is needed to show the complex deformation process, thereby providing a new idea for efficient and high-quality simulation. BRIEF DESCRIPTION OF DRAWINGS

[0062] Figure 1 a and Figure 1 b are different perspectives of the bridge hexahedral mesh model, respectively;

[0063] Figure 2 a is Figure 2 b are different perspectives of the bridge spline body model, respectively;

[0064] Figure 3 is the distribution of each spline body block of the bridge model;

[0065] Figure 4 is the result of the isogeometric analysis of the model by the application, Figure 3 different gray scales represent displacement solution sizes;

[0066] Figure 5 is the front result graph of the dynamic deformation of the bridge model under stress by the method of the application, Figure 5 a, 5b, 5c are three time points in the deformation process;

[0067] Figure 6 is the side result graph of the dynamic deformation of the bridge model under stress by the method of the application, Figure 6 a, 6b, 6c are three time points in the deformation process;

[0068] Figure 7 is the result graph of the dynamic deformation of the model C under stress by the method of the application, Figure 7 a, 7b, 7c are three time points in the deformation process;

[0069] Figure 8 is the result graph of the dynamic deformation of the model tentacle under stress by the method of the application, Figure 8 a, 8b, 8c are three time points in the deformation process;

[0070] Figure 9 is the overall flowchart. DETAILED DESCRIPTION

[0071] The application will be further described below with reference to the drawings.

[0072] A method for simulating isogeometric deformation of CC subdivision body, the specific steps are as follows:

[0073] 1. Establish the mapping relationship between the control points of the subdivision body and the spline body, specifically as follows

[0074] Given a hexahedral mesh model (such as Figure 1 a and 1b), record its corresponding control mesh as H, and generate an approximated spline body (such as Figure 2 a and 2b) for each hexahedral element by three-variable spline approximation method, and the spline body control points are H', and the calculation method is as follows:

[0075]

[0076] Where υ is any point in the hexahedral mesh, and the domain of the hexahedral mesh element where υ is located is L=v, e1, e2, e3, f1, f2, f3, c. Where e t , f t , c represent the coordinates of the adjacent edge points, adjacent face points and adjacent block points of v. The adjacent edge points are the points on the same edge as υ, and there are 3 adjacent edge points in a hexahedral mesh element; the adjacent face points are the opposite vertices on the same face as u, and there are 3 adjacent face points in a hexahedral; the adjacent block points are the opposite vertices in the same hexahedral mesh element as υ, and there is 1 adjacent block point in a hexahedral.

[0077] According to the calculated inner point coordinates, the new edge point, face point or corner point coordinates are obtained by averaging the adjacent inner point coordinates, that is

[0078]

[0079] Where n is the number of hexahedral mesh elements containing the to-be-solved edge point, face point or corner point, is the coordinates of the inner point closest to the to-be-solved edge point, face point or corner point in the tth hexahedral mesh element containing the to-be-solved edge point, face point or corner point.

[0080] The above calculation formula is expressed by a mapping matrix V to obtain the mapping relationship between the control points of the two:

[0081] H=VH′ (3)

[0082] 2. Input material properties and boundary conditions, and establish the constitutive model of hyperelastic material:

[0083] Input the hyper-elastic material properties, and establish the stress-strain relationship according to the material model. Unlike the traditional linear elastic problem, the stress-strain of hyper-elastic material is nonlinear, and there is an elastic potential energy function related to the strain, and the model stress is obtained by deriving the strain from the potential energy function. Take the St. Venant-Kirchhoff hyper-elastic material as an example, and the corresponding elastic potential energy function is:

[0084] Ψ(E) = λtr(E) + μtr(E 2 ) (4)

[0085] In the formula, Ψ(E) is the elastic potential energy size when the model is deformed, E represents the relative deformation size, that is, the strain, and λ and μ are the input material coefficients. The stress component is obtained by taking the partial derivative of the elastic potential energy function (1) with respect to the strain component, as follows:

[0086]

[0087] In the formula, S is the Kirchhoff stress tensor, and I is a unit matrix. According to the constitutive relationship of stress and strain, the derivative of stress with respect to strain is calculated to obtain the tangent material matrix D T :

[0088]

[0089] The strain indicates the relative deformation size of the object node, which can be expressed as:

[0090]

[0091] Where represents the displacement gradient of the object.

[0092] The input boundary conditions are generally two kinds, displacement boundary condition (Dirichlet condition) and external force condition (Neumann condition), which constrain the node displacement values according to the two conditions to ensure the uniqueness of the solution.

[0093] 3. Establish discrete equilibrium equations using the isogeometric method

[0094] The isogeometric method is used to input the spline basis function in step 1 as the basis function of the solution space, and the basis function is used to interpolate the discrete points in the spline element to obtain the numerical solution under the accurate geometric modeling expression. First, establish the element balance equation. According to the representation of the spline basis function, the displacement of the element discrete point is linearly represented by the displacement of each control point, as follows:

[0095] u(ξ) = ∑ i N i u i = Nu e (8)

[0096] where u(ξ) represents the displacement of any point of the model, u i represents the displacement of control point i, N i represents the basis function of control point i, and N is the matrix of basis function combination. According to the displacement of control points, the strain coefficient matrix B of the object can be derived from equation (7) as follows:

[0097] Bu e = E = LN u e (9)

[0098] L is a differential operator, and the above equation rewrites the strain into a tensor form. According to the Lagrange's virtual work principle, the total virtual work done by the system in the equilibrium state is zero, and the following element balance equation is established:

[0099] d e ∫ Ω B T S dΩ = d e ∫ Ω N T f dΩ + ∫ Γ N T t dΓ (10)

[0100] where d e represents the virtual displacement, f represents the volume force, and t represents the area force.

[0101] The left side of the equation represents the virtual work done by the internal force of the model on the virtual displacement, and the right side represents the virtual work done by the external force on the virtual displacement. The object remains in equilibrium, so the two kinds of virtual work cancel each other out, and equation (10) is obtained. The system balance equation is obtained by assembling the balance equations of each element into the overall system balance equation and eliminating the virtual displacement. Since numerical calculation needs to be performed on the parameter domain, the above equation is subjected to Jacobian transformation, i.e., Ω→Ω', to obtain the calculable system balance equation:

[0102] ∫ Ω′ B T S |J| dΩ' = ∫ Ω′ N T f |J| dΩ' + ∫ Γ′ N T t |J| dΓ' (11)

[0103] where |J| is the Jacobian determinant of the transformation between the physical domain and the parameter domain, and the above equation can also be simply written as:

[0104] Ku = f (12)

[0105] K is the stiffness matrix, f is the load vector, and u is the displacement solution to be solved. The displacement solution of the static analysis of the model can be obtained by solving the above equation, as shown in equation (13). Figure 4 ​

[0106] 4. Establish the model motion equation and use implicit time integration to establish the nonlinear solution system, as follows

[0107] When processing the model deformation process, the changes of state variables and deformation over time must be considered. According to Newton's law, the general form of the object motion equation can be expressed as follows:

[0108]

[0109] where M represents the mass matrix, which is related to the volume density of the model; C is the damping matrix, which is related to the mass matrix and the stiffness of the model; R is the initial load vector, which is related to the initial deformation size of the model; and P is the size of the external force received by the model. respectively represent the second-order derivative and the first-order derivative of displacement, i.e., acceleration and velocity.

[0110] After establishing the object motion equation, the Newmark implicit time integration is used for time discretization, and the velocity-displacement recursive relationship between time steps is established, which can be expressed as follows:

[0111]

[0112]

[0113] where Δt is the discretized time step, and the smaller the step size, the higher the solution accuracy; β and γ are self-defined parameters that determine the form of integral convergence. Equations (14) and (15) are derived from the n+1 step of velocity and acceleration to the n step of displacement. In the solution, the n step of displacement, velocity, and acceleration are known, so the above equations need to be moved and combined to obtain a new derivation:

[0114]

[0115]

[0116] where Δu n+1 = u n+1 -u n represents the displacement difference between the previous and next two time steps, α 1～6 is the transformed parameter, which is a constant term. Thus, the solution of the next time step is determined by the solution of the previous time step and the displacement change, and the only unknown quantity is Δu n+1 . By bringing (16) and (17) into the object motion equation (13), the nonlinear algebraic equation of unknown displacement is obtained:

[0117]

[0118] Equation (18) is a set of equations composed of N equations at N time steps, by solving the equation set, the displacement solution of each time step is obtained. Then each set of displacement solution is mapped back to the control points of the model, to obtain the motion state and deformation degree of the spline body model at each time step, as shown in Figure 5 a, 5b, 5c.

[0119] 5, solving the nonlinear equation set and mapping the displacement solution of each time step to the corresponding nodes of the CC subdivision body

[0120] For equation (18) of a time step, the displacement change Δu is solved by using Newton iteration method n+1 The tangent stiffness matrix needed to be solved in each sub-iteration is as follows:

[0121]

[0122] The matrix can be obtained by deriving the equilibrium equation obtained in step 3. Then, the following scheme is executed for each Newmark time step iteration:

[0123]

[0124]

[0125] Wherein, the initial value of iteration is the convergence value of the last iteration: After the iteration converges, the displacement, velocity and acceleration of the time step are calculated, which are substituted into the next time step to continue to solve, to obtain the displacement solution of each time step. At this time, the displacement of each control point of the spline body is obtained, and then the subdivision body and spline body control point mapping formula (3) in step 1 is used, to finally obtain the motion state and deformation size of the subdivision body, and the motion simulation is completed.

[0126] The implementation is as follows: Figure 5 The result graph of dynamic deformation of the bridge model under external force is shown in Figure 7 The result graph of dynamic deformation of the letter C is shown in Figure 8 The result graph of dynamic deformation of the tentacle model is shown in. The present application can use equal geometry to numerically simulate and simulate complex multi-piece models, to solve the motion state of the model, and the present method can be solved faster than the present method under the same precision.

Claims

1. A geometric analysis simulation method based on a hyperelastic material model and volume subdivision, characterized in that, Includes the following steps: S1 converts the CC subdivision mesh model into a spline model and establishes the control point mapping relationship between the CC subdivision and the spline. S2 establishes a constitutive model reflecting the stress-strain relationship of hyperelastic materials based on geometric boundary conditions and material parameters; S3 uses the constitutive model described in step 2 to establish discrete equilibrium equations using the isogeometric method. By solving the discrete equilibrium equations, the displacement solution of the spline model static analysis is obtained. S4 establishes the model's equations of motion, discretizes time using Newmark implicit time integration, establishes the velocity-displacement recursive relationship between time steps, and substitutes the velocity-displacement recursive relationship between time steps into the model's equations of motion to obtain the nonlinear solution system. S5 uses Newton's iteration method to solve the nonlinear system, obtains the displacement solution at each time step, maps the displacement solution at each time step back to the control points corresponding to the CC subdivision, obtains the motion state and deformation magnitude of the subdivision, and completes the motion simulation. In step S3, the establishment of the discrete equilibrium equations using the isogeometric method includes the following sub-steps: Use the spline basis functions of the spline model generated in step 1 as the basis functions of the solution space. The element equilibrium equations are established, and the displacements of the discrete points of the element are linearly represented by the displacements at each control point, as follows: u(ξ)=∑ i N i you i =Now e (8) In the formula, u(ξ) represents the magnitude of the displacement at any point in the model, u i N represents the magnitude of the displacement of control point i. i Let N represent the magnitude of the basis functions at control point i, and N be the matrix of the basis function combination. Based on the displacement of the control points, derive the representation of the strain coefficient matrix B of the object: Bu e =E=LNu e (9) Where L is the differential operator, the above equation rewrites the strain in tensor form; according to the Lagrange principle of virtual work, the total virtual work performed on the system in equilibrium is zero, and the following unit equilibrium equation is established: d e ∫ Ω B T SdΩ=d e ∫ Ω N T fdΩ+∫ Γ N T tdΓ (10) Where, d e Let f represent the virtual displacement, f represent the volume force, and t represent the area force. The left side of the equation represents the virtual work done by the internal forces of the model on the virtual displacement, and the right side represents the virtual work done by the external forces on the virtual displacement. The object remains in equilibrium, and the two virtual works cancel each other out. Assemble the equilibrium equations of each unit into the equilibrium equations of the overall system. After eliminating virtual displacements, the system equilibrium equations are obtained. Then, through the Jacobian transformation, i.e., Ω→Ω ′ The computable system equilibrium equations are obtained as discrete equilibrium equations: ∫ Ω′ BTS|J|dΩ′=∫ Ω′ N T f|J|dΩ′+∫ Γ′ N T t|J|dΓ′ (11) In the formula, |J| is the Jacobian determinant of the transformation between the physical domain and the parameter domain, and the above formula can be simplified as follows: Ku = f (12) K is the stiffness matrix, f is the load vector, and u is the displacement solution to be found. Solving the discrete equilibrium equations includes the following sub-steps: The discrete points within the element are interpolated using spline basis functions to obtain a numerical solution under the accurate geometric modeling expression, which serves as the displacement solution for the static analysis of the spline model.

2. The geometric analysis and simulation method for hyperelastic materials based on volume subdivision as described in claim 1, characterized in that, In S1, The mesh model of the CC subdivision is a complex, multi-layered, hexahedral CC subdivision mesh model; The spline model is a spline model of hexahedral elements generated by the three-variable spline approximation method; The control point mapping relationship is formed in the following way: Given a complex, multi-layered hexahedral mesh model, denoted as H, the spline body of each hexahedral element is generated using a three-variable spline approximation method. The spline body control points are H. ′ The calculation method is as follows: Where υ is any point within the hexahedral mesh, and the neighborhood of the hexahedral mesh element containing υ is L = υ, e1, e2, e3, f1, f2, f3, c; where e t f t , c represents the coordinates of adjacent edge points, adjacent face points, and adjacent block points of v, respectively; an adjacent edge point is a point on the same edge as v, and there are 3 adjacent edge points in a hexahedral mesh cell; an adjacent face point is a vertex opposite v on the same face, and there are 3 adjacent face points in a hexahedron; an adjacent block point is a diagonal point of v in the same hexahedral mesh cell, and there is 1 adjacent block point in a hexahedron. Based on the calculated interior point coordinates, new edge point, face point, or corner point coordinates are obtained by averaging the coordinates of adjacent interior points. Where n is the number of hexahedral mesh elements containing the edge points, face points, or corner points to be determined. Let V be the coordinates of the interior point on the t-th hexahedral mesh element containing the edge, face, or corner point to be determined, which is closest to that edge, face, or corner point. Representing the above calculation formula using a mapping matrix V, we obtain the mapping relationship between two sets of control points: H=VH′ (3).

3. The geometric analysis and simulation method for hyperelastic materials based on volume subdivision as described in claim 1, characterized in that, In S2, The superelastic material is: St. Venant-Kirchhoff superelastic material. The material parameters include the elastic potential energy function: Ψ(E)=λtr(E)+μtr(E 2 ) (4) In the formula, Ψ(E) is the magnitude of the elastic potential energy when the model undergoes deformation, E represents the magnitude of the relative deformation, i.e. strain, and λ and μ are constants obtained from the input material coefficients. The stress components are obtained by taking the partial derivative of the elastic potential energy function with respect to the strain components, as follows: In the formula, S is the Kirchhoff stress tensor, and I is an identity matrix; Based on the constitutive relation of stress and strain, the tangent material matrix D is obtained by calculating the derivative of stress with respect to strain. T : Strain refers to the relative deformation of the nodules of an object, and can be expressed as: in Represents the displacement gradient of an object; The geometric boundary conditions include: Boundary conditions reflecting displacement: Dirichlet conditions Boundary conditions reflecting external forces: Neumann conditions. By constraining the nodal displacement values ​​using the two boundary conditions described above, a unique solution is obtained.

4. The geometric analysis and simulation method for hyperelastic materials based on volume subdivision as described in claim 3, characterized in that, In S4, The model's equations of motion are in the following form: Where M represents the mass matrix, which is related to the volume density of the model; C is the damping matrix, which is related to the mass matrix and the stiffness of the model; R is the initial load vector, which is related to the initial deformation of the model; and P is the magnitude of the external force acting on the model. These represent the second derivative and first derivative of the displacement, respectively, and represent acceleration and velocity, respectively. The method for discretizing time using the Newmark implicit time integral and establishing the velocity-displacement recursive relationship between time steps specifically includes the following steps: Using Newmark implicit time integration to discretize time, a recursive relationship between velocity and displacement is established between time steps, expressed in the following form: Where Δt is the discretized time step, and the smaller the step, the higher the accuracy of the solution; β and γ are user-defined parameters that determine the form of integral convergence; Formulas (14) and (15) derive the displacement of the nth step from the (n+1)th step of velocity and acceleration. In the solution, the displacement, velocity, and acceleration of the nth step are known, so it is necessary to rearrange and combine the terms of the above formula to obtain a new derivation: Where, Δu n+1 =u n+1 -u n α represents the displacement difference between two consecutive time steps. 1～6 These are the transformed parameters, which are constant terms; The method for obtaining the nonlinear solution system by substituting the velocity-displacement recursive relationship between time steps into the model's equations of motion specifically includes the following steps: The solution for the next time step is determined by the solution for the previous time step and the change in displacement. By substituting (16) and (17) into the equation of motion of the object (13), the nonlinear algebraic equation for the unknown displacement is obtained: Formula (18) is a system of N equations consisting of N time steps and N equations.

5. The geometric analysis and simulation method for hyperelastic materials based on volume subdivision as described in claim 4, characterized in that, S5 specifically includes the following steps: Solving for the displacement change Δu using Newton's iterative method n+1 The tangent stiffness matrix to be solved in each sub-iteration is as follows: This matrix can be obtained by differentiating the equilibrium equations obtained in step S3; For each Newmark time step, the following scheme is executed iteratively: The initial value of the iteration is the convergence value of the previous iteration: After the iteration converges, the magnitudes of displacement, velocity, and acceleration at that time step are calculated and substituted into the next time step to continue solving, obtaining the displacement solution for each time step; At this point, the displacement magnitudes at each control point of the spline are obtained. Then, using the subdivision and spline control point mapping formula (3) in step S1, the motion state and deformation magnitude of the subdivision are finally obtained, and the motion simulation is completed.

Images