Blast furnace tuyere rotation area temperature field calculation method based on gas-solid two-phase flow

Abstract

The invention provides a blast furnace tuyere convolute area temperature field calculation method based on gas-solid two-phase flow, and relates to the technical field of complex industrial system modeling, and the blast furnace tuyere convolute area temperature field calculation method comprises the steps that an infrared camera is used for reading temperature of a tuyere convolute area to replace blast temperature to serve as model inlet temperature input; two-phase flow is adopted to conduct modeling research on a blast-furnace tuyere convolute area, and the influence of mass and heat transfer between gas and solid on multiple physical fields is considered; a discrete control equation of a gas-phase volume fraction field, a gas-phase velocity field and a gas-phase temperature field based on finite difference is constructed, meanwhile, high-order terms of viscous stress and Reynolds stress are ignored, and a regular grid is constructed, so that a discrete control equation of the gas-phase volume fraction field, the gas-phase velocity field and the gas-phase temperature field is constructed. The complexity and iteration cost of finite element and finite volume discrete formats are avoided, and the iteration speed is increased.

Description

Technical Field

[0001] The present invention relates to the technical field of complex industrial system modeling, and in particular to a method for calculating the temperature field in a blast furnace tuyere raceway based on gas-solid two-phase flow. Background Art

[0002] Traditional methods for studying the vortex zone of a tuyere are primarily based on single-fluid assumptions or simplified empirical models, ignoring gas-solid coupling effects. However, actual operating conditions involve complex gas-solid interactions. To accurately reveal the movement mechanisms and heat and mass transfer patterns of multiphase materials, a model that considers gas-solid coupling is necessary.

[0003] Currently, most models use commercial simulation software such as Fluent to perform iterative calculations using discrete methods. The calculation steps are complex and the amount of calculation is large, the calculation is slow, and there is a large lag. Summary of the Invention

[0004] In view of the shortcomings of the prior art, the present invention aims to propose a method for calculating the temperature field in the blast furnace tuyere raceway based on gas-solid two-phase flow, comprising:

[0005] Step 1: Establish the modeling area boundary of the blast furnace tuyere raceway, and generate the target mesh of the modeling area based on the modeling area boundary through mesh generation technology;

[0006] Step 2: For each grid in the target grid, construct the governing equation for the discretized gas phase volume fraction based on the gas-solid two-phase flow;

[0007] Step 3: For each grid in the target grid, construct an expression for the discretized velocity component based on the gas-solid two-phase flow;

[0008] Step 4: Construct the turbulent kinetic energy equation of the airflow phase, and then process the turbulent kinetic energy equation of the airflow phase to obtain a discrete turbulent kinetic energy equation;

[0009] Step 5: Construct the airflow phase turbulence dissipation rate equation, process the airflow phase turbulence dissipation rate equation, and obtain the discrete airflow phase turbulence dissipation rate equation;

[0010] Step 6: Construct a discrete temperature field mathematical model by combining fluid flow, chemical reaction exothermicity, chemical reaction endothermicity, chemical reaction heat transfer, and phase change processes;

[0011] Step 7: Based on the control equation of the discretized gas phase volume fraction, the expression of the discretized velocity component, the discrete turbulent kinetic energy equation, and the discrete gas phase turbulent dissipation rate equation, iteratively solve to obtain the final volume fraction field and the final velocity field of the blast furnace tuyere raceway;

[0012] Step 8: Based on the final velocity field, the discrete temperature field mathematical model is iteratively solved to obtain the final temperature field of the blast furnace tuyere raceway.

[0013] Optionally, step 2 specifically includes:

[0014] Step 2.1: For each grid in the target grid, construct a control equation for the volume fraction of the gas phase based on the gas-solid two-phase flow. The control equation for the volume fraction of the gas phase is expressed as:

[0015]

[0016] Among them, α g is the volume fraction of the gas phase, which represents the volume fraction of the gas in a grid, t represents time, and u x is the air flow phase velocity The velocity component in the horizontal direction x, u y is the air flow phase velocity The velocity component in the vertical direction y;

[0017] Step 2.2: Discretize the governing equation of the airflow phase volume fraction to obtain the governing equation of the discretized airflow phase volume fraction, which is expressed as:

[0018]

[0019] Among them, α g n+1 [i,j] is the volume fraction of the airflow phase on the grid located in the i-th row and the j-th column in the n+1th iteration, Δt is the time difference, is the horizontal velocity component on the grid at row i and column j in the nth iteration, Δx is the horizontal difference, and Δy is the vertical difference.

[0020] Optionally, step 3 specifically includes:

[0021] Step 3.1: For each grid in the target grid, the momentum conservation equation is adjusted based on the gas-solid two-phase flow to obtain an adjusted momentum conservation equation, which is expressed as:

[0022]

[0023] Among them, ρ g is the air flow phase density, is the viscous stress, is the Reynolds stress, p is the pressure, g is the acceleration of gravity, F g is the interphase force, which includes the drag force F drag , and the lift force F liftand virtual body force F vm ;

[0024] Step 3.2: Based on viscous stress Formula, Reynolds stress Formula, drag force F drag Formula, lift F lift Formula and virtual body force F vm Formula (3) is expanded and adjusted to obtain the expression of velocity component:

[0025]

[0026] Where μ is the molecular viscosity of the fluid, μ g,t is the turbulent kinematic viscosity of the air phase, F x is the interphase force F g The horizontal component, F y is the interphase force F g The vertical component of

[0027] Step 3.3: Discretize the expression of the velocity component to obtain the expression of the discretized velocity component, which is expressed as:

[0028]

[0029] Among them, u x n+1 [i,j] is the horizontal velocity component on the grid at row i and column j in the n+1th iteration. is the vertical velocity component on the grid at row i+1 and column j in the nth iteration; p n [i+1,j] is the pressure on the grid at row i+1 and column j in the nth iteration, F x [i,j] is the F on the grid located in row i and column j x , F y [i,j] is the F on the grid located in row i and column j y .

[0030] Optionally, step 4 specifically includes:

[0031] Step 4.1: Construct the turbulent kinetic energy equation of the air flow phase, which is expressed as:

[0032]

[0033] Among them, k is the turbulent kinetic energy of the air flow phase, ∈ is the turbulent dissipation rate of the air flow phase, μ k ' is the intermediate variable, G k is the generation term of turbulent kinetic energy, μ k' is expressed as:

[0034]

[0035]

[0036] Among them, μ t is the turbulent viscosity, σ k is the turbulent kinetic energy Prandtl number, C μ is a constant;

[0037] Among them, G k Expressed as:

[0038] G k =μ t S i,j S i,j (11)

[0039]

[0040] Among them, S i,j is the intermediate variable, u i represents the i-th component of velocity u, u j represents the jth component of velocity u, x i It represents the i-th component of the coordinate basis vector, x j It represents the jth component of the coordinate basis vector;

[0041] Step 4.2: Expand the turbulent kinetic energy equation of the airflow phase to obtain the expanded turbulent kinetic energy equation, which is expressed as:

[0042]

[0043] Step 4.3: Discretize the expanded turbulent kinetic energy equation to obtain the discrete turbulent kinetic energy equation, which is expressed as:

[0044]

[0045] Among them, k n+1 [i,j] is the turbulent kinetic energy of the airflow phase on the grid located in the i-th row and the j-th column in the n+1th iteration, μ k ' n [i,j] is μ on the grid at row i and column j in the nth iteration. k '; is the generated term of turbulent kinetic energy on the grid at row i and column j in the nth iteration, ∈ n [i, j] is the turbulent dissipation rate of the airflow phase on the grid located in the i-th row and j-th column in the n-th iteration.

[0046] Optionally, step 5 specifically includes:

[0047] Step 5.1: Construct the turbulent dissipation rate equation for the airflow phase, expressed as:

[0048]

[0049] Among them, σ ∈ is the dissipation rate Prandtl number, C ∈1 and C ∈2 is the value calibrated by experiment;

[0050] Step 5.2: Expand the airflow phase turbulence dissipation rate equation to obtain the expanded airflow phase turbulence dissipation rate equation, which is expressed as:

[0051]

[0052] Step 5.3: Discretize the expanded airflow phase turbulence dissipation rate equation to obtain the discrete airflow phase turbulence dissipation rate equation, which is expressed as:

[0053]

[0054] Among them, ∈ n+1 [i,j] is the turbulent dissipation rate of the airflow phase on the grid located in the i-th row and the j-th column in the n+1th iteration, is the turbulent viscosity on the grid at row i and column j in the nth iteration; is the turbulent viscosity on the grid at row i+1 and column j in the nth iteration.

[0055] Optionally, step 6 specifically includes:

[0056] Step 6.1: Combine fluid flow, chemical reaction exothermicity, chemical reaction endothermicity, chemical reaction heat transfer, and phase change processes to construct the energy conservation equation, which can be expressed as:

[0057]

[0058] Among them, E g is the energy of the air flow phase, k g is the thermal conductivity of the gas, Q1 is the heat generated by the chemical reaction, Q2 is the heat transferred between the gas and the solid, T g is the air flow phase temperature;

[0059] Among them, the air flow phase energy E g and air flow temperature T g The relationship between them is expressed as:

[0060] E g =ρ g C p Tg (19)

[0061] Among them, C p is the specific heat capacity at constant pressure;

[0062] Step 6.2: Based on the airflow phase energy E g and air flow temperature T g The relationship between and is used to expand the energy conservation equation and obtain the temperature field model of the tuyere raceway, which is expressed as:

[0063]

[0064] Step 6.3: According to the divergence product rule, the product of velocity and temperature in formula (20) is diverged and the expanded temperature field model of the tuyere raceway is obtained, which is expressed as:

[0065]

[0066] Step 6.4: When the fluid in the blast furnace tuyere raceway is an incompressible gas, adjust formula (21) to obtain the mathematical model of the temperature field, which is expressed as:

[0067]

[0068] Step 6.5: Expand the temperature field mathematical model in two-dimensional rectangular coordinates using the gradient and Laplace operators to obtain the expanded temperature field mathematical model, which is expressed as:

[0069]

[0070] Step 6.6: Discretize the expanded temperature field mathematical model to obtain a discrete temperature field mathematical model, which is expressed as:

[0071]

[0072] in, is the airflow phase temperature on the grid at row i and column j in the n+1th iteration.

[0073] Optionally, step 7 specifically includes:

[0074] Step 7.1: For each grid in the target grid, calculate the volume fraction α of the airflow phase of the grid g , the velocity component u in the horizontal direction x x , the velocity component u in the vertical direction y y , the airflow phase turbulent kinetic energy k and the airflow phase turbulent dissipation rate ∈ are initialized, and the initial iteration number n is set to 0, and the initial airflow phase volume fraction α gAs the volume fraction of the airflow phase in the nth iteration, the initial horizontal velocity component u x The velocity component u in the horizontal direction x of the nth iteration x , the velocity component u in the vertical direction y y As the velocity component in the vertical direction y of the nth iteration, the initial airflow phase turbulent kinetic energy k is used as the airflow phase turbulent kinetic energy of the nth iteration, and the initial airflow phase turbulent dissipation rate ∈ is used as the airflow phase turbulent dissipation rate of the nth iteration;

[0075] Step 7.2: Based on the n-th iteration airflow phase volume fraction, the n-th iteration horizontal x velocity component, the n-th iteration vertical y velocity component, the n-th iteration airflow phase turbulent kinetic energy, and the n-th iteration airflow phase turbulent dissipation rate of all grids, substitute the n-th iteration airflow phase volume fraction, the n-th iteration horizontal x velocity component, and the n-th iteration vertical y velocity component into the governing equation of the discretized airflow phase volume fraction to obtain the n+1-th iteration airflow phase volume fraction. Then, calculate the difference between the n-th iteration airflow phase volume fraction and the n+1-th iteration airflow phase volume fraction and use it as the residual of the airflow phase volume fraction.

[0076] Step 7.3: Substitute the airflow-phase turbulent dissipation rate of the nth iteration, the airflow-phase turbulent kinetic energy of the nth iteration, the velocity component in the horizontal direction x of the nth iteration, and the velocity component in the vertical direction y of the nth iteration into the discrete turbulent kinetic energy equation to obtain the airflow-phase turbulent kinetic energy of the n+1th iteration. Calculate the difference between the airflow-phase turbulent kinetic energy of the nth iteration and the airflow-phase turbulent kinetic energy of the n+1th iteration and use it as the residual of the airflow-phase turbulent kinetic energy.

[0077] Step 7.4: Substitute the airflow-phase turbulent kinetic energy of the n+1th iteration, the airflow-phase turbulent dissipation rate of the nth iteration, the velocity component in the horizontal direction x of the nth iteration, and the velocity component in the vertical direction y of the nth iteration into the discrete airflow-phase turbulent dissipation rate equation to obtain the airflow-phase turbulent dissipation rate of the n+1th iteration. Calculate the difference between the airflow-phase turbulent dissipation rate of the nth iteration and the airflow-phase turbulent dissipation rate of the n+1th iteration and use it as the residual of the airflow-phase turbulent dissipation rate.

[0078] Step 7.5: Substitute the horizontal x velocity component of the nth iteration, the vertical y velocity component of the nth iteration, the airflow phase volume fraction of the n+1th iteration, the airflow phase turbulent kinetic energy of the n+1th iteration, and the airflow phase turbulent dissipation rate of the n+1th iteration into the expression of the discretized velocity component to obtain the horizontal x velocity component of the n+1th iteration and the vertical y velocity component of the n+1th iteration. Then, calculate the difference between the horizontal x velocity component of the nth iteration and the horizontal x velocity component of the n+1th iteration and use it as the residual of the horizontal x velocity component. Calculate the difference between the vertical y velocity component of the nth iteration and the vertical y velocity component of the n+1th iteration and use it as the residual of the vertical y velocity component.

[0079] The residual of the airflow phase volume fraction, the residual of the airflow phase turbulent kinetic energy, the residual of the airflow phase turbulent dissipation rate, the residual of the velocity component in the horizontal direction x and the residual of the velocity component in the vertical direction y constitute a residual set;

[0080] Step 7.6: Determine whether the residual set satisfies the preset conditions. If the residual set satisfies the preset conditions, use the volume fraction of the gas phase in the n+1th iteration as the final volume fraction field of the blast furnace tuyere raceway, and use the horizontal velocity component x in the n+1th iteration and the vertical velocity component y in the n+1th iteration as the final velocity field of the blast furnace tuyere raceway. If the residual set does not satisfy the preset conditions, set n = n+1 and return to step 7.2.

[0081] Among them, the residual set satisfies the preset conditions, indicating that the residual of the volume fraction of the airflow phase is less than the preset threshold, and the residual of the turbulent kinetic energy of the airflow phase is less than the preset threshold, and the residual of the turbulent dissipation rate of the airflow phase is less than the preset threshold, and the residual of the velocity component in the horizontal direction x is less than the preset threshold, and the residual of the velocity component in the vertical direction y is less than the preset threshold.

[0082] Optionally, step 8 specifically includes:

[0083] Step 8.1: For each grid in the target grid, calculate the airflow phase temperature T of the grid g Initialize and set the initial iteration number n to 0, and set the initial airflow phase temperature T g As the airflow phase temperature of the nth iteration;

[0084] Step 8.2: Substitute the airflow phase temperature of the nth iteration based on all grids into the discrete temperature field mathematical model to obtain the airflow phase temperature of the n+1th iteration, and calculate the difference between the airflow phase temperature of the nth iteration and the airflow phase temperature of the n+1th iteration as the residual of the airflow phase temperature;

[0085] Step 8.3: Determine whether the residual of the airflow phase temperature is less than the preset threshold. If the residual of the airflow phase temperature is less than the preset threshold, use the airflow phase temperature of the n+1th iteration as the final temperature field of the blast furnace tuyere vortex zone. If the residual of the airflow phase temperature is not less than the preset threshold, set n=n+1 and return to execute step 8.2.

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

[0087] An infrared camera is used to read the temperature of the tuyere vortex zone instead of the blast temperature as the model inlet temperature input; two-phase flow is used to model the blast furnace tuyere vortex zone, considering the influence of mass transfer and heat transfer between gas and solid on multi-physical fields; the code iteration process is written by ourselves to optimize the discrete control equations of the gas volume fraction field, gas velocity field and gas temperature field based on finite differences, while ignoring the high-order terms of viscous stress and Reynolds stress. By constructing a regular grid, the complexity and iteration cost of the finite element and finite volume discrete formats are avoided, thereby improving the iteration speed. BRIEF DESCRIPTION OF THE DRAWINGS

[0088] Figure 1 Schematic diagram of a flow chart of a method for calculating the temperature field in a blast furnace tuyere raceway based on gas-solid two-phase flow in an embodiment of the present invention;

[0089] Figure 2 A schematic diagram of the structure of the target grid of the modeling area provided by the present invention;

[0090] Figure 3 Schematic diagram of the iterative algorithm in an embodiment of the present invention;

[0091] Figure 4 Schematic diagram of the volume fraction field of the final blast furnace tuyere raceway in an embodiment of the present invention;

[0092] Figure 5 Schematic diagram of the velocity field of the final blast furnace tuyere raceway in an embodiment of the present invention;

[0093] Figure 6 Schematic diagram of the temperature field of the final blast furnace tuyere raceway in an embodiment of the present invention. DETAILED DESCRIPTION

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

[0095] In response to the problems existing in the prior art, the present invention uses an infrared camera to read the temperature of the vortex zone of the tuyere instead of the blast temperature as the model inlet temperature input, and establishes a gas-solid coupling model under the Euler-Lagrangian framework based on the gas-solid two-phase flow characteristics of the vortex zone of the tuyere. The body-fitting coordinate system is used to achieve grid adaptation of complex geometric boundaries, and the staggered grid technology is used to eliminate the odd-even oscillation phenomenon in the calculation of the pressure gradient of the same grid. The SIMPLE algorithm is combined to realize the iterative correction of the velocity-pressure field to meet the continuity equation constraints. Based on the principle of conservation of energy, a two-way coupling solution system for the temperature field, velocity field, and volume fraction field is constructed to complete the multi-physics field joint simulation.

[0096] Specifically, the present invention provides a method for calculating the temperature field of the blast furnace tuyere raceway based on gas-solid two-phase flow, combined with Figure 1 , which may include the following steps:

[0097] Step 1: Establish the modeling area boundary of the blast furnace tuyere raceway. Based on the modeling area boundary, generate the target grid of the modeling area through mesh generation technology. Figure 2 ;

[0098] Step 2: For each grid in the target grid, construct the governing equation for the discretized gas phase volume fraction based on the gas-solid two-phase flow;

[0099] Step 2.1: For each grid in the target grid, construct a control equation for the volume fraction of the gas phase based on the gas-solid two-phase flow. The control equation for the volume fraction of the gas phase is expressed as:

[0100]

[0101] Among them, α g is the volume fraction of the gas phase, which represents the volume fraction of the gas in a grid, t represents time, and u x is the air flow phase velocity The velocity component in the horizontal direction x, u y is the air flow phase velocity The velocity component in the vertical direction y;

[0102] For the multi-physics field control equations of the gas phase field, in theory, these equations have analytical solutions when specific boundary conditions are met. However, due to the complexity of boundary conditions, the irregularity of the geometric domain, and the nonlinear characteristics of the equations themselves in actual engineering problems, it is usually difficult to obtain closed-form analytical solutions. Therefore, it is usually necessary to use numerical methods to discretize the continuous solution domain into a finite number of spatial points, and use the values of the dependent variables at these discrete positions as basic unknowns for discrete solution. By constructing a discrete format for the control equations, a closed algebraic equation group corresponding to the spatial discretization can be established. The linear algebra numerical solution method is then used to solve the equation group, and finally the numerical solution at the discrete grid nodes is obtained. For the areas between discrete points, appropriate numerical interpolation functions or extrapolation methods can be constructed to reconstruct the full-field physical quantities based on the numerical solutions of adjacent nodes.

[0103] The finite difference method discretizes partial differential equations directly by replacing the differential form of the derivative. Its essence is to convert the continuous partial differential equation into a difference equation based on the Taylor expansion, thereby obtaining an approximate solution. The finite difference method is simple and easy to implement, and is applicable to the case of regular grids. In particular, the derivative terms in the partial differential equation can be directly discretized. Now take the physical quantity φ as an example to illustrate. The physical quantity can be the volume fraction α, temperature T, velocity component u x , velocity component u y Common differential forms include forward differential and central differential, so the present invention can calculate the volume fraction α, temperature T, velocity component u x , velocity component u y Discretize.

[0104] Step 2.2: Discretize the governing equation of the airflow phase volume fraction. Specifically, for the transient term, use forward difference for the time gradient and use central difference for the spatial gradient to discretize the model. The governing equation of the discretized airflow phase volume fraction is obtained, which is expressed as:

[0105]

[0106] Among them, α g n+1 [i,j] is the volume fraction of the airflow phase on the grid located in the i-th row and the j-th column in the n+1th iteration, Δt is the time difference, is the horizontal velocity component on the grid at row i and column j in the nth iteration, Δx is the horizontal difference, and Δy is the vertical difference.

[0107] Step 3: For each grid in the target grid, construct an expression for the discretized velocity component based on the gas-solid two-phase flow;

[0108] Step 3.1: For each grid in the target grid, the momentum conservation equation is adjusted based on the gas-solid two-phase flow to obtain an adjusted momentum conservation equation, which is expressed as:

[0109]

[0110] Among them, ρ g is the air flow phase density, is the viscous stress, is the Reynolds stress, p is the pressure, g is the acceleration of gravity, F g is the interphase force, which includes the drag force F drag , and the lift force F lift and virtual body force F vm ;

[0111] Now we simplify each term in formula (3) in turn. First, we consider the transient term and the convection term. Under the assumption that the air flow phase is an incompressible gas, the transient term and the convection term can be simplified to the equation shown in (25).

[0112]

[0113] Considering the viscous stress term, viscous stress The expression of is shown in formula (26).

[0114]

[0115] Among them, μ g is the kinematic viscosity of gas phase fluid molecules. The kinematic viscosity of air at 1800°C is 3.8×10 -4 m 2 / s, I is the unit matrix.

[0116] Reynolds stress The expression of is shown in formula (27).

[0117]

[0118] Among them, μ g,t is the turbulent kinematic viscosity of the air flow phase, and k is the turbulent kinetic energy of the air flow phase.

[0119] Turbulent kinematic viscosity μ of the air phase g,t The expression of is shown in formula (28).

[0120]

[0121] Among them, ∈ is the turbulent dissipation rate of the air phase, C μ is a constant, and its empirical value is 0.09.

[0122] Now consider the interphase force F g In the gas-solid two-phase flow model, the dominant external forces include the drag force F drag , and the lift force F lift and virtual body force F vm Drag is the primary resistance exerted by the fluid on particle motion. It is generated by the combined action of viscous shear force and pressure differential force between the fluid and particle surfaces. Its direction is opposite to the relative velocity of the particle and fluid and is the dominant factor in the momentum exchange between the phases. Its expression is shown in Equation (29).

[0123]

[0124] Where β is the drag coefficient, is the flow velocity of the discrete phase (coke), R ep is the Reynolds stress coefficient, d p,b is the coke particle size.

[0125] The lift force is mainly caused by the velocity gradient of the flow field or the rotation of the particles, and is perpendicular to the direction of the relative velocity. Its expression is shown in formula (31).

[0126]

[0127] Among them, α coke is the volume fraction of the discrete phase (solid), ρ coke is the coke density, C l is a constant, and its empirical value is 0.5.

[0128] The virtual body force refers to the force that, when a particle accelerates in an unsteady flow field, it must drive the surrounding fluid to move together, which is equivalent to an increase in the particle's mass. This effect is characterized by the virtual body force, which is in the opposite direction to the relative acceleration of the particle. Its expression is shown in Equation (33).

[0129]

[0130] Among them, C vm is the virtual mass coefficient, which is 0.5 for spherical particles, and D / Dt is the reciprocal operator of the material.

[0131] The material reciprocal operator is defined as shown in formula (34):

[0132]

[0133] Since the velocity field model is relatively complex, each part is simplified and expanded here. When expanding, the original mathematical model is expanded into an expression in two-dimensional rectangular coordinates. For non-constant scalars, such as the volume fraction (field) α, it is a binary function about (x, y). For vectors, such as the gas phase velocity field There are two components (u x ,u y ), which represent the horizontal velocity component and the vertical velocity component respectively.

[0134] First, the transient term and the convection term. Formula (25) gives the rearranged expression, which is expanded as shown in formula (35).

[0135]

[0136] Formula (26) gives the expression of viscous stress. Here, the viscous stress term in the momentum equation is expanded, and the expanded expression is shown in (36).

[0137]

[0138] Formula (27) gives the expression of Reynolds stress. Here, the Reynolds stress term in the momentum equation is expanded, and the expanded expression is shown in (37).

[0139]

[0140] Step 3.2: Based on viscous stress Formula, Reynolds stress Formula, drag force F drag Formula, lift F lift Formula and virtual body force F vm Formula, expand and adjust formula (3), and the expression of velocity component is expressed as:

[0141]

[0142] Where μ is the molecular viscosity of the fluid, μ g,t is the turbulent kinematic viscosity of the air phase, F x is the interphase force F g The horizontal component, F y is the interphase force F g The vertical component of

[0143] Step 3.3: Discretize the expression of the velocity component. Specifically, use forward difference for the time derivative and central difference for the spatial derivative to obtain the expression of the discretized velocity component, which is expressed as:

[0144]

[0145]

[0146] Among them, u x n+1[i,j] is the horizontal velocity component on the grid at row i and column j in the n+1th iteration. is the vertical velocity component on the grid at row i+1 and column j in the nth iteration; p n [i+1,j] is the pressure on the grid at row i+1 and column j in the nth iteration, and is the F on the grid at row i and column j. x , F y [i,j] is the F on the grid located in row i and column j y .

[0147] Step 4: Construct the turbulent kinetic energy equation of the airflow phase, and then process the turbulent kinetic energy equation of the airflow phase to obtain a discrete turbulent kinetic energy equation;

[0148] Step 4.1: Construct the turbulent kinetic energy equation of the air flow phase, which is expressed as:

[0149]

[0150] Among them, k is the turbulent kinetic energy of the air flow phase, ∈ is the turbulent dissipation rate of the air flow phase, μ k ' is the intermediate variable, G k is the generation term of turbulent kinetic energy, μ k ' is expressed as:

[0151]

[0152] Among them, μ t is the turbulent viscosity, σ k is the turbulent kinetic energy Prandtl number, C μ is a constant;

[0153] Among them, G k Expressed as:

[0154] G k =μ t S i,j S i,j (11)

[0155]

[0156] Among them, S i,j is the intermediate variable, u i represents the i-th component of velocity u, u j represents the jth component of velocity u, x i It represents the i-th component of the coordinate basis vector, x j It represents the jth component of the coordinate basis vector;

[0157] Step 4.2: Expand the turbulent kinetic energy equation of the airflow phase to obtain the expanded turbulent kinetic energy equation, which is expressed as:

[0158]

[0159] Step 4.3: Discretize the expanded turbulent kinetic energy equation to obtain the discrete turbulent kinetic energy equation, which is expressed as:

[0160]

[0161] Among them, k n+1 [i,j] is the turbulent kinetic energy of the airflow phase on the grid located in the i-th row and the j-th column in the n+1th iteration, μ k ' n [i,j] is μ on the grid at row i and column j in the nth iteration. k '; is the generated term of turbulent kinetic energy on the grid at row i and column j in the nth iteration, ∈ n [i, j] is the turbulent dissipation rate of the airflow phase on the grid located in the i-th row and j-th column in the n-th iteration.

[0162] Step 5: Construct the airflow phase turbulence dissipation rate equation, process the airflow phase turbulence dissipation rate equation, and obtain the discrete airflow phase turbulence dissipation rate equation;

[0163] Step 5.1: Construct the turbulent dissipation rate equation for the airflow phase, expressed as:

[0164]

[0165] Among them, σ ∈ is the dissipation rate Prandtl number, C ∈1 and C ∈2 is the value calibrated by experiment;

[0166] Step 5.2: Expand the airflow phase turbulence dissipation rate equation to obtain the expanded airflow phase turbulence dissipation rate equation, which is expressed as:

[0167]

[0168] Step 5.3: Discretize the expanded airflow phase turbulence dissipation rate equation to obtain the discrete airflow phase turbulence dissipation rate equation, which is expressed as:

[0169]

[0170] Among them, ∈ n+1 [i,j] is the turbulent dissipation rate of the airflow phase on the grid located in the i-th row and the j-th column in the n+1th iteration, is the turbulent viscosity on the grid at row i and column j in the nth iteration; is the turbulent viscosity on the grid at row i+1 and column j in the nth iteration.

[0171] Step 6: Construct a discrete temperature field mathematical model by combining fluid flow, chemical reaction exothermicity, chemical reaction endothermicity, chemical reaction heat transfer, and phase change processes;

[0172] Specifically, we first analyze the chemical reactions in the blast furnace tuyere raceway:

[0173] The main chemical reactions occurring in the vortex zone of the blast furnace tuyere are the pyrolysis and combustion of pulverized coal. The pyrolysis of coal is a complex heterogeneous reaction involving the cracking of macromolecular organic matter, generating gases (such as CO, H2, CH4, etc.), liquids (tar), and solids (semi-coke). Because the pyrolysis process is extremely short, usually completed in milliseconds, and is affected by the combined effects of high temperature and rapid flow conditions, when establishing mathematical models, it is usually assumed that the pyrolysis process is completed instantaneously, thereby focusing on the impact of the pulverized coal combustion process on material balance and energy transfer. The chemical reaction equation for pulverized coal combustion is shown in (38).

[0174]

[0175] At the boundary of the cavity in the blast furnace tuyere raceway and the adjacent coke layer outside, the main chemical reactions are coke combustion and coke reduction, and the chemical reaction equations are shown in (39).

[0176]

[0177] At the same time, the chemical reaction rate in the blast furnace tuyere raceway is analyzed:

[0178] The chemical reaction rate of the gas-solid two-phase is affected not only by the gas concentration, but also by the solid surface area. For the chemical reaction rate of the gas-solid two-phase, the simplified reaction rate expression is shown as (40).

[0179] r=kNS (40)

[0180] Where r represents the chemical reaction rate, N represents the gas concentration, and S is the surface area of the solid participating in the reaction.

[0181] k is calculated by the Larenius equation, and its expression is shown in formula (41), the expression of the gas concentration N is shown in formula (42), and the expression of the surface area S of the reacting solid is shown in formula (43).

[0182]

[0183] At the same time, the heat in the blast furnace tuyere raceway is analyzed:

[0184] When calculating the heat generated by chemical reactions and the heat transfer between gas and solid, consider dividing the modeling area into several microcubes and studying the heat generation and transfer in each microcube separately. Then, the heat generated by the chemical reaction in a single cube is calculated as shown in (44).

