Thermodynamic modeling and simulation method of multi-component glass raw material melting process

Abstract

The application provides a thermodynamic modeling and simulation method for a multi-component glass raw material melting process, comprising: obtaining multi-component glass raw material parameters, glass melting furnace structure parameters and melting operation parameters, establishing a global calculation domain and a global low-order thermodynamic model; solving the global low-order thermodynamic model by time advancing to obtain a global low-order state field, calculating mass conservation residuals, energy conservation residuals, component conservation residuals and phase consistency residuals to form residual-state correlation data; identifying a continuous abnormal cluster according to the residual-state correlation data, determining a local upgrading area, a dominant residual type and an upgrading level, and establishing a local high-order thermodynamic model; solving the local high-order thermodynamic model, conservatively mapping local high-order solving results to corresponding basic space units, updating coarse-fine interface flux and equivalent model parameters, and obtaining a convergent global thermodynamic model and a local upgrading evolution record.

Description

Technical Field

[0001] This invention relates to the field of glass melting process simulation technology, and in particular to a thermodynamic modeling and simulation method for the melting process of multi-component glass raw materials. Background Technology

[0002] During the melting process of multi-component glass raw materials, several processes occur, including moisture evaporation, component reactions, phase transformation, formation of the primary liquid phase, dissolution of residual particles, and gas release. These processes are influenced by factors such as raw material composition, particle size distribution, moisture content, feeding rate, heating conditions, and glass melting furnace structure, exhibiting characteristics of multi-component, multi-phase, and strong coupling. Obtaining temperature, chemical composition, liquid phase ratio, melting front, and unmelted material distribution through thermodynamic modeling and simulation is crucial for analyzing the melting state and adjusting operating parameters.

[0003] Existing glass melting simulations typically employ a uniform spatial grid, time step, and reaction model across the entire domain. Using a coarse discretization scale makes it difficult to accurately describe the drastic temperature changes, component transformations, and phase transitions at the feeding port, melting front, and local heating areas. Conversely, using a globally refined model significantly increases computational load and simulation time. While some methods set up locally refined regions, their extent is often pre-defined and difficult to dynamically adjust as abnormal regions arise, spread, and disappear, easily leading to insufficient local accuracy or over-computation in unnecessary areas.

[0004] Furthermore, if the results transfer between the local fine model and the global model is achieved by directly averaging state variables, it may disrupt the conservation relationships of material mass, total enthalpy, chemical composition, and phase mass. Existing methods often use temperature or liquid phase ratio as the basis for determining the completion of melting, but they do not adequately consider residual particles, the degree of reaction completion, and the state of gas retention, making it difficult to accurately determine the complete melting time and evaluate operating strategies.

[0005] Therefore, this invention proposes a thermodynamic modeling and simulation method for the melting process of multi-component glass raw materials. The information disclosed in the background section is only for enhancing understanding of the background of this disclosure and may therefore contain prior art information that is not common knowledge to those skilled in the art. Summary of the Invention

[0006] The purpose of this invention is to address the shortcomings of existing technologies by providing a thermodynamic modeling and simulation method for the melting process of multi-component glass raw materials, thereby solving the technical problems mentioned in the background section.

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

[0008] Thermodynamic modeling and simulation methods for the melting process of multi-component glass raw materials include the following steps:

[0009] S1. Obtain parameters of multi-component glass raw materials, glass melting furnace structure parameters and melting operation parameters, and establish a global computational domain and a global low-order thermodynamic model.

[0010] S2. Perform time-progression solution on the global low-order thermodynamic model to obtain the global low-order state field, calculate the mass conservation residual, energy conservation residual, composition conservation residual and phase consistency residual of each basic space unit, and form residual and state correlation data.

[0011] S3. Identify persistent anomaly clusters based on residual and state correlation data, determine local upgrade regions, dominant residual types and upgrade levels by combining buffer space units, and establish local high-order thermodynamic models based on local states, boundary states and mapping relationships between coarse and fine space units.

[0012] S4. Solve the local higher-order thermodynamic model, conserve and map the local higher-order solution results to the corresponding basic space unit, update the flux at the coarse and fine interface and the equivalent model parameters, obtain the global dynamic model, and adjust the local order-up region according to the joint residual field until a converged global thermodynamic model and local order-up evolution record are obtained.

[0013] S5. Based on the convergent global thermodynamic model and local ascending evolution record, extract the melting front, spatial distribution of unmelted material and complete melting time, determine the sensitive kiln area and sensitive factor set, generate and simulate and evaluate candidate operating strategies, determine the target operating strategy from the candidate operating strategies that meet the process constraints and model convergence conditions, and output the operating setpoint and melting simulation results.

[0014] S1 specifically includes: acquiring multi-component glass raw material parameters, glass melting furnace structural parameters, and melting operation parameters; performing time alignment, unit conversion, and missing data processing on feeding, composition, thermophysical properties, reaction, and operation monitoring data to form a basic parameter set for melting simulation; establishing a three-dimensional computational domain based on the glass melting furnace structural parameters, dividing basic spatial units and determining surface adjacency relationships, setting initial states, dynamic function labels, and feeding, heating, discharging, and heat dissipation boundaries to form a global computational domain; and establishing and coupling sub-models for raw material heating, component reaction, phase transformation, gas generation, and material migration based on the basic parameter set for melting simulation and the global computational domain to form a global low-order thermodynamic model.

[0015] S2 specifically includes: reading the global low-order thermodynamic model, the melting simulation basic parameter set, and the global computational domain; performing time-progressive solutions for raw material heating, component reaction, phase transformation, gas generation, and material migration to form a global low-order state field; based on the global low-order state field, calculating the mass conservation residual, energy conservation residual, component conservation residual, and phase consistency residual of each basic spatial unit to form multiple types of residual subfields; standardizing the multiple types of residual subfields, calculating the joint residual intensity in combination with residual weights, determining the number of residual co-occurrences, spatial continuity, residual duration, and residual propagation direction, and associating them with the global low-order state field and dynamic function labels to form residual and state-related data.

[0016] S3 specifically includes: identifying persistent anomaly clusters based on residual and state correlation data, setting buffer space units around them to form a set of local upgrade regions; calculating the residual contribution ratio in each local upgrade region, determining the dominant residual type and upgrade level, selecting the process to be refined and the reaction stage to be split in combination with dynamic function tags, extracting the initial state, boundary state and coarse-fine spatial unit mapping relationship to form local upgrade configuration data; refining the spatial units and time steps based on the local upgrade configuration data, conserving and allocating the extensive state variables, splitting the actual reaction stages, establishing the relationships of heat transfer, component reaction, phase transformation, gas generation and material migration to form a local high-order thermodynamic model.

[0017] S4 specifically includes: solving the local higher-order thermodynamic model, executing multiple refined time steps within a preset global time step to obtain the local higher-order state field and the local higher-order solution results including the endpoint state, cumulative income and expenditure, cumulative flux at the coarse-fine interface, and local non-convergence markers; performing a conservation mapping on the local higher-order solution results based on the mapping relationship between coarse and fine spatial units, summarizing the extensive state variables and inversely calculating the temperature, updating the flux at the coarse-fine interface and the equivalent model parameters to form a global dynamic model; recalculating the joint residual field of the global dynamic model, adjusting the local ascending region and ascending level according to the residual changes, and repeating the local higher-order solution and conservation mapping until a converged global thermodynamic model and local ascending evolution record are obtained.

[0018] S5 specifically includes: reading the convergent global thermodynamic model and local ascending evolution records, extracting the melting front, spatial distribution of unmelted material, and complete melting time, and determining the sensitive kiln area and the set of sensitive factors; generating candidate operating strategies that meet process constraints based on the set of sensitive factors, re-executing each candidate operating strategy as a dynamic modeling driving force, recording melting indices and local ascending computational costs, and forming a candidate strategy simulation evaluation set; screening the candidate strategy simulation evaluation set based on process hard constraints, comprehensively evaluating the complete melting time and local ascending computational costs, verifying and determining the target operating strategy, and outputting the operating setpoints, melting simulation results, and constraint verification results.

[0019] The beneficial effects of this invention are as follows:

[0020] This invention establishes a global low-order thermodynamic model encompassing raw material heating, component reactions, phase transitions, gas generation, and material migration. This model continuously characterizes the coupled melting process of multi-component glass raw materials, improving the completeness of melting state simulation. By jointly analyzing mass conservation residuals, energy conservation residuals, component conservation residuals, and phase consistency residuals, it accurately identifies local anomaly regions, avoiding missed detections caused by relying solely on temperature or empirical location to determine sensitive areas.

[0021] This invention dynamically adjusts the local upgrade region and upgrade level based on the residual duration, spatial continuity, and propagation direction, allowing the fine calculation to change as the anomalous region arises, migrates, and disappears, thus balancing local computational accuracy with global computational efficiency. By refining the spatial units, time steps, and reaction stages within the local upgrade region, the computational resolution of the melting front, the concentrated reaction region of components, and the residual particle region can be improved, while reducing the computational resource waste caused by uniform refinement across the entire domain.

[0022] This invention ensures the consistency of total mass, total enthalpy, chemical composition mass, and mass of each phase by performing a conservation mapping on the local high-order solution results and updating the flux at the coarse-fine interface and equivalent model parameters, thereby reducing numerical errors caused by model integration. By comprehensively evaluating the melting state based on reaction progress, unmelted matter, residual gas mass fraction, and complete melting time, and by constraining and evaluating candidate operating strategies through simulation, this invention provides a basis for optimizing the feeding rate, heating load, and feeding interval. Attached Figure Description

[0023] Figure 1 This is a flowchart of the thermodynamic modeling and simulation method for the melting process of multi-component glass raw materials according to the present invention. Detailed Implementation

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

[0025] Example: Figure 1 As shown, this embodiment provides a thermodynamic modeling and simulation method for the melting process of multi-component glass raw materials, including the following steps:

[0026] S1. Obtain parameters of multi-component glass raw materials, glass melting furnace structure parameters and melting operation parameters, and establish a global computational domain and a global low-order thermodynamic model.

[0027] S2. Perform time-progression solution on the global low-order thermodynamic model to obtain the global low-order state field, calculate the mass conservation residual, energy conservation residual, composition conservation residual and phase consistency residual of each basic space unit, and form residual-state correlation data.

[0028] S3. Identify persistent anomaly clusters based on residual-state correlation data, determine local upgrade regions, dominant residual types and upgrade levels by combining buffer space units, and establish local high-order thermodynamic models based on local states, boundary states and mapping relationships between coarse and fine space units.

[0029] S4. Solve the local higher-order thermodynamic model, conserve and map the local higher-order solution results to the corresponding basic space unit, update the flux at the coarse and fine interface and the equivalent model parameters, obtain the global dynamic model, and adjust the local order-up region according to the joint residual field until a converged global thermodynamic model and local order-up evolution record are obtained.

[0030] S5. Based on the convergent global thermodynamic model and local ascending evolution record, extract the melting front, spatial distribution of unmelted material and complete melting time, determine the sensitive kiln area and sensitive factor set, generate and simulate and evaluate candidate operating strategies, determine the target operating strategy from the candidate operating strategies that meet the process constraints and model convergence conditions, and output the operating setpoint and melting simulation results.

[0031] S1 specifically includes the following sub-steps:

[0032] S110. Obtain parameters of multi-component glass raw materials, glass melting furnace structure parameters, and melting operation parameters. Multi-component glass raw material parameters include the wet-basis feeding rate, moisture content, particle size distribution, chemical component content, density, specific heat capacity, thermal conductivity, emissivity, viscosity, chemical component diffusion coefficient, reaction temperature range, reaction enthalpy, gas production per unit mass, and residual particle solubility coefficient for each raw material type.

[0033] The wet-basis feeding rate was obtained by a weighing feeder or belt scale; the moisture content was obtained by an online moisture meter or drying and weighing test; the particle size distribution was obtained by sieving test or laser particle size analysis; and the chemical component content was obtained by raw material batch inspection report or X-ray fluorescence spectroscopy. Density, specific heat capacity, thermal conductivity, emissivity, viscosity, and chemical component diffusion coefficient were obtained by density test, differential scanning calorimetry (DSC), thermal conductivity test, radiation calibration test, viscosity test, and isothermal diffusion test, respectively. The reaction temperature range, reaction enthalpy, and gas production per unit mass were obtained by differential scanning calorimetry-thermogravimetric-mass spectrometry (DSC-MS / MS). The residual particle solubility coefficient was obtained by particle dissolution tests under different temperature and particle size conditions. The chemical component content is the mass fraction of the corresponding chemical component relative to the dry basis mass of the raw material.

[0034] Convert the wet basis feeding rate of the i-th raw material to the dry basis feeding rate:

[0035]

[0036] In the formula, i represents the raw material type number; t represents the collection time; and These represent the dry basis feeding rate and the wet basis feeding rate, respectively, both in kg / h. The value represents the moisture content and is dimensionless. The difference between the wet basis feeding rate and the dry basis feeding rate is taken as the moisture mass flow rate; before entering the conservation calculation with seconds as the time unit, kg / h is divided by 3600 to convert it to kg / s.

[0037] The structural parameters of the glass melting furnace include the internal cavity dimensions and the location and dimensions of the feeding port, combustion zone, electric melting zone, and discharge port, which are derived from the as-built drawings and on-site re-measurement results; the heating area includes the combustion zone and the electric melting zone.

[0038] The melting operation parameters include temperature, fuel flow rate, lower heating value of fuel, electric melting power, kiln pressure, glass liquid level, discharge rate, and conveying speed. Temperature is obtained from thermocouples or infrared thermometers; fuel flow rate is obtained from fuel flow meters; lower heating value of fuel is obtained from fuel inspection reports; electric melting power is obtained from voltage, current, and power factor measurement modules; and kiln pressure, glass liquid level, discharge rate, and conveying speed are obtained from pressure sensors, liquid level detection devices, discharge metering devices, and conveying equipment encoders, respectively. Raw temperature values ​​are recorded in °C and converted to K by adding 273.15 before inputting into thermodynamic calculations.

[0039] All data are aligned according to the clock of the glass melting furnace control system. For data missing for no more than two consecutive acquisition cycles, linear interpolation is used. If data is missing for more than two consecutive acquisition cycles, the most recent valid value is retained to the current time step, and a data anomaly flag is set. Dynamic updates of the corresponding parameters cease when the preset failure time is reached. For example, when the wet basis feeding rate of quartz sand is 1000 kg / h and the moisture content is 0.03%, the dry basis feeding rate is 970 kg / h, with an additional 30 kg / h used as the moisture mass flow rate.

[0040] The standardized data is used to form a set of basic parameters for melting simulation and then input into S120.

[0041] S120. Using the reference point on the inner surface of the kiln bottom at the feeding end as the origin of the coordinate system, establish a 3D coordinate axis along the length, width, and height of the glass melting furnace. Establish a 3D calculation domain that includes the raw material layer, molten glass, and the upper gas phase space that participates in radiative heat exchange. When the refractory material is not included in the 3D calculation domain, the inner surface of the refractory material is used as the heat boundary.

[0042] Basic spatial units are generated based on the internal dimensions of the glass melting furnace. Local densification is applied to the areas around the feeding port, heating zone, discharge port, and concentrated melting front regions determined by historical stable operating data. If historical stable operating data is unavailable, an initialization simulation is performed first, and regions with concentrated temperature and liquid phase ratio gradients are designated as locally densified areas. When two basic spatial units share a common surface, they are identified as face-adjacent spatial units, and the area of ​​the common surface, center distance, and normal direction of the common surface are recorded. Basic spatial units that only contact each other through edges or vertices are not considered face-adjacent.

[0043] The feeding area and the discharging area are set as fixed functional areas, while the raw material heating area, component reaction area, primary liquid phase area and glass melt migration area are set as dynamic functional labels. The dynamic functional labels are updated step by step according to the reaction progress and phase state in the basic space unit. Only one main dynamic functional label is set in the same basic space unit at the same time step.

[0044] Each basic spatial unit is configured with an initial temperature, initial chemical composition content, initial reaction progress for each reaction stage, initial liquid phase ratio, initial retained gas mass fraction, and initial phase state. The liquid phase ratio is the ratio of the sum of the mass of the primary liquid phase and the mass of the molten glass to the total mass of the material; the retained gas mass fraction is the ratio of the mass of the gas that has not yet escaped to the total mass of the material; the phase state includes solid raw materials, reaction transition products, primary liquid phase, and molten glass.

[0045] The initial temperature is obtained by spatial interpolation of the measured temperature from the temperature measuring device; the initial chemical composition content is obtained by mapping the residence time distribution determined by the feeding record, conveying speed, and tracer test or historical stable operation data. When starting up for the first time and there is no historical stable operation data, the initial simulation results corresponding to the design conditions are used, and the initial simulation is continued until the state changes of each basic spatial unit meet the stability conditions.

[0046] The feeding boundary is set based on the dry basis feeding rate and moisture mass flow rate of S110; the heating boundary is set based on the fuel flow rate, lower heating value of the fuel, effective combustion feed coefficient, electric fluxing power, and effective electrothermal utilization coefficient; the discharge boundary is set based on the discharge rate; and the heat dissipation boundary is set based on the equivalent thermal resistance determined by the refractory material inner surface temperature, ambient temperature, refractory material thickness, and thermal conductivity. The effective combustion feed coefficient and effective electrothermal utilization coefficient are obtained from the heat balance test of the glass melting furnace.

[0047] For example, when the feed inlet width is 0.6m, basic spatial units with side lengths of 0.1 to 0.2m are used within its coverage area, while basic spatial units with side lengths of 0.3 to 0.5m are used in other areas to avoid merging and assigning values ​​to different feed boundaries. A global computational domain is formed and input into S130.

[0048] S130. Establish a raw material heating sub-model, a component reaction sub-model, a phase transformation sub-model, a gas generation sub-model, and a material migration sub-model within the global computational domain, and couple them to form a global low-order thermodynamic model.

[0049] A global low-order thermodynamic model refers to a thermodynamic numerical model that uses basic spatial units and a global preset time step to reduce continuous reactions to finite reaction stages. It is used to calculate temperature, chemical component content, reaction progress of each reaction stage, mass of each phase, liquid phase ratio, mass fraction of stagnant gas, mass flux, heat flux, and chemical component flux. The term "low-order" is used to characterize its spatial discrete scale, temporal resolution, and reaction stage resolution, which are lower than those of local high-order thermodynamic models.

[0050] The raw material heating sub-model updates the temperature based on the equivalent heat capacity of the basic spatial unit, the thermal power between adjacent spatial units, external heating power, heat dissipation power, reaction heat consumption power, and latent heat of phase transformation. The component reaction sub-model, based on the temperature-conversion rate curve obtained from S110, uses piecewise linear interpolation to determine the reaction progress increment for each merged reaction stage, and calculates the consumption, production, heat consumption, and gas production mass of chemical components according to stoichiometry. For example, if the conversion rates for a certain merged reaction stage are 0.55 and 0.95 at 750℃ and 850℃ respectively, then the conversion rate at 800℃ is 0.75.

[0051] The phase transformation sub-model determines the transformation variables between solid feedstock, reaction transition material, primary liquid phase, and molten glass based on temperature, chemical composition content, and reaction progress. The gas generation sub-model calculates the gas generation mass based on the reaction progress increment and gas production per unit mass, and calculates the gas escape mass and retained gas mass fraction based on the gas escape boundary. The material migration sub-model calculates the mass flux and chemical composition flux based on the mass, density, viscosity, chemical composition diffusion coefficient, surface adjacency relationship, and boundary conditions of each phase.

[0052] Each time step sequentially performs calculations for raw material heating, component reaction, phase transformation, gas generation, and material migration, and uses the output of this time step as the initial state for the next time step. Changes in raw material input, heating boundary conditions, and the state of the basic spatial units are defined as dynamic modeling drivers. These drivers are used to update the model state, equivalent model parameters, and boundary conditions without re-dividing the global computational domain. A global low-order thermodynamic model is formed and input into S210.

[0053] S2 specifically includes the following sub-steps:

[0054] S210: Read the global low-order thermodynamic model, melting simulation basic parameter set, and global computational domain. The preset time step is the minimum value among the unified acquisition cycle of the equipment, the time required for the material to cross the smallest basic spatial unit, and the shortest characteristic time of the merging reaction stage; the unified acquisition cycle of the equipment comes from the data time alignment result of S110, the material crossing time is obtained by dividing the distance between the centers of the basic spatial units by the material migration speed, and the shortest characteristic time is obtained by the time taken for the reaction conversion rate to change from 0.10 to 0.90 in the synchronous thermal analysis.

[0055] At the beginning of each time step, the dynamic modeling driving quantities are read, and the calculations of raw material heating, component reaction, phase transformation, gas generation and material migration are executed in sequence to obtain the temperature, chemical component content, reaction progress of each reaction stage, mass of each phase, liquid phase ratio, mass fraction of stagnant gas, mass flux, heat flux and chemical component flux of each basic space unit.

[0056] After completing one coupled calculation, the temperature, liquid phase ratio, and stagnant gas mass fraction are compared between two adjacent internal iterations. The current time step ends when the changes in all three do not exceed the corresponding iteration tolerance. Otherwise, the intermediate state is updated using under-relaxation, and iteration continues. This involves weighting the state variables calculated in the current internal iteration with those from the previous iteration using a preset under-relaxation factor, and using this weighted average as the input value for the next iteration. The iteration tolerance is determined by a time-step independence test, and the maximum number of iterations is determined by adding a safety margin to the upper limit of the number of iterations required for stable convergence in the time-step independence test. If convergence is not achieved after reaching the maximum number of iterations, a global non-convergence flag is set.

[0057] The state quantities of each basic spatial unit at the beginning and end of the time step, as well as the external feed mass, discharge mass, gas escape mass, migration mass, cumulative heat transfer, external heat supply, reaction heat consumption, phase transformation heat consumption, chemical component conversion amount, and gas generation mass within the time step, are stored according to the basic spatial unit number and time step number to form a global low-order state field and input into S220. The global low-order state field is a spatiotemporal data set used to recalculate the material, energy, chemical composition, and phase budget.

