Ethylene carbonate storage stability simulation method based on digital twinning

Abstract

The present application relates to the technical field of chemical material storage stability simulation, in particular to a carbonic acid ethylene storage stability simulation method based on digital twinning, aiming to solve the problem that traditional evaluation cannot deeply couple dynamic factors such as moisture penetration, temperature and humidity change and product accumulation, cannot accurately restore the moisture gradient balance and product accumulation law in the container, and is difficult to predict the evolution of key quality indicators and failure boundaries, by constructing a digital twinning body basic model containing molecular structure and hydrolysis kinetics parameters, embedding dynamic moisture adsorption and reaction kinetics modules, calculating parameters such as moisture diffusion and penetration rate, integrating temperature and humidity sensor data, forming feedback by coupling heat and mass transfer and reaction kinetics equations, outputting stability attenuation risk index, predicting acid value, moisture content, color evolution trend and identifying failure threshold, and iteratively correcting model parameters combined with actual detection data to improve simulation accuracy.

Description

Technical Field

[0001] This invention relates to the field of chemical material storage stability simulation technology, specifically to a method for simulating the storage stability of ethylene carbonate based on digital twins. Background Technology

[0002] Chemical material storage stability simulation is an important technology, specifically applied to the industrial storage stability assessment and prediction of ethylene carbonate. Its core lies in improving the accuracy of storage stability simulation through multi-field dynamic coupling, adapting to the core needs of chemical production for the safe and quality control of ethylene carbonate storage. Ethylene carbonate has an inherent characteristic of easy hydrolysis. During storage, external moisture diffuses and permeates through the walls of sealed containers. Dynamic changes in ambient temperature and relative humidity directly affect the hydrolysis reaction rate, while the accumulation of hydrolysis products in turn interferes with moisture diffusion and the reaction process. These factors interact to form a complex multi-physics coupling effect. Because traditional storage stability assessments do not deeply couple and provide real-time feedback on these dynamically related factors, they cannot accurately reproduce the dynamic balance of the moisture gradient within the container and the accumulation pattern of products. Consequently, it is difficult to accurately predict the evolution trend and stability failure boundary of key quality indicators such as acid value, moisture content, and color, failing to provide reliable support for the safe storage and quality assurance of ethylene carbonate. To solve this technical problem, we provide a digital twin-based method for ethylene carbonate storage stability simulation. Summary of the Invention

[0003] The purpose of this invention is to provide a method for simulating the storage stability of ethylene carbonate based on digital twins, so as to solve the problems mentioned in the background art.

[0004] To achieve the above objectives, this invention provides a method for simulating the storage stability of ethylene carbonate based on digital twins, comprising the following steps:

[0005] S1. Construct a digital twin model of the molecular structure and hydrolysis reaction kinetics of ethylene carbonate;

[0006] S2. In the digital twin basic model, a dynamic water adsorption and reaction kinetics module designed based on the easy hydrolysis characteristics of ethylene carbonate is embedded. The dynamic water adsorption and reaction kinetics module calculates the water molecule diffusion and permeation rate, hydrolysis reaction rate and product concentration distribution in real time.

[0007] S3. Based on the water molecule diffusion and permeation rate, hydrolysis reaction rate and product concentration distribution, a dynamic multiphysics field real-time feedback module is designed. The dynamic multiphysics field real-time feedback module integrates temperature and relative humidity sensor data streams, and adjusts the water molecule diffusion and permeation rate and hydrolysis reaction rate in real time by coupling heat transfer and mass transfer and chemical reaction kinetic equations to simulate the dynamic balance of the water gradient and product accumulation effect in a closed container, and outputs a stability decay risk index.

[0008] S4. Using the stability decay risk index, a time series is set to predict the evolution trend of key quality indicators of ethylene carbonate under different storage conditions. The key quality indicators include acid value, moisture content and color, and the stability failure threshold is identified.

[0009] S5. By applying machine learning algorithms to compare the detection data of the actual stored samples, the predicted values ​​of the digital twin basic model, and the stability failure threshold, the parameters of the digital twin basic model are iteratively corrected, and the visualized simulation results are output.

[0010] Preferably, when constructing the digital twin basic model of ethylene carbonate molecular structure and hydrolysis reaction kinetic parameters, quantum chemical calculations are used to determine the molecular conformation and charge distribution of ethylene carbonate, and molecular dynamics simulations are combined to obtain the transition state energy and reaction energy barrier in the hydrolysis reaction path. The transition state energy and reaction energy barrier are then input into a spatial discretization model based on a finite element mesh. The spatial discretization model uses adaptive mesh refinement technology to capture the influence of intermolecular forces on the hydrolysis reaction rate, thus forming the basic framework of the digital twin model.

[0011] Preferably, the dynamic water adsorption and reaction kinetics module designed based on the easy hydrolysis characteristics of ethylene carbonate specifically introduces the permeability tensor of porous media to describe the structural characteristics of the inner wall material of the sealed container, combines the diffusion barrier function of water molecules in the polymer material to establish a nonlinear mapping relationship between the water molecule adsorption isotherm and the local concentration gradient of ethylene carbonate, and dynamically updates the activation energy of the hydrolysis reaction based on the transition state theory, thus coupling the material interface characteristics with the chemical reaction.

[0012] Preferably, when the dynamic water adsorption and reaction kinetics module calculates the water molecule diffusion and permeation rate in real time, it uses the porous medium permeability tensor and the preset water molecule diffusion barrier function to solve the unsteady diffusion equation and obtain the instantaneous flux of water molecules in the container wall. Specifically, when calculating the hydrolysis reaction rate, the amount of reactant consumed per unit time is derived using a microscopic reaction probability model based on the local concentration gradient of ethylene carbonate and the updated activation energy of the hydrolysis reaction. When calculating the product concentration distribution, the product transport equation is constructed based on the reactant consumption and diffusion flux through the principle of mass conservation and the three-dimensional spatial concentration field is solved.

[0013] Preferably, when designing the dynamic multiphysics real-time feedback module, the water molecule diffusion and permeation rate is used as the mass transfer field input variable, the hydrolysis reaction rate is used as the chemical reaction field input variable, and the product concentration distribution is used as the convection diffusion field input variable. The temperature and relative humidity sensor data streams are integrated to drive the heat transfer field boundary conditions. The temporal synchronization of the four types of physical fields is coordinated by the field coupling controller, and an adaptive grid re-division mechanism with the product accumulation rate as the feedback signal is established.

[0014] Preferably, when constructing the heat transfer, mass transfer and chemical reaction kinetic equations, a water diffusion equation is established based on the principle of mass conservation, an unsteady-state heat conduction equation is established based on the principle of energy conservation, and a hydrolysis reaction rate equation is established based on the transition state theory. By introducing the inhibition coefficient of product concentration distribution on the reaction rate, the three are combined into a nonlinear partial differential equation system with cross terms, and the operator splitting method is used to decompose the equation system into independently solvable heat transfer sub-equations, proton transfer equations and reaction kinetic sub-equations.

[0015] Preferably, when adjusting the water molecule diffusion and hydrolysis reaction rate in real time, the water flux output by the proton transfer equation is used as the boundary condition update value of the reaction kinetic sub-equation by the solution obtained by the operator splitting method. At the same time, the product inhibition coefficient generated by the reaction kinetic sub-equation is fed back to the proton transfer equation to correct the diffusion coefficient, forming a two-way real-time feedback link. The solution step size of each sub-equation is dynamically adjusted by the field coupling controller to match the sensor data sampling frequency.

[0016] Preferably, when simulating the dynamic equilibrium of the moisture gradient and the product accumulation effect in a closed container, the three-dimensional moisture distribution field inside the container is reconstructed based on the output data of the two-way real-time feedback link using a spatial interpolation algorithm. The product accumulation area at different time points is predicted by combining the product transport equation. The extreme values ​​of the moisture gradient and the product accumulation concentration are mapped to a stability decay risk index by a normalization method. The stability decay risk index includes a moisture penetration risk component and a product corrosion risk component.

[0017] Preferably, when predicting the evolution trend of key quality indicators, the moisture penetration risk component in the stability decay risk index is associated with the change rate of acid value and moisture content, and the product corrosion risk component is associated with the change rate of color. The time series extrapolation method is used to generate the indicator evolution curve. When identifying the stability failure threshold, the acid value exceeding the standard critical point, the moisture saturation concentration point and the color change threshold are located based on the abrupt inflection point of the stability decay risk index. The acid value exceeding the standard critical point, the moisture saturation concentration point and the color change threshold are integrated into a composite failure criterion.

[0018] Preferably, when iteratively correcting the parameters of the digital twin basic model, the detection data sequence of the acid value, moisture content and color of the actual stored sample is time-series aligned with the predicted values ​​of the digital twin basic model at the corresponding time nodes. The prediction deviation weight is calculated through the composite failure criterion, and the correlation between the deviation weight and the water molecule diffusion and permeation rate coefficient and the hydrolysis reaction rate constant is analyzed using a long short-term memory neural network. The parameter correction amount is generated and backpropagated to the digital twin basic model to achieve closed-loop optimization.

[0019] Compared with the prior art, the beneficial effects of the present invention are:

[0020] This invention constructs a digital twin model that integrates quantum chemical calculations and molecular dynamics simulations to replicate the molecular characteristics and underlying laws of ethylene carbonate hydrolysis. It embeds a dynamic moisture adsorption and reaction kinetics module to achieve deep coupling between material interface properties and chemical reactions. This allows for real-time calculation of water molecule diffusion and permeation rates, hydrolysis reaction rates, and product concentration distribution. A dynamic multiphysics real-time feedback module integrates temperature and humidity sensor data streams. By coupling heat and mass transfer equations with reaction kinetics equations, a bidirectional feedback link is formed, accurately simulating the dynamic equilibrium of moisture gradients and product accumulation effects in a closed container. A quantified stability decay risk index is output. Based on this index, the evolution trends of key quality indicators such as acid value, moisture content, and color can be predicted. The failure threshold is identified by the inflection point of the stability decay risk index, and the model parameters are iteratively corrected using an LSTM neural network to continuously improve simulation accuracy. Attached Figure Description

[0021] Figure 1 This is a flowchart illustrating the overall workflow of the present invention. Detailed Implementation

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

[0023] Please see Figure 1 As shown, this invention provides a method for simulating the storage stability of ethylene carbonate based on digital twins, comprising the following steps:

[0024] S1. Construct a digital twin model of the molecular structure and hydrolysis reaction kinetics of ethylene carbonate;

[0025] S2. In the basic model of digital twin, a dynamic water adsorption and reaction kinetics module designed based on the easy hydrolysis characteristics of ethylene carbonate is embedded. The dynamic water adsorption and reaction kinetics module calculates the water molecule diffusion and permeation rate, hydrolysis reaction rate and product concentration distribution in real time.

[0026] S3. A dynamic multiphysics real-time feedback module is designed based on the water molecule diffusion and permeation rate, hydrolysis reaction rate and product concentration distribution. The dynamic multiphysics real-time feedback module integrates temperature and relative humidity sensor data streams. By coupling heat transfer and mass transfer with chemical reaction kinetic equations, the water molecule diffusion and permeation rate and hydrolysis reaction rate are adjusted in real time to simulate the dynamic balance of the water gradient and product accumulation effect in a closed container, and outputs the stability decay risk index.

[0027] S4. Use the stability decay risk index to set a time series and predict the evolution trend of key quality indicators of ethylene carbonate under different storage conditions. Key quality indicators include acid value, moisture content and color, and identify the stability failure threshold.

[0028] S5. By applying machine learning algorithms to compare the detection data of the actual stored samples, the predicted values ​​of the digital twin basic model, and the stability failure threshold, the parameters of the digital twin basic model are iteratively corrected, and the visualized simulation results are output.

[0029] In constructing the digital twin basic model of ethylene carbonate molecular structure and hydrolysis reaction kinetic parameters, quantum chemical calculations were used to determine the molecular conformation and charge distribution of ethylene carbonate. Molecular dynamics simulations were combined to obtain the transition state energy and reaction barrier in the hydrolysis reaction path. The transition state energy and reaction barrier were then input into a spatial discretization model based on a finite element mesh. The spatial discretization model used adaptive mesh refinement technology to capture the influence of intermolecular forces on the hydrolysis reaction rate, thus forming the basic framework of the digital twin model.

[0030] The dynamic water adsorption and reaction kinetics module designed based on the easy hydrolysis characteristics of ethylene carbonate specifically introduces the permeability tensor of porous media to describe the structural characteristics of the inner wall material of the sealed container, combines the diffusion barrier function of water molecules in the polymer material to establish a nonlinear mapping relationship between the water molecule adsorption isotherm and the local concentration gradient of ethylene carbonate, and dynamically updates the activation energy of the hydrolysis reaction based on the transition state theory, thus coupling the material interface characteristics with the chemical reaction.

[0031] When the dynamic water adsorption and reaction kinetics module calculates the water molecule diffusion and permeation rate in real time, it uses the porous medium permeability tensor and the preset water molecule diffusion barrier function to solve the unsteady diffusion equation and obtain the instantaneous flux of water molecules in the container wall. Specifically, when calculating the hydrolysis reaction rate, it uses the local concentration gradient of ethylene carbonate and the updated activation energy of the hydrolysis reaction to derive the reactant consumption per unit time using a microscopic reaction probability model. When calculating the product concentration distribution, it constructs the product transport equation based on the reactant consumption and diffusion flux through the principle of mass conservation and solves the three-dimensional spatial concentration field.

[0032] When designing the dynamic multiphysics real-time feedback module, the water molecule diffusion and permeation rate is used as the input variable of the mass transfer field, the hydrolysis reaction rate is used as the input variable of the chemical reaction field, and the product concentration distribution is used as the input variable of the convection and diffusion field. The data streams of temperature and relative humidity sensors are integrated to drive the boundary conditions of the heat transfer field. The temporal synchronization of the four types of physical fields is coordinated by the field coupling controller, and an adaptive mesh re-division mechanism with the product accumulation rate as the feedback signal is established.

[0033] When constructing the heat transfer, mass transfer, and chemical reaction kinetic equations, a water diffusion equation is established based on the principle of mass conservation, an unsteady-state heat conduction equation is established based on the principle of energy conservation, and a hydrolysis reaction rate equation is established based on the transition state theory. By introducing the inhibition coefficient of product concentration distribution on the reaction rate, the three are combined into a nonlinear partial differential equation system with cross terms. The operator splitting method is used to decompose the equation system into independently solvable heat transfer sub-equations, proton transfer equations, and reaction kinetic sub-equations.

[0034] When adjusting the water molecule diffusion and hydrolysis reaction rate in real time, the operator splitting method is used to solve the problem. The water flux output by the proton transfer equation is used as the boundary condition update value of the reaction kinetic sub-equation. At the same time, the product inhibition coefficient generated by the reaction kinetic sub-equation is fed back to the proton transfer equation to correct the diffusion coefficient, forming a two-way real-time feedback link. The solution step size of each sub-equation is dynamically adjusted by the field coupling controller to match the sensor data sampling frequency.

[0035] When simulating the dynamic equilibrium of moisture gradient and product accumulation effect in a closed container, the spatial interpolation algorithm is used to reconstruct the three-dimensional moisture distribution field inside the container based on the output data of the two-way real-time feedback link. Combined with the product transport equation, the product accumulation area at different time points is predicted. The extreme value of moisture gradient and product accumulation concentration are mapped to the stability decay risk index by the normalization method. The stability decay risk index includes the moisture penetration risk component and the product corrosion risk component.

[0036] When predicting the evolution trend of key quality indicators, the moisture penetration risk component in the stability decay risk index is correlated with the change rate of acid value and moisture content, and the product corrosion risk component is correlated with the change rate of color. The time series extrapolation method is used to generate the indicator evolution curve. When identifying the stability failure threshold, the abrupt change inflection point of the stability decay risk index is used to locate the acid value exceeding the standard critical point, the moisture saturation concentration point, and the color change threshold. The acid value exceeding the standard critical point, the moisture saturation concentration point, and the color change threshold are integrated into a composite failure criterion.

[0037] When iteratively correcting the parameters of the digital twin basic model, the detection data sequences of acid value, moisture content and color of the actual stored samples are time-series aligned with the predicted values ​​of the digital twin basic model at the corresponding time nodes. The prediction deviation weight is calculated through composite failure criteria, and the correlation between the deviation weight and the water molecule diffusion and permeation rate coefficient and the hydrolysis reaction rate constant is analyzed using a long short-term memory neural network. The parameter correction amount is generated and backpropagated to the digital twin basic model to achieve closed-loop optimization.

[0038] It is necessary to further explain that, as the core foundation of the entire digital twin-based method for simulating the storage stability of ethylene carbonate, the construction of the digital twin basic model is a crucial prerequisite for subsequent dynamic simulation of hydrolysis reactions, multiphysics coupling calculations, and storage stability predictions. Its core lies in the progressive integration of quantum chemical calculations, molecular dynamics simulations, and finite element spatial discretization. This integration reveals the intrinsic physicochemical properties of reduced ethylene carbonate and the underlying laws governing hydrolysis reactions across multiple scales, from electronic structure and molecular reaction characteristics to macroscopic spatial distribution. The specific implementation method is as follows:

[0039] First, quantum chemical calculations are performed to determine the molecular conformation and charge distribution of ethylene carbonate. Quantum chemical calculations, based on fundamental principles of quantum mechanics, are a computational method that precisely describes the intrinsic properties of molecules, such as geometric configuration, electron cloud distribution, and energy state, at the electronic structure level by solving the Schrödinger equation. This is a core method for obtaining precise physicochemical parameters of small organic molecules and can reveal the reactive sites and structural stability of ethylene carbonate from the most fundamental level. The specific calculation process is as follows:

[0040] Based on the five-membered ring chemical structure of ethylene carbonate (C3H4O3), the atomic connection between the two ester oxygen groups and one epoxy group was clarified. An initial three-dimensional configuration was constructed using molecular visualization modeling software. Empirical values ​​for bond lengths and bond angles of organic molecules were loaded as the starting point for calculations to ensure the initial configuration conforms to fundamental chemical principles. For small oxygen-containing five-membered ring organic molecules like ethylene carbonate, a balance between computational accuracy and efficiency was struck. The hybrid functional B3LYP paired with a 6-311G(d,p) basis set was selected. This combination has been proven to have excellent accuracy and stability in the geometric optimization and electronic structure calculations of organic ester molecules. A polarized continuum model (PCM) was introduced to simulate the bulk liquid environment of ethylene carbonate, correcting the discrepancy between gas-phase calculations and actual liquid storage environments, making the calculation results more consistent with real-world storage scenarios. Minimizing the energy of the molecular system was the convergence objective. Force convergence thresholds were set to 1×10^-5 Hartree / Bohr, displacement convergence thresholds to 1×10^-4 Bohr, and the maximum number of iterations to 100. Each iteration recalculates the electronic structure and atomic forces of the molecule, continuously adjusting atomic coordinates until the molecular configuration reaches the potential energy minimum, obtaining the optimal molecular conformation of ethylene carbonate that is most stable in the actual environment. Based on the optimal molecular conformation, the natural atomic charge of each atom in the molecule is calculated using the natural bond orbital analysis (NBO) method, obtaining a complete intramolecular charge distribution. The charge density of hydrolysis reaction active sites such as carbonyl carbon and epoxy groups is accurately located, clarifying the core sites of nucleophilic reactions. At the same time, the energy level difference between the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO) is calculated, providing basic data at the electronic structure level for subsequent hydrolysis reaction pathway analysis. The fifth step completes the configuration rationality verification. Frequency analysis is performed on the optimized molecular conformation to confirm that there are no imaginary frequencies, proving that the configuration is a stable ground state configuration at the potential energy minimum, avoiding misjudging the transition state configuration as a stable molecular conformation, and ensuring that the calculated molecular conformation and charge distribution are chemically reasonable.

[0041] After completing the quantum chemical calculations of the intrinsic properties of ethylene carbonate molecules, molecular dynamics simulations were further used to obtain the transition state energy and reaction barrier in the hydrolysis reaction pathway. The transition state energy refers to the system energy corresponding to the highest potential intermediate configuration encountered during the conversion of reactants into products in the ethylene carbonate hydrolysis reaction pathway; it is the energy peak that the reaction must overcome. The reaction barrier, also called the activation barrier, is the difference between the transition state energy and the ground state energy of the reactants. Its magnitude directly determines the ease with which the hydrolysis reaction occurs; the lower the barrier, the easier the hydrolysis reaction occurs, and the faster the reaction rate. Molecular dynamics simulations, based on mechanical principles, simulate the motion, collisions, and interactions of molecules over a certain timescale, capturing the configurational changes and energy evolution during the reaction process. This method can accurately reconstruct the complete hydrolysis reaction pathway and energy changes. The specific implementation process is as follows:

[0042] Based on the optimal molecular conformation of ethylene carbonate obtained through quantum chemical calculations, a reaction system conforming to practical storage scenarios was constructed. The core of the system contains one ethylene carbonate reactant molecule and one water molecule as a nucleophile. An appropriate amount of bulk ethylene carbonate molecules are added to construct a solvation environment, simulating the bulk effect of liquid storage. The total number of atoms in the system is controlled within the range of 100-200 to balance computational accuracy and efficiency. Based on the core mechanism of the ethylene carbonate hydrolysis reaction, the water molecule nucleophilically attacks the carbonyl carbon, initiating a ring-opening reaction, ultimately producing ethylene glycol and carbon dioxide, accompanied by acidic products. The initial state (the reactant system of ethylene carbonate and water molecules) and the final state (the hydrolysis product system) of the predefined reaction are used. A series of intermediate configurations along the reaction path are generated using linear interpolation as initial guesses for subsequent precise transition state searches. The CI-NEB (Climbing Elastic Band) method is used to locate the transition state. The intermediate configurations along the reaction path are used as image points. The distance between image points is constrained by spring force, and the image point with the highest energy is allowed to climb towards the peak potential energy surface to accurately lock the transition state configuration along the reaction path. The B3LYP / 6 method, consistent with quantum chemical calculations, is used in the search process. Using the -311G(d,p) method and basis set, the system energy and atomic forces at each pixel are calculated in each iteration. The convergence threshold is set to a maximum atomic force of less than 0.05 eV / Å. After convergence, the unique and reasonable transition state configuration of the hydrolysis reaction is obtained, along with the corresponding transition state energy. The ground state energy of the initial reactants is calculated first, and then the core reaction energy barrier of the ethylene carbonate hydrolysis reaction is obtained by subtracting the ground state energy of the reactants from the transition state energy. At the same time, frequency analysis is performed on the obtained transition state configuration to confirm that it has one and only one imaginary frequency, and the vibrational direction corresponding to the imaginary frequency is... By aligning the bond breaking and bond formation directions of the hydrolysis reaction with perfect matching, and then calculating using intrinsic reaction coordinates (IRC), path integrals were performed from the transition state towards both reactants and products. This confirmed that the transition state can smoothly connect reactants and products, thoroughly verifying the chemical rationality of the transition state configuration and reaction pathway. Through molecular dynamics simulations with controlled variables, the transition state energy and reaction barrier under different temperatures and moisture contents were calculated, establishing a correspondence library of temperature-moisture content-reaction barrier relationships. This provides a full-scenario kinetic basis for stability simulations under different storage conditions.

[0043] After obtaining the core kinetic parameters of the hydrolysis reaction, the transition state energy and reaction energy barrier are input into a spatial discretization model based on a finite element mesh. Adaptive mesh refinement technology is used to capture the influence of intermolecular forces on the hydrolysis reaction rate, ultimately forming a complete digital twin basic model framework. The spatial discretization model divides the continuous three-dimensional space of the sealed storage container into several non-overlapping discrete units using a finite element mesh, transforming the solution of the continuous physical field into a model for numerical solutions on each discrete unit. This is the core carrier for coupling the microscopic hydrolysis reaction with the macroscopic storage environment. The adaptive mesh refinement technology is a numerical processing technique that automatically adjusts the mesh density based on the gradient changes of physical quantities within the solution region. The mesh is refined in regions with drastic changes in physical quantities and coarsened in regions with gradual changes, significantly improving numerical computation efficiency while maintaining computational accuracy. The specific implementation process is as follows:

[0044] Based on the three-dimensional dimensions of the sealed container used for actual storage of ethylene carbonate, a 1:1 scale geometric model was constructed. Two core computational domains were clearly defined: the container wall (polymer material such as HDPE) and the internal ethylene carbonate liquid phase storage region. Topology repair was performed on the geometric model, removing minor chamfers, redundant structures, and other details that do not affect computational accuracy, ensuring the feasibility and quality of subsequent mesh generation. Initial tetrahedral unstructured mesh generation was performed on the geometric model, setting the initial global mesh size to 1 mm. At the interface between the container wall and the liquid phase, the initial mesh size was refined to 0.5 mm to ensure computational accuracy at the interface. The transition state energy and reaction barrier obtained from previous calculations were used as core kinetic parameters and embedded into the properties of each mesh element in the liquid phase region. Based on the Arrhenius equation, a constitutive relation for the reaction kinetics of each mesh element was established, converting the reaction barrier into a hydrolysis reaction rate constant at different temperatures. This enabled each mesh element to have numerical computational capabilities for microscopic hydrolysis reactions. Two core encryption triggering conditions were defined: one is the intermolecular force gradient threshold, where intermolecular forces, including van der Waals forces and hydrogen bonds, affect the moisture content. The core factors for the diffusion and collision probability of ethylene carbonate molecules are: first, mesh refinement is triggered when the intermolecular force gradient of adjacent mesh cells exceeds a preset threshold of 5%; second, the reaction rate gradient threshold is triggered when the hydrolysis reaction rate gradient of adjacent mesh cells exceeds a preset threshold of 10%. Both thresholds are calibrated through multiple sets of numerical experiments, achieving an optimal balance between computational accuracy and efficiency. First, molecular dynamics simulations are used to calculate the distribution of intermolecular forces between ethylene carbonate molecules and between ethylene carbonate and water molecules at different concentrations and temperatures, mapping this distribution to each element of the finite element mesh. Then, all mesh elements are traversed, detecting the force gradient and reaction rate gradient of each element and its adjacent elements. For regions meeting the refinement trigger conditions, mesh subdivision is performed, dividing the original mesh element into eight sub-mesh units half the original size. Simultaneously, parameters such as the reaction energy barrier and reaction rate constant of the original element are accurately mapped to the sub-mesh. For regions where the gradient is below the threshold and the physical quantity changes smoothly, mesh coarsening is performed, merging adjacent low-gradient meshes to reduce unnecessary computation. After each mesh update, the mesh quality is re-verified to ensure that the Jacobian determinant is greater than 0.7. To avoid computational errors caused by mesh distortion, this dynamic mesh density adjustment can accurately capture regions with drastic changes in intermolecular forces, such as the interface between the container wall and the liquid phase, and regions with large water concentration gradients, and their impact on the hydrolysis reaction rate. Changes in intermolecular forces in these regions directly alter the diffusion efficiency and collision probability of water molecules, thus affecting the hydrolysis reaction rate. The finer mesh can accurately reproduce the reaction rate changes in these regions, avoiding numerical dissipation caused by coarse meshes, ensuring the accuracy of simulation results, and realizing the intrinsic molecular properties such as molecular conformation and charge distribution obtained from quantum chemical calculations. The reaction kinetic parameters, such as transition state energy and reaction barrier, obtained from kinetic simulations are deeply coupled with a finite element spatial discretization model equipped with adaptive mesh refinement technology. This establishes a multi-scale mapping relationship between "microscopic molecular properties - mesoscopic reaction kinetics - macroscopic spatial distribution," forming a complete digital twin basic model framework. This framework can both reconstruct the essential laws of the ethylene carbonate hydrolysis reaction at the electronic and molecular levels and simulate the spatial distribution characteristics within the storage container at the macroscopic level. It provides a scalable and computable foundation for subsequent embedding of dynamic moisture adsorption and reaction kinetics modules and conducting multiphysics coupled simulations.

[0045] After completing the quantitative calculation of the intrinsic properties of ethylene carbonate molecules and the core kinetic parameters of the hydrolysis reaction, and constructing a finite element spatial discretization model equipped with adaptive mesh refinement technology to form a complete digital twin basic model framework, in order to accurately reproduce the entire process of quality degradation caused by the easy hydrolysis of ethylene carbonate in actual closed storage scenarios, a dynamic water adsorption and reaction kinetics module designed based on the easy hydrolysis characteristics of ethylene carbonate needs to be embedded in the digital twin basic model. This module takes the water permeation-adsorption process at the interface of the closed container material as the starting point and the hydrolysis reaction kinetics of the liquid phase as the core, realizing the full-chain coupling from the container wall material properties to the microscopic chemical reaction. The specific implementation method is as follows:

[0046] First, the structural characteristics of the inner wall material of the sealed container are described by introducing the porous media permeability tensor. The porous media permeability tensor is a second-order symmetric tensor describing the ability of the polymer material's pore structure to allow water molecules to permeate, unlike a single permeability scalar. It accurately characterizes the anisotropic porosity of the polymer material due to differences in molecular chain arrangement and crystallinity, comprehensively describing the ease with which water molecules permeate in different spatial directions of the material. Each element in the tensor corresponds to a permeability component in a specific direction; a larger value indicates easier water molecule permeation in that direction. The specific introduction process is as follows:

[0047] First, the polymer material on the inner wall of the target container was characterized in three dimensions using scanning electron microscopy (SEM) and X-ray micro-CT to obtain core structural parameters such as porosity, average pore size, pore connectivity, and crystallinity distribution, thus clarifying the heterogeneous pore structure characteristics of the material. Then, based on the characterization results, permeability tensors were constructed, corresponding one-to-one with the grid cells of the container wall in the finite element spatial discretization model. For each grid cell, three principal values ​​of the tensor were calibrated according to the material crystallinity and pore structure at its location, corresponding to the axial, radial, and circumferential permeability of the material, respectively. For highly crystalline regions with a crystallinity higher than 60%, the principal value is lowered to the order of 1×10^-18 m² due to the tight molecular chain arrangement and low effective porosity. For amorphous regions with low crystallinity, the principal value is raised to the order of 1×10^-16 m² due to the larger molecular chain gaps and easier diffusion of water molecules. Finally, the permeability tensor of each grid cell is embedded into the properties of the container wall computational domain, providing basic parameters at the material structure level for the subsequent numerical solution of the water molecule diffusion and permeation process, and realizing a precise mapping from the actual microstructure of the material to the macroscopic permeability characteristics.

[0048] After constructing the material permeability tensor, a nonlinear mapping relationship between the water molecule adsorption isotherm and the local concentration gradient of ethylene carbonate was established by combining the diffusion barrier function of water molecules in the polymer material. The diffusion barrier function is a quantitative function describing the energy barrier that water molecules need to overcome during diffusion between polymer chains, formed by intermolecular forces and steric hindrance. Its value is directly related to ambient temperature, polymer crystallinity, and local water molecule concentration. The higher the barrier value, the more energy is required for water molecule diffusion, and the slower the diffusion rate. This function is the core bridge connecting the material's adsorption characteristics and water molecule diffusion behavior. The specific construction and mapping process is as follows:

[0049] A water molecule adsorption isotherm model adapted to polymer materials was constructed. Adsorption isotherms describe the relationship between the equilibrium content of adsorbed water molecules on the surface and within the pores of a polymer material at a constant temperature and the relative humidity of the surrounding environment. For polymer materials commonly used in storage containers, the GAB adsorption isotherm model was selected. This model can accurately cover the water vapor adsorption behavior of polymers over a wide humidity range. Through pre-conducted water vapor adsorption experiments at different temperatures and relative humidities, the model's core parameters, such as the adsorption equilibrium constant and monolayer adsorption capacity, were fitted to establish the correspondence between "temperature-relative humidity-material equilibrium adsorbed water amount," providing a basis for subsequent adsorption capacity calculations. Based on this foundation, a diffusion barrier function related to temperature and concentration is constructed. Using the activation energy of water molecule diffusion in the amorphous and crystalline regions of the polymer obtained from molecular dynamics simulations, and combining the Arrhenius relation, a diffusion barrier function is constructed. A local water molecule concentration correction term is introduced into the function. When the local water molecule concentration within the material increases, water molecules plasticize the polymer molecular chains, reducing the intermolecular forces and thus lowering the diffusion barrier. This achieves dynamic adjustment of the barrier value with local concentration, rather than a fixed constant. A nonlinear mapping relationship between the adsorption isotherm and the local concentration gradient of ethylene carbonate is established. Firstly, based on the adsorption isotherm model, and combined with the real-time temperature outside the container... Relative humidity data was used to calculate the equilibrium adsorbed water volume on the outer surface of the container wall material, which was then used as the outer boundary condition for the mass transfer process. Subsequently, the diffusion barrier function was introduced into the unsteady-state diffusion equation to solve for the instantaneous concentration distribution of water molecules at different locations within the container wall material. The instantaneous flux of water molecules diffusing from the inner surface of the container wall material to the bulk ethylene carbonate liquid phase was then calculated. This water molecule flux was then used as the water source term for the bulk liquid phase. Combined with the liquid phase mesh elements of the finite element spatial discretization model, the local water concentration within each element was calculated. Then, through the material conservation relationship, the local concentration gradient formed by the consumption of ethylene carbonate in the hydrolysis reaction was obtained. Finally, the non-uniformity of the two was clarified. The linear mapping logic states that as the water molecule content adsorbed on the container wall increases, the diffusion barrier decreases, the water molecule permeation flux increases, and the local water concentration in the liquid phase increases. This accelerates the hydrolysis reaction rate of ethylene carbonate, leading to a further increase in the local concentration gradient of ethylene carbonate. Conversely, the decrease in ethylene carbonate concentration alters the solvation effect of the liquid phase, which in turn affects the diffusion rate of water in the liquid phase and the adsorption equilibrium at the container wall interface. This forms a strong nonlinear coupling between adsorption behavior and the reaction process, ultimately establishing a complete nonlinear mapping relationship from the water molecule adsorption isotherm to the local concentration gradient of ethylene carbonate, thus realizing the linkage between the adsorption behavior at the material interface and the main reaction process in the liquid phase.

[0050] Based on this, the activation energy of the hydrolysis reaction is dynamically updated according to the transition state theory. The transition state theory, also known as the activated complex theory, is the core theory describing the rate of elementary chemical reactions. Its core logic is as follows:

[0051] For a chemical reaction to occur, reactant molecules need to collide and form activated complexes, or transition states, with energy higher than that of the reactants and products. Only by crossing the reaction energy barrier corresponding to the transition state, i.e., the activation energy, can the reactants be converted into products. The reaction rate is determined by the Gibbs free energy difference between the transition state and the reactants. This theory provides a theoretical basis for the dynamic correction of the activation energy of hydrolysis reactions and forms a fundamental connection with the transition state energy and reaction energy barrier obtained from quantum chemical calculations. The specific dynamic update process is as follows:

[0052] Using the energy barrier of the ethylene carbonate hydrolysis reaction obtained from quantum chemical calculations and molecular dynamics simulations as the basic value of the activation energy, which corresponds to the activation energy of the hydrolysis reaction at 25℃, atmospheric pressure, and a reactant concentration of 1 mol / L under standard conditions, as a benchmark for dynamic updates, a multi-factor dynamic correction model for the activation energy is established. Based on the core formula of transition state theory, three core correction terms are introduced: temperature, local moisture concentration, and product concentration, to achieve real-time dynamic updates of the activation energy. First, the temperature correction term, based on the thermodynamic relationship between the Arrhenius equation and transition state theory, shows that when the storage environment temperature changes, the enthalpy change and entropy change of the reaction system change accordingly, leading to a change in the free energy difference between the transition state and the reactants. The activation energy adjustment value at different temperatures is calculated through the temperature correction term; for every 10℃ increase in temperature, the activation energy decreases by 3%-5%, matching the actual law that increasing temperature accelerates the reaction rate. Second, the moisture concentration correction term addresses the issue that water molecules are a core reactant in the ethylene carbonate hydrolysis reaction. When the local moisture concentration in the liquid phase increases, the collision probability of reactant molecules significantly increases, and the probability of transition state formation increases accordingly. Based on the activity coefficient of transition state theory... The correction method involves adjusting the activation energy by 2% for every 0.5% increase in local moisture concentration exceeding 1% of the ethylene carbonate molar concentration. This reflects the equivalent reduction in the reaction energy barrier due to reactant concentration. Finally, there is a product inhibition correction term. Ethylene glycol, carbon dioxide, and acidic products generated from the hydrolysis of ethylene carbonate accumulate in the liquid phase. Increased product concentration shifts the reaction equilibrium towards the reverse reaction, increasing the apparent activation energy of the forward reaction. Based on the transition state theory and the reverse reaction energy barrier relationship, when the total product concentration exceeds 0.1 mol / L, the activation energy is reduced by 2% for every 0.05 mol / L increase. The activation energy was increased by 4%, reflecting the inhibitory effect of product accumulation on the hydrolysis reaction. The third step is to achieve real-time grid-level updates of the activation energy. The modified activation energy is synchronously updated to each liquid phase grid cell of the finite element spatial discretization model. Each grid cell independently calculates and updates the corresponding hydrolysis reaction activation energy based on its real-time temperature, moisture concentration, ethylene carbonate concentration, and product concentration, rather than using a uniform fixed value across the entire region. This accurately captures the differences in reactivity at different spatial locations, providing dynamic kinetic core parameters for the accurate calculation of the subsequent hydrolysis reaction rate.

[0053] Ultimately, through the parameter construction and dynamic calculation of the entire process described above, a deep coupling between material interface properties and chemical reactions is achieved. The material interface properties, represented by the permeability tensor of the porous medium in the container wall and the adsorption-diffusion behavior of water molecules, are linked in a closed loop with the chemical reactions, represented by the hydrolysis reaction kinetics with dynamically updated activation energy. The diffusion and permeation process of water molecules within the container wall material is used as the reactant input boundary for the hydrolysis reaction, and the water consumption by the hydrolysis reaction is used as the source and sink terms of the liquid phase water concentration field. This inversely affects the water concentration difference between the inner surface of the container wall and the bulk liquid phase, thereby altering the adsorption equilibrium and diffusion flux of water molecules. Simultaneously, the hydrolysis reaction rate of each grid cell is measured... The calculation results are fed back to the local concentration gradient calculation of ethylene carbonate in real time. Then, the mapping relationship of the water molecule adsorption isotherm is corrected in reverse by the concentration gradient, forming a complete coupled link of "water adsorption at the material interface - diffusion - liquid phase hydrolysis reaction - concentration gradient change - reverse correction of adsorption and diffusion behavior". This allows the module to accurately describe the obstruction effect of the container material on water permeation and fully restore the dynamic change process of the ethylene carbonate hydrolysis reaction. It provides a complete coupled calculation framework for the real-time calculation of water molecule diffusion permeation rate, hydrolysis reaction rate and product concentration distribution. It also lays the core calculation foundation for the construction of the subsequent dynamic multiphysics real-time feedback module.

[0054] After completing the core framework construction of the dynamic water adsorption and reaction kinetics module and achieving deep coupling between the container material interface properties and hydrolysis reaction kinetics, the module will enter the core numerical calculation stage. It will sequentially solve for the real-time water molecule diffusion and permeation rate, hydrolysis reaction rate, and product concentration distribution, providing accurate dynamic data support for subsequent multiphysics coupling and stability risk assessment. The specific implementation method is as follows:

[0055] First, real-time calculation of the water molecule diffusion and permeation rate is carried out. The core is to solve the unsteady-state diffusion equation using the porous medium permeability tensor and a preset water molecule diffusion barrier function, ultimately obtaining the instantaneous flux of water molecules within the container wall. The unsteady-state diffusion equation is a partial differential equation describing the dynamic changes in the concentration of the diffusing substance with time and spatial location. Unlike the steady-state diffusion equation where the concentration does not change with time, it can accurately reproduce the evolution of water molecule concentration within the sealed container wall over storage time, fully capturing the dynamic process of water permeating from the external environment to the internal ethylene carbonate liquid phase. It is the core governing equation for calculating the water molecule diffusion and permeation rate. The specific solution process is executed in three progressive steps:

[0056] In the finite element spatial discretization model, all mesh elements corresponding to the container wall are used as the solution domain. Each mesh element is bound to a pre-calibrated porous media permeability tensor to ensure that the solution process can completely reproduce the anisotropic permeability characteristics of the material pores. Three types of boundary conditions are set to perfectly match the actual storage scenario: the outer boundary where the container wall contacts the external environment, using a constructed GAB adsorption isotherm model combined with real-time data streams from temperature and relative humidity sensors in the storage environment, calculates the equilibrium water molecule concentration on the outer surface of the container wall, which serves as a fixed-concentration first-type Dirichlet boundary condition; and the inner boundary where the container wall contacts the ethylene carbonate liquid phase, with the concentration determined by the liquid phase... The real-time moisture concentration of the main body is dynamically updated as a movable boundary condition that changes with the reaction process. The side edges of the container wall are set as flux-free Neumann boundary conditions to eliminate invalid water permeation from the side of the container wall, ensuring that the solution domain is completely consistent with the actual container structure. Based on the classic Fick diffusion law, the permeability tensor of the porous medium is first integrated into the calculation of the diffusion coefficient. Based on the coupling relationship between Darcy's law and the diffusion law, a diffusion coefficient tensor corresponding one-to-one with the permeability tensor is constructed, allowing the diffusion coefficient to fully reflect the differences in pore permeation in different directions of the polymer material. Then, a preset water molecule diffusion barrier function is introduced to dynamically correct the diffusion coefficient tensor. The Nieus relation correlates the diffusion barrier function with ambient temperature and local water molecule concentration. When the local water molecule concentration within the material increases, triggering polymer plasticization, the diffusion barrier decreases, and the diffusion coefficient automatically increases; conversely, it decreases. This ultimately forms an unsteady-state diffusion equation that can match material properties, environmental conditions, and local concentration changes in real time, avoiding calculation errors caused by a fixed diffusion coefficient. Based on the finite element method, the constructed unsteady-state diffusion equation is spatially discretized. The Galerkin weighted residual method is used to transform the partial differential equation into a solvable system of linear algebraic equations. In the time dimension, the unconditionally stable implicit Euler method is used for discretization, and a 1-hour time step is set to match the sensor sampling frequency. To ensure real-time and temporal synchronization of the calculations, within each time step, the system first synchronously updates the current ambient temperature and local water molecule concentration, then recalculates the diffusion barrier function and the corrected diffusion coefficient tensor, and then solves the linear algebraic equations to obtain the water molecule concentration distribution of each grid cell within the container wall at the current time step. Finally, based on Fick's first law, the instantaneous flux of water molecules on the inner surface of the container wall is calculated, which is the mass of water molecules that permeate into the liquid phase per unit area per unit time. This flux is the water molecule diffusion and permeation rate at the current moment. At the same time, the concentration distribution at this time step is used as the initial condition for solving the next time step, realizing the dynamic iterative calculation of the diffusion process.After obtaining the water molecule diffusion and permeation rate and determining the real-time water input of each grid cell in the liquid phase, the module will simultaneously calculate the hydrolysis reaction rate. The core is to deduce the reactant consumption per unit time using a microscopic reaction probability model based on the local concentration gradient of ethylene carbonate and the updated activation energy of the hydrolysis reaction.

[0057] Among them, the microscopic reaction probability model starts from the microscopic perspective of molecular collisions and effective reactions. Based on the local concentration of reactant molecules and the dynamic reaction rate constant, it calculates the number of molecules that undergo effective hydrolysis reaction per unit time and unit volume, and then derives the calculation model for reactant consumption. Unlike the macroscopic total reaction rate equation, it can accurately capture the influence of the local concentration gradient in the liquid phase on the probability of reaction, perfectly adapting to the fine calculation needs of finite element mesh micro-elements. The specific calculation process is executed according to the following logic:

[0058] From each grid cell in the liquid phase region of the finite element spatial discretization model, the local concentrations of ethylene carbonate and water molecules at the current moment, as well as the activation energy of the hydrolysis reaction dynamically updated according to transition state theory, are extracted. This ensures that each grid cell has independent reaction calculation parameters, rather than using a uniform fixed value across the entire region, accurately reproducing the differences in reactivity at different spatial locations. Simultaneously, based on the Arrhenius equation, the reaction rate constant for each grid cell is calculated using the dynamically updated activation energy of the hydrolysis reaction. The pre-exponential factor is determined by the transition state partition function of the hydrolysis reaction obtained from quantum chemical calculations and serves as a fixed fundamental constant. This ultimately achieves real-time dynamic calibration of the reaction rate constant for each grid cell. Considering the actual situation of excess water molecules in the ethylene carbonate storage scenario, the hydrolysis reaction is defined as a pseudo-first-order reaction, and the probability of the reaction occurring is mainly determined by the local concentration of ethylene carbonate. Based on the reaction rate constant and the calculation time step, the probability of a single ethylene carbonate molecule undergoing effective hydrolysis within a single time step is calculated. This probability increases with the increase of the reaction rate constant and decreases with the decrease of the time step, perfectly matching the time-cumulative effect of the reaction. First, based on the local concentration of ethylene carbonate and the volume of the grid cell, the total number of ethylene carbonate molecules in the grid is calculated. Then, combined with the effective reaction probability, the number of molecules undergoing hydrolysis in a single time step is calculated. This is converted into the mass consumption of ethylene carbonate based on the stoichiometric ratio. Finally, this is divided by the time step to obtain the ethylene carbonate reactant consumption per unit time in that grid cell, which is the hydrolysis reaction rate at the current moment. Simultaneously, based on the stoichiometric ratio of the hydrolysis reaction, the water molecule consumption and product generation per unit time are calculated to provide basic data for subsequent product concentration distribution calculations. During the calculation process, a concentration gradient correction term for adjacent grids is also introduced. When there is a ethylene carbonate concentration difference between adjacent grids, the mass transport flux is calculated through the convection-diffusion equation, the local concentration of the current grid is updated, and the reaction consumption is recalculated to ensure that the calculation results completely match the actual concentration distribution in the liquid phase. After completing the calculation of the hydrolysis reaction rate and reactant consumption, the module will construct the product transport equation based on the reactant consumption and water molecule diffusion flux, and solve the three-dimensional spatial concentration field through the principle of mass conservation, thus fully reconstructing the accumulation and diffusion process of hydrolysis products.

