CN113627064B - Roller kiln firing zone temperature prediction method based on mechanism and data mixed driving - Google Patents
Roller kiln firing zone temperature prediction method based on mechanism and data mixed driving Download PDFInfo
- Publication number
- CN113627064B CN113627064B CN202111033301.0A CN202111033301A CN113627064B CN 113627064 B CN113627064 B CN 113627064B CN 202111033301 A CN202111033301 A CN 202111033301A CN 113627064 B CN113627064 B CN 113627064B
- Authority
- CN
- China
- Prior art keywords
- temperature
- kiln
- heat transfer
- heat
- flue gas
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 97
- 230000007246 mechanism Effects 0.000 title claims abstract description 75
- 238000010304 firing Methods 0.000 title claims abstract description 70
- 239000000919 ceramic Substances 0.000 claims abstract description 63
- 230000008569 process Effects 0.000 claims abstract description 33
- 238000012628 principal component regression Methods 0.000 claims abstract description 19
- 238000004134 energy conservation Methods 0.000 claims abstract description 14
- 238000012546 transfer Methods 0.000 claims description 98
- UGFAIRIUMAVXCW-UHFFFAOYSA-N Carbon monoxide Chemical compound [O+]#[C-] UGFAIRIUMAVXCW-UHFFFAOYSA-N 0.000 claims description 80
- 239000003546 flue gas Substances 0.000 claims description 80
- VNWKTOKETHGBQD-UHFFFAOYSA-N methane Chemical compound C VNWKTOKETHGBQD-UHFFFAOYSA-N 0.000 claims description 60
- 239000011159 matrix material Substances 0.000 claims description 42
- 238000002485 combustion reaction Methods 0.000 claims description 30
- 239000003345 natural gas Substances 0.000 claims description 30
- 238000005245 sintering Methods 0.000 claims description 24
- 230000008859 change Effects 0.000 claims description 23
- 238000004364 calculation method Methods 0.000 claims description 20
- 230000005855 radiation Effects 0.000 claims description 20
- 239000007789 gas Substances 0.000 claims description 17
- 230000006870 function Effects 0.000 claims description 14
- 239000000779 smoke Substances 0.000 claims description 9
- 230000000694 effects Effects 0.000 claims description 4
- 238000010438 heat treatment Methods 0.000 claims description 4
- 238000013507 mapping Methods 0.000 claims description 4
- 230000005457 Black-body radiation Effects 0.000 claims description 3
- 238000000354 decomposition reaction Methods 0.000 claims description 3
- 238000009792 diffusion process Methods 0.000 claims description 3
- 230000004907 flux Effects 0.000 claims description 3
- 238000010606 normalization Methods 0.000 claims description 3
- 238000012545 processing Methods 0.000 claims description 3
- 239000013598 vector Substances 0.000 claims description 3
- 239000000523 sample Substances 0.000 description 26
- 238000004422 calculation algorithm Methods 0.000 description 16
- 238000010586 diagram Methods 0.000 description 10
- 238000004519 manufacturing process Methods 0.000 description 9
- 238000012549 training Methods 0.000 description 6
- 238000013459 approach Methods 0.000 description 3
- 239000000446 fuel Substances 0.000 description 3
- 238000011160 research Methods 0.000 description 3
- 230000008901 benefit Effects 0.000 description 2
- 238000006243 chemical reaction Methods 0.000 description 2
- 239000012774 insulation material Substances 0.000 description 2
- 230000007774 longterm Effects 0.000 description 2
- 238000013178 mathematical model Methods 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 244000179970 Monarda didyma Species 0.000 description 1
- 235000010672 Monarda didyma Nutrition 0.000 description 1
- 230000002159 abnormal effect Effects 0.000 description 1
- 230000005856 abnormality Effects 0.000 description 1
- 230000003044 adaptive effect Effects 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 239000011449 brick Substances 0.000 description 1
- 239000004568 cement Substances 0.000 description 1
- 230000008878 coupling Effects 0.000 description 1
- 238000010168 coupling process Methods 0.000 description 1
- 238000005859 coupling reaction Methods 0.000 description 1
- 238000013135 deep learning Methods 0.000 description 1
- 239000006185 dispersion Substances 0.000 description 1
- 238000005265 energy consumption Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 239000002737 fuel gas Substances 0.000 description 1
- 230000004927 fusion Effects 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000002156 mixing Methods 0.000 description 1
- 238000003062 neural network model Methods 0.000 description 1
- 238000005192 partition Methods 0.000 description 1
- 238000000513 principal component analysis Methods 0.000 description 1
- 239000002994 raw material Substances 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 230000001105 regulatory effect Effects 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/08—Thermal analysis or thermal optimisation
Abstract
The invention discloses a method for predicting the firing zone temperature of a roller kiln based on mechanism and data mixed driving, which comprises the following steps: s1: modeling a mechanism model of the mechanism process of the roller kiln firing zone by using a mass conservation law and an energy conservation law and solving by using a finite difference method; s2: based on the instant moving window and the local weighted kernel principal component regression, establishing an error compensation model based on data driving; s3: combining the mechanism model obtained in the step S1 with the error compensation model obtained in the step S2 to obtain a roller kiln firing zone temperature mixed prediction model; s4: and predicting the temperature of the firing zone of the roller kiln by using a mixed prediction model of the temperature of the firing zone of the roller kiln. The invention establishes a mechanism model of the firing zone temperature of the ceramic roller kiln, simultaneously uses a data modeling method to carry out temperature compensation on the error of the mechanism model result, and combines a mixed modeling strategy of mechanism modeling and data modeling to ensure the validity of the model and the high precision of the output result when the working condition changes.
Description
Technical Field
The invention relates to the field of temperature prediction in a ceramic roller kiln firing zone, in particular to a roller kiln firing zone temperature prediction method based on mechanism and data mixed driving.
Background
The firing of the ceramic is a process that a green brick is fired into a finished product, a process of multi-physical field coupling, nonlinearity, multiple variables, pure hysteresis, time variation and multiple disturbance is arranged in a kiln, the whole firing process of the ceramic is researched, a plurality of factors are involved, and a complex roller kiln global system is difficult to accurately establish. Therefore, it is scientific to grasp the most important steps in the ceramic production process and study the key firing process. The temperature of the firing zone is effectively predicted, the running condition in the kiln can be accurately mastered, and effective guidance is provided for operators, so that the products can be efficiently produced, and a certain guiding significance is provided for energy conservation and consumption reduction of enterprises.
Mechanism modeling is a mathematical model built up of aggregate objects, internal operating characteristics, and energy flow mechanism processes. The method is mainly based on a mathematical model established by an object-oriented process based on the law of conservation of mass, the law of conservation of energy, the law of conservation of momentum, certain physical equations, chemical reaction processes and the like. Although scholars at home and abroad do a great deal of research and have achieved a great deal of results in the aspect of mechanism modeling of the kiln, the research is mainly focused on tunnel kiln, rotary cement kiln and the like, and the research on mechanism modeling inside the ceramic roller kiln is less at present.
For the temperature prediction of the roller kiln, some students solve the problem by using an intelligent modeling method based on data driving. The method mainly comprises a multi-model approximation linear or nonlinear model, a neural network model, a multi-model fusion soft measurement model and the like.
The single internal mechanism modeling method can describe the working mechanism and the heat transfer mechanism in the kiln, but certain simplification and assumption of ideal state are carried out in the calculation process, so that the internal temperature change cannot be accurately described; the data-driven method can effectively predict the internal temperature, but cannot describe the heat transfer process and working mechanism of each part. So, a learner proposed to combine the internal mechanism modeling with the data-driven approach, and the combined approach is called the hybrid modeling approach.
Compared with a single mechanism model or a data driving model, the hybrid model has the advantages that the single mechanism model or the data driving model is adopted in the temperature prediction application of the ceramic roller kiln, the high precision and the robustness cannot be considered, and fewer cases are adopted.
The Chinese patent with publication number CN113190974A discloses a multipoint prediction method for a roller kiln temperature field for deep learning, which comprises the following steps: collecting each input parameter in the roller kiln by using a roller kiln sensor and the real sintering temperature in the roller kiln; simulating the temperature of the roller kiln by using simulation software according to the acquired input parameters, and outputting the simulated sintering temperature of the position points which cannot be measured by the sensor; defining an input data set according to input parameters and simulated sintering temperature, wherein the input data set comprises an input parameter matrix and an input temperature matrix, establishing a roller kiln temperature field multipoint prediction model based on a transducer, encoding the input parameter matrix by an encoder, converting an output matrix of the last encoder into a query matrix key matrix and taking the query matrix key matrix as input of each decoder, decoding the input temperature matrix by a plurality of decoders in a decoding stage, and finally outputting the multipoint temperature of a roller kiln sintering section at the current time point. The technical solution of this patent uses a single data-driven model, and cannot describe the heat transfer process and the working mechanism of each part.
Disclosure of Invention
The invention provides a method for predicting the firing zone temperature of a roller kiln based on mechanism and data mixed driving, which predicts the temperature in a production zone of a ceramic roller kiln and solves the problems that a single data driving modeling cannot describe the mechanism process in the kiln, the effectiveness of a model is prominent when the operation condition of the kiln changes and the robustness is not high.
In order to solve the technical problems, the technical scheme of the invention is as follows:
a method for predicting the firing zone temperature of a roller kiln based on mechanism and data mixed driving comprises the following steps:
s1: modeling a mechanism model of the mechanism process of the roller kiln firing zone by using a mass conservation law and an energy conservation law and solving by using a finite difference method;
s2: based on the instant moving window and the local weighted kernel principal component regression, establishing an error compensation model based on data driving;
s3: combining the mechanism model obtained in the step S1 with the error compensation model obtained in the step S2 to obtain a roller kiln firing zone temperature mixed prediction model;
s4: and predicting the temperature of the firing zone of the roller kiln by using a mixed prediction model of the temperature of the firing zone of the roller kiln.
Preferably, the step S1 specifically includes the following steps:
S1.1: analyzing the heat transfer process in the firing zone of the roller kiln, and establishing a roller kiln heat transfer process model;
s1.2: estimating state variables of the roller kiln firing belt, and setting assumption conditions of the state variables;
s1.3: based on mass conservation and energy conservation, establishing a flue gas temperature mechanism model of the roller kiln sintering belt;
s1.4: and carrying out numerical calculation and solving on the flue gas temperature mechanism model by adopting a finite difference method.
Preferably, step S1.1 is specifically:
the roller kiln sintering in-band heat transfer mode comprises convection heat transfer, radiation heat transfer and heat transfer, wherein:
the convective heat transfer comprises the convective heat transfer quantity between the flue gas and the product and the convective heat transfer quantity between the flue gas and the kiln wall, and the formula of the convective heat transfer is as follows:
Q convection current ＝h(t-t m )F
In which Q Convection current The total heat of convection heat transfer is expressed as W; h is the convection heat transfer coefficient, and the unit is W/(m) 2 K); t is the temperature of the surface of the product, and the unit is K; t is t m The unit is K, which is the temperature of the flue gas; f is the heating surface area of the product, and the unit is m 2 ；
The radiant heat transfer includes the radiant heat transfer between the article and the kiln wall and the radiant heat transfer between the transfer roller and the kiln wall, the radiant heat transfer having the formula:
in which Q Radiation of Heat radiated to the surface of the product by the gas is expressed as W; epsilon c Surface blackness of the product; f is the heated area of the product, and the unit is m 2 The method comprises the steps of carrying out a first treatment on the surface of the Sigma is the blackness of the gas; t is the absolute temperature of the gas, and the unit is K; t (T) p The absolute temperature of the surface of the product is expressed as K;
the heat conduction comprises heat transfer from a high-temperature region to a low-temperature region, heat absorbed by the product and the kiln wall and heat diffusion of the kiln wall, and the formula of the heat conduction is as follows:
in which Q Conduction of Is the heat flux density; k is the thermal conductivity; t is the temperature; x is the length in the heat transfer direction.
Preferably, step S1.2 is specifically:
the assumption conditions for setting the state variables include:
the system is in a stable running state;
the change of the temperature and the gas flow in the vertical direction is not considered, and only the change of the temperature and the gas flow in the length direction of the kiln is considered;
heat is transferred between the flue gas, the ceramic product and the kiln wall by means of convection, radiation and heat conduction;
the natural gas reacts with the combustion air instantaneously and then is evenly mixed with the flue gas;
the natural gas generates complete combustion reaction inside the burner;
irrespective of the effect of air humidity on the heat balance;
the kiln body is a closed container;
neglecting heat transfer of the roller way, and integrally observing the roller way and the product.
Preferably, step S1.3 is specifically:
Performing infinitesimal division on a firing interval, and establishing a flue gas temperature mechanism model in the roller kiln based on mass conservation and energy conservation:
conservation of mass:
assuming that the moisture of the ceramic is completely evaporated in a preheating zone, the micro-element division shows that the quality change of the flue gas on the micro-element delta x is derived from the inflow of natural gas and combustion air
In the method, in the process of the invention,the mass flow of the flue gas at the position x is expressed in kg/s; />The mass flow of the flue gas at the position x+Deltax is expressed in kg/s; />The input mass flow of combustion air at x is expressed in kg/s; />The input mass flow of natural gas at x is expressed in kg/s;
the conservation of energy includes conservation of flue gas energy, conservation of ceramic energy and conservation of wall energy, wherein:
the conservation of the energy of the flue gas is that the heat absorbed or lost by the flue gas is equal to the heat of convective heat transfer between the ceramic and the flue gas, the heat of convective heat transfer between the flue gas and the wall surface and the heat brought by natural gas and combustion air:
wherein, c pa The specific heat capacity of the flue gas is expressed as J/(kg.k); t (T) a The temperature of the flue gas is represented by k; t (T) c The temperature of the ceramic tile is expressed in k; t (T) w The temperature of the kiln wall surface is expressed as k; h is a c 、h w The convection heat transfer coefficients between the ceramics and the flue gas and between the wall surface and the flue gas are respectively expressed as w/(m) 2 ·k)；a c 、a w Respectively represent the heat transfer areas between ceramics and flue gas and between the wall surface and the flue gas,the unit is m 2 ；Q fg The heat brought by the natural gas and the combustion air comprises the heat release amount of the natural gas combustion and the sensible heat brought by the natural gas and the combustion air, and the unit is J;
the conservation of ceramic energy is that the heat absorbed or lost by the ceramic is equal to the heat of convective heat transfer between the ceramic and the flue gas and the radiant heat of the ceramic and the wall surface:
wherein:the mass flow rate of the ceramic is expressed in kg/s; c pc The specific heat capacity of the ceramic is expressed as J/(kg.k); epsilon c Emissivity of ceramic epsilon c =0.8 to 0.9; σ is the blackbody radiation constant, σ=5.67×10 -8 W/(m 2 ·K 4 )；
The conservation of wall energy comprises that under the steady state condition, the kiln wall is in a heat balance state, and the heat of convection heat transfer between smoke and the wall and the heat of radiation heat transfer between ceramics and the wall are equal to the heat dissipated by the wall and the outside air:
wherein k is w The thermal conductivity of the wall surface is w/(m.k); a, a w Is the area of the wall surface, m 2 ；
Wherein T is w | x K is the temperature of the inner surface of the kiln wall;the temperature K is the temperature of the outer surface of the kiln wall; Δy is the infinitesimal variation along the width of the kiln.
Preferably, step S1.4 is specifically:
dispersing the kiln length L of the roller kiln sintering section into a plurality of intervals which are mutually connected, and dispersing mass conservation and energy conservation equations at the nodes;
For heat conduction of the roller kiln sintering belt, the temperature changes along with time and space, and when the equation is discrete, not only space nodes but also time nodes are needed to be divided; dividing grids along the kiln length direction X according to the space step length delta X; starting from τ=0, the grid is divided by a time step Δτ, the coordinates of each dimension are as follows:
x i ＝i·Δx i＝1,2,…,I
τ＝j·Δτj＝1,2,…,J
Δx＝1/I
the differential equation of the internal node heat conduction of the roller kiln sintering belt can be expressed by a differential equation, the equation is discrete by adopting a Crank-Nicolson format, and the differential equation is solved by using a Gauss Seidel iteration method.
Preferably, the specific steps of the step S2 are as follows:
obtaining an error prediction model by utilizing local weighted kernel principal component regression, wherein the error prediction model outputs an error between an output result and an actual result of the mechanism model;
and utilizing the instant moving window judgment to adaptively update the error prediction model.
Preferably, the error prediction model is obtained by using local weighted kernel principal component regression, and the error between the output result and the actual result of the mechanism model is output by the error prediction model, specifically:
calculating the correlation weights of the input variable and the output variable by calculating Pearson coefficients:
Wherein ρ is x,y A Pearson correlation coefficient representing the input variable x and the output variable y;as the average value of the input variables,is the mean value of the output variables;
the correlation weight is:
correlation weight ρ m Representing the ratio of the weight of each input variable to the weight of the total input variable ρ x,y The larger the absolute value of (2) is, the stronger the correlation between the variables is, and ρ is m The larger the value of (c) is, the stronger the correlation of the input variable with the output variable among all variables is;
assuming that the ith input sample is X, its internal data is X 1 ,x 2 ,...,x n Written in matrix form as x= [ X 1 ,x 2 ,…,x n ] T The method comprises the steps of carrying out a first treatment on the surface of the Output samples y=y 1 ,y 2 ,…,y n Written in matrix form as y= [ Y ] 1 ,y 2 ,…,y n ] T The method comprises the steps of carrying out a first treatment on the surface of the Phi is the mapping function of the input variable projected into the high-dimensional space, and the corresponding sample of the high-dimensional space is phi (x i )；
The local weighting strategy is:
wherein x is i Represents the ith sample data, x ik Data representing a kth variable in the ith sample data;
assuming that phi (x) is the normalized data, the covariance matrix of the feature space:
λV＝CV
wherein lambda is the eigenvalue matrix of covariance matrix C, and V is the eigenvector matrix of covariance matrix C;
calculating a weighted kernel function, the kernel function should satisfyAnd (3) selecting a Gaussian kernel function for calculation:
then calculate byWeighted projection t of (2) i ：
The calculation process is as follows:
(1) Solving a eigenvalue eigenvector problem by eigenvalue decomposition, and arranging the obtained eigenvalues lambda in a descending order, wherein the corresponding eigenvectors are arranged similarly;
(2) The original data sample is described by taking the first d eigenvalues according to the contribution rate of the principal component, the corresponding data is the first d columns of matrix lambda and is marked as U d ＝[u 1 ,u 2 ,…,u d ]The first d eigenvalues U d The corresponding feature vectors are noted as
Projecting the dataset into a high-dimensional space:
wherein T is called a kernel principal component matrix, also called a score matrix, and the projected data is subjected to least squares regression, and a regression coefficient matrix theta is calculated as follows:
θ＝(T T T) -1 T T Y
wherein Y represents a related output data set and is subjected to normalization processing;
the predicted output of the data samples, i.e. the value of the error compensation, is therefore:
wherein the method comprises the steps ofIs the mean of the relevant dataset.
Preferably, the adaptive updating of the error prediction model is performed by utilizing the instant moving window judgment, specifically:
setting window length, upper threshold a 2 And a lower threshold value a 1 ；
Calculating the mean value mu and standard deviation sigma of samples in a window at the current moment;
calculating the mean value mu of samples in a next window at the next moment i And standard deviation sigma i ；
Judgment a 1 μ<μ i <a 2 And a 1 σ<σ i <a 2 Sigma, if yes, updating the error prediction model by using a sample in a next window at the next moment; if not, updating the step length, and recalculating the mean value mu and the standard deviation sigma of the samples in the window.
Preferably, when the step size is greater than or equal to the preset maximum step size, the error prediction model is updated by using samples in a next window at the next moment.
Compared with the prior art, the technical scheme of the invention has the beneficial effects that:
according to the invention, a mechanism model of the firing zone temperature of the ceramic roller kiln is established, meanwhile, a data modeling method is used for carrying out temperature compensation on errors of mechanism model results, and a mixed modeling strategy of mechanism modeling and data modeling is combined, so that a roller kiln firing zone temperature mixed prediction model is established, and the effectiveness of the model and the high precision of output results in the process of working condition change can be ensured.
Drawings
FIG. 1 is a schematic flow chart of the method of the present invention.
Fig. 2 is a schematic diagram of the overall structure of the roller kiln according to the embodiment.
Fig. 3 is a schematic diagram of a three-dimensional model of a roller kiln firing zone unit provided in the embodiment.
Fig. 4 is a schematic top view of hot air flow of a roller kiln firing belt.
Fig. 5 is a schematic diagram of the heat transfer pattern within the roller kiln firing zone.
Fig. 6 is a schematic diagram of the infinitesimal division of the firing window.
Fig. 7 is a one-dimensional unsteady conductive mesh partition diagram.
Fig. 8 is a numerical calculation flowchart.
Fig. 9 is a schematic diagram of the operation of the fixed-step moving window.
FIG. 10 is a flow chart diagram of an instant mobile window provided by an embodiment.
Fig. 11 is a schematic diagram of a roller kiln firing zone temperature mixing prediction model provided in the embodiment.
FIG. 12 is a schematic diagram of the flue gas mass flow of the firing zone.
Fig. 13 is a schematic diagram of a flue gas temperature variation curve.
FIG. 14 is a graph showing the error results output by the MWJIT-LWKPCR model.
FIG. 15 is a graph showing the comparison of the predicted value and the actual temperature value of the MWJIT-LWKPCR model.
FIG. 16 is a graph comparing predicted and actual values of the LWKPCR model.
FIG. 17 is a graph comparing predicted and actual values of MW-LWKPCR model.
Detailed Description
The drawings are for illustrative purposes only and are not to be construed as limiting the present patent;
for the purpose of better illustrating the embodiments, certain elements of the drawings may be omitted, enlarged or reduced and do not represent the actual product dimensions;
it will be appreciated by those skilled in the art that certain well-known structures in the drawings and descriptions thereof may be omitted.
The technical scheme of the invention is further described below with reference to the accompanying drawings and examples.
Example 1
The embodiment provides a method for predicting the firing zone temperature of a roller kiln based on mechanism and data mixed driving, which is shown in figure 1 and comprises the following steps:
S1: modeling a mechanism model of the mechanism process of the roller kiln firing zone by using a mass conservation law and an energy conservation law and solving by using a finite difference method;
s2: based on the instant moving window and the local weighted kernel principal component regression, establishing an error compensation model based on data driving;
s3: combining the mechanism model obtained in the step S1 with the error compensation model obtained in the step S2 to obtain a roller kiln firing zone temperature mixed prediction model;
s4: and predicting the temperature of the firing zone of the roller kiln by using a mixed prediction model of the temperature of the firing zone of the roller kiln.
Preferably, the step S1 specifically includes the following steps:
s1.1: analyzing the heat transfer process in the firing zone of the roller kiln, and establishing a roller kiln heat transfer process model;
s1.2: estimating state variables of the roller kiln firing belt, and setting assumption conditions of the state variables;
s1.3: based on mass conservation and energy conservation, establishing a flue gas temperature mechanism model of the roller kiln sintering belt;
s1.4: and carrying out numerical calculation and solving on the flue gas temperature mechanism model by adopting a finite difference method.
The overall structure layout of the roller kiln is shown in fig. 2, the firing zone of the roller kiln is assembled by a plurality of unit kilns, the three-dimensional structure of the unit kilns is shown in fig. 3, and the roller kiln is easy for ceramic enterprises to transport and install. 8 burners which are bilaterally symmetrical are respectively arranged outside the unit kiln, each burner is connected with a fuel and combustion-supporting air control system and uniformly regulated and controlled, and the fuel and the combustion-supporting air are fully mixed and combusted in the burners and then are sprayed into the kiln. The complex physical and chemical reaction of the ceramic product in the firing process mainly occurs in a firing zone, and the combustion reaction of natural gas and combustion supporting air also occurs in the firing zone, wherein the flow direction of high-temperature flue gas in a single kiln body of the ceramic roller kiln firing zone (from right to left) is the direction from the kiln outlet to the kiln inlet and is just opposite to the movement direction of the ceramic product (from left to right) as shown in figure 4.
Heat transfer within the firing zone occurs primarily between the ceramic article, flue gas, kiln walls and roller tables. In order to facilitate the calculation, the roller way and the ceramic product are regarded as a whole when the temperature distribution of the flue gas in the sintering zone is calculated.
The step S1.1 specifically comprises the following steps:
the roller kiln sintering in-band heat transfer mode comprises convection heat transfer, radiation heat transfer and heat transfer, wherein:
the convection heat transfer is the main heat transfer mode in the kiln, most of heat is transferred among different objects in the mode, the convection heat transfer comprises the convection heat transfer quantity between smoke and products and the convection heat transfer quantity between smoke and kiln walls, and the formula of the convection heat transfer is as follows:
Q convection current ＝h(t-t m )F
In which Q Convection current The total heat of convection heat transfer is expressed as W; h is the convection heat transfer coefficient, and the unit is W/(m) 2 K); t is the temperature of the surface of the product, and the unit is K; t is t m The unit is K, which is the temperature of the flue gas; f is the heating surface area of the product, and the unit is m 2 ；
The radiation heat transfer quantity inside the sintering belt is large, the radiation heat transfer is divided into gas radiation and solid radiation, the radiation heat transfer comprises radiation heat transfer between a product and a kiln wall and radiation heat transfer between a transmission roller and the kiln wall, and the radiation heat transfer has the following formula:
In which Q Radiation of Is a gasHeat radiated to the surface of the article in W; epsilon c Surface blackness of the product; f is the heated area of the product, and the unit is m 2 The method comprises the steps of carrying out a first treatment on the surface of the Sigma is the blackness of the gas; t is the absolute temperature of the gas, and the unit is K; t (T) p The absolute temperature of the surface of the product is expressed as K;
the heat conduction is also a common heat transfer mode in the sintering zone, and comprises the heat transfer from a temperature zone with higher temperature to a low temperature zone, the heat absorbed by a product and a kiln wall and the heat diffusion of the kiln wall, wherein the heat conduction formula is as follows:
in which Q Conduction of Is the heat flux density; k is the thermal conductivity; t is the temperature; x is the length in the heat transfer direction.
During the working process of the sintering belt, the three heat transfer modes exist at any moment, so that the heat transfer modes are taken into consideration when the roller kiln heat transfer process model is built. The firing belt heat transfer mode of the roller kiln is shown in figure 5.
The step S1.2 specifically comprises the following steps:
the assumption conditions for setting the state variables include:
the system is in a stable running state;
the change of the temperature and the gas flow in the vertical direction is not considered, and only the change of the temperature and the gas flow in the length direction of the kiln is considered;
heat is transferred between the flue gas, the ceramic product and the kiln wall by means of convection, radiation and heat conduction;
The natural gas reacts with the combustion air instantaneously and then is evenly mixed with the flue gas;
the natural gas generates complete combustion reaction inside the burner;
irrespective of the effect of air humidity on the heat balance;
the kiln body is a closed container;
neglecting heat transfer of the roller way, and integrally observing the roller way and the product.
The step S1.3 specifically comprises the following steps:
burners are regularly distributed on two sides of a firing zone, natural gas and combustion air are introduced into the burners and are sprayed out after being combusted in the burners, and the flow of high-temperature flue gas can cause the change of flue gas flow, energy and the like in the firing zone, so that the firing zone needs to be subjected to micro-element division when the flue gas temperature in the firing zone is studied, as shown in fig. 6, the firing zone is subjected to micro-element division, and a flue gas temperature mechanism model in a roller kiln is established based on mass conservation and energy conservation:
conservation of mass:
assuming that the moisture of the ceramic is completely evaporated in a preheating zone, the micro-element division shows that the quality change of the flue gas on the micro-element delta x is derived from the inflow of natural gas and combustion air
In the method, in the process of the invention,the mass flow of the flue gas at the position x is expressed in kg/s; />The mass flow of the flue gas at the position x+Deltax is expressed in kg/s; / >The input mass flow of combustion air at x is expressed in kg/s; />Represents the input mass flow of natural gas at x in units ofkg/s；
The conservation of energy includes conservation of flue gas energy, conservation of ceramic energy and conservation of wall energy, wherein:
the conservation of the energy of the flue gas is that the heat absorbed or lost by the flue gas is equal to the heat of convective heat transfer between the ceramic and the flue gas, the heat of convective heat transfer between the flue gas and the wall surface and the heat brought by natural gas and combustion air:
wherein, c pa The specific heat capacity of the flue gas is expressed as J/(kg.k); t (T) a The temperature of the flue gas is represented by k; t (T) c The temperature of the ceramic tile is expressed in k; t (T) w The temperature of the kiln wall surface is expressed as k; h is a c 、h w The convection heat transfer coefficients between the ceramics and the flue gas and between the wall surface and the flue gas are respectively expressed as w/(m) 2 ·k)；a c 、a w Respectively represent the heat transfer area between ceramics and smoke and between wall surface and smoke, the unit is m 2 ；Q fg The heat brought by the natural gas and the combustion air comprises the heat release amount of the natural gas combustion and the sensible heat brought by the natural gas and the combustion air, and the unit is J;
the conservation of ceramic energy is that the heat absorbed or lost by the ceramic is equal to the heat of convective heat transfer between the ceramic and the flue gas and the radiant heat of the ceramic and the wall surface:
wherein: The mass flow rate of the ceramic is expressed in kg/s; c pc The specific heat capacity of the ceramic is expressed as J/(kg.k); epsilon c Emissivity of ceramic epsilon c =0.8 to 0.9; σ is the blackbody radiation constant, σ=5.67×10 -8 W/(m 2 ·K 4 )；
The conservation of wall energy comprises that under the steady state condition, the kiln wall is in a heat balance state, and the heat of convection heat transfer between smoke and the wall and the heat of radiation heat transfer between ceramics and the wall are equal to the heat dissipated by the wall and the outside air:
wherein k is w The thermal conductivity of the wall surface is w/(m.k); a, a w Is the area of the wall surface, m 2 ；
Wherein T is w | x K is the temperature of the inner surface of the kiln wall;the temperature K is the temperature of the outer surface of the kiln wall; Δy is the infinitesimal variation along the width of the kiln.
The step S1.4 specifically comprises the following steps:
dispersing the kiln length L of the roller kiln sintering section into a plurality of intervals which are mutually connected, and dispersing mass conservation and energy conservation equations at the nodes;
for heat conduction of the roller kiln sintering belt, the temperature changes along with time and space, and when the equation is discrete, not only space nodes but also time nodes are needed to be divided; the coordinate dispersion of the time dimension and the space dimension is shown in fig. 7.
Dividing grids along the kiln length direction X according to the space step length delta X; starting from τ=0, the grid is divided by a time step Δτ, the coordinates of each dimension are as follows:
x i ＝i·Δx i＝1,2,…,I
τ＝j·Δτj＝1,2,…,J
Δx＝1aI
The basic idea of finite difference is based on the taylor series, the taylor expansion being as follows:
therefore, the differential equation of the internal node heat conduction of the roller kiln sintering belt can be expressed by a differential equation, the equation is discrete by adopting a Crank-Nicolson format, and the differential equation is solved by using a Gaussidel iteration method.
The partial differential format selects the Crank-Nicolson format:
a Gaussian Seidel iterative method is used to solve the discrete system of equations of the Crank-Nicolson format. The numerical calculation flow is shown in fig. 8.
The specific steps of the step S2 are as follows:
obtaining an error prediction model by utilizing local weighted kernel principal component regression, wherein the error prediction model outputs an error between an output result and an actual result of the mechanism model;
and utilizing the instant moving window judgment to adaptively update the error prediction model.
Kernel Principal Component Regression (KPCR) is an effective nonlinear data prediction method whose basic idea is to map data from a low-dimensional space into a high-dimensional feature space, i.e., (where Rn represents the original dataset and represents nonlinear mapping), perform principal component analysis in the high-dimensional space, and then build a least squares regression model to predict based on the extracted feature and quality variable data. The Local Weighted Kernel Principal Component Regression (LWKPCR) considers the influence of the correlation of the input variable and the output variable of the sample on the feature extraction on the basis of KPCR, and the higher the correlation of the input variable and the output variable is, the larger the influence on the quality prediction is, so that the correlation weight of the input variable and the output variable is added into the feature learning of KPCR.
The error prediction model is obtained by utilizing local weighted kernel principal component regression, and the error between the output result and the actual result of the mechanism model is output by the error prediction model, specifically:
calculating the correlation weights of the input variable and the output variable by calculating Pearson coefficients:
wherein ρ is x,y A Pearson correlation coefficient representing the input variable x and the output variable y;as the average value of the input variables,is the mean value of the output variables;
the correlation weight is:
correlation weight ρ m Representing the ratio of the weight of each input variable to the weight of the total input variable ρ x,y The larger the absolute value of (2) is, the stronger the correlation between the variables is, and ρ is m The larger the value of (c) is, the stronger the correlation of the input variable with the output variable among all variables is;
assuming that the ith input sample is X, its internal data is X 1 ,x 2 ,...,x n Written in matrix form as x= [ X 1 ,x 2 ,…,x n ] T The method comprises the steps of carrying out a first treatment on the surface of the Output samples y=y 1 ,y 2 ,…,y n Written in matrix formY= [ Y ] 1 ,y 2 ,…,y n ] T The method comprises the steps of carrying out a first treatment on the surface of the Phi is the mapping function of the input variable projected into the high-dimensional space, and the corresponding sample of the high-dimensional space is phi (x i )；
The local weighting strategy is:
wherein x is i Represents the ith sample data, x ik Data representing a kth variable in the ith sample data; assuming that phi (x) is the normalized data, the covariance matrix of the feature space:
λV＝CV
Wherein lambda is the eigenvalue matrix of covariance matrix C, and V is the eigenvector matrix of covariance matrix C;
calculating a weighted kernel function, the kernel function should satisfyAnd (3) selecting a Gaussian kernel function for calculation:
then calculate byWeighted projection t of (2) i ：
The calculation process is as follows:
(1) Solving a eigenvalue eigenvector problem by eigenvalue decomposition, and arranging the obtained eigenvalues lambda in a descending order, wherein the corresponding eigenvectors are arranged similarly;
(2) The original data sample is described by taking the first d eigenvalues according to the contribution rate of the principal component, the corresponding data is the first d columns of matrix lambda and is marked as U d ＝[u 1 ,u 2 ,…,u d ]The first d eigenvalues U d The corresponding feature vectors are noted as
Projecting the dataset into a high-dimensional space:
wherein T is called a kernel principal component matrix, also called a score matrix, and the projected data is subjected to least squares regression, and a regression coefficient matrix theta is calculated as follows:
θ＝(T T T) -1 T T Y
wherein Y represents a related output data set and is subjected to normalization processing;
the predicted output of the data samples, i.e. the value of the error compensation, is therefore:
/>
wherein the method comprises the steps ofIs the mean of the relevant dataset.
And utilizing an immediate Moving Window (MWJIT) to judge to adaptively update the error prediction model, wherein the method specifically comprises the following steps:
setting window length, upper threshold a 2 And a lower threshold value a 1 ；
Calculating the mean value mu and standard deviation sigma of samples in a window at the current moment;
calculating the mean value mu of samples in a next window at the next moment i And standard deviation sigma i ；
Judgment a 1 μ<μ i <a 2 And a 1 σ<σ i <a 2 Sigma, if yes, updating the error prediction model by using a sample in a next window at the next moment; if not, updating the step length, and recalculating the mean value mu and the standard deviation sigma of the samples in the window.
In the actual production process, the working conditions of the roller kiln, such as raw material batch, air pressure, gas flow and the like, tend to change slowly, and the change of the new working conditions can cause a certain deviation between real-time data and historical data, so that the accuracy of a prediction result can be reduced by continuously using a historical model for calculation. Therefore, a predictive model is required to be updated immediately, ensuring long-term validity of the model. When the self-adaptive updating is carried out, the moving window technology can continuously track the updating state of the process, and the change of the actual working condition is adapted by continuously updating the error compensation model.
The working principle of the traditional fixed step length moving window is as follows: the data samples in the first time window are { X ] w1 ,Y w1 The window length is L, and assuming the step size is D, the data sample of the second time window is { X } by moving forward by D in the second time window w2 ,Y w2 It follows that the second window discards data in the length D forward compared to the first window and then incorporates data of the same length D later to form a new data sample, as shown in fig. 9.
The instant Moving Window (MWJIT) works similarly to the fixed-step moving window, but the step size of the instant moving window is time-varying, assuming that the step size of the instant moving window is D i At this time D i The value of (2) is determined by the change condition of the data sample. Assume that the data sample of the ith time window is { X ] wi ,Y wi The mean mu and standard deviation sigma corresponding to the data sample can be calculated, and then whether to update the model is selected according to whether the mean and standard deviation of the subsequent sample data are consistent with the mean and variance of the samples in the window. Assuming that the upper and lower limits of the mean mu and the standard deviation sigma are 1.2 and 0.8, the algorithm flow of the corresponding instant moving windowAs shown in fig. 10
And when the step length is larger than or equal to the preset maximum step length, updating the error prediction model by using a sample in a next window at the next moment.
The running flow of the instant Moving Window (MWJIT) is as follows:
training a model of sample data of a first window, and under the condition that the data sample is not changed greatly, predicting by using a previous training model, namely, not updating the model; and when the statistics of the data samples change greatly, the window moves, and the Local Weighted Kernel Principal Component Regression (LWKPCR) model is updated. The amount of data samples between these two moving window local weighted kernel principal component regression model updates is then referred to as the step size (D) of the update, and the number of samples of the step size is used to calculate the predictions for the last updated model.
When the LWKPCR model is updated, whether the previously trained model is suitable for subsequent temperature prediction is judged, and whether the statistic of temperature error data is greatly changed is mainly judged, wherein the mean value mu and the standard deviation sigma are used as judgment bases. Taking the mean value mu and the standard deviation sigma of the trained LWKPCR model data as threshold values, and setting a certain upper limit and a certain lower limit (a 1 、a 2 ). If the new sample meets the condition:
a 1 μ<μ i <a 2 μand a 1 σ<σ i <a 2 σ
then it is indicated that the subsequent data samples are within acceptable variation and can be predicted using the previously trained LWKPCR model, where a 1 、a 2 Is a preset proportionality coefficient. And when the condition is not satisfied, adjusting the position of the moving window, and updating the LWKPCR model. In addition, sometimes the LWKPCR model is not updated all the time, possibly because the industrial process is relatively stable, and long-term data samples meet the conditions, thus requiring setting of an upper step size limit D max The LWKPCR model can be updated in time, so that timeliness of the model is guaranteed.
According to the temperature mechanism model established in the step S1, the error of the mechanism model is modeled by utilizing the instantaneous moving window local weighted kernel principal component regression model established in the step S2, and the two components form a hybrid modeling framework for temperature prediction based on a mechanism and data-driven hybrid model, wherein the hybrid modeling framework for temperature prediction is shown in the figure 11.
(1) Calculating to obtain data of temperature distribution in a sintering zone through a mechanism model of S1 by utilizing data acquired in real time in production;
(2) And modeling the error of the mechanism model of the S1 by using the error compensation model established by the S2.
S2, predicting value output by error compensation modelFor the error prediction value of the mechanism model and the actual temperature of the S1, the temperature of the mechanism is added to accurately represent the flue gas temperature of the firing zone, and the calculation formula is as follows:
wherein the method comprises the steps ofIs a predicted value of temperature; y is m Solving the mechanism model; />Predicted values for the algorithm solutions.
Taking a ceramic roller kiln as an example, implementation of the scheme is performed.
(1) Establishing a mechanism model
The actual production data of a ceramic enterprise in the bergamot is adopted as comparison, and in the data acquisition process, the enterprise places 11 temperature probes on the inner wall of a kiln body, and the probes are positioned on a plane 50mm below a nozzle. The temperature probes are uniformly arranged along the length direction of the kiln, and temperature data are acquired every 5 minutes. Equation solving is carried out by taking the data of one kiln section, so that the distribution of the temperature of high-temperature flue gas in the roller kiln firing zone along the length direction of the kiln can be effectively calculated.
The relevant parameters related to calculation are obtained by enterprise providing and consulting literature, the time step of solving is set to be 0.05, and the space step is set to be 0.01. Wherein the specification of the ceramic tile is 600mm, 11.31mm, the mass is 25.5kg/m < 2 >, the fuel consumption is 330kg/h, and the fuel gas is natural gas. The input and output parameters of the model are shown in table 1, the performance parameters of the kiln wall heat insulation material are shown in table 2, and the heat transfer coefficients of the outer surface of the kiln wall and air are shown in table 3:
Table 1 kiln firing zone input and output parameters
TABLE 2 kiln wall insulation material performance parameters
TABLE 3 heat transfer coefficient of kiln wall outer surface and air
The mass flow of the flue gas in the kiln is determined by the mass flow of the inlet flue gas, the mass flow of the natural gas and the mass flow of the combustion air together, as shown in fig. 12;
inputting parameters into a model, solving a differential equation through MATLAB programming, and calculating a result shown in FIG. 13;
it can be seen from fig. 13 that the flue gas temperature slowly rises in the kiln length direction in a section of the kiln body, the actual temperature profile will fluctuate at different points and the fluctuation range is not large, whereas the mechanism temperature profile does not appear as it is assumed that the natural gas reacts instantaneously with the combustion air and is uniformly mixed with the flue gas before solving, and the heat transfer is performed under ideal conditions, so that no abnormal fluctuation of the flue gas temperature occurs at a certain point.
The fluctuation in actual production may be due to abnormality in the firing process or that natural gas and combustion air are not uniformly mixed with flue gas in time after being combusted. The mechanism temperature is slowly increased and is similar to the temperature trend of an actual kiln in operation, and the error range is controlled within 0-30 ℃, so that the solving method and the result of the differential equation are consistent with the trend of the actual temperature change.
(2) Fitting error by instantaneous moving window local weighting kernel principal component regression model
It can be seen from fig. 13 that the calculation result of the mechanism model has a certain error with the actual temperature, so that further analysis of the error is required.
The method for the regression of the principal components of the instant moving window local weighting kernel (MWJIT-LWKPCR) adopted by the invention needs to determine the proportionality coefficient a 1 、a 2 And step upper limit D max Set to 0.8, 1.2, 200, respectively; the window length L is set to 200; setting the width parameter of the kernel function to 3 according to a trial-and-error method; the score threshold was set at 80% when the principal component was selected by CPV method at modeling.
The collected 1100 sets of input data are imported into the model, wherein the first 200 sets of data samples are selected as a training set for training of the model, and then the remaining 900 sets of samples are selected as a test set. The error fit results are shown in fig. 14.
(3) Creation of hybrid models
The mixed model consists of a mechanism model and a temperature error value, and outputs the resultWherein->Final temperature data obtained for the hybrid model; y is m Temperature data output for the mechanism model, +.>And a temperature error compensation value predicted for the error compensation model. The comparison of the temperature data obtained by the hybrid model with the actual workpiece temperature is shown in fig. 15.
For convenience of comparison, fig. 15 shows the coordinates of the collection point as follows: the center point of the kiln in the longitudinal direction (X axis), the center of the kiln in the width direction (Y axis) and the center point of the height (Z axis) are 10CM above.
To verify the advantages of the instant moving window kernel principal component regression (MWJIT-LWKPCR) proposed by the present invention, it was compared to the effects of Local Weighted Kernel Principal Component Regression (LWKPCR) and moving window local weighted kernel principal component regression (MW-LWKPCR).
The LWKPCR method is an offline algorithm, the first 200 data samples are selected as a training set for training of the model, and then the remaining 900 sets of samples are used as a test set. The window lengths L of MW-LWKPCR and MWJIT-LWKPCR are also set to 200. And setting the width parameter of the kernel function to 3 according to a trial and error method, and setting the score threshold value to 80% when the principal component is selected through a CPV method in modeling. In the MW-LWKPCR method, the moving window step length D is set to 100, and the proportionality coefficient a 1 、a 2 And step upper limit D max Set to 0.8, 1.2, 200, respectively.
The parameter settings are shown in the following table:
table 4 preset parameters for three algorithms
In performance evaluation of the three algorithms LWKPCR, MW-LWKPCR, and MWJIT-LWKPCR, root Mean Square Error (RMSE), mean Absolute Error (MAE), and maximum absolute error (MAX) were used as performance indicators. The specific calculation formula is as follows:
Root mean square error:
average absolute error:
/>
maximum absolute error:
wherein y is i As an actual value of the kiln temperature,is a predictive value of the hybrid model.
The average absolute error and the maximum absolute error represent the deviation between the predicted temperature and the actual temperature, and the maximum absolute error represents the maximum deviation between the predicted temperature and the actual temperature, and the two values can reflect the performance of the algorithm to a certain extent. The root mean square error characterizes the predictive power of the regression model, and the stronger the predictive power is, the closer the value is to zero.
Fig. 16 is a graph comparing predicted values and actual temperature values obtained using the LWKPCR algorithm, and fig. 17 is a graph comparing results of the MW-LWKPCR algorithm.
The model performance results for the three methods are shown in the following table:
table 5 comparison of predicted results
As can be seen from comparison of results, the LWKPCR belongs to an offline learning algorithm, so that the prediction performance is continuously reduced along with the change of working conditions when the LWKPCR is used for temperature prediction, and the prediction result obtained along with the change of the conditions is inevitably not in accordance with actual production requirements.
The MW-LWKPCR algorithm updates the prediction model at intervals, so that the effectiveness of the model can be maintained, but because the updating is timed, a great deal of calculation resources are wasted when the model is updated under the condition of similar working conditions, and when the working conditions change frequently, the MW-LWKPCR algorithm can reduce the effectiveness of the model because the updating cannot be responded in time.
The MWJIT-LWKPCR algorithm provided by the method selects an updated model according to the change of working conditions and data, the effectiveness of the model can be kept at all times, meanwhile, the result comparison can show that the prediction performance of the MWJIT-LWKPCR algorithm is obviously superior to that of the LWKPCR algorithm, and the calculation speed is improved by about 40% compared with that of the MW-LWKPCR algorithm although the prediction result is not quite different, so that the MWJIT-LWKPCR algorithm can better reflect the change of the temperature of a firing zone of the roller kiln, and can meet the actual production requirements of the roller kiln system.
The same or similar reference numerals correspond to the same or similar components;
the terms describing the positional relationship in the drawings are merely illustrative, and are not to be construed as limiting the present patent;
it is to be understood that the above examples of the present invention are provided by way of illustration only and not by way of limitation of the embodiments of the present invention. Other variations or modifications of the above teachings will be apparent to those of ordinary skill in the art. It is not necessary here nor is it exhaustive of all embodiments. Any modification, equivalent replacement, improvement, etc. which come within the spirit and principles of the invention are desired to be protected by the following claims.
Claims (8)
1. The method for predicting the firing zone temperature of the roller kiln based on the mechanism and data mixed driving is characterized by comprising the following steps of:
s1: modeling a mechanism model of the mechanism process of the roller kiln firing zone by using a mass conservation law and an energy conservation law and solving by using a finite difference method;
s2: based on the instant moving window and the local weighted kernel principal component regression, establishing an error compensation model based on data driving;
s3: combining the mechanism model obtained in the step S1 with the error compensation model obtained in the step S2 to obtain a roller kiln firing zone temperature mixed prediction model;
s4: predicting the temperature of the firing zone of the roller kiln by using a mixed prediction model of the temperature of the firing zone of the roller kiln;
the step S1 specifically comprises the following steps:
s1.1: analyzing the heat transfer process in the firing zone of the roller kiln, and establishing a roller kiln heat transfer process model;
s1.2: estimating state variables of the roller kiln firing belt, and setting assumption conditions of the state variables;
s1.3: based on mass conservation and energy conservation, establishing a flue gas temperature mechanism model of the roller kiln sintering belt;
s1.4: carrying out numerical calculation and solving on the flue gas temperature mechanism model by adopting a finite difference method;
the step S1.1 specifically comprises the following steps:
The roller kiln sintering in-band heat transfer mode comprises convection heat transfer, radiation heat transfer and heat transfer, wherein:
the convective heat transfer comprises the convective heat transfer quantity between the flue gas and the product and the convective heat transfer quantity between the flue gas and the kiln wall, and the formula of the convective heat transfer is as follows:
Q convection current ＝h(t-t m )F
In which Q Convection current The total heat of convection heat transfer is expressed as W; h is the convection heat transfer coefficient, and the unit is W/(m) 2 K); t is the temperature of the surface of the product, and the unit is K; t is t n The unit is K, which is the temperature of the flue gas; f is the heating surface area of the product, and the unit is m 2 ；
The radiant heat transfer includes the radiant heat transfer between the article and the kiln wall and the radiant heat transfer between the transfer roller and the kiln wall, the radiant heat transfer having the formula:
in which Q Radiation of Heat radiated to the surface of the product by the gas is expressed as W; epsilon c Surface blackness of the product; f is the heating surface area of the product, and the unit is m 2 The method comprises the steps of carrying out a first treatment on the surface of the Sigma is the blackness of the gas; t is the absolute temperature of the gas, singlyThe bit is K; t (T) ρ The absolute temperature of the surface of the product is expressed as K;
the heat conduction comprises heat transfer from a high-temperature region to a low-temperature region, heat absorbed by the product and the kiln wall and heat diffusion of the kiln wall, and the formula of the heat conduction is as follows:
in which Q Conduction of Is the heat flux density; k is the thermal conductivity; t is the absolute temperature of the gas; x is the length in the heat transfer direction.
2. The method for predicting the firing zone temperature of the roller kiln based on the mechanism and data mixed driving of claim 1, wherein the step S1.2 is specifically:
the assumption conditions for setting the state variables include:
the system is in a stable running state;
the change of the temperature and the gas flow in the vertical direction is not considered, and only the change of the temperature and the gas flow in the length direction of the kiln is considered;
heat is transferred between the flue gas, the ceramic product and the kiln wall by means of convection, radiation and heat conduction;
the natural gas reacts with the combustion air instantaneously and then is evenly mixed with the flue gas;
the natural gas generates complete combustion reaction inside the burner;
irrespective of the effect of air humidity on the heat balance;
the kiln body is a closed container;
and neglecting heat transfer of the roller way, and considering the roller way and the product as a whole.
3. The method for predicting the firing zone temperature of the roller kiln based on the mixed driving of the mechanism and the data according to claim 1 or 2, wherein the step S1.3 is specifically:
performing infinitesimal division on a firing interval, and establishing a flue gas temperature mechanism model in the roller kiln based on mass conservation and energy conservation:
conservation of mass:
assuming that the moisture of the ceramic is completely evaporated in a preheating zone, the micro-element division shows that the quality change of the flue gas on the micro-element delta x is derived from the inflow of natural gas and combustion air
In the method, in the process of the invention,the mass flow of the flue gas at the position x is expressed in kg/s; />The mass flow of the flue gas at the position x+Deltax is expressed in kg/s; />Representing the input mass flow of combustion air in kg/s at position x; />Represents the input mass flow of natural gas in kg/s at location x;
the conservation of energy includes conservation of flue gas energy, conservation of ceramic energy and conservation of wall energy, wherein:
the conservation of the energy of the flue gas is that the heat absorbed or lost by the flue gas is equal to the heat of convective heat transfer between the ceramic and the flue gas, the heat of convective heat transfer between the flue gas and the wall surface and the heat brought by natural gas and combustion air:
wherein, c pa The specific heat capacity of the flue gas is expressed as J/(kg.k); t (T) a The temperature of the flue gas is represented by k; t (T) c The temperature of the ceramic tile is expressed in k; t (T) w The temperature of the kiln wall surface is expressed as k; h is a c 、h w The convection heat transfer coefficients between the ceramics and the flue gas and between the wall surface and the flue gas are respectively expressed as w/(m) 2 •k)；a c 、a w Respectively represent the heat transfer area between ceramics and smoke and between wall surface and smoke, the unit is m 2 ；Q fg The heat brought by the natural gas and the combustion air comprises the heat release amount of the natural gas combustion and the sensible heat brought by the natural gas and the combustion air, and the unit is J;
The conservation of ceramic energy is that the heat absorbed or lost by the ceramic is equal to the heat of convective heat transfer between the ceramic and the flue gas and the radiant heat of the ceramic and the wall surface:
wherein:the mass flow rate of the ceramic is expressed in kg/s; c pc The specific heat capacity of the ceramic is expressed as J/(kg.k); epsilon c Emissivity of ceramic epsilon c ＝0.8～0.9；σ 0 Is a blackbody radiation constant, sigma 0 ＝5.67×10 -8 W/(m 2 ·K 4 )；
The conservation of wall energy comprises that under the steady state condition, the kiln wall is in a heat balance state, and the heat of convection heat transfer between smoke and the wall and the heat of radiation heat transfer between ceramics and the wall are equal to the heat dissipated by the wall and the outside air:
wherein k is w The unit is w/(m.k) for wall surface heat conductivity; a, a w Is the heat transfer area between the wall surface and the flue gas, and the unit is m 2 ；
Wherein T is w | x The unit is K, which is the temperature of the inner surface of the kiln wall;the unit is K, which is the temperature of the outer surface of the kiln wall; Δy is the infinitesimal variation along the width of the kiln.
4. The method for predicting the firing zone temperature of the roller kiln based on the mixed driving of the mechanism and the data according to claim 3, wherein the step S1.4 is specifically:
dispersing the kiln length L of the roller kiln sintering section into a plurality of intervals which are mutually connected, and dispersing mass conservation and energy conservation equations at the nodes;
For heat conduction of the roller kiln sintering belt, the temperature changes along with time and space, and when the equation is discrete, not only space nodes but also time nodes are needed to be divided; dividing grids along the kiln length direction X according to the space step length delta X; starting from τ=0, the grid is divided by a time step Δτ, the coordinates of each dimension are as follows:
x i ＝i·Δx,i＝1,2,…,I
τ＝j·Δτ,j＝1,2,…,J
Δx＝1/I
the differential equation of the internal node heat conduction of the roller kiln sintering belt can be expressed by a differential equation, the equation is discrete by adopting a Crank-Nicolson format, and the differential equation is solved by using a Gauss Seidel iteration method.
5. The method for predicting the firing zone temperature of the roller kiln based on the mechanism and data mixed driving according to claim 1, wherein the specific steps of the step S2 are as follows:
obtaining an error prediction model by utilizing local weighted kernel principal component regression, wherein the error prediction model outputs an error between an output result and an actual result of the mechanism model;
and utilizing the instant moving window judgment to adaptively update the error prediction model.
6. The method for predicting the firing zone temperature of the roller kiln based on mechanism and data mixed driving according to claim 5, wherein the error prediction model is obtained by utilizing local weighted kernel principal component regression, and the error between the output result and the actual result of the mechanism model is output by the error prediction model, specifically:
Calculating the correlation weights of the input variable and the output variable by calculating Pearson coefficients:
wherein ρ is X',Y' Pearson correlation coefficients representing the input variable X 'and the output variable Y';as the average value of the input variables,is the mean value of the output variables;
the correlation weight is:
correlation weight ρ m Representing the weight of each input variable and the total input variationRatio of the quantitative weights ρ X',Y' The larger the absolute value of (2) is, the stronger the correlation between the variables is, and ρ is m The larger the value of (c) is, the stronger the correlation of the input variable with the output variable among all variables is;
assuming that the ith input sample is X, its internal data is X 1 ,x 2 ,...,x n Written in matrix form as x= [ X 1 ,x 2 ,…,x n ] T The method comprises the steps of carrying out a first treatment on the surface of the Output samples y=y 1 ,y 2 ,…,y n Written in matrix form as y= [ Y ] 1 ,y 2 ,…,y n ] T The method comprises the steps of carrying out a first treatment on the surface of the Phi is the mapping function of the input variable projected into the high-dimensional space, and the corresponding sample of the high-dimensional space is phi (x i )；
The local weighting strategy is:
wherein x is i Represents the ith sample data, x ik Data representing a kth variable in the ith sample data;
assuming that phi (x) is the normalized data, the covariance matrix of the feature space:
λV＝CV
wherein lambda is the eigenvalue matrix of covariance matrix C, and V is the eigenvector matrix of covariance matrix C;
calculating a weighted kernel function, the kernel function should satisfy And (3) selecting a Gaussian kernel function for calculation:
then calculate byWeighted projection t of (2) i ：
The calculation process is as follows:
(1) Solving a eigenvalue eigenvector problem by eigenvalue decomposition, and arranging the obtained eigenvalues lambda in a descending order, wherein the corresponding eigenvectors are arranged similarly;
(2) The original data sample is described by taking the first d eigenvalues according to the contribution rate of the principal component, the corresponding data is the first d columns of matrix lambda and is marked as U d ＝[u 1 ,u 2 ,…,u d ]The first d eigenvalues U d The corresponding feature vectors are noted as
Projecting the dataset into a high-dimensional space:
wherein T is 1 The method is called a kernel principal component matrix, also called a score matrix, and performs least squares regression on the projected data to calculate a regression coefficient matrix theta as follows:
θ＝(T 1 T T 1 ) -1 T 1 T Y
wherein Y represents a related output data set and is subjected to normalization processing;
the predicted output of the data samples, i.e. the value of the error compensation, is therefore:
wherein the method comprises the steps ofIs the mean of the relevant dataset.
7. The mechanism and data hybrid driving-based roller kiln firing zone temperature prediction method according to claim 6, wherein the error prediction model is adaptively updated by utilizing instant moving window judgment, and specifically comprises the following steps:
setting window length, upper threshold a 2 And a lower threshold value a 1 ；
Calculating the mean value mu and standard deviation sigma of samples in a window at the current moment;
calculating the mean value mu of samples in a next window at the next moment i And standard deviation sigma i ；
Judgment a 1 μ＜μ i ＜a 2 And a 1 σ＜σ i ＜a 2 Sigma, if yes, updating the error prediction model by using a sample in a next window at the next moment; if not, updating the step length, and recalculating the mean value mu and the standard deviation sigma of the samples in the window.
8. The mechanism and data hybrid drive based method for predicting the firing zone temperature of a roller kiln according to claim 7, wherein when the step length is equal to or greater than a preset maximum step length, the error prediction model is updated by using samples in a next window at the next time.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111033301.0A CN113627064B (en) | 2021-09-03 | 2021-09-03 | Roller kiln firing zone temperature prediction method based on mechanism and data mixed driving |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111033301.0A CN113627064B (en) | 2021-09-03 | 2021-09-03 | Roller kiln firing zone temperature prediction method based on mechanism and data mixed driving |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113627064A CN113627064A (en) | 2021-11-09 |
| CN113627064B true CN113627064B (en) | 2023-11-21 |
Family
ID=78389103
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202111033301.0A Active CN113627064B (en) | 2021-09-03 | 2021-09-03 | Roller kiln firing zone temperature prediction method based on mechanism and data mixed driving |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113627064B (en) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114234634B (en) * | 2022-01-13 | 2022-06-03 | 中南大学 | Roller kiln temperature control system based on heat map processing |
Citations (12)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107038307A (en) * | 2017-04-18 | 2017-08-11 | 中南大学 | Mechanism predicts integrated modelling approach with the Roller Conveying Kiln for Temperature that data are combined |
| CN108549732A (en) * | 2017-12-19 | 2018-09-18 | 中南大学 | Roller Conveying Kiln for Temperature soft-measuring modeling method based on local secondary Weighted Kernel principal component regression |
| CN111128312A (en) * | 2019-12-19 | 2020-05-08 | 湖南工业大学 | Zinc oxide volatilization kiln mixed modeling method based on mechanism and support vector machine |
| CN111665809A (en) * | 2020-06-16 | 2020-09-15 | 济南大学 | Modeling method and system for segmentation mechanism of rotary cement kiln |
| CN111931417A (en) * | 2020-07-21 | 2020-11-13 | 广东工业大学 | Method for analyzing and optimizing technological parameters of roller kiln |
| CN111950191A (en) * | 2020-07-07 | 2020-11-17 | 湖南大学 | Rotary kiln sintering temperature prediction method based on hybrid deep neural network |
| CN112380738A (en) * | 2020-10-20 | 2021-02-19 | 济南大学 | Rotary cement kiln combustion field reconstruction error compensation and optimization method, storage medium and system |
| CN112668196A (en) * | 2021-01-04 | 2021-04-16 | 西安理工大学 | Mechanism and data hybrid driven generation type countermeasure network soft measurement modeling method |
| CN112685927A (en) * | 2020-12-16 | 2021-04-20 | 广东工业大学 | Soft measurement method for temperature of workpiece in vacuum sintering furnace |
| CN112949132A (en) * | 2021-01-22 | 2021-06-11 | 浙江大学 | Method for acquiring internal temperature distribution of ceramic body |
| CN113177332A (en) * | 2021-03-09 | 2021-07-27 | 广东工业大学 | Rotary kiln sintering temperature prediction method based on combination of mechanism and data |
| CN113190974A (en) * | 2021-04-09 | 2021-07-30 | 广东工业大学 | Roller kiln temperature field multipoint prediction method based on deep learning |
-
2021
- 2021-09-03 CN CN202111033301.0A patent/CN113627064B/en active Active
Patent Citations (12)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107038307A (en) * | 2017-04-18 | 2017-08-11 | 中南大学 | Mechanism predicts integrated modelling approach with the Roller Conveying Kiln for Temperature that data are combined |
| CN108549732A (en) * | 2017-12-19 | 2018-09-18 | 中南大学 | Roller Conveying Kiln for Temperature soft-measuring modeling method based on local secondary Weighted Kernel principal component regression |
| CN111128312A (en) * | 2019-12-19 | 2020-05-08 | 湖南工业大学 | Zinc oxide volatilization kiln mixed modeling method based on mechanism and support vector machine |
| CN111665809A (en) * | 2020-06-16 | 2020-09-15 | 济南大学 | Modeling method and system for segmentation mechanism of rotary cement kiln |
| CN111950191A (en) * | 2020-07-07 | 2020-11-17 | 湖南大学 | Rotary kiln sintering temperature prediction method based on hybrid deep neural network |
| CN111931417A (en) * | 2020-07-21 | 2020-11-13 | 广东工业大学 | Method for analyzing and optimizing technological parameters of roller kiln |
| CN112380738A (en) * | 2020-10-20 | 2021-02-19 | 济南大学 | Rotary cement kiln combustion field reconstruction error compensation and optimization method, storage medium and system |
| CN112685927A (en) * | 2020-12-16 | 2021-04-20 | 广东工业大学 | Soft measurement method for temperature of workpiece in vacuum sintering furnace |
| CN112668196A (en) * | 2021-01-04 | 2021-04-16 | 西安理工大学 | Mechanism and data hybrid driven generation type countermeasure network soft measurement modeling method |
| CN112949132A (en) * | 2021-01-22 | 2021-06-11 | 浙江大学 | Method for acquiring internal temperature distribution of ceramic body |
| CN113177332A (en) * | 2021-03-09 | 2021-07-27 | 广东工业大学 | Rotary kiln sintering temperature prediction method based on combination of mechanism and data |
| CN113190974A (en) * | 2021-04-09 | 2021-07-30 | 广东工业大学 | Roller kiln temperature field multipoint prediction method based on deep learning |
Non-Patent Citations (1)
| Title |
|---|
| 融合机理与数据的篦冷机温度软测量模型;王美琪 等;《仪器仪表学报》;第185-191页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113627064A (en) | 2021-11-09 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106920006B (en) | Subway station air conditioning system energy consumption prediction method based on ISOA-LSSVM | |
| CN104765350B (en) | Cement decomposing furnace control method and system based on Combined model forecast control technology | |
| CN113627064B (en) | Roller kiln firing zone temperature prediction method based on mechanism and data mixed driving | |
| CN105241239A (en) | Intelligent optimal control method and device for sintered brick tunnel kiln roasting process | |
| CN106327004A (en) | Cement firing process optimizing method based on clinker quality index | |
| CN105911868A (en) | Multi-batch intermittent reactor two-dimension iterative learning feedback control method | |
| CN107038307A (en) | Mechanism predicts integrated modelling approach with the Roller Conveying Kiln for Temperature that data are combined | |
| Hao et al. | Multi-objective collaborative optimization in cement calcination process: A time domain rolling optimization method based on Jaya algorithm | |
| Ding et al. | Application of fuzzy immune algorithm and soft computing in the design of 2-DOF PID controller | |
| CN108776713B (en) | Method for decomposing and modeling temperature field area of grate bed of grate | |
| CN114334025A (en) | Construction and verification method of cement clinker calcination environment variable | |
| Ji et al. | Application of the improved the ELM algorithm for prediction of blast furnace gas utilization rate | |
| Zhang et al. | Building energy performance optimization: a new multi-objective particle swarm method | |
| CN109631607B (en) | Intelligent sintering ignition temperature control method considering gas pressure fluctuation | |
| Wei et al. | Multi-objective adaptive optimization model predictive control: Decreasing carbon emissions from a zinc oxide rotary kiln | |
| CN116085823A (en) | Boiler combustion control method and system | |
| CN116127760A (en) | Heating furnace temperature optimization method based on hybrid modeling and improved sand dune cat swarm algorithm | |
| CN113176797B (en) | Kiln temperature automatic control method for celadon biscuit firing process | |
| CN113488111B (en) | Ammonia injection amount optimization modeling method for SCR denitration system | |
| Ning et al. | Integrated modeling for roller kiln temperature prediction | |
| CN107066659A (en) | A kind of method that limit of utilization learning machine predicts cement decomposing furnace temperature | |
| Ito et al. | Multipoint-measurement multipoint-heating greenhouse temperature control with wooden pellet fuel using an adaptive model predictive control approach with a genetic algorithm | |
| Fan et al. | Ultrashort-term heat load prediction for heat exchange stations based on dual-source data feature fusion | |
| Qu et al. | Establishment and Optimization of Gas Flow Prediction Model for Annealing Furnace Based on GA-SVM | |
| Luo et al. | Temperature prediction of roller kiln based on mechanism and SSA-ELM data-driven integrated model |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |