An isogeometric display dynamics calculation method for large deformation analysis of thin shell structures with complex non-conforming CAD models

Abstract

The present invention provides an isogeometric display dynamics calculation method for solving large deformation analysis of thin-shell structures in non-conforming complex CAD models, belonging to the field of display dynamics technology. First, a CAD design model is input and its NURBS parameters are obtained. Second, the CAD design model is subjected to isogeometric display dynamics analysis of large deformation thin-shell structures. Finally, the isogeometric display dynamics calculation results are output. The present invention effectively solves the processing difficulties of complex non-conforming CAD models in isogeometric analysis and display dynamics solutions for large deformation thin-shell structures.

Description

Technical Field

[0001] The invention belongs to the technical field of display dynamics and relates to an isogeometric display dynamics analysis method suitable for large deformation of thin shell structures of complex non-coordinated CAD models. Background Art

[0002] In the Finite Element Method (FEM), CAD (Auto Computer Aided Design) design models and CAE (Computer-Aided Engineering) engineering models are described using different mathematical languages. The transition from CAD design models to CAE engineering models requires time-consuming meshing, which is a discrete approximation and has certain errors. In addition, in order to obtain high-quality mesh units, the meshing process requires local refinement, which can also lead to parameter inconsistencies between units. If the CAD model can be used directly for computational solutions without the need for time-consuming additional operations such as meshing, the existing problems of the finite element analysis method will be eliminated.

[0003] Non-Uniform Rational B-Splines (NURBS) are a widely used modeling and interpolation tool in CAD design. They are used as a surface modeling method in various CAD design software. Surfaces modeled based on NURBS are called NURBS surfaces. Hughes et al. proposed isogeometric analysis (IGA) based on NURBS theory. IGA uses NURBS (non-uniform rational B-splines) curves and surfaces as the geometric basis for analyzing engineering models. This is the same representation used in CAD design models. This avoids the complex meshing process and allows CAD design models to be used directly as CAE engineering models. IGA has been widely applied to various areas of classical finite element analysis. Complex CAD design models often contain incompatible NURBS surfaces, making it particularly difficult to handle these inconsistent geometric information. This incompatibility exacerbates the solution difficulty, as the alignment and connection of these incompatible NURBS surfaces requires additional coupling algorithms to ensure model continuity and smoothness. This additional alignment and connection requires very complex coupling algorithms to implement, which is currently a difficult problem in isogeometric display dynamics analysis.

[0004] The Material Point Method (MPM) is a numerical computational method that combines a Lagrangian-Eulerian method (Lagrangian particles and an Euler background mesh) with a dual description. It effectively handles large deformation problems. Because the momentum equation is solved on an Euler background mesh, it can project incompatible NURBS surfaces onto a compatible background mesh and then solve them. This dual description of the Lagrangian and Euler methods allows the MPM framework to handle design models composed of discontinuous NURBS surfaces without requiring additional coupling algorithms. Standard MPMs often use linear interpolation functions for the background mesh, while Lagrangian particles are represented as point masses. This introduces numerical noise when material points cross the boundary between one background mesh and another. Numerous efforts have been made to address this issue. In 2004, Bardenhagen et al. proposed the Generalized Interpolating Material Point Method (GIMP). This method uses characteristic functions to describe Lagrangian particles, thereby improving the shape function and ensuring continuity at boundaries, eliminating the numerical noise inherent in the standard MPM. Based on GIMP, Sadeghirad et al. proposed the Convection Particle Domain Interpolation Material Point Method (CPDI1), which uses two vectors to represent the material point region and allows shear deformation of the material point region, solving the problem that GIMP has difficulty in handling parallelogram material point regions. Subsequently, Sadeghirad et al. proposed an improved Convection Particle Domain Interpolation Material Point Method (CPDI2) using linear Lagrangian basis functions. This method allows material points to deform arbitrarily, eliminating the problem of gaps in the material region due to deformation of material points in the CPDI1 method. Wang et al. proposed the GPDM method, which can perform large deformation analysis directly based on CAD data. The present invention studies the display dynamics problem of large deformation of thin shell structures of non-coordinated complex CAD models within the GPDM framework.

[0005] The development of shell theory and the numerical simulation of shell structures are the central topics of related research. The complexity and diversity of plate and shell theory also bring challenges to numerical methods. Shell theory essentially uses a mathematical two-dimensional configuration model to describe a thin-walled structure in three-dimensional space. Thin shell structures will produce significant geometric nonlinearity during large deformation, that is, the geometric shape of the structure after deformation is quite different from the initial state. Under large deformation, it is prone to buckling or instability, and its stress state and deformation mode will change rapidly. Buckling leads to sudden changes in stress distribution, which makes it difficult to ensure calculation accuracy and convergence, thereby increasing the complexity of problem solving. An effective solution algorithm for large deformation of thin shell structures is an important problem in the study of display dynamics.

[0006] The present invention proposes a new isogeometric display dynamics calculation method for solving large deformation analysis of thin shell structures in non-coordinated complex CAD models. Material points are generated from the NURBS surface information of the CAD design model. Within the framework of the GPDM material point method, a material point region is formed by mapping the material points to the control points of each NURBS unit. Each material point region represents a shell unit. Basis functions are calculated through a dual mapping from material points to NURBS unit control points and then to background grid nodes. The momentum equation is solved on the background grid, the background grid node positions are updated, and then mapped back to the material points, and the material point positions and velocities are updated. Internal forces are updated within the material point region based on the isogeometric shell unit theory, thereby realizing an isogeometric display dynamics analysis method suitable for large deformation of thin shell structures in complex non-coordinated CAD models. Summary of the Invention

[0007] In view of the defects of the existing technology in that it has poor adaptability to complex non-conforming CAD models and is complex in processing non-conforming NURBS coupling, resulting in inaccurate analysis results and low computational efficiency, the present invention provides an isogeometric display dynamic analysis method for solving large deformation analysis of thin shell structures of non-conforming complex CAD models.

[0008] In order to achieve the above object, the technical solution adopted by the present invention is:

[0009] A geometric display dynamics calculation method for solving large deformation analysis of thin shell structures of non-coordinated complex CAD models. The method is divided into three stages: model reading, algorithm solution and result processing. The specific steps are as follows: S1: Input the CAD design model and obtain the NURBS parameters of the CAD design model, including the number of NURBS surfaces, the degree of surfaces, the object parameter space grid topology information, the geometric positions of the NURBS surface control vertices and their weights, NURBS unit information, and the control point positions corresponding to each NURBS unit. S2: Perform geometric display dynamics solution analysis of large deformation of thin shell structures of the CAD design model. S3: Output the geometric display dynamics calculation results. Specifically including the following steps:

[0010] S1: Input the thin shell structure CAD design model and read the geometric dimensions of the CAD design model. The NURBS parameters of the CAD design model include NURBS surfaces, the number of NURBS surfaces, the degree of NURBS surfaces, the object parameter space grid topology information, NURBS surface control vertices and their weight factors, NURBS unit information, and NURBS control point basis functions.

[0011] S2: Perform iso-geometric display dynamic analysis of thin shell structure with large deformation on the CAD design model.

[0012] S21: Adding geometric display dynamics analysis calculation parameters such as large deformation of thin shell structure to the CAD design model described in S1, wherein the display dynamics analysis calculation parameters include: material information of each NURBS unit obtained in S1, initial boundary conditions of each NURBS unit obtained in S1, thickness of thin shell, etc., total time t for geometric display dynamics analysis solution end , etc. The geometric display dynamic analysis solves the initial time t0.

[0013] Furthermore, in S21, the material information of the NURBS unit includes detailed properties such as material density, mass, elastic modulus, Poisson's ratio, volume, bulk modulus, and shear modulus of each NURBS unit.

[0014] Furthermore, in S21, the initial boundary conditions of the NURBS unit include: the initial velocity, external force, and constraint conditions of each NURBS unit.

[0015] S22: Generate material points from the NURBS element information obtained in S1, with each NURBS element corresponding to a material point. The material information and boundary conditions for each NURBS element obtained in S21 are then assigned to the corresponding material point, along with the shell thickness parameter. This establishes a close connection between NURBS elements, material points, and related parameters, laying the foundation for subsequent calculations.

[0016] S23: Generate a background mesh according to the geometric dimensions of the CAD design model read in S1, and set the background mesh unit size according to the NURBS unit size.

[0017] S231: The geometric size of the background mesh is not less than the geometric size of the CAD design model described in S1, and the unit size of the background mesh is not less than the average size of the NURBS units described in S1. Setting this size relationship ensures that the background mesh completely covers the CAD design model and better matches the NURBS units in terms of accuracy, thereby improving the accuracy and reliability of the entire solution analysis.

[0018] S24: Input the current time t, and calculate the time step Δt based on the geometry display dynamic analysis based on the background grid size and the velocity calculation of the material point at the current time.

[0019] S241: The time step for the isogeometric display dynamics analysis calculation uses an adaptive time step, that is, the time step is not fixed, but is adjusted according to the velocity of the material point. This allows for more precise adjustment of the calculation time step, thereby improving calculation efficiency and accuracy. The calculation method is as follows:

[0020]

[0021] Where Δt is the time step of the isogeometric dynamics analysis; (h x ,h y ,h z ) is the cell size of the background grid; α is the time step coefficient, which ranges from 0 to 1; min indicates the minimum value; c x is the velocity coefficient of the material point in the x direction, c y is the velocity coefficient of the material point in the y direction, c z is the velocity coefficient of the material point in the z direction, c i The formula is as follows, where the subscript i = x, y, z:

[0022]

[0023] Where p represents the material point, v t xp is the x component of the velocity of the material point p at time t, v t yp is the y component of the velocity of the material point p at time t, max means the maximum value, v t zp is the z component of the velocity of the material point p, c dil is a coefficient.

[0024] c dil The formula is as follows:

[0025]

[0026] where ρ, K, and G are the particle density, bulk modulus, and shear modulus, respectively.

[0027] S25: Calculate the interpolation basis function between the material points generated by S22 and the background grid generated by S23, and the interpolation basis function gradient.

[0028] S251: The present invention performs calculations within the framework of the GPDM material point method. The mapping basis functions involved in this framework consist of two parts: one is the mapping function from the material point to the NURBS unit control point, and the other is the mapping function from the NURBS unit control point to the background grid. In the GPDM material point method, multiplying the NURBS basis function with the node interpolation function yields the weight mapping function and its gradient for the material point, as follows:

[0029] N I (x)=∑ c R c (ξ,η)N I (x c )

[0030]

[0031] Among them, x c is the position of the NURBS unit control point; Rc(ξ,η) is the NURBS basis function corresponding to the control point c; N I (xc) is the background grid interpolation function, N I (x c ) is generally calculated using the hat function. The hat function is a piecewise function defined on an interval. Intuitively, its image shape is similar to a "hat". The middle part is a linear rise and then a linear fall process, and the two ends are zero, hence the name hat function. I (x) is the weight function from the material point to the background grid;

[0032] is the gradient of the weight function from the material point to the background grid with respect to the spatial coordinate x; is the gradient of the spatial coordinate x; (ξ,η) is the parameter coordinate of the control point.

[0033] Integrate the weight mapping function and its gradient to obtain the interpolation basis function Φ between the material point and the background grid Ip and interpolation basis function gradients

[0034]

[0035] Where, is the weight of the interpolation basis function of the background mesh node corresponding to the NURBS control point, is the weight of the interpolation basis function derivative of the NURBS control point to the corresponding background mesh node, Ω represents the NURBS unit area, and (ξ,η) represents the local coordinates of the current NURBS unit control point. V p is the volume of the material point; Ω p is the area of ​​the material point domain.

[0036] Furthermore, in the framework of the material point method of GPDM, the material point domain is constructed by mapping each material point to the NURBS unit control point with the help of the mapping function, and its area is equivalent to the NURBS unit area.

[0037] S26: Map the physical information contained in each material point generated in S22 to the background grid generated in S23, and obtain the mass, momentum, internal force, and external force of each background grid node.

[0038] S261: The physical information mapped in S24 includes mass, velocity, momentum, node internal force, and node external force. The physical information can be expressed as:

[0039] Mapping 1: in, and m p are the position and mass of the material point P at time t, φ I is the mapping basis function. Represents the quality of background grid node I at time t.

[0040] Mapping 2:

[0041] Among them, v p t is the velocity and mass of the material point P at time t, Represents the velocity of background grid node I at time t.

[0042] Through these two mapping relationships, the velocity and mass information of the material point are first mapped to the background grid node, and then the momentum of the background grid node can be calculated, thereby establishing the momentum information corresponding to the material point on the background grid. The momentum information can be expressed as:

[0043]

[0044] Among them, P I t represents the momentum of background grid node I at time t.

[0045] Mapping 3: Among them, b(x p ) is the external force on the material point p described in S22.

[0046] Mapping 4: in, and are the volume and Cauchy stress of the material point P, is the gradient of the mapping basis function.

[0047] The Cauchy stresses are calculated based on the isogeometric shell stress update theory. Within the framework of the GPDM material point method, isogeometric shell elements have specific properties. Their area is equivalent to the area of ​​the material point domain constructed by each material point using the mapping function, which is equivalent to the area of ​​the NURBS element. Their thickness is equivalent to the shell thickness assigned to the material point by S22.

[0048] S27: Solve the momentum equation on the background grid and update the background grid momentum. The momentum equation can be expressed as:

[0049]

[0050] Where Δt is the S232 time step, is the momentum of the updated background grid node I, is the resultant force acting on background grid node I at time t, which can be expressed as:

[0051]

[0052] S28: Map the updated momentum of the background grid obtained in S26 back to the material point, and update the velocity, position, and stress of the material point. Reset the background grid and set all the background grid information obtained by S25 to 0.

[0053] This step maps the updated momentum information on the background grid back to the material point, thereby updating the material point's velocity, position, stress, and other parameters, achieving a two-way flow of information between the material point and the background grid. Resetting the background grid information to 0 allows for recalculation in the next iteration to ensure accuracy and continuity.

[0054] S281: The updating of the material point velocity can be expressed as:

[0055]

[0056] in, represents the velocity of the material point P after the update; represents the velocity of the material point P at time t; Δt represents the time step calculated by S241; represents the position of the material point P at time t; represents the resultant force acting on background grid node I at time t; represents the node quality of background grid node I at time t.

[0057] The position update of the material point can be expressed as:

[0058]

[0059] in, represents the position of the material point P after the update; represents the position of the material point P at time t; m represents the mass of the background grid node I after the update; v represents the velocity of the background grid node I after the update; Represents the momentum of background grid node I after update.

[0060] The stress update of the material point can be expressed as:

[0061]

[0062] Among them, Δσ p is the Cauchy stress increment at the material point

[0063] The Cauchy stress increment is calculated based on the isogeometric shell stress update theory.

[0064] S29: Update the current time t and determine whether the total time t of the dynamic analysis set in S21 is reached end If the current time t is greater than or equal to t end , the output isogeometric display dynamic analysis results.

[0065] Furthermore, the updating of the current time t described in S29 can be expressed by the formula: t=t+Δt

[0066] By continuously updating the current time t and comparing it with the set total time, the loop and termination of the entire solution analysis process are controlled to ensure that the entire calculation process can proceed according to the predetermined plan. When the updated current time t is greater than or equal to the total time t end , the dynamic analysis ends and the output geometry shows the dynamic analysis results.

[0067] Furthermore, the iso-geometric display of the dynamic analysis results described in S29 includes: solving the geometric changes of the CAD design model.

[0068] S3: Use VTK (Visualization Toolkit) format to perform output operations and convert the solution results into text files. Paraview can then read the output results to achieve result visualization.

[0069] The beneficial effects of the present invention are:

[0070] The present invention uses NURBS curves and surfaces as the geometric basis for analyzing engineering models, which is the same as the representation used in CAD design models. The complex meshing process is avoided during the analysis process, reducing errors caused by meshing and a series of problems such as parameter inconsistency between units due to local refinement. At the same time, the present invention addresses the problems that the existing technology has poor adaptability to complex non-conforming CAD models and is extremely complex in processing non-conforming NURBS surface coupling, resulting in inaccurate analysis results and low computational efficiency. Through a specific technical solution, the present invention uses the GPDM material point method framework with dual description of Lagrangian and Euler methods to analyze design models composed of discontinuous NURBS surfaces, establishes a close connection between material points, NURBS unit control points, and background grid node related parameters, eliminating the need to add additional coupling algorithms to ensure the continuity and smoothness of the model as in the past, thereby reducing the processing difficulty. At the same time, the isogeometric display dynamics calculation method proposed in combination with the isogeometric shell stress update theory for solving large deformation analysis of thin shell structures in non-conforming complex CAD models is specifically designed for such models, effectively overcoming the shortcomings of the existing technology in this regard and being able to better handle the analysis and calculation of such complex models. BRIEF DESCRIPTION OF THE DRAWINGS

[0071] Figure 1 is a flow chart of the method of the present invention;

[0072] Figure 2 CAD model of the car's B-pillar shell and rigid column;

[0073] Figure 3 Schematic diagram of the mapping relationship between material point, NURBS unit control point and background grid node;

[0074] Figure 4 The geometric deformation diagram of the B-pillar at each moment. DETAILED DESCRIPTION

[0075] In order to make the objectives, technical solutions and advantages of the present invention more clearly understood, the present invention is further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only intended to illustrate the present invention and are not intended to limit the present invention. In addition, the technical features involved in the various embodiments of the present invention described below may be combined with each other as long as they do not conflict with each other.

[0076] Example: Numerical example of automobile B-pillar shell colliding with rigid column;

[0077] In the field of automobile collision analysis, the B-pillar is a key load-bearing component in the vehicle body structure, and its deformation during a collision is crucial for vehicle safety assessment. During a collision, the vehicle's B-shell collides first, while the rigid pillar simulates an obstacle in an actual collision scenario. Analyzing the collision between the two effectively tests the present invention's ability to handle such complex structural collision deformation analysis. Therefore, these two models were selected as examples. Traditional isogeometric display dynamics algorithms are difficult to calculate because the CAD design model of the vehicle's B-pillar shell is often constructed from multiple non-coordinated NURBS surfaces. This example further illustrates the superiority of the present invention's technology. The vehicle's B-pillar shell used in this example is constructed from 65 non-coordinated NURBS surfaces, while the rigid pillar is constructed from two non-coordinated NURBS surfaces.

[0078] S1: Input Model 1: Car B-pillar and Model 2: Rigid column, the CAD design models of two thin shell structural parts are as follows Figure 2 As shown, the geometric dimensions of the two CAD design models and the NURBS parameters of the two CAD design models are read: including NURBS surfaces, the number of NURBS surfaces, the degree of NURBS surfaces, the object parameter space grid topology information, NURBS surface control vertices and their weight factors, NURBS unit information, and NURBS control point basis functions.

[0079] S2: Isogeometric display dynamics analysis of the large deformation of the thin shell structure of the rigid column in the collision of the automobile B-pillar.

[0080] S21: For the two thin shell CAD design models described in S1, add geometric display dynamic analysis calculation parameters such as large deformation of the thin shell structure.

[0081] Add rigid material to each NURBS unit in NURBS surface 1 and NURBS surface 2 on model 1, add plastic material to all NURBS units on the remaining NURBS surfaces and apply an initial velocity of 30m / s to all NURBS units. Figure 2 shown.

[0082] Add rigid materials to the two NURBS surfaces of model 2 and apply fixed constraints to the two NURBS surfaces of model 2. The shell thickness is set to 1mm; the total time for the dynamic analysis solution is t end =0.006s; the isogeometric display dynamic analysis solves the initial time t0=0s.

[0083] The basic material parameters of the plastic material are: (1) density: 7.8e -3 kg / m 3 (2) Poisson's ratio: 0.3 (3) Elastic modulus: 2e 5 Pa (4) shear modulus: 2.5e 5 Pa (4) yield strength: 2.2e 5 pa;

[0084] The basic material parameters of the rigid material are: (1) density: 7.8e 3 kg / m 3 (2) Poisson's ratio: 0.3 (3) Elastic modulus: 2e 5 ;

[0085] S22: Generate material points based on the NURBS element information of the two design models obtained in S1. Then, assign the material information of each NURBS element and the boundary conditions to the material points using S21, and also assign the shell thickness parameter to the material points.

[0086] S23: Generate a background grid according to the geometric dimensions of the two CAD design models read in S1, and calculate the time step Δt for isogeometric display dynamics analysis according to the background grid dimensions.

[0087] The geometric dimensions of the two design models are obtained by S1 reading to obtain the maximum and minimum values ​​of their x, y, and z coordinates in the Cartesian coordinate system, and a background grid is generated. The minimum size of the background grid is:

[0088] X min =1200mm; Y min =-1000mm; Z min =-200mm;

[0089] X max =1900mm; Y max =600mm; Z max =1400mm;

[0090] Among them, X min ,Y min ,Z min Indicates the coordinates of the minimum X, Y, and Z points of the background grid in the Cartesian coordinate system. max ,Y max ,Z max Indicates the coordinates of the maximum X, Y, and Z points of the background grid in the Cartesian coordinate system. Since the background grid is a standard cuboid, the background grid outline can be constructed based on these two points.

[0091] The average size of the NURBS unit obtained by S1 reading is 4.9mm*4.8mm. In order to match the average size of the NURBS unit, the background grid size is set to (h x ,h y ,h z )=(5mm, 5mm, 5mm).

[0092] where h x ,h y ,h z Represents the size of the background grid in the x, y, and z directions respectively.

[0093] S24: Enter the current time t and calculate the time step Δt based on the background grid size and the velocity calculation of the material point at the current time. The time step is calculated as described above.

[0094] S25: Calculate the mapping basis functions and gradients between the material points generated in S22 and the background mesh generated in S23. The NURBS degree information obtained from S1 indicates that the NURBS degree of both CAD design models is 2. Therefore, each NURBS unit is composed of 9 control points, which are interconnected to form a NURBS control network. The present invention is implemented within the framework of the GPDM material point method. Therefore, the mapping basis function between the material point and the generated background mesh can be expressed as:

[0095]

[0096]

[0097] Among them, xc is the position of the NURBS unit control point, the quadratic NURBS basis function corresponding to the Rc(ξ,η) control point, N I (x) is the weight function from the material point to the background grid; is the gradient of the weight function from the material point to the background grid with respect to the spatial coordinate x; is the gradient of the spatial coordinate x...; (ξ,η) are the parameter coordinates of the control point. (ξ,η) are the local coordinates of the NURBS unit control point, Ω represents the NURBS unit area, V p is the volume of the material point; Ω p is the area of ​​the material point domain.

[0098] Figure 3 It can help to better understand the mapping relationship between material points and background grids.

[0099] S26: Map the physical information contained in each material point of the two models generated in S22 to the background grid generated in S23, and obtain the mass, momentum, internal force, and external force of each background grid node.

[0100] Mapping 1 Material point mass mapping: in, and m p are the position and mass of the material point P at time t, φ I is the mapping basis function obtained by S25. Represents the quality of background grid node I at time t.

[0101] Mapping 2:

[0102] Among them, v p t is the velocity and mass of the material point P at time t, Represents the velocity of background grid node I at time t.

[0103] Through these two mapping relationships, the velocity and mass information of the material point are first mapped to the background grid node, and then the momentum of the background grid node can be calculated, thereby establishing the momentum information corresponding to the material point on the background grid. The momentum information can be expressed as:

[0104]

[0105] Among them, P I t represents the momentum of background grid node I at time t. represents the momentum of background grid node I at time t

[0106] Mapping 3: in, represents the external force on the background grid node I at time t, b(x p ) is the external force on the material point p described in S22.

[0107] Mapping 4: in, represents the internal force on the background grid node I at time t, and are the volume and Cauchy stress of the material point P at time t, is the gradient of the mapping basis function.

[0108] The Cauchy stresses are calculated using the isogeometric shell stress update theory.

[0109] S27: Solve the momentum equation on the background grid and update the background grid momentum.

[0110] S271: The momentum equation can be expressed as:

[0111]

[0112] Wherein, Δt is the time step described in S232, is the momentum of the updated background grid node I, is the resultant force acting on background grid node I at time t, which can be expressed as:

[0113]

[0114] S28: Map the updated momentum of the background grid obtained in S26 back to the material point, and update the velocity, position, and stress of the material point. Reset the background grid and set all the background grid information obtained by S25 to 0. S281: The update of the material point velocity can be expressed as:

[0115]

[0116] in, Represents the velocity of the material point P after the update.

[0117] The position update of the material point can be expressed as:

[0118]

[0119] in Indicates the updated material point position.

[0120] The stress update of the material point can be expressed as:

[0121]

[0122] in, is the stress of the updated material point, Δσ p is the Cauchy stress increment at the material point. The Cauchy stress increment is calculated using the isogeometric shell stress update theory. Resetting the background mesh information to 0 allows for a fresh calculation in the next loop to ensure accuracy and continuity.

[0123] S29: Update the current time t and determine whether the total time t of the dynamic analysis set in S21 is reached end If the current time t is greater than or equal to t end , then the output isogeometric display dynamic analysis solution results.

[0124] By continuously updating the current moment and comparing it with the set total time, the loop and termination of the entire solution analysis process are controlled to ensure that the entire calculation process can be carried out according to the predetermined plan. During the collision of the car B-pillar shell, as time goes by, by updating the current moment t and judging whether the total time has been reached, the entire dynamic process of the collision process from the beginning to the end can be fully simulated. When the current moment t is greater than or equal to the total time t after the update end , it indicates that the dynamic analysis is finished and the output shows the dynamic analysis results.

[0125] Furthermore, the isogeometric display dynamic analysis results described in S29 include: geometric changes of the B-pillar shell during the collision with the rigid column.

