Calculation Procedure of SCOZA
In the previous chapter, we have explained the equations required to numerically calculate the excess internal energy u(p, of a fluid.
In this chapter, we explain the calculation procedure for actually finding и using these equations.
Method to Determine the Potential Tail
It was shown in Chapter 4 that there are many intermolecular interactions that reproduce the same equation of state. Here, taking фС1 (r) as an example among them, the procedure for calculating the thermodynamic properties of the fluid is specifically described. First, how to determine the function form of ci(r) is described. In this book, we consider that water molecules are made of hard spheres with potential tails consisting of soft repulsive force and attractive force, and explore what shape tails reproduce well the thermodynamic properties of water. It is impossible to know the shape of such tails in advance. Therefore, it is required that work is done on several guessed tails to find the optimal interaction that reproduces the experimental results well while adjusting the shape. This is similar to exploring an ideal airplane shape that minimizes air resistance.
The Physics of Liquid Water Makoto Yasutomi
Copyright © 2021 Jenny Stanford Publishing Pte. Ltd.
ISBN 978-981-4877-25-1 (Hardcover), 978-1-003-05616-4 (eBook) www.jennystanford.com
Figure 6.1 A method to express the interaction between water molecules with multi-Yukawa potential. The guessed shape of the interaction between water molecules is shown by a broken line. Under the condition that the potential tail is made up of five Yukawa terms, the shape that best reproduces the dashed line is 0cl(r), shown in red.
As an example, the guessed tail shape is shown as a broken line in Fig. 6.1. Under the condition that tail 0cl is made of five Yukawa terms, the function фсi(r) that best fits the broken line is shown in red.
To confirm whether the potential tail ci(r) obtained in this way can reproduce the thermodynamic properties of water with high accuracy, it is necessary to calculate all the required physical quantities using SCOZA and thermodynamic related equations. Therefore, it is necessary to repeat the trial calculation several times while adjusting the shape of the potential tail that generates the soft repulsive force and attractive force to find the optimum shape. The potential cl(r) is a tail obtained as a result of repeating such a trial calculation several times.
Κ(1)(ρ) and z1(ρ)
The direct correlation function given by Eq. (5.79) of the hard sphere fluid is obtained by [Pini (1998)]. The quantities K^lp] and zx{p)
included in the formula are given by the following equation:
where,
Here /red is the reduced isothermal compressibility of the hard sphere fluid and is given by the following equation (see Fig. 6.2):
The quantity pp/р is given by Carnahan-Starling's empirical formula (5.76) (see blue line in Fig. 5.3). Also,
Figure 6.3 z^p).
The coefficients Zi(p) and /CW(p) calculated from the above equations are numerically unstable in the region where p is sufficiently small, and give occasionally incorrect values. Therefore, in such a region, the quantities Zi(p) and KPp) should be expanded to a power series of= лр/6 as shown below.
In the region of >7 < 10-2, we use the above power series expansion by >1 for zx and /CM. The obtained Zi{p) and Klp} are shown in Figs. 6.3 and 6.4, respectively. The density increment was set to Ар = 10-4.
Figure 6.5 @{p, 0).
Calculation of u(ρ, β = 0)
When p = 0 (Г ->• oo), the intermolecular interaction cl(r) is negligible compared to the thermal energy per molecule. In this case, the fluid can be regarded as a hard sphere fluid. As explained in Section 5.7.2, for a hard sphere fluid, 3>„ = 0 (2 < n < N = 6). Therefore, the values of ($?i, 3>{) can be obtained from Eqs. (5.128) and (5.129) using the Newton-Raphson method. The values of (2 < n < 6) can be obtained from Eq. (5.130). The value of ^(p, 0) obtained in this way is shown in Fig. 6.5.
The solutions (1 < n < 6) are shown in Fig. 6.6. Substituting %(p, 0) (2 < n < 6) into Eq. (5.123) gives Uhs = w(p, 0). The result is shown in Fig. 6.7. Substituting uHs into the diffusion equation
Figure 6.6 &„(p, 0) (1 < n < 6). Red line: ^(p, 0), Black line: &2(P, 0), Blue line: &3(p, 0), Green line: 0), Dashed line: ^(p, 0), Dotted line:
%(P, 0).
Figure 6.7 uHS = u(p, 0).
(5.124), the quantity u(p, f}) at the next temperature step = Д/Ь can be obtained.
Diffusion coefficient B{p, uhs) is obtained from Eq. (5.125). The resulting £?(p, uHs) is shown in Fig. 6.8. The methods to find A, dA/dx„, and dxn/du included in that equation are explained in Appendix D and Appendix E. In the density region of p > pj, there is no stable solution of the diffusion equation (5.124) because B{p, uHs) is negative. So we set the boundary density to p = pj (see Table 5.1). Therefore, the density region that should be considered is given by 0 < p < pj. In this case, the boundary condition (5.136) is given by the following equation:
Figure 6.8 S(p, uHs). #(Pi, uHs) = 176176.0 and S(p|, uHs) = 5.77 x 10 s.
Calculation up to Critical Temperature Tc
You can use 9n{p, p = 0), %{p, p = 0) (1 < n < N), uHs(p) = u[p, P = 0), and B(p, uhs) as initial values. If you solve the diffusion equation (5.124) using the method described in Section 5.8, the excess internal energy u{p, pk) for any temperature step pk can be obtained. The behavior of B{p, и) obtained in this way near the critical point (pc, /6C) is shown in Fig. 6.9. You can see how B(p, u) suddenly approaches 0 around p = pc = 0.0841 as p approaches pc.
Figure 6.9 The behavior of В(р, и) near the critical point (pc, pc). The curves represent the isotherms of $(p, u) for p = 0.1025,0.10776,0.10787, 0.107886, 0.107891 and pc = 0.107891313 from the top to the bottom, respectively.
Figure 6.10 Up, pc). Red line: [p, pc), Black line: $f2(p, Pc), Blue line: 'Up, pc), Green line: 'Up, Pc), Dashed line: ^(p, Pc), Dotted line: %(P, Pc)-
Figure 6.11 %{p, Д). Red line: &i{p, pc), Black line: &г{Р, Pc), Blue
line: ShSjp, Pc), Green line: &^{p, pc), Dashed line: @s(p, pc), Dotted line:
9ь{р. Pc)-
The solutions Уп[р, Pc) and 9>„{p, pc) (1 < n < 6) at the critical temperature P = pc are shown in Figs. 6.10 and 6.11.
Calculation below Critical Temperature
As explained in Section 5.7.3, below the critical temperature Tc, an unstable gas-liquid coexistence region appears in the middle density region where B[p, и) < 0. Gas is stable on the lower density side and liquid is stable on the higher density side than that. In order to obtain a solution in the stable region, a spinodal curve satisfying the boundary condition A{ps, P) = A[ps2, P) = 0 [see Eq. (5.131)] must be solved. The values of u(ps,) (i = 1, 2.) on the spinodal
Figure 6.13 9„(psi) (1 < n < 6). Red line: i^i(psi), Black line: S>2{Psi).
Blue line: ^(Psi), Green line: &^{psi), Dashed line: @s(psi), Dotted line:
^6(Psi)<
curve are shown in Fig. 6.12. The blue line on the lower density side than pc = 0.0841 represents u(psi), and the red line on the higher density side represents u(pS2). The density pS2 has an upper limit value of 0.3356, and the value of u(p52) diverges in the density region pS2 > 0.3356.
The values of 3>n[Psi) on the spinodal curve are shown in Fig. 6.13, @n(pS2) in Fig. 6.14, &„(pSi) in Fig. 6.15, and ^„(pS2) in Fig. 6.16, where (1 < n < 6).
Figure 6.15 The same as Fig. 6.13, but for%(psi) (1 < n < 6).
Figure 6.16 The same as Fig. 6.13, but for ^п(р52) (1 5 n < 6).
Calculations of Gas Phase Water below the Critical Temperature
Calculations for steam below the critical temperature can be made using the method described in Section 5.8.2. First, ^{p, ft), £*(/>, ft)] and y£ (7>si), are use<^ as approximate solutions
for the Newton-Raphson method to find solutions for the next temperature step = ft + Aft Similarly, if the solutions for this temperature are used as approximate solutions for the next temperature step, solutions for each temperature step can be
Figure 6.17 S(p, u) of gas phase water at typical temperatures fi = 0.1154 (orange line), 0.1414 (purple line), 0.2684 (green line), and 0.5954 (black line).
obtained one after another. The solutions B(p, u) and u(p, p) obtained in this way are shown for representative temperatures p = 0.1154, 0.1414, 0.2684, and 0.5954 in Figs. 6.17 and 6.18, respectively.
For the same representative temperatures as above, the solutions S>n[.P, P) and $„{p, fi) (1 < n < 6) are also exhibited in Figs. G.l, G.2, G.3, G.4, G.5, G.6, G.7, and G.8, respectively, in Section G.l in Appendix G.
Figure 6.18 The same as Fig. 6.17, but for u(p, P).
Calculations of Liquid Phase Water below the Critical Temperature
Liquid water below the critical temperature can be calculated by the method described in Section 5.8.3. The solutions of Eq. (5.124) are obtained in the same way as the calculation method for vapor phase water, using the solutions of the previous temperature step as approximate solutions.
The solutions B[p, и) and u(p, p) are shown for the same representative temperatures as previous subsection, respectively, in Figs. 6.19 and 6.20.
Figure 6.19 The same as Fig. 6.17, but for liquid water.
Figure 6.20 The same as Fig. 6.18, but for liquid water.
The solutions 3>„{p, p) and^„(p, p) (1 < n < 6) are exhibited in Figs. G.9, G.10, G.ll, G.12, G.13, G.14, G.15,and G. 16, respectively, in Section G.2 in Appendix G.
