Gas-liquid flash evaporation calculation method, equipment, medium and product

Abstract

The invention relates to the technical field of chemical process simulation, in particular to a gas-liquid flash calculation method and device, a medium and a product. The method comprises the following steps: acquiring initial parameters of a multi-component system, performing composition estimation on the multi-component system based on a preset thermodynamic equation, determining a liquid phase fugacity coefficient and a gas phase fugacity coefficient of each component, performing lower limit constraint on the fugacity coefficients, and determining a phase equilibrium constant of each component based on the constrained fugacity coefficients. The method comprises the following steps: constructing a gas phase fraction solving equation through a phase equilibrium constant and a feeding molar fraction, and determining a target expression corresponding to the gas phase fraction solving equation based on a relative size between a current gas phase fraction and a preset threshold value, so as to iteratively carry out mixed solving on the target expression to update the gas phase fraction until the gas phase fraction meets a preset convergence condition, the gas-liquid flash calculation result representing that the current phase composition of the multi-component system meets the gas-liquid two-phase equilibrium condition is obtained, so that the efficiency and the accuracy of gas-liquid flash calculation are improved.

Description

Technical Field

[0001] This application relates to the field of chemical process simulation technology, and in particular to a gas-liquid flash evaporation calculation method, equipment, medium, and product. Background Technology

[0002] Gas-liquid flash evaporation calculation, as a fundamental computational unit in simulations of petrochemical, natural gas processing, coal chemical, and liquefied chemical processes, has the core task of solving the gas-liquid two-phase equilibrium state of multi-component systems based on thermodynamic equations, determining the gas phase fraction and the distribution of each component between the two phases. An efficient flash evaporation algorithm is crucial for ensuring the computational accuracy, convergence speed, and overall stability of process simulators, and is of great significance for the simulation and optimization of complex processes.

[0003] Currently, common calculations of gas-liquid flash evaporation mainly rely on second-order iterative methods such as Newton's method to solve for the phase equilibrium constant K and the Rachford-Rice equation (RR equation) to predict gas-liquid distribution. However, traditional methods suffer from serious numerical defects in practical engineering applications, especially under extreme conditions such as highly asymmetric component distribution (e.g., containing trace amounts of heavy or light components), near-critical conditions, or single-phase boundaries. First, because the thermodynamic equation of state may return extremely small fugacity coefficient values, the calculated phase equilibrium constant can span multiple orders of magnitude or even become ill-conditioned. This causes drastic changes or unbounded tendencies in the first and second derivatives of the Rachford-Rice equation, leading to instability in the iteration step size of Newton's method, which is prone to oscillations or divergence. Furthermore, the traditional direct introduction of higher-order iterative methods such as Halley's method in pursuit of faster convergence often backfires. Higher-order methods are more sensitive to derivatives, and under conditions of increased ill-conditioning in the aforementioned equations, the denominator of the iteration step size easily approaches zero, resulting in uncontrolled step size and numerical stability even worse than Newton's method. Furthermore, traditional techniques typically design independent solution processes for different flash evaporation modes such as constant pressure and temperature, and constant pressure and constant gas phase ratio. This results in fragmented code structures and a lack of a unified and robust interface for extending to more complex phase equilibrium systems of three phases and above, increasing the complexity of integrating, maintaining, and extending the algorithm in industrial software.

[0004] Therefore, there is an urgent need for a gas-liquid flash evaporation calculation method to improve the efficiency and accuracy of gas-liquid flash evaporation calculations. Summary of the Invention

[0005] This invention provides a gas-liquid flash evaporation calculation method, device, storage medium, and program product to improve the efficiency and accuracy of gas-liquid flash evaporation calculation.

[0006] In a first aspect, this application provides a method for calculating gas-liquid flash evaporation, including: Obtain initial parameters for the multi-component system, wherein the initial parameters include at least two of the following: pressure, temperature, gas phase fraction, and feed mole fraction of the multi-component system; Based on the preset thermodynamic equations and the initial parameters, the composition of the multi-component system is estimated, and the liquid phase fugacity coefficient and gas phase fugacity coefficient of each component are determined. Lower limits are imposed on the liquid phase fugacity coefficient and the gas phase fugacity coefficient, and the phase equilibrium constants of each component are determined based on the constrained liquid phase fugacity coefficient and gas phase fugacity coefficient. Based on the phase equilibrium constant and the feed mole fraction, a gas phase fraction solution equation is constructed. Based on the relative magnitude between the current gas phase fraction and the preset threshold, the target expression corresponding to the gas phase fraction solution equation is determined; The objective expression is iteratively solved in a hybrid manner to update the gas phase fraction until the gas phase fraction meets the preset convergence condition, thereby obtaining the corresponding gas-liquid flash calculation result; the gas-liquid flash calculation result indicates that the current phase composition of the multi-component system meets the gas-liquid two-phase equilibrium condition.

[0007] Optionally, the iterative solution of the target expression to update the gas phase fraction includes: Based on the current gas phase fraction, determine the function value, first derivative value, and second derivative value of the equation for solving the gas phase fraction; Based on the function value, the first derivative value, and the second derivative value, a first update step size is determined; the first update step size represents the gas phase fraction update amount based on the Halley method, and the first update step size has third-order convergence characteristics. When the first update step size meets the preset degradation condition, the first update step size is replaced with the second update step size; the second update step size represents the gas phase fraction update amount based on Newton's method, and the second update step size has second-order convergence characteristics; the degradation condition represents that the absolute value of the denominator of the first update step size is less than the first preset value, or the absolute value of the first update step size is greater than the second preset value. The gas phase fraction is updated based on the determined target update step size.

[0008] Optionally, the iterative solution of the target expression to update the gas phase fraction includes: The target update step size is pruned. The gas phase fraction is updated based on the target update step size after trimming; If the updated gas phase fraction exceeds the preset feasible range of gas phase fraction, the updated gas phase fraction is projected into the feasible range of gas phase fraction.

[0009] Optionally, the gas phase fraction satisfies a preset convergence condition, including: Based on the phase equilibrium constant and gas phase fraction in the current iteration and the previous iteration, respectively, determine the changes in the phase equilibrium constant and the gas phase fraction. When both the change in the phase equilibrium constant and the change in the gas phase fraction are less than the corresponding change threshold, it is determined that the gas phase fraction satisfies the preset convergence condition.

[0010] Optionally, the lower limit constraint on the liquid phase fugacity coefficient and the gas phase fugacity coefficient includes: The liquid phase fugacity coefficient and the gas phase fugacity coefficient are compared with preset lower limit values ​​respectively; The fugacity coefficient value that is less than the preset lower limit value is corrected to the preset lower limit value.

[0011] Optionally, determining the target expression form corresponding to the gas phase fraction solution equation based on the relative magnitude between the current gas phase fraction and a preset threshold includes: When the current gas phase fraction is less than a preset threshold, the gas phase fraction solution equation is constructed using the first expression; When the current gas phase fraction is not less than a preset threshold, the gas phase fraction solution equation is constructed using the second expression; wherein, the first expression is: The second expression is: .

[0012] Optionally, after obtaining the corresponding gas-liquid flash evaporation calculation results, the method further includes: Based on the gas-liquid flash evaporation calculation results, the current phase composition of the multi-component system is determined; Based on a preset tangent plane distance strategy, phase stability analysis is performed on the current phase composition to determine the tangent plane distance value corresponding to the current phase composition; When the tangential plane distance value is less than zero, a new phase is introduced into the multi-component system to extend the phase equilibrium calculation of the multi-component system.

[0013] In a second aspect, this application provides a computer device, including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement any one of the gas-liquid flash evaporation calculation methods described in the first aspect above.

[0014] Thirdly, this application provides a computer storage medium storing computer program instructions, which are executed by a processor using any one of the gas-liquid flash evaporation calculation methods described in the first aspect above.

[0015] Fourthly, an embodiment of this application provides a computer program product, including computer program instructions, which, when executed by a processor, implement any one of the gas-liquid flash evaporation calculation methods described in the first aspect above.

[0016] The beneficial effects of this invention are as follows: This application provides a gas-liquid flash evaporation calculation method, device, medium, and product. The method acquires initial parameters of a multi-component system, estimates the composition of the system based on a preset thermodynamic equation, determines the liquid-phase fugacity coefficient and gas-phase fugacity coefficient of each component, and imposes a lower limit constraint on the fugacity coefficient. Based on the constrained fugacity coefficient, it determines the phase equilibrium constant of each component. Using the phase equilibrium constant and the feed mole fraction, it constructs a gas phase fraction solution equation. Based on the relative magnitude between the current gas phase fraction and a preset threshold, it determines the target expression corresponding to the gas phase fraction solution equation. The target expression is then iteratively solved using a mixed solution to update the gas phase fraction until the gas phase fraction meets a preset convergence condition. This yields a gas-liquid flash evaporation calculation result representing the current phase composition of the multi-component system satisfying the gas-liquid two-phase equilibrium condition, thereby improving the efficiency and accuracy of the gas-liquid flash evaporation calculation. Attached Figure Description

[0017] To more clearly illustrate the technical solutions in the embodiments of this application or related technologies, the drawings used in the description of the embodiments or related technologies will be briefly introduced below. Obviously, the drawings described below are only embodiments of this application. For those skilled in the art, other drawings can be obtained based on the provided drawings without creative effort.

[0018] Figure 1 A schematic flowchart illustrating a gas-liquid flash evaporation calculation method provided in an embodiment of this application; Figure 2 A schematic diagram of a feasible domain protection process based on piecewise equations is provided for an embodiment of this application; Figure 3 A schematic diagram of an extreme system provided for an embodiment of this application; Figure 4 This application provides a schematic diagram of an automatic degradation process. Figure 5 A schematic diagram of an extreme working condition comparison experiment provided in this application embodiment; Figure 6 A schematic diagram illustrating the relationship between the lower bound constraint of fugacity and the boundedness of the second derivative is provided for an embodiment of this application. Figure 7 A schematic diagram of a unified flash evaporation solution framework provided for embodiments of this application; Figure 8 A schematic diagram of a phase stability analysis process provided in an embodiment of this application; Figure 9 This is a schematic diagram of the structure of a computer device provided in an embodiment of this application. Detailed Implementation

[0019] To make the objectives, technical solutions, and advantages of this application clearer, the technical solutions in the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of this application, and not all of them. Based on the embodiments of this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application. Unless otherwise specified, the embodiments and features in the embodiments of this application can be arbitrarily combined with each other. Furthermore, although a logical order is shown in the flowchart, in some cases, the steps shown or described may be performed in a different order than that shown here.

[0020] The terms "first" and "second" in the specification, claims, and accompanying drawings of this application are used to distinguish different objects, not to describe a specific order. Furthermore, the term "comprising" and any variations thereof are intended to cover non-exclusive protection. For example, a process, method, system, product, or device that includes a series of steps or units is not limited to the listed steps or units, but may optionally include steps or units not listed, or may optionally include other steps or units inherent to these processes, methods, products, or devices. The term "multiple" in this application can mean at least two, for example, two, three, or more, and this application does not impose limitations.

[0021] The term "and / or" in the embodiments of this application is merely a description of the association relationship between related objects, indicating that three relationships can exist. For example, A and / or B can represent: A existing alone, A and B existing simultaneously, and B existing alone. Additionally, the character " / " in this document generally indicates that the preceding and following related objects have an "or" relationship.

[0022] It is understood that the following specific embodiments of this application involve data related to chemical production processes, etc. When the various embodiments of this application are applied to specific products or technologies, relevant licenses or consents are required, and the collection, use, and processing of related data must comply with the relevant laws, regulations, and standards of the relevant countries and regions. For example, relevant volunteers can be recruited and agreements can be signed to authorize their data, thereby enabling the implementation using the data of these volunteers; or, implementation can be carried out within an authorized organization, using data from members of the organization to implement the following implementation methods for data management; or, the relevant data used in the specific implementation may be simulated data, such as simulated data generated in a virtual scene.

[0023] The design concept of the embodiments of this application is briefly introduced below: Gas-liquid flash evaporation calculation, as a fundamental computational unit in simulations of petrochemical, natural gas processing, coal chemical, and liquefied chemical processes, has the core task of solving the gas-liquid two-phase equilibrium state of multi-component systems based on thermodynamic equations, determining the gas phase fraction and the distribution of each component between the two phases. An efficient flash evaporation algorithm is crucial for ensuring the computational accuracy, convergence speed, and overall stability of process simulators, and is of great significance for the simulation and optimization of complex processes.

[0024] Currently, common calculations of gas-liquid flash evaporation mainly rely on second-order iterative methods such as Newton's method to solve for the phase equilibrium constant K and the Rachford-Rice equation (RR equation) to predict gas-liquid distribution. However, traditional methods suffer from serious numerical defects in practical engineering applications, especially under extreme conditions such as highly asymmetric component distribution (e.g., containing trace amounts of heavy or light components), near-critical conditions, or single-phase boundaries. First, because the thermodynamic equation of state may return extremely small fugacity coefficient values, the calculated phase equilibrium constant can span multiple orders of magnitude or even become ill-conditioned. This causes drastic changes or unbounded tendencies in the first and second derivatives of the Rachford-Rice equation, leading to instability in the iteration step size of Newton's method, which is prone to oscillations or divergence. Furthermore, the traditional direct introduction of higher-order iterative methods such as Halley's method in pursuit of faster convergence often backfires. Higher-order methods are more sensitive to derivatives, and under conditions of increased ill-conditioning in the aforementioned equations, the denominator of the iteration step size easily approaches zero, resulting in uncontrolled step size and numerical stability even worse than Newton's method. Furthermore, traditional techniques typically design independent solution processes for different flash evaporation modes such as constant pressure and temperature, and constant pressure and constant gas phase ratio. This results in fragmented code structures and a lack of a unified and robust interface for extending to more complex phase equilibrium systems of three phases and above, increasing the complexity of integrating, maintaining, and extending the algorithm in industrial software.

[0025] In view of the above problems, this application provides a gas-liquid flash evaporation calculation method. This method obtains the initial parameters of a multi-component system and estimates the composition of the system based on a preset thermodynamic equation, determining the liquid phase fugacity coefficient and gas phase fugacity coefficient of each component. Then, it applies lower bound constraints to the liquid and gas phase fugacity coefficients and determines the phase equilibrium constant of each component based on these constrained coefficients. Using the phase equilibrium constant and the feed mole fraction, a gas phase fraction solution equation is constructed. Based on the relative magnitude between the current gas phase fraction and a preset threshold, the target expression corresponding to the gas phase fraction solution equation is determined. The target expression is iteratively solved using a mixed solution to update the gas phase fraction until the gas phase fraction meets a preset convergence condition. This yields a gas-liquid flash evaporation calculation result representing the current phase composition of the multi-component system satisfying the gas-liquid two-phase equilibrium condition, thereby improving the efficiency and accuracy of the gas-liquid flash evaporation calculation.

[0026] Furthermore, this application's embodiments effectively address the numerical instability problem of traditional flash evaporation calculation methods under extreme conditions through a systematic numerical protection mechanism and an intelligent iterative strategy. The lower limit constraint on the fugacity coefficient avoids abnormal fluctuations in the phase equilibrium constant, the selection of piecewise expressions improves the condition number of the equations, and the hybrid iterative method balances convergence speed and stability. Thus, this application demonstrates excellent numerical stability and computational efficiency when dealing with complex multi-component systems, providing reliable technical support for industrial process simulation and optimization.

[0027] Furthermore, this application constrains the condition number of the RR equation through piecewise Rachford–Rice expressions (i.e., RR expressions), avoiding near-zero derivatives or numerical cancellation when the gas phase fraction β is close to 0 or 1, thus significantly improving the solvability of the equation. Moreover, this application also ensures the Hessian boundedness of the RR equation through interval constraints of the fugacity lower limit and the gas phase fraction β, providing a mathematical guarantee for the stability of the Halley method, rather than just an empirical safeguard.

[0028] Furthermore, this application achieves third-order convergence with feasible region protection through a hybrid automatic degradation and step-size pruning strategy combining the Halley and Newton methods. Specifically, it fully utilizes the third-order convergence advantage of Halley's method within regions with good condition numbers, and automatically reverts to Newton's method in ill-conditioned regions, balancing convergence speed and stability. Thus, in complex multi-component systems, compared to traditional single Newton iteration, the method provided in this application significantly reduces the number of iterations and improves the convergence success rate, especially exhibiting good robustness under extreme conditions such as heavy-tailed components, near-critical conditions, or near single-phase boundaries.

[0029] Furthermore, this application simultaneously supports three modes—pressure-temperature (PT), pressure-vapor fraction (PV), and temperature-vapor fraction (TV)—within a unified flash evaporation framework. It also reserves a multiphase extension interface based on tangent plane distance (TPD), making the algorithm kernel and interface structurally unified, which facilitates integration and expansion in actual process simulation software.

[0030] Please refer to Figure 1 The following is a flowchart illustrating a gas-liquid flash evaporation calculation method provided in an embodiment of this application. The specific implementation process of this method is as follows: Step 101: Obtain the initial parameters of the multi-component system.

[0031] In the embodiments of this application, the initial parameters of the multi-component system include at least two of the known variable parameters such as the pressure, temperature, gas phase fraction, and feed mole fraction of the multi-component system.

[0032] Specifically, in this embodiment, pressure refers to the ambient pressure of the system, temperature is a thermodynamic state parameter of the system, gas phase fraction is the mole fraction of the gas phase in the total material, and feed mole fraction is the molar proportion distribution of each component in the total feed. After obtaining these initial known variable parameters, this application will normalize and verify the feed mole fraction to ensure that the sum of the mole fractions of all components is one, thereby avoiding calculation deviations caused by material non-conservation. Next, according to different flash evaporation calculation modes, the corresponding known parameters and variables to be determined are identified. For example, in the constant pressure and constant temperature flash evaporation mode, pressure and temperature are known inputs, and gas phase fraction is the variable to be determined. In the constant pressure and constant gas phase fraction flash evaporation mode, pressure and gas phase fraction are known inputs, and temperature is the variable to be determined, thus providing complete initial conditions for subsequent thermodynamic calculations.

[0033] In one possible implementation scheme, embodiments of this application may first construct a unified flash evaporation framework, transforming the solution variables in the PT, PV, and TV modes into a unified form (denoted as gas phase fraction). (Temperature T, pressure P): For the PT flash mode: (P, T) are known, and unknown... ; For PV flash mode: it is known that (P, ), unknown T; For TV flash mode: it is known that (T, ), unknown P.

[0034] And define a unified state vector:

[0035] Next, implement the following uniformly in the base class Flash: Feed composition Verification and normalization; Thermodynamics package interface based on Equation of State (EOS): Calculation (L) , (V) And the K value; General iterative control (residual evaluation, step size trimming, convergence judgment, etc.).

[0036] Thus, PTFlash, PVFlash, and TVFlash, as derived classes, only represent different unknown variables by overloading the "master equation residual F(u)", thereby forming a unified code framework and state space in terms of structure, but the numerical kernel is completely reused.

[0037] Step 102: Based on the preset thermodynamic equations and initial parameters, estimate the composition of the multi-component system and determine the liquid phase fugacity coefficient and gas phase fugacity coefficient of each component.

[0038] In this embodiment, a preset thermodynamic equation is invoked to estimate the composition of the multi-component system, thereby determining the liquid phase fugacity coefficient and gas phase fugacity coefficient of each component.

[0039] In one possible implementation, the thermodynamic equations in this application are typically in the form of equations of state, capable of describing the thermodynamic behavior of a fluid under non-ideal conditions. This application will obtain the fugacity coefficient of each component in the liquid and gas phases through iterative calculations based on the current pressure, temperature, and the gas-liquid phase composition estimated from the feed mole fraction and gas phase fraction. The fugacity coefficient is an important parameter measuring the degree of deviation between a real fluid and an ideal gas; its calculation requires consideration of factors such as intermolecular forces, system pressure, and temperature.

[0040] Step 103: Apply lower limit constraints to the liquid phase fugacity coefficient and the gas phase fugacity coefficient, and determine the phase equilibrium constant of each component based on the constrained liquid phase fugacity coefficient and gas phase fugacity coefficient.

[0041] In this embodiment, the calculated liquid-phase fugacity coefficient and gas-phase fugacity coefficient are subject to a lower limit constraint. Each fugacity coefficient is compared to a preset minimum positive value, and the larger value is taken as the constrained fugacity coefficient. This effectively prevents numerical calculation problems caused by excessively small fugacity coefficients. Based on the constrained fugacity coefficients, the phase equilibrium constants of each component are further calculated. The phase equilibrium constant represents the ratio of the liquid-phase fugacity coefficient to the gas-phase fugacity coefficient, reflecting the component's distribution tendency between the liquid and gas phases. Thus, by introducing a lower limit constraint on the fugacity coefficients, this application ensures that the phase equilibrium constants are within a reasonable numerical range, avoiding extreme cases where the phase equilibrium constants tend to infinity or zero due to excessively small fugacity coefficients, and providing numerical stability assurance for subsequent equation solving.

[0042] In one possible implementation, the embodiments of this application compare the liquid phase fugacity coefficient and the gas phase fugacity coefficient with a preset lower limit value, thereby correcting the fugacity coefficient value that is less than the preset lower limit value to the preset lower limit value, thus achieving a lower limit constraint on the fugacity coefficient.

[0043] Specifically, in the embodiments of this application, the K value is determined by the fugacity coefficients of the liquid and gas phases. , The decision is based on the following calculations:

[0044] like or If K_i is too small, it may tend to infinity or 0, causing extreme values ​​for f' and f'' in the RR equation, leading to... Unstable. Therefore, this application will set a lower limit for the fugacity coefficient at the EOS interface level:

[0045] in, For example, a lower limit for the value of 1e-12.

[0046] Corollary 1: The value of K is bounded on both sides.

[0047] Because of the EOS Given that the upper bound of the equation is finite within the engineering scope, and considering the lower bound constraint, we can conclude that there exists a constant M > 0 such that for all i, we have:

[0048] Corollary 2: The Hessian term f''(β) is bounded.

[0049] Recalling f'':

[0050] because Bounded, ; and Therefore, for any i, we have It is a positive number, and from this we can obtain

[0051] In summary, by introducing a lower bound on fugacity and a β interval restriction, this application can theoretically guarantee that the RR equation has a bounded second derivative (Hessian bounded) during the iteration process, thus providing a denominator for the Halley method. The stability of [the system] provides a mathematical foundation.

[0052] Step 104: Based on the phase equilibrium constant and feed mole fraction, construct the equation for solving the gas phase fraction.

[0053] In this embodiment, a gas phase fraction solution equation is constructed based on the obtained phase equilibrium constant and feed mole fraction. This equation, based on the principle of material conservation, describes the distribution relationship of each component between the two phases under gas-liquid two-phase equilibrium conditions. Specifically, it is expressed as a nonlinear function of the gas phase fraction, the root of which is the desired gas phase fraction value. In constructing this equation, this application ensures that the material conservation relationship for all components is satisfied, i.e., the sum of the mole fractions of each component in the gas and liquid phases equals the mole fraction in the feed. This transforms the physical phase equilibrium problem into a mathematical equation-solving problem, laying the foundation for subsequent numerical solutions.

[0054] Specifically, in this application embodiment, the phase equilibrium constant will be calculated based on the fugacity coefficient. And solve the gas phase fraction β equation, i.e., the Rachford–Rice equation (RR equation):

[0055] Specifically, in the embodiments of this application, for a given K value and composition... The RR equation can be written as:

[0056] Its first and second derivatives are:

[0057] The above derivation is based on each term The derivatives can be obtained separately, and they have clear physical meanings: the first derivative reflects the sensitivity of the gas phase fraction to the material balance, while the second derivative reflects the rate of change of this sensitivity.

[0058] Furthermore, this application also performs condition number analysis on the RR equation: For scalar nonlinear equations , at the root Nearby, allow for a small perturbation in the value of f. The perturbation of the solution The first-order approximation is:

[0059] Therefore, the condition number relative to f can be defined as:

[0060] Note the root Therefore, when specifically used for RR analysis, a neighborhood approximation can be used to approximate the value of the RR function. Treating it as a linear function, the following equivalent indices can be used to characterize the ease of solving the RR equation:

[0061] Intuitive meaning: When When the value is very small (the derivative approaches zero), any perturbation δf will be amplified into a huge δβ, making the equation extremely difficult to solve; The larger the value, the more "ill-conditioned" the equation becomes. Substituting, we get:

[0062] When there exists a component j such that Larger or Since β is extremely small, close to 0 or 1, the summation terms in the denominator may be very close to 0 or numerically cancel each other out, making... The magnitude is very large, which makes the RR equation highly sensitive to disturbances and difficult to converge iteratively. This is the mathematical root cause of why the RR equation is "difficult to solve" under extreme conditions.

[0063] Therefore, to improve the aforementioned ill-conditioning, this application will express the RR equation in a piecewise form. By performing algebraic transformations on the RR equation, an equivalent form can be obtained: When β is close to 0, use the regular form:

[0064] When β is close to 1, the RR equation can be written as:

[0065] It is easy to deduce that:

[0066] The solution when the RR equations of both are equal to zero They are completely identical, but their numerical properties differ, especially the weighting method in the derivative expression: when β < 0.5, the focus is mainly on the convergence process starting from β = 0. To get closer to 1, using the f1 form avoids unnecessary subtraction between 1 and small quantities. When β≥0.5, starting from β=1 is more reasonable, and using the f2 form avoids derivative cancellation caused by the extreme value of the denominator when β→1.

[0067] By taking the derivative of f2 and performing a similar analysis, a new derivative can be obtained:

[0068] The main difference from f1' lies in the weighting factor, which, in some extreme cases, can significantly reduce the mutual cancellation in the derivative summation, such that:

[0069] That is, in By switching the equation expression at both ends, we can give... Provide a positive lower bound C, thus allowing the condition number to be... A finite upper bound is given. Therefore, this application indirectly constrains the condition number of the RR equation through a piecewise RR expression, maintaining acceptable sensitivity even when approaching the single-phase boundary.

[0070] Step 105: Based on the relative magnitude between the current gas phase fraction and the preset threshold, determine the target expression corresponding to the gas phase fraction solution equation.

[0071] In this embodiment, the specific expression form of the gas phase fraction solution equation will be selected based on the relative magnitude between the current gas phase fraction and the preset threshold.

[0072] In one possible implementation, when the current gas phase fraction is less than a preset threshold, the present application embodiment will use a first expression to construct a gas phase fraction solution equation.

[0073] When the current gas phase fraction is not less than the preset threshold, the second expression is used to construct the gas phase fraction solution equation.

[0074] The first expression is: ; The second expression is: .

[0075] Specifically, when the gas phase fraction is less than a preset threshold, the standard form of the equation will be used; when the gas phase fraction is greater than or equal to the preset threshold, an equivalent transformed form will be used. Thus, this piecewise expression strategy prevents numerical computational difficulties caused by the denominator approaching zero when the gas phase fraction is close to zero or one. Furthermore, by selecting an appropriate target expression form, this application ensures that the equation exhibits good numerical behavior throughout its entire domain, thereby improving the stability and efficiency of iterative solutions.

[0076] Specifically, taking a preset threshold β of 0.5 as an example, when β < 0.5, the aforementioned RR equation is used for solving; when β ≥ 0.5, the above equation is equivalently rewritten as:

[0077] In one possible implementation, refer to Figure 2 The diagram shown is a schematic of a feasible region protection process based on piecewise equations provided in this application embodiment. This application embodiment will input the current gas phase fraction β and feed mole fraction z. i and the phase equilibrium constant Ki Then, proceed to segmented selection and judgment, that is, determine whether β is less than 0.5. If so, use RR form 1, that is, calculate using the following formula:

[0078] If not (i.e., β≥0.5), then we switch to RR form 2, which is calculated using the following formula:

[0079] After determining the objective expression, the function value f(β), the first derivative f'(β), and the second derivative f''(β) of the equation under the current β are calculated. Finally, the feasible region protection iteration module is entered, which uses the previously calculated f(β), f'(β), and f''(β) to perform subsequent mixed iterations to update the β value until convergence. Thus, the... Figure 2 This application demonstrates how it intelligently selects equation forms based on β values ​​to improve numerical conditions and prepares the necessary data for higher-order iterations.

[0080] For details, please refer to Figure 3 The diagram shown is a schematic representation of an extreme system provided in an embodiment of this application, which can represent a shale gas or condensate gas system. =1.5 MPa, When K = 329K, the value of K given by EOS in a certain initial iteration is approximately:

[0081] Around β=0.98:

[0082] After switching using the segmented RR form of the embodiment of this application:

[0083] The absolute value of the derivative no longer approaches zero, and the condition number is significantly reduced. The experiment shows that it is the piecewise segmentation that prevents the numerical disappearance of the derivative term, rather than a problem caused by EOS precision or initial values.

[0084] Step 106: Iterate the objective expression using a hybrid solution to update the gas phase fraction until the gas phase fraction meets the preset convergence condition, and obtain the corresponding gas-liquid flash evaporation calculation results.

[0085] In this embodiment of the application, the determined objective expression will be solved by mixed iteration to update the gas phase fraction value until the gas phase fraction meets the convergence condition. The flash evaporation calculation result obtained after convergence will accurately represent the gas-liquid two-phase equilibrium state of the multi-component system under given conditions.

[0086] Specifically, in the iterative process of this application, higher-order iterative methods are preferentially used to calculate the update step size, while the numerical stability of the iterative process is monitored in real time. When it is detected that the iteration step size may cause numerical divergence, the system automatically switches to a lower-order iterative method. After each iteration update, the gas phase fraction is protected within a feasible region to ensure that it always remains within the physically meaningful range. The iterative process continues until the change in gas phase fraction and the change in phase equilibrium constant are both less than the preset convergence threshold, at which point the calculation is considered to have converged. The flash evaporation calculation results obtained after convergence accurately characterize the gas-liquid two-phase equilibrium state of the multi-component system under given conditions.

[0087] In one possible implementation, this embodiment of the application will determine the function value, first derivative value, and second derivative value of the gas phase fraction solving equation based on the current gas phase fraction; based on the function value, first derivative value, and second derivative value, determine a first update step size, which represents the gas phase fraction update amount based on the Halley method, and the first update step size has third-order convergence characteristics; when the first update step size meets a preset degradation condition, the first update step size will be replaced with a second update step size, which represents the gas phase fraction update amount based on the Newton method, and the second update step size has second-order convergence characteristics. The degradation condition indicates that the absolute value of the denominator of the first update step size is less than a first preset value, or the absolute value of the first update step size is greater than a second preset value, thereby updating the gas phase fraction based on the determined target update step size.

[0088] Specifically, in this application embodiment, the Rachford–Rice equation will be solved by a hybrid iteration of the Halley method and the Newton method with feasible region protection.

[0089] The first derivative of Halley's method is:

[0090] The second derivative of Halley's method is:

[0091] When f(β) = 0 and f'(β) ≠ 0, the embodiments of this application use a piecewise RR expression, which can guarantee that the absolute value of f' has a positive lower bound, thereby restricting the condition number:

[0092] Furthermore, this application also ensures that the second derivative of the RR equation satisfies the following through fugacity lower bound and β feasible region constraints:

[0093] in, It is a positive number.

[0094] Furthermore, this application will prioritize calculating the first update step size corresponding to the Halley method.

[0095] Where f'(β) and f''(β) are the first and second derivatives of the above equation, respectively, if they satisfy...

[0096] Then, the second update step size corresponding to Newton's method is used instead:

[0097] Specifically, in this application embodiment, given the RR equation f(β) and its derivative, a hybrid iterative method combining the Halley method and the Newton method is employed:

[0098] Newton's method has a convergence order of 2 and is relatively robust to initial conditions. Halley's method converges near the root. Under ideal conditions, the convergence order is 3, but for stable operation, two conditions must be met: the denominator... Not close to 0; The value of must be controllable; otherwise, the order of magnitude of molecules will seriously affect stability.

[0099] Therefore, the following feasible domain protection strategy is designed in the algorithm of this application embodiment: Condition detection and automatic degradation Each iteration calculates the Halley denominator ,like

[0100] In this step, the Halley method is not used for updating; instead, it automatically degenerates into a Newtonian step.

[0101] in This is a preset safety threshold.

[0102] Specifically, the Halley method autodegradation in this application embodiment is as follows: For the same working condition, if Halley's method is used directly, the denominator term is:

[0103] At β≈0.98, |D|<1×10 -9If the Halley step size exceeds 1.5, it will directly jump out of the feasible region and diverge.

[0104] After using the automatic degradation strategy in the embodiments of this application:

[0105] The algorithm maintains a range of 0.0 < β < 1.0, with no divergence throughout. Once β enters the stable convergence region, this application will reactivate Halley, and from step 6 onwards, only 3 iterations are needed to reach the convergence threshold.

[0106] In one possible implementation, the embodiments of this application will determine the change in phase equilibrium constant and the change in gas phase fraction based on the phase equilibrium constant and gas phase fraction in the current iteration and the previous iteration, respectively. When both the change in phase equilibrium constant and the change in gas phase fraction are less than the corresponding change threshold, it is determined that the gas phase fraction meets the preset convergence condition.

[0107] Specifically, the convergence criteria in this application include:

[0108] Thus, this application will simultaneously judge the changes in K and β based on the convergence criterion, and output the result when both are less than the set threshold.

[0109] For details, please refer to Figure 4 The diagram illustrates an automatic degradation process provided in this application. First, the application calculates the update step size of the Halley method and determines whether the absolute value of its denominator is less than a preset threshold. If so, it indicates that directly using the Halley step size may lead to numerical instability, therefore automatic degradation occurs, and the Newton method step size is adopted instead. If not, the calculated Halley step size continues to be used. Regardless of the step size used, it is recorded as a candidate step size Δβ, and it is checked whether the updated β value (i.e., β + Δβ) exceeds the preset feasible region [β_min, β_max]. If it is determined to exceed, the step size is truncated, and the updated β value is forcibly projected back into the feasible region. After step size truncation and projection, the β value is formally updated. Finally, a convergence check is performed, checking whether the maximum change in the phase equilibrium constant K is less than the corresponding threshold, and whether the change in the gas phase fraction is less than the corresponding threshold. Only when both conditions are met simultaneously will the iteration terminate and the final flash evaporation calculation result be output; otherwise, the next iteration will begin based on the new β value.

[0110] For details, please refer to Figure 5The diagram illustrates an extreme condition comparison experiment provided in this application embodiment. Under the same EOS and initial value (Wilson), this application compares three methods: the traditional Newton's method, the Halley method without infeasibility region protection, and the hybrid iterative method provided in this application embodiment. It is evident that the method provided in this application embodiment can iterate stably and quickly to convergence. The fundamental reason for this is that it switches to the higher-order Halley method at appropriate times, and automatically degenerates to the Newton's method when the conditions are not met. Furthermore, because this application limits the fugacity, it directly restricts the range of the phase equilibrium constant, indirectly ensuring that the ranges of the Hessian value and the Halley denominator are reasonable, avoiding numerical divergence, and increasing the applicability of the Halley method.

[0111] For details, please refer to Figure 6 The diagram illustrates the relationship between the lower bound constraint of fugacity and the boundedness of the second derivative, as provided in an embodiment of this application. This application applies a preset positive lower bound constraint to the liquid and gas phase fugacity coefficients obtained from thermodynamic calculations, thereby limiting the phase equilibrium constant within a reasonable finite range and avoiding extreme values. Based on the boundedness of Ki, the derivatives of the Rachford-Rice equation are bounded, i.e., the first derivative has a positive lower bound and the second derivative has an upper bound. This improves the condition number of the equation and lays the foundation for the application of higher-order methods, ensuring that the denominator of the Halley iteration formula is far from zero and avoiding singularities in step size calculations. Thus, the aforementioned series of synergistic measures achieve controlled iteration step size, avoiding numerical divergence, enabling the algorithm to stably and efficiently solve for the gas phase fraction even under extreme conditions.

[0112] In one possible implementation, the present application embodiment will trim the target update step size and update the gas phase fraction based on the trimmed target update step size, so that when the updated gas phase fraction exceeds the preset feasible range of gas phase fraction, the updated gas phase fraction will be projected into the feasible range of gas phase fraction.

[0113] Specifically, this application will trim the obtained update step size and project β onto the feasible region:

[0114] In one possible implementation, the feasible region constraint of β in this application is: ,in: That is, the candidate step size δ obtained in this application. If the updated Not here (For example, within [1e-4, 0.9999]), the step size is linearly scaled, and β is strictly projected onto the feasible interval:

[0115] Thus, based on the aforementioned improvement of condition number (piecewise RR) and step size protection, the Halley method only works in local regions where the derivative condition number is good and the Hessian is bounded. This maintains the advantage of third-order convergence while avoiding the problem of divergence in the traditional Halley method in RR equations.

[0116] In one possible implementation, after obtaining the gas-liquid flash calculation results, the current phase composition of the multi-component system is determined based on the gas-liquid flash calculation results. Based on a preset tangential plane distance strategy, a phase stability analysis is performed on the current phase composition to determine the tangential plane distance value corresponding to the current phase composition. Thus, when the tangential plane distance value is less than zero, a new phase is introduced into the multi-component system to extend the phase equilibrium calculation of the multi-component system.

[0117] Specifically, in this application embodiment, under a unified framework, the three flash evaporation problems of PT, PV and TV are all reduced to "solving a coupled nonlinear system of one-dimensional main variable u and K value", with the RR equation or its variant as the core.

[0118] Based on this, in order to support multiphase flash evaporation, this application will reserve the phase stability analysis interface proposeNewPhaseByStabilityAnalysis in the base class and adopt the tangent plane distance (TPD) criterion. That is, when the tangent plane distance TPD(y) < 0 is satisfied, a new phase is introduced and the components are redistributed into the multiphase solution.

[0119] Specifically, this application assumes the current phase composition is x (liquid or gas) and the total Gibbs free energy is G. For any experimental phase composition y, the TPD function is defined as follows:

[0120] in, It represents the chemical potential.

