A cross-scale numerical simulation method based on ABAQUS explicit finite element and LIGGGHTS discrete element
Patent Information
- Application Number
- CN202610775805.6
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2026-06-01
- Publication Date
- 2026-10-13
AI Technical Summary
离散元根据接触模型和粒间作用求解颗粒运动，无需定义宏观本构模型，但因计算量大难以应用于边界条件问题的模拟
[0029]本发明的计算无需依赖经验性的宏观本构模型，避免因本构不同导致的结果差异，解决有限元在复杂工况下本构不完善的问题；
Smart Images
Figure CN122907188A_ABST
Abstract
Description
Technical Field
[0001] This invention belongs to the field of computational mechanics and numerical analysis of geotechnical engineering, and particularly relates to a cross-scale numerical simulation method based on ABAQUS explicit finite element and LIGGGHTS discrete element. Background Technology
[0002] Geotechnical materials exhibit significant discreteness and multi-scale characteristics. The finite element method (FEM) requires solving boundary condition problems based on material constitutive relations. The discrete element method (DIM) solves particle motion based on contact models and interparticle interactions, without requiring the definition of a macroscopic constitutive model, but its computational complexity makes it difficult to apply to the simulation of boundary condition problems. Combining ABAQUS and LIGGGHTS for multi-scale numerical computation can fully leverage the advantages of both methods and has good development potential.
[0003] Explicit analysis is a dynamic solution process that uses the central difference method for time integration. Displacement and stress updates depend only on the state at the previous time step, avoiding the construction and inversion of the global tangent stiffness matrix, resulting in extremely high computational efficiency in large-scale calculations. Explicit solution using the central difference integration method eliminates the need to assemble the global stiffness matrix, leading to a high degree of computational decoupling between the ABAQUS elements and the embedded LIGGGHTS particle swarm. Furthermore, the localization characteristics of the explicit algorithm are well-suited for MPI parallelism, significantly improving computational efficiency. Therefore, this study develops a multi-scale numerical simulation method based on explicit finite element method (FEM) and discrete element method (DIM). The FEM uses ABAQUS, a fully functional finite element program with powerful secondary development capabilities. The DEM uses the widely adopted open-source LIGGGHTS software, which offers advantages such as algorithmic transparency, robustness of explicit integration, and flexibility in multi-scale expansion. Multi-scale coupling between ABAQUS and LIGGGHTS is implemented using Python, thus addressing the issues of constitutive incompleteness in ABAQUS under complex conditions and computational limitations in LIGGGHTS. Summary of the Invention
[0004] Objective: Traditional ABAQUS requires pre-defined macroscopic constitutive relations (such as elastoplastic models); while LIGGGHTS can accurately simulate the discrete characteristics of materials, its computational cost increases non-linearly with the number of particles, limiting its engineering applicability. This invention proposes a cross-scale numerical simulation method based on ABAQUS explicit finite element method and LIGGGHTS discrete element method. This method allows the equivalent mechanical response of certain regions in ABAQUS to be obtained through LIGGGHTS simulation of the representative volume element (RVE) corresponding to that region, rather than relying on empirical macroscopic constitutive models. This significantly improves overall computational efficiency and makes the method applicable to solving engineering-scale problems.
[0005] Technical solution: The present invention provides a cross-scale numerical simulation method based on ABAQUS explicit finite element and LIGGGHTS discrete element method, comprising the following steps:
[0006] S1. Establish a finite element numerical model in the finite element software ABAQUS based on the information of the research object; divide the model according to the mesh generation scheme to achieve a balance between the accuracy of the results and the required computing resources; reasonably set the contact model and model parameters, analysis steps, constraints and load conditions according to the characteristics of the research object and the model features; set the soil and rock material properties of the Gaussian points in the cross-scale calculation region in ABAQUS.
[0007] S2. Specify a custom VUMAT path in the Job module of ABAQUS, trigger the subroutine interface in each calculation step of the finite element software ABAQUS, and associate the Gaussian integration point of each cross-scale element in the model with the discrete element software LIGGGHTS through Gaussian point numbering.
[0008] S3. Generate a set of LIGGGHTS particles in LIGGGHTS with the same number of Gaussian points. Set the consolidation pressure of the LIGGGHTS unit according to the stress conditions in ABAQUS. Achieve the target consolidation pressure in LIGGGHTS through servo control for subsequent calling and loading.
[0009] S4. Based on the Gaussian point number in ABAQUS, write the corresponding Gaussian point information into a .txt file in the VUMAT subroutine, read the strain information in the .txt file through Python, and convert it into the boundary control information of the LIGGGHTS model, which is then applied in LIGGGHTS as the boundary conditions of the corresponding Gaussian point in the LIGGGHTS model.
[0010] S5. After confirming that the actual strain increment in LIGGGHTS should be equal to the strain previously set in ABAQUS, calculate the stress state of the element in combination with the boundary force and the contact point location, and write it into a .txt file with a specific naming rule. The stress response of the Gaussian point representative element is fed back to the ABAQUS main program by the VUMAT subroutine according to the file name.
[0011] S6. Following the central difference algorithm, perform the next cycle of calculation based on the updated Gaussian point information from the current calculation results. Align the two in physical time by controlling the speed of boundary loading. After the macroscopic calculation step is completed, cross-scale data interaction is completed through a txt file. The calculation is repeated until the preset load conditions are met.
[0012] Furthermore, in step S1, a user-defined material is selected in ABAQUS to set the material properties of the Gaussian points in the cross-scale computational domain.
[0013] Furthermore, in step S3, the LIGGGHTS particle set includes particle gradation, contact models between particles and between particles and boundaries, and contact parameters such as elastic modulus, Poisson's ratio, friction coefficient, and damping. Loads are applied and displacement loading or load loading is achieved through boundary constraint conditions.
[0014] Further, step S4 specifically involves: obtaining the Gaussian point number in ABAQUS through VUMAT, and naming the corresponding data file (txt file) with the number; writing the incremental time, total time, and strain increment matrix of the current incremental step into the corresponding .txt file named according to the Gaussian point using VUMAT.
[0015] After writing is complete, VUMAT triggers the execution of an external Python script via a command. The Python script extracts the strain increment from the .txt file. Data; based on the time step required for the current computation step in ABAQUS. Combined with the number of time steps required to achieve this macroscopic step size as set in LIGGGHTS Calculate the time step size for each small step in LIGGGHTS ,Right now:
[0016]
[0017] By combining the strain increment information of the current calculation step, the speed of the current calculation step is calculated, so as to ensure that the LIGGGHTS program is consistent with the macroscopic increment state of ABAQUS after it is completed.
[0018] Using the formula:
[0019]
[0020] Calculate the corresponding deformation gradient tensor F, where I is the unit tensor. For strain; the rotation axis is calculated by the cross product of the wall normal vectors before and after deformation, and the angle between the two normal vectors is calculated by the magnitude of the normal vector, and then the rotation speed is calculated; then the translation speed of the center point of the LIGGGHTS boundary surface and the rotation angular velocity of the boundary surface are determined within the corresponding time step;
[0021] Subsequently, the Python script sends commands to the LIGGGHTS solver via an interface, applying the calculated kinematic parameters to the boundary. Once set up, LIGGGHTS runs for the specified number of steps, completing the LIGGGHTS calculations required for one incremental step in ABAQUS.
[0022] Further, step S5 specifically involves: confirming that the actual strain increment in LIGGGHTS should be equal to the strain previously set in ABAQUS; calculating the boundary stress through the force and area of the boundary wall; obtaining each stress component in Python through the LIGGGHTS interface; and writing it to the .txt file of the corresponding Gaussian point.
[0023] In the VUMAT subroutine, the corresponding txt file is located by looping through the element number, and the data is read into the stress-related variables. Then, the read values are assigned to the stressNew function to update the stress state of the Gaussian point in ABAQUS, which serves as the constitutive response input for macroscopic finite element analysis.
[0024] Furthermore, step S6 specifically involves: in the subsequent time step iteration phase, through explicit integration of the central difference, the motion information of the next instantaneous step is directly extrapolated and calculated based on the current acceleration, velocity, and displacement state in each incremental step; based on the VUMAT subroutine interface, real-time and efficient exchange of macro and micro information is realized; the above process is executed repeatedly until the preset external loading termination criterion is reached.
[0025] The present invention also discloses a computer device, including a memory, a processor, and a computer program stored in the memory, wherein the processor executes the computer program to implement the steps of the method of the present invention.
[0026] The present invention also discloses a computer-readable storage medium having a computer program / instructions stored thereon, which, when executed by a processor, implements the steps of the method of the present invention.
[0027] The present invention also discloses a computer program product, including a computer program / instructions that, when executed by a processor, implement the steps of the method of the present invention.
[0028] Beneficial effects: Compared with the prior art, the present invention has the following significant advantages:
[0029] The calculations of this invention do not rely on empirical macroscopic constitutive models, thus avoiding differences in results due to different constitutive models and solving the problem of constitutive imperfections in finite element methods under complex working conditions.
[0030] An explicit algorithm is adopted, which avoids the construction and inversion of the global tangent stiffness matrix, improves the computational efficiency of large-scale operations, and solves the problem of discrete element computation being limited by the amount of computation.
[0031] The asynchronous time-stepping algorithm only performs data exchange at key nodes, reducing the exchange frequency and thus improving computational efficiency. Moreover, compared with the results of fine simulation of the entire discrete element method, the relative error of the global macroscopic stress-strain curve is less than 5%.
[0032] By leveraging the powerful features and secondary development capabilities of ABAQUS, as well as the transparency, robustness of explicit integration, and high-performance parallel computing advantages of the LIGGGHTS open-source algorithm, both R&D efficiency and computational scale can be balanced. Attached Figure Description
[0033] Figure 1 This is a schematic diagram of the ABAQUS-LIGGGHTS multi-scale calculation provided by the present invention;
[0034] Figure 2 The multi-scale calculation flowchart provided for this invention;
[0035] Figure 3 This is the axial strain-stress ratio curve in Embodiment 1 of the present invention. Detailed Implementation
[0036] The technical solution of the present invention will be further described below with reference to the accompanying drawings.
[0037] like Figure 1 , Figure 2 As shown, the present invention provides a cross-scale numerical simulation method based on ABAQUS explicit finite element and LIGGGHTS discrete element methods, comprising the following steps:
[0038] S1: Establish a numerical model in ABAQUS based on the geometric information of the research object, and perform mesh generation, and assign the "user material" to the Gaussian points corresponding to the mesh;
[0039] S2: Specify a custom VUMAT path in the User subroutine file of the Job module in ABAQUS, trigger the subroutine interface in each calculation step of ABAQUS, and associate the Gaussian integration point of each cross-scale unit in the model with LIGGGHTS through Gaussian point numbering.
[0040] S3: Generate a set of LIGGGHTS particles in LIGGGHTS with the same number of Gaussian points, specify the particle size distribution, contact model and contact parameters between particles and between particles and boundaries, load control conditions, boundary constraint conditions, set the consolidation pressure of LIGGGHTS unit according to the stress conditions in ABAQUS, and achieve the target consolidation pressure in LIGGGHTS through servo control for subsequent calling and loading;
[0041] S4: Based on the Gaussian point number in ABAQUS, write the corresponding Gaussian point information into the ".txt" file in the VUMAT subroutine, read the strain information in the ".txt" file through Python, and convert it into the boundary control information of the LIGGGHTS model, which is then applied in LIGGGHTS as the boundary conditions of the corresponding Gaussian point in the LIGGGHTS model.
[0042] S5: After confirming that the actual strain increment in LIGGGHTS should be equal to the strain previously set in ABAQUS, calculate the stress state of the element in combination with the boundary forces and contact point locations, and write it to a ".txt" file with a specific naming rule. The stress response of the Gaussian point representative element is fed back to the ABAQUS main program by the VUMAT subroutine according to the file name.
[0043] S6: Following the central difference algorithm, the next cycle of calculation is performed based on the updated Gaussian point information according to the current calculation results. The speed of boundary loading is controlled to align the two in physical time. After the macroscopic calculation step is completed, cross-scale data interaction is completed through the ".txt" file. The calculation is repeated until the preset termination criterion is reached.
[0044] Furthermore, a finite element numerical model is established in ABAQUS based on the information of the research object; a suitable mesh generation scheme is selected to divide the model so as to balance the accuracy of the results with the required computational resources; the contact model and model parameters, analysis steps, constraints and load conditions are reasonably set according to the characteristics of the research object and the model features; in ABAQUS, "user material" is set as the material property of the Gaussian points in the region that needs to be calculated across scales.
[0045] Furthermore, a custom VUMAT path is specified in the User subroutine file of the Job module in ABAQUS, and the subroutine interface is triggered in each calculation step of ABAQUS. The Gaussian integration point of each cross-scale unit in the model is associated with LIGGGHTS through Gaussian point numbering.
[0046] Furthermore, in LIGGGHTS, the particles and boundaries of the unit cell sample are generated, and the contact models and contact parameters such as elastic modulus, Poisson's ratio, friction coefficient, and damping are set between particles and between particles and boundaries, thus completing the preparation of the unit cell in LIGGGHTS. Specifically, this includes:
[0047] In LIGGGHTS, a set of particles with a specific gradation is generated. Considering the computational accuracy of LIGGGHTS and the impact of the number of Gaussian points on the computation time during the coupling process, the model side length should be set reasonably. It is recommended that the ratio of the model side length to the particle diameter should be greater than 10:1, that is, the number of particles should be at least greater than 1000.
[0048] Set up contact models and particle parameters between particles and between particles and boundaries, as well as parameters required for the contact models such as elastic modulus, Poisson's ratio, friction coefficient, and damping; to ensure computational stability, control the time step within a certain range:
[0049]
[0050] Where ρ is the density of the particulate material, r is the particle radius, G is the shear modulus, and v is Poisson's ratio.
[0051] Furthermore, rigid walls are used as constraints on the unit body to control boundary conditions. Loads are applied and displacement loading or load loading is achieved through boundary constraints. Furthermore, the consolidation pressure of the LIGGGHTS unit body is set according to the stress conditions in ABAQUS. The target consolidation pressure is achieved in LIGGGHTS through servo control, so that the initial stress states of the two are consistent.
[0052] Furthermore, based on the Gaussian point numbers in ABAQUS, the corresponding Gaussian point information is written to a ".txt" file in the VUMAT subroutine. Python reads the strain information from the ".txt" file and converts it into boundary control information for the LIGGGHTS model. This information is then applied as the boundary conditions of the corresponding Gaussian point in the LIGGGHTS model. Specifically, this includes:
[0053] Based on the Gaussian point number in ABAQUS, the VUMAT subroutine writes the element number, increment step number, total time, and strain increment component information of each Gaussian point into the corresponding ".txt" file.
[0054] To complete the loading within LIGGGHTS physical time, the current increment step duration and strain increment are used. Calculate the translation or rotation speed of the wall, and align the two in physical time by controlling the speed of boundary loading;
[0055] The strain increment information at each Gaussian point is read through the VUMAT subroutine and converted into the translational velocity of the boundary surface center point and the rotational angular velocity of the boundary surface, which are then used as the boundary conditions for the element in LIGGGHTS.
[0056] Translational velocity of the center point of the boundary surface:
[0057]
[0058]
[0059] Angular velocity of boundary surface rotation:
[0060]
[0061] In the formula, Let be the coordinate vector of the center point of the boundary before deformation. The center point of the deformed boundary surface, The vector representing the difference between the centers of the boundary surfaces before and after deformation. In response, To calculate the step size, and F and I are the normal vectors of the boundary surfaces before and after deformation, respectively. F is the deformation gradient tensor, and I is the unit tensor.
[0062] The LIGGGHTS.py file provided in LIGGGHTS extends the Python interface to call the LIGGGHTS command line for calculation. The boundary condition commands obtained from the above calculations are sent to the LIGGGHTS command line, and then the discrete element RVE operation is run to obtain the material stress-strain relationship in the current calculation step.
[0063] Furthermore, after confirming that the actual strain increment in LIGGGHTS should be equal to the strain previously set in ABAQUS, the stress state of the element is calculated by combining the boundary forces and contact point locations, and written to a ".txt" file with a specific naming convention. The stress response of the Gaussian point representative element is then fed back to the ABAQUS main program by the VUMAT subroutine according to the filename, specifically including:
[0064] After LIGGGHTS completes the translational and rotational loading of the aforementioned boundary surface, the load passes through the center point of the deformed boundary surface. The area S and the normal vector of the boundary surface before deformation. Determine the deformation gradient tensor F of the RVE model:
[0065]
[0066] Furthermore, the actual strain during the discrete element method calculation process is determined:
[0067]
[0068] In the formula: V represents the volume of the particle assembly in the discrete element method (DEM). Except for the initial state where the particle assembly is a regular hexahedron, it is generally an irregular hexahedron at other times. Its volume calculation formula is related to the effective boundary area. Distance between origin and boundary surface Related:
[0069]
[0070] Furthermore, after confirming that the actual strain increment in LIGGGHTS should be equal to the strain previously set in ABAQUS, the stress state of the unit is calculated by extracting the wall force and geometric position through variables in LIGGGHTS, and written to a ".txt" file named with Gauss point number.
[0071] Furthermore, in VUMAT, the READ statement is used to read the data into stressNew, so that the stress calculated by LIGGGHTS is returned to the ABAQUS main program as the constitutive response input for macroscopic finite element analysis, and participates in the global nodal force balance calculation and acceleration update of the next time step.
[0072] Furthermore, a Python script is used to calculate the reaction force of the LIGGGHTS wall under this state, convert it into stress and write it into a ".txt" file. VUMAT reads the stress from the ".txt" file and updates the stressNew array to complete the calculation of the current Gaussian point and the interaction of multi-scale information.
[0073] Furthermore, through central difference explicit integration, based on the VUMAT subroutine interface, after obtaining the current stress state, strain information within this macroscopic time step is calculated based on the current acceleration, velocity, and displacement state in each increment step. The physical time alignment of the two is achieved by controlling the boundary loading speed. After the macroscopic calculation step is completed, cross-scale data interaction is completed through a ".txt" file. The calculation is repeated until the preset termination criterion is reached.
[0074] Example:
[0075] According to an embodiment of the present invention, a cross-scale numerical simulation technique based on ABAQUS explicit finite element method and LIGGGHTS discrete element method is provided.
[0076] S1: Taking the calculation of a single element as an example, a finite element-discrete element coupled calculation model of a single element is established with a size of 1m×1m×1m; the initial stress is set to 100kPa. In order to balance the internal and external forces at the nodes of the model, a confining pressure of 100kPa is applied around the perimeter to balance the initial stress, and the loading strain is set to 20%; when setting the material parameters, “usermaterial” is selected as the material property of the finite element-discrete element coupled calculation model of the single element.
[0077] S2: In the ABAQUS job module, set the working path of the VUMAT subroutine, call the self-written Python program, and start a discrete element model through VUMAT according to the number of Gaussian points (1 in this example);
[0078] S3: The discrete element RVE unit cell is a hexahedron. The boundary of the RVE model is formed by the intersection of 6 faces with side lengths greater than the side length of the unit cell. This ensures that the boundary can still constrain the particles during translation and rotation. The particles are generated using the expansion method. The actual particle radius is uniformly distributed in the range of 0.4-0.8 mm. In this example, approximately 4000 particles are used. The contact model between particles is set to linear contact. The particle density is 2500 kg / m3, the Young's modulus is 6.8e8 Pa, the Poisson's ratio is 0.25, and the coefficient of friction between particles is 0.5. Consolidation is achieved by controlling the movement of the boundary wall. The velocity of the boundary wall is 0.005 m / s. Finally, the consolidation pressure of the sample is controlled to be 100 kPa to keep it consistent with the initial state in ABAQUS.
[0079] S4: Write a VUMAT subroutine in Fortran to write the element number, integration point number, increment step number, iteration number, and strain increment array elements into a file named "1.txt". Name the data file using the element number and integration point number. Run a Python script based on the VUMAT subroutine to read the information from the corresponding "1.txt" file and store the strain increment information. Based on the strain increments output in Abaqus, obtain the deformation gradient tensor F representing the volume element (RVE). According to the notation convention of Abaqus, the strain output by Abaqus needs to be negativeized. Set the number of steps required in LIGGGHTS for each calculation step in ABAQUS to 2000, and the LIGGGHTS time step to 0.000005. Calculate the translational velocity of the center point of the LIGGGHTS boundary surface and the rotational angular velocity of the boundary surface within a macroscopic calculation step. Apply these as boundary conditions to the RVE model in LIGGGHTS within this calculation step. Use the liggghts.py file provided in LIGGGHTS to extend the Python script. The interface calls the LIGGGHTS command line to perform the loading and solving within the current time step. After confirming that the actual strain increment in LIGGGHTS should be equal to the strain previously set in ABAQUS, the calculation for the current time step is completed.
[0080] S5: After the calculation in LIGGGHTS is completed, update the actual strain of the LIGGGHTS element in this macroscopic time step. After confirming that it is equal to the strain set in ABAQUS, calculate the stress state of the element according to the boundary forces and contact point positions. Store the components of the stress tensor in the corresponding "1.txt" file according to the element number and integration point number. Read the file by file name in the VUNAT subroutine and return the stress tensor calculated above to the ABAQUS model to realize the stress update at the Gaussian integration point in the finite element model.
[0081] S6: Through explicit integration with central difference, based on the VUMAT subroutine interface, after obtaining the current stress state, the strain information within each incremental step is calculated based on the current acceleration, velocity and displacement state. The physical time is aligned between the two by controlling the boundary loading speed. After the macroscopic calculation step is completed, cross-scale data interaction is completed through the ".txt" file. The calculation is repeated until the preset strain condition (ε=20%) is reached.
[0082] like Figure 3 As shown, the stress-strain curves output by ABAQUS and LIGGGHTS when using coupled calculations, as well as the results of using LIGGGHTS alone, show that the output curves are completely consistent before reaching the peak value, and are also relatively consistent even under large deformations, indicating the effectiveness of this coupling result.
[0083] The above embodiments are merely preferred embodiments of the present invention. It should be noted that those skilled in the art can make several improvements and equivalent substitutions without departing from the principle of the present invention. All such improvements and equivalent substitutions to the claims of the present invention fall within the protection scope of the present invention.
Claims
1. A cross-scale numerical simulation method based on ABAQUS explicit finite element method and LIGGGHTS discrete element method, characterized in that, Includes the following steps: S1. Establish a finite element numerical model in the finite element software ABAQUS based on the information of the research object; divide the model according to the mesh generation scheme to achieve a balance between the accuracy of the results and the required computing resources; reasonably set the contact model and model parameters, analysis steps, constraints and load conditions according to the characteristics of the research object and the model features; set the soil and rock material properties of the Gaussian points in the cross-scale calculation region in ABAQUS. S2. Specify a custom VUMAT path in the Job module of ABAQUS, trigger the subroutine interface in each calculation step of the discrete element software LIGGGHTS, and associate the Gaussian integration point of each cross-scale unit in the model with LIGGGHTS through Gaussian point numbering. S3. Generate a discrete particle set in LIGGGHTS with the same number of Gaussian points. Set the consolidation pressure of the LIGGGHTS unit according to the stress conditions in ABAQUS. Achieve the target consolidation pressure in LIGGGHTS through servo control for subsequent calling and loading. S4. Based on the Gaussian point number in ABAQUS, write the corresponding Gaussian point information into a .txt file in the VUMAT subroutine, read the strain information in the .txt file through Python, and convert it into the boundary control information of the LIGGGHTS model, which is then applied in LIGGGHTS as the boundary conditions of the corresponding Gaussian point in the LIGGGHTS model. S5. After confirming that the actual strain increment in LIGGGHTS should be equal to the strain previously set in ABAQUS, calculate the stress state of the element in combination with the boundary force and the contact point location, and write it into a .txt file with a specific naming rule. The stress response of the Gaussian point representative element is fed back to the ABAQUS main program by the VUMAT subroutine according to the file name. S6. Following the central difference algorithm, the next cycle of calculation is performed based on the updated Gaussian point information according to the current calculation results. The speed of boundary loading is controlled to align the two in physical time. After the macroscopic calculation step is completed, cross-scale data interaction is completed through a .txt file. The calculation is repeated until the preset load conditions are reached.
2. The cross-scale numerical simulation method based on ABAQUS explicit finite element and LIGGGHTS discrete element method according to claim 1, characterized in that, In step S1, user-defined material is selected in ABAQUS to set the material properties of Gaussian points in the cross-scale calculation region.
3. The cross-scale numerical simulation method based on ABAQUS explicit finite element and LIGGGHTS discrete element method according to claim 1, characterized in that, In step S3, the LIGGGHTS particle set includes particle gradation, contact models between particles and between particles and boundaries, and contact parameters such as elastic modulus, Poisson's ratio, friction coefficient, and damping. Loads are applied and displacement loading or load loading is achieved through boundary constraint conditions.
4. The cross-scale numerical simulation method based on ABAQUS explicit finite element and LIGGGHTS discrete element method according to claim 1, characterized in that, Step S4 specifically involves: obtaining the Gaussian point number in ABAQUS using VUMAT, and naming the corresponding data file (i.e., .txt file) with the number; and writing the incremental time, total time, and strain increment matrix of the current incremental step into the corresponding .txt file named according to the Gaussian point using VUMAT. After writing is complete, VUMAT triggers the execution of an external Python script via a command. The Python script extracts the strain increment from the .txt file. Data; based on the time step required for the current computation step in ABAQUS. Combined with the number of time steps required to achieve this macroscopic step size as set in LIGGGHTS Calculate the time step size for each small step in LIGGGHTS ,Right now: ； By combining the strain increment information of the current calculation step, the speed of the current calculation step is calculated, so as to ensure that the LIGGGHTS program is consistent with the macroscopic increment state of ABAQUS after it is completed. Using the formula: ； Calculate the corresponding deformation gradient tensor F, where I is the unit tensor. For strain; the rotation axis is calculated by the cross product of the wall normal vectors before and after deformation, and the angle between the two normal vectors is calculated by the magnitude of the normal vector, and then the rotation speed is calculated; then the translation speed of the center point of the LIGGGHTS boundary surface and the rotation angular velocity of the boundary surface are determined in the corresponding time step. Subsequently, the Python script sends commands to the LIGGGHTS solver via an interface, applying the calculated kinematic parameters to the boundary. Once set up, LIGGGHTS runs for the specified number of steps, completing the LIGGGHTS calculations required for one incremental step in ABAQUS.
5. The cross-scale numerical simulation method based on ABAQUS explicit finite element and LIGGGHTS discrete element method according to claim 1, characterized in that, Step S5 specifically involves: confirming that the actual strain increment in LIGGGHTS should be equal to the strain previously set in ABAQUS; calculating the boundary stress through the force and area of the boundary wall; obtaining each stress component through the LIGGGHTS interface in Python; and writing it to the .txt file of the corresponding Gaussian point. In the VUMAT subroutine, the corresponding .txt file is located by looping through the element number, and the data is read into the stress-related variables. Then, the read values are assigned to the stressNew function to complete the update of the stress state of the Gaussian point in ABAQUS, which serves as the constitutive response input for macroscopic finite element analysis.
6. The cross-scale numerical simulation method based on ABAQUS explicit finite element and LIGGGHTS discrete element method according to claim 1, characterized in that, Step S6 specifically involves: in the subsequent time step iteration phase, by using explicit integration of the central difference, the motion information of the next transient step is directly extrapolated and calculated based on the current acceleration, velocity, and displacement state in each incremental step; Based on the VUMAT subroutine interface, the system enables real-time and efficient exchange of macro and micro information; the above process is executed cyclically until the preset external loading termination criteria are met.
7. A computer device comprising a memory, a processor, and a computer program stored in the memory, characterized in that, The processor executes the computer program to implement the steps of the method of claim 1.
8. A computer-readable storage medium having a computer program / instructions stored thereon, characterized in that, When the computer program / instructions are executed by the processor, they implement the steps of the method of claim 1.
9. A computer program product comprising a computer program / instructions, characterized in that, When the computer program / instructions are executed by the processor, they implement the steps of the method of claim 1.