A multi-physics coupling digital twin simulation and prediction method for industrial equipment

Abstract

The application discloses a kind of multi-physics field coupling digital twin simulation and prediction method of industrial equipment, comprising: obtaining the multi-physics field full-order numerical solution dataset of complex industrial equipment under multiple working conditions;After average field decomposition, low-dimensional modal base and modal coefficient are extracted using intrinsic orthogonal decomposition;Select representative physical quantity as the basis of first-level partition, then divide the three-dimensional calculation domain into several first-level subdomains;Construct secondary clustering feature vector in each first-level subdomain and cluster again;Adaptively determine the optimal clustering number of each secondary subdomain based on elbow method and improvement rate threshold method;In each secondary subdomain, different physical quantities and different modal orders are respectively established proxy model;Input new working condition parameters to each proxy model, predict modal coefficients and reconstruct full-order physical field distribution.The application significantly improves the accuracy and efficiency of multi-physics field prediction by hierarchical partition modeling driven by complex industrial equipment, meets the needs of real-time simulation of digital twin.

Description

Technical Field

[0001] This invention relates to the field of multiphysics simulation and digital twin technology, and more specifically to a multiphysics coupled digital twin simulation and prediction method for industrial equipment. Background Technology

[0002] Complex industrial equipment, such as large boilers, heavy-duty gas turbines, high-energy-density battery packs, and process systems, typically involve multiple coupled physical processes during operation, including fluid flow, heat and mass transfer, component transport, chemical reactions, and structural responses. To obtain the spatial distribution of state variables such as temperature, velocity, pressure, component concentration, and stress within the equipment, current technologies typically employ numerical methods such as computational fluid dynamics and finite element analysis to establish high-fidelity multiphysics coupled models. While these models can accurately describe the internal state of the equipment, they are usually characterized by a large number of meshes, large computational scale, and long solution times. When conducting multi-condition analysis, parameter studies, and state prediction, it is necessary to repeatedly solve the high-fidelity model, resulting in high computational costs.

[0003] Digital twins establish virtual models corresponding to physical entities in an information space and receive the physical entities' operating parameters or measurement data to achieve equipment state mapping, performance prediction, and operational status analysis. For industrial equipment with complex internal flow, heat transfer, and reaction processes, it is difficult to fully reflect the internal state of the equipment using only a small amount of measurement data. Therefore, it is necessary to introduce multiphysics models into the digital twin system. However, the computation time of high-fidelity multiphysics models is often insufficient to meet the real-time or near-real-time response requirements of digital twins; if only an empirical mapping between operating parameters and a small number of overall performance indicators is established, it is difficult to recover the complete physical field and accurately characterize local high-gradient regions and the spatial correlation between different physical fields.

[0004] Model reduction and surrogate models offer feasible approaches to improving the computational efficiency of digital twins. Proper Orthogonal Decomposition (POD) can extract key features from high-fidelity samples, transforming high-dimensional physical fields into low-dimensional representations; surrogate models can establish mapping relationships between operating parameters and low-dimensional features, enabling prediction of the physical field. However, the physical characteristics and variation patterns of different regions within complex industrial equipment vary significantly, making it difficult for reduced-order models to simultaneously guarantee the representation accuracy of each local region. For complex multiphysics fields with significant spatial differences, local reduction methods can enhance the feature representation capabilities of different regions, but the method of partitioning the computational domain, the number of local models, and the processing of results between subdomains all significantly impact model performance.

[0005] Therefore, there is a need to provide a multi-physics digital twin simulation method for complex industrial equipment. While preserving high-fidelity multi-physics information, this method combines physical mechanism partitioning with data-driven partitioning, rationally determines different local modeling regions during the model reduction stage, and restores the complete physical field through hierarchical splicing. Furthermore, it combines a surrogate model to achieve state prediction of complex industrial equipment under different operating conditions. Summary of the Invention

[0006] The purpose of this invention is to provide a multi-physics coupled digital twin simulation and prediction method for industrial equipment, in order to solve the problems of high computational cost of existing high-fidelity multi-physics models for complex industrial equipment, insufficient local representation capability of reduced-order models, difficulty in adapting fixed partitions to different regions of complexity, and difficulty in ensuring the prediction accuracy of complete physical fields in digital twins.

[0007] To achieve the above objectives, the present invention adopts the following technical solution:

[0008] This invention provides a multi-physics coupled digital twin simulation and prediction method for industrial equipment, comprising: Step S1: Determine the operating parameters and their value range of the target complex industrial equipment, generate multi-condition samples using experimental design methods, and divide the multi-condition samples into training conditions and independent test conditions. Step S2: Based on the geometric structure, material properties, boundary conditions and operating mechanism of the target complex industrial equipment, establish a high-fidelity multiphysics coupling model, perform numerical solutions for each sample working condition, and obtain a multiphysics snapshot database. Step S3: Based on the multi-physics snapshot database of the training conditions, perform global intrinsic orthogonal decomposition on each physical quantity, determine the number of modes to be retained according to the preset cumulative energy contribution rate, and establish a global order reduction model. Step S4: Determine at least one primary partitioning direction based on the geometric structure and physical field distribution characteristics of the target complex industrial equipment, and divide the computational domain into multiple primary core subdomains according to the gradient change of the selected representative physical quantity along the direction. Step S5: Set up a first-level overlapping region between adjacent first-level core subdomains, and construct complementary continuous weights based on the normalized spatial position of grid points or grid cells in the first-level overlapping region. Step S6: Construct a first-level local training snapshot matrix for each first-level subdomain and each physical quantity, perform intrinsic orthogonal decomposition, determine the number of retained modes for each first-level subdomain according to the preset cumulative energy contribution rate, and use the complementary continuous weights to fuse the local physical fields of each first-level subdomain to establish an intrinsic orthogonal decomposition model for the first-level subdomain. Step S7: Within each primary core subdomain, a secondary partitioning is performed using a multiphysics data-driven method. Joint clustering features are constructed using spatial coordinates and the statistical characteristics of physical quantities under training conditions. Physical quantities with large magnitude differences are transformed and standardized. A clustering method with spatial neighborhood constraints is used to partition the secondary subdomains. Step S8: Set multiple candidate second-level subdomains for each first-level core subdomain, and adaptively determine the number of second-level subdomains for each first-level core subdomain based on the local order reduction error and partitioning results under internal verification conditions. Step S9: Construct a second-level local training snapshot matrix for each second-level subdomain and each physical quantity, calculate the local training mean field, perform eigenorthogonal decomposition on the snapshot matrix after removing the local training mean field, and determine the number of modes to be retained based on the preset cumulative energy contribution rate. Step S10: The second-level subdomains are assigned non-overlapping weights. The weights of the second-level subdomains are multiplied by the complementary continuous weights of the corresponding first-level subdomains to form hierarchical weights. The hierarchical weights are then used to accumulate the local physical fields of each second-level subdomain to recover the complete multiphysics field and establish the intrinsic orthogonal decomposition model of the second-level subdomains. Step S11: The intrinsic orthogonal decomposition models of the whole domain, first-level subdomain, and second-level subdomain are compared using unified evaluation conditions, and the target reduced-order model for prediction is determined by combining model error and scale. Step S12: Using the operating parameters of the target complex industrial equipment as input and the modal coefficients of the target reduced-order model as output, establish multiple candidate surrogate models respectively, and use cross-validation method to evaluate the predictive performance of each candidate surrogate model to determine the mapping relationship between operating parameters and modal coefficients. Step S13: Input the operating parameters to be predicted into the surrogate model to obtain the corresponding modal coefficients, and reconstruct the complete multiphysics field by combining the mean field and the reduced-order basis to build a multiphysics digital twin simulation model for predicting the operating status of complex industrial equipment.

[0009] In one embodiment, in step S4, the first-level partition direction is determined based on the geometry and main physical processes of the target complex industrial equipment, and the representative physical quantity is a physical quantity closely related to the main physical process; the gradient change is determined by calculating the spatial gradient of the representative physical quantity along the direction under each training condition and taking the aggregated gradient, and the aggregated gradient is the maximum value or statistical quantile of the absolute value of the gradient under each training condition.

