Non-intrusive transient contact numerical calculation method based on penalty function and conjugate gradient
By combining the penalty function and conjugate gradient methods into a non-intrusive transient contact numerical calculation method, the problem of non-physical penetration in transient contact calculation using the penalty function method is solved, achieving higher simulation accuracy and stability, and making it applicable to a wide range of collision problems.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2025-12-30
- Publication Date
- 2026-04-07
AI Technical Summary
Existing penalty function methods cannot strictly satisfy contact boundary conditions in transient contact calculations, leading to non-physical penetration and divergent calculation results. In particular, errors accumulate severely in high-speed collision problems, affecting simulation accuracy and stability.
A non-intrusive transient contact numerical calculation method based on penalty function and conjugate gradient is adopted. By estimating the contact gap and displacement, and combining the penalty function method and conjugate gradient method for iterative solution, the contact calculation is ensured to meet the non-intrusive condition. The penalty function method is modified by the preprocessed conjugate gradient method to improve convergence and accuracy.
It effectively avoids non-physical penetration and computational divergence, improves the simulation accuracy of transient collision problems and the prediction accuracy and stability of structural impact dynamics simulation, and is applicable to high-speed collision problems with rigidity and large plastic deformation.
Smart Images
Figure CN121809168A_ABST
Abstract
Description
Technical Field
[0001] This application relates to the field of numerical simulation of structural impact dynamics, and in particular to a non-intrusive transient contact numerical calculation method based on penalty functions and conjugate gradients. Background Technology
[0002] The transient contact algorithm is one of the most core algorithms in the numerical simulation technology of structural impact dynamics. It directly affects the accuracy of the collision impact load obtained by simulation, and is also related to the accuracy and stability of the entire structural impact event simulation.
[0003] For the calculation of finite deformation contact collision problems, two main algorithms are the Lagrange multiplier method and the penalty function method. The Lagrange multiplier method introduces constraints using Lagrange multipliers, enabling a rigorous solution for the contact boundary. However, its implicit iteration often fails to converge in strongly nonlinear problems such as transient collisions, and even in large deformation contact analysis of simple structures, making it unsuitable for impact collision analysis of engineering structures like metals. The penalty function method applies a penalty force to the already intruded contact boundary to approximate the non-intrusion condition. This method is widely used in collision engineering due to its high computational efficiency and stability. However, its inability to strictly satisfy the contact boundary condition leads to non-physical penetration, significantly impacting the accuracy of the calculation results. Especially in high-speed collision problems, computational errors accumulate gradually, and penetration increases rapidly over time, causing the calculation to diverge quickly. Summary of the Invention
[0004] To address the aforementioned problems and technical requirements, this application proposes a non-intrusive transient contact numerical calculation method based on penalty functions and conjugate gradients. The technical solution of this application is as follows: A non-intrusive transient contact numerical calculation method based on penalty functions and conjugate gradients includes the following steps: For any (n+1)th calculation time step, based on the motion state parameters of the first and second colliding objects at the nth calculation time step, the predicted displacement values of the first and second colliding objects at the (n+1)th calculation time step are estimated. ; Determine the contact gap between the first and second colliders at the (n+1)th computation time step. And based on the contact gap and displacement prediction value Determine the predicted contact gap value at the (n+2)th computation time step. ; Based on the predicted contact gap value and displacement prediction value The penalty function method is used to calculate the contact force between the first and second colliding objects, and the estimated contact force is obtained. Based on the estimated contact force and displacement prediction value The displacement at the (n+1)th calculation time step is updated. and the contact gap at the (n+2)th calculation time step ; When the updated contact gap If the non-intrusion condition is not met, use the updated displacement. and contact gap Using the non-intrusion condition as the initial value for iteration and the pre-processed conjugate gradient method as the constraint, the normal gap equation is solved iteratively to obtain the contact force correction value. And update the displacement at the (n+1)th calculation time step. The normal gap equation characterizes the state of zero normal contact between two colliding objects. Based on the estimated contact force and contact force correction value Determine the contact force at the (n+1)th calculation time step. And use the displacement of the (n+1)th calculation time step Update the motion state parameters of the first and second colliders, and proceed to the next calculation time step.
[0005] A further technical solution involves specifying that the motion state parameters include velocity and acceleration, and calculating the velocity and acceleration at the nth calculation time step in time increments. Numerical integration yields the predicted displacement values of the first and second colliding objects. Time increment It is the time interval between the nth computation time step and the (n+1)th computation time step.
[0006] A further technical solution involves determining the predicted contact gap value at the (n+2)th calculation time step. ; The penalty function method is used to calculate the contact force between the first and second colliding objects, and the estimated contact force is obtained. ,in, It is the outward normal vector of the main facet. It is the elastic force coefficient. It is the damping force coefficient. It is a penalty factor. It is the spatial gradient of the contact gap and , express transpose, Indicates contact gap For position coordinates Perform differentiation. This indicates taking the maximum value.
[0007] A further technical solution involves updating the displacement at the (n+1)th calculation time step. The contact gap at the (n+2)th calculation time step is updated. ; in, , It is a time increment. Represents the inverse of a matrix; It is the spatial gradient of the contact gap and , express transpose, Indicates contact gap For position coordinates Perform differentiation.
[0008] A further technical solution involves using the finite element method to discretize the first and second colliding objects into multiple elements, each element including multiple nodes. The surface element of the outer surface of any one of the first and second colliding objects is taken as the principal surface element, and the nodes on the outer surface of the other colliding object are taken as slave nodes. The contact gap includes the contact gap of all contact pairs between the first and second colliding objects, with each contact pair including one principal surface element and one slave node. The normal gap equation is:
[0009] The non-invasive condition is:
[0010] in, , It is a time increment. Represents the inverse of a matrix; It is the spatial gradient of the contact gap and , express transpose, Indicates contact gap For position coordinates Perform differentiation; It is the normal gap factor. It is the outward normal vector of the main face.
[0011] A further technical solution involves using the pre-processed conjugate gradient method to iteratively solve the normal gap equation, thereby obtaining the contact force correction value. include: Constructing the preprocessing projection matrix , , , It is the spatial gradient of the contact gap. express The transpose of M, where M is the mass matrix. It is a time increment. Let A be the diagonal matrix of matrix A. The inverse of the matrix is represented; the preprocessed projection matrix is used to project the normal gap equation onto the constraint space for solution. Initialize integer parameter k=1, initialize contact gap Initialize contact force correction value Initialize the search step size constraint coefficient ; Utilizing preprocessed projection matrices Handling the initial gap The initial search direction components are obtained. ; For any k-th iteration, the contact gap Contact pairs that do not meet the non-intrusive condition are included in the set of iterative basic variables; Iteration step size for calculating the current contact force ,Will and The iteration step size corresponding to any i-th contact pair Updated to And remove the contact pairs from the set of iterative basic variables. It is the search direction component corresponding to the i-th contact pair in the (k-1)-th iteration. It is the contact force correction value corresponding to the i-th contact pair in the (k-1)-th iteration; Update the contact force correction value according to the updated iteration step size. Update contact gap , It is the contact force correction value for the (k-1)th iteration. It is the search direction component of the (k-1)th iteration. This is the contact gap in the (k-1)th iteration; the updated contact gap will be... Contact pairs that satisfy the non-intrusive condition are removed from the set of iterative basic variables; Calculate the iteration step size for the current search direction. And update the search direction component. And update the search step size constraint coefficient. , It is the search step size constraint coefficient for the (k-1)th iteration; Let k = k + 1 to enter the next iteration, until the set of iterative basic variables is empty, or the contact gap... and contact force correction value The preset iterative convergence conditions are met.
[0012] Its further technical solution is, when and ,or At that time, contact gap and contact force correction value The iterative convergence condition is met. It is the error threshold. It is the outward normal vector of the main facet. This indicates a search for the norm.
[0013] The further technical solution is that the motion state parameters include velocity and acceleration; The update yields the displacement at the (n+1)th computation time step. And update the speed at the (n+1)th computation time step. The acceleration at the (n+1)th calculation time step is updated. ; in, M is the mass matrix. It is a time increment. It is the spatial gradient of the contact gap. express transpose, It is the speed of the nth computation time step. It is the acceleration at the nth calculation time step.
[0014] The beneficial technical effects of this application are: This application discloses a non-intrusive transient contact numerical calculation method based on penalty functions and conjugate gradients. By modifying the penalty function contact algorithm through a pre-processed conjugate gradient method, the contact calculation strictly satisfies the non-intrusive condition, effectively avoiding problems such as divergence in calculation results and element distortion caused by non-physical penetration in the penalty function method. Moreover, compared to the traditional conjugate gradient contact algorithm, this application's method pre-applies partial contact force through the penalty function method, ensuring that the initial iteration of the conjugate gradient method is based on the modified calculation results. This significantly reduces the iterative convergence cost of the large deformation contact calculation process, giving it good convergence and engineering practicality. It has a wider range of applications, handling not only rigid collision problems but also high-speed collision problems with large plastic deformation. This method improves the simulation accuracy of transient collision problems, thereby further enhancing the prediction accuracy and stability of structural impact dynamics simulation. Attached Figure Description
[0015] Figure 1 This is a schematic diagram of the non-intrusive and intrusive states of the contact discrete model.
[0016] Figure 2 This is a flowchart of a non-invasive transient contact numerical calculation method.
[0017] Figure 3This is a flowchart of the preprocessing conjugate gradient method iterative solution process.
[0018] Figure 4 This is a comparison chart of calculation results for a thin-walled cylinder collision problem as an example.
[0019] Figure 5 This is a comparison chart of calculation results for a solid rod undergoing large deformation and collision. Detailed Implementation
[0020] The specific embodiments of this application will be further described below with reference to the accompanying drawings.
[0021] This application discloses a non-intrusive transient contact numerical calculation method based on penalty functions and conjugate gradients, used to simulate the transient collision process between any two colliding objects. Before performing simulation calculations, a contact discrete model is first constructed. The finite element method is used to discretize the first and second colliding objects into multiple elements, each element including multiple nodes. The surface element of the outer surface of any one of the first and second colliding objects is taken as the master surface, and the node of the outer surface of the other colliding object is taken as the slave node. The contact gap includes the contact gap of all contact pairs between the first and second colliding objects. Each contact pair includes a master surface and a slave node.
[0022] Please refer to Figure 1 The diagram shows the intrusion state in the contact region of the discrete contact model of the two colliding objects. Figure 1 (a) indicates a non-intrusive state, where the gap between the two colliding objects is greater than zero, and the contact gap is... Outer normal vector pointing to the main face Direction, no contact force between the two colliding objects, contact gap It is the contact gap between the main face and the slave node, which is the minimum distance from the slave node to the main face. Figure 1 (b) indicates the intrusion state, the contact gap between the two colliding objects. The outward normal vector pointing between the first and second colliders In opposite directions, and the contact force between the two colliding objects Outward normal vector direction.
[0023] Based on this contact discrete model, please refer to the numerical calculation method for non-invasive transient contact in this application. Figure 2 The flowchart shown illustrates that the method includes the following steps: Step 1: For any (n+1)th computation time step in the simulation process, there are two cases: When the (n+1)th calculation time step is the initial calculation time step (n=0), the initial state of the simulation calculation is initialized, including the mass matrix M, initial motion state parameters, initial position coordinates, and the time increment from the initial time to the first calculation time step. The mass matrix M includes the mass of all nodes of the first and second colliding objects, which is predetermined before the simulation calculation and remains unchanged throughout the entire calculation process. The time increment for each calculation time step can be preset to a fixed value or calculated based on structural deformation and material parameters. When the element deformation is large, the calculation time step will change accordingly. The motion state parameters include velocity and acceleration, which are updated in each calculation time step, and the position coordinates are also updated with the displacement changes in each calculation time step.
[0024] When the (n+1)th calculation time step is not the initial calculation time step (i.e., n≥1), the output of the previous calculation time step is directly used as the initial value of the current calculation time step, that is, the motion state parameters, position coordinates, displacement, etc. of the nth calculation time step are determined.
[0025] Step 2: Based on the motion state parameters of the first and second colliding objects at the nth calculation time step, estimate the predicted displacement values of the first and second colliding objects at the (n+1)th calculation time step. .
[0026] The method in this application requires calculating a preliminary step that does not consider contact before performing contact calculations at each computation time step. This step plays a crucial role in subsequent contact calculations. Specifically, for the velocity and acceleration in the nth computation time step, the calculation is performed on the time increment... Numerical integration yields the predicted displacement values of the first and second colliding objects. Time increment It is the time interval between the nth computation time step and the (n+1)th computation time step.
[0027] Step 3: Perform a contact search on the first and second colliders to identify all contact pairs between them and determine the contact gap between them at the (n+1)th calculation time step. And based on the contact gap and displacement prediction value Determine the predicted contact gap value at the (n+2)th computation time step. ,in, It is the spatial gradient of the contact gap and , express transpose, Indicates contact gap For position coordinates Differentiate, position coordinates The position coordinates of each node in the contact discrete model of the first and second colliders at the (n+1)th calculation time step are included. These position coordinates have been calculated based on the displacement of the previous calculation time step before the contact solution of the current calculation time step and belong to the known input of the current calculation time step.
[0028] The contact search method and the time difference scheme for differentiating the contact gap can be any of the existing technologies, such as bucket sort search or envelope box search. This application does not impose any restrictions, as long as it can be ensured that the calculated values do not affect the accuracy of the finite element calculation.
[0029] Step 4, based on the predicted contact gap value and displacement prediction value The penalty function method is used to calculate the contact force between the first and second colliding objects, and the estimated contact force is obtained. .
[0030] After obtaining the predicted displacement and contact gap values through the prediction step and contact search, further precise contact calculations are needed. First, a penalty function method is used to apply preliminary contact forces. This method, by applying a penalty contact force, allows the contact boundary to approximately satisfy the physical constraint of non-penetration, and is relatively simple and cost-effective. While applying the penalty contact force cannot guarantee the absence of non-physical penetration between structures, it significantly reduces the penetration depth from the node to the main surface after the prediction step. This will improve the convergence rate of subsequent conjugate gradient calculations, reducing convergence difficulty and computational cost.
[0031] Specifically, intrusion checks are performed on each contact pair individually, and in the case of boundary intrusion, i.e. At that time, a penalty contact force is applied, and the estimated value of the contact force is obtained. ,in, It is the outward normal vector of the main facet. It is the elastic force coefficient. It is the damping force coefficient. It is a penalty factor. This indicates taking the maximum value. The first term in the formula for calculating the estimated contact force. The penalty force directly caused by the amount of intrusion is similar to an elastic force; the second term The penalty force caused by the intrusion rate is similar to a damping force. , , The specific value can be customized according to the actual application requirements. Typically, =0.25, =2.0, =10 9 .
[0032] After applying the penalty contact force, the displacement and contact gap need to be recalculated based on the estimated contact force. and displacement prediction value The displacement at the (n+1)th calculation time step is updated. and the contact gap at the (n+2)th calculation time step .
[0033] The displacement and contact gap can be obtained from the governing equations of the contact solution step. Neglecting friction and damping, the governing equations of the contact solution step are constructed as follows: ,in, , It is the time increment. In the governing equations This can be viewed as the resistance to the inertial force generated by a unit displacement change. Essentially, this governing equation states that the inertial force caused by the displacement change due to the contact force should be balanced with the contact force.
[0034] Based on this governing equation, the displacement at the (n+1)th calculation time step can be solved and updated. ; and then update the contact gap at the (n+2)th calculation time step. , This represents the inverse of a matrix.
[0035] Step 5, when the updated contact gap If the non-intrusion condition is not met, use the updated displacement. and contact gap Using the non-intrusion condition as the initial value for iteration and the pre-processed conjugate gradient method as the constraint, the normal gap equation is solved iteratively to obtain the contact force correction value. And update the displacement at the (n+1)th calculation time step. The normal gap equation characterizes the zero-gap state of normal contact between two colliding objects, that is, the state in which the two colliding objects are in normal contact with each other without intrusion.
[0036] Due to the penalty factor in the contact force prediction formula of the penalty function method Elastic force coefficient Damping force coefficient These are all parameters that are set manually, so the contact force calculated by this method cannot be accurate. When the two colliding objects continue to make contact, it will be difficult to meet the non-invasive contact conditions, and further corrections are needed.
[0037] Considering that the conjugate gradient method can adapt to the strictly non-intrusive nature of explicit integration calculations, and that preprocessing can effectively improve the convergence speed of the conjugate gradient method, this application utilizes the preprocessed conjugate gradient method to modify the traditional penalty function contact algorithm to improve the computational accuracy of collision problems. When the updated contact gap... If the non-intrusion condition is not met, a second calculation using the pre-processed conjugate gradient method is required to update the displacements and contact forces that did not strictly meet the constraints in the initial contact calculation.
[0038] The expression for the normal gap equation is: , It is the normal gap factor, which is usually a small positive value.
[0039] The non-intrusive condition is a fundamental constraint ensuring the accuracy of conjugate gradient calculation. This application ensures that no contact intrusion occurs at the start of the next calculation time step by constraining the contact gap at time n+2, thus achieving true non-intrusive contact. This method achieves decoupled solution of internal and external forces and contact forces during the solution of the normal gap equation, and provides strict constraints to ensure non-intrusive contact solution. In one embodiment, the non-intrusive condition is:
[0040] in, It is the normal gap factor. It is the outward normal vector of the main face.
[0041] For details, please refer to Figure 3 The flowchart shown uses the pre-processed conjugate gradient method to iteratively solve the normal gap equation and obtain the contact force correction value. include: (1) Constructing the preprocessing projection matrix , Due to the symmetry of matrix A, the computation is transformed into solving the optimal solution problem satisfying the constraints using the conjugate gradient method. Furthermore, to accelerate convergence, a pre-projected conjugate gradient method is employed, projecting the equations onto the constraint space for solution. The pre-processed projection matrix is used to project the normal gap equations onto the constraint space for solution. The iterative basic variable set contains all contact pairs that do not satisfy the non-intrusion condition. Initialize integer parameter k=1, initialize contact gap Initialize contact force correction value Initialize the search step size constraint coefficient ; Utilizing preprocessed projection matrices Handling the initial gap The initial search direction components are obtained. ; (2) For any k-th iteration, the contact gap is... Contact pairs that do not meet the non-intrusive condition are included in the iterative basic variable set; that is, each contact pair is determined individually, and for any i-th contact pair, its contact gap is... In the case of base removal, the corresponding contact force correction value Otherwise, they will be admitted to the foundation.
[0042] (3) Calculate the iteration step size of the current contact force. ,Will and The iteration step size corresponding to any i-th contact pair Updated to And remove the contact pair from the set of iterative basic variables. It is the search direction component corresponding to the i-th contact pair in the (k-1)-th iteration. It is the contact force correction value corresponding to the i-th contact pair in the (k-1)-th iteration.
[0043] (4) Update the contact force correction value according to the updated iteration step size. Update contact gap , It is the contact force correction value for the (k-1)th iteration. It is the search direction component of the (k-1)th iteration. It is the contact gap in the (k-1)th iteration.
[0044] Further analysis is performed on each contact pair to update the contact gap. Contact pairs that satisfy the non-intrusive condition are removed from the set of iterative basic variables.
[0045] (5) Calculate the iteration step size for the current search direction. And update the search direction component. And update the search step size constraint coefficient. , It is the search step size constraint coefficient for the (k-1)th iteration; (6) Let k = k + 1 to enter the next iteration, until the set of iterative basic variables is empty, or the contact gap. and contact force correction value The preset iterative convergence condition is met. The preset iterative convergence condition is: when... and ,or At that time, contact gap and contact force correction value The iterative convergence condition is met. This is the error threshold, which can be customized to an acceptable error value, for example, by setting... =10-6, This indicates a search for the norm.
[0046] Step 6, based on the estimated contact force and contact force correction value Determine the contact force at the (n+1)th calculation time step. And use the displacement of the (n+1)th calculation time step Update the motion state parameters of the first and second colliders.
[0047] Contact force based on the (n+1)th calculation time step and displacement prediction value The displacement at the (n+1)th calculation time step is updated. Since the motion state parameters include velocity and acceleration, the velocity at the (n+1)th calculation time step is updated. The acceleration at the (n+1)th calculation time step is updated. ; It is the speed of the nth computation time step. It is the acceleration at the nth calculation time step.
[0048] Let n = n + 1, and proceed to the next calculation time step. Repeat the above contact force calculation process to update the motion state parameters for the next calculation time step until the maximum calculation time step or the total calculation time is reached. The maximum calculation time step and the total calculation time are predetermined according to the actual simulation requirements.
[0049] The effectiveness of the method in this application is further verified through simulation experiments using the following two examples: (1) Collision problem of thin-walled cylinder In this example, two identical thin-walled cylindrical shells collide with each other at the same velocity of 35 m / s. The shells are 0.46 m long, 0.2 m in diameter, and 0.005 m thick. Their elastic modulus is 25 GPa, Poisson's ratio is 0.3, and density is 7640 kg / m³. 3 The material is isotropic and linearly hardening plastic with an initial yield strength of 100 MPa and a hardening modulus of 230 MPa. The initial distance between the two cylinders is 0.01 m, and each thin cylinder is uniformly discretized into 2944 reduced integral shell elements. The problem is simulated using the method of this application and the penalty function method, respectively.
[0050] Figure 4 The calculation results of two contact algorithms were compared at 3ms and 6ms. Figure 4 (a) and (b) are the simulation results of the method of this application at 3ms and 6ms, respectively. Figure 4 (c) and (d) show the simulation results of the penalty function method at 3ms and 6ms, respectively. As can be seen from the figures, the method proposed in this application can simulate the contact collision process between shells relatively well. The simulation results of the penalty function method show non-physical penetration, and the penetration depth increases with time. This example further demonstrates the superiority of the method proposed in non-intrusive contact calculation.
[0051] (2) Large deformation and collision problem of solid rod In this example, a metal cylindrical rod impacts a target at a velocity of 400 m / s. The rod is 10 cm long and has a radius of 1 cm. The target has a radius of 4.5 cm and a thickness of 6 cm. The initial distance between the rod and the target is 0.5 mm. The material has an elastic modulus of 210 GPa, a Poisson's ratio of 0.3, and a density of 7800 kg / m³. 3 The rod material is isotropically power-hardening plastic with an initial yield strength of 249 MPa, a hardening modulus of 889 MPa, and a hardening index of 0.746. The target material density is 7800 kg / m³. 3 The elastic modulus is 210 GPa, and Poisson's ratio is 0.3. The rod and target are uniformly discretized into 38,400 and 53,760 reduced integral solid elements, respectively. The problem is simulated using the method of this application and the preprocessed conjugate gradient method.
[0052] like Figure 5 Simulation results were compared using two methods at the instant of contact and a period of time after the contact collision. The instant of contact was 0.015 ms. Figure 5 (a) and (b) are the simulation results of the method of this application at 0.015ms and 0.025ms. Figure 5 (c) and (d) show the simulation results of the simple conjugate gradient method at 0.015 ms and 0.025 ms, respectively. As can be seen from the figures, the simple conjugate gradient method struggles to converge at the contact instant to obtain an accurate solution to this strongly nonlinear problem, exhibiting element distortion and intrusion in the collision region. As shown in Table 1, it consumes significantly more computation time in this process compared to the method in this application, mainly due to the large number of iterations. Once a step fails to converge, subsequent steps become even more difficult to converge, resulting in a much longer computation time and the inability to obtain correct results. This example demonstrates that the method in this application has a lower convergence cost than the simple conjugate gradient method and has advantages in large deformation transient contact simulation.
[0053] Table 1 Comparison of calculation time for collision problems involving large deformation of solid rods
[0054] The above descriptions are merely preferred embodiments of this application, and this application is not limited to the above embodiments. It is understood that other improvements and variations that can be directly derived or conceived by those skilled in the art without departing from the spirit and concept of this application should be considered to be included within the protection scope of this application.
Claims
1. A non-intrusive transient contact numerical calculation method based on penalty functions and conjugate gradients, characterized in that, The non-invasive transient contact numerical calculation method includes: For any (n+1)th calculation time step, based on the motion state parameters of the first and second colliding objects at the nth calculation time step, the predicted displacement values of the first and second colliding objects at the (n+1)th calculation time step are estimated. ; Determine the contact gap between the first and second colliders at the (n+1)th computation time step. And based on the contact gap and the predicted displacement value Determine the predicted contact gap value at the (n+2)th computation time step. ; Based on the predicted contact gap value and the predicted displacement value The penalty function method is used to calculate the contact force between the first and second colliding objects, and the estimated contact force is obtained. Based on the estimated contact force and the predicted displacement value The displacement at the (n+1)th calculation time step is updated. and the contact gap at the (n+2)th calculation time step ; When the updated contact gap If the non-intrusion condition is not met, use the updated displacement. and contact gap Using the non-intrusion condition as the initial value for iteration and the pre-processed conjugate gradient method as the constraint, the normal gap equation is solved iteratively to obtain the contact force correction value. And update the displacement at the (n+1)th calculation time step. The normal gap equation characterizes the zero-gap state of normal contact between two colliding objects. Based on the estimated contact force and the contact force correction value Determine the contact force at the (n+1)th calculation time step. And use the displacement of the (n+1)th calculation time step Update the motion state parameters of the first and second colliders, and proceed to the next calculation time step.
2. The non-invasive transient contact numerical calculation method according to claim 1, characterized in that, The motion state parameters include velocity and acceleration, and the velocity and acceleration at the nth calculation time step are incremented by the time increment. Numerical integration yields the predicted displacement values of the first and second colliding objects. Time increment It is the time interval between the nth computation time step and the (n+1)th computation time step.
3. The non-invasive transient contact numerical calculation method according to claim 2, characterized in that, Determine the predicted contact gap value at the (n+2)th computation time step. ; The penalty function method is used to calculate the contact force between the first and second colliding objects, and the estimated contact force is obtained. ,in, It is the outward normal vector of the main facet. It is the elastic force coefficient. It is the damping force coefficient. It is a penalty factor. It is the spatial gradient of the contact gap and , express transpose, Indicates contact gap For position coordinates Perform differentiation. This indicates taking the maximum value.
4. The non-invasive transient contact numerical calculation method according to claim 2, characterized in that, The update yields the displacement at the (n+1)th computation time step. ; The contact gap is updated to the (n+2)th calculation time step. ; in, , It is a time increment. Represents the inverse of a matrix; It is the spatial gradient of the contact gap and , express transpose, Indicates contact gap For position coordinates Perform differentiation.
5. The non-invasive transient contact numerical calculation method according to claim 1, characterized in that, The finite element method is used to discretize the first and second colliding objects into multiple elements, each element including multiple nodes. The surface element of the outer surface of any one of the first and second colliding objects is taken as the principal surface element, and the nodes on the outer surface of the other colliding object are taken as slave nodes. The contact gap includes the contact gap of all contact pairs between the first and second colliding objects, and each contact pair includes one principal surface element and one slave node. The normal gap equation is: The non-invasive condition is: in, , It is a time increment. Represents the inverse of a matrix; It is the spatial gradient of the contact gap and , express transpose, Indicates contact gap For position coordinates Perform differentiation; It is the normal gap factor. It is the outward normal vector of the main face.
6. The non-invasive transient contact numerical calculation method according to claim 5, characterized in that, The pre-processed conjugate gradient method is used to iteratively solve the normal gap equation to obtain the contact force correction value. include: Constructing the preprocessing projection matrix , , , It is the spatial gradient of the contact gap. express The transpose of M, where M is the mass matrix. It is a time increment. Let A be the diagonal matrix of matrix A. The inverse of the matrix is represented; the preprocessed projection matrix is used to project the normal gap equation onto the constraint space for solution. Initialize integer parameter k=1, initialize contact gap Initialize contact force correction value Initialize the search step size constraint coefficient ; using the preprocessed projection matrix Handling the initial gap The initial search direction components are obtained. ; For any k-th iteration, the contact gap Contact pairs that do not meet the non-intrusive condition are included in the set of iterative basic variables; Iteration step size for calculating the current contact force ,Will and The iteration step size corresponding to any i-th contact pair Updated to And remove the contact pair from the set of iterative basic variables. It is the search direction component corresponding to the i-th contact pair in the (k-1)-th iteration. It is the contact force correction value corresponding to the i-th contact pair in the (k-1)-th iteration; Update the contact force correction value according to the updated iteration step size. Update contact gap , It is the contact force correction value for the (k-1)th iteration. It is the search direction component of the (k-1)th iteration. This is the contact gap in the (k-1)th iteration; the updated contact gap will be... Contact pairs that satisfy the non-intrusive condition are removed from the set of iterative basic variables; Calculate the iteration step size for the current search direction. And update the search direction component. And update the search step size constraint coefficient. , It is the search step size constraint coefficient for the (k-1)th iteration; Let k = k + 1 to enter the next iteration, until the set of iterative basic variables is empty, or, the contact gap. and contact force correction value The preset iterative convergence conditions are met.
7. The non-invasive transient contact numerical calculation method according to claim 6, characterized in that, when and ,or At that time, contact gap and contact force correction value The iterative convergence condition is met. It is the error threshold. It is the outward normal vector of the main facet. This indicates a search for the norm.
8. The non-invasive transient contact numerical calculation method according to claim 1, characterized in that, Motion state parameters include velocity and acceleration; The update yields the displacement at the (n+1)th computation time step. And update the speed at the (n+1)th computation time step. The acceleration at the (n+1)th calculation time step is updated. ; in, M is the mass matrix. It is a time increment. It is the spatial gradient of the contact gap. express transpose, It is the speed of the nth computation time step. It is the acceleration at the nth calculation time step.