Modeling method for velocity field of tuyere convolute area of blast furnace

Abstract

The invention provides a blast-furnace tuyere convolute zone velocity field modeling method, and relates to the technical field of complex industrial system.The blast-furnace tuyere convolute zone velocity field modeling method comprises the steps that on the basis that it is determined that a blast-furnace tuyere convolute zone is in a turbulent flow state, a two-dimensional volume fraction model, a two-dimensional momentum conservation equation and a mass conservation equation under turbulent flow are obtained through derivation on the basis of a Reynolds time-averaging method; meanwhile, a reformed group k-epsilon model is obtained, then a finite difference method is adopted to discretize the equation, iterative solution is conducted on the basis of the discretized equation to obtain the velocity field, the gas volume fraction, the turbulent flow kinetic energy and the turbulent flow dissipation rate, and through the solution method, the gas volume fraction and the turbulent flow kinetic energy are obtained. The calculation loss caused by excessive generalization in commercial fluid mechanics simulation software is avoided, the simulation efficiency under complex working conditions is improved, and transplantation and expansion of subsequent research are facilitated.

Description

Technical Field

[0001] The present invention relates to the technical field of complex industrial system modeling, in particular to a method for modeling the velocity field of a blast furnace tuyere raceway. Background Art

[0002] The blast furnace's tuyere raceway, the core region for energy transfer and reducing gas generation, requires precise characterization of its dynamic characteristics to understand its underlying operating mechanisms. However, traditional research methods have significant limitations, often relying on empirical formulas or simple mathematical models to roughly estimate parameters in the tuyere raceway. This approach fails to accurately reflect the complex physical and chemical processes occurring in this region. While some studies have employed simulation software, these approaches suffer from slow results and difficulty in intervention and transplantation. Summary of the Invention

[0003] In view of the shortcomings of the prior art, the present invention aims to propose a method for modeling the velocity field in the raceway of a blast furnace tuyere, comprising:

[0004] Step 1: Determine the modeling area boundary of the blast furnace tuyere raceway, divide the modeling area boundary into grids, and obtain the grid after the modeling area division;

[0005] Step 2: For the fluid in each grid in the grid after the modeling area is divided, the fluid in each grid is regarded as a microelement and a discrete gas volume fraction equation is constructed;

[0006] Step 3: Based on the gas volume fraction and the interphase force of the coke bed obstruction on the microelement, the discrete forms of the velocity component in the x-axis direction and the discrete forms of the velocity component in the y-axis direction are constructed;

[0007] Step 4: Based on the mass conservation equation under turbulence, construct a discrete form of the mass conservation equation;

[0008] Step 5: Based on the renormalization group k-ε model, construct the discrete forms of the turbulent kinetic energy k equation and the turbulent dissipation rate ε equation;

[0009] Step 6: For the grid after the modeling area is divided, the velocity field, pressure field, gas volume fraction, turbulent kinetic energy and turbulent dissipation rate are initialized. Based on the discrete gas volume fraction equation, the discrete form of the velocity component in the x-axis direction, the discrete form of the velocity component in the y-axis direction, the discrete form of the mass conservation equation, the discrete form of the turbulent kinetic energy k equation and the discrete form of the turbulent dissipation rate ε equation, the initialized velocity field, pressure field, gas volume fraction, turbulent kinetic energy and turbulent dissipation rate are iterated to obtain the final velocity field, final gas volume fraction, final turbulent kinetic energy and final turbulent dissipation rate.

[0010] Optionally, step 2 specifically includes:

[0011] Step 2.1: Consider the gas blown into the blast furnace tuyere raceway as an incompressible fluid and construct a volume fraction model, which can be expressed as:

[0012]

[0013] Among them, α g is the gas volume fraction, D is the diffusion coefficient, t is the time, u is the velocity component of the microelement in the x-axis direction, v is the velocity component of the microelement in the y-axis direction, and w is the velocity component of the microelement in the z-axis direction;

[0014] Step 2.2: Take the Reynolds-time average of both sides of the volume fraction model equation to obtain the two-dimensional volume fraction model under turbulence, which is expressed as:

[0015]

[0016] in, is the time-averaged value of the velocity component in the x-axis direction, is the time-averaged value of the velocity component in the y-axis direction;

[0017] Step 2.3: Use the finite difference method to discretize the two-dimensional volume fraction model under turbulence and obtain the discrete gas volume fraction equation, which is expressed as:

[0018]

[0019] in, is the gas volume fraction at the grid in the i-th row and j-th column of the grid after the modeling area is divided, and Δt is the time difference. is the velocity component in the x-axis direction at time step n at the grid in the i+1 / 2th row and jth column of the grid after the modeling area is divided, and Δy is the difference in the vertical direction. is the velocity component in the y-axis direction at the grid in the i-th row and j+1 / 2 column of the grid after the modeling area is divided, and Δx is the difference in the horizontal direction.

[0020] Optionally, step 3 specifically includes:

[0021] Step 3.1: Based on the gas volume fraction, construct the interphase force of the microelement due to the obstruction of the coke bed, which can be expressed as:

[0022]

[0023] Among them, F x1 is the interphase force in the x-axis direction after considering the coke volume fraction, F y1 is the interphase force in the y-axis direction after considering the coke volume fraction, F z1is the interphase force in the z-axis direction after considering the coke volume fraction, f1 and f2 are the source term coefficients;

[0024] Step 3.2: Construct the momentum conservation equation of the fluid based on the interphase force, which is expressed as:

[0025]

[0026] Among them, f x is the volume force density of the microelement in the x-axis direction, f y is the volume force density of the microelement in the y-axis direction, f z is the volume force density of the microelement in the z-axis direction, p is the pressure, μ is the dynamic viscosity coefficient, ρ g is the gas density;

[0027] Step 3.3: Take the Reynolds average of both sides of the momentum conservation equation for the fluid and obtain the two-dimensional momentum conservation equation for turbulence, which is expressed as:

[0028]

[0029] in, f x Time averaging processing, f y Time averaging processing, is the time-averaged pressure, u′ is the pulsating value of u, v′ is the pulsating value of v, is the time-averaged product of u′ and v′, F x1 Time averaging processing, F y1 The time-averaged processing of and Expressed as:

[0030]

[0031] Among them, C1 and C2 are the correction coefficients during time averaging processing;

[0032] Step 3.4: Use the finite difference method to discretize the two-dimensional momentum conservation equation in the x-axis direction under turbulence, and obtain the discrete form of the velocity component in the x-axis direction, which is expressed as:

[0033]

[0034] in, is the velocity component in the x-axis direction at the time step n+1 in the grid of the i+1 / 2th row and jth column in the modeling area. is the velocity component in the y-axis direction at the grid in the i+1 / 2th row and j+1 / 2th column of the grid after the modeling area is divided, It represents the pressure at time step n at the grid of row i+1 and column j in the grid after the modeling area is divided;

[0035] Among them, μ t is the turbulent viscosity coefficient, expressed as:

[0036]

[0037] Among them, C μ is an empirical constant, ε is the turbulent dissipation rate of turbulent kinetic energy, and k represents the turbulent kinetic energy;

[0038] Step 3.5: Use the finite difference method to discretize the two-dimensional momentum conservation equation in the y-axis direction under turbulence. The discrete form of the velocity component in the y-axis direction is obtained, which is expressed as:

[0039]

[0040] in, is the velocity component in the y-axis direction at the time step n+1 at the grid in the i-th row and j+1 / 2 column of the grid after the modeling area is divided, It is the velocity component in the x-axis direction at the grid in the i+1 / 2th row and j+1 / 2th column of the grid after the modeling area is divided.

[0041] Optionally, step 4 specifically includes:

[0042] Step 4.1: Obtain the mass conservation equation under turbulent flow, which is expressed as:

[0043]

[0044] Step 4.2: Use the finite difference method to discretize the mass conservation equation under turbulence and obtain the discrete form of the mass conservation equation, which is expressed as:

[0045]

[0046] in, is the velocity component in the x-axis direction at time step n in the grid of the i+1 / 2th row and jth column in the grid after the modeling area is divided, It is the velocity component in the y-axis direction at the grid in the i-th row and j+1 / 2 column of the grid after the modeling area is divided, at time step n.

[0047] Optionally, step 5 specifically includes:

[0048] The finite difference method is used to discretize the renormalization group k-ε model, and the discrete form of the turbulent kinetic energy k equation is obtained, which is expressed as formula (13). At the same time, the discrete form of the turbulent dissipation rate ε equation is obtained, which is expressed as formula (14). Formulas (13) and (14) are expressed as follows:

[0049]

[0050]

[0051] in, is the turbulent kinetic energy at the time step n+1 at the grid in the i-th row and j-th column of the grid after the modeling area is divided, is the turbulence dissipation rate at time step n+1 at the grid in the i-th row and j-th column of the grid after the modeling area is divided, is the dimensionless strain parameter at time step n at the grid of row i and column j in the grid after the modeling area is divided, η0 is the dimensionless strain threshold, It is expressed as the turbulent kinetic energy generation term at time step n at the grid in the i-th row and j-th column of the grid after the modeling area is divided.

[0052] Optionally, step 6 specifically includes:

[0053] Step 6.1: For the grid after the modeling area is divided, set the initial time step n = 0 and initialize the velocity field, that is, initialize u and v at the initial time step to obtain u n and v n Initialize the pressure field, gas volume fraction, turbulent kinetic energy and turbulent dissipation rate to obtain the pressure field p n , gas volume fraction Turbulent kinetic energy k n and the turbulent dissipation rate ε n ;

[0054] Step 6.2: Based on the turbulent kinetic energy k n and the turbulent dissipation rate ε n , calculate the turbulent viscosity μ t ;

[0055] Step 6.3: Set the correction value Δp. If the value of the correction value Δp is unknown, set μ t 、u n 、v n 、p n +Δp and Substitute the discrete form of the velocity component in the x-axis direction and the discrete form of the velocity component in the y-axis direction to obtain the initial velocity field u* and v* of time step n+1. Substitute the initial velocity field u* and v* into the discrete form of the mass conservation equation to calculate the value of the correction amount Δp. Substitute the value of the correction amount Δp into the initial velocity field u* and v* to obtain the velocity field u of time step n+1 n+1 and v n+1 ;

[0056] Step 6.4: Put u n+1 and v n+1 Substituting the discrete form of the turbulent kinetic energy k equation into the equation, we get the turbulent kinetic energy k at time step n+1: n+1 ,u n+1 and v n+1 Substituting the discrete form of the turbulence dissipation rate ε equation into the equation, we can obtain the turbulence dissipation rate ε at time step n+1: n+1 ;

[0057] Step 6.5: Put u n+1 、v n+1 and Substitute into the discrete gas volume fraction equation to obtain the gas volume fraction at time step n+1

[0058] Step 6.6: Calculate k n and k n+1 The difference between the two values is used as the residual of turbulent kinetic energy to calculate ε n and ε n+1 The difference between and The difference between the two values is used as the residual of gas volume fraction to calculate u n and u n+1 The difference is used as the velocity component residual in the x-axis direction to calculate v n and v n+1 The difference is taken as the residual of the velocity component in the y-axis direction;

[0059] Step 6.7: Determine whether the turbulent kinetic energy residual, turbulent dissipation rate residual, gas volume fraction residual, velocity component residual in the x-axis direction, and velocity component residual in the y-axis direction meet the preset conditions. If the turbulent kinetic energy residual, turbulent dissipation rate residual, gas volume fraction residual, velocity component residual in the x-axis direction, and velocity component residual in the y-axis direction meet the preset conditions, set k n+1 As the final turbulent kinetic energy, ε n+1 As the final turbulent dissipation rate, As the final gas volume fraction, u n+1 and v n+1as the final velocity field; if the turbulent kinetic energy residual, turbulent dissipation rate residual, gas volume fraction residual, velocity component residual in the x-axis direction, and velocity component residual in the y-axis direction do not meet the preset conditions, set n = n + 1 and return to step 6.2.

