Merge pull request #130 from ipqa-research/dev

v2.1

app/critical.f90

@@ -26,28 +26,40 @@ program main
    real(pr) :: F(3), X(3)
    integer :: i, j
 
+   integer :: unit_cl
+   integer :: unit_pt
+   integer :: unit_pt_cp, unit_pt_hpl
+
    model = get_model()
 
-   a = real(1, pr)/100._pr
+   open(newunit=unit_cl, file="tmp/cl.dat")
+   open(newunit=unit_pt, file="tmp/pt.dat")
+   open(newunit=unit_pt_cp, file="tmp/pt_cp.dat")
+   open(newunit=unit_pt_hpl, file="tmp/pt_hpl.dat")
+
    print *, "1stCL"
-   cl = critical_line(model, a, z0, zi, 0.1_pr)
+   a = real(1, pr)/100._pr
+   cl = critical_line(model, a0=a, z0=z0, zi=zi, ns=spec_CP%a, S=a, dS0=0.1_pr)
+   print *, size(cl%a)
    do i=1, size(cl%a)
-      write(2, *) cl%a(i), cl%V(i), cl%T(i), cl%P(i)
+      write(unit_cl, *) cl%a(i), cl%V(i), cl%T(i), cl%P(i)
    end do
-   write (2, *)
-   write (2, *)
+   write (unit_cl, *)
+   write (unit_cl, *)
 
    print *, "2ndCL"
-   a = 0.001
-   cl = critical_line(model, a0=a, z0=zi, zi=z0, dS0=0.01_pr)
+   a = 1-epsilon(1._pr)
+   cl = critical_line(model, a0=a, z0=z0, zi=zi, ns=spec_CP%a, S=a, dS0=-0.01_pr)
+   print *, size(cl%a)
    do i=1, size(cl%a)
-      write(2, *) 1-cl%a(i), cl%V(i), cl%T(i), cl%P(i)
+      write(unit_cl, *) cl%a(i), cl%V(i), cl%T(i), cl%P(i)
    end do
 
    z = a*zi + (1-a)*z0
    T = sum(model%components%Tc * z)
    P = sum(model%components%Pc * z)
    call model%volume(n=z, P=P, T=T, V=V, root_type="stable")
+
    X = [a, log(V), log(T)]
    ns = 1
    S = X(ns)
@@ -60,39 +72,34 @@ program main
    write(20, *) env
    write(21, *) env
 
-
-
-   open(unit=4, file="pt")
-   open(unit=60, file="pt_cp")
-   open(unit=61, file="pt_hpl")
-
-!   !$OMP PARALLEL DO PRIVATE(j, a, z, sat, env, i) shared(model, z0, zi)
-   do j=9999999, 999999999, 1000000
-      print *, j
-      a = real(j, pr)/1000000000
+   !$OMP PARALLEL DO PRIVATE(j, a, z, sat, env, i) shared(model, z0, zi)
+   do j=1, 99, 3
+      a = real(j, pr)/100
       z = a*zi + (1-a)*z0
       sat = saturation_temperature(model, z, P=0.01_pr, kind="dew")
       env = pt_envelope_2ph(model, z, sat)
 
       do i=1,size(env%points)
-         write(4, *) a, env%points(i)%T, env%points(i)%P
+         write(unit_pt, *) a, env%points(i)%T, env%points(i)%P
       end do
-      write(4, *)
-      write(4, *)
+      write(unit_pt, *)
+      write(unit_pt, *)
 
-      write(60, *) a, env%cps
+      write(unit_pt_cp, *) a, env%cps
 
       env = find_hpl(model, z, t0=500._pr, P0=1000._pr)
       do i=1,size(env%points)
-         write(61, *) a, env%points(i)%T, env%points(i)%P
+         write(unit_pt_hpl, *) a, env%points(i)%T, env%points(i)%P
       end do
-      write(61, *)
-      write(61, *)
+      write(unit_pt_hpl, *)
+      write(unit_pt_hpl, *)
    end do
-!   !$OMP END PARALLEL DO
-   close(4)
-   close(60)
-   close(61)
+   !$OMP END PARALLEL DO
+
+   close(unit_cl)
+   close(unit_pt)
+   close(unit_pt_cp)
+   close(unit_pt_hpl)
 contains
 
    type(CubicEoS) function get_model()

app/gpec.f90

@@ -10,13 +10,15 @@ program gpec
    class(ArModel), allocatable :: model !! Thermodynamic model to use
 
    type(EquilibriumState) :: sat_point !! Saturation point
