CFD-DEM coupling acceleration method

Abstract

The invention provides a CFD-DEM coupling acceleration method, which is applied to the field of computational fluid dynamics and discrete element method coupling simulation and comprises a CFD module, a DEM module, a coupling controller and an optimization module. The DEM module predicts particle distribution based on LSTM to generate a probability distribution diagram; a gradient is calculated in combination with flow field data in a CFD module, and an irregular dynamic coupling region is divided by using a geometric deep learning model; constructing a GNN-based proxy model, and adaptively switching the proxy model and a traditional drag force model based on Reynolds number; the coupling controller monitors a calculation state, and triggers data synchronization of a hierarchical compression strategy when an asynchronous time step threshold value is reached; the optimization module monitors a momentum error after data synchronization, and if the momentum error exceeds a threshold value, a neural network parameter and a threshold value coefficient are fed back and adjusted; according to the method, the problems that a traditional method is low in calculation efficiency, rigid in region division and lack of error control are solved, the calculation efficiency is remarkably improved by 5-10 times, the data synchronization amount is reduced by 58-65%, and it is guaranteed that the precision error is smaller than or equal to 5%.

Description

Technical Field

[0001] The present invention relates to the technical field of coupled simulations of computational fluid dynamics (CFD) and discrete element method (DEM), and more particularly to a CFD-DEM coupling acceleration method based on intelligent partitioning of dynamic coupling regions, adaptive switching of proxy models, and asynchronous time-stepping collaborative optimization. This method significantly improves computational efficiency by 5-10 times and ensures an accuracy error of ≤5% by dynamically optimizing the coupling mechanism, thus breaking through the computational bottleneck of traditional synchronous coupling. This method is suitable for high-precision and efficient simulation of particle-fluid interactions, with specific applications in chemical reactor optimization, fluidized bed design, and particle transport system analysis. Background Art

[0002] Particle-fluid multiphase flow systems are widely present in industrial scenarios such as energy, chemical industry, and pharmaceutical industry, and their accurate simulation is crucial for equipment design and process optimization. The CFD-DEM coupling method can restore complex multiphase flow behavior with high fidelity by solving the dynamic equations of the fluid phase (CFD) and the discrete particle phase (DEM) separately and exchanging the interphase forces (such as drag and lift) in real time. However, the traditional coupling method has the following significant bottlenecks: (1) Low computational efficiency: DEM needs to track the microscopic motion of millions of particles (time step of about 10 -6 s), while the CFD time step is usually larger than 10 2 ~10 4 times. Forced synchronous coupling causes the DEM module to be in a waiting state for more than 90% of the time, and the computing resource utilization rate is less than 40%, especially when the local particle volume fraction is higher than 30%; (2) Rigid regional division: Fixed coupling regions or uniform partitioning strategies cannot adapt to the dynamic aggregation / diffusion behavior of particles, resulting in redundant calculations in low-interaction regions; (3) Insufficient model adaptability: Traditional drag models (such as Gidaspow and DiFelice) have reduced accuracy in transient high Reynolds number flows, while pure data-driven proxy models are prone to non-physical solutions due to the lack of conservation law embedding; (4) Lack of error control: Asynchronous calculations are prone to accumulate phase errors, and existing methods lack an adaptive feedback mechanism for the coupling state, resulting in divergent simulation results.

[0003] Although there have been some attempts to improve existing technologies, there are still limitations. Asynchronous coupling strategies (such as MPC algorithms) allow CFD / DEM to adopt different time steps, but the threshold setting relies on empirical parameters. When the flow field suddenly changes (such as when the particle surge velocity changes by more than 20%), it is easy to cause coupling instability, and the cumulative error of the interphase force exceeds 15%. Secondly, the application of surrogate models is mostly limited to steady-state scenarios, and the model switching problem of transient high Reynolds number flow is not solved. In addition, the optimization of the coupling region is usually based on regular geometric divisions, which makes it difficult to identify unstructured particle aggregation patterns (such as chain clusters and sparse dispersions).

[0004] At the same time, with the simultaneous improvement of the demand for multiphase flow simulation accuracy and efficiency in industrial scenarios (such as lithium battery electrode slurry mixing and fluidized bed combustion optimization), higher requirements are placed on advanced technologies that can support efficient and accurate simulation of engineering-level large-scale particle systems.

[0005] In summary, there is an urgent need to develop a coupling acceleration framework that integrates intelligent partitioning of dynamic coupling regions, adaptive switching of proxy models, asynchronous time step collaborative optimization, and physical constraint embedding to achieve high-precision and high-efficiency CFD-DEM coupling simulation, which is necessary to break through the existing computing bottleneck and meet the needs of complex industrial applications. Summary of the Invention

[0006] The embodiment of the present invention provides a CFD-DEM coupling acceleration method, which aims to solve the problems in traditional CFD-DEM coupling technology, such as the large coupling area and large amount of data synchronization caused by the use of regular boundary division, the slow calculation speed and low efficiency of the traditional drag model, and the lack of error control mechanism.

[0007] In order to solve the above-mentioned technical problems, the present invention is achieved as follows:

[0008] An embodiment of the present invention provides a CFD-DEM coupling acceleration method, which includes at least a CFD module, a DEM module, a coupling controller, and an optimization module; a dynamic coupling region is divided between the CFD module and the DEM module; within the dynamic coupling region, the coupling controller monitors the calculation state and determines whether an asynchronous time step threshold is reached, and performs data synchronization after the asynchronous time step threshold is reached; the optimization module is used to monitor the coupling error after the data synchronization, determine whether the coupling error reaches an error threshold, and feedback-adjust parameters of a neural network model and threshold parameters if the error threshold is reached; if the asynchronous time step threshold or the error threshold is not reached, the calculation of the original CFD module and DEM module continues to be executed.

[0009] The steps for dividing the dynamic coupling area are as follows: first, input the initial flow field data of the CFD module and the particle properties (position, velocity) of the DEM module; then, in the DEM module, predict the particle distribution based on the long short-term memory neural network (LSTM), generate a probability distribution map of the particle distribution and output it to the CFD module; then, in the CFD module, calculate the velocity gradient, pressure gradient and turbulence intensity in combination with the flow field data, and determine the area where the particle concentration exceeds the dynamic coupling area threshold based on the probability distribution map; finally, use the geometric deep learning model to identify the particle aggregation pattern, divide and output the irregularly shaped dynamic coupling area.

[0010] The neural network model includes at least one of the LSTM, the geometric deep learning model, and a dynamic graph neural network (GNN).

[0011] The thresholds include but are not limited to: the dynamic coupling region threshold, the asynchronous time step threshold, and the error threshold.

[0012] The data synchronization adopts a hierarchical compression strategy; the hierarchical compression strategy clusters particles through a spatial hash algorithm, uses ZFP lossless compression for key data, uses FP8 quantization for flow velocity in low gradient areas, and uses run-length encoding for static particle positions.

[0013] The flow field data at least includes fluid velocity vector, pressure, and turbulence intensity.

[0014] The training data of the LSTM, the geometric deep learning model, and the GNN all come from high-fidelity CFD-DEM coupled simulation, and data enhancement technology is used to improve generalization.