[0126] S3: The geometric display dynamic analysis results of the collision between the car B-pillar shell and the rigid column are output using VTK (V format, and the solution results are converted into text files. The output results are read by Paraview to realize the visualization of the results. It shows the geometric changes of the car B-pillar at 0.003s and 0.006s during the collision process obtained by the dynamic analysis, such as Figure 4 As shown. Figure 4 It can be seen that during the collision process, the geometric deformation of the B-pillar where the collision occurs gradually increases, which is in line with the expected effect.

[0127] The above-described embodiments merely express the implementation methods of the present invention, but should not be understood as limiting the scope of the patent of the present invention. It should be pointed out that for those skilled in the art, several variations and improvements can be made without departing from the concept of the present invention, and these all fall within the scope of protection of the present invention.

Claims

1. An isogeometric display dynamics calculation method for solving large deformation analysis of thin shell structures of non-coordinated complex CAD models, characterized by: The method described includes three stages: model reading, algorithm solving, and result processing, as follows: S1: Input the thin shell structure CAD design model and obtain the NURBS parameters of the CAD design model. The NURBS parameters of the CAD design model include NURBS surface, number of NURBS surfaces, degree of NURBS surface, object parameter space grid topology information, geometric positions of NURBS surface control vertices and their weight factors, NURBS unit information, control point positions corresponding to each NURBS unit, and NURBS control point basis functions; S2: Perform iso-geometric display dynamics analysis of thin shell structure with large deformation on the CAD design model, as follows: S21: Adding geometric display dynamic analysis calculation parameters such as large deformation of thin shell structure to the CAD design model described in S1, wherein the display dynamic analysis calculation parameters include: material information of each NURBS unit obtained in S1, initial boundary conditions of each NURBS unit obtained in S1, thickness of thin shell, total time t for geometric display dynamic analysis solution, etc. end , etc. Geometry shows that the dynamic analysis solves the initial time t0; S22: Generate material points based on the NURBS unit information obtained in S1, with each NURBS unit corresponding to a material point one-to-one. The material information and boundary conditions of each NURBS unit obtained in S21 are assigned to the corresponding material point, and the shell thickness parameter is also assigned to each material point. Through this operation, the relationship between NURBS units, material points, and related parameters is established. S23: Generate a background mesh according to the geometric dimensions of the CAD design model read in S1, and set the background mesh unit size according to the NURBS unit size; S24: Input the current time t, and calculate the time step Δt based on the geometry display dynamic analysis according to the background grid size and the velocity calculation of the material point at the current time; S25: Calculate the interpolation basis function between the material points generated by S22 and the background grid generated by S23, and the interpolation basis function gradient; S26: Map the physical information contained in each material point generated by S22 to the background grid generated by S23; obtain the mass, momentum, internal force and external force of each background grid node; S27: Solve the momentum equation on the background grid and update the background grid momentum; S28: Map the updated momentum of the background grid obtained in S26 back to the material point, and update the velocity, position, and stress of the material point; reset the background grid and set all the background grid information obtained by S25 to 0; S29: Update the current time t and determine whether the total time t of the dynamic analysis set in S21 is reached end If the time is not reached, return to S23; if the current time t is greater than or equal to t end , then the output isogeometric display dynamic analysis results; S3: Output isogeometric display of dynamic calculation results.

2. The isogeometric display dynamics calculation method for solving large deformation analysis of thin shell structures of non-coordinated complex CAD models according to claim 1 is characterized in that: In S21, the material information of the NURBS unit includes: the material density, mass, elastic modulus, Poisson's ratio, volume, bulk modulus, and shear modulus of each NURBS unit; in S21, the initial boundary conditions of the NURBS unit include: the initial velocity, external force, and constraint conditions of each NURBS unit.

3. The isogeometric display dynamics calculation method for solving large deformation analysis of thin shell structures of non-coordinated complex CAD models according to claim 2 is characterized in that: In said S23, the geometric size of the background grid is not less than the geometric size of the CAD design model described in S1, and the unit size of the background grid is not less than the average size of the NURBS unit described in S1.

4. According to the isogeometric display dynamics calculation method for solving large deformation analysis of thin shell structures of non-coordinated complex CAD models according to claim 3, the S4 is specifically as follows: The calculation time step of the isogeometric display dynamics analysis adopts the adaptive time step, which is calculated as follows: in, Δt is the time step for isogeometric dynamics analysis; (h x ,h y ,h z ) is the cell size of the background grid; α is the time step coefficient, which ranges from 0 to 1; min indicates the minimum value; c x is the velocity coefficient of the material point in the x direction, c y is the velocity coefficient in the y direction of the material point, c z is the velocity coefficient of the material point in the z direction, c i The formula is as follows, where the subscript i = x, y, z: Where p represents the material point, v t xp is the x component of the velocity of the material point p at time t, v t yp is the y component of the velocity of the material point p at time t, max means the maximum value, v t zp is the z component of the velocity of the material point p, c dil is a coefficient; c dil The formula is as follows: where ρ, K, and G are the particle density, bulk modulus, and shear modulus, respectively.

5. According to the isogeometric display dynamics calculation method for solving large deformation analysis of thin shell structures of non-coordinated complex CAD models according to claim 4, the S5 is specifically as follows: The calculation is performed under the framework of the GPDM material point method. The mapping basis function involved consists of two parts: one is the mapping function from the material point to the NURBS unit control point, and the other is the mapping function from the NURBS unit control point to the background grid. In the GPDM material point method, the NURBS basis function is multiplied by the node interpolation function to obtain the weight mapping function of the material point and its gradient. The formula is as follows: N I (x)=Σ c R c (ξ,η)N I (x c ) ▽ x N I (x)=Σ c ▽ x R c (ξ,η)N I (x c ) in, x c is the position of the NURBS unit control point; R c (ξ,η) NURBS basis function corresponding to control point c; N I (x c ) is the background grid interpolation function, N I (x c ) is calculated using the hat function; N I (x) is the weight function from the material point to the background grid; ▽xN I (x) is the gradient of the weight function from the material point to the background grid with respect to the spatial coordinate x; x is the gradient of the spatial coordinate x; (ξ,η) is the parameter coordinate of the control point; Integrate the weight mapping function and its gradient to obtain the interpolation basis function Φ between the material point and the background grid Ip and the interpolation basis function gradient ▽Φ Ip : Where, is the weight of the interpolation basis function of the background mesh node corresponding to the NURBS control point, is the weight of the derivative of the interpolation basis function of the NURBS control point to the corresponding background grid node, Ω represents the area of ​​the NURBS unit, (ξ,η) represents the local coordinates of the current NURBS unit control point; V p is the volume of the material point; Ω p is the area of ​​the material point domain; In the framework of the material point method of GPDM, the material point domain is constructed by mapping each material point to the NURBS unit control point with the help of the mapping function, and its area is equal to the NURBS unit area.

6. According to the isogeometric display dynamics calculation method for solving large deformation analysis of thin shell structures of non-harmonious complex CAD models according to claim 5, the specific S6 is as follows: The mapped physical information of S24 includes mass, velocity, momentum, node internal force, and node external force; the physical information is represented as: Mapping 1: in, and m p are the position and mass of the material point P at time t, φ I is the mapping basis function; represents the quality of background grid node I at time t; Mapping 2: Among them, v p t is the velocity and mass of the material point P at time t, represents the velocity of background grid node I at time t; Through these two mapping relationships, the velocity and mass information of the material point are first mapped to the background grid node, and then the momentum of the background grid node can be calculated, thereby establishing the momentum information corresponding to the material point on the background grid. The momentum information can be expressed as: Among them, P I t represents the momentum of background grid node I at time t; Mapping 3: Among them, b(x p ) is the external force at the material point p described in S22; Mapping 4: in, and are the volume and Cauchy stress of the material point P, respectively, I is the mapping basis function gradient; The Cauchy stress is calculated based on the isogeometric shell stress update theory; under the framework of the GPDM material point method, the area of ​​the isogeometric shell unit is equivalent to the area of ​​the material point domain constructed by each material point with the help of the mapping function, and this situation is equivalent to the NURBS unit area, and its thickness is equivalent to the thin shell thickness assigned to the material point by S22.

7. The isogeometric display dynamics calculation method for solving large deformation analysis of thin shell structures of non-harmonious complex CAD models according to claim 6, wherein said S7 is specifically as follows: The momentum equation is expressed as: in, Δt is the S232 time step, is the momentum of the updated background grid node I, is the resultant force acting on background grid node I at time t, which can be expressed as: The S8 is specifically as follows: The update of material point velocity can be expressed as: in, represents the velocity of the material point P after the update; represents the velocity of the material point P at time t; Δt represents the time step calculated by S241; represents the position of the material point P at time t; represents the resultant force acting on background grid node I at time t; represents the node quality of background grid node I at time t; The position update of the material point can be expressed as: in, represents the position of the material point P after the update; represents the position of the material point P at time t; m represents the mass of the background grid node I after the update; v represents the velocity of the background grid node I after the update; represents the momentum of background grid node I after update; The stress update of material points can be expressed as: Among them, Δσ p is the Cauchy stress increment at the material point The Cauchy stress increment is calculated based on the isogeometric shell stress update theory.

8. According to the isogeometric display dynamics calculation method for solving large deformation analysis of thin shell structures of non-harmonious complex CAD models according to claim 7, the S9 is specifically as follows: The updating current time t can be expressed by the formula: t = t + Δt; by continuously updating the current time t and comparing it with the set total time, the cycle and termination of the entire solution analysis process are controlled; when the updated current time t is greater than or equal to the total time t end , the dynamic analysis ends and the output geometry shows the dynamic analysis results.

Images