+   type(EquilibriumState) :: cp !! Individual critical point
    type(PTEnvel2) :: psats(2) !! Saturation curves
    type(CriticalLine) :: cl
 
    real(pr) :: z(nc) !! Molar fractions
    real(pr) :: a !! Fraction between component 1 and 2
-   real(pr), parameter :: z0(nc) = [1, 0] !! Component 1 molar fractions
-   real(pr), parameter :: zi(nc) = [0, 1] !! Component 2 molar fractions
+   real(pr), parameter :: eps = 1e-300
+   real(pr), parameter :: z0(nc) = [1-eps, eps] !! Component 1 molar fractions
+   real(pr), parameter :: zi(nc) = [eps, 1-eps] !! Component 2 molar fractions
    real(pr) :: P !! Pressure [bar]
    real(pr) :: T !! Temperature [K]
    real(pr) :: V !! Volume [L/mol]
@@ -33,7 +35,7 @@ program gpec
    ! ---------------------------------------------------------------------------
    ! model = get_model_nrtl_mhv()
    sus(1) = Substance("methane")
-   sus(2) = Substance("n-butane")
+   sus(2) = Substance("carbon dioxide")
    model = PengRobinson76(&
       Tc=sus%critical%critical_temperature%value, &
       Pc=sus%critical%critical_pressure%value/1e5, &
@@ -46,41 +48,71 @@ program gpec
    print *, model%components%Tc
    print *, model%components%Pc
    do i=0,1
-      sat_point = critical_point(model, z0, zi, S=S, spec=CPSpec%a, max_iters=1000, a0=S)
-
-      print *, sat_point%iters, sat_point
-      psats(i+1) = pt_envelope_2ph(model, z, sat_point, delta_0=-0.1_pr)
+      S = i
+      z = zi*S + z0*(1-S)
+      sat_point = critical_point(model, z0, zi, S=S, spec=spec_CP%a, max_iters=1000, a0=S)
+
+      select case(i)
+       case(0)
+         sat_point = saturation_pressure(model, z, T=100._pr, kind="bubble")
+       case(1)
+         sat_point = saturation_pressure(model, z, T=100._pr, kind="bubble")
+      end select
+
+      psats(i+1) = pt_envelope_2ph(model, z, sat_point, delta_0=0.001_pr)
    end do
-   call exit
 
    print *, "Psat 1"
+   open(unit=1, file="gpec_psat1.dat")
    do i=1,size(psats(1)%points)
-      print *, psats(1)%points(i)%T, psats(1)%points(i)%P
+      write(1, *) psats(1)%points(i)%T, psats(1)%points(i)%P
    end do
+   close(1)
 
    print *, ""
    print *, ""
-   
+
    print *, "Psat 2"
+   open(unit=1, file="gpec_psat2.dat")
    do i=1,size(psats(2)%points)
       print *, psats(2)%points(i)%T, psats(2)%points(i)%P
    end do
+   close(2)
 
    ! ===========================================================================
    ! Calculate the first critical line (2 -> 1)
    ! ---------------------------------------------------------------------------
-   cl = critical_line(model, a0=0.99_pr, z0=z0, zi=zi, dS0=-0.01_pr)
+   cl = critical_line(model, a0=0.99_pr, z0=z0, zi=zi, ns=spec_CP%a, S=0.99_pr, dS0=-0.01_pr)
+   open(unit=1, file="gpec_cl2.dat")
    do i=1,size(cl%a)
       write(1, *) cl%a(i), cl%T(i), cl%P(i), cl%V(i)
    end do
+   close(1)
 
-   call exit
+
+   cl = critical_line(model, a0=0.001_pr, z0=z0, zi=zi, ns=spec_CP%a, S=0.001_pr, dS0=0.001_pr)
+   open(unit=1, file="gpec_cl1.dat")
+   do i=1,size(cl%a)
+      write(1, *) cl%a(i), cl%T(i), cl%P(i), cl%V(i)
+   end do
+   close(1)
 
-   ! cl = critical_line(model, a0=0.001_pr, z0=z0, zi=zi, dS0=0.001_pr)
-   ! do i=1,size(cl%a)
-   !    write(2, *) cl%a(i), cl%T(i), cl%P(i), cl%V(i)
-   ! end do
+   cp = critical_point(model, z0, zi, S=log(200._pr), spec=spec_CP%P, max_iters=1000, a0=0.5_pr)
+
+   !TODO: Si inicializo con S=200 converge a otro lado, ver por qué pasa!
+   cl = critical_line(model, a0=cp%x(2), z0=z0, zi=zi, ns=spec_CP%P, S=log(cp%P), dS0=-0.01_pr)
+   open(unit=1, file="gpec_cl3.dat")
+   do i=1,size(cl%a)
+      write(1, *) cl%a(i), cl%T(i), cl%P(i), cl%V(i)
+   end do
+   close(1)
+   print *, cp%iters, cp
+
+
+
 
+   call exit
+
    call plot_pts([(real(i,pr)/100, i=1,99,10)])
 
    if (cl%a(size(cl%a)) < 1e-3) then
@@ -135,10 +167,10 @@ subroutine plot_pts(zs)
          sat = saturation_pressure(model, z, T=200._pr, kind="bubble")
          env = pt_envelope_2ph(model, z, sat)
          write(i+10, *) env
-         
+
          sat = saturation_temperature(model, z, P=0.01_pr, kind="dew")
          env = pt_envelope_2ph(model, z, sat)
          write(i+10, *) env
       end do
-   end subroutine
+   end subroutine plot_pts
 end program gpec

app/tpd.f90

@@ -1,39 +1,65 @@
 program phase_diagram
    use forsus, only: Substance, forsus_dir
    use yaeos
-   use yaeos__phase_equilibria_stability, only: tm, min_tpd
-   use yaeos, only: flash
    implicit none
 
-   integer, parameter :: nc=2
+   integer, parameter :: nc=12
+
 
    class(ArModel), allocatable :: model
    type(Substance) :: sus(nc)
    real(pr) :: tc(nc), pc(nc), ac(nc)
    real(pr) :: z(nc), T, P
-   real(pr) :: w(nc), mintpd
-
-   forsus_dir = "build/dependencies/forsus/data/json"
-   sus(1) = Substance("methane")
-   sus(2) = Substance("hydrogen sulfide")
-
-   z = [0.13, 1-0.13]
-   z = z/sum(z)
-
-   P = 20.0_pr
-   T = 190._pr
-
-   tc = sus%critical%critical_temperature%value
-   pc = sus%critical%critical_pressure%value/1e5_pr
-   ac = sus%critical%acentric_factor%value
-
-   model = SoaveRedlichKwong(tc, pc, ac)
-
-   call min_tpd(model, z, P, T, mintpd, w)
-   print *, mintpd, w/sum(w)
-
-   P = 15
-   call min_tpd(model, z, P, T, mintpd, w)
-   print *, mintpd, w/sum(w)
+   real(pr) :: w(nc), mintpd, allmin(nc, nc), di(nc), dw(nc)
+   real(pr) :: val, lnphi_w(nc)
+   type(PTEnvel2) :: env
+   type(EquilibriumState) :: sat
+   integer :: i, j
+
+   P = 3
+   T = 140
+
+   model = get_model()
+
+   call min_tpd(model, z, P, T, mintpd, w, allmin)
+   ! print *, mintpd, w/sum(w)
+
+   do i=1,nc
+      print *, allmin(i,:)/sum(allmin(i,:))
+   end do
+
+   print *, "SS"
+
+   call model%lnphi_pt(z, P, T, root_type="stable", lnPhi=di)
+   di = log(z) + di
+
+   sat = saturation_temperature(model, z, P=0.001_pr, kind="dew", T0=500._pr)
+   env = pt_envelope_2ph(model, z, sat)
+   write(1, *) env
+   env = find_hpl(model, z, T0=300._pr, p0=maxval(env%points%P))
+   write(1, *) env
+
+   ! w =[0.001, 0.999]
+
+contains
+
+   type(CubicEoS) function get_model()
+      real(pr) :: tc(nc), pc(nc), w(nc), kij(nc,nc)
+      z=[0.0656,0.3711,0.0538,0.0373,0.0261,0.0187,0.0218,0.1791,0.091,0.0605,0.0447,0.0302]
+
+      tc=[304.088888888889,190.6,305.4,369.8,425.2,469.6,507.4,616.2,698.9,770.4,853.1,1001.2]
+      pc=[73.7343491450634,45.9196083838941,48.7516547159404,42.3795504688362, &
+          37.9291919470491,33.6811224489796,29.6353419746277,28.8261858797573,&
+          19.3186017650303,16.5876999448428,15.2728212906784,14.6659542195256]
+      w= [0.228,0.008,0.098,0.152,0.193,0.251,0.296,0.454,0.787,1.048,1.276,1.299]
+
+      kij = 0
+      kij(1, 2) = 0.12
+      kij(1, 3:) = 0.15
+      kij(:, 1) = kij(1, :)
+
+      get_model = PengRobinson78(tc, pc, w, kij=kij)
+   end function get_model
+
 
 end program phase_diagram