CN115862771A - Equal geometric analysis simulation method of super-elastic material model based on volume subdivision - Google Patents

Equal geometric analysis simulation method of super-elastic material model based on volume subdivision Download PDF

Info

Publication number
CN115862771A
CN115862771A CN202211154628.8A CN202211154628A CN115862771A CN 115862771 A CN115862771 A CN 115862771A CN 202211154628 A CN202211154628 A CN 202211154628A CN 115862771 A CN115862771 A CN 115862771A
Authority
CN
China
Prior art keywords
displacement
model
equation
point
subdivision
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202211154628.8A
Other languages
Chinese (zh)
Inventor
徐岗
王宁
黄浩
许金兰
吴海燕
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hangzhou Dianzi University
Original Assignee
Hangzhou Dianzi University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Hangzhou Dianzi University filed Critical Hangzhou Dianzi University
Priority to CN202211154628.8A priority Critical patent/CN115862771A/en
Publication of CN115862771A publication Critical patent/CN115862771A/en
Pending legal-status Critical Current

Links

Images

Landscapes

  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

The invention discloses a hyperelastic material model isogeometric analysis simulation method based on subdivision. The isogeometric analysis adopts the division of a spline model parameter domain and the mapping from the parameter domain to a physical model, avoids the time consumption of grid division in a finite element, and combines a CC subdivision method to realize an integrated process from modeling to simulation. The method comprises the following steps: converting the CC subdivision into a spline model, and establishing a control point mapping relation between the subdivision and the spline; establishing a superelasticity material constitutive model according to the geometric boundary conditions and the material parameters; establishing a discrete balance equation of the model according to an isogeometric method; combining with a balance equation, establishing a model motion equation and a nonlinear solving system; and solving the nonlinear system by using a Newton iteration method to obtain displacement solutions of each time step, and mapping each displacement solution back to the subdivision body to realize the deformation simulation of the subdivision body.

Description

