A CFD-DEM coupling acceleration method

Abstract

The application provides a CFD-DEM coupling acceleration method 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 generates a probability distribution graph based on LSTM prediction of particle distribution; the CFD module combines flow field data to calculate a gradient, and utilizes a geometric deep learning model to divide an irregular dynamic coupling region; a proxy model based on GNN is constructed, and the proxy model and a traditional drag force model are adaptively switched based on a Reynolds number; the coupling controller monitors a calculation state, and triggers data synchronization of a hierarchical compression strategy when an asynchronous time step threshold is reached; the optimization module monitors momentum error after data synchronization, and feeds back adjustment of neural network parameters and a threshold coefficient if the threshold is exceeded; the application solves the problems of low calculation efficiency, rigid region division and lack of error control of traditional methods, significantly improves the calculation efficiency by 5-10 times, reduces the data synchronization amount by 58-65%, and ensures that the accuracy error is less than or equal to 5%.

Description

Technical Field

[0001] This invention relates to the field of coupled simulation technology of computational fluid dynamics (CFD) and discrete element method (DEM), and particularly to a CFD-DEM coupling acceleration method based on intelligent dynamic coupling region partitioning, adaptive surrogate model switching, and asynchronous time step collaborative optimization. This method significantly improves computational efficiency by 5-10 times while ensuring an accuracy error ≤5% through dynamic optimization of the coupling mechanism, breaking through the computational bottleneck of traditional synchronous coupling. This method is suitable for high-precision and efficient simulation of particle-fluid interactions, specifically applied to fields such as chemical reactor optimization, fluidized bed design, and particle transport system analysis. Background Technology

[0002] Particle-fluid multiphase flow systems are widely used in industrial scenarios such as energy, chemical, and pharmaceutical industries, and their accurate simulation is crucial for equipment design and process optimization. The CFD-DEM coupling method can reproduce 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 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 typically larger than 10 s. 2 ~10 4 (1) Forced synchronous coupling causes the DEM module to be in a waiting state for more than 90% of the time, with a computational resource utilization rate of less than 40%, especially when the local particle volume fraction is higher than 30%; (2) Rigid region partitioning: Fixed coupling region or uniform partitioning strategy cannot adapt to the dynamic aggregation / diffusion behavior of particles, resulting in redundant computation in low interaction regions; (3) Insufficient model adaptability: Traditional drag force models (such as Gidaspow, DiFelice) have decreased accuracy under transient high Reynolds number flow conditions, 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 computation is prone to accumulating phase error, and existing methods lack an adaptive feedback mechanism for coupling state, resulting in divergent simulation results.

[0003] While some improvements have been attempted in existing technologies, limitations remain. Asynchronous coupling strategies (such as the MPC algorithm) allow CFD / DEM to use different time steps, but the threshold setting relies on empirical parameters, which can easily lead to coupling instability during abrupt changes in the flow field (such as changes in particle inrush velocity >20%), resulting in accumulated errors in interphase forces exceeding 15%. Secondly, surrogate model applications are mostly limited to steady-state scenarios, failing to address the model switching problem in transient high Reynolds number flow states. Furthermore, coupling region optimization is typically based on regular geometric partitioning, making it difficult to identify unstructured particle aggregation patterns (such as chain clusters and sparse dispersions).

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

[0005] In summary, developing a coupling acceleration framework that integrates intelligent partitioning of dynamically coupled regions, adaptive switching of proxy models, asynchronous time-step collaborative optimization, and physical constraint embedding to achieve high-precision and high-efficiency CFD-DEM coupled simulation is urgently needed to overcome existing computational bottlenecks and meet the requirements of complex industrial applications. Summary of the Invention

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

[0007] To solve the above-mentioned technical problems, the present invention is implemented as follows:

[0008] This invention provides a CFD-DEM coupling acceleration method, comprising at least a CFD module, a DEM module, a coupling controller, and an optimization module. A dynamic coupling region is defined between the CFD module and the DEM module. Within the dynamic coupling region, the coupling controller monitors the computational state and determines whether an asynchronous time step threshold has been reached. If the asynchronous time step threshold is reached, data synchronization is performed. The optimization module monitors the coupling error after data synchronization and determines whether the coupling error has reached an error threshold. If the error threshold is reached, the optimization module adjusts the parameters of the neural network model and the threshold parameters. If neither the asynchronous time step threshold nor the error threshold is reached, the computation of the original CFD module and DEM module continues.

