Background technology
Granule fluid two-phase stream extensively is present in the process industrial system.With the variation of rerum natura and operating parameter, diphasic flow presents complicated structure and parameter changes in distribution, and this has directly for pressure drop distribution, mass transfer, heat transfer efficiency and reaction conversion ratio in the industrial reactor, sometimes or even the influence of the order of magnitude.Therefore, how fluidal texture in the accurate description granule fluid two-phase stream reactor and parameter distribution are important process to design, operation, control and the amplification of process industrial.The alternate drag force Forecasting Methodology of two-phase flow then is one of core content of realizing this target.
Alternate drag coefficient commonly used is derived from effective drag coefficient that the experimental data association under fixed bed pressure drop extrapolation relation or the particulate fluidization condition forms in the two-phase flow at present, and it is the average disturbance influence of the individual particle drag coefficient having been introduced imaginary equally distributed particle swarm in essence.This account form is adopted by most researchers and practitioner, and the Fluid Mechanics Computation of being passed through in the world (CFD) software CFX, Fluent, Phoenics, STAR-CD etc. recommend (referring to list of references 1:Gidaspow, D., Jung, J., Singh, R.K., Powder Technology, Vol.148, P.123 (2004)).This method be applicable to flow to change relatively evenly, system slowly.For flowing fast of complexity, as fluid catalytic cracking common in the petroleum industry (FCC) reactor, flowing in the core cells such as (CFBB) of the Circulating Fluidized Bed Boiler in the power industry, because the agglomerating structure size of particle is indefinite, microcosmic and meso-scale boundary are fuzzy, can't find and satisfy continuous medium statistical requirements and inner even structureless infinitesimal space, so structure is inevitable to the influence of calculating.In the case, do not consider that the calculated value of structure influence compares with actual value, widely different, the result is in the different orders of magnitude, flow field granule density and reactor internal drop distribution distortion when its direct result is design or simulation calculation.Must in two-phase flow computation process, introduce the complex effects of structure, could accurately predict the flow distribution in the two-phase flow reactor, and really grasp the inside reactor rule, fundamentally ensure design, operation, control and amplify.
The variation behavior and the drag force of two-phase flow labyrinth intercouple, and are difficult to experiment measuring and model prediction, are mainly reflected in:
● non-homogeneous aggregate structure changes rapidly, does not have the geometric shape of determining, its geometry is difficult to characterize with mechanics parameter, is difficult to directly measure the drag force variation under the non-homogeneous dynamic structure;
● the rerum natura under the alternate exchange intensity of drag coefficient representative and a plurality of yardsticks, operating parameter are relevant: it is not only relevant with average two-phase superficial velocity, the concentration of reactor, and is also closely related with the variation of local velocity and structure.This association of striding yardstick is difficult to dwindling survey region earlier, extends to macroscopical method from the microcosmic localized variation again and obtains.
Previous work pays attention to distinguish the force difference of aggregate particle and individual particle, and the sight aggregate structure employing that is situated between is simplified, and static often, the bulky grain of shape invariance replaces.Existing a small amount of exploration work adopts the method for empirical correlation or experience to obtain the characteristic dimension of aggregate, the dynamic structure of particle cluster is expressed as the potpourri of a kind of static bulky grain and primary granule, and calculate alternate drag force based on this (referring to list of references 2:Heynderickx, G.J., Das, A.K., De Wilde, J., Marin, G.B., Industrial ﹠amp; EngineeringChemistry Research, Vol.43, P.4635, (2004)).This static state, experimental method of estimation can't be applied to the industrial reactor from original experiment condition extrapolation in essence.
Existing work only noted the substructure micro-scale to the influence of mesoscopic structure (referring to list of references 3:Agrawal, K., Loezos, P.N., Syamlal, M., Sundaresan, S.Journal of Fluid Mechanics.Vol.445, P.151. (2001)), attempt to obtain of the influence of complicated mesoscopic structure, and do not have to consider to stride the influence of yardstick association mesoscopic structure to drag force from the result of calculation direct correlation in sample infinitesimal space.The sample space of real process is difficult to choose, and is infeasible in the application of industrial reactor.
Summary of the invention
The objective of the invention is does not have to consider to stride the influence of yardstick association to mesoscopic structure in order to overcome prior art, result of calculation and actual conditions that flow parameter distributes differ shortcoming far away, thus the method that provides a kind of flow parameter that can accurately measure in the particle two-phase flow reactor to distribute.
To achieve these goals, the invention provides a kind of method that flow parameter distributes in the granule fluid two-phase stream reactor of measuring, comprising:
1), according to the rerum natura of reality and the operating conditions situation of reactor, initialization flow field and boundary condition;
2), utilize two-fluid model Accounting Legend Code or commercialization cfdrc, calculate mass-conservation equation, momentum conservation equation in each space infinitesimal, obtain speed, CONCENTRATION DISTRIBUTION (u
g, u
p, ε
g), described u
gRepresent the gas true velocity in the infinitesimal, described u
pRepresent the solid true velocity in the infinitesimal, described ε
gRepresent average void fraction in the infinitesimal;
3), judge whether current speed and CONCENTRATION DISTRIBUTION satisfy the convergence of whole flow field, if satisfy, export current speed and CONCENTRATION DISTRIBUTION, and end operation, otherwise, carry out next step;
4), according to known reactor monolith operating conditions (U
g, U
p), the combination of calculating and being met all variable roots of overall dynamic (dynamical) unstable state, Nonlinear System of Equations; Described U
gRepresent the gas phase superficial velocity; Described U
pRepresent particle phase superficial velocity;
5), from the combination of all variable roots, seek the optimum root satisfy extremum conditions, and preserve mesoscopic structure parameter ε wherein
c, ε
f, d
Cl, described ε
cThe close phase voidage of expression aggregate, described ε
fExpression dilute phase voidage, described d
ClExpression aggregate diameter;
6), by step 2) speed, the CONCENTRATION DISTRIBUTION (u that obtain
g, u
p, ε
g) and the mesoscopic structure parameter ε that obtains of step 5)
c, ε
f, d
Cl, in each infinitesimal space of reactor, find the solution the combination of all variable roots that satisfy the local dynamics unstable state of microcosmic, Nonlinear System of Equations; Described variable root comprises U
f, U
Pf, U
c, U
Pc, f, a
i, a
c, a
f, U wherein
fRepresent the dilute phase fluid velocity, U
PfRepresent the dilute phase particle speed, U
cDelegation's polymers dense-phase fluid speed, U
PcThe close phase particle speed of delegation's polymers, f represents close phase volume fraction, a
iRepresentative is situated between and sees interface acceleration, a
cRepresent close phase acceleration, a
fRepresent the dilute phase acceleration;
7), the combination of the variable root that obtains by step 6), therefrom seek the optimum root that satisfies extremum conditions;
8), by all variate-value (U in step 5) and the resulting infinitesimal of step 7)
f, U
Pf, U
c, U
Pc, ε
c, ε
f, f, d
Cl, a
i, a
c, a
f) combination, calculate drag force and drag coefficient;
9), drag coefficient substitution step 2 that step 8) is obtained) in the two-fluid model Accounting Legend Code or commercialization cfdrc that are adopted, replace original drag coefficient, forward step 2 then to), do iterative computation.
In the technique scheme, described unstable state, Nonlinear System of Equations comprise:
In the unit volume infinitesimal space, concentrated phase aggregate particle drag force F
D, cBalance equation:
In the unit volume infinitesimal space, the dilute phase particle is subjected to drag force F
D, fBalance equation:
Be situated between in the unit volume infinitesimal space and see the drag force F of phase interface
D, iBalance equation:
Fluid and particle mass-conservation equation mutually:
U
g=fU
c+(1-f)U
f
U
p=fU
pc+(1-f)U
pf
ε
g=fε
c+(1-f)ε
f
The aggregate diameter d
ClFormula:
Wherein, ρ represents density; d
pExpression individual particle diameter; G represents acceleration of gravity; Subscript g represents gas phase; Subscript p represents solid phase; Subscript c represents concentrated phase; Subscript f represents dilute phase; C
D, cRepresent close phase average individual particle drag coefficient; C
D, fThe average individual particle drag coefficient of expression dilute phase; C
D, iExpression is situated between and sees alternate average single aggregate drag coefficient;
The definition of above-mentioned parameter is as follows:
Re
f=d
pρ
gU
sf/μ
g
Re
c=d
pρ
gU
sc/μ
g
Re
i=d
clρ
gU
si/μ
g
Each apparent sliding velocity is defined as:
U
sc=U
c-ε
cU
pc/(1-ε
c)
U
sf=U
f-ε
fU
pf/(1-ε
f)
In the technique scheme, in described step 5) and the step 7), described extremum conditions is meant that the particle suspending conveying capacity of unit mass is tending towards minimum, the particle suspending conveying capacity N of described unit mass
StExpression, it is defined as:
In the technique scheme, in the described step 8), described drag force F
dExpression, its computing formula is: F
d=fF
D, c+ (1-f) F
D, f+ F
D, iDescribed drag coefficient is represented with β, can be obtained by drag force, and the computing formula of described drag coefficient is:
In the technique scheme, in described step 2) in, described commercialization cfdrc adopts CFX or Fluent or Phoenics or STAR-CD.
In the technique scheme, described granule fluid two-phase stream reactor is fluidized calcinator or coal-powder boiler or circulating fluid bed reactor or gas-solid and flows to upward riser reactor.
The invention has the advantages that:
1, structural parameters really are coupled at infinitesimal level and Fluid Mechanics Computation, have obtained the alternate momentum exchange coefficient or the drag coefficient of reflection actual conditions, have corrected the higher situation of drag coefficient in the conventional method;
2, obtained the interior structure distribution of any infinitesimal, reaction, heat transfer, mass transport process on this basis can correct calculation granule fluid two-phase stream reactor, accurately realize the hot calculating of reaction of two-phase flow reactor under the industrial condition, thereby instruct correlated process control, optimal design and even technology to amplify.
Embodiment
Below in conjunction with the drawings and specific embodiments the method that flow parameter in the measurement granule fluid two-phase stream reactor of the present invention distributes is elaborated.
The granule fluid two-phase stream reactor is fluidized calcinator, coal-powder boiler, circulating fluid bed reactor and gas-solid and flows to upward riser reactor.Method of the present invention is a kind of computing method in conjunction with extreme value multi-Scale Structural Model and two-fluid model, before the method for the invention is described, at first extreme value multi-Scale Structural Model and two-fluid model is made an explanation.
Extreme value multi-Scale Structural Model in the inventive method is meant: model comes computation structure parameter and flow parameter by the physical quantity on a plurality of characteristic time scales and feature space yardstick, carries out related by conservation equation with extremum conditions between the physical quantity.Wherein, extreme value multi-Scale Structural Model related dynamics conservation equation on all yardsticks adopts two-fluid model again.
All have non-homogeneous multiple dimensioned structure in any infinitesimal control volume of granule fluid two-phase stream, this structure can equivalence be aggregate concentrated phase and two kinds of two-fluid media of discrete particle dilute phase, and the interfacial structure and the form of depositing.Momentum conservation equation separately can be write as the form of dynamic balance through transforming.Being without loss of generality, is example with the gas particles two-phase flow:
In the unit volume infinitesimal space, concentrated phase aggregate particle drag force F
DcBalance equation:
In the unit volume infinitesimal space, the dilute phase particle is subjected to drag force F
DfBalance equation:
Particle is stressed different with its inside at the interface owing to the concentrated phase aggregate, can obtain being situated between in the unit volume infinitesimal space seeing the drag force F of phase interface
DiBalance equation:
Fluid and particle mass-conservation equation mutually:
U
g=fU
c+(1-f)U
f (4)
U
p=fU
pc+(1-f)U
pf (5)
ε
g=fε
c+(1-f)ε
f (6)
Wherein, ρ, density; U
g, the gas phase superficial velocity; U
p, particle phase superficial velocity; U
f, the dilute phase fluid velocity; U
Pf, the dilute phase particle speed; U
c, aggregate dense-phase fluid speed; U
Pc, the close phase particle speed of aggregate; ε
g, average void fraction in the infinitesimal; ε
c, the close phase voidage of aggregate; ε
f, the dilute phase voidage; F, close phase volume fraction; d
Cl, the aggregate diameter; d
p, the individual particle diameter; a
i, be situated between and see the interface acceleration; a
c, close phase acceleration; a
f, the dilute phase acceleration; G, acceleration of gravity; Subscript g, gas phase; Subscript p, solid phase; Subscript c, concentrated phase; Subscript f, dilute phase.C
Dc, close phase average individual particle drag coefficient; C
Df, the average individual particle drag coefficient of dilute phase; C
Di, be situated between and see alternate average single aggregate drag coefficient.Definition is as follows:
Re
f=d
pρ
gU
sf/μ
g
Re
c=d
pρ
gU
sc/μ
g
Each apparent sliding velocity is defined as:
U
sc=U
c-ε
cU
pc/(1-ε
c)
U
sf=U
f-ε
fU
pf/(1-ε
f)
And the aggregate diameter d
ClAs give a definition:
Wherein, subscript m f represents ε
g=ε
MfThe time the minimum fluidized state of broad sense, ε
MfVoidage during for minimum fluidized state, subscript m ax represents ε
g=ε
MaxThe time sparse ultimate limit state, ε
MaxFor having the maximum interspace rate of particle agglomeration, the particle suspending conveying capacity N of unit mass
StBe defined as follows:
Above-mentioned formula (1)-(7) are the dynamics conservation equation, include 11 undetermined parameter: U in these equations
f, U
Pf, U
c, U
Pc, ε
c, ε
f, f, d
Cl, a
i, a
c, a
fSimple dynamics conservation equation is not enough to the structural factor in definite complicated diphasic flow, also needs the extremum conditions in the extreme value multi-Scale Structural Model.
By the coordinating analysis to particle and fluid two-phase movement mechanism, the trend of finding fluid motion is to select the path of resistance minimum, and the motion of particle is to be tending towards potential energy minimum or poly-group in the scope of Boundary Condition Control.Two kinds of locomotory mechanisms have coordinated to produce the extremum conditions of extreme value multi-Scale Structural Model mutually, that is: the particle suspending conveying capacity of unit mass is tending towards minimum (N
St→ min).
Like this, in the extreme value multi-Scale Structural Model 11 flow and structural parameters can be able to mathematics by the extreme-value problem that 7 conservation equations and 1 extremum conditions are formed and find the solution.That is: at given rerum natura and flow parameter (U
g, U
p) etc. under the condition, find the solution satisfy kinetics equation (1)-(7) satisfy the root of Nst → min in might root simultaneously, mathematical expression is
Min?N
st(X)s.t.F
i(X)=0,i=1,2,……7. (9)
As shown in Figure 1, the method that flow parameter distributes in the measurement granule fluid two-phase stream reactor of the present invention may further comprise the steps:
Step 10, according to the rerum natura of reality and the operating conditions situation of reactor, initialization flow field and boundary condition.Described initialization flow field and boundary condition are exactly will be according to the rerum natura of reality and physical dimension, grain diameter, particle density, physical properties of fluids, temperature of reaction, apparent fluid velocity and the particle flow flux etc. of the given granule fluid two-phase stream reactor of operating conditions.
Step 20, utilize two-fluid model Accounting Legend Code or commercialization cfdrc, calculate mass-conservation equation, momentum conservation equation in each space infinitesimal, obtain speed, CONCENTRATION DISTRIBUTION.The described cfdrc of this step can adopt CFX, F1uent, Phoenics, STAR-CD etc.Speed and CONCENTRATION DISTRIBUTION can be used (u
g, u
p, ε
g) expression.u
gAnd u
pRepresent gas and solid true velocity in the infinitesimal.
Step 30, judge whether current speed and CONCENTRATION DISTRIBUTION satisfy the convergence of the whole audience,, export current speed and CONCENTRATION DISTRIBUTION if satisfy the convergence of the whole audience, and end operation; Otherwise, carry out next step.
Step 40, according to known reactor monolith operating conditions (U
g, U
p), the combination of calculating and being met all variable roots of overall dynamic (dynamical) unstable state, Nonlinear System of Equations.In this step, unstable state, Nonlinear System of Equations are meant kinetics equation (1)-(7), in these 7 equations, 11 variablees are arranged, and obtain variable root (U by 7 equations
f, U
Pf, U
c, U
Pc, ε
c, ε
f, f, d
Cl, a
i, a
c, a
f) combination.
Step 50, from the combination of all variable roots, seek the optimum root satisfy extremum conditions, and preserve mesoscopic structure parameter ε wherein
c, ε
f, d
ClIn the present invention, extremum conditions is meant that the particle suspending conveying capacity of unit mass is tending towards minimum, i.e. N
St→ min.
Step 60, the speed, the CONCENTRATION DISTRIBUTION (u that obtain by step 20
g, u
p, ε
g) and the mesoscopic structure parameter ε that obtains of step 50
c, ε
f, d
Cl, in each infinitesimal space of reactor, find the solution all variable root (U that satisfy the local dynamics unstable state of microcosmic, Nonlinear System of Equations (1)-(7)
f, U
Pf, U
c, U
Pc, f, a
i, a
c, a
f) combination.
The optimum root that satisfies extremum conditions is therefrom sought in the combination of step 70, the variable root that obtained by step 60.Described extremum conditions is meant that the particle suspending conveying capacity of unit mass is tending towards minimum, i.e. N
St→ min
Step 80, by all variate-value (U in step 50 and the resulting infinitesimal of step 70
f, U
Pf, U
c, U
Pc, ε
c, ε
f, f, d
Cl, a
i, a
c, a
f) combination, calculate drag force and drag coefficient.The computing formula of drag force as shown in Equation (10).
F
d=fF
d,c+(1-f)F
d,f+F
d,i (10)
According to resulting drag force, also can further obtain drag coefficient β, the relation between drag force and drag coefficient is as shown in Equation (11).
The drag coefficient β of non-structure influence
0Expression
In the two-fluid model Accounting Legend Code or commercialization cfdrc that is adopted in step 90, the drag coefficient substitution step 20 that step 80 is obtained, replace original drag coefficient, forward step 20 then to, do iterative computation.
In the infinitesimal that the inventive method calculates, on the basis of dense dilute phase structure and flow parameter, can calculate dense, dilute phase rate of heat transfer h, mass transfer velocity k and reaction rate r separately in the infinitesimal, and then sum up.Because heat transfer, mass transfer and reaction are non-linear mostly, the result of calculation with structural parameters is realistic, and homogeneous model does not conform to the actual conditions.Such as for reaction,
r
Actual=r (concentrated phase parameter)+r (dilute phase parameter)+r (interface parameter) ≠ r (mean parameter)
With the inventive method be coupled to quality, energy conservation equation mutually in, can further obtain the mass concentration in the granule fluid two-phase stream reactor, the information of heat distribution, thus the more accurate hot calculating of reaction that must realize the two-phase flow reactor.
As shown in Figure 2, the drag coefficient in the true drag coefficient that has comprised structural factor and the commercial software under the homogeneous condition is different, and heterogeneous texture is drag forces when evenly flowing all in most of the cases.
The method that flow parameter distributes in the measurement granule fluid two-phase stream reactor of the present invention has a good application industrial.
Embodiment 1 is applied to produce per year 1000 tons of mobile predictions in the Coaseries kaolin fast fluidization calcinator with the inventive method.Feed kaolin mean grain size 0.002mm, particle density 800kg/m
3, calcining burner hearth internal diameter 0.45m, the high 10.5m of burner hearth.Fire box temperature is controlled at about 900 ℃.The key reaction that takes place is: C+O
2→ CO
2, Al
2O
32SiO
2H
2O → Al
2O
32SiO
2On the general-purpose platform that commercialization CFD software provides, the drag coefficient that the inventive method is calculated according to computing block diagram is compiled into interface routine to replace the method for software inhouse during calculating.At first reaction unit is carried out how much structure bodies of 1:1, the selection account form is the two-fluid model under two Eulerian coordinates; Secondly design temperature and set particle, gas property according to actual rerum natura system; Material is initially set to the nature stacking states in the burner hearth, and piling up voidage is 0.6, and gas enters from the burner hearth bottom inlet, superficial gas velocity 1.5m/s, and burner hearth middle part feedback outlet is discharged and be circulated to granule materials by the top-side outlet.To calculate that interface document is called in and compilation run.Calculating is calculated 20 seconds actual flow time according to carrying out with the dynamic similation mode of time correlation.The voidage distribution that this method obtains, solid Flux Distribution, solid speed distribute and have dynamic heterogeneous texture, realistic mode of operation; And the even flow field that commercial software calculates does not have tangible heterogeneous texture, and this and actual conditions are not inconsistent.For actual bed mean concentration is 5% situation, and the burner hearth total pressure drop (3557Pa) that the inventive method is calculated is compared more approaching actual test value (3058Pa) with the result of calculation (5608Pa) of commercial software.
Embodiment 2 is applied to the inventive method per hour mobile prediction in the wuhan petrochemical industry factory Circulating Fluidized Bed Boiler of 75 tons of producing steams with the auxiliary detection fault.As shown in Figure 3, boiler furnace principal section 5.4 * 3.99m, the high about 22.5m of burner hearth, circulating ash particle mean grain size 0.21mm in the burner hearth, particle density 2000kg/m
3, fire box temperature is controlled at about 900 ℃, design primary air flow 49701Nm
3/ h, secondary air flow 33134Nm
3/ h.Main combustion reaction: the C+O that takes place in the burner hearth
2→ CO
2On the general-purpose platform that commercialization CFD software provides, write the drag coefficient that the inventive method is calculated the external interface program of User Defined form as to replace the method for software inhouse.At first boiler furnace is carried out how much structure bodies of 1:1 during calculating, the selection account form is the two-fluid model under two Eulerian coordinates; Secondly design temperature and set particle, gas property according to actual rerum natura system; Material is initially set to the nature stacking states in the burner hearth, piling up voidage is 0.6, the burner hearth bottom superficial gas velocity calculates according to the primary air flow of design, secondary air realizes that by add the quality source item in the gas phase conservation equation of respective heights burner hearth side feedback outlet is discharged and be circulated to granule materials by the top-side outlet.Call in the calculating interface document and the compilation run of drag coefficient.Calculate 60 seconds actual flow time.Fig. 4 has shown the pressure in the burner hearth distribution plan that the inventive method calculates, each point is the force value of 40 equally distributed radial positions on the cross section of burner hearth corresponding height among the figure, all point overlaps to any illustrates and does not exist radial pressure to distribute on this height, 40 radial position data are overlapping to be that one section line illustrates and radially has pressure distribution on this height, and flowing, it is violent, unstable to mix.The burner hearth ash circulating ratio that the present invention dopes is 11, and factory's operation site is tested the empirical value of summary for a long time about 10, and commercial software result of calculation is 46.The inventive method is better than present commercial software.
Embodiment 3 is applied to the prediction of China Petrochemical Industry's novel oil refining multistage reactor flow field distribution with auxiliary optimal design with the inventive method.Reactor body section diameter 3.4m, inlet gas speed 15m/s, import amount of solid 1500t/h, particles used mean grain size 0.075mm, particle density 1300kg/m
3, about 400 ℃ of burner hearth operating temperature.Write the drag coefficient that the inventive method is calculated as the external interface program to replace the method for commercial software inside.At first reactor is carried out how much structure bodies of 1:1 during calculating, the selection account form is the two-fluid model under two Eulerian coordinates; Next design temperature, and according to actual rerum natura system setting particle, gas property; Bottom superficial gas velocity and solid speed are given according to design load, and granule materials is discharged by top exit.Calculating is calculated 120 seconds actual flow time according to carrying out with the dynamic similation mode of time correlation.Fig. 5 shown the influence that distribution grid that the inventive method calculates distributes to granule density (▲-distribution grid arranged; Zero-distribution-free plate).The result shows that its interior distribution grid of reflection influences material distribution form in the reactor, but little to extension diameter section overall reactor mean concentration (remaining on about 2.4%) influence.This conforms to practice in factory.
Embodiment 4 is applied to pilot-scale circulating fluidization reducing reactor flow field with the inventive method and distributes with the amplification design industrial-scale reactor.Reactor epimere diameter 0.09m, hypomere diameter 0.045m, the high about 6.2m of epimere, particle mean grain size 0.082mm, particle density 1450kg/m
3, reactor bottom inlet gas speed 10m/s, solid flux 95kg/m
2S, side import solid flux 12.6kg/m
2S.Write the drag coefficient that the inventive method is calculated the method for external interface program as with replacement commercial software inside, and compilation run.Fig. 6 shows that the predicting the outcome of the inventive method (●) and experiment value (■) coincide finely, and that existing business software calculated value (▲) predicts the outcome is rare partially, departs from bigger with experiment value.