CN101716427B - System and method for simulating dynamic flow of internal thermally coupled air separation column - Google Patents

Abstract

A system for simulating dynamic flow of an internal thermally coupled air separation column comprises a field intelligent instrument connected with the internal thermally coupled air separation column, a control station, a database and an upper computer. The upper computer comprises a signal acquisition module used for acquiring the current production condition data, and a main solving and computing module. The main solving and computing module adopts the following processes: setting the structural parameter and the operation parameter of the column and setting the start time tstart and the end time tend; specifying the liquid composition and liquid flow rate of each column plate at initial time and supposing the current iteration time t=tstart; respectively computing the balance temperature and vapor composition of each column plate by a bubble point method and computing vapor-liquid enthalpy of each column plate; computing the liquid composition and liquid flow rate of each column plate at (t+deltat); and supposing t=t+deltat, using new liquid composition and liquid flow rate of each column plate to return for iteration until t is not less than tend, completing iteration and outputting the result. The invention also provides a method for simulating dynamic flow of the internal thermally coupled air separation column. The system and the method can quickly and accurately simulate the dynamic flow of the internal thermally coupled air separation column.

Description

Dynamic flow of internal thermally coupled air separation column simulation system and method

Technical field

The present invention relates to empty branch field, especially, relate to a kind of dynamic flow of internal thermally coupled air separation column simulation system and method.

Background technology

The application of oxygen, nitrogen and argon gas is very extensive.Cryogenic air separation process is a current domestic and international air separation sector application method the most widely.In the air separation industry, energy cost has accounted for 75% of air products price.Under the situation that energy crisis is constantly deepened, the energy efficiency that improves air separation technology has important society and economic implications.

Internal thermal coupled distillation technology is a highest power-saving technology of energy-conservation usefulness in up to now the four big energy-saving rectifying technology at first, has worldwide obtained huge attention.In the cryogenic air separation process, air separation column is an important operating unit, also is most important power consumption unit.By internal thermal coupled Study on Technology is found: the realization of this technology need to have in the process apparatus the different just zones of pressure, and the technological process of conventional air separation unit itself has just provided such zone; Studies show that internal thermal coupled technology will have bigger economic benefit on the high purity product production process, and air separation process is a high process of product purity requirement just; Internal thermal coupled distillation technology helps bringing into play its energy-saving potential more under the low operating temperature situation, and conventional air separation unit itself is exactly a cryogenic rectification; Internal thermal coupled distillation technology can embody very high energy-saving efficiency, and nitrogen in the air--argon--oxygen ternary system has satisfied this requirement just when handling the approaching system of boiling point, also improved highly beneficial condition for the exploitation of thermal coupling The Application of Technology.Therefore, internal thermal coupled technology can be applied to air separation process, changes traditional air separating tower structure, reaches good energy-saving effect.

Traditional research method of Chemical Engineering is to reduce the master with experience, promptly uses the relation between the laboratory facilities searching system variable.Along with the maximization of production process and the raising of automatization level, for the research of complicated chemical process, dimensional analysis and similarity method often can not satisfy the needs of research.The internal thermally coupled air separation column process is to include heat transfer, mass transfer, the mobile complex process of fluid, present the strong nonlinear relation between many variablees, rule over time such as CONCENTRATION DISTRIBUTION, VELOCITY DISTRIBUTION in will the understanding process, an effective method utilizes computer to carry out the dynamic flow simulation exactly.By the dynamic characteristic and the response of dynamic flow of internal thermally coupled air separation column sunykatuib analysis air separation column, be the prerequisite of flow scheme design and control.

Summary of the invention

For the deficiency of the accurately simulated interior thermally coupled air separation column dynamic flow of the internal thermal coupled distillation process that overcomes existing air separation column, the invention provides a kind of dynamic flow of internal thermally coupled air separation column simulation system and method for accurately simulated interior thermally coupled air separation column dynamic flow.

The technical solution adopted for the present invention to solve the technical problems is:

A kind of dynamic flow of internal thermally coupled air separation column simulation system, comprise the field intelligent instrument and control station, database and the host computer that are connected with internal thermally coupled air separation column, intelligence instrument is connected with control station, database and host computer, described host computer comprises: signal acquisition module, in order to gather current production status data; Find the solution the calculating primary module,, adopt following process to finish in order to find the solution calculating:

1) structural parameters and the operating parameter of setting tower are set initial moment tstart, stop tend constantly;

2) specify each column plate liquid phase of initial time to form and the liquid phase flow, make current iteration time t=tstart;

3), calculate its equilibrium temperature by the bubble point method respectively and vapour phase is formed to each column plate;

4), calculate the enthalpy of vapour-liquid phase respectively to each column plate;

5) calculate the vapour-liquid phase flow rate of each column plate by formula (1) (2):

$V_{j+1}{H}_{j+1}^G+U_{j-1}{H}_{j-1}^L+F_j{H}_j^F-(V_j+S_j^G)H_j^G-(U_j+S_j^L)H_j^L-Q_j=0---\left(1\right)$

$V_{j+1}+U_{j-1}+F_j^G+F_j^L-(V_j+S_j^G)-(U_j+S_j^L)=0---\left(2\right)$

Wherein, V represents the vapour phase flow, and U represents the liquid phase flow, and F represents feed rate, H ^FExpression charging enthalpy, S represent that side carries flow, H ^GAnd H ^LBe respectively vapour-liquid phase enthalpy, subscript j-1, j, j+1 represent j-1, j, j+1 piece plate respectively, and subscript L represents liquid phase, and subscript G represents vapour phase, and Q represents the thermal coupling amount, is calculated by following formula:

Q＝UAΔT (3)

Wherein, UA represents heat couple coefficient, the Δ T temperature difference between column plate of representing to be coupled;

6) calculating (t+ Δ t) each column plate liquid phase constantly forms and the liquid phase flow:

$\frac{{\mathrm{dx}}_{i,j}}{\mathrm{dt}}{M}_j=U_{j-1}{x}_{i,j-1}+V_{j+1}{y}_{i,j+1}+F_j{z}_{i,j}-(U_j+S_j^L)x_{i,j}-(V_j+S_j^G)y_{i,j}----\left(4\right)$

$[U_{j-1}+V_{j+1}+F_j-(U_j+S_j^L)-(V_j+S_j^G)]x_{i,j}$

$\frac{d{U}_j}{\mathrm{dt}}=518.13{\mathrm{\ρ}}_j^{-1/3}{A}_a^{-1}{l}^{2/3}{U}_j^{1/3}[U_{j-1}+V_{j+1}+F_j-(U_j+S_j^L)-(V_j+S_j^G)]---\left(5\right)$

$x_{i,j}^{(t+\mathrm{\Δt})}=x_{i,j}^{\left(t\right)}+\frac{d{x}_{i,j}}{\mathrm{dt}}\mathrm{\Δt}---\left(6\right)$

$U_j^{(t+\mathrm{\Δt})}=U_j^{\left(t\right)}+\frac{d{U}_j}{\mathrm{dt}}\mathrm{\Δt}---\left(7\right)$

Wherein, Δ t is an iteration step length, and x is that liquid phase is formed, and y is that vapour phase is formed, and z is a feed composition, and ρ is a density of liquid phase, A _aBe the column plate effective area, l is that the column plate weir is long, subscript i=1,2,3 expression components, and corresponding successively nitrogen, argon, oxygen, subscript (t) and (t+ Δ t) are represented t and t+ Δ t constantly respectively, M represents the column plate liquid holdup, is calculated by following formula:

$M_j={\mathrm{\ρ}}_j{A}_{\mathrm{aj}}[h_{\mathrm{wj}}+0.00284{\left(\frac{U_j}{{\mathrm{\ρ}}_j{l}_j}\right)}^{2/3}]---\left(8\right)$

H wherein _wIt is the overflow height of weir;

7) make t=t+ Δ t, return 3 with new each column plate liquid phase composition and liquid phase flow) iteration, up to t 〉=tend, finishing iteration, output result.

As preferred a kind of scheme: described host computer also comprises: bubble point method module, and in order to calculate its equilibrium temperature by the bubble point method and vapour phase is formed, its process is as follows:

3.1) supposition column plate equilibrium temperature;

3.2) calculate the VLE constant, adopt following process to finish:

$\ln {\mathrm{\Φ}}_i^L=\ln \frac{\mathrm{RT}}{P(v^L-b^L)}-\frac{b_i}{b^L}(1-Z^L)+{\mathrm{\ξ}}^L{a}^L(\frac{b_i}{b^L}-\frac{2\underset{m}{\mathrm{\Σ}}{x}_m{a}_{i,m}}{a^L})/b^L\mathrm{RT}---\left(9\right)$

$\ln {\mathrm{\Φ}}_i^G=\ln \frac{\mathrm{RT}}{P(v^G-b^G)}-\frac{b_i}{b^G}(1-Z^G)+{\mathrm{\ξ}}^G{a}^G(\frac{b_i}{b^G}-\frac{2\underset{m}{\mathrm{\Σ}}{x}_m{a}_{i,m}}{a^G})/b^G\mathrm{RT}---\left(10\right)$

$K_i={\mathrm{\Φ}}_i^L/{\mathrm{\Φ}}_i^G---\left(11\right)$

y _i＝K _ix _i (12)

Wherein, Φ represents fugacity coefficient, and subscript G represents vapour phase, and Q represents the thermal coupling amount, and R is a gas constant, and T is a temperature, and P is a column plate pressure, subscript m=1,2,3 expression components, corresponding successively nitrogen, argon, oxygen, molal volume v, parameter b ^G, b ^L, b _i, a ^G, a ^L, a _{I, m}, ξ ^G, ξ ^L, vapour phase compressibility factor Z ^G, liquid phase compressibility factor Z ^LCalculate by the rerum natura module;

3.3) check

Whether set up, set up then finishing iteration, return result of calculation, otherwise, upgrade the column plate equilibrium temperature, return 3.2) the continuation iteration.

As preferred another kind of scheme: described host computer also comprises: the enthalpy module, and in order to calculate vapour-liquid phase enthalpy of mixing, its process is as follows:

$H_i^{*}=c_i+d_{i}T+e_i{T}^2+f_i{T}^3+h_i{T}^4---\left(13\right)$

$H^{*}=\underset{i}{\mathrm{\Σ}}{y}_i{H}_i^{*}---\left(14\right)$

$H^G=H^{*}-\mathrm{RT}(1-Z^G)-{\mathrm{\ξ}}^G(a^G-T\frac{d{a}^G}{\mathrm{dT}})/b^G---\left(15\right)$

$H^L=H^{*}-\mathrm{RT}(1-Z^L)-{\mathrm{\ξ}}^L(a^L-T\frac{d{a}^L}{\mathrm{dT}})/b^L---\left(16\right)$

H wherein _i ^*The enthalpy of representing i pure component perfect gas, H ^*Be mixture perfect gas enthalpy, c, d, e, f, h are constant.

As preferred another scheme: described host computer also comprises: the rerum natura module, and in order to calculate physical parameter, its process is as follows:

$a_{i,m}={\mathrm{\Ω}}_{\mathrm{ai},m}{R}^2{T}_{\mathrm{ci},m}^2/P_{\mathrm{ci},m}---\left(17\right)$

b _i＝Ω _bRT _ci/P _cia (18)

$T_{\mathrm{ci},m}=\sqrt{T_{\mathrm{ci}}{T}_{\mathrm{cj}}}(1-k_{i,m})---\left(19\right)$

$V_{\mathrm{ci},m}=0.125{(V_{\mathrm{ci},m}^{1/3}+V_{\mathrm{ci},m}^{1/3})}^3---\left(20\right)$

Z _ci，m＝0.5(Z _ci+Z _cm) (21)

P _ci，m＝RT _ci，mZ _ci，m/V _ci，m (22)

Ω _ai，m＝0.5(Ω _ai+Ω _am) (23)

To vapour phase:

$a^G=\underset{i}{\mathrm{\Σ}}\underset{m}{\mathrm{\Σ}}{y}_i{y}_m{a}_{i,m}---\left(24\right)$

$b^G=\underset{i}{\mathrm{\Σ}}{y}_i{b}_i---\left(25\right)$

Order

A ^G＝a ^GP/R ²T ² (26)

B ^G＝b ^GP/RT (27)

α ^G＝2B ^G-1 (28)

${\mathrm{\β}}^G=A^G-3B^G-5B^{G^2}---\left(29\right)$

${\mathrm{\γ}}^G=2(B^{G^3}+B^{G^2})-A^G{B}^G---\left(30\right)$

Getting initial value is 1-0.6P _r, separate following equation with Newton method, promptly obtain vapour phase compressibility factor Z ^G

$Z^{G^3}+{\mathrm{\α}}^G{Z}^{G^2}+\mathrm{\β}{Z}^G+{\mathrm{\γ}}^G=0---\left(31\right)$

Then,

v ^G＝RT/PZ ^G (32)

${\mathrm{\ξ}}^G=0.242536\ln \frac{v^G+3.561553b^G}{v^G-0.561553b^G}---\left(33\right)$

To liquid phase:

$a^L=\underset{i}{\mathrm{\Σ}}\underset{m}{\mathrm{\Σ}}{x}_i{x}_m{a}_{i,m}---\left(34\right)$

$b^L=\underset{i}{\mathrm{\Σ}}{x}_i{b}_i---\left(35\right)$

Order

A ^L＝a ^LP/R ²T ² (36)

B ^L＝b ^LP/RT (37)

α ^L＝2B ^L-1 (38)

${\mathrm{\β}}^L=A^L-3B^L-5B^{L^2}---\left(39\right)$

${\mathrm{\γ}}^L=2(B^{L^3}+B^{L^2})-A^L{B}^L---\left(40\right)$

Getting initial value is P _r(0.106+0.078P _r), separate following equation with Newton method, promptly obtain liquid phase compressibility factor Z ^L

$Z^{L^3}+{\mathrm{\α}}^L{Z}^{L^2}+\mathrm{\β}{Z}^L+{\mathrm{\γ}}^L=0---\left(41\right)$

Then,

v ^L＝RT/PZ ^L (42)

${\mathrm{\ξ}}^L=0.242536\ln \frac{v^L+3.561553b^L}{v^L-0.561553b^L}---\left(43\right)$

Ω _ai＝C _i-D _iτ+E _iτ ²-W _iτ ³ (44)

Ω _b＝0.070721 (45)

τ＝0.01T (46)

Wherein, A, B, α, β, γ, τ are intermediate variables, and C, D, E, W are constants, T _c, P _c, V _c, Z _cBe respectively critical-temperature, pressure, volume and compressibility factor, P _rBe reduced pressure, R is a gas constant, k _{I, m}The mutual coefficient of binary of representing i component and m component, k _{I, m}Be constant, subscript c represents the character of critical point, and subscript r represents reduced state, subscript i, and m represents the binary mixture of i component and m component, Ω _a, Ω _bIt is intermediate variable.

Further, described host computer also comprises: display module as a result is used for that result of calculation is passed to control station and shows, and by fieldbus result of calculation is delivered to operator station and shows.

A kind of dynamic flow of internal thermally coupled air separation column analogy method, described flowsheeting method may further comprise the steps:

1) structural parameters of setting tower are gathered the production status data, set initial moment tstart, stop tend constantly;

2) specify each column plate liquid phase of initial time to form and the liquid phase flow, make current iteration time t=tstart;

3), calculate its equilibrium temperature by the bubble point method respectively and vapour phase is formed to each column plate;

4), calculate the enthalpy of its vapour-liquid phase to each column plate;

5) simultaneous formula (1) (2) is calculated the vapour-liquid phase flow rate of each column plate:

$V_{j+1}{H}_{j+1}^G+U_{j-1}{H}_{j-1}^L+F_j{H}_j^F-(V_j+S_j^G)H_j^G-(U_j+S_j^L)H_j^L-Q_j=0---\left(1\right)$

$V_{j+1}+U_{j-1}+F_j^G+F_j^L-(V_j+S_j^G)-(U_j+S_j^L)=0---\left(2\right)$

Wherein, V represents the vapour phase flow, and U represents the liquid phase flow, and F represents feed rate, H ^FExpression charging enthalpy, S represent that side carries flow, H ^GAnd H ^LBe respectively vapour-liquid phase enthalpy, subscript j-1, j, j+1 represent j-1, j, j+1 piece plate respectively, and subscript L represents liquid phase, and subscript G represents vapour phase, and Q represents the thermal coupling amount, is calculated by following formula:

Q＝UAΔT (3)；

Wherein, UA represents heat couple coefficient, the Δ T temperature difference between column plate of representing to be coupled;

6) calculating (t+ Δ t) each column plate liquid phase constantly forms and the liquid phase flow:

$\frac{{\mathrm{dx}}_{i,j}}{\mathrm{dt}}{M}_j=U_{j-1}{x}_{i,j-1}+V_{j+1}{y}_{i,j+1}+F_j{z}_{i,j}-(U_j+S_j^L)x_{i,j}-(V_j+S_j^G)y_{i,j}----\left(4\right)$

$[U_{j-1}+V_{j+1}+F_j-(U_j+S_j^L)-(V_j+S_j^G)]x_{i,j}$

$\frac{d{U}_j}{\mathrm{dt}}=518.13{\mathrm{\ρ}}_j^{-1/3}{A}_a^{-1}{l}^{2/3}{U}_j^{1/3}[U_{j-1}+V_{j+1}+F_j-(U_j+S_j^L)-(V_j+S_j^G)]---\left(5\right)$

$x_{i,j}^{(t+\mathrm{\Δt})}=x_{i,j}^{\left(t\right)}+\frac{d{x}_{i,j}}{\mathrm{dt}}\mathrm{\Δt}---\left(6\right)$

$U_j^{(t+\mathrm{\Δt})}=U_j^{\left(t\right)}+\frac{d{U}_j}{\mathrm{dt}}\mathrm{\Δt}---\left(7\right)$

Wherein, Δ t is an iteration step length, and x is that liquid phase is formed, and y is that vapour phase is formed, and z is a feed composition, and ρ is a density of liquid phase, A _aBe the column plate effective area, l is that the column plate weir is long, subscript i=1,2,3 expression components, and corresponding successively nitrogen, argon, oxygen, subscript (t) and (t+ Δ t) are represented t and t+ Δ t constantly respectively, M represents the column plate liquid holdup, is calculated by following formula:

$M_j={\mathrm{\ρ}}_j{A}_{\mathrm{aj}}[h_{\mathrm{wj}}+0.00284{\left(\frac{U_j}{{\mathrm{\ρ}}_j{l}_j}\right)}^{2/3}]---\left(8\right)$

H wherein _wIt is the overflow height of weir;

7) make t=t+ Δ t, return the step 3) iteration with new each column plate liquid phase composition and liquid phase flow, up to t 〉=tend, finishing iteration, output result.

As preferred a kind of scheme: in the described step 3), the bubble point method is calculated its equilibrium temperature and vapour phase, adopts following process to finish:

3.1) supposition column plate equilibrium temperature;

3.2) calculate the VLE constant, adopt following process to finish:

$\ln {\mathrm{\Φ}}_i^L=\ln \frac{\mathrm{RT}}{P(v^L-b^L)}-\frac{b_i}{b^L}(1-Z^L)+{\mathrm{\ξ}}^L{a}^L(\frac{b_i}{b^L}-\frac{2\underset{m}{\mathrm{\Σ}}{x}_m{a}_{i,m}}{a^L})/b^L\mathrm{RT}---\left(9\right)$

$\ln {\mathrm{\Φ}}_i^G=\ln \frac{\mathrm{RT}}{P(v^G-b^G)}-\frac{b_i}{b^G}(1-Z^G)+{\mathrm{\ξ}}^G{a}^G(\frac{b_i}{b^G}-\frac{2\underset{m}{\mathrm{\Σ}}{x}_m{a}_{i,m}}{a^G})/b^G\mathrm{RT}---\left(10\right)$

$K_i={\mathrm{\Φ}}_i^L/{\mathrm{\Φ}}_i^G---\left(11\right)$

y _i＝K _ix _i (12)

Wherein, Φ represents fugacity coefficient, and subscript L represents liquid phase, and subscript G represents vapour phase, and R is a gas constant, and T is a temperature, and P is a column plate pressure, subscript m=1,2,3 expression components, corresponding successively nitrogen, argon, oxygen, molal volume v, parameter b ^G, b ^L, b _i, a ^G, a ^L, a _{I, m}, ξ ^G, ξ ^L, vapour phase compressibility factor Z ^G, liquid phase compressibility factor Z ^LCalculate by the physical parameter computational methods;

3.3) check

Whether set up, set up then finishing iteration, return result of calculation, otherwise, upgrade the column plate equilibrium temperature, return 3.2) the continuation iteration.

As preferred another kind of scheme: in the described step 4), described enthalpy computational methods process is as follows:

$H_i^{*}=c_i+d_{i}T+e_i{T}^2+f_i{T}^3+h_i{T}^4---\left(13\right)$

$H^{*}=\underset{i}{\mathrm{\Σ}}{y}_i{H}_i^{*}---\left(14\right)$

$H^G=H^{*}-\mathrm{RT}(1-Z^G)-{\mathrm{\ξ}}^G(a^G-T\frac{d{a}^G}{\mathrm{dT}})/b^G---\left(15\right)$

$H^L=H^{*}-\mathrm{RT}(1-Z^L)-{\mathrm{\ξ}}^L(a^L-T\frac{d{a}^L}{\mathrm{dT}})/b^L---\left(16\right)$

H wherein _i ^*The enthalpy of representing i pure component perfect gas, H ^*Be mixture perfect gas enthalpy, c, d, e, f, h are constant.

As preferred another scheme: the process of described physical parameter computational methods is as follows:

$a_{i,m}={\mathrm{\Ω}}_{\mathrm{ai},m}{R}^2{T}_{\mathrm{ci},m}^2/P_{\mathrm{ci},m}---\left(17\right)$

b _i＝Ω _bRT _ci/P _cia (18)

$T_{\mathrm{ci},m}=\sqrt{T_{\mathrm{ci}}{T}_{\mathrm{cj}}}(1-k_{i,m})---\left(19\right)$

$V_{\mathrm{ci},m}=0.125{(V_{\mathrm{ci},m}^{1/3}+V_{\mathrm{ci},m}^{1/3})}^3---\left(20\right)$

Z _ci，m＝0.5(Z _ci+Z _cm) (21)

P _ci，m＝RT _ci，mZ _ci，m/V _ci，m (22)

Ω _ai，m＝0.5(Ω _ai+Ω _am) (23)

To vapour phase:

$a^G=\underset{i}{\mathrm{\Σ}}\underset{m}{\mathrm{\Σ}}{y}_i{y}_m{a}_{i,m}---\left(24\right)$

$b^G=\underset{i}{\mathrm{\Σ}}{y}_i{b}_i---\left(25\right)$

Order

A ^G＝a ^GP/R ²T ² (26)

B ^G＝b ^GP/RT (27)

α ^G＝2B ^G-1 (28)

${\mathrm{\β}}^G=A^G-3B^G-5B^{G^2}---\left(29\right)$

${\mathrm{\γ}}^G=2(B^{G^3}+B^{G^2})-A^G{B}^G---\left(30\right)$

Getting initial value is 1-0.6P _r, separate following equation with Newton method, promptly obtain vapour phase compressibility factor Z ^G

$Z^{G^3}+{\mathrm{\α}}^G{Z}^{G^2}+\mathrm{\β}{Z}^G+{\mathrm{\γ}}^G=0---\left(31\right)$

Then,

v ^G＝RT/PZ ^G (32)

${\mathrm{\ξ}}^G=0.242536\ln \frac{v^G+3.561553b^G}{v^G-0.561553b^G}---\left(33\right)$

To liquid phase:

$a^L=\underset{i}{\mathrm{\Σ}}\underset{m}{\mathrm{\Σ}}{x}_i{x}_m{a}_{i,m}---\left(34\right)$

$b^L=\underset{i}{\mathrm{\Σ}}{x}_i{b}_i---\left(35\right)$

Order

A ^L＝a ^LP/R ²T ² (36)

B ^L＝b ^LP/RT (37)

α ^L＝2B ^L-1 (38)

${\mathrm{\β}}^L=A^L-3B^L-5B^{L^2}---\left(39\right)$

${\mathrm{\γ}}^L=2(B^{L^3}+B^{L^2})-A^L{B}^L---\left(40\right)$

Getting initial value is P _r(0.106+0.078P _r), separate following equation with Newton method, promptly obtain liquid phase compressibility factor Z ^L

$Z^{L^3}+{\mathrm{\α}}^L{Z}^{L^2}+\mathrm{\β}{Z}^L+{\mathrm{\γ}}^L=0---\left(41\right)$

Then,

v ^L＝RT/PZ ^L (42)

${\mathrm{\ξ}}^L=0.242536\ln \frac{v^L+3.561553b^L}{v^L-0.561553b^L}---\left(43\right)$

Ω _ai＝C _i-D _iτ+E _iτ ²-W _iτ ³ (44)

Ω _b＝0.070721 (45)

τ＝0.01T (46)

Wherein, A, B, α, β, γ, τ are intermediate variables, and C, D, E, W are constants, T _c, P _c, V _c, Z _cBe respectively critical-temperature, pressure, volume and compressibility factor, P _rBe reduced pressure, R is a gas constant, k _{I, m}The mutual coefficient of binary of representing i component and m component, k _{I, m}Be constant, subscript c represents the character of critical point, and subscript r represents reduced state, subscript i, and m represents the binary mixture of i component and m component, Ω _a, Ω _bIt is intermediate variable.

Further, in described step 7), host computer is passed to control station with result of calculation and is shown, and by fieldbus result of calculation is delivered to operator station and shows.

Beneficial effect of the present invention mainly shows: internal thermally coupled air separation column is carried out the dynamic flow simulation, and computational speed is fast, and analog result is accurate, and can be used for guidance production and further optimize, control and study, thus the raising productivity effect.

Description of drawings

Fig. 1 is the hardware structure diagram of dynamic flow simulation system proposed by the invention.

Fig. 2 is an internal thermally coupled air separation column structural representation of the present invention.

Fig. 3 is the functional block diagram of host computer of the present invention.

The specific embodiment

Below in conjunction with accompanying drawing the present invention is further described.

Embodiment 1

With reference to Fig. 1, Fig. 2, Fig. 3, a kind of dynamic flow of internal thermally coupled air separation column simulation system, comprise field intelligent instrument 2, data-interface 3, control station 4, database 5 and host computer 6 that internal thermally coupled air separation column 1 connects, intelligence instrument 2 connects fieldbus, described fieldbus connects data-interface 3, described data-interface 3 is connected with control station 4, database 5 and host computer 6, and described host computer 6 comprises:

Signal acquisition module 7 is in order to gather current production status data;

Find the solution and calculate primary module 8,, adopt following process to finish in order to find the solution calculating:

1) structural parameters and the operating parameter of setting tower are set initial moment tstart, stop tend constantly;

2) specify each column plate liquid phase of initial time to form and the liquid phase flow, make current iteration time t=tstart;

3), calculate its equilibrium temperature by the bubble point method respectively and vapour phase is formed to each column plate;

4), calculate the enthalpy of vapour-liquid phase respectively to each column plate;

5) calculate the vapour-liquid phase flow rate of each column plate by formula (1) (2):

$V_{j+1}{H}_{j+1}^G+U_{j-1}{H}_{j-1}^L+F_j{H}_j^F-(V_j+S_j^G)H_j^G-(U_j+S_j^L)H_j^L-Q_j=0---\left(1\right)$

$V_{j+1}+U_{j-1}+F_j^G+F_j^L-(V_j+S_j^G)-(U_j+S_j^L)=0---\left(2\right)$

Wherein, V represents the vapour phase flow, and U represents the liquid phase flow, and F represents feed rate, H ^FExpression charging enthalpy, S represent that side carries flow, H ^GAnd H ^LBe respectively vapour-liquid phase enthalpy, subscript j-1, j, j+1 represent j-1, j, j+1 piece plate respectively, and subscript L represents liquid phase, and subscript G represents vapour phase, and Q represents the thermal coupling amount, is calculated by following formula:

Q＝UAΔT (3)

Wherein, UA represents heat couple coefficient, the Δ T temperature difference between column plate of representing to be coupled;

6) calculating (t+ Δ t) each column plate liquid phase constantly forms and the liquid phase flow:

$\frac{{\mathrm{dx}}_{i,j}}{\mathrm{dt}}{M}_j=U_{j-1}{x}_{i,j-1}+V_{j+1}{y}_{i,j+1}+F_j{z}_{i,j}-(U_j+S_j^L)x_{i,j}-(V_j+S_j^G)y_{i,j}----\left(4\right)$

$[U_{j-1}+V_{j+1}+F_j-(U_j+S_j^L)-(V_j+S_j^G)]x_{i,j}$

$\frac{d{U}_j}{\mathrm{dt}}=518.13{\mathrm{\ρ}}_j^{-1/3}{A}_a^{-1}{l}^{2/3}{U}_j^{1/3}[U_{j-1}+V_{j+1}+F_j-(U_j+S_j^L)-(V_j+S_j^G)]---\left(5\right)$

$x_{i,j}^{(t+\mathrm{\Δt})}=x_{i,j}^{\left(t\right)}+\frac{d{x}_{i,j}}{\mathrm{dt}}\mathrm{\Δt}---\left(6\right)$

$U_j^{(t+\mathrm{\Δt})}=U_j^{\left(t\right)}+\frac{d{U}_j}{\mathrm{dt}}\mathrm{\Δt}---\left(7\right)$

Wherein, Δ t is an iteration step length, and x is that liquid phase is formed, and y is that vapour phase is formed, and z is a feed composition, and ρ is a density of liquid phase, A _aBe the column plate effective area, l is that the column plate weir is long, subscript i=1,2,3 expression components, and corresponding successively nitrogen, argon, oxygen, subscript (t) and (t+ Δ t) are represented t and t+ Δ t constantly respectively, M represents the column plate liquid holdup, is calculated by following formula:

$M_j={\mathrm{\ρ}}_j{A}_{\mathrm{aj}}[h_{\mathrm{wj}}+0.00284{\left(\frac{U_j}{{\mathrm{\ρ}}_j{l}_j}\right)}^{2/3}]---\left(8\right)$

H wherein _wIt is the overflow height of weir;

7) make t=t+ Δ t, return 3 with new each column plate liquid phase composition and liquid phase flow) iteration, up to t 〉=tend, finishing iteration, output result.

Described host computer 6 also comprises: bubble point method module 9, and in order to calculate its equilibrium temperature by the bubble point method and vapour phase is formed, its process is as follows:

3.1) supposition column plate equilibrium temperature;

3.2) calculate the VLE constant, adopt following process to finish:

$\ln {\mathrm{\Φ}}_i^L=\ln \frac{\mathrm{RT}}{P(v^L-b^L)}-\frac{b_i}{b^L}(1-Z^L)+{\mathrm{\ξ}}^L{a}^L(\frac{b_i}{b^L}-\frac{2\underset{m}{\mathrm{\Σ}}{x}_m{a}_{i,m}}{a^L})/b^L\mathrm{RT}---\left(9\right)$

$\ln {\mathrm{\Φ}}_i^G=\ln \frac{\mathrm{RT}}{P(v^G-b^G)}-\frac{b_i}{b^G}(1-Z^G)+{\mathrm{\ξ}}^G{a}^G(\frac{b_i}{b^G}-\frac{2\underset{m}{\mathrm{\Σ}}{x}_m{a}_{i,m}}{a^G})/b^G\mathrm{RT}---\left(10\right)$

$K_i={\mathrm{\Φ}}_i^L/{\mathrm{\Φ}}_i^G---\left(11\right)$

y _i＝K _ix _i (12)

Wherein, Φ represents fugacity coefficient, and subscript L represents liquid phase, and subscript G represents vapour phase, and R is a gas constant, and T is a temperature, and P is a column plate pressure, subscript m=1,2,3 expression components, corresponding successively nitrogen, argon, oxygen, molal volume v, physical parameter b ^G, b ^L, b _i, a ^G, a ^L, a _{I, m}, ξ ^G, ξ ^L, vapour phase compressibility factor Z ^G, liquid phase compressibility factor Z ^LCalculate by the rerum natura module;

3.3) check

Whether set up, set up then finishing iteration, return result of calculation, otherwise, upgrade the column plate equilibrium temperature, return 3.2) the continuation iteration;

