Intelligent simulation optimization method for extraction separation process based on digital twinning

Abstract

The application provides an extraction separation process intelligent simulation optimization method based on digital twinning, relates to the technical field of chemical process control, and comprises the following steps: establishing a digital twinning model and a numerical calculation model of an extraction separation physical device, and generating a simulation benchmark state; collecting real-time operation data and comparing the real-time operation data with the simulation benchmark state, identifying a deviation area, adjusting a local mass transfer coefficient, and solving to obtain a correction state field; constructing a coupling matrix of a concentration gradient and a velocity distribution based on the correction state field, calibrating a mass transfer dominant flow area through characteristic decomposition; setting a virtual disturbance source in the area to perform multi-working-condition simulation, and screening an optimal disturbance frequency; extracting a column segment residence time distribution variance, calibrating a flow deviation column segment, and calculating the spatial overlap degree of the flow deviation column segment and the mass transfer dominant flow area; and selecting and executing distributor opening rate or feed pulsation frequency adjustment according to the overlap degree. The application realizes accurate simulation and dynamic optimization of the extraction separation process, and improves the mass transfer efficiency and separation performance.

Description

Intelligent Simulation Optimization Method for Extraction and Separation Processes Based on Digital Twins Technical Field

[0001] This invention relates to chemical process control technology, and more particularly to an intelligent simulation optimization method for extraction and separation processes based on digital twins. Background Technology

[0002] Extraction separation technology is a core unit operation for material separation and purification in fields such as chemical engineering, pharmaceuticals, and nuclear fuel reprocessing. It achieves selective transfer of target components between two phases through two-phase contact. The optimization of traditional extraction separation processes mainly relies on empirical adjustments and offline experiments. Existing digital twin models mostly employ static or quasi-static modeling methods, lacking dynamic correction mechanisms based on real-time operational data. The flow field and concentration field inside the extraction equipment are affected by various factors such as operating conditions, physical property parameters, and equipment aging, exhibiting significant time-varying and spatial non-uniformity.

[0003] Traditional simulation models, once established, have fixed parameters and cannot adaptively adjust to actual operating conditions. This leads to a gradual decline in model prediction accuracy over time, resulting in significant deviations between simulation results and actual operating conditions, making it difficult to provide a reliable basis for optimization decisions. Traditional optimization methods often employ global adjustment strategies, failing to identify and implement local enhancement measures for mass transfer-dominant regions. This results in limited optimization effects and may even disrupt existing favorable flow structures due to blind adjustments, leading to decreased mass transfer efficiency. Different types of deviations require different control measures, such as adjusting structural parameters, optimizing operating parameters, or changing process conditions. Existing technologies often rely on trial-and-error adjustments based on experience, lacking intelligent decision-making mechanisms based on deviation characteristics. This results in long optimization cycles, significant resource waste, and difficulty in achieving rapid and accurate process control. Summary of the Invention

[0004] This invention provides an intelligent simulation optimization method for extraction and separation processes based on digital twins, which can solve the problems in the prior art.

[0005] A first aspect of this invention provides an intelligent simulation optimization method for extraction and separation processes based on digital twins, comprising:

[0006] Establish digital twin models and numerical calculation models of the extraction and separation physical equipment to generate simulation baseline states;

[0007] Real-time operating data of the extraction and separation physical equipment is collected, and the real-time operating data is spatially compared with the simulation baseline state. Spatial regions where the deviation exceeds the set deviation threshold are identified. The local mass transfer coefficient of the spatial region is adjusted and the numerical calculation model is updated to obtain the corrected state field.

[0008] Based on the corrected state field, a coupling matrix of concentration gradient and velocity distribution is constructed. The coupling matrix is ​​decomposed into features, the dominant feature vector is extracted and mapped to the device space, and the dominant flow region of mass transfer is calibrated.

[0009] A virtual disturbance source is set up in the mass transfer-dominant flow region. Multi-condition simulation is carried out by changing the disturbance frequency of the virtual disturbance source, and the disturbance frequency that maximizes the mass transfer flux is selected.

[0010] Based on the correction state field, the residence time distribution variance of the tower section is extracted. The tower section with variance exceeding the preset variance threshold is labeled as the flow offset tower section. The spatial overlap between the flow offset tower section and the mass transfer dominant flow region is calculated.

[0011] Based on the spatial overlap, the distributor opening ratio or feed pulsation frequency is adjusted, and the adjusted flow rate distribution and concentration distribution data are fed back to the digital twin model.

[0012] Establish digital twin models and numerical calculation models of the extraction and separation physical equipment, and generate simulation baseline states including:

[0013] Obtain the geometric structure parameters and process operation parameters of the extraction and separation physical equipment, construct a three-dimensional digital space model as a digital twin model, and generate a spatial mesh based on the digital twin model;

[0014] Multiple virtual measurement points are constructed in the spatial grid, and the virtual measurement points are mapped to the corresponding physical measurement point locations in the extraction and separation physical equipment;

[0015] Pressure and concentration data are collected at the physical measurement points, and the pressure and concentration data are then substituted back into the virtual measurement points. The pressure field distribution and concentration field distribution are constructed using a grid interpolation method.

[0016] The flow field velocity distribution and flow field direction distribution are calculated based on the pressure field distribution and substituted into the mass transfer calculation equation. The concentration field distribution and the mass transfer calculation equation are combined to construct a numerical calculation model. The boundary constraints of the numerical calculation model are set based on the process operating parameters.

[0017] The numerical calculation model is iteratively calculated, and the iteration stops when the calculation residual is less than the convergence criterion. The calculation result is output as the simulation baseline state.

[0018] By spatially comparing real-time operating data with the simulation baseline state, spatial regions where deviations exceed a set deviation threshold are identified. The local mass transfer coefficients of these spatial regions are adjusted, and the numerical calculation model is updated. The corrected state field is then obtained by solving for:

[0019] The deviation values ​​between the real-time running data and the simulation baseline state are calculated point by point on the spatial grid to generate a spatial deviation field. The spatial deviation field is then analyzed to extract connected regions whose deviation values ​​exceed a set deviation threshold and are marked as spatial regions.

[0020] Extract the spatial coordinate range of the spatial region, and obtain the local mass transfer coefficient of the spatial region from the solution domain of the numerical calculation model corresponding to the simulation baseline state based on the spatial coordinate range;

[0021] Calculate the difference in concentration distribution between real-time running data and simulation baseline state within the spatial region, integrate the concentration distribution difference along the flow direction, and calculate the mass transfer coefficient correction.

[0022] The mass transfer coefficient correction is superimposed with the local mass transfer coefficient to obtain the adjusted local mass transfer coefficient. The adjusted local mass transfer coefficient is then substituted into the mass transfer source term corresponding to the spatial region in the numerical calculation model, while keeping the mass transfer coefficient outside the spatial region unchanged, to obtain the updated numerical calculation model.

[0023] Apply boundary conditions corresponding to real-time running data to the updated numerical calculation model and perform numerical solutions. Terminate the iterative calculation when the solution residual is less than the set convergence threshold, and output the corrected state field.

[0024] A coupling matrix between concentration gradient and velocity distribution is constructed based on the corrected state field. The coupling matrix is ​​then subjected to eigenvalue decomposition to extract the dominant feature vector, which is mapped to the device space. The dominant mass transfer flow region is calibrated as follows:

[0025] The concentration field in the correction state field is differentiated according to spatial coordinates to generate a concentration gradient field. The velocity field in the correction state field is orthogonally decomposed to generate a velocity distribution field. Multi-scale tensor product operation is performed on the concentration gradient field and the velocity distribution field to generate a multi-level coupling matrix.

[0026] The multi-level coupling matrix is ​​block diagonalized to generate a block matrix group. Singular value decomposition is performed on each block matrix in the block matrix group and regularization constraints are applied to obtain a local dominant feature vector group. The local dominant feature vector group is then combined with adaptive weighting to generate a global dominant feature vector.

[0027] A nonlinear spatial mapping function is constructed to transform the global dominant feature vector into a feature distribution field in the physical coordinate system of the device through the nonlinear spatial mapping function, and the gradient value and curvature value of each spatial point in the feature distribution field are calculated.

[0028] The characteristic distribution field is spatially partitioned, and the connected regions where both the gradient value and curvature value are greater than the partition mean are designated as the mass transfer-dominant flow regions.

[0029] A virtual disturbance source is set up in the mass transfer-dominant flow region. Multi-condition simulations are performed by changing the disturbance frequency of the virtual disturbance source. The disturbance frequencies that maximize the mass transfer flux are selected from the following:

[0030] A spatial discrete grid is constructed in the mass transfer-dominant flow region, and the nodes of the spatial discrete grid are set as virtual disturbance sources. The velocity field distribution and concentration field distribution of the mass transfer-dominant flow region are extracted, and the mass transfer channel distribution is obtained through gradient calculation.

[0031] The location of the vortex center is identified based on the mass transfer channel distribution. A spatial weight function is constructed with the location of the vortex center as a reference. The spatial weight function is convolved with the virtual disturbance source to generate a non-uniform intensity distribution. The non-uniform intensity distribution is normalized to obtain the disturbance intensity distribution field.

