Solution thermodynamic equilibrium calculation system and method based on DFP gradient algorithm

Abstract

The invention belongs to a solution system, particularly relates to the technical field of thermodynamic equilibrium calculation of a complex solution system containing a plurality of chemical (such as polycondensation clustering) reactions, and discloses a solution thermodynamic equilibrium calculation system and method based on a DFP gradient algorithm. According to the method, constraint variables are converted into independent variables, on the basis of the Gibbs free energy minimization thought, the problem of solving multiple chemical equilibrium is converted into the problem of unconstrained optimization, and the DFP gradient algorithm is properly selected to solve the problem. The invention provides a solution thermodynamic equilibrium calculation method and system based on a DFP gradient algorithm. The system is composed of a preprocessing module, a solution thermodynamic equilibrium calculation module and a DFP optimization calculation module, and solution multiple chemical equilibrium calculation is carried out. Compared with other algorithms, the DFP gradient optimization algorithm has the advantages that the calculation process is simplified, the convergence and calculation time is effectively shortened, the calculation efficiency and the calculation flux are improved, and higher calculation precision can be obtained as required.

Description

Technical Field

[0001] The present invention relates to the technical field of thermodynamic equilibrium calculation of multi-species (cluster) solution systems with multiple chemical reactions, and in particular to a system and method for calculating thermodynamic equilibrium of a solution based on a DFP gradient algorithm. Background Art

[0002] Thermodynamic models of chemical reactions are suitable for describing complex solution systems of multiple species (clusters) generated by multiple chemical reactions. Especially for salt or acid solution systems containing poly (tungsten, molybdenum, vanadium, chromium, zirconium, hafnium, aluminum, phosphorus, silicon, boron, arsenic, etc.) high-valent oxygen-containing acid radical ion clusters, it has a wide range of applications in catalysis, life sciences, material synthesis and hydrometallurgy. Due to the complex structure of such solutions, in-depth research on such solutions has both theoretical significance and practical application value. Thermodynamic description of such solutions is the premise and basis for in-depth understanding and research of such solutions.

[0003] For this type of solution system, it is often difficult to effectively handle or extremely inefficient to use a variety of state-type thermodynamic models including the Pitzer model. The chemical reaction thermodynamic model has become a powerful description and calculation tool that is difficult to match with other thermodynamic models because it can handle chemical reactions very well. However, solving multiple chemical equilibria when applying this model is an algorithmic problem, and it is also a bottleneck for whether the model can smoothly describe such complex solutions. Although there are currently a variety of algorithms available to solve this problem, there are still some problems to varying degrees in terms of algorithm complexity, computational efficiency, and computational accuracy. In the prior art, the chemical reaction thermodynamic model in the Chinese invention patent with application number ZL201410151779.7 calls optimization programs such as donlp2 or SNOPT to solve multiple chemical equilibrium algorithms, which are not specifically designed for chemical reaction thermodynamic models, and their algorithms are relatively complex and their computational efficiency is relatively low. Summary of the invention

[0004] The present invention and the Chinese invention patent with application number ZL 201410151779.7 are technologies of the same line, and the purpose is to overcome the defects of the prior art, combine the chemical reaction thermodynamic model, and select the DFP unconstrained gradient algorithm to solve multiple chemical equilibrium problems, so as to help the chemical reaction thermodynamic model to quickly and efficiently describe the complex solution system, and is particularly suitable for solving the thermodynamic properties of complex solutions containing condensation clustering reactions, where the system is not easy to converge and the amount of calculation is large.

[0005] The technical solution adopted by the present invention to achieve the above-mentioned purpose is: for a complex solution system with multiple species (clusters) formed by multiple chemical reactions, a Gibbs free energy equation (i.e., a chemical reaction thermodynamic model) capable of describing the thermodynamic properties of the solution is constructed, and the DFP unconstrained gradient algorithm is selected to calculate the minimum Gibbs free energy to achieve the purpose of solving multiple chemical equilibria. When the Gibbs free energy of the solution system reaches the minimum value, the solution reaches thermodynamic equilibrium. At this time, according to the amount of substance X of each species (cluster) in the equilibrium state, i , and then the properties of the series equilibrium states of the mixed system can be calculated.

[0006] The improved technical points and beneficial effects of the present invention are as follows

[0007] (1) In the prior art, the chemical reaction thermodynamic model program externally calls optimization programs such as donlp2 or SNOPT to solve multiple chemical equilibrium algorithms; in comparison, the external optimization program of the present invention is compared with the directly written DFP native program (the algorithm is relatively simple and easy to program), and the calculation program of the present invention is more direct, more efficient, and easier to maintain and expand.

[0008] (2) The present invention is based on a chemical reaction thermodynamic model and selects to adopt the DFP gradient unconstrained optimization algorithm. In order to adapt to this algorithm, the present invention uses independent variables instead of constrained variables, which reduces the number of variables, simplifies the algorithm, greatly reduces the complexity of the algorithm and the program, and improves the calculation efficiency.

[0009] (3) Compared with other gradient algorithms, the correction matrix H in the DFP gradient algorithm of the present invention is continuously iterated. k Finally, it converges to the second-order partial derivative at the extreme point, that is, This ensures that the correction matrix H k Positive definiteness, and can converge stably and quickly to the minimum value g of the objective function (i.e., the Gibbs free energy of the solution) min . BRIEF DESCRIPTION OF THE DRAWINGS

[0010] In order to make the content of the present invention more clearly understood, the present invention is further described in detail below according to specific schemes and in conjunction with the accompanying drawings, wherein

[0011] Figure 1 .DFP optimization calculation process diagram;

[0012] Figure 2 .Calculation flow chart and logical relationship between modules;

[0013] In the figure: I is the pretreatment program module; II is the solution thermodynamic equilibrium calculation program module; III is the DFP optimization calculation program module. DETAILED DESCRIPTION

[0014] 1. Thermodynamic description of complex solution system

[0015] A type of complex solution system, such as an acid solution system containing phosphate, tungstate, molybdate, silicate, etc., is characterized by the self-association of oxygen-containing acid anions into chains to form polynucleic acid radical polyanions (clusters) in different forms. This type of solution has a wide range of applications in biochemistry, materials science, nanotechnology, hydrometallurgy and other fields. Due to the complex structure of this type of solution, in-depth research on this type of solution has both theoretical significance and practical application value. Thermodynamic description of this type of solution is the premise and basis for in-depth understanding and research of this type of solution. In order to describe this type of solution, taking the simplest single electrolyte aqueous solution system as an example, the solution is composed of the simplest three components: anions, cations and water molecules. All species (clusters) present in the solution can be regarded as being generated by reactions with different stoichiometric coefficients. Taking concentrated phosphoric acid solution as an example, The species (clusters) can be regarded as composed of components H + , H2O is generated according to a certain reaction coefficient, such as It can be seen as a combination of 2 H + , 3 Reacts with 3 H2O to generate. Taking tungstic acid solution as an example, different tungsten-containing species (clusters) such as WO3(aq), It can be seen that H + , H2O is generated according to a certain reaction coefficient, such as It can be seen as a combination of 9 H + , 7 It is generated by reacting with -4 H2O. In particular, in the polyacid solution system, the mononuclear oxygen acid anion will lose a certain amount of water during the condensation process of the polyanion (cluster), and the generated species (cluster) may also carry a certain number of water molecules, so the reaction coefficient of the water molecule can be positive or negative.

[0016] In a more general way, in a single electrolyte-aqueous solution complex system, the reaction of the basic components to generate species (clusters) can be expressed as formula (1)

[0017]

[0018] It is represented as species (clusters) It can be regarded as X n- and H2O is generated by reaction. Where i = 1, 2, ... n, n is the number of H2O molecules that do not contain the basic component M m+ , X n-The number of species (clusters) with H2O.

[0019] After the reaction is balanced, the species (cluster) The amount of substance is X i , then the amount of each species (cluster) generated is X1, X2, …, X i , X n These species (clusters) are related to the basic component M m+ , X n- There is a material conservation relationship between the amount of H2O and the amount of substance.

[0020] X a +∑r ai ·X i =N a X b +∑r bi ·X i =N b X c +∑r ci ·X i =N c Formula (2)

[0021] In the above formula, N a 、N b With N c The basic components M m+ , X n- Before the reaction with H2O (the amount of each species (cluster) at this time is X i The total amount of substance is 0, and It is not difficult to calculate the basic component M after the reaction m+ , X n- The amount of substances remaining with H2O are and

[0022] Therefore, according to the thermodynamic model of electrolyte solution chemical reaction, the Gibbs free energy g of the solution system can be expressed as

[0023] g=g s +g DH Formula (3)

[0024]

[0025] In the above formula, g s is the Gibbs free energy of the solution when the solution is ideally mixed, g DH is the long-range electrostatic interaction energy of the solution. n t 、V t, I are the standard Gibbs free energy of species (clusters) and basic components (cations and water molecules), the total amount of substance, the total volume (in L) and the ionic strength. R is the gas constant, with a value of 8.31. T is the Kelvin temperature, and the room temperature is 298.15K. A and b are constants, with values ​​of 1.178 and 1.579, respectively. j = a, b, c, i corresponds to the components cations, anions, water molecules and species (cluster) i subscripts, species (cluster) i refers to the basic components according to a certain stoichiometric coefficient and The generated particles.

[0026] In particular, where V i With z j is the molar volume (in L / mol) and charge of the particle, n t 、V t and I are respectively represented as follows

[0027]

[0028] According to the above formula, since the basic components can be expressed as X j To represent the cation Anions Water molecules X i Here, it is the independent variable, and the subscript i traverses all species (clusters) generated by the basic component reactions. The use of independent variables can reduce the use of variables and bring simplicity to the calculation. i Re-express n as the independent variable t 、V t and I are as follows

[0029] n t =(N a -∑r ai X i )+(N b -∑r bi X i )+(N c -∑r ci X i )+∑X i Formula (8)

[0030] V t =(N a -∑r ai X i )V a +(N b -∑r bi X i )V b +(N c -∑rci X i )V c +∑X i V i Formula (9)

[0031]

[0032] Similarly, g s It can also be further expressed as

[0033]

[0034] Based on formulas (3), (5), and (11), the combined Gibbs free energy of the solution is expressed as:

[0035]

[0036]

[0037] Formula (3)g s The amount of substance X for all particles j Differentiate (j=a, b, c, i) and we get

[0038]

[0039] Formula (5) g DH The amount of substance X for all particles j Differentiate (j=a, b, c, i) and we get

[0040]

[0041] Considering the use of independent variables, n t 、V t The expression of I is re-expressed from formula (7) as formula (8-10), where n t 、V t and I respectively to X i (i is the species (cluster) generated by the reaction of the basic components, the same below) partial differentiation to obtain

[0042]

[0043] Therefore, g expressed by the independent variable s (See formula (11)) for X i Partial differentiation gives

[0044]

[0045] Similarly, the long-range electrostatic interaction energy g DH X i Partial differentiation gives

[0046]

[0047]

[0048] Also because

[0049]

[0050] so

[0051]

[0052] So the total Gibbs free energy of the solution is

[0053]

[0054] The above formula a i 、a a 、a b with a c are the activities of species (clusters) and basic components, cations, anions and water molecules

[0055]

[0056] G i , G a , G b With G c The chemical potential or Gibbs free energy per unit mole of the species (cluster) and the basic components cations, anions and water molecules It is the potential difference between the two ends of the equal sign of the reaction equation during the reaction process.

[0057] When the solution system reaches thermodynamic equilibrium, g = g min , the partial derivative of each species (cluster) is 0, which means

[0058]

[0059] Therefore, there is

[0060]

[0061] If lnK r ,lnK c It is expressed as follows

[0062]

[0063] Then further

[0064]

[0065] Combining the reaction formula (1) of species (cluster) generation with equation (25), it is not difficult to see that the potential difference on both sides of the reaction equal sign is In the reaction process, the amount of substance X is not equal. For example, in the process of species generation, the amount of substance X is i is gradually increased, and the former chemical potential Increase, while the chemical potential The two are equal when equilibrium is reached. Naturally, when the Gibbs free energy of the system reaches a minimum, g = g min At this time, the solution system reaches thermodynamic equilibrium and the reactions between species (clusters) are in a state of chemical equilibrium.

[0066] 2. DFP algorithm calculates thermodynamic equilibrium of complex solution system

[0067] From an algorithmic point of view, calculating the thermodynamic equilibrium of a complex solution system is essentially to minimize the total Gibbs free energy of the solution system and calculate the minimum value of the objective function, that is, formula (3). i ) is a nonlinear complex multivariate function. Since its variables are independent variables and the solution has the minimum Gibbs free energy, the function has the characteristics of being unconstrained and having a single extreme value. Its partial derivative is an n-dimensional vector, which has been obtained in the previous analysis process. Its second-order partial derivative is It is an n×n matrix, and its analytical expression is still difficult to find.

[0068] According to these characteristics, the Gibbs free energy of the solution system, i.e. the objective function, is suitable for calculating the minimum value of the function using a gradient-related optimization algorithm. Comparing several different gradient algorithms, the steepest descent method only converges to the first order, and its search paths advance in a perpendicular zigzag path, so its convergence speed is relatively slow. Other gradient methods require the calculation of second-order derivatives. Should It is difficult to calculate in this calculation, so it is not suitable. By comparison, the present invention uses the DFP algorithm in the quasi-Newton method combined with the chemical reaction model to calculate the thermodynamic equilibrium of the complex solution system. This algorithm does not require the calculation of the second-order derivative However, it has a second-order convergence characteristic and thus has better search efficiency. The DFP algorithm was first proposed by Davidon (1959) and later improved by Fletcher and Powell (1963). It is one of the most effective methods in unconstrained optimization. A brief introduction and derivation of the DFP algorithm is as follows.

[0069] Omitting higher-order infinitesimals, g(X) in X (k+1) The Taylor expansion of

[0070]

[0071] make And differentiate with respect to X to get

[0072]

[0073] Let X = X (k) , then the above formula becomes

[0074]

[0075] That is, Newton's iteration formula. Among them, Θ(X) is the second-order derivative of the function, but there are two problems with Θ(X). First, for complex target functions, its second-order derivative is generally difficult to obtain, and second, it needs to be ensured to be positive definite. In order to ensure the positive definiteness of Θ(X), it is often replaced by Θ(X)+I(X), where I(X) is the unit matrix. Since it is generally difficult to obtain Θ(X), the commonly used method is to use the correction matrix H k To replace Θ(X), take H k =Θ -1 (X)

[0076] Let δ k =X (k+1) -X (k) , Then δ k =H k γ k

[0077] Let the directional derivative

[0078] In order to further achieve the optimal search effect, it is further modified to X (k+1) -X (k) =t k P k , where t k The step parameter can be expressed by g(X (k) +t k P k ) Differentiate to obtain

[0079]

[0080] When g(X (k) +t k P k ) for t k When the differential is 0, solve the k The expression of t can be obtained k However, due to g(X (k) +t k P k ) is a complex expression and is difficult to obtain by differentiation. kIt can also be obtained by an approximate method, for g(X k +t k P k ) is expanded by Taylor second-order approximation to obtain

[0081]

[0082] Since dg(X (k) +t k P k ) / dt k =0, so the above formula is for t k Differential

[0083]

[0084] Substitution Further obtain t k expression

[0085]

[0086] However, the above k There are still some problems in the algorithm process of the approximate calculation formula. In the actual operation process, the approximation method is often used to more accurately obtain t k value

[0087] t k =argming(X (k) +t k P k ) Formula (37)

[0088] Structure H k Iteration formula H k =H k-1 +E k , and let So there is

[0089]

[0090] where u k 、v k They are the undetermined n-dimensional column vectors, a k ≠0, β k ≠0 is the unknown coefficient. All the terms in the above formula are multiplied by γ k have to

[0091]

[0092] Substitute δ k =H k γ k ,have to

[0093]

[0094] In the above formula, u k 、v k There are many ways to get Get

[0095]

[0096] Since matrices and vectors cannot be divided directly, the above formula is further changed to

[0097]

[0098] Substituting the above formula (38) into (37) yields

[0099]

[0100] In particular, H0 is initialized to be equal to the identity matrix I. Since H0 is positive definite, each subsequent H k must be positive; because H k Positive definite, the direction vectors generated during the iteration are different, so P i T H(X)P j =0, (i≠j), when X (k) Converges to X min At the extreme point, H n =G -1 (X) With the correction matrix, we can further find that when k = 0 and k = 1, we take t0 = 1 and t k The expression

[0101] t k =argming(X (k) +t k P k )

[0102] In particular,

[0103] 3. Solution system calculation process and system combined with DFP algorithm

[0104] 3.1 Calculation process and pseudo code

[0105] ①(S1-1) Define variables:

[0106] The total amount of each component N a 、N b 、N c and the amount of substance of each species X{X i}; Standard Gibbs free energy of each component And the standard Gibbs free energy of each species (cluster) Activity coefficient γ of each component a , γ b , γ c and the activity coefficients of each species (cluster) γ{γ i}; Activity of each component a a 、a b 、a c and the activity of each species (cluster) a{a i}; Molar volume of each component V a 、V b 、V c and the molar volume of each species (cluster) V{V i}; Charge z of each component a 、z b 、z c and the charge of each species (cluster) z{z i}; Generate the reaction coefficient r of each species (cluster) a {r ai}、 r c {r ci} and the reaction equilibrium constant K of each species (cluster) r {K r，i} and K c {K c，i}; gas constant R, Kelvin temperature T and Debye-Huckel coefficients A and b;

[0107] ②(S1-2) Declare sub-function:

[0108] The objective function g(X) is calculated using formula (12):

[0109]

[0110] Where X = {X i}The vector is the amount series of each species (cluster);

[0111] The activity coefficient is calculated using formula (20-23) as follows:

[0112]

[0113] The activity is calculated using formula (26):

[0114]

[0115] The partial derivative of the objective function is calculated using formula (25):

[0116]

[0117] Where X = {X i} vector is the amount series of each species (cluster); use formula (7) to calculate n t 、V t ,I

[0118]

[0119] ③(S1-3) Input data initialization list:

[0120] N a 、N b 、N c , {X i};

[0121]

[0122] V a 、V b 、V c , {V i}; z a 、z b 、z c 、{z i};

[0123] {r ai}、 {r ci};

[0124] R, T, A, b;

[0125] ④(S2-1) Call the DFP program to calculate the minimum value of the solution Gibbs free energy (i.e., the objective function) g(X) ⑤(S3-1) Define the DFP program variables:

[0126] Correction matrix H; direction vector P, vector δ, γ; step size adjustment parameter t;

[0127] Calculation accuracy ε, loop variable k, maximum number of iterations I max ;

[0128] ⑥(S3-2) Set ε(0<ε<<1) and I max , the optimal process initialization:

[0129] k=0, calculate g(X),

[0130] According to {X i}Assignment, determine X (0) starting point;

[0131] Set H0 = I. t0=1;

[0132] ⑦

[0133]

[0134] ⑧(S3-4) Output the minimum value of the objective function g(X), X={X i}Optimal value;

[0135] ⑨(S2-2) Calculate the properties of the solution system when the solution is in thermodynamic equilibrium:

[0136] Solution Gibbs free energy g(X) and its

[0137] γ a , γ b , γ c , {γ i}; a a 、a b 、a c 、{a i};

[0138] Equilibrium constant value K r , K c

[0139]

[0140] 3.2 Computing system composition

[0141] like Figure 2 As shown, the system includes the entire computing system including: S1. pre-processing program module, S2. solution thermodynamic equilibrium calculation program module, S3. DFP optimization calculation program module. The contents of each program module are as follows

[0142] S1. Preprocessing program module: It consists of a variable definition program segment (S1-1), a sub-function calculation program segment (S1-2) and a data initialization list (S1-3).

[0143] In the variable definition block (S1-1):

[0144] The variables defined include the solution Gibbs free energy g and its partial derivatives The total amount of each component N a 、N b 、N c and the amount of each species (cluster) X{X i}; Standard Gibbs free energy of each component And the standard Gibbs free energy of each species (cluster) Activity coefficient γ of each component a , γ b , γ c and the activity coefficients of each species (cluster) γ{γi}; Activity of each component a a 、a b 、a c and the activity of each species (cluster) a{a i}; Molar volume of each component V a 、V b 、V c and the molar volume of each species (cluster) V{V i}; Charge z of each component a 、z b 、z c and the charge of each species (cluster) z{z i}; Generate the reaction coefficient r of each species (cluster) a {r ai}、 r c {r ci} and the reaction equilibrium constant K of each species (cluster) r {K r，i} and K c {K c，i}; gas constant R, temperature T in Kelvin and Debye-Huckel coefficients A and b.

[0145] Sub-function declaration program segment (S1-2): including the calculation of the objective function g(X); the calculation of the activity coefficient lnγ i , lnγ a , lnγ b , lnγ c ; Activity calculation lna i , lna a , lna b , lna c ; Partial derivative of the objective function calculate

[0146] Enter the initialization data list (S1-3):

[0147] N a 、N b 、N c , {X i};

[0148]

[0149] V a 、V b 、V c , {V i}; z a 、z b 、z c 、{z i};

[0150] {rai}、 {r ci};

[0151] R, T, A, b;