[0015] The dynamic coupling region is based on the particle aggregation point set output by the geometric deep learning, and the α-shape algorithm (α = 0.5d p ~2d p , where d p A tightly packed smooth convex hull boundary (where φ is the particle diameter) is constructed, which is used to define the closed computational domain for CFD-DEM coupled data exchange.

[0016] The geometric deep learning model includes at least a graph attention network, PointNet++, and a graph isomorphism network.

[0017] The CFD-DEM coupling acceleration method is used to calculate two-phase flow, including gas-solid two-phase flow and liquid-solid two-phase flow; the gas phase and the liquid phase are compressible or incompressible fluids, and the solid phase is particles of regular or irregular shape.

[0018] The CFD-DEM coupling acceleration method is applicable to gas-solid fluidized beds, lithium battery slurry stirring tanks, and coal powder pneumatic conveying systems, and is particularly applicable to strongly non-steady-state conditions where the particle concentration dynamically changes by more than 40%, such as the preparation process of carbon nanotubes in a fluidized bed.

[0019] The CFD module includes at least a proxy model.

[0020] The construction of the proxy model at least includes: in the CFD module, using particles as node features and particle-flow field grid interaction information as edge features, inputting the dynamic coupling region, the particle properties of the DEM module, and the flow field data; constructing the proxy model based on GNN to replace the drag model based on physical equations, and calculating the particle-flow field interaction force F coupling: the local Reynolds number of the surrogate model is in the low Reynolds number interval, the interphase force F is calculated by the drag force model instead coupling GNN : when the local Reynolds number is in the high Reynolds number interval, the interphase force F is automatically switched back to be calculated by the drag force model based on physical equations from the surrogate model coupling Drag : in the transition Reynolds number interval, the interphase force F is calculated by the weighted fusion of the surrogate model and the drag force model in parallel coupling GNN Drag , wherein the weight ; wherein is a Sigmoid activation function, and PMRM is a reliability index of the surrogate model; and the calculation result of the surrogate model or the drag force model is output to the DEM module; the DEM module continuously calculates the motion of the particles according to the calculation result of the surrogate model or the drag force model until the coupling controller monitors that the asynchronous time step threshold is reached and then performs the data synchronization.

[0021] The node features at least include particle position particle velocity particle size d p , particle temperature T p , forming a 6-dimensional feature vector; the edge features at least include local flow velocity difference grid pressure gradient , relative distance , constituting a 4-dimensional edge feature; the particle temperature T p is an optional input, which is enabled when the simulation involves heat exchange.

[0022] The LSTM at least includes a spatio-temporal attention mechanism and a loss function; the loss function at least includes Newton's motion equation; through the spatio-temporal attention mechanism, it can focus on the high gradient area and avoid resource waste caused by uniform processing of all data; the Newton's motion equation is embedded into the loss function as a physical constraint to improve the credibility of single particle motion.

[0023] The high gradient area refers to the area with high velocity gradient, pressure gradient and turbulent intensity in the flow field, and the interaction between fluid and particles in these areas is strong, at least including turbulent boundary layer and particle aggregation area.

[0024] The formula of the spatio-temporal attention mechanism at least includes attention weight vector, adjusted input vector and LSTM state update; the formula of the attention weight vector is Wherein, a t ​​​​is the attention weight vector at time step t; σ is the Sigmoid activation function, and the output range is (0-1); W a is a learnable parameter matrix; h t-1 is the historical hidden state of the LSTM unit at time step t-1; is the position coordinate vector of the particle at time step t-1; is the velocity vector of the particle at time step t-1; is the velocity vector of the fluid at the position of the particle at time step t-1; the formula of the adjusted input vector is wherein, is the attention weight a t weighted adjusted particle position input vector; ⊙ represents the Hadamard product; the formula of the LSTM state update is wherein, h t is the updated LSTM hidden state.

[0025] The Hadamard product represents element-wise multiplication, which adaptively weights the particle position coordinates and suppresses redundant data calculation in the low gradient area.

[0026] The low gradient area refers to an area with low velocity gradient, pressure gradient and turbulence intensity in the flow field, and the interaction between the fluid and the particle in the area is weak.

[0027] The spatio-temporal attention mechanism dynamically focuses on the high gradient area by adjusting the attention weight a t , the historical hidden state h t-1 , the particle position particle velocity and local flow field velocity .

[0028] The formula of the Newtonian motion equation includes at least the calculation method of the LSTM predicted displacement increment, the Newtonian law displacement increment, the fused position prediction and the dynamic weight; the formula of the LSTM predicted displacement increment is wherein, is the predicted displacement increment of the particle at time step t based on the historical state prediction by the LSTM unit; the formula of the Newtonian law displacement increment is wherein, is the predicted displacement increment of the particle at time step t based on the Newton's second law at time step t-1; Δt dem is the time step of the DEM module; is the acceleration calculated by the total force acting on the particle at time step t-1 The formula of the fused position prediction is wherein, is the particle position at time step t; η is a dynamic weight (0-1) to balance the contribution of LSTM predicted displacement and Newton's law predicted displacement; the formula for calculating the dynamic weight is wherein, is the modulus of velocity gradient tensor; is the maximum value of velocity gradient modulus in the current flow field; P density is the normalized particle distribution probability (value range 0-1).

[0029] The LSTM learns the nonlinear motion pattern through the historical state to predict the path; the Newton's second law calculates the displacement increment of the particle based on force analysis; the dynamic weight η increases the weight of the LSTM in the high gradient area, representing that the contribution proportion of the LSTM predicted displacement is large, indicating that the nonlinear flow field effect is more significant, and vice versa, in the low gradient area, the Newton's mechanics is relied on for calculation.

[0030] The specific way of the geometric deep learning model to divide the dynamic coupling area is: input the probability distribution map, velocity gradient, pressure gradient and turbulent intensity; analyze the local features and global features of particle aggregation based on graph convolution network (GCN); generate the dynamic coupling area according to dense cluster, chain distribution and sparse dispersion mode.

[0031] The GCN extracts spatial features by aggregating particle neighborhood information.

[0032] The neighborhood is defined as the particles in the spherical region with the target particle as the center and the radius of 2-3 times the particle size.

[0033] Wherein, the radius of 2-3 times the particle size is set based on the physical scale of particle-fluid interaction: this range can cover the particle wake effect and near-field drag interference area, while avoiding too large radius to introduce irrelevant noise.

[0034] The local feature is the aggregation of the hidden state and spatial distance information of directly adjacent particles; the global feature is the screening of a set of potentially remote acting particles, and then capturing the non-local topological structure.

[0035] The surrogate model at least includes multi-scale spatio-temporal aggregation and loss function; the loss function at least includes mass-momentum conservation equation; the multi-scale spatio-temporal aggregation is realized by GCN operation, which can realize cross-scale coupling modeling and comprehensively cover the particle-flow field coupling effect; the mass-momentum conservation equation is embedded in the loss function as a physical constraint to improve the rationality and stability of gas-solid flow field simulation.

[0036] The formula of the multi-scale spatio-temporal aggregation includes message passing function and node state update; the formula of the message passing function is wherein, is the message transmitted from node j to node i at time step t; φ is the MLP function for message generation and state update; and is the hidden state of nodes i and j at time step t; is the Euclidean distance between nodes i and j; the formula for updating the node status is in, is the updated hidden state of node i at time step t+1; ψ is the MLP function for message generation and state update; is the set of directly adjacent nodes of node i; is the set of nodes (not directly adjacent) that have potential remote interactions with node i, selected by the spatial hashing algorithm; is a message delivered from remote node k to node i.

[0037] The φ includes two fully connected layers (the number of neurons is 64→32), and the activation function is ReLU.

[0038] The formula of the mass-momentum conservation equation includes a total loss function, an adaptive weight, and a turbulence intensity; the formula of the total loss function is in, is the total loss function; is the mean square error loss, which measures the deviation between the predicted value of the proxy model and the true value; λ is the weight coefficient of the physical constraint term; N cell is the number of CFD grid cells; is the mass conservation equation (continuity equation). For incompressible flow, ρ is a constant and the equation is simplified to is the fluid inertia term (material derivative); is the pressure gradient term; μ is the dynamic viscosity, is the viscous force term; F coupling is the interphase force coupling force, calculated by the agent model or the drag model; ρg is the gravity term; the formula of the adaptive weight is λ=0.1+0.9·I t , where I t is the turbulence intensity index; the formula for the turbulence intensity is Where k is the turbulent kinetic energy; ||u|| is the local average flow velocity; ε is a minimum constant that is always zero.

[0039] Among them, ε is a very small constant (typical value is 10 -6 ) is used to prevent numerical errors caused by zero denominators. When the local flow velocity ||u|| approaches zero (such as at stagnation points or in low Reynolds number regions), ε ensures numerical stability and prevents abnormal program termination.

[0040] The formula for the dynamic coupling region threshold is: Among them, TH region is the dynamic coupling region threshold; is the modulus of the velocity gradient tensor; is the modulus of the pressure gradient vector; α, β, γ are weight coefficients and satisfy α+β+γ=1.

[0041] The formula of the reliability index PMRM of the proxy model is: in, L is the characteristic length, which is set as the particle size or grid size according to the simulation object; is the modulus of the time derivative of the pressure gradient; ||u|| is the modulus of the velocity vector; f1(Re), f2(Re), f3(Re) are the adaptive coefficients of dynamic adjustment,

[0042] The reliability index PMRM of the proxy model is dynamically calculated; the low Reynolds number interval of the proxy model is Re<500, the high Reynolds number interval of the proxy model is Re>2000, and the transition Reynolds number interval of the proxy model is 500≤Re≤2000.

[0043] The asynchronous time step threshold is TH async =k1·CFL+k2·N coll ·τ coll , where TH async is the asynchronous time step threshold; CFL is the Courant-Friedrichs-Lewy number, CFL<1, u max is the maximum flow velocity, Δt is the current time step, Δx is the characteristic size of the CFD grid; N coll is the total number of particle collisions that occurred within the most recent DEM coupling step; τ coll is the characteristic time of a single particle collision; k1 and k2 are empirical coefficients determined by pre-calibration experiments (k1∈[0.1,0.3], k2∈[0.05,0.2]).

[0044] The parameters of the thresholds, including but not limited to α, β, γ of the dynamic coupling region thresholds and k1 and k2 of the asynchronous time step thresholds, are determined by parameter scanning or optimization algorithms based on sensitivity analysis of typical working conditions with the goal of minimizing the coupling error.

[0045] The layered compression strategy includes: performing lossless compression on key data and lossy compression on secondary data; clustering particles into clusters using the spatial hash algorithm, with the cluster side length being 2 to 3 times the particle diameter.

[0046] The cube cluster size of 2 to 3 times the particle diameter is a dedicated parameter for data compression, which matches the average spacing of a typical particle system, ensures a high spatial correlation between the positions of particles within the cluster, and avoids the loss of compression accuracy due to excessively large clusters.

[0047] The lossless compression includes at least the ZFP algorithm, which retains full precision for key data (such as particle force, flow field pressure, and coupling interface velocity); the lossy compression includes at least FP16 half-precision, SZ compression, FP8 precision quantization, and run-length coding compression. The algorithm is dynamically selected according to the data type, and the FP16 half-precision or SZ compression or FP8 precision quantization is used for secondary data (such as the flow field velocity and flow field temperature in the low gradient area). The run-length coding compression is used for the particle position outside the dynamic coupling area and the stationary particle velocity.

[0048] The method for determining the error threshold is to calculate the momentum error function; the momentum error function is Among them, E m is the momentum error function; N is the total number of particles; F DEM,i F is the interphase force actually used by the DEM module in the previous calculation step and calculated by the CFD module; CFD,i The force acting on the position of particle i is obtained by interpolating the phase forces on the CFD grid nodes around particle i using the inverse distance weighted interpolation method.

[0049] The inverse distance weighted interpolation is a local interpolation method based on the inverse of distance weighting. The attribute value of a point in space is obtained by weighted averaging the attribute values ​​of the surrounding known points. The closer the point is, the greater the influence on the interpolation result.

[0050] In an embodiment of the invention, by using LSTM to predict particle distribution and combining it with a geometric deep learning model to identify particle aggregation patterns, it is possible to adaptively identify unstructured particle dynamic behavior, significantly reducing redundant calculation areas by 30% to 50%, and effectively overcoming the problems of excessive coupling areas and high data synchronization caused by the use of regular boundary divisions; a proxy model constructed using GNN replaces the traditional drag model for calculations under low Reynolds number flow states, and automatically switches back to the drag model based on physical equations at high Reynolds numbers, so that the simulation accuracy error under transient high Reynolds number conditions is controlled to ≤5%, realizing adaptive intelligent switching of the proxy model and solving the bottleneck of low computational efficiency of the traditional drag model; the data synchronization amount is significantly reduced through a hierarchical compression strategy, and the coupling error feedback mechanism is used to dynamically adjust key threshold parameters (such as asynchronous time step threshold and dynamic coupling area threshold), achieving a 40% to 60% reduction in synchronization waiting time, while stably controlling the coupling momentum error within a preset threshold (≤3%), effectively improving overall computing efficiency and enhancing error control capabilities. The embodiments of the present invention solve the problems in traditional CFD-DEM coupling technology, such as excessively large coupling area and large amount of data synchronization caused by the use of regular boundary division, slow and inefficient calculation of traditional drag models, and lack of error control mechanism. BRIEF DESCRIPTION OF THE DRAWINGS

[0051] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the following briefly introduces the drawings required for use in the description of the embodiments of the present invention. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative labor.

[0052] Figure 1 A schematic diagram showing the structure of a CFD-DEM coupling acceleration method provided by an embodiment of the present invention;

[0053] Figure 2 A schematic diagram showing a basic flow chart of a CFD-DEM coupling acceleration method provided by an embodiment of the present invention;

[0054] Figure 3 A schematic diagram showing a flow chart of a proxy model in a CFD-DEM coupling acceleration method provided by an embodiment of the present invention;

[0055] Description of reference numerals:

[0056] 10. CFD module; 20. DEM module; 30. Coupling controller; 40. Optimization module. DETAILED DESCRIPTION

[0057] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of them. All other embodiments obtained by ordinary technicians in this field based on the embodiments of the present invention without making any creative efforts shall fall within the scope of protection of the present invention.

[0058] It should be understood that references throughout this specification to "one embodiment" or "an embodiment" mean that a particular feature, structure, or characteristic associated with the embodiment is included in at least one embodiment of the present invention. Therefore, the appearances of "in one embodiment" or "in an embodiment" throughout this specification do not necessarily refer to the same embodiment. Furthermore, these particular features, structures, or characteristics may be combined in any suitable manner in one or more embodiments.

[0059] See also Figures 1 to 3 The embodiment of the present invention provides a CFD-DEM coupling acceleration method for the field of computational fluid dynamics and discrete element method coupling simulation, comprising at least a CFD module (10), a DEM module (20), a coupling controller (30), and an optimization module (40); a dynamic coupling region is divided in the CFD module (10) and the DEM module (20); within the dynamic coupling region, the coupling controller (30) monitors the calculation state and determines whether an asynchronous time step threshold is reached, and performs data synchronization after reaching the asynchronous time step threshold; the optimization module (40) is used to monitor the coupling error after the data synchronization, determine whether the coupling error reaches an error threshold, and if the error threshold is reached, feedback is used to adjust the parameters of the neural network model and the coefficient of the threshold; if the asynchronous time step threshold or the error threshold is not reached, the calculation of the original CFD module (10) and the DEM module (20) is continued.

[0060] The steps for dividing the dynamic coupling region are as follows: first, inputting the initial flow field data of the CFD module (10) and the particle properties (position, velocity) of the DEM module (20); then, in the DEM module (20), predicting the particle distribution based on a long short-term memory neural network (LSTM), generating a probability distribution map of the particle distribution and outputting it to the CFD module (10); then, in the CFD module (10), calculating the velocity gradient, pressure gradient and turbulence intensity in combination with the flow field data, and judging the particle region exceeding the threshold of the dynamic coupling region based on the probability distribution map; finally, using a geometric deep learning model to identify the particle aggregation pattern, dividing and outputting the irregularly shaped dynamic coupling region.

[0061] The neural network model includes at least one of the LSTM, the geometric deep learning model, and a dynamic graph neural network (GNN).

[0062] The thresholds include but are not limited to: the dynamic coupling region threshold, the asynchronous time step threshold, and the error threshold.

[0063] The data synchronization adopts a layered compression strategy.

[0064] The flow field data at least includes fluid velocity vector, pressure, and turbulence intensity.

[0065] The CFD module (10) is used to solve the fluid control equation and calculate the interphase force; the DEM module (20) tracks the particle motion trajectory based on Newton's law of motion; the coupling controller (30) is used to dynamically monitor the calculation state, execute the asynchronous time step threshold judgment and the data synchronization; the optimization module (40) is used to analyze the coupling error in real time and feedback-adjust the neural network parameters and the threshold coefficient.

[0066] In an embodiment of the present invention, by adopting LSTM to predict particle distribution and combining with geometric deep learning model to identify particle aggregation pattern, it is possible to adaptively identify unstructured particle dynamic behavior, significantly reduce redundant calculation area by 30% to 50%, and effectively overcome the problems of excessive coupling area and high data synchronization caused by the use of regular boundary division; the proxy model constructed by GNN replaces the traditional drag model for calculation under low Reynolds number flow state, and automatically switches back to the drag model based on physical equations at high Reynolds number, so that the simulation accuracy error under transient high Reynolds number condition is controlled at ≤5%, realizing adaptive intelligent switching of proxy model and solving the problem of traditional drag model. The invention overcomes the bottleneck of low computational efficiency of the CFD-DEM coupling technology; significantly reduces the amount of data synchronization through a layered compression strategy, and dynamically adjusts key threshold parameters (such as the asynchronous time step threshold and the dynamic coupling area threshold) using a coupling error feedback mechanism, thereby shortening the synchronization waiting time by 40% to 60%. At the same time, the coupling momentum error is stably controlled within a preset threshold (≤3%), effectively improving the overall computational efficiency and strengthening the error control capability; the embodiments of the invention solve the problems of the conventional CFD-DEM coupling technology in the prior art, such as the excessively large coupling area and the large amount of data synchronization caused by the use of regular boundary division, the slow and inefficient calculation of the conventional drag model, and the lack of an error control mechanism.

[0067] Preferably, the CFD module (10) comprises at least a proxy model.

[0068] Preferably, the construction of the proxy model at least includes: in the CFD module (10), using particles as node features and particle-flow field grid interaction information as edge features, inputting the dynamic coupling region, the particle properties (position, velocity) of the DEM module (20), and the flow field data; constructing the proxy model based on GNN to replace the drag model based on physical equations, and calculating the particle-flow field interaction force F coupling The local Reynolds number of the proxy model replaces the drag model to calculate the interphase force F when the local Reynolds number is in the low Reynolds number range coupling =F GNN :When the local Reynolds number is in the high Reynolds number range, the proxy model is automatically switched back to the drag model based on the physical equation to calculate the interphase force F coupling =F Drag In the transition Reynolds number range, the proxy model and the drag model are used to calculate the weighted fusion interphase force in parallel: F coupling =ω·F GNN +(1-ω)·F Drag , where the weight ;in, is a Sigmoid activation function, PMRM is a reliability index of the proxy model; and the calculation results of the proxy model or the drag model are output to the DEM module (20); the DEM module (20) continuously calculates the movement of the particles according to the calculation results of the proxy model or the drag model until the coupling controller (30) monitors that the asynchronous time step threshold is reached and then performs the data synchronization.

[0069] It should be noted that the proxy model is deeply integrated with the solver of the CFD module (10) through the FluentUDF interface; the FluentUDF interface includes but is not limited to: model mounting points, data channel optimization, and real-time hot switching; the proxy model is mounted on the DEFINE_EXCHANGE_PROPERTY macro to realize the real-time hot switching of the interphase force calculation.

[0070] It should be noted that a multi-physics data bridging engine needs to be implemented in the coupling controller (30) to solve the problem of heterogeneous data compatibility between CFD and DEM software, including at least grid-particle mapping optimization and data type conversion; the grid-particle mapping optimization uses a spatial hash grid to accelerate neighbor search, and the complexity is reduced from O(N 2) is reduced to O(N); the data type conversion converts the Euler field data (including but not limited to velocity / pressure tensors) of the CFD module (10) into Lagrangian particle force vectors that can be recognized by the DEM module (20); the data bridging engine uses FLUENT_RP_GET_REAL to read the flow field gradient in real time; the coupling controller (30) is an independent C++ module.

[0071] Preferably, the LSTM includes at least a spatiotemporal attention mechanism and a loss function; the loss function includes at least Newton's equation of motion; through the spatiotemporal attention mechanism, it is possible to focus on areas with high velocity gradients, high pressure gradients, and high turbulence intensity, avoiding the waste of resources caused by uniformly processing all data; the Newton's equation of motion is embedded in the spatiotemporal attention mechanism as a physical constraint to enhance the credibility of single-particle motion.

[0072] It should be noted that the DEM module (20) injects the particle motion trajectory into the LSTM predictor through the API.

[0073] Preferably, the formula of the spatiotemporal attention mechanism includes at least an attention weight vector, an adjusted input vector, and an LSTM state update; the formula of the attention weight vector is Among them, a t is the attention weight vector at time step t; σ is the Sigmoid activation function with an output range of (0 to 1); W a is the learnable parameter matrix; h t-1 is the historical hidden state of the LSTM unit at time step t-1; is the position coordinate vector of the particle at time step t-1; is the velocity vector of the particle at time step t-1; is the velocity vector of the fluid at the particle position at time step t-1; the formula of the adjusted input vector is in, is the attention weight a t The weighted particle position input vector; ⊙ represents the Hadamard product; the formula for updating the LSTM state is: Among them, h t is the updated LSTM hidden state.

[0074] It should be noted that directly inputting high-dimensional flow field data into the LSTM leads to computational redundancy. The flow field gradient filter is introduced only when The spatiotemporal attention mechanism is calculated while

[0075] It should be noted that the LSTM requires the state of particles and flow fields, including but not limited to position speed Local flow field velocity Obtained through FluentUDF interpolation.

[0076] Preferably, the formula of the Newton's motion equation includes at least LSTM predicted displacement increment, Newton's law displacement increment, fusion position prediction, and dynamic weight calculation method; the formula of the LSTM predicted displacement increment is in, The particle displacement increment at time step t is predicted by the LSTM unit based on the historical state; the formula for the Newton law displacement increment is in, is the displacement increment at time step t predicted based on Newton's second law at time step t-1; Δt dem is the time step of the DEM module (20); is the acceleration calculated from the sum of the forces acting on the particle at time step t-1 The formula for fusion position prediction is: in, is the particle position at time step t; η is a dynamic weight (0-1) used to balance the contribution of LSTM predicted displacement and Newton's law predicted displacement; the formula for calculating the dynamic weight is η = in, is the modulus of the velocity gradient tensor; is the maximum value of the velocity gradient modulus in the current flow field; P density is the normalized particle distribution probability (range 0 to 1).

[0077] It should be noted that the historical step length of the LSTM is 5 to 15 steps; preferably, the historical step length of the LSTM is 10 steps.

[0078] Preferably, the specific way of dividing the dynamic coupling area by the geometric deep learning model is: inputting the probability distribution map, velocity gradient, pressure gradient and turbulence intensity; analyzing the local and global characteristics of particle aggregation based on the graph convolutional network (GCN); generating the dynamic coupling area according to dense clusters, chain distribution and sparse dispersion patterns.

[0079] It should be noted that the boundary of the dynamic coupling region adopts a non-convex envelope algorithm, including but not limited to NURBS curve fitting.

[0080] It should be noted that the geometric deep learning model accurately identifies irregular structures such as the fluidized bed bubbling area and bubble tail vortex chain.

[0081] Preferably, the GCN extracts spatial features by aggregating particle neighborhood information; the neighborhood is defined as particles within a spherical region centered on the target particle and with a radius of 2 to 3 times the particle diameter.

[0082] Preferably, a particle size of 3 times the radius is selected.

[0083] Preferably, the proxy model includes at least multi-scale spatiotemporal aggregation and a loss function; the loss function includes at least the mass-momentum conservation equation; the multi-scale spatiotemporal aggregation is realized through GCN operation, which can enable cross-scale coupling modeling and comprehensively cover the particle-flow field coupling effect; the mass-momentum conservation equation is embedded in the multi-scale spatiotemporal aggregation as a physical constraint to improve the rationality and stability of gas-solid flow field simulation.

[0084] Preferably, the formula of the multi-scale spatiotemporal aggregation includes a message passing function and a node state update; the formula of the message passing function is: in, is the message transmitted from node j to node i at time step t; φ is the MLP function for message generation and state update; and is the hidden state of nodes i and j at time step t; is the Euclidean distance between nodes i and j; the formula for updating the node status is in, is the updated hidden state of node i at time step t+1; ψ is the MLP function for message generation and state update; is the set of directly adjacent nodes of node i; is the set of nodes (not directly adjacent) that have potential remote interactions with node i, selected by the spatial hashing algorithm; is a message delivered from remote node k to node i.

[0085] Preferably, the formula of the mass-momentum conservation equation includes a total loss function, an adaptive weight, and a turbulence intensity; the formula of the total loss function is in, is the total loss function; is the mean square error loss, which measures the deviation between the predicted value of the proxy model and the true value; λ is the weight coefficient of the physical constraint term; N cell is the number of CFD grid cells; is the mass conservation equation (continuity equation). For incompressible flow, ρ is a constant and the equation is simplified to is the fluid inertia term (material derivative); is the pressure gradient term; μ is the dynamic viscosity, is the viscous force term; Fcoupling is the interphase force coupling force, calculated by the agent model or the drag model; ρg is the gravity term; the formula of the adaptive weight is λ=0.1+0.9·I t , where I t is the turbulence intensity index; the formula for the turbulence intensity is Where k is the turbulent kinetic energy; ||u|| is the local average flow velocity; ε is a minimum constant that is always zero.

[0086] Preferably, the formula for the dynamic coupling region threshold is: Among them, TH region is the dynamic coupling region threshold; is the modulus of the velocity gradient tensor; is the modulus of the pressure gradient vector; α, β, γ are weight coefficients and satisfy α+β+γ=1.

[0087] It should be noted that the grid unit TH region When the value is >0.4, it is marked as a candidate for the dynamic coupling region.

[0088] Preferably, the formula of the reliability index PMRM of the proxy model is: in, L is the characteristic length, which is set as the particle size or grid size according to the simulation object; is the modulus of the time derivative of the pressure gradient; ||u|| is the modulus of the velocity vector; f1(Re), f2(Re), f3(Re) are the adaptive coefficients of dynamic adjustment,

[0089] Preferably, the asynchronous time step threshold is TH async =k1·CFL+k2·N coll ·τ coll , where TH async is the asynchronous time step threshold; CFL is the Courant-Friedrichs-Lewy number, CFL<1, u max is the maximum flow velocity, Δt is the current time step, Δx is the characteristic size of the CFD grid; N coll is the total number of particle collisions that occurred within the most recent DEM coupling step; τ coll is the characteristic time of a single particle collision; k1 and k2 are empirical coefficients determined by pre-calibration experiments (k1∈[0.1,0.3], k2∈[0.05,0.2]).

[0090] It should be noted that TH async>0.35, the data synchronization is performed.

[0091] Preferably, the layered compression strategy includes: performing lossless compression on key data and lossy compression on secondary data; clustering particles into clusters using the spatial hash algorithm, with the cluster side length being 2 to 3 times the particle size.

[0092] The method for determining the error threshold is to calculate the momentum error function; the momentum error function is Among them, E m is the momentum error function; N is the total number of particles; F DEM,i F is the interphase force acting on particle i calculated by the DEM module (20); CFD,i The force acting on the position of particle i is obtained by interpolating the phase forces on the CFD grid nodes around particle i using the inverse distance weighted interpolation method.

[0093] It should be noted that F DEM,i F is the force actually used by the DEM module (20) in the previous calculation step (calculated and transferred by the CFD module (10) at an earlier time). CFD,i The interphase force is recalculated (using the same rules of the proxy model or the drag model) based on the current latest flow field state solved by the CFD module (10) and the current latest position of the particle i at the time of data synchronization; the momentum error function E m It reflects the cumulative degree of deviation in calculation of interphase forces caused by asynchronous calculation of the CFD module (10) and the DEM module (20) and the delay in data synchronization.

[0094] It should be noted that when E m >0.03, starts the optimization process of the parameters and threshold coefficients of the feedback adjustment neural network model.

[0095] It should be noted that the CFD module (10) of the present invention is compatible with a variety of fluid dynamics solvers, including but not limited to ANSYS Fluent and OpenFOAM; the DEM module (20) is compatible with a variety of discrete element solvers, including but not limited to EDEM and LIGGGHTS; and the coupling controller (30) realizes cross-platform data interaction through a standardized interface (such as UDF, API).

[0096] It should be noted that, in this document, the terms "comprises," "includes," or any other variations thereof are intended to encompass non-exclusive inclusion, such that a process, method, article, or apparatus comprising a series of elements includes not only those elements but also other elements not explicitly listed, or elements inherent to such process, method, article, or apparatus. In the absence of further limitations, an element defined by the phrase "comprising a ..." does not exclude the presence of other identical elements in the process, method, article, or apparatus comprising the element.

[0097] Through the description of the above embodiments, those skilled in the art can clearly understand that the above-mentioned embodiment methods can be implemented by means of software plus the necessary general hardware platform, or of course by hardware, but in many cases the former is a better embodiment. Based on this understanding, the technical solution of the present invention, or the part that contributes to the prior art, can be embodied in the form of a software product. The computer software product is stored in a storage medium (such as ROM / RAM, magnetic disk, optical disk), and includes a number of instructions for enabling a terminal (which can be a mobile phone, computer, server, air conditioner, or network device, etc.) to execute the methods described in each embodiment of the present invention.

[0098] The embodiments of the present invention are described above in conjunction with the accompanying drawings, but the present invention is not limited to the above-mentioned specific implementation methods. The above-mentioned specific implementation methods are merely illustrative and not restrictive. Under the guidance of the present invention, ordinary technicians in this field can also make many forms without departing from the scope of protection of the present invention and the claims, all of which are protected by the present invention. Example 1

[0099] In order to more clearly illustrate the technical implementation of the present invention, the following Figures 1 to 3 The specific implementation steps of the present invention are described in detail through a typical application scenario, such as the simulation of the fluidization process of carbon nanotube particles in a gas-solid fluidized bed reactor.

[0100] First, data preparation and system initialization were performed. A conical fluidized bed (height 1.2m, bottom diameter 0.02m, top diameter 0.05m) was constructed in Fluent 2019, using an unstructured hexahedral mesh (boundary layer encryption) with a total mesh size of 2.15 million. The gas was set to incompressible nitrogen with a density of 1.225kg / m 3 , the dynamic viscosity is 1.7894×10 -5kg / (m·s), the velocity inlet is 6m / s, the gauge pressure at the pressure outlet is 0Pa, and the wall is a no-slip boundary. EDEM2022 is used to simulate the particle phase. The particle system contains 500,000 spherical particles to simulate carbon nanotube agglomerates. The mixed particles have a dual particle size distribution, with small particles having a diameter of 150μm, accounting for 60%, and large particles having a diameter of 300μm, accounting for 40%. The particle density is 150kg / m 3 , Poisson's ratio is 0.25, and the shear modulus is 1×10 10 Pa, the contact model between particles and between particles and the wall adopts Hertz-Mindlin with JKR. Based on the analysis of the time scale of particle collision, the time step of DEM is set to 1×10 -6 s, based on the CFL condition and flow field stability requirements, the CFD time step is set to 1×10 -5 s.

[0101] Secondly, the dynamic coupling area is intelligently divided. A two-layer LSTM network (each layer has a hidden state dimension of 128) is deployed in FluentUDF, and the particle positions of the first 10 DEM steps are input. and speed And the interpolated flow field velocity at the particle position at the corresponding moment Applying spatiotemporal attention mechanism The network output is a particle existence probability distribution map P on a 800×800 two-dimensional grid covering the entire CFD calculation domain. density (Each grid value range is [0, 1]). Calculate the global velocity gradient tensor at each time step in the Fluent UDF and pressure gradient Then calculate the dynamic coupling region threshold After parameter scanning optimization, we first fixed β = 0.25, γ = 0.15, and scanned in the range of α∈[0.5,0.7] with a step size of 0.01, and determined that the error was the smallest when α = 0.6; then fixed α = 0.6, γ = 0.15, and scanned in the range of β∈[0.2,0.3], and determined that β = 0.25 was the optimal; finally fixed α = 0.6, β = 0.25, and scanned in the range of γ∈[0.1,0.2], and determined that γ = 0.15 was the optimal. In this case, α = 0.6, β = 0.25, and γ = 0.15 were finally used to calculate TH for each grid cell. region . Mark the previous step more than TH regionIn this study, a geometric deep learning model, such as the PointNet++ model (containing four set abstraction layers and two feature propagation layers), was used to process the point cloud. The model outputs point-level clustering pattern labels (0: sparse dispersion, 1: chain distribution, 2: dense cluster) and feature vectors. Based on the predicted labels, all particle point sets labeled as dense clusters (label 2) and chain distributions (label 1) were extracted. The α-shape algorithm was applied to generate an irregular, smooth convex hull boundary, which defined the dynamic coupling region. The sparsely dispersed region (label 0) was excluded from the coupling calculation. During simulations, this strategy accurately identified the dense bubbling region at the bottom of the fluidized bed, the chain-like particle belt in the rising bubble wake, and the dilute phase region at the top. The dynamic coupling region occupied an average of only 28-35% of the original computational domain, reducing data synchronization by over 40%.

[0102] Then, a proxy model is constructed in the dynamic coupling region and the model is switched adaptively. In the dynamic coupling region, each particle is regarded as a graph node, and its node feature is the particle position. Particle speed Particle size d p , particle temperature T p (6-dimensional), whose edge feature is the local flow velocity difference Grid pressure gradient Relative distance (4D). Construct a 4-layer dynamic graph neural network, each layer contains: message passing function (φ: 2-layer MLP, 64→32 neurons, ReLU), node update function (ψ: 2-layer MLP, 64→64 neurons, ReLU), the output layer is MLP (64->3 neurons), predicting the interphase force F acting on particle i coupling Calculate the local Reynolds number of particle i. When Re<500, use the proxy model to calculate F GNN As F coupling When Re>2000, switch to Gidaspow drag model to calculate F Drag As F coupling , when 500≤Re≤2000, the proxy model and drag model are calculated in parallel, and the final force is the weighted fusion F coupling =ω·F GNN +(1-ω)·F Gidaspow , where the weight in, is the Sigmoid activation function, PMRM is the reliability index of the proxy model; the formula of the reliability index of the proxy model is in, L is the characteristic length, which is set as the particle size or grid size according to the simulation object; is the modulus of the time derivative of the pressure gradient; ||u|| is the modulus of the velocity vector; f1(Re), f2(Re), f3(Re) are the adaptive coefficients of dynamic adjustment, A nonlinear weight function with Sigmoid transition is used at the boundary of Re=500 / 2000. The calculated F coupling Passed to EDEM 2022 for particle motion solution.

[0103] Then, based on the calculation results of CFD and DEM, it is determined whether asynchronous data synchronization is required. At the end of each CFD time step, the coupling controller calculates N coll Count the total number of collisions among all particles in the most recent complete coupling step, τ coll The characteristic time of a single collision of a small particle of 150 μm is calculated based on Hertz theory. The typical value is about 1.2×10 -6 , and then calculate the asynchronous time step threshold TH async =k1·CFL+k2·N coll ·τ coll In this case, k1=0.2,k2=0.12 are used after calibration. async If the value is greater than 0.35 (the preset trigger threshold), data synchronization is triggered. Otherwise, EDEM2022 continues to run on its smaller time step and Fluent2019 enters the next time step. When executing data synchronization, the layered compression strategy is used. The total force on the particle, the flow field pressure, and the fluid velocity near the coupling interface (in the dynamic coupling region) are measured by the ZFP compression library. In ABS mode with a tolerance of 1×10 -8 Under these conditions, the compressibility ratio reaches 1.8:1; in the low turbulence intensity region within the dynamic coupling region (I t <0.05) fluid velocities are stored in FP16 half-precision, achieving a 50% compression rate. Particle positions outside the dynamic coupling region are compressed using run-length encoding (processed separately on each coordinate axis), achieving a compression rate of 65-75% during the stable simulation phase. By combining these strategies, the average amount of data transmitted per synchronization is reduced by approximately 58-65% compared to the full-precision, uncompressed raw data.

[0104] Finally, error feedback control is performed, and the neural network parameters and threshold coefficients are adjusted according to the results of data synchronization. After each data synchronization is completed, the optimization module calculates the momentum error (normalized relative error), F DEM,i The interphase force calculated and applied to particle i by the DEM module (i.e., F transferred from the CFD module in the previous step) coupling ), F CFD,iIn the CFD module, the interphase force is recalculated at the current position of particle i based on the current latest flow field using the surrogate model or drag model. The calculation uses inverse distance weighted interpolation to search for the interphase force field on the 8 nodes of the CFD grid cell where particle i is located and its adjacent nodes. If the normalized momentum error E m >0.03, the optimization module starts the adjustment process. First, based on the current error E m and the flow / particle field state (such as the average Reynolds number Re, the average particle concentration), a small step gradient descent or rule-based empirical adjustment is used, such as using the gradient-based Adam optimizer (learning rate 10 -4 ) to fine-tune the asynchronous threshold coefficients k1, k2 and the region threshold weights α, β, γ. For example: if the error is mainly due to the lack of synchronization in the high collision area, slightly increase k2; if the error is due to insufficient division in the high gradient area, slightly increase β. Second, under the premise of preserving the main network weight of LSTM, only the particle motion trajectory data of the recent few time steps is used, and the Newton motion equation is used as a physical constraint term to fine-tune the parameter matrix W a of the spatiotemporal attention layer in LSTM, with a learning rate set to one-tenth of the initial training. The training target is the position prediction error of 2 CFD steps, while satisfying the weak constraint of Newton's second law. The updated parameters (k1, k2, α, β, γ) and the fine-tuned LSTM attention layer weight W a are applied to the subsequent simulation process.