[0060] Optionally, the turbulent kinetic energy residual, turbulent dissipation rate residual, gas volume fraction residual, x-axis velocity component residual and y-axis velocity component residual described in step 6.7 meet the preset conditions, indicating that the turbulent kinetic energy residual is less than the preset threshold, the turbulent dissipation rate residual is less than the preset threshold, the gas volume fraction residual is less than the preset threshold, the x-axis velocity component residual is less than the preset threshold, and the y-axis velocity component residual is less than the preset threshold.

[0061] The beneficial effects of adopting the above technical solution are:

[0062] Based on the determination that the vortex zone of the blast furnace tuyere is in a turbulent state, the present invention derives a two-dimensional volume fraction model, a two-dimensional momentum conservation equation and a mass conservation equation under turbulence based on the Reynolds time-averaged method, and simultaneously obtains the renormalization group k-ε model, and then discretizes the above equations using the finite difference method. The velocity field, gas volume fraction, turbulent kinetic energy and turbulent dissipation rate are obtained by iterative solution based on the discretized equations. Through the above-mentioned solution method, the computational loss caused by "over-generalization" in commercial fluid mechanics simulation software is avoided, the simulation efficiency under complex working conditions is improved, and it is convenient for transplantation and expansion of subsequent research. BRIEF DESCRIPTION OF THE DRAWINGS

[0063] Figure 1 Schematic diagram of turbulent motion synthesized by superposition of laminar flow and pulsation in an embodiment of the present invention;

[0064] Figure 2 Schematic diagram of the flow of the velocity field modeling method of the blast furnace tuyere raceway in an embodiment of the present invention. DETAILED DESCRIPTION

[0065] The following embodiments of the present invention are described in further detail with reference to the accompanying drawings and examples. The following examples are used to illustrate the present invention but are not intended to limit the scope of the present invention.

[0066] The vortex zone of the blast furnace tuyere is a critical area within the blast furnace where combustion reactions occur and provide energy for the entire furnace. Establishing an appropriate physical model for this key area reflects the internal state of the blast furnace at a mechanistic level and is fundamental to reflecting the health of the blast furnace. This chapter describes its physical characteristics and models it from a fluid mechanics perspective. Fluid motion follows the laws of conservation of mass and momentum. First, a differential equation model for conservation of mass will be established. Then, based on this, a gas volume fraction model will be derived in combination with a two-phase flow model. Next, a momentum conservation model will be established for the fluid motion state, and the main forces acting on the fluid element will be explained. Finally, the momentum equation under the laminar flow form will be expanded to a turbulent form, making it more suitable for the fluid motion state in the blast furnace tuyere vortex zone.

[0067] From top to bottom, a blast furnace is divided into the throat, shaft, waist, bosh, and hearth. The tuyere (or whirlpool) is the area in front of the lower tuyere. Pulverized coal and hot air enter the tuyere through blast pipes, carrying kinetic energy that propels coke and iron oxides out of an ellipsoidal cavity. Here, the coke is oxidized to produce carbon dioxide, which provides energy for the blast furnace's smelting process. At high temperatures, the carbon dioxide and coke react to form carbon monoxide, which maintains the reduction of iron oxides.

[0068] When the flow rate is very low, the fluid flows in layers without mixing. As the flow rate gradually increases, the fluid streamlines begin to oscillate in a wave-like pattern. As the flow rate continues to increase, the streamlines become invisible and vortices appear, completely disrupting the laminar flow and transforming the fluid into a turbulent state. The state of the fluid can be determined by the Reynolds number (Re). When the Reynolds number is small, the flow is laminar, while when the Reynolds number is large, the flow is turbulent. In theory, turbulent states are difficult to solve due to the randomness of the particle motion parameters in time and space. In engineering applications, more attention is paid to the macroscopic effects caused by turbulent motion. In 1886, Reynolds decomposed turbulent motion into time-averaged motion and pulsating motion, and subsequently derived the Reynolds time-averaged equation. The emergence of the Reynolds time-averaged equation introduces Reynolds stress, which makes the time-averaged equation non-closed. To counteract the influence of Reynolds stress, various auxiliary models have emerged.

[0069] The present invention starts from the Navier-Stokes equation under laminar flow, explains the necessity of establishing a turbulence model for the blast furnace tuyere raceway, gradually obtains the Reynolds time-averaged equation, and adopts the RNG k-ε two-equation model to close the time-averaged equation.

[0070] Since the three-dimensional Navier-Stokes equations need to solve the velocity components (u, v, w) and pressure field (p) in three directions simultaneously and the turbulence model needs to solve the transport equations of k and ε, the number of grids increases exponentially with the dimension (the number of three-dimensional grids is on the order of N 3 , the number of two-dimensional grids is N 2), resulting in a significant increase in computational time and memory requirements. Considering the vertical cross-sectional symmetry of the blast furnace structure and tuyere location, this study subsequently simplified the blast furnace tuyere raceway modeling area into a two-dimensional model, encompassing both the radial direction from the tuyere to the furnace core and the axial direction from the furnace bottom to the furnace top.

[0071] First, perform a Reynolds number analysis:

[0072] The gas injection velocity in the vortex zone of the tuyere can reach 200 m / s, the depth of the vortex zone is about 1.5 m, and the temperature in the vortex zone can reach 2000 degrees Celsius. Combining these data, the Reynolds number of each gas component in the vortex zone of the tuyere can be calculated according to formula (15).

[0073]

[0074] Where u is the characteristic velocity of the fluid, l is the characteristic length, and μ is the dynamic viscosity of the fluid. Considering the reaction between the blown hot air and the coke, the main components of the coal gas formed are carbon monoxide, carbon dioxide, nitrogen, and hydrogen. When calculating the Reynolds number for these gas components, the following assumptions are made: To simplify the calculation, the pressure of the space in which the gas is located is taken as standard atmospheric pressure. The density formula for an ideal gas is analyzed as follows:

[0075]

[0076] Where P is the absolute pressure of the space in which the gas is located; M is the molar mass; R is the universal gas constant; and T is the thermodynamic temperature.

[0077] It can be found that the density of the gas decreases as the pressure decreases, and thus the gas Reynolds number decreases according to Equation (16). This shows that if the gas is in a turbulent state at standard atmospheric pressure and 2000 degrees Celsius, then the gas must be in a turbulent state under the high temperature environment of the blast furnace. The gas viscosity μ is mainly affected by temperature, and pressure has little effect on it. Use the Satterland formula to calculate the viscosity of the gas at 2000 degrees Celsius:

[0078]

[0079] Where μ0 is the gas viscosity at reference temperature T0 and standard atmospheric pressure, T is the target temperature at which the gas viscosity is required, and S is the Satterland coefficient, which is related to the type of gas. Taking the above into consideration, the Reynolds number for each gas component can be calculated as shown in Table 1:

[0080] Table 1 Reynolds numbers of main gas components

[0081]

[0082] When the Reynolds number Re < 2000, the fluid is considered to be in a laminar state, and when Re > 4000, the fluid is considered to be in a turbulent state. As shown in the table, the Reynolds number of each gas component in the fluid is much greater than 4000. Therefore, this section will improve on the previous section and establish a turbulent flow model for the fluid to better adapt to the production conditions in the blast furnace tuyere raceway.

[0083] Furthermore, the Reynolds time-averaged method is determined:

[0084] Turbulence contains vortices of various sizes. The size and rotation axis of the vortices are random. Large-scale vortices are mainly determined by boundary conditions, while small-scale vortices are mainly determined by viscous forces. Large-scale vortices are stretched and broken into small-scale vortices, and small-scale vortices are broken into even smaller-scale vortices. Large vortices mainly obtain energy from the mainstream and transmit it step by step through the interaction between vortices. Finally, due to the effect of viscosity, the small-scale vortices gradually disappear, and the mechanical energy is converted into heat energy. Reynolds regards turbulent motion as a motion composed of time-averaged values and pulsating values. Figure 1 As shown, from this point on, the mass conservation and momentum equations in laminar flow are converted to turbulent flow. The following shows the relationship between the instantaneous value and the time-averaged value and the pulsating value in the Reynolds time-averaged method.

[0085] Write the instantaneous variables u, v, w, and p in the laminar flow form as a combination of time-averaged values and pulsating values:

[0086]

[0087] in is the time-averaged value of the velocity component and pressure in the x and y directions of the microelement; u′, v′, and p′ are the corresponding pulsation values. Let s and t represent the aforementioned velocity components, respectively. The relationship between their time-averaged values and pulsation values is expressed as:

[0088]

[0089] Based on this, in order to solve the problems existing in the prior art, the present invention provides a method for modeling the velocity field of the blast furnace tuyere raceway. Figure 2 , which may include the following steps:

[0090] Step 1: Determine the modeling area boundary of the blast furnace tuyere raceway, divide the modeling area boundary into grids, and obtain the grid after the modeling area division;

[0091] Step 2: For the fluid in each grid in the grid after the modeling area is divided, the fluid in each grid is regarded as a microelement and a discrete gas volume fraction equation is constructed;

[0092] In the vortex zone of the blast furnace tuyere, substances are divided into two categories based on their dynamic behavior and physical scale differences: the gas phase of the blast hot air flow and the solid phase of the coke particle layer. The gas phase of the blast hot air flow that the gas is subjected to when flowing in the cavity is composed of high-speed blast air flow (containing coal powder) and combustion products (CO, CO2, etc.). Since the coal powder particle size is very small, its movement is dominated by gas turbulence and can be regarded as a continuous gas phase. The Euler-Euler model is used in the calculation, and the mixed phase of coal powder and gas is treated as a single continuous medium. Combined with the conservation of mass and taking into account the diffusion effect of gas, the gas volume fraction is expressed as:

[0093]

[0094] Step 2.1: To reduce the computational complexity, the gas blown into the blast furnace tuyere raceway is considered as an incompressible fluid and a volume fraction model is constructed, which is expressed as:

[0095]

[0096] Among them, α g is the gas volume fraction, D is the diffusion coefficient, t is the time, u is the velocity component of the microelement in the x-axis direction, v is the velocity component of the microelement in the y-axis direction, and w is the velocity component of the microelement in the z-axis direction;

[0097] Step 2.2: Perform Reynolds-time averaging on both sides of the volume fraction model equation. Specifically, substitute formula (18) into the volume fraction model, ignore the velocity and physical quantity changes in the z direction, and perform Reynolds-time averaging on both sides of the equation to obtain a two-dimensional volume fraction model under turbulence, which is expressed as:

[0098]

[0099] in, is the time-averaged value of the velocity component in the x-axis direction, is the time-averaged value of the velocity component in the y-axis direction;

[0100] Step 2.3: Use the finite difference method to discretize the two-dimensional volume fraction model under turbulence and obtain the discrete gas volume fraction equation, which is expressed as:

[0101]

[0102] in, is the gas volume fraction at the grid in the i-th row and j-th column of the grid after the modeling area is divided, and Δt is the time difference. is the velocity component in the x-axis direction at time step n at the grid in the i+1 / 2th row and jth column of the grid after the modeling area is divided, and Δy is the difference in the vertical direction. is the velocity component in the y-axis direction at the grid in the i-th row and j+1 / 2 column of the grid after the modeling area is divided, and Δx is the difference in the horizontal direction.

[0103] Step 3: Based on the gas volume fraction and the interphase force of the coke bed obstruction on the microelement, the discrete forms of the velocity component in the x-axis direction and the discrete forms of the velocity component in the y-axis direction are constructed;

[0104] The flow of fluid in the vortex zone of the blast furnace tuyere is hindered by the coke bed. Considering it as the source term of momentum conservation, it is expressed using the Ergen equation as follows:

[0105]

[0106] Among them, F x 、F y 、F z is the interphase force, which represents the resistance of the coke bed to the fluid flow, and its direction is opposite to the flow velocity. f1 and f2 are source term coefficients:

[0107]

[0108] Where η g is the gas viscosity; ε s is the porosity of the charge; ψ s is the shape factor of the solid particles; d s is the diameter of the solid particles.

[0109] Taking into account the obvious differences in gas volume fractions at different positions in the blast furnace tuyere raceway, when the position is close to the tuyere, the gas volume fraction is larger and the interphase force it is subjected to is smaller. When it is close to the edge of the tuyere raceway, the solid volume fraction increases and the interphase force is enhanced. Therefore, it is necessary to take the gas volume fraction into account in the interphase force.

[0110] Step 3.1: Based on the gas volume fraction, construct the interphase force of the microelement due to the obstruction of the coke bed, which can be expressed as:

[0111]

[0112] Among them, F x1 is the interphase force in the x-axis direction after considering the coke volume fraction, F y1 is the interphase force in the y-axis direction after considering the coke volume fraction, F z1 is the interphase force in the z-axis direction after considering the coke volume fraction, f1 and f2 are the source term coefficients;

[0113] Step 3.2: Construct the momentum conservation equation of the fluid based on the interphase force, which is expressed as:

[0114]

[0115] Among them, f x is the volume force density of the microelement in the x-axis direction, f y is the volume force density of the microelement in the y-axis direction, f z is the volume force density of the microelement in the z-axis direction, p is the pressure, μ is the dynamic viscosity coefficient, ρ g is the gas density;

[0116] Step 3.3: Perform Reynolds-averaged calculations on both sides of the momentum conservation equation. Specifically, substitute formula (18) into the momentum conservation equation of the fluid, ignore the velocity and physical quantity changes in the z direction, and perform Reynolds-averaged calculations on both sides of the equation to obtain the two-dimensional momentum conservation equation under turbulence, which is expressed as:

[0117]

[0118] in, f x Time averaging processing, f y Time averaging processing, is the time-averaged pressure, u′ is the pulsating value of u, v′ is the pulsating value of v, is the time-averaged product of u′ and v′, F x1 Time averaging processing, F y1 The time-averaged processing of and Expressed as:

[0119]

[0120] Among them, C1 and C2 are the correction coefficients during time averaging processing.

[0121] Step 3.4: Use the finite difference method to discretize the two-dimensional momentum conservation equation in the x-axis direction under turbulence, and obtain the discrete form of the velocity component in the x-axis direction, which is expressed as:

[0122]

[0123] in, is the velocity component in the x-axis direction at the time step n+1 in the grid of the i+1 / 2th row and jth column in the modeling area. is the velocity component in the y-axis direction at the grid in the i+1 / 2th row and j+1 / 2th column of the grid after the modeling area is divided, It represents the pressure at time step n at the grid of row i+1 and column j in the grid after the modeling area is divided;

[0124] Among them, μ t is the turbulent viscosity coefficient, expressed as:

[0125]

[0126] Among them, C μ is an empirical constant, ε is the turbulent dissipation rate of turbulent kinetic energy, and k represents the turbulent kinetic energy;

[0127] Step 3.5: Use the finite difference method to discretize the two-dimensional momentum conservation equation in the y-axis direction under turbulence. The discrete form of the velocity component in the y-axis direction is obtained, which is expressed as:

[0128]

[0129] in, is the velocity component in the y-axis direction at the time step n+1 at the grid in the i-th row and j+1 / 2 column of the grid after the modeling area is divided, It is the velocity component in the x-axis direction at the grid in the i+1 / 2th row and j+1 / 2th column of the grid after the modeling area is divided.

[0130] Step 4: Based on the mass conservation equation under turbulence, construct a discrete form of the mass conservation equation;

[0131] The mass conservation equation is expressed as:

[0132]

[0133] Substituting formula (18) into formula (23), we obtain the two-dimensional mass conservation equation under turbulence:

[0134]

[0135] Combining (23) with (19) we can get:

[0136]

[0137] in, They represent the time-averaged values of the pulsating terms of the velocity components of the microelement in the x and y directions respectively.

[0138] From formula (19), we know that the time-averaged result of the velocity pulsation value is 0, so we can get:

[0139]

[0140] Step 4.1: Obtain the mass conservation equation under turbulence. Specifically, combine Equation (24) and Equation (25) to obtain the mass conservation equation under turbulence, which is expressed as:

[0141]

[0142] Step 4.2: Use the finite difference method to discretize the mass conservation equation under turbulence and obtain the discrete form of the mass conservation equation, which is expressed as:

[0143]

[0144] in, is the velocity component in the x-axis direction at time step n in the grid of the i+1 / 2th row and jth column in the grid after the modeling area is divided, It is the velocity component in the y-axis direction at the grid in the i-th row and j+1 / 2 column of the grid after the modeling area is divided, at time step n.

[0145] Step 5: Based on the renormalization group k-ε model, construct the discrete forms of the turbulent kinetic energy k equation and the turbulent dissipation rate ε equation;

[0146] Among them, the renormalization group k-ε model is expressed as:

[0147]

[0148] Among them, r is the free index, r = 1, 2, U r Including U1 or U2, U1 is the velocity component of the microelement in the x-axis direction, that is, u, U2 is the velocity component of the microelement in the y-axis direction, that is, v; s is a dummy index, s = 1, 2, X s Including X1 or X2, X1 is the coordinate of the microelement in the x-axis direction, X2 is the coordinate of the microelement in the y-axis direction, α k The diffusion Prandtl number, μ, represents the turbulent kinetic energy. eff The effective viscosity is expressed as molecular viscosity μ and turbulent viscosity μ t The sum of G k represents the turbulent kinetic energy generation term, α ε is the diffusion Prandtl number of the dissipation rate, C * 1ε 、C 1ε and C 2ε is the corrected empirical coefficient. More detailed physical quantity expressions and empirical values are shown below.