[0185] Q1=∫∫∫qρdV (44)

[0186] Where q is the equivalent calorific value of a single cube (J / kg), which represents the heat generated per unit mass; ρ is the equivalent density of the cube, which is calculated by weighting the mass ratio of pulverized coal to coke.

[0187] Considering that the chemical reactions in a single cube mainly include the combustion of coal powder, the combustion of coke, and the reduction of coke, the expression of the equivalent calorific value q is shown in (45).

[0188]

[0189] Where Q 11 is the heat generated by the combustion of unit mass of pulverized coal (J / kg), M c is the molar mass of coal powder (kg / mol), r 11 is the chemical reaction rate of pulverized coal combustion (mol / (m 3 ·s)),ρ c is the density of pulverized coal, Q 12 is the heat generated by the combustion of unit mass of coke, M coke is the molar mass of coke, r 12 is the chemical reaction rate of coke combustion, Q 13 is the heat generated by the reduction reaction of unit mass of coke, r 13 is the chemical reaction rate of coke reduction, ρ coke is the density of coke.

[0190] Furthermore, a heat transfer model is constructed:

[0191] In the tuyere area, heat transfer occurs between the gas phase fluid and the solid coke, and the gas-solid two-phase heat transfer model is shown in (46).

[0192] Q2=hS(T coke -T g ) (46)

[0193] Where h is the heat transfer coefficient between gas and solid, S is the surface area of the solid, T coke is the temperature of coke, T g is the temperature of the air flow.

[0194] The expression of the heat transfer coefficient h between gas and solid is shown in (47).

[0195]

[0196] Where k g is the thermal conductivity of the gas, d coke is the particle size of the solid, N u is the Nusselt number.

[0197] Nusselt number N u The expression of is shown in (48).

[0198]

[0199] Where Re is the Reynolds number, Pr g is the Prandtl number for convective heat transfer of gas.

[0200] Reynolds number Re and Prandtl number Pr of gas convective heat transfer g The calculation formula is shown in formula (49).

[0201]

[0202] Where η g is the gas viscosity, ψ s is the shape factor of the solid particles, d s is the equivalent particle size, C p,g is the specific heat of gas at constant pressure.

[0203] Step 6.1: Based on the above formula, combined with fluid flow, chemical reaction exothermicity, chemical reaction endothermicity, chemical reaction heat transfer, and phase change processes, construct the energy conservation equation, which is expressed as:

[0204]

[0205] Among them, E g is the energy of the airflow phase, k g is the thermal conductivity of the gas, Q1 is the heat generated by the chemical reaction, Q2 is the heat transferred between the gas and the solid, T g is the air flow phase temperature;

[0206] Among them, the air flow phase energy E g and air flow temperature T g The relationship between them is expressed as:

[0207] E g =ρ g C p T g (19)

[0208] Among them, Cp is the specific heat capacity at constant pressure;

[0209] Step 6.2: Based on the airflow phase energy E g and air flow temperature T g The relationship between and is used to expand the energy conservation equation and obtain the temperature field model of the tuyere raceway, which is expressed as:

[0210]

[0211] Step 6.3: According to the divergence product rule, the product of velocity and temperature in formula (20) is diverged and the expanded temperature field model of the tuyere raceway is obtained, which is expressed as:

[0212]

[0213] Step 6.4: When the fluid in the blast furnace tuyere raceway is an incompressible gas, adjust formula (21) to obtain the mathematical model of the temperature field, which is expressed as:

[0214]

[0215] Step 6.5: Expand the temperature field mathematical model in two-dimensional rectangular coordinates using the gradient and Laplace operators to obtain the expanded temperature field mathematical model, which is expressed as:

[0216]

[0217] Step 6.6: Discretize the expanded temperature field mathematical model. Specifically, use forward difference for the time derivative and central difference for the spatial derivative to obtain the discrete temperature field mathematical model, which can be expressed as:

[0218]

[0219] in, is the airflow phase temperature on the grid at row i and column j in the n+1th iteration.

[0220] Step 7: Based on the control equation of the discretized gas phase volume fraction, the expression of the discretized velocity component, the discrete turbulent kinetic energy equation, and the discrete gas phase turbulent dissipation rate equation, iteratively solve to obtain the final volume fraction field and the final velocity field of the blast furnace tuyere raceway;

[0221] Step 7.1: Combine Figure 3 , for each grid in the target grid, the airflow phase volume fraction α of the grid g , the velocity component u in the horizontal direction x x , the velocity component u in the vertical direction y y, the airflow phase turbulent kinetic energy k and the airflow phase turbulent dissipation rate ∈ are initialized, and the initial iteration number n is set to 0, and the initial airflow phase volume fraction α g As the volume fraction of the airflow phase in the nth iteration, the initial horizontal velocity component u x The velocity component u in the horizontal direction x of the nth iteration x , the velocity component u in the vertical direction y y As the velocity component in the vertical direction y of the nth iteration, the initial airflow phase turbulent kinetic energy k is used as the airflow phase turbulent kinetic energy of the nth iteration, and the initial airflow phase turbulent dissipation rate ∈ is used as the airflow phase turbulent dissipation rate of the nth iteration;

[0222] In the specific implementation process, for the initial conditions of the gas phase volume fraction, at the tuyere position, the gas phase volume fraction α g Initialized to 1, and the rest of the positions are initialized to 0.1. It should be noted that in the expression of the discretized velocity component, the gas phase volume fraction appears in the denominator. If it is initialized to 0, it will cause a division by zero anomaly and destroy the balance of the equation. Therefore, after each update of the volume fraction, it will also be corrected to avoid the volume fraction having a value of 0.

[0223] For the boundary conditions of the gas velocity field, the left boundary of the modeling region is considered a symmetric boundary and its value is 0. The right and lower boundaries are considered walls, and the gas volume fraction at these boundaries is equal to the gas volume fraction close to the wall. The upper boundary is considered a free flow outlet, and the gradient of the gas volume fraction is 0.

[0224] For the initial conditions of the gas phase velocity field, at the tuyere position, the radial velocity u x Set to -258m / s (the symbol indicates the opposite direction of the positive direction of the x-axis), the axial speed u y It is set to 90 m / s, and the radial and axial velocities at other positions are initialized to 0.

[0225] For the boundary conditions, the left boundary is regarded as a symmetric boundary, requiring the radial velocity to be 0 and the radial gradient of the axial velocity to be 0, which is expressed as follows

[0226]

[0227] The right and lower boundaries are considered as walls. Assuming no slip, the velocity close to the wall is equal to the wall velocity. The upper boundary is a free outlet with a velocity gradient of 0.

[0228] The initial conditions for the gas temperature field are 1800°C at the center of the tuyere and 1500°C at all other locations. For the boundary conditions, the bottom and right boundaries are considered to be walls. Assuming no slip, the temperature close to the wall is equal to the wall temperature. The upper boundary is a free outlet with a temperature gradient of 0. The left boundary is a symmetric interface with a temperature gradient of 0.

[0229] The volume fraction, velocity, and temperature were all solved iteratively. The key control parameters included the iteration step size Δt and the unit grid lengths Δx and Δy. In subsequent experiments, the iteration step size Δt was 0.0001 s, and the unit grid sizes Δx and Δy were 0.05 m and 0.05 m, respectively.

[0230] Step 7.2: Based on the n-th iteration airflow phase volume fraction, the n-th iteration horizontal x velocity component, the n-th iteration vertical y velocity component, the n-th iteration airflow phase turbulent kinetic energy, and the n-th iteration airflow phase turbulent dissipation rate of all grids, substitute the n-th iteration airflow phase volume fraction, the n-th iteration horizontal x velocity component, and the n-th iteration vertical y velocity component into the governing equation of the discretized airflow phase volume fraction to obtain the n+1-th iteration airflow phase volume fraction. Then, calculate the difference between the n-th iteration airflow phase volume fraction and the n+1-th iteration airflow phase volume fraction and use it as the residual of the airflow phase volume fraction.

