Acceleration method of train-ballast track-roadbed system coupling simulation model

Abstract

The invention discloses a method and a system for accelerating a train-ballast track-roadbed system MBD-DEM-FDM coupling simulation model based on a self-adaptive PID (Proportion Integration Differentiation) algorithm. The method comprises the following steps: firstly, establishing a data interaction interface of MBD and DEM-FDM modules, and calculating and transmitting an acting force and a counter-acting force; a self-adaptive PID control algorithm is introduced, loading force is dynamically corrected according to interface errors, and fast approximation balance in a single time step is achieved; the convergence state of the system is monitored through a sliding window error judgment mechanism, a local fuzzy self-learning or global beetle antenna search optimization algorithm is intelligently triggered according to the convergence state, and PID control parameters are dynamically adjusted in real time and globally optimized. According to the method, a traditional time-consuming space iteration balance process is converted into dynamic feedback adjustment of a time domain, the problems that an existing coupling method is low in calculation efficiency, poor in numerical stability and difficult to apply to a large-scale model are effectively solved, and an efficient and reliable simulation tool is provided for design, operation and maintenance of a railway system.

Description

Technical Field

[0001] This invention relates to an acceleration method, and more particularly to an acceleration method for a coupled simulation model of a train-ballast track-subgrade system MBD-DEM-FDM based on an adaptive PID algorithm. Background Technology

[0002] Heavy-haul and some high-speed railway lines consist of ballasted track, roadbed, bridges, and tunnels, forming a multi-layered and heterogeneous longitudinal structure. To ensure the stable operation of the railway line, each structural layer needs to meet functional and technical requirements and work collaboratively. However, existing studies often separate the upper train, track structure, and substructure, neglecting their mutual influences and failing to fully characterize the overall mechanical properties of the train-track-roadbed coupled system.

[0003] Currently, some scholars have conducted extensive research on the dynamic response of train-ballast track-subgrade. Multibody dynamics (MBD) can describe the dynamic interaction between rails, sleepers, and trains well, but it usually simplifies the ballast layer as a spring-damped system, ignoring particle morphology and lateral shear effects, making it difficult to reflect the true microscopic mechanical behavior of ballast. Discrete element method (DEM) can simulate the contact and breakage behavior of single particles at the mesoscopic level. Finite difference method (FLAC) is suitable for representing the stress wave propagation and dynamic response of continuous media such as subgrade and foundation, and can be quickly coupled with DEM to characterize the interaction between ballast and subgrade, but it cannot independently realize the multibody interaction between train load and track superstructure. Therefore, some scholars have proposed a high-speed train-ballast track-subgrade system based on MBD-DEM-FDM coupling, pioneering the realization of three-dimensional fully coupled dynamic simulation of train-discrete medium ballast track-continuous medium subgrade system. However, these three types of models differ significantly in terms of time step, numerical stability, and data communication mechanism. The traditional iterative solution methods currently used often result in low computational efficiency, frequent numerical oscillations, large coupling errors, and poor parallel efficiency, making it difficult to extend to multi-sleeper and large-scale models.

[0004] Existing technologies, such as CN119862654A, disclose a simulation method for a train-ballast track-subgrade system based on MBD-DEM-FDM coupling, including: establishing a train model based on the multibody dynamics MBD method, including: establishing a dynamic model of the train system, simplifying the train's suspension system through springs and dampers, and considering the interaction between the train and the track; establishing a three-dimensional ballast track bed model based on the discrete element method DEM; establishing a subgrade finite difference model based on the finite difference method FDM, including: treating the subgrade as a continuous medium, simulating the mechanical behavior and deformation characteristics of the subgrade under train load based on the finite difference method FDM; setting a coupling medium to couple MBD-DEM-FDM, and performing real-time information exchange between the MBD, DEM, and FDM systems; establishing a coupled model based on the MBD-DEM-FDM coupling and verifying and applying the coupled model.