Described host computer 6 also comprises: enthalpy module 10, and in order to calculate vapour-liquid phase enthalpy of mixing, its process is as follows:

$H_i^{*}=c_i+d_{i}T+e_i{T}^2+f_i{T}^3+h_i{T}^4---\left(13\right)$

$H^{*}=\underset{i}{\mathrm{\Σ}}{y}_i{H}_i^{*}---\left(14\right)$

$H^G=H^{*}-\mathrm{RT}(1-Z^G)-{\mathrm{\ξ}}^G(a^G-T\frac{d{a}^G}{\mathrm{dT}})/b^G---\left(15\right)$

$H^L=H^{*}-\mathrm{RT}(1-Z^L)-{\mathrm{\ξ}}^L(a^L-T\frac{d{a}^L}{\mathrm{dT}})/b^L---\left(16\right)$

H wherein _i ^*The enthalpy of representing i pure component perfect gas, H ^*Be mixture perfect gas enthalpy, c, d, e, f, h are constant.

Described host computer 6 also comprises: rerum natura module 11, and in order to calculate physical parameter, its process is as follows:

$a_{i,m}={\mathrm{\Ω}}_{\mathrm{ai},m}{R}^2{T}_{\mathrm{ci},m}^2/P_{\mathrm{ci},m}---\left(17\right)$

b _i＝Ω _bRT _ci/P _cia (18)

$T_{\mathrm{ci},m}=\sqrt{T_{\mathrm{ci}}{T}_{\mathrm{cj}}}(1-k_{i,m})---\left(19\right)$

$V_{\mathrm{ci},m}=0.125{(V_{\mathrm{ci},m}^{1/3}+V_{\mathrm{ci},m}^{1/3})}^3---\left(20\right)$

Z _ci，m＝0.5(Z _ci+Z _cm) (21)

P _ci，m＝RT _ci，mZ _ci，m/V _ci，m (22)

Ω _ai，m＝0.5(Ω _ai+Ω _am) (23)

To vapour phase:

$a^G=\underset{i}{\mathrm{\Σ}}\underset{m}{\mathrm{\Σ}}{y}_i{y}_m{a}_{i,m}---\left(24\right)$

$b^G=\underset{i}{\mathrm{\Σ}}{y}_i{b}_i---\left(25\right)$

Order

A ^G＝a ^GP/R ²T ² (26)

B ^G＝b ^GP/RT (27)

α ^G＝2B ^G-1 (28)

${\mathrm{\β}}^G=A^G-3B^G-5B^{G^2}---\left(29\right)$

${\mathrm{\γ}}^G=2(B^{G^3}+B^{G^2})-A^G{B}^G---\left(30\right)$

Getting initial value is 1-0.6P _r, separate following equation with Newton method, promptly obtain vapour phase compressibility factor Z ^G

$Z^{G^3}+{\mathrm{\α}}^G{Z}^{G^2}+\mathrm{\β}{Z}^G+{\mathrm{\γ}}^G=0---\left(31\right)$

Then,

v ^G＝RT/PZ ^G (32)

${\mathrm{\ξ}}^G=0.242536\ln \frac{v^G+3.561553b^G}{v^G-0.561553b^G}---\left(33\right)$

To liquid phase:

$a^L=\underset{i}{\mathrm{\Σ}}\underset{m}{\mathrm{\Σ}}{x}_i{x}_m{a}_{i,m}---\left(34\right)$

$b^L=\underset{i}{\mathrm{\Σ}}{x}_i{b}_i---\left(35\right)$

Order

A ^L＝a ^LP/R ²T ² (36)

B ^L＝b ^LP/RT (37)

α ^L＝2B ^L-1 (38)

${\mathrm{\β}}^L=A^L-3B^L-5B^{L^2}---\left(39\right)$

${\mathrm{\γ}}^L=2(B^{L^3}+B^{L^2})-A^L{B}^L---\left(40\right)$

Getting initial value is P _r(0.106+0.078P _r), separate following equation with Newton method, promptly obtain liquid phase compressibility factor Z ^L

$Z^{L^3}+{\mathrm{\α}}^L{Z}^{L^2}+\mathrm{\β}{Z}^L+{\mathrm{\γ}}^L=0---\left(41\right)$

Then,

v ^L＝RT/PZ ^L (42)

${\mathrm{\ξ}}^L=0.242536\ln \frac{v^L+3.561553b^L}{v^L-0.561553b^L}---\left(43\right)$

Ω _ai＝C _i-D _iτ+E _iτ ²-W _iτ ³ (44)

Ω _b＝0.070721 (45)

τ＝0.01T (46)

Wherein, A, B, α, β, γ, τ are intermediate variables, and C, D, E, W are constants, T _c, P _c, V _c, Z _cBe respectively critical-temperature, pressure, volume and compressibility factor, P _rBe reduced pressure, R is a gas constant, k _{I, m}The mutual coefficient of binary of representing i component and m component, k _{I, m}Be constant, subscript c represents the character of critical point, and subscript r represents reduced state, subscript i, and m represents the binary mixture of i component and m component, Ω _a, Ω _bIt is intermediate variable.

Described host computer 6 also comprises: display module 12 as a result, are used for that result of calculation is passed to control station and show, and by fieldbus result of calculation is delivered to operator station and shows.

The hardware structure diagram of the dynamic flow of internal thermally coupled air separation column simulation system of present embodiment as shown in Figure 1, described process simulation system core by comprise signal acquisition module 7, find the solution calculate primary module 8, bubble point method module 9, enthalpy module 10, rerum natura module 11, the host computer 6 of display module 12 and man-machine interface constitutes as a result, comprise in addition: field intelligent instrument 2, data-interface 3, control station 4, database 5 and fieldbus.Internal thermally coupled air separation column 1, intelligence instrument 2, data-interface 3, control station 4, database 5, host computer 6 link to each other successively by fieldbus, realize uploading and assigning of information flow.Process simulation system moves on host computer 6, can carry out information exchange with first floor system easily.

The functional block diagram of the optimization system of present embodiment mainly comprises signal acquisition module 7, finds the solution and calculate primary module 8, bubble point method module 9, enthalpy module 10, rerum natura module 11, display module 12 etc. as a result as shown in Figure 3.

Described dynamic flow analogy method is implemented according to following steps:

1) structural parameters of setting tower are gathered the production status data, set initial moment tstart, stop tend constantly:

2) specify each column plate liquid phase of initial time to form and the liquid phase flow, make current iteration time t=tstart;

3), calculate its equilibrium temperature by the bubble point method respectively and vapour phase is formed to each column plate;

4), calculate the enthalpy of its vapour-liquid phase to each column plate;

5) simultaneous formula (1) (2) is calculated the vapour-liquid phase flow rate of each column plate:

$V_{j+1}{H}_{j+1}^G+U_{j-1}{H}_{j-1}^L+F_j{H}_j^F-(V_j+S_j^G)H_j^G-(U_j+S_j^L)H_j^L-Q_j=0---\left(1\right)$

$V_{j+1}+U_{j-1}+F_j^G+F_j^L-(V_j+S_j^G)-(U_j+S_j^L)=0---\left(2\right)$

Wherein, V represents the vapour phase flow, and U represents the liquid phase flow, and F represents feed rate, H ^FExpression charging enthalpy, S represent that side carries flow, H ^GAnd H ^LBe respectively vapour-liquid phase enthalpy, subscript j-1, j, j+1 represent j-1, j, j+1 piece plate respectively, and subscript L represents liquid phase, and subscript G represents vapour phase, and Q represents the thermal coupling amount, is calculated by following formula:

Q＝UAΔT (3)；

Wherein, UA represents heat couple coefficient, the Δ T temperature difference between column plate of representing to be coupled;

6) calculating (t+ Δ t) each column plate liquid phase constantly forms and the liquid phase flow:

$\frac{{\mathrm{dx}}_{i,j}}{\mathrm{dt}}{M}_j=U_{j-1}{x}_{i,j-1}+V_{j+1}{y}_{i,j+1}+F_j{z}_{i,j}-(U_j+S_j^L)x_{i,j}-(V_j+S_j^G)y_{i,j}----\left(4\right)$

$[U_{j-1}+V_{j+1}+F_j-(U_j+S_j^L)-(V_j+S_j^G)]x_{i,j}$

$\frac{d{U}_j}{\mathrm{dt}}=518.13{\mathrm{\ρ}}_j^{-1/3}{A}_a^{-1}{l}^{2/3}{U}_j^{1/3}[U_{j-1}+V_{j+1}+F_j-(U_j+S_j^L)-(V_j+S_j^G)]---\left(5\right)$

$x_{i,j}^{(t+\mathrm{\Δt})}=x_{i,j}^{\left(t\right)}+\frac{d{x}_{i,j}}{\mathrm{dt}}\mathrm{\Δt}---\left(6\right)$

$U_j^{(t+\mathrm{\Δt})}=U_j^{\left(t\right)}+\frac{d{U}_j}{\mathrm{dt}}\mathrm{\Δt}---\left(7\right)$

Wherein, Δ t is an iteration step length, and x is that liquid phase is formed, and y is that vapour phase is formed, and z is a feed composition, and ρ is a density of liquid phase, A _aBe the column plate effective area, l is that the column plate weir is long, subscript i=1,2,3 expression components, and corresponding successively nitrogen, argon, oxygen, subscript (t) and (t+ Δ t) are represented t and t+ Δ t constantly respectively, M represents the column plate liquid holdup, is calculated by following formula:

$M_j={\mathrm{\ρ}}_j{A}_{\mathrm{aj}}[h_{\mathrm{wj}}+0.00284{\left(\frac{U_j}{{\mathrm{\ρ}}_j{l}_j}\right)}^{2/3}]---\left(8\right)$

H wherein _wIt is the overflow height of weir;

7) make t=t+ Δ t, return 3 with new each column plate liquid phase composition and liquid phase flow) iteration, up to t 〉=tend, finishing iteration, output result.

Embodiment 2

With reference to Fig. 1, Fig. 2, Fig. 3, a kind of dynamic flow of internal thermally coupled air separation column analogy method, described flowsheeting method may further comprise the steps:

1) structural parameters of setting tower are gathered the production status data, set initial moment tstart, stop tend constantly;

2) specify each column plate liquid phase of initial time to form and the liquid phase flow, make current iteration time t=tstart;

3), calculate its equilibrium temperature by the bubble point method respectively and vapour phase is formed to each column plate;

4), calculate the enthalpy of its vapour-liquid phase to each column plate;

5) simultaneous formula (1) (2) is calculated the vapour-liquid phase flow rate of each column plate:

$V_{j+1}{H}_{j+1}^G+U_{j-1}{H}_{j-1}^L+F_j{H}_j^F-(V_j+S_j^G)H_j^G-(U_j+S_j^L)H_j^L-Q_j=0---\left(1\right)$

$V_{j+1}+U_{j-1}+F_j^G+F_j^L-(V_j+S_j^G)-(U_j+S_j^L)=0---\left(2\right)$

Wherein, V represents the vapour phase flow, and U represents the liquid phase flow, and F represents feed rate, H ^FExpression charging enthalpy, S represent that side carries flow, H ^GAnd H ^LBe respectively vapour-liquid phase enthalpy, subscript j-1, j, j+1 represent j-1, j, j+1 piece plate respectively, and subscript L represents liquid phase, and subscript G represents vapour phase, and Q represents the thermal coupling amount, is calculated by following formula:

Q＝UAΔT (3)；

Wherein, UA represents heat couple coefficient, the Δ T temperature difference between column plate of representing to be coupled;

6) calculating (t+ Δ t) each column plate liquid phase constantly forms and the liquid phase flow:

$\frac{{\mathrm{dx}}_{i,j}}{\mathrm{dt}}{M}_j=U_{j-1}{x}_{i,j-1}+V_{j+1}{y}_{i,j+1}+F_j{z}_{i,j}-(U_j+S_j^L)x_{i,j}-(V_j+S_j^G)y_{i,j}----\left(4\right)$

$[U_{j-1}+V_{j+1}+F_j-(U_j+S_j^L)-(V_j+S_j^G)]x_{i,j}$

$\frac{d{U}_j}{\mathrm{dt}}=518.13{\mathrm{\ρ}}_j^{-1/3}{A}_a^{-1}{l}^{2/3}{U}_j^{1/3}[U_{j-1}+V_{j+1}+F_j-(U_j+S_j^L)-(V_j+S_j^G)]---\left(5\right)$

$x_{i,j}^{(t+\mathrm{\Δt})}=x_{i,j}^{\left(t\right)}+\frac{d{x}_{i,j}}{\mathrm{dt}}\mathrm{\Δt}---\left(6\right)$

$U_j^{(t+\mathrm{\Δt})}=U_j^{\left(t\right)}+\frac{d{U}_j}{\mathrm{dt}}\mathrm{\Δt}---\left(7\right)$

Wherein, Δ t is an iteration step length, and x is that liquid phase is formed, and y is that vapour phase is formed, and z is a feed composition, and ρ is a density of liquid phase, A _aBe the column plate effective area, l is that the column plate weir is long, subscript i=1,2,3 expression components, and corresponding successively nitrogen, argon, oxygen, subscript (t) and (t+ Δ t) are represented t and t+ Δ t constantly respectively, M represents the column plate liquid holdup, is calculated by following formula:

$M_j={\mathrm{\ρ}}_j{A}_{\mathrm{aj}}[h_{\mathrm{wj}}+0.00284{\left(\frac{U_j}{{\mathrm{\ρ}}_j{l}_j}\right)}^{2/3}]---\left(8\right)$

H wherein _wIt is the overflow height of weir;

7) make t=t+ Δ t, return the step 3) iteration with new each column plate liquid phase composition and liquid phase flow, up to t 〉=tend, finishing iteration, output result.

In the described step 3), the bubble point method is calculated its equilibrium temperature and vapour phase, adopts following process to finish:

3.1) supposition column plate equilibrium temperature;

3.2) calculate the VLE constant, adopt following process to finish:

$\ln {\mathrm{\Φ}}_i^L=\ln \frac{\mathrm{RT}}{P(v^L-b^L)}-\frac{b_i}{b^L}(1-Z^L)+{\mathrm{\ξ}}^L{a}^L(\frac{b_i}{b^L}-\frac{2\underset{m}{\mathrm{\Σ}}{x}_m{a}_{i,m}}{a^L})/b^L\mathrm{RT}---\left(9\right)$

$\ln {\mathrm{\Φ}}_i^G=\ln \frac{\mathrm{RT}}{P(v^G-b^G)}-\frac{b_i}{b^G}(1-Z^G)+{\mathrm{\ξ}}^G{a}^G(\frac{b_i}{b^G}-\frac{2\underset{m}{\mathrm{\Σ}}{x}_m{a}_{i,m}}{a^G})/b^G\mathrm{RT}---\left(10\right)$

$K_i={\mathrm{\Φ}}_i^L/{\mathrm{\Φ}}_i^G---\left(11\right)$

y _i＝K _ix _i (12)

Wherein, Φ represents fugacity coefficient, and subscript L represents liquid phase, and subscript G represents vapour phase, and R is a gas constant, and T is a temperature, and P is a column plate pressure, subscript m=1,2,3 expression components, corresponding successively nitrogen, argon, oxygen, molal volume v, parameter b ^G, b ^L, b _i, a ^G, a ^L, a _{I, m}, ξ ^G, ξ ^L, vapour phase compressibility factor Z ^G, liquid phase compressibility factor Z ^LCalculate by the physical parameter computational methods;

3.3) check

Whether set up, set up then finishing iteration, return result of calculation, otherwise, upgrade the column plate equilibrium temperature, return 3.2) the continuation iteration.

In the described step 4), described enthalpy computational methods process is as follows:

$H_i^{*}=c_i+d_{i}T+e_i{T}^2+f_i{T}^3+h_i{T}^4---\left(13\right)$

$H^{*}=\underset{i}{\mathrm{\Σ}}{y}_i{H}_i^{*}---\left(14\right)$

$H^G=H^{*}-\mathrm{RT}(1-Z^G)-{\mathrm{\ξ}}^G(a^G-T\frac{d{a}^G}{\mathrm{dT}})/b^G---\left(15\right)$

$H^L=H^{*}-\mathrm{RT}(1-Z^L)-{\mathrm{\ξ}}^L(a^L-T\frac{d{a}^L}{\mathrm{dT}})/b^L---\left(16\right)$

H wherein _i ^*The enthalpy of representing i pure component perfect gas, H ^*Be mixture perfect gas enthalpy, c, d, e, f, h are constant.

The process of described physical parameter computational methods is as follows:

$a_{i,m}={\mathrm{\Ω}}_{\mathrm{ai},m}{R}^2{T}_{\mathrm{ci},m}^2/P_{\mathrm{ci},m}---\left(17\right)$

b _i＝Ω _bRT _ci/P _cia (18)

$T_{\mathrm{ci},m}=\sqrt{T_{\mathrm{ci}}{T}_{\mathrm{cj}}}(1-k_{i,m})---\left(19\right)$

$V_{\mathrm{ci},m}=0.125{(V_{\mathrm{ci},m}^{1/3}+V_{\mathrm{ci},m}^{1/3})}^3---\left(20\right)$

Z _ci，m＝0.5(Z _ci+Z _cm) (21)

P _ci，m＝RT _ci，mZ _ci，m/V _ci，m (22)

Ω _ai，m＝0.5(Ω _ai+Ω _am) (23)

To vapour phase:

$a^G=\underset{i}{\mathrm{\Σ}}\underset{m}{\mathrm{\Σ}}{y}_i{y}_m{a}_{i,m}---\left(24\right)$

$b^G=\underset{i}{\mathrm{\Σ}}{y}_i{b}_i---\left(25\right)$

Order

A ^G＝a ^GP/R ²T ² (26)

B ^G＝b ^GP/RT (27)

α ^G＝2B ^G-1 (28)

${\mathrm{\β}}^G=A^G-3B^G-5B^{G^2}---\left(29\right)$

${\mathrm{\γ}}^G=2(B^{G^3}+B^{G^2})-A^G{B}^G---\left(30\right)$

Getting initial value is 1-0.6P _r, separate following equation with Newton method, promptly obtain vapour phase compressibility factor Z ^G

$Z^{G^3}+{\mathrm{\α}}^G{Z}^{G^2}+\mathrm{\β}{Z}^G+{\mathrm{\γ}}^G=0---\left(31\right)$

Then,

v ^G＝RT/PZ ^G (32)

${\mathrm{\ξ}}^G=0.242536\ln \frac{v^G+3.561553b^G}{v^G-0.561553b^G}---\left(33\right)$

To liquid phase:

$a^L=\underset{i}{\mathrm{\Σ}}\underset{m}{\mathrm{\Σ}}{x}_i{x}_m{a}_{i,m}---\left(34\right)$

$b^L=\underset{i}{\mathrm{\Σ}}{x}_i{b}_i---\left(35\right)$

Order

A ^L＝a ^LP/R ²T ² (36)

B ^L＝b ^LP/RT (37)

α ^L＝2B ^L-1 (38)

${\mathrm{\β}}^L=A^L-3B^L-5B^{L^2}---\left(39\right)$

${\mathrm{\γ}}^L=2(B^{L^3}+B^{L^2})-A^L{B}^L---\left(40\right)$

Getting initial value is P _r(0.106+0.078P _r), separate following equation with Newton method, promptly obtain liquid phase compressibility factor Z ^L

$Z^{L^3}+{\mathrm{\α}}^L{Z}^{L^2}+\mathrm{\β}{Z}^L+{\mathrm{\γ}}^L=0---\left(41\right)$

Then,

v ^L＝RT/PZ ^L (42)

${\mathrm{\ξ}}^L=0.242536\ln \frac{v^L+3.561553b^L}{v^L-0.561553b^L}---\left(43\right)$

Ω _ai＝C _i-D _iτ+E _iτ ²-W _iτ ³ (44)

Ω _b＝0.070721 (45)

τ＝0.01T (46)

Wherein, A, B, α, β, γ, τ are intermediate variables, and C, D, E, W are constants, T _c, P _c, V _c, Z _cBe respectively critical-temperature, pressure, volume and compressibility factor, P _rBe reduced pressure, R is a gas constant, k _{I, m}The mutual coefficient of binary of representing i component and m component, k _{I, m}Be constant, subscript c represents the character of critical point, and subscript r represents reduced state, subscript i, and m represents the binary mixture of i component and m component, Ω _a, Ω _bIt is intermediate variable.

In described step 7), host computer is passed to control station with result of calculation and is shown, and by fieldbus result of calculation is delivered to operator station and shows.

Dynamic flow of internal thermally coupled air separation column simulation system and method proposed by the invention, be described by above-mentioned concrete implementation step, person skilled obviously can be in not breaking away from content of the present invention, spirit and scope to device as herein described with method of operating is changed or suitably change and combination, realize the technology of the present invention.Special needs to be pointed out is, the replacement that all are similar and change apparent to one skilled in the artly, they all can be regarded as being included in spirit of the present invention, scope and the content.

Claims (5)

1. dynamic flow of internal thermally coupled air separation column simulation system, comprise the field intelligent instrument and control station, database and the host computer that are connected with internal thermally coupled air separation column, intelligence instrument is connected with control station, database and host computer, it is characterized in that: described host computer comprises:

Signal acquisition module is in order to gather current production status data;

Find the solution the calculating primary module,, adopt following process to finish in order to find the solution calculating:

1) structural parameters and the operating parameter of setting tower are set initial moment tstart, stop tend constantly;

2) specify each column plate liquid phase of initial time to form and the liquid phase flow, make current iteration time t=tstart;

3), calculate its equilibrium temperature by the bubble point method respectively and vapour phase is formed to each column plate;

4), calculate the enthalpy of vapour-liquid phase respectively to each column plate;

5) calculate the vapour-liquid phase flow rate of each column plate by formula (1) (2):

$V_{j+1}{H}_{j+1}^G+U_{j-1}{H}_{j-1}^L+F_j{H}_j^F-(V_j+S_j^G)H_j^G-(U_j+S_j^L)H_j^L-Q_j=0---\left(1\right)$

$V_{j+1}+U_{j-1}+F_j^G+F_j^L-(V_j+S_j^G)-(U_j+S_j^L)=0---\left(2\right)$

Wherein, V represents the vapour phase flow, and U represents the liquid phase flow, and F represents feed rate, H ^FExpression charging enthalpy, S represent that side carries flow, H ^GAnd H ^LBe respectively vapour-liquid phase enthalpy, subscript j-1, j, j+1 represent j-1, j, j+1 piece plate respectively, and subscript L represents liquid phase, and subscript G represents vapour phase, and Q represents the thermal coupling amount, is calculated by following formula:

Q＝UAΔT (3)

Wherein, UA represents heat couple coefficient, the Δ T temperature difference between column plate of representing to be coupled;

6) calculating (t+ Δ t) each column plate liquid phase constantly forms and the liquid phase flow:

$\frac{{\mathrm{dx}}_{i,j}}{\mathrm{dt}}{M}_j=U_{j-1}{x}_{i,j-1}+V_{j+1}{y}_{i,j+1}+F_j{z}_{i,j}-(U_j+S_j^L)x_{i,j}-(V_j+S_j^G)y_{i,j}----\left(4\right)$

$[U_{j-1}+V_{j+1}+F_j-(U_j+S_j^L)-(V_j+S_j^G)]x_{i,j}$

$\frac{{\mathrm{dU}}_j}{\mathrm{dt}}=518.13{\mathrm{\ρ}}_j^{-1/3}{A}_a^{-1}{l}^{2/3}{U}_j^{1/3}[U_{j-1}+V_{j+1}+F_j-(U_j+S_j^L)-(V_j+S_j^G)]---\left(5\right)$

$x_{i,j}^{(t+\mathrm{\Δt})}=x_{i,j}^{\left(t\right)}+\frac{{\mathrm{dx}}_{i,j}}{\mathrm{dt}}\mathrm{\Δt}---\left(6\right)$

$U_j^{(t+\mathrm{\Δt})}=U_j^{\left(t\right)}+\frac{d{U}_j}{\mathrm{dt}}\mathrm{\Δt}---\left(7\right)$

Wherein, Δ t is an iteration step length, and x is that liquid phase is formed, and y is that vapour phase is formed, and z is a feed composition, and ρ is a density of liquid phase, A _aBe the column plate effective area, l is that the column plate weir is long, subscript i=1,2,3 expression components, and corresponding successively nitrogen, argon, oxygen, subscript (t) and (t+ Δ t) are represented t and t+ Δ t constantly respectively, M represents the column plate liquid holdup, is calculated by following formula:

$M_j={\mathrm{\ρ}}_j{A}_{\mathrm{aj}}[h_{\mathrm{wj}}+0.00284{\left(\frac{U_j}{{\mathrm{\ρ}}_j{l}_j}\right)}^{2/3}]---\left(8\right)$

H wherein _wIt is the overflow height of weir;

7) make t=t+ Δ t, return 3 with new each column plate liquid phase composition and liquid phase flow) iteration, up to t 〉=tend, finishing iteration, output result.

2. dynamic flow of internal thermally coupled air separation column simulation system as claimed in claim 1 is characterized in that: described host computer also comprises: bubble point method module, and in order to calculate its equilibrium temperature by the bubble point method and vapour phase is formed, its process is as follows:

3.1) supposition column plate equilibrium temperature;

3.2) calculate the VLE constant, adopt following process to finish:

$\ln {\mathrm{\Φ}}_i^L=\ln \frac{\mathrm{RT}}{P(v^L-b^L)}-\frac{b_i}{b^L}(1-Z^L)+{\mathrm{\ξ}}^L{a}^L(\frac{b_i}{b^L}-\frac{2\underset{m}{\mathrm{\Σ}}{x}_m{a}_{i,m}}{a^L})/b^L\mathrm{RT}---\left(9\right)$

$\ln {\mathrm{\Φ}}_i^G=\ln \frac{\mathrm{RT}}{P(v^G-b^G)}-\frac{b_i}{b^G}(1-Z^G)+{\mathrm{\ξ}}^G{a}^G(\frac{b_i}{b^G}-\frac{2\underset{m}{\mathrm{\Σ}}{x}_m{a}_{i,m}}{a^G})/b^G\mathrm{RT}---\left(10\right)$

$K_i={\mathrm{\Φ}}_i^L/{\mathrm{\Φ}}_i^G---\left(11\right)$

y _i＝K _ix _i (12)

Wherein, Φ represents fugacity coefficient, and subscript L represents liquid phase, and subscript G represents vapour phase, and R is a gas constant, and T is a temperature, and P is a column plate pressure, subscript m=1,2,3 expression components, corresponding successively nitrogen, argon, oxygen, molal volume v, parameter b ^G, b ^L, b _i, a ^G, a ^L, q _{I, m}, ξ ^G, ξ ^L, vapour phase compressibility factor Z ^G, liquid phase compressibility factor Z ^LCalculate by the rerum natura module;

3.3) check

Whether set up, set up then finishing iteration, return result of calculation, otherwise, upgrade the column plate equilibrium temperature, return 3.2) the continuation iteration.

3. dynamic flow of internal thermally coupled air separation column simulation system as claimed in claim 1 or 2 is characterized in that: described host computer also comprises: the enthalpy module, and in order to calculate vapour-liquid phase enthalpy of mixing, its process is as follows:

$H_i^{*}=c_i+d_{i}T+e_i{T}^2+f_i{T}^3+h_i{T}^4---\left(13\right)$

$H^{*}=\underset{i}{\mathrm{\Σ}}{y}_i{H}_i^{*}---\left(14\right)$

$H^G=H^{*}-\mathrm{RT}(1-Z^G)-{\mathrm{\ξ}}^G(a^G-T\frac{d{a}^G}{\mathrm{dt}})/b^G---\left(15\right)$

$H^L=H^{*}-\mathrm{RT}(1-Z^L)-{\mathrm{\ξ}}^L(a^L-T\frac{d{a}^L}{\mathrm{dt}})/b^L---\left(16\right)$

Wherein

The enthalpy of representing i pure component perfect gas, H ^*Be mixture perfect gas enthalpy, c, d, e, f, h are constant.

4. dynamic flow of internal thermally coupled air separation column simulation system as claimed in claim 2 is characterized in that: described host computer also comprises: the rerum natura module, and in order to calculate physical parameter, its process is as follows:

$a_{i,m}={\mathrm{\Ω}}_{\mathrm{ai},m}{R}^2{T}_{\mathrm{ci},m}^2/P_{\mathrm{ci},m}---\left(17\right)$

b _i＝Ω _bRT _ci/P _cia (18)

$T_{\mathrm{ci},m}=\sqrt{T_{\mathrm{ci}}{T}_{\mathrm{cj}}}(1-k_{i,m})---\left(19\right)$

$V_{\mathrm{ci},m}=0.125{(V_{\mathrm{ci},m}^{1/3}+V_{\mathrm{ci},m}^{1/3})}^3---\left(20\right)$

Z _ci，m＝0.5(Z _ci+Z _cm) (21)

P _ci，m＝RT _ci，mZ _ci，m/V _ci，m (22)

Ω _ai，m＝0.5(Ω _ai+Ω _am) (23)

To vapour phase:

$a^G=\underset{i}{\mathrm{\Σ}}\underset{m}{\mathrm{\Σ}}{y}_i{y}_m{a}_{i,m}---\left(24\right)$

$b^G=\underset{i}{\mathrm{\Σ}}{y}_i{b}_i---\left(25\right)$

Order

A ^G＝a ^GP/R ²T ² (26)

B ^G＝b ^GP/RT (27)

α ^G＝2B ^G-1 (28)

${\mathrm{\β}}^G=A^G-3B^G-5B^{G^2}---\left(29\right)$

${\mathrm{\γ}}^G=2(B^{G^3}+B^{G^2})-A^G{B}^G---\left(30\right)$

Getting initial value is 1-0.6P _r, separate following equation with Newton method, promptly obtain vapour phase compressibility factor Z ^G

$Z^{G^3}+{\mathrm{\α}}^G{Z}^{G^2}+\mathrm{\β}{Z}^G+{\mathrm{\γ}}^G=0---\left(31\right)$

Then,

v ^G＝RT/PZ ^G (32)

${\mathrm{\ξ}}^G=0.242536\ln \frac{v^G+3.561553b^G}{v^G-0.561553b^G}---\left(33\right)$

To liquid phase:

$a^L=\underset{i}{\mathrm{\Σ}}\underset{m}{\mathrm{\Σ}}{x}_i{x}_m{a}_{i,m}---\left(34\right)$

$b^L=\underset{i}{\mathrm{\Σ}}{x}_i{b}_i---\left(35\right)$

Order

A ^L＝a ^LP/R ²T ² (36)

B ^L＝b ^LP/RT (37)

α ^L＝2B ^L-1 (38)

${\mathrm{\β}}^L=A^L-3B^L-5B^{L^2}---\left(39\right)$

${\mathrm{\γ}}^L=2(B^{L^3}+B^{L^2})-A^L{B}^L---\left(40\right)$

Getting initial value is P _r(0.106+0.078P _r), separate following equation with Newton method, promptly obtain liquid phase compressibility factor Z ^L

$Z^{L^3}+{\mathrm{\α}}^L{Z}^{L^2}+\mathrm{\β}{Z}^L+{\mathrm{\γ}}^L=0---\left(41\right)$

Then,

v ^L＝RT/PZ ^L (42)

${\mathrm{\ξ}}^L=0.242536\ln \frac{v^L+3.561553b^L}{v^L-0.561553b^L}---\left(43\right)$

Ω _ai＝C _i-D _iτ+E _iτ ²-W _iτ ³ (44)

Ω _b＝0.070721 (45)

τ＝0.01T (46)

Wherein, A, B, α, β, γ, carry intermediate variable, C, D, E, W are constants, T _c, P _c, V _c, Z _cBe respectively critical-temperature, pressure, volume and compressibility factor, P _rBe reduced pressure, R is a gas constant, k _{I, m}The mutual coefficient of binary of representing i component and m component, k _{I, m}Be constant, subscript c represents the character of critical point, and subscript r represents reduced state, subscript i, and m represents the binary mixture of i component and m component, Ω _a, Ω _bIt is intermediate variable.

5. dynamic flow of internal thermally coupled air separation column simulation system as claimed in claim 1 or 2 is characterized in that: described host computer also comprises:

Display module is used for that result of calculation is passed to control station and shows as a result, and by fieldbus result of calculation is delivered to operator station and shows.

Images