[0105] The above embodiment 1 is quantitatively analyzed and compared with the traditional CFD-DEM coupling method on a single node computing platform AMD EPYC 9754 (128 cores), 512GB DDR5 RAM, NVIDIA Tesla A100 computing under the same conditions. The calculation time, CPU utilization and data synchronization volume of the present embodiment 1 are significantly better than those of the traditional CFD-DEM coupling method. Specifically, the calculation time of the present embodiment 1 is 8 times that of the traditional CFD-DEM coupling method, which indicates that the present method significantly improves the calculation efficiency and fully embodies the effectiveness of the present application in solving the problems of low efficiency, rigid region division, insufficient model adaptability and lack of error control of the traditional coupling method. The performance comparison between the traditional CFD-DEM coupling method and the method of the present application is shown in the following table:

Claims

1. A CFD-DEM coupling acceleration method, applied to the field of computational fluid dynamics and discrete element method coupling simulation, characterized by: The invention is used to realize the accelerated calculation of a CFD-DEM coupling model, and at least comprises a CFD module (10), a DEM module (20), a coupling controller (30), and an optimization module (40); a dynamic coupling region is divided in the CFD module (10) and the DEM module (20); in the dynamic coupling region, the coupling controller (30) monitors the calculation state and determines whether an asynchronous time step threshold is reached, and performs data synchronization after the asynchronous time step threshold is reached; the optimization module (40) is used to monitor the coupling error after the data synchronization, determine whether the coupling error reaches an error threshold, and feedback-adjust the parameters of the neural network model and the parameters of the threshold if the error threshold is reached; if the asynchronous time step threshold or the error threshold is not reached, the calculation of the original CFD module (10) and the DEM module (20) is continued; The steps of dividing the dynamic coupling region are as follows: first, inputting the initial flow field data of the CFD module (10) and the particle attributes (position, velocity) of the DEM module (20); then, in the DEM module (20), predicting the particle distribution based on a long short-term memory neural network (LSTM), generating a probability distribution map of the particle distribution and outputting it to the CFD module (10); then, in the CFD module (10), calculating the velocity gradient, pressure gradient and turbulence intensity in combination with the flow field data, and judging the particle region exceeding the threshold of the dynamic coupling region based on the probability distribution map; finally, using a geometric deep learning model to identify the particle aggregation pattern, dividing and outputting the irregularly shaped dynamic coupling region; The neural network model includes at least one of the LSTM, the geometric deep learning model, and a dynamic graph neural network (GNN); The thresholds include but are not limited to: the dynamic coupling region threshold, the asynchronous time step threshold, and the error threshold; The data synchronization adopts a layered compression strategy; The flow field data at least includes fluid velocity vector, pressure, and turbulence intensity.

2. A CFD-DEM coupling acceleration method according to claim 1, characterized in that: The CFD module (10) includes at least a proxy model; The construction of the proxy model at least includes: in the CFD module (10), using particles as node features and particle-flow field grid interaction information as edge features, inputting the dynamic coupling area, the particle properties (position, velocity) of the DEM module (20), and the flow field data; constructing the proxy model based on GNN to replace the drag model based on physical equations, and calculating the particle-flow field phase interaction force F coupling :The local Reynolds number of the proxy model is in the low Reynolds number range (Re≤Re low ) is used to replace the drag model to calculate the interphase force F coupling =F GNN ; When the local Reynolds number is in the high Reynolds number range (Re≥Re high ), the proxy model is automatically switched back to the drag model based on the physical equation to calculate the interphase force F coupling =F Drag ; In the transition Reynolds number range (Re low ≤Re≤Re high ), the proxy model and the drag model are used to calculate the weighted fusion in parallel to calculate the interphase force: F coupling =ω·F GNN +(1-ω)·F Drag , where the weight in, is a Sigmoid activation function, PMRM is a reliability index of the proxy model; the calculation result of the proxy model or the drag model is output to the DEM module (20); the DEM module (20) continuously calculates the movement of the particles according to the calculation result of the proxy model or the drag model until the coupling controller (30) monitors that the asynchronous time step threshold is reached and then performs the data synchronization.

3. A CFD-DEM coupling acceleration method according to claim 1, characterized in that: The LSTM includes at least a spatiotemporal attention mechanism and a loss function; the loss function includes at least Newton's equations of motion; the spatiotemporal attention mechanism can focus on areas with high velocity gradients, high pressure gradients, and high turbulence intensity, avoiding the waste of resources caused by uniformly processing all data; the Newton's equations of motion are embedded in the spatiotemporal attention mechanism as physical constraints to improve the credibility of single-particle motion; The formula of the spatiotemporal attention mechanism includes at least the attention weight vector, the adjusted input vector, and the LSTM state update; the formula of the attention weight vector is Among them, a t is the attention weight vector at time step t; σ is the Sigmoid activation function with an output range of (0 to 1); W a is the learnable parameter matrix; h t-1 is the historical hidden state of the LSTM unit at time step t-1; is the position coordinate vector of the particle at time step t-1; is the velocity vector of the particle at time step t-1; is the velocity vector of the fluid at the particle position at time step t-1; the formula of the adjusted input vector is in, is the attention weight a t The weighted particle position input vector; ⊙ represents the Hadamard product; the formula for updating the LSTM state is: Among them, h t is the updated LSTM hidden state; The formula of the Newton motion equation includes at least LSTM predicted displacement increment, Newton's law displacement increment, fusion position prediction, and dynamic weight calculation method; the formula of the LSTM predicted displacement increment is in, The particle displacement increment at time step t is predicted by the LSTM unit based on the historical state; the formula for the Newton law displacement increment is in, is the displacement increment at time step t predicted based on Newton's second law at time step t-1; Δt dem is the time step of the DEM module (20); is the acceleration calculated from the sum of the forces acting on the particle at time step t-1 The formula for fusion position prediction is: in, is the particle position at time step t; η is the dynamic weight (0-1), which is used to balance the contribution of LSTM predicted displacement and Newton's law predicted displacement; the formula for calculating the dynamic weight is in, is the modulus of the velocity gradient tensor; is the maximum value of the velocity gradient modulus in the current flow field; P density is the normalized particle distribution probability (range 0 to 1).