[0005] The data interaction interface between the MBD and DEM-FDM parts is as follows: Figure 1 As shown, the train motion response updated at this time step is transmitted to the corresponding position on the top of the sleeper in the DEM-FDM model through the fastener, and the reaction force of the track bed and subgrade system is calculated. If the difference between the action force and the reaction force meets the set threshold, the equilibrium is calculated for the dynamic analysis of the next time step. If this requirement is not met, the input force needs to be corrected. The correction formula is shown below:

[0006]

[0007] In the formula: F ai F is the initial iterative force at the current time step. r1 The reaction force is calculated at the current time step, and μ is the relaxation factor in the iteration process, which is used to ensure the computational efficiency of the correction process.

[0008] Although the above technologies have established an MBD-DEM-FDM coupling framework, their coupling mechanism still relies on traditional spatial domain iterative strategies, which has obvious limitations:

[0009] 1. Low computational efficiency: This method requires iterative updates of the forces and displacements at the coupling interface within each time step, and the convergence process relies on a fixed relaxation factor or a preset number of iterations. As the number of coupling interfaces increases, the number of iterations grows exponentially, leading to frequent data interactions, time-consuming convergence judgment, and a significant decrease in overall computational efficiency, making it impossible to realize multi-sleeper and long-distance track line models.

[0010] 2. Poor numerical stability: Due to the lack of a dynamic feedback adjustment mechanism, the force-displacement transmission between the DEM and FDM / MBD is prone to oscillations under conditions of high coupling stiffness or severe train excitation. Traditional methods rely on manually set iterative parameters to suppress oscillations, but the effect is unstable and iterative divergence easily occurs after error amplification, leading to unstable simulation operation.

[0011] 3. Model Scale Limitations: The convergence characteristics of existing methods are highly sensitive to computational scale. As the coupling region expands, data synchronization, iteration coordination, and memory consumption increase dramatically, making it difficult to support dynamic response studies of multi-sleeper or track-level structures. Therefore, CN119862654A is practically only applicable to coupling analysis of local structures or single-sleeper areas and cannot meet the engineering calculation needs of large-scale railway systems. Summary of the Invention

[0012] To overcome the shortcomings of existing MBD-DEM-FDM coupled simulation methods, such as CN119862654A, in terms of computational efficiency, numerical stability, and model scalability, this invention proposes an accelerated method for coupled simulation of train-ballast track-subgrade systems based on an adaptive PID control algorithm. This method improves the computational efficiency of traditional MBD-DEM-FDM coupled simulations, enabling dynamic response simulation of multi-sleeper, large-scale models of train-track-ballast bed systems. Compared with existing technologies, this invention has the following significant differences and advantages:

[0013] 1. Significantly Improved Computational Efficiency: Existing technologies employ traditional spatial domain iterative methods, requiring repeated updates of forces and displacements at the coupling interface within each time step. As the number of coupling elements increases, the number of iterations grows exponentially, severely limiting computational efficiency. This invention, however, transforms the coupling process into a dynamic feedback control process in the time domain. By introducing an adaptive PID error adjustment mechanism, it achieves single-step approximate equilibrium of interface forces, significantly reducing the number of iterations and fundamentally improving computational efficiency.

[0014] 2. Significantly Improved Numerical Stability: Existing technologies lack automatic parameter tuning capabilities and are prone to interface oscillations or even divergence when faced with high-stiffness coupling or strong excitation loads. This invention constructs a sliding window error judgment mechanism and, based on this, introduces a local fuzzy self-learning and global beetle antenna search optimization algorithm to achieve real-time dynamic adjustment and global optimization of PID parameters, effectively suppressing numerical oscillations and significantly improving coupling stability.

[0015] 3. Support for large-scale simulation: Existing technologies struggle to scale to large-scale models at the multi-sleeper or track level, and the complexity of data synchronization and iterative coordination increases dramatically with the computational domain. This invention reduces the computational burden of the spatial iteration process and lowers the data interaction density between subdomains through dynamic parameter optimization, enabling the MBD-DEM-FDM coupled simulation to maintain efficient and stable operation even under large-scale conditions. This numerical simulation method can reproduce actual on-site operating conditions, reducing the burden of on-site measurements and thus ensuring the safe operation and long-term service of railways.

[0016] This invention discloses an acceleration method based on a coupled simulation model of a train-ballast track-subgrade system, the technical solution of which is as follows:

[0017] An acceleration method for the MBD-DEM-FDM coupled simulation model of the train-ballast track-subgrade system based on an adaptive PID algorithm is characterized by the following time-series coupling steps:

[0018] Step 1: Calculate the equivalent force at the fastener in the three-dimensional direction in the upper MBD calculation module through the data interaction interface, obtain the force in the MBD calculation module of the previous time step, and apply the force to the sleeper at the corresponding coupling position in the DEM-FDM model;

[0019] Step 2: After calculating the same time step in the DEM-FDM system, extract the displacement and velocity of the sleeper in three dimensions, and calculate the reaction force. Among them, the reaction force in the Z direction is corrected based on the sleeper's rotational motion around the y-axis.

[0020] Step 3: Using the force from the previous time step as the initial iterative force for the current time step, introduce the proportional-integral-derivative PID control module and calculate the correction amount for the force based on the interface error of the previous time step.

[0021] Step 4: Correct the force applied at the current time step according to the correction amount, and input it into the DEM-FDM model to calculate the dynamic response of the updated rail and sleeper;

[0022] Step 5: Calculate the difference between the MBD action force and the DEM-FDM reaction force at the current time step to obtain the interface error at the current time step;

[0023] Step 6: Introduce a moving average stationarity determination mechanism, calculate the root mean square value of the interface error within a preset time window, and determine whether the convergence condition is met: the root mean square value of the interface error is less than the preset tolerance threshold within multiple consecutive time windows.

[0024] Step 7: Based on the judgment results of Step 6, perform parameter optimization:

[0025] If the convergence condition is not met, the global optimization mechanism based on the beetle antenna search BAS algorithm is triggered to perform global optimization correction on the proportional, integral and derivative gain coefficients of the PID control module.

[0026] If the convergence condition is met, the fuzzy self-learning module is triggered to perform local adaptive correction on the proportional, integral, and derivative gain coefficients of the PID control module.

[0027] Step 8: Record the updated dynamic response of the rails and sleepers at the current time step;

[0028] Step 9: Transfer all dynamic response data of the lower DEM-FDM structure to the upper MBD calculation module to calculate the fastener force at this step and use it as input for the next step.

[0029] This invention also discloses a train-track-subgrade coupling simulation acceleration system based on an adaptive PID algorithm for implementing the above method, characterized in that it includes:

[0030] The data interaction and force initialization module is used to execute step one;

[0031] The reaction force calculation and correction module is used to execute step two.

[0032] The PID control and force correction module is used to execute steps three and four.

[0033] The interface error calculation module is used to execute step five.

[0034] The convergence determination module is used to execute step six;

[0035] The intelligent parameter optimization module is used to perform step seven, and it includes a BAS global optimization unit and a fuzzy self-learning unit.

[0036] The data recording and transmission module is used to execute steps eight and nine.

[0037] Preferably, the data interaction and force initialization module is implemented using a Python program.

[0038] Preferably, the intelligent parameter optimization module is configured to receive the determination result from the convergence determination module and selectively activate the BAS global optimization unit or the fuzzy self-learning unit accordingly.

[0039] The present invention also discloses a non-volatile storage medium, characterized in that the non-volatile storage medium includes a stored program, wherein the program, when running, controls the device where the non-volatile storage medium is located to execute the above-described method.

[0040] This invention also discloses a terminal device, characterized in that the terminal device includes: a processor, a memory, a communication interface, and a bus; the processor, the memory, and the communication interface are connected via the bus and communicate with each other; the memory stores executable program code; the processor reads the executable program code stored in the memory to run a program corresponding to the executable program code, for executing the method described above. Beneficial Effects

[0041] 1. Significantly improves convergence efficiency

[0042] Existing technical limitations: Traditional iterative balancing methods rely on spatial iterative solutions, requiring multiple iterations for each step to approximate equilibrium, resulting in high computational costs, slow convergence speed, and susceptibility to initial conditions.

[0043] Technical advantages: This invention transforms the spatial iterative process into dynamic feedback regulation in the time domain, achieving single-step approximation equilibrium of interface forces through error-driven PID control. The system automatically corrects errors within each time step, reaching a stable state without multiple iterations, thereby significantly reducing computational load and improving convergence speed and numerical stability.

[0044] 2. Adaptive parameter tuning

[0045] Existing technical defects: Traditional PID controllers rely on manual parameter tuning, and the parameters are fixed and cannot be adjusted in real time according to system errors or changes in operating conditions, which often leads to overshoot, lag and instability under nonlinear or time-varying loads.

[0046] Technical advantages: This invention achieves K-axis parameter tuning by constructing a self-learning parameter tuning mechanism based on error amplitude, direction, and rate of change. P K I K D The online dynamic updating of parameters enables the control system to automatically identify the current error state and make immediate corrections. Compared with fixed-parameter PID, this method can adaptively cope with different loads, boundary conditions, and coupling conditions, achieving true "self-learning with error" control and significantly improving control accuracy and response speed.

[0047] 3. Global-local dual-layer optimization collaboratively enhances system intelligence.

[0048] Existing technology limitations: Traditional control systems only have single-layer control or single-time optimization capabilities, lacking multi-scale and multi-level adaptive adjustment logic, making it difficult to continuously optimize performance in complex coupled systems.

[0049] Technical advantages: This invention forms a two-layer collaborative control framework of "local self-learning + global BAS optimization": the local module is responsible for rapid response and real-time correction, while the global module performs cross-time step parameter optimization and long-term trend adjustment. This framework combines speed and globality, enabling the control system to maintain efficient and stable dynamic equilibrium even in nonlinear, multi-coupled, and multi-condition environments. Attached Figure Description

[0050] Figure 1 This is a schematic diagram of the data interaction interface between the MBD and DEM-FDM parts in the existing MBD-DEM-FDM coupling method;

[0051] Figure 2 This is a flowchart of the acceleration method for the MBD-DEM-FDM coupled simulation model of the train-ballast track-subgrade system based on the adaptive PID algorithm of the present invention. Detailed Implementation

[0052] To overcome the shortcomings of existing technologies in MBD-DEM-FDM train-ballast track-subgrade coupled systems, this invention proposes an accelerated simulation method for the MBD-DEM-FDM coupled simulation model of the train-ballast track-subgrade system based on an adaptive PID algorithm. This method optimizes the cross-software coupling mechanism using the adaptive PID algorithm, eliminating the need for step-by-step iteration, thus improving model running efficiency and establishing a large-scale coupled model that simultaneously reflects the dynamics of the train system and the microscopic behavior of ballast particles. The technical solution includes the following detailed steps:

[0053] Step 1: Establish a data interaction interface between the MBD and DEM-FDM parts using a Python program, and calculate the equivalent force F at the fastener in the three-dimensional direction in the MBD calculation module of the previous step. a (i-1) ={F ax1 , F ay1 , F az1 , F ax2 , F ay2 , F az2 The calculation formula is as follows, and it is applied to the sleepers at the corresponding coupling positions.

[0054] (1)

[0055] (2)

[0056] In the formula: F x1 F y1 F z1 These represent the reaction forces exerted by the left rail on the fastener in the x, y, and z directions, respectively; F x2 F y2 Fz2 These represent the reaction forces acting on the right-side rail in the x, y, and z directions, respectively, returning to the fastener; k x k y k z These represent the spring stiffness coefficients of the fastener in the x, y, and z directions, respectively; c x c y c z These represent the damping coefficients of the fastener in the x, y, and z directions, respectively; x rx1 x ry1 x rz1 These represent the displacements of the left rail in the x, y, and z directions, respectively; x rx2 x ry2 x rz2 These represent the displacements of the right-hand rail in the x, y, and z directions, respectively; v rx1 v ry1 v rz1 These represent the displacements of the left rail in the x, y, and z directions, respectively; v rx2 v ry2 v rz2 These represent the displacements of the right-hand rail in the x, y, and z directions, respectively; x sx x sy x sz These represent the displacements of the sleeper in the x, y, and z directions, respectively; v sx v sy v sz These represent the velocities of the sleeper in the x, y, and z directions, respectively.

[0057] Step 2: After calculating the same time step within the DEM system, extract the displacement and velocity of the sleeper in three dimensions at the last time step, and calculate the reaction force F using the following formula. r (i) ={F x1 , F y1 , F z1_new , F x2 , F y2 ,F z2_new It should be noted that, due to the significant difference in the forces on both sides of the sleeper in the three-dimensional state of this study caused by the introduction of DEM ballast particles, this study corrects the formulas as shown in (4) to (8) by considering the fixed-axis rotation of the rigid body and the three-dimensional forces. The force F in the Z direction is... z1_new and F z2_new The corrected formula is as follows:

[0058] (3)

[0059] (4)

[0060] (5)

[0061] (6)

[0062] (7)

[0063] In the formula: F z1_new and F z2_new For the corrected reaction force in the Z direction, θ y ω y α y To represent the rotation angle, angular velocity, and angular acceleration of the sleeper about the y-axis; l is half the distance between the left and right fasteners of the same sleeper, x sz1 x sz2 These represent the displacements of the left and right ends of the sleeper in the z-direction, respectively; v sz1 v sz2 These represent the velocities of the left and right ends of the sleeper in the z-direction.

[0064] Step 3: The initial iterative force is the fastening force F from the previous time step. a (i) =F a (i-1) A proportional-integral-derivative (PID) control module is introduced to realize real-time error calculation and dynamic correction of interface forces, so as to ensure the mechanical equilibrium accuracy and calculation convergence efficiency of the multi-physics coupling interface. The correction formula is as follows:

[0065] (8)

[0066] In the formula, ΔF a (i) Let e(i) be the fastener force correction amount, e(i) be the error at step i, and K be the fastener force correction amount. P (i), K I (i), K D (i) represents the proportional, integral, and derivative gain coefficients at time step i, initially set to 0.5, 0.1, and 0, respectively, with Δt being the time step. This control law achieves stable control and dynamic adjustment of the loading process by simultaneously considering the current value, cumulative value, and rate of change of the error.

[0067] Step 4: Adjust the applied force F at the current time step based on the calculated correction amount. a (i) The input to the DEM model is used to calculate and update the rail and sleeper responses. The calculation formula is as follows:

[0068] (9)

[0069] Step 5: The initial iterative force is the fastening force F from the previous time step. a (i) =F a (i-1) Calculate the force F applied by the MBD calculation module at the current time step. a (i) The reaction force F of the same DEM calculation module r (i) The difference is used to obtain the interface error e(i) at the current time step. The error calculation formula is as follows:

[0070] e(i)=F a (i) -F r (i) (10)

[0071] The interface error e(i) is defined as the norm of the error vector (such as the 2-norm) to ensure that it is treated as a scalar in subsequent steps.

[0072] Step 6: To prevent misjudgment and oscillation during sudden load changes or interface nonlinearity stages, this invention introduces a moving average stationarity determination mechanism in the interface error control stage. By performing dynamic smoothing judgment in the time domain, instantaneous noise interference is eliminated, achieving stable convergence detection of the error. Specifically, a time window ω is defined (5-10 time steps, preferably 8 time steps, with a tolerance threshold ε of 0.01, adjusted according to the specific model stiffness), and the root mean square error within the window is calculated:

[0073] (11)

[0074] If the following conditions are met within multiple consecutive time windows:

[0075] (12)

[0076] In the formula, RMSE ω If ω is the root mean square error within the time window, and ε is the preset tolerance threshold, then the system is considered to have reached convergence.

[0077] Step 7: To further enhance the system's adaptive capability under complex boundary conditions and multi-field coupling scenarios, this invention designs a self-learning module to perform local adaptive adjustment and global intelligent optimization of the PID control parameters, achieving dynamic self-evolution of the control law. During system operation, the self-learning module dynamically corrects the proportional, integral, and derivative parameters based on the amplitude, direction, and rate of change of the interface error.

[0078] (13)

[0079] In the formula, α P αI α D The adaptive step size coefficient for each parameter is α, whose value is empirically set based on the system stiffness and damping characteristics. P =10 -3 ~10 -2 α I =10 -4 ~10 -3 α D =10 -4 ~10 -2 By introducing a fuzzy logic adjustment mechanism to achieve nonlinear adaptive correction, when the error is large and changes rapidly, the proportional coefficient K is increased. P To improve response sensitivity; when the error is small but there is a steady-state deviation, the integral coefficient K should be appropriately increased. I To eliminate steady-state error; when error oscillations are frequent or the system has an overshoot tendency, increase the differential coefficient K. D This is to improve damping and suppress oscillations.

[0080] To prevent local adaptive adjustment from getting stuck in local optima, a global optimization mechanism based on the Beetle Antennae Search (BAS) algorithm is integrated into the self-learning module. This mechanism is triggered after running several coupled time steps or when the root mean square error fails to reach the convergence criterion within a set time.

[0081] The optimization process first defines the overall fitness function, with the root mean square error of the time window as the optimization objective. The fitness function is as follows, with the current parameter vector K = [K...]. P , K I , K D [Starting from], execute a global optimization strategy based on BAS.

[0082] (14)

[0083] In the formula, J=[K P , K I , K D ] represents the selected fitness function.

[0084] Establish a unit random vector b to represent the random position and orientation of the longhorn beetle's two antennae, since its movement direction is random during foraging.

[0085] (15)

[0086] In the formula, rmd() is a function that generates uniformly distributed random numbers in the range [-1, 1], and || and ||2 are the L2 norms.

[0087] Determine the coordinates of the left and right tentacles of the longhorn beetle at time t.

[0088] (16)

[0089] (17)

[0090] In the formula, p r p l Let d be the coordinates of the left and right antennae. t Let p be the distance from the beetle's center of mass to its antenna at time t. t This represents the final position of the longhorn beetle after flight time t.

[0091] At time t, by comparing the fitness function J values ​​of the left and right antennae, if the J value of the left antenna is less than that of the right antenna, it indicates that the food odor on the right is stronger, and the longhorn beetle will fly to the right; otherwise, it will fly to the left. The longhorn beetle's position after flight time t is:

[0092] (18)

[0093] In the formula, s t The search step size factor at time t, and sign() is the sign function.

[0094] The location of the longhorn beetle at position p is calculated using formula (15). t+l The fitness function value J(p) at time 1 t+l ), and with optimal adaptation J best In comparison, if J(p) t+l ) <J best Then we get J best =J(p t+l ), p best =p t+l。 Step 8: The specific coupling process is as follows: the initial iterative force is the fastening force under the previous time step, and the input force is corrected based on adaptive PID in this time step to update the dynamic response of the rail and sleeper. The reaction force is updated according to Equations 4 and 5. At this time, it is calculated and judged whether the root mean square error judgment condition of the time window is met. If it is not met, the BAS corrects the proportional, integral and derivative gain coefficients. If it is met, the fuzzy self-learning module corrects the proportional, integral and derivative gain coefficients. The rail and sleeper response at this time is recorded.

[0095] Step 9: Transfer all dynamic response data of the substructure to the upper train-track calculation module, calculate the fastener force of this step and input it to the next step.

[0096] This invention, based on a three-dimensional numerical calculation model of MBD-DEM-FDM coupling, can overcome the shortcomings of existing studies in analyzing the overall mechanical behavior of coupled systems, especially considering the micromechanical properties of particles within the track bed and the comprehensive impact of train dynamics on the track and subgrade. However, existing MBD-DEM-FDM coupling methods typically rely on time-step iterative solutions, resulting in extremely high computational costs and limiting the application of large-scale coupled models. To address these issues, this invention proposes an accelerated method and system for the MBD-DEM-FDM coupled simulation model of the train-ballast track-subgrade system based on an adaptive PID algorithm. This method combines the proportional-integral-derivative (PID) control principle with an intelligent parameter optimization mechanism to real-time correct the coupling force, displacement, and error signals at the train-sleeper-ballast interface, and dynamically corrects the PID parameters (K) using an adaptive adjustment strategy. p K I K D The adaptive PID algorithm achieves rapid convergence of input and reaction forces, avoiding the bottleneck of extensive repetitive calculations in traditional iterative algorithms. Utilizing comprehensive feedback from current error, historical error, and changing trends, the PID algorithm enables dynamic equilibrium force prediction and correction without iteration, significantly improving computational efficiency and reducing memory consumption. Furthermore, the adaptive PID algorithm avoids the poor adaptability of fixed-parameter PID controllers to model uncertainties, which necessitates frequent manual parameter tuning and limits the dynamic simulation and prediction capabilities of multi-field dynamic responses between train, track bed, and roadbed.

[0097] The foregoing has shown and described the basic principles, main features, and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The embodiments and descriptions in the specification are merely principles of the invention. Various changes and modifications can be made to the invention without departing from its spirit and scope, and all such changes and modifications fall within the scope of the claimed invention. The scope of protection claimed by the appended claims and their equivalents is defined.

Claims

1. An acceleration method for a coupled simulation model of a train-ballast track-subgrade system based on an adaptive PID algorithm (MBD-DEM-FDM), characterized in that... This includes the following timing coupling steps: Step 1: Calculate the equivalent force at the fastener in the three-dimensional direction in the upper MBD calculation module through the data interaction interface, obtain the force in the MBD calculation module of the previous time step, and apply the force to the sleeper at the corresponding coupling position in the DEM-FDM model; Step 2: After calculating the same time step in the DEM-FDM system, extract the displacement and velocity of the sleeper in three dimensions, and calculate the reaction force. Among them, the reaction force in the Z direction is corrected based on the sleeper's rotational motion around the y-axis. Step 3: Using the force from the previous time step as the initial iterative force for the current time step, introduce the proportional-integral-derivative PID control module and calculate the correction amount for the force based on the interface error of the previous time step. Step 4: Correct the force applied at the current time step according to the correction amount, and input it into the DEM-FDM model to calculate the dynamic response of the updated rail and sleeper; Step 5: Calculate the difference between the MBD action force and the DEM-FDM reaction force at the current time step to obtain the interface error at the current time step; Step 6: Introduce a moving average stationarity determination mechanism, calculate the root mean square value of the interface error within a preset time window, and determine whether the convergence condition is met: the root mean square value of the interface error is less than the preset tolerance threshold within multiple consecutive time windows. Step 7: Based on the judgment results of Step 6, perform parameter optimization: If the convergence condition is not met, the global optimization mechanism based on the beetle antenna search BAS algorithm is triggered to perform global optimization correction on the proportional, integral and derivative gain coefficients of the PID control module. If the convergence condition is met, the fuzzy self-learning module is triggered to perform local adaptive correction on the proportional, integral, and derivative gain coefficients of the PID control module. Step 8: Record the updated dynamic response of the rails and sleepers at the current time step; Step 9: Transfer all dynamic response data of the lower DEM-FDM structure to the upper MBD calculation module to calculate the fastener force at this step and use it as input for the next step.

2. The method according to claim 1, characterized in that, In step one, the calculation of the equivalent force is based on the spring stiffness coefficient and damping coefficient of the fastener in three dimensions, as well as the relative displacement and relative velocity between the rail and the sleeper: （1） （2） In the formula: F x1 F y1 F z1 These represent the reaction forces exerted by the left rail on the fastener in the x, y, and z directions, respectively; F x2 F y2 F z2 These represent the reaction forces acting on the right-side rail in the x, y, and z directions, respectively, returning to the fastener; k x k y k z These represent the spring stiffness coefficients of the fastener in the x, y, and z directions, respectively; c x c y c z These represent the damping coefficients of the fastener in the x, y, and z directions, respectively; x rx1 x ry1 x rz1 These represent the displacements of the left rail in the x, y, and z directions, respectively; x rx2 x ry2 x rz2 These represent the displacements of the right-hand rail in the x, y, and z directions, respectively; v rx1 v ry1 v rz1 These represent the displacements of the left rail in the x, y, and z directions, respectively; v rx2 v ry2 v rz2 These represent the displacements of the right-hand rail in the x, y, and z directions, respectively; x sx x sy x sz These represent the displacements of the sleeper in the x, y, and z directions, respectively; v sx v sy v sz These represent the velocities of the sleeper in the x, y, and z directions, respectively.

3. The method according to claim 2, characterized in that, In step two, the formula for correcting the reaction force in the Z direction is: （3） （4） （5） （6） （7） In the formula: F z1_new and F z2_new For the corrected reaction force in the Z direction, θ y ω y α y To represent the rotation angle, angular velocity, and angular acceleration of the sleeper about the y-axis; l is half the distance between the left and right fasteners of the same sleeper, x sz1 x sz2 These represent the displacements of the left and right ends of the sleeper in the z-direction, respectively; v sz1 v sz2 These represent the velocities of the left and right ends of the sleeper in the z-direction, respectively.

4. The method according to claim 1, characterized in that, In step three, the formula for calculating the correction amount of the PID control module is as follows: （8） In the formula, ΔF a (i) Let e(i) be the fastener force correction amount, e(i) be the error value at step i, and K be the fastener force correction amount. P (i), K I (i), K D (i) are the proportional, integral, and differential gain coefficients at time step i, respectively, and Δt is the time step.

5. The method according to claim 1, characterized in that, In step four, the initial iterative force is the fastening force F from the previous time step. a (i) =F a (i-1) Based on the calculated correction amount, the applied force F at the current time step is corrected. a (i) The input to the DEM model is used to calculate and update the rail and sleeper responses. The calculation formula is as follows: （9）。 6. The method according to claim 1, characterized in that, In step five, the force F exerted by the MBD calculation module at the current time step is calculated. a (i) The reaction force F of the same DEM calculation module r (i) The difference is used to obtain the interface error e(i) at the current time step. The error calculation formula is as follows: e(i)=F a (i) -F r (i) (10); In the formula, the interface error e(i) is defined as the norm of the error vector (such as the 2-norm) to ensure that it is treated as a scalar in subsequent steps.

7. The method according to claim 1, characterized in that, In step six, a time window ω is defined, and the root mean square error within the window is calculated: （11） If the following conditions are met within multiple consecutive time windows: （12） In the formula, RMSE ω If ω is the root mean square error within the time window, and ε is the preset tolerance threshold, then the system is considered to have reached convergence.

8. The method according to claim 1, characterized in that, In step seven, the self-learning module dynamically corrects the proportional, integral, and differential parameters based on the amplitude, direction, and rate of change of the interface error. （13） In the formula, α P α I α D These are the adaptive step size coefficients for each parameter.

9. A train-track-subgrade coupling simulation acceleration system based on an adaptive PID algorithm for implementing the method of any one of claims 1 to 8, characterized in that, include: The data interaction and force initialization module is used to execute step one; The reaction force calculation and correction module is used to perform steps two and three. The PID control and force correction module is used to execute step four. The interface error calculation module is used to execute step five. The convergence determination module is used to execute step six; The intelligent parameter optimization module, used to perform step seven, includes a BAS global optimization unit and a fuzzy self-learning unit. The data recording and transmission module is used to execute steps eight and nine.

10. The system according to claim 7, characterized in that, The data interaction and force initialization module is implemented using a Python program.

11. The system according to claim 8, characterized in that, The intelligent parameter optimization module is configured to receive the determination result from the convergence determination module and selectively activate the BAS global optimization unit or the fuzzy self-learning unit accordingly.

12. A non-volatile storage medium, characterized in that, The non-volatile storage medium includes a stored program, wherein the program, when executed, controls the device where the non-volatile storage medium is located to perform the method described in any one of claims 1 to 6.

13. A terminal device, characterized in that, The terminal device includes: a processor, a memory, a communication interface, and a bus; the processor, the memory, and the communication interface are connected through the bus and communicate with each other; the memory stores executable program code; the processor reads the executable program code stored in the memory to run a program corresponding to the executable program code, so as to perform the method as described in any one of claims 1-6 above.

Images