Equal geometric analysis simulation method of super-elastic material model based on volume subdivision
Technical Field
The invention relates to the field of simulation of hyperelastic materials, in particular to a hyperelastic material model isogeometric analysis simulation method based on volume subdivision
Background
The invention relates to a deformation simulation process of a superelasticity material model in the field of numerical simulation, in particular to a motion deformation process of a Catmull-Clark (CC) subdivision model under stress solved by using an isogeometric analysis method.
Simulation of superelastic models has been a research hotspot of computer graphics, and in industrial manufacturing, numerical simulation calculation is usually performed on the models by using a finite element method, and the simulation precision is determined by the number and quality of grid divisions. However, the divided grids are approximate to the original precise geometry, so that the model has a certain degree of data distortion after grid division, and for improving the precision of a complex model, the grids are often refined to ensure the simulation quality, which also leads to greatly improved simulation time. The isogeometric analysis is based on the idea of finite element dispersion, non-uniform rational B splines (NURBS) adopted by a CAD model are used as a shape function in the analysis process, the traditional grid division process is skipped, the accurate geometry is directly subjected to simulation analysis, and the simulation precision is greatly improved. The NURBS function used by the isogeometric analysis has high-order continuity and can meet the requirement of high-order differentiation in engineering analysis, and the accurate geometric description enables the coarsest discrete grid to have better simulation analysis results than finite elements.
Body sample bar modeling of complex geometric models is a key problem in isogeometric analysis, and the corresponding spline bodies can be defined by a set of control grids and a set of refinement rules. These rules are applied to a given control grid to generate a new refined grid. Catmull-Clark (CC) volume subdivision for an input hexahedral mesh, an unstructured three-variable sample volume can be directly generated through two simple steps to serve as an approximation of a subdivision limit volume. Compared with the traditional approach subdivision, the method is simple and convenient to implement, the required resource consumption is low, and the refined grids with any precision can be quickly generated.
The research on geometric simulation analysis mainly focuses on statics, hydrodynamics, plate and shell problems, electromagnetism and the like. Deformation simulation of a superelasticity material model belongs to nonlinear large deformation simulation analysis, and isogeometry is less researched in the field.
Disclosure of Invention
In order to provide a rapid and effective superelasticity material deformation simulation method facing CC subdivision bodies, the invention uses isogeometric analysis to replace finite element analysis, carries out simulation analysis on a spline model converted from the CC subdivision bodies, can calculate the displacement of each time step aiming at a complex three-dimensional geometric model, and realizes the deformation simulation of the superelasticity material model.
The technical scheme adopted by the invention is as follows: converting the CC subdivision into a spline model, and establishing a control point mapping relation between the subdivision and the spline; establishing a constitutive model of the superelastic material according to the geometric boundary conditions and the material parameters; establishing a discrete balance equation of an isogeometric method; establishing a model motion equation and a nonlinear solving system; solving the nonlinear system by using a Newton iteration method, and demapping the displacement back to the subdivision body, wherein the specific steps are as follows:
step 1, converting the CC subdivision body into a B spline body model, and establishing a mapping relation between the subdivision body and each control point of the spline body, wherein the mapping relation is as follows:
giving a complex multi-sheet hexahedron mesh model, recording the corresponding control mesh as H, generating a spline body of each hexahedron unit by a three-variable spline approximation method, wherein the control point of the spline body is H', and firstly calculating the coordinate v of an internal point of the mesh in The formula is as follows:
Figure BDA0003857970520000021
wherein upsilon is any point in the hexahedral mesh, and the field of the hexahedral mesh unit in which upsilon is located is L = v, e 1 ,e 2 ,e 3 ,f 1 ,f 2 ,f 3 And c, performing the following steps. Wherein e t ,f t And c represents the coordinates of the adjacent side point, the adjacent surface point and the adjacent block point of v. The adjacent edge points are points on the same edge with the upsilon, and 3 adjacent edge points are arranged in one hexahedron grid unit; the abutment points are opposite vertices on the same plane as v, and there are 3 abutment points in a hexahedron(ii) a The adjacent block points are diagonal points which are on the same hexahedral grid unit with upsilon, and 1 adjacent block point is arranged in a hexahedron.
From the calculated coordinates of the interior points, new coordinates of edge, face or corner points, i.e. coordinates of adjacent interior points, are obtained by averaging
Figure BDA0003857970520000022
Wherein n is the number of hexahedral mesh units containing edge points, face points or corner points to be solved,
Figure BDA0003857970520000023
the coordinate of the internal point which is closest to the edge point, the face point or the angular point to be solved on the t hexahedron grid unit containing the edge point, the face point or the angular point to be solved.
The calculation formula is represented by a mapping matrix V, and the mapping relation of two groups of control points is obtained:
H=VH′ (3)
step 2, establishing a stress-strain relation according to the input material attribute and the boundary condition, and selecting a calculation form corresponding to the superelasticity model, wherein the calculation form is as follows:
inputting the attribute of the super-elastic material, and establishing a stress-strain relation according to a material model. Different from the traditional linear elasticity problem, the stress strain of the super-elastic material is in a nonlinear relation, an elastic potential energy function related to the strain exists, and the model stress is obtained by derivation of the potential energy function on the strain. Using st.venant-Kirchhoff superelastic material as an example, the corresponding elastic potential energy function is:
Ψ(E)=λtr(E)+μtr(E 2 ) (4)
in the formula, Ψ (E) is the magnitude of the elastic potential energy when the model deforms, E represents the relative deformation magnitude, i.e., strain, and λ and μ are the input material coefficients. Obtaining a stress component by solving the partial derivative of the strain component through an elastic potential energy function, wherein the stress component comprises the following steps:
Figure BDA0003857970520000031
wherein S is the Schiff stress tensor, and I is an identity matrix. According to the constitutive relation of stress and strain, calculating the derivative of stress and strain to obtain a tangent material matrix D T
Figure BDA0003857970520000032
Strain refers to the relative deformation of the object node, and can be expressed as:
Figure BDA0003857970520000033
wherein
Figure BDA0003857970520000034
Representing the displacement gradient of the object.
The input boundary conditions are generally two, namely a displacement boundary condition (Dirichlet condition) and an external force condition (Noeman condition), and the node displacement value is constrained according to the two conditions, so that the solution is unique.
Step 3, establishing a discrete balance equation by using an isogeometric method
And (3) taking the spline basis function input in the step (1) as a basis function of a solution space by using an isogeometric method, and performing interpolation solution on discrete points in the spline unit by using the basis function to obtain a numerical solution under accurate geometric modeling expression. Firstly, establishing a cell balance equation, and linearly representing the displacement of the cell discrete point by the displacement on each control point according to the representation of the spline basis function, wherein the displacement is as follows:
u(ξ)=∑ i N i u i =Nu e (8)
in the formula, u (ξ) represents the magnitude of displacement at any point of the model, and u (ξ) represents the magnitude of displacement at any point of the model i Indicating the magnitude of the displacement of the control point i, N i The size of the basis function of the control point i is shown, and N is a matrix of the combination of the basis functions. From the control point displacement, a representation of the strain coefficient matrix B of the object can be derived from equation (7):
Bu e =E=LNu e (9)
l is a differential operator, and the above equation rewrites the strain into a tensor form. According to the Lagrange virtual work principle, the sum of virtual work suffered by the system in a balanced state is zero, and the following unit balance equation is established:
d eΩ B T SdΩ=d eΩ N T fdΩ+∫ Γ N T tdΓ (10)
in the formula, d e Representing the virtual displacement, f the received volumetric force, and t the received area force. The left side of the equation shows that the internal force of the model does virtual work on virtual displacement, the right side shows that the external force does virtual work on virtual displacement, and the object keeps balance, so that the two virtual works are mutually counteracted to obtain (10). And assembling the balance equations of the units into a balance equation of the whole system, and eliminating the virtual displacement to obtain the balance equation of the system. Since numerical calculations need to be performed in the parameter domain, the above equation is subjected to the jacobian transformation, i.e., Ω → Ω', to obtain a calculable system balance equation:
Ω′ B T S|J|dΩ′=∫ Ω′ N T f|J|dΩ′+∫ Γ′ N T t|J|dΓ′ (11)
where | J | is the jacobian determinant of the physical domain to parametric domain transform, the above equation can also be abbreviated as:
Ku=f (12)
k is a rigidity matrix, f is a load vector, u is the solved displacement solution, and the displacement solution of the model static analysis can be obtained by solving the equation.
Step 4, establishing a model motion equation and establishing a nonlinear solving system by using implicit time integration, wherein the method specifically comprises the following steps:
when dealing with the process of model deformation, the state variables and the change of deformation over time must be taken into account, and according to newton's law, the general form of the equation of motion of an object can be expressed as follows:
Figure BDA0003857970520000041
wherein M represents a mass matrix, related to the bulk density of the model; c is a damping matrix, related to the mass matrix and the model stiffness; r is an initial load vector and is related to the initial deformation of the model, and P is the external force applied to the model.
Figure BDA0003857970520000042
Figure BDA0003857970520000043
Representing the second derivative of displacement, the first derivative, i.e. acceleration and velocity, respectively.
After an object motion equation is established, a velocity displacement recurrence relation between time steps is established through Newmark implicit time integration to time dispersion, and the expression form is as follows:
Figure BDA0003857970520000044
Figure BDA0003857970520000045
wherein, Δ t is a time step after dispersion, and the smaller the step is, the higher the solving precision is; beta and gamma are self-defined parameters, which determine the form of integral convergence. Equations (14) and (15) are the n-th step displacement derived from the n + 1-th steps of velocity and acceleration, and the displacement, velocity and acceleration of the n-th step are known in the solution, so the above equations need to be combined to obtain a new derivation:
Figure BDA0003857970520000051
Figure BDA0003857970520000052
wherein, Δ u n+1 =u n+1 -u n When two times are shownStep-by-step displacement difference, alpha 1~6 Is a transformed parameter, which is a constant term. The solution of the next time step is then determined by the solution of the previous time step and the variation of the displacement, the unknowns being only Δ u n+1 By substituting (16) (17) into the object motion equation (13), a nonlinear algebraic equation of unknown displacement is obtained:
Figure BDA0003857970520000053
equation (18) is a system of N equations at N time steps, and by solving the system of equations, a displacement solution is obtained for each time step. Then, each group of displacement is demapped back to the control points of the model, and the motion state and the deformation degree of each time-step spline model are obtained.
Step 5, solving a nonlinear equation system and mapping the displacement solution of each time step back to the control point corresponding to the CC subdivision body:
for equation (18) at a certain time step, the displacement variation delta u is solved by using a Newton iteration method n+1 The tangential stiffness matrix to be solved for each sub-iteration is as follows:
Figure BDA0003857970520000054
this matrix can be obtained by derivation of the equilibrium equation (11) obtained in step 3. Then, the following scheme is performed for each Newmark time step iteration:
Figure BDA0003857970520000055
Figure BDA0003857970520000056
wherein the initial value of the iteration is the convergence value of the last iteration:
Figure BDA0003857970520000057
and after iterative convergence, calculating the displacement, the speed and the acceleration of the time step, substituting the values into the next time step to continuously solve to obtain the displacement solution of each time step. And (3) obtaining the displacement on each control point of the spline body, and finally obtaining the motion state and the deformation size of the subdivision by using the subdivision and spline body control point mapping formula (3) in the step (1) to finish motion simulation.
The invention has the advantages that the deformation simulation analysis of the superelasticity model is realized by using the isogeometric analysis method, and compared with the traditional finite element analysis, the method has higher simulation precision and efficiency on the model with the same fineness. The method is directly butted with the model after CC subdivision, thereby realizing the deformation simulation of the complex CC subdivision model and perfecting the integrated process from modeling to simulation. The invention constructs the isogeometric framework for numerical simulation of the CC subdivision model, and the framework is also suitable for structural optimization of the model.
The invention converts the CC subdivision model into a spline model with accurate geometric representation, then carries out analog simulation by using isogeometric analysis, directly connects the modeling and simulation processes, and realizes the integrated flow of modeling and simulation. Moreover, for the simulation of large deformation of the superelastic material, the finite element analysis often causes the precision to be reduced due to the distortion of the elements or the solving time to be too long due to the dense mesh. The special high-precision high-continuity characteristics of the geometric model can well capture a deformation area, and a complex deformation process can be represented only by a few degrees of freedom, so that a new idea is provided for high-efficiency high-quality analog simulation.
Drawings
FIG. 1a and FIG. 1b are different views of the hexahedral mesh model of the bridge, respectively;
FIG. 2a is a view of FIG. 2b showing different views of the bridge spline model;
FIG. 3 is a distribution of spline blocks of a bridge model;
FIG. 4 is the result of the present invention after performing isogeometric analysis on the model of FIG. 3, wherein different gray scales represent the magnitude of the displacement solution;
FIG. 5 is a front result diagram of the dynamic deformation of the bridge model under stress according to the method of the present invention, and FIGS. 5a,5b, and 5c are three moments intercepted during the deformation process;
FIG. 6 is a side result graph of the dynamic deformation of the bridge model under stress in the method of the present invention, and FIGS. 6a,6b, and 6c are three times intercepted during the deformation process;
FIG. 7 is a graph of the results of the dynamic deformation of the model C under force according to the method of the invention, FIG. 7a,7b,7C being three moments intercepted by the deformation process;
FIG. 8 is a result diagram of dynamic deformation of a model tentacle under stress in the method of the present invention, and FIG. 8a,8b and 8c are three moments intercepted in the deformation process;
fig. 9 is an overall flowchart.
Detailed Description
The invention will be further explained with reference to the drawings.
An equal geometric deformation simulation method facing CC subdivision comprises the following specific steps:
1. establishing a mapping relation between the subdivision body and each control point of the sample strip body as follows
Giving a hexahedral mesh model (as shown in fig. 1a and 1 b), recording the corresponding control mesh as H, and generating an approximated spline body (as shown in fig. 2a and 2 b) for each hexahedral unit by a three-variable spline approximation method, wherein the control points of the spline body are H', and the calculation method is as follows:
Figure BDA0003857970520000071
wherein upsilon is any point in the hexahedral mesh, and the field of the hexahedral mesh unit in which upsilon is located is L = v, e 1 ,e 2 ,e 3 ,f 1 ,f 2 ,f 3 And c, c. Wherein e t ,f t And c represents the coordinates of the adjacent side point, the adjacent surface point and the adjacent block point of v. The adjacent edge points are points on the same edge with the upsilon, and 3 adjacent edge points are arranged in one hexahedron grid unit; the adjacent surface points are opposite vertex points on the same surface with the adjacent surface points, and 3 adjacent surface points are arranged in a hexahedron; point of adjacent blockIs at the diagonal point of the same hexahedral mesh unit as υ, and there are 1 contiguous block points in one hexahedron.
From the calculated coordinates of the interior points, new coordinates of edge, face or corner points, i.e. coordinates of adjacent interior points, are obtained by averaging
Figure BDA0003857970520000072
Wherein n is the number of hexahedral mesh units containing edge points, face points or corner points to be solved,
Figure BDA0003857970520000073
the coordinate of the internal point which is closest to the edge point, the face point or the angular point to be solved on the t hexahedron grid unit containing the edge point, the face point or the angular point to be solved. />
And expressing the calculation formula by a mapping matrix V to obtain the mapping relation of the control points of the two control points:
H=VH′ (3)
2. inputting material properties and boundary conditions, and establishing a constitutive model of the super-elastic material:
inputting the attribute of the super-elastic material, and establishing a stress-strain relation according to a material model. Different from the traditional linear elasticity problem, the stress strain of the super-elastic material is in a nonlinear relation, an elastic potential energy function related to the strain exists, and the model stress is obtained by derivation of the potential energy function on the strain. Using st.venant-Kirchhoff superelastic material as an example, the corresponding elastic potential energy function is:
Ψ(E)=λtr(E)+μtr(E 2 ) (4)
in the formula, Ψ (E) is the magnitude of the elastic potential energy when the model deforms, E represents the relative deformation magnitude, i.e., strain, and λ and μ are the input material coefficients. Obtaining a stress component by solving a partial derivative of the strain component through an elastic potential energy function (1), wherein the stress component comprises the following steps:
Figure BDA0003857970520000074
where S is the kSchiff stress tensor and I is an identity matrix. According to the constitutive relation of stress and strain, calculating the derivative of stress and strain to obtain a tangent material matrix D T
Figure BDA0003857970520000081
Strain refers to the magnitude of the relative deformation of the object's node, and can be expressed as:
Figure BDA0003857970520000082
wherein
Figure BDA0003857970520000083
Representing the displacement gradient of the object.
The input boundary conditions are generally two, namely a displacement boundary condition (Dirichlet condition) and an external force condition (Noeman condition), and the node displacement value is constrained according to the two conditions, so that the solution is unique.
3. Establishing discrete equilibrium equation by isogeometric method
And (3) taking the spline basis function input in the step (1) as a basis function of a solution space by using an isogeometric method, and performing interpolation solution on discrete points in the spline unit by using the basis function to obtain a numerical solution under accurate geometric modeling expression. Firstly, establishing a cell balance equation, and linearly representing the displacement of the cell discrete point by the displacement on each control point according to the representation of the spline basis function, wherein the displacement is as follows:
u(ξ)=∑ i N i u i =Nu e (8)
in the formula, u (xi) represents the displacement of any point of the model, and u (xi) represents i Indicating the magnitude of the displacement of control point i, N i The size of the basis function of the control point i is shown, and N is a matrix of the combination of the basis functions. From the control point displacement, a representation of the strain coefficient matrix B of the object can be derived from equation (7):
Bu e =E=LNu e (9)
l is a differential operator, and the above equation rewrites the strain into a tensor form. According to the Lagrange virtual work principle, the sum of virtual work suffered by the system in the equilibrium state is zero, and the following unit equilibrium equation is established:
d eΩ B T SdΩ=d eΩ N T fdΩ+∫ Γ N T tdΓ (10)
in the formula (d) e Representing the virtual displacement, f the received volumetric force, and t the received area force.
The left side of the equation shows that the internal force of the model does virtual work on virtual displacement, the right side shows that the external force does virtual work on virtual displacement, and the object keeps balance, so that the two virtual works are mutually counteracted to obtain (10). And assembling the balance equations of the units into an integral system balance equation, and eliminating the virtual displacement to obtain the system balance equation. Since numerical calculations need to be performed in the parameter domain, the above equation is subjected to the jacobian transformation, i.e., Ω → Ω', resulting in a calculable system equilibrium equation:
Ω′ B T S|J|dΩ′=∫ Ω′ N T f|J|dΩ′+∫ Γ′ N T t|J|dΓ′ (11)
where | J | is the jacobian determinant of the physical domain to parametric domain transform, the above equation can also be abbreviated as:
Ku=f (12)
k is a rigidity matrix, f is a load vector, u is the solved displacement solution, and the displacement solution of the model static analysis can be obtained by solving the equation, as shown in FIG. 4.
4. Establishing a model equation of motion and establishing a nonlinear solving system by using implicit time integral, which is as follows
When dealing with the process of model deformation, the state variables and the change of deformation over time must be taken into account, and according to newton's law, the general form of the equation of motion of an object can be expressed as follows:
Figure BDA0003857970520000091
wherein M represents a mass matrix, related to the bulk density of the model; c is a damping matrix, related to the mass matrix and the model stiffness; r is an initial load vector and is related to the initial deformation of the model, and P is the external force applied to the model.
Figure BDA0003857970520000092
Representing the second derivative of displacement, the first derivative, i.e. acceleration and velocity, respectively.
After an object motion equation is established, a velocity displacement recurrence relation between time steps is established through Newmark implicit time integration to time dispersion, and the expression form is as follows:
Figure BDA0003857970520000093
Figure BDA0003857970520000094
wherein, Δ t is a time step after dispersion, and the smaller the step is, the higher the solving precision is; beta and gamma are self-defining parameters and determine the form of integral convergence. Equations (14) and (15) are the displacement of the nth step derived from the (n + 1) th steps of velocity and acceleration, and the displacement, velocity and acceleration of the nth step are known in the solution, so the above equation shift terms need to be combined to obtain a new derived equation:
Figure BDA0003857970520000095
Figure BDA0003857970520000096
wherein, Δ u n+1 =u n+1 -u n Representing the difference in displacement, alpha, of two time steps before and after 1~6 Is a transformed parameter and is a constant term. The solution of the next time step is then determined by the solution of the previous time step and the amount of change in the displacementThe unknown quantity is only Deltau n+1 By substituting (16) (17) into the object motion equation (13), a nonlinear algebraic equation of unknown displacement is obtained:
Figure BDA0003857970520000097
equation (18) is a system of N equations at N time steps, and by solving the system of equations, a displacement solution is obtained for each time step. Each set of displacements is then de-mapped back into the control points of the model, resulting in the state of motion and the degree of deformation for each time step spline model, as shown in figures 5a,5b, 5c.
5. Solving a nonlinear equation system and mapping the displacement solution of each time step back to the corresponding node of the CC subdivision body
For equation (18) at a certain time step, the displacement variation delta u is solved by using a Newton iteration method n+1 The tangential stiffness matrix to be solved for each sub-iteration is as follows:
Figure BDA0003857970520000101
this matrix can be obtained by deriving the equilibrium equation obtained in step 3. Then, the following scheme is performed for each Newmark time step iteration:
Figure BDA0003857970520000102
Figure BDA0003857970520000103
wherein the initial value of the iteration is the convergence value of the last iteration:
Figure BDA0003857970520000104
calculating the displacement, the speed and the acceleration of the time step after iterative convergence, substituting the values into the next time step to continuously solve to obtainThe shift solution at each time step. And (4) obtaining the displacement size of each control point of the spline body, and finally obtaining the motion state and the deformation size of the subdivision body by using the subdivision body and the spline body control point mapping formula (3) in the step 1 to finish motion simulation.
The examples are as follows: fig. 5 is a diagram showing a result of dynamic deformation of the bridge model under an external force, fig. 7 is a diagram showing a result of dynamic deformation of a letter C, and fig. 8 is a diagram showing a result of dynamic deformation of the tentacle model. The method can use isogeometry to carry out numerical simulation on the complex multi-chip model and solve the motion state of the model, and can solve the problem faster than the method under the same precision.

Claims (6)

1. A geometric analysis simulation method based on a subdivision hyperelastic material model and the like is characterized by comprising the following steps:
s1, converting a grid model of the CC subdivision body into a spline body model, and establishing a control point mapping relation between the CC subdivision body and the spline body;
s2, establishing a constitutive model reflecting the stress-strain relation of the super-elastic material according to the geometric boundary conditions and the material parameters;
s3, establishing a discrete balance equation of an isogeometric method by using the constitutive model in the step 2, and solving the discrete balance equation to obtain a displacement solution of the spline model static analysis;
s4, establishing a model motion equation, establishing a speed and displacement recurrence relation between time steps through Newmark implicit time integration to time dispersion, and substituting the speed and displacement recurrence relation between the time steps into the model motion equation to obtain a nonlinear solving system;
and S5, solving the nonlinear system by using a Newton iteration method to obtain a displacement solution of each time step, mapping the displacement solution of each time step back to the control point corresponding to the CC subdivision body, obtaining the motion state and the deformation size of the subdivision body, and completing motion simulation.
2. The method for simulating geometric analysis based on a subdivided superelastic material model according to claim 1, wherein in S1,
the grid model of the CC subdivision is a complex multi-slice hexahedron grid model of the CC subdivision;
the spline body model is a spline body model of a hexahedral unit generated by a three-variable spline approximation method;
the control point mapping relation is formed by the following steps:
giving a complex multi-sheet hexahedral mesh model, recording a corresponding control mesh as H, and generating a spline body of each hexahedral unit by a three-variable spline approximation method, wherein the control point of the spline body is H', and the calculation mode is as follows:
Figure FDA0003857970510000011
wherein upsilon is any point in the hexahedral mesh, and the field of the hexahedral mesh unit in which upsilon is located is L = v, e 1 ,e 2 ,e 3 ,f 1 ,f 2 ,f 3 C, c; wherein e t ,f t C represents the coordinates of the adjacent side point, the adjacent surface point and the adjacent block point of v respectively; the adjacent edge points are points on the same edge with the upsilon, and 3 adjacent edge points are arranged in one hexahedron grid unit; the adjacent surface points are opposite vertexes on the same surface with the upsilon, and 3 adjacent surface points are arranged in a hexahedron; the adjacent block points are diagonal points which are on the same hexahedral grid unit with upsilon, and 1 adjacent block point is arranged in a hexahedron;
from the calculated coordinates of the interior points, new coordinates of edge, face or corner points, i.e. coordinates of adjacent interior points, are obtained by averaging
Figure FDA0003857970510000021
Wherein n is the number of hexahedral mesh units containing edge points, face points or corner points to be solved,
Figure FDA0003857970510000022
to compriseCoordinates of an inner point which is closest to the edge point, the face point or the angular point to be solved on the t hexahedron grid unit of the edge point, the face point or the angular point to be solved; and expressing the calculation formula by using a mapping matrix V to obtain the mapping relation of two groups of control points:
M=VM′ (3)。
3. the method for simulating geometric analysis based on subdivided superelastic material model according to claim 1, wherein in S2,
the super-elastic material is: venant-Kirchhoff superelastic material,
the material parameters include elastic potential energy function:
Ψ(E)=λtr(E)+μtr(E 2 ) (4)
in the formula, psi (E) is the elastic potential energy when the model deforms, E represents the relative deformation magnitude, namely strain, and lambda and mu are constants obtained by input material coefficients;
obtaining a stress component by solving the partial derivative of the strain component through an elastic potential energy function, wherein the stress component comprises the following steps:
Figure FDA0003857970510000023
wherein S is a kSchiff stress tensor, and I is an identity matrix;
according to the constitutive relation of stress and strain, calculating the derivative of stress and strain to obtain a tangent material matrix D T
Figure FDA0003857970510000024
Strain refers to the magnitude of the relative deformation of the object's node, and can be expressed as:
Figure FDA0003857970510000025
wherein
Figure FDA0003857970510000026
Representing a displacement gradient of the object;
the geometric boundary conditions include:
boundary conditions reflecting displacement: dirichlet conditions
Boundary conditions reflecting external forces: under the conditions of the Noelman (R) state,
and constraining the node displacement value according to the two boundary conditions to obtain a unique solution.
4. The method for simulating geometric analysis based on subdivided superelastic material model according to claim 3, wherein in S3,
the establishment of the discrete balance equation of the isogeometric method comprises the following substeps:
taking the spline basis function of the spline model generated in the step 1 as a basis function of a solution space,
establishing a cell balance equation, wherein the displacement of the cell discrete point is linearly represented by the displacement on each control point, and the following steps are carried out:
u(ξ)=∑ i N i u i =Nu e (8)
in the formula, u (ξ) represents the magnitude of displacement at any point of the model, and u (ξ) represents the magnitude of displacement at any point of the model i Indicating the magnitude of the displacement of the control point i, N i Representing the size of a basis function of the control point i, wherein N is a matrix of basis function combination; from the control point displacement, a representation of the strain coefficient matrix B of the object is derived from equation (7):
Bu e =E=LNu e (9)
wherein L is a differential operator, the above formula rewrites the strain into a tensor form, and according to the Lagrange virtual work principle, the virtual work sum suffered by the system in an equilibrium state is zero, and the following unit balance equation is established:
d eΩ B T SdΩ=d eΩ N T fdΩ+∫ Γ N T tdΓ (10)
wherein d is e Representing a virtual displacement, f-tableThe volume force is shown, t is the area force, the left side of the equation represents the virtual work of the internal force of the model on the virtual displacement, the right side represents the virtual work of the external force on the virtual displacement, the object is kept balanced, and the two virtual works are mutually counteracted;
assembling the unit balance equations into an integral system balance equation, eliminating the virtual displacement to obtain a system balance equation, and obtaining a computable system balance equation as a discrete balance equation through Jacobian transformation, namely omega → omega':
Ω′ B T S|J|dΩ′=∫ Ω′ N T f|J|dΩ′+∫ Γ′ N T t|J|dΓ′ (11)
in the formula, | J | is jacobian of the transformation between the physical domain and the parameter domain, and the above formula is abbreviated as:
Ku=f (12)
k is a stiffness matrix, f is a load vector, and u is the solved displacement solution;
the solution of the discrete balance equation comprises the following sub-steps:
and carrying out interpolation solution on discrete points in the unit by using the spline basis function to obtain a numerical solution under accurate geometric modeling expression as the displacement solution of the static analysis of the spline model.
5. The method for simulating geometric analysis based on subdivided superelastic material model according to claim 4, wherein in S4,
the model equation of motion is of the form:
Figure FDA0003857970510000031
wherein M represents a mass matrix, related to the bulk density of the model; c is a damping matrix, related to the mass matrix and the model stiffness; r is an initial load vector and is related to the initial deformation of the model, and P is the external force applied to the model;
Figure FDA0003857970510000041
a second derivative and a first derivative representing displacement respectively, and representing acceleration and velocity respectively;
the method for establishing the velocity displacement recurrence relation between time steps through the dispersion of Newmark implicit time integration to time specifically comprises the following steps of:
and (3) establishing a speed displacement recurrence relation between time steps by using Newmark implicit time integration to time dispersion, wherein the expression form is as follows:
Figure FDA0003857970510000042
Figure FDA0003857970510000043
wherein, Δ t is a time step after dispersion, and the smaller the step is, the higher the solving precision is; beta and gamma are self-defined parameters and determine the form of integral convergence; equations (14) and (15) are the displacement of the nth step derived from the (n + 1) th steps of velocity and acceleration, and the displacement, velocity and acceleration of the nth step are known in the solution, so the above equation shift terms need to be combined to obtain a new derived equation:
Figure FDA0003857970510000044
Figure FDA0003857970510000045
wherein, Δ u n+1 =u n+1 -u n Representing the difference in displacement, alpha, of two time steps before and after 1~6 Is a transformed parameter, is a constant term,
the method for obtaining the nonlinear solving system by substituting the velocity displacement recurrence relation between the time steps into a model motion equation specifically comprises the following steps:
the solution of the next time step is determined by the solution of the previous time step and the displacement variation, and a nonlinear algebraic equation of unknown displacement is obtained by substituting (16) and (17) into the object motion equation (13):
Figure FDA0003857970510000046
equation (18) is a system of N equations at N time steps.
6. The method for geometric analysis simulation based on a subdivided superelastic material model according to claim 5, wherein S5 specifically comprises the following steps:
method for solving displacement variation delta u by utilizing Newton iteration method n+1 The tangential stiffness matrix to be solved for each sub-iteration is as follows:
Figure FDA0003857970510000051
the matrix may be obtained by deriving the equilibrium equation obtained in step S3;
the following scheme is performed for each Newmark time step iteration:
Figure FDA0003857970510000052
Figure FDA0003857970510000053
wherein the initial value of the iteration is the convergence value of the last iteration:
Figure FDA0003857970510000054
calculating the displacement, speed and acceleration of the time step after iterative convergence, substituting the calculated values into the next time step to continue solving to obtain the position of each time stepPerforming dismissal;
and (3) obtaining the displacement on each control point of the spline body, and finally obtaining the motion state and the deformation size of the subdivision by using the subdivision and spline body control point mapping formula (3) in the step S1 to finish motion simulation.
CN202211154628.8A 2022-09-22 2022-09-22 Equal geometric analysis simulation method of super-elastic material model based on volume subdivision Pending CN115862771A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211154628.8A CN115862771A (en) 2022-09-22 2022-09-22 Equal geometric analysis simulation method of super-elastic material model based on volume subdivision

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211154628.8A CN115862771A (en) 2022-09-22 2022-09-22 Equal geometric analysis simulation method of super-elastic material model based on volume subdivision

