Digital Twin Method for the Production Operation Status of Blast Furnaces

Abstract

The present invention provides a digital twin method for the production operation status of a blast furnace, which relates to the technical field of blast furnace production. Starting from the blast furnace top image, cross thermocouples, wall thermocouples, tuyeres, and measurable data at the top of the blast furnace, the present invention uses the boundary sub-model (Steps 1 and 2) as the boundary of the two-dimensional mathematical model of the blast furnace, analyzes the mechanism model of the blast furnace based on computational fluid dynamics, and the solved parameters (velocity, pressure, temperature, and composition) reflect the production operation status of the blast furnace, and the production operation status of the blast furnace can be updated according to the update of the measured data. Based on the mechanism and data, the physical model and dynamic changes are realized, and the reduced-order method is used to accelerate the solution speed, realizing the dynamic update of the blast furnace operation parameters to a certain extent, and proposing a digital twin method for establishing the production operation status of the blast furnace.

Description

Technical Field

[0001] The present invention relates to the technical field of blast furnace production, and particularly to a digital twin method for the operating state of blast furnace production. Background Art

[0002] Blast furnace digital twin is an important part of digital transformation. It combines the physical model of the blast furnace with computer simulation technology to establish a digital simulation model of the blast furnace, so as to realize real-time monitoring and optimal control of the operating state of the blast furnace. Traditional blast furnace production faces many problems, such as low production efficiency, difficult control of product quality, high energy consumption and emissions. These problems are closely related to the physical and chemical processes inside the blast furnace, and in-depth understanding of the working principle and operating rules inside the blast furnace is required to solve these problems. With the continuous development of information technology, digital transformation has become an important way to promote industrial upgrading and development. As an important part of digital transformation, digital twin technology can more accurately simulate and predict the physical and chemical processes inside the blast furnace by combining the establishment of physical models and computer simulation models, thereby improving the production efficiency and quality of the blast furnace.

[0003] The application of digital twin technology in blast furnaces has broad application prospects in improving blast furnace production efficiency and quality, reducing production costs, promoting blast furnace industrial upgrading and realizing sustainable development, etc., which are specifically reflected in: 1) Improving blast furnace production efficiency and quality. Blast furnace digital twin technology can real-time monitor the physical state and chemical reactions inside the blast furnace, and optimize the operation parameters and control strategies of the blast furnace by simulating and predicting the operation of the blast furnace, thereby improving the production efficiency and product quality of the blast furnace. 2) Reducing production costs. Through blast furnace digital twin technology, intelligent control of blast furnace operation can be realized, production costs can be reduced, and production benefits can be improved. At the same time, digital twin technology can also predict blast furnace failures, perform maintenance and repairs in a timely manner, and reduce downtime and maintenance costs. 3) Promoting blast furnace industrial upgrading. Blast furnace digital twin technology not only improves blast furnace production efficiency and quality, but also promotes blast furnace industrial upgrading and technological innovation. Through digital twin technology, the physical and chemical processes inside the blast furnace can be better understood, blast furnace design and operation control strategies can be optimized, and the technical level and core competitiveness of the blast furnace industry can be further improved. 4) Realizing sustainable development. Blast furnace digital twin technology can reduce the energy consumption and emissions of blast furnaces, improve the energy utilization rate and environmental protection level of blast furnaces, and realize sustainable development. Summary of the Invention

[0004] Aiming at the deficiencies of the prior art, the present invention provides a digital twin method for the operating state of blast furnace production.

[0005] A digital twin method for the operating state of blast furnace production includes the following steps:

[0006] Step 1: Calculate the temperature field of the blast furnace top based on the infrared image, perform gray correction on the infrared image, and use the lower surface of the temperature field as the boundary condition at the model furnace top;

[0007] The gray correction is as follows: Determine several gray levels (G0, G1, …, G n ) in the infrared image, where G0 is the lowest gray value and G n is the highest gray value. Use the image segmentation algorithm to extract the all-black and all-white regions of the image. Draw isotherms in the all-white and all-black regions according to the temperatures calculated by the cross thermometry. The number of isotherms m b in the all-black region and the number of isotherms m w in the all-white region are both related to the gray level difference, as shown below:

[0008]

[0009] In the formula, k is the gray difference of the isotherm;

[0010] Adjust the image gray level. For the all-white region, first set the gray level of the highest temperature region to G n , the lowest temperature region to G n-1 , and the gray levels of the intermediate temperature regions are gradually adjusted according to the temperature from high to low; for the all-black region, first set the gray level of the lowest temperature region to G0, the highest temperature region to G1, and the gray levels of the intermediate temperature regions are gradually adjusted according to the temperature from low to high; according to the position symmetry of the cross thermometry device and the quantity level of the reference temperature, divide the gray levels into corresponding several levels; count the gray values near the thermocouple in the image to determine the reference gray value corresponding to the thermocouple reference temperature; the temperature values between the reference grays are linearly interpolated between the corresponding reference temperatures by the two-point calibration method; take the average of the results of multiple calibrations to improve the calibration accuracy;

[0011] The description of the two-point calibration method is as follows. Assume the characteristics of temperature calibration:

[0012] y = ax + b (2)

[0013] where x is the gray input signal and y is the temperature value output signal. By finding the gain factor a and the offset factor b, the corresponding relationship between gray level and temperature can be obtained; let G min be the lowest gray value, and the corresponding temperature at this time is T min , G max be the highest gray value, and its corresponding temperature is T max . Substitute into the equation to get:

[0014]

[0015] Step 2: Calculate the inner surface temperature of the blast furnace using a mechanism model as the boundary condition at the model furnace wall;

[0016] First, establish the heat conduction differential equation of the blast furnace wall, determine the boundary conditions, establish a heat transfer model for the blast furnace wall and solve it to obtain the temperature distribution inside the wall;

[0017] First, establish the heat conduction differential equation of the blast furnace wall.

[0018] For the differential control volume, according to the law of conservation of energy:

[0019] Input - Output + R = S (4)

[0020] In the formula, Input is the sum of the heat input into the control volume per unit time, Input = Input x + Input y + Input z , where Input x = q x dydz, Input y = q y dxdz, Input z = q z dxdy are the heat input into the control volume from the x, y, and z directions per unit time respectively; Output is the sum of the heat output from the control volume per unit time, Output = Output x + Output y + Output z , where are the heat output from the control volume in the x, y, and z directions per unit time respectively; R is the heat generated inside the control volume, that is, the internal heat source. S is the change in internal energy of the control volume over time after absorbing or releasing heat;

[0021] Therefore, there is:

[0022]

[0023] There is Fourier's law:

[0024]

[0025] where q is the heat conduction flux, is the temperature gradient, and λ is the thermal conductivity;

[0026] Since the furnace wall is a solid and the only heat transfer method inside it is conduction, substituting Fourier's law (6) into the heat conduction differential equation (5) gives:

[0027]

[0028] Therefore, it is transformed into the heat conduction differential equation in the cylindrical coordinate system. From the heat conduction differential equation in the rectangular coordinate system, through the mathematical coordinate transformation \(x = r\cos\theta\), \(y = r\sin\theta\), \(z = h\), we get

[0029]

[0030] where \(r\) is the cross-sectional radius, \(\theta\) is the cross-sectional angle, \(h\) is the longitudinal axis height, and \(T\) is the temperature distribution of the furnace wall;

[0031] In order to obtain the circumferential temperature distribution at different heights on the inner surface of the blast furnace, first establish a heat transfer mechanism model for the cross-section of the blast furnace wall. Without considering the axial distribution, we get:

[0032]

[0033] The above heat conduction differential equation is discretized by central difference in the radial and circumferential directions, and the discretized difference equation is as follows:

[0034]

[0035] where \(r_0\) is the radius of the inner ring of the ring model, \(r_2\) is the radius of the outer ring of the ring model, \(\Delta r\) is the radial step size, \(\Delta\theta\) is the circumferential step size, \(r\) i is the radius at the \(i\)th position, \(M\) is the number of radial differential steps, and \(N\) is the number of circumferential differential steps;

[0036] Inner surface boundary of the furnace wall: \(T\) 1,j =\(T\) in,j , \(r = r_0\), \(0\leq j\leq N\) (11)

[0037] Boundary at the junction of the refractory brick and the cooling stave:

[0038]

[0039] After discretization:

[0040] Outer surface boundary of the cooling stave:

[0041] After discretization:

[0042] where: \(T\) in , \(T\) out are the temperatures of the inner surface of the blast furnace and the outer surface of the cooling stave; \(T\) CW is the average temperature of the cooling water around the cooling stave, \(T\) CW =( \(T\) CW,in +\(T\) CW,out ) / 2, \(T\) CW,in is the inlet temperature of the cooling water, \(T\) CW,outis the outlet temperature of the cooling water; r0, r1, and r2 are the radii at the inner surface of the blast furnace refractory brick, the junction of the refractory brick and the cooling stave, and the cooling pipe; λ1 and λ2 are the thermal conductivities of the refractory brick and the cooling stave respectively; h x is the convective heat transfer coefficient between the cooling stave and the cooling water, and its value is calculated according to the following empirical formula, where v 水 is the flow velocity of the cooling water:

[0043] h x = 208 + 47.5v 水 (16)

[0044] Taking the known temperature of the thermocouple on the furnace wall as the comparison value, using the above equations (9) and (10) and the boundary condition equations (11) to (15) as the calculation model, with the cooling water pipe as the outer boundary and the inner surface as the inner boundary, the temperature field of the furnace wall is solved; assuming the inner surface temperature to solve the thermocouple temperature, calculating the error between the true value and the calculated value of the thermocouple temperature, when the error is satisfied, the assumed inner surface temperature is the true inner surface temperature;