[0010] In one embodiment, in step S5, the complementary continuous weights are constructed using a cosine function: let the overlapping interval be [ , The normalized position corresponding to grid position y is:

[0011] In the formula, The left boundary of the overlapping interval. The right boundary of the overlapping interval; the complementary continuous weights of the left and right first-level subdomains are as follows:

[0012]

[0013] In the formula, The normalized coordinates of the grid positions within the overlapping interval. The complementary continuous weights of the left-hand first-order subdomain, Let be the complementary continuous weights of the right-hand first-order subdomain; where satisfies .

[0014] In one embodiment, in step S7, the joint clustering features include the three-dimensional spatial coordinates of the grid cells and the mean and standard deviation of each physical quantity under the internal fitting condition; for physical quantities with skewed distribution characteristics, a logarithmic transformation is performed on them before constructing the clustering features; the spatial neighborhood constraint is achieved by calculating the proportion of each candidate cluster in the spatial nearest neighbors of each grid cell and introducing a regularization term to correct the clustering cost function.

[0015] In one embodiment, the cost function is:

[0016] In the formula, Number the global grid cells. Candidate tags The corresponding evaluation cost, For the first Clustering feature vectors of individual grid cells; For the first Cluster centers, For the feature dimension, For grid cells j Spatial neighbors are tagged as c proportion, is the space regularity coefficient.

[0017] In one embodiment, step S8, the adaptive determination of the number of second-level subdomains, includes: For each primary core subdomain, the number K of candidate secondary subdomains is set to a value within a preset range; The training scenarios are divided into internal fitting scenarios and internal validation scenarios. For each candidate K, a local intrinsic orthogonal decomposition is performed and the order reduction error of the internal validation scenario is calculated. Based on the curve of the reduction error as a function of K, the optimal K value is automatically selected using the elbow rule or the improvement rate threshold method.

[0018] In one embodiment, the improvement rate threshold method includes: calculating the total improvement rate R from the smallest candidate feasible K to the largest candidate feasible K; selecting the smallest feasible K when R is less than a first threshold; otherwise, calculating the elbow score of the normalized error curve; selecting the K with the largest elbow score when the maximum elbow score is greater than or equal to a second threshold; otherwise, selecting the K that first reaches the preset error improvement ratio.

[0019] In one embodiment, in step S10, the hierarchical weight is:

[0020] in, For hierarchical weights, For the first Each first-level subdomain in the grid cell complementary continuous weights at the point, This is the attribution indicator function for the second-level subdomain, and it satisfies the condition for all grid cells. have , This represents the total number of first-level subdomains. This represents the number of second-level subdomains within the first-level subdomain p.

[0021] In one embodiment, in step S12, the candidate surrogate model includes one or more of Gaussian process regression, radial basis function, support vector regression, random forest, and polynomial response surface; different surrogate models are selected for different subdomains, different physical quantities, or different modal orders.

[0022] In one embodiment, in step S2, the multiphysics field includes one or more of a velocity field, a temperature field, a component concentration field, and a pressure field.

[0023] As can be seen from the above technical solution, compared with the prior art, the present invention has the following technical advantages: 1. Two-level partitioning strategy improves local representation accuracy: The macroscopic first-level partitioning driven by physics is based on gradient change, while the microscopic second-level partitioning driven by data is based on spatial-physical property joint clustering. The two-level partitioning makes the physical field change more uniform in each local subdomain. The intrinsic orthogonal decomposition can capture a higher proportion of energy with fewer modes, thereby significantly improving the prediction accuracy of the reduced-order model in the whole domain.

[0024] 2. Adaptive mechanism for determining the number of secondary subdomains: An adaptive method combining internal validation with the elbow rule and improvement rate threshold is adopted to automatically select the optimal number of secondary subdomains for each primary core subdomain. This avoids the subjectivity of manual trial and error and the insufficient adaptability of fixed partitions in different regions, achieving a reasonable balance between model accuracy and complexity.

[0025] 3. Hierarchical weight splicing ensures the continuity of the entire field: The first-level overlapping area adopts cosine complementary weight to achieve a smooth transition, and the second-level subdomain adopts hard subdomain and non-overlapping division. The hierarchical weight formed by the product of the two ensures that the weight sum of any grid point in the entire field is 1, thereby making the reconstructed physical field naturally continuous and eliminating artificial discontinuities at the partition boundary.

[0026] 4. Highly efficient digital twin prediction capability: Through a phased strategy of offline training and online prediction, complete multiphysics prediction at the millisecond to second level is achieved, meeting the real-time or near-real-time requirements of digital twins.

[0027] 5. Good versatility and scalability: This method does not depend on the physical model of specific industrial equipment and can be extended to various complex systems with multi-physics coupling, such as boilers, gas turbines, chemical reactors, and battery thermal management systems. Attached Figure Description

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

[0029] Figure 1 A schematic diagram of the overall process of the multi-physics field coupled digital twin simulation and prediction method for industrial equipment provided in the embodiments of the present invention.

[0030] Figure 2 This is a three-dimensional model diagram of the boiler in an embodiment of the present invention.

[0031] Figure 3 This is a schematic diagram of dividing a first-level subdomain based on the change in the sensitive direction of a representative physical quantity in an embodiment of the present invention.

[0032] Figure 4 This is a schematic diagram of dividing secondary subdomains based on spatial joint clustering features in an embodiment of the present invention.

[0033] Figure 5a This is the original simulation calculation cloud map in the embodiment of the present invention.

[0034] Figure 5b This is a global POD reconstruction cloud map in an embodiment of the present invention.

[0035] Figure 5c This is a reconstruction cloud map of the first-level subdomain POD in an embodiment of the present invention.

[0036] Figure 5d This is a reconstruction cloud map of the second-level adaptive subdomain POD in an embodiment of the present invention. Detailed Implementation

[0037] 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 embodiments of the present invention, and not all embodiments. 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.

[0038] Example 1: This invention discloses a multi-physics coupled digital twin simulation and prediction method for industrial equipment, referring to... Figure 1 As shown, steps S1 to S13 are included: Step S1: Determine the operating parameters and their value range of the target complex industrial equipment, generate multi-condition samples using experimental design methods, and divide the multi-condition samples into training conditions and independent test conditions. This embodiment takes a 330MW coal-fired boiler as the research object, and uses boiler load, excess air coefficient, and burnout air rate as operating parameters. The parameter vector for each training scenario is represented as follows:

[0039] In the formula, For the first The parameter vector for each training scenario. For boiler load, The excess air coefficient, To burn out the wind rate, This indicates transpose.

[0040] The boiler load varied from 30% to 100% of the rated load, the excess air coefficient varied from 1.10 to 1.40, and the burnout air rate varied from 0.15 to 0.25. Twenty-five operating conditions were generated using the Latin hypercube experimental design method, of which 21 were used as training conditions and 4 as independent test conditions. The training conditions were used for subsequent first-level subdomain partitioning, second-level subdomain partitioning, intrinsic orthogonal decomposition mode extraction, and subsequent surrogate model training; the independent test conditions did not participate in determining the mean field, modal basis, and the number of adaptive second-level subdomains, but were only used to evaluate the model's predictive ability for unknown operating conditions.

[0041] Step S2: Based on the geometric structure, material properties, boundary conditions and operating mechanism of the target complex industrial equipment, establish a high-fidelity multiphysics coupling model, perform numerical solutions for each sample working condition, and obtain a multiphysics snapshot database. Based on the structure of the boiler furnace, burners, burnout air nozzles, and heating surfaces, a three-dimensional computational domain is established using three-dimensional modeling software, referencing... Figure 2 The air inlets for the main burner zone and the burnout air zone are symmetrically distributed at the four corners, including five layers of primary air inlets (A~E) and seven layers of secondary air inlets (AA~FF) interleaved, while four layers of burnout air inlets (SOFA1~SOFA4) are separately arranged in the upper part of the furnace. Local structures that do not affect the main flow and combustion processes within the furnace are appropriately simplified, and the computational domain is meshed. The same mesh is used for solving all operating conditions.

[0042] The coal consumption and air supply for each operating condition were determined based on boiler load, excess air coefficient, and burnout air rate. Boundary conditions were set for fuel inlet, primary air inlet, secondary air inlet, burnout air inlet, furnace outlet, and wall surface. Computational fluid dynamics software was used to numerically solve boiler combustion, establishing mass conservation, momentum conservation, energy conservation, and component transport equations. Three-dimensional steady-state calculations were performed using a pressure solver, with the SIMPLE algorithm used for pressure-velocity coupling. The standard k-ε model was used for gas-phase turbulence, the P-1 model for radiative heat transfer, and the finite-rate / eddy dissipation model for gas-phase chemical reactions. A stochastic trajectory model was used for pulverized coal particle tracking, a two-step competing rate model was used for volatile matter analysis, and a dynamic / diffusion control model was used for coke combustion. x The generation simulation uses a post-processing method, in which thermal NO... x and fuel-type NO x The extended Zeldovich mechanism and the de-Soete model are used for description, respectively. The same physical model, numerical discretization scheme, and control parameters are used to solve for each operating condition.

[0043] During the calculation, the residuals of the continuity equation, momentum equation, energy equation, and component transport equation are monitored, along with the furnace outlet temperature, outlet component concentration, and mass and energy balance. Convergence is considered achieved when the residuals of each equation reach the preset convergence criteria and the main monitored quantities remain stable throughout the iteration process.

[0044] For each operating condition, eight physical quantities are derived: velocity components in three directions, temperature, oxygen mass fraction, carbon dioxide mass fraction, carbon monoxide mass fraction, and nitrogen-containing pollutant mass fraction. These are integrated into a multi-condition multi-physics field snapshot database. For any physical quantity... The snapshot matrix is ​​represented as:

[0045] In the formula, physical quantity Snapshot matrix, For the first Physical quantities under training conditions In all The physical quantity column vector on each grid cell, where ,for example This represents 21 training scenarios. Snapshot matrices are established for each of the eight physical quantities to avoid scale bias caused by direct combination of different units.

[0046] Step S3: Based on the multi-physics snapshot database of the training conditions, perform global intrinsic orthogonal decomposition on each physical quantity, determine the number of modes to be retained according to the preset cumulative energy contribution rate, and establish a global order reduction model. Calculate the training mean field for the snapshot matrix of each physical quantity and then center it:

[0047]

[0048] In the formula, For a centralized snapshot matrix, This is the original physical field column vector of the simulation calculation results. To train the average field, For length is A column vector of all 1s. This indicates transpose.

[0049] Economic singular value decomposition of the centralized snapshot matrix:

[0050] In the formula, It is a left singular vector matrix, whose column vectors correspond to the spatial POD modes of each order; Let be a diagonal matrix composed of singular values, where Indicates the first Energy contribution of the first-order POD mode; Let be a right singular vector matrix, representing the evolution characteristics of each POD mode over time; This indicates transpose.

[0051] The cumulative energy contribution rate is calculated based on the squared singular values. The minimum number of modes that satisfies the following formula is selected. :

[0052] In the formula, To retain the number of modes for candidates, physical quantity The final number of POD modes retained. The cumulative energy threshold is set to 0.999 in this embodiment; Centralized snapshot matrix rank, physical quantity The There are 10 singular values. The mode matrix is ​​preserved as follows: The modal coefficients and reconstructed field for any given operating condition are expressed as:

[0053]

[0054] In the formula, For any given working condition, the vector parameters are... For the physical quantity under this working condition In all Column vectors of physical quantities on each grid cell For POD modal coefficients, Reconstructing the physics field for POD This indicates transpose. Independent test cases are projected only onto the modal space obtained from the 21 training cases and do not participate in the mean field and modal basis calculations. Step S4: Determine at least one primary partitioning direction based on the geometric structure and physical field distribution characteristics of the target complex industrial equipment, and divide the computational domain into multiple primary core subdomains according to the gradient change of the selected representative physical quantity along the direction. First, based on the boiler's geometry, boundary conditions, and major physical processes such as flow, combustion, and heat transfer, the furnace height direction is determined as the primary partition direction to reflect the main variation patterns of the multi-physical fields within the furnace along the height direction. Then, temperature, closely related to the aforementioned physical processes, is selected as a representative physical quantity, and the boundaries of the primary subdomains are determined based on the spatial gradient variation of temperature along the furnace height direction.

[0055] A height-divided grid was established along the furnace height, ranging from 0 m to 55.25 m, with a fixed width of 0.1 m. For each training condition, the volume-weighted average temperature was calculated within each height-divided grid based on the grid cell volume. The training scenario is in the first... The average temperature inside each height compartment is:

[0056] In the formula, For the first The training scenario is in the first... Average temperature within each height compartment For the first A collection of grid cells within a height-divided bin. For the first Volume of each grid cell For the first The training scenario is in the first... Temperature on each grid cell.

[0057] Calculate the absolute value of the height gradient of the temperature profiles under each of the 21 operating conditions, and take the maximum value from the training conditions at each height position as the aggregate gradient:

[0058] In the formula, For height position The polymerization temperature gradient at that location, For the first The volume-weighted average temperature profile for each training condition is obtained. The aggregated gradient is smoothed using a Gaussian filter; in this embodiment, the window width is 5. The 92nd percentile of the smoothed gradient is used as the high gradient threshold, and candidate boundaries are extracted based on the location of high gradient state changes. Intervals shorter than 8% of the total furnace height are merged with adjacent intervals, forming three primary core subdomains. Figure 3 The first primary subdomain is the cold ash hopper area, with a height range of approximately 0. 11.99 m; the second primary subdomain is the main combustion zone, with a height range of approximately 11.99 m. 26.28 m; the third primary subdomain is the heat exchange and flow development zone, with a height range of approximately 26.28 m. 55.25 m.

[0059] Step S5: Set up a first-level overlapping region between adjacent first-level core subdomains, and construct complementary continuous weights based on the normalized spatial position of grid points or grid cells in the first-level overlapping region. A symmetrical overlapping region is set at the interface of adjacent first-level core subdomains. The total width of the overlapping region is 5% of the shorter length of the two adjacent core subdomains, and is evenly distributed on both sides of the interface. Let the overlapping interval be [ , ], grid position The corresponding normalized position is:

[0060] In the formula, The left boundary of the overlapping interval. This represents the right boundary of the overlapping interval. The complementary continuous weights of the left and right first-level subdomains are as follows:

[0061]

[0062] In the formula, The normalized coordinates of the grid positions within the overlapping interval. The complementary continuous weights of the left-hand first-order subdomain, These are the complementary continuous weights of the right-hand first-order subdomain. Therefore, they satisfy... The weight of the first-level subdomain to which the non-overlapping core region belongs is 1. Based on this, a first-level core grid mapping, a first-level extended unit, and a first-level cosine splicing weight table are established.

[0063] Step S6: Construct a first-level local training snapshot matrix for each first-level subdomain and each physical quantity, perform intrinsic orthogonal decomposition, determine the number of retained modes for each first-level subdomain according to the preset cumulative energy contribution rate, and use the complementary continuous weights to fuse the local physical fields of each first-level subdomain to establish an intrinsic orthogonal decomposition model for the first-level subdomain. Local snapshot matrices are established for eight physical quantities in three first-order subdomains. Each local model employs training mean field centering and singular value decomposition, and the number of retained modes is determined based on a cumulative energy threshold of 0.999. The local eigenorthogonal decomposition method is consistent with that of the global model.

[0064] For independent test cases, orthogonal projection is performed using the trained local modes. During full field recovery, the local projected fields of each first-level subdomain are first calculated, and then the product of the first-level cosine weights and the local projected values ​​is accumulated to the global domain based on the global grid index. If there are uncovered grids, or if the first-level weights corresponding to any global grid deviate from 1 by more than the preset tolerance, the stitching is deemed invalid.

[0065] Step S7: Within each primary core subdomain, a secondary partitioning is performed using a multiphysics data-driven method. Joint clustering features are constructed using spatial coordinates and the statistical characteristics of physical quantities under training conditions. Physical quantities with large magnitude differences are transformed and standardized. A clustering method with spatial neighborhood constraints is used to partition the secondary subdomains. Within each primary core subdomain, multiphysics data obtained from high-fidelity training conditions are used as the basis for secondary partitioning, with grid cells as the clustering objects. Joint clustering features are constructed using the three-dimensional spatial coordinates of the grid cells, as well as the mean and standard deviation of each physical quantity under the internal fitting conditions. The mean reflects the average physical state of the grid cells, and the standard deviation reflects the sensitivity of the grid cells to changes in operating conditions.

[0066] The mass fractions of carbon monoxide and nitrogenous pollutants can vary significantly in magnitude across different spatial regions, and their distributions exhibit a marked skew. To reduce the excessive influence of local extreme values ​​on clustering distances and enhance the ability to identify differences between low-concentration regions, a logarithmic transformation is performed on the aforementioned physical quantities before constructing the clustering features:

[0067] In the formula, For quality fraction, It is a very small positive number, in this embodiment . No. On the _ grid cell _ The mean and standard deviation of each physical quantity are as follows:

[0068]

[0069] In the formula, F This is the set of internal fitting conditions; ; Number the physical quantities; Number the grid cells; physical quantity In the Under the internal fitting training condition, the first The values ​​taken on each grid cell; physical quantity In the The mean value over each grid cell; physical quantity In the The standard deviation over each grid cell.

[0070] No. The clustering feature vector of each grid cell is represented as:

[0071] In the formula, coordinate weights , , In this embodiment, the values ​​are 1, 0.5, and 1, respectively. and The first The physical quantity in the first... Mean and standard deviation of the internal fitting condition at each grid cell; Number the physical quantities. ; The total number of physical quantities in this implementation case. ; This indicates transpose. Since the first-level subdomains have already been partitioned along the height direction, the second-level clustering reduces the weight of the height coordinates to enhance the study of the lateral space, but spatial adjacency relationships still use the original three-dimensional coordinates.

[0072] After standardizing each feature dimension, MiniBatch K-means++ is used to obtain initial labels, and spatial regularization is performed based on the 8 spatial nearest neighbors of each grid cell. Candidate labels. The evaluation cost is expressed as:

[0073] In the formula, Number the global grid cells. Candidate tags The corresponding evaluation cost, For the first Clustering feature vectors of grid cells, For the first Cluster centers, For the feature dimension, The label of the spatial nearest neighbor of grid cell j is c proportion, For spatial regularization coefficients, in this embodiment... Set the value to 0.35. Spatial regularization iterates a maximum of 8 times, when the label change rate is no higher than... Stop when the time comes.

[0074] The minimum allowable number of grid cells for the second-level core subdomain takes into account the absolute size, the relative size of the first-level subdomain, and the ideal average cluster size:

[0075] In the formula, For the first The number of candidate second-level subdomains for each first-level subdomain is The minimum number of subdomain grid cells allowed at that time. This represents the number of core grids in the current first-level subdomain. This represents the number of candidate second-level subdomains. The absolute minimum number of grid cells; and This is a proportionality coefficient. In this embodiment, it can be taken as... , and After clustering, overly small subdomains are repaired by region growing, and connected components are identified using a grid nearest neighbor graph. Non-primary connected components are redistributed according to boundary contact degree and physical feature distance, so that the final secondary subdomains meet the spatial connectivity requirements.

[0076] Step S8: Set multiple candidate second-level subdomains for each first-level core subdomain, and adaptively determine the number of second-level subdomains for each first-level core subdomain based on the local order reduction error and partitioning results under internal verification conditions. Each first-level subdomain determines the number of candidate second-level subdomains. Four out of the 21 training scenarios (2 to 15) were selected as internal validation scenarios, and the remaining 17 were selected as internal fitting scenarios; the four official independent test scenarios did not participate in any testing. Value selection. All candidates. All cases undergo internal local POD verification, and only infeasible situations such as empty clusters, insufficient number of valid categories, or subdomain size below the preset minimum value are considered.

[0077] For candidates and the For each physical quantity, the sum of squared errors from all internal verification conditions is combined to calculate the relative L2 error, and then the eight physical quantities are averaged with equal weights to obtain the result. Select the principal error:

[0078] In the formula, For the first Each first-level subdomain is divided into The overall internal validation error when considering two second-level subdomains This is a set of internal verification operating conditions. The total number of physical quantities. Numbering for first-level subdomains. For the first Physical quantities under training conditions In the Vectors obtained from simulation within each first-order subdomain The number of candidate second-level subdomains is Time Local POD reconstruction results for a verification working condition. To prevent extremely small positive numbers with a denominator of zero.

[0079] Among the candidates that have completed internal POD validation, only those that meet the following criteria are retained: minimum subdomain size, worst-case maximum connected component ratio not less than 0.95, and fragmented mesh ratio not greater than 0.02. Feasible candidates are then categorized as follows: Arrange the data in ascending order and construct the monotonic historical optimal envelope of the verification error:

[0080] In the formula, This represents the permutation number of the currently feasible candidates. Not greater than Candidate permutation number, For the first The first-level subdomain of the Number of feasible candidate second-level subdomains For the first The first-level subdomain of the Number of historically feasible candidate second-level subdomains For the first Each first-level subdomain uses the number of candidates The overall internal verification error at that time As of the current candidate The historical minimum internal validation error was obtained. From the minimum feasible candidate... To the largest feasible candidate The total relative error improvement rate obtained across the entire candidate range is:

[0081] In the formula, For the first The overall relative improvement rate from the smallest feasible candidate to the candidate with the lowest historical error in each first-level subdomain. To prevent extremely small positive numbers with a denominator of zero. When When choosing the least feasible option, select the least feasible option. ;when At that time, The normalized coordinates of the sum of the error envelope are and :

[0082]

[0083] In the formula At the minimum feasible value Corresponding value to the historical lowest error The normalization results between them For the corresponding error envelope The normalized result between the first-end error and the minimum error It is a very small positive number.

[0084] Calculate the vertical deviation of each candidate point in the normalized error curve relative to the first and last chords, and use this deviation as the elbow score:

[0085] In the formula, For the first The first-level subdomain of the The error elbow score corresponding to each feasible candidate. When the maximum elbow score is not less than 0.05, the candidate corresponding to the maximum elbow score is selected. When the elbow is unclear, use the following formula to select the percentage of error improvement achieved for the first time. Minimum candidate :

[0086] In the formula, For the first The optimal number of second-level subdomains for each first-level parent domain; To achieve the target improvement ratio, this embodiment takes... =0.5. According to the above rules, the optimal number of the three first-level subdomains in this embodiment is... A total of 13 second-level subdomains are formed, referring to Figure 4 The first-level subdomain is divided into two second-level subdomains to satisfy the minimum subdomain size. The second primary subdomain was divided into 5 secondary subdomains, and the verification error curve obtained the maximum elbow score of 0.2294; the third primary subdomain was divided into 6 secondary subdomains, and the verification error curve obtained the maximum elbow score of 0.3331, both of which are higher than the elbow judgment threshold of 0.05.

[0087] Step S9: Construct a second-level local training snapshot matrix for each second-level subdomain and each physical quantity, calculate the local training mean field, perform eigenorthogonal decomposition on the snapshot matrix after removing the local training mean field, and determine the number of modes to be retained based on the preset cumulative energy contribution rate. To analyze the impact of the number of second-level subdomains on the reduced-order model, while maintaining consistency in training conditions, clustering features, standardization methods, spatial regularization coefficients, number of nearest neighbors, minimum region size constraints, connectivity restoration methods, POD inner product forms, and cumulative energy thresholds, the three first-level subdomains were uniformly fixed. , and .

[0088] Thus, the second-level subdomain part forms four schemes: adaptive There are a total of 13 second-level subdomains; fixed There are a total of 12 second-level subdomains; fixed There are a total of 15 second-level subdomains; fixed There are a total of 18 second-level subdomains. (Fixed) The proposed scheme serves as a control to evaluate whether the adaptive method achieves a reasonable balance between projection accuracy and model complexity under unknown operating conditions.

[0089] Step S10: The second-level subdomains are assigned non-overlapping weights. The weights of the second-level subdomains are multiplied by the complementary continuous weights of the corresponding first-level subdomains to form hierarchical weights. The hierarchical weights are then used to accumulate the local physical fields of each second-level subdomain to recover the complete multiphysics field and establish the intrinsic orthogonal decomposition model of the second-level subdomains. Local POD models were established for eight physical quantities in each second-level subdomain of the four second-level schemes. The training mean field, singular value decomposition, cumulative energy threshold of 0.999, and orthogonal projection method for independent test conditions were consistent with the global model. The first-level cosine weights and second-level attribution weights were not involved in mode extraction.

[0090] No second-level overlapping region is set between second-level subdomains. For global grid cells... First-level subdomain and its internal second-level subdomains The second-level subdomain affiliation indicator function is defined as follows:

[0091] In the formula, For the first Within the first-level subdomain The set of grid cells contained in each second-level subdomain. For the first The attribution indicator function for each grid cell to a second-level subdomain. The hierarchical weights are formed by multiplying the second-level attribution indicator value by the complementary continuous weight of the corresponding first-level subdomain:

[0092] In the formula, For hierarchical weights, For the first Each first-level subdomain in the grid cell complementary continuous weights at the point, For the attribution indicator function of the second-level subdomain, when the first The grid cell belongs to the first Within the first-level subdomain When there are two second-level subdomains, Otherwise, it is 0. The complete physical field is reconstructed as follows:

[0093] In the formula, For any working condition The next physical quantity in the first Reconstructed values ​​at each grid cell, For the second-level subfield in the th Reconstructed values ​​at each grid cell, This represents the total number of first-level subdomains. This represents the number of second-level subdomains within the first-level subdomain p.

[0094]

[0095] Since the second-level subdomains within the same first-level subdomain are hard partitioned, and the cosine weights of adjacent first-level subdomains are complementary, the sum of the hierarchical weights of any grid cell in the entire domain is 1.

[0096] Step S11: The intrinsic orthogonal decomposition models of the whole domain, first-level subdomain, and second-level subdomain are compared using unified evaluation conditions, and the target reduced-order model for prediction is determined by combining model error and scale. This embodiment forms three types of order reduction methods: global POD, first-level subdomain POD, and second-level subdomain POD; among which second-level subdomain POD includes adaptive methods. ,fixed ,fixed and fixed Four specific schemes. For the above six models, the reconstruction error for 21 training scenarios and the projection error for 4 independent test scenarios are calculated respectively, referring to... Figures 5a-5d As shown, Figure 5a5b is the original simulation calculation cloud map, 5c is the global POD reconstruction cloud map, 5d is the first-level subdomain POD reconstruction cloud map, and 5d is the second-level adaptive subdomain POD reconstruction cloud map.

[0097] Evaluation metrics include relative L2 error, root mean square error, and mean absolute error. The number of retained modes, basis matrix size, decomposition time, and complete field reconstruction and stitching time for each physical quantity are also statistically analyzed. The mean, standard deviation, and maximum values ​​of temperature, velocity modulus, carbon dioxide, carbon monoxide, and nitrogen-containing pollutants are calculated in the furnace outlet grid region.

[0098] For a physical field containing M grid cells, the relative L2 error, root mean square error, and mean absolute error are respectively expressed as:

[0099]

[0100]

[0101] In the formula, To simulate physical field vectors, To reconstruct the physical field vectors, The total number of grid cells participating in the evaluation; For the first Simulation calculation values ​​at each grid cell; For the first Reduced-order model reconstruction value at each grid cell; It is a very small positive number.

[0102] Step S12: Using the operating parameters of the target complex industrial equipment as input and the modal coefficients of the target reduced-order model as output, establish multiple candidate surrogate models respectively, and use cross-validation method to evaluate the predictive performance of each candidate surrogate model to determine the mapping relationship between operating parameters and modal coefficients. With running parameters As input, and using the locally preserved modal coefficients in the selected reduced-order kernel as output, a mapping from running parameters to local modal coefficients is established:

[0103]

[0104] In the formula, Numbering for first-level subdomains. For the first Numbering of second-level subdomains within a first-level subdomain. Number the physical quantities. For modal order, The number of modes to retain for the corresponding local physical quantity. For the corresponding proxy model, For the predicted individual mode coefficients, The prediction vector consists of all retained modal coefficients. Candidate surrogate models include Gaussian process regression, radial basis function, support vector regression, random forest, and multinomial response surface.

[0105] Leave-one-out cross-validation is used to evaluate the modal coefficient prediction error of candidate surrogate models. Furthermore, the predicted modal coefficients are substituted into local order reduction models and hierarchical splicing models to evaluate the complete physical field recovery error. Based on the cross-validation results, surrogate models can be selected for different subdomains, different physical quantities, or different modal orders.

[0106] Substituting the modal coefficients predicted by the surrogate model into the local reduced-order model, the local physics field is reconstructed as follows:

[0107] In the formula, To locally reconstruct the physical field, For local training average field, To preserve the POD mode matrix locally, This is the vector of local modal coefficients predicted by the surrogate model.

[0108] Step S13: Input the operating parameters to be predicted into the surrogate model to obtain the corresponding modal coefficients, and reconstruct the complete multiphysics field by combining the mean field and the reduced-order basis to build a multiphysics digital twin simulation model for predicting the operating status of complex industrial equipment.

[0109] The results of first- and second-level subdomain partitioning, local mean fields, local POD modes, surrogate model parameters, and hierarchical splicing relationships are uniformly stored to construct a multiphysics digital twin rapid simulation model. Boiler operating parameters for the predicted operating conditions are input into the surrogate model to obtain the modal coefficients corresponding to each local physical field. Based on the modal coefficients, local mean fields, and local POD modes, the physical fields of each subdomain are reconstructed, and the complete temperature field, velocity field, and component fields are obtained through hierarchical splicing, enabling the prediction of boiler multiphysics fields under different operating conditions.

[0110] Example 2: This embodiment is basically the same as Embodiment 1, focusing on the combustion chamber of a gas turbine. Operating parameters include fuel flow rate, air flow rate, and inlet temperature. Representative physical quantities are selected from temperature or OH radical concentration. The first-level zoning direction is determined based on the axial structure of the flame tube within the combustion chamber. The coordinate weights in the second-level clustering features are adjusted according to the geometric characteristics of the combustion chamber. The remaining steps are the same as in Embodiment 1, enabling efficient and accurate prediction of multiphysics fields.

[0111] Example 3: This embodiment is basically the same as Embodiment 1, focusing on the thermal management system of a lithium-ion battery. Operating parameters include charge / discharge rate, ambient temperature, and coolant flow rate. Multiple physical fields are used, including temperature, potential, and current density fields. The primary zoning direction is determined based on the electrode stacking direction within the battery cell. Temperature is selected as the representative physical quantity. The remaining steps are the same as in Embodiment 1, enabling rapid thermal field prediction for the battery thermal management system.

[0112] The various embodiments in this specification are described in a progressive manner, with each embodiment focusing on its differences from other embodiments. Similar or identical parts between embodiments can be referred to interchangeably. For the apparatus disclosed in the embodiments, since they correspond to the methods disclosed in the embodiments, the description is relatively simple; relevant parts can be referred to the method section.

[0113] The above description of the disclosed embodiments enables those skilled in the art to make or use the invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the invention. Therefore, the invention is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A multi-physics coupled digital twin simulation and prediction method for industrial equipment, characterized in that, include: Step S1: Determine the operating parameters and their value range of the target complex industrial equipment, generate multi-condition samples using experimental design methods, and divide the multi-condition samples into training conditions and independent test conditions. Step S2: Based on the geometric structure, material properties, boundary conditions and operating mechanism of the target complex industrial equipment, establish a high-fidelity multiphysics coupling model, perform numerical solutions for each sample working condition, and obtain a multiphysics snapshot database. Step S3: Based on the multi-physics snapshot database of the training conditions, perform global intrinsic orthogonal decomposition on each physical quantity, determine the number of modes to be retained according to the preset cumulative energy contribution rate, and establish a global order reduction model. Step S4: Determine at least one primary partitioning direction based on the geometric structure and physical field distribution characteristics of the target complex industrial equipment, and divide the computational domain into multiple primary core subdomains according to the gradient change of the selected representative physical quantity along the direction. Step S5: Set up a first-level overlapping region between adjacent first-level core subdomains, and construct complementary continuous weights based on the normalized spatial position of grid points or grid cells in the first-level overlapping region. Step S6: Construct a first-level local training snapshot matrix for each first-level subdomain and each physical quantity, perform intrinsic orthogonal decomposition, determine the number of retained modes for each first-level subdomain according to the preset cumulative energy contribution rate, and use the complementary continuous weights to fuse the local physical fields of each first-level subdomain to establish an intrinsic orthogonal decomposition model for the first-level subdomain. Step S7: Within each primary core subdomain, a secondary partitioning is performed using a multiphysics data-driven method. Joint clustering features are constructed using spatial coordinates and the statistical characteristics of physical quantities under training conditions. Physical quantities with large magnitude differences are transformed and standardized. A clustering method with spatial neighborhood constraints is used to partition the secondary subdomains. Step S8: Set multiple candidate second-level subdomains for each first-level core subdomain, and adaptively determine the number of second-level subdomains for each first-level core subdomain based on the local order reduction error and partitioning results under internal verification conditions. Step S9: Construct a second-level local training snapshot matrix for each second-level subdomain and each physical quantity, calculate the local training mean field, perform eigenorthogonal decomposition on the snapshot matrix after removing the local training mean field, and determine the number of modes to be retained based on the preset cumulative energy contribution rate. Step S10: The second-level subdomains adopt non-overlapping assignment, and the assignment weight of the second-level subdomains is multiplied with the complementary continuous weight of the corresponding first-level subdomains to form the hierarchical weight. The second-level local physical fields are accumulated according to the hierarchical weights to recover the complete multiphysics field and establish a second-level subdomain eigenorthogonal decomposition model. Step S11: The intrinsic orthogonal decomposition models of the whole domain, first-level subdomain and second-level subdomain are compared using unified evaluation conditions, and the target reduced-order model for prediction is determined by combining model error and scale. Step S12: Using the operating parameters of the target complex industrial equipment as input and the modal coefficients of the target reduced-order model as output, establish multiple candidate surrogate models respectively, and use cross-validation method to evaluate the predictive performance of each candidate surrogate model to determine the mapping relationship between operating parameters and modal coefficients. Step S13: Input the operating parameters to be predicted into the surrogate model to obtain the corresponding modal coefficients, and reconstruct the complete multiphysics field by combining the mean field and the reduced-order basis to build a multiphysics digital twin simulation model for predicting the operating status of complex industrial equipment.

2. The method according to claim 1, characterized in that, In step S4, the first-level partition direction is determined based on the geometric structure and main physical processes of the target complex industrial equipment, and the representative physical quantity is a physical quantity closely related to the main physical process; the gradient change is determined by calculating the spatial gradient of the representative physical quantity along the direction under each training condition and taking the aggregated gradient, and the aggregated gradient is the maximum value or statistical quantile of the absolute value of the gradient under each training condition.

3. The method according to claim 1, characterized in that, In step S5, the complementary continuous weights are constructed using a cosine function: let the overlapping interval be [ , The normalized position corresponding to grid position y is: In the formula, The left boundary of the overlapping interval. The right boundary of the overlapping interval; the complementary continuous weights of the left and right first-level subdomains are as follows: In the formula, The normalized coordinates of the grid positions within the overlapping interval. The complementary continuous weights of the left-hand first-order subdomain, Let be the complementary continuous weights of the right-hand first-order subdomain; where satisfies .

4. The method according to claim 1, characterized in that, In step S7, the joint clustering features include the three-dimensional spatial coordinates of the grid cells and the mean and standard deviation of each physical quantity under the internal fitting conditions; for physical quantities with skewed distributions, a logarithmic transformation is performed on them before constructing the clustering features; the spatial neighborhood constraint is achieved by calculating the proportion of each candidate cluster in the spatial nearest neighbors of each grid cell and introducing a regularization term to correct the clustering cost function.

5. The method according to claim 4, characterized in that, The cost function is: In the formula, Number the global grid cells. Candidate tags The corresponding evaluation cost, For the first Clustering feature vectors of individual grid cells; For the first Cluster centers, For the feature dimension, For grid cells Spatial neighbors are tagged as c proportion, is the space regularity coefficient.

6. The method according to claim 1, characterized in that, In step S8, the adaptive determination of the number of second-level subdomains includes: For each primary core subdomain, the number K of candidate secondary subdomains is set to a value within a preset range; The training scenarios are divided into internal fitting scenarios and internal validation scenarios. For each candidate K, a local intrinsic orthogonal decomposition is performed and the order reduction error of the internal validation scenario is calculated. Based on the curve of the reduction error as a function of K, the optimal K value is automatically selected using the elbow rule or the improvement rate threshold method.

7. The method according to claim 6, characterized in that, The improvement rate threshold method includes: calculating the total improvement rate R from the smallest candidate feasible K to the largest candidate feasible K; selecting the smallest feasible K when R is less than the first threshold; otherwise, calculating the elbow score of the normalized error curve; selecting the K with the largest elbow score when the maximum elbow score is greater than or equal to the second threshold; otherwise, selecting the K that first reaches the preset error improvement ratio.

8. The method according to claim 1, characterized in that, In step S10, the hierarchical weights are: in, For hierarchical weights, For the first Each first-level subdomain in the grid cell complementary continuous weights at the point, This is the attribution indicator function for the second-level subdomain, and it satisfies the condition for all grid cells. have , This represents the total number of first-level subdomains. This represents the number of second-level subdomains within the first-level subdomain p.

9. The method according to claim 1, characterized in that, In step S12, the candidate surrogate model includes one or more of Gaussian process regression, radial basis function, support vector regression, random forest and polynomial response surface; different surrogate models are selected for different subdomains, different physical quantities or different modal orders.

10. The method according to claim 1, characterized in that, In step S2, the multiphysics field includes one or more of the following: velocity field, temperature field, component concentration field, and pressure field.

Images