[0032] A perturbation frequency sequence is constructed according to a preset sampling interval. The perturbation intensity distribution field and the perturbation frequency sequence are orthogonally decomposed to generate a multi-condition perturbation field. A Fourier transform is performed on the multi-condition perturbation field to construct a dynamic response matrix. The dynamic response matrix is ​​substituted into a numerical calculation model to solve the problem.

[0033] The transient mass transfer flux is extracted from the solution results, and the overall mass transfer flux distribution is obtained by integrating along the spatial direction. The overall mass transfer flux distribution is then subjected to eigenvalue decomposition to extract the flow modes.

[0034] The flow modes and the overall mass flux distribution are combined to construct a frequency response feature spectrum, and the perturbation frequency corresponding to the maximum mass flux value is extracted from the frequency response feature spectrum as the output of the optimization result.

[0035] Based on the variance of the residence time distribution in the tower section extracted from the corrected state field, tower sections with variances exceeding a preset variance threshold are labeled as flow-shifted tower sections. The spatial overlap between the flow-shifted tower sections and the mass transfer-dominant flow region is calculated, including:

[0036] The correction state field is divided into tower segment calculation units. The velocity field distribution is extracted from the tower segment calculation units. The motion trajectory of fluid particles is calculated based on the velocity field distribution. The motion trajectory of fluid particles is sampled according to time intervals to obtain a motion trajectory sequence.

[0037] The residence time of fluid particles is obtained by performing time statistics on the motion trajectory sequence. A probability distribution matrix is ​​constructed based on the residence time of fluid particles. Eigenvalue decomposition is performed on the probability distribution matrix to extract the variance component.

[0038] The variance components are combined according to the tower segment calculation units to construct a variance distribution field. The variance distribution field is compared with a preset variance threshold, and the tower segment calculation units that exceed the preset variance threshold are marked as flow offset tower segments.

[0039] A multi-level mesh structure is generated by performing mesh mapping on the flow offset tower section, and an overlapping mesh is obtained by performing projection operation on the multi-level mesh structure and the mass transfer dominant flow region.

[0040] The area normalization calculation is performed on the overlapping region, and the calculation result is output as the spatial overlap degree.

[0041] Based on the spatial overlap, the system selects to adjust either the distributor opening ratio or the feed pulsation frequency, and feeds the adjusted flow rate and concentration distribution data back to the digital twin model, including:

[0042] Calculate the difference between the spatial overlap and the preset reference threshold to generate a spatial overlap deviation value;

[0043] When the spatial overlap deviation is positive, the pressure gradient distribution in the mass transfer-dominant flow region is extracted, and the pressure gradient distribution is multiplied with the spatial overlap deviation to obtain the orifice ratio correction matrix. Based on the orifice ratio correction matrix, the orifice ratio of the distributor is adjusted to generate the adjusted orifice ratio of the distributor.

[0044] When the spatial overlap deviation is negative, the velocity field distribution in the mass transfer-dominant flow region is extracted, and the velocity field distribution is multiplied with the spatial overlap deviation to obtain the frequency correction matrix. The feed pulsation frequency is adjusted based on the frequency correction matrix to generate the adjusted feed pulsation frequency.

[0045] The flow rate distribution data and concentration distribution data were calculated based on the adjusted distributor opening ratio and the adjusted feed pulsation frequency, respectively.

[0046] The velocity distribution data and concentration distribution data are written into the digital twin model.

[0047] A second aspect of the present invention provides an electronic device comprising:

[0048] processor;

[0049] Memory used to store processor-executable instructions;

[0050] The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0051] A third aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0052] In this embodiment, a joint simulation framework combining a digital twin model and a numerical calculation model was established to achieve real-time status tracking and dynamic correction of the extraction and separation process. By employing a spatial comparison method between real-time operating data and the simulation baseline state, deviation areas can be accurately identified and local mass transfer coefficients can be adjusted accordingly. This significantly improves the simulation model's adaptability to actual operating conditions and its predictive accuracy, solving the technical challenge of traditional static models struggling to track changes in the actual equipment's operating state. Setting virtual disturbance sources in calibrated key areas and conducting multi-condition simulation optimization allows for the precise selection of the optimal disturbance frequency that maximizes mass transfer flux. Simultaneously, by combining residence time distribution variance analysis to identify flow deviation tower sections, and by calculating spatial overlap, a quantitative assessment of key influencing factors is achieved, enabling adaptive control of operating parameters. This effectively reduces energy consumption and operating costs, providing a complete technical solution for the intelligent operation and maintenance of the extraction and separation process. Attached Figure Description

[0053] Figure 1 is a flowchart illustrating the intelligent simulation optimization method for extraction and separation processes based on digital twins according to an embodiment of the present invention.

[0054] Figure 2 is a flowchart of flow field feature analysis and region calibration according to an embodiment of the present invention. Detailed Implementation

[0055] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, 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.

[0056] The technical solution of the present invention will be described in detail below with reference to specific embodiments. These specific embodiments can be combined with each other, and the same or similar concepts or processes may not be described again in some embodiments.

[0057] Figure 1 is a flowchart illustrating the intelligent simulation optimization method for extraction and separation processes based on digital twins according to an embodiment of the present invention. As shown in Figure 1, the method includes:

[0058] Establish digital twin models and numerical calculation models of the extraction and separation physical equipment to generate simulation baseline states;

[0059] Real-time operating data of the extraction and separation physical equipment is collected, and the real-time operating data is spatially compared with the simulation baseline state. Spatial regions where the deviation exceeds the set deviation threshold are identified. The local mass transfer coefficient of the spatial region is adjusted and the numerical calculation model is updated to obtain the corrected state field.

[0060] Based on the corrected state field, a coupling matrix of concentration gradient and velocity distribution is constructed. The coupling matrix is ​​decomposed into features, the dominant feature vector is extracted and mapped to the device space, and the dominant flow region of mass transfer is calibrated.

[0061] A virtual disturbance source is set up in the mass transfer-dominant flow region. Multi-condition simulation is carried out by changing the disturbance frequency of the virtual disturbance source, and the disturbance frequency that maximizes the mass transfer flux is selected.

[0062] Based on the correction state field, the residence time distribution variance of the tower section is extracted. The tower section with variance exceeding the preset variance threshold is labeled as the flow offset tower section. The spatial overlap between the flow offset tower section and the mass transfer dominant flow region is calculated.

[0063] Based on the spatial overlap, the distributor opening ratio or feed pulsation frequency is adjusted, and the adjusted flow rate distribution and concentration distribution data are fed back to the digital twin model.

[0064] In one optional implementation, a digital twin model and a numerical calculation model of the extraction and separation physical equipment are established to generate a simulation baseline state, including:

[0065] Obtain the geometric structure parameters and process operation parameters of the extraction and separation physical equipment, construct a three-dimensional digital space model as a digital twin model, and generate a spatial mesh based on the digital twin model;

[0066] Multiple virtual measurement points are constructed in the spatial grid, and the virtual measurement points are mapped to the corresponding physical measurement point locations in the extraction and separation physical equipment;

[0067] Pressure and concentration data are collected at the physical measurement points, and the pressure and concentration data are then substituted back into the virtual measurement points. The pressure field distribution and concentration field distribution are constructed using a grid interpolation method.

[0068] The flow field velocity distribution and flow field direction distribution are calculated based on the pressure field distribution and substituted into the mass transfer calculation equation. The concentration field distribution and the mass transfer calculation equation are combined to construct a numerical calculation model. The boundary constraints of the numerical calculation model are set based on the process operating parameters.

[0069] The numerical calculation model is iteratively calculated, and the iteration stops when the calculation residual is less than the convergence criterion. The calculation result is output as the simulation baseline state.

[0070] When establishing the digital twin model and numerical calculation model of the extraction and separation physical equipment, the geometric structural parameters of the equipment were obtained through on-site measuring instruments. Specifically, these parameters included a tower height of 12 meters, a tower diameter of 2.4 meters, a tray spacing of 0.5 meters, a tray opening ratio of 15%, and an opening diameter of 8 millimeters per tray. Simultaneously, process operating parameters were acquired, including a feed flow rate of 8 tons per hour, an extractant flow rate of 10 tons per hour, an operating temperature of 45 degrees Celsius, and an operating pressure of 0.2 MPa. Based on these parameters, a three-dimensional digital space model was constructed in the digital twin software platform according to the actual dimensions. This model fully reproduces the internal structure of the extraction tower, including the tray arrangement, the location of the inlet and outlet pipes, the distributor structure, and the spatial distribution characteristics inside the tower.

[0071] After constructing the 3D digital space model, the model space was meshed. An unstructured meshing method was used, with a mesh size of 50 mm in the main tower area, 10 mm in critical areas near the tower plates, and 5 mm at the opening locations. Approximately 4.5 million mesh cells were generated in the entire model space. The minimum mesh quality was 0.35, and the average mesh quality was 0.78, meeting the accuracy requirements for subsequent numerical calculations. The total number of mesh nodes reached approximately 1.8 million, with each mesh cell having clearly defined spatial coordinates and topological relationships.

[0072] Virtual measuring points are arranged in the generated spatial grid. A measurement section is set every 1 meter along the axial direction of the tower, with 8 measuring points evenly distributed radially in each section, for a total of 12 measurement sections along the tower height, resulting in 96 virtual measuring points. Each virtual measuring point has a unique spatial coordinate identifier; for example, the coordinates of measuring point No. 3 in the 5th layer are 0.9 meters radially from the tower center, 135 degrees circumferentially, and 5 meters axially. The virtual measuring points are mapped one-to-one with the physical measuring point positions of the pressure sensors and concentration analyzers actually installed in the extraction and separation physical equipment. The physical measuring points are installed on the tower wall via flange interfaces, with the measuring probe penetrating 0.2 meters into the tower body, ensuring that the spatial deviation between the measurement position and the virtual measuring point position does not exceed 0.05 meters.

[0073] Pressure and concentration data at physical measurement points are acquired in real time through a data acquisition system. Pressure data is acquired 10 times per second, and concentration data is acquired every 30 seconds via an online analyzer. Data is continuously collected for 2 hours under stable operating conditions. After data filtering and averaging, representative values ​​for each physical measurement point are obtained. For example, the pressure value at measurement point 1 in layer 3 is 0.185 MPa, and the target component concentration is 23.5%; the pressure value at measurement point 5 in layer 7 is 0.192 MPa, and the target component concentration is 45.8%. The acquired pressure and concentration data are then fed back to the corresponding virtual measurement points in the digital twin model via a data interface.

[0074] The pressure and concentration field distributions were constructed across the entire spatial grid using radial basis function grid interpolation. During interpolation, the value of each grid node was obtained by a weighted average of the values ​​from its six nearest virtual measuring points, with the weighting coefficient inversely proportional to the square of the distance. After interpolation, the pressure and concentration values ​​for all 1.8 million nodes in the spatial grid were obtained. The pressure field distribution showed a pressure of 0.205 MPa at the bottom of the column, gradually decreasing to 0.178 MPa at the top. The concentration field distribution showed that the target component concentration was 18% at the bottom feed position, reaching 62% at the extractant outlet at the top after the extraction process, and decreasing to 5% at the raffinate outlet at the bottom.

[0075] Based on the established pressure field distribution, the velocity distribution and direction distribution of the flow field are calculated using the correlation between the pressure gradient and the flow velocity. Within each grid cell, the pressure difference between adjacent nodes is calculated; dividing the pressure difference by the distance between nodes yields the pressure gradient, and the direction of the pressure gradient is the driving direction of the fluid flow. The continuous phase density is 850 kg / m³, which is considered in the physical property parameters. 3 The viscosity is 0.00, and the dispersed phase density is 920 kg / m³. 3 With a viscosity of 0.0015, the calculated peak velocity of the continuous phase at the tray opening was 0.35 m / s, and the average rising velocity of the dispersed phase droplets was 0.12 m / s. The flow field distribution showed that the continuous phase flowed downwards as a whole, while the dispersed phase moved upwards as a whole, forming obvious cross-flow characteristics in the tray region.

[0076] The obtained velocity and direction distributions of the flow field were substituted into the mass transfer calculation equations. These equations describe the concentration variation over time and space, involving convective and diffusion mass transfer terms. The convective mass transfer term is determined by the product of the flow field velocity and the concentration gradient, while the diffusion mass transfer term is determined by the product of the diffusion coefficient and the second-order concentration gradient. The diffusion coefficient of the target component in the continuous phase is 1.2 × 10⁻⁶. -9 m 2 / s, with a diffusion coefficient of 0.8×10 in the dispersed phase. -9 m 2 / s. The interfacial mass transfer coefficient was calculated using droplet diameter and contact time, with an average value of 0.0015 meters per second. The established concentration field distribution was used as the initial condition for the mass transfer calculation equation, and the flow field information was used as known coefficients in the equation, combining them to form a complete numerical calculation model.

[0077] Boundary constraints for the numerical calculation model are set based on process operating parameters. At the column inlet boundary, a continuous phase feed flow rate of 8 tons per hour corresponds to an inlet velocity of 0.25 m / s, and a feed concentration of 18%. At the extractant inlet boundary, a dispersed phase flow rate of 10 tons per hour corresponds to an inlet velocity of 0.18 m / s, and the target component concentration in the extractant is 0%. At the top of the column, the extractant phase outlet boundary is set with a pressure boundary condition of 0.178 MPa, allowing free flow of fluid. At the bottom of the column, the raffinate phase outlet boundary is set with a pressure boundary condition of 0.205 MPa. At the column wall boundary, a no-slip boundary condition is set, with zero fluid velocity, and an adiabatic boundary condition is also set, eliminating heat exchange.

[0078] The established numerical calculation model was solved iteratively. A non-steady-state calculation method with a time step of 0.01 seconds was used, with the pressure, velocity, and concentration fields solved in a coupled manner within each time step. During the iterative calculation, the changes in the calculation residuals of each physical quantity were monitored in real time. The pressure calculation residual was defined as the maximum absolute value of the difference between the pressure value of the current iteration step and the previous iteration step across all grid nodes. The velocity and concentration calculation residuals were defined in the same way. The convergence criteria were set as pressure residual less than 100 Pa, velocity residual less than 0.001 m / s, and concentration residual less than 0.1%. After 18,000 time steps (corresponding to 180 seconds of actual time), all calculation residuals met the convergence criteria requirements: the pressure residual decreased to 65 Pa, the velocity residual decreased to 0.0006 m / s, and the concentration residual decreased to 0.07%. The pressure field distribution, velocity field distribution, flow direction distribution, and concentration field distribution calculated at this time are used as the simulation baseline state output and stored in the state database of the digital twin system. This simulation baseline state fully describes the steady-state operating characteristics of the extraction and separation physical equipment under the current process operating parameters.

[0079] In one optional implementation, real-time running data is spatially compared with the simulation baseline state to identify spatial regions where the deviation exceeds a set deviation threshold. The local mass transfer coefficient of the spatial region is adjusted and the numerical calculation model is updated. The corrected state field is then obtained by solving for:

[0080] The deviation values ​​between the real-time running data and the simulation baseline state are calculated point by point on the spatial grid to generate a spatial deviation field. The spatial deviation field is then analyzed to extract connected regions whose deviation values ​​exceed a set deviation threshold and are marked as spatial regions.

[0081] Extract the spatial coordinate range of the spatial region, and obtain the local mass transfer coefficient of the spatial region from the solution domain of the numerical calculation model corresponding to the simulation baseline state based on the spatial coordinate range;

[0082] Calculate the difference in concentration distribution between real-time running data and simulation baseline state within the spatial region, integrate the concentration distribution difference along the flow direction, and calculate the mass transfer coefficient correction.

[0083] The mass transfer coefficient correction is superimposed with the local mass transfer coefficient to obtain the adjusted local mass transfer coefficient. The adjusted local mass transfer coefficient is then substituted into the mass transfer source term corresponding to the spatial region in the numerical calculation model, while keeping the mass transfer coefficient outside the spatial region unchanged, to obtain the updated numerical calculation model.

[0084] Apply boundary conditions corresponding to real-time running data to the updated numerical calculation model and perform numerical solutions. Terminate the iterative calculation when the solution residual is less than the set convergence threshold, and output the corrected state field.

[0085] First, the deviation between the real-time operating data and the simulation baseline state is calculated point-by-point on a spatial grid. Specifically, real-time operating data, including physical quantities such as pressure, concentration, and flow rate, is acquired from the extraction and separation physical equipment every 10 seconds. The acquired real-time operating data is mapped onto the same spatial grid as the simulation baseline state using a grid interpolation method, forming a real-time state field. The point differences between the real-time state field and the simulation baseline state are calculated to generate a spatial deviation field. For the concentration field, a deviation threshold of 5% relative error is set; for the pressure field, a deviation threshold of 3% relative error is set; and for the velocity field, a deviation threshold of 8% relative error is set.

[0086] Connectivity analysis was performed on the generated spatial deviation field, and a region growing algorithm was used to identify adjacent grid cells exceeding the deviation threshold. Starting with the grid cell with the largest deviation value, the grid cells in its six adjacent directions were checked. If the deviation values ​​of adjacent cells also exceeded the set threshold, they were merged into the current connected region. This process was repeated until no more grid cells could be merged, forming a complete connected region. The above operation was performed on all grid cells in the spatial deviation field that exceeded the threshold, ultimately obtaining multiple independent connected regions, each of which was labeled as a spatial region. A significant spatial region was identified between sections 5 and 7 of the extraction column, containing approximately 8500 grid cells, accounting for 1.9% of the total number of grid cells. The maximum deviation value occurred in section 6 near the column wall, with a concentration deviation reaching 7.8%.

[0087] For each marked spatial region, its spatial coordinate range is extracted, and the minimum and maximum coordinate values ​​of the region are recorded to determine its boundary in three-dimensional space. Based on this coordinate range, the corresponding solution domain is located from the numerical calculation model, and the local mass transfer coefficient within that region is extracted. For the liquid-liquid two-phase system in the extraction tower, the mass transfer coefficient is mainly controlled by the interfacial mass transfer coefficient, which is related to the phase interface area, the degree of fluid turbulence, and the interface renewal rate. Within the identified spatial regions, the initial local mass transfer coefficient is 0.0015 m / s, a value obtained based on theoretical calculations and initial calibration.