[0045] Step 3: Taking the blast furnace shaft mass point as the research object, analyze and establish the basic model of the blast furnace shaft temperature field, use the N-S equation described by Euler to solve the physical quantities of each fluid particle at different positions and different times in the entire flow field domain, including pressure, velocity, and temperature, and describe the physical state of the entire flow field domain through these physical quantities; divide the inside of the blast furnace into countless fluid micro-elements, analyze a single fluid micro-element, and obtain the basic control equation set of the blast furnace model;

[0046] Based on the continuous medium hypothesis, according to the three conservation laws of mass, momentum, and energy, taking the control volume described by Euler as the research object, without considering the body force and source term, analyze the control volume; the increase in the momentum of the control volume is equal to the momentum flowing in through the control surface, for the control volume there is:

[0047]

[0048] In the formula, ρ is the gas density; △x, △y, and △z are the sizes of the differential control volume in three directions respectively; △t is the time step; are the momentum changes in three directions.

[0049] Transform the momentum conservation equation into:

[0050]

[0051] In the formula: μ is the dynamic viscosity coefficient;

[0052] Similarly, the mass conservation equation is:

[0053]

[0054] Similarly, the energy conservation equation is:

[0055] where: k is the thermal conductivity of Fourier's law; c v is the specific heat;

[0056] Step 4: Analyze the transfer phenomenon between the heterogeneous and homogeneous phases existing in the blast furnace as the source term; specifically including gas-solid momentum transfer and gas-solid energy transfer;

[0057] The gas-solid momentum transfer: Analyze the relationship between the multiphase flows in the blast furnace, and uniformly describe the momentum transfer amount per unit time between each phase state in the form of the following equation:

[0058] F i-j = f i-j (U i - U j ) (21)

[0059] In the formula, F i-j is the momentum transfer amount per unit time between the i-th phase and the j-th phase; f i-j is the momentum transfer coefficient between the i-th phase and the j-th phase; U i , U j are the velocities of the i-th phase and the j-th phase respectively;

[0060] Assume that the inclination angle of each layer of solid phase is θ, and the momentum transfer phenomenon between the solid phase and the gas phase is described as follows:

[0061]

[0062] In the formula: R perp is the flow resistance coefficient in the direction perpendicular to the burden surface, and is expressed by the following formula:

[0063]

[0064] R para is the flow resistance coefficient in the direction parallel to the burden surface, and is expressed by the following formula:

[0065]

[0066] R k is the resistance coefficient of different components of the material, and is expressed by the Ergun equation:

[0067]

[0068] f k is the volume fraction of different components of the material;

[0069] Substitute Eqs. (23)-(25) into (22), assuming the inclination angle is 0, we get:

[0070]

[0071] where F g-s is the momentum transfer amount per unit time between the gas and the solid; η g is the gas viscosity; u g is the gas velocity vector; u s is the solid velocity vector; ρ g is the gas density; ε s is the porosity of the burden; ψ s is the shape factor of the solid particles; d s is the solid particle diameter.

[0072] The gas-solid energy transfer: The energy transfer amount between phases is uniformly described in the form of the following equation:

[0073] F i-j = A i-j h i-j (T i - T j ) (27)

[0074] where F i-j is the energy transfer amount between phase i and phase j; A i-j is the contact area; h i-j is the energy transfer coefficient between phase i and phase j; T i and T j are the temperatures of phase i and phase j, respectively;

[0075] Calculate the energy transfer coefficient according to the modified Ranz-Marshall formula, and the energy transfer phenomenon between the gas phase and the solid phase is described as follows:

[0076]

[0077] where k g is the thermal conductivity of the gas; d s is the average particle size; Re is the Reynolds number, which is calculated according to the following formula:

[0078]

[0079] where η g is the gas viscosity; u g is the gas velocity vector; u s is the solid velocity vector; ρ g is the gas density; ψ s is the shape factor of the solid particles; d s is the solid particle diameter,

[0080] Prandtl number Pr of convective heat transfer of gas g A dimensionless combination number representing the mutual influence of the energy and momentum transfer processes in a fluid. The relationship between the surface temperature boundary layer and the flow boundary layer is calculated according to the following formula, where C g is the specific heat capacity of the gas:

[0081]

[0082] Step 5: Since the blast furnace is a symmetric model, a two-dimensional mathematical model of the blast furnace is established;

[0083] By introducing the source term in Step 4 into the basic control equations analyzed in Step 3, a two-dimensional mathematical model of the blast furnace is established, and the specific continuity equation, momentum conservation equation, and energy conservation equation of the blast furnace are established to specify the model, and a two-dimensional mathematical model of the blast furnace for the blast furnace process is obtained. The two-dimensional mathematical model of the blast furnace includes a gas flow model, a gas temperature model, and a solid temperature model; the magnitude of the body force acting on the fluid volume element in the gas flow model is proportional to the volume of the fluid element, and f x 、f y 、f z are used to represent the body force per unit mass of the fluid. After considering the body force and converting the three-dimensional to two-dimensional, the momentum conservation equation (18) is reduced to:

[0084]

[0085] The momentum transfer formula (26) from the gas phase to the solid phase is reduced to:

[0086]

[0087] Considering the momentum transfer formula (32) between the gas and the solid in the source term of the gas momentum conservation control equation, formula (31) is reduced to:

[0088]

[0089] In the formula:

[0090] μ is the dynamic viscosity coefficient.

[0091]

[0092] At the same time, the gas flow model also includes a continuity equation. Reducing the continuity equation (19) to two dimensions gives:

[0093]

[0094] The gas temperature model is established according to the energy conservation equation (20), and at the same time, the heat transfer between phases and the heat of chemical reaction are considered in the source term;

[0095] Considering the heat transfer formula (28) between gas and solid in the source term of the energy conservation control equation, Equation (20) is transformed into:

[0096]

[0097] Where:

[0098] k is the thermal conductivity of Fourier's law, c v is the specific heat:

[0099]

[0100] k g is the thermal conductivity of the gas, d s is the average particle size, Re is the Reynolds number, Pr g is the convective heat transfer Prandtl coefficient of the gas, representing the dimensionless combination number of the mutual influence of the energy and momentum transfer processes in the fluid, and the relationship between the surface temperature boundary layer and the flow boundary layer.

[0101] The heat transfer between solids in the solid temperature model is described by the heat conduction equation:

[0102]

[0103] Where: k is the thermal conductivity, which is determined by the thermal conductivity, density and heat capacity of the material;

[0104] Therefore, the solid temperature model is as follows:

[0105]

[0106] Where:

[0107]

[0108] k g is the thermal conductivity of the gas, d s is the average particle size, Re is the Reynolds number, Pr g is the convective heat transfer Prandtl coefficient of the gas, representing the dimensionless combination number of the mutual influence of the energy and momentum transfer processes in the fluid, and the relationship between the surface temperature boundary layer and the flow boundary layer.

[0109] The reaction heat of the burden itself is obtained by substituting the empirical formula:

[0110]

[0111] Where:

[0112] D pis the diameter of the ore balls in the furnace. Q1, Q2, E, R, and K1 are all constants related to the compositions of the ore and coke.

[0113] Step 6: Solve the 2D mathematical model of the blast furnace using the finite difference method. Since the finite difference method solves on a grid, it is necessary to first divide the modeling area into grids and then discretely solve the model on the grid.

[0114] Use the BFC grid generation technique to obtain the body-fitted coordinates. Through coordinate transformation, the irregular region in the physical space is mapped one-to-one with the regular region in the computational space, and the discretization of the partial differential equations in the model is transformed from the original irregular region to the regular region in the computational space.

[0115] Step 6.1: Generate the grid for the 2D blast furnace modeling area.

[0116] Step 6.1.1: First, set the number of grids and the boundary conditions of the physical region. Use the linear interpolation method to divide the grids in the physical region, and the division result is used as the initial value for the calculation iteration.

[0117] Step 6.1.2: Discretize the elliptic differential equations (45) to (47). According to the coordinate distribution at the current time, obtain the values of the coefficients J, α, β, and γ.

[0118]

[0119]

[0120] P and Q are source terms, specifically:

[0121]

[0122] Step 6.1.3: Calculate the source terms: First, use Equation (48) to calculate the values of φ(i,j) and ψ(i,j) on the boundary conditions:

[0123]

[0124] Then, use linear interpolation to obtain φ(i,j) and ψ(i,j) at all points. Finally, obtain P(i,j) and Q(i,j) according to Equation (47).

[0125] Step 6.1.4: After obtaining all the parameters, calculate the grid coordinate distribution in the physical plane according to Equation (45) and determine whether it meets the set accuracy requirements. If it meets, generate the grid diagram; otherwise, return to Step 6.1.2.

[0126] Step 6.2: Solve the 2D mathematical model of the blast furnace.

[0127] The finite difference method is selected to discretize the model, and the staggered grid method is selected. The time-direction iterative solution is carried out on the grid generated by the BFC grid generation technology in step 6.1;

[0128] Step 6.2.1: Determine that the blast furnace tuyere is the gas inlet and the top is the gas outlet. According to the parameters of the gas inlet and the gas outlet collected on-site, these are the boundary conditions at these two locations; the furnace wall and the bottom are set with no-slip boundary conditions, and the furnace wall temperature is obtained in step 2; the symmetry axis is set with the boundary condition that the gradients of all parameters are 0, and the internal parameters of the blast furnace in the steady state are used as the initial conditions;

[0129] Step 6.2.2: Calculate the velocity at the next moment through the gas flow model and the gas velocity and pressure distributions at the current moment; and

[0130] Step 6.2.3: Calculate the gas temperature at the next moment through the gas temperature model and the gas velocity, gas temperature, and solid temperature distributions at the current moment;

[0131] Step 6.2.4: Calculate the solid temperature at the next moment through the solid temperature model and the gas temperature and solid temperature distributions at the current moment;

[0132] Step 6.2.5: Judge whether the change in the temperature value at the current moment and the temperature value at the next moment meets the accuracy requirement. If it meets, go to step 7; if it does not meet, return to step 6.2.2.

[0133] Step 7: POD reduced-order optimization;

[0134] Step 7.1: First, use step 6.2 to solve for the L time period and then stop, and give the approximate solutions of the first L time periods to form a snapshot matrix; and These five snapshot matrices respectively represent the radial velocity, circumferential velocity, pressure, gas temperature, and solid temperature in the two-dimensional modeling area of the blast furnace in the first L time periods. Where m = MN is the total number of grid points in the modeling area, 1 ≤ i ≤ m, 1 ≤ l ≤ L.

[0135] Step 7.2: Solve the following linear equations according to the above five snapshot matrices:

[0136]

[0137] Obtain five sets of eigenvalues of the above linear equations;

[0138]

[0139] and the corresponding five sets of eigenvectors

[0140]

[0141] Step 7.3: On the premise of satisfying the POD basis error estimation condition , determine the acceptable error e = O(△t, △x 2 , △y 2 ) of the 2D mathematical model of the blast furnace in the POD optimization stage. The order of the POD basis is the model reduction order and constitute the initial POD basis:

[0142]

[0143] Step 7.4: Solve the POD eigenvalues; assume the solutions of the entire blast furnace shaft temperature model are and Then, they converge to

[0144]

[0145] According to the 2D mathematical model of the blast furnace, we have

[0146]

[0147] where is obtained according to the 2D mathematical model of the blast furnace;

[0148] Let

[0149] Calculate the POD eigenvalues of each of the first L segments based on the data of the first L time periods and the POD basis

[0150]

[0151] Calculate the POD eigenvalues of each segment from the L - th time period to the (T - 1)-th time period

[0152]

[0153] where

[0154]

[0155] Calculate the POD eigenvalues of each segment from the L - th time period to the convergence end time period using equations (56) and (57).

[0156] Step 7.5: Solve the model distribution;

[0157] Based on the POD basis and POD eigenvalues obtained in the above steps 7.3 and 7.4, we get from equation (58)

[0158]

[0159] Use to replace and the target value distribution of the internal temperature model of the blast furnace can be obtained and First, calculate the differences between the pressure, velocity, and temperature at the current moment and the previous moment to determine whether the iteration is completed. If the termination accuracy is satisfied, the iteration is completed and the program ends. If not, first, according to Step 7.6, determine whether the POD basis obtained at the current moment meets the POD error requirement. If not, perform Step 7.6 to update the POD basis. If it meets the requirement, calculate the distribution at the next moment based on the POD eigenvalues and POD basis at the next moment until the termination accuracy is satisfied

[0160] Step 7.6: Use to replace is the main error source of this POD method. Analyze the error and determine the update of the POD basis

[0161] Let the POD error be

[0162]

[0163] According to Equation (52), we have

[0164]

[0165]

[0166] According to Equation (52), there is

[0167]

[0168] Let On the premise of ensuring the stability of the difference equation, according to Equation (61), we have

[0169]

[0170] Under the condition of ensuring the stability of the differential equation, there is the following error estimate between the accurate solution of the blast furnace temperature distribution model and the approximate solution of the POD reduction algorithm

[0171]

[0172] where

[0173]

[0174]

[0175] Therefore, in step 7.5, for each time period, let Judge Whether it holds. If it holds, then And Are the solutions that meet the accuracy requirements for the current time period. Otherwise, repeat steps 7.1 to 7.4 until the And That meet the above accuracy are obtained, and then continue with step 7.5;

[0176] Step 8: Use the results obtained in steps 6 and 7 as the calculation basis for the reaction rates in all chemical reactions in the blast furnace. Calculate the enthalpy change during the chemical reaction process based on the chemical reaction rate and the composition distribution generated by the flow, and perform temperature field correction and coupling modeling based on the enthalpy change;

[0177] Step 9: Temperature field optimization and composition concentration. Calculate the chemical reaction rates of each reaction based on the temperature field, and then calculate the enthalpy change △H and the change in CO gas flow concentration based on the chemical reaction rate. Use the enthalpy change to correct the temperature distribution and composition distribution, and finally realize the digital twin of the temperature field and composition field during the operation of the blast furnace.

[0178] Step 9.1: Assume the internal composition distribution of the blast furnace without considering chemical reactions based on the velocity field in the blast furnace modeling area obtained in steps 6 and 7.

[0179] Step 9.2: Calculate the chemical reaction rates of each chemical reaction in step 8 based on the temperature field in the blast furnace modeling area obtained in steps 6 and 7.

[0180] Step 9.3: Calculate the enthalpy change and composition change of each chemical reaction based on the chemical rate;

[0181] Step 9.4: Update the temperature field and composition field based on the enthalpy change and composition change.

[0182] Step 9.5: Judge whether the coupling degree error between the temperature field and the composition field meets the accuracy. If it meets, end. Otherwise, continue with step 9.1.

[0183] The beneficial effects produced by adopting the above technical solutions are as follows:

[0184] The present invention provides a digital twin method for the production operation state of a blast furnace. Starting from the top-of-furnace image, cross thermocouples, wall thermocouples, tuyeres, and measurable data at the top of the blast furnace, using the boundary sub-model (Steps 1, 2) as the boundary of the two-dimensional mathematical model of the blast furnace, analyzing the mechanism model of the blast furnace based on computational fluid dynamics, the solved parameters (velocity, pressure, temperature, and composition) reflect the production operation state of the blast furnace, and the production operation state of the blast furnace can be updated according to the update of the measured data. Based on mechanism and data, the physical model and dynamic changes are realized, and a reduced-order method is used to accelerate the solution speed, realizing the dynamic update of the blast furnace operation parameters to a certain extent, and a digital twin method for establishing the production operation state of the blast furnace is proposed. BRIEF DESCRIPTION OF THE DRAWINGS

[0185] Figure 1 It is a schematic cross-sectional view of the blast furnace wall of the present invention;

[0186] Figure 2 It is the mass conservation differential control volume of the present invention;

[0187] Figure 3 It is a schematic diagram of the sectionalization of the blast furnace body of the present invention;

[0188] Figure 4 It is a flow chart of network generation of the present invention;

[0189] Figure 5 It is a flow chart of the temperature field accuracy of the present invention;

[0190] Figure 6 It is a schematic diagram of the chemical reaction in the blast furnace of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0191] The following will further describe in detail the specific embodiments of the present invention in conjunction with the drawings and embodiments. The following embodiments are used to illustrate the present invention, but are not used to limit the scope of the present invention.

[0192] A digital twin method for the production operation state of a blast furnace includes the following steps:

[0193] Step 1: Calculate the temperature field of the blast furnace top based on the infrared image, perform gray correction on the infrared image, and use the lower surface of the temperature field as the boundary condition at the model furnace top; it can be known from the characteristics of the blast furnace infrared image that the gray level of the outer image is positively correlated with the temperature of the burden surface. The gray level of the image at the place where the burden surface temperature is high is relatively high, and the gray level of the image at the place where the temperature is low is relatively low. The measurement range of the infrared camera system is between 200°C and 600°C. The gray level corresponding to the position of the image where the burden surface is higher than 600°C is 255, and the gray level of the image where it is lower than the Celsius degree is 0. The bright spot at the center of the blast furnace corresponds to a gray level of 255, which is called the saturation area; among them, the black area near the furnace wall corresponds to a gray level of 0, which is called the detection dead area. Due to the limitation of the measurement range of the camera, the gray level information in the saturation area and the dead area in the image cannot effectively reflect the temperature of the burden surface field. Moreover, considering the influence of the cross temperature measurement beam on the infrared image, the gray level value in the area blocked by the cross temperature measurement beam cannot reflect the temperature information of the burden surface field either. Therefore, in this embodiment, the gray level of the infrared image is corrected and recalibrated.

[0194] Since the blast furnace burden distribution adopts circular burden distribution, generally, the temperature at the center of the burden surface image is high, and the temperature at the periphery is low, and it has the characteristic of central symmetry. Further considering that the gray level value of the burden surface infrared image is another measure of the burden surface temperature field, and its gray level value distribution has the same change trend as the temperature field distribution. Therefore, the gray level value of the infrared image is corrected and recalibrated by using the obtained cross temperature measurement interpolation curve.

[0195] The gray correction is as follows: Determine several gray levels (G0, G1,..., G n ) in the infrared image, where G0 is the lowest gray level value, and G n is the highest gray level value. Use the image segmentation algorithm to extract the all-black area and the all-white area of the image. Draw isothermal lines in the all-white and all-black areas according to the temperature calculated by the cross temperature measurement. The number of isothermal lines m b in the all-black area and the number of isothermal lines m w in the all-white area are both related to the gray level difference; adjust the gray level of the image. For the all-white area, first set the gray level of the highest temperature area to G n , the lowest temperature area to G n-1 , and the gray level of the intermediate temperature area is adjusted gradually according to the temperature from high to low; for the all-black area, first set the gray level of the lowest temperature area to G0, the highest temperature area to G1, and the gray level of the intermediate temperature area is adjusted gradually according to the temperature from low to high; according to the position symmetry of the cross temperature measurement device and the number level of the reference temperature, divide the gray level into corresponding several levels; count the gray level values near the thermocouple in the image to determine the reference gray level value corresponding to the reference temperature of the thermocouple; the temperature values between the reference gray levels are linearly interpolated between the corresponding reference temperatures by using the two-point calibration method; take the average value of the results of multiple calibrations to improve the accuracy of the calibration;