[0058] Mass flux is the mass passing through a common surface per unit area per unit time, and its unit is kg / (m²). 2 Mass flow rate (·s) is the product of mass flux and common surface area, expressed in kg / s; transported mass is the product of mass flow rate and time step, expressed in kg. Heat flux, heat power, and cumulative heat transfer are expressed in W / m². 2W and J. For example, if the initial mass of a basic space unit is 500 kg, 15 kg flows in and 14 kg flows out within 5 seconds without any external material input, then the theoretical final mass is 501 kg. Therefore, both state and process quantities should be preserved.

[0059] S220. For each basic spatial unit and each time step, calculate the mass conservation residual, energy conservation residual, composition conservation residual, and phase consistency residual respectively.

[0060] The residual of mass conservation is determined according to the following formula:

[0061]

[0062] In the formula, This represents the mass conservation residual of the basic spatial unit c at the nth time step, in kg. and These represent the total mass of material at the beginning and end of the time step, respectively, in kg; This indicates the preset time step, in seconds. This represents the set of face-adjacent spatial units of the basic spatial unit c; This represents the net mass flow rate from adjacent spatial unit p to the base spatial unit c, expressed in kg / s, with inflow being positive and outflow being negative. This represents the external net mass flow rate resulting from material feeding, discharging, and gas escape, expressed in kg / s. Before gas escape, the total mass of the material is included; after gas escape, only the external net mass flow rate is considered.

[0063] The energy conservation residual is determined according to the following formula:

[0064]

[0065] In the formula, This represents the energy conservation residual of the basic spatial unit c at the nth time step, in J. and These represent the total enthalpy at the beginning and end of the time step, respectively, in J; and These represent the cumulative energy of the input and output basic spatial unit c within this time step, respectively, in J, including thermal conduction, radiation, external heating, and material-carried enthalpy; This represents the cumulative heat consumed in a chemical reaction, expressed in J. This represents the cumulative heat consumption during phase transition, expressed in J. Energy already included in the heat consumption of the reaction or phase transition is not included again in the total enthalpy increment.

[0066] For each chemical component j defined in S110, calculate its independent component mass conservation residual:

[0067]

[0068] In the formula, Represents the residual of chemical component j; and These represent the remaining quantities of component j at the end and beginning of the time step, respectively. The net mass flow rate of component j, which is input to the adjacent spatial cell; and These represent the reaction formation rate and the consumption rate, respectively. After calculating the residuals of all chemical components, the maximum value among the standardized absolute residuals of each component is taken as the final component conservation residual for that basic spatial unit.

[0069] The phase consistency residual is the larger of the absolute difference between the sum of the mass ratios of the four phases and 1, and the difference between the current liquid phase ratio and the liquid phase ratio recalculated based on temperature, chemical composition content, and reaction progress. For example, when the mass ratios of the four phases are 0.40, 0.30, 0.25, and 0.15, the sum of the ratios is 1.10, and the phase consistency residual is at least 0.10.

[0070] Record the values, locations, time step numbers, durations, and non-convergence flags of various residuals to form multiple residual subfields, and input them into S230.

[0071] S230. Based on the grid independence test, time step independence test, and stable operation samples of the glass melting furnace, four types of residual tolerance values ​​are determined respectively. The residual tolerance value shall not be lower than the upper limit of the corresponding numerical dispersion error, and the 95th quantile value of the standardized residual corresponding to the stable operation sample shall be used as the initial tolerance value.

[0072] When the system is deployed for the first time and there is a lack of historical stable operating samples, the initial residual tolerance values ​​for mass conservation, energy conservation, composition conservation, and phase consistency are preset to 10 by default. -4 10 -3 10 -4 and 10 -4 Zero reference quantity The value is uniformly set to a fixed minimum value (e.g., 10). -6 This ensures that the system can perform cold start operations when there is no historical data.

[0073] The residuals of different dimensions are standardized, and the joint residual strength is calculated:

[0074]

[0075] In the formula, k represents the residual type number; This represents the standardized residual, which is dimensionless. The basis for the income and expenditure of the k-th type of residual is represented by the larger value between the corresponding inventory and the input and output quantities for the mass, energy and component residuals, respectively, and the basis for the income and expenditure of the phase consistency residual is taken as 1. This represents the zero-prevention reference quantity corresponding to the k-th type of residual, and its unit is . The same applies, and the corresponding measurement resolution or lower limit of numerical calculation accuracy is taken; Let represent the weight of the k-th residual, which is dimensionless and the sum of the weights of the four residuals is 1. This represents the joint residual strength, which is dimensionless.

[0076] Each residual weight is calibrated based on stable operating samples and known abnormal operating condition samples: the recognition accuracy of the four types of standardized residuals for known abnormal spatial units is calculated, and the corresponding residual weight is determined according to the proportion of each recognition accuracy to the total recognition accuracy of the four types; when there are no known abnormal operating condition samples, the weights of the four types of residuals are all 0.25, and remain unchanged in this simulation.

[0077] The number of residual types exceeding the corresponding allowable value within the same basic spatial unit and the same time step is defined as the residual co-occurrence quantity; the proportion of the number of face-adjacent spatial units with the same type of out-of-limit residuals to the total number of face-adjacent spatial units is defined as the spatial continuity; the product of the number of time steps in which the joint residual intensity continuously exceeds the limit and the preset time step length is defined as the residual duration; the displacement direction of the weighted center of the residual intensity of adjacent time steps is defined as the residual propagation direction, and when the displacement is less than the distance between the centers of one basic spatial unit, it is recorded as no clear propagation direction.

[0078] The joint residual state is composed of four types of standardized residuals, joint residual strength, number of residual co-occurrences, spatial continuity, residual duration, and residual propagation direction, and is arranged in space and time to form a joint residual field.

[0079] Using the basic spatial unit number and time step number as the association key, the joint residual state is mapped one by one to the state variables, process variables and dynamic function labels in the global low-order state field to form residual-state association data, which is then input into S310.

[0080] S3 specifically includes the following sub-steps:

[0081] S310: Read the residual-state correlation data output by S230, and use the surface adjacency relationship formed by S120, the residual tolerance value determined by S230, and the non-convergence flag. Define the basic spatial unit that has at least one type of standardized residual exceeding the corresponding tolerance value in the same time step, or that has a global non-convergence flag, as the initial abnormal spatial unit; define the basic spatial unit that continuously exceeds the limit in consecutive time steps, or that has a spatial continuity degree reaching a preset continuity ratio, as the effective abnormal spatial unit. The preset continuity ratio and preset duration are respectively adopted from the 95th percentile value of the corresponding index in the stable running samples.

[0082] Within each time step, effective anomaly spatial units that are face-adjacent and have the same type of out-of-limit residuals are aggregated into instantaneous anomaly clusters. Then, instantaneous anomaly clusters in adjacent time steps are tracked. When two instantaneous anomaly clusters have spatial overlap, face adjacency, or the subsequent instantaneous anomaly cluster is located in the residual propagation direction of the previous instantaneous anomaly cluster, they are merged into candidate anomaly clusters. When both instantaneous anomaly clusters have a definite residual propagation direction, the angle between the residual propagation direction and the direction of the line connecting the cluster centers should not exceed 45°; basic spatial units that only contact each other through edges or vertices do not establish connectivity.

[0083] Candidate anomalous clusters are identified as transient anomalous clusters and excluded if they persist for only one time step, all residuals recover to the acceptable range in the next time step, they do not contain non-converged markers, and no new face-adjacent anomalous spatial elements are added; other candidate anomalous clusters are identified as persistent anomalous clusters. One face-adjacent normal spatial element is added around the persistent anomalous cluster as a buffer spatial element. The persistent anomalous cluster and the buffer spatial element together constitute a local ascending-order region. The local ascending-order region stops expanding outward from the physical boundary when it reaches the feed port, discharge port, or the inner surface boundary of the refractory material. The buffer spatial element is used to provide local higher-order solution boundaries and does not participate in the dominant residual statistics.

[0084] For example, if four adjacent basic spatial units simultaneously exceed the limits of energy conservation residuals and phase consistency residuals for three consecutive time steps, and the residual weighting center moves along the glass melt migration direction, it is identified as a continuous anomalous cluster; if a single basic spatial unit recovers only after exceeding the limit in one time step, it does not trigger local order increase.

[0085] The list of abnormal spatial cells, the list of buffer spatial cells, the residual duration, the residual propagation direction, and the dynamic function label are combined to form a local ascending region set, and then input into S320.

[0086] S320. For each local up-order region, establish a set of basic spatial unit-time step data for that region during the tracking period. And calculate the contribution ratio of the k-th type residual in the local up-order region g:

[0087]

[0088] In the formula, g represents the local upgrade region index; k represents the residual type index to be calculated; and u represents the residual type summation index. This represents the contribution ratio of the k-th type of residual, and is dimensionless. The definition in S230 is adopted.

[0089] The residual type with the largest contribution ratio is determined as the dominant residual type; when the difference between the largest contribution ratio and the second largest contribution ratio is less than the preset contribution difference value, the corresponding two types of residuals are determined as parallel dominant residual types. The preset contribution difference value is the 95th percentile of the contribution ratio difference value in the stable running samples.

[0090] When mass conservation residuals are dominant, the mass migration of solid raw materials, reaction transition products, primary liquid phase, molten glass, and stagnant gases is refined; when energy conservation residuals are dominant, heat conduction, radiative heat transfer, material enthalpy, reaction heat consumption, and latent heat of phase transformation are refined; when composition conservation residuals are dominant, merged chemical components are broken down into specific reactants, intermediate reactants, and products; when phase consistency residuals are dominant, the phase transformation processes are refined. The upgraded content corresponding to parallel dominant residuals is activated simultaneously.

[0091] Based on the dynamic function labels, the content to be refined is as follows: in the raw material heating area, the water evaporation and heat transfer should be refined first; in the component reaction area, the specific reaction stages should be refined first; in the primary liquid phase area, the phase transformation should be refined first; and in the glass melt migration area, mass migration and chemical component diffusion should be refined first.

[0092] The escalation level is determined based on the joint residual strength, residual duration, and spatial continuity: the first escalation level is used when only one item exceeds the corresponding escalation threshold, the second escalation level is used when any two items exceed the corresponding escalation threshold, the third escalation level is used when all three items exceed the corresponding escalation threshold, or when any one item reaches twice the corresponding escalation threshold; each escalation threshold is the 95th percentile of the corresponding indicator in the stable operating sample.

[0093] The temperature, chemical composition content, reaction progress, mass of each phase, liquid phase ratio, and mass fraction of retained gas in the local ascending region are extracted, as well as the heat flux, mass flux, and chemical composition flux on the common surface between the buffer space unit and the surrounding basic space unit.

[0094] The local upgrade configuration data is composed of the region number, dominant residual type, upgrade level, spatial refinement factor, temporal refinement factor, reaction stage to be split, local initial state, local boundary state, and mapping relationship between coarse and fine spatial units, and then input into S330.

[0095] S330, a local higher-order thermodynamic model, refers to a thermodynamic numerical model that improves spatial resolution, temporal resolution, or reaction stage resolution within a local higher-order region, and is used to calculate the state of refined spatial units and cross-surface fluxes.

[0096] Spatial and temporal refinement is performed based on local upgrade configuration data:

[0097]

[0098] In the formula, and These represent the feature sizes of the refined spatial unit and the corresponding basic spatial unit within the local upgraded region g, respectively, in meters; Indicates the space refinement factor; This indicates the refined time step, in seconds (s). The definition in S210 is adopted; This represents the temporal refinement factor; both the spatial refinement factor and the temporal refinement factor are integers not less than 2, and their specific values ​​are determined by the mesh independence test and the local time step independence test.

[0099] The total mass of materials, total enthalpy, mass of chemical components, mass of each phase, and mass of stagnant gas in the basic spatial unit are initially allocated according to the volume of the refined spatial unit; when there are temperature gradients, chemical component gradients, or phase interfaces, spatial interpolation is performed based on the state of the adjacent basic spatial units.

[0100] Let the sum of a certain extensive state variable after spatial interpolation in all refined spatial cells be... The extended state variable corresponding to the basic spatial unit is Then the interpolation results of each refined spatial unit are multiplied uniformly by . Make the corrected sum equal to The extended state quantities include the total mass of the material, total enthalpy, mass of each chemical component, mass of each phase, and mass of the retained gas.

[0101] Based on the types of raw materials, chemical component contents, and simultaneous thermal analysis results in S110, the actual stages of water evaporation, carbonate decomposition, sulfate reaction, primary liquid phase formation, residual particle dissolution, and gas release in the current formulation are determined. If the corresponding raw material is absent or the corresponding reaction peak is not detected, the corresponding reaction branch is not established. For example, the global low-order thermodynamic model groupes the decomposition of limestone and dolomite into one reaction stage, while the thermogravimetric curve shows two independent weight loss intervals. Therefore, separate limestone decomposition and dolomite decomposition stages are established, and the reaction progress, enthalpy, and carbon dioxide generation are configured for each stage.

[0102] Establish relationships for heat transfer, chemical component reactions, phase transitions, gas generation, and material migration between refined spatial units, ensuring that the sum of the cumulative transfer amounts at each refined sub-interface of the coarse-fine interface equals the cumulative transfer amount at the corresponding basic common surface. Form a local high-order thermodynamic model and input it into S410.

[0103] S4 specifically includes the following sub-steps:

[0104] S410 reads the local high-order thermodynamic model output by S330, along with its refined spatial units, refined time steps, local initial states, local boundary conditions, specific reaction stages, and mapping relationships between coarse and fine spatial units.

[0105] For a local upgrade region g, within a global preset time step... Continuous execution within A refined time step, satisfying:

[0106]

[0107] In the formula, This represents the time refinement factor of the local up-order region g; This indicates the refined time step, in seconds (s). This indicates the preset time step for the entire domain, in seconds.

[0108] The boundary temperature, chemical component content, phase ratio, and stagnation gas mass fraction of each refined time step are obtained by piecewise linear interpolation of the boundary states of adjacent global time nodes; the heat flux, mass flux, and chemical component flux at the coarse-fine interface are applied according to the distribution relationship established in S330.

[0109] Each refined time step sequentially calculates heat conduction and radiation heat transfer, chemical component consumption and generation, reaction heat consumption, phase transformation, gas generation and escape, and material migration, and then updates the temperature, chemical component content, reaction progress of each specific reaction stage, mass of each phase, liquid phase ratio, and mass fraction of retained gas.

[0110] In two consecutive internal iterations, if the changes in temperature, liquid phase ratio, and stagnant gas mass fraction of all refined spatial units do not exceed the corresponding local iteration tolerance value, the current refinement time step ends. The local iteration tolerance value is determined by the local time step independence test and is not greater than the corresponding iteration tolerance value in S210. If convergence is not achieved after reaching the maximum number of local iterations, a local non-convergence flag is set.

[0111] The state variables, process variables, and cross-plane fluxes arranged according to the refined spatial unit number and refined time step number are defined as local higher-order state fields. From the local higher-order state fields, the state variables of the global preset time step endpoint, the cumulative reaction amount, cumulative heat balance, cumulative material balance, cumulative flux of the coarse and fine interface, and local non-convergence markers are extracted to form local higher-order solution results and input into S420.

[0112] For example, if the global preset time step is 5s and the time refinement factor is 5, five 1s refinement time steps should be completed before the local result is matched with the global 5s endpoint state. The 1s state should not be used to replace the 5s state.

[0113] S420, conservation mapping, refers to converting local high-order solution results into basic spatial unit states based on the mapping relationship between coarse and fine spatial units, while keeping the total mass of materials, total enthalpy, mass of each chemical component, mass of each phase, and mass of stagnant gas consistent before and after the mapping.

[0114] For a basic spatial unit c, the set of refined spatial units it contains is denoted as c. Summarize the total mass and total enthalpy of the materials according to the following formula:

[0115]

[0116] In the formula, q represents the refined spatial unit number; and Let represent the total mass and total enthalpy of the material in the basic spatial unit c, obtained by mapping from the local higher-order solution results, in kg and J, respectively; and Let H and L represent the total mass and total enthalpy of the material in the refined spatial unit q, respectively; the superscripts H and L represent the higher-order model and the lower-order model, respectively. The masses of each chemical component, each phase, and the mass of the retained gas are summed using the same method.

[0117] The mapped temperature is not calculated by arithmetic averaging of the temperatures of the refined spatial units, but rather by inversely calculating the total enthalpy and material composition after mapping. The total enthalpy of the mapped basic spatial unit satisfies:

[0118]

[0119] In the formula, v represents the category number of the phase or chemical component; Indicates the number of phases or chemical components involved in the total enthalpy calculation; This represents the mass of the v-th type of matter within the basic space unit c, in kg. This indicates the reference temperature, expressed in Kelvin (K). This represents the temperature to be determined, in Kelvin (K). This represents the specific heat capacity of the vth substance as a function of temperature, expressed in J / (kg·K). This represents the enthalpy increase of the phase transition of this type of substance, which has been factored into the total enthalpy, expressed in J / kg. The equation is determined using the interval bisection method within the allowable temperature range of the glass melting furnace. .

[0120] In computer program implementation, the specific heat capacity function It is expressed in polynomial function form based on the fitting of differential scanning calorimetry experimental data, i.e. coefficients of each item , , Based on different material properties, the parameters are pre-stored in the melting simulation basic parameter set.

[0121] The liquid phase ratio is determined by the ratio of the sum of the mass of the initial liquid phase and the mass of the glass melt to the total mass of the material; the proportion of the mass of the reacted chemical components in each reaction stage to the mass of the initial reactive chemical components in that reaction stage is calculated to obtain the corresponding reaction progress.

[0122] Sum the cumulative transfer amounts of each refined sub-interface of the coarse-fine interface within a preset global time step, and then divide by . Based on the corresponding basic common surface area, obtain the heat flux, mass flux, and chemical component flux of the basic common surface, and replace the original flux of the same type to prevent double counting.

[0123] The equivalent reaction rate, equivalent heat transfer parameters, and equivalent migration parameters are calculated back from the cumulative reaction amount, cumulative heat transfer, and cumulative migration mass, and updated using the following formula:

[0124]

[0125] In the formula, 'a' represents the equivalent model parameter type number; , and These represent the updated parameters, the original parameters, and the parameters calculated from the results of local higher-order solutions, respectively. This represents the parameter update coefficient, with a value greater than 0 and not greater than 1.

[0126] The parameter update coefficients were calibrated using stable running samples: multiple candidate update coefficients were selected, and the local high-order solution results were mapped to the global low-order thermodynamic model. The joint residual strength of consecutive time steps after the update was compared, and the largest candidate update coefficient that could continuously reduce the joint residual strength without producing oscillations was determined. .

[0127] The global dynamic model refers to a global thermodynamic model formed by mapping local high-order solution results to corresponding basic spatial units based on a global low-order thermodynamic model, and updating the flux at the coarse-fine interface and the equivalent model parameters; the global dynamic model uses the global computational domain formed by S120. The formed global dynamic model is then input into S430.

[0128] S430. Continue to solve the global dynamic model for at least one global preset time step. According to the same residual definition, standardization benchmark, residual weight and residual tolerance value as S220 and S230, recalculate the joint residual field of the local upgrade region and its surrounding 1-level adjacent basic spatial units, and compare it with the result before the local upgrade item by item.

[0129] If the cumulative standardized residual of the dominant residual type in the local upgrade region is lower than the corresponding cumulative standardized residual before the upgrade, and no new local non-convergence marker is added, the local upgrade is considered valid. If the residual does not decrease and the current upgrade level is lower than the 3rd upgrade level, the upgrade level is increased by 1, the local upgrade configuration data is updated, and the process jumps to S330 to re-establish the local high-order thermodynamic model. If the residual still does not decrease after reaching the 3rd upgrade level, the integrity of the boundary data, the configuration of the reaction branches, and the local convergence status are checked. If there is still no improvement after 3 consecutive corrections, the simulation is marked as non-converged and does not proceed to S510.

[0130] When the adjacent basic space units on the outer side of a local upgraded region continuously exceed the limits, and their positions are consistent with the residual propagation direction, the local upgraded region is expanded; when the edge basic space units are continuously within the residual tolerance range, have no local non-convergence markers, and are not buffer space units necessary to maintain coarse and fine boundaries, shrinkage is performed; when the original abnormal space units, buffer space units, and peripheral associated space units continuously meet the residual tolerance values, the corresponding local high-order thermodynamic model is canceled, and the region is restored to the global low-order thermodynamic model calculation.

[0131] The number of continuous time steps is the value obtained by dividing the preset duration in S310 by the preset time step size of the entire domain and rounding it up. When the original abnormal spatial unit returns to normal but a new abnormality appears on the periphery along the material migration direction, the local upgrade region is not canceled, but rather moved or expanded.

[0132] Repeat steps S320 to S430 until the mass conservation residuals, energy conservation residuals, composition conservation residuals, and phase consistency residuals of all basic space units are within their respective allowable ranges, there are no global non-convergence markers or local non-convergence markers, and the preset duration is maintained continuously.

[0133] The global dynamic model that meets the above conditions is defined as a convergent global thermodynamic model. The region number, dominant residual type, upgrade level, expansion range, reduction range, cancellation time, residual changes before and after upgrade, and local upgrade computational cost of each local upgrade region are recorded in chronological order to form a local upgrade evolution record. The convergent global thermodynamic model and the local upgrade evolution record are then input into S510.

[0134] S5 specifically includes the following sub-steps:

[0135] S510: Read the convergent global thermodynamic model and local ascending evolution records, and extract temperature, chemical component content, reaction progress of each reaction stage, mass of each phase, liquid phase ratio, mass fraction of stagnant gas, mass flux and heat flux.

[0136] The melting front is defined as the position between a base space unit dominated by solid raw materials or reaction transition products and a base space unit dominated by primary liquid phase or molten glass, where the liquid phase ratio crosses the front liquid phase ratio. The front liquid phase ratio is calibrated by corresponding the liquid phase ratio of the model with the microscopic observation results of samples after interrupted melting and rapid cooling at different melting stages. When the liquid phase ratios of adjacent base space units are located on both sides of the front liquid phase ratio, linear interpolation is performed along the line connecting their centers to determine the melting front position.

[0137] Solid raw materials, reaction transition products, and residual solid particles with a D90 particle size larger than the determination size that have not completed the necessary reaction stages are defined as unmelted materials. The necessary reaction stages refer to the reaction stages that must be experienced in the process of converting raw materials into target glass melt, based on the types of raw materials, chemical composition content, and simultaneous thermal analysis results obtained from S110 in the current glass formulation; stages where corresponding raw materials are not present or where corresponding reaction peaks are not detected are not considered necessary reaction stages.

[0138] The D90 particle size indicates that 90% of the residual particle volume corresponds to a particle size no larger than this value. The residual particle size and the allowable proportion of unmelted material are derived from the glass product defect acceptance standards, enterprise melting process specifications, and laboratory melting tests and finished product stone detection results. The unmelted material mass, unmelted material mass fraction, and residual particle D90 particle size of each basic spatial unit are statistically analyzed to form the spatial distribution of unmelted material.

[0139] The state of complete melting is defined as the proportion of the total mass of unmelted material in the total dry basis of the input raw materials not exceeding the allowable proportion of unmelted material, the reaction progress of each necessary reaction stage not being lower than the corresponding reaction completion threshold, and the state being continuously maintained for a preset confirmation time. The time from the start of the simulation to the first entry into and continuous maintenance of the state of complete melting is defined as the complete melting time.

[0140] The completion threshold for each reaction is determined by the reaction endpoint of the corresponding reaction stage in the simultaneous thermal analysis, and the preset confirmation time is the 95th percentile of the duration of short-term fluctuations in the stable production sample. For example, if the overall liquid phase ratio reaches 0.98, but the mass fraction of high-silica residual particles is still 0.8%, and the allowable proportion of unmelted matter is 0.3%, it is not considered as complete melting.

[0141] The baseline operating strategy is composed of the dry basis feeding rate of each raw material, the heat load of each heating zone, and the feeding time interval corresponding to the convergence of S430. Based on the local order-up evolution record, the number of order-up times, the cumulative order-up time, the dominant residual type, and the amount of residual reduction before and after order-up are statistically analyzed for each basic spatial unit. Continuous spatial regions that repeatedly trigger local order-up or whose cumulative order-up time exceeds the set value are defined as sensitive kiln areas.

[0142] Positive and negative single-factor perturbations are applied to the total feeding rate, dry basis feeding rate of each raw material, heat load of the heating area and feeding time interval in the baseline operating strategy, respectively, with the perturbation amplitude taken as the corresponding equipment control resolution; the adjustable operating parameters that cause the peak changes of the combined residual strength of the sensitive kiln area, unmelted material mass or stagnant gas mass fraction to exceed the sensitivity tolerance value are used to form a set of sensitive factors.

[0143] The sensitivity tolerance value is taken as the 95th percentile of the corresponding change in the stable operating sample. The local order-up computational cost of the baseline operating strategy is statistically calculated based on the local order-up evolution record and denoted as... The melting simulation results and a set of sensitive factors are generated and input into S520.

[0144] S520. Generate candidate operating strategies based on the set of sensitive factors. A candidate operating strategy is a set of operating settings that can be issued to the raw material feeding equipment, combustion equipment, and electric melting equipment, including the dry basis feeding rate of each raw material, the heat load of each heating zone, and the feeding time interval. The upper and lower limits of each operating parameter are taken from the rated capacity of the equipment, the historical stable operating range, the allowable deviation of the target glass formulation, the safety operating specifications of the glass melting furnace, and the operating cycle of the feeding equipment; the adjustment step size is taken as the minimum control resolution of the corresponding equipment, and the parameter change rate does not exceed the ramp setting value allowed by the equipment controller.

[0145] Centered on the baseline operating strategy, sensitive factors are discretely increased, decreased, and combined according to preset adjustment step sizes to form initial candidate operating strategies; parameter combinations that violate constraints on equipment load, target formula ratio, total heat supply, or parameter change rate are deleted. When the total feed rate changes, the dry basis feed rate of each raw material is adjusted in conjunction with the target formula ratio. For example, when the total dry basis feed rate increases by 3%, the feed rates of quartz sand, soda ash, limestone, and dolomite are increased synchronously according to the original dry basis ratio; the target glass composition must not be changed by only increasing the soda ash feed rate.

[0146] Each executable candidate running strategy is written as a new dynamic modeling driver into the convergent global thermodynamic model, and the final state formed in S430 is used as the common initial state to re-execute S210 to S430. Candidate running strategies that do not meet the preset convergence conditions or whose dominant residuals have not decreased at the third ascending order are marked as invalid candidate running strategies and excluded.

[0147] For valid candidate operating strategies, record the complete melting time, total mass of unmelted material at the simulation endpoint, maximum local mass fraction of unmelted material, peak mass fraction of stagnant gas during the simulation period, maximum temperature standard deviation during the simulation period, total heat supply per unit mass of raw material, and local order-up calculations.

[0148] Temperature standard deviation is calculated in the nascent liquid phase region and the glass melt migration region; total heat supply per unit mass of raw material is the cumulative heat supply during the simulation period divided by the total dry mass of the input raw material; local order-scaling computation is the sum of the products of the number of refined spatial elements and the number of refined time steps in each local higher-order solution, and the local order-scaling computation corresponding to candidate running strategy s is denoted as... The values ​​are dimensionless. The above evaluation results are used to form a candidate strategy simulation evaluation set, which is then input into S530.

[0149] S530. Based on the glass product defect acceptance standards, the calibration results between bubble defects and the mass fraction of retained gas, the design temperature range of the glass melting furnace, the allowable temperature of refractory materials, the rated capacity of the equipment, and the allowable deviation of the target formula, determine the hard constraints for the total mass of unmelted material, the maximum local mass fraction of unmelted material, the peak value of the mass fraction of retained gas, the maximum temperature standard deviation, the heating load, the feeding capacity, and the formula deviation. Delete candidate operating strategies that violate any hard constraint or fail to meet the model convergence conditions.

[0150] Calculate the optimization evaluation value of the remaining candidate running strategies:

[0151]

[0152] In the formula, s represents the candidate running strategy number; This represents the optimization evaluation value of the candidate running strategy s, which is dimensionless. and These represent the complete melting time of the candidate running strategy s and the baseline running strategy, respectively, in seconds; and represents the local order-up computational cost of the candidate running strategy s and the baseline running strategy, respectively, both of which are dimensionless; and These represent the weights for complete melting time and local order-up computation, respectively, both of which are not less than 0 and their sum is 1.

[0153] The weights for complete melting time and local order-scaling computation are determined before strategy optimization. Let the allowable increase in complete melting time be the proportion of the baseline complete melting time. The allowable increase in local upgrade computational cost relative to the baseline local upgrade computational cost is: Then respectively for and Normalization is performed to obtain and The increment is allowed to be determined by production cycle time requirements and online computing resource limits. Therefore, indicators with smaller changes are allowed to receive larger weights, and these weights remain constant throughout this round of strategy optimization.

[0154] Among all feasible candidate operating strategies, the one with the smallest optimized evaluation value is selected as the initial operating strategy. When the difference in optimized evaluation values ​​among multiple candidate operating strategies does not exceed the preset evaluation tolerance, the strategy selected in the following order is chosen: smaller total mass of unmelted material and maximum local mass fraction; lower peak mass fraction of stagnant gas; smaller maximum temperature standard deviation; smaller total heat supply per unit mass of raw material; and smaller adjustment range relative to the baseline operating strategy. The preset evaluation tolerance is the 95th percentile of the difference in optimized evaluation values ​​obtained from repeated simulations of the same operating strategy.

[0155] The initial running strategy is subjected to complete and repeated simulations with the same initial state, equivalent model parameters, residual tolerance, and upgrade level as S520 to confirm that it meets the preset convergence conditions and all hard constraints, and that the deviation between the complete melting time and the local upgrade calculation amount and the first simulation result does not exceed the verification tolerance. If the verification fails, the second-best candidate running strategies are verified in turn.

[0156] The candidate operating strategies that have passed the review are determined as the target operating strategies. The output includes the dry basis feeding rate of each raw material, the heat load of each heating area, the feeding time interval, the execution time, the duration and the allowable rate of change. The output also includes the corresponding temperature field, chemical component content field, reaction progress field of each reaction stage, liquid phase ratio field, residual gas mass fraction field, melting front, spatial distribution of unmelted matter, complete melting time and constraint verification results.

[0157] When the actual weighing feeder uses the wet basis set value, the dry basis feeding rate is converted into the wet basis feeding rate based on the moisture content conversion relationship in S110 and then output, forming the simulation optimization results of the multi-component glass raw material melting process.

[0158] All the above formulas are performed using dimensionless numerical calculations; the relevant formulas are based on empirical models that approximate the real situation, obtained through extensive data collection and software simulation fitting. The preset parameters and thresholds involved in the formulas can be conventionally set and adjusted by those skilled in the art according to the physical constraints of the actual application scenario.

[0159] Those skilled in the art will recognize that the modules and algorithm steps of the various examples described in conjunction with the embodiments disclosed herein can be implemented in electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are implemented in hardware or software depends on the specific application and design constraints of the technical solution. Those skilled in the art can use different methods to implement the described functions for each specific application, but such implementation should not be considered beyond the scope of this application.

[0160] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

[0161] In conclusion, the above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A thermodynamic modeling and simulation method for the melting process of multi-component glass raw materials, characterized in that, Includes the following steps: S1. Obtain parameters of multi-component glass raw materials, glass melting furnace structure parameters and melting operation parameters, and establish a global computational domain and a global low-order thermodynamic model. S2. Perform time-progression solution on the global low-order thermodynamic model to obtain the global low-order state field, calculate the mass conservation residual, energy conservation residual, composition conservation residual and phase consistency residual of each basic space unit, and form residual and state correlation data. S3. Identify persistent anomaly clusters based on residual and state correlation data, determine local upgrade regions, dominant residual types and upgrade levels by combining buffer space units, and establish local high-order thermodynamic models based on local states, boundary states and mapping relationships between coarse and fine space units. S4. Solve the local higher-order thermodynamic model, conserve and map the local higher-order solution results to the corresponding basic space unit, update the flux at the coarse and fine interface and the equivalent model parameters, obtain the global dynamic model, and adjust the local order-up region according to the joint residual field until a convergent global thermodynamic model and local order-up evolution record are obtained.

2. The thermodynamic modeling and simulation method for the melting process of multi-component glass raw materials according to claim 1, characterized in that, Also includes: S5. Based on the convergent global thermodynamic model and local ascending evolution record, extract the melting front, spatial distribution of unmelted material and complete melting time, determine the sensitive kiln area and sensitive factor set, generate and simulate and evaluate candidate operating strategies, determine the target operating strategy from the candidate operating strategies that meet the process constraints and model convergence conditions, and output the operating setpoint and melting simulation results.

3. The thermodynamic modeling and simulation method for the melting process of multi-component glass raw materials according to claim 1, characterized in that, S1 specifically includes: Acquire parameters of multi-component glass raw materials, structural parameters of glass melting furnace, and melting operation parameters. Perform time alignment, unit conversion, and missing data processing on data of feeding, composition, thermal properties, reaction, and operation monitoring to form a basic parameter set for melting simulation. A three-dimensional computational domain is established based on the structural parameters of the glass melting furnace. Basic spatial units are divided and surface adjacency relationships are determined. Initial states, dynamic function labels, and boundaries for feeding, heating, discharging, and heat dissipation are set to form a global computational domain.

4. The thermodynamic modeling and simulation method for the melting process of multi-component glass raw materials according to claim 3, characterized in that, Also includes: Based on the basic parameter set of melting simulation and the global computational domain, sub-models of raw material heating, component reaction, phase transformation, gas generation and material migration are established and coupled to form a global low-order thermodynamic model.

5. The thermodynamic modeling and simulation method for the melting process of multi-component glass raw materials according to claim 1, characterized in that, S2 specifically includes: Read the global low-order thermodynamic model, the basic parameter set of melting simulation and the global computational domain, and perform time-progression solution of raw material heating, component reaction, phase transformation, gas generation and material migration to form a global low-order state field; Based on the global low-order state field, the mass conservation residual, energy conservation residual, composition conservation residual and phase consistency residual of each basic spatial unit are calculated to form multiple types of residual subfields; The residual subfields of various types are standardized, and the joint residual intensity is calculated by combining the residual weights. The number of residual co-occurrences, spatial continuity, residual duration and residual propagation direction are determined, and associated with the global low-order state field and dynamic function labels to form residual and state-related data.

6. The thermodynamic modeling and simulation method for the melting process of multi-component glass raw materials according to claim 1, characterized in that, S3 specifically includes: Based on residual and state correlation data, identify persistent abnormal clusters, set up buffer space units around them, and form a set of local upgrade regions. Calculate the residual contribution ratio in each local upgrade region, determine the dominant residual type and upgrade level, select the process to be refined and the reaction stage to be split in combination with dynamic function labels, extract the initial state, boundary state and coarse and fine spatial unit mapping relationship, and form local upgrade configuration data.

7. The thermodynamic modeling and simulation method for the melting process of multi-component glass raw materials according to claim 6, characterized in that, Also includes: Based on the local higher-order configuration data, the spatial units and time steps are refined, the extensive state variables are conserved and allocated, the actual reaction stages are broken down, and the relationships of heat transfer, component reaction, phase transformation, gas generation and material migration are established to form a local higher-order thermodynamic model.

8. The thermodynamic modeling and simulation method for the melting process of multi-component glass raw materials according to claim 1, characterized in that, S4 specifically includes: Solve the local high-order thermodynamic model by executing multiple refined time steps within a preset time step in the global domain, and obtain the local high-order state field and the local high-order solution results including the endpoint state, cumulative income and expenditure, cumulative flux at the coarse and fine interface, and local non-convergence markers; Based on the mapping relationship between coarse and fine spatial units, the local high-order solution results are conserved and mapped, the extensive state variables are summarized and the temperature is calculated inversely, the flux at the coarse and fine interface and the equivalent model parameters are updated, and a global dynamic model is formed. The joint residual field of the global dynamic model is recalculated. The local order-up region and order-up level are adjusted according to the residual changes. The local higher-order solution and conservation mapping are repeated until a convergent global thermodynamic model and local order-up evolution record are obtained.

9. The thermodynamic modeling and simulation method for the melting process of multi-component glass raw materials according to claim 2, characterized in that, S5 specifically includes: Read the convergent global thermodynamic model and local ascending evolution records, extract the melting front, spatial distribution of unmelted matter and complete melting time, and determine the sensitive kiln area and the set of sensitive factors; Based on the set of sensitive factors, candidate operating strategies that meet the process constraints are generated. Each candidate operating strategy is then re-executed as a dynamic modeling driving force. The melting index and local order-up computational cost are recorded to form a simulation evaluation set of candidate strategies.

10. The thermodynamic modeling and simulation method for the melting process of multi-component glass raw materials according to claim 9, characterized in that, Also includes: Based on the process hard constraints, a simulation evaluation set of candidate strategies is selected. The complete melting time and local order increase computation are comprehensively evaluated. The target operation strategy is verified and determined. The operation set value, melting simulation results and constraint verification results are output.

Images