Isogeometric display dynamics calculation method for solving large deformation analysis of thin shell structure of non-coordinated complex CAD (computer-aided design) model

Abstract

The invention provides an isogeometric display dynamics calculation method for solving large deformation analysis of a thin shell structure of a non-coordinated complex CAD model, and belongs to the technical field of display dynamics. The method comprises the steps that firstly, a CAD design model is input, and NURBS parameters of the CAD design model are obtained; secondly, performing isogeometric display dynamics solution analysis on the large deformation of the thin shell structure on the CAD design model; and finally, outputting an isogeometric display dynamics calculation result. According to the method, the processing problem of the complex non-coordinated CAD model in geometric analysis display dynamics solution of large deformation of the thin shell structure and the like can be effectively solved.

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 a thin shell structure of a complex non-coordinated CAD model. 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 idea with 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 modeling method for surface modeling by various CAD design software. Surfaces modeled based on NURBS are called NURBS surfaces. Hughes et al. proposed isogeometric analysis (IGA) based on non-uniform rational B-spline theory. IGA uses NURBS (non-uniform rational B-spline) curves and surfaces as the geometric basis for analyzing engineering models. This is the same as the representation used in CAD design models. The complex meshing process is avoided during the analysis process, and the CAD design model can be used directly as a CAE engineering model. It has been extended to various fields of classical finite element analysis. Complex CAD design models usually contain incompatible NURBS surfaces, and it is particularly difficult to process these inconsistent geometric information. This incompatibility will increase the difficulty of solving, because additional coupling algorithms are required in the process of aligning and connecting incompatible NURBS surfaces to ensure the continuity and smoothness of the model. This additional alignment and connection requires very complex coupling algorithms to implement, which is currently a challenge in isogeometric display dynamics analysis.

[0004] The Material Point Method (MPM) is a numerical calculation method with dual description of Lagrangian-Euler method (Lagrangian particles and Euler background grid), which can effectively deal with large deformation problems. Because its momentum equation is solved on the Euler background grid, the originally inconsistent NURBS surface can be projected onto the coordinated background grid and then solved. The material point method framework with dual description of Lagrangian and Euler methods can handle the design model composed of discontinuous NURBS surfaces in the analysis without adding additional coupling algorithms. The background grid of the standard material point method mostly uses linear interpolation functions, and the Lagrangian particles are represented by particles. In this way, numerical noise will be generated when the material point passes through the boundary from the current background grid to another background grid. In order to solve this problem, scholars have made a lot of efforts. In 2004, Bardenhagen et al. proposed the generalized interpolation material point method (GIMP), which uses characteristic functions to describe Lagrangian particles, thereby improving the shape function, making the shape function continuous at the boundary, and eliminating the numerical noise of the standard material point method. Based on GIMP, Sadeghirad et al. proposed the Convective 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 Convective 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 a three-dimensional space. Thin shell structures will produce significant geometric nonlinearity during large deformations, that is, the geometric shape of the structure after deformation is quite different from the initial state, and it is prone to buckling or instability under large deformations, 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 solving the problem. An effective solution algorithm for large deformations 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 of non-coordinated complex CAD models. Material points are generated through the NURBS surface information of the CAD design model, and a material point area is formed by mapping the material points with the control points of each NURBS unit under the framework of the GPDM material point method. Each material point area represents a shell unit. The basis function is calculated by double mapping from material points to NURBS unit control points and then to background grid nodes, and the momentum equation is solved on the background grid, the background grid node position is updated, and then mapped back to the material point and the material point position and velocity are updated. The internal force is updated in the material point area based on the isogeometric shell unit theory, thereby realizing an isogeometric display dynamics analysis method suitable for large deformation of thin shell structures of complex non-coordinated CAD models. Summary of the invention

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

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

[0009] A method for isogeometric display dynamics calculation 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 to 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 isogeometric display dynamics solution analysis of large deformation of thin-shell structures of the CAD design model. S3: Output the isogeometric display dynamics calculation results. Specifically including the following steps:

[0010] S1: Input the thin shell structure CAD design model, read the geometric dimensions of the CAD design model, and 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, NURBS surface control vertices and their weight factors, NURBS unit information, and NURBS control point basis functions.

[0011] S2: Perform iso-geometric display dynamic solution analysis of thin shell structure large deformation on the thin shell structure 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 by S1, initial boundary conditions of each NURBS unit obtained by S1, thickness of thin shell, and total time t for solving geometric display dynamics analysis. 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 of the material density, mass, elastic modulus, Poisson's ratio, volume, bulk modulus, shear modulus, etc. of each NURBS unit.

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

[0015] S22: Generate material points from the NURBS unit information obtained in S1, and each NURBS unit corresponds to a material point one by one; assign the material information and boundary conditions of each NURBS unit obtained in S21 to the corresponding material point, and also assign the shell thickness parameter to each material point. Through such operations, a close connection between NURBS units, material points and related parameters is established, laying the foundation for subsequent calculations.

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

[0017] S231: 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. Setting such a size relationship can ensure that the background grid can completely cover the CAD design model and can better match the NURBS unit 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 such as the background grid size and the velocity calculation of the material point at the current time.

[0019] S241: The calculation time step of the isogeometric display dynamics analysis adopts an adaptive time step, that is, the time step is not fixed, but adjusted according to the speed of the material point, which can more accurately adjust the calculation time step, thereby improving the calculation efficiency and accuracy. The calculation method is as follows:

[0020]

[0021] Where Δt is the time step of 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 means the minimum value; c x is the velocity coefficient in the x direction of the material point, 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 a 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 taking 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 and the interpolation basis function gradient between the material points generated by S22 and the background grid generated by S23.

[0028] S251: The present invention is a calculation performed in the framework of the GPDM material point method. The mapping basis function involved in this framework 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 weight mapping function of the material point and its gradient can be obtained by multiplying the NURBS basis function with the node interpolation function. The formula is 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 a 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 material point to 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 the interpolation basis function gradient

[0034]

[0035] In the formula, 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 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 equal to the NURBS unit area.

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

[0038] S261: The mapped physical information described in S24 includes mass, velocity, momentum, node internal force, and node external force. The physical information described 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] Map 3: Among them, b(x p ) is the external force on the material point p described in S22.

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

[0047] The Cauchy stress is calculated based on the isogeometric shell stress update theory. Under the framework of the material point method of GPDM, the isogeometric shell element has specific properties. Its area is equivalent to the area of ​​the material point domain constructed by each material point with the help of the mapping function, which is equivalent to the area of ​​the NURBS unit, and its thickness is equivalent to the thin 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 on the background grid node I at time t, which can be expressed by the formula:

[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 mapping to 0.

[0053] This step is to map the momentum information updated on the background grid back to the material point, thereby updating the velocity, position, stress and other parameters of the material point, and realizing the two-way flow of information between the material point and the background grid. Resetting the background grid information to 0 is to recalculate in the next cycle to ensure the accuracy and continuity of the calculation.

[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 material point P at time t; represents the resultant force on the 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 updated position of the material point P; 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 dynamics analysis set in S21 is reached end If the current time t is greater than or equal to t end , the output isogeometric display of dynamic analysis results.

[0065] Furthermore, the updating 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 be carried out according to the predetermined plan. When the current time t is greater than or equal to the total time t after the update end , the dynamic analysis ends and the output isogeometric display 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 output the solution results into text files, which can then be read by Paraview to visualize the results.

[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 method used in CAD design models. The complex meshing process is avoided during the analysis process, and the errors caused by meshing and a series of problems such as incoordination of parameters between units due to local refinement processing are reduced. At the same time, the present invention is aimed at the poor adaptability of the prior art to complex non-coordinated CAD models, and the extremely complex processing of non-coordinated NURBS surface coupling, resulting in inaccurate analysis results and low calculation efficiency. Through a specific technical solution, the present invention analyzes the design model composed of discontinuous NURBS surfaces. The present invention adopts the GPDM material point method framework of the dual description of Lagrangian and Euler methods to establish a close connection between material points, NURBS unit control points, and background grid node related parameters. There is no need to add additional coupling algorithms as in the past to ensure the continuity and smoothness of the model, which reduces the processing difficulty; at the same time, the isogeometric display dynamics calculation method for solving the large deformation analysis of thin shell structures of non-coordinated complex CAD models proposed by combining the isogeometric shell stress update theory is specially designed for such models, effectively overcoming the defects of the prior art in this regard, and can 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 B-pillar shell and rigid column;

[0073] Figure 3 It is a 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 time. DETAILED DESCRIPTION

[0075] In order to make the purpose, technical solutions and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain 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 can 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 automobile B-pillar is a key load-bearing component in the body structure, and its deformation during the collision process is crucial to the vehicle safety assessment; during the collision process, the automobile B shell collides first, while the rigid column can simulate the obstacles in the actual collision scene. The analysis of the collision between the two can effectively test the ability of the present invention in dealing with such complex structure collision deformation analysis, so these two models are selected as examples. For traditional isogeometric display dynamics algorithms, since the CAD design model of the automobile B-pillar shell is often constructed by multiple non-coordinated NURBS surfaces, it is difficult to calculate. This example can better illustrate the superiority of the technology of the present invention. The automobile B-pillar shell used in this example is constructed by 65 non-coordinated NURBS surfaces, and the rigid column is constructed by two non-coordinated NURBS surfaces.

[0078] S1: Input Model 1: Automobile 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 solution analysis of the large deformation of the thin shell structure of the rigid column when the B-pillar of the car collides.

[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 thickness of the thin shell is set to 1mm; the total time for solving the dynamic analysis 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 from the NURBS unit information of the two design models obtained in S1. And assign the material information of each NURBS unit in S21, assign the boundary conditions to the material points, 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 a 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 Represents 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, 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 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 Respectively represent the size of the background grid in the x, y, and z directions.

[0093] S24: Input the current time t, and calculate the time step Δt based on the geometry display dynamic analysis such as 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 function and mapping basis function gradient between the material point generated by S22 and the background grid generated by S23. The NURBS degree information read by S1 is obtained to obtain that the NURBS degrees of the two CAD design models are both 2. Therefore, each NURBS unit is composed of 9 control points, which are interconnected to form a NURBS control network. The present invention is carried out under the framework of the GPDM material point method, so the mapping basis function between the material point and the generated background grid can be expressed as:

[0095]

[0096]

[0097] Among them, xc is the position of the NURBS unit control point, Rc(ξ,η) is the quadratic NURBS basis function corresponding to the control point, N I (x) is the weight function from material point to 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, node internal force, and node 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] Map 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] Map 4: in, represents the internal force of 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 mapping basis function gradient.

[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 on the background grid node I at time t, which can be expressed by the formula:

[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 velocity of the material point can be expressed as:

[0115]

[0116] in, Represents the velocity of the material point P after 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 material point after update, Δσ p is the Cauchy stress increment of the material point. The Cauchy stress increment is calculated by the isogeometric shell stress update theory. Resetting the background grid information to 0 is to recalculate in the next cycle to ensure the accuracy and continuity of the calculation.

[0123] S29: Update the current time t, and determine whether the total time t of the dynamics 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 solution results.

[0124] By continuously updating the current time and comparing it with the set total duration, 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 B-pillar shell of the car, as time goes by, by updating the current time t and judging whether the total duration has been reached, the entire dynamic process from the beginning to the end of the collision process can be fully simulated. When the current time t after the update is greater than or equal to the total duration t end , the kinetic analysis ends and the output shows the kinetic analysis results.

[0125] Furthermore, the iso-geometry 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 iso-geometric display dynamic analysis solution of the car B-pillar shell colliding with the rigid column is output using VTK (V format, and the solution result is converted into a text file. The output result is read by Paraview to realize the visualization of the result. It displays 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 they cannot 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 modifications and improvements can be made without departing from the concept of the present invention, which all belong to the protection scope 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 in that: 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, the number of NURBS surfaces, the degree of NURBS surface, the object parameter space grid topology information, the geometric position of the NURBS surface control vertex and its weight factor, NURBS unit information, the control point position corresponding to each NURBS unit, and the NURBS control point basis function; 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 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 by S1, initial boundary conditions of each NURBS unit obtained by S1, thickness of thin shell, and total time t for solving geometric display dynamics analysis. end , etc. Geometry shows that the dynamic analysis solves the initial time t0; S22: Generate material points from the NURBS unit information obtained in S1, and each NURBS unit corresponds to a material point one by one; assign the material information and boundary conditions of each NURBS unit obtained in S21 to the corresponding material point, and also assign the shell thickness parameter to each material point; through such operations, establish the connection between NURBS units, material points and related parameters; S23: Generate a background grid according to the geometric dimensions of the CAD design model read in S1, and set the size of the background grid unit 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 calculation such as the background grid size and the velocity calculation of the material point at the current time; S25: Calculate the interpolation basis function and the gradient of the interpolation basis function between the material points generated by S22 and the background grid generated by S23; S26: Map the physical information contained in each material point generated by S22 to the background grid generated by S23; obtain the mass, momentum, node internal force, and node external force of each background grid node; S27: solving the momentum equation on the background grid and updating the momentum of the background grid; 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 mapping to 0; S29: Update the current time t, and determine whether the total time t of the dynamics 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 results; S3: Output isogeometric display of dynamics 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 the S21, the material information of the NURBS unit includes: material density, mass, elastic modulus, Poisson's ratio, volume, bulk modulus, and shear modulus of each NURBS unit; in the S21, the initial boundary conditions of the NURBS unit include: 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 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 isogeometric display dynamics analysis adopts adaptive time step, which is calculated as follows: in, Δt is the time step of isogeometric dynamics analysis calculation; (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 means the minimum value; c x is the velocity coefficient in the x direction of the material point, 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: Where p represents a 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 taking 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, and 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 : In the formula, 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 NURBS unit area, (ξ,η) 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-coordinated complex CAD models according to claim 5, the S6 is specifically as follows: The mapped physical information described in 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, so as to establish 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; Map 3: Among them, b(x p ) is the external force described in S22 for the material point p; Map 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 material point method of GPDM, 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. 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 6, the 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 on the background grid node I at time t, which can be expressed by the formula: 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 material point P at time t; represents the resultant force on the 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 updated position of the material point P; 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-coordinated 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 isogeometric display dynamic analysis results.

Images