[0196] The two-point calibration method is described as follows. Assume the characteristics of temperature calibration:

[0197] y = ax + b (2)

[0198] where x is the grayscale input signal, y is the temperature value output signal. By obtaining the gain factor a and the offset factor b, the corresponding relationship between grayscale and temperature can be obtained. Let G min be the lowest grayscale value, and the corresponding temperature at this time is T min , G max be the highest grayscale value, and its corresponding temperature is T max . Substituting into the equation, we get:

[0199]

[0200] Step 2: Use the mechanism model to calculate the inner surface temperature of the blast furnace as the boundary condition at the model furnace wall; heat transfer occurs from the inner surface of the blast furnace wall to the furnace shell. Therefore, to obtain the inner surface temperature field, it is necessary to establish a heat conduction differential equation for the furnace wall structure. First, establish the heat conduction differential equation of the blast furnace wall, determine the boundary conditions, establish a heat transfer model for the blast furnace wall and solve it to obtain the temperature distribution inside the furnace wall;

[0201] First, establish the heat conduction differential equation of the blast furnace wall.

[0202] For the differential control volume, according to the law of conservation of energy, we have:

[0203] Input - Output + R = S (4)

[0204] In the formula, Input is the sum of the heat input into the control volume per unit time, Input = Input x + Input y + Input z , where Input x = q x dydz, Input y = q y dxdz, Input z = q z dxdy are the heat input into the control volume from the x, y, and z directions per unit time respectively; Output is the sum of the heat output from the control volume per unit time, Output = Output x + Output y + Output z , where are the heat output from the control volume in the x, y, and z directions per unit time, respectively; R is the heat generated inside the control volume, i.e., the internal heat source. For the furnace wall model in this embodiment, this value is 0; S is the change in internal energy with time after the control volume absorbs or releases heat. In the blast furnace wall model of this embodiment, the environment is after the blast furnace operates stably, its state is stable and does not fluctuate with time, and this value is 0.

[0205] Therefore, we have:

[0206]

[0207] There is Fourier's law:

[0208]

[0209] where q is the heat conduction flux, is the temperature gradient, and λ is the thermal conductivity; since the furnace wall is solid, the heat transfer mode inside it is only heat conduction. Substituting Fourier's law (6) into the heat conduction differential equation (5), we get:

[0210]

[0211] Therefore, converting it into the heat conduction differential equation in the cylindrical coordinate system, through the mathematical coordinate transformation x = r cosθ, y = r sinθ, z = h from the heat conduction differential equation in the rectangular coordinate system, we get

[0212]

[0213] In the formula, r is the cross-sectional radius, θ is the cross-sectional angle, h is the longitudinal axis height, and T is the temperature distribution of the furnace wall; in order to obtain the circumferential temperature distribution at different heights on the inner surface of the blast furnace, first establish a heat transfer mechanism model of the cross-section of the blast furnace wall, without considering the axial distribution, we get:

[0214]

[0215] Discretize the above heat conduction differential equation by central difference in the radial and circumferential directions, and the discretized difference equation is as follows:

[0216]

[0217] In the formula, r0 is the radius of the inner ring of the ring model, r2 is the radius of the outer ring of the ring model, △r is the radial step, △θ is the circumferential step, r i is the radius at the i-th position, M is the number of radial differential steps, and N is the number of circumferential differential steps;

[0218] For the furnace wall, the temperature is uniform near the junction of the cooling stave and the cooling pipe. Forced convective heat transfer occurs with the cooling water, which can be considered as the third type of boundary condition. The inner surface of the furnace wall is in contact with the inside of the blast furnace, and the temperature at this point can be considered the same as the temperature of the gas inside the blast furnace, which is treated as the first type of boundary condition. Therefore, the boundary conditions of this model are determined as follows: Inner surface boundary of the furnace wall:

[0219] T 1,j =T in,j , r=r0, 0≤j≤N (11)

[0220] Boundary at the junction of the refractory brick and the cooling stave:

[0221]

[0222] After discretization:

[0223]

[0224] Outer surface boundary of the cooling stave:

[0225] After discretization:

[0226] In the formula: T in 、T out are the temperatures of the inner surface of the blast furnace and the outer surface of the cooling stave; T CW is the average temperature of the cooling water around the cooling stave, T CW =(T CW,in +T CW,out ) / 2, T CW,in is the inlet temperature of the cooling water, T CW,out is the outlet temperature of the cooling water; r0, r1, r2 are the radii of the inner surface of the blast furnace refractory brick, the junction of the refractory brick and the cooling stave, and the cooling pipe; λ1 and λ2 are the thermal conductivities of the refractory brick and the cooling stave respectively; h x is the convective heat transfer coefficient between the cooling stave and the cooling water, and its value is calculated according to the following empirical formula, where v 水 is the flow velocity of the cooling water:

[0227] h x =208 + 47.5v 水 (16)

[0228] With the known thermocouple on the furnace wall (inside the furnace wall, see Figure 1)Taking temperature as the comparison value, using the above equations (9) and (10) and the boundary condition equations (11) to (15) as the calculation model, with the cooling water pipe as the outer boundary and the inner surface as the inner boundary, the temperature field of the furnace wall is solved; since what needs to be solved in this step is the inner surface temperature field, it is a classical thermodynamic inverse problem. Assume the inner surface temperature to solve the thermocouple temperature, calculate the error between the true value and the calculated value of the thermocouple temperature, and when the error is satisfied, the assumed inner surface temperature is the true inner surface temperature. The PSO algorithm is used for optimization here, and the initial population and population update are set using the PSO idea, the error is calculated, and the optimal population (i.e., the true inner surface temperature) is selected. The PSO will not be introduced further here.

[0229] Step 3: Taking the mass points of the blast furnace shaft as the research object, analyze and establish the basic model of the blast furnace shaft temperature field, use the N - S equations described by Euler to solve the physical quantities (including pressure, velocity, temperature) of each fluid particle at different positions and different times in the entire flow field domain, and describe the physical state of the entire flow field domain through these physical quantities; divide the inside of the blast furnace into countless fluid micro - groups, starting from the most basic conservation laws in nature (mass conservation, momentum conservation, energy conservation), analyze a single fluid micro - group, and obtain the basic control equation set of the blast furnace model as equations (18) to (20); the mass conservation differential control volume is as Figure 2 shown;

[0230] Based on the continuous medium hypothesis, according to the three major conservation laws of mass, momentum, and energy, taking the control volume described by Euler as the research object, without considering the body force and source term, analyze the control volume; the increase in the momentum of the control volume is equal to the momentum flowing in through the control surface. According to Figure 1 the conservation differential control volume, for the control volume, there is:

[0231]

[0232] In the formula, ρ is the gas density; △x, △y, and △z are the sizes of the differential control volume in three directions respectively; △t is the time step; are the momentum change amounts in three directions. The momentum conservation equation is transformed into:

[0233]

[0234] In the formula: μ is the dynamic viscosity coefficient. Similarly, analyzing mass according to the above - mentioned momentum analysis, the mass conservation equation and the energy conservation equation are obtained as:

[0235]

[0236]

[0237] Step 4: Analyze the transfer phenomena between the heterogeneous and homogeneous phases existing in the blast furnace as the source term; use it as the main basis for establishing the blast furnace multiphase flow model in Step 5. Specifically, it includes gas-solid momentum transfer and gas-solid energy transfer;

[0238] The gas-solid momentum transfer: Based on the overall framework of the three-dimensional mathematical model of the entire blast furnace, analyze the relationship between the multiphase flows in the blast furnace, and uniformly describe the momentum transfer amount per unit time between each phase state in the form of the following equation:

[0239] F i-j =f i-j (U i -U j ) (21)

[0240] In the formula, F i-j is the momentum transfer amount per unit time between the i-phase and the j-phase; f i-j is the momentum transfer coefficient between the i-phase and the j-phase; U i and U j are the velocities of the i-phase and the j-phase respectively;

[0241] Assume that the inclination angle of each layer of solid phase is θ, and the momentum transfer phenomena between the solid phase and the gas phase are described as follows:

[0242]

[0243] In the formula: R perp is the flow resistance coefficient in the direction perpendicular to the burden surface, and is expressed by the following formula:

[0244]

[0245] R para is the flow resistance coefficient in the direction parallel to the burden surface, and is expressed by the following formula:

[0246]

[0247] R k is the resistance coefficient of different components of the material, and is expressed by the Ergun equation, where f k is the volume fraction of different components of the material:

[0248]

[0249] Substitute equations (23)-(25) into (22), assuming the inclination angle is 0, we get:

[0250]

[0251] In the formula, F g-s is the momentum transfer amount per unit time between the gas and the solid; η g is the gas viscosity; ug is the gas velocity vector; u s is the solid velocity vector; ρ g is the gas density; ε s is the porosity of the burden; ψ s is the shape factor of solid particles; d s is the solid particle diameter. The gas-solid energy transfer: Similar to the above gas-solid momentum transfer, the amount of energy transfer between phases is uniformly described in the form of the following equation:

[0252] F i-j = A i-j h i-j (T i - T j ) (27)

[0253] where F i-j is the amount of energy transfer between phase i and phase j; A i-j is the contact area; h i-j is the energy transfer coefficient between phase i and phase j; T i , T j are the temperatures of phase i and phase j respectively; The energy transfer coefficient is calculated according to the modified Ranz-Marshall formula, and the energy transfer phenomenon between the gas phase and the solid phase is described as follows:

[0254]

[0255] where k g is the thermal conductivity of the gas; d s is the average particle size; Re is the Reynolds number, which is calculated according to the following formula:

[0256]

[0257] where η g is the gas viscosity; u g is the gas velocity vector; u s is the solid velocity vector; ρ g is the gas density; ψ s is the shape factor of solid particles; d s is the solid particle diameter. The convective heat transfer Prandtl number Pr of the gas g represents a dimensionless combination number that reflects the mutual influence of the energy and momentum transfer processes in the fluid. The relationship between the surface temperature boundary layer and the flow boundary layer is calculated according to the following formula, where C g is the specific heat capacity of the gas:

[0258]

[0259] Step 5: Since the blast furnace is a symmetric model, a two-dimensional mathematical model of the blast furnace is established. By introducing the source term of the control equation analyzed in Step 4 into the basic control equation set analyzed in Step 3, a two-dimensional mathematical model of the blast furnace is established, and specific continuity equation, momentum conservation equation and energy conservation equation of the blast furnace are established to specify the model, and a two-dimensional mathematical model of the blast furnace for the blast furnace process is obtained. The two-dimensional mathematical model of the blast furnace includes a gas flow model, a gas temperature model and a solid temperature model. Next, the modeling region and the two-dimensional mathematical model are introduced. The variables and the modeling region are as Figure 3 shown, from the dead stock column to the cross temperature measurement, and this part is modeled and solved.

[0260] In this embodiment, the following variables are set for the two-dimensional mathematical model of the whole blast furnace:

[0261]

[0262]

[0263] The magnitude of the body force acting on the fluid volume element in the gas flow model is proportional to the volume of the fluid element, including gravity, inertial force, electromagnetic force, etc., and acts at the centroid of the control volume. Use f x 、f y 、f z to represent the body force per unit mass of the fluid. After considering the body force and converting the three-dimensional to two-dimensional, the momentum conservation equation (18) is reduced to:

[0264]

[0265] The source term includes different parts in various equations. For example, all terms such as chemical reaction rate, interaction between phases, heat exchange between phases and phase change that cannot be classified into other terms need to be considered in the source term.

[0266] The descending speed of the blast furnace burden is about 3 - 4 m / h, and the ascending speed of the gas is about 2.5 - 6.8 m / s. For the convenience of solving and the descending speed of the burden can be ignored relative to the ascending speed of the gas. Therefore, the momentum transfer equation (26) from the gas phase to the solid phase is reduced to:

[0267]

[0268] Considering the momentum transfer equation (32) between the gas and the solid in the source term of the gas momentum conservation control equation, equation (31) is reduced to:

[0269]

[0270]

[0271] In the formula:

[0272]

[0273] μ is the dynamic viscosity coefficient.

[0274]

[0275] Meanwhile, the gas flow model also includes the continuity equation. Reducing the continuity equation (19) to two dimensions gives:

[0276]

[0277] The gas temperature model is established based on the energy conservation equation (20). Meanwhile, the heat transfer between phases and the heat of chemical reaction are considered in the source term; the heat transfer equation between gas and solid (28) is considered in the source term of the energy conservation control equation, and equation (20) is transformed into:

[0278]

[0279] Where:

[0280]

[0281] k is the thermal conductivity of Fourier's law, c v is the specific heat:

[0282]

[0283] k g is the thermal conductivity of the gas, d s is the average particle size, Re is the Reynolds number, Pr g is the convective heat transfer Prandtl coefficient of the gas, representing a dimensionless combination number that reflects the mutual influence of the energy and momentum transfer processes in the fluid, and the relationship between the surface temperature boundary layer and the flow boundary layer.

