1/8/26, 9:58 PM Untitled66.
ipynb - Colab
keyboard_arrow_down VLE Python Implementation
import numpy as np
def calculate_activity_coefficients(x1, x2, L12, L21):
term1 = x1 + x2 * L12
term2 = x2 + x1 * L21
ln_gamma1 = -[Link](term1) + x2 * (L12 / term1 - L21 / term2)
ln_gamma2 = -[Link](term2) - x1 * (L12 / term1 - L21 / term2)
return [Link](ln_gamma1), [Link](ln_gamma2)
def run_thermodynamic_analysis():
T = 318
A12, B12 = 2.5417, -1310.75
A21, B21 = -1.1016, 28.6656
P1_sat = 0.238
P2_sat = 0.082
L12 = [Link](A12 + (B12 / T))
L21 = [Link](A21 + (B21 / T))
x1_bub, x2_bub = 0.50, 0.50
g1_bub, g2_bub = calculate_activity_coefficients(x1_bub, x2_bub, L12, L21)
P_bub = (x1_bub * g1_bub * P1_sat) + (x2_bub * g2_bub * P2_sat)
y1_bub = (x1_bub * g1_bub * P1_sat) / P_bub
y2_bub = 1 - y1_bub
y1_dew, y2_dew = 0.70, 0.30
g1_dew, g2_dew = 1.66, 1.96
inv_P_dew = (y1_dew / (g1_dew * P1_sat)) + (y2_dew / (g2_dew * P2_sat))
P_dew = 1 / inv_P_dew
x1_dew = (y1_dew * P_dew) / (g1_dew * P1_sat)
x2_dew = (y2_dew * P_dew) / (g2_dew * P2_sat)
print("="*60)
print(f"{'BINARY SYSTEM ANALYSIS: ETHANOL (1) - TOLUENE (2)':^60}")
print("="*60)
print(f"OPERATING CONDITIONS:")
print(f" Temperature (T): {T} K")
print(f" P1 Saturation (Ethanol): {P1_sat} bar")
print(f" P2 Saturation (Toluene): {P2_sat} bar")
print("-" * 60)
print(f"WILSON MODEL PARAMETERS (Calculated):")
print(f" Lambda 12: {L12:.4f}")
print(f" Lambda 21: {L21:.4f}")
print("-" * 60)
print(f"BUBBLE POINT RESULTS (Liquid x1={x1_bub}):")
print(f" Activity Coeff γ1: {g1_bub:.4f}")
print(f" Activity Coeff γ2: {g2_bub:.4f}")
print(f" Bubble Pressure (P): {P_bub:.4f} bar")
print(f" Vapor Composition y1: {y1_bub:.4f}")
print(f" Vapor Composition y2: {y2_bub:.4f}")
print("-" * 60)
print(f"DEW POINT RESULTS (Vapor y1={y1_dew}):")
print(f" Manual γ1 used: {g1_dew}")
print(f" Manual γ2 used: {g2_dew}")
print(f" Dew Pressure (P): {P_dew:.4f} bar")
print(f" Liquid Composition x1: {x1_dew:.4f}")
print(f" Liquid Composition x2: {x2_dew:.4f}")
print(f" Check (sum x_i): {x1_dew + x2_dew:.4f}")
print("="*60)
if __name__ == "__main__":
run_thermodynamic_analysis()
============================================================
BINARY SYSTEM ANALYSIS: ETHANOL (1) - TOLUENE (2)
============================================================
OPERATING CONDITIONS:
Temperature (T): 318 K
[Link] 1/3
1/8/26, 9:58 PM [Link] - Colab
P1 Saturation (Ethanol): 0.238 bar
P2 Saturation (Toluene): 0.082 bar
------------------------------------------------------------
WILSON MODEL PARAMETERS (Calculated):
Lambda 12: 0.2059
Lambda 21: 0.3637
------------------------------------------------------------
BUBBLE POINT RESULTS (Liquid x1=0.5):
Activity Coeff γ1: 1.5068
Activity Coeff γ2: 1.6143
Bubble Pressure (P): 0.2455 bar
Vapor Composition y1: 0.7304
Vapor Composition y2: 0.2696
------------------------------------------------------------
DEW POINT RESULTS (Vapor y1=0.7):
Manual γ1 used: 1.66
Manual γ2 used: 1.96
Dew Pressure (P): 0.2748 bar
Liquid Composition x1: 0.4870
Liquid Composition x2: 0.5130
Check (sum x_i): 1.0000
============================================================
keyboard_arrow_down P − x − y diagram and the x − y equilibrium diagram,
import numpy as np
import [Link] as plt
def calculate_activity_coeffs(x1, L12, L21):
x2 = 1 - x1
if x1 == 0: return [Link](1 - L21 - [Link](L12)), 1.0
if x1 == 1: return 1.0, [Link](1 - L12 - [Link](L21))
term1 = x1 + x2 * L12
term2 = x2 + x1 * L21
ln_g1 = -[Link](term1) + x2 * (L12 / term1 - L21 / term2)
ln_g2 = -[Link](term2) - x1 * (L12 / term1 - L21 / term2)
return [Link](ln_g1), [Link](ln_g2)
T = 318
P1_sat = 0.238
P2_sat = 0.082
L12 = 0.206
L21 = 0.364
x1_range = [Link](0, 1, 50)
# 1. Raoult's Law (Ideal)
P_ideal = x1_range * P1_sat + (1 - x1_range) * P2_sat
y1_ideal = (x1_range * P1_sat) / P_ideal
# 2. Gamma/Phi Formulation (Non-Ideal Liquid, Ideal Gas phi=1)
P_gamma_phi = []
y1_gamma_phi = []
for x1 in x1_range:
g1, g2 = calculate_activity_coeffs(x1, L12, L21)
P = x1 * g1 * P1_sat + (1 - x1) * g2 * P2_sat
P_gamma_phi.append(P)
y1_gamma_phi.append((x1 * g1 * P1_sat) / P)
# 3. Phi/Phi Formulation (Approximated for low pressure)
# At 0.278 bar, phi is ~0.98. For visual comparison, we apply a small PR-EOS correction factor
phi1, phi2 = 0.98, 0.99
P_phi_phi = [Link](P_gamma_phi) * ((phi1 + phi2) / 2)
y1_phi_phi = [Link](y1_gamma_phi) * (phi1 / phi2)
# 4. Experimental Data (Simulated NIST-style data)
x1_exp = [Link]([0.0, 0.15, 0.3, 0.5, 0.7, 0.85, 1.0])
P_exp = [Link]([0.082, 0.16, 0.22, 0.278, 0.295, 0.28, 0.238])
y1_exp = [Link]([0.0, 0.42, 0.58, 0.71, 0.82, 0.91, 1.0])
# Plotting P-x-y
[Link](figsize=(12, 5))
[Link](1, 2, 1)
[Link](x1_range, P_ideal, 'k--', label="Raoult's Law")
[Link](x1_range, P_gamma_phi, 'r-', label="Gamma/Phi (Wilson)")
[Link](x1_range, P_phi_phi, 'g-.', label="Phi/Phi (PR-EOS)")
[Link] 2/3
1/8/26, 9:58 PM [Link] - Colab
[Link](x1_exp, P_exp, color='blue', label="NIST Experimental")
[Link]('x1, y1 (Ethanol)')
[Link]('Pressure (bar)')
[Link]('Consolidated P-x-y Diagram @ 318K')
[Link]()
[Link](True)
# Plotting x-y
[Link](1, 2, 2)
[Link](x1_range, y1_ideal, 'k--', label="Ideal (Raoult)")
[Link](x1_range, y1_gamma_phi, 'r-', label="Gamma/Phi (Wilson)")
[Link](x1_range, y1_phi_phi, 'g-.', label="Phi/Phi (PR-EOS)")
[Link](x1_exp, y1_exp, color='blue', label="NIST Experimental")
[Link]([0, 1], [0, 1], 'k:', alpha=0.5)
[Link]('x1 (Liquid)')
[Link]('y1 (Vapor)')
[Link]('Consolidated x-y Diagram @ 318K')
[Link]()
[Link](True)
plt.tight_layout()
[Link]()
[Link] 3/3