Publications (1)

Publication Number Publication Date
CN115862771A true CN115862771A (en) 2023-03-28

Family

ID=85661079

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211154628.8A Pending CN115862771A (en) 2022-09-22 2022-09-22 Equal geometric analysis simulation method of super-elastic material model based on volume subdivision

Country Status (1)

Country Link
CN (1) CN115862771A (en)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116757026A (en) * 2023-06-08 2023-09-15 上海交通大学 Plate frame structure ultimate strength analysis implementation method based on isogeometric analysis

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116757026A (en) * 2023-06-08 2023-09-15 上海交通大学 Plate frame structure ultimate strength analysis implementation method based on isogeometric analysis
CN116757026B (en) * 2023-06-08 2024-04-05 上海交通大学 Plate frame structure ultimate strength analysis implementation method based on isogeometric analysis

Similar Documents

Publication Publication Date Title
Cirak et al. Integrated modeling, finite-element analysis, and engineering design for thin-shell structures using subdivision
Dhatt et al. Finite element method
Benson et al. Isogeometric shell analysis: the Reissner–Mindlin shell
Liew et al. A review of meshless methods for laminated and functionally graded plates and shells
Schramm et al. The coupling of geometric descriptions and finite elements using NURBs—A study in shape optimization
CN110543654B (en) Method for determining distributed process parameters of laser shot blasting forming complex curved surface
Illoul et al. On some aspects of the CNEM implementation in 3D in order to simulate high speed machining or shearing
Somireddy et al. Meshless natural neighbor Galerkin method for the bending and vibration analysis of composite plates
Ding et al. An efficient dynamic mesh generation method for complex multi-block structured grid
CN107038270B (en) Method for calculating machining deformation caused by surface machining residual stress field
Rendall et al. Improved radial basis function fluid–structure coupling via efficient localized implementation
CN115862771A (en) Equal geometric analysis simulation method of super-elastic material model based on volume subdivision
Pagani et al. Cross-sectional mapping for refined beam elements with applications to shell-like structures
Zhan et al. A high efficient surface-based method for predicting part distortions in machining and shot peening
CN112035980A (en) Construction method of isogeometric mixed Kirchhoff-Love shell unit
Cui et al. A high-order edge-based smoothed finite element (ES-FEM) method with four-node triangular element for solid mechanics problems
Jonsson et al. Development of flutter constraints for high-fidelity aerostructural optimization
CN107145630B (en) Plate shell structure design and analysis integrated method based on CAD (computer-aided design) cutting curved surface
CN111177960B (en) Thin shell collision contact calculation method based on ANCF
Karimian et al. Immersed boundary method for the solution of 2d inviscid compressible flow using finite volume approach on moving cartesian grid
Ding et al. Isogeometric topology optimization of compliant mechanisms using transformable triangular mesh (TTM) algorithm
van Schie et al. Solver-independent aeroelastic coupling for large-scale multidisciplinary design optimization
EP4246362A1 (en) Modelling an object, determination of load capacity, improvement of the design and generating component, system
Neumann et al. Coupling strategies for large industrial models
Recio et al. Locking and hourglass phenomena in an element‐free Galerkin context: the B‐bar method with stabilization and an enhanced strain method

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination