Skip to content

Commit

Permalink
HSG
Browse files Browse the repository at this point in the history
  • Loading branch information
@SalvadorBrandolin
SalvadorBrandolin committed Feb 5, 2025
1 parent 0b904cc commit 88a87aa
Show file tree
Hide file tree
Showing 2 changed files with 250 additions and 37 deletions.
174 changes: 174 additions & 0 deletions doc/page/usage/eos/index.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -316,4 +316,178 @@ fugacity_vt: block
)
end block fugacity_vt
```

## Residual entropy (V,T): \(S^r(n,V,T)\)

We start explaining the residual entropy because it is the easiest property to
calculate. Also, helps us to understand how to calculate the next properties.

The residual entropy is calculated as follows:

$$
S^r = - \left(\frac{\partial A^r}{\partial T} \right)_{V,n}
$$

And its derivatives:

$$
\left(\frac{\partial S^r}{\partial T} \right)_{V,n} =
- \left(\frac{\partial^2 A^r}{\partial T^2} \right)_{V,n}
$$


$$
\left(\frac{\partial S^r}{\partial V} \right)_{T,n} =
- \left(\frac{\partial^2 A^r}{\partial V \partial T} \right)_n
$$

$$
\left(\frac{\partial S^r}{\partial n_i} \right)_{V,T} =
- \left(\frac{\partial^2 A^r}{\partial n_i \partial T} \right)_V
$$

```fortran
residual_entropy: block
real(pr) :: T, V, S, dSdV, dSdT, dSdn(2)
T = 300.0_pr ! Set temperature to 300 K
V = 1.0_pr ! Set volume to 1 liter
! Calculate residual entropy and its derivatives
call eos%residual_entropy(n, V, T, Sr=Sr, SrT=SrT, SrV=SrV, Srn=Srn)
print *, "Residual entropy: ", Sr
print *, "SrV: ", SrV
print *, "SrT: ", SrT
print *, "Srn: ", Srn
end block residual_entropy
```

## Residual enthalpy (V,T): \(H^r(n,V,T)\)

The residual enthalpy is calculated as follows:

$$
H^r(n,V,T) = H^r(n,P,T) = A^r(n,V,T) + T S^r(n,V,T) + PV - nRT
$$

We know that pressure can be calculated as:

$$
P = - \left(\frac{\partial A^r}{\partial V} \right)_{T,n} + \frac{nRT}{V}
$$

Then, we can obtain:

$$
P - \frac{nRT}{V} = - \left(\frac{\partial A^r}{\partial V} \right)_{T,n}
$$

$$
PV - nRT = - V \left(\frac{\partial A^r}{\partial V} \right)_{T,n}
$$

And we also know how to calculate residual entropy as:

$$
S^r = - \left(\frac{\partial A^r}{\partial T} \right)_{V,n}
$$

Then, we can obtain the residual enthalpy as:

$$
H^r = A^r - T \left(\frac{\partial A^r}{\partial T} \right)_{V,n}
- V \left(\frac{\partial A^r}{\partial V} \right)_{T,n}
$$

And its derivatives:

$$
\left(\frac{\partial H^r}{\partial T} \right)_{V,n} =
- T \left(\frac{\partial^2 A^r}{\partial T^2} \right)_{V,n}
- V \left(\frac{\partial^2 A^r}{\partial V \partial T} \right)_n
$$

$$
\left(\frac{\partial H^r}{\partial V} \right)_{T,n} =
- T \left(\frac{\partial^2 A^r}{\partial V \partial T} \right)_n
- V \left(\frac{\partial^2 A^r}{\partial V^2} \right)_{T,n}
$$

$$
\left(\frac{\partial H^r}{\partial n_i} \right)_{V,T} =
\left(\frac{\partial A^r}{\partial n_i} \right)_{V,T}
- T \left(\frac{\partial^2 A^r}{\partial T \partial n_i} \right)_V
- V \left(\frac{\partial^2 A^r}{\partial V \partial n_i} \right)_T
$$

```fortran
residual_enthalpy: block
real(pr) :: T, V, Hr, HrV, HrT, Hrn(2)
T = 300.0_pr ! Set temperature to 300 K
V = 1.0_pr ! Set volume to 1 liter
! Calculate residual enthalpy and its derivatives
call eos%residual_enthalpy(n, V, T, Hr=Hr, HrV=HrV, HrT=HrT, Hrn=Hrn)
print *, "Residual enthalpy: ", Hr
print *, "HrV: ", HrV
print *, "HrT: ", HrT
print *, "Hrn: ", Hrn
end block residual_enthalpy
```

## Residual Gibbs free energy (V,T): \(G^r(n,V,T)\)

The residual Gibbs free energy is calculated as follows:

$$
G^r(n,V,T) = A^r(n,V,T) + PV - nRT
$$

As with residual enthalpy we can easily deduce:

$$
G^r = A^r - V \left(\frac{\partial A^r}{\partial V} \right)_{T,n}
$$

And its derivatives:

$$
\left(\frac{\partial G^r}{\partial T} \right)_{V,n} =
\left(\frac{\partial A^r}{\partial T} \right)_{V,n}
- V \left(\frac{\partial^2 A^r}{\partial T \partial V} \right)_n
$$

$$
\left(\frac{\partial G^r}{\partial V} \right)_{T,n} =
- V \left(\frac{\partial^2 A^r}{\partial V^2} \right)_{T,n}
$$

$$
\left(\frac{\partial G^r}{\partial n_i} \right)_{V,T} =
\left(\frac{\partial A^r}{\partial n_i} \right)_{V,T}
- V \left(\frac{\partial^2 A^r}{\partial V \partial n_i} \right)_T
$$

```fortran
residual_gibbs: block
real(pr) :: T, V, Gr, GrV, GrT, Grn(2)
T = 300.0_pr ! Set temperature to 300 K
V = 1.0_pr ! Set volume to 1 liter
! Calculate residual Gibbs free energy and its derivatives
call eos%residual_gibbs(n, V, T, Gr=Gr, GrV=GrV, GrT=GrT, Grn=Grn)
print *, "Residual Gibbs free energy: ", Gr
print *, "GrV: ", GrV
print *, "GrT: ", GrT
print *, "Grn: ", Grn
end block residual_gibbs
```
113 changes: 76 additions & 37 deletions src/models/residual_helmholtz/ar_models.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -265,9 +265,9 @@ subroutine lnphi_pt(eos, &
!! Calculate natural logarithm of fugacity given pressure and temperature.
!!
!! Calculate the natural logarithm of the fugacity coefficient and its
!! derivatives given pressure and temperature. This routine will obtain
!! the desired volume root at the specified pressure and calculate
!! fugacity at that point.The routine gives the possibility to calculate
!! derivatives given pressure and temperature. This routine will obtain
!! the desired volume root at the specified pressure and calculate
!! fugacity at that point.The routine gives the possibility to calculate
!! the pressure derivatives and volume.
!!
!! # Examples
Expand Down Expand Up @@ -410,11 +410,11 @@ subroutine lnphi_vt(eos, &
dPdn_in = RT/V - ArVn

if (present(lnPhi)) lnPhi = Arn(:)/RT - log(Z)

if (present(dlnPhidP)) then
dlnPhidP(:) = -dPdn_in(:)/dPdV_in/RT - 1._pr/P_in
end if

if (present(dlnPhidT)) then
dlnPhidT(:) = (ArTn(:) - Arn(:)/T)/RT + dPdn_in(:)*dPdT_in/dPdV_in/RT + 1._pr/T
end if
Expand Down Expand Up @@ -558,68 +558,108 @@ subroutine lnfug_vt(eos, &
if (present(dPdn)) dPdn = dPdn_in
end subroutine lnfug_vt

subroutine enthalpy_residual_vt(eos, n, V, T, Hr, HrT, HrV, Hrn)
subroutine enthalpy_residual_vt(eos, n, V, T, Hr, HrV, HrT, Hrn)
!! Calculate residual enthalpy given volume and temperature.
!!
!! # Examples
!!
!! ```fortran
!! eos = PengRobinson76(Tc, Pc, w)
!!
!! n = [1.0_pr, 1.0_pr]
!! T = 300.0_pr
!! V = 1.0_pr
!!
!! call eos%enthalpy_residual_vt(&
!! n, V, T, Hr=Hr, HrV=HrV, HrT=HrT, Hrn=Hrn &
!! )
!! ```
class(ArModel), intent(in) :: eos !! Model
real(pr), intent(in) :: n(:) !! Moles number vector
real(pr), intent(in) :: T !! Temperature [K]
real(pr), intent(in) :: V !! Volume [L]
real(pr), intent(out) :: Hr !! Residual enthalpy [bar L]
real(pr), optional, intent(out) :: HrT !! \(\frac{dH^r}}{dT}\)
real(pr), optional, intent(out) :: Hr !! Residual enthalpy [bar L]
real(pr), optional, intent(out) :: HrV !! \(\frac{dH^r}}{dV}\)
real(pr), optional, intent(out) :: HrT !! \(\frac{dH^r}}{dT}\)
real(pr), optional, intent(out) :: Hrn(size(n)) !! \(\frac{dH^r}}{dn}\)

real(pr) :: Ar, ArV, ArT, Arn(size(n))
real(pr) :: ArV2, ArT2, ArTV, ArVn(size(n)), ArTn(size(n))

call eos%residual_helmholtz(&
n, V, T, Ar=Ar, ArV=ArV, ArT=ArT, ArTV=ArTV, ArV2=ArV2, ArT2=ArT2, Arn=Arn, ArVn=ArVn, ArTn=ArTn &
n, V, T, Ar=Ar, ArV=ArV, ArT=ArT, ArTV=ArTV, &
ArV2=ArV2, ArT2=ArT2, Arn=Arn, ArVn=ArVn, ArTn=ArTn &
)

Hr = Ar - T*ArT - V*ArV

if (present(HrT)) HrT = - t*ArT2 - v*ArTV
if (present(HrV)) HrV = - t*ArTV - v*ArV2
if (present(HrN)) HrN(:) = Arn(:) - t*ArTn(:) - v*ArVn(:)
if (present(Hr)) Hr = Ar - T*ArT - V*ArV
if (present(HrT)) HrT = - T*ArT2 - V*ArTV
if (present(HrV)) HrV = - T*ArTV - V*ArV2
if (present(Hrn)) Hrn(:) = Arn(:) - T*ArTn(:) - V*ArVn(:)
end subroutine enthalpy_residual_vt

subroutine gibbs_residual_VT(eos, n, V, T, Gr, GrT, GrV, Grn)
subroutine gibbs_residual_vt(eos, n, V, T, Gr, GrV, GrT, Grn)
!! Calculate residual Gibbs energy given volume and temperature.
!!
!! # Examples
!!
!! ```fortran
!! eos = PengRobinson76(Tc, Pc, w)
!!
!! n = [1.0_pr, 1.0_pr]
!! T = 300.0_pr
!! V = 1.0_pr
!!
!! call eos%gibbs_residual_vt(&
!! n, V, T, Gr=Gr, GrV=GrV, GrT=GrT, Grn=Grn &
!! )
!! ```
use yaeos__constants, only: R
class(ArModel), intent(in) :: eos !! Model
real(pr), intent(in) :: n(:) !! Moles number vector
real(pr), intent(in) :: V !! Volume [L]
real(pr), intent(in) :: T !! Temperature [K]
real(pr), intent(out) :: Gr !! Gibbs energy [bar L]
real(pr), optional, intent(out) :: GrT !! \(\frac{dG^r}}{dT}\)
real(pr), optional, intent(out) :: Gr !! Gibbs energy [bar L]
real(pr), optional, intent(out) :: GrV !! \(\frac{dG^r}}{dV}\)
real(pr), optional, intent(out) :: GrT !! \(\frac{dG^r}}{dT}\)
real(pr), optional, intent(out) :: Grn(size(n)) !! \(\frac{dG^r}}{dn}\)

real(pr) :: Ar, ArV, ArT, Arn(size(n))
real(pr) :: Ar, ArV, ArT, Arn(size(n)), ArTV, ArV2, ArVn(size(n))
real(pr) :: p, dPdV, dPdT, dPdn(size(n)), z, totn

totn = sum(n)
call pressure(eos, n, V, T, P, dPdV=dPdV, dPdT=dPdT, dPdn=dPdn)
z = P*V/(totn*R*T)

call eos%residual_helmholtz(n, v, t, Ar=Ar, ArV=ArV, ArT=ArT, Arn=Arn)

Gr = Ar + P*V - totn*R*T
call eos%residual_helmholtz(&
n, v, t, Ar=Ar, ArV=ArV, &
ArT=ArT, Arn=Arn, ArTV=ArTV, ArV2=ArV2, ArVn=ArVn &
)

if (present(GrT)) GrT = ArT + V*dPdT - totn*R
if (present(GrV)) GrV = ArV + V*dPdV + P
if (present(GrN)) GrN(:) = Arn(:) + V*dPdn(:) - R*T
end subroutine gibbs_residual_VT
if (present(Gr)) Gr = Ar - V*ArV
if (present(GrT)) GrT = ArT - V*ArTV
if (present(GrV)) GrV = -V*ArV2
if (present(Grn)) Grn(:) = Arn(:) - V*ArVn(:)
end subroutine gibbs_residual_vt

subroutine entropy_residual_vt(eos, n, V, T, Sr, SrT, SrV, Srn)
subroutine entropy_residual_vt(eos, n, V, T, Sr, SrV, SrT, Srn)
!! Calculate residual entropy given volume and temperature.
!!
!! # Examples
!!
!! ```fortran
!! eos = PengRobinson76(Tc, Pc, w)
!!
!! n = [1.0_pr, 1.0_pr]
!! T = 300.0_pr
!! V = 1.0_pr
!!
!! call eos%entropy_residual_vt(&
!! n, V, T, Sr=Sr, SrV=SrV, SrT=SrT, Srn=Srn &
!! )
!! ```
class(ArModel), intent(in) :: eos !! Model
real(pr), intent(in) :: n(:) !! Moles number vector
real(pr), intent(in) :: V !! Volume [L]
real(pr), intent(in) :: T !! Temperature [K]
real(pr), intent(out) :: Sr !! Entropy [bar L / K]
real(pr), optional, intent(out) :: SrT !! \(\frac{dS^r}}{dT}\)
real(pr), optional, intent(out) :: Sr !! Entropy [bar L / K]
real(pr), optional, intent(out) :: SrV !! \(\frac{dS^r}}{dV}\)
real(pr), optional, intent(out) :: SrT !! \(\frac{dS^r}}{dT}\)
real(pr), optional, intent(out) :: Srn(size(n)) !! \(\frac{dS^r}}{dn}\)

real(pr) :: Ar, ArT, ArT2, ArTV, ArTn(size(n))
Expand All @@ -628,11 +668,10 @@ subroutine entropy_residual_vt(eos, n, V, T, Sr, SrT, SrV, Srn)
n, v, t, Ar=Ar, ArT=ArT, ArTV=ArTV, ArT2=ArT2, ArTn=ArTn &
)

Sr = - ArT

if (present(SrT)) SrT = - ArT2
if (present(SrV)) SrV = - ArTV
if (present(SrN)) SrN = - ArTn
if (present(Sr)) Sr = -ArT
if (present(SrT)) SrT = -ArT2
if (present(SrV)) SrV = -ArTV
if (present(SrN)) Srn = -ArTn
end subroutine entropy_residual_vt

subroutine Cv_residual_vt(eos, n, V, T, Cv)
Expand Down

0 comments on commit 88a87aa

Please sign in to comment.