CN109783904B - Wide parameter range carbon dioxide physical property solving method - Google Patents

Wide parameter range carbon dioxide physical property solving method Download PDF

Info

Publication number
CN109783904B
CN109783904B CN201811626329.3A CN201811626329A CN109783904B CN 109783904 B CN109783904 B CN 109783904B CN 201811626329 A CN201811626329 A CN 201811626329A CN 109783904 B CN109783904 B CN 109783904B
Authority
CN
China
Prior art keywords
zone
pressure
region
carbon dioxide
critical
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201811626329.3A
Other languages
Chinese (zh)
Other versions
CN109783904A (en
Inventor
吴攀
单建强
高春天
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xian Jiaotong University
Original Assignee
Xian Jiaotong University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Xian Jiaotong University filed Critical Xian Jiaotong University
Priority to CN201811626329.3A priority Critical patent/CN109783904B/en
Publication of CN109783904A publication Critical patent/CN109783904A/en
Application granted granted Critical
Publication of CN109783904B publication Critical patent/CN109783904B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Landscapes

  • Organic Low-Molecular-Weight Compounds And Preparation Thereof (AREA)

Abstract

The invention discloses a wide parameter range carbon dioxide physical property solving method, which adopts a real gas state equation to carry out iterative solution in a region of a carbon dioxide critical point; directly solving the physical properties in a region far away from the critical point of the carbon dioxide by adopting a fitting relation; and (4) performing interpolation solution in the boundary area by adopting a fitting relation and a real gas state equation, and calculating to obtain saturated liquid enthalpy value, saturated gas enthalpy value, temperature, specific volume, thermal conductivity and dynamic viscosity physical property parameters. The method not only ensures the calculation precision of the carbon dioxide, but also improves the solving speed of the reactor numerical calculation program.

Description

Wide parameter range carbon dioxide physical property solving method
Technical Field
The invention belongs to the technical field of nuclear reactor safety analysis and calculation, and particularly relates to a method for solving physical properties of carbon dioxide in a wide parameter range.
Background
The supercritical carbon dioxide Brayton cycle has the advantages of high thermal efficiency, simplified and compact system and the like, and is a promising fourth-generation nuclear reactor energy conversion system. When the supercritical carbon dioxide is in Brayton cycle steady-state operation, the supercritical carbon dioxide at the inlet of the compressor is near a critical point, the carbon dioxide undergoes violent and rapid physical property changes, and key physical property parameters such as specific volume, specific heat, thermal conductivity and the like of the carbon dioxide are all distorted. In addition, in a transient state or an accident, the operating state point of the supercritical carbon dioxide Brayton cycle will experience a large range of physical property changes, so the accurate carbon dioxide physical property calculation method has very important use value for realizing the steady-state design, the transient control strategy and the accident analysis of the supercritical carbon dioxide Brayton cycle reactor.
At present, the thermodynamic physical property of fluid is calculated mainly by three methods, namely a graph method, a state equation method and a fitting relational expression method. The graphical method is to calculate the thermodynamic properties of a fluid by means of a prepared graph and table, and has the advantages of simplicity and easy understanding, but has the disadvantages of low calculation efficiency and is not suitable for large-scale thermal analysis programs. The equation of state method is a fundamental method for fluid thermodynamic property calculation, and the method is established based on strict theoretical and experimental researches, so that the method has the advantage of high precision, but has the defects of complex model, long calculation time and complex programming. The fitting relation method is a simplified calculation method of fluid thermodynamic physical properties, and obtains a mathematical relation meeting certain precision by fitting the existing thermodynamic physical property data. Because the adopted relational expressions are simple mathematical relational expressions such as polynomial expressions, fraction expressions and the like, the method has the advantages of simple calculation, small calculated amount, higher calculation precision and suitability for computer programming solution, and is widely adopted in thermodynamic calculation programs.
The physical properties of the carbon dioxide with higher precision can be calculated by the method. However, the equation of state method requires a large number of iterations to obtain the physical property result, and the calculation time is long. The physical property calculation is very important in the thermal calculation program, and particularly in the transient thermal calculation program, the physical property calculation occupies about 20% of the program calculation time. If the physical property calculation speed is slow, the calculation speed of the transient thermodynamic calculation program is reduced by times. The physical property is calculated by using the fitting relation, so that the calculation efficiency can be effectively improved. However, the physical property distortion of carbon dioxide in the vicinity of the critical point results in insufficient calculation accuracy of the fitting relational expression.
The physical property characteristics of the carbon dioxide cause that the existing transient and accident numerical analysis program which can be used for the Brayton cycle of the supercritical carbon dioxide reactor does not meet the requirements on the calculation precision or the calculation speed. The method for rapidly and accurately solving the physical property of the carbon dioxide in the wide parameter range is urgently needed to be developed, and the calculation precision and the calculation speed of the transient thermodynamic calculation program are improved.
Disclosure of Invention
The technical problem to be solved by the present invention is to provide a method for solving the physical properties of carbon dioxide in a wide parameter range, which can ensure the calculation accuracy of carbon dioxide and improve the solving speed of a reactor numerical calculation program, aiming at the defects in the prior art.
The invention adopts the following technical scheme:
a wide parameter range carbon dioxide physical property solving method adopts a real gas state equation to carry out iterative solution in a region of a carbon dioxide critical point; solving the physical property in the region far away from the critical point of the carbon dioxide by adopting a fitting relation; and (4) performing interpolation solution in the boundary area by adopting a fitting relation and a real gas state equation, and calculating to obtain saturated liquid enthalpy value, saturated gas enthalpy value, temperature, specific volume, thermal conductivity and dynamic viscosity physical property parameters.
Specifically, solving the saturated liquid phase enthalpy value by adopting a polynomial fitting relation related to pressure as follows:
Figure BDA0001928065530000031
wherein h isl1(p) saturated liquid phase enthalpy of the low pressure zone, hl2(p) is the saturated liquid phase enthalpy of the high pressure zone, hl3And (p) is the saturated liquid phase enthalpy value of the interpolation zone, p is the pressure, MPa, a, b and i are all constant coefficients.
Specifically, solving the saturated vapor phase enthalpy value by adopting a polynomial fitting relation related to pressure as follows:
Figure BDA0001928065530000032
wherein h isg1(p) the saturated vapor phase enthalpy value g, h of the low-pressure zoneg2(p) is the saturated vapor phase enthalpy of the high pressure zone, hg3And (p) is the saturated vapor phase enthalpy value of the interpolation zone, p is the pressure, MPa, a, b and i are all constant coefficients.
Specifically, the temperature calculation area is divided into an overcooling area, an overheating area I and an overheating area II, the pressure is 0.1-20 MPa, and the specific enthalpy is 200-1600 kJ-kg-1。
Further, the range of the supercooling region is the lower boundary (312.25 kJ. kg) of the saturated liquid specific enthalpy line and the critical specific enthalpy line-1) A region under a pressure of from 0.1MPa to 20 MPa; the range of the first superheat zone is the saturation gas specific enthalpy line and the upper boundary of the critical specific enthalpy line (352.25 kJ-kg)-1) Above, 600 kJ.kg-1A region below the specific enthalpy line and at a pressure of from 0.1MPa to 20 MPa; the range of the second superheat zone is 600 kJ-kg -11600 kJ.kg above the specific enthalpy line-1A region where the pressure is from 0.1MPa to 20MPa below the specific enthalpy line.
Further, the temperature T of the supercooling region1(p, h) is calculated as follows:
Figure BDA0001928065530000033
wherein, aijIs constant coefficient, p is pressure, MPa, h is enthalpy, kJ/kg, and p is not more than pcritical，h＜hl(p) or p > pcritical，h≤hset10，hset10＝hcritical-20.0＝312.25kJ/kg；
Temperature T of superheat zone2(p, h) is calculated as follows:
Figure BDA0001928065530000041
wherein, bijIs a constant coefficient, p is less than or equal to pcritical，hg(p)＜h≤hset20Or p > pcritical，hset11＜h≤hset20，hset11＝hcritical+20.0＝352.25kJ/kg，hset20Is constant, 580 kJ/kg;
temperature T of the superheat zone two3(p, h) is calculated as follows:
Figure BDA0001928065530000042
wherein, cijIs a constant coefficient, h > hset21，hset21Is constant, 620 kJ/kg;
the interpolation zone is the transition zone between the supercooling zone and the superheating zone, and has a temperature T4(p, h) is calculated as follows:
Figure BDA0001928065530000043
wherein h is the enthalpy value, p > pcritical，hset10＜h≤hset11；
The second interpolation zone is a transition zone between the first superheating zone and the second superheating zone, and the temperature T of the second interpolation zone is5(p, h) is calculated as follows:
Figure BDA0001928065530000044
specifically, the calculation of specific volume is divided intoThe area is divided into an overcooling area, an equation of state area, an overheating area I and an overheating area II, the pressure is 0.1-20 MPa, and the specific enthalpy is 200-1600 kJ.kg-1。
Further, the range of the supercooling zone is a saturated liquid specific enthalpy line and 260 kJ.kg-1A region below the specific enthalpy line at a pressure of from 0.1MPa to 20 MPa; the state equation region is solved by using a real gas state equation, and the range comprises two parts, wherein the first part is below a saturated liquid specific enthalpy line and 296 kJ.kg-1Above specific enthalpy line, above saturated gas line and 380 kJ.kg-1A pressure of 7MPa to 7.38MPa below the specific enthalpy line, and a second fraction of 296 kJ-kg-1Above the specific enthalpy line and 380kJ/kg-1Below the specific enthalpy line, the pressure is from 7.38MPa to 20 MPa; the range of the first superheating area is more than the saturated gas line and 580 kJ.kg-1A region below the specific enthalpy line and at a pressure of from 0.1MPa to 7 MPa; the enthalpy value is 410kJ kg-1Above the specific enthalpy line and 580kJ/kg-1A region having a pressure of 7MPa to 20MPa below a specific enthalpy line; the range of the second superheat zone is 600 kJ-kg -11600 kJ.kg above the specific enthalpy line-1A region where the pressure is from 0.1MPa to 20MPa below the specific enthalpy line.
Further, the specific volume v of the supercooling region1(p, h) is calculated as follows:
Figure BDA0001928065530000051
wherein, aijA constant coefficient, p is less than or equal to 7MPa, h is less than hl(p) or p is more than 7MPa, h is less than or equal to hset30，hset30Is constant, 260 kJ/kg;
the equation of state region is calculated using the Hall Motz function:
φ(,τ)＝φo(,τ)+φr(,τ)
where φ is the Hall Moz energy of the actual gas, φoIs the ideal gas portion of the Hall Motz energy, phirFor the remainder of the Hall Motz energy, ═ ρ/ρc，τ＝T/Tc，ρc＝467.6kg/m3，Tc＝31.1℃，p≥pcritical，hset00≤h≤hset01Or p is more than or equal to 7MPa and less than or equal to pcritical，hset00≤h≤hl(p) or p is not less than 7MPa and not more than pcritical，hg(p)≤h≤hset01，hset00Is constant, 296kJ/kg, hset01Is constant, 380kJ/kg, pcriticalCritical pressure of carbon dioxide, 7.31 MPa;
specific volume v of superheat region3(p, h) is calculated as follows:
Figure BDA0001928065530000052
wherein, bijIs a constant coefficient, p is less than or equal to 7MPa, hg(p)＜h≤hset20Or p > 7MPa, hset31＜h≤hset20，hset31Is constant, 410 kJ/kg;
specific volume v of the overheated second zone4(p, h) is calculated as follows:
Figure BDA0001928065530000053
wherein, cijIs a constant coefficient, h > hset21；
The interpolation zone is a transition zone of the supercooling zone and the equation of state zone and has a specific volume v5(p, h) is calculated as follows:
Figure BDA0001928065530000061
wherein h isset30＝260kJ/kg，hset00＝296kJ/kg，p＞7MPa，hset30＜h≤hset00；
The interpolation two zone is a transition zone of the equation of state zone and the overheating one zone and has a specific volume v6(p, h) is calculated as follows:
Figure BDA0001928065530000062
wherein h isset31＝410kJ/kg，hset01＝380kJ/kg，p＞7MPa，hset01＜h≤hset31；
The three interpolation zones are transition zones of the first overheating zone and the second overheating zone and have specific volume v7(p, h) is calculated as follows:
Figure BDA0001928065530000063
wherein h isset20＝580kJ/kg，hset21＝620kJ/kg，hset20＜h≤hset21。
Specifically, the calculation range of the thermal conductivity is 216-1100K, and the calculation range of the thermal conductivity is as follows:
λ(ρ,T)＝λ0(T)+Δλ(ρ,T)+Δλc(ρ,T)
wherein λ is0(T) represents the thermal conductivity at zero density of the lean gas, Δ λ (ρ, T) represents the increase in thermal conductivity due to the increase in density of the lean gas at the same temperature, Δ λc(ρ, T) represents the increase in thermal conductivity in the critical point peak region;
the dynamic viscosity is calculated in the temperature range of 200-1500K and the density range of 1400kg/m3The method comprises the following steps:
η(ρ,T)＝η0(T)+Δη(ρ,T)+Δηc(ρ,T)
wherein eta is0(T) represents the viscosity at zero density of the lean gas, Δ η (ρ, T) represents the increase in viscosity due to the increase in density of the lean gas at the same temperature, Δ ηc(ρ, T) represents the increase in viscosity in the peak region at the critical point.
Compared with the prior art, the invention has at least the following beneficial effects:
the invention relates to a wide parameter range carbon dioxide physical property solving method, which combines a fitting polynomial and a state equation and can quickly and accurately calculate physical property parameters such as saturated liquid enthalpy value, saturated gas enthalpy value, temperature, specific volume, thermal conductivity, dynamic viscosity and the like.
Further, for saturated liquid specific enthalpy calculation, 6777 pressure points are taken by taking 1kPa as a step length to perform saturated liquid specific enthalpy calculation under 0.6-7.377 MPa, the maximum relative error of the calculation of the method is 0.08%, and the calculation speed is 118 times of that of commercial software NIST PREPROP.
Further, for saturated gas specific enthalpy calculation, 6777 pressure points are taken by taking 1kPa as a step length to perform saturated liquid specific enthalpy calculation under 0.6-7.377 MPa, the maximum relative error of the calculation of the method is 0.009%, and the calculation speed is 230 times of that of commercial software NIST PREPROP.
Further, for the temperature calculation, the pressure range is 0.1-20 MPa, the step length is 0.1MPa, and the specific enthalpy range is 200--1Step length of 2 kJ/kg-1And 131914 state points are taken for calculation, the maximum relative error calculated by the method is 0.1%, and the calculation speed is 15 times of that of commercial software NIST PREPROP.
Further, for specific volume calculation, the pressure range is 0.1-20 MPa, the step length is 0.1MPa, and the specific enthalpy range is 200--1Step length of 2 kJ/kg-1131914 state points are taken for calculation, the maximum relative error of the calculation of the method is within 1 percent, and the calculation speed is 28 times of that of commercial software NIST PREPROP.
Further, for the calculation of the thermal conductivity, the maximum relative error is within 5%, and the calculation speed is equivalent to that of the commercial software NIST PREPROP.
Further, for the calculation of dynamic viscosity, the maximum relative error is within 5%, and the calculation speed is equivalent to commercial software NIST PREPROP.
In summary, the invention develops a carbon dioxide physical property solving method based on an autonomously developed fitting polynomial and a carbon dioxide physical property calculation relational expression and a real gas state equation reported in a published literature, and the carbon dioxide physical property solving method is used for solving the physical property of a reactor numerical analysis program related to carbon dioxide. The physical property solving method combining the fitting polynomial and the state equation is provided, the fitting relational expression which is independently researched and developed is adopted to directly solve the parameters such as saturated liquid specific enthalpy, saturated gas specific enthalpy and temperature, the calculation formula published by the open literature is adopted to directly solve the parameters such as heat conductivity and dynamic viscosity, and the calculation speed of the physical property bag is improved. In the specific volume parameter, near the critical point, the carbon dioxide state equation is mainly adopted for iterative solution, so that the specific volume calculation precision is ensured, and when the specific volume parameter is far away from the critical point, the self-developed fitting relation is adopted for direct solution.
The technical solution of the present invention is further described in detail by the accompanying drawings and embodiments.
Drawings
FIG. 1 is a schematic diagram of temperature calculation;
FIG. 2 is a schematic diagram of specific volume calculation;
FIG. 3 is a diagram showing the calculation results of saturated liquid specific enthalpy, wherein (a) is the relative error of the saturated liquid specific enthalpy low-pressure region relation, and (b) is the relative error of the saturated liquid specific enthalpy high-pressure region relation;
FIG. 4 is a comparison of saturated liquid specific enthalpy fit relation to commercial software NIST REFPROP calculated velocity;
FIG. 5 is a schematic diagram of saturated specific enthalpy, wherein (a) is the relative error of the relationship between the low pressure regions of saturated specific enthalpy and (b) is the relative error of the relationship between the high pressure regions of saturated specific enthalpy;
FIG. 6 is a comparison of a saturated air specific enthalpy fit relation to a commercial software NIST REFPROP calculated velocity;
FIG. 7 is a graph showing the temperature calculation results, wherein (a) is the relative error of the temperature relationship of the supercooling region, (b) is the relative error of the temperature relationship of the first superheating region, and (c) is the relative error of the temperature relationship of the second superheating region;
FIG. 8 is a comparison of the temperature fit relationship with the commercial software NIST REFPROP calculated speed;
FIG. 9 is a diagram showing the calculated specific volume, wherein (a) is the relative error of the specific volume relation of the over-cooling area, and (b) is the relative error of the temperature relation of the over-heating area;
FIG. 10 is a comparison of the specific volume fit relationship with the commercial software NIST REFPROP calculation speed.
Detailed Description
The invention provides a wide parameter range carbon dioxide physical property solving method, which adopts a fitting relation formula to solve the physical property in the region of a carbon dioxide critical point; solving by adopting a real gas state equation in a region far away from the critical point; in the boundary area, interpolation solution is carried out by adopting a fitting relation and a real gas state equation; can solve the physical parameters of saturated liquid enthalpy value, saturated gas enthalpy value, temperature, specific volume, thermal conductivity and dynamic viscosity.
And substituting the pressure and enthalpy values provided by the user into a saturated liquid and saturated gas enthalpy value calculation formula of the carbon dioxide, and calculating to obtain the saturated liquid and saturated gas enthalpy values. Substituting the pressure and the enthalpy value into a carbon dioxide temperature calculation formula to obtain the temperature of the carbon dioxide; substituting into the specific volume calculation formula to obtain the specific volume (or density) of the carbon dioxide. And substituting the temperature and the density of the carbon dioxide into a calculation formula of the thermal conductivity and the dynamic viscosity, so that the thermal conductivity and the dynamic viscosity of the carbon dioxide can be calculated. Through the physical property solving method provided by the invention, under the condition of known pressure and enthalpy, parameters such as saturated liquid enthalpy, saturated gas enthalpy, temperature, specific volume, thermal conductivity and dynamic viscosity of carbon dioxide can be obtained, and the parameters can be used for a reactor numerical calculation program.
The invention discloses a wide parameter range carbon dioxide physical property solving method, which comprises the following steps:
s1, solving a saturated liquid phase enthalpy value and a saturated gas phase enthalpy value by adopting a pressure-related polynomial fitting relation;
the input parameter of the enthalpy value of the saturated liquid is pressure
Figure BDA0001928065530000091
Wherein h isl1(p) is the saturated liquid phase enthalpy of the low pressure zone, kJ/kg, hl2(p) is the saturated liquid phase enthalpy value of the high pressure zone, kJ/kg, hl3(p) is the saturated liquid phase enthalpy value of the interpolation zone, kJ/kg, p is pressure, MPa, a, b and i are all constant coefficients, and the specific values are shown in Table 1.
TABLE 1 coefficient of saturated liquid specific enthalpy fitting relation
Figure BDA0001928065530000101
The input parameter of the saturated gas enthalpy value is pressure
Figure BDA0001928065530000102
Wherein h isg1(p) is the saturated vapor phase enthalpy of the low pressure zone, kJ/kg, hg2(p) is the saturated vapor phase enthalpy of the high pressure zone, kJ/kg, hg3(p) is the saturated vapor phase enthalpy of the interpolation zone, kJ/kg, p is the pressure, MPa, a, b and i are all constant coefficients, and the specific values are shown in Table 2.
TABLE 2 coefficient of fitting relationship of saturated gas specific enthalpy
Figure BDA0001928065530000103
Figure BDA0001928065530000111
S2, dividing the temperature calculation area into an overcooling area, an overheating one area and an overheating two area, wherein the pressure is 0.1-20 MPa, and the specific enthalpy is 200-1600 kJ-kg-1The specific temperature calculation zones are shown in fig. 1.
The range of the supercooling zone is a saturated liquid specific enthalpy line and a lower boundary of a critical specific enthalpy line (h ═ h)set10312.25kJ/kg), the pressure is calculated as follows in the region from 0.1MPa to 20 MPa:
Figure BDA0001928065530000112
wherein, T1(p, h) the temperature of the supercooling zone, DEG C, p the pressure, MPa, h the enthalpy, kJ/kg, aijIs constant, and the specific values are shown in Table 3; (ii) a
TABLE 3 supercooling zone temperature fitting relationship coefficient aij
Figure BDA0001928065530000113
Figure BDA0001928065530000121
The range of the first overheating zone is a saturation gas specific enthalpy line and an upper boundary of a critical specific enthalpy line (h ═ h%set11352.25kJ/kg) or more, 600 kJ/kg-1A region below the specific enthalpy line and at a pressure of from 0.1MPa to 20 MPa;
Figure BDA0001928065530000122
wherein, T2(p, h) is the temperature of the first superheating area, DEG C, p is less than or equal to pcritical，hg(p)＜h≤hset20Or p > pcritical，hset11＜h≤hset20，hset20Is constant, 580kJ/kg, bijThe coefficient is constant, and specific values are shown in Table 4;
TABLE 4 superheat-zone temperature fitting relation coefficient bij
Figure BDA0001928065530000123
The range of the second superheated zone is 620 kJ-kg -11600 kJ.kg above the specific enthalpy line-1The pressure in the region from 0.1MPa to 20MPa below the specific enthalpy line is calculated as follows:
Figure BDA0001928065530000124
wherein, T3(p, h) is the temperature of the second superheating zone, DEG C, h & gt hset21，hset21Is constant, ═ 620kJ/kg, cijThe specific values are shown in Table 5;
TABLE 5 superheat two zone temperature fitting relation coefficient cij
Figure BDA0001928065530000125
Figure BDA0001928065530000131
The calculation formula for interpolating a region is as follows:
Figure BDA0001928065530000132
wherein, T4(p, h) is the temperature of the interpolated region, DEG C, p > pcritical，hset10＜h≤hset11，hset10Is constant, ═ 312.25kJ/kg, hset11Constant, 352.25 kJ/kg;
the calculation formula of the interpolation zone two is as follows:
Figure BDA0001928065530000133
wherein, T5(p, h) is the temperature of the interpolation zone two, DEG C, hset20＜h≤hset21，hset20Is constant, 580kJ/kg, hset21Is constant, 620 kJ/kg.
S3, dividing the specific volume calculation area into an over-cooling area, a state equation area, an over-heating first area, an over-heating second area and a corresponding interpolation area;
the specific volume is calculated within the range of 0.1-20 MPa of pressure and 200-1600 kJ.kg of specific enthalpy-1As shown in fig. 3.
The range of the supercooling zone is a saturated liquid specific enthalpy line and 260 kJ.kg-1A region below the specific enthalpy line at a pressure of from 0.1MPa to 20 MPa;
the state equation region adopts a real gas state equation to carry out iterative solution, and the range comprises two parts, wherein the first part is below a saturated liquid specific enthalpy line and 296 kJ.kg-1Above specific enthalpy line, above saturated gas line and 380 kJ.kg-1A pressure of 7MPa to 7.38MPa below the specific enthalpy line, and a second fraction of 296 kJ-kg-1Above the specific enthalpy line and 380kJ/kg-1Below the specific enthalpy line, the pressure is from 7.38MPa to 20 MPa;
the range of the first superheating area is more than the saturated gas line and 580 kJ.kg-1A region below the specific enthalpy line and at a pressure of from 0.1MPa to 7 MPa; the enthalpy value is 410kJ kg-1Above the specific enthalpy line and 580kJ/kg-1A pressure in a region of 7MPa to 20MPa below the specific enthalpy line.
The range of the second superheat zone is 600 kJ-kg -11600 kJ.kg above the specific enthalpy line-1A region where the pressure is from 0.1MPa to 20MPa below the specific enthalpy line.
The formula for the supercooling region is as follows:
Figure BDA0001928065530000141
wherein v is1(p, h) is the specific volume of the supercooling region, m3/kg，p≤7MPa，h＜hl(p) or p is more than 7MPa, h is less than or equal to hset30，hset30Is constant, 260 kJ/kg; a isijThe specific values are shown in Table 6.
TABLE 6 fitting of specific volume of supercooling region coefficient aij
Figure BDA0001928065530000142
The state equation region is solved by adopting a Hall Motz energy equation:
φ(,τ)＝φo(,τ)+φr(,τ)
where φ is the Hall Moz energy of the actual gas, φoIs the ideal gas portion of the Hall Motz energy, phirFor the remainder of the Hall Motz energy, ═ ρ/ρc，τ＝T/Tc，ρc＝467.6kg/m3，Tc＝31.1℃，p≥pcritical，hset00≤h≤hset01Or p is more than or equal to 7MPa and less than or equal to pcritical，hset00≤h≤hl(p) or p is not less than 7MPa and not more than pcritical，hg(p)≤h≤hset01，hset00Is constant, 296kJ/kg, hset01Is constant, 380kJ/kg, pcriticalCritical pressure of carbon dioxide, 7.31 MPa.
The superheat region is calculated using the following relationship:
Figure BDA0001928065530000151
wherein v is3(p, h) is the specific volume of the overheating zone, p is less than or equal to 7MPa, hg(p)＜h≤hset20Or p > 7MPa, hset31＜h≤hset20，hset31Is constant, 410 kJ/kg; bijThe specific values are shown in Table 7.
TABLE 7 fitting relation coefficient b of specific volume of overheat-zoneij
Figure BDA0001928065530000152
The calculated relationship for superheat region two is as follows:
Figure BDA0001928065530000153
wherein v is4(p, h) is the specific volume of the superheat region two, m3/kg，h＞hset21；cijThe specific values are shown in Table 8.
TABLE 8 fitting relation coefficient a of specific volume of two-zone over-heatij
Figure BDA0001928065530000154
Figure BDA0001928065530000161
An interpolation area:
Figure BDA0001928065530000162
wherein v is5(p, h) is the specific volume of the interpolation area, m3/kg，hset30＝260kJ/kg，hset00＝296kJ/kg，p＞7MPa，hset10＜h≤hset00；
Figure BDA0001928065530000163
Wherein v is6(p, h) isSpecific volume of interpolated two zone, m3/kg，hset31410kJ/kg ofset01＝380kJ/kg，p＞7MPa，hset01＜h≤hset31；
Figure BDA0001928065530000164
Wherein v is7(p, h) is the specific volume of the interpolated three regions, m3/kg，hset20Is constant, 580kJ/kg, hset21Is constant, 620kJ/kg, hset20＜h≤hset21。
S4, calculating the thermal conductivity, wherein the calculation range is 216-1100K, and the pressure is 200 MPa;
λ(ρ,T)＝λ0(T)+Δλ(ρ,T)+Δλc(ρ,T)
wherein λ is0(T) represents the thermal conductivity at zero density of the lean gas, Δ λ (ρ, T) represents the increase in thermal conductivity due to the increase in density of the lean gas at the same temperature, Δ λc(ρ, T) represents the increase in thermal conductivity in the peak region at the critical point.
And S5, calculating the dynamic viscosity.
The calculation range is 200-1500K, the density is 1400kg/m3
η(ρ,T)＝η0(T)+Δη(ρ,T)+Δηc(ρ,T)
Wherein eta is0(T) represents the viscosity at zero density of the lean gas, Δ η (ρ, T) represents the increase in viscosity due to the increase in density of the lean gas at the same temperature, Δ ηc(ρ, T) represents the increase in viscosity in the peak region at the critical point.
In order to make the objects, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the drawings in the embodiments of the present invention, and it is obvious that the described embodiments are some, but not all, embodiments of the present invention. The components of the embodiments of the present invention generally described and illustrated in the figures herein may be arranged and designed in a wide variety of different configurations. Thus, the following detailed description of the embodiments of the present invention, presented in the figures, is not intended to limit the scope of the invention, as claimed, but is merely representative of selected embodiments of the invention. All other embodiments, which can be derived by a person skilled in the art from the embodiments given herein without making any creative effort, shall fall within the protection scope of the present invention.
For saturated liquid specific enthalpy calculation, 6777 pressure points are taken by taking 1kPa as a step length to perform saturated liquid specific enthalpy calculation under 0.6-7.377 MPa, and the maximum relative error calculated by the method is 0.08%, as shown in FIG. 3; the calculation speed is 118 times that of commercial software NIST PREPROP as shown in fig. 4.
For saturated gas specific enthalpy calculation, 6777 pressure points are taken by taking 1kPa as a step length to perform saturated liquid specific enthalpy calculation under 0.6-7.377 MPa, and the maximum relative error calculated by the method is 0.009%, as shown in FIG. 5; the calculation speed is 230 times that of commercial software NIST PREPROP as shown in fig. 6.
For the temperature calculation, the pressure range is 0.1-20 MPa, the step length is 0.1MPa, and the specific enthalpy range is 200--1Step length of 2 kJ/kg-1131914 state points are taken for calculation, and the maximum relative error calculated by the method is 0.25 percent, as shown in FIG. 7; the calculation speed is 15 times that of commercial software NIST PREPROP as shown in fig. 8.
For specific volume calculation, the pressure range is 0.1-20 MPa, the step length is 0.1MPa, and the specific enthalpy range is 200--1Step length of 2 kJ/kg-1And 131914 state points are taken for calculation, the maximum relative error of the calculation of the method is within 1 percent, as shown in FIG. 9, and the calculation speed is 28 times of that of commercial software NIST PREPROP, as shown in FIG. 10.
For the calculation of the thermal conductivity, the maximum relative error is within 5%, and the calculation speed is equivalent to the commercial software NIST PREPROP.
For the calculation of dynamic viscosity, the maximum relative error is within 5%, and the calculation speed is equivalent to commercial software NISTPREPROP.
The above-mentioned contents are only for illustrating the technical idea of the present invention, and the protection scope of the present invention is not limited thereby, and any modification made on the basis of the technical idea of the present invention falls within the protection scope of the claims of the present invention.

Claims (10)

1. A wide parameter range carbon dioxide physical property solving method is characterized in that iteration solving is carried out in a region of a carbon dioxide critical point by adopting a real gas state equation; solving the physical property in the region far away from the critical point of the carbon dioxide by adopting a fitting relation; and (4) performing interpolation solution in the boundary area by adopting a fitting relation and a real gas state equation, and calculating to obtain saturated liquid enthalpy value, saturated gas enthalpy value, temperature, specific volume, thermal conductivity and dynamic viscosity physical property parameters.
2. The wide parameter range carbon dioxide solving method of claim 1, wherein the saturated liquid phase enthalpy value is solved using a pressure-dependent polynomial fit relation as follows:
Figure FDA0002573610710000011
wherein h isl1(p) saturated liquid phase enthalpy of the low pressure zone, hl2(p) is the saturated liquid phase enthalpy of the high pressure zone, hl3And (p) is the saturated liquid phase enthalpy value of the interpolation zone, p is the pressure, MPa, a, b and i are all constant coefficients.
3. The wide parameter range carbon dioxide solving method of claim 1, wherein the saturated vapor phase enthalpy is solved using a pressure-dependent polynomial fit relation as follows:
Figure FDA0002573610710000012
wherein h isg1(p) the saturated vapor phase enthalpy value g, h of the low-pressure zoneg2(p) is the saturated vapor phase enthalpy of the high pressure zone, hg3And (p) is the saturated vapor phase enthalpy value of the interpolation zone, p is the pressure, MPa, a, b and i are all constant coefficients.
4. The method for solving the physical property of the carbon dioxide in the wide parameter range according to claim 1, wherein the temperature calculation area is divided into an overcooling area, an overheating one area and an overheating two area, the pressure is 0.1-20 MPa, and the specific enthalpy is 200-1600 kJ-kg-1。
5. The method for solving the physical property of carbon dioxide in the wide parameter range according to claim 4, wherein the range of the supercooling region is 312.25 kJ-kg under a saturated liquid specific enthalpy line and a critical specific enthalpy line-1A region under a pressure of from 0.1MPa to 20 MPa;
the range of the first overheating zone is saturation gas specific enthalpy line and upper boundary 352.25 kJ.kg of critical specific enthalpy line-1Above, 600 kJ.kg-1A region below the specific enthalpy line and at a pressure of from 0.1MPa to 20 MPa;
the range of the second superheat zone is 600 kJ-kg-11600 kJ.kg above the specific enthalpy line-1A region where the pressure is from 0.1MPa to 20MPa below the specific enthalpy line.
6. The method for solving the wide parameter range carbon dioxide physical property as claimed in claim 5, wherein the temperature T of the supercooling region1(p, h) is calculated as follows:
Figure FDA0002573610710000021
wherein, aijIs constant coefficient, p is pressure, MPa, h is enthalpy, kJ/kg, and p is not more than pcritical，h＜hl(p) or p > pcritical，h≤hset10，hset10＝hcritical-20.0＝312.25kJ/kg；
Temperature T of superheat zone2(p, h) is calculated as follows:
Figure FDA0002573610710000022
wherein, bijIs a constant coefficient, p is less than or equal to pcritical，hg(p)＜h≤hset20Or p > pcritical，hset11＜h≤hset20，hset11＝hcritical+20.0＝352.25kJ/kg，hset20Is a constant number hset20＝580kJ/kg；
Temperature T of the superheat zone two3(p, h) is calculated as follows:
Figure FDA0002573610710000023
wherein, cijIs a constant coefficient, h > hset21，hset21Is a constant number hset21＝620kJ/kg；
The interpolation zone is the transition zone between the supercooling zone and the superheating zone, and has a temperature T4(p, h) is calculated as follows:
Figure FDA0002573610710000024
wherein h is the enthalpy value, p > pcritical，hset10＜h≤hset11；
The second interpolation zone is a transition zone between the first superheating zone and the second superheating zone, and the temperature T of the second interpolation zone is5(p, h) is calculated as follows:
Figure FDA0002573610710000031
7. the method for solving the physical property of the carbon dioxide in the wide parameter range according to claim 1, wherein the calculation area of the specific volume is divided into an overcooling area, an equation of state area, an overheating one area and an overheating two area, the pressure is 0.1-20 MPa, and the specific enthalpy is 200-1600 kJ-kg-1。
8. The method of claim 7, wherein the range of the supercooling region is a saturated liquid specific enthalpy line and 260 kJ-kg-1A region below the specific enthalpy line at a pressure of from 0.1MPa to 20 MPa;
the state equation region is solved by using a real gas state equation, and the range comprises two parts, wherein the first part is below a saturated liquid specific enthalpy line and 296 kJ.kg-1Above specific enthalpy line, above saturated gas line and 380 kJ.kg-1A pressure of 7MPa to 7.38MPa below the specific enthalpy line, and a second fraction of 296 kJ-kg-1Above the specific enthalpy line and 380kJ/kg-1Below the specific enthalpy line, the pressure is from 7.38MPa to 20 MPa;
the range of the first superheating area is more than the saturated gas line and 580 kJ.kg-1A region below the specific enthalpy line and at a pressure of from 0.1MPa to 7 MPa; the enthalpy value is 410kJ kg-1Above the specific enthalpy line and 580kJ/kg-1A region having a pressure of 7MPa to 20MPa below a specific enthalpy line;
the range of the second superheat zone is 600 kJ-kg-11600 kJ.kg above the specific enthalpy line-1A region where the pressure is from 0.1MPa to 20MPa below the specific enthalpy line.
9. The method of claim 8, wherein the specific volume v of the supercooling region is set as the specific volume v1(p, h) is calculated as follows:
Figure FDA0002573610710000032
wherein, aijA constant coefficient, p is less than or equal to 7MPa, h is less than hl(p) or p is more than 7MPa, h is less than or equal to hset30，hset30Is a constant number hset30＝260kJ/kg；
The equation of state region is calculated using the Hall Motz function:
φ(,τ)＝φo(,τ)+φr(,τ)
where φ is the Hall Moz energy of the actual gas, φoIs the ideal gas portion of the Hall Motz energy, phirFor the remainder of the Hall Motz energy, ═ ρ/ρc，τ＝T/Tc，ρc＝467.6kg/m3，Tc＝31.1℃，p≥pcritical，hset00≤h≤hset01Or p is more than or equal to 7MPa and less than or equal to pcritical，hset00≤h≤hl(p) or p is not less than 7MPa and not more than pcritical，hg(p)≤h≤hset01，hset00Is a constant number hset00＝296kJ/kg，hset01Is a constant number hset01＝380kJ/kg，pcriticalIs the critical pressure of carbon dioxide, pcritical＝7.31MPa；
Specific volume v of superheat region3(p, h) is calculated as follows:
Figure FDA0002573610710000041
wherein, bijIs a constant coefficient, p is less than or equal to 7MPa, hg(p)＜h≤hset20Or p > 7MPa, hset31＜h≤hset20，hset31Is a constant number hset31＝410kJ/kg；
Specific volume v of the overheated second zone4(p, h) is calculated as follows:
Figure FDA0002573610710000042
wherein, cijIs a constant coefficient, h > hset21；
The interpolation zone is a transition zone of the supercooling zone and the equation of state zone and has a specific volume v5(p, h) is calculated as follows:
Figure FDA0002573610710000043
wherein h isset30＝260kJ/kg，hset00＝296kJ/kg，p＞7MPa，hset30＜h≤hset00；
The interpolation two zone is a transition zone of the equation of state zone and the overheating one zone and has a specific volume v6(p, h) is calculated as follows:
Figure FDA0002573610710000044
wherein h isset31＝410kJ/kg，hset01＝380kJ/kg，p＞7MPa，hset01＜h≤hset31；
The three interpolation zones are transition zones of the first overheating zone and the second overheating zone and have specific volume v7(p, h) is calculated as follows:
Figure FDA0002573610710000051
wherein h isset20＝580kJ/kg，hset21＝620kJ/kg，hset20＜h≤hset21。
10. The wide parameter range carbon dioxide physical property solving method according to claim 1, wherein the thermal conductivity is calculated in the temperature range of 216-1100K and the pressure range of 200MPa as follows:
λ(ρ,T)＝λ0(T)+Δλ(ρ,T)+Δλc(ρ,T)
wherein λ is0(T) represents the thermal conductivity at zero density of the lean gas, Δ λ (ρ, T) represents the increase in thermal conductivity due to the increase in density of the lean gas at the same temperature, Δ λc(ρ, T) represents the increase in thermal conductivity in the critical point peak region;
the dynamic viscosity is calculated in the temperature range of 200-1500K and the density range of 1400kg/m3The method comprises the following steps:
η(ρ,T)＝η0(T)+Δη(ρ,T)+Δηc(ρ,T)
wherein eta is0(T) represents the viscosity at zero density of the lean gas, Δ η (ρ, T) represents the increase in viscosity due to the increase in density of the lean gas at the same temperature, Δ ηc(ρ, T) represents the increase in viscosity in the peak region at the critical point.
CN201811626329.3A 2018-12-28 2018-12-28 Wide parameter range carbon dioxide physical property solving method Active CN109783904B (en)

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 CN109783904A (en) 2019-05-21
CN109783904B true 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)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113490391B (en) * 2021-05-27 2024-06-21 合肥通用机械研究院有限公司 Rectangular micro-channel unit one-dimensional temperature distribution calculation method and system for electronic cooling
CN116030906B (en) * 2023-03-28 2023-06-27 西安交通大学 Physical property parameter calculation method for lead bismuth fast reactor multicomponent fluid

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104777183A (en) * 2014-10-31 2015-07-15 北京卫星环境工程研究所 Satellite electric propulsion system xenon filling thermodynamic characteristic numerical simulation method
CN108134156A (en) * 2017-11-28 2018-06-08 江西爱驰亿维实业有限公司 The calculation method of parameters of tube refrigerant fluid interchange, system, medium, terminal, battery pack

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP4980198B2 (en) * 2007-10-26 2012-07-18 日東電工株式会社 Foam material with slit for mounting step structure
CN107506534B (en) * 2017-08-04 2021-02-19 陕西延长石油（集团）有限责任公司 Method for evaluating closure of cover layer in carbon dioxide flooding and sequestration
CN108763670B (en) * 2018-05-15 2021-04-20 西安交通大学 Method for solving Brayton cycle transient process of supercritical carbon dioxide reactor

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104777183A (en) * 2014-10-31 2015-07-15 北京卫星环境工程研究所 Satellite electric propulsion system xenon filling thermodynamic characteristic numerical simulation method
CN108134156A (en) * 2017-11-28 2018-06-08 江西爱驰亿维实业有限公司 The calculation method of parameters of tube refrigerant fluid interchange, system, medium, terminal, battery pack