[0231] Step 7.3: Substitute the airflow-phase turbulent dissipation rate of the nth iteration, the airflow-phase turbulent kinetic energy of the nth iteration, the velocity component in the horizontal direction x of the nth iteration, and the velocity component in the vertical direction y of the nth iteration into the discrete turbulent kinetic energy equation to obtain the airflow-phase turbulent kinetic energy of the n+1th iteration. Calculate the difference between the airflow-phase turbulent kinetic energy of the nth iteration and the airflow-phase turbulent kinetic energy of the n+1th iteration and use it as the residual of the airflow-phase turbulent kinetic energy.

[0232] Step 7.4: Substitute the airflow-phase turbulent kinetic energy of the n+1th iteration, the airflow-phase turbulent dissipation rate of the nth iteration, the velocity component in the horizontal direction x of the nth iteration, and the velocity component in the vertical direction y of the nth iteration into the discrete airflow-phase turbulent dissipation rate equation to obtain the airflow-phase turbulent dissipation rate of the n+1th iteration. Calculate the difference between the airflow-phase turbulent dissipation rate of the nth iteration and the airflow-phase turbulent dissipation rate of the n+1th iteration and use it as the residual of the airflow-phase turbulent dissipation rate.

[0233] Step 7.5: Substitute the horizontal x velocity component of the nth iteration, the vertical y velocity component of the nth iteration, the airflow phase volume fraction of the n+1th iteration, the airflow phase turbulent kinetic energy of the n+1th iteration, and the airflow phase turbulent dissipation rate of the n+1th iteration into the expression of the discretized velocity component to obtain the horizontal x velocity component of the n+1th iteration and the vertical y velocity component of the n+1th iteration. Then, calculate the difference between the horizontal x velocity component of the nth iteration and the horizontal x velocity component of the n+1th iteration and use it as the residual of the horizontal x velocity component. Calculate the difference between the vertical y velocity component of the nth iteration and the vertical y velocity component of the n+1th iteration and use it as the residual of the vertical y velocity component.

[0234] The residual of the airflow phase volume fraction, the residual of the airflow phase turbulent kinetic energy, the residual of the airflow phase turbulent dissipation rate, the residual of the velocity component in the horizontal direction x and the residual of the velocity component in the vertical direction y constitute a residual set;

[0235] Step 7.6: Determine whether the residual set satisfies the preset conditions. If the residual set satisfies the preset conditions, use the volume fraction of the gas phase in the n+1th iteration as the final volume fraction field of the blast furnace tuyere raceway, and use the horizontal velocity component x in the n+1th iteration and the vertical velocity component y in the n+1th iteration as the final velocity field of the blast furnace tuyere raceway. If the residual set does not satisfy the preset conditions, set n = n+1 and return to step 7.2.

[0236] Among them, the residual set satisfies the preset conditions, indicating that the residual of the volume fraction of the airflow phase is less than the preset threshold, and the residual of the turbulent kinetic energy of the airflow phase is less than the preset threshold, and the residual of the turbulent dissipation rate of the airflow phase is less than the preset threshold, and the residual of the velocity component in the horizontal direction x is less than the preset threshold, and the residual of the velocity component in the vertical direction y is less than the preset threshold.

[0237] Among them, the final volume fraction field of the blast furnace tuyere raceway is as follows: Figure 4 , where the origin represents the center of the blast furnace, the horizontal axis represents the radius of the blast furnace, and the vertical axis represents the height of the blast furnace. There are three key factors that affect the gas phase volume fraction. The first is the change in gas composition caused by chemical reactions, such as the combustion reaction consuming O2 and generating products such as CO2. The second is the change in the volume fraction of discrete phase particles (coke) in the microelement control volume, whose spatial distribution directly affects the effective gas phase volume. The third factor is the diffusion and mass transfer effect between multi-component gases, which leads to changes in the gas concentration gradient in the local area.

[0238] The analysis of the modeling results is as follows: At the blast inlet, assuming that there is no particle phase, since the injection medium is oxygen-enriched air, the initial gas phase volume fraction can be set to 1. As the gas moves in the radial / axial direction, the pulverized coal undergoes pyrolysis and combustion reactions: O2 is continuously consumed and CO2 is generated. At the same time, the concentration of coke particles increases gradually in the radial and axial directions, resulting in a decrease in the gas phase volume fraction due to dual effects, which is affected by both the change in the number of moles of the reaction gas and the volume exclusion effect of the particle phase. This numerical distribution characteristic is consistent with Figure 4 The simulation results are consistent with those of , which verifies the rationality of the discrete phase-continuous phase coupling model.

[0239] Among them, the velocity field of the final blast furnace tuyere raceway is as follows: Figure 5 , where the origin represents the blast furnace center, the abscissa represents the blast furnace radius, and the ordinate represents the blast furnace height. The main factors influencing the gas phase velocity field distribution are the direction of the initial velocity and the direction of the net external force. Based on the vector decomposition of the initial velocity, the initial velocity of the injected gas is the superposition of the radial and axial velocities, with the radial velocity being approximately twice the axial velocity. Assuming the components of the net external force are equal in both directions, i.e., the acceleration is the same, and according to kinematic formulas, displacement is proportional to the square of the velocity. Ideally, the radial displacement should be four times the axial displacement. However, in actual operation, the presence of a dead material column in the tuyere cavity, as shown by the lower left slope of the model, significantly alters the net external force distribution: the force exerted by the dead material column on the gas is perpendicular to the slope and can be decomposed into components in the positive axial direction (upward) and the positive radial direction (rightward). The radial velocity is in the opposite direction to the radial force component of the dead material column, causing it to decay faster due to resistance. The axial velocity is in the same direction as the axial force component of the dead material column, providing dynamic support and slowing its decay.

[0240] Therefore, although the coke layer's macroscopic resistance in both directions is similar, the mechanical decomposition effect of the dead stock column causes the radial velocity to decay faster than the axial velocity, resulting in the ratio of the actual radial displacement to the axial displacement being four times less than the ideal value. Furthermore, the farther from the blast inlet, the more significant the dead stock column's interference with the local flow field, further reducing the ratio of radial to axial displacement at the same terminal velocity. Figure 5 The results show that the radial displacement at the typical position is about twice that of the axial displacement, and the ratio decreases with increasing distance, which is consistent with the above-mentioned mechanical mechanism and model prediction.

[0241] Step 8: Based on the final velocity field, the discrete temperature field mathematical model is iteratively solved to obtain the final temperature field of the blast furnace tuyere raceway.

[0242] Step 8.1: For each grid in the target grid, calculate the airflow phase temperature T of the grid g Initialize and set the initial iteration number n to 0, and set the initial airflow phase temperature T g As the airflow phase temperature of the nth iteration;

[0243] Step 8.2: Substitute the airflow phase temperature of the nth iteration based on all grids into the discrete temperature field mathematical model to obtain the airflow phase temperature of the n+1th iteration, and calculate the difference between the airflow phase temperature of the nth iteration and the airflow phase temperature of the n+1th iteration as the residual of the airflow phase temperature;

[0244] Step 8.3: Determine whether the residual of the airflow phase temperature is less than the preset threshold. If the residual of the airflow phase temperature is less than the preset threshold, use the airflow phase temperature of the n+1th iteration as the final temperature field of the blast furnace tuyere vortex zone. If the residual of the airflow phase temperature is not less than the preset threshold, set n=n+1 and return to execute step 8.2.