[0088] The concentration distribution difference between real-time operating data and the simulation baseline state within the spatial region is calculated. The concentration difference is obtained by subtracting the baseline concentration value from the real-time concentration value in each grid cell. The flow field direction distribution within the spatial region is acquired, and the main flow direction vector for each grid cell is determined. The spatial region is discretized into multiple flow tubes parallel to the main flow direction, each containing several grid cells arranged along the flow direction. The concentration difference is spatially integrated along the flow direction of each flow tube, and the integration result is divided by the flow tube length to obtain the average concentration change rate. The average concentration change rate is divided by the local fluid velocity to obtain the mass transfer coefficient correction. For the identified spatial regions from the 5th to the 7th segment of the extraction tower, the calculated mass transfer coefficient correction is -0.0004 m / s, indicating that the actual mass transfer efficiency is lower than the model prediction.

[0089] The adjusted local mass transfer coefficient is obtained by superimposing the mass transfer coefficient correction amount with the local mass transfer coefficient. In the current case, the adjusted local mass transfer coefficient is 0.0011 m / s, which is about 26.7% lower than the original value. The adjusted mass transfer coefficient is then updated in the mass transfer source term of the corresponding spatial region in the numerical calculation model. The mass transfer source term describes the mass transfer rate per unit volume, and its value is equal to the mass transfer coefficient multiplied by the interfacial area density and then multiplied by the interphase concentration difference. The mass transfer coefficient outside the spatial region is kept unchanged to ensure the locality and specificity of the model correction and to prevent over-correction from causing the overall model to become inaccurate. The boundary conditions are updated to the operating conditions corresponding to the real-time operating data, including adjusting the feed flow rate to 7.8 tons per hour, maintaining the extractant flow rate at 10 tons per hour, and increasing the operating temperature to 47 degrees Celsius.

[0090] The updated numerical model was iteratively solved, with convergence thresholds set as follows: pressure residual less than 80 Pa, velocity residual less than 0.0008 m / s, and concentration residual less than 0.08%. An explicit solution method with a time step of 0.01 seconds was used, and residual changes were dynamically monitored during the iteration. The calculation terminated when all residuals were less than the set thresholds, and the corrected state field was output. The correction calculation was performed for 12,000 time steps, corresponding to 120 seconds of actual time. Finally, the pressure residual decreased to 72 Pa, the velocity residual decreased to 0.0007 m / s, and the concentration residual decreased to 0.06%, meeting the convergence requirements.

[0091] In this embodiment, dynamic correction of the digital twin model is achieved through real-time data feedback and local model adjustment, effectively improving the model's prediction accuracy for the actual extraction and separation process. Spatial region identification and local mass transfer coefficient adjustment avoid the model instability issues that may arise from global parameter adjustments, thus improving computational efficiency. Based on the corrected state field, the impact of operating condition changes on the separation effect can be predicted more accurately, providing a reliable basis for process parameter optimization. This method can also promptly identify changes in equipment performance or abnormal conditions, supporting predictive maintenance, reducing equipment failure risks, extending equipment lifespan, and simultaneously improving product quality and production efficiency in the extraction and separation process.

[0092] Figure 2 shows the flowchart of flow field feature analysis and region calibration in this embodiment.

[0093] In one optional implementation, a coupling matrix of concentration gradient and velocity distribution is constructed based on the corrected state field. The coupling matrix is ​​then subjected to eigenvalue decomposition to extract the dominant eigenvector and map it to the device space. The calibrated mass transfer dominant flow region includes:

[0094] The concentration field in the correction state field is differentiated according to spatial coordinates to generate a concentration gradient field. The velocity field in the correction state field is orthogonally decomposed to generate a velocity distribution field. Multi-scale tensor product operation is performed on the concentration gradient field and the velocity distribution field to generate a multi-level coupling matrix.

[0095] The multi-level coupling matrix is ​​block diagonalized to generate a block matrix group. Singular value decomposition is performed on each block matrix in the block matrix group and regularization constraints are applied to obtain a local dominant feature vector group. The local dominant feature vector group is then combined with adaptive weighting to generate a global dominant feature vector.

[0096] A nonlinear spatial mapping function is constructed to transform the global dominant feature vector into a feature distribution field in the physical coordinate system of the device through the nonlinear spatial mapping function, and the gradient value and curvature value of each spatial point in the feature distribution field are calculated.

[0097] The characteristic distribution field is spatially partitioned, and the connected regions where both the gradient value and curvature value are greater than the partition mean are designated as the mass transfer-dominant flow regions.

[0098] When constructing the coupling matrix of concentration gradient and velocity distribution based on the calibration state field, it is necessary to process the concentration field data in the calibration state field. In the digital twin model of extraction and separation, the concentration field data is stored as scalar values ​​on three-dimensional spatial grid nodes, covering the entire extraction tower space. The concentration field data is differentiated according to spatial coordinates to calculate the concentration gradient at each grid node. In specific implementation, the central difference scheme is used to calculate the concentration gradient. For internal nodes, the concentration gradient components along the three coordinate axes are equal to the concentration difference between adjacent nodes divided by the distance between the two nodes; for boundary nodes, one-sided difference is used to calculate the gradient. The concentration gradient detected at the 10th layer position on the central axis of the extraction tower is 0.08 radially, 0.02 circumferentially, and 0.35 axially, indicating that the axial concentration change is the most significant.

[0099] The velocity field in the corrected state field contains information about the direction and speed of fluid motion, requiring orthogonal decomposition to extract its main features. Orthogonal decomposition of the velocity field employs a decomposition method based on the characteristic flow field, representing the velocity field as a linear combination of a series of orthogonal basis functions. The velocity autocorrelation matrix of the entire computational domain is calculated, and its eigenvectors are extracted as orthogonal basis functions. Typically, the eigenvectors corresponding to the first 10 eigenvalues ​​are selected, which can reconstruct more than 85% of the energy of the original velocity field. The resulting velocity distribution field contains the projection coefficients of the velocity amplitude onto each orthogonal component. In the region near the 6th tray of the extraction tower, the dominant component of the velocity distribution field exhibits obvious swirling characteristics, with the swirling center located 0.8 meters from the tower center, and the circumferential velocity component reaching 0.15 m / s.

[0100] A multi-scale tensor product operation is performed on the concentration gradient field and the velocity distribution field to generate a multi-level coupling matrix. The tensor product operation correlates the concentration gradient and velocity components at different spatial scales. In implementation, wavelet transform is used to decompose the concentration gradient field and velocity distribution field into different frequency components, and the tensor product of the two fields is calculated at each frequency scale. For the extraction tower model, four analysis scales are set, corresponding to feature sizes of 0.05 m, 0.2 m, 0.5 m, and 1 m. At each scale, the tensor product result forms a matrix, where the matrix elements represent the coupling strength between the concentration gradient and velocity components at that scale. Combining the results from each scale yields the multi-level coupling matrix, whose size is the product of the dimensions of the concentration gradient field and the velocity distribution field.

[0101] The multi-level coupling matrix is ​​block diagonalized, dividing it into several sub-blocks based on spatial relationships. Each sub-block corresponds to a local region within the extraction tower, such as a section of the tower or the area near a tray. Block diagonalization employs a graph partitioning algorithm, dividing the computational domain into non-overlapping sub-regions. Each grid node within a sub-region corresponds to a sub-block in the coupling matrix. For a 12-meter-high extraction tower, it is divided into 24 sub-regions along its height, each approximately 0.5 meters high. Singular value decomposition (SVD) is performed on each sub-block matrix to extract its main feature patterns. SVD decomposes the sub-block matrix into the product of three matrices, where the elements of the diagonal matrix are singular values, representing the importance of the corresponding feature patterns.

[0102] Regularization constraints are applied to the singular value decomposition results to suppress the influence of noise and minor features. Regularization employs a truncated singular value decomposition method, retaining only eigenvectors with singular values ​​greater than 20% of the average singular value, while setting the remaining eigenvectors to zero. This method typically retains 2 to 4 principal eigenvectors for each sub-block matrix, forming a locally dominant eigenvector group. In the 8th layer region of the extraction column, the first eigenvector of the locally dominant eigenvector group captures the strong coupling mode between radial flow and axial concentration gradient, with the maximum value of the eigenvector element occurring at a distance of 0.3 meters from the column wall.

[0103] The global dominant feature vector is generated by adaptively weighting and combining the local dominant feature vectors of each sub-region. During weighting, the weight coefficients are proportional to the singular values ​​corresponding to the local feature vectors and inversely proportional to the spatial distance between sub-regions, ensuring a continuous feature distribution in space. For feature vectors of adjacent sub-regions, a smooth transition is performed at the boundary to avoid discontinuous jumps in feature distribution. The global dominant feature vector has the same dimension as the original spatial grid nodes and contains the main features coupling the concentration gradient and velocity distribution throughout the entire computational domain.

[0104] A nonlinear spatial mapping function is constructed to map vectors in the feature space back to the device's physical coordinate system. The spatial mapping is implemented using a radial basis function network, selecting 200 key nodes as the basis function centers. The basis functions are Gaussian functions, and the width parameter is automatically adjusted based on the average distance between nodes. The mapping function transforms the globally dominant feature vectors into a feature distribution field in the device's physical coordinate system. The value at each spatial point in the feature distribution field represents the strength of the coupling between the concentration gradient and velocity at that location. The gradient and curvature values ​​of the feature distribution field at each spatial point are calculated. The gradient value is calculated using central difference, and the curvature value is obtained using the Laplacian operator.

[0105] Spatial partitioning analysis was performed on the characteristic distribution field, calculating the average gradient and curvature values ​​across the entire computational domain. The mean gradient value was 0.42, and the mean curvature value was 1.85. Spatial points in the characteristic distribution field where both gradient and curvature values ​​were greater than their respective mean values ​​were marked as candidate points, and connected component analysis was performed on these candidate points. A region growing algorithm was used to group spatially adjacent candidate points into the same connected region, ultimately resulting in several independent connected regions. These connected regions were designated as mass transfer-dominant flow regions, representing locations within the extraction column where the concentration gradient and fluid motion are strongly coupled, and are the most active regions in the mass transfer process. In the actual extraction column model, three main mass transfer-dominant flow regions were identified, located in the 3rd-4th, 7th-9th, and 11th-12th layers, respectively. The 7th-9th layer region had the largest area, accounting for approximately 35% of the column's cross-sectional area.

[0106] In this embodiment, the mass transfer-dominant flow region was accurately calibrated by constructing a coupling matrix between the concentration gradient and velocity distribution. Based on multi-scale analysis and nonlinear mapping techniques, the mass transfer characteristics under complex flow conditions can be effectively captured, overcoming the limitations of traditional methods in handling spatiotemporal multi-scale coupling problems. The accurately calibrated mass transfer-dominant flow region provides a clear target for subsequent optimization design, enabling control measures to be applied specifically to key areas and improving control efficiency. It can also adapt to different operating conditions and equipment geometry changes, exhibiting good versatility and providing strong support for digital twin simulation and intelligent optimization of extraction and separation processes, significantly improving separation efficiency and product quality.

[0107] In one optional implementation, a virtual disturbance source is set in the mass transfer-dominant flow region. Multi-condition simulations are performed by changing the disturbance frequency of the virtual disturbance source, and the disturbance frequencies that maximize the mass transfer flux are selected, including:

[0108] A spatial discrete grid is constructed in the mass transfer-dominant flow region, and the nodes of the spatial discrete grid are set as virtual disturbance sources. The velocity field distribution and concentration field distribution of the mass transfer-dominant flow region are extracted, and the mass transfer channel distribution is obtained through gradient calculation.

[0109] The location of the vortex center is identified based on the mass transfer channel distribution. A spatial weight function is constructed with the location of the vortex center as a reference. The spatial weight function is convolved with the virtual disturbance source to generate a non-uniform intensity distribution. The non-uniform intensity distribution is normalized to obtain the disturbance intensity distribution field.

[0110] A perturbation frequency sequence is constructed according to a preset sampling interval. The perturbation intensity distribution field and the perturbation frequency sequence are orthogonally decomposed to generate a multi-condition perturbation field. A Fourier transform is performed on the multi-condition perturbation field to construct a dynamic response matrix. The dynamic response matrix is ​​substituted into a numerical calculation model to solve the problem.

[0111] The transient mass transfer flux is extracted from the solution results, and the overall mass transfer flux distribution is obtained by integrating along the spatial direction. The overall mass transfer flux distribution is then subjected to eigenvalue decomposition to extract the flow modes.

[0112] The flow modes and the overall mass flux distribution are combined to construct a frequency response feature spectrum, and the perturbation frequency corresponding to the maximum mass flux value is extracted from the frequency response feature spectrum as the output of the optimization result.

[0113] When setting up virtual disturbance sources in the mass transfer-dominant flow region, this region needs to be refined. A spatial discrete grid is constructed for the calibrated mass transfer-dominant flow region. Unstructured grid generation technology is used to refine the grid in the active mass transfer region to improve computational accuracy. For the mass transfer-dominant flow region of the 7th-9th layers of the extraction tower, approximately 42,000 grid cells are constructed. The minimum grid cell size is 5 mm, located in areas of drastic flow field gradient changes; the maximum grid cell size is 25 mm, located in relatively uniform flow field regions. The total number of nodes in the discrete grid is approximately 35,000, and each grid node is set as a virtual disturbance source for subsequent disturbance energy injection. Velocity and concentration field distribution data for this mass transfer-dominant flow region are extracted from the digital twin model. The velocity field data contains three directional components, and the concentration field data represents the mole fraction distribution of the target component. Gradient calculations are performed on the velocity and concentration field data to calculate the velocity gradient tensor and concentration gradient vector at each grid node. The velocity gradient tensor characterizes the local flow field deformation characteristics, and the concentration gradient vector indicates the direction of the fastest concentration change. The inner product operation of the velocity gradient tensor and the concentration gradient vector results in a scalar field that characterizes the degree of coordination between velocity and concentration changes. This scalar field is defined as the mass transfer channel distribution.

[0114] Based on the mass transfer channel distribution, vortex structures in the flow field were identified, and the Q-criterion was used to identify the location of the vortex center. The Q-criterion is calculated based on the invariant of the velocity gradient tensor. When the Q value is greater than zero and reaches a local maximum, the corresponding location is identified as the vortex center. A significant vortex center was identified 50 mm above the 8th tray of the extraction tower, with coordinates of 0.75 m radially and 135 degrees circumferentially from the tower center. A spatial weighting function was constructed based on this vortex center, using a Gaussian distribution. The weight is maximum at the vortex center (1.0) and decreases exponentially with increasing distance. The standard deviation parameter of the Gaussian distribution was set to 0.3 m to ensure that the weighting function covers the vortex influence range. The spatial weighting function was convolved with the virtual disturbance source grid, assigning different intensity values ​​to each virtual disturbance source node to generate a non-uniform intensity distribution. Virtual disturbance sources near the vortex center received higher intensity values, with a maximum of 0.95; virtual disturbance sources farther from the vortex center received lower intensity values, with a minimum of 0.15. The non-uniform intensity distribution is normalized so that the sum of the intensities is 1.0, thus obtaining the disturbance intensity distribution field.

[0115] A perturbation frequency sequence is constructed according to a preset sampling interval. Considering the characteristics of the extraction and separation process and the equipment response characteristics, the frequency range is set from 0.1 Hz to 2.0 Hz, with a sampling interval of 0.1 Hz, generating a total of 20 frequency points. The perturbation intensity distribution field and the perturbation frequency sequence are orthogonally decomposed to generate a multi-condition perturbation field. For each frequency point, the periodic change of the perturbation intensity over time is calculated to form a time-domain perturbation signal. The time-domain perturbation signal is applied to the virtual perturbation source node to form a dynamic perturbation field. The multi-condition perturbation field contains 20 sets of dynamic perturbation fields with different frequencies, each set corresponding to a perturbation frequency. A Fourier transform is performed on the multi-condition perturbation field to convert the time-domain signal into a frequency-domain representation, constructing a dynamic response matrix. The rows of the dynamic response matrix correspond to spatial grid nodes, the columns correspond to frequency components, and the matrix elements represent the response amplitude and phase information of a specific node at a specific frequency. The dynamic response matrix is ​​substituted as a boundary condition into the numerical calculation model, and the solution is performed separately for each perturbation frequency condition.

[0116] Transient mass transfer flux data were extracted from the solution results. Mass transfer flux is defined as the amount of mass transferred per unit area per unit time. For a liquid-liquid two-phase system, the mass transfer flux equals the mass transfer coefficient multiplied by the interfacial area density and then multiplied by the interphase concentration difference. Multiple monitoring sections were set up along the flow direction within the computational domain, with 72 monitoring points evenly distributed on each section. The transient mass transfer flux value at each monitoring point was extracted, and the mass transfer flux on the monitoring section was integrated along the spatial direction to obtain the total mass transfer flux of the section. The total mass transfer flux of all monitoring sections was averaged over time to obtain the overall mass transfer flux distribution. For the 0.5 Hz disturbance frequency condition, the total mass transfer flux at the 8th tray reached its maximum value of 8.2 × 10⁻⁶. -3 This represents an improvement of 18.3% compared to the unperturbed baseline state.

[0117] Eigenvalue decomposition was performed on the overall mass flux distribution, and the dominant flow modes were extracted using an intrinsically orthogonal decomposition method. Mass flux distribution data at different times were collected, a spatiotemporal correlation matrix was constructed, and its eigenvalues ​​and eigenvectors were solved. The eigenvectors represent the spatial modes of the flow, and the eigenvalues ​​represent the energy contribution of each mode. Typically, the top five modes with the largest energy contributions are selected, as these modes can reconstruct more than 85% of the energy of the original signal. For the 0.8 Hz perturbation frequency condition, the first mode exhibits axial mass transfer enhancement, contributing 56%; the second mode exhibits radial mass transfer enhancement, contributing 22%. The extracted flow modes were combined with the overall mass flux distribution to construct a frequency response characteristic spectrum. The frequency response characteristic spectrum is a two-dimensional matrix, where rows represent spatial locations, columns represent perturbation frequencies, and matrix elements represent the mass flux values ​​under corresponding conditions. By analyzing the frequency response characteristic spectrum, the maximum mass flux and its corresponding perturbation frequency were identified. In all test conditions, the overall mass flux reached its maximum value of 8.5 × 10⁻⁶ under the 0.7 Hz perturbation frequency condition. -3 This represents a 23.2% improvement over the unperturbed baseline. 0.7 Hz was determined as the optimal perturbation frequency for the output results.

[0118] In this embodiment, by setting a virtual disturbance source in the mass transfer-dominant flow region and combining multi-condition simulation with frequency optimization analysis, the mass transfer efficiency of the extraction and separation process is improved. Digital twin technology allows for virtual testing and optimization without interfering with actual production, avoiding the high cost and low efficiency of traditional trial-and-error methods. The non-uniform intensity distribution design based on mass transfer channel distribution and vortex center identification ensures that disturbance energy is precisely applied to key flow regions, maximizing energy utilization efficiency. This not only improves the mass transfer efficiency of the extraction and separation process but also reduces energy consumption and the amount of separating agent used.

[0119] In one optional implementation, the variance of the residence time distribution of the tower section is extracted based on the corrected state field, and the tower section with a variance exceeding a preset variance threshold is labeled as a flow-off tower section. The spatial overlap between the flow-off tower section and the mass transfer-dominant flow region is calculated, including:

[0120] The correction state field is divided into tower segment calculation units. The velocity field distribution is extracted from the tower segment calculation units. The motion trajectory of fluid particles is calculated based on the velocity field distribution. The motion trajectory of fluid particles is sampled according to time intervals to obtain a motion trajectory sequence.

[0121] The residence time of fluid particles is obtained by performing time statistics on the motion trajectory sequence. A probability distribution matrix is ​​constructed based on the residence time of fluid particles. Eigenvalue decomposition is performed on the probability distribution matrix to extract the variance component.

[0122] The variance components are combined according to the tower segment calculation units to construct a variance distribution field. The variance distribution field is compared with a preset variance threshold, and the tower segment calculation units that exceed the preset variance threshold are marked as flow offset tower segments.

[0123] A multi-level mesh structure is generated by performing mesh mapping on the flow offset tower section, and an overlapping mesh is obtained by performing projection operation on the multi-level mesh structure and the mass transfer dominant flow region.

[0124] The area normalization calculation is performed on the overlapping region, and the calculation result is output as the spatial overlap degree.

[0125] When extracting the variance of residence time distribution in the extraction column based on the corrected state field, the corrected state field needs to be reasonably divided. The corrected state field of the extraction column is divided into column segment calculation units according to the tray position. Each column segment calculation unit contains one tray and a 25 cm space above and below it. For an extraction column with a height of 12 meters and 24 trays, 24 column segment calculation units are formed, each unit approximately 50 cm high. Velocity field distribution data is extracted from the column segment calculation units. The velocity field contains components in three directions: radial velocity, circumferential velocity, and axial velocity. The velocity field data is stored using a three-dimensional array, with the array dimensions consistent with the grid division of the calculation units. The average axial velocity in the 8th column segment calculation unit of the extraction column is 0.15 m / s, the maximum radial velocity is 0.08 m / s, and the maximum circumferential velocity is 0.05 m / s.

[0126] The motion trajectory of fluid particles was calculated based on the extracted velocity field distribution. 1000 virtual fluid particles were uniformly distributed in each column segment's computational unit, with the initial position of each particle marked using three-dimensional coordinates. The fourth-order Runge-Kutta integral method was used to solve the particle motion equations, with a time step of 0.01 seconds. For each particle, its local velocity was calculated by interpolation based on its current position, and then the particle position was updated based on the velocity. If a particle moved to the boundary of the computational unit, its departure time was recorded and marked as out of bounds. The trajectory was continuously calculated for all particles that had not yet out of bounds until the simulation time reached the preset maximum time of 300 seconds or all particles out of bounds. In the 8th column segment's computational unit, the average residence time of the fluid particles was 32.5 seconds, the shortest residence time was 15.8 seconds, and the longest residence time was 89.3 seconds. The motion trajectory of the fluid particles was sampled at 0.5-second intervals, and the particle position at each sampling time was recorded to form a motion trajectory sequence. The motion trajectory sequence of each particle contains multiple three-dimensional coordinate points, which are arranged in chronological order to fully describe the motion history of the particle within the computing unit.

[0127] The generated motion trajectory sequence is statistically analyzed to calculate the residence time of each fluid particle. Residence time is defined as the time interval from when a particle enters the computational unit to when it leaves. For particles that have not left the computational unit by the end of the simulation, their residence time is counted as the total simulation time. A probability distribution matrix is ​​constructed based on the residence time data of all fluid particles. The residence time range is divided into 30 equally spaced time intervals, and the number of particles falling into each time interval is counted. The probability value is obtained by dividing the number of particles by the total number of particles. The rows of the probability distribution matrix represent the computational units of the tower segment, the columns represent the time intervals, and the matrix element values ​​are the probability values ​​under the corresponding conditions. Eigenvalue decomposition is performed on the probability distribution matrix, and principal component analysis is used to extract the main distribution features. In the eigenvalue decomposition results, the eigenvalues ​​represent the variance contribution of each distribution pattern, and the eigenvectors represent the shape of the distribution patterns. The variance component is extracted as an indicator to measure the uniformity of the residence time distribution. The variance component is equal to the weighted sum of the eigenvalues, and the weights are the second moments of the corresponding eigenvectors.

[0128] The extracted variance components are combined according to the column segment calculation units to construct a variance distribution field. The variance distribution field is a one-dimensional array, with each element corresponding to the variance value of a column segment calculation unit. The variance values ​​of the 24 column segment calculation units in the extraction column are distributed as follows: the variance values ​​of columns 1-5 are between 0.25 and 0.35, indicating relatively uniform flow; the variance values ​​of columns 6-9 suddenly increase to 0.58-0.72, indicating significant flow non-uniformity; the variance values ​​of columns 10-16 drop back to 0.3-0.4; the variance values ​​of columns 17-20 rise again to 0.45-0.55; and the variance values ​​of columns 21-24 decrease to 0.2-0.3. A preset variance threshold of 0.45 is set, and column segment calculation units with variance values ​​exceeding this threshold are marked as flow-off columns. Based on the variance comparison results, sections 6-9 and 17-20 were marked as flow-displaced sections. The fluid residence time distribution in these sections was uneven, and there were local short circuits or dead zones.

[0129] Mesh mapping was performed on the marked flow displacement tower segment to generate a multi-level mesh structure. An octree mesh subdivision method was used to progressively subdivide the computational mesh of the flow displacement tower segment. The original mesh cell size was 25 mm; after the first level of subdivision, the mesh cell size was 12.5 mm; and after the second level of subdivision, the mesh cell size was 6.25 mm. The mesh subdivision level was adaptively adjusted according to the flow field characteristics, with higher levels of subdivision used in regions of drastic flow field changes. The generated multi-level mesh structure contained approximately 280,000 mesh cells, accurately describing the geometric features and flow field distribution of the flow displacement tower segment. Projection operations were performed on the multi-level mesh structure and the mass transfer-dominant flow region to identify the spatial overlap. The projection operation used a spatial intersection algorithm to determine whether each mesh cell belonged to both the flow displacement tower segment and the mass transfer-dominant flow region. If it belonged to both, the mesh cell was marked as part of the overlapping region. The overlapping region consisted of a series of mesh cells that simultaneously exhibited flow inhomogeneity and active mass transfer.

[0130] Area normalization calculations were performed on the identified overlapping areas, and the total volume of the overlapping areas was divided by the total volume of the flow displacement tower segments to obtain the spatial overlap. For the flow displacement areas of tower segments 6-9, the overlapping area volume was 0.36 m³. 3 The total volume of the flow offset tower section is 2.26m³. 3 The calculated spatial overlap is 0.159. For the flow displacement region of tower segments 17-20, the overlapping region volume is 0.58 m³. 3 The total volume of the flow offset tower section is 2.26m³. 3 The calculated spatial overlap was 0.257. The spatial overlap value represents the degree of spatial correlation between flow inhomogeneity and mass transfer activity; a larger value indicates a more significant spatial overlap, and control measures need to consider both flow uniformity and mass transfer efficiency. The calculated spatial overlap result was stored in the state database of the digital twin system as a basis for subsequent optimization decisions.

[0131] In this embodiment, by extracting the variance of the residence time distribution in the tower section and calculating its spatial overlap with the dominant mass transfer flow region, accurate identification and quantitative assessment of flow nonuniformity during the extraction and separation process are achieved. Based on digital twin technology and fluid trajectory analysis, detailed information on the internal flow state of the equipment can be obtained without interfering with actual production, overcoming the limitations of traditional methods in observing the internal flow field. Through variance analysis and spatial overlap calculation, a quantitative correlation between flow nonuniformity and mass transfer efficiency is established, significantly improving production efficiency and product quality while reducing energy and material consumption.

[0132] In one optional implementation, selecting to adjust the distributor orifice ratio or the feed pulsation frequency based on the spatial overlap, and feeding the adjusted flow rate distribution and concentration distribution data back to the digital twin model includes:

[0133] Calculate the difference between the spatial overlap and the preset reference threshold to generate a spatial overlap deviation value;

[0134] When the spatial overlap deviation is positive, the pressure gradient distribution in the mass transfer-dominant flow region is extracted, and the pressure gradient distribution is multiplied with the spatial overlap deviation to obtain the orifice ratio correction matrix. Based on the orifice ratio correction matrix, the orifice ratio of the distributor is adjusted to generate the adjusted orifice ratio of the distributor.

[0135] When the spatial overlap deviation is negative, the velocity field distribution in the mass transfer-dominant flow region is extracted, and the velocity field distribution is multiplied with the spatial overlap deviation to obtain the frequency correction matrix. The feed pulsation frequency is adjusted based on the frequency correction matrix to generate the adjusted feed pulsation frequency.

[0136] The flow rate distribution data and concentration distribution data were calculated based on the adjusted distributor opening ratio and the adjusted feed pulsation frequency, respectively.

[0137] The velocity distribution data and concentration distribution data are written into the digital twin model.

[0138] When selecting to adjust the distributor opening ratio or feed pulsation frequency based on spatial overlap, it is necessary to calculate the difference between the spatial overlap and a preset reference threshold. The preset reference threshold is determined based on the design specifications and process requirements of the extraction and separation equipment. For an industrial extraction tower with a height of 12 meters and a diameter of 2.4 meters, the preset reference threshold is set to 0.2. The previously calculated spatial overlap value is obtained; for the flow offset region of sections 17-20, the spatial overlap is 0.257. The difference between the spatial overlap and the preset reference threshold is calculated, i.e., 0.257 minus 0.2, resulting in a spatial overlap deviation value of 0.057, which is positive. The sign of the spatial overlap deviation value indicates the adjustment direction: a positive value indicates a high degree of overlap between the flow offset and the mass transfer dominant region, requiring adjustment of the distributor opening ratio; a negative value indicates a low degree of overlap, requiring adjustment of the feed pulsation frequency.

[0139] When the spatial overlap deviation is positive, the pressure gradient distribution of the mass transfer-dominant flow region is extracted from the digital twin model. The pressure gradient distribution is stored as a three-dimensional array, representing the rate of change of pressure in three spatial directions. In the mass transfer-dominant flow region of tower sections 17-20, the average pressure gradient is 1200 Pa / m, and the maximum is 3500 Pa / m, occurring in the region near the tower wall of section 18. The extracted pressure gradient distribution is multiplied by the spatial overlap deviation, with each component of the pressure gradient multiplied by the deviation value of 0.057, to obtain the opening ratio correction matrix. The opening ratio correction matrix has the same spatial dimension as the pressure gradient distribution, and the matrix element values ​​represent the amount of opening ratio adjustment required at each location. The opening ratio adjustment is larger in the high-pressure gradient region and smaller in the low-pressure gradient region.

[0140] The distributor's aperture ratio is adjusted based on an aperture ratio correction matrix. The distributor in the extraction column is located on each tray, with an original designed aperture ratio of 15% and an aperture diameter of 8 mm. The aperture ratio is adjusted by changing the aperture diameter or the number of apertures. For the four trays of sections 17-20, the aperture ratio adjustment for each region is calculated. Each tray plane is divided into five annular regions with radii from the inside to the outside: 0-0.4 m, 0.4-0.8 m, 0.8-1.2 m, 1.2-1.6 m, and 1.6-2.0 m. For the 17th tray, the opening ratio adjustments for the five zones from the inside out are: increase by 2.1%, increase by 1.5%, decrease by 0.8%, decrease by 1.2%, and increase by 0.5%, respectively; for the 18th tray, the adjustments are: increase by 1.8%, decrease by 0.9%, decrease by 1.5%, increase by 2.3%, and increase by 1.2%, respectively; for the 19th tray, the adjustments are: decrease by 0.6%, decrease by 1.1%, increase by 1.9%, increase by 1.4%, and decrease by 0.8%, respectively; and for the 20th tray, the adjustments are: increase by 0.7%, increase by 1.3%, increase by 0.4%, decrease by 1.7%, and decrease by 1.0%, respectively.

[0141] The calculated adjustment amount of the opening ratio is superimposed on the original opening ratio to obtain the adjusted distributor opening ratio. For areas where the opening ratio increases, the opening diameter or the number of openings is appropriately increased; for areas where the opening ratio decreases, the opening diameter or the number of openings is appropriately decreased. After adjustment, the opening ratios of the 18th tray distributor, from the inside out, are 16.8%, 14.1%, 13.5%, 17.3%, and 16.2% for the five regions, with an average opening ratio of 15.58%. By adjusting the opening distribution, the fluid distribution pattern on the tray is changed, allowing more fluid to flow through the mass transfer-dominant region and reducing flow drift.

[0142] When the spatial overlap deviation is negative, the velocity field distribution of the mass transfer-dominant flow region is extracted from the digital twin model. Taking the flow offset region of tower sections 6-9 as an example, its spatial overlap is 0.159, and compared with the preset reference threshold of 0.2, the deviation is -0.041, which is negative. The velocity field distribution of this region is extracted, containing velocity components in three directions. The average axial velocity at the 7th tower plate is 0.13 m / s, the average radial velocity is 0.06 m / s, and the average circumferential velocity is 0.04 m / s. The velocity field distribution is multiplied by the spatial overlap deviation, and each velocity component is multiplied by the absolute value of the deviation, 0.041, to obtain the frequency correction matrix. The frequency correction matrix contains three components, corresponding to the frequency adjustment amounts in the axial, radial, and circumferential directions, respectively.

[0143] The feed pulsation frequency is adjusted based on a frequency correction matrix. The original feed condition is a steady-state feed with no pulsation frequency. The optimal pulsation frequency is calculated based on the frequency correction matrix, and the frequency adjustment is calculated for the axial, radial, and circumferential directions. The frequency adjustment is directly proportional to the local velocity and inversely proportional to the characteristic length. For the 6th-9th tower section, the calculated axial pulsation frequency is 0.65 Hz, the radial pulsation frequency is 0.38 Hz, and the circumferential pulsation frequency is 0.25 Hz. Considering the frequencies in the three directions, the axial pulsation frequency is selected as the main adjustment parameter, and the feed pulsation frequency is set to 0.65 Hz. Feed pulsation is achieved by installing a pulsation generator before the feed pump. This device can periodically change the feed flow rate at a set frequency to generate pulsating flow.

[0144] Based on the adjusted distributor orifice ratio and the adjusted feed pulsation frequency, velocity distribution data and concentration distribution data were calculated, respectively. For sections 17-20 of the column with adjusted distributor orifice ratio, a new velocity distribution was calculated using numerical simulation. Simulation results show that the adjusted distributor orifice ratio resulted in a more uniform fluid distribution on the trays, reducing short-circuit flow and dead zones. A comparison of the velocity distribution before and after adjustment showed that the standard deviation of the velocity decreased from 0.038 m / s to 0.025 m / s, improving flow uniformity by 34%. For sections 6-9 of the column with the feed pulsation frequency set, simulation results showed that pulsating feed enhanced fluid mixing and mass transfer. The pulsation frequency of 0.65 Hz is close to the characteristic frequency of the fluid in this region, creating a resonance effect that significantly enhanced interface renewal and mass transfer rates. Concentration distribution data showed that the mass transfer efficiency of the target component in this region increased by 21%, and the steeper concentration gradient indicated a reduction in mass transfer resistance.

[0145] The calculated velocity and concentration distribution data are written into a digital twin model. The velocity distribution data contains velocity components in three directions, formatted as a three-dimensional array corresponding to the grid division of the computational domain. The concentration distribution data is a scalar field, storing the target component concentration value at each grid node. An incremental update strategy is used during data writing, updating only the changed regions to reduce data transfer volume and storage requirements. The updated digital twin model reflects the extraction and separation process status after adjustments and can be used for subsequent monitoring and optimization. Simultaneously, the adjustment parameters (distributor orifice distribution or feed pulsation frequency) are recorded in the process parameter database as a basis for actual equipment adjustments.

[0146] In this embodiment, two different mechanisms—distributor orifice ratio adjustment and feed pulsation frequency adjustment—provide effective solutions for different types of flow deviation problems. Based on virtual simulation and feedback optimization using digital twin technology, the orifice ratio adjustment mechanism alters the spatial distribution of the fluid, enhancing flow uniformity; the feed pulsation frequency adjustment mechanism introduces time-varying disturbances, strengthening mixing and mass transfer effects, achieving synergistic optimization of the flow field and mass transfer field, and overcoming the limitations of single control methods. This not only improves extraction and separation efficiency but also reduces energy consumption and separator dosage, and enhances product quality stability.

[0147] A second aspect of the present invention provides an electronic device, comprising:

[0148] processor;

[0149] Memory used to store processor-executable instructions;

[0150] The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0151] A third aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0152] This invention can be a method, apparatus, system, and / or computer program product. The computer program product may include a computer-readable storage medium having computer-readable program instructions loaded thereon for performing various aspects of the invention.