Also Published As

Publication number Publication date
CN109783904A (en) 2019-05-21

Similar Documents

Publication Publication Date Title
CN109783904B (en) Wide parameter range carbon dioxide physical property solving method
Scaccabarozzi et al. Thermodynamic optimization and part-load analysis of the NET power cycle
CN107368680A (en) A kind of steam turbine optimum vacuum real-time computing technique
CN112464472B (en) Method for improving heat exchange calculation performance of steam generator in sodium-cooled fast reactor system program
CN110925741B (en) Method for acquiring inlet and outlet temperatures of main steam superheaters based on desuperheating water quantity
CN111651851A (en) Containment solving method and containment solver
CN109033724B (en) Main steam temperature consumption difference correction curve optimization method applied to sliding pressure operation of steam turbine
Li et al. Triple-objective optimization of He Brayton cycles for ultra-high-temperature solar power tower
Chan et al. Defining the thermodynamic efficiency in a wave rotor
Alsagri et al. Performance comparison and parametric analysis of sCO2 power cycles configurations
CN104483152B (en) The heat consumption rate assay method of non-reheat backheat combined-circulation unit
CN117725787A (en) Calculation method for non-fluid phase-change molten salt heat exchange based on fin heat exchanger
Li et al. Technical benefits of the subcritical inlet condition for high-speed CO2 centrifugal compressor in the advanced power-generation cycle
CN105911088A (en) Comparative method for condenser saturation pressures under different thermal loads
CN113153570B (en) Pulse detonation tube performance calculation method and device
Mohammed et al. THERMOECONOMIC OPTIMIZATION OF TRIPLE PRESSURE HEAT RECOVERY STEAM GENERATOR OPERATING PARAMETERS FOR COMBINED CYCLE PLANTS.
CN113868887A (en) Simulation calculation method for water washing tower
Armijo et al. System Design of a 2.0 MWth Sodium/Molten salt pilot system
JP6067594B2 (en) Steam turbine plant
CN106610917B (en) The novel calculation method that coolant temperature changes in reactor core runner when heap is opened in simulation
Aktershev et al. Numerical modeling of the vapor bubble growth in a homogenously superheated liquid (thermal energy scheme)
CN113821998B (en) Method for solving shell side pressure of condenser real-time dynamic simulation model by Newton iteration method
Salah et al. Off-design performance assessment of an axial turbine for a 100 MWe concentrated solar power plant operating with CO2 mixtures
CN112360810B (en) Impeller inlet design method of supercritical carbon dioxide centrifugal compressor
CN118211385A (en) Phase change coupling method for condensing and evaporating vapor containing non-condensable gas

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