Collision simulation-oriented co-rotation embedded domain isogeometric explicit dynamics method
By introducing co-rotation embedding domain and co-rotation method in isogeometric analysis methods, the problems of low calculation accuracy and difficult parameterization of complex three-dimensional models in collision simulation are solved, and efficient nonlinear dynamic analysis is achieved.
Patent Information
- Application Number
- CN202510042237.4
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2025-01-10
- Publication Date
- 2025-05-09
AI Technical Summary
Existing geometric analysis methods have problems such as low calculation accuracy, prone to locking and difficulty in geometric parameterization when dealing with collision simulation of complex three-dimensional models.
A geometric explicit dynamics method such as co-rotating embedding domain is proposed. By embedding complex three-dimensional models into regular background embedding domain geometry, ray tracing and collision detection are used to determine relative positions, nonlinear explicit dynamics are directly used to solve the background embedding domain basis function, and a co-rotating coordinate system is established at the Gaussian point of the background embedding domain unit, and a co-rotating method is used to perform accurate kinematic descriptions.
This method can effectively deal with the nonlinear dynamics problem of any complex three-dimensional model, significantly improves the calculation accuracy and efficiency, slows down the occurrence of locking phenomena, and reduces the difficulty of solving.
Smart Images
Figure CN119962301A_ABST
Abstract
Description
Technical Field
[0001] The invention belongs to the field of computational mechanics and relates to a geometric explicit dynamics method such as co-rotation embedded domain for collision simulation. Background Art
[0002] Vehicle collision simulation, as an important application field of Computer-Aided Engineering (CAE), has been widely used in the automotive design stage. Its core task is to solve large-scale transient strong nonlinear dynamic problems through numerical methods, especially in the collision process involving complex conditions such as geometric nonlinearity, material nonlinearity, and contact nonlinearity. The explicit finite element analysis method has become an important tool in this field due to its superiority in dealing with strong nonlinear problems. However, traditional explicit finite element analysis methods face many challenges in solving these nonlinear dynamic problems. For example, the single-point quadrature method of low-order elements often leads to the "hourglass phenomenon", which distorts the calculation results; high-order elements have significant high-order modal errors and poor calculation efficiency. In addition, the geometric errors introduced in the meshing process seriously affect the stability and accuracy of nonlinear solutions.
[0003] Isogeometric Analysis (IGA) avoids the meshing process by directly using non-uniform rational B-splines (NURBS) widely used in computer-aided design (CAD), realizes the seamless integration of CAD and CAE, and provides higher geometric accuracy and solution accuracy. Compared with the traditional finite element method, IGA can effectively reduce high-order modal errors, solve the problem of negative mass vertices, and has higher computational efficiency and robustness in explicit dynamics solutions. The full Lagrangian and updated Lagrangian methods have low computational accuracy and are prone to locking when dealing with large geometric nonlinear large deformations dominated by large rotations. For complex three-dimensional solid structures, the common B-Rep representation method makes it difficult to represent the internal geometric information of the solid, and only the surface information is available; the geometric parameterization of complex three-dimensional models is extremely difficult, and it is difficult to directly analyze complex models using traditional isogeometric analysis.
[0004] On the basis of the above-mentioned background technology, the present invention proposes a new co-rotating embedded domain isogeometric explicit dynamics method for collision simulation. This method is based on the idea of geometric analysis such as embedded domain, embeds a complex three-dimensional model into a regular background embedded domain geometry, uses ray tracing and collision detection schemes to determine the relative position of the background embedded domain unit and the model, and then directly uses the background embedded domain basis function to directly perform the corresponding nonlinear explicit dynamics solution, avoiding the grid division and complicated geometric parameterization process. This method innovatively establishes a co-rotating coordinate system at the Gaussian point of the background embedded domain unit, and uses the co-rotation method for accurate kinematic description, which greatly reduces the occurrence of locking phenomenon and improves calculation accuracy and efficiency. In addition, this method proposes a new processing solution for special cases such as collision with rigid wall surfaces, which reduces the difficulty of solution and improves calculation efficiency. Summary of the invention
[0005] The purpose of the present invention is to overcome the shortcomings of existing isogeometric analysis methods in processing collision simulations and to provide a co-rotating embedded domain isogeometric explicit dynamics analysis method that can process any complex three-dimensional models.
[0006] The technical solution of the present invention:
[0007] A co-rotating embedded domain isogeometric explicit dynamics method for collision simulation includes the following steps:
[0008] Step 1: Read the model file:
[0009] Read the STEP file of the 3D geometric model and the BREP file of the boundary condition application surface;
[0010] Step 2: Define the background embedding domain:
[0011] The background embedding domain geometry is usually defined by three-parameter NURBS, which is a regular hexahedron with a size that can accommodate the model. The H refinement and P refinement methods are used to perform relevant refinement and order-raising operations on the NURBS basis functions of the background embedding domain geometry to meet the accuracy requirements of dynamic finite element simulation and obtain a number of background embedding domain units.
[0012] Step 3: Determine the relationship between the background embedding domain unit and the 3D geometric model:
[0013] Ray tracing and collision detection schemes are used to interpret three types of positional relationships between background embedded domain units and 3D geometric models: intersection, interior, and exterior, and the corresponding intersection units, interior domain units, and exterior domain units are obtained;
[0014] Step 4: Derive the equations of motion for the embedded domain IGA:
[0015] The corresponding embedded domain IGA motion equations are derived. When the basis function order of the background embedded domain geometry is low, the mass matrix and internal force matrix of the out-of-domain unit are 0, the in-domain unit is solved by standard Gaussian integral, and the intersecting unit is processed by combined Gaussian integral. If the basis function order of the background embedded domain geometry is high, the mass matrix and control point force matrix of the out-of-domain unit need to be multiplied by a decimal to impose a penalty, reducing the contribution of the out-of-domain unit to the embedded domain IGA motion equations. The processing of the in-domain and intersecting units is consistent with that when the basis function order of the background embedded domain geometry is low.
[0016] Step 5: Stress update based on co-rotation method
[0017] By establishing a co-rotating coordinate system that rotates with the Gaussian points of the background embedded domain unit to accurately describe the nonlinear motion, the rotated Cauchy stress and deformation rate tensors defined in the co-rotating coordinate system are used as stress and strain metrics, and the corresponding stress update is performed under this theoretical framework to study the nonlinear motion of linear elastic materials and elastoplastic materials based on the von Mises yield criterion.
[0018] Step 6: Contact Nonlinear Formulation
[0019] In each time step, the contact of the objects is not considered first, and the trial configuration after the initial movement is obtained. If the non-embedded condition is met, the trial configuration is the real configuration, and there is no contact between the objects; if the non-embedded condition is not met, the penalty function method is used to calculate the normal contact force to make corresponding corrections to the force at the contact point to meet the non-embedded condition. For special cases such as rigid wall contact, the normal contact acceleration and velocity are directly processed to meet the contact condition.
[0020] Step 7: Discretization in time domain
[0021] The central difference method is used to explicitly discretize the IGA motion equations in time domain.
[0022] Beneficial effects of the invention: The invention provides a new type of co-rotation embedded domain explicit dynamics analysis method that can perform isogeometric analysis on any complex model without geometric reconstruction. Compared with the traditional isogeometric analysis explicit dynamics method, the proposed method can directly analyze the nonlinear dynamics problem of any complex model. Compared with the full Lagrangian and updated Lagrangian methods, the co-rotation method can more accurately describe the nonlinear kinematic model, significantly improving the calculation accuracy and engineering efficiency. BRIEF DESCRIPTION OF THE DRAWINGS
[0023] Figure 1 is a schematic diagram of a geometric explicit dynamics method such as co-rotational embedded domain for collision simulation according to the present invention;
[0024] Figure 2It is the result of the collision detection + ray tracing method described in the present invention;
[0025] Figure 3 is the initialization flow chart of the present invention;
[0026] Figure 4 It is the three-dimensional B-Rep model and the generated octree integration unit; wherein, (a) is the B-Rep model of three-dimensional geometry, and (b) is the schematic diagram of octree space partition;
[0027] Figure 5 It is a schematic diagram of the impact on a finite rigid wall surface described in the present invention. DETAILED DESCRIPTION
[0028] The specific implementation of the present invention is further described below in conjunction with the accompanying drawings and technical solutions.
[0029] A co-rotating embedded domain isogeometric explicit dynamics method for collision simulation, the flow chart is as follows Figure 1 As shown, the following steps are included:
[0030] (1) Read model geometry information:
[0031] Read the STEP file of the 3D geometric model to obtain the minimum and maximum point coordinates of the axis-aligned bounding box of the 3D geometric model; read the BREP file of the boundary condition application surface to obtain the geometric parameter information of the boundary condition application surface, such as node vectors, control points, and control point weights;
[0032] (2) Define the background embedding domain:
[0033] According to the coordinates of the minimum and maximum points of the axis-aligned bounding box of the three-dimensional geometric model, the background embedding domain geometry is defined by using the three-parameter non-uniform rational B-spline (NURBS) in the form of tensor product; the shape of the background embedding domain geometry is a regular hexahedron, and its size can accommodate the entire three-dimensional geometric model; the NURBS basis function of the background embedding domain geometry (hereinafter referred to as the background embedding domain basis function) is refined and upgraded by H refinement and P refinement to meet the accuracy requirements of dynamic finite element simulation, and a number of background embedding domain units are obtained;
[0034] (3) Determine the relationship between the background embedding domain unit and the 3D geometric model:
[0035] There are three types of positional relationships between background embedded domain units and three-dimensional geometric models: intersection, interior, and exterior. The positional relationship can be determined by using ray tracing and collision detection: collision detection determines whether the background embedded domain unit intersects with the three-dimensional geometric model. If not, the point position judgment method is used to determine whether the maximum physical coordinate point and the minimum physical coordinate point of the background embedded domain unit are within the three-dimensional geometric model, thereby determining the internal and external relationship between the embedded domain unit and the model. Among them, the unit that intersects with the three-dimensional geometric model is hereinafter referred to as the intersection unit, the background embedded domain unit located inside the three-dimensional geometric model is hereinafter referred to as the in-domain unit, and the background embedded domain unit located outside the three-dimensional geometric model is hereinafter referred to as the out-of-domain unit. Ray tracing is used to determine the internal and external relationship between the background embedded domain unit and the B-Rep model, and the corresponding judgment results are shown as follows. Figure 2 The corresponding initialization flowchart is as follows Figure 3 shown.
[0036] (4) Derivation of IGA motion equations
[0037] Starting from the weak form of the momentum conservation equation and based on the Galerkin principle, the imaginary power equation is obtained:
[0038]
[0039] In the formula, δv is the virtual velocity, t is the surface force acting on the boundary Γ, b is the body force acting on the unit mass of the object, ρ is the density, is the second-order derivative of displacement with respect to time, σ is the Cauchy stress. δD is the virtual deformation rate tensor, defined as:
[0040]
[0041] The IGA motion equation is spatially discretized using the background embedding domain basis function, and the virtual power equation is written as:
[0042]
[0043] Where M is the mass matrix, f int and f ext are the internal force and external force of the control point; the specific expressions are as follows:
[0044] M=∫ Ω ρN T NdΩ
[0045] f int =∫ V B T ·σdV
[0046] f ext =∫ V ρN·bdV+∫ Γ N·tdΓ
[0047] Where N represents the background embedding domain basis function, and B is the shape function matrix;
[0048] B=[B1...B A ...B n ]
[0049]
[0050] When the order of the background embedding domain basis function is quadratic or cubic, the cells outside the domain do not participate in the calculation directly, and only the cells inside the domain and the intersecting cells participate in the calculation. Among them, the cells inside the domain are normally Gaussian integrated, and the intersecting cells are integrated and solved by the combined Gaussian integral scheme. Using the octree space partitioning method, an intersecting cell is subdivided into eight integral sub-cells. For the background embedding domain basis functions of B of p, q, and s times in three directions, we first start from the 0th layer integration unit. If the integration unit is not clipped and is located inside the integration domain, (p+1)×(q+1)×(s+1) Gaussian points are selected in the unit for integration. If the integration unit is clipped by the boundary of the three-dimensional geometric model, it is evenly divided along the three directions of its parameter space. The subunit is divided into eight integration subunits. For each subunit completely inside the model, (p+1)×(q+1)×(s+1) Gaussian points are selected for integration. This process is recursively performed on all clipped units until the number of recursions reaches the target value to achieve relatively accurate Gaussian integration for the clipped units. The relevant process results are shown in Figure 4 As shown in the figure; when the order of the background embedding domain basis function is fourth order or above, the out-of-domain unit needs to multiply by 10 when solving the mass matrix and the internal and external forces of the control point. -5 ≤α≤10 -10 Penalize the domain to reduce the contribution of the out-of-domain units to the virtual power equation, and keep the intersection units, the in-domain units and the background embedded domain basis functions consistent when the order is low;
[0051] (5) Stress update based on co-rotation method
[0052] The co-rotation method requires corresponding stress and strain measurements in a co-rotation coordinate system; the co-rotation coordinate system is a local Cartesian coordinate system defined on the Gaussian points of the background embedded domain unit and rotates with the Gaussian points on the background embedded domain unit. The definition process is as follows:
[0053] The covariant basis vectors g1 and g2 are expressed as:
[0054]
[0055] Where x = (x, y, z) Tare the spatial coordinates of the material point measured in the global coordinate system, (ξ,η,ζ) are the parameter coordinates of the Gaussian points of the background embedding domain unit involved in the calculation; the covariant basis vectors g1 and g2 define a surface tangent to the layer ζ=d,d∈(-1,1), and the normal vector perpendicular to the surface is:
[0056]
[0057] Define a set of auxiliary vectors:
[0058]
[0059] Get the other two mutually perpendicular basis vectors of the co-rotational coordinate system:
[0060]
[0061] From this, we can define the rotation matrix R, whose matrix components R ij It can be expressed as:
[0062]
[0063] In the formula, the value ranges of the right subscripts i and j are 0, 1, and 2, corresponding to x, y, and z, respectively, that is, e0 = e x 、e1=e y 、e2=e z , e x Represents the unit vector in the x direction of the global coordinate system, e y Represents the unit vector in the y direction of the global coordinate system, e z Represents the unit vector in the z direction of the global coordinate system.
[0064] The transformation relationship between velocity, displacement, Cauchy stress and deformation rate tensor in the co-rotation coordinate system and the global coordinate system is:
[0065]
[0066] The deformation rate tensor defined in the corotational coordinate system and the rotational Cauchy stress As an objective measure of strain rate and stress, it does not change with the rotation or translation of the background embedded domain unit; therefore, in the co-rotating coordinate system, the rotational Cauchy stress update rule is:
[0067]
[0068] Where Δt is the time step, is the Cauchy stress rate of rotation.
[0069] For linear elastic materials, the rotational Cauchy stress rate is It is expressed as:
[0070]
[0071] In the formula, is the strain matrix in the co-rotational coordinate system, is the local constitutive matrix defined in the corotational coordinate system.
[0072] For the nonlinear plastic deformation of isotropic metal materials, the radial return algorithm based on the von Mises yield criterion can be used to update the corresponding stress, which includes two steps: first, according to t n The converged solution at the moment is used to calculate the trial stress. If the trial stress is within the linear elastic region of the material deformation, it is regarded as the true stress; otherwise, it is projected to the nearest point on the yield surface of the material.
[0073] (6) Contact nonlinearity formula
[0074] Objects A and B interact in the framework of the principles of continuum mechanics and must comply with the non-embedding condition, namely:
[0075] Ω A ∩Ω B =0
[0076] Where Ω A represents the space occupied by object A, Ω B Represents the space occupied by object B.
[0077] When dealing with nonlinear contact problems, the explicit dynamics analysis method first ignores the contact between objects in each time step, and independently calculates the IGA motion equation of each object to obtain the trial configuration of each object after preliminary motion; if the trial configuration satisfies the non-embedded condition, it indicates that there is no contact between the objects, the motion between the objects does not affect each other, and the trial configuration is the real configuration; if the trial configuration does not meet the non-embedded condition, it indicates that contact occurs between the objects, and the contact force needs to intervene in the calculation and correct the force at the contact point to meet the contact condition.
[0078] The contact force is calculated using the penalty function method: in each time step, it is first detected whether the point on object A is embedded in the surface of object B. If there is no penetration, no processing is required; otherwise, a normal contact force similar to the "normal spring" effect is applied between the contact point of object A and the contact surface of object B to limit penetration.
[0079] The normal contact force is:
[0080] f c =p c ×l×N
[0081] Where l is the penetration distance and N is the normal direction of the contact point at the contact surface.c represents the normal stiffness, and its calculation formula is:
[0082]
[0083] Where m c is the mass at the point, Δt is the time step;
[0084] For simple nonlinear contact problems such as impacting a rigid wall, separate treatment is performed accordingly: Figure 5 As shown, if the rigid wall is a finite rectangle with side lengths L and M, l and m are the unit vectors of the two adjacent sides of the rigid wall, It is represented as the vector of the point on the contact detection surface and the origin; after each time step Δt of the global update of velocity and acceleration, all points k on the surface of each object need to be checked to ensure that all points on the surface of the object are located inside the rigid wall by checking whether two inequalities are satisfied:
[0085]
[0086] For the case of an infinite rigid wall, skip the above steps and only need to check the penetration condition to see whether the contact detection point penetrates the rigid wall surface, that is, you need to check the following formula:
[0087]
[0088] Where n is the normal vector of the contact point. If penetration occurs, corresponding processing is required to meet the non-embedded condition: the normal contact acceleration of the contact point where penetration has occurred is set to the negative value of the original normal acceleration, and the normal velocity of the contact point is set to 0 to correct the force at the contact point;
[0089] (7) Time domain discretization
[0090] The central difference method is used to explicitly discretize the IGA motion equations in the time domain. First, the relevant parameters are initialized, namely the initial velocity, initial stress, and total simulation time; the mass matrix is solved; the Dirichlet boundary conditions are applied; and then the central difference method is used for time iteration. Finally, the corresponding result configuration is output.
[0091] The present invention provides a novel co-rotation embedded domain explicit dynamic analysis method that does not require geometric reconstruction and can perform iso-geometric analysis on any complex model. Compared with the traditional iso-geometric nonlinear explicit dynamic analysis method, the proposed method can directly perform nonlinear dynamic analysis on any complex three-dimensional B-Rep model, and the corresponding solution can be performed without meshing and complex geometric parameterization. For large-scale transient strong nonlinear problems such as collisions, the newly proposed co-rotation method for three-dimensional solid elements can effectively deal with geometric nonlinear problems dominated by large rotations and large deformations, and proposes a stress update scheme for elastoplastic materials under the theoretical framework of the co-rotation method. The penalty function method is used to deal with contact nonlinear problems, and a new processing method is proposed for special contact situations such as impacting rigid walls, which greatly improves the calculation accuracy and efficiency of large-scale transient strong nonlinear simulation problems such as collision simulations.
[0092] Although the embodiments of the present invention have been described for the purpose of illustration, it will be understood by those skilled in the art that various modifications, additions and substitutions may be made in form and detail without departing from the scope and spirit of the invention disclosed in the appended claims, and all these changes shall fall within the scope of protection of the appended claims of the present invention, and the various steps in the various departments and methods of the products claimed by the present invention may be combined together in any combination. Therefore, the description of the embodiments disclosed in the present invention is not intended to limit the scope of the present invention, but is used to describe the present invention. Accordingly, the scope of the present invention is not limited by the above embodiments, but is defined by the claims or their equivalents.
Claims
1. A co-rotating embedded domain isogeometric explicit dynamics method for collision simulation, characterized in that: The following steps are involved: (1) Read the model file: (2) Define the background embedding domain: (3) Determine the relationship between the background embedding domain unit and the 3D geometric model: (4) Derivation of the equations of motion for the embedded domain IGA; (5) Accurate kinematic description based on the co-rotation method and the corresponding stress update algorithm; (6) Derivation of nonlinear contact formulas for embedded domain and other geometric contacts and contact algorithms for rigid wall collisions; (7) Time domain discretization.
2. The co-rotating embedded domain isogeometric explicit dynamics method for collision simulation according to claim 1, characterized in that: The specific steps of step (1) are as follows: Read the STEP file of the three-dimensional geometric model to obtain the minimum and maximum point coordinates of the axis-aligned bounding box of the three-dimensional geometric model; read the BREP file of the boundary condition application surface to obtain the geometric parameter information of the boundary condition application surface, including node vectors, control points, and control point weights.
3. The co-rotating embedded domain isogeometric explicit dynamics method for collision simulation according to claim 1, characterized in that: The specific steps of step (2) are as follows: According to the minimum and maximum point coordinates of the axis-aligned bounding box of the three-dimensional geometric model, the background embedding domain geometry is defined using three-parameter non-uniform rational B-splines in the form of tensor product; the shape of the background embedding domain geometry is a regular hexahedron, and its size can accommodate the entire three-dimensional geometric model; the NURBS basis function of the background embedding domain geometry, i.e., the background embedding domain basis function, is refined and upgraded using the H refinement and P refinement methods to meet the accuracy requirements of dynamic finite element simulation, and a number of background embedding domain units are obtained.
4. The co-rotating embedded domain isogeometric explicit dynamics method for collision simulation according to claim 1, characterized in that: The specific steps of step (3) are as follows: There are three types of positional relationships between background embedded domain units and three-dimensional geometric models: intersecting units, in-domain units, and out-of-domain units. Ray tracing and collision detection can be used to determine the positional relationship: collision detection determines whether the background embedded domain units intersect with the three-dimensional geometric model. If not, the point position determination method is used to determine whether the maximum physical coordinate point and the minimum physical coordinate point of the background embedded domain unit are within the three-dimensional geometric model, thereby determining the relationship between the embedded domain unit and the inside and outside of the model. Among them, the unit that intersects with the three-dimensional geometric model is called an intersecting unit, the background embedded domain unit located inside the three-dimensional geometric model is called an in-domain unit, and the background embedded domain unit located outside the three-dimensional geometric model is called an out-of-domain unit.
5. The co-rotating embedded domain isogeometric explicit dynamics method for collision simulation according to claim 1, characterized in that: The specific steps of step (4) are as follows: Starting from the weak form of the momentum conservation equation and based on the Galerkin principle, the imaginary power equation is obtained: In the formula, δv is the virtual velocity, t is the surface force acting on the boundary Γ, b is the body force acting on the unit mass of the object, ρ is the density, is the second-order derivative of displacement with respect to time, σ is the Cauchy stress; δD is the virtual deformation rate tensor, defined as: The IGA motion equation is spatially discretized using the background embedding domain basis function, and the virtual power equation is written as: Where M is the mass matrix, f int and f ext are the internal force and external force of the control point; the specific expressions are as follows: M=∫ Ω ρN T NdΩ in int =∫ V B T ·σdV f ext =∫ V ρN·bdV+∫ Γ N·tdΓ Where N represents the background embedding domain basis function, and B is the shape function matrix; B=[B1...B A ...B n ] When the order of the background embedding domain basis function is quadratic or cubic, the units outside the domain do not participate in the calculation directly, and only the units inside the domain and the intersection units participate in the calculation; among them, the units inside the domain are normally Gaussian integrated, and the intersection units are integrated and solved by the combined Gaussian integration scheme; using the octree space partitioning method, an intersection unit is subdivided into eight integral sub-units; for the background embedding domain basis functions of p, q, and s B in three directions, first start from the 0th layer integration unit. If the integration unit is not clipped and is located inside the integration domain, then select (p+1)×(q+1)×( s+1) Gaussian points are integrated. If the integration unit is clipped by the boundary of the three-dimensional geometric model, it is evenly divided along the three directions of its parameter space. The subunit is divided into eight integration subunits. For each subunit completely inside the model, (p+1)×(q+1)×(s+1) Gaussian points are selected for integration. This process is recursively performed on all clipped integration units until the number of recursions reaches the target value to achieve relatively accurate Gaussian integration for the clipped units. When the order of the background embedding domain basis function is fourth order or above, the units outside the domain need to be multiplied by 10 when solving the mass matrix and the internal and external forces of the control points. -5 ≤α≤10 -10 In order to impose a penalty on the out-of-domain units and reduce the contribution of the out-of-domain units to the virtual power equation, the intersection units, the in-domain units and the background embedded domain basis functions are kept consistent when the order is low.
6. The co-rotating embedded domain isogeometric explicit dynamics method for collision simulation according to claim 1, characterized in that: The specific steps of step (5) are as follows: The co-rotation method requires corresponding stress and strain measurements in a co-rotation coordinate system; the co-rotation coordinate system is a local Cartesian coordinate system defined on the Gaussian points of the background embedded domain unit and rotates with the Gaussian points on the background embedded domain unit. The definition process is as follows: The covariant basis vectors g1 and g2 are expressed as: Where x = (x, y, z) T are the spatial coordinates of the material point measured in the global coordinate system, (ξ,η,ζ) are the parameter coordinates of the Gaussian points of the background embedding domain unit involved in the calculation; the covariant basis vectors g1 and g2 define a surface tangent to the layer ζ=d,d∈(-1,1), and the normal vector perpendicular to the surface is: Define a set of auxiliary vectors: Get the other two mutually perpendicular basis vectors of the co-rotational coordinate system: From this, we can define the rotation matrix R, whose matrix components R ij It can be expressed as: In the formula, the value ranges of the right subscripts i and j are 0, 1, and 2, corresponding to x, y, and z, respectively, that is, e0 = e x 、e1=e y 、e2=e z , e x Represents the unit vector in the x direction of the global coordinate system, e y Represents the unit vector in the y direction of the global coordinate system, e z Represents the unit vector in the z direction of the global coordinate system; The transformation relationship between velocity, displacement, Cauchy stress and deformation rate tensor in the co-rotation coordinate system and the global coordinate system is: The deformation rate tensor defined in the corotational coordinate system and the rotational Cauchy stress As an objective measure of strain rate and stress, it does not change with the rotation or translation of the background embedded domain unit; therefore, in the co-rotating coordinate system, the rotational Cauchy stress update rule is: Where Δt is the time step, is the rotating Cauchy stress rate, n is the number of time iterations; For linear elastic materials, the rotational Cauchy stress rate is It is expressed as: In the formula, is the strain matrix in the co-rotational coordinate system, is the local constitutive matrix defined in the corotational coordinate system; For the nonlinear plastic deformation of isotropic metal materials, the radial return algorithm based on the von Mises yield criterion is used to update the corresponding stress, which includes two steps: first, according to t n The converged solution at time t is used to calculate the trial stress. If the trial stress is within the linear elastic region of material deformation, it is considered as the true stress. Otherwise, it is projected to the nearest point on the yield surface of the material.
7. The co-rotating embedded domain isogeometric explicit dynamics method for collision simulation according to claim 1, characterized in that: The specific steps of step (6) are as follows: Objects A and B interact in the framework of the principles of continuum mechanics and must comply with the non-embedding condition, namely: Oh A ∩Oh B =0 Where Ω A represents the space occupied by object A, Ω B represents the space occupied by object B; When dealing with contact nonlinear problems, the explicit dynamics analysis method first ignores the contact between objects in each time step, and independently calculates the IGA motion equation of each object to obtain the trial configuration of each object after preliminary motion. If the trial configuration meets the non-embedded condition, it means that there is no contact between the objects, the motion between the objects does not affect each other, and the trial configuration is the real configuration. If the trial configuration does not satisfy the non-embedded condition, it indicates that the objects are in contact, and the contact force needs to intervene in the calculation and correct the force at the contact point to meet the contact condition; The contact force is calculated using the penalty function method: in each time step, first check whether the point on object A is embedded in the surface of object B. If there is no penetration, no processing is required; otherwise, a normal contact force similar to the "normal spring" effect is applied between the contact point of object A and the contact surface of object B to limit penetration; The normal contact force is: f c =p c ×l×N Where l is the penetration distance, N is the normal direction of the contact point at the contact surface; p c represents the normal stiffness, and its calculation formula is: In the formula, m c is the mass of the point; For the nonlinear contact problem of impacting the rigid wall, a separate treatment is performed accordingly: if the rigid wall is a finite rectangle with side lengths L and M, l and m are the unit vectors of the two adjacent sides of the rigid wall, r k n It is represented as the vector between the point on the contact detection surface and the origin; after each time step Δt of the global update of velocity and acceleration, all points k on the surface of each object need to be checked to ensure that all points on the surface of the object are located inside the rigid wall by checking whether two inequalities are satisfied: For the case of an infinite rigid wall, skip the above steps and only need to check the penetration condition to see whether the contact detection point penetrates the rigid wall surface, that is, you need to check the following formula: Where n is the normal vector of the contact point; If penetration occurs, corresponding processing is required to meet the non-embedded condition: the normal contact acceleration of the contact point where contact penetration has occurred is set to the negative value of the original normal acceleration, and the normal velocity of the contact point is set to 0 to correct the force at the contact point.
8. The co-rotating embedded domain isogeometric explicit dynamics method for collision simulation according to claim 1, characterized in that: The specific steps of step (7) are as follows: The central difference method is used to explicitly discretize the IGA motion equation in the time domain. First, the relevant parameters are initialized, namely the initial velocity, initial stress and total time, etc.; the mass matrix is solved. Dirichlet boundary conditions are applied; then the central difference method is used for time iteration; finally, the corresponding result configuration is output.