[0153] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A digital twin-based intelligent simulation optimization method for extraction and separation processes, characterized in that, include: Establish digital twin models and numerical calculation models of the extraction and separation physical equipment to generate simulation baseline states; Real-time operating data of the extraction and separation physical equipment is collected, and the real-time operating data is spatially compared with the simulation baseline state. Spatial regions with deviations exceeding a set deviation threshold are identified. The local mass transfer coefficient of the spatial region is adjusted, and the numerical calculation model is updated to obtain the corrected state field. This includes: calculating the deviation value of the real-time operating data and the simulation baseline state point by point on the spatial grid to generate a spatial deviation field; performing connectivity analysis on the spatial deviation field to extract connected regions with deviation values ​​exceeding the set deviation threshold and marking them as spatial regions; extracting the spatial coordinate range of the spatial region; obtaining the local mass transfer coefficient of the spatial region from the solution domain of the numerical calculation model corresponding to the simulation baseline state based on the spatial coordinate range; calculating the concentration distribution difference between the real-time operating data and the simulation baseline state in the spatial region; integrating the concentration distribution difference along the flow direction to calculate the mass transfer coefficient correction; superimposing the mass transfer coefficient correction with the local mass transfer coefficient to obtain the adjusted local mass transfer coefficient; and substituting the adjusted local mass transfer coefficient into the spatial region in the numerical calculation model. The corresponding mass transfer source term is used to keep the mass transfer coefficient outside the spatial region unchanged, resulting in an updated numerical calculation model. Boundary conditions corresponding to real-time running data are applied to the updated numerical calculation model and numerical solutions are performed. The iterative calculation is terminated when the solution residual is less than a set convergence threshold, and a corrected state field is output. A coupling matrix of concentration gradient and velocity distribution is constructed based on the corrected state field. The coupling matrix is ​​subjected to eigenvalue decomposition to extract the dominant feature vector and map it to the equipment space to calibrate the dominant mass transfer flow region. A virtual disturbance source is set in the dominant mass transfer flow region. Multi-condition simulation is performed by changing the disturbance frequency of the virtual disturbance source to select the disturbance frequency that maximizes the mass transfer flux. The variance of the residence time distribution of the tower section is extracted based on the corrected state field. The tower section with variance exceeding a preset variance threshold is calibrated as a flow offset tower section, and the spatial overlap between the flow offset tower section and the dominant mass transfer flow region is calculated. Based on the spatial overlap, the distributor opening ratio adjustment or feed pulsation frequency adjustment is selected, and the adjusted velocity distribution and concentration distribution data are fed back to the digital twin model.

2. The method according to claim 1, characterized in that, The process of establishing a digital twin model and a numerical calculation model of the extraction and separation physical equipment to generate a simulation baseline state includes: acquiring the geometric structural parameters and process operating parameters of the extraction and separation physical equipment; constructing a three-dimensional digital space model as the digital twin model; generating a spatial mesh based on the digital twin model; constructing multiple virtual measuring points in the spatial mesh and mapping the virtual measuring points to the corresponding physical measuring point locations in the extraction and separation physical equipment; collecting pressure and concentration data at the physical measuring point locations; substituting the pressure and concentration data back into the virtual measuring points; constructing pressure field distribution and concentration field distribution using a mesh interpolation method; calculating the flow field velocity distribution and flow field direction distribution based on the pressure field distribution and substituting them into the mass transfer calculation equation; combining the concentration field distribution with the mass transfer calculation equation to construct a numerical calculation model; setting boundary constraints for the numerical calculation model based on the process operating parameters; performing iterative calculations on the numerical calculation model; stopping the iteration when the calculation residual is less than the convergence criterion; and outputting the calculation results as the simulation baseline state.

3. The method according to claim 1, characterized in that, Based on the corrected state field, a coupling matrix of concentration gradient and velocity distribution is constructed. The coupling matrix is ​​then subjected to eigenvalue decomposition to extract dominant feature vectors and map them to the device space. The calibration of the mass transfer dominant flow region includes: differentiating the concentration field in the corrected state field according to spatial coordinates to generate a concentration gradient field; orthogonally decomposing the velocity field in the corrected state field to generate a velocity distribution field; performing multi-scale tensor product operations on the concentration gradient field and velocity distribution field to generate a multi-level coupling matrix; performing block diagonalization on the multi-level coupling matrix to generate a block matrix group; performing singular value decomposition on each block matrix in the block matrix group and applying regularization constraints to obtain a local dominant feature vector group; and generating a global dominant feature vector group through adaptive weighted combination of the local dominant feature vector group; constructing a nonlinear spatial mapping function; transforming the global dominant feature vector through the nonlinear spatial mapping function to generate a feature distribution field in the device physical coordinate system; calculating the gradient value and curvature value of each spatial point in the feature distribution field; spatially partitioning the feature distribution field; and calibrating the connected regions where both the gradient value and curvature value are greater than the partition mean as the mass transfer dominant flow region.

4. The method according to claim 1, characterized in that, Virtual disturbance sources are set up in the mass transfer-dominant flow region. Multi-condition simulations are conducted by changing the disturbance frequency of these virtual sources. The disturbance frequency that maximizes the mass transfer flux is selected by: constructing a spatial discrete grid in the mass transfer-dominant flow region, setting the nodes of the spatial discrete grid as virtual disturbance sources, extracting the velocity and concentration field distributions of the mass transfer-dominant flow region, and obtaining the mass transfer channel distribution through gradient calculation; identifying the vortex center position based on the mass transfer channel distribution, constructing a spatial weight function based on the vortex center position, performing a convolution operation between the spatial weight function and the virtual disturbance source to generate a non-uniform intensity distribution, and normalizing the non-uniform intensity distribution. The disturbance intensity distribution field is obtained; a disturbance frequency sequence is constructed according to a preset sampling interval; the disturbance intensity distribution field and the disturbance frequency sequence are orthogonally decomposed to generate a multi-condition disturbance field; a Fourier transform is performed on the multi-condition disturbance field to construct a dynamic response matrix; the dynamic response matrix is ​​substituted into a numerical calculation model for solution; transient mass transfer flux is extracted from the solution result; the overall mass transfer flux distribution is obtained by integration along the spatial direction; eigenvalue decomposition is performed on the overall mass transfer flux distribution to extract flow modes; the flow modes and the overall mass transfer flux distribution are combined to construct a frequency response feature spectrum; the disturbance frequency corresponding to the maximum mass transfer flux is extracted from the frequency response feature spectrum as the optimization result output.

5. The method according to claim 1, characterized in that, Based on the variance of the residence time distribution of the tower segment extracted from the corrected state field, tower segments with variances exceeding a preset variance threshold are labeled as flow-off tower segments. The spatial overlap between the flow-off tower segments and the mass transfer-dominant flow region is calculated as follows: the corrected state field is divided into tower segment calculation units; velocity field distribution is extracted from the tower segment calculation units; fluid particle trajectories are calculated based on the velocity field distribution; the fluid particle trajectories are sampled at time intervals to obtain a trajectory sequence; the residence time of fluid particles is obtained by performing time statistics on the trajectory sequence; a probability distribution matrix is ​​constructed based on the residence time of fluid particles; eigenvalue decomposition is performed on the probability distribution matrix to extract variance components; the variance components are combined according to the tower segment calculation units to construct a variance distribution field; the variance distribution field is compared with a preset variance threshold, and tower segment calculation units exceeding the preset variance threshold are labeled as flow-off tower segments; mesh mapping is performed on the flow-off tower segments to generate a multi-level mesh structure; projection operation is performed on the multi-level mesh structure and the mass transfer-dominant flow region to obtain an overlapping mesh; area normalization calculation is performed on the overlapping region, and the calculation result is output as the spatial overlap.

6. The method according to claim 1, characterized in that, Based on the spatial overlap, the system selects to adjust either the distributor opening ratio or the feed pulsation frequency, and feeds the adjusted velocity and concentration distribution data back to the digital twin model. This includes: calculating the difference between the spatial overlap and a preset reference threshold to generate a spatial overlap deviation value; when the spatial overlap deviation is positive, extracting the pressure gradient distribution in the mass transfer-dominant flow region, multiplying the pressure gradient distribution with the spatial overlap deviation value to obtain an opening ratio correction matrix, and adjusting the distributor opening ratio based on this matrix to generate the adjusted distributor opening ratio; when the spatial overlap deviation is negative, extracting the velocity field distribution in the mass transfer-dominant flow region, multiplying the velocity field distribution with the spatial overlap deviation value to obtain a frequency correction matrix, and adjusting the feed pulsation frequency based on this matrix to generate the adjusted feed pulsation frequency; calculating the velocity and concentration distribution data based on the adjusted distributor opening ratio and the adjusted feed pulsation frequency; and writing the velocity and concentration distribution data into the digital twin model.

7. An electronic device, characterized in that, include: processor; A memory for storing processor-executable instructions; wherein the processor is configured to invoke the instructions stored in the memory to perform the method according to any one of claims 1 to 6.

8. A computer-readable storage medium having computer program instructions stored thereon, characterized in that, When the computer program instructions are executed by the processor, they implement the method described in any one of claims 1 to 6.

Images