[0149]

[0150] Among them, E rsis the strain rate tensor; η is the dimensionless strain parameter; η0 is the dimensionless strain threshold; β is a damping coefficient that controls the sensitivity of the correction term to high strain rate flow.

[0151] The finite difference method is used to discretize the renormalization group k-ε model, and the discrete form of the turbulent kinetic energy k equation is obtained, which is expressed as formula (13). At the same time, the discrete form of the turbulent dissipation rate ε equation is obtained, which is expressed as formula (14). Formulas (13) and (14) are expressed as follows:

[0152]

[0153] in, is the turbulent kinetic energy at the time step n+1 at the grid in the i-th row and j-th column of the grid after the modeling area is divided, is the turbulence dissipation rate at time step n+1 at the grid in the i-th row and j-th column of the grid after the modeling area is divided, is the dimensionless strain parameter at time step n at the grid of row i and column j in the grid after the modeling area is divided, η0 is the dimensionless strain threshold, It is expressed as the turbulent kinetic energy generation term at time step n at the grid in the i-th row and j-th column of the grid after the modeling area is divided.

[0154] Step 6: For the grid after the modeling area is divided, the velocity field, pressure field, gas volume fraction, turbulent kinetic energy and turbulent dissipation rate are initialized. Based on the discrete gas volume fraction equation, the discrete form of the velocity component in the x-axis direction, the discrete form of the velocity component in the y-axis direction, the discrete form of the mass conservation equation, the discrete form of the turbulent kinetic energy k equation and the discrete form of the turbulent dissipation rate ε equation, the initialized velocity field, pressure field, gas volume fraction, turbulent kinetic energy and turbulent dissipation rate are iterated to obtain the final velocity field, final gas volume fraction, final turbulent kinetic energy and final turbulent dissipation rate.

[0155] Step 6.1: For the grid after the modeling area is divided, set the initial time step n = 0 and initialize the velocity field, that is, initialize u and v at the initial time step to obtain u n and v n Initialize the pressure field, gas volume fraction, turbulent kinetic energy and turbulent dissipation rate to obtain the pressure field p n , gas volume fraction Turbulent kinetic energy k n and the turbulent dissipation rate ε n ;

[0156] Step 6.2: Based on the turbulent kinetic energy k n and the turbulent dissipation rate ε n , calculate the turbulent viscosity μ t ;

[0157] Step 6.3: Set the correction value Δp. If the value of the correction value Δp is unknown, set μ t 、u n 、v n 、p n +Δp and Substitute the discrete form of the velocity component in the x-axis direction and the discrete form of the velocity component in the y-axis direction to obtain the initial velocity field u* and v* of time step n+1. Substitute the initial velocity field u* and v* into the discrete form of the mass conservation equation to calculate the value of the correction amount Δp. Substitute the value of the correction amount Δp into the initial velocity field u* and v* to obtain the velocity field u of time step n+1 n+1 and v n+1 ;

[0158] Step 6.4: Put u n+1 and v n+1 Substituting the discrete form of the turbulent kinetic energy k equation into the equation, we get the turbulent kinetic energy k at time step n+1: n+1 ,u n+1 and v n+1 Substituting the discrete form of the turbulence dissipation rate ε equation into the equation, we can obtain the turbulence dissipation rate ε at time step n+1: n+1 ;

[0159] Step 6.5: Put u n+1 、v n+1 and Substitute into the discrete gas volume fraction equation to obtain the gas volume fraction at time step n+1

[0160] Step 6.6: Calculate k n and k n+1 The difference between the two values is used as the residual of turbulent kinetic energy to calculate ε n and ε n+1 The difference between and The difference between the two values is used as the residual of gas volume fraction to calculate u n and u n+1 The difference is used as the velocity component residual in the x-axis direction to calculate v n and v n+1 The difference is taken as the residual of the velocity component in the y-axis direction;

[0161] Step 6.7: Determine whether the turbulent kinetic energy residual, turbulent dissipation rate residual, gas volume fraction residual, velocity component residual in the x-axis direction, and velocity component residual in the y-axis direction meet the preset conditions. If the turbulent kinetic energy residual, turbulent dissipation rate residual, gas volume fraction residual, velocity component residual in the x-axis direction, and velocity component residual in the y-axis direction meet the preset conditions, set k n+1As the final turbulent kinetic energy, ε n+1 As the final turbulent dissipation rate, As the final gas volume fraction, u n+1 and v n+1 as the final velocity field; if the turbulent kinetic energy residual, turbulent dissipation rate residual, gas volume fraction residual, velocity component residual in the x-axis direction, and velocity component residual in the y-axis direction do not meet the preset conditions, set n = n + 1 and return to step 6.2.

[0162] Among them, the turbulent kinetic energy residual, the turbulent dissipation rate residual, the gas volume fraction residual, the velocity component residual in the x-axis direction, and the velocity component residual in the y-axis direction meet the preset conditions, indicating that the turbulent kinetic energy residual is less than the preset threshold, the turbulent dissipation rate residual is less than the preset threshold, the gas volume fraction residual is less than the preset threshold, the velocity component residual in the x-axis direction is less than the preset threshold, and the velocity component residual in the y-axis direction is less than the preset threshold.

[0163] This paper examines the key physical properties of the blast furnace tuyere raceway and their impact on the smelting process. It constructs a coupled mathematical model based on the conservation of mass and momentum, and systematically explains the correction method for the governing equations under turbulent flow conditions. This foundation, combined with auxiliary differential equations and a pressure-velocity correction algorithm, provides a complete solution process.

[0164] The present invention first establishes a set of governing equations for mass-momentum conservation in the tuyere vortex zone. By introducing a gas volume fraction model, the momentum transfer effect between the gas and solid phases is quantified, providing theoretical support for the correction of the interphase force term. In the derivation of the momentum conservation equations, differential control volume analysis combined with Newton's second law is employed to systematically characterize the dynamic equilibrium relationships between the pressure gradient force, viscous stress, and other factors acting on the fluid element.

[0165] Since the blast furnace tuyere vortex zone is in a high temperature and high pressure environment, and the inlet airflow has a high speed characteristic (Re>1×10 4 ), this study uses the Reynolds time-averaged method (RANS) to transform the laminar flow control equations into conservation equations suitable for turbulent conditions. By introducing the RNG k-ε auxiliary differential equation model closed equation system, the problem of the equation not being closed caused by the Reynolds stress term in the time-averaged process is effectively solved.

[0166] In terms of numerical solution, a discretization strategy based on the finite difference method is proposed, and the RNG k-ε auxiliary differential equation is coupled with the SIMPLE algorithm to realize the iterative correction of the velocity-pressure field.

[0167] The above description is merely a preferred embodiment of the present disclosure and an explanation of the technical principles employed. 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 a specific combination of the above-mentioned technical features, but should also encompass other technical solutions formed by any combination of the above-mentioned technical features or their equivalents without departing from the above-mentioned inventive concept. For example, a technical solution formed by mutually replacing the above-mentioned features with (but not limited to) technical features with similar functions disclosed in the embodiments of the present disclosure.

Claims

1. A velocity field modeling method for a blast furnace tuyere raceway, characterized in that: include: Step 1: Determine the modeling area boundary of the blast furnace tuyere raceway, divide the modeling area boundary into grids, and obtain the grid after the modeling area division; Step 2: For the fluid in each grid in the grid after the modeling area is divided, the fluid in each grid is regarded as a microelement and a discrete gas volume fraction equation is constructed; Step 3: Based on the gas volume fraction and the interphase force of the coke bed obstruction on the microelement, the discrete forms of the velocity component in the x-axis direction and the discrete forms of the velocity component in the y-axis direction are constructed; Step 4: Based on the mass conservation equation under turbulence, construct a discrete form of the mass conservation equation; Step 5: Based on the renormalization group k-ε model, construct the discrete forms of the turbulent kinetic energy k equation and the turbulent dissipation rate ε equation; Step 6: For the grid after the modeling area is divided, the velocity field, pressure field, gas volume fraction, turbulent kinetic energy and turbulent dissipation rate are initialized. Based on the discrete gas volume fraction equation, the discrete form of the velocity component in the x-axis direction, the discrete form of the velocity component in the y-axis direction, the discrete form of the mass conservation equation, the discrete form of the turbulent kinetic energy k equation and the discrete form of the turbulent dissipation rate ε equation, the initialized velocity field, pressure field, gas volume fraction, turbulent kinetic energy and turbulent dissipation rate are iterated to obtain the final velocity field, final gas volume fraction, final turbulent kinetic energy and final turbulent dissipation rate.

2. The blast furnace tuyere raceway velocity field modeling method according to claim 1, characterized in that: Step 2 specifically includes: Step 2.1: Consider the gas blown into the blast furnace tuyere raceway as an incompressible fluid and construct a volume fraction model, which can be expressed as: Among them, α g is the gas volume fraction, D is the diffusion coefficient, t is the time, u is the velocity component of the microelement in the x-axis direction, v is the velocity component of the microelement in the y-axis direction, and w is the velocity component of the microelement in the z-axis direction; Step 2.2: Take the Reynolds-time average of both sides of the volume fraction model equation to obtain the two-dimensional volume fraction model under turbulence, which is expressed as: in, is the time-averaged value of the velocity component in the x-axis direction, is the time-averaged value of the velocity component in the y-axis direction; Step 2.3: Use the finite difference method to discretize the two-dimensional volume fraction model under turbulence and obtain the discrete gas volume fraction equation, which is expressed as: in, is the gas volume fraction at the grid in the i-th row and j-th column of the grid after the modeling area is divided, and Δt is the time difference. is the velocity component in the x-axis direction at time step n at the grid in the i+1 / 2th row and jth column of the grid after the modeling area is divided, and Δy is the difference in the vertical direction. is the velocity component in the y-axis direction at the grid in the i-th row and j+1 / 2 column of the grid after the modeling area is divided, and Δx is the difference in the horizontal direction.

3. The blast furnace tuyere raceway velocity field modeling method according to claim 2, characterized in that: Step 3 specifically includes: Step 3.1: Based on the gas volume fraction, construct the interphase force of the microelement due to the obstruction of the coke bed, which can be expressed as: Among them, F x1 is the interphase force in the x-axis direction after considering the coke volume fraction, F y1 is the interphase force in the y-axis direction after considering the coke volume fraction, F z1 is the interphase force in the z-axis direction after considering the coke volume fraction, f1 and f2 are the source term coefficients; Step 3.2: Construct the momentum conservation equation of the fluid based on the interphase force, which is expressed as: Among them, f x is the volume force density of the microelement in the x-axis direction, f y is the volume force density of the microelement in the y-axis direction, f z is the volume force density of the microelement in the z-axis direction, p is the pressure, μ is the dynamic viscosity coefficient, ρ g is the gas density; Step 3.3: Take the Reynolds average of both sides of the momentum conservation equation for the fluid and obtain the two-dimensional momentum conservation equation for turbulence, which is expressed as: in, f x Time averaging processing, f y Time averaging processing, is the time-averaged pressure, u′ is the pulsating value of u, v′ is the pulsating value of v, is the time-averaged product of u′ and v′, F x1 Time averaging processing, F y1 The time-averaged processing of and Expressed as: Among them, C1 and C2 are the correction coefficients during time averaging processing; Step 3.4: Use the finite difference method to discretize the two-dimensional momentum conservation equation in the x-axis direction under turbulence, and obtain the discrete form of the velocity component in the x-axis direction, which is expressed as: in, is the velocity component in the x-axis direction at the time step n+1 in the grid of the i+1 / 2th row and jth column in the modeling area. is the velocity component in the y-axis direction at the grid in the i+1 / 2th row and j+1 / 2th column of the grid after the modeling area is divided, It represents the pressure at time step n at the grid of row i+1 and column j in the grid after the modeling area is divided; Among them, μ t is the turbulent viscosity coefficient, expressed as: Among them, C μ is an empirical constant, ε is the turbulent dissipation rate of turbulent kinetic energy, and k represents the turbulent kinetic energy; Step 3.5: Use the finite difference method to discretize the two-dimensional momentum conservation equation in the y-axis direction under turbulence. The discrete form of the velocity component in the y-axis direction is obtained, which is expressed as: in, is the velocity component in the y-axis direction at the time step n+1 at the grid in the i-th row and j+1 / 2 column of the grid after the modeling area is divided, It is the velocity component in the x-axis direction at the grid in the i+1 / 2th row and j+1 / 2th column of the grid after the modeling area is divided.

4. The blast furnace tuyere raceway velocity field modeling method according to claim 3, characterized in that: Step 4 specifically includes: Step 4.1: Obtain the mass conservation equation under turbulent flow, which is expressed as: Step 4.2: Use the finite difference method to discretize the mass conservation equation under turbulence and obtain the discrete form of the mass conservation equation, which is expressed as: in, is the velocity component in the x-axis direction at time step n in the grid of the i+1 / 2th row and jth column in the grid after the modeling area is divided, It is the velocity component in the y-axis direction at the grid in the i-th row and j+1 / 2 column of the grid after the modeling area is divided, at time step n.

5. The blast furnace tuyere raceway velocity field modeling method according to claim 4, characterized in that: Step 5 specifically includes: The finite difference method is used to discretize the renormalization group k-ε model, and the discrete form of the turbulent kinetic energy k equation is obtained, which is expressed as formula (13). At the same time, the discrete form of the turbulent dissipation rate ε equation is obtained, which is expressed as formula (14). Formulas (13) and (14) are expressed as follows: in, is the turbulent kinetic energy at the time step n+1 at the grid in the i-th row and j-th column of the grid after the modeling area is divided, is the turbulence dissipation rate at time step n+1 at the grid in the i-th row and j-th column of the grid after the modeling area is divided, is the dimensionless strain parameter at time step n at the grid of row i and column j in the grid after the modeling area is divided, η0 is the dimensionless strain threshold, It is expressed as the turbulent kinetic energy generation term at time step n at the grid in the i-th row and j-th column of the grid after the modeling area is divided.

6. The blast furnace tuyere raceway velocity field modeling method according to claim 5, characterized in that: Step 6 specifically includes: Step 6.1: For the grid after the modeling area is divided, set the initial time step n = 0 and initialize the velocity field, that is, initialize u and v at the initial time step to obtain u n and v n Initialize the pressure field, gas volume fraction, turbulent kinetic energy and turbulent dissipation rate to obtain the pressure field p n , gas volume fraction Turbulent kinetic energy k n and the turbulent dissipation rate ε n ; Step 6.2: Based on the turbulent kinetic energy k n and the turbulent dissipation rate ε n , calculate the turbulent viscosity μ t ; Step 6.3: Set the correction value Δp. If the value of the correction value Δp is unknown, set μ t 、u n 、v n 、p n +Δp and Substitute the discrete form of the velocity component in the x-axis direction and the discrete form of the velocity component in the y-axis direction to obtain the initial velocity field u* and v* of time step n+1. Substitute the initial velocity field u* and v* into the discrete form of the mass conservation equation to calculate the value of the correction amount Δp. Substitute the value of the correction amount Δp into the initial velocity field u* and v* to obtain the velocity field u of time step n+1 n+1 and v n+1 ; Step 6.4: Put u n+1 and v n+1 Substituting the discrete form of the turbulent kinetic energy k equation into the equation, we get the turbulent kinetic energy k at time step n+1: n+1 ,u n+1 and v n+1 Substituting the discrete form of the turbulence dissipation rate ε equation into the equation, we can obtain the turbulence dissipation rate ε at time step n+1: n+1 ; Step 6.5: Put u n+1 、v n+1 and Substitute into the discrete gas volume fraction equation to obtain the gas volume fraction at time step n+1 Step 6.6: Calculate k n and k n+1 The difference between the two values is used as the residual of turbulent kinetic energy to calculate ε n and ε n+1 The difference between and The difference between the two values is used as the gas volume fraction residual to calculate u n and u n+1 The difference is used as the velocity component residual in the x-axis direction to calculate v n and v n+1 The difference is taken as the residual of the velocity component in the y-axis direction; Step 6.7: Determine whether the turbulent kinetic energy residual, turbulent dissipation rate residual, gas volume fraction residual, velocity component residual in the x-axis direction, and velocity component residual in the y-axis direction meet the preset conditions. If the turbulent kinetic energy residual, turbulent dissipation rate residual, gas volume fraction residual, velocity component residual in the x-axis direction, and velocity component residual in the y-axis direction meet the preset conditions, set k n+1 As the final turbulent kinetic energy, ε n+1 As the final turbulent dissipation rate, As the final gas volume fraction, u n+1 and v n+1 as the final velocity field; if the turbulent kinetic energy residual, turbulent dissipation rate residual, gas volume fraction residual, velocity component residual in the x-axis direction, and velocity component residual in the y-axis direction do not meet the preset conditions, set n = n + 1 and return to step 6.

2.

7. The blast furnace tuyere raceway velocity field modeling method according to claim 5, characterized in that: The turbulent kinetic energy residual, turbulent dissipation rate residual, gas volume fraction residual, x-axis velocity component residual and y-axis velocity component residual described in step 6.7 meet the preset conditions, indicating that the turbulent kinetic energy residual is less than the preset threshold, the turbulent dissipation rate residual is less than the preset threshold, the gas volume fraction residual is less than the preset threshold, the x-axis velocity component residual is less than the preset threshold, and the y-axis velocity component residual is less than the preset threshold.

Images