Merge pull request #96 from ipqa-research/dev

Set of examples

example/README.md

@@ -0,0 +1,14 @@
+# yaeos examples
+Here are located a set of examples of how the library could be used.
+
+We show them in three different categories:
+
+- basics: The basics of what can be done with the library
+    1. Thermodynamic properties
+    2. Saturation points
+    3. Phase split (FlashPT and FlashVT)
+    4. Phase envelopes tracing
+    5. Implementation of a simple algorithm to calculate pure component
+      saturation pressure.
+- Advanced: How a more advanced user can take benefit of the library.
+- extra: Other set of examples with no particular structure.

example/tutorials/new_alpha_function.f90 → example/advanced/new_alpha_function.f90

@@ -44,18 +44,21 @@ end subroutine alpha
 end module alpha_mathias_copeman
 
 program new_alpha_example
-   use yaeos__example_tools, only: methane_butane_pr76
    use yaeos, only: pr, PengRobinson76, CubicEoS, QMR
    use alpha_mathias_copeman, only: MathiasCopeman
    type(CubicEoS) :: eos
    type(QMR) :: mr
    type(MathiasCopeman) :: alpha
 
-   real(pr) :: n(2), v, t, P
-   integer :: i
+   real(pr) :: n(2), v, t, P, tc(2), pc(2), w(2)
+
+   ! Methane/ Butane mixture
+   tc = [190.564, 425.12]     ! Critical temperatures
+   pc = [45.99, 37.96]        ! Critical pressures
+   w = [0.0115478, 0.200164]  ! Acentric factors
 
    ! Get the example PR76 binary model
-   eos = methane_butane_pr76()
+   eos = PengRobinson76(tc, pc, w)
 
    ! Define the new alpha function parameters
    alpha%c1 = [0.49258, 0.84209]
@@ -75,6 +78,4 @@ program new_alpha_example
 
    call eos%pressure(n, V, T, P=P)
    print *, "Peng-Robinson76-MC:", P
-
-
 end program new_alpha_example

example/basics/1_basics.f90

@@ -0,0 +1,71 @@
+program basics
+   use yaeos ! First the module is imported (used). This will bring all the
+             ! types and functions defined in the module to this program.
+   implicit none
+   ! ===========================================================================
+
+   ! In yaeos, all residual Helmholtz models are based on the structure
+   !  `ArModel` and are derivatives from it.
+   ! Defining the model variable as an allocatable ArModel means that this
+   ! variable can be redefined as any other kind of variable that derives from
+   ! `ArModel`, for example `CubicEoS`.
+   class(ArModel), allocatable :: model
+
+   integer, parameter :: nc=2 !! In this case we define a constant number of
+                              !! components.
+
+   ! Used input variables
+   real(pr) :: n(nc), Tc(nc), Pc(nc), w(nc)
+
+   ! Output variables
+   real(pr) :: P, dPdV
+   ! ---------------------------------------------------------------------------
+
+   ! Here the executed code starts
+
+   ! A classic Cubic model can be defined with each component critial constants
+   tc = [190._pr,  310._pr]
+   pc = [14._pr,   30._pr ]
+   w  = [0.001_pr, 0.03_pr]
+
+   ! Here we chose to model with the PengRobinson EoS
+   model = PengRobinson78(tc, pc, w)
+
+   ! Number of moles vector
+   n = [0.3, 0.7]
+
+   ! Pressure calculation
+   call model%pressure(n, V=2.5_pr, T=150._pr, P=P)
+   print *, "P: ", P
+
+   ! Derivatives can also be calculated when included as optional arguments!
+   call model%pressure(n, V=2.5_pr, T=150._pr, P=P, dPdV=dPdV)
+   print *, "dPdV: ", dPdV
+
+   thermoproperties: block
+      !! Al the bulk thermodynamic properties currently available
+      real(pr) :: Cpr, Cvr
+      real(pr) :: Gr, GrT, GrV, Grn(nc)
+      real(pr) :: Hr, HrT, HrV, Hrn(nc)
+      real(pr) :: lnPhi(nc), dlnPhidT(nc), dlnPhidP(nc), dlnPhidn(nc, nc)
+
+      real(pr) :: n(nc), P, V, T
+
+      n = [0.3_pr, 0.7_pr]
+      V = 3.5_pr
+      P = 5.0_pr
+      call model%Cp_residual_vt(n, V=2.5_pr, T=150._pr, Cp=Cpr)
+      call model%Cv_residual_vt(n, V=2.5_pr, T=150._pr, Cv=Cvr)
+      call model%gibbs_residual_vt(n, V=2.5_pr, T=150._pr, Gr=Gr, GrT=GrT, GrV=GrV, Grn=Grn)
+      call model%enthalpy_residual_vt(n, V=2.5_pr, T=150._pr, Hr=Hr, HrT=HrT, HrV=HrV, Hrn=Hrn)
+      call model%lnphi_pt(&
+         n, P=P, T=T, root_type="stable", &
+         lnPhi=lnPhi, dlnPhidT=dlnphidT, dlnPhidP=dlnPhidP, dlnPhidn=dlnPhidn &
+      )
+      call model%lnphi_vt(&
+         n, V=V, T=T, &
+         lnPhi=lnPhi, dlnPhidT=dlnphidT, dlnPhidP=dlnPhidP, dlnPhidn=dlnPhidn &
+      )
+      call model%volume(n, P=P, T=T, root_type="vapor", V=V)
+   end block thermoproperties
+end program basics

example/tutorials/saturation_points.f90 → example/basics/2_saturation_points.f90

@@ -7,14 +7,20 @@
 ! All the outputs of this functions use the `EquilibriumState` type
 program saturation
     use yaeos
-    use yaeos__example_tools, only: methane_butane_pr76
 
     class(ArModel), allocatable :: model
     type(EquilibriumState) :: sat_point
 
     real(pr) :: n(2), T
+    real(pr) :: Tc(2), Pc(2), w(2)
 
-    model = methane_butane_pr76()
+    ! Methane/ Butane mixture
+    tc = [190.564, 425.12]     ! Critical temperatures
+    pc = [45.99, 37.96]        ! Critical pressures
+    w = [0.0115478, 0.200164]  ! Acentric factors
+
+    ! Get the example PR76 binary model
+    model = PengRobinson76(tc, pc, w)
 
     n = [2.5, 6.7]
 
@@ -24,7 +30,7 @@ program saturation
     write(*, *) "Bubble pressure:"
     T = 150
     sat_point = saturation_pressure(model, n, T=T, kind="bubble")
-    write (*, *) "kind, T, P: ", sat_point
+    write (*, *) "kind, T, P: ", sat_point%kind, sat_point%T, sat_point%P
     write (*, *) "x: ", sat_point%x
     write (*, *) "y: ", sat_point%y
 
@@ -34,7 +40,7 @@ program saturation
     write(*, *) ""
     write (*, *) "Bubble temperature:"
     sat_point = saturation_temperature(model, n, P=15._pr, kind="bubble")
-    write (*, *) "kind, T, P: ", sat_point
+    write (*, *) "kind, T, P: ", sat_point%kind, sat_point%T, sat_point%P
     write (*, *) "x: ", sat_point%x
     write (*, *) "y: ", sat_point%y
 
@@ -44,7 +50,7 @@ program saturation
     write(*, *) ""
     write(*, *) "Dew pressure:"
     sat_point = saturation_pressure(model, n, T=150._pr, kind="dew")
-    write (*, *) "kind, T, P: ", sat_point
+    write (*, *) "kind, T, P: ", sat_point%kind, sat_point%T, sat_point%P
     write (*, *) "x: ", sat_point%x
     write (*, *) "y: ", sat_point%y
 
@@ -53,8 +59,8 @@ program saturation
     ! --------------------------------------------------------------------------
     write(*, *) ""
     write (*, *) "Dew temperature:"
-    sat_point = saturation_temperature(model, n, P=15._pr, kind="dew")
-    write (*, *) "kind, T, P: ", sat_point
+    sat_point = saturation_temperature(model, n, P=15._pr, kind="dew", t0=330._pr)
+    write (*, *) "kind, T, P: ", sat_point%kind, sat_point%T, sat_point%P
     write (*, *) "x: ", sat_point%x
     write (*, *) "y: ", sat_point%y
 end program

example/basics/3_phase_split.f90

@@ -0,0 +1,46 @@
+program flasher
+    !! Example program to make flash calculations with a PengRobinson76 
+    !! EoS model.
+
+    ! Import the relevant 
+    use yaeos, only: pr, EquilibriumState, flash, PengRobinson76, ArModel, k_wilson
+    implicit none
+
+    ! Variables definition:
+    class(ArModel), allocatable :: model !! Model to use
+    type(EquilibriumState) :: flash_result !! Result of Flash calculation
+
+    real(pr) :: tc(2), pc(2), w(2)
+    real(pr) :: n(2), t, p, k0(2)
+    integer :: iter
+
+    ! Code starts here:
+    print *, "FLASH EXAMPLE:"
+    print *, "Methane/Butane mixture"
+
+    ! Methane/ Butane mixture
+    n = [0.4, 0.6]                      ! Composition
+    tc = [190.564, 425.12]              ! Critical temperatures
+    pc = [45.99, 37.96]                 ! Critical pressures
+    w = [0.0115478, 0.200164]           ! Acentric factors
+
+    ! Use the PengRobinson76 model
+    model = PengRobinson76(tc, pc, w)
+
+    ! Set pressure and temperatures
+    P = 60
+    T = 294
+
+    ! Calculate flashes
+    flash_result = flash(model, n, t=t, p_spec=p, iters=iter)
+    write(*, *) flash_result
+
+    ! K-wilson factors for initialization
+    k0 = k_wilson(model, T=T, P=P)
+    flash_result = flash(model, n, t=t, p_spec=p, k0=k0, iters=iter)
+    write(*, *) flash_result
+
+    ! Specify volume
+    flash_result = flash(model, n, t=t, V_spec=1.0_pr, iters=iter)
+    write(*, *) flash_result
+end program

example/basics/4_phase_envelope.f90

@@ -0,0 +1,54 @@
+program phase_envelope
+   use yaeos
+   implicit none
+
+   class(ArModel), allocatable :: model ! Model
+   type(EquilibriumState) :: bubble ! Saturation point
+   type(PTEnvel2) :: envelope
+
+   integer, parameter :: nc=2
+   real(pr) :: z(nc), tc(nc), pc(nc), w(nc)
+
+   tc = [190.564, 425.12]     ! Critical temperatures
+   pc = [45.99, 37.96]        ! Critical pressures
+   w = [0.0115478, 0.200164]  ! Acentric factors
+
+   ! Get the example PR76 binary model
+   model = PengRobinson76(tc, pc, w)
+
+   ! Composition
+   z = [0.1_pr, 0.9_pr]
+
+   ! Calculate a bubble point at 100K to initialize the rest of the phase
+   ! envelope calculation
+   bubble = saturation_pressure(model, z, T=150._pr, kind="bubble")
+
+   ! Calculate the whole phase envelope using the information from the converged
+   ! dew point
+   envelope = pt_envelope_2ph(model, z, bubble)
+
+   ! Write the resulting phase envelope to the screen
+   write(*, *) envelope
+
+   ! Write the resulting phase envelope to a file
+   open(1, file="phase_envelope.dat")
+      write(1, *) envelope
+   close(1)
+
+   ! It is also possible to initialize the phase_envelope with a manually
+   ! defined point
+   manual: block
+      real(pr) :: k(nc), T, P
+      T = 150
+      P = 1
+      k = k_wilson(model, T=T, P=P)
+
+      ! Define manually an initial point
+      bubble%kind = "bubble"
+      bubble%x = z
+      bubble%y = k*z
+      bubble%beta = 0
+
+      envelope = pt_envelope_2ph(model, z, bubble)
+   end block manual
+end program phase_envelope

example/extra/pure_psat.f90 → example/basics/5_pure_psat.f90

@@ -1,16 +1,17 @@
 !! Program to calculate the vapor pressure of pure components
-!!
-
 module pure_psat
    !! Module used to calculate the saturation pressure of pure components at
    !! a given temperature.
    use yaeos
 contains
    real(pr) function Psat(eos, ncomp, T)
-      use yaeos__math, only: solve_system
-      class(ArModel), intent(in) :: eos
-      integer, intent(in) :: ncomp
-      real(pr), intent(in) :: T
+      !! Calculation of saturation pressure of a pure component using the
+      !! secant method.
+      class(ArModel), intent(in) :: eos !! Model that will be used
+      integer, intent(in) :: ncomp 
+         !! Number of component in the mixture from which the saturation pressure
+         !! will be calculated
+      real(pr), intent(in) :: T !! Temperature [K]
 
       real(pr) :: P1, P2
       real(pr) :: f1, f2
@@ -30,7 +31,6 @@ real(pr) function Psat(eos, ncomp, T)
          P1 = P2
          P2 = Psat
       end do
-
    contains
       real(pr) function diff(P)
          real(pr), intent(in) :: P
@@ -47,11 +47,13 @@ program main
    use yaeos
    use forsus, only: Substance, forsus_default_dir, forsus_dir
    use pure_psat, only: Psat
+
    implicit none
+
    integer, parameter :: nc=2
    type(CubicEoS) :: eos
    type(Substance) :: sus(nc)
-   real(pr) :: n(nc), T,  f
+   real(pr) :: n(nc), T
    integer :: i, j
 
    forsus_dir = "build/dependencies/forsus/" // forsus_default_dir
@@ -62,6 +64,8 @@ program main
       sus%critical%critical_pressure%value/1e5,&
       sus%critical%acentric_factor%value &
       )
+
+
    T = 273.15_pr + 50
    do i=273+90, nint(maxval(sus%critical%critical_temperature%value))
       T = real(i,pr)

example/extra/benchmark.f90

@@ -2,7 +2,7 @@ module bench
     use yaeos, only: pr, R, Substances, AlphaSoave, CubicEoS, &
         QMR, PengRobinson76, ArModel
     use hyperdual_pr76, only: PR76, setup_adiff_pr76 => setup
-    use TapeRobinson, only: setup_taperobinson => setup_model
+    ! use TapeRobinson, only: setup_taperobinson => setup_model
     implicit none
 
     real(pr), allocatable :: z(:), tc(:), pc(:), w(:), kij(:,:), lij(:,:)
@@ -40,7 +40,7 @@ subroutine yaeos__run(n, dn, f_p, model_name)
         case ("Adiff PR76")
             model = setup_adiff_pr76(tc, pc, w, kij, lij)
         case ("Tape PR76")
-            model = setup_taperobinson(tc, pc, w, kij, lij)
+            ! model = setup_taperobinson(tc, pc, w, kij, lij)
         end select
 
         v = 1.0_pr

example/extra/demo.f90

@@ -1,18 +1,18 @@
 program examples
-    use bench, only: benchmarks => main
+    ! use bench, only: benchmarks => main
     use hyperdual_pr76, only: run_hyperdual_pr76 => main
     use flashing, only: run_flashes => main
-    use TapeRobinson, only: run_tape_pr76 => main
     use tape_nrtl, only: run_tape_nrtl => main
 
-    ! print *, "Running Tapenade generated NRTL model"
-    ! call run_tape_nrtl
-    ! print *, "Running Tapenade generated PR76"
-    ! call run_tape_pr76
-    ! print *, "Running Hyperdual generated PR76"
-    ! call run_hyperdual_pr76
-    print *, "Running bencharks O(f(N))"
-    call benchmarks
-    ! print *, "Flash example"
-    ! call run_flashes
+    print *, "========================================="
+    print *, "YAEOS DEMO"
+    print *, "Running Tapenade generated NRTL model"
+    call run_tape_nrtl
+    print *, "Running Hyperdual generated PR76"
+    call run_hyperdual_pr76
+    ! print *, "Running bencharks O(f(N))"
+    ! call benchmarks
+    print *, "Flash example"
+    call run_flashes
+    print *, "========================================="
 end program