Multiphysics Simulation Method for Rock Wool Integrated Panel Composite Production Line Based on Digital Twin

Abstract

This invention relates to the field of digital twin technology, specifically a multiphysics simulation method for a rock wool integrated panel composite production line based on digital twins. The method includes: first, constructing a representative volume element and parameterized model reflecting the microstructure of rock wool, and then building a thermal-fluid coupling parameter response library using the lattice Boltzmann method. Subsequently, full-order simulation is performed using the SST turbulence model, the medium radiation model, and the component migration equation, combined with dynamic boundary compensation to map panel displacement and obtain the full-order distribution field of the production line. Further, spatial basis functions are extracted using intrinsic orthogonal decomposition, and a low-dimensional feature space is constructed using discrete empirical interpolation and Galerkin projection. Finally, using a parameterized background data weak observer and Gaussian process regression, the measured data is projected onto the low-dimensional space for coefficient reconstruction and fitting, achieving self-correcting twin field reconstruction. This invention solves the problem of balancing simulation accuracy and real-time performance under complex production line conditions, improving the predictive reliability of the simulation.

Description

Technical Field

[0001] This invention relates to the field of digital twin technology, specifically to a multiphysics simulation method for a rock wool integrated panel composite production line based on digital twins. Background Technology

[0002] In the production process of integrated rock wool panels, traditional numerical simulation methods are mostly based on static physical parameters for coarse calculations, which are difficult to accurately describe the dynamic evolution of the internal microporous structure of the panels under the high-temperature wind field of the curing furnace. Existing business scheduling systems and simulation platforms lack the ability to characterize the multi-field coupling process of heat, flow, mass, and radiation in porous media in real time, resulting in significant deviations between simulation results and actual production line conditions, making it difficult to meet the real-time requirements of sub-second self-correction and refined control in high-speed production. Traditional simulation methods often leave key physical fields such as heat transfer, component migration, and media radiation in a "logical island" state during numerical calculations, making it difficult to achieve high-fidelity collaborative solutions for multiple physical fields. Moreover, the simulation architecture based on full-order full-mesh has an excessively large degree of freedom. When dealing with iterative calculations involving turbulent source terms and radiation nonlinear residuals, the computational complexity increases exponentially, and the solution time far exceeds the sub-second response requirements of the controlled objects on the production line, making real-time numerical extrapolation infeasible in engineering practice.

[0003] To address the challenge of performing dynamic simulations and self-correction techniques with high physical reliability and real-time requirements on the spatial morphology evolution and internal state of rock wool throughout its production lifecycle in complex multi-physics coupled environments, a multi-physics simulation method for integrated rock wool panel composite production lines based on digital twins is proposed. Summary of the Invention

[0004] This invention aims to provide a multi-physics simulation method for rock wool integrated board composite production lines based on digital twins. It constructs a multi-scale simulation framework that deeply integrates microscopic property characterization and macroscopic order reduction model, so as to realize real-time accurate prediction and virtual-real dynamic correction of the multi-physics state of rock wool production lines.

[0005] To achieve the above objectives, the present invention provides the following technical solution: A multi-physics simulation method for integrated rock wool panel composite production lines based on digital twins includes: Using a random fiber filling algorithm and fractal geometry theory, a representative volume element reflecting the microstructure of rock wool fibers is constructed, and a parameterized structural model based on the porosity evolution law is established. The heat flow coupling field in the micropores is simulated by combining the lattice Boltzmann method and the volume is averaged to construct a parameter response library. The SST turbulence model is used to simulate the high-speed circulating air jet in the curing oven. The parameter response library is retrieved as the variable parameter source term of the porous medium region. It is integrated to participate in the medium radiation model and the migration of porous medium components. Combined with dynamic boundary compensation, the displacement of the board on the production line is mapped to obtain the full-order distribution field of the production line. The spatial basis functions of the full-order distribution field of the production line are extracted by eigenorthogonal decomposition. The control equations are then projected onto the feature subspace by discrete empirical interpolation and Galerkin projection to obtain a low-dimensional feature space. By using a parameterized background data weak observer and Gaussian process regression, the sensor data is projected onto the low-dimensional feature space for coefficient reconstruction and least squares fitting to obtain a self-correcting twin field.

[0006] Preferably, the step of constructing a representative volume element reflecting the microstructure of rock wool fibers using a random fiber filling algorithm and fractal geometry theory, and establishing a parameterized structural model based on the porosity evolution law, includes the following steps: Based on the average diameter and aspect ratio distribution of rock wool fibers, an initial set of fiber central axes is generated in a set of three-dimensional Euclidean spaces. The anisotropic distribution of the central axis set is controlled by a fiber orientation probability density function set based on the frequency of the production line pendulum. Combined with fractal geometry theory, multi-level bifurcation and surface roughness iteration are performed on individual fibers to generate fiber geometric entities with interpenetrating structures, constituting the initial representative volume element. Using the real-time pressure distribution of the production line pressure rollers as the driving variable, a geometric compression transformation operator describing the spatial deformation of the fiber skeleton is established. The initial representative volume element is then spatially scaled and mapped using the geometric compression transformation operator, dynamically adjusting the number of contact points and the porosity between fibers to obtain a parameterized structural model that evolves in real-time with the porosity under production line conditions.

[0007] Preferably, the step of constructing the parameter response library includes: The representative volume element is discretized into a three-dimensional grid space, and velocity distribution functions describing fluid motion and temperature distribution functions describing energy transfer are constructed respectively. Pressure gradients and temperature difference boundaries are applied to the grid space, and iterative calculations are performed using collision and migration operators until the velocity and temperature fields within the micropores converge. The local forces and heat fluxes at the fluid-solid interface are calculated using the momentum exchange method. The converged microscale field is spatially integrated and averaged. The anisotropic permeability tensor is obtained based on the linear mapping relationship between the average filtration velocity and the pressure gradient. The equivalent thermal conductivity tensor is obtained based on the mapping relationship between the heat flux density distribution and the macroscopic temperature gradient. The fiber-fluid heat transfer coefficient is obtained based on the ratio of the total heat transfer at the fluid-solid interface to the average temperature difference. By changing the porosity of the parameterized structural model, the above steps are repeated to establish a correspondence table with porosity as the independent variable and the anisotropic permeability tensor, the equivalent thermal conductivity tensor, and the fiber-fluid heat transfer coefficient as dependent variables, thus forming the parameter response library.

[0008] Preferably, the step of obtaining the full-order distribution field of the production line includes: A continuity equation and momentum equation based on the SST turbulence model are established to describe the wind speed field inside the curing oven. Simultaneously, the energy relationship between the radiation term and component migration term of the integrated participating medium is established to describe the heat and mass exchange process inside and on the surface of the board. The board region is defined as a porous medium region. According to the current porosity of the parameterized structural model, the corresponding values ​​in the parameter response library are retrieved in real time as the resistance source term in the momentum equation and the heat conduction and heat exchange correction terms in the energy equation. According to the actual conveying speed of the production line, the geometric position of the board is updated using a coordinate transformation operator to compensate the flow field information in the curing oven to the moving board surface boundary. Then, full-order numerical solutions are performed, and the output is aggregated to obtain the full-order distribution field of the production line, which includes the surface heat flux spectrum, the internal moisture content gradient, and the topological space wind pressure gradient line.

[0009] Preferably, the step of obtaining the low-dimensional feature space includes: The full-order distribution fields of the production line at multiple time steps are combined into a snapshot matrix. Singular value decomposition is then performed on the snapshot matrix using eigenorthogonal decomposition to extract spatial basis functions reflecting the principal components of the production line's physical evolution, constructing a low-dimensional linear subspace. Nonlinear terms in the control equations are identified, and sparse sampling points are selected in the physical space using discrete empirical interpolation. An approximate projection basis for the nonlinear operator is constructed based on the spatial basis functions. Using the spatial basis functions and the approximate projection basis of the nonlinear operator, the full-order momentum and energy equations are mapped to the low-dimensional linear subspace via Galerkin projection, eliminating the dependency of computational complexity on the full-order mesh degrees of freedom, thus obtaining the low-dimensional feature space.

[0010] Preferably, the step of obtaining the self-correcting twin field includes: Using a parameterized background data weak observer, discrete measurement point data collected by sensors are used as constraints, and variational assimilation calculations are performed in the low-dimensional feature space to extract initial modal coefficients that match the current measured conditions. A Gaussian process regression model is used to predict the nonlinear residuals between the low-dimensional feature space and the full-order distribution field, generating a calibration operator to compensate for unmodeled physical effects of the system. The initial modal coefficients and the calibration operator are linearly combined using a least-squares fitting operator to reconstruct the corrected dynamic modal coefficients. The corrected dynamic modal coefficients are projected back to the spatial basis functions to synthesize a self-correcting twin field with physical conservation properties and real-time bias correction.

[0011] Preferably, the self-correcting twin field is a real-time dynamically distributed field with physical conservation constraints, and its synthesis process includes: The corrected dynamic mode coefficients are projected back to the characteristic subspace defined by the spatial basis functions, and the set of grid node values ​​under the full-order computational domain is obtained through matrix linear reconstruction. The temperature field evolution process and moisture content gradient distribution inside the rock wool board are presented in real time on the monitoring terminal in the form of a three-dimensional voxel cloud map.

[0012] Compared with the prior art, the beneficial effects of the present invention are as follows: 1. Achieve precise spatiotemporal correction: By introducing a parameterized background data weak observer and a variational assimilation mechanism, a technological leap from static simulation to dynamic spatiotemporal calibration is achieved, effectively compensating for the unmodeled effects of the system and ensuring that the twin field boundary does not cause topological interference with the actual scheduling.

[0013] 2. Constructing a highly reliable multi-scale correlation: By utilizing the lattice Boltzmann method and fractal geometry theory, a deep alignment between the evolution of microscopic pore structures and macroscopic physical property parameters was achieved, accurately capturing the property decay laws under the influence of pressure and heat flow, and providing physically reliable evolution data for the refined management of three-dimensional storage yards.

[0014] 3. Breakthrough in real-time computing efficiency bottleneck: By using the POD-DEIM order reduction algorithm and Galerkin projection technology, the degree of freedom dimension of the control equations is greatly reduced, achieving sub-second response while preserving nonlinear characteristics, thus solving the technical limitation of excessive computational load for complex multiphysics models in high-speed production line environments. Attached Figure Description

[0015] Figure 1 The figure shows the multiphysics simulation method for the integrated rock wool board composite production line based on digital twins according to the present invention. Figure 2 A schematic diagram illustrating the process of constructing the parameter response library for this invention; Figure 3This is a schematic diagram of the process for solving the full-order distribution field and generating the self-correcting twin field in this invention. Detailed Implementation

[0016] 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.

[0017] Please see Figures 1 to 3 This invention provides a multiphysics simulation method for integrated rock wool panel composite production lines based on digital twins. The technical solution is as follows: Using a random fiber filling algorithm and fractal geometry theory, a representative volume element reflecting the microstructure of rock wool fibers is constructed, and a parameterized structural model based on the porosity evolution law is established. The heat flow coupling field in the micropores is simulated by combining the lattice Boltzmann method and the volume is averaged to construct a parameter response library. The SST turbulence model is used to simulate the high-speed circulating air jet in the curing oven. The parameter response library is retrieved as the variable parameter source term of the porous medium region. It is integrated to participate in the medium radiation model and the migration of porous medium components. Combined with dynamic boundary compensation, the displacement of the board on the production line is mapped to obtain the full-order distribution field of the production line. The spatial basis functions of the full-order distribution field of the production line are extracted by eigenorthogonal decomposition. The control equations are then projected onto the feature subspace by discrete empirical interpolation and Galerkin projection to obtain a low-dimensional feature space. By using a parameterized background data weak observer and Gaussian process regression, the sensor data is projected onto the low-dimensional feature space for coefficient reconstruction and least squares fitting to obtain a self-correcting twin field.

[0018] Example 1: This embodiment is set in the scenario of real-time quality monitoring and digital operation and maintenance of building insulation rock wool board production line. In response to the problem that the micro-topology of rock wool fibers is complex and dynamically changes with the pressure of the production line during the processes of blowing, pendulum laying and pressure roller collecting, the system performs the reconstruction and evolution of micro-representative volume elements through a structural simulation engine integrated into the edge computing node to support the subsequent online prediction of thermal conductivity.

[0019] Furthermore, by utilizing the random fiber filling algorithm and fractal geometry theory, a representative volume element reflecting the microstructure of rock wool fibers is constructed, and a parameterized structural model based on the porosity evolution law is established. Specifically, the system first acquires the statistical characteristics of rock wool fibers obtained offline by the production line vision sensor, sets the average diameter value to 5μm, the average length-to-diameter ratio distribution to 200, and initializes a discrete geometric grid coordinate system in a cubic three-dimensional Euclidean space with a side length of 1mm.

[0020] Specifically, the random fiber filling algorithm adopts a Monte Carlo sampling-based generation strategy. Its input receives the aforementioned diameter and aspect ratio parameters, and in memory, it determines the coordinates (x, y, z) of the starting point of each fiber and its direction vector by generating uniform random numbers in three-dimensional space.

[0021] Specifically, the algorithm logic generates a sequence of line segments along the direction vector by setting a linear iterative operator with a step size of 0.1 mm. When the cumulative length of a single fiber reaches the physical upper limit determined by the product of the diameter and the aspect ratio, growth stops, thereby outputting an initial set of fiber central axes containing 5000 discrete line segments in the virtual computing domain.

[0022] Summary of results: This step achieved a preliminary mapping from macroscopic statistical features to microscopic topological structures through random sampling, laying a high-precision skeletal foundation for the subsequent evolution of complex geometric forms.

[0023] Specifically, the system acquires the operating parameters of the pendulum on the production line in real time through the industrial Ethernet interface, extracts the pendulum frequency of 0.5Hz, and inputs it as an independent variable into the orientation probability distribution function based on the Gaussian mixture model. The distribution bias of the fiber central axis in the X-Y plane and the Z-axis direction is adjusted by calculating the weight coefficient.

[0024] Specifically, the bifurcation iteration algorithm in fractal geometry theory is used to perform topological expansion on the central axis. This algorithm adopts the fractal rules based on the Lindenmayer system. By setting the probability of second-order bifurcation to 0.15 and the probability of third-order bifurcation to 0.05, lateral branch nodes are automatically generated on each main central axis according to the preset offset angle, simulating the fine fiber branches generated by uneven stretching of rock wool in the high-temperature melting state.

[0025] Specifically, the system calls a surface reconstruction operator based on fractal Brownian motion. By adjusting the Hurst exponent component to 0.7, micron-level random perturbation displacements are superimposed on the surface of the fiber cylinder to generate a non-smooth geometric surface that can simulate the distribution of real slag. Finally, a fiber geometric entity with a complex interpenetrating topology is generated in the graphics processor memory through Boolean merging operations. The output is an initial representative volume element containing one million mesh vertices.

[0026] Specifically, the system collects physical load data of the pressure roller section of the production line in real time through a pressure transmitter, obtains a real-time pressure value of 50 kPa, and converts it into a compression drive variable represented in 32-bit floating-point format.

[0027] Specifically, the geometric compression transformation operator adopts non-uniform scaling logic based on the affine transformation matrix, which internally nests a deformation coefficient prediction model based on a multilayer perceptron (MLP) architecture. This model architecture includes an input layer that receives pressure values, two hidden layers each containing 32 neurons and using the ReLU activation function, and an output layer that outputs a three-dimensional spatial deformation tensor.

[0028] Specifically, the model is trained using the Adam optimizer on a sample set containing 3000 sets of stress-strain experimental data. Its algorithm processes data by mapping the 50 kPa pressure component to a contraction scaling factor of 0.85 in the Z-axis direction and X−Y. The lateral expansion coefficient of the plane is 1.05, thus generating a geometric compression transformation matrix that describes the spatial deformation of the fiber skeleton.

[0029] Specifically, the computing unit performs matrix multiplication, applies the generated geometric compression transformation matrix to the coordinates of all grid vertices of the initial representative volume element, and compacts the fiber geometry in the Z-axis direction by performing nonlinear reprojection of spatial coordinates.

[0030] Specifically, the algorithm logic uses a KD-tree-based spatial retrieval operator to re-search and count the physical overlap points between compressed fiber entities every 10ms, increasing the number of physical contact points between fibers from the initial 500 to 1200.

[0031] Specifically, the Monte Carlo ray tracer operator is used to calculate the porosity characteristics after compaction. By emitting 1,000 virtual probe rays into the volume element and calculating the curvature of the ray path through the void, the effective porosity value of the volume element is updated in real time, dynamically evolving from 95% to 88%.

[0032] Specifically, the system encapsulates the above-mentioned multidimensional features, including dynamic contact weights and transient porosity distribution, into a 32-bit binary sequence to obtain a parameterized structural model that evolves in real time with the porosity of the production line under different operating conditions.

[0033] The geometric feature space of the representative volume element includes: acquiring the real-time oscillation frequency trajectory of the production line pendulum; calculating the angular distribution weight of the fiber in the central axis set using the trajectory; and determining the initial spatial vector of each fiber using a probability sampling algorithm based on the angular distribution weight, so as to characterize the quasi-layered anisotropic distribution of the fiber formed by the reciprocating motion of the mechanical pendulum.

[0034] The real-time oscillation frequency of the pendulum on the production line is 1.2Hz. Combined with the rotation speed of the cotton-laying mesh belt at 0.5m / s, the fiber interlacing trajectory within a unit area is calculated. Based on this trajectory, the algorithm establishes an angular probability density distribution function in polar coordinates and performs spatial orientation mapping on the central axes of 5000 initial fibers to ensure that the permeability of the generated representative volume element in the main oscillation direction is higher than that in the vertical direction. By directly mapping the macroscopic pendulum mechanical motion parameters to the probability weights of microscopic fiber orientation, the problem of physical anisotropy error caused by the inability of traditional random filling models to reflect the real production process is solved, improving the prediction accuracy of the microscopic permeability tensor by about 15.4%.

[0035] By constructing a random fiber filling algorithm and a fractal geometric model, the problem of quantifying the dynamic evolution of microstructure in rock wool manufacturing was solved. By using a geometric compression transformation operator and a closed-loop mapping of real-time pressure data, sub-second updates of porosity and the number of contact points were achieved, ensuring the consistency of the digital twin system's representation under varying working conditions.

[0036] Reference Figure 2 The flowchart of the present invention for constructing a parameter response library is shown in the figure. Further, the thermal flux coupling field in the micropores is simulated by the lattice Boltzmann method and the volume is averaged to construct the parameter response library. Specifically, the system receives a cubic parametric structural model with a side length of 1 mm output by the preceding module, and discretizes it into a three-dimensional grid space with a resolution of 256×256×256 using voxel mesh generation technology. The grid spacing is set to 3.9 μm to ensure that the flow boundary layer characteristics of a single fiber surface can be captured.

[0037] Specifically, the system initializes the velocity distribution function describing fluid motion at each grid point, and uses the D3Q19 model with 19 discrete velocity vectors to characterize the air flow state inside the rock wool. At the same time, it constructs the temperature distribution function describing energy transfer, and uses the D3Q7 model with 7 discrete directions to characterize the heat migration behavior between the fiber and the fluid.

[0038] Specifically, the algorithm logic applies a driving pressure gradient of 50 Pa along the Z-axis in the grid space and sets the temperature deviation between the top and bottom surfaces to 50 K, thereby establishing the initial physical environment boundary for simulating heat flow penetration under pressure on the rock wool board.

[0039] Specifically, the system performs collision operator operations based on a single relaxation time model. It simulates the collision dissipation process between molecules by calculating the deviation between the distribution function in each direction and the local equilibrium distribution function. Then, it uses a migration operator to translate the post-collision distribution function along the corresponding 19 discrete velocity vectors to adjacent grid points.

[0040] Specifically, the algorithm performs more than 10,000 spatiotemporal iterations on the computing server, monitors the rate of change of kinetic energy and internal energy across the entire field in real time, and determines the model convergence when the relative error of the velocity field and the relative error of the temperature field are both below 0.000001 in three consecutive time steps.

[0041] Specifically, the momentum exchange method is used to process the distribution function at the fluid-solid interface. By calculating the momentum difference between the distribution function rebounding back to the fluid region and the distribution function flowing towards the solid surface, the frictional resistance on the fiber surface and the local heat flux of 5.2 W per square meter are accurately extracted.

[0042] Specifically, the system performs spatial integral averaging on the converged microscale velocity field and temperature field, and obtains the macroscopic average filtration velocity by accumulating the velocity vectors of 256×256×256 grid points in the entire computational domain and dividing by the total volume.

[0043] Specifically, the algorithm maps the average filtration velocity to the applied 50Pa pressure gradient using a tensor division based on Darcy's law, and calculates and outputs a 3×3 anisotropic permeability tensor that reflects the gas permeability inside the rock wool.

[0044] Specifically, the system calculates the volume average heat flux density within the spatial grid and divides it by the macroscopic temperature gradient generated by the 50K temperature difference, thereby obtaining the equivalent thermal conductivity tensor characterizing the thermal conductivity properties of the rock wool microstructure.

[0045] Specifically, the total heat transfer at the fluid-solid interface is statistically analyzed using the area integral operator. By calculating the ratio of this heat transfer to the difference between the average temperature of the fiber and the average temperature of the fluid, the fiber-fluid heat transfer coefficient, which characterizes the intensity of microscale convective heat transfer, is obtained.

[0046] Specifically, in order to realize the intelligent calling and prediction of the parameter library, a deep residual shrinking network is adopted as the specific implementation method of the parameter response library. A proxy model based on the deep residual shrinking network architecture is introduced. The model architecture includes one input layer for receiving porosity parameters, three consecutive residual blocks, and one output layer containing 15 neurons. The training details of the deep residual shrinking network are as follows: the loss function is the mean squared error; 1000 samples are divided into 800 training sets, 100 validation sets, and 100 test sets in a ratio of 8:1:1; the Adam optimizer is used during training, with an initial learning rate of 0.001 and a learning rate decay strategy of decreasing to 0.5 times every 500 rounds, and a weight decay coefficient (L2 regularization) of 0.0001; the input porosity is normalized to the [0, 1] interval, and the output physical property parameters are normalized to their respective [0, 1] intervals, and inverse normalization is performed during inference to restore physical units; the training termination condition is that the validation set loss does not decrease for 50 consecutive rounds, or the maximum number of iterations of 3000 rounds is reached.

[0047] Specifically, the residual block integrates a sub-network for automatically extracting feature importance. This sub-network contains two fully connected layers and uses the Sigmoid activation function to output weight coefficients. It is trained for 3000 rounds on 1000 sample sets obtained by iteratively calculating by changing the porosity from 80% to 98% using the Adam optimizer with a learning rate of 0.001.

[0048] Specifically, the model generates a set of continuous response sequences covering permeability components, thermal conductivity components, and heat transfer coefficients through the output layer, and establishes a correspondence with porosity as the independent variable in the form of a binary mapping table. Finally, the parameter response library is constructed in non-volatile memory.

[0049] The specific steps of the lattice Boltzmann method simulation include: setting the relaxation time τ to a value between 0.6 and 1.0, so that the calculated kinematic viscosity is similar to the actual kinematic viscosity of the air inside the rock wool. The matching rate is 1 square meter per second; a pressure gradient is implemented at the boundary of the grid space using the Zou-He boundary condition, with the top surface set as the pressure inlet and the bottom surface as the pressure outlet, and the pressure difference is 50 Pa; the temperature difference boundary uses an isothermal boundary condition, with the top surface temperature set at 323 K and the bottom surface temperature set at 273 K; the ratio of the time step Δt to the spatial step Δx is set so that the grid velocity satisfies the low Mach number approximation; the convergence criterion is that within three consecutive time steps, the L2 norm change rate of the entire velocity field is less than 0.000001, and the L2 norm change rate of the entire temperature field is less than 0.000001.

[0050] By performing in-depth deconstruction of the physical field state of rock wool microstructure using the lattice Boltzmann method, and transforming the solution of complex partial differential equations into parameterized response relationships using the momentum exchange method and volume averaging logic, the problem of difficulty in obtaining thermal indicators in real time under dynamic pressure conditions of production lines is solved, providing accurate physical mechanism support for the construction of high-fidelity digital twin systems.

[0051] The parameter response library uses a deep residual shrinkage network as a surrogate model. The input is a porosity scalar, and the output is a 15-dimensional vector containing 6 independent components of the permeability tensor, 6 independent components of the thermal conductivity tensor, and 3 components of the heat transfer coefficient. When calling the macroscopic simulation, the current porosity is input into the neural network for forward inference to directly obtain the corresponding physical property parameters without the need for discrete interpolation.

[0052] Furthermore, the SST turbulence model is used to simulate the high-speed circulating air jet in the curing oven. The parameter response library is retrieved as the variable parameter source term of the porous medium region, which is integrated into the medium radiation model and the migration of porous medium components. Combined with dynamic boundary compensation, the displacement of the board on the production line is mapped to obtain the full-order distribution field of the production line. Specifically, the system first acquires the real-time rotational speed signal fed back by the curing oven fan, extracts the initial boundary parameters with a value of 1500 revolutions per minute, and transforms them into a high-speed circulating air jet boundary with an air volume of 15 m / s. In the three-dimensional computational domain, the system establishes the continuity equation and momentum equation based on the SST turbulence model. Using the mixing function in this model, the system adopts the logic of using the k-omega model near the wall and the k-epsilon model in the mainstream region to accurately solve the complex vortex structure inside the curing oven.

