|
6 | 6 | "source": [ |
7 | 7 | "# Cubic Equations of State\n", |
8 | 8 | "\n", |
9 | | - "Computing compressibility factors and molar volumes using the Peng-Robinson and Soave-Redlich-Kwong cubic equations of state for pure methane and a methane-ethane mixture." |
| 9 | + "Simulating the compressibility factor of methane across a pressure sweep using the Peng-Robinson and Soave-Redlich-Kwong equations of state wired into a PathSim simulation." |
10 | 10 | ] |
11 | 11 | }, |
12 | 12 | { |
13 | 13 | "cell_type": "markdown", |
14 | 14 | "metadata": {}, |
15 | 15 | "source": [ |
16 | | - "## Setup\n", |
| 16 | + "Cubic equations of state compute the compressibility factor $Z = Pv/(RT)$ by solving a cubic polynomial at each $(T, P)$ condition. The EoS blocks take two inputs (temperature and pressure) and produce two outputs (molar volume and compressibility factor).\n", |
17 | 17 | "\n", |
18 | | - "Import the equation of state blocks." |
| 18 | + "At low pressure $Z \\to 1$ (ideal gas). At moderate pressure attractive forces cause $Z < 1$, and at high pressure excluded-volume effects push $Z > 1$." |
19 | 19 | ] |
20 | 20 | }, |
21 | 21 | { |
|
24 | 24 | "metadata": {}, |
25 | 25 | "outputs": [], |
26 | 26 | "source": [ |
27 | | - "import numpy as np\n", |
28 | 27 | "import matplotlib.pyplot as plt\n", |
29 | 28 | "\n", |
| 29 | + "from pathsim import Simulation, Connection\n", |
| 30 | + "from pathsim.blocks import Source, Constant, Scope\n", |
| 31 | + "\n", |
30 | 32 | "from pathsim_chem.thermodynamics import PengRobinson, RedlichKwongSoave" |
31 | 33 | ] |
32 | 34 | }, |
33 | 35 | { |
34 | 36 | "cell_type": "markdown", |
35 | 37 | "metadata": {}, |
36 | 38 | "source": [ |
37 | | - "## Pure Methane\n", |
38 | | - "\n", |
39 | | - "The Peng-Robinson EoS solves a cubic equation in the compressibility factor $Z = Pv/(RT)$:\n", |
| 39 | + "## System Definition\n", |
40 | 40 | "\n", |
41 | | - "$$Z^3 - (1-B)Z^2 + (A - 3B^2 - 2B)Z - (AB - B^2 - B^3) = 0$$\n", |
42 | | - "\n", |
43 | | - "where $A$ and $B$ are dimensionless parameters computed from the critical properties. Each block takes $(T, P)$ as inputs and returns $(v, Z)$." |
| 41 | + "Create Peng-Robinson and SRK blocks for pure methane. A `Constant` block supplies a fixed temperature while a `Source` sweeps pressure from 0.1 MPa to 30 MPa." |
44 | 42 | ] |
45 | 43 | }, |
46 | 44 | { |
|
50 | 48 | "outputs": [], |
51 | 49 | "source": [ |
52 | 50 | "# Critical properties of methane\n", |
53 | | - "Tc_CH4 = 190.6 # K\n", |
54 | | - "Pc_CH4 = 4.6e6 # Pa\n", |
55 | | - "omega_CH4 = 0.011 # acentric factor\n", |
56 | | - "\n", |
57 | | - "# Create EoS blocks\n", |
58 | | - "pr = PengRobinson(Tc=Tc_CH4, Pc=Pc_CH4, omega=omega_CH4)\n", |
59 | | - "rks = RedlichKwongSoave(Tc=Tc_CH4, Pc=Pc_CH4, omega=omega_CH4)\n", |
60 | | - "\n", |
61 | | - "def eval_eos(block, T, P):\n", |
62 | | - " \"\"\"Evaluate an EoS block and return (v, Z).\"\"\"\n", |
63 | | - " block.inputs[0] = T\n", |
64 | | - " block.inputs[1] = P\n", |
65 | | - " block.update(None)\n", |
66 | | - " return block.outputs[0], block.outputs[1]\n", |
67 | | - "\n", |
68 | | - "# Evaluate at 300 K, 1 atm\n", |
69 | | - "v_pr, Z_pr = eval_eos(pr, 300, 101325)\n", |
70 | | - "v_rks, Z_rks = eval_eos(rks, 300, 101325)\n", |
71 | | - "\n", |
72 | | - "print(f\"Methane at 300 K, 1 atm:\")\n", |
73 | | - "print(f\" PR: Z = {Z_pr:.6f}, v = {v_pr:.6e} m³/mol\")\n", |
74 | | - "print(f\" RKS: Z = {Z_rks:.6f}, v = {v_rks:.6e} m³/mol\")\n", |
75 | | - "print(f\" Ideal gas: Z = 1.0, v = {8.314462 * 300 / 101325:.6e} m³/mol\")" |
| 51 | + "Tc, Pc, omega = 190.6, 4.6e6, 0.011\n", |
| 52 | + "\n", |
| 53 | + "# EoS blocks\n", |
| 54 | + "pr = PengRobinson(Tc=Tc, Pc=Pc, omega=omega)\n", |
| 55 | + "rks = RedlichKwongSoave(Tc=Tc, Pc=Pc, omega=omega)\n", |
| 56 | + "\n", |
| 57 | + "# Fixed temperature, logarithmic pressure sweep\n", |
| 58 | + "import numpy as np\n", |
| 59 | + "T_const = Constant(250) # 250 K (above Tc, supercritical)\n", |
| 60 | + "P_src = Source(func=lambda t: 10**(4 + t * 0.035)) # 10 kPa to ~30 MPa over 100s\n", |
| 61 | + "\n", |
| 62 | + "# Scopes: record Z from both EoS (output port 1)\n", |
| 63 | + "scp_pr = Scope(labels=[\"Z_PR\"])\n", |
| 64 | + "scp_rks = Scope(labels=[\"Z_RKS\"])" |
76 | 65 | ] |
77 | 66 | }, |
78 | 67 | { |
79 | 68 | "cell_type": "markdown", |
80 | 69 | "metadata": {}, |
81 | 70 | "source": [ |
82 | | - "## Compressibility Factor vs Pressure\n", |
| 71 | + "## Wiring\n", |
83 | 72 | "\n", |
84 | | - "At low pressure, $Z \\to 1$ (ideal gas). As pressure increases, intermolecular forces cause $Z$ to deviate. At very high pressures, repulsive forces dominate and $Z > 1$." |
| 73 | + "Both EoS blocks receive the same $(T, P)$ inputs. We record the compressibility factor (output port 1) from each." |
85 | 74 | ] |
86 | 75 | }, |
87 | 76 | { |
|
90 | 79 | "metadata": {}, |
91 | 80 | "outputs": [], |
92 | 81 | "source": [ |
93 | | - "T_fixed = 250 # K (above Tc for methane, supercritical)\n", |
94 | | - "P_range = np.logspace(4, 7.5, 100) # 10 kPa to ~30 MPa\n", |
95 | | - "\n", |
96 | | - "Z_pr_arr = []\n", |
97 | | - "Z_rks_arr = []\n", |
98 | | - "\n", |
99 | | - "for P in P_range:\n", |
100 | | - " _, Z = eval_eos(pr, T_fixed, P)\n", |
101 | | - " Z_pr_arr.append(Z)\n", |
102 | | - " _, Z = eval_eos(rks, T_fixed, P)\n", |
103 | | - " Z_rks_arr.append(Z)\n", |
| 82 | + "sim = Simulation(\n", |
| 83 | + " blocks=[T_const, P_src, pr, rks, scp_pr, scp_rks],\n", |
| 84 | + " connections=[\n", |
| 85 | + " # Temperature -> both EoS (input port 0)\n", |
| 86 | + " Connection(T_const, pr, rks),\n", |
| 87 | + " # Pressure -> both EoS (input port 1)\n", |
| 88 | + " Connection(P_src, pr[1], rks[1]),\n", |
| 89 | + " # Z output (port 1) -> scopes\n", |
| 90 | + " Connection(pr[1], scp_pr),\n", |
| 91 | + " Connection(rks[1], scp_rks),\n", |
| 92 | + " ],\n", |
| 93 | + " dt=1.0,\n", |
| 94 | + ")\n", |
| 95 | + "\n", |
| 96 | + "sim.run(100)" |
| 97 | + ] |
| 98 | + }, |
| 99 | + { |
| 100 | + "cell_type": "code", |
| 101 | + "execution_count": null, |
| 102 | + "metadata": {}, |
| 103 | + "outputs": [], |
| 104 | + "source": [ |
| 105 | + "time, Z_pr = scp_pr.read()\n", |
| 106 | + "_, Z_rks = scp_rks.read()\n", |
| 107 | + "P_vals = 10**(4 + time * 0.035) / 1e6 # MPa\n", |
104 | 108 | "\n", |
105 | 109 | "fig, ax = plt.subplots(figsize=(7, 5))\n", |
106 | | - "ax.semilogx(P_range / 1e6, Z_pr_arr, label=\"Peng-Robinson\")\n", |
107 | | - "ax.semilogx(P_range / 1e6, Z_rks_arr, \"--\", label=\"Soave-Redlich-Kwong\")\n", |
| 110 | + "ax.semilogx(P_vals, Z_pr[0], label=\"Peng-Robinson\")\n", |
| 111 | + "ax.semilogx(P_vals, Z_rks[0], \"--\", label=\"Soave-Redlich-Kwong\")\n", |
108 | 112 | "ax.axhline(1.0, color=\"gray\", linestyle=\"-.\", alpha=0.5, label=\"Ideal gas\")\n", |
109 | 113 | "ax.set_xlabel(\"Pressure [MPa]\")\n", |
110 | 114 | "ax.set_ylabel(\"Compressibility Factor Z\")\n", |
111 | | - "ax.set_title(f\"Methane at T = {T_fixed} K\")\n", |
| 115 | + "ax.set_title(\"Methane at T = 250 K\")\n", |
112 | 116 | "ax.legend()\n", |
113 | 117 | "ax.grid(True, alpha=0.3)\n", |
114 | 118 | "plt.tight_layout()\n", |
|
126 | 130 | "cell_type": "markdown", |
127 | 131 | "metadata": {}, |
128 | 132 | "source": [ |
129 | | - "## Z vs Reduced Temperature\n", |
| 133 | + "## Mixture\n", |
130 | 134 | "\n", |
131 | | - "Vary the temperature at a fixed pressure to see how the compressibility factor changes with reduced temperature $T_r = T / T_c$." |
| 135 | + "The EoS blocks also support mixtures through van der Waals one-fluid mixing rules. Here we set up a methane-ethane mixture and sweep pressure." |
132 | 136 | ] |
133 | 137 | }, |
134 | 138 | { |
|
137 | 141 | "metadata": {}, |
138 | 142 | "outputs": [], |
139 | 143 | "source": [ |
140 | | - "P_fixed = 5e6 # 5 MPa\n", |
141 | | - "T_range = np.linspace(200, 600, 100)\n", |
142 | | - "Tr_range = T_range / Tc_CH4\n", |
143 | | - "\n", |
144 | | - "Z_vs_T = []\n", |
145 | | - "for T in T_range:\n", |
146 | | - " _, Z = eval_eos(pr, T, P_fixed)\n", |
147 | | - " Z_vs_T.append(Z)\n", |
148 | | - "\n", |
149 | | - "fig, ax = plt.subplots(figsize=(7, 5))\n", |
150 | | - "ax.plot(Tr_range, Z_vs_T)\n", |
151 | | - "ax.axhline(1.0, color=\"gray\", linestyle=\"-.\", alpha=0.5)\n", |
152 | | - "ax.axvline(1.0, color=\"red\", linestyle=\":\", alpha=0.5, label=r\"$T_r = 1$ (critical)\")\n", |
153 | | - "ax.set_xlabel(r\"Reduced Temperature $T_r = T / T_c$\")\n", |
154 | | - "ax.set_ylabel(\"Compressibility Factor Z\")\n", |
155 | | - "ax.set_title(f\"Peng-Robinson: Methane at P = {P_fixed/1e6:.0f} MPa\")\n", |
156 | | - "ax.legend()\n", |
157 | | - "ax.grid(True, alpha=0.3)\n", |
158 | | - "plt.tight_layout()\n", |
159 | | - "plt.show()" |
160 | | - ] |
161 | | - }, |
162 | | - { |
163 | | - "cell_type": "markdown", |
164 | | - "metadata": {}, |
165 | | - "source": [ |
166 | | - "## Methane-Ethane Mixture\n", |
167 | | - "\n", |
168 | | - "The EoS blocks support mixtures through standard van der Waals one-fluid mixing rules. Supply arrays of critical properties and mole fractions." |
| 144 | + "pr_mix = PengRobinson(\n", |
| 145 | + " Tc=[190.6, 305.3],\n", |
| 146 | + " Pc=[4.6e6, 4.872e6],\n", |
| 147 | + " omega=[0.011, 0.099],\n", |
| 148 | + " x=[0.7, 0.3],\n", |
| 149 | + ")\n", |
| 150 | + "\n", |
| 151 | + "T_const2 = Constant(300)\n", |
| 152 | + "P_src2 = Source(func=lambda t: 10**(4 + t * 0.035))\n", |
| 153 | + "scp_mix = Scope(labels=[\"Z_mixture\"])\n", |
| 154 | + "\n", |
| 155 | + "sim_mix = Simulation(\n", |
| 156 | + " blocks=[T_const2, P_src2, pr_mix, scp_mix],\n", |
| 157 | + " connections=[\n", |
| 158 | + " Connection(T_const2, pr_mix),\n", |
| 159 | + " Connection(P_src2, pr_mix[1]),\n", |
| 160 | + " Connection(pr_mix[1], scp_mix),\n", |
| 161 | + " ],\n", |
| 162 | + " dt=1.0,\n", |
| 163 | + ")\n", |
| 164 | + "\n", |
| 165 | + "sim_mix.run(100)" |
169 | 166 | ] |
170 | 167 | }, |
171 | 168 | { |
|
174 | 171 | "metadata": {}, |
175 | 172 | "outputs": [], |
176 | 173 | "source": [ |
177 | | - "# Critical properties: methane, ethane\n", |
178 | | - "Tc = [190.6, 305.3] # K\n", |
179 | | - "Pc = [4.6e6, 4.872e6] # Pa\n", |
180 | | - "omega = [0.011, 0.099] # acentric factors\n", |
181 | | - "\n", |
182 | | - "# Compare Z at different compositions\n", |
183 | | - "x1_range = np.linspace(0, 1, 20) # methane mole fraction\n", |
184 | | - "T_mix, P_mix = 300, 3e6 # 300 K, 3 MPa\n", |
185 | | - "\n", |
186 | | - "Z_mix = []\n", |
187 | | - "for x1 in x1_range:\n", |
188 | | - " pr_mix = PengRobinson(\n", |
189 | | - " Tc=Tc, Pc=Pc, omega=omega,\n", |
190 | | - " x=[x1, 1 - x1],\n", |
191 | | - " )\n", |
192 | | - " _, Z = eval_eos(pr_mix, T_mix, P_mix)\n", |
193 | | - " Z_mix.append(Z)\n", |
| 174 | + "time_m, Z_mix = scp_mix.read()\n", |
| 175 | + "P_mix = 10**(4 + time_m * 0.035) / 1e6\n", |
194 | 176 | "\n", |
195 | 177 | "fig, ax = plt.subplots(figsize=(7, 5))\n", |
196 | | - "ax.plot(x1_range, Z_mix, \"o-\")\n", |
197 | | - "ax.set_xlabel(r\"$x_{\\mathrm{CH_4}}$\")\n", |
| 178 | + "ax.semilogx(P_vals, Z_pr[0], label=\"Pure CH₄ (250 K)\")\n", |
| 179 | + "ax.semilogx(P_mix, Z_mix[0], \"--\", label=\"70/30 CH₄-C₂H₆ (300 K)\")\n", |
| 180 | + "ax.axhline(1.0, color=\"gray\", linestyle=\"-.\", alpha=0.5)\n", |
| 181 | + "ax.set_xlabel(\"Pressure [MPa]\")\n", |
198 | 182 | "ax.set_ylabel(\"Compressibility Factor Z\")\n", |
199 | | - "ax.set_title(f\"PR: Methane-Ethane at T = {T_mix} K, P = {P_mix/1e6:.0f} MPa\")\n", |
| 183 | + "ax.set_title(\"Peng-Robinson: Pure vs Mixture\")\n", |
| 184 | + "ax.legend()\n", |
200 | 185 | "ax.grid(True, alpha=0.3)\n", |
201 | 186 | "plt.tight_layout()\n", |
202 | 187 | "plt.show()" |
|
206 | 191 | "cell_type": "markdown", |
207 | 192 | "metadata": {}, |
208 | 193 | "source": [ |
209 | | - "The mixture compressibility factor varies smoothly with composition. Pure ethane (right side) has a lower $Z$ at these conditions because it is closer to its critical point ($T_c = 305.3$ K) and therefore deviates more from ideal gas behavior." |
| 194 | + "The mixture shows a deeper dip because ethane ($T_c = 305.3$ K) is near its critical temperature at 300 K, leading to stronger non-ideal behavior." |
210 | 195 | ] |
211 | 196 | } |
212 | 197 | ], |
|
0 commit comments