A method of simulating the conversion of a powder into a continuum during a compaction process

Abstract

The application provides a method for simulating powder transformation into a continuum in a pressing process, which comprises the following steps: discretizing a powder pressing area into two areas of discrete elements and finite elements which are completely overlapped; selecting a space mapping mode to map discrete element particles in the discrete element area to finite element units in the corresponding finite element area; introducing a contact constitutive relation between the discrete element particles and introducing a viscoelastic plastic constitutive relation on the finite element units in the finite element area; selecting a transformation model of the discrete element into the finite element and setting transformation parameters; creating a mold and upper and lower pressing plates around the powder; selecting a quasi-static velocity load applied on the set pressing plate; calculating the coupling transformation of the discrete element and the finite element, calculating a transformation factor of the finite element unit according to a transformation factor of the discrete element particle in the overlapped area; and calculating the overall stress of the pressing area. The application can truly reproduce the discrete-continuum gradual transformation process in the powder pressing process and has the advantages of good physical consistency and high simulation precision.

Description

Technical Field

[0001] This invention relates to the field of powder forming and numerical simulation technology, specifically to a method for simulating the transformation of powder into a continuum during the pressing process. Background Technology

[0002] The powder pressing process involves pressing molding powder into a continuous body under high pressure. The size of the molding powder is in the range of tens of micrometers to millimeters, the size of the continuous body is in the range of centimeters to meters, the pressing pressure is in the range of megapascals to hundreds of megapascals, and the pressing time is in the range of minutes to hours.

[0003] During the loading process, the discrete particle system of molding powder gradually exhibits motion, deformation, compaction, and adhesion phenomena, thus achieving a gradual transformation from a discrete body to a continuous body. The finite element method (FEM) and the discrete element method (DEM) are two main numerical simulation methods for molding. The FEM focuses on the changes in macroscopic mechanical behavior during molding, such as the evolution of the load-displacement curve and the evolution of spatial density. The computational accuracy of this method mainly depends on the accuracy of the macroscopic constitutive model. However, the establishment of the macroscopic constitutive model and the acquisition of constitutive parameters mainly rely on experiments. Since macroscopic experimental data are often limited, it is impossible to accurately characterize the macroscopic evolution of molding powder during molding, let alone directly predict the macroscopic mechanical behavior during molding based on the microscopic parameters of the molding powder. The DEM focuses on the microscopic behavior between molding powder particles. By meticulously establishing a microscopic model of the molding powder and introducing microscopic parameters, it can accurately simulate the motion, deformation, and compaction behavior of the molding powder during molding. However, limited by the computational complexity of the DEM, this type of algorithm is only suitable for small-scale molding mechanism analysis.

[0004] In general, current numerical calculation methods either focus only on the calculation of the continuous part (such as the finite element method) or only describe discrete particles (such as the discrete element method), and lack a calculation method for quantitatively analyzing the transformation of molding powder from a discrete state to a continuous state during the pressing process. Summary of the Invention

[0005] To address the technical problems existing in the background art, the present invention proposes a method for simulating the transformation of powder into a continuum during the pressing process. This method can realistically reproduce the discrete-to-continuous gradual transformation process during powder pressing, and has advantages such as good physical consistency and high simulation accuracy. It can solve the technical problems in the prior art, such as low computational efficiency of microscopic simulation and inability of macroscopic simulation to consider the influence of microstructure.

[0006] To solve the above-mentioned technical problems, the present invention provides a method for simulating the transformation of powder into a continuous body during the pressing process, which includes the following steps:

[0007] (1) Discretize the powder pressing area into two regions: discrete element and finite element, with the two regions completely overlapping;

[0008] (2) Select the spatial mapping mode to map the discrete element particles in the discrete element region to the finite element elements in the corresponding finite element region;

[0009] (3) Introduce contact constitutive models between particles in the discrete element and introduce viscoelastic-plastic constitutive models on the elements of the finite element in the finite element region;

[0010] (4) Select the conversion model for converting discrete elements to finite elements and set the conversion parameters;

[0011] (5) Create a mold and upper and lower pressure plates around the powder;

[0012] (6) According to the set loading mode, select to apply a quasi-static velocity load on the set pressure plate;

[0013] (7) Calculate the coupling transformation between discrete element and finite element, and calculate the transformation factor of finite element based on the transformation factor of discrete element particles in the superimposed region;

[0014] (8) Take the powder pressing area as a model and calculate the overall force of the pressing area based on the conversion factor of discrete element particles and the conversion factor of finite element elements.

[0015] As a preferred embodiment of the present invention, the minimum mesh size of the finite element unit is greater than or equal to the maximum particle diameter of the discrete element particles, specifically expressed as:

[0016] ;

[0017] In the above formula, For the first The feature size of each unit; For the first The diameter of each particle.

[0018] As a preferred embodiment of the present invention: the spatial mapping mode in step (2) includes two types: body-centered mapping mode and area mapping mode;

[0019] When the spatial mapping mode is a body-centered mapping mode, if the particle body-centered Located inside a finite element, the particle forms a bonded relationship with that element; a finite element is composed of... Composed of 1 facet, the 1st facet The equations for the planes are:

[0020] ;

[0021] In the above formula, Let be the unit outward normal vector of the plane. It is a constant;

[0022] If the following equation is satisfied, it means that the particle's body center is located inside a certain finite element:

[0023] ;

[0024] When the spatial mapping mode is an area mapping mode, the intersection of the discrete element particle and each finite element element needs to be calculated. If the intersection is not empty, then the discrete element particle and the finite element element form a binding relationship. Let the spatial region occupied by the discrete element particle be... The spatial region occupied by the finite element element is ;like Then the intersection volume between the discrete element particle and each finite element is defined as: ,Right now:

[0025] ;

[0026] In the above formula, For discrete element particles, , This is the indicator function for the finite element; when the intersection is empty, the integral value is 0.

[0027] As a preferred embodiment of the present invention: the contact constitutive model in step (3) is any one of the linear contact model, the Hertz contact model and the DPC contact model; the viscoelastic-plastic constitutive model in step (3) is any one of the DPC model, the MC model and the Burgers model.

[0028] As a preferred embodiment of the present invention, the conversion model for converting discrete elements into finite elements in step (4) is either the volume stress conversion model or the bulk modulus conversion model.

[0029] In the volumetric stress transformation model, the discrete element particle transformation factor The calculation formula is expressed as:

[0030] ;

[0031] In the above formula, The body stress at the start of the transformation. The body stress at the point of termination of the transformation;

[0032] When the volume stress experienced by a discrete element particle is ≥ The transformation begins; when the volume stress experienced by a discrete element particle is ≥ The conversion is complete;

[0033] In the bulk modulus conversion model, the discrete element particle conversion factor The calculation formula is expressed as:

[0034] ;

[0035] In the above formula, The equivalent bulk modulus of the region where the discrete element particles are located at the start of the transformation. The equivalent bulk modulus of the region where the discrete element particles are located at the time of termination of the transformation;

[0036] When the equivalent bulk modulus of the region where a discrete element particle is located is ≥ The transformation begins; when the equivalent bulk modulus of the region where a discrete element particle is located is ≥ The conversion is complete.

[0037] As a preferred embodiment of the present invention, the discrete element calculation in step (7) includes the calculation of the contact force of discrete element particles, the force calculation of discrete element particles, and the motion calculation of discrete element particles.

[0038] The specific process for calculating the contact force of the discrete element particles is as follows:

[0039] The calculation of interparticle contact force includes two parts: normal force calculation and tangential friction calculation. The formula for calculating the normal force when using a linear contact model is as follows:

[0040] ;

[0041] In the above formula, , These are the normal contact stiffness and normal damping coefficient. , For embedding quantity and embedding rate;

[0042] The tangential force calculation uses an incremental-update mode, with each time step... First, calculate the test tangential force, then correct it according to Coulomb's law of friction; the formula for calculating the test stress is:

[0043] ;

[0044] In the above formula, To test the tangential force at the current time step, The tangential force of the previous time step, For tangential stiffness, This represents the increment of the tangential displacement difference;

[0045] Next, calculate the maximum static friction force at the current time step, which is:

[0046] ;

[0047] In the above formula, For the maximum static friction force, The coefficient of friction;

[0048] Next, a correction for the tangential force is performed, calculated using the following formula:

[0049] ;

[0050] The specific process for calculating the force on the discrete element particles is as follows:

[0051] The formulas for calculating the resultant force and resultant torque on discrete element particles are as follows:

[0052] ;

[0053] In the above formula, and The net force and net torque acting on the particle are shown in order, respectively. and The forces applied to the particles are gravity and torque, respectively. and These are, in order, the contact force and the connecting torque between particles; This refers to the contact force between the particle and the rigid surface. This is the position vector of the point of contact between particles relative to the center of mass of the particle; is the position vector of the point of action between the particle and the rigid surface relative to the particle's center of mass;

[0054] The motion calculation of the discrete element particles includes particle translation calculation and particle rotation calculation:

[0055] The formula for calculating particle translational motion is:

[0056] ;

[0057] In the above formula, , , The vectors are acceleration, velocity, and displacement, respectively. For particle mass, This is the local translational damping coefficient;

[0058] The formula for calculating particle rotation is:

[0059] ;

[0060] In the above formula, , , The values ​​are, in order, rotational angular acceleration, angular velocity, and angular displacement; This is the local rotational damping coefficient; It is the moment of inertia; The quantity is represented by a quaternion of rotational speed.

[0061] As a preferred embodiment of the present invention, the finite element calculation in step (7) includes the calculation of the deformation and force of the finite element element and the calculation of the force and motion of the finite element node.

[0062] The calculation process for the deformation and stress of the finite element is as follows:

[0063] Calculate a Gaussian point in a given element based on the nodal displacement increment. The strain increment is calculated using the following formula:

[0064] ;

[0065] In the above formula, Let be the strain increment matrix at a certain Gaussian point in a certain element. The strain matrix at the Gaussian point. This is the element strain increment matrix.

[0066] Calculate the stress increment at the current moment based on the strain increment. The calculation formula is:

[0067] ;

[0068] In the above formula, The material matrix is ​​determined by the constitutive model used by the solid element.

[0069] Calculate the total stress based on the stress increment at the current moment. The calculation formula is:

[0070] ;

[0071] In the above formula, This represents the total stress from the previous time step;

[0072] The calculation process for the force-induced motion of the finite element nodes is as follows:

[0073] First, calculate the resultant force at each node of the finite element element. The calculation formula is as follows:

[0074] ;

[0075] In the above formula, The net external force acting on a finite element node. The force generated by the deformation of the finite element on this node. The external force applied to this node;

[0076] Then, based on the resultant force of each node Calculate the acceleration of the node; if the node's current velocity is greater than zero, then the acceleration... If the node's current velocity is less than or equal to zero, then the acceleration... ;in, This is the local damping coefficient. For node quality;

[0077] Next, the velocity is updated based on the acceleration, and then the displacement increment and displacement of this step are updated based on the velocity. The calculation formula is as follows:

[0078] ;

[0079] In the above formula, The nodal acceleration vector; The node velocity vector; This is the vector of nodal displacement increments; The nodal displacement vector; This is for calculating the time step.

[0080] As a preferred embodiment of the present invention: the quasi-static velocity load can be applied to the upper pressure plate or simultaneously to the upper and lower pressure plates, and the range of the quasi-static load is set to 10. -7 ~10 -9 m / step.

[0081] As a preferred embodiment of the present invention, the formula for calculating the transformation factor of the finite element element is expressed as follows:

[0082] ;

[0083] In the above formula, The total number of discrete element particles that are bound to a certain finite element element. For the first The volume of a discrete element particle. For the first The transformation factor of each discrete element particle.

[0084] As a preferred embodiment of the present invention, the formula for the overall force of the calculation model is:

[0085] ;

[0086] In the above formula, For a certain finite element node The combined force on the area These are the nodal forces obtained directly from finite element analysis. The nodal forces are directly calculated from discrete element particles. Let be the weight coefficient of the i-th discrete element particle at the j-th finite element node;

[0087] calculate First, the nodal force matrix of each node in the finite element is calculated based on the stress at each Gaussian point in the finite element. The calculation formula is:

[0088] ;

[0089] In the above formula, Gaussian point The strain matrix, This represents the total stress at the current time step. and Gaussian point The weight coefficients and the values ​​of the Jacobian determinant;

[0090] Next, in the nodal force matrix Find the node The corresponding nodal force vector is calculated using the following formula:

[0091] ;

[0092] Particles The contact force exerted on it by other particles The calculation formula is:

[0093] ;

[0094] In the above formula, To be with particles The number of other particles in contact, For the first Individual and Particle Contacting particles act on particles The force on top.

[0095] As a preferred embodiment of the present invention: the method of converting powder into a continuum can be used for two-dimensional problem simulation or three-dimensional problem simulation; when the numerical method performs two-dimensional simulation, the discrete element particles are circles, ellipses or polygons, and the finite element units are triangles or quadrilaterals; when the numerical method performs three-dimensional simulation, the discrete element particles are spheres, ellipsoids or polyhedra, and the finite element units are tetrahedrons, pyramids, triangular prisms or hexahedrons.

[0096] By adopting the above technical solution, the present invention has the following beneficial effects:

[0097] This invention couples the Discrete Element Method (DEM) with the Finite Element Method (FEM), effectively overcoming the limitation of traditional numerical simulations that cannot achieve integrated calculations of fine and macroscopic dimensions. By employing a progressive transformation strategy from DEM to FEM, it ensures a natural, oscillating, and robust transformation process, realistically reproducing the entire mechanical response from loose particles to a dense continuum. This invention constructs a spatially superimposed overlay model of DEM and FEM, using volume stress and bulk modulus as transformation criteria, and dynamically allocating the contributions of the DEM and FEM to the spatial region using a transformation factor, thus achieving a progressively smooth transformation from loose powder media to dense continuum media.

[0098] The main features and advantages of this invention are reflected in the following aspects:

[0099] (1) Achieved high-fidelity simulation integrating micro- and macro-scale: This invention innovatively superimposes the discrete element method (DEM) and the finite element method (FEM) in space and achieves dynamic coupling between the two through a transformation factor. This method can accurately capture the micro-mechanical behavior of powder particles in the early stage of pressing, such as movement, rearrangement, and crushing, and can also efficiently simulate the macro-stress-strain response of the dense continuum in the later stage of pressing by progressively transforming it into the finite element method. This breaks through the limitation of traditional single numerical methods that cannot simultaneously take into account both micro-mechanisms and macro-scale performance.

[0100] (2) The transformation process has clear physical meaning and high numerical stability: The transformation model based on volume stress or bulk modulus proposed in this invention has clear physical meaning. The transformation factor smoothly transitions between 0 and 1, avoiding non-physical oscillations or numerical mutations when switching between discrete and finite element methods. At the same time, by weighting the resultant force of the nodes through the transformation factor, the continuity of the system's mechanical response during the transformation process is ensured. The algorithm is robust and can realistically reproduce the entire mechanical evolution from loose powder to dense continuum.

[0101] (3) Strong applicability, capable of simulating complex geometries and multi-material systems: This invention provides two flexible spatial mapping modes and supports combinations of various discrete element contact constitutive and finite element constitutive methods. This method supports two-dimensional and three-dimensional problem analysis, allows powder particles to be spherical, ellipsoidal, polygonal or polyhedral, and can flexibly adapt to the simulation of pressing processes of powders with different shapes, complex mold geometry and different material systems. Attached Figure Description

[0102] To more clearly illustrate the specific embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the specific embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained from these drawings without creative effort.

[0103] Figure 1 This is a schematic diagram of the process for simulating the conversion of powder into a continuum during the pressing process of the present invention;

[0104] Figure 2 This is a schematic diagram illustrating the transformation process from discrete element to finite element method.

[0105] Figure 3 This is a vertical displacement cloud diagram during the two-dimensional pressing process in Embodiment 1 of the present invention;

[0106] Figure 4 This is a graph showing the changes in finite element transformation factor and top plate stress with compression displacement during the two-dimensional pressing process in Embodiment 1 of the present invention.

[0107] Figure 5 This is a vertical stress cloud diagram during the three-dimensional pressing process in Embodiment 2 of the present invention;

[0108] Figure 6 This is a schematic diagram of the finite element mesh and discrete element particles in the computational region involved in the method of simulating the transformation of powder into a continuum during the pressing process of this invention. Detailed Implementation

[0109] The technical solution of the present invention will now be clearly and completely described with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0110] The present invention will be further explained below with reference to specific embodiments.

[0111] like Figure 1 As shown in the figure, this embodiment provides a method for simulating the conversion of powder into a continuous body during the pressing process, which includes the following steps:

[0112] Step 100: Discretize the powder pressing area into two regions: a discrete element method and a finite element method, with the two regions completely overlapping; in specific implementation, such as... Figure 6 As shown, the powder region is discretized using two sets of meshes: discrete element and finite element. The discrete element mesh consists of a series of spherical or non-spherical particles that are in contact with each other and meet certain porosity and loose packing density requirements. The finite element mesh consists of regular or irregular elements that cover the entire powder region, with each element tightly connected without gaps. The mold and upper and lower pressure plates around the powder are described using finite element elements.

[0113] Step 200: Select the spatial mapping mode to map the discrete element particles in the discrete element region onto the finite element elements in the corresponding finite element region.

[0114] Step 300: Introduce contact constitutive models between particles in the discrete element and introduce viscoelastic-plastic constitutive models on the elements of the finite element in the finite element region.

[0115] Step 400: Select the conversion model for converting discrete element method to finite element method and set the conversion parameters;

[0116] Step 500: Create a mold and upper and lower pressure plates around the powder;

[0117] Step 600: According to the set loading mode, select the setting platen to apply a quasi-static velocity load. The set loading mode includes three modes: one is that only the upper platen loads downward at a certain quasi-static loading rate; the second is that only the lower platen loads upward at a certain quasi-static velocity; and the third is that both the upper and lower platens load towards the center of the powder at a certain quasi-static velocity simultaneously.

[0118] Step 700: Calculate the coupling transformation between discrete element and finite element, and calculate the transformation factor of finite element based on the transformation factor of discrete element particles in the superimposed region.

[0119] Step 800: Treat the powder pressing area as a model, and calculate the overall force of the pressing area based on the conversion factors of discrete element particles and finite element elements.

[0120] The numerical method of this invention can perform two-dimensional problem simulation as well as three-dimensional problem simulation. In two-dimensional simulation, the discrete element particles can be circles, ellipses or polygons; the finite element elements can be triangles or quadrilaterals. In three-dimensional simulation, the discrete element particles can be spheres, ellipsoids or polyhedra, and the finite element elements can be tetrahedrons, pyramids, triangular prisms or hexahedrons.

[0121] The minimum mesh size of the finite element element mentioned above is greater than or equal to the maximum particle diameter of the discrete element particles; the calculation formula can be expressed as:

[0122] ;

[0123] in, For the first The feature size of each unit For the first The diameter of each particle.

[0124] The aforementioned spatial mapping modes include two categories: body-centered mapping mode and area mapping mode.

[0125] When the space mapping mode is selected as body-centered mapping mode, if the particle body-centered Located inside a finite element, the particle forms a bonded relationship with that element; a finite element is composed of... Composed of 1 facet, the 1st facet The equations for the planes are:

[0126] ;

[0127] in, Let be the unit outward normal vector of the plane. It is a constant.

[0128] If the following equation is satisfied, it means that the particle's body center is located inside a certain finite element:

[0129] ;

[0130] When the area mapping mode is selected in the spatial mapping mode, the intersection of the discrete element particle and each finite element element needs to be calculated. If the intersection is not empty, then the discrete element particle and the finite element element are bound together. Let the spatial region occupied by the discrete element particle be... The spatial region occupied by the finite element element is ;like (Indicating that the intersection of two spatial regions is non-empty), then the intersection volume of the discrete element particle and each finite element is defined as: ,Right now:

[0131] ;

[0132] In the above formula, For discrete element particles, , This is the indicator function for the finite element; when the intersection is empty, the integral value is 0.

[0133] The contact constitutive model of the discrete element described above can be any one of the linear contact model, Hertz contact model, or DPC contact model.

[0134] The element constitutive model of the finite element method described above can be any one of the DPC model, MC model, or Burgers model.

[0135] The transformation model for converting discrete elements to finite elements can be either the volume stress transformation model or the bulk modulus transformation model.

[0136] In the volumetric stress transformation model, the discrete element particle transformation factor The calculation formula can be expressed as:

[0137] ;

[0138] In the above formula, The body stress at the start of the transformation. The body stress at the point of termination of the transformation.

[0139] When the volume stress experienced by a discrete element particle is ≥ The transformation begins; when the volume stress experienced by a discrete element particle is ≥ The conversion is complete.

[0140] In the bulk modulus conversion model, the discrete element particle conversion factor The calculation formula can be expressed as:

[0141] ;

[0142] In the above formula, The equivalent bulk modulus of the region where the discrete element particles are located at the start of the transformation. It is the equivalent bulk modulus of the region where the discrete element particles are located at the time of termination of the transformation.

[0143] When the equivalent bulk modulus of the region where a discrete element particle is located is ≥ The transformation begins; when the equivalent bulk modulus of the region where a discrete element particle is located is ≥ The conversion is complete.

[0144] The aforementioned mold and upper and lower pressure plates can be described using rigid surface element units or finite element units.

[0145] Step 700 involves calculating the discrete element, including the contact force of discrete element particles and the force-induced motion of discrete element particles.

[0146] The specific process for calculating the contact force of the discrete element particles is as follows:

[0147] The calculation of interparticle contact force includes two parts: normal force calculation and tangential friction calculation. The formula for calculating the normal force when using a linear contact model is as follows:

[0148] ;

[0149] In the above formula, , These are the normal contact stiffness and normal damping coefficient. , For embedding quantity and embedding rate;

[0150] The tangential force calculation uses an incremental-update mode, with each time step... First, calculate the test tangential force, then correct it according to Coulomb's law of friction; the formula for calculating the test stress is:

[0151] ;

[0152] In the above formula, To test the tangential force at the current time step, The previous step tangential force, For tangential stiffness, This represents the increment of the tangential displacement difference;

[0153] Next, calculate the maximum static friction force at the current time step, which is:

[0154] ;

[0155] In the above formula, For the maximum static friction force, The coefficient of friction;

[0156] Next, a correction for the tangential force is performed, calculated using the following formula:

[0157] ;

[0158] The specific process for calculating the force and motion of the discrete element particles is as follows:

[0159] The formula for calculating the resultant force and resultant moment acting on discrete element particles is as follows:

[0160] ;

[0161] In the above formula, and The net force and net torque acting on the particle are given by the particle. and The force and torque applied to the particles, and The contact force and connecting torque between particles. The contact force between the particle and the rigid surface. and For the contact points between particles and for the point of interaction between a particle and a rigid surface, the position vectors relative to the particle's center of mass are denoted as .

[0162] Discrete element particle motion calculations include particle translation and particle rotation calculations:

[0163] The formula for calculating particle translational motion is:

[0164] ;

[0165] In the above formula, , , For acceleration, velocity, and displacement vectors, For particle mass, This is the local translational damping coefficient;

[0166] The formula for calculating particle rotation is:

[0167] ;

[0168] In the above formula, , , For rotational angular acceleration, angular velocity, and angular displacement, The rotational local damping coefficient is... For rotational inertia, The quantity is represented by a quaternion of rotational speed.

[0169] In step 700, the finite element calculation is performed using the incremental finite element method, which includes the calculation of the deformation and stress of finite element elements and the calculation of the force and motion of finite element nodes.

[0170] The calculation process for the deformation and stress of the finite element is as follows:

[0171] Calculate a Gaussian point in a given element based on the nodal displacement increment. The strain increment is calculated using the following formula:

[0172] ;

[0173] in, Let be the strain increment matrix at a certain Gaussian point in a certain element. The strain matrix at the Gaussian point. This is the element strain increment matrix.

[0174] Calculate the stress increment at the current moment based on the strain increment. The calculation formula is:

[0175] ;

[0176] in, The material matrix is ​​determined by the constitutive model used by the solid element.

[0177] Calculate the total stress based on the stress increment at the current moment. The calculation formula is:

[0178] ;

[0179] in, This represents the total stress from the previous time step;

[0180] The calculation process for the force-induced motion of the finite element nodes is as follows:

[0181] First, calculate the resultant force at each node of the finite element element. The calculation formula is as follows:

[0182] ;

[0183] In the above formula, The net external force acting on a finite element node. The force generated by the deformation of the finite element on this node. The external force applied to this node;

[0184] Then, based on the resultant force of each node Calculate the acceleration of the node; if the node's current velocity is greater than zero, then the acceleration... If the node's current velocity is less than or equal to zero, then the acceleration... ;in, This is the local damping coefficient. For node quality;

[0185] Next, the velocity is updated based on the acceleration, and then the displacement increment and displacement of this step are updated based on the velocity. The calculation formula is as follows:

[0186] ;

[0187] In the above formula, For nodal acceleration vectors, The node velocity vector; This is the vector of nodal displacement increments; The nodal displacement vector; This is for calculating the time step.

[0188] The coupling and transformation between the finite element method (FEM) and the discrete element method (DIM) are mainly achieved through a transformation factor. After each iteration step, the contact force obtained by each discrete particle is interpolated to the nodes of the finite element method (FEM), and then the resultant force at the nodes of the finite element method (FEM) is calculated using the transformation factor between the FEM and DEM.

[0189] The aforementioned quasi-static velocity load can be applied to the upper pressure plate or simultaneously to both pressure plates. The range of the quasi-static load can be set to 10. -7 ~10 -9 m / step.

[0190] The formula for calculating the transformation factor of the finite element element described above can be expressed as follows:

[0191] ;

[0192] In the above formula, The total number of discrete element particles that are bound to a certain finite element element. For the first The volume of a discrete element particle. For the first The transformation factor of each discrete element particle.

[0193] The formula for the overall force in the above calculation model can be expressed as:

[0194] ;

[0195] In the above formula, For a certain finite element node The combined force on the area These are the nodal forces obtained directly from finite element analysis. The nodal forces are directly calculated from discrete element particles. Let be the weight coefficient of the i-th discrete element particle at the j-th finite element node;

[0196] calculate First, the nodal force matrix of each node in the element is calculated based on the stress at each Gaussian point in the finite element. The calculation formula can be expressed as:

[0197] ;

[0198] in, Gaussian point The strain matrix, This represents the total stress at the current time step. and Gaussian point The weight coefficients and the values ​​of the Jacobian determinant;

[0199] Next, in the nodal force matrix Find the node The corresponding nodal force vector is calculated using the following formula:

[0200] ;

[0201] Particles The contact force exerted on it by other particles The calculation formula is:

[0202] ;

[0203] In the above formula, To be with particles The number of other particles in contact, For the first Individual and Particle Contacting particles act on particles The force on top.

[0204] Example 1:

[0205] This embodiment is a two-dimensional compression simulation case, employing a volumetric stress transformation mode to achieve the transition from discrete element method (DEM) to finite element method (FEM). Specifically, a rectangular model with a width of 0.2m and a height of 1m is constructed. It is divided into 2726 finite element triangular elements and 5210 discrete element particles. Normal constraints are applied to the bottom and left and right sides of the model, and a 10-degree stress is applied to the top of the model.-8 Quasi-static vertically downward velocity load at m / step.

[0206] The contact model between discrete element particles is the Hertz contact model, and the constitutive model for the finite element method is the DPC constitutive model. The density of the molding powder particles is 1800 kg / m³. 3 The model has an elastic modulus of 1 GPa, a Poisson's ratio of 0.25, a cohesive strength and tensile strength of 0.8 MPa, an internal friction angle of 35°, and a local damping of 0.1. Volume stress is used as the criterion for the transition from a discrete medium to a continuous medium, with an initial transition stress of 5 MPa and a final transition stress of 15 MPa.

[0207] The evolution of the vertical displacement of the model during the calculation process is as follows: Figure 3 As shown in the figure, under the vertical displacement load at the top of the model, the modeling powder particles exhibit downward movement and compaction. When the volume stress of the particles exceeds the lower limit of transformation, the contribution of discrete element particles to the model gradually decreases, while the contribution of finite element elements gradually increases. When the transformation factor of the discrete element is equal to 1, the discrete element particles are hidden; when the transformation factor of the finite element is equal to 1, the finite element elements are displayed. Therefore, the process of discrete element particles gradually evolving into finite element elements can be seen from the figure.

[0208] The variation of the finite element average transformation factor and the top vertical stress with the loading displacement of the model are as follows: Figure 4 As shown in the figure, as the top displacement load increases, the average conversion factor of the finite element model gradually increases, and the rate of increase gradually slows down as the displacement increases. As the displacement load increases, the vertical stress at the top of the model first increases linearly. When the lower limit of conversion is reached, the rate of increase of the vertical stress gradually increases. After all discrete elements are converted into finite elements, the stress increases linearly again, but the slope is significantly greater than that of the initial linear segment.

[0209] Example 2:

[0210] This embodiment is a three-dimensional compression simulation case, which uses the bulk modulus conversion mode to realize the transition from discrete element to finite element.

[0211] A cylindrical model with a diameter of 0.2m and a height of 1m was constructed, consisting of 138,590 finite element tetrahedral elements and 26,294 discrete element particles. Normal constraints were applied to the bottom and sides of the model, and a 10-degree constraint was applied to the top. -9 m / step quasi-static vertical downward velocity.

[0212] The contact between discrete element particles was modeled using a linear contact model, and the finite element model used a Monte Carlo (MC) constitutive model. The density of the molding powder particles was 1800 kg / m³. 3The model has an elastic modulus of 1 GPa, a Poisson's ratio of 0.25, a cohesive strength and tensile strength of 0.8 MPa, an internal friction angle of 35°, and a local damping of 0.1. The bulk modulus is used as the criterion for the transformation from a discrete medium to a continuous medium, with an initial transformation stress of 10 MPa and a final transformation stress of 500 MPa.

[0213] During loading, the gradual transformation of discrete element particles into finite element elements is as follows: Figure 5 As shown in the figure, as the vertical displacement increases, the vertical stress of the model gradually increases, and the discrete element particles gradually evolve into finite element elements. When the vertical compressive stress of the model reaches 60~70MPa, the transformation is basically completed.

[0214] This invention can realistically reproduce the discrete-to-continuous gradual transition process during powder pressing, and has the advantages of good physical consistency and high simulation accuracy. It can solve the technical problems in the prior art, such as low computational efficiency of microscopic simulation and inability of macroscopic simulation to consider the influence of microstructure.

[0215] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features therein. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for simulating the transformation of powder into a continuous body during the pressing process, characterized in that, Includes the following steps: (1) Discretize the powder pressing area into two regions: discrete element and finite element, with the two regions completely overlapping; (2) Select the spatial mapping mode to map the discrete element particles in the discrete element region to the finite element elements in the corresponding finite element region; (3) Introduce contact constitutive models between particles in the discrete element and introduce viscoelastic-plastic constitutive models on the elements of the finite element in the finite element region; (4) Select the conversion model for converting discrete elements to finite elements and set the conversion parameters; (5) Create a mold and upper and lower pressure plates around the powder; (6) According to the set loading mode, select to apply a quasi-static velocity load on the set pressure plate; (7) Calculate the coupling transformation between discrete element and finite element, and calculate the transformation factor of finite element based on the transformation factor of discrete element particles in the superimposed region; (8) Take the powder pressing area as a model and calculate the overall force of the pressing area based on the conversion factor of discrete element particles and the conversion factor of finite element elements.

2. The method for converting powder into a continuous body during the simulated pressing process as described in claim 1, characterized in that, The minimum mesh size of the finite element element is greater than or equal to the maximum particle diameter of the discrete element particles, specifically expressed as: ； In the above formula, For the first The feature size of each unit; For the first The diameter of each particle.

3. The method for converting powder into a continuous body during the simulated pressing process as described in claim 1, characterized in that: The spatial mapping modes in step (2) include two types: body-centered mapping mode and area mapping mode; When the spatial mapping mode is a body-centered mapping mode, if the particle body-centered Located inside a finite element, the particle forms a bonded relationship with that element; a finite element is composed of... Composed of 1 facet, the 1st facet The equations of the planes are: ； In the above formula, Let be the unit outward normal vector of the plane. It is a constant; If the following equation is satisfied, it means that the particle's body center is located inside a certain finite element: ； When the spatial mapping mode is an area mapping mode, the intersection of the discrete element particle and each finite element element needs to be calculated. If the intersection is not empty, then the discrete element particle and the finite element element form a binding relationship. Let the spatial region occupied by the discrete element particle be... The spatial region occupied by the finite element element is ;like Then the intersection volume between the discrete element particle and each finite element is defined as: ,Right now: ； In the above formula, For discrete element particles, , This is the indicator function for the finite element; when the intersection is empty, the integral value is 0.

4. The method for converting powder into a continuous body during the simulated pressing process as described in claim 1, characterized in that: The contact constitutive model in step (3) is any one of the linear contact model, Hertz contact model and DPC contact model; the viscoelastic-plastic constitutive model in step (3) is any one of the DPC model, MC model and Burgers model.

5. The method for converting powder into a continuous body during the simulated pressing process as described in claim 1, characterized in that, In step (4), the conversion model for converting discrete elements into finite elements can be either the volume stress conversion model or the bulk modulus conversion model. In the volumetric stress transformation model, the discrete element particle transformation factor The calculation formula is expressed as: ； In the above formula, The body stress at the start of the transformation. The body stress at the point of termination of the transformation; When the volume stress experienced by a discrete element particle is ≥ The transformation begins; when the volume stress experienced by a discrete element particle is ≥ The conversion is complete; In the bulk modulus conversion model, the discrete element particle conversion factor The calculation formula is expressed as: ； In the above formula, The equivalent bulk modulus of the region where the discrete element particles are located at the start of the transformation. The equivalent bulk modulus of the region where the discrete element particles are located at the time of termination of the transformation; When the equivalent bulk modulus of the region where a discrete element particle is located is ≥ The transformation begins; when the equivalent bulk modulus of the region where a discrete element particle is located is ≥ The conversion is complete.