[0245] Combine Figure 6 , Figure 6 The figure is a schematic diagram of the temperature field of the final blast furnace tuyere raceway, where the origin represents the center of the blast furnace, the horizontal axis represents the radius of the blast furnace, and the vertical axis represents the height of the blast furnace. The core factors affecting the temperature distribution include the thermal effect of the chemical reaction (heat absorption or heat release) and the heat transfer process inside and outside the system. Figure 6 The rationality of the model results is analyzed from the perspective of thermodynamics and reaction kinetics: in the blast inlet cavity area, the rapid combustion and cracking reaction of pulverized coal dominate (ΔH < 0, strong exothermic), resulting in local high temperature concentration; the reaction type in the area outside the cavity changes - coke combustion (ΔH < 0, exothermic) and the reduction reaction of coke and CO2 (ΔH > 0, endothermic) coexist, and at the same time, the direct reduction reaction of coke and FeO in the bosh area (ΔH > 0, endothermic) further intensifies the endothermic effect.

[0246] Therefore, as the spatial position moves away from the cavity, the intensity of the exothermic reaction decreases while the proportion of the endothermic reaction increases, forming a decreasing trend in the temperature gradient. Figure 6 The results show that the highest temperature zone is concentrated in the core of the cyclotron cavity, with the temperature decreasing gradually radially and axially outward, consistent with the aforementioned reaction thermal coupling mechanism. It is important to note that the competitive effect of endothermic reactions outside the cavity (such as C + CO2 → 2CO and C + FeO → Fe + CO) significantly weakens the temperature rise potential, while convective heat transfer and radiative heat dissipation between the gas and solid phases further enhance the spatial attenuation of the temperature field. The temperature distribution morphology and gradient variation in the simulation results are consistent with the thermal equilibrium characteristics of a multi-reaction coupled system, verifying the physical rationality of the model.

[0247] The simulation results show that the total absolute error at the six measurement points is 168.178°C, the average absolute error at the measurement points is 28.03°C, the standard deviation of the absolute error is 16.14°C, the average relative error is 1.644%, and the model error is controlled within the engineering threshold of ±2.5%.

[0248] 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 method for calculating the temperature field in the blast furnace tuyere raceway based on gas-solid two-phase flow, characterized in that: include: Step 1: Establish the modeling area boundary of the blast furnace tuyere raceway, and generate the target mesh of the modeling area based on the modeling area boundary through mesh generation technology; Step 2: For each grid in the target grid, construct the governing equation for the discretized gas phase volume fraction based on the gas-solid two-phase flow; Step 3: For each grid in the target grid, construct an expression for the discretized velocity component based on the gas-solid two-phase flow; Step 4: Construct the turbulent kinetic energy equation of the airflow phase, and then process the turbulent kinetic energy equation of the airflow phase to obtain a discrete turbulent kinetic energy equation; Step 5: Construct the airflow phase turbulence dissipation rate equation, process the airflow phase turbulence dissipation rate equation, and obtain the discrete airflow phase turbulence dissipation rate equation; Step 6: Construct a discrete temperature field mathematical model by combining fluid flow, chemical reaction exothermicity, chemical reaction endothermicity, chemical reaction heat transfer, and phase change processes; Step 7: Based on the control equation of the discretized gas phase volume fraction, the expression of the discretized velocity component, the discrete turbulent kinetic energy equation, and the discrete gas phase turbulent dissipation rate equation, iteratively solve to obtain the final volume fraction field and the final velocity field of the blast furnace tuyere raceway; Step 8: Based on the final velocity field, the discrete temperature field mathematical model is iteratively solved to obtain the final temperature field of the blast furnace tuyere raceway.

2. The method for calculating the temperature field of the blast furnace tuyere raceway based on gas-solid two-phase flow according to claim 1, characterized in that: Step 2 specifically includes: Step 2.1: For each grid in the target grid, construct a control equation for the volume fraction of the gas phase based on the gas-solid two-phase flow. The control equation for the volume fraction of the gas phase is expressed as: Among them, α g is the volume fraction of the gas phase, which represents the volume fraction of the gas in a grid, t represents time, and u x is the air flow phase velocity The velocity component in the horizontal direction x, u y is the air flow phase velocity The velocity component in the vertical direction y; Step 2.2: Discretize the governing equation of the airflow phase volume fraction to obtain the governing equation of the discretized airflow phase volume fraction, which is expressed as: Among them, α g n+1 [i,j] is the volume fraction of the airflow phase on the grid located in the i-th row and the j-th column in the n+1th iteration, Δt is the time difference, is the horizontal velocity component on the grid at row i and column j in the nth iteration, Δx is the horizontal difference, and Δy is the vertical difference.

3. The method for calculating the temperature field of the blast furnace tuyere raceway based on gas-solid two-phase flow according to claim 2, characterized in that: Step 3 specifically includes: Step 3.1: For each grid in the target grid, the momentum conservation equation is adjusted based on the gas-solid two-phase flow to obtain an adjusted momentum conservation equation, which is expressed as: Among them, ρ g is the air flow phase density, is the viscous stress, is the Reynolds stress, p is the pressure, g is the acceleration of gravity, F g is the interphase force, which includes the drag force F drag , and the lift force F lift and virtual body force F vm ; Step 3.2: Based on viscous stress Formula, Reynolds stress Formula, drag force F drag Formula, lift F lift Formula and virtual body force F vm Formula (3) is expanded and adjusted to obtain the expression of velocity component: Where μ is the molecular viscosity of the fluid, μ g,t is the turbulent kinematic viscosity of the air phase, F x is the interphase force F g The horizontal component, F y is the interphase force F g The vertical component of Step 3.3: Discretize the expression of the velocity component to obtain the expression of the discretized velocity component, which is expressed as: Among them, u x n+1 [i,j] is the horizontal velocity component on the grid at row i and column j in the n+1th iteration. is the vertical velocity component on the grid at row i+1 and column j in the nth iteration; p n [i+1,j] is the pressure on the grid at row i+1 and column j in the nth iteration, F x [i,j] is the F on the grid located in row i and column j x , F y [i,j] is the F on the grid located in row i and column j y .

4. The method for calculating the temperature field of the blast furnace tuyere raceway based on gas-solid two-phase flow according to claim 3, characterized in that: Step 4 specifically includes: Step 4.1: Construct the turbulent kinetic energy equation of the air flow phase, which is expressed as: Among them, k is the turbulent kinetic energy of the air flow phase, ∈ is the turbulent dissipation rate of the air flow phase, μ k ' is the intermediate variable, G k is the generation term of turbulent kinetic energy, μ k ' is expressed as: Among them, μ t is the turbulent viscosity, σ k is the turbulent kinetic energy Prandtl number, C μ is a constant; Among them, G k Expressed as: G k =μ t S i,j S i,j (11) Among them, S i,j is the intermediate variable, u i represents the i-th component of velocity u, u j represents the jth component of velocity u, x i It represents the i-th component of the coordinate basis vector, x j It represents the jth component of the coordinate basis vector; Step 4.2: Expand the turbulent kinetic energy equation of the airflow phase to obtain the expanded turbulent kinetic energy equation, which is expressed as: Step 4.3: Discretize the expanded turbulent kinetic energy equation to obtain the discrete turbulent kinetic energy equation, which is expressed as: Among them, k n+1 [i,j] is the turbulent kinetic energy of the airflow phase on the grid located in the i-th row and the j-th column in the n+1th iteration, μ k ' n [i,j] is μ on the grid at row i and column j in the nth iteration. k '; is the generated term of turbulent kinetic energy on the grid at row i and column j in the nth iteration, ∈ n [i, j] is the turbulent dissipation rate of the airflow phase on the grid located in the i-th row and j-th column in the n-th iteration.

5. The method for calculating the temperature field of the blast furnace tuyere raceway based on gas-solid two-phase flow according to claim 4, characterized in that: Step 5 specifically includes: Step 5.1: Construct the turbulent dissipation rate equation for the airflow phase, expressed as: Among them, σ ∈ is the dissipation rate Prandtl number, C ∈1 and C ∈2 is the value calibrated by experiment; Step 5.2: Expand the airflow phase turbulence dissipation rate equation to obtain the expanded airflow phase turbulence dissipation rate equation, which is expressed as: Step 5.3: Discretize the expanded airflow phase turbulence dissipation rate equation to obtain the discrete airflow phase turbulence dissipation rate equation, which is expressed as: Among them, ∈ n+1 [i,j] is the turbulent dissipation rate of the airflow phase on the grid located in the i-th row and the j-th column in the n+1th iteration, is the turbulent viscosity on the grid at row i and column j in the nth iteration; is the turbulent viscosity on the grid at row i+1 and column j in the nth iteration.