[0059] First, it is clear that the core products of the hydrolysis reaction of ethylene carbonate are ethylene glycol, carbon dioxide, and organic acids. These products are the key factors affecting the acid value, color, and other critical quality indicators of ethylene carbonate, and are also the key basis for subsequent stability degradation risk assessment. Therefore, the product transport equation mainly focuses on solving the two core components: ethylene glycol and acidic products. The specific solution process is as follows:

[0060] First, a product transport equation is constructed based on the principle of mass conservation. The core logic of the mass conservation principle is that the change in the mass of the product within a grid cell per unit time is equal to the product flux entering the cell through diffusion and convection, plus the amount of product generated by the hydrolysis reaction within the cell, minus the product flux leaving the cell. The product transport equation constructed based on this includes a diffusion term describing the spatial diffusion of the product, a source term describing the product formation rate, and a natural convection correction term induced by the temperature gradient, which closely matches the product transport characteristics within the liquid phase in actual storage scenarios. The source term is directly obtained from the calculated reactant consumption per unit time, converted according to the stoichiometric ratio of the hydrolysis reaction, ensuring complete synchronization between the reaction process and the product formation process. Then, the initial and boundary conditions of the transport equation are determined. The initial condition is set as the product concentration in the container at the start of storage, which can be either 0% of the new material or the initial product concentration of the actually detected material, fully adapting to different simulation scenarios. The boundary condition sets the contact surface between the container wall and the liquid phase as a zero-flux boundary, because the hydrolysis products cannot be... Permeation through the polymer container wall allows diffusion only within the liquid phase, ensuring boundary conditions are completely consistent with the actual storage scenario. Then, equation solving and 3D concentration field reconstruction are performed. Similarly, based on the finite element method, all grid cells in the liquid phase region of the product transport equation are spatially discretized, transforming it into a system of linear algebraic equations. The time dimension uses the 1-hour implicit Euler method, consistent with diffusion and reaction calculations, ensuring complete synchronization of the entire module's computational timing. Within each time step, the system first updates the product generation rate source term based on the hydrolysis reaction rate of the previous step, then corrects the product diffusion coefficient based on temperature and concentration changes, solving the linear equations to obtain the product concentration value of each liquid phase grid cell at the current time step. Finally, a spatial interpolation algorithm is used to interpolate the concentration data of all grid cells, reconstructing the 3D product concentration field of the entire liquid phase region. This visually presents the spatial distribution characteristics and aggregation areas of the product within the container, while simultaneously recording the concentration field data for each time step, providing complete data support for subsequent product accumulation effect simulation and stability decay risk index calculation.

[0061] The entire calculation process forms a complete dynamic closed loop. The calculation results of the water molecule diffusion and permeation rate provide the boundary conditions of liquid phase water concentration for the calculation of the hydrolysis reaction rate. The consumption of reactants in the hydrolysis reaction will change the local concentration gradient between ethylene carbonate and water molecules, and inversely correct the inner boundary conditions of the unsteady diffusion equation of water molecules. The calculation results of product concentration distribution will dynamically update the activation energy of the hydrolysis reaction through the product inhibition coefficient, thereby adjusting the hydrolysis reaction rate and realizing dynamic linkage of the whole process. This ensures that each set of data output by the module can accurately reproduce the dynamic changes of the hydrolysis reaction during the storage of ethylene carbonate, and provides the core dynamic calculation foundation for the subsequent construction of the dynamic multiphysics real-time feedback module.

[0062] After completing the core numerical calculations of water molecule diffusion and permeation rate, hydrolysis reaction rate, and product concentration distribution in the dynamic moisture adsorption and reaction kinetics module, and establishing a coupled calculation framework between container material interface properties and hydrolysis reaction kinetics, in order to further realize the real-time linkage between dynamic changes in the storage environment and the hydrolysis reaction process of ethylene carbonate, and accurately reproduce the entire process of material quality degradation under the interaction of multi-physics fields in a closed container, a dynamic multi-physics field real-time feedback module needs to be designed based on the core calculation results output by the aforementioned module. This module takes the full coupling of multi-physics fields as its core logic, real-time sensor data as its driving boundary, and time-series synchronous control and adaptive grid optimization as its computational guarantee, ultimately realizing closed-loop feedback and dynamic adjustment from changes in environmental parameters to the hydrolysis reaction process. The specific implementation method is as follows:

[0063] The first step in module construction is to define four types of physical fields and accurately map them to core input variables. This establishes clear physical boundaries and data links for subsequent multi-field coupled calculations. These four physical fields correspond to the mass transfer field, chemical reaction field, convection-diffusion field, and heat transfer field, respectively. Each field is deeply integrated with the calculation results to ensure the continuity and accuracy of the module input. The core physical meaning of the mass transfer field is to describe the entire process of water molecules diffusing and transporting from the external environment through the container wall material to the bulk ethylene carbonate liquid phase. Its core input variable is the water molecule diffusion and permeation rate calculated in real time by the dynamic water adsorption and reaction kinetics module. This rate directly serves as the core source term of the mass transfer field's governing equation, while simultaneously relating it to the porous medium permeability tensor of the container wall material and the water molecule diffusion barrier function. The constructed unsteady diffusion equation is used as the core governing equation for the mass transfer field to ensure the consistency of the water diffusion process calculation. The core physical meaning of the chemical reaction field is to describe the dynamic process of the hydrolysis reaction between ethylene carbonate and water molecules. Its core input variable is the hydrolysis reaction rate calculated by the module, i.e., the reaction rate per unit time derived from the microscopic reaction probability model. The system considers the consumption of materials and the activation energy of the hydrolysis reaction, which is dynamically updated based on the transition state theory. The hydrolysis reaction rate equation, constructed based on the Arrhenius formula, serves as the core governing equation, fully reproducing the dynamic changes in reaction rate with environmental and material states. The core physical meaning of the convection-diffusion field is to describe the transport, diffusion, and aggregation processes of products such as ethylene glycol and organic acids generated by the hydrolysis reaction in the ethylene carbonate liquid phase. Its core input variable is the product concentration distribution obtained by the module solution, i.e., the three-dimensional spatial concentration field calculated through the product transport equation. The convection-diffusion equation with a reaction source term serves as the core governing equation, where the source term is directly updated in real-time by the product generation rate output from the chemical reaction field, accurately reproducing the spatial distribution evolution of products in the liquid phase. The core physical meaning of the heat transfer field is to describe the transmission process of temperature changes in the storage environment between the container wall and the liquid phase, as well as the influence of temperature on diffusion and reaction processes. Its boundary conditions are directly driven in real-time by integrated temperature and relative humidity sensor data streams. The unsteady-state heat conduction equation, constructed based on the principle of energy conservation, serves as the core governing equation, acting as the core bridge connecting the actual storage environment and the digital twin simulation system.

[0064] After completing the framework construction of the four types of physical fields, the module will integrate the data streams from temperature and relative humidity sensors to achieve dynamic driving of the boundary conditions of the heat transfer field. This will enable the digital twin simulation system to respond in real time to changes in the actual storage environment. The specific implementation process is as follows:

[0065] First, in the actual warehouse environment where ethylene carbonate is stored, high-precision temperature and humidity sensors are deployed along the arrangement of sealed storage containers at a density of one set every 5 meters. The sampling frequency of the sensors is set to 1 hour / time, perfectly matching the main time step of the unsteady-state equation solution, ensuring that the sensor data and simulation calculations are synchronized. The raw temperature and humidity data collected by the sensors are first transmitted to the module's preprocessing unit. Abnormal jump values ​​caused by environmental disturbances are removed by a 3-point moving average filter, and missing data caused by signal interruption are filled by linear interpolation, finally obtaining a smooth and continuous effective environmental data stream. Subsequently, the preprocessed data stream is divided according to function. The relative humidity data is synchronously transmitted to the mass transfer field solver. Through the constructed GAB adsorption isotherm model, the equilibrium water molecule concentration on the outer surface of the container wall is updated in real time to correct the outer boundary conditions of the mass transfer field. The temperature data is directly used as the first type of Diri for the heat transfer field. The chlet boundary conditions drive the heat transfer field solver to calculate the three-dimensional temperature distribution within the entire container computational domain. Simultaneously, temperature data is transmitted to the mass transfer field and chemical reaction field. The diffusion barrier function corrects the water molecule diffusion coefficient, and the Arrhenius relation corrects the activation energy and rate constant of the hydrolysis reaction. This achieves a full-chain linkage effect of ambient temperature on diffusion and reaction processes, ensuring the simulation system perfectly reflects the dynamic changes of the actual storage environment. After integrating sensor data streams and setting physical field boundaries, the module coordinates the temporal synchronization of the four types of physical fields through a field coupling controller. This is the core element ensuring the convergence and accuracy of multi-physics coupling calculations. The field coupling controller, as the temporal hub of the entire module, is crucial for resolving the temporal matching of the four types of physical field solutions and real-time data interaction across fields. This avoids computational deviations and result distortions caused by inconsistent solution step sizes and asynchronous data transmission between different fields. The specific execution logic is as follows:

[0066] First, the timing initialization settings are completed, with a main calculation step of 1 hour, completely consistent with the sensor sampling frequency and the time step of single-step solution. At the beginning of each main step, the field coupling controller first synchronously acquires the preprocessed temperature and humidity sensor data stream at the current moment and sends it to the solvers of the heat transfer field and mass transfer field respectively, as the initial values ​​of the boundary conditions for the calculation of this step. At the same time, it retrieves the global water molecule concentration distribution output by the mass transfer field, the hydrolysis reaction rate distribution output by the chemical reaction field, and the product concentration distribution output by the convection diffusion field after the end of the previous main step, as the initial conditions for solving each field in this step, ensuring the temporal continuity of the entire simulation process. The data breakpoints and time were misaligned. Subsequently, the field coupling controller set a progressive solution sequence of "heat transfer field → mass transfer field → chemical reaction field → convection diffusion field." This sequence perfectly matches the logical relationship of the actual physical process—ambient temperature is the core prerequisite factor affecting the diffusion and hydrolysis reaction rate of water molecules. Therefore, the heat transfer field is solved first to obtain the three-dimensional temperature distribution in the entire computational domain. The temperature distribution is then written into the module's shared data pool in real time. The mass transfer field solver synchronously retrieves the temperature distribution from the shared data pool, updates the water molecule diffusion barrier function and diffusion coefficient tensor, completes the solution of the mass transfer field for this step, and obtains the updated water molecule diffusion permeation rate and liquid phase main... The global moisture concentration distribution is also written into the shared data pool. Next, the chemical reaction field solver retrieves the latest moisture concentration and temperature distributions from the shared data pool, dynamically updates the activation energy and rate constant of the hydrolysis reaction, and completes the chemical reaction field solution for this step. This yields the updated hydrolysis reaction rate, reactant consumption per unit time, and product formation, which are then written into the shared data pool. Finally, the convection-diffusion field solver retrieves the product formation data, updates the source terms of the product transport equation, and completes the convection-diffusion field solution for this step, obtaining the latest three-dimensional product concentration distribution. During the solution process for each field, the field coupling controller verifies the solution results in real time. Convergence is ensured by the fact that the solution for the current field only starts when the solution meets the convergence criterion (residual less than 1×10^-6), thus avoiding subsequent calculation deviations caused by non-convergent data transfer. Based on this, the field coupling controller dynamically adjusts the sub-solution step size of each field according to the degree of nonlinearity and convergence speed of each field. For heat transfer fields with high linearity and fast convergence speed, the sub-step size is appropriately shortened to improve computational efficiency. For chemical reaction fields with strong nonlinearity and slow convergence speed, the sub-step size is appropriately extended to ensure computational convergence, ensuring that the solutions for all four fields are completed before the end of the main step, achieving complete timing synchronization of the entire process.

[0067] After each main step is solved, the field coupling controller uses the final solution results of the four fields in this step as the initial conditions for the next main step, and simultaneously updates the sensor data stream at the next moment, starting a new iteration cycle. This ensures the time-series closed loop and continuous iteration of the entire simulation process. After completing the time-series synchronous solution of the multiphysics field, the module will establish an adaptive mesh re-division mechanism with the product accumulation rate as the feedback signal. This achieves an optimal balance between ensuring the computational accuracy of the core region and the overall computational efficiency. This mechanism is a deepening and dynamic optimization of the adaptive mesh refinement technology in the digital twin basic model, using the dynamic changes of the actual reaction process as the core basis for mesh adjustment. The specific implementation process is as follows:

[0068] First, the core definition of the feedback signal is clarified. The product accumulation rate refers to the increase in the concentration of hydrolysis products per unit time and per unit volume. It is calculated from the global product concentration distribution obtained by solving the convection-diffusion field, using the concentration difference between the current step size and the previous step size. This value directly reflects the intensity of the hydrolysis reaction within the corresponding grid cell and the magnitude of the physical quantity gradient change, and is the core trigger signal for grid remapping. Subsequently, a three-level grid remapping triggering criterion is established. The first level is the threshold for high-reaction, high-gradient regions. When the product accumulation rate of a single grid cell exceeds three times the average value of the entire computational domain, or the product accumulation rate gradient between adjacent grid cells exceeds 10%, it is determined to be a core region with intense hydrolysis and rapid physical quantity changes, triggering grid refinement. The second level... Level 1 is the threshold for the medium-reaction gradient region. When the product accumulation rate of a grid cell is between 1 and 3 times the average value of the entire region, or the rate gradient between adjacent grid cells is between 5% and 10%, it is determined to be a region with moderate reaction and gradual changes in physical quantities, and the current grid size remains unchanged. Level 3 is the threshold for the low-reaction, low-gradient region. When the product accumulation rate of a grid cell is lower than the average value of the entire region, and the rate gradient between adjacent grid cells is lower than 5%, it is determined to be a region with extremely low reaction and almost no change in physical quantities, triggering grid coarsening. At the same time, the two-phase interface region between the container wall and the liquid bulk, and the permeation boundary region with a large water concentration gradient are set as forced densification zones. Regardless of the product accumulation rate, a fine grid is always maintained to ensure the calculation accuracy of the interface mass transfer process. Based on this triggering criterion, the system will dynamically re-mesh the finite element mesh across the entire computational domain. For mesh elements that trigger refinement, they will be subdivided into eight sub-meshes, each half the size of the original mesh. For mesh elements that trigger coarsening, eight adjacent low-gradient meshes will be merged into one large mesh. During the re-meshing process, the topological continuity of the mesh will be strictly guaranteed, with the size ratio of adjacent meshes not exceeding 2:1, to avoid computational errors and numerical dissipation caused by abrupt changes in mesh size. Simultaneously, a conformal interpolation algorithm will be used to accurately map all physical quantities, such as temperature, concentration, and reaction rate, within the original mesh elements to the newly re-meshed mesh elements, ensuring the spatial continuity of physical quantities without data loss or bias. After the re-meshing is completed, the new mesh model will be synchronously transmitted to the solvers of the four physical fields. This forms the computational grid for the next main step, creating a closed-loop feedback loop of "product accumulation rate calculation → adaptive grid re-partitioning → improved multiphysics solution accuracy → more accurate product accumulation rate calculation." This ensures computational accuracy in high-reaction, high-gradient core regions while significantly reducing computational load in low-gradient regions, achieving an optimal balance between computational accuracy and efficiency for the entire simulation system. Ultimately, through fully coupled multiphysics solution and dynamic feedback, the module adjusts the water molecule diffusion and hydrolysis reaction rates in real time, accurately simulating the dynamic balance of the water gradient and product accumulation effects in a closed container, and outputting a stability decay risk index. Specifically, through synchronous solution by the field coupling controller, the real-time temperature distribution of the heat transfer field is continuously updated to the mass transfer field.By correcting the water molecule diffusion coefficient and diffusion barrier, the real-time dynamic adjustment of the water molecule diffusion and permeation rate is achieved. The water concentration distribution in the mass transfer field, the temperature distribution in the heat transfer field, and the product concentration distribution in the convection-diffusion field are synchronously updated to the chemical reaction field, dynamically correcting the activation energy and reaction rate constant of the hydrolysis reaction, and achieving real-time dynamic adjustment of the hydrolysis reaction rate. At the same time, the product concentration output from the chemical reaction field is fed back to the mass transfer field through the product inhibition coefficient to correct the diffusion coefficient, and back to the chemical reaction field to correct the reaction activation energy, forming a two-way real-time feedback link between multiple physical fields. This fully recreates the entire dynamic process of water in a closed container from permeation and adsorption to participation in the reaction, and then to product accumulation and reaction inhibition, accurately simulating the water gradient. The module considers the dynamic equilibrium and product accumulation effect. Based on this, it uses a spatial interpolation algorithm to reconstruct the three-dimensional moisture distribution field inside the container using the final results of multi-physics coupled solutions. Combined with the output results of the convection-diffusion field, it locates the product accumulation region. Through normalization, the extreme values ​​of the moisture gradient are mapped to a moisture permeation risk component, and the product accumulation concentration is mapped to a product corrosion risk component. The two components are weighted and summed in a 4:6 ratio to obtain a stability degradation risk index ranging from 0 to 1. A higher value indicates a greater risk of storage stability degradation for ethylene carbonate. This index will serve as the core output, providing core data support for subsequent prediction of key quality indicator evolution trends and failure threshold identification.

[0069] After completing the physical field framework construction, sensor data stream integration, and timing synchronization mechanism design of the dynamic multiphysics real-time feedback module, the core task of the module is to construct a set of heat transfer, mass transfer, and chemical reaction kinetic equations that can fully describe the intrinsic relationship of each physical process and fit the actual storage scenario of ethylene carbonate in order to achieve deep coupling and accurate numerical solution of heat transfer, mass transfer, and hydrolysis chemical reaction processes. The entire construction process follows the progressive logic of "accurate establishment of single-field equations - introduction of coupling terms to realize multi-field simultaneous equations - decoupling and efficient solution using operator splitting method". This not only fully restores the essential laws of each physical process, but also solves the problems of high solution complexity and poor convergence caused by strong coupling of multiple fields. The specific implementation method is as follows:

[0070] First, the independent construction of the single-physics field governing equations is carried out to lay a precise physical foundation for subsequent multi-field simultaneous equations. The establishment of the three core equations strictly follows the corresponding physical conservation laws and reaction kinetics theories, and is deeply integrated with the core parameters of the constructed digital twin basic model and the dynamic water adsorption and reaction kinetics module. The first equation is a water diffusion equation based on the principle of mass conservation. The principle of mass conservation states that within a fixed control volume, the change in mass of a substance per unit time is equal to the difference between the mass flux flowing in and out through the boundary of the control volume, plus the difference between the amount of the substance generated and consumed within the control volume. This is the core criterion for describing the diffusion and transport process of substances. Here, the control volume is each grid element in the finite element spatial discretization model, ensuring that the equation solution is fully adapted to the spatial discretization model. The specific establishment process is as follows:

[0071] First, based on Fick's diffusion law, the molecular diffusion process of water molecules in the polymer material of the container wall and the liquid phase of ethylene carbonate is described. The diffusion coefficient is directly adopted from the constructed dynamic diffusion coefficient tensor, which is related to the permeability tensor of porous media and the diffusion barrier function of water molecules. This fully reflects the dynamic influence of material pore anisotropy, temperature, and local moisture concentration on diffusion capacity. Then, considering the actual characteristics of ethylene carbonate storage scenarios, a water consumption term due to hydrolysis is introduced into the equation. This consumption term is directly related to the calculation results of the subsequent hydrolysis reaction rate equation, reflecting the direct influence of chemical reaction on mass transfer process. Finally, an unsteady-state moisture diffusion equation is constructed. The left side of the equation is the partial derivative of the water concentration in the control volume with respect to time, representing the change in water mass in the control volume per unit time. The first part of the right side of the equation is the divergence term of diffusion flux, describing the diffusion and transport process of water molecules in space. The first equation consists of negative reaction consumption terms, representing the real-time water consumption of the hydrolysis reaction. It comprehensively covers the entire process of water permeating from the container wall, diffusing in the liquid phase, and being consumed in the hydrolysis reaction. Simultaneously, the outer boundary conditions of the equation are tied to the equilibrium water molecule concentration on the outer surface of the container wall calculated by the GAB adsorption isotherm model, while the inner boundary conditions are linked to the real-time water concentration in the bulk liquid phase, ensuring that the equation boundaries perfectly match the actual storage scenario. The second equation is an unsteady-state heat conduction equation based on the principle of energy conservation. The principle of energy conservation states that within a fixed control volume, the change in internal energy per unit time is equal to the heat flow through the control volume boundary, plus the heat generated by the heat source inside the control volume. This is the core criterion for describing the heat transfer process. This equation is used to reconstruct the transfer process of temperature changes in the storage environment within the container, as well as the key influence of temperature on diffusion and reaction processes. The specific establishment process is as follows:

[0072] First, temperature-dependent thermophysical properties, including thermal conductivity, density, and isobaric specific heat capacity, were calibrated for both the container wall polymer material and the ethylene carbonate liquid phase to ensure dynamic updates with temperature changes. Then, based on Fourier's law of heat conduction, the heat conduction process in the container wall and liquid phase was described, while an internal heat source term was introduced. The intensity of the internal heat source is directly proportional to the hydrolysis reaction rate, thus reducing the influence of the thermal effect of the hydrolysis reaction on the system temperature. Finally, an unsteady-state heat conduction equation was constructed. The left side of the equation represents the partial derivative of the internal energy per unit volume of the control volume with respect to time, obtained by multiplying the material density, isobaric specific heat capacity, and the partial derivative of temperature with respect to time, characterizing the change in internal energy of the control volume per unit time. The right side of the equation consists of two parts: the first part is the divergence term of the heat conduction flux, describing the heat conduction process in space; the second part is the internal heat source term, characterizing the heat change caused by the hydrolysis reaction. The outer boundary conditions of the process are directly driven by real-time data streams from integrated temperature sensors. The interface between the container wall and the liquid phase is set as a thermally continuous boundary to ensure the continuity of temperature and heat flux at the interface, fully reproducing the transfer process of ambient temperature changes to the interior of the container. Simultaneously, it provides global temperature field data for the subsequent dynamic updates of diffusion coefficients and reaction rate constants. The third equation is the hydrolysis reaction rate equation based on transition state theory. The core logic of transition state theory is clear: the occurrence of a chemical reaction requires reactant molecules to form a higher-energy transition state activated complex. Only by crossing the reaction energy barrier corresponding to the transition state can the reaction be completed. The reaction rate is directly determined by the Gibbs free energy difference between the transition state and the reactants. This equation is the core link connecting the mass transfer and heat transfer processes, and it is also the core governing equation describing the dynamics of the ethylene carbonate hydrolysis reaction. The specific establishment process is as follows:

[0073] First, based on the transition state energy and reaction energy barrier of the hydrolysis reaction obtained from quantum chemical calculations and molecular dynamics simulations, the pre-exponential factor and activation energy of the reaction were determined. The pre-exponential factor, calculated from the partition function of the transition state, is a fixed fundamental constant. The activation energy follows a dynamically corrected model constructed based on transition state theory, which can be updated in real time with temperature, local moisture concentration, and product concentration. Then, considering the actual storage scenario of ethylene carbonate, where ethylene carbonate is a large excess of the main material and moisture is a limiting reactant with trace permeation, the hydrolysis reaction is defined as a pseudo-first-order reaction. The reaction rate is mainly determined by the moisture concentration and the reaction rate constant. Finally, a hydrolysis reaction... The core expression of the rate equation is the amount of reactant consumed per unit volume of ethylene carbonate per unit time, which is equal to the product of the reaction rate constant, the concentration of ethylene carbonate, and the water concentration. The reaction rate constant is directly related to the system temperature and the dynamically updated activation energy through the Arrhenius equation, perfectly connecting the global temperature field output by the heat transfer equation and the water concentration field output by the mass transfer equation. The calculation result of the equation is the hydrolysis reaction rate in the dynamic water adsorption and reaction kinetics module. At the same time, the amount of water molecules consumed and the amount of products generated can be obtained synchronously through stoichiometry, providing accurate input data for the reaction consumption term of the mass transfer equation and the internal heat source term of the heat conduction equation.

[0074] After constructing the three single-physics field governing equations, the inhibition coefficient of product concentration distribution on the reaction rate is introduced, and the three are combined into a system of nonlinear partial differential equations with cross terms. The inhibition coefficient of product concentration distribution on the reaction rate is a dimensionless correction coefficient, ranging from 0 to 1. Its core physical meaning is to quantify the inhibitory effect of products such as ethylene glycol and organic acids generated in the hydrolysis reaction on the forward reaction—the accumulation of products in the liquid phase changes the polarity of the system, reduces the reactivity of water molecules, and simultaneously pushes the reaction equilibrium towards the reverse reaction direction. Apparently, this manifests as an increase in the apparent activation energy of the hydrolysis reaction and a decrease in the reaction rate. The higher the total product concentration, the smaller the inhibition coefficient, and the stronger the inhibitory effect on the reaction rate. The expression for this coefficient is obtained by fitting the hydrolysis reaction equilibrium experiment and shows a nonlinear relationship with a negative correlation to the total product concentration. The specific process of combining the equations is as follows:

[0075] First, a product inhibition coefficient is introduced into the hydrolysis reaction rate equation. Multiplying the original reaction rate by this inhibition coefficient establishes a direct correlation between the reaction rate and the product concentration distribution. The product concentration distribution is determined by the time integral of the reaction rate and the product transport process, forming an internal feedback loop in the reaction process. Then, the cross-coupling terms among the three equations are analyzed to clarify their strong correlation: the global temperature field output by the unsteady-state heat conduction equation appears simultaneously in the diffusion coefficient term of the moisture diffusion equation and the reaction rate constant term of the hydrolysis reaction rate equation, serving as the core coupling variable connecting the heat transfer process with the mass transfer and reaction processes; the global moisture concentration field output by the moisture diffusion equation is the core input term of the hydrolysis reaction rate equation, directly determining the reaction rate. The magnitude of the rate, and the amount of water consumed calculated by the hydrolysis reaction rate equation, in turn, serves as the source and sink terms in the water diffusion equation, thus inversely altering the distribution of the water concentration field; the amount of products generated calculated by the hydrolysis reaction rate equation, through the product inhibition coefficient, inversely corrects the reaction rate, while the reaction rate determines the magnitude of the heat source term in the heat conduction equation, thereby changing the temperature field distribution of the system; finally, the three equations are integrated into a complete set of nonlinear partial differential equations, in which the solution result of each equation is the input variable of the other two equations, forming a strongly coupled nonlinear equation set with multiple sets of cross-coupling terms, fully restoring the inherent physical logic of mutual influence and bidirectional linkage between the three processes of heat transfer, mass transfer, and hydrolysis reaction, without omitting any key processes. To address the challenges of solving strongly coupled nonlinear equation systems with poor convergence, an operator splitting method is employed. This method decomposes the equation system into independently solvable heat transfer, proton transfer, and reaction kinetics sub-equations. The operator splitting method, also known as the fractional-step method, is a classic numerical approach for solving multi-physics coupled nonlinear partial differential equation systems. Its core idea is to break down the complex multi-physics coupled process into multiple independent sub-processes, each corresponding to an independent sub-operator and sub-equation. Within the same time step, each sub-equation is solved sequentially according to the causal logic of the physical processes. The solution of the previous sub-equation serves as the input boundary condition for the next sub-equation. This approach not only preserves the coupling relationships between multiple fields but also decomposes the complex coupled equation system into multiple simple, independently solvable sub-equations, significantly reducing the difficulty of solving the system while ensuring computational efficiency and convergence.The specific decomposition and solution process strictly follows the physical causal logic of "temperature affects diffusion, diffusion determines reaction, and reaction in turn affects temperature and diffusion," perfectly matching the physical field solution sequence set by the field coupling controller. Within a 1-hour time step [t_n, t_{n+1}] that is completely consistent with the sensor sampling frequency and the main step size of the field coupling controller, it is executed in three progressive steps: The first step solves the heat transfer equation, corresponding to the independent operators of the heat transfer process. In this step, the variables of the mass transfer and reaction processes are temporarily frozen. The water concentration, product concentration, and reaction rate all adopt the convergence result of the previous time step t_n, retaining only the temperature-related terms in the heat conduction equation. By fixing the internal heat source term, the unsteady heat conduction equation in the original coupled equation set is decomposed into independent heat transfer sub-equations. The external boundary conditions of the sub-equations are based on the pre-processed real-time temperature sensor data stream at the current moment. The sub-equations are spatially discretized and solved using the finite element method to obtain the intermediate temperature field T^* of the entire computational domain at the current time step, thus completing the independent solution of the heat transfer process. Correspondingly, the independent operators for the mass transfer process are obtained. This step uses the intermediate temperature field T^* obtained in the first step to update the moisture diffusion coefficient tensor and diffusion barrier function in real time. At the same time, the reaction kinetics process is frozen, and the reaction consumption term uses the result of the previous time step t_n. The original coupled equations are then decomposed into independent heat transfer sub-equations. The moisture diffusion equation in the equation set is decomposed into independent proton transfer equations. The outer boundary condition of the sub-equations is the equilibrium water molecule concentration on the outer surface of the container wall, calculated from the GAB adsorption isotherm using relative humidity sensor data at the current time. This is also solved using the finite element method to obtain the intermediate moisture concentration field C_H2O^* in the entire computational domain at the current time step, thus completing the independent solution for the mass transfer process. Corresponding to the independent operators for the hydrolysis reaction process, this step uses the intermediate temperature field T^* and intermediate moisture concentration field C_H2O^* obtained in the first two steps. First, the activation energy and rate constant of the hydrolysis reaction are dynamically updated based on transition state theory, and then the current time step's... The product concentration distribution is used to calculate the product inhibition coefficient. The hydrolysis reaction rate equation in the original coupled equation set is decomposed into an independent reaction kinetic sub-equation. This sub-equation is an ordinary differential equation, which can be solved quickly using the fourth-order Runge-Kutta method to obtain the hydrolysis reaction rate r_{n+1} at the current time step. Then, based on the conservation of mass and stoichiometry, the total product concentration C_p^{n+1} at the current time step and the final water concentration field C_H2O^{n+1} after deducting the reaction consumption are calculated. At the same time, the temperature correction is obtained based on the reaction heat, and the final global temperature field T^{n+1} is updated to complete the independent solution of the reaction process.

[0076] After the three-step solution is completed in one time step, the convergence results of the temperature field, moisture concentration field, product concentration, hydrolysis reaction rate, etc., will be used as the initial conditions for the solution in the next time step. At the same time, the field coupling controller synchronously updates the sensor data stream at the next moment, starts a new solution cycle, and forms a complete time-series closed loop. The entire operator splitting process adopts an implicit solution format. The solution of each sub-equation is unconditionally stable and will not cause the problem of solution divergence due to the choice of time step. After the solution of each sub-equation is completed, convergence verification is performed. Only when the solution residual is less than 1×10^-6 will the next solution be entered, ensuring that the calculation accuracy after decoupling is completely consistent with that of the coupled solution. Through the entire process of constructing single-field equations, combining multiple fields, and decoupling solutions using the operator splitting method, the strong coupling inherent laws of heat transfer, mass transfer, and hydrolysis reaction during the storage of ethylene carbonate were fully restored. Furthermore, efficient and stable solutions to complex nonlinear equations were achieved. This provides a core numerical solution framework for the real-time dynamic adjustment of water molecule diffusion and hydrolysis reaction rates in the dynamic multiphysics real-time feedback module, as well as for the construction of a two-way real-time feedback link. It also lays a solid computational foundation for subsequent simulations of dynamic equilibrium of moisture gradients within sealed containers, analysis of product accumulation effects, and accurate output of the stability decay risk index.

[0077] After independently decoupling and solving the heat transfer, mass transfer, and reaction kinetic sub-equations using the operator splitting method, to achieve dynamic linkage and real-time adaptation between multiple physics fields and ensure that the simulation process accurately follows changes in the actual storage environment and the evolution of material reaction states, it is necessary to perform closed-loop adjustments to the water molecule diffusion and permeation rate and the hydrolysis reaction rate based on the solution results of each sub-equation. The core is to form a bidirectional real-time feedback link through "forward transmission of boundary conditions - reverse transmission of correction parameters," while using a field coupling controller to dynamically adjust the solution step size to ensure complete synchronization between the simulation timing and the sensor data sampling frequency. The specific implementation method is as follows:

[0078] First, the core parameters required for real-time adjustment are extracted using the solution results of the operator splitting method. This is a fundamental prerequisite for achieving dynamic rate adjustment. The decoupling characteristic of the operator splitting method ensures that the solution results of each sub-equation are both independent and complete, yet interconnected, providing a clear data source for parameter extraction. The specific process is as follows:

[0079] Within each main computation step [t_n, t_{n+1}], after completing independent calculations in the order of solving the "heat transfer sub-equation → proton transfer equation → reaction kinetic sub-equation", the system extracts three types of key results from the shared data pool of the field coupling controller: First, the global three-dimensional temperature field data output from solving the heat transfer sub-equation, including the real-time temperature values ​​of each grid cell of the container wall and the liquid phase bulk. This data will serve as the environmental basis parameters for subsequent dynamic correction of the diffusion coefficient and reaction rate constant. Second, the key data on water molecule diffusion and permeation output from solving the proton transfer equation. The core data is the instantaneous water flux (i.e., water molecule diffusion and permeation rate) of each grid cell on the inner surface of the container wall, and the water concentration distribution field of the entire liquid phase bulk. The former is the core parameter for reactant supply in the reaction process, and the latter is the direct input for reaction rate calculation. Third, the hydrolysis reaction-related data output from solving the reaction kinetic sub-equation, including the real-time hydrolysis reaction rate of each grid cell, the total product concentration distribution, and the product inhibition coefficient calculated based on the product concentration. These data are the core basis for the reverse correction of the mass transfer process.

[0080] During the extraction process, the system verifies the validity of all data, removing outliers caused by grid distortion and numerical oscillations to ensure the accuracy and reliability of the extracted parameters. Simultaneously, the verified valid data is indexed and bound according to grid cell coordinates, forming a one-to-one correspondence between coordinates and parameters, providing a spatial positioning basis for subsequent forward and reverse transmission. After completing the core parameter extraction, forward parameter transmission is first executed, using the water flux output from the proton transfer equation as the boundary condition update value for the reaction kinetics sub-equation, enabling the mass transfer process to drive the chemical reaction process in real time. The specific process is as follows: First, clarify the type of water boundary condition in the reaction kinetics sub-equation. Since water molecules are the core reactant in the hydrolysis reaction, their supply rate directly determines the extent of the reaction. Therefore, the water input boundary of the reaction sub-equation is set as a second-type Neumann boundary condition (flux boundary). The core of this boundary condition is to dynamically control the supply intensity of reactants by specifying the water flux entering the reaction system per unit time. Second, perform spatial matching and transfer of water flux. The system transfers the instantaneous water flux of each grid cell on the inner surface of the container wall extracted by the proton transfer equation to the boundary of the corresponding liquid phase grid cell in the reaction kinetics sub-equation one-to-one according to the coordinate index, ensuring that the water supply flux of each reaction region is completely consistent with the permeation result of the mass transfer process, avoiding spatial mismatch. The third step involves adjusting the hydrolysis reaction rate based on the updated boundary conditions. After receiving water flux data, the reaction kinetics sub-equation first calculates the total amount of water molecules entering the grid cell per unit time based on the grid cell volume. Then, combined with the current water concentration in the liquid phase, it dynamically adjusts the reaction rate constant through a microscopic reaction probability model. When the water flux increases, the number of water molecules participating in the reaction per unit time increases, and the probability of the reaction increases. The system will synchronously increase the reaction rate constant at a 1:1 ratio to the flux increase. When the water flux decreases, the water molecule supply is insufficient, and the reaction probability decreases. The system will synchronously decrease the reaction rate constant, ultimately achieving dynamic adaptation of the hydrolysis reaction rate to the water molecule diffusion and permeation rate, allowing the reaction process to truly reflect changes in water permeation supply. After the reaction rate adjustment is completed in the forward transfer, the reverse parameter backhaul is executed synchronously. The product inhibition coefficient generated by the reaction kinetics sub-equation is backhauled to the proton transfer equation to correct the water molecule diffusion coefficient, realizing the reverse constraint of the chemical reaction process on the mass transfer process. The specific process is as follows:

[0081] The product inhibition coefficient is a dimensionless parameter calculated based on the total product concentration distribution obtained from the reactant equation. Its value ranges from 0.1 to 1.0; the higher the product concentration, the smaller the inhibition coefficient. Its core function is to quantify the hindering effect of product accumulation on water diffusion. Product aggregation in the liquid phase alters the system's viscosity and polarity, increasing the diffusion resistance of water molecules in the liquid phase. Simultaneously, product adsorption on the inner surface of the container wall blocks some pores, reducing water permeability. The product inhibition coefficient of each grid cell output from the reactant equation is fed back to the corresponding grid cell of the proton transfer equation according to its coordinate index. The inhibition coefficient in the liquid phase region is directly related to the water diffusion coefficient within the liquid phase, while the inhibition coefficient in the container wall region is related to the water diffusion coefficient within the wall material, ensuring accurate coverage of the inhibition effect across the entire mass transfer region. The proton transfer equation receives the inhibition coefficient... After the initial diffusion coefficient is calculated, it is coupled with the original diffusion coefficient tensor to form the corrected dynamic diffusion coefficient. The correction logic is "corrected diffusion coefficient = original diffusion coefficient × product inhibition coefficient". The original diffusion coefficient is a basic value calculated based on the porous medium permeability tensor, diffusion barrier function and current temperature field. By multiplying it with the product inhibition coefficient, the dynamic change of the diffusion coefficient with the product concentration is realized. When the product concentration increases and the inhibition coefficient drops to 0.5, the diffusion coefficient drops to 50% of the original value, and the water diffusion and permeation rate slows down. When the product concentration is low and the inhibition coefficient is close to 1.0, the diffusion coefficient basically remains at the original value, and water diffusion is not significantly affected. Through this correction process, the product accumulation is used to reverse the constraint on water diffusion, so that the mass transfer process truly reflects the aggregation effect of the reaction products.

[0082] Through forward transmission and reverse feedback, a bidirectional real-time feedback loop has been formed, consisting of "moisture diffusion and permeation → hydrolysis reaction → product accumulation → diffusion coefficient correction → moisture diffusion and permeation adjustment". To ensure the stable operation and simulation accuracy of this loop, the solution step size of each sub-equation needs to be dynamically adjusted by the field coupling controller to perfectly match the sensor data sampling frequency. The specific process is as follows: First, clarify the benchmark for step size adjustment. The sensor data sampling frequency is set to 1 hour / time, which is the basic reference for the main simulation step size. The step size adjustment of the field coupling controller must not deviate from this benchmark frequency, ensuring that the latest sensor data is synchronized within each main step size to avoid decoupling between the simulation and the actual environment. Second, ensure... The core judgment indicators for fixed-step size adjustment are monitored in real time by the controller in two key dimensions: First, the convergence of solving each sub-equation. The convergence speed is judged by calculating the rate of change of the residual of each sub-equation (the ratio of the current residual to the residual of the previous step). A rate of change of the residual below 1×10^-3 indicates fast convergence, while a rate of change above 1×10^-1 indicates slow convergence. Second, the gradient of physical quantity changes, including the spatial gradient and time rate of change of the temperature field, moisture concentration field, and product concentration field. A gradient value higher than the set threshold (temperature gradient > 2℃ / cm, concentration gradient > 0.1mol / (L·cm)) indicates drastic changes in physical quantities, while a value lower than the threshold indicates gradual changes. The third step is to perform dynamic step size adjustment, based on the controller's... The above-mentioned criteria allow for independent adjustment of the sub-step size for each sub-equation: For sub-equations with fast convergence and gradual changes in physical quantities (such as heat transfer sub-equations without significant temperature fluctuations), the sub-step size is appropriately increased within the range of the main step size (maximum not exceeding 1 / 2 of the main step size) to reduce the number of iterations and improve computational efficiency; for sub-equations with slow convergence and drastic changes in physical quantities (such as reaction kinetics sub-equations when products accumulate rapidly, and proton transfer sub-equations when water flux changes abruptly), the sub-step size is reduced (minimum not less than 1 / 10 of the main step size) to increase the number of iterations and ensure that subtle changes in physical quantities are captured; during the adjustment process, the controller strictly ensures that the sum of the sub-step sizes of the three sub-equations equals the main step size. The first step is to avoid timing misalignment by setting a maximum step size adjustment limit (no more than 50% of the original step size in a single adjustment) to prevent numerical oscillations and calculation divergence caused by sudden changes in step size. The fourth step is to achieve synchronous calibration with the sensor sampling frequency. After each main step, the controller will compare the current solution step size configuration with the sensor sampling period. If the simulation timing deviates from the sensor data sampling time due to step size adjustment (the deviation exceeds 1 minute), the sub-step size allocation will be fine-tuned in the next main step to ensure that the simulation time node is completely aligned with the sensor data acquisition time, so that each rate adjustment can be based on the latest actual environmental data and achieve real-time synchronization between the simulation and the actual scene.

[0083] Through the aforementioned forward transfer, reverse feedback, and step size adjustment throughout the entire process, a bidirectional real-time feedback link between heat transfer, mass transfer, and chemical reaction is formed. This ensures that the water molecule diffusion and permeation rate and the hydrolysis reaction rate can be dynamically adjusted according to environmental changes and material states, accurately simulating the dynamic balance of the moisture gradient and the product accumulation effect in a closed container. Furthermore, the dynamic optimization of the step size ensures the computational efficiency and convergence stability of the simulation. This allows the entire simulation process to not only conform to actual physical laws but also efficiently respond to real-time environmental data, providing a reliable dynamic adjustment mechanism for the accurate output of the subsequent stability decay risk index.

[0084] After dynamically adjusting the water molecule diffusion and permeation rate and hydrolysis reaction rate through a bidirectional real-time feedback link, and achieving precise synchronization between the solution step size and sensor data sampling frequency via a field coupling controller, the dynamic multiphysics real-time feedback module focuses on the visualization and quantitative evaluation of the core physicochemical processes inside the sealed container. Based on the global grid-level data output by the link, through the progressive logic of three-dimensional moisture distribution field reconstruction, product aggregation region prediction, and risk index mapping, it fully simulates the dynamic balance of the moisture gradient and the product accumulation effect, providing an intuitive and accurate quantitative basis for subsequent stability assessment. The specific implementation method is as follows:

[0085] First, based on the output data of the bidirectional real-time feedback link, a spatial interpolation algorithm is used to reconstruct the three-dimensional moisture distribution field inside the container. The output data of the bidirectional real-time feedback link includes the instantaneous moisture concentration value of each finite element mesh element of the container wall and the liquid phase body. However, due to the limitation of mesh density (especially the relatively sparse mesh in the central region of the liquid phase), the directly output data is discretely distributed and cannot fully present the continuous migration law and gradient change of moisture in the container. The spatial interpolation algorithm is a numerical method that calculates the values ​​of adjacent unknown points by knowing the values ​​of discrete points. The radial basis function interpolation method is selected here because it has a good fitting ability for nonlinear distribution data, can accurately restore the smooth distribution characteristics of moisture concentration in three-dimensional space, and the interpolation error is controllable. The specific reconstruction process involves four core steps: data preprocessing, interpolation model construction, global field generation, and accuracy verification. The first step is data preprocessing. The system extracts water concentration data from each grid cell from a shared data pool, sorts it by container 3D coordinates, and removes outliers caused by numerical oscillations (such as concentrations exceeding the 0-saturation range). Simultaneously, boundary conditions are calibrated—the outer surface of the container wall uses equilibrium concentration correction calculated from GAB adsorption isotherms, while the inner surface uses instantaneous flux from the proton transfer equation to back-calculate concentration values. At the interface between the liquid phase and the container wall, concentration continuity is ensured to avoid interpolation distortion caused by abrupt boundary changes. The second step is interpolation model construction. The interpolation node density is determined to be one 3D interpolation point every 0.5 cm (this density has been experimentally verified to ensure detail restoration without excessive computational burden). A Gaussian radial basis function is selected as the interpolation kernel function, and its shape parameters are determined through cross-validation (using known grid data). The grid points are randomly divided into training and validation sets. The shape parameters are iteratively adjusted until the interpolation error of the validation set is minimized, ultimately constructing an interpolation model that fits the ethylene carbonate storage scenario. The third step is global field generation. The preprocessed known grid point concentration data is input into the interpolation model. The moisture concentration values ​​of all unknown interpolation points are calculated through model fitting. Then, according to the three-dimensional physical dimensions of the container, the concentration data of all known points and interpolation points are integrated to generate a continuous and smooth three-dimensional moisture distribution field. This field can clearly show the concentration gradient change of moisture permeating from the container wall to the center of the liquid phase, as well as the hindering effect of product aggregation areas on moisture diffusion (manifested as an increase in the moisture concentration gradient around the aggregation area). The fourth step is interpolation accuracy verification. The relative error between the interpolation results of known grid points and the original data is calculated. The average relative error is required to be less than 5%. If the error in a certain area exceeds the standard, the grid in that area is densified and re-interpolated to ensure that the reconstructed three-dimensional moisture distribution field is accurate and reliable.

[0086] Based on this, using the reconstructed continuous moisture distribution field, the moisture gradient (i.e., the rate of change of concentration along the three spatial coordinate axes) at each interpolation point is calculated by solving the spatial partial derivatives, resulting in a three-dimensional moisture gradient distribution field. This allows for the identification of extreme moisture gradient regions—typically concentrated on the inner surface of the container wall and at the boundaries of product aggregation regions. These areas represent the core regions where moisture permeation and reaction coupling are most intense and are the focus of subsequent risk assessments. After reconstructing the three-dimensional moisture distribution field, the product aggregation regions at different time points are predicted using the product transport equation. The product transport equation is well-defined; its core function is to describe the diffusion, convection, and aggregation processes of products in the liquid phase. By combining the reconstructed moisture distribution field with real-time updated reaction rate data, the spatial distribution patterns of products at different storage stages can be accurately predicted. The specific prediction process is detailed as follows: First, key time nodes are determined. Based on the typical storage period of ethylene carbonate (usually 30 days), 1 day, 3 days, 7 days, 15 days, and 30 days are selected as core time nodes. User-defined special time nodes are also supported to ensure coverage of both short-term rapid reactions and long-term slow accumulation. Second, for each time node, the three-dimensional moisture distribution field, global temperature field, and hydrolysis reaction rate distribution data at that moment are input. The diffusion coefficient (which varies with temperature and product concentration), convection term (natural convection caused by temperature gradient), and reaction source term (product generation amount determined by the hydrolysis reaction rate) in the product transport equation are dynamically updated. The transport equation is solved using the finite volume method to obtain the product concentration value of each grid cell at that time node. Then, the product aggregation threshold is determined. This threshold is calibrated experimentally—based on ethylene carbonate storage failure case data, when the product concentration exceeds 0.5 mol / L... When the concentration of product exceeds the aggregation threshold, it will significantly accelerate material corrosion and quality degradation. Therefore, 0.5 mol / L is set as the aggregation threshold, and the threshold can be adjusted according to the product type (such as ethylene glycol, organic acid). Then, aggregation region identification is performed, traversing all grid cells and marking cells with product concentrations exceeding the aggregation threshold. Then, through spatial connectivity analysis (using the 8-neighborhood connectivity criterion), adjacent aggregation grid cells are merged into continuous product aggregation regions, and the spatial coordinates, volume, average concentration and maximum concentration of each aggregation region are recorded. Finally, the prediction results are output. For each time node, a three-dimensional spatial distribution map of the product aggregation region is generated, which clarifies the location (such as near the container wall, liquid phase center or local hot spot), morphology and concentration change trend of the aggregation region. At the same time, the spatial correlation between the aggregation region and the extreme value region of the moisture gradient is analyzed. Usually, the product aggregation region and the extreme value region of the moisture gradient highly overlap, because the reaction is more active in the region with intense water penetration, and the product is more likely to accumulate.After completing the prediction of moisture gradient distribution and product aggregation region, the extreme values ​​of moisture gradient and product aggregation concentration are mapped to a stability decay risk index by normalization method. This index is a comprehensive quantitative assessment of the storage stability of ethylene carbonate, including two core dimensions: moisture penetration risk component and product corrosion risk component. The value range is 0-1, and the larger the index, the higher the corresponding risk. The specific mapping process is as follows: First, quantitative indicators are extracted. The maximum value of the global moisture gradient (i.e., the extreme value of the moisture gradient) is extracted from the three-dimensional moisture gradient distribution field. This value directly reflects the severity of moisture penetration. The maximum product concentration value of all aggregation regions is extracted from the product aggregation region prediction results. This value reflects the severity of product corrosion. Second, normalized benchmark values ​​are determined. The normalized benchmark for moisture penetration risk is the critical value of the moisture gradient that leads to moisture exceeding the standard and failure in historical stored data (experimentally calibrated to 5 mol / (L·cm)). The normalized benchmark for product corrosion risk is the product failure concentration that leads to excessive acid value and color change of ethylene carbonate (experimentally calibrated to 1.0 mol / L). Both benchmark values ​​are dynamically adjusted according to changes in container material and storage conditions. Third, a single risk component is calculated: Moisture penetration risk component = current extreme value of moisture gradient / critical benchmark value of moisture gradient. If the calculated result is greater than 1, it is taken as 1 (indicating that the critical state of permeation failure has been reached), and if it is less than 0, it is taken as 0 (indicating no permeation risk). The product corrosion risk component = the current maximum concentration of the product / the baseline value of the product failure concentration. Similarly, when the result exceeds the 0-1 range, it is treated as a boundary value. The fourth step is to calculate the comprehensive risk index. The two components are integrated by weighted summation. The weights are calibrated by multiple sets of storage experiments. The weight of the moisture permeation risk component is 0.4, and the weight of the product corrosion risk component is 0.6. Since the product accumulation has a more direct and significant impact on the quality degradation of ethylene carbonate, the final stability degradation risk index = moisture permeation risk component × 0.4 + product corrosion risk component × 0.6. The value range of this index is 0-1, where 0-0.3 indicates low risk (stable storage), 0.3-0.7 indicates medium risk (changes need to be monitored), and 0.7-1.0 indicates high risk (approaching or reaching the failure state). Through the entire process of reconstructing the three-dimensional moisture distribution field, predicting the product aggregation region, and mapping the risk index, we not only fully simulated the entire process of moisture in a closed container from infiltration and diffusion to dynamic equilibrium, as well as the entire process of product generation, transport, and aggregation, but also transformed the abstract physicochemical process into an intuitive spatial distribution and a quantified risk index. This provides core input data for predicting the evolution trend of subsequent key quality indicators and provides accurate decision-making basis for optimizing storage strategies (such as adjusting storage temperature and changing container materials), thus realizing a closed-loop connection from multi-physics field coupling simulation to stability quantitative assessment.

[0087] After the dynamic multiphysics real-time feedback module outputs a stability decay risk index containing both moisture penetration risk and product corrosion risk components, in order to transform the abstract risk quantification value into an intuitive quality control basis, it is necessary to predict the evolution trend of three key quality indicators—acid value, moisture content, and color—based on this index, and accurately identify the failure boundary of storage stability. The entire process revolves around the correlation logic between risk components and quality indicators, uses time series extrapolation to achieve trend prediction, relies on the abrupt change characteristics of the risk index to locate the failure threshold, and finally forms a feasible composite failure criterion. The specific implementation method is as follows:

[0088] First, a correlation model is established between the stability degradation risk component and the rate of change of key quality indicators. This is the core premise for trend prediction. The correlation logic is based on the hydrolysis reaction mechanism and quality degradation law of ethylene carbonate: the moisture penetration risk component directly reflects the intensity of external moisture penetration into the container. Moisture is the core reactant of the hydrolysis reaction. The more intense the penetration, the more acidic products are generated by the hydrolysis reaction, and the moisture content of the liquid phase will also increase simultaneously. Therefore, a positive correlation is established between the moisture penetration risk component and the rate of change of acid value and the rate of change of moisture content. The product corrosion risk component characterizes the accumulation degree of hydrolysis products (ethylene glycol, organic acids, etc.). The higher the product concentration, the easier it is to trigger side reactions such as oxidation and polymerization of ethylene carbonate, resulting in a darker color of the material (e.g., from colorless and transparent to pale yellow or yellow). Therefore, a positive correlation is established between the product corrosion risk component and the rate of change of color. The specific correlation process is as follows:

[0089] By collecting historical data under different storage conditions over the past 12 months, including the stability degradation risk index, corresponding acid value, moisture content, and color detection values ​​at each time point, a correlation model was constructed using a multiple linear regression algorithm. The calculation expressions for the change rates of the three quality indicators were obtained as follows: Acid value change rate = 0.8 × moisture penetration risk component + 0.05 (where the coefficients were obtained by fitting historical data, R² ≥ 0.92, ensuring the reliability of the correlation); Moisture content change rate = 1.2 × moisture penetration risk component - 0.03 (after coefficient fitting, R² ≥ 0.95, considering the balance between moisture evaporation and penetration); Color change rate = 0.9 × product corrosion risk component + 0.02 (after coefficient fitting, R² ≥ 0.95).88. (Covering the color change patterns under different product concentrations). After the correlation model is established, iterative optimization will be performed quarterly using newly added stored experimental data to correct the regression coefficients and ensure that the correlation always closely matches the actual decay pattern. After the correlation model is completed, the evolution curves of key quality indicators are generated using time series extrapolation. Time series extrapolation refers to a prediction method that analyzes the trend, periodicity, and stability of changes in historical data sequences of things over time, fits historical patterns through mathematical models, and then extrapolates the trend of data changes in the future. This method does not rely on additional external variables and can achieve prediction solely through the inherent patterns of its own time series, making it particularly suitable for ethylene carbonate storage. In scenarios involving gradual changes and relatively stable trends in storage stability, the system can accurately capture the slow evolution of quality indicators. The specific implementation process follows these steps: First, a historical time-series dataset is constructed. The system extracts existing time-node data within the current storage period from the digital twin's simulated database, arranging them chronologically to form a complete historical sequence—the time nodes align with the key nodes predicted for product aggregation areas, including stability decay risk indices (containing two components) for 1 day, 3 days, 7 days, 15 days, and 30 days. Then, the acid value change rate, moisture content change rate, and color change rate corresponding to each node are calculated using the aforementioned correlation model. The rate of change in color, combined with the initial detection values ​​of each indicator (measured values ​​at the start of storage), is used to derive the actual values ​​of the quality indicators at each node, ultimately forming three historical time series: "time-acid value," "time-moisture content," and "time-color." The second step involves time series preprocessing. Each series undergoes a stationarity test (using the ADF test). If a series exhibits non-stationarity (e.g., trend shifts due to minor fluctuations in storage conditions), it is transformed into a stationary series using a first-order difference transformation. Simultaneously, outliers caused by sensor malfunctions or simulation errors are removed (using the 3σ criterion, i.e., removing data points deviating from the series mean by three times the standard deviation), ensuring... The reliability of the sequence data; the third step is to select a suitable extrapolation model and calibrate the parameters. Considering that the evolution trend of the storage stability index of ethylene carbonate is mainly gradual and without obvious periodicity, the quadratic exponential smoothing method is selected as the core extrapolation model. This model can effectively capture the linear trend of the sequence and balance prediction accuracy and computational efficiency. The model parameter calibration adopts the historical data back-substitution verification method, using the first 80% of the historical data as the training set and the last 20% as the validation set. The smoothing coefficient α (within the range of 0.1-0.3) is optimized by the least squares method so that the average relative error between the predicted value and the actual value of the validation set is less than 5%. Finally, the value of α is determined (usually calibrated to 0).2) The fourth step is to perform future trend extrapolation prediction. Based on the calibrated quadratic exponential smoothing model, and using existing historical time series as a basis, the key time node quality indicator values ​​for the next 60 days are extrapolated. The extrapolation time nodes are set at 5 days, covering short-term (1-30 days) and medium-to-long-term (31-60 days) storage periods to ensure complete trend coverage. The fifth step is to generate indicator evolution curves. Historical actual values ​​and future predicted values ​​are integrated in chronological order. Each quality indicator corresponds to an evolution curve. Different colors are used to distinguish historical data segments from predicted data segments in the curves. At the same time, the prediction interval (based on the 95% confidence interval of the model prediction error) is marked to intuitively present the changing trend of the indicators. For example, the acid value evolution curve gradually increases over time, the moisture content curve rises rapidly and then flattens out, and the color curve deepens slowly. This allows staff to clearly understand the quality change patterns at different stages. After obtaining the key quality indicator evolution curves, the stability failure threshold is identified based on the abrupt inflection point of the stability decay risk index, and then integrated into a composite failure criterion. The specific implementation process is as follows:

[0090] The first step is to locate the abrupt inflection point of the stability decay risk index. An abrupt inflection point refers to a significant, non-gradual change in the risk index within a continuous time period. This reflects the hydrolysis reaction transitioning from a gradual to a rapid phase, or the accumulation of products reaching a critical state. The specific identification method is as follows: using a sliding window method (with a window size set to three consecutive time points), calculate the average rate of change of the risk index within each window, and set a mutation threshold (calibrated to 0.2 / day based on historical failure cases). When the average rate of change in a certain window exceeds this threshold, the time point in the middle of the window is determined to be the abrupt inflection point. Simultaneously, the evolution curves of three quality indicators are used to verify whether the rate of change of the indicators at the inflection point changes synchronously, ensuring the accuracy of inflection point identification and avoiding misjudgments due to fluctuations in a single data point. The second step is to locate the three major failure thresholds based on the abrupt inflection point. For the acid value exceeding the standard threshold, the predicted acid value corresponding to the abrupt inflection point is queried, while also referring to the upper limit of the acid value of ethylene carbonate in the industry standard. For example, according to GB / T19281-2014, the acid value of superior grade products should be ≤0.01mgKOH / g. The smaller of the two values ​​is taken as the critical point for exceeding the acid value standard. That is, when the acid value reaches this value, it is judged as a quality failure. For the moisture saturation concentration point, the rate of change of moisture content will decrease significantly at the abrupt inflection point (because moisture penetration and hydrolysis consumption reach equilibrium). The moisture content corresponding to this point is the moisture saturation concentration point. This value is the upper limit threshold of moisture content. Exceeding it will accelerate the hydrolysis side reaction. For the color change threshold, the color change is quantified by the color difference meter detection standard (such as using the ΔE value of the CIELab color space). The ΔE value corresponding to the abrupt inflection point is the color change threshold. When the ΔE value exceeds this threshold, the appearance quality of the material decreases significantly and is considered a failure. The third step is to integrate the composite failure criteria. Considering that the failure of any one of the three quality indicators will affect the performance of ethylene carbonate, the composite failure criteria are integrated using "OR logic".

[0091] When the acid value reaches the critical point of exceeding the standard, the moisture content reaches the saturation concentration point, or the color reaches the abrupt change threshold, the storage stability of ethylene carbonate is determined to have failed, and an early warning mechanism must be activated. To ensure the flexibility of the criteria, the weights can be adjusted according to the application scenario. For example, the acid value exceeding the standard critical point can be set to a more stringent value for high-purity electronic-grade ethylene carbonate, while the moisture saturation concentration point can be appropriately relaxed for industrial-grade products, allowing the criteria to adapt to different quality requirements. Through the correlation between the aforementioned risk components and quality indicators, trend prediction using time series extrapolation, and failure threshold identification driven by abrupt change inflection points, not only is the abstract risk index transformed into an intuitive quality evolution curve, providing a visual basis for quality control in the storage process, but the composite failure criteria also clarify the boundaries of stability, solving the problem of "judgment based on experience and delayed early warning" in traditional storage management. This achieves accurate prediction and early control of the storage stability of ethylene carbonate, providing a core failure judgment standard for the subsequent parameter iteration and correction of the digital twin basic model.

[0092] After completing the prediction of the evolution trend of key quality indicators and the identification of stability failure thresholds, in order to further improve the simulation accuracy of the digital twin basic model and make it more consistent with the quality degradation law of actual storage scenarios, it is necessary to iteratively correct the core parameters of the model by comparing and analyzing the detection data of actual storage samples with the model's predicted values, forming a closed loop of "simulation-verification-correction-optimization". The specific implementation method is as follows:

[0093] First, time-series alignment of actual test data with model predictions is crucial for ensuring the effectiveness of deviation calculations. The core of time-series alignment is to ensure a one-to-one correspondence between actual test data and model predictions across time, eliminating misjudgments caused by time misalignment. The specific process is as follows: First, collect the test data sequence of the actual stored samples. Under the same storage conditions as the digital twin simulation, conduct laboratory tests on the stored samples for acid value, moisture content, and color at set key time nodes (1 day, 3 days, 7 days, 15 days, 30 days) and newly added failure threshold proximity nodes (e.g., the three time points before the acid value approaches the critical point of exceeding the standard). The testing methods strictly adhere to industry standards (acid value is measured using potentiometric titration, moisture content using the Karl Fischer method, and color using a CIELab colorimeter) to ensure the accuracy of the test data. Second, organize the model prediction value sequence from the digital twin... The first step involves extracting predicted acid value, moisture content, and color values ​​for corresponding time points from the simulation database of the basic model, forming a predicted value sequence consistent with the actual detection data structure. The second step is to perform time-series alignment, using the actual detection timestamp as a benchmark to calibrate the model's predicted value sequence. If there is a slight deviation (≤2 hours) between a certain actual detection time point and the model's predicted time point, linear interpolation is used to correct the model's predicted value, ensuring that each actual detection data point can find a unique corresponding model predicted value. If actual detection data is missing (e.g., a time point was not detected), it is filled by trend fitting between two adjacent valid detection data points to avoid alignment interruptions due to data loss. Finally, an aligned dataset of "timestamp - actual acid value - predicted acid value - actual moisture content - predicted moisture content - actual color - predicted color" is formed, providing a unified time-series basis for subsequent deviation analysis.After time-series alignment, prediction deviation weights are calculated using integrated composite failure criteria. The core function of these weights is to quantify the impact of prediction deviations at different time points on model accuracy, giving higher weights to critical deviations approaching the failure threshold. This ensures the correction process focuses on core influencing factors. The specific calculation process is as follows: First, calculate the prediction deviation for each single indicator. For acid value, moisture content, and color at each time point, calculate the absolute deviation: "Single indicator deviation = |Actual detected value - Model predicted value|". Simultaneously, calculate the relative deviation: "Relative deviation = Single indicator deviation / Actual detected value", avoiding distortion in deviation assessment due to differences in indicator magnitude. Second, introduce composite failure criteria for weight allocation. These criteria include the acid value exceeding the critical point, the moisture saturation concentration point, and the color change threshold. The system first determines the distance between the actual detected value at each time point and the corresponding failure threshold, defining "Safety margin = (Failure threshold - Actual detected value) / Failure threshold". (If the actual value has exceeded the failure threshold, the safety margin is set to 0); then, a weight coefficient is assigned to each single indicator deviation. The smaller the safety margin (i.e., the closer to the failure threshold), the larger the weight coefficient. Specifically, the weight coefficient is set as follows: when the safety margin ≤ 0.2, the weight coefficient = 1.5; when 0.2 < safety margin ≤ 0.5, the weight coefficient = 1.2; when the safety margin > 0.5, the weight coefficient = 0.8. This setting can highlight the prediction deviation in the critical failure region and allow the model to prioritize the correction of parameters that are crucial to the stability judgment; the third step is to calculate the comprehensive prediction deviation weight. For the relative deviation of the three single indicators at each time point, multiply them by the corresponding weight coefficient, and then sum them according to the weighted ratio of "acid value deviation weight 40%, moisture content deviation weight 35%, color deviation weight 25%" to obtain the comprehensive prediction deviation weight of that time point. The weight value ranges from 0 to 1.5. The larger the value, the more significant the impact of the prediction deviation of that point on the model accuracy.Next, a Long Short-Term Memory (LSTM) neural network is used to analyze the correlation between the comprehensive prediction bias weights and the water molecule diffusion and permeation rate coefficient and the hydrolysis reaction rate constant. LSTM, due to its core advantage of capturing long-term dependencies in time-series data, can accurately uncover the intrinsic correlation between the changes in bias weights over time and the core parameters of the model, avoiding nonlinear relationships that traditional linear models cannot fit. The specific implementation process is as follows: First, construct the input and output datasets of the LSTM model. The input data is the time-aligned sequence of comprehensive prediction bias weights, divided into sliding windows of length 5 (covering 5 consecutive time nodes) according to time order. Each window corresponds to one... The first step involves inputting samples to ensure the model can capture the temporal trend of bias changes. The output data consists of the actual effective values ​​of the water molecule diffusion and permeation rate coefficient and the hydrolysis reaction rate constant within the corresponding time window (obtained by back-calculation of the output data from the bidirectional real-time feedback link), forming a one-to-one correspondence dataset of "input bias weight window - output core parameters". The second step involves dividing the dataset into training and validation sets, splitting it in a 7:3 ratio, with 70% used as the training set for model parameter training and 30% as the validation set for model accuracy verification. Simultaneously, the input and output data are normalized (mapped to the 0-1 range) to avoid model training errors due to differences in numerical magnitude. Imbalance; the third step is to configure the LSTM model structure and training parameters. The model adopts a three-layer architecture of "input layer - hidden layer - output layer". The number of neurons in the input layer is the same as the sliding window length (5). There are 2 hidden layers, each containing 64 neurons. The activation function is the ReLU function. The output layer contains 2 neurons (corresponding to the diffusion rate coefficient and the correlation coefficient of the reaction rate constant, respectively). The training process uses the Adam optimizer, with an initial learning rate of 0.001, which decays by 50% every 20 iterations. The maximum number of iterations is 200. The loss function is the mean square error between the predicted and actual values. When the loss changes after 10 consecutive iterations... The value is considered convergent when it is less than 1×10^-4; the fourth step is model training and correlation analysis. The training set is input into the LSTM model for iterative training. The model weight parameters are adjusted through backpropagation. After training, the model accuracy is verified using the validation set. The average relative error of the validation set is required to be less than 8%; finally, the new comprehensive prediction bias weight sequence is input into the trained LSTM model, and the correlation coefficient between the two core parameters and the bias weight is output (the value range is -1 to 1, a positive value indicates that the parameter is too high and causes the bias, a negative value indicates that the parameter is too low and causes the bias, and the larger the absolute value, the stronger the correlation). The core parameters that cause the prediction bias and the direction of influence are accurately located.After clarifying the correlation between core parameters and prediction bias, parameter correction amounts are generated and backpropagated to the digital twin base model to achieve closed-loop optimization. The specific process is as follows: First, calculate the parameter correction amount. Using the water molecule diffusion and permeation rate coefficient and hydrolysis reaction rate constant in the current digital twin base model as benchmark values, and combining the correlation coefficient output by the LSTM with the average of the comprehensive prediction bias weights, calculate according to the logic of "correction amount = benchmark value × correlation coefficient × average bias weight". For example, if the correlation coefficient of the diffusion rate coefficient is 0.3 (indicating that the parameter is too high), the average bias weight will be adjusted accordingly. If the average deviation weight is 0.6, then the correction amount = baseline value × 0.3 × 0.6 = baseline value × 0.18, meaning the diffusion rate coefficient needs to be reduced by 18%. If the correlation coefficient of the reaction rate constant is -0.2 (indicating the parameter is too low), and the average deviation weight is 0.5, then the correction amount = baseline value × (-0.2) × 0.5 = - baseline value × 0.1, meaning the reaction rate constant needs to be increased by 10%. Simultaneously, upper and lower limit thresholds are set for the correction amount (the correction range should not exceed 30% of the baseline value) to avoid excessive parameter correction due to misjudgment of correlation. The second step is to execute parameter reverse transmission. The first step involves updating the model by feeding the calculated correction back into the corresponding parameter module of the digital twin's base model. This updates the diffusion coefficient tensor in the finite element space discretization model and the rate constant in the hydrolysis reaction rate equation. Simultaneously, it updates the calculation logic in the dynamic water adsorption and reaction kinetics module and the dynamic multiphysics real-time feedback module that depend on these parameters, ensuring the consistency of parameters across the entire model. The second step verifies the correction effect by reusing the updated model for simulation calculations. This generates new predicted values ​​for key quality indicators, which are then time-aligned and their deviations calculated again with the actual detection data. If the average relative deviation after correction is reduced by more than 30% compared to before correction, the correction is considered effective, and the updated parameters are retained. If the deviation is not significantly reduced, the LSTM model is retrained, and the model structure (such as increasing the number of hidden layer neurons) or weight allocation ratio is adjusted. The correction is then generated again for iteration until the average relative deviation of the model prediction is below 5% (meeting the accuracy requirements for engineering applications). This complete closed-loop optimization process ensures that the simulation capability of the digital twin's base model continuously improves with the accumulation of actual data, guaranteeing the accuracy and reliability of the simulation results for subsequent storage stability.

[0094] The foregoing has shown and described the basic principles, main features, and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The embodiments and descriptions in the specification are merely preferred examples and are not intended to limit the invention. Various changes and modifications can be made to the invention without departing from its spirit and scope, and all such changes and modifications fall within the scope of the present invention as claimed. The scope of protection of the present invention is defined by the appended claims and their equivalents.

Claims

1. A method for simulating the storage stability of ethylene carbonate based on digital twins, characterized in that: Includes the following steps: S1. Construct a digital twin model of the molecular structure and hydrolysis reaction kinetics of ethylene carbonate; S2. In the digital twin basic model, a dynamic water adsorption and reaction kinetics module designed based on the easy hydrolysis characteristics of ethylene carbonate is embedded. The dynamic water adsorption and reaction kinetics module calculates the water molecule diffusion and permeation rate, hydrolysis reaction rate and product concentration distribution in real time. S3. Based on the water molecule diffusion and permeation rate, hydrolysis reaction rate and product concentration distribution, a dynamic multiphysics field real-time feedback module is designed. The dynamic multiphysics field real-time feedback module integrates temperature and relative humidity sensor data streams, and adjusts the water molecule diffusion and permeation rate and hydrolysis reaction rate in real time by coupling heat transfer and mass transfer and chemical reaction kinetic equations to simulate the dynamic balance of the water gradient and product accumulation effect in a closed container, and outputs a stability decay risk index. S4. Using the stability decay risk index, a time series is set to predict the evolution trend of key quality indicators of ethylene carbonate under different storage conditions. The key quality indicators include acid value, moisture content and color, and the stability failure threshold is identified. S5. By applying machine learning algorithms to compare the detection data of the actual stored samples, the predicted values ​​of the digital twin basic model, and the stability failure threshold, the parameters of the digital twin basic model are iteratively corrected, and the visualized simulation results are output.

2. The method for simulating the storage stability of ethylene carbonate based on digital twins according to claim 1, characterized in that: In constructing the digital twin basic model of ethylene carbonate molecular structure and hydrolysis reaction kinetic parameters, quantum chemical calculations were used to determine the molecular conformation and charge distribution of ethylene carbonate. Molecular dynamics simulations were combined to obtain the transition state energy and reaction barrier in the hydrolysis reaction path. The transition state energy and reaction barrier were then input into a spatial discretization model based on a finite element mesh. The spatial discretization model used adaptive mesh refinement technology to capture the influence of intermolecular forces on the hydrolysis reaction rate, thus forming the basic framework of the digital twin model.

3. The method for simulating the storage stability of ethylene carbonate based on digital twins according to claim 1, characterized in that: The dynamic water adsorption and reaction kinetics module designed based on the easy hydrolysis characteristics of ethylene carbonate specifically introduces the permeability tensor of porous media to describe the structural characteristics of the inner wall material of the sealed container, combines the diffusion barrier function of water molecules in the polymer material to establish a nonlinear mapping relationship between the water molecule adsorption isotherm and the local concentration gradient of ethylene carbonate, and dynamically updates the activation energy of the hydrolysis reaction based on the transition state theory, thus coupling the material interface characteristics with the chemical reaction.

4. The method for simulating the storage stability of ethylene carbonate based on digital twins according to claim 1, characterized in that: When the dynamic water adsorption and reaction kinetics module calculates the water molecule diffusion and permeation rate in real time, it uses the porous medium permeability tensor and the preset water molecule diffusion barrier function to solve the unsteady diffusion equation and obtain the instantaneous flux of water molecules in the container wall. Specifically, when calculating the hydrolysis reaction rate, the amount of reactant consumed per unit time is derived based on the local concentration gradient of ethylene carbonate and the updated activation energy of the hydrolysis reaction using a microscopic reaction probability model. When calculating the product concentration distribution, the product transport equation is constructed based on the reactant consumption and diffusion flux, and the three-dimensional spatial concentration field is solved by applying the principle of mass conservation.

5. The method for simulating the storage stability of ethylene carbonate based on digital twins according to claim 4, characterized in that: When designing the dynamic multiphysics real-time feedback module, the water molecule diffusion and permeation rate is used as the mass transfer field input variable, the hydrolysis reaction rate is used as the chemical reaction field input variable, and the product concentration distribution is used as the convection diffusion field input variable. The data streams from temperature and relative humidity sensors are integrated to drive the boundary conditions of the heat transfer field. The temporal synchronization of the four types of physical fields is coordinated through a field coupling controller, and an adaptive mesh re-division mechanism with the product accumulation rate as the feedback signal is established.

6. The method for simulating the storage stability of ethylene carbonate based on digital twins according to claim 5, characterized in that: When constructing the heat transfer, mass transfer, and chemical reaction kinetic equations, a water diffusion equation is established based on the principle of mass conservation, an unsteady-state heat conduction equation is established based on the principle of energy conservation, and a hydrolysis reaction rate equation is established based on the transition state theory. By introducing the inhibition coefficient of product concentration distribution on the reaction rate, the three are combined into a nonlinear partial differential equation system with cross terms. The operator splitting method is used to decompose the equation system into independently solvable heat transfer sub-equations, proton transfer equations, and reaction kinetic sub-equations.

7. The method for simulating the storage stability of ethylene carbonate based on digital twins according to claim 6, characterized in that: When adjusting the water molecule diffusion and hydrolysis reaction rate in real time, the water flux output by the proton transfer equation is used as the boundary condition update value of the reaction kinetic sub-equation, and the product inhibition coefficient generated by the reaction kinetic sub-equation is fed back to the proton transfer equation to correct the diffusion coefficient, forming a two-way real-time feedback link. The solution step size of each sub-equation is dynamically adjusted by the field coupling controller to match the sensor data sampling frequency.

8. The method for simulating the storage stability of ethylene carbonate based on digital twins according to claim 7, characterized in that: When simulating the dynamic equilibrium of moisture gradient and product accumulation effect in a closed container, the spatial interpolation algorithm is used to reconstruct the three-dimensional moisture distribution field inside the container based on the output data of the two-way real-time feedback link. Combined with the product transport equation, the product accumulation area at different time points is predicted. The extreme value of moisture gradient and product accumulation concentration are mapped to the stability decay risk index by the normalization method. The stability decay risk index includes the moisture penetration risk component and the product corrosion risk component.

9. The method for simulating the storage stability of ethylene carbonate based on digital twins according to claim 8, characterized in that: When predicting the evolution trend of key quality indicators, the moisture penetration risk component in the stability decay risk index is associated with the change rate of acid value and moisture content, and the product corrosion risk component is associated with the change rate of color. The time series extrapolation method is used to generate the indicator evolution curve. When identifying the stability failure threshold, the acid value exceeding the standard critical point, the moisture saturation concentration point and the color change threshold are located based on the abrupt change inflection point of the stability decay risk index. The acid value exceeding the standard critical point, the moisture saturation concentration point and the color change threshold are integrated into a composite failure criterion.

10. The method for simulating the storage stability of ethylene carbonate based on digital twins according to claim 8, characterized in that: When iteratively correcting the parameters of the digital twin basic model, the detection data sequence of acid value, moisture content and color of the actual stored sample is time-series aligned with the predicted values ​​of the digital twin basic model at the corresponding time nodes. The prediction deviation weight is calculated through the composite failure criterion, and the correlation between the deviation weight and the water molecule diffusion and permeation rate coefficient and the hydrolysis reaction rate constant is analyzed using a long short-term memory neural network. The parameter correction amount is generated and backpropagated to the digital twin basic model to achieve closed-loop optimization.

Images