A cold briquetting elastic stress distribution prediction method based on finite element method

Abstract

The application relates to a cold-pressed briquette elastic stress distribution prediction method based on a finite element method, and belongs to the technical field of material structure mechanics simulation. The method comprises the following steps: obtaining a microstructure image of a cold-pressed briquette, performing image segmentation to identify three-phase regions of a mineral powder particle phase, a binder film phase and a pore phase; extracting a phase boundary from the microstructure image to generate a finite element grid; giving material properties to each phase unit; adopting diffusion interface modeling at the interface between the binder and the particle to smoothly transition the elastic modulus at the interface between the binder and the particle; introducing mesoscopic parameters to perform degradation correction on the yield strength; applying boundary conditions to solve linear elastic finite element equations to obtain stress tensors of nodes; and based on the yield strength degradation correction result and the stress tensors of the nodes, calculating and outputting parameters such as a von Mises equivalent stress, a stress concentration coefficient and a high-risk area volume fraction. The application can avoid artificial simplification and make the stress distribution more physically reasonable.

Description

Technical Field

[0001] This invention belongs to the field of material structure mechanics simulation technology, specifically involving a method for predicting the elastic stress distribution of cold-pressed agglomerates based on the finite element method (FEM). Background Technology

[0002] Cold-pressed agglomerates (made by pressing iron ore powder with sodium silicate binder) contain a complex multiphase structure: ore particles (high elastic modulus, approximately 50–70 GPa), a binder film (medium modulus, approximately 5–15 GPa), and pores (modulus approaching 0). Under external compressive loads, due to the significant differences in the material properties of each phase and the discontinuities in geometry (such as particle edges, abrupt changes in binder film thickness, and sharp edges of pores), severe stress concentrations occur in localized areas. These high-stress zones often become the starting points for microcrack initiation and propagation, ultimately leading to the macroscopic failure of the agglomerates.

[0003] Currently, the main methods for analyzing the internal stress distribution of cold-pressed blocks fall into two categories: experimental and numerical simulation. Experimental methods: The strain distribution on the surface of the mass can be obtained by compression test combined with digital image correlation (DIC) technology, but the internal three-dimensional stress field cannot be obtained; moreover, the experiment is costly and time-consuming, making it difficult to conduct a large number of parameterized studies.

[0004] Numerical simulation: The Discrete Element Method (DEM) can simulate the contact and bond breaking between particles, but its output is a discrete force chain distribution, making it difficult to directly obtain the stress field in the sense of a continuous medium. The Finite Element Method (FEM), as a mature tool in continuum mechanics, can calculate detailed stress distributions based on microstructural models.

[0005] However, existing FEM applications have the following shortcomings: Obtaining microscopic geometry is challenging: generating meshes directly from Micro-CT images of real samples involves enormous computational costs and is difficult to perform parametric studies. While existing simulation methods (such as DEM) can provide particle information, they cannot directly generate continuous meshes suitable for FEM.

[0006] Simplified material property assignment: The adhesive film is usually regarded as a homogeneous material, ignoring the change of its elastic modulus with the degree of polymerization (Q4 fraction) and the gradual change of properties in the transition zone of the adhesive-particle interface.

[0007] The influence of mesoscopic defects was not considered: mesoscopic parameters such as crack density and porosity were not used to reduce local strength, which may lead to the predicted stress concentration being too conservative or deviating from reality.

[0008] Lack of coupling with multi-scale framework: It is impossible to directly map the degree of polymerization information output at the atomic scale to the elastic modulus of the binder, nor can it use the crack density output at the mesoscale as the initial defect field input.

[0009] Therefore, there is an urgent need to develop a method that can seamlessly integrate with atomic / mesoscopic multiscale simulations, accurately identify hazardous areas, and provide reliable mechanical basis for the formulation design and process optimization of cold-pressed agglomerates. Summary of the Invention

[0010] In view of this, the purpose of this invention is to provide a method for predicting the elastic stress distribution of cold-pressed agglomerates based on the finite element method, so as to solve the problems existing in the current application of the finite element method, such as difficulty in obtaining microscopic geometry, simplification of material property assignment, failure to consider the influence of mesoscopic defects, and lack of coupling with multi-scale frameworks.

[0011] This method automatically generates high-quality three-phase finite element meshes based on binarized images output by cellular automata (CA) or phase-field method (PF); and calculates the average degree of aggregation and Q based on atomic scale. 4 This method fractionally calibrates the elastic modulus of the binder film, achieving a quantitative mapping from microscopic chemical structure to macroscopic elastic properties. A diffusion interface model is employed at the binder-particle interface to ensure a smooth modulus transition and avoid spurious stress singularities. Crack density output at the mesoscale is introduced to downgrade the yield strength and predict potential failure initiation locations. Parameters such as von Mises stress contour maps, maximum / minimum principal stresses, stress concentration factors, and volume fraction of high-risk areas are output. This method can seamlessly integrate with atomic / mesoscale multi-scale simulations, accurately identifying critical areas and providing reliable mechanical basis for the formulation design and process optimization of cold-pressed agglomerates.

[0012] To achieve the above objectives, the present invention provides the following technical solution: A method for predicting the elastic stress distribution of cold-pressed agglomerates based on the finite element method, comprising: 1. Obtain microstructure images of cold-pressed lumps.

[0013] The microstructure image can be a binarized image of the cross section of a cold-pressed agglomerate output from a cellular automata (CA) model, or an order parameter cloud map output from a phase field (PF) model.

[0014] The image should clearly distinguish three phases: the mineral particle phase (grayscale value 2), the binder film phase (grayscale value 1), and the porous phase (grayscale value 0). The image resolution should be at least 600×600 pixels, and the physical size of each pixel should be consistent with the mesh size of the mesoscopic simulation (typically 1~4 μm).

[0015] 2. Perform image segmentation on the microstructure image to identify the three-phase regions: mineral powder particle phase, binder film phase, and porous phase.

[0016] The Otsu threshold segmentation algorithm is used to determine the optimal threshold for distinguishing the three phases. If it is a sequence parameter output by the PF... Cloud map, then It is a silicon-rich phase (binder). It is a silicon-poor phase (pores), and the area between the two is the interface transition zone.

[0017] 3. Extract phase boundaries from the microstructure image to generate an unstructured finite element mesh.

[0018] Phase boundaries are extracted from the binarized image using the Marching Squares algorithm (2D) or the Marching Cubes algorithm (3D), and unstructured triangular (2D) or tetrahedral (3D) meshes are generated based on the extracted phase boundaries. In the generated mesh, the smallest element size is less than 1 / 3 of the adhesive film thickness (typically ≤2 μm) to ensure the analytical accuracy of the stress gradient.

[0019] 4. Assign material properties to each phase unit.

[0020] For the mineral powder particle phase, its elastic modulus (Determined based on the type of iron ore, e.g., 65 GPa for hematite), Poisson's ratio .

[0021] Porous phase, elastic modulus (Taking the minimum value to avoid numerical oddities), Poisson's ratio .

[0022] The binder film phase, with its elastic modulus calculated from the atomic scale, represents the average degree of polymerization. and Q 4 Fraction Calibration: (1) in, This is the reference modulus at low polymerization degree (approximately 1~5 GPa for sodium silicate gel systems). The coupling coefficient (experimentally calibrated, typical value 2.5~10). Poisson's ratio of the binder film phase. .

[0023] 5. A diffusion interface model is used at the interface between the binder film phase and the mineral powder particle phase to ensure a smooth transition of the elastic modulus at the interface between the binder film phase and the mineral powder particle phase.

[0024] To avoid spurious stress singularities caused by abrupt changes in material properties, a smooth transition is employed at the interface between the binder film phase and the mineral powder particle phase. For each finite element mesh node located within the transition region at the binder film phase-mineral powder particle phase interface, the distance from each node to its nearest mineral powder particle phase boundary is calculated. And based on this distance, the elastic modulus at each node is determined: (2) (3) in, Let be the elastic modulus at the node. The width of the interface is 1 to 2 times the unit size.

[0025] A smooth transition ensures a continuous and smooth change in the elastic modulus from the binder film phase to the mineral powder particle phase, making the stress distribution more in line with physical laws.

[0026] 6. Consider the strength degradation caused by mesoscopic defects.

[0027] Obtain mesoscopic parameters—crack density—from CA or PF models. (Dimensionless) Yield strength of binder film phase material Perform a downgrade correction: (4) in, This is the yield sensitivity coefficient (2~5 for brittle materials).

[0028] The downgrading of yield strength reflects the physical fact that the strength of the concentrated microcrack area is reduced, making the prediction of high-risk areas more accurate.

[0029] 7. Apply boundary conditions, solve the linear elastic finite element equation, and obtain the stress tensor of each node.

[0030] Boundary conditions are applied to the bottom, top, and side surfaces of the generated unstructured finite element mesh.

[0031] Among them, the bottom surface ( The compression direction is Apply a fixed displacement constraint to the shaft: , , , The nodes are respectively , , Displacement components along the three coordinate axes.

[0032] Top surface ( , Apply a uniform vertical displacement (to the height of the mass) ,in Take 0.1% of the mass height to ensure that the deformation is within the linear elastic range.

[0033] The sides are configured as free boundaries.

[0034] Using linear elastic constitutive equations Solve for, where The elastic stiffness matrix is ​​uniquely determined by Young's modulus and Poisson's ratio. For the stress tensor (a second-order symmetric tensor, components) express Directional plane Directional stress). For the strain tensor (a second-order symmetric tensor, components) , for Directional displacement).

[0035] 8. Calculate and output key stress parameters, including von Mises equivalent stress distribution cloud map, maximum principal stress, minimum principal stress, stress concentration factor, volume average stress, stress distribution histogram, and volume fraction of high-risk area.

[0036] von Mises equivalent stress : (5) In the two-dimensional plane strain problem, it simplifies to: (6) in, , They are respectively , Normal stress components in the direction, for In-plane shear stress components; for Normal stress components in the direction; , They are respectively flat, In-plane shear stress components .

[0037] Volume-average von Mises equivalent stress : (7) in, The total number of grid cells. For the first von Mises equivalent stress per unit For the first Volume of each unit For the entire computational domain, This is the total volume.

[0038] Maximum principal stress and minimum principal stress By solving the characteristic equation The three principal stresses (eigenvalues) are obtained and arranged in descending order of their algebraic values ​​as follows: .in, These are the eigenvalues ​​(i.e., principal stresses). It is a second-order unit tensor. Let be the stress tensor. The maximum principal stress (the largest algebraic value, a positive value indicates tensile stress) indicates the location where tensile failure is most likely to occur; It represents the minimum principal stress (the smallest algebraic value; a negative value indicates compressive stress).

[0039] Stress concentration factor : (8) in, This is the nominal stress (far-field average stress). This represents the local stress peak. This reflects the degree of stress amplification caused by geometric discontinuities.

[0040] Volume fraction of high-risk areas : (9) in, This is an indicator function (1 if the condition is true, 0 otherwise). For the total volume, The entire computational domain is represented. This metric is used to quantitatively assess the potential extent of plastic yielding or failure within the material.

[0041] Stress distribution histogram: Divides the stress range into several intervals and counts the proportion of unit volume or number falling into each interval, reflecting the degree of stress non-uniformity.

[0042] 9. Marking and outputting localized stress concentration areas: Marking The system identifies regions and outputs their coordinates and corresponding principal stress directions to guide experimental verification.

[0043] Furthermore, in equation (1), and Different Q values ​​can be measured through nanoindentation experiments. 4 The modulus of the fractional samples is determined by linear regression.

[0044] Furthermore, the interface width Use twice the unit size to ensure a smooth transition between at least two units.

[0045] Furthermore, yield strength Yield sensitivity coefficient can be obtained through uniaxial tensile testing. The strength decay curves of specimens with different crack densities can be determined by fitting them.

[0046] Furthermore, when generating the mesh based on the extracted phase boundaries, a structured mesh can also be generated. If a structured mesh scheme is adopted, quadrilateral elements (2D) or hexahedral elements (3D) are preferred to reduce numerical artifacts. The smallest element size in the structured mesh is less than 1 / 3 of the adhesive film thickness.

[0047] The beneficial effects of this invention are as follows: (1) Microscopic geometric reality: The mesh is constructed directly from the image output by the mesoscopic simulation (CA / PF), which ensures that the distribution of particles, binders and pores is consistent with the simulation results and avoids artificial simplification.

[0048] (2) Material properties have a physical basis: the phase modulus of the binder film is derived from Q at the atomic scale. 4 Fractional calibration established a quantitative correlation between chemical structure and macroscopic elastic properties.

[0049] (3) Eliminating spurious stress singularities: The diffusion interface method (tanh smoothing function) is adopted to avoid abrupt changes in modulus at phase boundaries, making the stress distribution more physically reasonable.

[0050] (4) Consider the impact of defects on strength: introduce crack density to reduce yield strength, making the prediction of high-risk areas more accurate and in line with the principles of damage mechanics.

[0051] (5) Abundant output parameters: It not only provides stress cloud diagrams, but also quantitatively calculates practical engineering indicators such as stress concentration factor and volume fraction of high-risk area, which can be used to compare the advantages and disadvantages of different process schemes.

[0052] (6) Seamless integration with multi-scale frameworks: The input comes from atomic / mesoscopic simulations, and the output can be directly used to guide process optimization, forming a closed loop.

[0053] Other advantages, objectives, and features of the invention will be set forth in part in the description which follows, and in part will be apparent to those skilled in the art from the following examination, or may be learned from practice of the invention. The objectives and other advantages of the invention can be realized and obtained through the following description. Attached Figure Description

[0054] To make the objectives, technical solutions, and advantages of the present invention clearer, the preferred embodiments of the present invention will be described in detail below with reference to the accompanying drawings, wherein: Figure 1 This is a schematic diagram of the process for predicting the elastic stress distribution of cold-pressed lumps based on the finite element method in an embodiment of the present invention.

[0055] Figure 2 Here is an example of a von Mises stress distribution cloud map, where (a) shows the overall distribution, indicating that the high-stress areas are mainly located near the surface of the agglomerates; and (b) shows a magnified view of a local area, indicating stress concentration at the corner of the binder-particle contact (red area).

[0056] Figure 3 This is a histogram of stress distribution and a schematic diagram of the volume fraction of high-risk areas. (a) is the stress distribution histogram, with the horizontal axis representing the stress value (MPa) and the vertical axis representing the percentage of elements; (b) is a schematic diagram of the mass cross-section, with the red area marking the high-risk area. ).

[0057] Figure 4 A bar chart comparing the volume fraction of high-risk areas under different CO2 treatment conditions (0s, 15s, 30s, 60s). Detailed Implementation

[0058] The following specific examples illustrate the implementation of the present invention. Those skilled in the art can easily understand other advantages and effects of the present invention from the content disclosed in this specification. The present invention can also be implemented or applied through other different specific embodiments, and various details in this specification can be modified or changed based on different viewpoints and applications without departing from the spirit of the present invention. It should be noted that the illustrations provided in the following embodiments are only schematic representations of the basic concept of the present invention. Unless otherwise specified, the following embodiments and features can be combined with each other.

[0059] The accompanying drawings are for illustrative purposes only and are schematic diagrams, not actual pictures, and should not be construed as limiting the invention. To better illustrate the embodiments of the invention, some parts in the drawings may be omitted, enlarged, or reduced, and do not represent the actual product dimensions. It is understandable to those skilled in the art that some well-known structures and their descriptions may be omitted in the drawings.

[0060] In the accompanying drawings of the embodiments of the present invention, the same or similar reference numerals correspond to the same or similar components. In the description of the present invention, it should be understood that if terms such as "upper," "lower," "left," "right," "front," and "rear" indicate the orientation or positional relationship based on the orientation or positional relationship shown in the drawings, they are only for the convenience of describing the present invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation. Therefore, the terms used to describe positional relationships in the drawings are only for illustrative purposes and should not be construed as limiting the present invention. For those skilled in the art, the specific meaning of the above terms can be understood according to the specific circumstances.

[0061] Example 1: 2D plane strain analysis based on CA microstructure images Data Preparation: A 600×600 pixel binarized image was exported from the output of the CA model (CNT-enhanced system with CO2 processing), with a physical size of 1.2 mm × 1.2 mm (2 μm per pixel). In the image, gray represents mineral powder particles (65%), blue represents binder film (25%), and white represents pores (10%). After extracting phase boundaries, a triangular mesh was generated with a total of approximately 150,000 elements and a minimum element size of 1.5 μm.

[0062] Material parameters: Mineral powder particle phase: , .

[0063] Porous phase: , .

[0064] Binder thin film phase: output from atomic scale ( , ),Pick , ,have to , .

[0065] Diffuse Interface: Sets the interface width (2 units). Calculate the distance from each node in the transition region of the interface to the nearest mineral powder phase boundary, substitute it into the tanh function, and obtain a smooth modulus field.

[0066] Strength degradation: crack density Yield strength (Taken from nanoindentation), take ,but (Downgraded by approximately 0.75%).

[0067] Boundary and Solution: Bottom surface fixed, displacement applied to top surface (0.1% height). Solved using ABAQUS Standard, linear statics analysis.

[0068] Output results: von Mises stress cloud diagram ( Figure 2 (a) shows the high-stress area (red, The main concentrations are 1 to 2 particle depths below the surface of the agglomerate, especially at the corners where the adhesive film contacts the particles.

[0069] Local magnification ( Figure 2 (b) shows that the maximum local stress occurs at the sharp corner of a particle. Nominal stress (Estimated from total load / cross-sectional area), then the stress concentration factor .

[0070] Volume fraction of high-risk areas (That is, approximately 12.3% of the volumetric internal stress exceeds 0.8 × 11.91 = 9.53 MPa).

[0071] The stress distribution histogram shows that the stress in most elements is concentrated in the range of 2 to 6 MPa, with the maximum stress value reaching 58 MPa (corresponding to a very small number of elements in the histogram).

[0072] Example 2: Comparison under different CO2 treatment conditions Using the same microstructure image from Example 1, but modifying the crack density parameter to simulate CO2-free treatment ( ) and excessive CO2 treatment ( The three cases are calculated separately, as shown in Table 1: Table 1

[0073] The results show that CO2 treatment reduces crack density, effectively increases yield strength (slightly), and reduces the volume fraction of the high-risk zone. However, overtreatment does not significantly improve the crack density, indicating that 30 s of CO2 treatment is sufficient.

[0074] Example 3: The Influence of Material Type (CNT vs SiO2) Keeping all other conditions the same as in Example 1, only changing the phase modulus of the adhesive film: SiO2 system (correspond CNT system (correspond ).

[0075] The calculated maximum von Mises stress of the CNT system is slightly higher (62 MPa), but its yield strength is also higher. (Relatively high), but the actual volume fraction of high-risk areas is even lower (CNT is 10.8%, SiO2 is 13.5%), meaning that CNT-reinforced agglomerates have higher crack resistance.

[0076] Example 4: Mesh Sensitivity Verification To verify the influence of mesh size, based on Example 1, the same microstructure was divided using unit sizes of 1.5 μm, 2.0 μm, and 3.0 μm, respectively.

[0077] Calculated The values ​​of 7.25, 7.18, and 6.95 respectively indicate that the results converge when the cell size is less than 1 / 3 of the adhesive film thickness (approximately 8 μm). Therefore, a minimum cell size of no more than 2 μm is recommended.

[0078] Example 5: Verification by comparison with experiments Micro-CT images were obtained from actual compressed samples (SiO2 system, optimal CO2 treatment), and a 3D finite element model was reconstructed. The locations of high-stress zones predicted by simulation were compared with the locations of surface crack initiation observed after compression experiments.

[0079] The comparison results show that surface cracks first appear in the high-stress corner area predicted by the simulation (with an error within 2 particle diameters), verifying the reliability of the present invention.

[0080] Example 6: Application of Parametric Research Using the method of the present invention to study different sodium silicate moduli Stress analysis was performed on the aggregates with moduli of (2.15, 2.6, 3.0). The elastic modulus of the adhesive film increased with increasing modulus. The maximum von Mises stress increased (from 7 GPa to 11 GPa), but... The increase was even greater, leading to The percentage decreased from 15% to 9%. Therefore, increasing the modulus helps reduce the volume of high-risk areas and enhance crack resistance.

[0081] In summary, this invention provides a method for predicting the elastic stress distribution of cold-pressed blocks based on the finite element method. It can provide stress cloud maps, quantitatively calculate practical engineering indicators such as stress concentration factor and volume fraction of high-risk areas, and can be used to compare the advantages and disadvantages of different process schemes and guide process optimization.

[0082] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of the present invention without departing from the spirit and scope of the present invention, and all such modifications or substitutions should be covered within the scope of the claims of the present invention.

Claims

1. A method for predicting the elastic stress distribution of cold-pressed agglomerates based on the finite element method, characterized in that, The method includes: Obtain microstructure images of cold-pressed lumps; Image segmentation is performed on the microstructure image to identify the three-phase regions of mineral powder particles, binder film phase, and porous phase; Phase boundaries are extracted from microstructure images to generate finite element meshes; Assign material properties to each phase unit; A diffusion interface model is used at the interface between the binder and the particles to ensure a smooth transition in the elastic modulus at the interface. A mesoscopic parameter is introduced to downgrade the yield strength; By applying boundary conditions, the linear elastic finite element equations are solved to obtain the stress tensors at each node. Based on the yield strength downgrade correction results and the stress tensor of each node, key stress parameters are calculated and output, including the von Mises equivalent stress distribution cloud map, maximum principal stress, minimum principal stress, stress concentration factor, volume average von Mises equivalent stress, stress distribution histogram, and volume fraction of high-risk area.

2. The method according to claim 1, characterized in that, The generated finite element mesh can be a structured finite element mesh or an unstructured finite element mesh; in the finite element mesh, the minimum element size is set to be less than 1 / 3 of the adhesive film thickness.

3. The method according to claim 1, characterized in that, Material properties are assigned to each phase unit, including elastic modulus and Poisson's ratio for the mineral powder particle phase, binder film phase and porous phase.

4. The method according to claim 1, characterized in that, Diffusion interface modeling at the interface between the binder film phase and the mineral powder particle phase includes calculating the distance from each node to its nearest mineral powder particle phase boundary for each node within the interface region, and determining the elastic modulus at each node based on this distance.

5. The method according to claim 4, characterized in that, The elastic modulus at each node is determined by the following formula: In the formula, Let be the elastic modulus at the node. The elastic modulus of the adhesive film phase. This represents the elastic modulus of the mineral powder particle phase. is the distance from the node to its nearest mineral particle phase boundary. This represents the width of the interface.

6. The method according to claim 5, characterized in that, Interface width Use 1 to 2 times the unit size.

7. The method according to claim 1, characterized in that, The degradation correction of yield strength by introducing mesoscopic parameters includes: introducing crack density. Yield strength of binder film phase material Perform a downgrade correction: In the formula, This represents the yield sensitivity coefficient.

8. The method according to claim 1, characterized in that, Applying boundary conditions includes applying boundary conditions to the bottom, top, and side surfaces of the unstructured finite element mesh as a whole; wherein, a fixed displacement constraint is applied to the bottom surface; a uniform vertical displacement is applied to the top surface; and the side surfaces are configured as free boundaries.

9. The method according to claim 8, characterized in that, Solving the linear elastic finite element equations using the linear elastic constitutive equations To solve, where, Here is the elastic stiffness matrix. For stress tensor, For strain tensor.

10. The method according to claim 1, characterized in that, von Mises equivalent stress is calculated using the following formula: In the formula, , , They are respectively , , Normal stress components in the direction, , , They are respectively flat, flat, In-plane shear stress components; The stress concentration factor is calculated using the following formula: In the formula, This is a local peak value. Nominal stress; The high-risk volume fraction refers to a von Mises equivalent stress exceeding 0.8 times the material yield strength. The proportion of volume occupied by the region; The volume-average von Mises equivalent stress is calculated using the following formula: In the formula, The total number of grid cells. For the first von Mises equivalent stress per unit For the first Volume of each unit For the entire computational domain, This is the total volume.

Images