[0284] The solid temperature model has nothing to do with the basic equation in step 3 and is analyzed separately. The heat transfer between solids is described using the heat conduction equation (derived form of Fourier's law):

[0285]

[0286] Where: k is the thermal conductivity, which is determined by the thermal conductivity, density, and heat capacity of the material;

[0287] In the blast furnace, in addition to the heat transfer between solids, the influence of the gas on the solids and the heat of the reaction of the burden itself on the solids also exist. Therefore, the solid temperature model is as follows.

[0288]

[0289] Where:

[0290]

[0291] k g is the thermal conductivity of the gas, d s is the average particle size, Re is the Reynolds number, Pr g is the convective heat transfer Prandtl coefficient of the gas, representing the dimensionless combination number that reflects the mutual influence of the energy and momentum transfer processes in the fluid, and the relationship between the surface temperature boundary layer and the flow boundary layer. The reaction heat of the burden itself is obtained by substituting into the empirical formula:

[0292]

[0293] In the formula:

[0294]

[0295] D p is the diameter of the ore balls in the furnace, and Q1, Q2, E, R, K1 are all constants, which are related to the compositions of the ore and coke;

[0296] Step 6: Solve the 2D mathematical model of the blast furnace using the finite difference method. The finite difference method is solved on the grid, so it is necessary to first divide the grid of the modeling area and then discretely solve the model on the grid;

[0297] Use the BFC grid generation technology to solve and obtain the body-fitted coordinates. Through coordinate transformation, the irregular area in the physical space is mapped one-to-one with the regular area in the computational space, and the discretization of the partial differential equations in the model is transformed from the original irregular area to the regular area in the computational space. The present invention uses the elliptic partial differential equation in the differential equation transformation method to generate the BFC grid.

[0298] Step 6.1: Generate the grid for the 2D blast furnace modeling area;

[0299] As Figure 4 shown, the grid generation process is divided into the following 4 steps

[0300] Step 6.1.1: First, set the number of grids and the boundary conditions of the physical area, and use the linear interpolation method to divide the grid in the physical area. The division result is used as the initial value for the calculation iteration;

[0301] Step 6.1.2: Discretize the elliptic differential equations (45) to (47), and obtain the values of J, α, β, γ according to the coordinate distribution at the current moment;

[0302]

[0303] In the formula, J, α, β, γ are coefficients, specifically

[0304]

[0305] P and Q are source terms, and the calculation form established by Wei Wenli et al. is used, specifically as follows:

[0306]

[0307] Step 6.1.3: Calculate the source terms: First, use Equation (48) to calculate the values of φ(i,j) and ψ(i,j) on the boundary conditions:

[0308]

[0309] Then, use linear interpolation to obtain φ(i,j) and ψ(i,j) at all points. Finally, obtain P(i,j) and Q(i,j) according to Equation (47);

[0310] Step 6.1.4: After obtaining all parameters, calculate the grid coordinate distribution in the physical plane according to Equation (45), and determine whether the set accuracy requirements are met; if so, generate a grid diagram, otherwise, return to Step 6.1.2.

[0311] Step 6.2: Solve the two-dimensional mathematical model of the blast furnace; select the finite difference method to discretize the model, select the staggered grid method, and perform iterative solution in the time direction on the grid generated by the BFC grid generation technology in Step 6.1;

[0312] Step 6.2.1: Determine that the tuyere of the blast furnace is the gas inlet and the top is the gas outlet. According to the parameters (velocity, pressure, temperature) of the gas inlet and gas outlet collected on site, these are the boundary conditions at these two locations; the furnace wall and the bottom are set with no-slip boundary conditions, and the furnace wall temperature is obtained in Step 2; the symmetry axis adopts the boundary condition with the gradient of each parameter being 0, and the internal parameters of the steady state of the blast furnace are used as the initial conditions;

[0313] Step 6.2.2: Calculate the velocity at the next moment through the gas flow model and the gas velocity and pressure distributions at the current moment and pressure

[0314] Step 6.2.3: Calculate the gas temperature at the next moment through the gas temperature model and the gas velocity, gas temperature, and solid temperature distributions at the current moment

[0315] Step 6.2.4: Calculate the solid temperature at the next moment through the solid temperature model and the gas temperature and solid temperature distributions at the current moment

[0316] Step 6.2.5: Determine whether the change in the temperature value at the current moment and the temperature value at the next moment meets the accuracy requirements. If it meets, go to Step 7; if not, return to Step 6.2.2.

[0317] Step 7: POD order reduction optimization; Since the number of iterations in Step 6.2 is too large and the solution speed is slow, the present invention uses the POD order reduction method to optimize Step 6.2. All the time periods described below are the time periods during the iteration process in Step 6.2.

[0318] Step 7.1: First, use Step 6.2 to solve for the L time period and then stop, and give the approximate solutions of the first L time periods to form a snapshot matrix and These five snapshot matrices respectively represent the radial velocity, circumferential velocity, pressure, gas temperature, and solid temperature in the two-dimensional modeling area of the blast furnace in the first L time periods. Where m = MN is the total number of grid points in the modeling area, 1 ≤ i ≤ m, 1 ≤ l ≤ L.

[0319] Step 7.2: Solve the following linear equations according to the above five snapshot matrices:

[0320]

[0321] Obtain five groups of eigenvalues of the above linear equations

[0322]

[0323] and the corresponding five groups of eigenvectors

[0324]

[0325] Step 7.3: On the premise of meeting the POD basis error estimation condition , determine the acceptable error e = O(△t, △x 2 , △y 2 ) of the two-dimensional mathematical model of the blast furnace in the POD optimization stage. The order of the POD basis is the model order reduction order and constitute the initial POD basis:

[0326]

[0327] Step 7.4: Solve the POD eigenvalues; Assume that the solution of the entire blast furnace shaft temperature model is and Then, they converge to

[0328]

[0329] Obtained according to the two-dimensional mathematical model of the blast furnace

[0330]

[0331] wherein Calculated according to the two-dimensional mathematical model of the blast furnace;

[0332] Let

[0333] Based on the data of the previous L time periods and the POD basis, calculate the POD eigenvalue of each of the previous L segments

[0334]

[0335] Calculate the POD eigenvalue of each segment from the L-th time period to the (T - 1)-th time period

[0336]

[0337] wherein

[0338]

[0339] Use equations (56) and (57) to calculate the POD eigenvalue of each segment from the L-th time period to the convergence end time period.

[0340] Step 7.5: Solve the model distribution; based on the POD basis and POD eigenvalues obtained in the above steps 7.3 and 7.4, obtain from equation (58)

[0341]

[0342] Use to replace to obtain the target value distribution of the internal temperature model of the blast furnace and In the solution of this step, due to errors, the POD basis and eigenvalues need to be updated. Therefore, after solving the parameter distribution for each time step, first calculate the differences in pressure, velocity, and temperature from the previous moment to determine whether the iteration is completed. If the termination accuracy is met, the iteration is completed and the program ends. If not, first, according to what is described in step 7.6, determine whether the result of obtaining the POD basis at the current moment meets the POD error requirement. If not, perform step 7.6 to update the POD basis. If it meets the requirement, calculate the distribution at the next moment based on the POD eigenvalue and POD basis at the next moment until the termination accuracy is met.

[0343] Step 7.6: Use to replace It is the main error source of the POD method. This step analyzes the error and determines the update of the POD basis.

[0344] Let the POD error be

[0345]

[0346] According to Equation (52), we get (60)

[0347]

[0348] According to Equation (52), we have

[0349]

[0350] Let On the premise of ensuring the stability of the difference equation, according to Equation (61), we get

[0351]

[0352] Under the conditional condition of ensuring the stability of the differential equation, there is the following error estimation between the accurate solution of the blast furnace temperature distribution model and the approximate solution of the POD reduced-order algorithm:

[0353]

[0354] Where

[0355]

[0356]

[0357] Therefore, in step 7.5, let Judge Whether it holds. If it holds, then And Are the solutions that meet the accuracy requirements for the current time period. Otherwise, repeat steps 7.1 to 7.4 until the And That meet the above accuracy are obtained, and then continue with step 7.5;

[0358] Step 8: Use the two-dimensional distribution results (velocity, temperature) obtained in steps 6 and 7 as the calculation basis for the reaction rate in all chemical reactions in the blast furnace. Calculate the enthalpy change during the chemical reaction according to the chemical reaction rate and the composition distribution generated by the flow, and perform temperature field correction and coupling modeling based on the enthalpy change, as Figure 5 Shown;

[0359] All chemical reactions in the blast furnace are listed in this embodiment, specifically including:

[0360] 1. Combustion reaction: The reaction equation is as follows:

[0361] 2C(s) + O2(g) = 2CO(g) (67)

[0362] Due to the very high gas velocity in the raceway of the tuyere, the reaction is intense, which is very different from the gas velocity and reaction efficiency in the shaft furnace, and is not conducive to the simulation of the model. Therefore, in the present invention, the coke combustion reaction is implicitly processed, and the combustion reaction is directly processed at the tuyere, and the substances and energy generated by the combustion reaction are directly added to the blast tuyere.

[0363] 2. Direct reduction reaction: The reaction equation is as follows:

[0364] FeO(l) + C(s) = Fe(l) + CO(g) (68)

[0365] The present invention believes that the direct reduction reaction occurs on the surface of the coke, reducing the molten FeO to Fe. The reaction rate is proportional to the particle size of the coke and the activity of the molten FeO. The reaction rate of the direct reduction reaction is:

[0366] R2 = k2(A c / V B )a FeO (69)

[0367] In the formula, R2 is the calculated reaction rate of the direct reduction reaction; k2 is the reaction rate constant of the direct reduction reaction, obtained by fitting experimental data:

[0368]

[0369] where T s is the solid temperature, A c / V B is the effective specific surface area of the direct reduction reaction, approximated from experimental data:

[0370]

[0371] where ε s is the porosity of the burden; ψ s is the shape factor of the solid particles; d s is the diameter of the solid particles, a FeO is the activity factor of FeO, calculated according to the following formula:

[0372] a FeO =(N FeO ) 0.55 (72)

[0373] where N FeO is the mole fraction of FeO;

[0374] 3. Indirect reduction reaction: The reaction equation is as follows:

[0375]

[0376] The unreacted shrinking core model is used to simulate the indirect reduction reaction. The overall chemical reaction rate of the indirect reduction reaction is

[0377]

[0378] where A s is the diffusion area; P g is the gas pressure; R is the ideal gas constant; T s is the solid temperature; is the difference in the CO mole fractions inside and outside the iron ore; k3′ is the gas film diffusion coefficient of CO on the surface of the iron ore and is solved according to the following formula:

[0379]

[0380] where Sh CO is the Sherwood number of CO; Sc CO is the Schmidt number of CO; D CO is the diffusion coefficient of CO; η g is the gas viscosity, ρ g is the gas density; Re g-s is the modified Reynolds number; d s is the solid particle diameter; P g is the gas pressure; T g is the gas temperature, f s is the reduction degree of the iron ore:

[0381]

[0382] where r0 is the radius of the iron ore particle; r1 is the radius of the reaction spherical surface, k3″ is the diffusion coefficient of CO inside the iron ore and is solved according to the following formula:

[0383]

[0384] where ε Fe is the porosity of Fe in the iron ore particle; ζ io is the tortuosity factor of the iron ore particle; ε io is the porosity of the iron ore particle, k3″′ is the chemical reaction rate constant of the indirect reduction reaction:

[0385]

[0386] K3 is the equilibrium constant of the indirect reduction reaction. It is calculated according to the following formula based on the different solid phase temperatures and reduction degrees:

[0387]

[0388] 4. Carbon loss reaction: The chemical reaction equation is:

[0389] C(s) + CO2(g) = 2CO(g) (80)

[0390] The unreacted shrinking core model is used to simulate the indirect reduction reaction. The overall chemical reaction rate of the carbon loss reaction is

[0391]

[0392] where k5′ is the mass transfer coefficient of CO2 in the coke gas film diffusion stage. The calculation method of k5′ is similar to that of k3′ and is shown in the following formula.

[0393]

[0394] E5 is the effective factor of coke loss:

[0395]

[0396] where m is the Thiele modulus; k5″′ is the chemical reaction rate constant of carbon loss:

[0397]

[0398] ε dreg is the porosity of the reaction slag on the coke particle surface; ζ C is the tortuosity factor of the coke particle.

[0399]

[0400] where ε C is the porosity of the coke.

[0401] Step 9: Temperature field optimization and component concentration. Calculate the chemical reaction rate of each reaction according to the temperature field, and then calculate the enthalpy change △H and the change in CO gas flow concentration according to the chemical reaction rate. Use the enthalpy change to correct the temperature distribution and component distribution, and finally realize the digital twin of the temperature field and component field during the operation of the blast furnace. As Figure 4 shown.

[0402] Step 9.1: Assume the internal component distribution of the blast furnace without considering chemical reactions based on the velocity field in the blast furnace modeling area obtained in Steps 6 and 7.

[0403] Step 9.2: Calculate the chemical reaction rate of each chemical reaction in Step 8 according to the temperature field in the blast furnace modeling area obtained in Steps 6 and 7.

[0404] Step 9.3: Calculate the enthalpy change and composition change of each chemical reaction according to the chemical rate. The chemical reactions are as follows Figure 6 as shown.

[0405] Step 9.4: Update the temperature field and composition field according to the enthalpy change and composition change.

[0406] Step 9.5: Determine whether the coupling degree error between the temperature field and the composition field meets the accuracy. If it meets, end; otherwise, continue with Step 9.1.

[0407] The above description is only the preferred embodiments of the present disclosure and the explanation of the applied technical principles. Those skilled in the art should understand that the scope of the invention involved in the embodiments of the present disclosure is not limited to the technical solutions formed by the specific combination of the above technical features, but should also cover other technical solutions formed by any combination of the above technical features or their equivalent features without departing from the above inventive concept. For example, the technical solutions formed by mutually replacing the above features with the technical features (but not limited to) disclosed in the embodiments of the present disclosure that have similar functions.

Claims

1. A digital twin method for the production operation state of a blast furnace, characterized in that, It includes the following steps: Step 1: Calculate the temperature field of the blast furnace top based on the infrared image, perform gray correction on the infrared image, and use the lower surface of the temperature field as the boundary condition at the model furnace top; Step 2: Use the mechanism model to calculate the inner surface temperature of the blast furnace as the boundary condition at the model furnace wall; Step 3: Take the blast furnace shaft mass points as the research object, analyze and establish the basic model of the blast furnace shaft temperature field, use the N - S equation described by Euler to solve the physical quantities of each fluid mass point at different positions and different times in the entire flow field domain, including pressure, velocity, and temperature, and describe the physical state of the entire flow field domain through these physical quantities; Divide the inside of the blast furnace into countless fluid micro - groups, analyze a single fluid micro - group, and obtain the basic control equation set of the blast furnace model; Step 4: Analyze the transfer phenomena between heterogeneous and homogeneous phases existing in the blast furnace as the source term; Specifically include gas - solid momentum transfer and gas - solid energy transfer; Step 5: Since the blast furnace is a symmetric model, a two - dimensional mathematical model of the blast furnace is established; Step 6: Solve the two - dimensional mathematical model of the blast furnace using the finite - difference method. The finite - difference method is solved on the grid, so it is necessary to first divide the grid of the modeling area, and then discretely solve the model on the grid; Use the BFC grid generation technology to solve and obtain the body - fitted coordinates. Through coordinate transformation, the irregular area in the physical space is put into one - to - one correspondence with the regular area in the computational space, and the discretization of the partial differential equation in the model is transformed from the original irregular area to the regular area in the computational space; Step 8: POD reduction optimization; Step 9: Use the results obtained in Step 6 and Step 7 as the calculation basis for the reaction rate in all chemical reactions in the blast furnace. Calculate the enthalpy change during the chemical reaction according to the chemical reaction rate and the composition distribution generated by the flow, and perform temperature field correction and coupled modeling according to the enthalpy change; Step 10: Temperature field optimization and composition concentration; Calculate the chemical reaction rate of each reaction according to the temperature field, then calculate the enthalpy change △H and the change in the concentration of the CO gas flow according to the chemical reaction rate, and use the enthalpy change to correct the temperature distribution and composition distribution, and finally realize the digital twin of the temperature field and composition field during the operation of the blast furnace.

2. The digital twin method for the production operation state of a blast furnace according to claim 1, wherein The gray-scale correction described in Step 1 is as follows: Determine several gray-scale levels (G0, G1, …, G n ) in the infrared image, where G0 is the lowest gray-scale value and G n is the highest gray-scale value. Use an image segmentation algorithm to extract the completely black and completely white regions of the image. Draw isotherms in the completely white and completely black regions according to the temperatures calculated by cross thermometry. The number of isotherms m b in the completely black region and the number of isotherms m w in the completely white region are both related to the gray-scale level difference, as shown below: In the formula, k is the gray - level difference of the isotherm; Adjust the image grayscale. For the all-white area, first set the grayscale of the highest-temperature area to G n , and set the grayscale of the lowest-temperature area to G n-1 . The grayscale of the intermediate-temperature area is adjusted gradually according to the temperature from high to low; for the all-black area, first set the grayscale of the lowest-temperature area to G0, the grayscale of the highest-temperature area to G1, and the grayscale of the intermediate-temperature area is adjusted gradually according to the temperature from low to high; according to the position symmetry of the cross thermometry device and the order of magnitude of the reference temperature, the grayscale is divided into corresponding levels; count the grayscale values near the thermocouple in the image to determine the reference grayscale value corresponding to the thermocouple reference temperature; the temperature values between the reference grayscales are realized by linear interpolation between the corresponding reference temperatures using the two-point calibration method; take the average of the results of multiple calibrations to improve the accuracy of calibration; The two - point calibration method is described as follows. Assume the characteristics of temperature calibration: y = ax + b (2) Among them, x is the grayscale input signal, y is the temperature value output signal. By obtaining the gain factor a and the offset factor b, the corresponding relationship between grayscale and temperature can be obtained. Let G min be the lowest grayscale value, and the corresponding temperature at this time is T min , G max be the highest grayscale value, and its corresponding temperature is T max . Substituting into the equation gives:

3. A digital twin method for the production operation state of a blast furnace according to claim 1, characterized in that, The specific content of Step 2 is: First, establish the heat conduction differential equation of the blast furnace wall, determine the boundary conditions, establish a heat transfer model for the blast furnace wall and solve it to obtain the temperature distribution inside the furnace wall; First, establish the heat conduction differential equation of the blast furnace wall; For the differential control volume, according to the law of conservation of energy: Input - Output + R = S (4) where Input is the sum of the heat input into the control volume per unit time, Input = Input x + Input y + Input z , where Input x = q x dydz, Input y = q y dxdz, Input z = q z dxdy are the heat input into the control volume from the x, y, and z directions per unit time respectively; Output is the sum of the heat output from the control volume per unit time, Output = Output x + Output y + Output z , where are the heat output from the control volume from the x, y, and z directions per unit time respectively; R is the heat generated inside the control volume, i.e., the internal heat source; S is the change in internal energy of the control volume over time after absorbing or releasing heat; Therefore: There is Fourier's law: where q is the heat conduction flux, is the temperature gradient, and λ is the thermal conductivity; Since the furnace wall is solid and the heat transfer mode inside it is only heat conduction, substitute Fourier's law (6) into the heat conduction differential equation (5) to get: Therefore, transform it into the heat conduction differential equation in the cylindrical coordinate system. Obtain it from the heat conduction differential equation in the rectangular coordinate system through the mathematical coordinate transformation x = rcosθ, y = rsinθ, z = h In the formula, r is the cross - sectional radius, θ is the cross - sectional angle, h is the longitudinal axis height, and T is the temperature distribution of the furnace wall; In order to obtain the circumferential temperature distribution at different heights on the inner surface of the blast furnace, first, a heat transfer mechanism model of the cross-section of the blast furnace wall is established. Without considering the axial distribution, we have: The above heat conduction differential equation is discretized by central difference in the radial and circumferential directions, and the discretized difference equation is as follows: where \(r_0\) is the radius of the inner ring of the circular ring model, \(r_2\) is the radius of the outer ring of the circular ring model, \(\Delta r\) is the radial step size, \(\Delta\theta\) is the circumferential step size, \(r\) i is the radius at the \(i\)th position, \(M\) is the number of radial differential steps, and \(N\) is the number of circumferential differential steps; Inner surface boundary of the furnace wall: T 1,j = T in,j , r = r0, 0 ≤ j ≤ N (11) Boundary at the junction of refractory bricks and cooling stave: After discretization: Outer surface boundary of the cooling stave: After discretization: Where: T in , T out are the temperatures of the inner surface of the blast furnace and the outer surface of the cooling stave; T CW is the average temperature of the cooling water around the cooling stave, and T CW =(T CW,in +T CW,out ) / 2, T CW,in is the inlet temperature of the cooling water, and T CW,out is the outlet temperature of the cooling water; r0, r1, and r2 are the radii of the inner surface of the blast furnace refractory brick, the junction between the refractory brick and the cooling stave, and the cooling pipe; λ1 and λ2 are the thermal conductivities of the refractory brick and the cooling stave, respectively; h x is the convective heat transfer coefficient between the cooling stave and the cooling water, and its value is calculated according to the following empirical formula, where v 水 is the flow velocity of the cooling water: h x = 208 + 47.5v 水 (16) Taking the known temperature of the thermocouple on the furnace wall as the comparison value, using the above equations (9) and (10) and the boundary conditions (11) to (15) as the calculation model, with the cooling water pipe as the outer boundary and the inner surface as the inner boundary, the temperature field of the furnace wall is solved; assuming the inner surface temperature to solve the thermocouple temperature, calculating the error between the true value and the calculated value of the thermocouple temperature. When the error is satisfied, the assumed inner surface temperature is the true inner surface temperature.

4. A digital twin method for the production operation state of a blast furnace according to claim 1, characterized in that Step 3 is specifically as follows: Based on the continuum hypothesis, according to the three conservation laws of mass, momentum, and energy, taking the control volume described by Euler as the research object, without considering the body force and source term, the control volume is analyzed; the increase in the momentum of the control volume is equal to the momentum flowing in through the control surface. For the control volume, we have: where ρ is the gas density; Δx, Δy, and Δz are the sizes of the differential control volume in three directions respectively; Δt is the time step; are the momentum changes in three directions; The momentum conservation equation is transformed into: In the formula: μ is the dynamic viscosity coefficient; Similarly, the mass conservation equation is: Similarly, the energy conservation equation is: Wherein: k is the thermal conductivity of Fourier's law; c v is the specific heat.

5. A digital twin method for the production operation state of a blast furnace according to claim 1, characterized in that, The gas-solid momentum transfer in Step 4: Analyze the relationship between the multiphase flows in the blast furnace, and uniformly describe the momentum transfer amount per unit time between each phase state in the form of the following equation: F i-j = f i-j (U i - U j )(21) where F i-j is the amount of momentum transfer per unit time between the i-th and j-th phases; f i-j is the momentum transfer coefficient between the i-phase and the j-phase; U i and U j are the velocities of the i-phase and the j-phase respectively; Assuming the inclination angle of each layer of solid phase is θ, the momentum transfer phenomena between the solid phase and the gas phase are described as follows: Where: R perp is the flow resistance coefficient in the direction perpendicular to the material surface and is expressed by the following formula: R para is the flow resistance coefficient in the direction parallel to the material surface and is expressed by the following formula: R k is the resistance coefficient of different components of the material and is expressed using the Ergun equation: f k is the volume fraction of different components of the material; Substituting equations (23)-(25) into (22) and assuming the inclination angle is 0, we get: where F g-s is the amount of momentum transfer per unit time between the gas and the solid; η g is the gas viscosity; u g is the gas velocity vector; u s is the solid velocity vector; ρ g is the gas density; ε s is the porosity of the burden; ψ s is the shape factor of the solid particles; d s is the solid particle diameter; The gas-solid energy transfer: The energy transfer amount between phase states is uniformly described in the form of the following equation: F i-j = A i-j h i-j (T i - T j ) (27) where F i-j is the energy transfer amount between the i-th phase and the j-th phase; A i-j is the contact area; h i-j is the energy transfer coefficient between the i-th phase and the j-th phase; T i and T j are the temperatures of the i-th phase and the j-th phase, respectively; Calculating the energy transfer coefficient according to the modified Ranz-Marshall formula, the energy transfer phenomena between the gas phase and the solid phase are described as follows: where k g is the thermal conductivity of the gas; d s is the average particle size; Re is the Reynolds number, which is calculated according to the following formula: where η g is the gas viscosity; u g is the gas velocity vector; u s is the solid velocity vector; ρ g is the gas density; ψ s is the shape factor of solid particles; d s is the solid particle diameter, Prandtl number Pr of convective heat transfer of gas g A dimensionless combination number representing the mutual influence of the energy and momentum transfer processes in a fluid. The relationship between the surface temperature boundary layer and the flow boundary layer is calculated according to the following formula, where C g is the specific heat capacity of the gas:

6. The digital twin method for the production operation state of a blast furnace according to claim 1, characterized in that Step 5 is specifically as follows: By introducing the source term in Step 4 into the basic control equation set analyzed in Step 3, a two-dimensional mathematical model of the blast furnace is established, establishing the specific continuity equation, momentum conservation equation, and energy conservation equation of the blast furnace, and specifying the model to obtain the two-dimensional mathematical model of the blast furnace process. The two-dimensional mathematical model of the blast furnace includes a gas flow model, a gas temperature model, and a solid temperature model; The magnitude of the body force exerted by the gas flow model on the fluid volume element is proportional to the volume of the fluid element, and f x , f y , f z are used to represent the body force per unit mass of the fluid. After considering the body force and converting from three dimensions to two dimensions, the momentum conservation equation (18) is transformed into: The momentum transfer equation (26) from the gas phase to the solid phase is transformed into: Considering the momentum transfer equation (32) between the gas and the solid in the source term of the gas momentum conservation control equation, equation (31) is transformed into: In the formula: μ is the dynamic viscosity coefficient; At the same time, the gas flow model also includes a continuity equation. Reducing the continuity equation (19) to two dimensions, we have: The gas temperature model is established according to the energy conservation equation (20), and at the same time, the heat transfer between phases and the chemical reaction heat are considered in the source term; Considering the heat transfer equation (28) between the gas and the solid in the source term of the energy conservation control equation, equation (20) is transformed into: In the formula: k is the thermal conductivity of Fourier's law, c v is the specific heat capacity: k g is the thermal conductivity of the gas, d s is the average particle size, Re is the Reynolds number, Pr g is the convective heat transfer Prandtl coefficient of the gas, which represents the dimensionless combination number that reflects the mutual influence of the energy and momentum transfer processes in the fluid, and the relationship between the surface temperature boundary layer and the flow boundary layer; The heat transfer between solids in the solid temperature model is described by the heat conduction equation: In the formula: k is the thermal conductivity, which is determined by the thermal conductivity, density, and heat capacity of the material; Therefore, the solid temperature model is as follows: In the formula: k g is the thermal conductivity of the gas, d s is the average particle size, Re is the Reynolds number, Pr g is the convective heat transfer Prandtl coefficient of the gas, which represents the dimensionless combination number of the mutual influence of the energy and momentum transfer processes in the fluid, and the relationship between the surface temperature boundary layer and the flow boundary layer; The reaction heat of the burden itself is obtained by substituting the empirical formula: In the formula: D p is the diameter of the ore balls in the furnace, and Q1, Q2, E, R, and K1 are all constants related to the compositions of the ore and coke.

7. A digital twin method for the production operation status of a blast furnace according to claim 1, characterized in that The specific steps of step 6 are as follows: Step 6.1: Generation of the grid for the two-dimensional blast furnace modeling area; Step 6.1.1: First, set the number of grids and the boundary conditions of the physical area, and use the linear interpolation method to divide the grids in the physical area. The division result is used as the initial value for the calculation iteration; Step 6.1.2: Discretize the elliptic differential equations (45) to (47), and obtain the values of the coefficients J, α, β, and γ according to the coordinate distribution at the current moment; P and Q are source terms, specifically: Step 6.1.3: Calculate the source terms: First, use formula (48) to calculate the values of φ(i,j) and ψ(i,j) on the boundary conditions: Then use linear interpolation to obtain φ(i,j) and ψ(i,j) at all points, and finally obtain P(i,j) and Q(i,j) according to formula (47); After obtaining all the parameters, calculate the grid coordinate distribution in the physical plane according to formula (45), and determine whether it meets the set accuracy requirements; if it meets, generate a grid diagram, otherwise, return to step 6.1.2; Step 6.2: Solve the two-dimensional mathematical model of the blast furnace; Select the finite difference method to discretize the model, select the staggered grid method, and perform iterative solution in the time direction on the grid generated by the BFC grid generation technology in step 6.1; Step 6.2.1: Determine that the tuyere of the blast furnace is the air inlet and the top is the air outlet. According to the parameters of the air inlet and the air outlet collected on-site, they are the boundary conditions at these two places; the furnace wall and the bottom are set with no-slip boundary conditions, and the furnace wall temperature is obtained in step 2; the symmetry axis adopts the boundary condition with a gradient of 0 for each parameter, and the internal parameters in the steady state of the blast furnace are used as the initial conditions; Step 6.2.2: Calculate the velocity and pressure at the next moment based on the gas flow model and the gas velocity and pressure distribution at the current moment and pressure Step 6.2.3: Calculate the gas temperature at the next moment through the gas temperature model and the gas velocity, gas temperature, and solid temperature distributions at the current moment Step 6.2.4: Calculate the solid temperature at the next moment based on the solid temperature model, the gas temperature at the current moment, and the solid temperature distribution Step 6.2.5: Judge whether the change in the temperature value at the current moment and the temperature value at the next moment meets the accuracy requirements. If it meets, go to step 7; if it does not meet, return to step 6.2.

2.

8. A digital twin method for the production operation state of a blast furnace according to claim 1, characterized in that, The specific steps of step 7 are as follows: Step 7.1: First, use Step 6.2 to solve for the L time periods and then stop, giving the approximate solutions for the first L time periods to form a snapshot matrix and These five snapshot matrices respectively represent the radial velocity, circumferential velocity, pressure, gas temperature, and solid temperature in the blast furnace two-dimensional modeling region for the first L time periods; where m = MN is the total number of grid points in the modeling region, 1 ≤ i ≤ m, 1 ≤ l ≤ L; Step 7.2: Solve the following linear equations according to the above five snapshot matrices: Obtain five groups of eigenvalues of the above linear equations And the corresponding five groups of eigenvectors Step 7.3: On the premise of satisfying the POD basis error estimation condition , determine the acceptable error e = O(△t, △x 2 , △y 2 ) of the 2D mathematical model of the blast furnace in the POD optimization stage. The order of the POD basis is the model reduction order and constitute the initial POD basis: Step 7.4: Solve the POD eigenvalues; Assume that the solution of the entire blast furnace shaft temperature model is and Then, they converge to According to the two-dimensional mathematical model of the blast furnace Among them Obtained according to the 2D mathematical model of the blast furnace; Let According to the data of the first L time periods and the POD basis, calculate the POD eigenvalues of each of the first L segments Calculate the POD eigenvalues of each segment from the L time period to the T-1 time period Where Use formulas (56) and (57) to calculate the POD eigenvalues of each segment from the L time period to the convergence end time period; Step 7.5: Solve the model distribution; The POD basis and POD eigenvalues obtained according to the above steps 7.3 and 7.4 are obtained from Equation (58). Use instead of The target value distribution of the internal temperature model of the blast furnace is obtained and First, calculate the differences between the pressure, velocity, and temperature and those at the previous moment to determine whether the iteration is completed. If the termination accuracy is met, the iteration is completed and the program ends; if not, first, as described in step 7.6, determine whether the POD basis calculation result at the current moment meets the POD error requirement. If not, perform step 7.6 to update the POD basis; if it meets the requirement, calculate the distribution at the next moment based on the POD eigenvalues and POD basis at the next moment until the termination accuracy is met; Step 7.6: Use to replace which is the main error source of the POD method; analyze the error and determine the update of the POD basis; Let the POD error be According to formula (52) According to formula (52), there is Let On the premise of ensuring the stability of the difference equation, according to Equation (61), we have Under the conditional condition of ensuring the stability of the differential equation, there is the following error estimate between the accurate solution of the blast furnace temperature distribution model and the approximate solution of the POD reduction algorithm: Where Therefore, in step 7.5, for each time period, let Judge Whether it holds. If it holds, then And Are solutions that meet the accuracy requirements for the current time period. Otherwise, repeat steps 7.1 to 7.4 until the And That meet the above accuracy are obtained, and then continue with step 7.

5.

9. A digital twin method for the production operation state of a blast furnace according to claim 1, characterized in that, The specific steps of step 9 are as follows: Step 9.1: Assume the internal component distribution of the blast furnace without considering chemical reactions according to the velocity field of the blast furnace modeling area obtained in step 6 and step 7; Step 9.2: Calculate the chemical reaction rates of each chemical reaction in Step 8 based on the temperature field of the blast furnace modeling area obtained in Step 6 and Step 7; Step 9.3: Calculate the enthalpy change and composition change of each chemical reaction according to the chemical rate; Step 9.4: Update the temperature field and composition field according to the enthalpy change and composition change; Step 9.5: Determine whether the coupling degree error between the temperature field and the composition field meets the accuracy; if it meets, end, otherwise continue with Step 9.1.

Images