CN109783904A - A kind of width parameter area carbon dioxide physical property method for solving - Google Patents
A kind of width parameter area carbon dioxide physical property method for solving Download PDFInfo
- Publication number
- CN109783904A CN109783904A CN201811626329.3A CN201811626329A CN109783904A CN 109783904 A CN109783904 A CN 109783904A CN 201811626329 A CN201811626329 A CN 201811626329A CN 109783904 A CN109783904 A CN 109783904A
- Authority
- CN
- China
- Prior art keywords
- area
- enthalpy
- pressure
- region
- carbon dioxide
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
Abstract
The invention discloses a kind of wide parameter area carbon dioxide physical property method for solving, are iterated solution using actual gas state equation in the region of carbon dioxide critical point；Physical property is carried out using fit correlation formula in the region far from carbon dioxide critical point and carries out direct solution；Interpolation solution is carried out using fit correlation formula and actual gas state equation in borderline region, saturated solution enthalpy, saturated air enthalpy, temperature, specific volume, thermal conductivity and dynamic viscosity physical parameter is calculated.The present invention had not only guaranteed the computational accuracy of carbon dioxide, but also promoted the solving speed of reactor program of numerical calculation.
Description
Technical field
The invention belongs to nuclear reactor safety analytical calculation technical fields, and in particular to a kind of width parameter area carbon dioxide
Physical property method for solving.
Background technique
Supercritical carbon dioxide Brayton cycle has many advantages, such as that the thermal efficiency is high, system simplification is compact, is to have very much prospect
Forth generation nuclear reactor energy conversion system.When supercritical carbon dioxide Brayton cycle steady-state operation, at suction port of compressor
Supercritical carbon dioxide be in Near The Critical Point, carbon dioxide will undergo violent and rapid physical property to change, carbon dioxide
Specific volume, specific heat, the key physical parameter such as thermal conductivity can be distorted.In addition, in transient state or accident, supercritical carbon dioxide
The operating status point of Brayton cycle will undergo large-scale physical property to change, therefore, accurate carbon dioxide Calculation of Physical Properties method,
For realizing bistable design, transient control strategy and the crash analysis of supercritical carbon dioxide Brayton cycle reactor,
With very important use value.
The calculation method of thermodynamic fluid physical property mainly has diagram method, state Equation Method and fit correlation formula method three at present
Class.Diagram method is the thermal transport cofficients for calculating fluid manufactured figure and table, and this method advantage is to be easily understood but lack
Point is that computational efficiency is very low, is not suitable for large-scale thermodynamic analysis program.State Equation Method is the root of thermodynamic fluid Calculation of Physical Properties
This method, this method be established based on stringent theory and experimental study, therefore such methods have the advantages that it is with high accuracy, but
It is the disadvantage is that model is complicated, calculating the time, long and programming is complicated.Fit correlation formula method is that simplifying for thermodynamic fluid physical property calculates
Method, this method are satisfied the relationship of certain precision by being fitted to existing thermal transport cofficients data.
Since the relational expression of use is generally the simple mathematicals relational expression such as multinomial, fraction, therefore this method has and calculates simple, calculation amount
It is small, computational accuracy is higher, be suitable for computer programming solve the advantages of, largely used in Thermodynamic Calculation Program.
The current main adoption status equation method of carbon dioxide thermodynamic analysis program carries out Calculation of Physical Properties, and this method can be counted
Calculation obtains the higher carbon dioxide physical property of precision.But state Equation Method needs a large amount of iteration to obtain physical property as a result, when calculating
Between it is long.Calculation of Physical Properties is extremely important in Thermodynamic Calculation Program, and especially in transient state Thermodynamic Calculation Program, Calculation of Physical Properties is occupied
Program calculates 20% or so of time.If Calculation of Physical Properties speed is slow, will lead to the calculating speed of transient state Thermodynamic Calculation Program at
It reduces again.Physical property is calculated with fit correlation formula can effectively improve computational efficiency.But carbon dioxide is in the object of Near The Critical Point
Sex distortion causes the computational accuracy of fit correlation formula insufficient.
The physical property feature of carbon dioxide leads to the existing transient state that can be used for supercritical carbon dioxide reactor Brayton cycle
It is unsatisfactory for requiring in computational accuracy or calculating speed with accident numerical analysis programs.It is wide urgently to develop a kind of quick and precisely solution
The method of parameter area carbon dioxide physical property promotes the computational accuracy and calculating speed of transient state Thermodynamic Calculation Program.
Summary of the invention
In view of the above-mentioned deficiencies in the prior art, the technical problem to be solved by the present invention is that providing a kind of wide parameter model
Carbon dioxide physical property method for solving is enclosed, had not only guaranteed the computational accuracy of carbon dioxide, but also promote asking for reactor program of numerical calculation
Solve speed.
The invention adopts the following technical scheme:
A kind of width parameter area carbon dioxide physical property method for solving uses actual gas in the region of carbon dioxide critical point
State equation is iterated solution；Physical property solution is carried out using fit correlation formula in the region far from carbon dioxide critical point；?
Borderline region carries out interpolation solution using fit correlation formula and actual gas state equation, and saturated solution enthalpy, saturation is calculated
Gas enthalpy, temperature, specific volume, thermal conductivity and dynamic viscosity physical parameter.
Specifically, it is as follows to solve saturation liquid phase enthalpy using fitting of a polynomial relational expression relevant to pressure:
Wherein, hl1It (p) is the saturation liquid phase enthalpy of low-pressure area, hl2It (p) is the saturation liquid phase enthalpy of higher-pressure region, hl3(p) it is
The saturation liquid phase enthalpy in interpolation area, p are pressure, and MPa, a, b, i are constant coefficient.
Specifically, it is as follows to solve saturation gas phase enthalpy using fitting of a polynomial relational expression relevant to pressure:
Wherein, hg1It (p) is saturation gas phase the enthalpy g, h of low-pressure areag2It (p) is the saturation gas phase enthalpy of higher-pressure region, hg3(p)
For the saturation gas phase enthalpy in interpolation area, p is pressure, and MPa, a, b, i are constant coefficient.
Specifically, by the zoning of temperature be divided into cold-zone, overheat one area and overheat 2nd area, pressure be 0.1~
20MPa, specific enthalpy are 200~1600kJkg-1。
Further, the range for crossing cold-zone is enthalpy of saturated liquid line and critical specific enthalpy line lower boundary (312.25kJkg-1)
Under, region of the pressure from 0.1MPa to 20MPa；The range for overheating an area is saturated air specific enthalpy line and critical specific enthalpy line coboundary
(352.25kJ·kg-1) more than, 600kJkg-1Specific enthalpy line is hereinafter, region of the pressure from 0.1MPa to 20MPa；Overheat 2nd area
Range be 600kJkg-1It is more than specific enthalpy line, 1600kJkg-1Specific enthalpy line is hereinafter, area of the pressure from 0.1MPa to 20MPa
Domain.
Further, the temperature T of cold-zone is crossed1(p, h) calculates as follows:
Wherein, aijFor constant coefficient, p is pressure, and MPa, h are enthalpy, kJ/kg, p≤pcritical, h < hl(p) or p >
pcritical, h≤hset10, hset10=hcritical- 20.0=312.25kJ/kg；
Overheat the temperature T in an area2(p, h) calculates as follows:
Wherein, bijFor constant coefficient, p≤pcritical, hg(p) < h≤hset20Or p > pcritical, hset11< h≤hset20,
hset11=hcritical+ 20.0=352.25kJ/kg, hset20For constant ,=580kJ/kg；
Overheat the temperature T in 2nd area3(p, h) calculates as follows:
Wherein, cijFor constant coefficient, h > hset21, hset21For constant ,=620kJ/kg；
One area of interpolation was the transitional region in one area of cold-zone and overheat, temperature T4(p, h) calculates as follows:
Wherein, h is enthalpy, p > pcritical, hset10< h≤hset11；
2nd area of interpolation is the transitional region for overheating an area and overheating 2nd area, temperature T5(p, h) calculates as follows:
Specifically, the zoning of specific volume was divided into cold-zone, 2nd area of state equation area, one area of overheat and overheat, pressure
It is 200~1600kJkg for 0.1~20MPa, specific enthalpy-1。
Further, the range for crossing cold-zone is enthalpy of saturated liquid line and 260kJkg-1Under specific enthalpy line, pressure from
The region of 0.1MPa to 20MPa；State equation area uses actual gas Solving Equation of State, and range includes two parts, and first
Dividing is enthalpy of saturated liquid line or less and 296kJkg-1More than specific enthalpy line, more than saturated air line and 380kJkg-1Specific enthalpy line with
Under, region of the pressure from 7MPa to 7.38MPa, second part is 296kJkg-1More than specific enthalpy line and 380kJkg-1Specific enthalpy
Line hereinafter, pressure from 7.38MPa to 20MPa；The range for overheating an area is saturated air line or more and 580kJkg-1Specific enthalpy line with
Under, region of the pressure from 0.1MPa to 7MPa；Enthalpy is in 410kJkg-1More than specific enthalpy line and 580kJkg-1Specific enthalpy line with
Under, region of the pressure from 7MPa to 20MPa；The range for overheating 2nd area is 600kJkg-1It is more than specific enthalpy line, 1600kJkg-1
Specific enthalpy line is hereinafter, region of the pressure from 0.1MPa to 20MPa.
Further, the specific volume v of cold-zone is crossed1(p, h) calculates as follows:
Wherein, aijFor constant coefficient, p≤7MPa, h < hl(p) or p > 7MPa, h≤hset30, hset30For constant ,=
260kJ/kg；
Hereby function does not calculate using Hall for state equation area:
φ (δ, τ)=φo(δ,τ)+φr(δ,τ)
Wherein, φ is the Hall not hereby energy of real gas, φoFor the Hall not hereby perfect gas part of energy, φrFor
The Hall not hereby remainder of energy, δ=ρ/ρc, τ=T/Tc, ρc=467.6kg/m3, Tc=31.1 DEG C, p >=pcritical,
hset00≤h≤hset01Or 7MPa≤p≤pcritical, hset00≤h≤hl(p) or 7MPa≤p≤pcritical, hg(p)≤h≤
hset01, hset00For constant ,=296kJ/kg, hset01For constant ,=380kJ/kg, pcriticalFor the critical pressure of carbon dioxide,
=7.31MPa；
Overheat the specific volume v in an area3(p, h) calculates as follows:
Wherein, bijFor constant coefficient, p≤7MPa, hg(p) < h≤hset20Or p > 7MPa, hset31< h≤hset20, hset31For
Constant ,=410kJ/kg；
Overheat the specific volume v in 2nd area4(p, h) calculates as follows:
Wherein, cijFor constant coefficient, h > hset21；
One area of interpolation was the transition region of cold-zone and state equation area, specific volume v5(p, h) calculates as follows:
Wherein, hset30=260kJ/kg, hset00=296kJ/kg, p > 7MPa, hset30< h≤hset00；
2nd area of interpolation is the transition region in one area of state equation area and overheat, specific volume v6(p, h) calculates as follows:
Wherein, hset31=410kJ/kg, hset01=380kJ/kg, p > 7MPa, hset01< h≤hset31；
3rd area of interpolation is the transition region for overheating an area and overheating 2nd area, specific volume v7(p, h) calculates as follows:
Wherein, hset20=580kJ/kg, hset21=620kJ/kg, hset20< h≤hset21。
Specifically, the computer capacity of thermal conductivity, in 216~1100K of temperature, pressure to 200MPa calculates as follows:
λ (ρ, T)=λ0(T)+Δλ(ρ,T)+Δλc(ρ,T)
Wherein, λ0(T) thermal conductivity in zero density low density gas is indicated, Δ λ (ρ, T) represents thin gas at identical temperature
Volume density increases bring thermal conductivity and increases, Δ λc(ρ, T) represents the increase in critical point peak region thermal conductivity；
The computer capacity of dynamic viscosity is in 200~1500K of temperature, density to 1400kg/m3, it is specific as follows:
η (ρ, T)=η0(T)+Δη(ρ,T)+Δηc(ρ,T)
Wherein, η0(T) viscosity in zero density low density gas is indicated, Δ η (ρ, T) represents low density gas at identical temperature
Density increases bring viscosity and increases, Δ ηc(ρ, T) represents the increase in critical point peak region viscosity.
Compared with prior art, the present invention at least has the advantages that
The one wide parameter area carbon dioxide physical property method for solving of the present invention, polynomial fitting is mutually tied with state equation
It closes, can rapidly and accurately calculate the physical property such as saturated solution enthalpy, saturated air enthalpy, temperature, specific volume, thermal conductivity and dynamic viscosity ginseng
Number.
Further, enthalpy of saturated liquid is calculated, in 0.6~7.377MPa, takes 6777 pressure by step-length of 1kPa
Point carries out enthalpy of saturated liquid calculating, and it is 0.08% that method of the invention, which calculates maximum relative error, and calculating speed is business software
118 times of NIST PREPROP.
Further, saturated air specific enthalpy is calculated, in 0.6~7.377MPa, takes 6777 pressure by step-length of 1kPa
Point carries out enthalpy of saturated liquid calculating, and it is 0.009% that method of the invention, which calculates maximum relative error, and calculating speed is business software
230 times of NIST PREPROP.
Further, for temperature computation, in 0.1~20MPa of pressure limit, step-length 0.1MPa, than enthalpy spectrum 200-
1600kJ·kg-1, step-length 2kJkg-1, take 131914 state points to be calculated, it is opposite accidentally that method of the invention calculates maximum
Difference is 0.1%, and calculating speed is 15 times of business software NIST PREPROP.
Further, specific volume is calculated, in 0.1~20MPa of pressure limit, step-length 0.1MPa, than enthalpy spectrum 200-
1600kJ·kg-1, step-length 2kJkg-1, take 131914 state points to be calculated, it is opposite accidentally that method of the invention calculates maximum
For difference within 1%, calculating speed is 28 times of business software NIST PREPROP.
Further, for the calculating of thermal conductivity, for maximum relative error within 5%, calculating speed is business software
NIST PREPROP is suitable.
Further, for the calculating of dynamic viscosity, for maximum relative error within 5%, calculating speed is business software
NIST PREPROP is suitable.
In conclusion the present invention is based on the carbon dioxide objects reported in the polynomial fitting of independent development and open source literature
Property calculation relational expression and actual gas state equation, develop a carbon dioxide physical property method for solving, are used for and carbon dioxide
The physical property of relevant reactor numerical analysis programs solves.Physical property solution proposes polynomial fitting and state equation combines
Physical property method for solving, in the parameters such as enthalpy of saturated liquid, saturated air specific enthalpy and temperature, using the fit correlation formula of independent research
Direct solution is carried out, in the parameters such as thermal conductivity and dynamic viscosity, the calculation formula delivered using open source literature is directly asked
Solution, improves the calculating speed of physical property packet.In specific volume parameter, in Near The Critical Point, carbon dioxide state equation is mainly used
It is iterated solution, ensure that the precision that specific volume calculates, when far from critical point, is carried out using the fit correlation formula of independent research
Direct solution.
Below by drawings and examples, technical scheme of the present invention will be described in further detail.
Detailed description of the invention
Fig. 1 is temperature computation schematic diagram；
Fig. 2 is that specific volume calculates schematic diagram；
Fig. 3 is enthalpy of saturated liquid schematic diagram of calculation result, wherein (a) is that enthalpy of saturated liquid low-pressure area relational expression is opposite accidentally
Difference (b) is enthalpy of saturated liquid higher-pressure region relational expression relative error；
Fig. 4 is that enthalpy of saturated liquid fit correlation formula and business software NIST REFPROP calculating speed compare；
Fig. 5 is saturated air specific enthalpy schematic diagram, wherein (a) is saturated air specific enthalpy low-pressure area relational expression relative error, (b) is
Saturated air specific enthalpy higher-pressure region relational expression relative error；
Fig. 6 is that saturated air specific enthalpy fit correlation formula and business software NIST REFPROP calculating speed compare；
Fig. 7 is temperature computation result schematic diagram, wherein (a) was cold-zone temperature dependence relative error, was (b) overheat
One area's temperature dependence relative error is (c) two area's temperature dependence relative errors of overheat；
Fig. 8 is that temperature foh relational expression and business software NIST REFPROP calculating speed compare；
Fig. 9 is specific volume schematic diagram of calculation result, wherein (a) was cold-zone specific volume relational expression relative error, was (b) overheat
Two area's temperature dependence relative errors；
Figure 10 is that specific volume fit correlation formula and business software NIST REFPROP calculating speed compare.
Specific embodiment
The present invention provides a kind of wide parameter area carbon dioxide physical property method for solving, in the region of carbon dioxide critical point
Physical property solution is carried out using fit correlation formula；In the region far from critical point, solved using actual gas state equation；?
Borderline region carries out interpolation solution using fit correlation formula and actual gas state equation；Solution saturated solution enthalpy can be solved, satisfied
Gentle enthalpy, temperature, specific volume, thermal conductivity and dynamic viscosity physical parameter.
The pressure and enthalpy provided according to user, substitutes into the saturated solution and saturated air enthalpy calculation formula of carbon dioxide, can
Saturated solution and saturated air enthalpy is calculated.Pressure and enthalpy substitute into the temperature computation formula of carbon dioxide, available two
The temperature of carbonoxide；Specific volume calculation formula is substituted into, the specific volume (or density) of carbon dioxide can be obtained.By the temperature of carbon dioxide and
Density substitutes into the calculation formula of thermal conductivity and dynamic viscosity, then the thermal conductivity and dynamic viscosity of carbon dioxide can be calculated.
Pass through physical property method for solving provided herein, in the case where known pressure and enthalpy, the saturated solution of available carbon dioxide
The parameters such as enthalpy, saturated air enthalpy, temperature, specific volume, thermal conductivity and dynamic viscosity, these parameters can be used for the calculating of reactor numerical value
Program.
A kind of wide parameter area carbon dioxide physical property method for solving of the present invention, comprising the following steps:
S1, saturation liquid phase enthalpy and saturation gas phase enthalpy are solved using fitting of a polynomial relational expression relevant to pressure；
The input parameter of saturated solution enthalpy is pressure
Wherein, hl1It (p) is the saturation liquid phase enthalpy of low-pressure area, kJ/kg, hl2It (p) is the saturation liquid phase enthalpy of higher-pressure region,
KJ/kg, hl3It (p) is the saturation liquid phase enthalpy in interpolation area, kJ/kg, p are pressure, and MPa, a, b, i are constant coefficient, and occurrence is shown in
Table 1.
1 enthalpy of saturated liquid fit correlation formula coefficient of table
The input parameter of saturated air enthalpy is pressure
Wherein, hg1It (p) is the saturation gas phase enthalpy of low-pressure area, kJ/kg, hg2It (p) is the saturation gas phase enthalpy of higher-pressure region,
KJ/kg, hg3It (p) is the saturation gas phase enthalpy in interpolation area, kJ/kg, p are pressure, and MPa, a, b, i are constant coefficient, and occurrence is shown in
Table 2.
2 saturated air specific enthalpy fit correlation formula coefficient of table
S2, temperature computation region was divided into 2nd area of cold-zone, one area of overheat and overheat, pressure is 0.1~20MPa, specific enthalpy
For 200~1600kJkg-1, it is as shown in Figure 1 that actual temp calculates subregion.
The range for crossing cold-zone is enthalpy of saturated liquid line and critical specific enthalpy line lower boundary (h=hset10=312.25kJ/kg) it
Under, region calculation formula of the pressure from 0.1MPa to 20MPa is as follows:
Wherein, T1(p, h) was the temperature of cold-zone, DEG C, p is pressure, and MPa, h are enthalpy, kJ/kg, aijIt is constant, specifically
Value is shown in Table 3；；
Table 3 crosses cold-zone temperature foh relational expression coefficient aij
The range for overheating an area is saturated air specific enthalpy line and critical specific enthalpy line coboundary (h=hset11=352.25kJ/kg)
Above, 600kJkg-1Specific enthalpy line is hereinafter, region of the pressure from 0.1MPa to 20MPa；
Wherein, T2(p, h) is the temperature for overheating an area, DEG C, p≤pcritical, hg(p) < h≤hset20Or p > pcritical,
hset11< h≤hset20, hset20For constant ,=580kJ/kg, bijFor constant coefficient, occurrence is shown in Table 4；
Table 4 overheats an area temperature foh relational expression coefficient bij
The range for overheating 2nd area is 620kJkg-1It is more than specific enthalpy line, 1600kJkg-1Specific enthalpy line hereinafter, pressure from
The region of 0.1MPa to 20MPa, calculation formula are as follows:
Wherein, T3(p, h) is the temperature for overheating 2nd area, DEG C, h > hset21, hset21For constant ,=620kJ/kg, cijIt is normal
Coefficient, occurrence are shown in Table 5；
Table 5 overheats two area temperature foh relational expression coefficient cij
The calculation formula in one area of interpolation is as follows:
Wherein, T4(p, h) is the temperature in one area of interpolation, DEG C, p > pcritical, hset10< h≤hset11, hset10For constant,
=312.25kJ/kg, hset11For constant ,=352.25kJ/kg；
The calculation formula in 2nd area of interpolation is as follows:
Wherein, T5(p, h) is the temperature in 2nd area of interpolation, DEG C, hset20< h≤hset21, hset20For constant ,=580kJ/kg,
hset21For constant ,=620kJ/kg.
S3, specific volume zoning was divided into cold-zone, state equation area, one area of overheat and 2nd area of overheat and was inserted accordingly
It is worth area；
Specific volume computer capacity is 0.1~20MPa of pressure, 200~1600kJkg of specific enthalpy-1, as shown in Figure 3.
The range for crossing cold-zone is enthalpy of saturated liquid line and 260kJkg-1Under specific enthalpy line, pressure is from 0.1MPa to 20MPa
Region；
State equation area is iteratively solved using actual gas state equation, and range includes two parts, and first part is
Below enthalpy of saturated liquid line and 296kJkg-1More than specific enthalpy line, more than saturated air line and 380kJkg-1Specific enthalpy line is hereinafter, pressure
Region of the power from 7MPa to 7.38MPa, the second part are 296kJkg-1More than specific enthalpy line and 380kJkg-1Specific enthalpy line with
Under, pressure is from 7.38MPa to 20MPa；
The range for overheating an area is saturated air line or more and 580kJkg-1Specific enthalpy line hereinafter, pressure from 0.1MPa to 7MPa
Region；Enthalpy is in 410kJkg-1More than specific enthalpy line and 580kJkg-1Specific enthalpy line hereinafter, pressure from 7MPa to 20MPa
Region.
The range for overheating 2nd area is 600kJkg-1It is more than specific enthalpy line, 1600kJkg-1Specific enthalpy line hereinafter, pressure from
The region of 0.1MPa to 20MPa.
The calculation formula for crossing cold-zone is as follows:
Wherein, v1(p, h) was the specific volume of cold-zone, m3/ kg, p≤7MPa, h < hl(p) or p > 7MPa, h≤hset30,
hset30For constant ,=260kJ/kg；aijSpecific value be shown in Table 6.
Table 6 crosses cold-zone specific volume fit correlation formula coefficient aij
Using Hall, hereby energy equation is not solved in state equation area:
φ (δ, τ)=φo(δ,τ)+φr(δ,τ)
Wherein, φ is the Hall not hereby energy of real gas, φoFor the Hall not hereby perfect gas part of energy, φrFor
The Hall not hereby remainder of energy, δ=ρ/ρc, τ=T/Tc, ρc=467.6kg/m3, Tc=31.1 DEG C, p >=pcritical,
hset00≤h≤hset01Or 7MPa≤p≤pcritical, hset00≤h≤hl(p) or 7MPa≤p≤pcritical, hg(p)≤h≤
hset01, hset00For constant ,=296kJ/kg, hset01For constant ,=380kJ/kg, pcriticalFor the critical pressure of carbon dioxide,
=7.31MPa.
An area is overheated to be calculated using following relational expression:
Wherein, v3(p, h) is the specific volume for overheating an area, p≤7MPa, hg(p) < h≤hset20Or p > 7MPa, hset31< h
≤hset20, hset31For constant ,=410kJ/kg；bijSpecific value be shown in Table 7.
Table 7 overheats an area specific volume fit correlation formula coefficient bij
The calculation relational expression for overheating 2nd area is as follows:
Wherein, v4(p, h) is the specific volume for overheating 2nd area, m3/ kg, h > hset21；cijSpecific value be shown in Table 8.
Table 8 overheats two area specific volume fit correlation formula coefficient aij
Interpolation area:
Wherein, v5(p, h) is the specific volume in one area of interpolation, m3/ kg, hset30=260kJ/kg, hset00=296kJ/kg, p >
7MPa, hset10< h≤hset00；
Wherein, v6(p, h) is the specific volume in 2nd area of interpolation, m3/ kg, hset31=410kJ/kg is hset01=380kJ/kg, p
> 7MPa, hset01< h≤hset31；
Wherein, v7(p, h) is the specific volume in 3rd area of interpolation, m3/ kg, hset20For constant ,=580kJ/kg, hset21For constant,
=620kJ/kg, hset20< h≤hset21。
S4, thermal conductivity is calculated, computer capacity is in 216~1100K of temperature, pressure to 200MPa；
λ (ρ, T)=λ0(T)+Δλ(ρ,T)+Δλc(ρ,T)
Wherein, λ0(T) thermal conductivity in zero density low density gas is indicated, Δ λ (ρ, T) represents thin gas at identical temperature
Volume density increases bring thermal conductivity and increases, Δ λc(ρ, T) represents the increase in critical point peak region thermal conductivity.
S5, dynamic viscosity is calculated.
Computer capacity is in 200~1500K of temperature, density to 1400kg/m3
η (ρ, T)=η0(T)+Δη(ρ,T)+Δηc(ρ,T)
Wherein, η0(T) viscosity in zero density low density gas is indicated, Δ η (ρ, T) represents low density gas at identical temperature
Density increases bring viscosity and increases, Δ ηc(ρ, T) represents the increase in critical point peak region viscosity.
In order to make the object, technical scheme and advantages of the embodiment of the invention clearer, below in conjunction with the embodiment of the present invention
In attached drawing, technical scheme in the embodiment of the invention is clearly and completely described, it is clear that described embodiment is
A part of the embodiment of the present invention, instead of all the embodiments.The present invention being described and shown in usually here in attached drawing is real
The component for applying example can be arranged and be designed by a variety of different configurations.Therefore, below to the present invention provided in the accompanying drawings
The detailed description of embodiment be not intended to limit the range of claimed invention, but be merely representative of of the invention selected
Embodiment.Based on the embodiments of the present invention, those of ordinary skill in the art are obtained without creative efforts
The every other embodiment obtained, shall fall within the protection scope of the present invention.
Enthalpy of saturated liquid is calculated, in 0.6~7.377MPa, takes 6777 pressure spots to be saturated by step-length of 1kPa
Liquor ratio enthalpy calculates, and it is 0.08% that method of the invention, which calculates maximum relative error, as shown in Figure 3；Calculating speed is business software
118 times of NIST PREPROP, as shown in Figure 4.
Saturated air specific enthalpy is calculated, in 0.6~7.377MPa, takes 6777 pressure spots to be saturated by step-length of 1kPa
Liquor ratio enthalpy calculates, and it is 0.009% that method of the invention, which calculates maximum relative error, as shown in Figure 5；Calculating speed is business software
230 times of NIST PREPROP, as shown in Figure 6.
For temperature computation, in 0.1~20MPa of pressure limit, step-length 0.1MPa, than enthalpy spectrum 200-1600kJkg-1、
Step-length 2kJkg-1, take 131914 state points to be calculated, it is 0.25% that method of the invention, which calculates maximum relative error, such as
Shown in Fig. 7；Calculating speed is 15 times of business software NIST PREPROP, as shown in Figure 8.
Specific volume is calculated, in 0.1~20MPa of pressure limit, step-length 0.1MPa, than enthalpy spectrum 200-1600kJkg-1、
Step-length 2kJkg-1, take 131914 state points to be calculated, method of the invention calculates maximum relative error within 1%,
As shown in figure 9, calculating speed is 28 times of business software NIST PREPROP, as shown in Figure 10.
Calculating for thermal conductivity, for maximum relative error within 5%, calculating speed is business software NIST PREPROP
Quite.
Calculating for dynamic viscosity, for maximum relative error within 5%, calculating speed is business software NIST
PREPROP is suitable.
The above content is merely illustrative of the invention's technical idea, and this does not limit the scope of protection of the present invention, all to press
According to technical idea proposed by the present invention, any changes made on the basis of the technical scheme each falls within claims of the present invention
Protection scope within.
Claims (10)
1. a kind of width parameter area carbon dioxide physical property method for solving, which is characterized in that adopted in the region of carbon dioxide critical point
Solution is iterated with actual gas state equation；Object is carried out using fit correlation formula in the region far from carbon dioxide critical point
Property solve；Interpolation solution is carried out using fit correlation formula and actual gas state equation in borderline region, saturated solution is calculated
Enthalpy, saturated air enthalpy, temperature, specific volume, thermal conductivity and dynamic viscosity physical parameter.
2. width parameter area carbon dioxide physical property method for solving according to claim 1, which is characterized in that use and pressure
It is as follows that relevant fitting of a polynomial relational expression solves saturation liquid phase enthalpy:
Wherein, hl1It (p) is the saturation liquid phase enthalpy of low-pressure area, hl2It (p) is the saturation liquid phase enthalpy of higher-pressure region, hl3It (p) is interpolation
The saturation liquid phase enthalpy in area, p are pressure, and MPa, a, b, i are constant coefficient.
3. width parameter area carbon dioxide physical property method for solving according to claim 1, which is characterized in that use and pressure
It is as follows that relevant fitting of a polynomial relational expression solves saturation gas phase enthalpy:
Wherein, hg1It (p) is saturation gas phase the enthalpy g, h of low-pressure areag2It (p) is the saturation gas phase enthalpy of higher-pressure region, hg3It (p) is slotting
It is worth the saturation gas phase enthalpy in area, p is pressure, and MPa, a, b, i are constant coefficient.
4. width parameter area carbon dioxide physical property method for solving according to claim 1, which is characterized in that by the meter of temperature
It calculates region and was divided into 2nd area of cold-zone, one area of overheat and overheat, pressure is 0.1~20MPa, specific enthalpy is 200~1600kJkg-1。
5. width parameter area carbon dioxide physical property method for solving according to claim 4, which is characterized in that cross the model of cold-zone
Enclose is enthalpy of saturated liquid line and critical specific enthalpy line lower boundary (312.25kJkg-1) under, area of the pressure from 0.1MPa to 20MPa
Domain；
The range for overheating an area is saturated air specific enthalpy line and critical specific enthalpy line coboundary (352.25kJkg-1) more than, 600kJ
kg-1Specific enthalpy line is hereinafter, region of the pressure from 0.1MPa to 20MPa；
The range for overheating 2nd area is 600kJkg-1It is more than specific enthalpy line, 1600kJkg-1Specific enthalpy line hereinafter, pressure from 0.1MPa to
The region of 20MPa.
6. width parameter area carbon dioxide physical property method for solving according to claim 5, which is characterized in that cross the temperature of cold-zone
Spend T1(p, h) calculates as follows:
Wherein, aijFor constant coefficient, p is pressure, and MPa, h are enthalpy, kJ/kg, p≤pcritical, h < hl(p) or p > pcritical, h
≤hset10, hset10=hcritical- 20.0=312.25kJ/kg；
Overheat the temperature T in an area2(p, h) calculates as follows:
Wherein, bijFor constant coefficient, p≤pcritical, hg(p) < h≤hset20Or p > pcritical, hset11< h≤hset20, hset11
=hcritical+ 20.0=352.25kJ/kg, hset20For constant, hset20=580kJ/kg；
Overheat the temperature T in 2nd area3(p, h) calculates as follows:
Wherein, cijFor constant coefficient, h > hset21, hset21For constant, hset21=620kJ/kg；
One area of interpolation was the transitional region in one area of cold-zone and overheat, temperature T4(p, h) calculates as follows:
Wherein, h is enthalpy, p > pcritical, hset10< h≤hset11；
2nd area of interpolation is the transitional region for overheating an area and overheating 2nd area, temperature T5(p, h) calculates as follows:
7. width parameter area carbon dioxide physical property method for solving according to claim 1, which is characterized in that by the meter of specific volume
Calculate region and be divided into cold-zone, 2nd area of state equation area, one area of overheat and overheat, pressure is 0.1~20MPa, specific enthalpy be 200~
1600kJ·kg-1。
8. width parameter area carbon dioxide physical property method for solving according to claim 7, which is characterized in that cross the model of cold-zone
Enclose is enthalpy of saturated liquid line and 260kJkg-1Under specific enthalpy line, region of the pressure from 0.1MPa to 20MPa；
State equation area uses actual gas Solving Equation of State, and range includes two parts, and first part is enthalpy of saturated liquid line
Below and 296kJkg-1More than specific enthalpy line, more than saturated air line and 380kJkg-1Specific enthalpy line hereinafter, pressure from 7MPa to
The region of 7.38MPa, second part are 296kJkg-1More than specific enthalpy line and 380kJkg-1Specific enthalpy line hereinafter, pressure from
7.38MPa arriving 20MPa；
The range for overheating an area is saturated air line or more and 580kJkg-1Specific enthalpy line is hereinafter, area of the pressure from 0.1MPa to 7MPa
Domain；Enthalpy is in 410kJkg-1More than specific enthalpy line and 580kJkg-1Specific enthalpy line is hereinafter, region of the pressure from 7MPa to 20MPa；
The range for overheating 2nd area is 600kJkg-1It is more than specific enthalpy line, 1600kJkg-1Specific enthalpy line hereinafter, pressure from 0.1MPa to
The region of 20MPa.
9. width parameter area carbon dioxide physical property method for solving according to claim 8, which is characterized in that cross the ratio of cold-zone
Hold v1(p, h) calculates as follows:
Wherein, aijFor constant coefficient, p≤7MPa, h < hl(p) or p > 7MPa, h≤hset30, hset30For constant ,=260kJ/kg；
Hereby function does not calculate using Hall for state equation area:
φ (δ, τ)=φo(δ,τ)+φr(δ,τ)
Wherein, φ is the Hall not hereby energy of real gas, φoFor the Hall not hereby perfect gas part of energy, φrFor Hall
The not hereby remainder of energy, δ=ρ/ρc, τ=T/Tc, ρc=467.6kg/m3, Tc=31.1 DEG C, p >=pcritical, hset00≤h
≤hset01Or 7MPa≤p≤pcritical, hset00≤h≤hl(p) or 7MPa≤p≤pcritical, hg(p)≤h≤hset01, hset00
For constant ,=296kJ/kg, hset01For constant ,=380kJ/kg, pcriticalFor the critical pressure of carbon dioxide ,=
7.31MPa；
Overheat the specific volume v in an area3(p, h) calculates as follows:
Wherein, bijFor constant coefficient, p≤7MPa, hg(p) < h≤hset20Or p > 7MPa, hset31< h≤hset20, hset31It is normal
Number ,=410kJ/kg；
Overheat the specific volume v in 2nd area4(p, h) calculates as follows:
Wherein, cijFor constant coefficient, h > hset21；
One area of interpolation was the transition region of cold-zone and state equation area, specific volume v5(p, h) calculates as follows:
Wherein, hset30=260kJ/kg, hset00=296kJ/kg, p > 7MPa, hset30< h≤hset00；
2nd area of interpolation is the transition region in one area of state equation area and overheat, specific volume v6(p, h) calculates as follows:
Wherein, hset31=410kJ/kg, hset01=380kJ/kg, p > 7MPa, hset01< h≤hset31；
3rd area of interpolation is the transition region for overheating an area and overheating 2nd area, specific volume v7(p, h) calculates as follows:
Wherein, hset20=580kJ/kg, hset21=620kJ/kg, hset20< h≤hset21。
10. width parameter area carbon dioxide physical property method for solving according to claim 1, which is characterized in that thermal conductivity
Computer capacity calculates as follows in 216~1100K of temperature, pressure to 200MPa:
λ (ρ, T)=λ0(T)+Δλ(ρ,T)+Δλc(ρ,T)
Wherein, λ0(T) thermal conductivity in zero density low density gas is indicated, it is close that Δ λ (ρ, T) represents low density gas at identical temperature
Degree increases bring thermal conductivity and increases, Δ λc(ρ, T) represents the increase in critical point peak region thermal conductivity；
The computer capacity of dynamic viscosity is in 200~1500K of temperature, density to 1400kg/m3, it is specific as follows:
η (ρ, T)=η0(T)+Δη(ρ,T)+Δηc(ρ,T)
Wherein, η0(T) viscosity in zero density low density gas is indicated, Δ η (ρ, T) represents low density gas density at identical temperature
It increases bring viscosity to increase, Δ ηc(ρ, T) represents the increase in critical point peak region viscosity.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811626329.3A CN109783904B (en) | 2018-12-28 | 2018-12-28 | Wide parameter range carbon dioxide physical property solving method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811626329.3A CN109783904B (en) | 2018-12-28 | 2018-12-28 | Wide parameter range carbon dioxide physical property solving method |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109783904A true CN109783904A (en) | 2019-05-21 |
| CN109783904B CN109783904B (en) | 2020-10-27 |
Family
ID=66498719
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201811626329.3A Active CN109783904B (en) | 2018-12-28 | 2018-12-28 | Wide parameter range carbon dioxide physical property solving method |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109783904B (en) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116030906A (en) * | 2023-03-28 | 2023-04-28 | 西安交通大学 | Physical property parameter calculation method for lead bismuth fast reactor multicomponent fluid |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP4980198B2 (en) * | 2007-10-26 | 2012-07-18 | 日東電工株式会社 | Foam material with slit for mounting step structure |
| CN104777183A (en) * | 2014-10-31 | 2015-07-15 | 北京卫星环境工程研究所 | Satellite electric propulsion system xenon filling thermodynamic characteristic numerical simulation method |
| CN107506534A (en) * | 2017-08-04 | 2017-12-22 | 陕西延长石油（集团）有限责任公司 | A kind of carbon dioxide drive seals middle cap rock sealed harmonic drive method up for safekeeping |
| CN108134156A (en) * | 2017-11-28 | 2018-06-08 | 江西爱驰亿维实业有限公司 | The calculation method of parameters of tube refrigerant fluid interchange, system, medium, terminal, battery pack |
| CN108763670A (en) * | 2018-05-15 | 2018-11-06 | 西安交通大学 | A kind of solution supercritical carbon dioxide reactor Brayton cycle transient process method |
-
2018
- 2018-12-28 CN CN201811626329.3A patent/CN109783904B/en active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP4980198B2 (en) * | 2007-10-26 | 2012-07-18 | 日東電工株式会社 | Foam material with slit for mounting step structure |
| CN104777183A (en) * | 2014-10-31 | 2015-07-15 | 北京卫星环境工程研究所 | Satellite electric propulsion system xenon filling thermodynamic characteristic numerical simulation method |
| CN107506534A (en) * | 2017-08-04 | 2017-12-22 | 陕西延长石油（集团）有限责任公司 | A kind of carbon dioxide drive seals middle cap rock sealed harmonic drive method up for safekeeping |
| CN108134156A (en) * | 2017-11-28 | 2018-06-08 | 江西爱驰亿维实业有限公司 | The calculation method of parameters of tube refrigerant fluid interchange, system, medium, terminal, battery pack |
| CN108763670A (en) * | 2018-05-15 | 2018-11-06 | 西安交通大学 | A kind of solution supercritical carbon dioxide reactor Brayton cycle transient process method |
Non-Patent Citations (2)
| Title |
|---|
| PAN WU等: "Development and Verification of a Transient Analysis Tool for Reactor System Using Supercritical CO2 Brayton Cycle as Power Conversion System", 《SCIENCE AND TECHNOLOGY OF NUCLEAR INSTALLATIONS》 * |
| 蒋雪峰: "超临界CO_2压缩机数值模拟与设计研究", 《中国优秀硕士学位论文全文数据库 工程科技Ⅱ辑》 * |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116030906A (en) * | 2023-03-28 | 2023-04-28 | 西安交通大学 | Physical property parameter calculation method for lead bismuth fast reactor multicomponent fluid |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109783904B (en) | 2020-10-27 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Wagner et al. | IAPWS industrial formulation 1997 for the thermodynamic properties of water and steam | |
| Woudstra et al. | Thermodynamic evaluation of combined cycle plants | |
| Wang et al. | A new complex variable meshless method for transient heat conduction problems | |
| Županović et al. | Kirchhoff’s loop law and the maximum entropy production principle | |
| CN109325283A (en) | A kind of cogeneration units amount of energy saving calculation method and device | |
| CN109783904A (en) | A kind of width parameter area carbon dioxide physical property method for solving | |
| CN103335538B (en) | The computational methods of a kind of condenser of power station pressure and terminal temperature difference | |
| Kunick et al. | CFD Analysis of steam turbines with the IAPWS standard on the Spline-Based Table Look-Up Method (SBTL) for the fast calculation of real fluid properties | |
| Mezei | Direct calculation of the excess free energy of the dense Lennard-Jones fluid | |
| CN111241728A (en) | Intermittent Galerkin finite element numerical solution method of Euler equation | |
| CN109344423A (en) | A kind of calculation method for closing the practical IP efficiency of cylinder steam turbine | |
| Wallace | Statistical entropy and a qualitative gas-liquid phase diagram | |
| Aktershev et al. | Numerical modeling of the vapor bubble growth in a homogenously superheated liquid (thermal energy scheme) | |
| Woods | On the role of the harmonic mean isentropic exponent in the analysis of the closed-cycle gas turbine | |
| Zou et al. | Numerical Uncertainties vs. Model Uncertainties in Two-Phase Flow Simulations | |
| CN113821998B (en) | Method for solving shell side pressure of condenser real-time dynamic simulation model by Newton iteration method | |
| CN105509522A (en) | Manufacturing method of sintered copper powder and high-porosity copper foam composited heat pipe | |
| CN108038224A (en) | A kind of new water physical property querying method of computer based | |
| Chen et al. | A Cartesian grid-based unified gas kinetic scheme | |
| CN104318092B (en) | SCV (submerged combustion vaporizer) design method | |
| Li et al. | Numerical investigation of equilibrium wet steam flow property based on S2 calculation code | |
| Awad | Comments on “Analysis of heat transfer and irreversibility of ORC evaporator for selecting working fluid and operating conditions” | |
| CN114415748B (en) | Sleeve regulating valve resistance coefficient estimation method and system considering cavitation influence | |
| Cai et al. | Application of EXCEL Univariate Solution Function in Thermal Calculation | |
| Chen | Short-time critical behavior of anisotropic cubic systems |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |