feos-org · prehner · Apr 14, 2025 · Apr 14, 2025
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -6,6 +6,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
+## [0.2.1] - 2025-04-14
+### Added
+- Added `StateAD::molar_isochoric_heat_capacity` and `StateAD::molar_isobaric_heat_capacity`. [#3](https://github.com/feos-org/feos-ad/pull/3)
+
 ## [0.2.0] - 2025-01-27
 ### Changed
 - Made `PcSaftBinary` generic for associating/non-associating systems. [#2](https://github.com/feos-org/feos-ad/pull/2)

diff --git a/Cargo.toml b/Cargo.toml
@@ -1,6 +1,6 @@
 [package]
 name = "feos-ad"
-version = "0.2.0"
+version = "0.2.1"
 authors = ["Philipp Rehner <prehner@ethz.ch"]
 edition = "2021"
 readme = "README.md"

diff --git a/src/core/state.rs b/src/core/state.rs
@@ -174,6 +174,24 @@ impl<'a, E: TotalHelmholtzEnergy<N>, D: DualNum<f64> + Copy, const N: usize> Sta
             &self.molefracs,
         ))
     }
+
+    pub fn molar_isochoric_heat_capacity(&self) -> MolarEntropy<D> {
+        MolarEntropy::from_reduced(E::molar_isochoric_heat_capacity(
+            &self.eos.parameters,
+            self.reduced_temperature,
+            self.reduced_molar_volume,
+            &self.molefracs,
+        ))
+    }
+
+    pub fn molar_isobaric_heat_capacity(&self) -> MolarEntropy<D> {
+        MolarEntropy::from_reduced(E::molar_isobaric_heat_capacity(
+            &self.eos.parameters,
+            self.reduced_temperature,
+            self.reduced_molar_volume,
+            &self.molefracs,
+        ))
+    }
 }
 
 impl<'a, E: ResidualHelmholtzEnergy<1>, D: DualNum<f64> + Copy> StateAD<'a, E, D, 1> {
@@ -569,4 +587,45 @@ mod test {
         );
         Ok(())
     }
+
+    #[test]
+    fn test_heat_capacities() -> EosResult<()> {
+        let (joback, ideal_gas) = joback()?;
+        let (pcsaft, residual) = pcsaft()?;
+        let eos = Arc::new(EquationOfState::new(ideal_gas, residual));
+        let eos_ad = EquationOfStateAD::new([joback], pcsaft)
+            .wrap()
+            .derivatives(([joback.params()], pcsaft.params()));
+
+        let temperature = 300.0 * KELVIN;
+        let pressure = 5.0 * BAR;
+
+        let state = State::new_npt(
+            &eos,
+            temperature,
+            pressure,
+            &(arr1(&[1.0]) * MOL),
+            DensityInitialization::None,
+        )?;
+        let state_ad = StateAD::new_tp(
+            &eos_ad,
+            temperature,
+            pressure,
+            SVector::from([1.0]),
+            DensityInitialization::None,
+        )?;
+
+        let c_v = state.molar_isochoric_heat_capacity(Contributions::Total);
+        let c_p = state.molar_isobaric_heat_capacity(Contributions::Total);
+        let c_v_ad = state_ad.molar_isochoric_heat_capacity();
+        let c_p_ad = state_ad.molar_isobaric_heat_capacity();
+
+        println!("{c_v} {c_p}");
+        println!("{c_v_ad} {c_p_ad}");
+
+        assert_relative_eq!(c_v, c_v_ad, max_relative = 1e-10);
+        assert_relative_eq!(c_p, c_p_ad, max_relative = 1e-10);
+
+        Ok(())
+    }
 }
diff --git a/src/core/total.rs b/src/core/total.rs
@@ -1,6 +1,8 @@
 use super::{ParametersAD, ResidualHelmholtzEnergy};
 use nalgebra::SVector;
-use num_dual::{first_derivative, gradient, Dual, DualNum, DualVec};
+use num_dual::{
+    first_derivative, gradient, hessian, second_derivative, Dual, Dual2, Dual2Vec, DualNum, DualVec,
+};
 
 /// Implementation of an ideal gas Helmholtz energy contribution.
 pub trait IdealGasAD: ParametersAD {
@@ -156,6 +158,43 @@ pub trait TotalHelmholtzEnergy<const N: usize>: ResidualHelmholtzEnergy<N> {
         a - temperature * da_dt - molar_volume * da_dv
     }
 
+    fn molar_isochoric_heat_capacity<D: DualNum<f64> + Copy>(
+        parameters: &Self::Parameters<D>,
+        temperature: D,
+        molar_volume: D,
+        molefracs: &SVector<D, N>,
+    ) -> D {
+        let params = Self::params_from_inner(parameters);
+        let molar_volume = Dual2::from_re(molar_volume);
+        let molefracs = molefracs.map(Dual2::from_re);
+        let (_, _, d2a) = second_derivative(
+            |temperature| {
+                Self::molar_helmholtz_energy(&params, temperature, molar_volume, &molefracs)
+            },
+            temperature,
+        );
+        -temperature * d2a
+    }
+
+    fn molar_isobaric_heat_capacity<D: DualNum<f64> + Copy>(
+        parameters: &Self::Parameters<D>,
+        temperature: D,
+        molar_volume: D,
+        molefracs: &SVector<D, N>,
+    ) -> D {
+        let params = Self::params_from_inner(parameters);
+        let molefracs = molefracs.map(Dual2Vec::from_re);
+        let (_, _, d2a) = hessian(
+            |x| {
+                let [temperature, molar_volume] = x.data.0[0];
+                Self::molar_helmholtz_energy(&params, temperature, molar_volume, &molefracs)
+            },
+            SVector::from([temperature, molar_volume]),
+        );
+        let [[a_tt, a_tv], [_, a_vv]] = d2a.data.0;
+        temperature * (a_tv * a_tv / a_vv - a_tt)
+    }
+
     fn pressure_entropy<D: DualNum<f64> + Copy>(
         parameters: &Self::Parameters<D>,
         temperature: D,