From 8a421c3cd7155c8ed68082862ab0fc61218b888f Mon Sep 17 00:00:00 2001
From: wengzf20 <67515007+wengzf20@users.noreply.github.com>
Date: Tue, 27 Jun 2023 14:13:48 +0800
Subject: [PATCH] Update compressibility calculation for real fluid in
 CanteraMixture.H

the compressibility for constant enthalpy was used
---
 src/dfCanteraMixture/CanteraMixture.H | 20 ++++++++++++++++++--
 1 file changed, 18 insertions(+), 2 deletions(-)

diff --git a/src/dfCanteraMixture/CanteraMixture.H b/src/dfCanteraMixture/CanteraMixture.H
index 80443318a..87e4b6165 100644
--- a/src/dfCanteraMixture/CanteraMixture.H
+++ b/src/dfCanteraMixture/CanteraMixture.H
@@ -119,8 +119,24 @@ public:
             const scalar& T
         ) const
         {
-            // for ideal gas, psi = CanteraGas_->density()/p = 1/RT
-            return CanteraGas_->density()/p;
+            const Cantera::AnyMap& map = CanteraGas_->input();
+            if(map["thermo"].asString()=="Peng-Robinson" || map["thermo"].asString()=="Redlich-Kwong")
+            {
+                // for non-ideal EoS, psi = \partial rho/ \partial p at constant enthalpy
+                // Currently, Peng-Robinson and Redlich-Kwong are available in Cantera
+                // The partial derivative was approximated using backward difference with a step size dx
+                double h0 = CanteraGas_->enthalpy_mass();
+                double dx = 0.0001;
+                double rho1 = CanteraGas_->density();
+                CanteraGas_->setState_HP(h0,p*(1-dx));
+                double rho2 = CanteraGas_->density();
+                return (rho1-rho2)/(dx*p);
+            }
+            else
+            {
+                // for ideal gas, psi = W/RT
+                return CanteraGas_->meanMolecularWeight()/CanteraGas_->RT();
+            }
         }
 
         scalar rho