Material Point Method Based on Improved Contact Algorithm for Bicolor Coin Embossing Forming Simulation

Abstract

The present invention provides a material point method based on an improved contact algorithm for bicolor coin embossing forming simulation. The method includes the following steps: establishing a dynamic physical model for the bicolor coin embossing forming process; performing material point integration on the dynamic model of the embossing forming process; discretizing the forming die with finite element triangular / quadrilateral elements; establishing a background grid according to the blank data model, discretizing the bicolor coin blank with tetrahedral / hexahedral elements, and establishing a set of mass points for the copper and aluminum material regions; calculating the strain increment and spin rate increment of the mass points, solving the stress and elastoplastic related variables of the mass points according to the constitutive model, and updating the density of all mass points. The present invention avoids the problem of material self-contact judgment in the finite element simulation of bicolor coin embossing forming through the material point method.

Description

Technical Field

[0001] The present invention relates to the technical field of stamping, and particularly to a material point method based on an improved contact algorithm applied to the simulation of two-color coin stamping forming Background Art

[0002] The production of two-color metal coins is a very important technology for coin anti-counterfeiting today, which is conducive to improving the anti-counterfeiting performance of residents' use of metal coins. At the same time, the beautifully designed two-color coins have strong ornamental value and great collection value. With the continuous innovation of science and technology, higher requirements are put forward for the manufacture of metal coins, especially for the production of two-color coins with obvious advantages and high technical content. During the actual stamping process of two-color coins, the mutual contact and extrusion between two different materials are likely to cause stress concentration problems, which in turn lead to the generation of microcracks or even cracks, seriously affecting the product quality. In the actual production process, continuous trial molds are required to overcome the stress concentration problem. This will greatly increase the production cycle and production cost. Analyzing the stamping process using a simulation algorithm can assist engineers in providing a theoretical basis for solving the stress concentration problem, and at the same time effectively reducing the production cycle and production cost. The finite element method has been widely applied and recognized in many engineering fields, but there are also mesh distortion problems in the large deformation problems of stamping forming, resulting in the termination of calculations or the deviation of simulation results from reality. At the same time, for the stamping simulation analysis of two-color coin materials, the finite element method needs to additionally judge the contact penetration problem at the interface of different materials, greatly increasing the simulation time and difficulty. The material point method is a numerical calculation method that uses the dual description of Lagrangian particles and Eulerian grids and is very suitable for analyzing problems involving large deformations and contact problems. It makes up for the problems of Lagrangian method in mesh division, such as mesh distortion leading to a decrease in the accuracy of element interpolation and integration and non-convergence of results, and has obvious advantages in problems involving large deformations such as contact and impact

[0003] Based on the obvious advantages of the material point method in the field of large deformation numerical simulation, this research invention intends to use the material point method to simulate and analyze the stamping forming of two-color coins. However, the traditional material point method has the defect of virtual contact, resulting in non-conservation of momentum and interface defects. The present invention proposes a material point method with an improved contact algorithm, which avoids virtual contact and improves the accuracy of contact judgment of the material point method. In addition, the material point method is used to avoid the problem of self-contact judgment of materials in the finite element simulation of two-color coin stamping forming, greatly reducing the difficulty and calculation amount of the contact algorithm in the simulation Summary of the Invention

[0004] The purpose of the present invention is to provide a material point method based on an improved contact algorithm applied to the simulation of two-color coin stamping forming for the problem of virtual contact caused by point-to-point contact in the material point method. The purpose of the present invention is achieved as follows: The method includes the following steps:

[0005] Step 1: Establish a dynamic physical model for the stamping forming process of the two-color coin;

[0006] Step 2: Perform material point integration on the dynamic model of the stamping forming process;

[0007] Step 3: Discretize the forming die with finite element triangular / quadrilateral elements;

[0008] Step 4: Establish a background grid according to the blank data model, discretize the two-color coin blank with tetrahedral / hexahedral elements, and establish a set of mass points for the copper and aluminum material regions;

[0009] Step 5: Calculate the node loads of the background grid;

[0010] Step 6: Map the physical quantities carried by the mass points to the background grid nodes through interpolation shape functions;

[0011] Step 7: Solve the system adaptive time step according to the different material properties of the two-color coin;

[0012] Step 8: Solve the momentum equation and apply boundary conditions;

[0013] Step 9: Update the physical quantities of the mass points;

[0014] Step 10: Update the physical quantities of the background grid points;

[0015] Step 11: Update the spatial position of the die;

[0016] Step 12: Judge the penetration state of the mass points and solve the penetration distance;

[0017] Step 13: Solve the contact force;

[0018] Step 14: Update the physical quantities of the mass points again;

[0019] Step 15: Apply boundary conditions;

[0020] Step 16: Calculate the strain increment and spin rate increment of the mass points, solve the stress and elastoplastic related variables of the mass points according to the constitutive model, and update the density of all mass points.

[0021] The present invention also includes the following features:

[0022] 1. Specifically, Step 1 is as follows:

[0023] From a mechanical perspective, the embossing forming process is a quasi-static low-speed motion process, which is generally solved by implicit integration methods. However, considering the complex surface pattern of the commemorative coin, millions of grids are required to describe it, resulting in a large amount of memory consumption for solving the linear equations. At the same time, the embossing process is accompanied by complex material nonlinearity. Using the Newton-Raphson iteration to solve may lead to non-convergence. Therefore, this study intends to use the dynamic explicit algorithm to describe and solve the embossing forming process. Let the blank be the object under study. When it undergoes plastic deformation, it should satisfy the basic equilibrium equation:

[0024]

[0025] where σ ij and b i are the Cauchy stress tensor and body force respectively, ρ is the material density, γ is the damping coefficient, and are the velocity and acceleration of any point in the object respectively, and i and j are the coordinate directions. According to the principle of virtual displacement, relevant stress boundary conditions, and Gauss divergence theorem, the virtual power equation of the system can be derived as:

[0026]

[0027] where V is the object region, is the virtual velocity, is the virtual strain rate corresponding to the Cauchy stress σ ij . Among them, S p is the surface of the known external force p i ; S c is the surface in contact with another object, and the contact surface force is denoted as q i ; the power terms of inertia force and damping force reflect the inertia effect and physical damping effect of the object system.

[0028] 2. Step 2 is specifically as follows:

[0029] In the material point method, each material point occupies a certain volume space and carries a certain mass and other related physical quantities. The material density at any spatial point x can be expressed as:

[0030]

[0031] where n p is the number of material points, m p is the mass of the material point, δ() is the Dirac function, and x p is the coordinate of the material point p. The displacement u ip of any material point p in the background grid can be obtained by interpolating the shape function of the background grid nodes:

[0032]

[0033] where N Ip is the shape function of node I related to the particle p, and u iI is the velocity in the i-th direction of the background grid node I. Substitute the above density and displacement formulas into the system virtual power equation in step 1, and at the same time consider that there is no body force b i and external force p i . We can get:

[0034]

[0035] where is the mass of node I, and are the velocity and acceleration in the i-th direction of node I respectively. is the external force in the i-th direction of node I, is the internal force in the i-th direction of node I.

[0036] 3. Step 4 is specifically as follows:

[0037] First, divide the region occupied by the blank into regular hexahedron background grids (the number of grids in the x, y, and z directions are Gx, Gy, and Gz respectively), and find the integration points of each background grid. Then, perform tetrahedral mesh division on the regions of the two different materials, copper and aluminum. As Figure 1 shown, it is a simple 2D example (the background grid is a rectangle, and the blank grid is a triangle.. The actual 3D blank size is as Figure 2 shown). Each small rectangle represents a background grid cell, the solid circle represents the grid node, the small solid triangle represents the copper particle, and the small solid square represents the aluminum particle. Next, find the set I Cu of background grid cells of each tetrahedral element in the copper material region. The specific method to find this set is as follows: Number the background grid cells from 1 in the x, y, and z directions to G = Gx * Gy * Gz. According to the four node coordinates of the tetrahedral element, a spatial cuboid can be found to just enclose the tetrahedron. Assume the coordinates of the left-front point (x min , y min , z min ) and the right-back point (x max , y max , z max ) of the background grid. Based on the background grid number p of any point x = (x, y, z) in space, find the set I Cu of all background grids contained in this cuboid. Define the set composed of the Gaussian points of all cells in this set as point set II CuAccording to the surface boundary elements of the copper material (automatically generated by the mesh generation system and programmed to read these surface elements), the Gauss points outside the copper material in the point set II are removed, and the set formed by the remaining points is the mass point set III of the copper material domain. Cu As shown by the solid small triangle in Figure 1 . Similarly, for the aluminum material region, similar steps are carried out to obtain the mass point set III of the aluminum material. cu As Figure 1 shown by the solid small triangle. Similarly, for the aluminum material region, similar steps are carried out to obtain the mass point set III of the aluminum material. Al As Figure 1 shown by the solid small square. Finally, the two sets III Cu and III Al are merged into the mass point set III of the blank cake, and corresponding material properties such as density, volume, Poisson's ratio, yield stress, hardening index, strengthening coefficient and other parameters are assigned to each mass point. In order to accurately describe the initial contour of the coin, the Gauss point sets of each triangular element on the surface of the blank cake are added to the mass point set of the blank cake to form a new mass point set of the blank cake. As Figure 3 shown, a quarter model is cut off to facilitate viewing of the internal mass points. The mass point set on the surface of the blank cake is divided into three parts: the upper surface, the side surface and the lower surface mass point sets, which form contact pairs with the upper die, the middle ring and the lower die respectively. During the simulation process, contact judgment processing is carried out for these three pairs of contact pairs respectively.

[0038] 4. Step 7 is specifically as follows:

[0039] Under normal circumstances, for the single-material embossing forming simulation problem, we use the central difference dynamic explicit algorithm for solution. This algorithm is conditionally stable, and the system time step Δt must be less than the critical step Δt cr That is

[0040] Δt = αΔt cr

[0041] In the formula, α is a constant, and its value range is: 0.8 ≤ α ≤ 0.98

[0042] Due to the CFL (Courant Friedrichs Lewy) condition, in order to ensure the stability of the explicit time integration, when using the finite element method for dynamic explicit analysis, the distance that the wave propagates within one time step must not exceed one element. And the sound speed of the elastic material is:

[0043]

[0044] In the formula, E is the elastic modulus, v is Poisson's ratio, and ρ is the material density. Assuming the nominal length of the element is l e Then the critical time step is: As the deformation increases, the nominal element size gradually decreases, and the critical time step also gradually decreases. When the mesh undergoes severe distortion, the minimum nominal length of the solution approaches 0, resulting in the time step also approaching 0, making the simulation analysis impossible to carry out.

[0045] The material point method is different from the finite element method, and all time integrations need to be carried out on a fixed background mesh. When determining the critical time step, not only the sound speed of the material needs to be considered, but also the velocity of the material points needs to be considered. In particular, in the problem of hypervelocity impact, the velocity of the material points is relatively large, on the same order of magnitude as the sound speed of the material or even exceeds the sound speed of the material. Therefore, the influence of the velocity of the material points on the simulation calculation cannot be ignored. In the process of embossing simulation, the material point method based on the dynamic explicit central difference algorithm is used for solution, and the critical time step needs to be calculated at each step to make the numerical solution stable:

[0046]

[0047] In the formula, c is the sound speed of the material point, and u is the velocity of the material point. When using a uniform background mesh, l g can also be taken as d c , that is, the background mesh size. In the present invention, the blank has two different materials, namely copper inside and aluminum outside. The size of the time step is determined by the properties of the system and the material parameters of the model. Since two different materials are involved, the time steps in the two material regions need to be calculated separately and the smaller value is taken. The time step in the copper material is:

[0048]

[0049] In the formula, and are the sound speed in the material point set III Cu and the velocity of any material point p1, respectively.

[0050] Similarly, the time step in the aluminum material is:

[0051]

[0052] In the formula, and are the sound speed in the material point III Al and the velocity of any material point p2, respectively.

[0053] We need to calculate the wave speeds of the two materials and the velocities of the material points of the two materials during the embossing process, and then compare the critical time steps of the two, and finally determine that the critical time step of the system is the minimum of the two. Therefore, the critical time step of the system is:

[0054]

[0055] In the formula, The time step for Cu The time step for Al

[0056] 5. Step 12 is to judge the penetration state of material points and solve the penetration distance, which mainly includes the following content: As Figure 4 shown, according to the traditional material point method theory, both the mold particle j′ and the blank cake particle j contribute to the velocity of the background grid node I (because the background grids they are located in share the node I), so it is determined that the mold is in contact with the blank cake particle. In fact, they are not in contact, thus resulting in false contact. The present invention uses a point - surface contact algorithm to replace the point - point contact algorithm of traditional material points to avoid virtual contact and improve the accuracy of simulation calculation. When judging the contact state between each particle on the surface of the blank cake and the upper die, lower die, and middle ring, and calculating the contact values (contact force, penetration depth, and friction force), it will cause a large amount of calculation in each time step. In order to reduce the calculation time of contact judgment, a strategy combining global search and local search is adopted. The global search is to find the candidate mold elements in contact with each particle, and the process is as follows.

[0057] Construct a block containing specific mold grids: The size of the block is defined according to the maximum and minimum spatial coordinates of each mold. The coordinates of the lower left corner and the lower right corner of the block are (x min , y min , z min ) and (x max , y max , z max ) respectively. Then divide the block area into sub - cells, and the number of cells in each direction is N x , N y , N z respectively, and number each cell along the x, y, z three directions. Note that this block moves with the mold grid. In order to distinguish this kind of cell from the background cell grid, we name this kind of cell the spatial cell.

[0058] Search for the candidate mold grid elements contained in each cell. As Figure 5 shown is the cell schematic diagram of a certain two - dimensional object. There are a total of N = N x ×N y = 8×5 = 40 spatial cells, and only the cell numbers at the four corners are marked in the figure. The cell numbered 21 contains five mold grid elements, and its nodes are represented by solid circles. Finding the candidate elements contained in each cell generally includes two steps. First, expand a certain mold spatial cell domain along three directions, and the expansion distance is half of the sub - spatial cell size, and loop through each mold element to obtain the cells containing it, and finally obtain the candidate mold elements contained in each cell. For a given particle x p=(x, y, z), and the unit number it locates is

[0059] where [] is a rounding operation. For the particle x p , the subspace lattice containing the particle can be located. According to the above global search results, the candidate die unit set of the particle can be obtained. Next, local search is performed in this candidate set and contact judgment is carried out. In order to judge the contact state between meshes and the particle x p , we only need to find whether the particle contacts the candidate unit or penetrates it. The following steps can be executed to achieve this goal.

[0060] Find potential contact units: As Figure 6 shown, for the particle x p , there are four candidate units S1, S2, S3, S4. Find the node m p nearest to x s , which does not coincide with x p . If the following inequalities are satisfied, the particle will contact the S i unit

[0061] (C i ×g2)·(C i ×C i+1 )>0, (C i ×g2)·(g2×C i+1 )>0

[0062] Here, i = 1, 2, 3, 4 is the i-th die unit to be judged. C i and C i+1 are two unit vectors with a common starting point, along two adjacent sides of the S i unit, g1 starts from m s and ends at x p . g2 is the projection of the unit S i on the vector g1, written as:

[0063] g2 = g1 - (g1·m)m

[0064] Here

[0065]

[0066] As Figure 7 shown, in order to calculate the spatial coordinates of x c , a parametric formula is established on the unit:

[0067] r(ξ, η) = f1(ξ, η)i1 + f2(ξ, η)i2 + f3(ξ, η)i3

[0068] Among them N j (ξ, η) is a bilinear function of the shell element, is the i-th coordinate of node j, (ξ, η) is the local coordinate of any point within the element, and i1, i2, i3 are the unit vectors in the spatial coordinates respectively. t is the vector starting from the original coordinate and ending at x p ending vector. Because the vector is perpendicular to the contact element, the local coordinates (ξ c , η c ) of the contact point x must satisfy:

[0069]

[0070]

[0071] Judgment of the contact state and calculation of the penetration depth d:

[0072]

[0073] where n c is the unit outward normal vector of the contact element at the point (ξ c , η c ), and the calculation method is: Description of the drawings

[0074] Figure 1 is a schematic diagram of the discretization of material points

[0075] Figure 2 is a schematic diagram of the cross-section of a two-color coin

[0076] Figure 3 is a schematic diagram of potential contact pairs of mass points

[0077] Figure 4 is a schematic diagram of virtual contact of material points

[0078] Figure 5 is a schematic diagram of the division of the background grid cell area

[0079] Figure 6 is a schematic diagram for judging the contact state between mass points

[0080] Figure 7 is a schematic diagram for calculating the contact position of mass points

[0081] Figure 8 is a schematic diagram for calculating the contact force of mass points Specific implementation manners

[0082] In order to clarify the purpose, technical solution and advantages of the present invention, the present invention will be described in more detail below in conjunction with the attached structures and embodiments. The present invention is a material point method based on an improved contact algorithm applied to the simulation of the embossing forming of two-color coins, including the following steps:

[0083] Step 1: Establish a dynamic physical model for the embossing forming process of two-color coins;

[0084] From a mechanical perspective, the embossing forming process is a quasi-static low-speed motion process. Generally, an implicit integration method is used to solve it. However, considering the complex surface pattern of commemorative coins, millions of meshes are required to describe it, resulting in a large amount of memory occupied by the solution of the linear equations. At the same time, the embossing process is accompanied by complex material nonlinearity. Using the Newton-Raphson iteration to solve may lead to non-convergence. Therefore, this study intends to use the dynamic explicit algorithm to describe and solve the embossing forming process. Let the blank be the object under study, and when it undergoes plastic deformation, it should satisfy the basic equilibrium equation:

[0085]

[0086] where σ ij，j are the Cauchy stress tensors respectively, b i is the body force, ρ is the material density, γ is the damping coefficient, and are the velocity and acceleration of any point in the object respectively, and i and j are the coordinate directions respectively. According to the principle of virtual displacement, relevant stress boundary conditions and Gauss divergence theorem, the virtual power equation of the system can be deduced as:

[0087] In the formula, V is the object region, is the virtual velocity, is the virtual strain rate corresponding to the Cauchy stress σ ij .

[0088] where S p is the surface of the known external force p i ; S c is the surface in contact with another object, and the contact surface force is denoted as q i ; The power terms of inertial force and damping force reflect the inertial effect and physical damping effect of the object system.

[0089] Step 2: Perform material point integration on the dynamic model of the embossing forming process;

[0090] In the material point method, each material point occupies a certain volume space and carries a certain mass and other related physical quantities. The material density at any spatial point x can be expressed as:

[0091]

[0092] where n p is the number of mass points, m p is the mass of a mass point, δ() is the Dirac function, and x p is the coordinate of mass point p. The displacement u ip of any mass point p in the background grid can be obtained by interpolating with the shape functions of the background grid nodes:

[0093]

[0094] where, N Ip is the shape function of node I related to mass point p, and u iI is the velocity of the background grid node I in the i-th direction. Substitute the above density and displacement formulas into the system virtual power equation in Step 1, and at the same time consider that there is no body force b i and external force p i during the embossing process, we can get:

[0095]

[0096] where, is the mass of node I, and are the velocity and acceleration of node I in the i-th direction respectively. is the external force on node I in the i-th direction, is the internal force on node I in the i-th direction.

[0097] Step 3: Discretize the forming die with finite element triangular / quadrilateral elements;

[0098] Use open-source software to perform tetrahedral mesh generation on the two-color coin and obtain the surface triangular mesh discretization, and establish a background grid that only contains the material region of the two-color coin.

[0099] Step 4: Establish a background grid according to the blank data model, and perform tetrahedral / hexahedral element discretization on the two-color coin blank, and establish a set of mass points for the copper and aluminum material regions;

[0100] After the blank is discretized into mass points, its mass point density can be expressed as:

[0101]

[0102] where n p is the number of mass points, m p is the mass of a mass point, δ() is the Dirac function, and x p is the coordinate of mass point p. In this paper, the Material Point Method (MPM) discretizes the material domain, i.e., the initial two-color coin blank model, into a set of particles, such as Figure 1As shown below. First, divide the hexahedron background grid (the number of grids in the x, y, and z directions are Gx, Gy, and Gz respectively) for the area occupied by the blank cake, and find the integration points of each background grid. Then, perform tetrahedral mesh division on the areas of the two different materials, copper and aluminum. As Figure 1 shown, it is a simple 2D example (the background grid is rectangular, and the blank grid is triangular. The actual 3D blank cake size is as Figure 2 shown). Each small rectangle represents a background grid unit, the solid circle represents the grid node, the small solid triangle represents the copper material point, and the small solid square represents the aluminum material point. Next, find the set I of background grid units for each tetrahedral unit in the copper material area Cu . The specific method to find this set is as follows: Number the background grid units from 1 in the x, y, and z directions to G = Gx * Gy * Gz. According to the four node coordinates of the tetrahedral unit, a spatial cuboid can be found to just enclose the tetrahedron. Assume the coordinates of the left-front most point of the background grid (x min , y min , z min ) and the coordinates of the right-back most point (x max , y max , z max ). Based on the background grid number p of any point x = (x, y, z) in space, find the set I of all background grids contained in this cuboid Cu . Define the set composed of the Gauss points of all units in this set as point set II Cu . According to the surface boundary units of the copper material (automatically generated by the mesh generation system and read by the program to obtain these surface units), remove the Gauss points outside the copper material from point set II Cu . The remaining set of points forms the set of material points III of the copper material domain cu , as shown by the solid small triangles in Figure 1 . Similarly, perform similar steps on the aluminum material area to obtain the set of material points III of the aluminum material Al , as shown by the solid small squares in Figure 1 . Finally, merge the two sets III Cu and III Al into the set of material points III of the blank cake, and assign corresponding material properties such as density, volume, Poisson's ratio, yield stress, hardening index, strengthening coefficient, etc. to each material point. To accurately describe the initial contour of the coin, add the Gauss point set of each triangular unit on the surface of the blank cake to the set of material points of the blank cake to form a new set of material points of the blank cake. As Figure 3As shown, one - quarter of the model is removed for easy viewing of internal mass points. The mass point set on the blank surface is divided into three parts: the upper - surface, side - surface, and lower - surface mass point sets, which form contact pairs with the upper die, middle ring, and lower die respectively. During the simulation process, contact judgment processing is performed on these three pairs of contact pairs respectively.

[0103] Step 5: Calculate the load of the background grid nodes;

[0104] The internal force vector of the node is:

[0105]

[0106] In the formula, is the i - th component of the internal force of node I, m p is the mass of mass point p, σ ij is the stress of mass point p, is the j - th partial derivative of the shape function of node I at mass point p, is the position of mass point p at time n, ρ p is the density of mass point p.

[0107] The external force vector of the node is:

[0108]

[0109] In the formula, is the external force of the node, b i is the body force, and h is the thickness of the boundary layer. For the embossing forming problem, without considering contact, this external force of the node is 0.

[0110] The resultant force at node I at time n is:

[0111]

[0112] In the formula, is the resultant force at node I at time n.

[0113] Step 6: Map the physical quantities carried by the mass points to the background grid nodes through the interpolation shape function;

[0114] The mass of node I is:

[0115]

[0116] The momentum and velocity of the node are:

[0117]

[0118]

[0119] In the formula, is the i - th component of the shape function of node I, is the i-th component of the velocity of the p-th particle at the moment of n - 1 / 2, is the i-th component of the momentum of node I at the moment of n + 1 / 2, is the i-th component of the velocity of node I at the moment of n + 1 / 2, represents the mass of node I at the moment of n.

[0120] Step 7: Solve the system adaptive time step according to the different material properties of the two-color coin;

[0121] Under normal circumstances, for the stamping forming simulation problem of a single material, we use the central difference dynamic explicit algorithm to solve it. This algorithm is conditionally stable, and the system time step Δt must be less than the critical step length Δt cr , that is

[0122] Δt = αΔt cr (13)

[0123] In the formula, α is a constant, and its value range is: 0.8 ≤ α ≤ 0.98

[0124] Due to the CFL (Courant Friedrichs Lewy) condition, in order to ensure the stability of explicit time integration, when using the finite element method for dynamic explicit analysis, the distance that the wave travels within a time step must not exceed one element. And the sound speed of the elastic material is:

[0125]

[0126] In the formula, E is the elastic modulus, v is the Poisson's ratio, and ρ is the material density. Assuming that the nominal length of the element is l e , then the critical time step is: As the deformation increases, the nominal size of the element gradually decreases, and the critical time step also gradually decreases. When the mesh is severely distorted, the minimum nominal length of the solution approaches 0, resulting in the time step also approaching 0, so that the simulation analysis cannot be carried out.

[0127] The material point method is different from the finite element method. All time integrations need to be carried out on a fixed background grid. When determining the critical time step, not only the sound speed of the material needs to be considered, but also the velocity of the particles needs to be considered. In particular, in the problem of hypervelocity impact, the particle velocity is relatively large, of the same order of magnitude as the material sound speed or even exceeds the material sound speed. Therefore, the influence of the particle velocity on the simulation calculation cannot be ignored. In the process of stamping forming simulation, based on the material point method solution of the dynamic explicit central difference algorithm, the critical time step needs to be calculated at each step to make the numerical solution stable:

[0128]

[0129] In the formula, c is the sound speed of the particle, and u is the particle velocity. When using a uniform background grid, l g can also be taken as d c , that is, the background grid size. In the present invention, the blank cake has two different materials, namely copper inside and aluminum outside. The size of the time step is determined by the properties of the system and the model material parameters. Since two different materials are involved, it is necessary to calculate the time steps in the two material regions separately and take the smaller value. The time step in the copper material is:

[0130]

[0131] In the formula, and are the sound speed and the particle velocity of any particle p1 in particle set III Cu respectively.

[0132] Similarly, the time step in the aluminum material is:

[0133]

[0134] In the formula, and are the sound speed and the particle velocity of any particle p2 in particle III Al respectively.

[0135] We need to calculate the wave speeds of the two materials and the velocities of the material points of the two materials during the imprinting process, and then compare the critical time steps of the two, and finally determine that the critical time step of the system is the minimum value of the two. Therefore, the critical time step of the system is:

[0136]

[0137] In the formula, is the time step of Cu, is the time step of Al.

[0138] Step 8: Solve the momentum equation and apply boundary conditions;

[0139] Integrate the momentum equation of the background grid nodes:

[0140]

[0141] In the formula, and are the momenta in the i-th direction of node I at n + 1 / 2 and n - 1 / 2 moments respectively.

[0142] Step 9: Update the physical quantities of the particles;

[0143] Update the particle velocity:

[0144]

[0145] It represents the i-th component of the velocity of the p material point at the n + 1 / 2 moment.

[0146] Update the position of the material point:

[0147]

[0148] In the formula, represents the i-th component of the coordinate of the p material point at the n + 1 moment, represents the i-th component of the coordinate of the p material point at the n moment.

[0149] Step 10: Update the physical quantities of the background grid points;

[0150] Update the velocity of the background grid nodes

[0151]

[0152] In the formula, represents the velocity of node I at the n + 1 / 2 moment.

[0153] Step 11: Update the spatial position of the mold;

[0154] Update the position of the mold:

[0155] d punch = d punch + Δd punch (23)

[0156] d punch represents the current total displacement of the upper mold, and Δd punch represents the displacement of the upper mold at the current moment.

[0157] Step 12: Judge the penetration state of the material point and solve the penetration distance;

[0158] Such as Figure 4As shown, according to the traditional material point method theory, both the die particle j' and the blank particle j contribute to the velocity of the background grid node I (because the background grids they are in share the node I), so it is determined that the die is in contact with the blank particle. In fact, they are not in contact, resulting in false contact. The present invention uses a point-plane contact algorithm to replace the point-point contact algorithm of the traditional material point to avoid virtual contact and improve the accuracy of simulation calculation. When judging the contact state between each particle on the blank surface and the upper die, lower die, and middle ring, and calculating the contact values (contact force, penetration depth, and friction force), a large amount of computational work is caused in each time step. To reduce the computational time of contact judgment, a strategy combining global search and local search is adopted. The global search is to find the candidate die elements in contact with each particle, and the process is as follows.

[0159] Construct a block containing specific die grids: The size of the block is defined according to the maximum and minimum spatial coordinates of each die. The coordinates of the lower left and lower right corners of the block are (x min , y min , z min ) and (x max , y max , z max ). Then divide the block area into sub-cells, and the number of cells in each direction is N x , N y , N z , and number each cell along the x, y, and z directions respectively. Note that this block moves with the die grid. To distinguish this kind of cell from the background cell grid, we name this kind of cell a spatial cell.

[0160] Search for the candidate die grid cells contained in each cell. As Figure 5 shown is the schematic diagram of the cells of a two-dimensional object. There are a total of N = N x ×N y = 8×5 = 40 spatial cells, and only the cell numbers at the four corners are marked in the figure. The cell numbered 21 contains five die grid cells, and their nodes are represented by solid circles. Finding the candidate cells contained in each cell generally includes two steps. First, expand a certain die spatial grid area along three directions, and the expansion distance is half of the size of the sub-spatial grid, and loop through each die cell to obtain the cells containing it, and finally obtain the candidate die cells contained in each cell. For a given particle x p = (x, y, z), the cell number where it is located is

[0161] where [] is a rounding operation. For the particle x p, the subspace lattice containing the mass point can be located. Based on the above global search results, the candidate die unit set of the mass point can be obtained. Next, local search is performed in this candidate set and contact judgment is carried out. To judge the contact state between meshes and the mass point x p , we only need to find whether the mass point contacts or penetrates the candidate unit. The following steps can be executed to achieve this goal.

[0162] Find potential contact units: As Figure 6 shown, for the mass point x p , there are four candidate units S1, S2, S3, S4. Find the node m p nearest to x s , which does not coincide with x p . If the following inequalities are satisfied, the mass point will contact the S i unit

[0163] (C i ×g2)·(C i ×C i+1 )>0, (C i ×g2)·(g2×C i+1 )>0 (25)

[0164] Here, i = 1, 2, 3, 4 is the i-th die unit to be judged. C i and C i+1 are two unit vectors with a common starting point, along two adjacent sides of the S i unit. g1 starts from m s and ends at x p . g2 is the projection of the unit S i on the vector g1, written as:

[0165] g2 = g1 - (g1·m)m (26)

[0166] Here

[0167]

[0168] As Figure 7 shown, to calculate the spatial coordinates of x c , a parametric formula is established on the unit:

[0169] r(ξ, η) = f1(ξ, η)i1 + f2(ξ, η)i2 + f3(ξ, η)i3 (28)

[0170] where N j (ξ, η) is the bilinear function of the shell element, is the i-th coordinate of node j, (ξ, η) is the local coordinate of any point within the element, and i1, i2, i3 are the unit vectors in the spatial coordinates. t is the vector starting from the original coordinates to x p ending. Since the vector is perpendicular to the contact element, the local coordinates (ξ c , η c ) of the contact point x must satisfy:

[0171]

[0172]

[0173] Judgment of the contact state and calculation of the penetration depth d:

[0174]

[0175] where n c is the unit outward normal vector of the contact element at the point (ξ c , η c ), and the calculation method is:

[0176]

[0177] Step 13: Solve the contact force;

[0178] As Figure 8 shown, the resistance F R is a function of d, and its expression is:

[0179] F R =-km p dn c / (△t n ) 2 (33)

[0180] where m p is the mass of the particle and k is the interface stiffness.

[0181] The expressions for the normal and tangential contact forces applied to the particle are:

[0182]

[0183]

[0184] In the formula, is the normal contact force of particle p, is the tangential contact force of particle p, is the resistance, are the normal vectors of the contact points. The total contact force can be expressed as:

[0185]

[0186] where τ ip is the unit tangential vector of the contact point, μ is the friction coefficient, is the total contact force of the particle p.

[0187] Step 14: Update the physical quantities of the particles again:

[0188] For particle x p the updated coordinates can be written as:

[0189]

[0190]

[0191] In the formula, is the velocity of particle p at time n + 1 / 2, is the total contact force of particle p, m p is the mass of particle p, is the coordinate of particle p at time n + 1 / 2.

[0192] Finally, update the background grid momentum:

[0193]

[0194] In the formula, is the momentum of node I at time n + 1 / 2, n p is the number of particles, N Ip the shape function of particle p, is the velocity of particle p at time n + 1 / 2.

[0195] Step 15: Apply the boundary conditions:

[0196] For the nodes on the fixed edge, let

[0197] Step 16: Calculate the strain increment and spin increment of the particles, solve the stress and elastoplastic related variables of the particles according to the constitutive model, and update the density of all particles. Calculate the strain increment and spin increment of the material points:

[0198]

[0199]

[0200] In the formula, is the strain increment, is the spin increment, is the velocity of the Ith node at time n + 1 / 2, is the velocity of the Ith node at time n + 1 / 2.

[0201] Update the density of the material point:

[0202]

[0203] In the formula, is the density of the p-th material point at the (n + 1)-th moment, is the density of the p-th material point at the n-th moment. In summary, the present invention discloses a material point method based on an improved contact algorithm for bicolor coin embossing forming simulation, and the method includes the following steps:

[0204] Step 1: Establish a dynamic physical model for the bicolor coin embossing forming process;

[0205] Step 2: Perform material point integration on the dynamic model of the embossing forming process;

[0206] Step 3: Discretize the forming die with finite element triangular / quadrilateral elements;

[0207] Step 4: Establish a background grid according to the blank data model, discretize the bicolor coin blank with tetrahedral / hexahedral elements, and establish a set of material points for the copper and aluminum material regions;

[0208] Step 5: Calculate the node loads of the background grid;

[0209] Step 6: Map the physical quantities carried by the material points to the background grid nodes through interpolation shape functions;

[0210] Step 7: Solve the system adaptive time step according to the different material properties of the bicolor coin;

[0211] Step 8: Solve the momentum equation and apply boundary conditions;

[0212] Step 9: Update the physical quantities of the material points;

[0213] Step 10: Update the physical quantities of the background grid points;

[0214] Step 11: Update the spatial position of the die;

[0215] Step 12: Judge the penetration state of the material points and solve the penetration distance;

[0216] Step 13: Solve the contact force;

[0217] Step 14: Update the physical quantities of the material points again;

[0218] Step 15: Apply boundary conditions;

[0219] Step 16: Calculate the strain increment and rotation rate increment of the material points, solve the stress and elastoplastic related variables of the material points according to the constitutive model, and update the density of all material points.

[0220] When the traditional finite element algorithm is used to handle actual embossing forming or similar machining problems with molds, the material interface contact algorithm needs to additionally establish a material self-contact judgment process. The present invention proposes a material point method that improves the point-to-face contact algorithm, which not only avoids virtual contact to improve the accuracy of the contact judgment of the material point method, but also does not require an additional material self-contact judgment algorithm, providing an effective way for such problems involving material self-contact.

Claims