[0053] Specifically, the algorithm logic synchronously establishes the energy relationship between the radiation term and the component migration term of the integrated participating medium. By introducing the discrete coordinate DO radiation model to process the thermal infrared radiation energy with a wavelength of 2.5μm, and combining Fick's diffusion law to establish a component migration operator to describe the evaporation of phenolic resin moisture in the board, the physicochemical evolution process of the board after being heated is described on a macroscopic scale.

[0054] Specifically, the system defines the space occupied by the board in the 10m curing oven calculation domain as a porous medium region. It extracts the current effective porosity value of 88% in real time from the previously generated parameterized structural model and uses this value as the index key to retrieve and match the corresponding 3×3 anisotropic permeability tensor, equivalent thermal conductivity tensor, and fiber-fluid heat transfer coefficient from the parameter response library.

[0055] Specifically, the algorithm logic performs a source term injection operation, inverts the retrieved permeability tensor and performs a tensor product operation with the velocity vector, which serves as the Darcy-Fochheimer drag source term in the momentum equation. At the same time, it uses the retrieved heat transfer coefficient to correct the fluid-solid heat exchange term in the energy equation, ensuring that the flow resistance and heat storage effect of the micro-fiber skeleton can be reflected in the macroscopic flow field solution.

[0056] Specifically, the system obtains real-time displacement data of the production line conveyor belt through the PLC communication interface, extracts the actual conveying speed of 5m per minute, and executes coordinate transformation operators within each simulation time step of 50ms.

[0057] Specifically, the operator uses Lagrange tracking logic to perform position translation mapping on the moving plate mesh, and compensates the flow field pressure vector and temperature gradient data output by the preset fixed jet nozzle in the curing oven to the moving boundary layer on the plate surface through spatial linear interpolation technology, thereby realizing the accurate restoration of the heating process of the plate during dynamic transportation.

[0058] Specifically, in order to achieve efficient aggregation and dimensionality reduction feature extraction of massive computational data, the system introduces a full-order field feature reconstruction model based on a convolutional neural network (CNN) architecture. This model architecture includes one input layer for receiving fluid node data, three consecutive convolutional layers configured with 32, 64, and 128 three-dimensional convolutional kernels of size 3×3×3 respectively, and using the ReLU activation function. After passing through a global average pooling layer containing 512 neurons, the Adam optimizer is used to output the reconstructed full-order distribution features after 1000 iterations.

[0059] Specifically, the algorithm logic performs pressure-velocity correction calculations of the coupled equations in parallel in the distributed memory of the high-performance computing node, monitors the overall quality residual in real time until it drops below the convergence threshold of 0.0001, and finally outputs and aggregates the surface heat flux map containing the distribution of the heated power on the plate surface, the internal moisture content gradient cloud map characterizing the moisture migration state, and the topological space wind pressure gradient line describing the airflow penetration resistance.

[0060] The specific implementation of the dynamic boundary compensation is as follows: a fixed Euler coordinate system is established to describe the airflow field inside the curing oven, and a Lagrange coordinate system that moves with the board is established to describe the temperature field and moisture content field of the board; within each time step Δt, the displacement of the board is calculated based on the board conveying speed v. The surface node coordinates of the plate in the Lagrange coordinate system are transformed to the Eulerian coordinate system, and the fluid element corresponding to the node in the Eulerian mesh is found. The fluid velocity, temperature and pressure at this location are extracted from the Eulerian mesh by trilinear interpolation and used as the convective heat transfer boundary condition of the plate surface. Through the above mapping, the continuous heating effect of the fixed nozzle jet on the moving plate surface is accurately simulated.

[0061] By bidirectionally decoupling and real-time nesting of the microscopic parameter response library and the macroscopic turbulence model, physical consistency modeling of the entire rock wool curing process was achieved. Coordinate transformation technology was used to solve the problem of boundary unsteady fluctuation caused by board displacement, providing a high-fidelity digital image of the production line for accurate tracing of subsequent quality defects.

[0062] Furthermore, the step of establishing the energy relationship between the integrated participating medium radiation term and the component migration term to describe the heat and mass transfer process inside and on the surface of the plate includes: A component migration term is introduced into the energy equation to characterize the effect of moisture content on apparent heat capacity, and a component migration and diffusion equation reflecting the moisture evaporation process is established. The energy equation and the component migration equation are coupled through the latent heat of phase change term. The absorption, emission, and scattering effects of radiation intensity in the fiber skeleton are calculated using the discrete coordinate method, and the calculated radiation source term is superimposed into the energy equation. Coupled boundary conditions based on convective heat transfer coefficient and convective mass transfer coefficient are established at the contact interface between the board and the circulating air to calculate the convective heat transfer and moisture evaporation rate on the board surface. The energy equation, component migration equation, and radiation transfer equation are solved alternately in each calculation step until the temperature field and moisture content field inside the board converge.

[0063] Specifically, the system first extracts the local moisture content parameter of 10% of the rock wool board collected in real time by the sensor, and determines the specific heat capacity of water (4200 joules per kilogram of Celsius) and the specific heat capacity of the fiber skeleton (800 joules per kilogram of Celsius) by retrieving the physical property database. The above components are synthesized using a volume weighted average operator and injected into the unsteady-state term of the energy equation to characterize the change in apparent heat capacity of the board caused by the dynamic decrease in moisture content in real time.

[0064] Specifically, the algorithm logic constructs a component migration and diffusion equation describing the movement of water within micropores. This equation calculates the water migration flux due to the concentration gradient by reading the temperature field data of the current computational grid and combining it with an effective diffusion coefficient of 0.00002 m² / s. It also uses a latent heat of phase change of 2,260,000 joules per kilogram as a negative source term for energy coupling to deduct the heat of vaporization absorbed by water evaporation in real time.

[0065] Specifically, the system applies the discrete coordinate method to perform spatial deconstruction of the radiation intensity in the three-dimensional computational domain, divides the 4π spatial solid angle into 8×8 discrete directions, and solves the radiation transfer equation in each direction using the finite volume method. By inputting the absorption coefficient value of 0.5 and the scattering coefficient value of 0.3 of the fiber skeleton, the energy deposition of radiation energy in the porous medium is quantified, thereby calculating and generating the radiation source scalar field superimposed on the energy equation.

[0066] Specifically, in order to accurately characterize the coupled boundary conditions on the surface of the plate, the system introduces a boundary mass transfer prediction model based on a multilayer perceptron (MLP) architecture. This model architecture includes a 3D input layer for receiving local wind speed of 15 m / s, surface temperature of 150 degrees Celsius, and moisture content gradient, three hidden layers each containing 128 neurons and using the ReLU activation function, and an output layer. The Adam optimizer is used to perform 2000 iterations on a training set containing 5000 sets of porous media convection heat transfer experimental samples to lock the weights.

[0067] Specifically, the boundary mass transfer prediction model generates a local convection mass transfer coefficient of 0.05 m / s in the output layer by aggregating the nonlinear features of the hidden layer neurons. Based on this, the system establishes coupled boundary conditions that include Newton's cooling law and mass diffusion flux, and calculates the sensible heat exchange and moisture evaporation rate between the plate surface and the circulating air in real time.

[0068] Specifically, the computing engine executes alternating solution logic within each preset 50-millisecond simulation step. First, it uses the physical field distribution of the previous moment as the initial condition to solve the component migration equation. Then, it corrects the physical property matrix in the energy equation according to the updated local water content and drives the radiation transfer equation to perform energy balance iteration in discrete directions.

[0069] Specifically, the algorithm monitors the L2 norm residual distribution of the temperature field and moisture content field between two consecutive inner iterations in real time. When the convergence criterion of the highest temperature rise deviation of the entire grid points being less than 0.01 degrees Celsius and the moisture content residual being less than 0.0001 is identified, the calculation of the current time step is stopped, thereby outputting the internal physical state distribution data of the board with thermal-mass balance consistency.

[0070] The specific implementation of the discrete coordinate method includes: using a specific directional discretization format to divide the radiation transmission path of the entire three-dimensional space into dozens of independent discrete directions; ensuring uniform sampling of radiation energy in three-dimensional space by uniformly dividing the zenith angle and azimuth angle into intervals; establishing a radiation intensity evolution logic in each discrete direction to describe the physical changes of radiation energy when passing through the rock wool medium. This logic covers the absorption effect of the medium on radiation energy, its own outward emission effect, and the energy extinction and scattering effects caused by the fiber skeleton; configuring preset absorption and scattering coefficients, and setting a scattering phase function based on the isotropic assumption, that is, assuming that radiation energy is uniformly scattered in all directions after colliding with the fiber skeleton; using gray body surface boundary conditions for the inner wall of the curing furnace and the surface of the board, setting uniform emissivity and reflectivity parameters, and describing the surface reflection behavior as a diffuse reflection process; calculating the local radiation vector divergence by spatially integrating the radiation intensity in all discrete directions, thereby characterizing the net contribution of radiative heat flow to the energy equation, and superimposing it as a heat source term into the macroscopic energy conservation equation, thus realizing the closed-loop coupling solution of the radiation field and temperature field.

[0071] By establishing a full-order distribution field of the production line based on the SST turbulence model and parameter response library, the simulation problem of unsteady heat and mass exchange during the movement of large-size plates was solved. The boundary mass transfer coefficient was dynamically corrected using a multilayer perceptron model, ensuring that the digital twin mirror image can maintain more than 98% synchronization accuracy with the physical production situation at different moisture content stages.

[0072] In describing the implementation steps of heat and mass exchange within the board material, the system establishes a dynamically coupled component migration mechanism. Specifically, based on the real-time porosity of the board material at different processing stages, the effective diffusion coefficient of moisture in the porous fiber skeleton is determined using Brügmann correction logic. For example, when the porosity fluctuates within the range of 0.85 to 0.95, the effective diffusion coefficient is dynamically adjusted according to the 1.5th power of the porosity to characterize the fiber's hindering effect on moisture diffusion.

[0073] Furthermore, the algorithm monitors the temperature status of each grid node inside the sheet in real time. When the local temperature reaches or exceeds the saturation temperature, the system automatically triggers the phase change evaporation model. At this time, the program calculates the evaporation rate using the difference between the current local vapor pressure and the ambient partial pressure, and introduces the latent heat of vaporization constant (valued at 2260 kJ / kg) to convert the mass flow into a negative source term (heat sink) in the energy equation. In this way, the system achieves the physical conservation constraints between moisture content change, moisture evaporation, and temperature field distribution without directly listing partial differential equations, ensuring that the error between the simulation results and actual production data of the drying section of the production line is controlled within 3%.

[0074] In the implementation of the full-order flow field solution, this embodiment made adaptive modifications to the SST k-omega model to address the interaction between the high-speed nozzle jet and the rock wool board inside the curing furnace, taking into account the porous media characteristics. Specifically, for the metal wall region of the curing furnace, a low Reynolds number mesh refinement technique was adopted to ensure that the dimensionless distance from the center of the first mesh layer to the wall is less than 1, thereby directly capturing the laminar sublayer characteristics of the boundary layer without relying on a general wall function.

[0075] For the flow of porous media within the plate, the system transforms the anisotropic permeability tensor retrieved from the parameter response library into a momentum dissipation source term. In the computational logic, this source term is composed of the superposition of a linear viscous drag term and a nonlinear inertial drag term: the linear part is determined by the product of the fluid kinematic viscosity and the inverse matrix of the permeability tensor, while the nonlinear part is determined based on the mapping relationship between the inertia factor and the square of the local flow velocity. Simultaneously, to compensate for the suppression of turbulent fluctuations by the fiber skeleton, the program introduces an additional dissipation operator based on porosity correction into the turbulent kinetic energy transport equation in the porous media region. This enforces constraints on the pressure drop characteristics and energy loss of the fluid as it passes through the plate, thereby ensuring the physical rigor of the multiphysics simulation under complex wind field environments.

[0076] Furthermore, the spatial basis functions of the full-order distribution field of the production line are extracted using intrinsic orthogonal decomposition, and the control equations are projected onto the feature subspace using discrete empirical interpolation and Galerkin projection to obtain a low-dimensional feature space. Specifically, the system first extracts the full-order distribution field data of the production line under the past 500 simulation time steps output by the preceding module, including the velocity, temperature and moisture content components of each grid node, and arranges the 100,000 degrees of freedom vectors of each time step into columns to construct a feature snapshot matrix with a dimension of 100,000×500.

[0077] Specifically, the algorithm logic calls the intrinsic orthogonal decomposition operator to perform singular value decomposition on the snapshot matrix. By calculating the energy ratio of the diagonal elements of the eigenvalue matrix, and according to the criterion that the cumulative contribution rate of the eigenvalues ​​exceeds 99.5%, the top 20 orthogonal eigenvectors from the left singular vector matrix are selected as spatial basis functions.

[0078] Specifically, in order to achieve discrete representation of the principal components of physical evolution, the system introduces a feature enhancement model based on a convolutional autoencoder. The model architecture consists of a decoder composed of one input layer, three convolutional layers, one bottleneck layer with 20 neurons, and three transposed convolutional layers.

[0079] Specifically, the convolutional layers are configured with 16, 32, and 64 convolutional kernels of size 3×3, respectively, and employ the ReLU activation function. They are trained on a historical snapshot set with an initial learning rate of 0.001 using the Adam optimizer, and the output spatial basis function is projected into a 20-dimensional low-dimensional linear subspace.

[0080] Specifically, the system automatically identifies the nonlinear convection terms in the governing equations of the SST turbulence model and the fourth-order radiation source terms in the radiation model of the participating medium, and calls the discrete empirical interpolation method to perform sampling point optimization in the full-order physical grid space.

[0081] Specifically, the algorithm logic adopts a greedy search strategy, selecting the grid index position that maximizes the second norm of the approximation error of the current nonlinear operator in each iteration step, and locking 100 key sparse sampling points from 100,000 original grids as DEIM metric points.

[0082] Specifically, the system constructs an approximate projection basis matrix for the nonlinear operator based on the physical field values ​​at the sparse sampling points and the aforementioned generated 20th-order spatial basis functions.

[0083] Specifically, the system uses the Galerkin projection method to map the partial differential operators in the full-order momentum equation and energy equation to a low-dimensional linear subspace defined by 20 basis functions. By performing orthogonal dot product operations between the full-order residual vector and the spatial basis functions, the original nonlinear algebraic equation system with 100,000 unknowns is simplified to an ordinary differential equation system containing only 20 generalized coordinate coefficients.

[0084] Specifically, the algorithm logic uses the output approximate projection basis to perform discrete evaluation of the nonlinear terms in the equation on a sample-by-sample basis, thereby eliminating the strong dependence of computational complexity on the number of full-order grids, and finally generating a low-dimensional feature space with significantly compressed dimensions in memory.

[0085] Specifically, the reduced-order model parameters in the low-dimensional feature space are sent to the real-time evolution prediction module in binary tensor format as the core dynamic operator driving the rapid iteration of the digital twin.

[0086] The snapshot acquisition strategy is as follows: During the complete cycle from cold start to steady-state operation of the production line, a full-order distribution field is acquired every 50ms, for a total of 500 snapshots, covering the entire operating conditions of the production line, including start-up, acceleration, steady-state, and deceleration. Each column of the snapshot matrix corresponds to the physical quantity of the entire field at a given moment. For each grid node, its three components of velocity vector, temperature scalar, and moisture content scalar, a total of five variables, are arranged sequentially to form a column vector containing 100,000 grid nodes × 5 physical quantities = 500,000 degrees of freedom. The number of basis functions is selected as follows: After performing singular value decomposition on the snapshot matrix, the energy percentage of each singular value is calculated. When the cumulative energy percentage first exceeds 99.5%, the corresponding order N is recorded, and the number of basis functions is set to N, thus determining it to be 20th order.

[0087] The method of selecting sparse sampling points in physical space using discrete empirical interpolation includes: calculating the energy gradient distribution of the nonlinear residuals of the control equations at the current time step; extracting the spatial region where the residual energy gradient exceeds a preset threshold as a high-sensitivity sensitive region; dynamically increasing the distribution density of sparse sampling points in the high-sensitivity sensitive region and synchronously updating the approximate projection basis of the nonlinear operator to compensate for the projection error caused by the nonlinear abrupt change of the radiation term.

[0088] In simulating the jet region of the curing oven nozzle, the nonlinear residual energy gradient surges due to severe pressure pulsations. The discrete empirical interpolation algorithm identifies the increased sensitivity weights in this region and automatically increases the sampling point ratio of this local area from the initial 5.0% to 12.5%, reconstructing the projection operator for the turbulence source term in real time to ensure that the reduced-order model does not become unstable under strong disturbance conditions. By introducing an adaptive sampling mechanism driven by the residual energy gradient, the problem of accuracy collapse of the traditional static DEIM algorithm when dealing with local strong nonlinear disturbances (such as jet impact and radiation abrupt changes) is solved, and the relative error of the reduced-order model under high disturbances is controlled within 3.0%.

[0089] By constructing a reduced-order model architecture based on POD-DEIM, the technical bottleneck of real-time simulation of complex porous media heat and mass exchange processes was solved. The sparse sampling and subspace projection techniques ensured the accurate inheritance of the physical laws of the full order by the low-dimensional model, providing efficient computing power support for online quality monitoring of large-scale rock wool production lines.

[0090] Furthermore, by utilizing a parameterized background data weak observer, the discrete measurement point data collected by the sensor is used as a constraint term, and variational assimilation calculation is performed in the low-dimensional feature space to extract the initial modal coefficients that match the current measured working conditions.

[0091] Specifically, the system obtains the temperature values ​​of 245.5 degrees Celsius from 24 thermocouple sensors distributed along the length of the curing oven in real time through an industrial data gateway, and encapsulates them into a discrete observation vector with a dimension of 24×1. At the same time, it retrieves the 20th-order spatial basis function matrix generated by the preceding module from the memory buffer.

[0092] Specifically, the parameterized background data weak observer adopts a cost function minimization logic based on the variational principle. Its algorithm process involves constructing an augmented cost function containing background field bias terms and observation residual terms within a low-dimensional linear subspace defined by 20 generalized coordinates. The gradient descent operator is used to search for 20 weighting coefficients with a step size of 0.05 that minimize the sum of the squared Euclidean distances between the predicted and measured values ​​at the 24 physical sampling points.

[0093] Specifically, the system performs the adjoint equation solution in the arithmetic logic unit of the high-performance central processing unit. By calculating the partial derivative of the cost function with respect to the generalized coordinates, it completes 50 iterations within 10ms to eliminate transient disturbances introduced by sudden wind speed fluctuations in the environment, thereby locking in a set of initial modal coefficients that characterize the optimal projection state of the current physical field.

[0094] Reference Figure 3 The present invention provides a schematic diagram of the process for solving the full-order distribution field and generating a self-correcting twin field. Furthermore, a Gaussian process regression model is used to predict the nonlinear residual between the low-dimensional feature space and the full-order distribution field, thereby generating a calibration operator for compensating for unmodeled physical effects of the system.

[0095] Specifically, the Gaussian process regression model adopts a stochastic process modeling architecture based on the Matern 5×2 kernel function. This model architecture includes an input layer for receiving 20-dimensional initial modal coefficients, a kernel function mapping layer for calculating the covariance matrix, and an output layer for outputting the residual prediction mean and variance.

[0096] Specifically, the system uses the Adam optimizer to perform maximum likelihood estimation training on a historical dataset containing 10,000 pairs of residuals from "full-order simulation field - reduced-order simulation field". Its internal logic is to calculate the distance measure of the input sample points in the 20-dimensional manifold space and use the inversion operation of the covariance matrix to achieve nonlinear mapping of unmodeled physical effects caused by subtle turbulent fluctuations or material constitutive approximations.

[0097] Specifically, the algorithm receives 20 initial modal coefficients as query variables in each simulation step and calculates and outputs a residual compensation vector with a dimension of 20×1 as a calibration operator, where the compensation component of the first-order mode has a value of 0.015, which is used to accurately correct the systematic loss of the reduced-order model in the energy dissipation term.

[0098] By leveraging the uncertainty quantification capability of Gaussian process regression, this process achieves intelligent extraction and probabilistic characterization of order reduction error, effectively compensating for the generalization failure problem of physical mechanism models under extremely complex working conditions.

[0099] Furthermore, the initial modal coefficients are linearly combined with the calibration operator using the least squares fitting operator to reconstruct the corrected dynamic modal coefficients.

[0100] Specifically, the system constructs a linear fusion operator based on weighted average logic. By reading the trust weight parameters preset in the configuration register, the weight of the initial modal coefficient is set to 0.85 and the weight of the calibration operator is set to 0.15.

[0101] Specifically, the computation unit performs scalar multiplication and vector addition operations on the weight coefficients and corresponding vectors, and then inputs the calculation results into a second-order least squares fitter. By performing parameter optimization based on minimizing the sum of squared residuals within a sliding window in the time domain, the numerical oscillations generated during mode switching are dynamically eliminated.

[0102] Specifically, the algorithm logic outputs a set of smoothed, 20-dimensional modified dynamic mode coefficients in real time and maps them to the real-time variable area in the video memory.

[0103] By combining and fitting the virtual and real features linearly, the system achieves smooth evolution of the dynamic trajectory while ensuring computational efficiency, thus ensuring the continuity of calibration parameters over time.

[0104] Furthermore, the corrected dynamic mode coefficients are projected back to the spatial basis functions to synthesize a self-correcting twin field with physical conservation properties and real-time bias correction.

[0105] Specifically, the system retrieves a 100000×20 dimension spatial basis function matrix from the graphics processor's video memory and performs a matrix multiplication operation between the basis function matrix and the corrected 20-dimensional dynamic mode coefficient vector.

[0106] Specifically, the algorithm restores the physical attribute values ​​(such as temperature and moisture content) of each grid node to specific values ​​in all dimensions, generating a digital cloud map containing 100,000 discrete nodes.

[0107] Specifically, the system uses the Laplace smoothing operator to perform boundary consistency correction on the synthesized distribution field, ensuring that the output self-correcting twin field satisfies physical conservation laws such as the Navier-Stokes equations, while keeping the error between the output and the measured values ​​at 24 physical measurement points within an extremely narrow envelope of less than 1.5%.

[0108] Specifically, the self-correcting twin field is pushed to the visualization monitoring terminal in three-dimensional voxel format to display the real physical distribution of the 4% moisture content gradient and the 150-degree Celsius temperature field evolution inside the rock wool board in real time.

[0109] Specifically, the specific steps of the variational assimilation calculation include: Initial values ​​for the modal coefficients to be optimized are set, and the background field error covariance matrix and the sensor observation error covariance matrix are preset based on system experience. The background field error is used to limit the model's confidence in the preset values, and the observation error is used to characterize the level of random noise during sensor acquisition. The cost function is solved using the gradient descent algorithm. In each iteration, the rate of change of the cost function relative to the modal coefficients (i.e., the gradient) is calculated. This gradient comprehensively reflects the background field deviation pressure and the observation field fitting pressure. The modal coefficient vector is continuously corrected along the opposite direction of the gradient according to a set step size. After each iteration, the second norm of the gradient vector is detected. When this value drops to a preset minimum threshold, the calculation is considered to have reached convergence. The iterative calculation is terminated, and the optimal modal coefficients that satisfy the error weighted sum minimization are output as key driving parameters for subsequent reconstruction of the self-correcting twin field.

[0110] The method of using a Gaussian process regression model to predict nonlinear residuals includes: injecting the local extrema of heat flux calculated based on the participating medium radiation model as a priori physical constraint into the kernel function of the Gaussian process regression model; using the priori physical constraint to limit the search space of the regression output, so that the generated calibration operator, under the premise of satisfying the fitting of discrete sensor data, is forced to follow the heat flow continuity constraint in the law of conservation of energy.

[0111] When predicting the residual heat on the surface of the plate, the Gaussian process regression model uses the heat density distribution calculated by the radiation model as the benchmark mean function. By setting a physical soft constraint operator, when the heat flux fluctuation predicted by the calibration operator exceeds 1.5 times the physical maximum heat flux, the kernel function automatically performs a penalty weight reduction, so that the output result returns to the physical solution space that conforms to the law of energy conservation. By embedding prior thermodynamic constraints into the kernel of the regression model, the problem of artifacts that violate physical laws are easily generated by the pure data-driven calibration model when the sensor data is sparse or the noise is large is solved, ensuring the robustness and physical reliability of the twin field correction process.

[0112] The training steps of the Gaussian process regression model include: In the offline stage, 10,000 full-order simulations are performed within a reasonable range by changing the process parameters of the production line (such as wind speed, pressure, and conveying speed). For each simulation, a reduced-order model is run simultaneously to obtain the reduced-order field. The difference between the physical quantities of the full-order field and the reduced-order field at each grid point is calculated, and this difference is projected onto a 20th-order spatial basis function to obtain a 20-dimensional residual mode coefficient vector. The input 20-dimensional initial mode coefficients and the output 20-dimensional residual mode coefficients are used to form training sample pairs, totaling 10,000 pairs. The Matrn 5 / 2 kernel function is used. By maximizing the marginal likelihood function, the hyperparameters of the kernel function are optimized using the Adam optimizer to obtain a length scale l=1.5 and a signal variance σ_f^2=0.01. In the online stage, given the current initial mode coefficients, the trained Gaussian process model is used to predict the residual mode coefficients, which are then used as calibration operators.

[0113] By constructing a collaborative correction mechanism between a parameterized background data weak observer and a Gaussian process regression model, deep assimilation of multi-source sensor data and low-dimensional dynamic models was achieved. Under the premise of ensuring a computational latency of less than 50ms, the problem of prediction drift that is easily generated by reduced-order models under dynamic disturbances was overcome, providing real-time full-field data support with physical conservation constraints for precise quality control in the rock wool production process.

[0114] Furthermore, the self-correcting twin field is a real-time dynamic distribution field with physical conservation constraints. Its synthesis process includes: projecting the corrected dynamic mode coefficients back to the characteristic subspace defined by the spatial basis function, restoring the set of grid node values ​​under the full-order computational domain through matrix linear reconstruction, and presenting the temperature field evolution process and moisture content gradient distribution inside the rock wool board in real time on the monitoring terminal in the form of a three-dimensional voxel cloud map.

[0115] Specifically, the system extracts the corrected dynamic mode coefficients from the preceding steps. These coefficients are represented as a 20×1 column vector, where each element precisely corresponds to the weight contribution value of the first-order intrinsic mode.

[0116] Specifically, the system retrieves the spatial basis function matrix generated by intrinsic orthogonal decomposition from the dedicated video memory partition of the graphics processor. This matrix has a dimension of 100000×20, and each column represents a physical distribution mode with a specific spatial frequency.

[0117] Specifically, the matrix linear reconstruction process is performed in the hardware multiplier array of the computation kernel. By performing matrix multiplication operations on a 100000×20 spatial basis function matrix and a 20×1 dynamic modal coefficient column vector, that is, by weighted summation of the 20 basis function components at each spatial location, the abstract coefficients in the feature subspace are restored to the values ​​of 100000 grid nodes in the full-order computation domain.

[0118] Specifically, in order to ensure that the reconstructed field data conforms to the physical conservation law, the system introduces a physical conservation verification model based on the convolutional neural network (CNN) architecture. This model includes one three-dimensional input layer, two hidden convolutional layers with 64 feature channels each, and one output layer. The convolutional kernel size is set to 3×3×3 with a stride of 1.

[0119] Specifically, the physical conservation verification model is trained on a set of 10,000 fluid snapshots that conform to the conservation of mass and energy using the Adam optimizer. The numerical gradient between grid nodes is extracted using a convolution operator, and the residual of the mass source term in the whole field is calculated through a custom physical loss function layer. This forces the outlier node values ​​that deviate from the physical evolution trajectory to be corrected to the theoretical envelope range.

[0120] Specifically, the system maps the verified set of grid node values ​​back to physical space according to their original three-dimensional coordinate indices, and constructs a structured three-dimensional voxel tensor with a resolution of 256×128×128.

[0121] Specifically, the self-correcting twin field synthesis engine receives the voxel tensor, performs spatial interpolation and pseudo-color mapping using a volume rendering algorithm, maps the temperature range of 25 degrees Celsius to 250 degrees Celsius into a chromatographic sequence from dark blue to bright red according to the numerical value, and simultaneously maps the moisture content gradient from 0% to 10% into a cloud layer with varying transparency.

[0122] Specifically, the synthesized self-correcting twin field data is pushed to the monitoring terminal in real time through a 10 Gigabit Ethernet interface, and the internal temperature rise evolution curve and the three-dimensional dynamic trajectory of moisture migration of the rock wool board as it moves in the curing furnace are presented on the graphics display unit at a refresh rate of 20 frames per second.

[0123] This embodiment constructs a complete digital twin evolution system, from microscopic property response to the full-order distribution field of the production line, by deeply coupling gridded Boltzmann microscopic simulation, SST turbulence macroscopic modeling, and POD-DEIM order reduction technology. A real-time feedback mechanism established using a parameterized background data weak observer and a Gaussian process regression model solves the problem of prediction trajectory drift that traditional order reduction models easily encounter when facing dynamic disturbances in production, ensuring a high degree of consistency in statistical distribution between the simulation baseline and the physically measured state. Through cross-scale correlation extraction and closed-loop correction, the system's real-time quantification capability of internal quality indicators in the rock wool solidification process is enhanced, and finally, physical conservation constraints ensure the authenticity and robustness of the twin field data.

[0124] Example 2: This embodiment is set in a large-scale flexible manufacturing scenario of rock wool boards with multiple production lines working together. In response to the microscopic fiber orientation mutations that occur during the continuous curing process of rock wool products of different densities and the risk of local thermal runaway caused by the exothermic curing of phenolic resin, a cross-scale feature deconstruction model based on variational information bottleneck is introduced to perform deep correlation analysis of physical patterns.

[0125] Furthermore, the dynamic time warping algorithm is used to correct the temporal phase difference between the spatial gradient of the environmental field and the micro-rheological features. Cross-scale manifold alignment is performed through the variational information bottleneck, and conjugate latent vectors are extracted.

[0126] Specifically, the system extracts the transient temperature gradient sequence from the three-dimensional voxel mesh from the full-order distribution field of the production line, and simultaneously extracts the shear stress evolution sequence of the fiber contact point from the parameterized structural model. Both sets of data are encapsulated with a sampling interval of 10ms.

[0127] Specifically, the system constructs a phase correction operator based on the dynamic time warping algorithm. By establishing a 500×500 dimension cost matrix in memory, the cosine distance between the feature vectors of the sampling points between the two sequences is used as the local cost. The dynamic programming search algorithm is applied to identify the nonlinear warping path with the minimum cumulative residual, thereby offsetting the approximately 250ms time-domain phase offset caused by thermal conduction hysteresis.

[0128] Specifically, a feature deconstruction model based on the physical constraint variational autoencoder (PC-VAE) architecture is constructed. This model architecture consists of an environment feature encoder with three convolutional layers, a dynamic feature encoder with two layers of GRU recurrent units, and a common reparameterization layer.

[0129] Specifically, the environmental feature encoder is equipped with 32, 64, and 128 three-dimensional convolutional kernels of size 3×3×3 to extract heat flux distribution texture from full-order field cloud maps, and the dynamic feature encoder is equipped with 256 neurons to capture the long-range dependence of stress fluctuations. The hybrid model is trained for 2000 rounds on a large dataset containing 100,000 production snapshots using the Adam optimizer until the mutual information term converges.

[0130] Specifically, the algorithm logic introduces a variational information bottleneck operator at the end of the encoder. By setting a KL divergence penalty term with a weight of 0.01 in the loss function, random background noise unrelated to the physical mechanism is forcibly filtered out, thereby extracting conjugate latent vectors that can characterize the essence of physical evolution in the 64-dimensional manifold space.

[0131] Furthermore, within the probabilistic latent space established with Fisher information metric, the geodesic displacement of the conjugate latent vector deviating from the standard physical steady-state distribution is calculated in real time, and non-steady-state modes are identified in combination with physical hard constraint boundaries.

[0132] Specifically, the system constructs a Riemannian metric operator based on the Fisher information matrix within the probabilistic latent space. Each element in this matrix is ​​calculated by the expected value of the second-order partial derivative of the log-likelihood function with respect to the latent variable, which is used to define the geometric properties of the anisotropy of the latent space.

[0133] Specifically, the algorithm logic extracts the coordinate components of the current conjugate latent vector in each calculation cycle, searches for the nearest standard physical steady-state centroid that conforms to the Hamiltonian mechanics preset, calculates the geodesic distance between the two points that satisfies the energy conservation constraint by performing path integral operation on the Riemannian manifold, and outputs a scalar displacement value of 0.75.

[0134] Specifically, the system retrieves a preset physical hard constraint boundary, which is defined by the maximum tensile modulus of rock wool fiber (1.5 GPa) and the upper limit of the resin curing exothermic rate (200 W / kg). When the rate of change of the first derivative of the geodesic displacement exceeds 0.5 within 5 consecutive sampling points and the latent vector component overflows the boundary range, the current working condition is determined to enter an unsteady state mode.

[0135] Specifically, by using the Jacobian matrix inverse mapping to map the associated features back to the physical dimension, and coupling the local stress concentration points and the circumferential heat release rate under the higher-order tensor, a latent image of quality risk with a spatiotemporal evolution probability distribution is generated.

[0136] Specifically, the system extracts the Jacobian matrix corresponding to the variational autoencoder decoder layer and transforms the unsteady offset in the latent space into parameter increment components at the physical level by performing matrix inversion.

[0137] Specifically, the algorithm logic identifies the core process cause of instability as a 15% deviation in wind pressure pulsation by calculating the sensitivity of the partial derivatives of each parameter increment to the curing oven wind speed and pressure roller pressure, and couples this feature into the three-dimensional voxel model.

[0138] Specifically, the image rendering engine receives the failure probability values ​​of each grid node from the calculation output, uses volume rendering technology to generate a quality risk latent image where the pixel brightness represents the risk intensity, marks the fiber fracture area inside the board caused by uneven density distribution in real time, and pushes closed-loop intervention instructions to the DCS system.

[0139] This embodiment achieves sub-millisecond-level deconstruction and tracing of microscopic physical patterns in ultra-large-scale rock wool production lines by introducing a deep learning architecture with dynamic time warping and Fisher information metric. It solves the information redundancy problem in cross-scale feature alignment by utilizing variational information bottlenecks, and achieves a closed-loop mapping from abstract latent features to specific process parameters through Jacobian matrix inverse mapping. This ensures that the system's advanced prediction accuracy for fiber breakage and lint risks reaches over 96%, realizing digital control of high-performance porous materials.

[0140] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A multi-physics simulation method for integrated rock wool panel composite production line based on digital twins, characterized in that, include: Using a random fiber filling algorithm and fractal geometry theory, a representative volume element reflecting the microstructure of rock wool fibers is constructed, and a parameterized structural model based on the porosity evolution law is established. A parameter response library was constructed by simulating the heat-fluid coupling field in micropores using the lattice Boltzmann method and performing volume averaging. use The SST turbulence model simulates the high-speed circulating air jet in the curing oven, retrieves the parameter response library as the variable parameter source term of the porous medium region, integrates and participates in the medium radiation model and the migration of porous medium components, and combines dynamic boundary compensation to map the displacement of the board on the production line to obtain the full-order distribution field of the production line. The spatial basis functions of the full-order distribution field of the production line are extracted using intrinsic orthogonal decomposition. When the cumulative energy ratio first exceeds 99.5%, the corresponding order N is recorded, and the number of basis functions is N. The energy gradient distribution of the nonlinear residual of the control equation at the current time step is calculated. When the energy gradient of the nonlinear residual surges, the discrete empirical interpolation algorithm identifies that the sensitivity weight of this region is increased, and automatically increases the sampling point ratio of this local region from the initial 5.0% to 12.5%. The control equation is projected to the feature subspace through discrete empirical interpolation and Galerkin projection to obtain a low-dimensional feature space. Using a parameterized background data weak observer and Gaussian process regression, sensor data is projected onto the low-dimensional feature space for coefficient reconstruction and least squares fitting to obtain a self-correcting twin field. The parameterized background data weak observer constructs an augmented cost function containing background field bias and observation residual terms within a low-dimensional linear subspace defined by 20 generalized coordinates. It uses a gradient descent operator to search with a step size of 0.05, completing 50 iterations within 10ms to detect the magnitude of the second norm of the gradient vector. When this value drops to the preset minimum threshold ε, the calculation is determined to have reached a convergence state; The step of obtaining the self-correcting twin field includes: using a parameterized background data weak observer, taking the discrete measurement point data collected by the sensor as a constraint term, and performing variational assimilation calculation in the low-dimensional feature space to extract the initial modal coefficients that match the current measured working conditions; The nonlinear residuals between the low-dimensional feature space and the full-order distribution field are predicted using a Gaussian process regression model to generate a calibration operator for compensating for unmodeled physical effects in the system. The initial modal coefficients are linearly combined with the calibration operator using a least-squares fitting operator to reconstruct the corrected dynamic modal coefficients. The corrected dynamic modal coefficients are projected back onto the spatial basis functions to synthesize a self-correcting twin field with physical conservation properties and real-time bias correction.

2. The multiphysics simulation method for integrated rock wool panel composite production line based on digital twins according to claim 1, characterized in that, The method of constructing a representative volume element reflecting the microstructure of rock wool fibers using a random fiber filling algorithm and fractal geometry theory, and establishing a parameterized structural model based on the porosity evolution law, includes the following steps: Based on the average diameter and aspect ratio distribution of rock wool fibers, an initial set of fiber central axes is generated in a set of three-dimensional Euclidean spaces. The anisotropic distribution of the central axis set is controlled by a fiber orientation probability density function set based on the frequency of the production line pendulum. Combined with fractal geometry theory, multi-level bifurcation and surface roughness iteration are performed on individual fibers to generate fiber geometric entities with interpenetrating structures, constituting the initial representative volume element. Using the real-time pressure distribution of the production line pressure rollers as the driving variable, a geometric compression transformation operator describing the spatial deformation of the fiber skeleton is established. The initial representative volume element is then spatially scaled and mapped using the geometric compression transformation operator, dynamically adjusting the number of contact points and the porosity between fibers to obtain a parameterized structural model that evolves in real-time with the porosity under production line conditions.

3. The multiphysics simulation method for integrated rock wool panel composite production line based on digital twins according to claim 1, characterized in that, The steps for constructing the parameter response library include: The representative volume element is discretized into a three-dimensional grid space, and velocity distribution functions describing fluid motion and temperature distribution functions describing energy transfer are constructed respectively. Pressure gradients and temperature difference boundaries are applied to the grid space, and iterative calculations are performed using collision and migration operators until the velocity and temperature fields within the micropores converge. The local forces and heat fluxes at the fluid-solid interface are calculated using the momentum exchange method. The converged microscale field is spatially integrated and averaged. The anisotropic permeability tensor is obtained based on the linear mapping relationship between the average filtration velocity and the pressure gradient. The equivalent thermal conductivity tensor is obtained based on the mapping relationship between the heat flux density distribution and the macroscopic temperature gradient. The fiber-fluid heat transfer coefficient is obtained based on the ratio of the total heat transfer at the fluid-solid interface to the average temperature difference. By changing the porosity of the parameterized structural model, the above steps are repeated to establish a correspondence table with porosity as the independent variable and the anisotropic permeability tensor, the equivalent thermal conductivity tensor, and the fiber-fluid heat transfer coefficient as dependent variables, thus forming the parameter response library.

4. The multiphysics simulation method for integrated rock wool panel composite production line based on digital twins according to claim 1, characterized in that, The steps for obtaining the full-order distribution field of the production line include: A continuity equation and momentum equation based on the SST turbulence model are established to describe the wind speed field inside the curing oven. Simultaneously, the energy relationship between the radiation term and component migration term of the integrated participating medium is established to describe the heat and mass exchange process inside and on the surface of the board. The board region is defined as a porous medium region. According to the current porosity of the parameterized structural model, the corresponding values ​​in the parameter response library are retrieved in real time as the resistance source term in the momentum equation and the heat conduction and heat exchange correction terms in the energy equation. According to the actual conveying speed of the production line, the geometric position of the board is updated using a coordinate transformation operator to compensate the flow field information in the curing oven to the moving board surface boundary. Then, full-order numerical solutions are performed, and the output is aggregated to obtain the full-order distribution field of the production line, which includes the surface heat flux spectrum, the internal moisture content gradient, and the topological space wind pressure gradient line.

5. The multiphysics simulation method for integrated rock wool panel composite production line based on digital twins according to claim 1, characterized in that, The steps for obtaining the low-dimensional feature space include: The full-order distribution fields of the production line at multiple time steps are combined into a snapshot matrix. Singular value decomposition is then performed on the snapshot matrix using eigenorthogonal decomposition to extract spatial basis functions reflecting the principal components of the production line's physical evolution, constructing a low-dimensional linear subspace. Nonlinear terms in the control equations are identified, and sparse sampling points are selected in the physical space using discrete empirical interpolation. An approximate projection basis for the nonlinear operator is constructed based on the spatial basis functions. Using the spatial basis functions and the approximate projection basis of the nonlinear operator, the full-order momentum and energy equations are mapped to the low-dimensional linear subspace via Galerkin projection, eliminating the dependency of computational complexity on the full-order mesh degrees of freedom, thus obtaining the low-dimensional feature space.

6. The multiphysics simulation method for integrated rock wool panel composite production line based on digital twins according to claim 1, characterized in that, The self-correcting twin field is a real-time dynamic distribution field with physical conservation constraints, and its synthesis process includes: The corrected dynamic mode coefficients are projected back to the characteristic subspace defined by the spatial basis functions, and the set of grid node values ​​under the full-order computational domain is obtained through matrix linear reconstruction. The temperature field evolution process and moisture content gradient distribution inside the rock wool board are presented in real time on the monitoring terminal in the form of a three-dimensional voxel cloud map.

Images