[0152] S2. Solution thermodynamic equilibrium calculation program module: It consists of the solution system thermodynamic equilibrium process calculation program segment (S2-1) and the equilibrium system property calculation program segment (S2-2). The details are as follows

[0153] Solution system thermodynamic equilibrium process calculation program segment (S2-1): mainly calls the DFP optimization program to calculate the minimum Gibbs free energy of the solution to solve multiple chemical equilibrium; according to the calculation results of the DFP optimization program, the minimum value and extreme point of g(X) are returned. i}.

[0154] Equilibrium system property calculation program segment (S2-2): From X = {X i}Calculate the concentration of each species (cluster), the activity coefficient γ of each component and each species (cluster) a , γ b , γ c , {γ i}Calculation; Activity of each component and species (cluster) a a 、a b 、a c 、{a i}Calculate; calculate the equilibrium constant K of each reaction c and other solution properties;

[0155] S3.DFP optimization calculation program module: It consists of the optimal process variable definition program segment (S3-1), parameter setting and optimal process initialization program segment (S3-2), optimization iteration cycle process program segment (S3-3) and calculation result output program segment (S3-4). The details are as follows

[0156] DFP program variable definition program segment (S3-1): correction matrix H; direction vector P, vector δ, γ, step adjustment parameter t; calculation accuracy ε, loop variable k, maximum number of iterations I max ;

[0157] Setting parameters and optimal process initialization program segment (S3-2): Setting ε and I max ; Set k = 0, calculate g(X), According to {X i}Assignment determines the starting point X (0) ; Set H0 = I, t0=1;

[0158] Optimization iterative cycle process program segment (S3-3): Perform the following calculation X through k=k+1 cycles k+1 =X k +t k P k , P k =-H k γ k

[0159]

[0160] t k =argming(X (k) +t k P k )

[0161]

[0162] Calculation result output (S3-4): Output the minimum value of the objective function g(X) and determine the optimal value X={X i};

[0163] In particular, the calculation process of the DFP optimization calculation program module is as follows: Figure 1 shown.

[0164] 4. Algorithm and system calculation examples

[0165] Referring to Table 1 and Table 2, based on the above embodiments, this embodiment takes the polytungstic acid solution as an example and calculates the thermodynamic equilibrium of the solution based on the DFP gradient algorithm.

[0166] Table 1. Reaction stoichiometric coefficients, equilibrium constants and standard reaction Gibbs free energy changes of tungstate species (clusters)

[0167]

[0168] Note: In the table, there are negative reaction coefficients for H and O, because the reaction of the generated species does not consume water but generates water.

[0169] Table 2. Theoretical values ​​of K for chemical equilibrium at different concentrations and pH calculation accuracy r and the calculated value K c (K c , K r See calculation formulas (29) and (30))

[0170]

[0171] Note: ΔK i =K c，i -K r，i , See formula (31)

[0172] Taking polytungstate solution as an example, according to the open literature of Xiao Chao et al. [Xiao Chao, Xiao Liansheng, Cao Zuoying. Study on the mechanism of separation of trace tungsten from molybdate solution [J]. China Molybdenum Industry, 2011, 35(02): 29-32], the main species (clusters) generated by the solution are set as WO3(aq) is considered to be composed of basic components H + , H2O reacts with a certain stoichiometric coefficient. The reaction stoichiometric coefficients, logarithmic values ​​of equilibrium constants, and standard Gibss free energy changes of each species (cluster) are generated. Each is listed in Table 1.

[0173] The chemical reaction thermodynamics model was used to calculate the DFP gradient optimization algorithm for the above-mentioned tungstate solution system. a The calculation results under the conditions of pH, acidity, and calculation accuracy Eps are shown in Table 2. It can be seen from Table 2 that the chemical reaction thermodynamic model combined with the DFP gradient algorithm can better calculate the set solution thermodynamic equilibrium. The higher the calculation accuracy requirement set by the program, the closer the chemical equilibrium constant Kc of each species (cluster) calculated is to the theoretical calculated value Kr, and the smaller the difference between the two. Even if the solution concentration ma is as low as 0.005 mol / kg, its calculation accuracy can easily reach 1.0e-10.

[0174] In summary, the present invention adopts the above technical solution and has the following technical advantages: the method and system for calculating thermodynamic equilibrium of complex solution system based on DFP gradient algorithm provided by the present invention converts constraint variables into independent variable method, converts the problem of solving multiple chemical equilibrium into optimization problem based on the idea of ​​minimizing Gibbs free energy, and appropriately selects DFP gradient algorithm to solve this problem. Compared with other algorithms, DFP gradient optimization algorithm simplifies the calculation process, effectively shortens the convergence and calculation time, improves the calculation efficiency, and can obtain high calculation accuracy as needed. Compared with the prior art, the method and system for calculating thermodynamic equilibrium of complex solution system based on DFP gradient algorithm provided by the present invention calls optimization programs such as donlp2 or SNOPT to solve multiple chemical equilibrium algorithms. Compared with the directly written DFP native program, this external optimization program is relatively simple and easy to program. The method provided by the present invention is more direct in calculation, more efficient in calculation, and easier to maintain and expand.

[0175] In addition, a class of salt or acid complex solution systems containing poly (tungsten, molybdenum, vanadium, chromium, zirconium, hafnium, aluminum, phosphorus, silicon, boron, arsenic, etc.) oxygen-containing acid anions clusters, characterized by self-association between oxygen-containing acid anions into chains, rings or stacks, forming poly acid anion clusters with different polymerization degrees and different structures, and complex solution systems with multiple chemical equilibria are suitable for the technology of the present invention when dealing with multiple chemical equilibria, solution properties, especially thermodynamic properties. By extension, multi-species solution systems with multiple chemical reactions are also suitable for the technology of the present invention when dealing with multiple chemical equilibria, solution properties, especially thermodynamic properties.

[0176] The above is a detailed description of the specific implementation of the present invention in conjunction with the accompanying drawings. The specific implementation of this article is only used to help understand the method of the present invention. For ordinary technicians in this technical field, according to the idea of ​​the present invention, changes and modifications can be made in the specific implementation and application scope, so the description should not be understood as limiting the present invention.

Claims

1. A solution thermodynamic equilibrium calculation system based on DFP gradient algorithm, characterized in that: The solution thermodynamic equilibrium calculation system based on the DFP gradient algorithm includes a preprocessing module S1, a solution thermodynamic equilibrium calculation module S2 and a DFP optimization calculation module S3.

2. The solution thermodynamic equilibrium calculation system based on the DFP gradient algorithm according to claim 1, characterized in that: The preprocessing module S1 performs preparation work for the thermodynamic equilibrium calculation of the DFP gradient algorithm solution, and is composed of a variable definition program segment S1-1, a sub-function calculation program segment S1-2, and a data initialization list S1-3; in S1-1, variable definition program segment: define variables, the variables include the total amount of each component N a 、N b 、N c And the amount of substance X of each species (cluster) i ; Standard Gibbs free energy of each component and species (cluster) Activity coefficient γ of each component and species (cluster) a , γ b , γ c , {γ i }; Activity of each component and species (cluster) a a 、a b 、a c 、{a i }; Generate the response coefficient of each species (cluster) {r ai }、 {r ci } and reaction equilibrium constant K r , K c ; Gas constant R; Kelvin temperature T and Debye-Huckel coefficients A and b; S1-2, sub-function calculation program segment: used to define sub-function calculation, the sub-function includes: using formula (12): To calculate the solution Gibbs free energy g(X), use formula (25): To calculate the partial derivative of Gibbs free energy Using formula (20): Calculate the activity coefficient {γ i } and each component γ a , γ b , γ c , use formula (26) to calculate the activity a of each component a 、a b 、a c and the activity of each species (cluster) {a i }; S1-3, data initialization list: the data in the list includes: The total amount of each component: N a 、N b 、N c and the initial amount of each species {X i }; Standard Gibbs free energy of each component: And the standard Gibbs free energy of each species (cluster) Molar volume of each component: V a 、V b 、V c and the molar volume of each species (cluster) {V i }; Charge of each component: z a 、z b 、z c and the charge of each species (cluster) {z i }; Generate the reaction coefficient of each species (cluster): {r ai }、 {r ci }.

3. The solution thermodynamic equilibrium calculation system based on the DFP gradient algorithm according to claim 1, characterized in that: The solution thermodynamic equilibrium calculation module S2 is used to calculate the thermodynamic properties of the solution system. The calculation module consists of a solution thermodynamic equilibrium process calculation program segment and an equilibrium system property calculation program segment; in S2-1, solution thermodynamic equilibrium process calculation program segment: call DFP optimization algorithm to calculate the minimum Gibbs free energy of solution to solve multiple chemical equilibrium; return the minimum value and extreme point of g(X) X={X i }; S2-2, Equilibrium system property calculation program segment: used to calculate the properties of the solution system when the solution is in thermodynamic equilibrium, by X = {X i }Calculate the concentration of each species (cluster) and the activity coefficient γ of each component and each species (cluster) a , γ b , γ c , {γ i }; Calculate the activity a of each component and each species (cluster) a 、a b 、a c 、{a i }; Use formula (29): Calculate the equilibrium constant K for each reaction c and K r .

4. The solution thermodynamic equilibrium calculation system based on the DFP gradient algorithm according to claim 1, characterized in that: The DFP optimization calculation module S3 calculates the minimum value of the objective function, i.e., the Gibbs free energy g(X) of the solution, based on the DFP gradient optimization algorithm, and returns the amount of each species (cluster) when the minimum value is reached, i.e., when the solution is in equilibrium. i }; The DFP optimization calculation module S3 includes a DFP gradient algorithm variable definition program segment, a parameter setting and initialization program segment, an optimization iterative cycle process program segment, an optimization iterative cycle process program segment and an optimization calculation result output program segment; wherein: S3-1, DFP gradient algorithm variable definition program segment: correction matrix H; direction vector P, vector δ, γ, step adjustment parameter t; calculation accuracy ε, loop variable k, maximum number of iterations I max ; S3-2, parameter setting and initialization program segment: setting ε and I max ; Set k = 0, calculate g(X), According to {X i }Assignment determines the starting point X (0) ; Set H0 = I, t0=1; S3-3, Optimization iterative cycle process program segment: through k<I max , k=k+1 and perform the following calculations in a loop: X k+1 =X k +t k P k ,P k =-H k γ k δ k =X k -X k-1 , t k =argminf(X k +t k P k ) When k>I max or When , the loop terminates the calculation; S3-4, optimization calculation result output program segment: output the minimum value of the objective function g(X) and determine the optimal value.

5. A method for calculating the thermodynamic equilibrium of a solution based on a DFP gradient algorithm, characterized in that: The calculation is performed using the preprocessing module S1, the solution thermodynamic equilibrium calculation module S2 and the DFP optimization calculation module S3 in the solution thermodynamic equilibrium calculation system based on the DFP gradient algorithm according to any one of claims 1 to 4, specifically including: Execute the preprocessing module S1. The preprocessing module first executes S1-1 to prepare variable definitions, then executes S1-2 to declare sub-functions, and then executes S1-3 to obtain initialization data by inputting or reading the initialization data list; After the preprocessing module S1 is executed, the solution thermodynamic equilibrium calculation module S2 is continued to run. S2 executes S2-1 to call the DFP optimization calculation module S3 to calculate the minimum Gibbs free energy g of the solution system. min (X), get g min The amount of each species (cluster) in thermodynamic equilibrium at (X) = {X i }; After the DFP optimization calculation module S3 is completed, the solution thermodynamic equilibrium calculation module S2 is returned to execute. The solution thermodynamic equilibrium calculation module S2 executes S2-2 and uses the amount of each species (cluster) returned X = {X i }Calculate the thermodynamic properties of solution equilibrium systems; The DFP optimization calculation module S3 is executed by first executing S3-1 to define variables, then executing S3-2 to set and initialize parameters, and then executing S3-3 to perform optimization iteration loop. When the number of iteration loops k < I max or Terminate the calculation and continue to execute S3-4 to output the minimum value of the objective function g(X) and determine the optimal value X={X i }.

6. The method for calculating solution thermodynamic equilibrium based on the DFP gradient algorithm according to claim 5, characterized in that: The Gibbs free energy g(X) of the solution is calculated using the following formula (12): The partial derivative of the objective function g(X) is given by the following formula (25): Calculation sub-function: The activity coefficients of each species (cluster) and each component are calculated using the following formulas (20), (21), (22), and (23): The following formula (26) is used to calculate the activity of each species (cluster) and each component: The following formulas (8), (9) and (10) are used to calculate the total amount of substance n t , total volume V measured in L t , I is the ionic strength: n t =(N a -∑r ai X i )+(N b -∑r bi X i )+(N c -∑r ci X i )+∑X i Formula (8) V t = (N a - ∑r ai X i )V a + (N b - ∑r bi X i )V b + (N c - ∑r ci X i )V c + ∑X i V i Formula (9)

Images