4. The CFD-DEM coupling acceleration method according to claim 1, characterized in that: The specific method of dividing the dynamic coupling region by the geometric deep learning model is as follows: inputting the probability distribution map, velocity gradient, pressure gradient and turbulence intensity; analyzing the local and global characteristics of particle aggregation based on a graph convolutional network (GCN); generating the dynamic coupling region according to dense cluster, chain distribution and sparse dispersion patterns; The GCN extracts spatial features by aggregating particle neighborhood information; the neighborhood is defined as particles within a spherical region centered on the target particle and with a radius of 2 to 3 times the particle diameter.

5. A CFD-DEM coupling acceleration method according to claim 2, characterized in that: The proxy model includes at least multi-scale spatiotemporal aggregation and a loss function; the loss function includes at least the mass-momentum conservation equation; the multi-scale spatiotemporal aggregation is achieved through GCN operation, which enables cross-scale coupled modeling and comprehensively covers the particle-flow field coupling effect; the mass-momentum conservation equation is embedded in the multi-scale spatiotemporal aggregation as a physical constraint to improve the rationality and stability of gas-solid flow field simulation; The formula of the multi-scale spatiotemporal aggregation includes a message passing function and a node state update; the formula of the message passing function is: in, is the message transmitted from node j to node i at time step t; φ is the MLP function for message generation and state update; and is the hidden state of nodes i and j at time step t; is the Euclidean distance between nodes i and j; the formula for updating the node status is in, is the updated hidden state of node i at time step t+1; ψ is the MLP function for message generation and state update; is the set of directly adjacent nodes of node i; is the set of nodes (not directly adjacent) that have potential remote interactions with node i, selected by the spatial hashing algorithm; is the message transmitted from remote node k to node i; The formula of the mass-momentum conservation equation includes a total loss function, an adaptive weight, and a turbulence intensity; the formula of the total loss function is in, is the total loss function; is the mean square error loss, which measures the deviation between the predicted value of the proxy model and the true value; λ is the weight coefficient of the physical constraint term; N cell is the number of CFD grid cells; is the mass conservation equation (continuity equation). For incompressible flow, ρ is a constant and the equation is simplified to is the fluid inertia term (material derivative); is the pressure gradient term; μ is the dynamic viscosity, is the viscous force term; F coupling is the interphase force coupling force, calculated by the agent model or the drag model; ρg is the gravity term; the formula of the adaptive weight is λ=0.1+0.9·I t , where I t is the turbulence intensity index; the formula for the turbulence intensity is Where k is the turbulent kinetic energy; ‖u‖ is the local average flow velocity; ε is a minimum constant that is always zero.