[0121] If there exists a y such that TPD(y) < 0, then the current phase is unstable and a new phase should be introduced. The specific steps are as follows: The composition of the new phase is initialized based on y where TPD < 0; The total component z is redistributed between the old and new phases to ensure mass conservation; The extended RR equations are a set of multiphase material balance equations. Under the same unified framework, the new β vector (mole fraction of each phase) and the K value of each phase queue are solved iteratively. Although the current implementation can activate only two-phase modes, the multiphase extension has a clear mathematical foundation and feasible path through the above definition, leaving a unified interface and algorithmic constraints for subsequent implementations.

[0122] In one possible implementation, refer to Figure 7 The diagram shown is a schematic of a unified flash evaporation solution framework provided in an embodiment of this application. It uses the base class Flash as the core numerical unit and includes: a feed verification and normalization module; and an EOS call interface for calculation. , The application includes a K-value calculation module; a main variable solution module (capable of calculating T / P / β); a feasible region-protected high-order iterative module; and a stability analysis interface. In this application, the PT, PV, and TV flash modes all employ a unified state space and iterative control, only changing the solution equations for the main variables β, T, or P. For example, when the main variable is T or P, the Wilson formula for inversely estimating the initial K-value is used to obtain the initial value. Thus, the PTFlash, PVFlash, and TVFlash flash calculation models in this application will all use Flash as the parent class, only overriding the main variable solution equations, allowing different flash modes to share a unified data structure and numerical flow without needing to switch between different algorithmic logics.

[0123] In one possible implementation, refer to Figure 8 The diagram illustrates a phase stability analysis process provided in this application embodiment. This application performs phase stability analysis to determine and expand multiphase equilibrium after obtaining preliminary convergence results from flash evaporation calculations. First, the application inputs the currently calculated phase state, including temperature T, pressure P, and composition x. Based on this current phase composition, a family of experimental compositions y is constructed through perturbation. According to the tangent plane distance criterion, the TPD(y) value corresponding to each experimental composition y is calculated, and it is determined whether there is an experimental phase with TPD(y) < 0. If so, the y that minimizes the TPD value is selected. The most unstable phase is used to initialize the new phase composition. Then, according to the principle of mass conservation, the total feed composition z is redistributed between the original phases and this new phase. Next, the extended multiphase material balance equations are constructed and solved to obtain the updated mole fractions βk and phase compositions xk of each phase. Finally, the updated number of phases and the multiphase equilibrium results are output. If the TPD(y) of all experimental compositions is not less than zero, it indicates that the current system is stable. The original number of phases will be maintained (e.g., two phases in this embodiment), and subsequent iterations will continue without introducing a new phase.

[0124] Specifically, in this embodiment of the application, for the above-mentioned multi-component system, proposeNewPhaseByStabilityAnalysis can be called, and the TPD function can be used:

[0125] During the iterative convergence phase, all perturbations conforming to y are detected, and TPD ≥ 0 is obtained. Therefore, the system can maintain two phases and does not require expansion to three phases.

[0126] Please see Figure 9 As shown, based on the same technical concept, this application also provides a computer device 90. In one embodiment, the computer device can be a device specifically for gas-liquid flash evaporation calculation, or it can be a device for overall control of industrial production. The computer device, as shown... Figure 9 As shown, it includes a memory 901, a communication module 903, and one or more processors 902.

[0127] The memory 901 is used to store computer programs executed by the processor 902. The memory 901 may mainly include a program storage area and a data storage area. The program storage area may store the operating system and programs required to run instant messaging functions, etc.; the data storage area may store various instant messaging information and operation instruction sets, etc.

[0128] Memory 901 may be volatile memory, such as random-access memory (RAM); memory 901 may also be non-volatile memory, such as read-only memory, flash memory, hard disk drive (HDD), or solid-state drive (SSD); or memory 901 may be any other medium capable of carrying or storing desired program code in the form of instructions or data structures and accessible by a computer, but is not limited thereto. Memory 901 may be a combination of the above-mentioned memories.

[0129] The processor 902 may include one or more central processing units (CPUs) or digital processing units, etc. The processor 902 is used to implement the above-described gas-liquid flash evaporation calculation method when it calls the computer program stored in the memory 901.

[0130] The communication module 903 is used to communicate with the industrial control system.

[0131] This application embodiment does not limit the specific connection medium between the memory 901, communication module 903, and processor 902 described above. This application embodiment... Figure 9 The memory 901 and the processor 902 are connected via a bus 904, which is in... Figure 9 The diagram uses thick lines to describe the connections between other components; these are for illustrative purposes only and should not be considered limiting. The 904 bus can be divided into address bus, data bus, control bus, etc. For ease of description, Figure 9 It is described using only a thick line, but does not indicate that there is only one bus or one type of bus.

[0132] The memory 901 stores a computer storage medium, which stores computer-executable instructions. The computer-executable instructions are used to implement the gas-liquid flash evaporation calculation method of the embodiments of this application, and the processor 902 is used to execute the gas-liquid flash evaporation calculation method of the above embodiments.

[0133] Based on the same inventive concept, embodiments of this application also provide a storage medium storing a computer program that, when run on a computer, causes the computer to execute the steps in the gas-liquid flash evaporation calculation method according to various exemplary embodiments of this application described above.

[0134] In some possible implementations, various aspects of the gas-liquid flash evaporation calculation method provided in this application can also be implemented in the form of a computer program product, which includes a computer program that, when run on a computer device, causes the computer device to perform the steps in the gas-liquid flash evaporation calculation method according to various exemplary embodiments of this application as described above. For example, the computer device can perform the steps of the various embodiments.

[0135] The program product may employ any combination of one or more readable media. A readable medium may be a readable signal medium or a readable storage medium. A readable storage medium may be, for example, but not limited to, an electrical, magnetic, optical, electromagnetic, infrared, or semiconductor system, apparatus, or device, or any combination thereof. More specific examples of readable storage media (a non-exhaustive list) include: electrical connections having one or more wires, portable disks, hard disks, random access memory (RAM), read-only memory (ROM), erasable programmable read-only memory (EPROM or flash memory), optical fiber, portable compact disk read-only memory (CD-ROM), optical storage devices, magnetic storage devices, or any suitable combination thereof.

[0136] The program product of the embodiments of this application may employ a portable compact disc read-only memory (CD-ROM) and include a computer program, and may run on a computer device. However, the program product of this application is not limited thereto. In this application, the readable storage medium may be any tangible medium that contains or stores a program, and the computer program included therein may be used by or in conjunction with a command execution system, apparatus, or device.

[0137] A readable signal medium may include a data signal propagated in baseband or as part of a carrier wave, carrying a readable computer program. This propagated data signal may take various forms, including but not limited to electromagnetic signals, optical signals, or any suitable combination thereof. A readable signal medium may also be any readable medium other than a readable storage medium, capable of sending, propagating, or transmitting a program for use by or in conjunction with a command execution system, apparatus, or device.

[0138] Computer programs contained on readable media may be transmitted using any suitable medium, including but not limited to wireless, wired, optical fiber, RF, etc., or any suitable combination thereof.

[0139] Computer programs for performing the operations of this application can be written in any combination of one or more programming languages, including object-oriented programming languages ​​such as Java and C++, as well as conventional procedural programming languages ​​such as the "C" language or similar programming languages.

[0140] It should be noted that although several units or sub-units of the device have been mentioned in the detailed description above, this division is merely exemplary and not mandatory. In fact, according to embodiments of this application, the features and functions of two or more units described above can be embodied in one unit. Conversely, the features and functions of one unit described above can be further divided and embodied by multiple units.

[0141] Furthermore, although the operations of the method of this application are described in a specific order in the accompanying drawings, this does not require or imply that these operations must be performed in that specific order, or that all the operations shown must be performed to achieve the desired result. Additionally or alternatively, certain steps may be omitted, multiple steps may be combined into one step, and / or one step may be broken down into multiple steps.

[0142] Those skilled in the art will understand that embodiments of this application can be provided as methods, systems, or computer program products. Therefore, this application can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, this application can take the form of a computer program product embodied on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0143] Although preferred embodiments of this application have been described, those skilled in the art, upon learning the basic inventive concept, can make other changes and modifications to these embodiments. Therefore, the appended claims are intended to be interpreted as including the preferred embodiments as well as all changes and modifications falling within the scope of this application.

[0144] Obviously, those skilled in the art can make various modifications and variations to this application without departing from the spirit and scope of this application. Therefore, if such modifications and variations fall within the scope of the claims of this application and their equivalents, this application also intends to include such modifications and variations.

Claims

1. A method for calculating gas-liquid flash evaporation, characterized in that, The method includes: Obtain initial parameters for the multi-component system, wherein the initial parameters include at least two of the following: pressure, temperature, gas phase fraction, and feed mole fraction of the multi-component system; Based on the preset thermodynamic equations and the initial parameters, the composition of the multi-component system is estimated, and the liquid phase fugacity coefficient and gas phase fugacity coefficient of each component are determined. Lower limits are imposed on the liquid phase fugacity coefficient and the gas phase fugacity coefficient, and the phase equilibrium constants of each component are determined based on the constrained liquid phase fugacity coefficient and gas phase fugacity coefficient. Based on the phase equilibrium constant and the feed mole fraction, a gas phase fraction solution equation is constructed. Based on the relative magnitude between the current gas phase fraction and the preset threshold, the target expression corresponding to the gas phase fraction solution equation is determined; The objective expression is iteratively solved in a hybrid manner to update the gas phase fraction until the gas phase fraction meets the preset convergence condition, thereby obtaining the corresponding gas-liquid flash calculation result; the gas-liquid flash calculation result indicates that the current phase composition of the multi-component system meets the gas-liquid two-phase equilibrium condition.

2. The method as described in claim 1, characterized in that, The iterative solution of the objective expression to update the gas phase fraction includes: Based on the current gas phase fraction, determine the function value, first derivative value, and second derivative value of the equation for solving the gas phase fraction; Based on the function value, the first derivative value, and the second derivative value, a first update step size is determined; the first update step size represents the gas phase fraction update amount based on the Halley method, and the first update step size has third-order convergence characteristics. When the first update step size meets the preset degradation condition, the first update step size is replaced with the second update step size; the second update step size represents the gas phase fraction update amount based on Newton's method, and the second update step size has second-order convergence characteristics; the degradation condition represents that the absolute value of the denominator of the first update step size is less than the first preset value, or the absolute value of the first update step size is greater than the second preset value. The gas phase fraction is updated based on the determined target update step size.

3. The method as described in claim 2, characterized in that, The iterative solution of the objective expression to update the gas phase fraction includes: The target update step size is pruned. The gas phase fraction is updated based on the target update step size after trimming; If the updated gas phase fraction exceeds the preset feasible range of gas phase fraction, the updated gas phase fraction is projected into the feasible range of gas phase fraction.

4. The method as described in claim 1, characterized in that, The gas phase fraction satisfies preset convergence conditions, including: Based on the phase equilibrium constant and gas phase fraction in the current iteration and the previous iteration, respectively, determine the changes in the phase equilibrium constant and the gas phase fraction. When both the change in the phase equilibrium constant and the change in the gas phase fraction are less than the corresponding change threshold, it is determined that the gas phase fraction satisfies the preset convergence condition.

5. The method as described in claim 1, characterized in that, The lower limit constraint on the liquid phase fugacity coefficient and the gas phase fugacity coefficient includes: The liquid phase fugacity coefficient and the gas phase fugacity coefficient are compared with preset lower limit values ​​respectively; The fugacity coefficient value that is less than the preset lower limit value is corrected to the preset lower limit value.

6. The method as described in claim 1, characterized in that, The step of determining the target expression form corresponding to the gas phase fraction solution equation based on the relative magnitude between the current gas phase fraction and a preset threshold includes: When the current gas phase fraction is less than a preset threshold, the gas phase fraction solution equation is constructed using the first expression; When the current gas phase fraction is not less than a preset threshold, the gas phase fraction solution equation is constructed using the second expression; wherein, the first expression is: The second expression is: .

7. The method as described in claim 1, characterized in that, After obtaining the corresponding gas-liquid flash evaporation calculation results, the method further includes: Based on the gas-liquid flash evaporation calculation results, the current phase composition of the multi-component system is determined; Based on a preset tangent plane distance strategy, phase stability analysis is performed on the current phase composition to determine the tangent plane distance value corresponding to the current phase composition; When the tangential plane distance value is less than zero, a new phase is introduced into the multi-component system to extend the phase equilibrium calculation of the multi-component system.

8. A computer device, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes the computer program, it implements the steps of the method according to any one of claims 1 to 7.

9. A computer storage medium storing computer program instructions thereon, characterized in that, When executed by a processor, the computer program instructions implement the steps of the method according to any one of claims 1 to 7.

10. A computer program product comprising computer program instructions, characterized in that, When executed by a processor, the computer program instructions implement the steps of the method according to any one of claims 1 to 7.

Images