6. The method for calculating the temperature field of the blast furnace tuyere raceway based on gas-solid two-phase flow according to claim 5, characterized in that: Step 6 specifically includes: Step 6.1: Combine fluid flow, chemical reaction exothermicity, chemical reaction endothermicity, chemical reaction heat transfer, and phase change processes to construct the energy conservation equation, which can be expressed as: Among them, E g is the energy of the air flow phase, k g is the thermal conductivity of the gas, q1 is the heat generated by the chemical reaction, q2 is the heat transferred between the gas and the solid, T g is the air flow phase temperature; Among them, the air flow phase energy E g and air flow temperature T g The relationship between them is expressed as: E g =ρ g C p T g (19) Among them, C p is the specific heat capacity at constant pressure; Step 6.2: Based on the airflow phase energy E g and air flow temperature T g The relationship between and is used to expand the energy conservation equation and obtain the temperature field model of the tuyere raceway, which is expressed as: Step 6.3: According to the divergence product rule, the product of velocity and temperature in formula (20) is diverged and the temperature field model of the tuyere raceway is obtained after the expansion, which is expressed as: Step 6.4: When the fluid in the blast furnace tuyere raceway is an incompressible gas, adjust formula (21) to obtain the mathematical model of the temperature field, which is expressed as: Step 6.5: Expand the temperature field mathematical model in two-dimensional rectangular coordinates using the gradient and Laplace operators to obtain the expanded temperature field mathematical model, which is expressed as: Step 6.6: Discretize the expanded temperature field mathematical model to obtain a discrete temperature field mathematical model, which is expressed as: in, is the airflow phase temperature on the grid at row i and column j in the n+1th iteration.

7. The method for calculating the temperature field of the blast furnace tuyere raceway based on gas-solid two-phase flow according to claim 6, characterized in that: Step 7 specifically includes: Step 7.1: For each grid in the target grid, calculate the volume fraction α of the airflow phase of the grid g , the velocity component u in the horizontal direction x x , the velocity component u in the vertical direction y y , the airflow phase turbulent kinetic energy k and the airflow phase turbulent dissipation rate ∈ are initialized, and the initial iteration number n is set to 0, and the initial airflow phase volume fraction α g As the volume fraction of the airflow phase in the nth iteration, the initial horizontal velocity component u x The velocity component u in the horizontal direction x of the nth iteration x , the velocity component u in the vertical direction y y As the velocity component in the vertical direction y of the nth iteration, the initial airflow phase turbulent kinetic energy k is used as the airflow phase turbulent kinetic energy of the nth iteration, and the initial airflow phase turbulent dissipation rate ∈ is used as the airflow phase turbulent dissipation rate of the nth iteration; Step 7.2: Based on the n-th iteration airflow phase volume fraction, the n-th iteration horizontal x velocity component, the n-th iteration vertical y velocity component, the n-th iteration airflow phase turbulent kinetic energy, and the n-th iteration airflow phase turbulent dissipation rate of all grids, substitute the n-th iteration airflow phase volume fraction, the n-th iteration horizontal x velocity component, and the n-th iteration vertical y velocity component into the governing equation of the discretized airflow phase volume fraction to obtain the n+1-th iteration airflow phase volume fraction. Then, calculate the difference between the n-th iteration airflow phase volume fraction and the n+1-th iteration airflow phase volume fraction and use it as the residual of the airflow phase volume fraction. Step 7.3: Substitute the airflow-phase turbulent dissipation rate of the nth iteration, the airflow-phase turbulent kinetic energy of the nth iteration, the velocity component in the horizontal direction x of the nth iteration, and the velocity component in the vertical direction y of the nth iteration into the discrete turbulent kinetic energy equation to obtain the airflow-phase turbulent kinetic energy of the n+1th iteration. Calculate the difference between the airflow-phase turbulent kinetic energy of the nth iteration and the airflow-phase turbulent kinetic energy of the n+1th iteration and use it as the residual of the airflow-phase turbulent kinetic energy. Step 7.4: Substitute the airflow-phase turbulent kinetic energy of the n+1th iteration, the airflow-phase turbulent dissipation rate of the nth iteration, the velocity component in the horizontal direction x of the nth iteration, and the velocity component in the vertical direction y of the nth iteration into the discrete airflow-phase turbulent dissipation rate equation to obtain the airflow-phase turbulent dissipation rate of the n+1th iteration. Calculate the difference between the airflow-phase turbulent dissipation rate of the nth iteration and the airflow-phase turbulent dissipation rate of the n+1th iteration and use it as the residual of the airflow-phase turbulent dissipation rate. Step 7.5: Substitute the horizontal x velocity component of the nth iteration, the vertical y velocity component of the nth iteration, the airflow phase volume fraction of the n+1th iteration, the airflow phase turbulent kinetic energy of the n+1th iteration, and the airflow phase turbulent dissipation rate of the n+1th iteration into the expression of the discretized velocity component to obtain the horizontal x velocity component of the n+1th iteration and the vertical y velocity component of the n+1th iteration. Then, calculate the difference between the horizontal x velocity component of the nth iteration and the horizontal x velocity component of the n+1th iteration and use it as the residual of the horizontal x velocity component. Calculate the difference between the vertical y velocity component of the nth iteration and the vertical y velocity component of the n+1th iteration and use it as the residual of the vertical y velocity component. The residual of the airflow phase volume fraction, the residual of the airflow phase turbulent kinetic energy, the residual of the airflow phase turbulent dissipation rate, the residual of the velocity component in the horizontal direction x and the residual of the velocity component in the vertical direction y constitute a residual set; Step 7.6: Determine whether the residual set satisfies the preset conditions. If the residual set satisfies the preset conditions, use the volume fraction of the gas phase in the n+1th iteration as the final volume fraction field of the blast furnace tuyere raceway, and use the horizontal velocity component x in the n+1th iteration and the vertical velocity component y in the n+1th iteration as the final velocity field of the blast furnace tuyere raceway. If the residual set does not satisfy the preset conditions, set n = n+1 and return to step 7.

2. Among them, the residual set satisfies the preset conditions, indicating that the residual of the volume fraction of the airflow phase is less than the preset threshold, and the residual of the turbulent kinetic energy of the airflow phase is less than the preset threshold, and the residual of the turbulent dissipation rate of the airflow phase is less than the preset threshold, and the residual of the velocity component in the horizontal direction x is less than the preset threshold, and the residual of the velocity component in the vertical direction y is less than the preset threshold.

8. The method for calculating the temperature field of the blast furnace tuyere raceway based on gas-solid two-phase flow according to claim 6, characterized in that: Step 8 specifically includes: Step 8.1: For each grid in the target grid, calculate the airflow phase temperature T of the grid g Initialize and set the initial iteration number n to 0, and set the initial airflow phase temperature T g As the airflow phase temperature of the nth iteration; Step 8.2: Substitute the airflow phase temperature of the nth iteration based on all grids into the discrete temperature field mathematical model to obtain the airflow phase temperature of the n+1th iteration, and calculate the difference between the airflow phase temperature of the nth iteration and the airflow phase temperature of the n+1th iteration as the residual of the airflow phase temperature; Step 8.3: Determine whether the residual of the airflow phase temperature is less than the preset threshold. If the residual of the airflow phase temperature is less than the preset threshold, use the airflow phase temperature of the n+1th iteration as the final temperature field of the blast furnace tuyere vortex zone. If the residual of the airflow phase temperature is not less than the preset threshold, set n=n+1 and return to execute step 8.2.

Images