6. The method for converting powder into a continuous body during the simulated pressing process as described in claim 1, characterized in that, The discrete element calculation in step (7) includes the calculation of the contact force of discrete element particles, the force calculation of discrete element particles, and the motion calculation of discrete element particles. The specific process for calculating the contact force of the discrete element particles is as follows: The calculation of interparticle contact force includes two parts: normal force calculation and tangential friction calculation. The formula for calculating the normal force when using a linear contact model is as follows: ； In the above formula, , These are the normal contact stiffness and normal damping coefficient. , For embedding quantity and embedding rate; The tangential force calculation uses an incremental-update mode, with each time step... First, calculate the test tangential force, then correct it according to Coulomb's law of friction; the formula for calculating the test stress is: ； In the above formula, To test the tangential force at the current time step, The tangential force of the previous step. For tangential stiffness, This represents the increment of the tangential displacement difference; Next, calculate the maximum static friction force at the current time step, which is: ； In the above formula, For the maximum static friction force, The coefficient of friction; Next, a correction for the tangential force is performed, calculated using the following formula: ； The specific process for calculating the force on the discrete element particles is as follows: The formulas for calculating the resultant force and resultant torque on discrete element particles are as follows: ； In the above formula, and The net force and net torque acting on the particle are shown in order, respectively. and The forces applied to the particles are gravity and torque, respectively. and These are, in order, the contact force and the connecting torque between particles; This refers to the contact force between the particle and the rigid surface. This is the position vector of the point of contact between particles relative to the center of mass of the particle; is the position vector of the point of action between the particle and the rigid surface relative to the particle's center of mass; The motion calculation of the discrete element particles includes particle translation calculation and particle rotation calculation: The formula for calculating particle translational motion is: ； In the above formula, , , The vectors are acceleration, velocity, and displacement, respectively. For particle mass, This is the local translational damping coefficient; The formula for calculating particle rotation is: ； In the above formula, , , The values ​​are, in order, rotational angular acceleration, angular velocity, and angular displacement; This is the local rotational damping coefficient; It is the moment of inertia; The quantity is characterized by a quaternion representing the rotational speed.

7. The method for converting powder into a continuous body during the simulated pressing process as described in claim 1, characterized in that, The finite element calculation in step (7) includes the calculation of the deformation and force of the finite element elements and the calculation of the force and motion of the finite element nodes. The calculation process for the deformation and stress of the finite element is as follows: Calculate a Gaussian point in a given element based on the nodal displacement increment. The strain increment is calculated using the following formula: ； In the above formula, Let be the strain increment matrix at a certain Gaussian point in a certain element. The strain matrix at the Gaussian point. This is the element strain increment matrix. Calculate the stress increment at the current moment based on the strain increment. The calculation formula is: ； In the above formula, The material matrix is ​​determined by the constitutive model used by the solid element. Calculate the total stress based on the stress increment at the current moment. The calculation formula is: ； In the above formula, This represents the total stress from the previous time step; The calculation process for the force-induced motion of the finite element nodes is as follows: First, calculate the resultant force at each node of the finite element element. The calculation formula is as follows: ； In the above formula, The net external force acting on a finite element node. The force generated by the deformation of the finite element on this node. The external force applied to this node; Then, based on the resultant force of each node Calculate the acceleration of the nodes; If the node's current velocity is greater than zero, then the acceleration... ; If the node's current velocity is less than or equal to zero, then the acceleration... ;in, This is the local damping coefficient. For node quality; Next, the velocity is updated based on the acceleration, and then the displacement increment and displacement of this step are updated based on the velocity. The calculation formula is as follows: ； In the above formula, The nodal acceleration vector; The node velocity vector; This is the vector of nodal displacement increments; The nodal displacement vector; This is for calculating the time step.

8. The method for converting powder into a continuous body during the simulated pressing process as described in claim 1, characterized in that: The quasi-static velocity load can be applied to the upper pressure plate or simultaneously to both the upper and lower pressure plates. The range of the quasi-static load is set to 10. -7 ~10 -9 m / step.

9. The method for converting powder into a continuous body during the simulated pressing process as described in claim 1, characterized in that, The formula for calculating the transformation factor of the finite element element is as follows: ； In the above formula, The total number of discrete element particles that are bound to a certain finite element element. For the first The volume of a discrete element particle. For the first The transformation factor of each discrete element particle.

10. The method for converting powder into a continuous body during the simulated pressing process as described in claim 1, characterized in that, The formula for the overall force in the calculation model is: ； In the above formula, For a certain finite element node The combined force on the area These are the nodal forces obtained directly from finite element analysis. The nodal forces are directly calculated from discrete element particles. Let be the weight coefficient of the i-th discrete element particle at the j-th finite element node; calculate First, the nodal force matrix of each node in the finite element is calculated based on the stress at each Gaussian point in the finite element. The calculation formula is: ； In the above formula, Gaussian point The strain matrix, This represents the total stress at the current time step. and Gaussian point The weight coefficients and the values ​​of the Jacobian determinant; Next, in the nodal force matrix Find the node The corresponding nodal force vector is calculated using the following formula: ； Particles The contact force exerted on it by other particles The calculation formula is: ； In the above formula, To be with particles The number of other particles in contact, For the first Individual and Particle Contacting particles act on particles The force on top.

11. The method for converting powder into a continuous body during the simulated pressing process as described in claim 1, characterized in that: The method for converting powder into a continuum can be used for two-dimensional or three-dimensional problem simulation. When performing two-dimensional simulation using the numerical method, the discrete element particles are circles, ellipses, or polygons, and the finite element units are triangles or quadrilaterals. When performing three-dimensional simulation using the numerical method, the discrete element particles are spheres, ellipsoids, or polyhedra, and the finite element units are tetrahedrons, pyramids, triangular prisms, or hexahedrons.

Images