[0009] The steps for delineating the dynamic coupling region are as follows: First, input the initial flow field data of the CFD module and the particle attributes (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; next, in the CFD module, calculate the velocity gradient, pressure gradient, and turbulence intensity based on the flow field data, and determine the region where the particle concentration exceeds the threshold of the dynamic coupling region based on the probability distribution map; finally, use a geometric deep learning model to identify the particle aggregation pattern, delineate and output the irregularly shaped dynamic coupling region.

[0010] The neural network model includes at least one of the LSTM, the geometric deep learning model, and the 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 uses spatial hashing algorithm to cluster particles, ZFP lossless compression is used for key data, FP8 quantization is used for flow velocity in low gradient regions, and run-length encoding is used for static particle positions.

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

[0014] The training data for the LSTM, the geometric deep learning model, and the GNN are all derived from high-fidelity CFD-DEM coupled simulations, and data augmentation techniques are used to improve generalization.

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

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

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

[0018] The CFD-DEM coupled acceleration method described above is applicable to gas-solid fluidized beds, lithium battery slurry mixing tanks, and coal powder pneumatic conveying systems. It is especially suitable for highly unsteady working conditions where the dynamic change in particle concentration exceeds 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 surrogate model includes at least the following steps: in the CFD module, using particles as node features and particle-flow field mesh interaction information as edge features, inputting the dynamic coupling region, the particle attributes of the DEM module, and the flow field data; constructing the surrogate model based on GNN to replace the drag force model based on physical equations, and calculating the particle-flow field interphase force F. couplingThe surrogate model, when its local Reynolds number is in the low Reynolds number range, replaces the drag force model in calculating the interphase force F. coupling =F GNN When the local Reynolds number is in the high Reynolds number range, the system automatically switches from the surrogate model back to the drag force model based on the physical equations to calculate the interphase force F. coupling =F Drag In the transition Reynolds number range, the surrogate model and the drag force model are used in parallel to calculate the weighted and fused interphase force: F coupling =ω·F GNN +(1-ω)·F Drag Among them, weight ;in, The sigmoid activation function is used, and PMRM is the reliability index of the surrogate model. The calculation results of the surrogate model or the drag model are output to the DEM module. The DEM module continuously calculates the particle motion based on the calculation results of the surrogate model or the drag model until the coupling controller detects that the asynchronous time step threshold has been reached, and then performs the data synchronization.

[0021] The node features include at least the particle location. Particle speed Particle size d p Particle temperature T p This forms a 6-dimensional feature vector; the edge features include at least the local velocity difference. Grid pressure gradient Relative distance Constitutes a 4-dimensional edge feature; the particle temperature T p This is an optional input and is enabled when the simulation involves heat exchange.

[0022] The LSTM includes at least a spatiotemporal attention mechanism and a loss function; the loss function includes at least Newton's equation of motion; the spatiotemporal attention mechanism can focus on high gradient regions, avoiding resource waste caused by uniformly processing all data; the Newton's equation of motion is embedded as a physical constraint into the loss function to improve the reliability of single-particle motion.

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

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

[0025] The Hadamard product represents element-wise multiplication, adaptively weighting the particle position coordinates to suppress redundant data calculations in low-gradient regions.

[0026] The low gradient region refers to the region in the flow field where the velocity gradient, pressure gradient, and turbulence intensity are low, and the interaction between fluid and particles is weak in these regions.

[0027] The spatiotemporal attention mechanism, through attention weight a t Hide the historical state h t-1 Particle position Particle speed and local flow field velocity Dynamically focus on the high gradient region.

[0028] The formula for the Newtonian equation of motion includes at least the calculation methods for LSTM predicted displacement increment, Newton's law displacement increment, fused position prediction, and dynamic weights; the formula for the LSTM predicted displacement increment is as follows: in, The particle displacement increment at time step t is predicted by the LSTM unit based on historical states; the formula for the displacement increment according to Newton's law is: in, Δt represents the displacement increment at time step t predicted at time step t-1 based on Newton's second law. dem It is the time step of the DEM module; The acceleration is calculated from the sum of forces acting on the particle at time step t-1. The formula for the fusion location prediction is as follows: in, It refers to the particle position at time step t; η is the dynamic weight (0~1), used to balance the contributions of LSTM predicted displacement and Newton's law predicted displacement; the formula for calculating the dynamic weight is... in, It is the magnitude of the velocity gradient tensor; P represents the maximum value of the velocity gradient magnitude in the current flow field. density It is the normalized particle distribution probability (range 0 to 1).

[0029] The LSTM learns nonlinear motion patterns from historical states to predict paths; Newton's second law determines the particle displacement increment based on force analysis; the dynamic weight η increases the weight of the LSTM in the high gradient region, indicating a larger contribution of the LSTM to the predicted displacement, which shows a more significant nonlinear flow field effect; conversely, in the low gradient region, calculations rely on Newtonian mechanics.

[0030] The specific method by which the geometric deep learning model divides the dynamic coupling region is as follows: input the probability distribution map, velocity gradient, pressure gradient and turbulence intensity; analyze the local and global features of particle aggregation based on graph convolutional network (GCN); generate the dynamic coupling region according to the dense cluster, chain distribution and sparse dispersion pattern.

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

[0032] The neighborhood is defined as the particles within a spherical region centered on the target particle and with a radius of 2 to 3 times the particle diameter.

[0033] The radius setting of 2 to 3 times the particle size is based on the physical scale of particle-fluid interaction: this range can cover the particle wake effect and near-field drag interference zone, while avoiding the introduction of irrelevant noise by an excessively large radius.

[0034] The local features are the aggregation of the hidden states and spatial distance information of directly adjacent particles; the global features are the screening of potential long-range particle sets to capture non-local topological structures.

[0035] The surrogate 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 computation, enabling cross-scale coupled modeling and fully covering particle-flow field coupling effects; the mass-momentum conservation equation is embedded as a physical constraint in the loss function to improve the rationality and stability of gas-solid flow field simulation.

[0036] The formula for the multi-scale spatiotemporal aggregation includes a message passing function and node state updates; the formula for the message passing function is as follows: in, φ is the message passed from node j to node i at time step t; φ is an MLP function used for message generation and state update. and It represents the hidden states of nodes i and j at time step t; It is the Euclidean distance between nodes i and j; the formula for updating the node state is: in, ψ is the hidden state of node i after the update at time step t+1; ψ is an MLP function used for message generation and state update. Let i be the set of the direct neighbors of node i. Let i be the set of nodes (not directly adjacent) that have potential long-range interactions with node i, selected by the spatial hashing algorithm. It is a message passed from remote node k to node i.

[0037] The φ layer contains two fully connected layers (the number of neurons changes from 64 to 32), and the activation function is ReLU.

[0038] The mass-momentum conservation equation includes the total loss function, adaptive weights, and turbulence intensity; the formula for the total loss function is as follows: in, It is the total loss function; Mean squared error loss measures the deviation between the predicted value and the true value of the surrogate model; λ is the weighting coefficient of the physical constraint term; N cell The number of CFD mesh cells; The mass conservation equation (continuity equation) is given. For incompressible flow, ρ is a constant, and the equation simplifies to: This is the fluid inertia term (mass derivative); The pressure gradient term is represented by μ, where μ is the dynamic viscosity. For viscous forces; F coupling The interphase coupling force is calculated by the surrogate model or the drag force model; ρg is the gravity term; the formula for the adaptive weight is λ = 0.1 + 0.9·I t , among which, I t This is a turbulence intensity index; the formula for the turbulence intensity is: Where k is the turbulent kinetic energy; ||u|| is the magnitude of the local average velocity; and ε is a local constant, which is divided by zero.

[0039] Where ε is a local constant (typically taking the value of 10). -6 ε is used to prevent numerical calculation errors caused by a denominator of zero. When the local flow velocity ||u|| is close to zero (such as at stagnation points or in low Reynolds number regions), ε ensures the numerical stability of the formula and avoids abnormal program termination.

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

[0041] The formula for the reliability index PMRM of the proxy model is as follows: in, L is the feature length, which is set to particle size or mesh size depending on the simulation object; Let |u|| be the magnitude of the time derivative of the pressure gradient; ||u|| be the magnitude of the velocity vector; and f1(Re), f2(Re), and f3(Re) be dynamically adjusted adaptive coefficients.

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

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

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

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

[0046] Among them, the cubic cluster size of 2 to 3 times the particle size is a dedicated parameter for data compression. It matches the average spacing of typical particle systems, ensuring the spatial correlation of particle positions within the cluster, while avoiding excessively large clusters that could lead to a decrease in compression accuracy.

[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 encoding compression. The algorithm is dynamically selected according to the data type. For secondary data (such as flow field velocity and flow field temperature in the low gradient region), the FP16 half-precision, SZ compression, or FP8 precision quantization is used. For particle positions and stationary particle velocities outside the dynamic coupling region, the run-length encoding compression is used.

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

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

[0050] In the embodiments of the invention, by using LSTM to predict particle distribution and combining it with a geometric deep learning model to identify particle aggregation patterns, unstructured particle dynamic behavior can be adaptively identified, significantly reducing redundant computational regions by 30% to 50%, effectively overcoming the problems of excessively large coupling regions and high data synchronization caused by the use of regular boundary partitioning; a surrogate model constructed using GNN replaces the traditional drag force model for calculation under low Reynolds number flow conditions, and automatically switches back to the drag force model based on physical equations at high Reynolds numbers, so that the simulation accuracy error under transient high Reynolds number conditions is controlled within ≤5%, realizing the adaptive intelligent switching of the surrogate model and solving the bottleneck of low computational efficiency of the traditional drag force model; by using a layered compression strategy to significantly reduce the amount of data synchronization, and by using a coupling error feedback mechanism to dynamically adjust key threshold parameters (such as asynchronous time step threshold and dynamic coupling region threshold), the synchronization waiting time is shortened by 40% to 60%, while 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 present invention solve the problems of excessively large coupling area and large amount of data synchronization caused by the use of regular boundary division in traditional CFD-DEM coupling technology, slow calculation speed and low efficiency of traditional traction model, and lack of error control mechanism. Attached Figure Description

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

[0052] Figure 1 This diagram illustrates the structure of a CFD-DEM coupling acceleration method provided in an embodiment of the present invention.

[0053] Figure 2 This is a schematic diagram illustrating the basic process of a CFD-DEM coupling acceleration method provided in an embodiment of the present invention.

[0054] Figure 3 This is a flowchart illustrating the proxy model in a CFD-DEM coupling acceleration method provided in an embodiment of the present invention.

[0055] Explanation of reference numerals in the attached figures:

[0056] 10. CFD module; 20. DEM module; 30. Coupling controller; 40. Optimization module. Detailed Implementation

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

[0058] It should be understood that the phrase "one embodiment" or "an embodiment" throughout the specification means that a specific feature, structure, or characteristic related to the embodiment is included in at least one embodiment of the invention. Therefore, "in one embodiment" or "in an embodiment" appearing throughout the specification do not necessarily refer to the same embodiment. Furthermore, these specific features, structures, or characteristics can be combined in any suitable manner in one or more embodiments.

[0059] See Figures 1 to 3 This invention provides a CFD-DEM coupling acceleration method for computational fluid dynamics and discrete element method coupled 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 between the CFD module (10) and the DEM module (20); within the dynamic coupling region, the coupling controller (30) monitors the computational state and determines whether an asynchronous time step threshold has been reached, and performs data synchronization after the asynchronous time step threshold is reached; the optimization module (40) monitors the coupling error after the data synchronization, determines whether the coupling error has reached an error threshold, and if the error threshold is reached, feeds back to adjust the parameters of the neural network model and the coefficients of the threshold; if the asynchronous time step threshold or the error threshold is not reached, the original calculations of the CFD module (10) and the DEM module (20) continue to be executed.

[0060] The steps for dividing the dynamic coupling region are as follows: First, input 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), 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 (10); subsequently, in the CFD module (10), calculate the velocity gradient, pressure gradient and turbulence intensity in combination with the flow field data, and determine the particle regions that exceed the threshold of the dynamic coupling region based on the probability distribution map; finally, use a geometric deep learning model to identify the particle aggregation pattern, divide and output the irregularly shaped dynamic coupling region.

[0061] The neural network model includes at least one of the LSTM, the geometric deep learning model, and the 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 employs a layered compression strategy.

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

[0065] The CFD module (10) is used to solve the fluid control equations and calculate the interphase forces; the DEM module (20) tracks the particle motion trajectory based on Newton's laws of motion; the coupling controller (30) is used to dynamically monitor the calculation status, perform 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 adjust the neural network parameters and the threshold coefficients.

[0066] In this embodiment of the invention, by employing LSTM to predict particle distribution and combining it with a geometric deep learning model to identify particle aggregation patterns, unstructured particle dynamic behavior can be adaptively identified, significantly reducing redundant computational regions by 30% to 50%. This effectively overcomes the problems of excessively large coupling regions and high data synchronization requirements caused by using regular boundary partitioning. A surrogate model constructed using GNN replaces the traditional drag force model for calculations under low Reynolds number flow conditions and automatically switches back to the physical equation-based drag force model at high Reynolds numbers, keeping the simulation accuracy error under transient high Reynolds number conditions within ≤5%. This achieves adaptive intelligent switching of the surrogate model and solves the problems associated with traditional drag force models. This invention addresses the bottleneck of low computational efficiency in conventional CFD-DEM coupling technology. By employing a layered compression strategy to significantly reduce data synchronization, and utilizing a coupling error feedback mechanism to dynamically adjust key threshold parameters (such as asynchronous time step threshold and dynamic coupling region threshold), the synchronization waiting time is reduced by 40% to 60%, while the coupling momentum error is stably controlled within a preset threshold (≤3%). This effectively improves overall computational efficiency and strengthens error control capabilities. The embodiments of this invention solve the problems in prior art where traditional CFD-DEM coupling technology suffers from excessively large coupling regions, large data synchronization volumes, slow computation speed and low efficiency of traditional drag models, and the lack of error control mechanisms due to the use of regular boundary division.

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

[0068] Preferably, the construction of the surrogate model includes at least the following steps: in the CFD module (10), using particles as node features and particle-flow field mesh interaction information as edge features, inputting the dynamic coupling region, the particle attributes (position, velocity) of the DEM module (20), and the flow field data; constructing the surrogate model based on GNN to replace the drag force model based on physical equations, and calculating the particle-flow field interphase force F. coupling The surrogate model, when its local Reynolds number is in the low Reynolds number range, replaces the drag force model in calculating the interphase force F. coupling =F GNN When the local Reynolds number is in the high Reynolds number range, the system automatically switches from the surrogate model back to the drag force model based on the physical equations to calculate the interphase force F. coupling =F Drag In the transition Reynolds number range, the surrogate model and the drag force model are used in parallel to calculate the weighted and fused interphase force: F coupling =ω·F GNN +(1-ω)·F Drag Among them, weight ;in, The Sigmoid activation function is used, and PMRM is the reliability index of the surrogate model. The calculation results of the surrogate model or the drag model are output to the DEM module (20). The DEM module (20) continuously calculates the motion of the particles based on the calculation results of the surrogate model or the drag model until the coupling controller (30) detects that the asynchronous time step threshold has been reached and performs the data synchronization.

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

[0070] It should be noted that a multiphysics data bridging engine needs to be implemented in the coupling controller (30) to solve the compatibility problem of heterogeneous data between CFD and DEM software, including at least grid-particle mapping optimization and data type conversion; the grid-particle mapping optimization uses spatial hash grid to accelerate neighbor search, reducing the complexity from O(N) 2The data type conversion reduces the Eulerian field data (including but not limited to velocity / pressure tensors) of the CFD module (10) to the Lagrange particle force vector that the DEM module (20) can recognize; 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 equations of motion; the spatiotemporal attention mechanism can focus on regions with high velocity gradients, high pressure gradients, and high turbulence intensity, avoiding resource waste caused by uniformly processing all data; the Newton's equations of motion are embedded as physical constraints into the spatiotemporal attention mechanism to improve the reliability of single-particle motion.

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

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

[0074] It should be noted that directly inputting high-dimensional flow field data into the LSTM results in computational redundancy. Therefore, flow field gradient filtering is introduced, and this filtering is only applied when... The spatiotemporal attention mechanism is calculated in real time.

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

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

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

[0078] Preferably, the specific method by which the geometric deep learning model divides the dynamic coupling region is as follows: inputting the probability distribution map, velocity gradient, pressure gradient, and turbulence intensity; analyzing the local and global features of particle aggregation based on graph convolutional networks (GCN); and generating the dynamic coupling region according to dense clusters, chain distributions, 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 bubbling zone and bubble tail chains in the fluidized bed.

[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 with a radius of 3 times is selected.

[0083] Preferably, the surrogate 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 computation, enabling cross-scale coupled modeling and fully covering particle-flow field coupling effects; the mass-momentum conservation equation is embedded as a physical constraint in the multi-scale spatiotemporal aggregation to improve the rationality and stability of gas-solid flow field simulation.

[0084] Preferably, the formula for the multi-scale spatiotemporal aggregation includes a message passing function and node state updates; the formula for the message passing function is as follows: in, φ is the message passed from node j to node i at time step t; φ is an MLP function used for message generation and state update. and It represents the hidden states of nodes i and j at time step t; It is the Euclidean distance between nodes i and j; the formula for updating the node state is: in, ψ is the hidden state of node i after the update at time step t+1; ψ is an MLP function used for message generation and state update. Let i be the set of the direct neighbors of node i. Let i be the set of nodes (not directly adjacent) that have potential long-range interactions with node i, selected by the spatial hashing algorithm. It is a message passed from remote node k to node i.

[0085] Preferably, the formula for the mass-momentum conservation equation includes a total loss function, adaptive weights, and turbulence intensity; the formula for the total loss function is as follows: in, It is the total loss function; Mean squared error loss measures the deviation between the predicted value and the true value of the surrogate model; λ is the weighting coefficient of the physical constraint term; N cell The number of CFD mesh cells; The mass conservation equation (continuity equation) is given. For incompressible flow, ρ is a constant, and the equation simplifies to: This is the fluid inertia term (mass derivative); The pressure gradient term is represented by μ, where μ is the dynamic viscosity. For viscous forces; Fcoupling The interphase coupling force is calculated by the surrogate model or the drag force model; ρg is the gravity term; the formula for the adaptive weight is λ = 0.1 + 0.9·I t , among which, I t This is a turbulence intensity index; the formula for the turbulence intensity is: Where k is the turbulent kinetic energy; ||u|| is the magnitude of the local average velocity; and ε is a local constant, which is divided by zero.

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

[0087] It should be noted that the grid cell TH region The dynamic coupling region is marked as a candidate when the value is greater than 0.4.

[0088] Preferably, the formula for the reliability index PMRM of the proxy model is as follows: in, L is the feature length, which is set to particle size or mesh size depending on the simulation object; Let |u|| be the magnitude of the time derivative of the pressure gradient; ||u|| be the magnitude of the velocity vector; and f1(Re), f2(Re), and f3(Re) be dynamically adjusted adaptive coefficients.

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

[0090] It should be noted that TH asyncThe data synchronization is performed when the value is greater than 0.35.

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

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

[0093] It should be noted that F DEM,i The DEM module (20) calculates the forces actually used in previous calculation steps (calculated and passed by the CFD module (10) at an earlier time). CFD,i At this data synchronization moment, based on the latest flow field state solved by the CFD module (10) and the latest position of particle i, the interphase force is recalculated (using the same rules of the surrogate model or the drag force model); the momentum error function E m This reflects the cumulative degree of deviation in the calculation of interphase forces caused by the asynchronous calculation between 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, initiate the optimization process of adjusting the parameters and threshold coefficients of the feedback 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 standardized interfaces (such as UDF and API).

[0096] It should be noted that, in this document, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or apparatus. Unless otherwise specified, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or apparatus that includes that element.

[0097] Through the above description of the embodiments, those skilled in the art can clearly understand that the methods of the above embodiments can be implemented by means of software plus necessary general-purpose hardware platforms. Of course, they can also be implemented by hardware, but in many cases the former is a better implementation method. 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. This computer software product is stored in a storage medium (such as ROM / RAM, magnetic disk, optical disk) and includes several instructions to cause a terminal (which may be a mobile phone, computer, server, air conditioner, or network device, etc.) to execute the methods described in the various embodiments of the present invention.

[0098] The embodiments of the present invention have been described above with reference to the accompanying drawings. However, the present invention is not limited to the specific embodiments described above. The specific embodiments described above are merely illustrative and not restrictive. Those skilled in the art can make many other forms under the guidance of the present invention without departing from the spirit and scope of the claims, and all of these forms are within the protection scope of the present invention. Example 1

[0099] To more clearly illustrate the technical implementation of this invention, the following description is in conjunction with the appendix. Figures 1-3 The specific implementation steps of the present invention will be explained 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 (1.2m height, 0.02m bottom diameter, 0.05m top diameter) was constructed in Fluent 2019 using an unstructured hexahedral mesh (with boundary layer refinement), totaling 2.15 million meshes. Incompressible nitrogen gas was used, with a density of 1.225 kg / m³. 3 The dynamic viscosity is 1.7894 × 10⁻⁶. -5The velocity inlet was 6 m / s, the pressure outlet gauge pressure was 0 Pa, and the wall surface was a no-slip boundary. EDEM2022 was used to simulate the particulate phase. The particulate system contained 500,000 spherical particles to simulate carbon nanotube aggregates. The mixed particles had a dual-size distribution: small particles with a diameter of 150 μm (60%) and large particles with a diameter of 300 μm (40%). The particle density was 150 kg / m³. 3 Poisson's ratio is 0.25, and shear modulus is 1×10⁻⁶. 10 The contact model between particles and between particles and the wall was modeled using Hertz-Mindlin with JKR. Based on the particle collision timescale analysis, the time step of the DEM was set to 1×10⁻⁶. -6 Based on CFL conditions and flow field stability requirements, the CFD time step is set to 1×10⁻⁶. -5 s.

[0101] Secondly, the dynamically coupled regions are intelligently partitioned. A two-layer LSTM network (each layer with 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 time. Applying spatiotemporal attention mechanisms The network output is a probability distribution map P of particle existence on an 800×800 two-dimensional grid covering the entire CFD computation domain. density (Each grid value ranges from [0, 1]). The global velocity gradient tensor is computed at each time step in the Fluent UDF. and pressure gradient Then calculate the threshold of the dynamic coupling region. After parameter scanning optimization, the following approach was determined: First, with β = 0.25 and γ = 0.15 fixed, a scan was performed within the range of α ∈ [0.5, 0.7] with a step size of 0.01. The minimum error was found when α = 0.6. Then, with α = 0.6 and γ = 0.15 fixed, a scan was performed within the range of β ∈ [0.2, 0.3], and β = 0.25 was found to be optimal. Finally, with α = 0.6 and β = 0.25 fixed, a scan was performed within the range of γ ∈ [0.1, 0.2], and γ = 0.15 was found to be optimal. In this case, α = 0.6, β = 0.25, and γ = 0.15 were ultimately used to calculate TH for each grid cell. region The previous step marked more than TH regionIn the region, a geometric deep learning model, such as the PointNet++ model (containing 4 ensemble abstraction layers and 2 feature propagation layers), is 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) are extracted respectively. The α-shape algorithm is applied to generate irregular, smooth convex hull boundaries, which are the defined dynamic coupling regions. Sparse dispersion regions (label 0) are excluded from the coupling calculation. In the simulation, this strategy accurately identifies the dense bubbling region at the bottom of the fluidized bed, the chain-like particle band of the rising bubble tail vortex, and the dilute phase region at the top. The dynamic coupling region accounts for only 28-35% of the original computational domain on average, and the data synchronization amount is reduced by more than 40%.

[0102] Subsequently, a proxy model is constructed within the dynamically coupled region, and the model is adaptively switched. Within the dynamically coupled region, each particle is treated as a graph node, and its node feature is the particle's position. Particle speed Particle size d p Particle temperature T p (6-dimensional), its edge feature is local velocity difference Grid pressure gradient relative distance (4-dimensional). Construct a 4-layer dynamic graph neural network, each layer containing: a message passing function. (φ: 2-layer MLP, 64→32 neurons, ReLU), node update function (ψ: 2-layer MLP, 64→64 neurons, ReLU), 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 a surrogate model to calculate F. GNN As F coupling When Re > 2000, switch to the Gidaspow traction model to calculate F. Drag As F coupling When 500≤Re≤2000, the proxy model and the drag model are computed in parallel, and the final force is the weighted fusion F. coupling =ω·F GNN +(1-ω)·F Gidaspow Among them, weight in, Let be the Sigmoid activation function, and PMRM be the reliability metric of the surrogate model; the formula for the reliability metric of the surrogate model is: in, L is the feature length, which is set to particle size or mesh size depending on the simulation object; Let |u|| be the magnitude of the time derivative of the pressure gradient; ||u|| be the magnitude of the velocity vector; and f1(Re), f2(Re), and f3(Re) be dynamically adjusted adaptive coefficients. A nonlinear weighting function with a Sigmoid transition is used at the boundary where Re = 500 / 2000. The calculated F... coupling Passed to EDEM 2022 for particle motion solving.

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

[0104] Finally, error feedback control is implemented, adjusting the neural network parameters and threshold coefficients based on the data synchronization results. After each data synchronization, the optimization module calculates the momentum error. (Normalized relative error), F DEM,i The interphase interaction 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,iTo recalculate the interphase force at the current position of particle i in the CFD module, based on the latest flow field, using either a surrogate model or a drag model. The calculation employs inverse distance weighted interpolation, searching for the interphase force field on the CFD mesh cell containing particle i and its eight adjacent nodes. If the normalized momentum error E... m When the value is greater than 0.03, the optimization module initiates the adjustment process. First, based on the current error E... m Given the flow / particle field conditions (e.g., mean Reynolds number Re, mean particle concentration), use small-step gradient descent or rule-based empirical tuning, such as employing a gradient-based Adam optimizer (learning rate 10). -4 Fine-tune the asynchronous threshold coefficients k1 and k2 and the region threshold weights α, β, and γ. For example, if the error mainly stems from untimely synchronization in high-collision regions, slightly increase k2; if the error stems from insufficient partitioning of high-gradient regions, slightly increase β. Secondly, while retaining the weights of the main LSTM network, utilize only the particle motion trajectory data from the most recent few time steps, and use Newton's equations of motion as physical constraints to adjust the parameter matrix W of the spatiotemporal attention layer in the LSTM. a Fine-tuning is performed, with the learning rate set to one-tenth of the initial training rate. The training objective is the position prediction error over two CFD steps, while satisfying the weak constraint of Newton's second law. The updated parameters (k1, k2, α, β, γ) and the fine-tuned LSTM attention layer weights W are then used. a It is applied to subsequent simulation processes.

[0105] The above Example 1 was computed on a single-node computing platform AMD EPYC 9754 (128 cores), 512GB DDR5 RAM, and NVIDIA Tesla A100, and quantitatively compared with the traditional CFD-DEM coupling method under the same conditions. Example 1 significantly outperforms the traditional CFD-DEM coupling method in terms of computation time, CPU utilization, and data synchronization. Specifically, Example 1's computation time is 8 times faster than the traditional CFD-DEM coupling method, indicating a significant improvement in computational efficiency and fully demonstrating the effectiveness of this invention in solving the problems of low efficiency, rigid region partitioning, insufficient model adaptability, and lack of error control in traditional coupling methods. The performance comparison between the traditional CFD-DEM coupling method and the method of this invention is shown in the table below:

Claims

1. A CFD-DEM coupled acceleration method, applied to the field of coupled simulation of computational fluid dynamics and discrete element method, characterized in that, The system is designed to accelerate the computation of a CFD-DEM coupled model, 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 defined between the CFD module (10) and the DEM module (20). Within the dynamic coupling region, the coupling controller (30) monitors the computation status and determines whether an asynchronous time step threshold has been reached. Data synchronization is performed after the asynchronous time step threshold is reached. The optimization module (40) monitors the coupling error after data synchronization and determines whether the coupling error has reached an error threshold. If the error threshold is reached, the system feeds back and adjusts the parameters of the neural network model and the threshold parameters. If the asynchronous time step threshold or the error threshold is not reached, the computation of the original CFD module (10) and DEM module (20) continues. The steps for dividing the dynamic coupling region are as follows: First, input the initial flow field data of the CFD module (10) and the particle position and velocity of the DEM module (20); then, in the DEM module (20), 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 (10); subsequently, in the CFD module (10), calculate the velocity gradient, pressure gradient and turbulence intensity in combination with the flow field data, and determine the particle regions that exceed the threshold of the dynamic coupling region based on the probability distribution map; finally, use a geometric deep learning model to identify the particle aggregation pattern, divide and output the irregularly shaped dynamic coupling region. The neural network model includes at least one of the LSTM, the geometric deep learning model, and the 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 employs a layered compression strategy; The flow field data includes at least the fluid velocity vector, pressure, and turbulence intensity.

2. The CFD-DEM coupling acceleration method according to claim 1, characterized in that, The CFD module (10) includes at least an agent model; The construction of the surrogate model includes at least the following: in the CFD module (10), using particles as node features and particle-flow field mesh interaction information as edge features, inputting the dynamic coupling region, the particle position and velocity of the DEM module (20), and the flow field data; constructing the surrogate model based on GNN to replace the drag force model based on physical equations, and calculating the particle-flow field interphase force. The local Reynolds number of the surrogate model is in the low Reynolds number range. The time is used to replace the traction force model to calculate the interphase force. When the local Reynolds number is in the high Reynolds number range When the time comes, the system automatically switches from the proxy model back to the drag force model based on physical equations to calculate the interphase forces. In the transition Reynolds number interval The surrogate model and the drag force model are used in parallel for weighted fusion calculation to calculate the interphase forces: Among them, weight ;in, It is the Sigmoid activation function. The reliability index of the proxy model is used; 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 motion of the particles based on the calculation results of the proxy model or the drag model until the coupling controller (30) detects that the asynchronous time step threshold has been reached and then performs the data synchronization.

3. The 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 regions with high velocity gradients, high pressure gradients, and high turbulence intensity, avoiding resource waste caused by uniformly processing all data; the Newton's equations of motion are embedded as physical constraints into the spatiotemporal attention mechanism to improve the reliability of single-particle motion; The formula for the spatiotemporal attention mechanism includes at least an attention weight vector, an adjusted input vector, and an LSTM state update; the formula for the attention weight vector is: ,in, In time step Attention weight vector; Use the Sigmoid activation function; The learnable parameter matrix; For LSTM cells at time steps The hidden state of history; For particles at time step The position coordinate vector; For particles at time step The velocity vector; In time step The velocity vector of the fluid at the location of the particle; the formula for the adjusted input vector is... ,in, Attention weight Weighted adjustment of particle position input vector; Represents the Hadamard product; the formula for the LSTM state update is: ,in, The updated LSTM hidden state; The formula for the Newtonian equation of motion includes at least the calculation methods for LSTM predicted displacement increment, Newton's law displacement increment, fused position prediction, and dynamic weights; the formula for the LSTM predicted displacement increment is as follows: ,in, The prediction at time step by the LSTM unit based on historical states The particle displacement increment; the formula for the displacement increment according to Newton's law is: ,in, Based on Newton's second law at time step Predicted at time step The displacement increment; It is the time step of the DEM module (20); It is the particles at time step The acceleration is calculated from the sum of the forces acting on it. The formula for the fusion location prediction is: ,in, In time step The position of the particles; The dynamic weights, ranging from 0 to 1, are used to balance the contributions of LSTM-predicted displacement and Newton's law-predicted displacement; the formula for calculating the dynamic weights is as follows: ,in, It is the magnitude of the velocity gradient tensor; This represents the maximum value of the velocity gradient magnitude in the current flow field; It is the normalized particle distribution probability.

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

5. The CFD-DEM coupling acceleration method according to claim 2, characterized in that, The surrogate model includes at least a multi-scale spatiotemporal aggregation and a loss function; the loss function includes at least a mass-momentum conservation equation; the multi-scale spatiotemporal aggregation is achieved through GCN computation, enabling cross-scale coupled modeling and fully covering particle-flow field coupling effects; the mass-momentum conservation equation is embedded as a physical constraint into the multi-scale spatiotemporal aggregation to improve the rationality and stability of gas-solid flow field simulation. The formula for the multi-scale spatiotemporal aggregation includes a message passing function and node state updates; the formula for the message passing function is as follows: ,in, In time step From node Passed to node The message; These are MLP functions used for message generation and state updates; and It is a node and At time step The hidden state; It is a node and The Euclidean distance between them; the formula for updating the node state is: ,in, It is a node At time step The updated state is now hidden. These are MLP functions used for message generation and state updates; For nodes The set of directly adjacent nodes; For nodes selected through spatial hashing algorithm A set of nodes with potential long-range interactions; From remote node Passed to node The message; The mass-momentum conservation equation includes the total loss function, adaptive weights, and turbulence intensity; the formula for the total loss function is as follows: ,in, It is the total loss function; Mean squared error loss measures the deviation between the predicted value and the true value of the surrogate model. The weighting coefficients for the physical constraint terms; The number of CFD mesh cells; For the mass conservation equation, for incompressible flow, Since is a constant, the equation simplifies to ; This is the fluid inertia term; The pressure gradient term is represented by μ, where μ is the dynamic viscosity. It is a viscous force term; The interphase coupling force is calculated by the surrogate model or the drag force model; For gravity; the formula for the adaptive weight is: ,in, This is a turbulence intensity index; the formula for the turbulence intensity is: ,in, It is turbulent kinetic energy; This represents the magnitude of the local average flow velocity. It is a very small constant, so we avoid dividing by zero.

6. The CFD-DEM coupling acceleration method according to claim 1, characterized in that, The formula for the threshold of the dynamic coupling region is: ,in, It is the threshold of the dynamic coupling region; It is the normalized particle distribution probability; The magnitude of the velocity gradient tensor; The magnitude of the pressure gradient vector; , , The weighting coefficients are and satisfy the following conditions: .

7. The CFD-DEM coupling acceleration method according to claim 2, characterized in that, The reliability index of the proxy model The formula is ,in, L is the feature length, which is set to particle size or mesh size depending on the simulation object; Let be the modulus of the time derivative of the pressure gradient; The magnitude of the flow velocity vector; , , For dynamically adjusted adaptive coefficients, .

8. The CFD-DEM coupling acceleration method according to claim 1, characterized in that, The asynchronous time step threshold is ,in, It is the threshold value for the asynchronous time step; For Courant-Friedrichs-Lewy numbers, , must meet , It is the maximum flow rate. For the current time step, The feature size of the CFD mesh; This represents the total number of particle collisions that occurred within the most recent DEM coupling step. The characteristic time of a single particle collision; and These are the empirical coefficients determined through pre-calibration experiments. .

9. The 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; and using a spatial hashing algorithm to cluster particles into clusters, with the cluster side length being 2 to 3 times the particle diameter.

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... ,in, It is the momentum error function; The total number of particles; The effect of the DEM module (20) on the particles is calculated. Interphase forces on; To use the inverse distance weighted interpolation method, based on the particle size distribution... The effect of interphase force interpolation on the surrounding CFD mesh nodes on the particles The force exerted by position.

Images