1. A material point method based on an improved contact algorithm for bicolor coin embossing forming simulation, characterized in that: It includes the following steps: Step 1: Establish a dynamic physical model for the embossing forming process of the two-color coin; Step 2: Perform material point integration on the dynamic model of the embossing forming process; Step 3: Discretize the forming die with finite element triangular / quadrilateral elements; Step 4: Establish a background grid according to the blank data model, discretize the two-color coin blank with tetrahedral / hexahedral elements, and establish a set of material points for the copper and aluminum material regions; Step 5: Calculate the node loads of the background grid; Step 6: Map the physical quantities carried by the material points to the background grid nodes through the interpolation shape function; Step 7: Solve the system adaptive time step according to the different material properties of the two-color coin; Step 8: Solve the momentum equation and apply boundary conditions; Step 9: Update the physical quantities of the material points; Step 10: Update the physical quantities of the background grid points; Step 11: Update the spatial position of the die; Step 12: Judge the penetration state of the material points and solve the penetration distance; Step 13: Solve the contact force; Step 14: Update the physical quantities of the material points again; Step 15: Apply boundary conditions; Step 16: Calculate the strain increment and spin rate increment of the material points, solve the stress and elastoplastic related variables of the material points according to the constitutive model, and update the density of all material points; The main content of Step 1 is to establish a dynamic physical model for the embossing forming process of the two-color coin, and the specific content is as follows: From a mechanical point of view, the embossing forming process is a quasi-static low-speed motion process. The dynamic explicit algorithm is used to describe and solve the embossing forming process. Assuming the blank is the object under study, when it undergoes plastic deformation, it should satisfy the basic equilibrium equation: Among them, σ ij,j is the Cauchy stress tensor respectively, b i is the body force, ρ is the material density, γ is the damping coefficient, and are the velocity and acceleration of any point in the object respectively, i and j are the coordinate directions respectively. According to the principle of virtual displacement, relevant stress boundary conditions and Gauss divergence theorem, the virtual power equation of the system can be derived as: where \(V\) is the object region, is the virtual velocity, is the virtual strain rate corresponding to the Cauchy stress \(\sigma\) ij ; where \(S\) p is the surface of the known external force \(p\) i ; \(S\) c is the surface in contact with another object, and the contact surface force is denoted as \(q\) i ; the power terms of the inertial force and the damping force reflect the inertial effect and the physical damping effect of the object system.

2. The material point method based on an improved contact algorithm applied to the simulation of bicolor coin embossing forming according to claim 1, wherein: Step 2 is the material point integration of the dynamic equation of the embossing forming process, and the specific steps are as follows: In the material point method, each material point occupies a certain volume space, carries a certain mass and other related physical quantities. The material density at any spatial point x can be expressed as: where n p is the number of mass points, m p is the mass of a mass point, δ() is the Dirac function, x p is the coordinate of mass point p, and the displacement u ip of any mass point p in the background grid can be obtained by interpolation using the shape functions of the background grid nodes: Among them, N Ip is the shape function of node I related to particle p, and u iI is the velocity of the background grid node I in the i-th direction. Substitute the above density and displacement formulas into the system virtual power equation in step 1, and at the same time consider that there is no body force b i and external force p i . We can get: Among them, is the mass of node I, and are the velocity and acceleration of node I in the i-th direction respectively, is the external force of node I in the i-th direction, is the internal force of node I in the i-th direction.

3. The material point method based on an improved contact algorithm applied to the simulation of bicolor coin embossing forming according to claim 1, characterized in that: The main content of Step 4 is to extract the set of material points in the two different material regions of copper and aluminum (copper is the inner core and aluminum is the outer ring), and the specific content is as follows: First, the area occupied by the cake is divided into a regular hexahedral background grid, and the number of grids along the x, y, and z directions are Gx, Gy, and Gz respectively, and the integral points of each background grid are calculated. Then, the copper and aluminum material areas are divided into tetrahedral grids respectively. Each small rectangle represents a background grid unit, the solid circle represents a grid node, the small solid triangle represents a copper material point, and the small solid square represents an aluminum material point; then, the background grid unit set I of each tetrahedral unit is found in the copper material area. Cu The specific method of finding this set is as follows: number the background grid units from 1 to G=Gx*Gy*Gz along the x, y, z directions. According to the coordinates of the four nodes of the tetrahedron unit, a space cuboid can be found to just surround the tetrahedron. Assuming that the coordinates of the leftmost front point of the background grid (x min ,y min , z min ) and the rightmost point coordinate (x max ,y max , z max ), according to any point x in space p = background grid number of (x, y, z) Find all background grid sets contained in the cuboid I Cu The set of Gaussian points of all units in this set is defined as point set II Cu According to the surface boundary unit of the copper material, the mesh generation system automatically generates and writes a program to read the surface unit, and the point set II Cu The Gaussian points outside the copper material are removed, and the set formed by the remaining points is the particle set III of the copper material domain. cu Similarly, similar steps are performed on the aluminum material area to obtain the particle set III of the aluminum material Al , finally, the two sets III Cu and III Al The material point set III of the blank is merged, and each material point is assigned corresponding material properties, such as density, volume, Poisson's ratio, yield stress, hardening index, strengthening coefficient and other parameters. In order to accurately describe the initial outline of the coin, the Gaussian point set of each triangular unit on the surface of the blank is added to the blank particle set to form a new blank particle set. The blank surface particle set is divided into three parts: upper surface, side surface and lower surface particle sets, which respectively form contact pairs with the upper die, middle circle and lower die; during the simulation process, these three pairs of contact pairs are respectively subjected to contact judgment processing.