6. A CFD-DEM coupling acceleration method according to claim 1, characterized in that: The formula for the dynamic coupling region threshold is: Among them, TH region is the dynamic coupling region threshold; is the modulus of the velocity gradient tensor; is the modulus of the pressure gradient vector; α, β, γ are weight coefficients and satisfy α+β+γ=1.

7. A CFD-DEM coupling acceleration method according to claim 2, characterized in that: The formula of the reliability index PMRM of the proxy model is: in, L is the characteristic length, which is set as the particle size or grid size according to the simulation object; is the modulus of the time derivative of the pressure gradient; ‖u‖ is the modulus of the velocity vector; f1(Re), f2(Re), f3(Re) are the adaptive coefficients of dynamic adjustment, 8. The CFD-DEM coupling acceleration method according to claim 1, characterized in that: The asynchronous time step threshold is TH async =k1·CFL+k2·N coll ·τ coll , where TH async is the asynchronous time step threshold; CFL is the Courant-Friedrichs-Lewy number, CFL<1, u max is the maximum flow velocity, Δt is the current time step, Δx is the characteristic size of the CFD grid; N coll is the total number of particle collisions that occurred within the most recent DEM coupling step; τ coll is the characteristic time of a single particle collision; k1 and k2 are empirical coefficients determined by pre-calibration experiments (k1∈[0.1,0.3], k2∈[0.05,0.2]).

9. A CFD-DEM coupling acceleration method according to claim 1, characterized in that: The hierarchical compression strategy includes: performing lossless compression on key data and lossy compression on secondary data; clustering particles into clusters using a spatial hashing algorithm, with the cluster side length being 2 to 3 times the particle size.

10. The CFD-DEM coupling acceleration method according to claim 1, characterized in that: The method for determining the error threshold is to calculate the momentum error function; the momentum error function is Among them, E m is the momentum error function; N is the total number of particles; F DEM,i F is the interphase force acting on particle i calculated by the DEM module (20); CFD,i The force acting on the position of particle i is obtained by interpolating the phase forces on the CFD grid nodes around particle i using the inverse distance weighted interpolation method.

Images