4. The material point method based on the improved contact algorithm applied to the simulation of bicolor coin embossing forming according to claim 1, characterized in that: Step 7 is the determination of the time integration step size of the dynamic explicit central difference algorithm, which specifically includes the following content: Under normal circumstances, for the stamping forming simulation problem of a single material, the central difference dynamic explicit algorithm is used for solution. This algorithm is conditionally stable, and the system time step Δt must be less than the critical step Δt cr , that is: Δt = αΔt cr In the formula, α is a constant, and its value range is: 0.8 ≤ α ≤ 0.98; Due to the CFL (Courant Friedrichs Lewy) condition, in order to ensure the stability of the explicit time integration, when using the finite element method for dynamic explicit analysis, the distance traveled by the wave within one time step must not exceed one element, and the sound speed of the elastic material is: In the formula, E is the elastic modulus, v is the Poisson's ratio, ρ is the material density, and it is assumed that the nominal length of the element is l e , then the critical time step is as follows: As the deformation increases, the nominal size of the element gradually decreases, and the critical time step also gradually decreases. When the mesh undergoes severe distortion, the minimum nominal length of the solution approaches 0, resulting in the time step also approaching 0, thus making the simulation analysis impossible to carry out; The material point method is different from the finite element method. All time integrations need to be carried out on a fixed background grid. When determining the critical time step size, not only the sound speed of the material needs to be considered, but also the speed of the material points needs to be considered. In particular, in the problem of ultra-high-speed impact, the speed of the material points is relatively large, on the same order of magnitude as the material sound speed or even exceeds the material sound speed; therefore, the influence of the material point speed on the simulation calculation cannot be ignored. In the embossing forming simulation process, based on the material point method of the dynamic explicit central difference algorithm, the critical time step size needs to be calculated at each step to make the numerical solution stable: In the formula, c is the sound speed of the particle, and u is the particle velocity. When using a uniform background grid, l g can also be taken as d c , that is, the background grid size; in the present invention, the blank cake has two different materials, namely copper inside and aluminum outside. The size of the time step is determined by the properties of the system and the model material parameters. Since two different materials are involved, it is necessary to calculate the time steps in the two material regions separately and take the smaller value. The time step in the copper material is: In the formula, and are the speed of sound and the velocity of any particle p1 in the particle set III Cu respectively; Similarly, the time step size in the aluminum material is: In the formula, and are the speed of sound and the particle velocity of any particle p2 in particle III Al respectively; It is necessary to calculate the wave velocities of the two materials and the velocities of the material points of the two materials during the embossing process, and then compare the critical time steps of the two to finally determine that the critical time step of the system is the minimum of the two. Therefore, the critical time step of the system is as follows: In the formula, is the time step of Cu, is the time step of Al.

5. The material point method based on the improved contact algorithm applied to the simulation of bicolor coin embossing and forming according to claim 1, characterized in that: Step 12 is to judge the penetration state of the material points and solve the contact force, which specifically includes the following content: According to the traditional material point method theory, the velocities of the die particle j′ and the blank particle j contribute to the background grid node I (because the background grids they are in share the node I). Therefore, it is determined that the die is in contact with the blank particle, but in fact they are not in contact, resulting in false contact. The present invention uses a point-plane contact algorithm to replace the point-point contact algorithm of the traditional material point method to avoid virtual contact and improve the accuracy of simulation calculation. When judging the contact state and calculating the contact values (contact force, penetration depth, and friction force) between each particle on the blank surface and the upper die, lower die, and middle ring, it will cause a large amount of calculation in each time step. To reduce the calculation time of contact judgment, a strategy combining global search and local search is adopted. The global search is to find the candidate die elements in contact with each particle, and the process is as follows: Construct a block containing a specific die grid: The size of the block is defined according to the maximum and minimum spatial coordinates of each die. The coordinates of the lower left corner and the lower right corner of the block are (x min , y min , z min ),) and (x max , y max , z max ), respectively. Then divide the block area into sub-cells. The number of cells in each direction is N x , N y , N z , and number each cell along the x, y, and z directions respectively. Note that this block moves with the die grid. To distinguish this type of cell from the background cell grid, this type of cell is named a spatial cell; Search for the candidate die grid cells contained in each cell. There are a total of N = N x × N y = 8 × 5 = 40 spatial cells. Only the cell numbers at the four corners are marked in the figure; the cell numbered 21 contains five die grid cells, and its nodes are represented by solid circles. Finding the candidate cells contained in each cell generally includes two steps; first, expand a certain die space cell domain in three directions, and the expansion distance is half of the subspace cell size, and loop through each die cell to obtain the cells containing it, and finally obtain the candidate die cells contained in each cell; for a given particle x p = (x, y, z), the cell number where it is located is where [] is a rounding operation for the particle x p , the subspace lattice containing the particle can be located, and the candidate die element set of the particle can be obtained according to the above global search results; next, local search is performed in the candidate set and contact judgment is carried out. In order to judge the contact state between grids and the particle x p , it is only necessary to find whether the particle contacts or penetrates the candidate cell. The following steps can be executed to achieve this goal; Find potential contact unit: particle x p , S1, S2, S3, S4 have four candidate units, find the one closest to x p The nearest node m s , which does not coincide with x p , if the following inequality is satisfied, the particle will contact the S i unit (C i × g2)·(C i × C i+1 ) > 0, (C i × g2)·(g2 × C i+1 ) > 0 Here, i = 1, 2, 3, 4 represents the i-th die unit to be judged, C i and C i+1 are two unit vectors that have a common starting point and are along S i Two adjacent sides of the unit, g1 starts from m s and ends at x p End, g2 is the projection of the unit S i On the vector g1, written as: g2 = g1 - (g1·m)m Here To calculate the spatial coordinates of x c A parametric formula is established on the element: r(ξ, η) = f1(ξ, η)i1 + f2(ξ, η)i2 + f3(ξ, η)i3 Among them N j (ξ, η) is a bilinear function of the shell element, is the i-th coordinate of node j, (ξ, η) is the local coordinate of any point within the element, and i1, i2, i3 are the unit vectors in the spatial coordinate system; t is the vector starting from the original coordinate to x p ending, because the vector is perpendicular to the contact element, and the local coordinates (ξ c , η c ) of the contact point x must satisfy: Judgment of the contact state and calculation of the penetration depth d: where n c is the unit outward normal vector of the contact element at the point (ξ c , η c ), and the calculation method is as follows:

Images