-
Notifications
You must be signed in to change notification settings - Fork 4
/
mod_pcsaft.f08
3894 lines (3365 loc) · 186 KB
/
mod_pcsaft.f08
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
!-----------------------------------------------------------------------------------------------------------------------------------
!
! -- CSTES, TYPE definitions for the PC-SAFT EOS --
!
!-----------------------------------------------------------------------------------------------------------------------------------
! D. Cordier, Institut UTINAM, Besançon, France (until september, 2015).
! GSMA, Reims, France (from september, 2015).
!
! - Research Gate ...: https://www.researchgate.net/profile/Daniel_Cordier
! - OrCID ...........: http://orcid.org/0000-0003-4515-6271
!
! - 3 novembre 2013 --: écriture première version.
! - 21 janvier 2014 --: introduction du fichier contenant les donnes pour les paramtres d'interation "k_ij".
! - 18 novembre 2019 --: passage des constantes de "1.3806488d-23" "1.3806488E-23_dp".
! - 7 mai 2020 -------: introduction d'une nouvelle subroutine 'extract_species_prop' qui permet d'extraire directement
! dans la base de données 'PC-SAFT' certaines propriétés des espèces étudiées. Cette souroutine permet d'utiliser
! ces données depuis un programme appelant, avant tout appel à PC-SAFT stricto-sensus.
! - 12 mai 2020 -------: - mise en place de 'ares' comme sortie optionelles de 'Zandmu_pcsaft_PT'.
! - mise en place de la sortie de 'h^res/RT' de 'pcsaft_PT'.
! - 30 juin 2020 ------: - mise en place du paramètre logique 'normized_xi_sum' qui contrôle la condition de normalisation des fractions molaires
! à l'entrée de certaines subroutines. Cela permet, par exemple, de calculer des dérivées numériques, par différences
! finie, sans se préoccuper de cette condition de normalisation.
! - 11 janvier 2021 ---: - introduction de la subroutine 'alphaV', qui provide the volumetric thermal expansion coefficient (K^-1).
!
module mod_pcsaft
use utils_DC
! =========================================================================================================================================
! Paramètre qui permet de contrôler la condition de normalisation des fractions molaires à l'entrée de certaines subroutines,
! Cela permet, par exemple, de calculer des dérivées numériques, par différences finie, sans se préoccuper de cette condition de normalisation.
logical, parameter :: normized_xi_sum = .true. ! In this case, we force the sum of mole fraction to equal the unity !
!logical, parameter :: normized_xi_sum = .false. ! In this case, we _DO NOT_ force the sum of mole fraction to equal the unity !
! =========================================================================================================================================
character(len=*), private, parameter :: mod_name = 'mod_pcsaft' ! Name of this module (PRIVATE).
character(len=*), parameter :: datafilename = './DATABASE_PCSAFT/COMPOUNDS_DATA_PC-SAFT_withASSO.data' ! PC-SAFT parameters database file name.
! in this file there are: Molar_mass, Num._seg., Seg._diam., Seg._engergy and References.
character(len=*), parameter :: datafilename_inter = './DATABASE_PCSAFT/COMPOUNDS_DATA_PC-SAFT-PARAM_INTERAC.data' ! PC-SAFT interaction
! parameters database file name, in this file one can find the k_ij's and/or the parameters "a" and "b"
! required to compute the "k_ij" at a given temperature.
! Physical constants used this implementation of PC-SAFT:
integer, parameter :: comp_name_ncmax = 15 ! Max. number of characters for compound names.
integer, parameter :: comp_ref_ncmax = 100 ! Max. number of characters for compound reference (for PC-SAFT parameters).
real(dp), parameter :: kB = 1.3806488E-23_dp ! Boltzmann's constant in J.K^-1
real(dp), parameter :: angstrom_meter = 1.E-10_dp ! one Angstrom in meter
real(dp), parameter :: Navo = 6.0221413E+23_dp ! Avogadro's number
real(dp), parameter :: onebar_in_Pa = 1.0000E+5_dp ! One bar in pascal
real(dp), parameter :: onekgm3_in_gcm3 = 1.0000E-3_dp ! One kg.m^-3 in g.cm^-3
real(dp), parameter :: R_gas = 8.3144621_dp ! The gas contant in J.K^-1.mol^-1
! Definitions of derived types associated to chemical species properties:
type compdata
character(len=comp_name_ncmax) :: name ! The name of the considered compound.
real(dp) :: molmass ! The molar mass (in kg.mol^-1) for the considered compound.
real(dp) :: numseg ! The number of segment for the considered compound.
real(dp) :: segdiam ! The segment diameter.
real(dp) :: segenerg ! The segment energy
integer :: assoflag ! Flag: 0=non associative species, 1= associative species
real(dp) :: assenerg ! The association energy
real(dp) :: assvol ! The association volume
character(len=comp_ref_ncmax) :: ref ! The reference in which these parameters have been taken.
end type compdata
type compdataINTER
character(len=comp_name_ncmax) :: name1 ! The name of the FIRST considered compound.
character(len=comp_name_ncmax) :: name2 ! The name of the SECOND considered compound.
real(dp) :: a ! The "a" coefficient in Eq. (8) of Tan et al. (2013).
real(dp) :: b ! The "b" coefficient in Eq. (8) of Tan et al. (2013).
character(len=comp_ref_ncmax) :: ref ! The reference in which these parameters have been taken.
end type compdataINTER
type compprop
character(len=comp_name_ncmax) :: name ! The name of the considered compound.
real(dp) :: x ! Mole fraction
real(dp) :: molmass ! The molar mass (in kg.mol^-1) for the considered compound.
real(dp) :: numseg ! The number of segment for the considered compound.
real(dp) :: segdiam ! The segment diameter.
real(dp) :: segenerg ! The segment energy
integer :: assoflag ! Flag: 0=non associative species, 1= associative species
real(dp) :: assenerg ! The association energy
real(dp) :: assvol ! The association volume
end type compprop
! --------------------------------------------------------------------------------------------------------------------------------
! Subroutines available in this module:
public :: read_database_pcsaft, & ! Read the 'PC-SAFT database', i.e. files containing the parameters for numerous species.
extract_species_prop, & ! Extract species properties from 'PC-SAFT database', useful if we need these propertie outside PC-SAFT.
P_pcsaft_rho, & ! Provides 'P' when temperature T and density rho are given.
pcsaft_rhoT, pcsaft_PT, pcsaft_Prho, &
alphaV, & ! Provide the volumetric thermal expansion coefficient (K^-1)
cp_solid, cpcv_ssound, cp0_idealgas_joback
! --- PC-SAFT internal subroutines:
private :: P_pcsaft_eta, Zandmu_pcsaft_PT, Zandmu_pcsaft_Prho, univ_cst
private :: association_terms
private :: intext_linear ! Subroutines utilitaires, également fournie par le module '../../prog_Fortran/INTERPOL/interpol_mod.f08',
! avec lequel elle ne rentre pas en conflit du fait de sa déclaration en 'private' ici.
!===================================================================================================================================
!===================================================================================================================================
contains
!===================================================================================================================================
subroutine extract_species_prop (path, ncall, compound, species_prop)
! Goal of this subroutine: extract the properties of the species specified in 'compound' and leave them in 'species_prop'.
! D. Cordier, 7 mai 2020.
use utils_dc
implicit none
character(len=*), parameter :: subname = 'extract_species_prop'
! Inputs/outputs of the program:
character(len=*), intent(in) :: path
integer, intent(in) :: ncall
character(len=len_trim(path)+len_trim(subname)+1) :: subtemp
character(len=comp_name_ncmax), dimension(:), intent(in) :: compound ! List of the names of the involved species.
type(compprop), dimension(:), intent(out) :: species_prop ! proprieties associated to involved species.
integer(si) :: i, j
logical :: first = .true.
logical :: present
character(len=6) :: temp_string
character(len=comp_name_ncmax) :: temp_string2
type(compdata), dimension(:), allocatable, save :: compdatabase
type(compdataINTER), dimension(:), allocatable, save :: compdatabase_interac
character(len=comp_name_ncmax) :: name1, name2, name3, name4
! --------------------------------------------------------------------------------------------------------------------------------
! We read the database of PC-SAFT parameters:
if (first) then
call read_database_pcsaft(subname,1,compdatabase,compdatabase_interac)
first= .false.
write(6,*) ''
write(temp_string,'(I4)') size(compdatabase)
call wf(' > Size of the database: ', 'bright red', &
temp_string(1:len_trim(temp_string)), 'bright black', &
' species:', 'bright red')
do i= 1, size(compdatabase)
temp_string2= compdatabase(i)%name
call wf(' - ', 'bright red', &
temp_string2(1:len_trim(temp_string2)), 'bright black' )
end do
write(6,*) ''
call wf(' > Species taken into account in this computation: ', 'bright red')
do i= 1, size(compound)
temp_string2= compound(i)
call wf(' - ', 'bright red', &
temp_string2(1:len_trim(temp_string2)), 'bright blue' )
end do
write(6,*) ''
end if
do i= 1, size(compound)
name1 = trim(ADJUSTL(compound(i))) ! Enlève les blancs en début et fin de chaîne.
present= .false.
do j= 1, size(compdatabase)
name2= trim(ADJUSTL(compdatabase(j)%name))
!write(6,*) name2
if ( name1(1:len_trim(name1)) == name2(1:len_trim(name2)) ) then
present= .true.
species_prop(i)%name = name1 ! Name
species_prop(i)%x = -99. ! Mole fraction, UNDEFINED HERE -> -99. magic number!
species_prop(i)%molmass = compdatabase(j)%molmass ! The molar mass (in kg.mol^-1) for the considered compound.
species_prop(i)%numseg = compdatabase(j)%numseg ! The number of segment for the considered compound.
species_prop(i)%segdiam = compdatabase(j)%segdiam ! The segment diameter.
species_prop(i)%segenerg = compdatabase(j)%segenerg ! The segment energy
! Next a series of parameters UNDEFINED here -> magic number.
species_prop(i)%assoflag = -99._dp ! Flag: 0=non associative species, 1= associative species
species_prop(i)%assenerg = -99._dp ! The association energy
species_prop(i)%assvol = -99._dp ! The association volume
exit
end if
end do
if ( .NOT. present ) then
write(6,*) ''
call wf(' > In module: ', 'bright red', mod_name(1:len_trim(mod_name)), 'bright yellow')
write(6,*) ' > In subroutine ''', subname(1:len_trim(subname)), ''' called by ''', subtemp(1:len_trim(subtemp)), ''''
write(6,*) ' the species "', name1(1:len_trim(name1)), '" is not present in the database!'
write(6,*) ''
stop
end if
end do
return
end subroutine extract_species_prop
!===================================================================================================================================
!-----------------------------------------------------------------------------------------------------------------------------------
!
! -- Computation of "Z" (compressibility factor) and "mu^res" in the frame of the PC-SAFT theory --
!
!-----------------------------------------------------------------------------------------------------------------------------------
! D. Cordier, Institut UTINAM, Besançon, France.
! 1er novembre 2013 -- [email protected]
!
! References: * Sunil Kumar Maity, "Modeling and simulation of solid-liquid equilibrium by
! perturbed-chain statistical associating fluid theory", Master of Technology in Chemical Engineering,
! 2003,
! Indian Institute of Technology, Kharagpur, India.
! * Abdelkrim Belkadi, Thse de l'Universit de Toulouse, 2008.
! * Gross & Sadowski, Ind. Eng. Chem. Res. 2001, 40, 1244-1260.
!
! - June 2sd 2014: change the name 'Zandmu_pcsaft' -> 'Zandmu_pcsaft_PT'.
! - 27 novembre 2014 : correction du nom de la routine dans 'subname'.
!
! - 15 avril 2020 : - passage des constantes numriques de l'criture '1.d0' vers l'criture '1._dp'
! - mise dans le module 'NEW_MOD_PCSAFT'.
!
! - 12 mai 2020 : - introduction d'une sortie optionelle : 'ares'.
!
subroutine Zandmu_pcsaft_PT(sub, ncall, state, species_prop, k_ij, P, T, Z, muskT_k, eta, rho, ares)
use utils_DC
implicit none
character(len=*), parameter :: subname = 'Zandmu_pcsaft_PT'
! Inputs/outputs of the program:
character(len=*), intent(in) :: sub
integer, intent(in) :: ncall
character(len=len_trim(sub)+len_trim(subname)+1) :: subtemp
character(len=1), intent(in) :: state ! the physical state of the phase: - state='L' = liquid phase
! - state='V' = vapor phase.
type(compprop), dimension(:), intent(in) :: species_prop ! The properties of studied species.
real(dp), dimension(:,:), intent(in) :: k_ij ! interaction parameters at the given temperature T
real(dp), intent(in) :: P, T ! Pressure (in Pa) and temperature (in K).
real(dp), intent(out) :: Z ! Compressibility factor (no unit)
real(dp), dimension(:), intent(out) :: muskT_k ! The residual chemical potential divided by kT (no unit)
real(dp), intent(out) :: eta ! the reduced density (no unit)
real(dp), intent(out) :: rho ! the density (kg.m^-3) of the system at P and T
real(dp), optional, intent(out) :: ares ! Residual (compare to ideal gas value) Helmoltz energy (J.mol^-1).
integer :: ntour_NR
integer, parameter :: ntour_NR_max = 100 ! Max. number of tours in the Newton-Raphson's algorithm.
real(dp), parameter :: delta_NR_max = 1.E-8_dp ! Convergence criterion for the Newton-Raphson's algorithm.
real(dp), parameter :: pd = 1.E-5_dp
real(dp) :: Psaft, Psaft1, Psaft2, Psaft_0, eta1, eta2, rho1, rho2, Z1, Z2, dP_pcsaftsdeta, Delta_P
real(dp) :: ln_eta, ln_eta1, ln_eta2, Delta_P1, Delta_P2
logical :: increase
!---------------------------------------------------------------------------------------------------------------------------------
! Call path of this subroutine:
subtemp= sub(1:len_trim(sub))//'/'//subname(1:len_trim(subname))
!---------------------------------------------------------------------------------------------------------------------------------
! We determine the value of "eta" that matches the pressure of the system P (i.e. to get Psaft=P):
!
if ( state == 'L' ) then
eta= 0.5_dp
end if
if ( state == 'V' ) then
eta= 1.E-10_dp
end if
!---------------------------------------------------------------------------------------------------------------------------------
! Recherche par Newton-Raphson :
! Initialisation de la routine 'P_pcsaft_eta' :
call P_pcsaft_eta(subtemp,1,species_prop,k_ij,T,eta,Psaft_0,rho,Z,muskT_k)
! The Newton-Raphson's algorithm itself:
Delta_P= Psaft_0/P-1._dp
ntour_NR= 1
do while ( abs(Delta_P) > delta_NR_max)
call P_pcsaft_eta(subtemp,2,species_prop,k_ij,T,eta,Psaft_0,rho,Z,muskT_k)
Delta_P= Psaft_0/P-1._dp
eta1= eta * (1._dp-pd)
call P_pcsaft_eta(subtemp,3,species_prop,k_ij,T,eta1,Psaft1,rho1,Z1,muskT_k)
eta2= eta * (1._dp+pd)
call P_pcsaft_eta(subtemp,4,species_prop,k_ij,T,eta2,Psaft2,rho2,Z2,muskT_k)
dP_pcsaftsdeta= (Psaft2-Psaft1)/(eta2-eta1)/P
eta= eta - Delta_P/dP_pcsaftsdeta
!if ( eta < 0.d0 ) then
! print*, ' eta < 0.d0 on impose eta= 1.d-10 ! (Dans Zandmu_pcsaft)'
! eta= 1.d-20
!end if
ntour_NR= ntour_NR + 1
if ( ntour_NR > ntour_NR_max ) then
write(6,*) ''
call wf(' > In module: ', 'bright red', mod_name(1:len_trim(mod_name)), 'bright yellow')
write(6,*) ' > In subroutine ''', subname(1:len_trim(subname)), ''' called by ''', subtemp(1:len_trim(subtemp)), ''''
write(6,*) ' the number of iteration in the Newton-Raphson''s algorithm is larger than the allowed maximum value!'
write(6,*) ' ntour_NR = ', ntour_NR
write(6,*) ' ntour_NR_max = ', ntour_NR_max
write(6,*) ''
stop
end if
end do
!
! Fin recherche par Newton-Raphson
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
! Recherche par dichotomie :
!!
!!call P_pcsaft_eta(subname,1,species_prop,T,eta1,Psaft1,rho,Z,muskT_k)
!!Delta_P1= Psaft1/P-1.d0
!!
!!call P_pcsaft_eta(subtemp,2,species_prop,T,eta2,Psaft2,rho,Z,muskT_k)
!!Delta_P2= Psaft2/P-1.d0
!!
!!if ( Delta_P1 * Delta_P2 >= 0.d0 ) then
!! write(6,*) ''
!! write(6,*) ' > In subroutine ''', subname(1:len_trim(subname)), ''' called by ''', subtemp(1:len_trim(subtemp)), ''''
!! write(6,*) ' Delta_P1 * Delta_P2 => 0.d0'
!! write(6,*) ' the bisection method can be applied here, please check the values of "eta1" and "eta2"'
!! write(6,*) ''
!! stop
!!end if
!!
!!if ( (Delta_P1 < 0.d0) .AND. (Delta_P2 > 0.d0) ) then
!! increase= .true.
!!else
!! increase= .false.
!!end if
!!
!!ntour_NR= 0
!!eta= (eta1+eta2)/2.d0
!!do while (abs(eta1-eta2)/eta > 1.d-6)
!! eta= (eta1+eta2)/2.d0
!! call P_pcsaft_eta(subtemp,3,species_prop,T,eta,Psaft_0,rho,Z,muskT_k)
!! Delta_P= Psaft_0/P-1.d0
!! if (increase) then
!! if (Delta_P > 0.d0) then
!! eta2= eta
!! else
!! eta1= eta
!! end if
!! else
!! if (Delta_P > 0.d0) then
!! eta1= eta
!! else
!! eta2= eta
!! end if
!! end if
!! !print*, ' eta= ', eta
!! !read(5,*)
!! ntour_NR= ntour_NR + 1
!! if ( ntour_NR > ntour_NR_max ) then
!! write(6,*) ''
!! write(6,*) ' > In subroutine ''', subname(1:len_trim(subname)), ''' called by ''', sub(1:len_trim(sub)), ''''
!! write(6,*) ' the number of iteration in the bisection method''s algorithm is larger than the allowed maximum value!'
!! write(6,*) ' ntour_NR = ', ntour_NR
!! write(6,*) ' ntour_NR_max = ', ntour_NR_max
!! write(6,*) ''
!! stop
!! end if
!!end do
! Fin recherche par dichotomie.
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
! Final values of "rho", "Z", "muskT_k":
if (present(ares)) then
call P_pcsaft_eta(subtemp, 5, species_prop, k_ij, T, eta, Psaft_0, rho, Z, muskT_k, ares)
else
call P_pcsaft_eta(subtemp, 6, species_prop, k_ij, T, eta, Psaft_0, rho, Z, muskT_k)
end if
!---------------------------------------------------------------------------------------------------------------------------------
return
end subroutine Zandmu_pcsaft_PT
!===================================================================================================================================
!-----------------------------------------------------------------------------------------------------------------------------------
!
! -- Computation of "Z" (compressibility factor) and "mu^res" in the frame of the PC-SAFT theory --
!
!-----------------------------------------------------------------------------------------------------------------------------------
! D. Cordier, Institut UTINAM, Besançon, France.
! June 2sd 2014 -- [email protected]
!
! References: * Sunil Kumar Maity, "Modeling and simulation of solid-liquid equilibrium by
! perturbed-chain statistical associating fluid theory", Master of Technology in Chemical Engineering,
! 2003,
! Indian Institute of Technology, Kharagpur, India.
! * Abdelkrim Belkadi, Thse de l'Universit de Toulouse, 2008.
! * Gross & Sadowski, Ind. Eng. Chem. Res. 2001, 40, 1244-1260.
!
! - 27 novembre 2014 : correction du nom de la routine dans 'subname'.
!
! - 15 avril 2020 : - passage des constantes numriques de l'criture '1.d0' vers l'criture '1._dp'
! - mise dans le module 'NEW_MOD_PCSAFT'.
!
subroutine Zandmu_pcsaft_Prho(sub,ncall,state,species_prop,k_ij,P,rho,Z,muskT_k,T,ares)
use utils_DC
implicit none
character(len=*), parameter :: subname = 'Zandmu_pcsaft_Prho'
! Inputs/outputs of the program:
character(len=*), intent(in) :: sub
integer, intent(in) :: ncall
character(len=len_trim(sub)+len_trim(subname)+1) :: subtemp
character(len=1), intent(in) :: state ! the physical state of the phase: - state='L' = liquid phase
! - state='V' = vapor phase.
type(compprop), dimension(:), intent(in) :: species_prop ! The properties of studied species.
real(dp), dimension(:,:), intent(in) :: k_ij ! interaction parameters at the given temperature T
real(dp), intent(in) :: P, rho ! Pressure (in Pa) and density (in kg.m^-3).
real(dp), intent(out) :: Z ! Compressibility factor (no unit)
real(dp), dimension(:), intent(out) :: muskT_k ! The residual chemical potential divided by kT (no unit)
real(dp), intent(out) :: T ! the temperature (in K)
real(dp), intent(out), optional :: ares ! Residual (compare to ideal gas value) Helmoltz energy (in J.mol^-1)
integer :: ntour_NR
integer, parameter :: ntour_NR_max = 100000 ! Max. number of tours in the Newton-Raphson's algorithm.
real(dp), parameter :: delta_NR_max = 1.E-6_dp ! Convergence criterion for the Newton-Raphson's algorithm.
real(dp), parameter :: pd = 1.E-5_dp
real(dp) :: Psaft, Psaft1, Psaft2, Psaft_0, temp, temp1, temp2, Z1, Z2, dP_pcsaftsdtemp, Delta_P
!---------------------------------------------------------------------------------------------------------------------------------
! Call path of this subroutine:
subtemp= sub(1:len_trim(sub))//'/'//subname(1:len_trim(subname))
!---------------------------------------------------------------------------------------------------------------------------------
! We determine the value of "eta" that matches the pressure of the system P (i.e. to get Psaft=P):
!
!if ( state == 'L' ) then
! eta= 0.5
!end if
!if ( state == 'V' ) then
! eta= 1.d-10
!end if
!---------------------------------------------------------------------------------------------------------------------------------
! Recherche par Newton-Raphson :
! Initialisation de la routine 'P_pcsaft_rho' :
temp= 273._dp ! On initialise la valeur de la température.
call P_pcsaft_rho(subtemp,1,species_prop,k_ij,temp,rho,Psaft_0,Z,muskT_k)
! The Newton-Raphson's algorithm itself:
Delta_P= Psaft_0/P-1._dp
ntour_NR= 1
do while ( abs(Delta_P) > delta_NR_max)
call P_pcsaft_rho(subtemp,2,species_prop,k_ij,temp,rho,Psaft_0,Z,muskT_k)
Delta_P= Psaft_0/P-1._dp
temp1= temp * (1._dp-pd)
call P_pcsaft_rho(subtemp,3,species_prop,k_ij,temp1,rho,Psaft1,Z,muskT_k)
temp2= temp * (1._dp+pd)
call P_pcsaft_rho(subtemp,4,species_prop,k_ij,temp2,rho,Psaft2,Z,muskT_k)
dP_pcsaftsdtemp= (Psaft2-Psaft1)/(temp2-temp1)/temp
temp= temp - Delta_P/dP_pcsaftsdtemp
!print*, ' temp= ', temp
!read(5,*)
ntour_NR= ntour_NR + 1
if ( ntour_NR > ntour_NR_max ) then
write(6,*) ''
call wf(' > In module: ', 'bright red', mod_name(1:len_trim(mod_name)), 'bright yellow')
write(6,*) ' > In subroutine ''', subname(1:len_trim(subname)), ''' called by ''', subtemp(1:len_trim(subtemp)), ''''
write(6,*) ' the number of iteration in the Newton-Raphson''s algorithm is larger than the allowed maximum value!'
write(6,*) ' ntour_NR = ', ntour_NR
write(6,*) ' ntour_NR_max = ', ntour_NR_max
write(6,*) ''
stop
end if
end do
!
! Fin recherche par Newton-Raphson
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
! Final values of "rho", "Z", "muskT_k":
if (present(ares)) then
call P_pcsaft_rho(subtemp,5,species_prop,k_ij,temp,rho,Psaft2,Z,muskT_k,ares)
else
call P_pcsaft_rho(subtemp,6,species_prop,k_ij,temp,rho,Psaft2,Z,muskT_k)
end if
T= temp
!---------------------------------------------------------------------------------------------------------------------------------
return
end subroutine Zandmu_pcsaft_Prho
!===================================================================================================================================
!-----------------------------------------------------------------------------------------------------------------------------------
!
! -- PC-SAFT parameters: reading of data files --
!
!-----------------------------------------------------------------------------------------------------------------------------------
! D. Cordier, Institut UTINAM, Besançon, France.
!
! - 3 novembre 2013 : criture de la premire version.
! - 21 janvier 2014 : mise en place de la lecture du deuxme fichier (celui contenant les donnes des paramtres d'interaction
! "k_ij".
! - 18 septembre 2014 : scurit : dtection des doublons dans le fichier de 'k_ij'.
! - 27 septembre 2016 : rcriture des "do while" pour les lectures de fichiers de longueurs indtermines.
!
! - 15 avril 2020 : - mise dans le module 'NEW_MOD_PCSAFT'.
!
subroutine read_database_pcsaft(sub, ncall, compdatabase, compdatabase_interac)
use utils_DC
implicit none
character(len=*), parameter :: subname = 'read_database_pcsaft'
! Inputs/outputs of the program:
character(len=*), intent(in) :: sub
integer, intent(in) :: ncall
type(compdata), dimension(:), allocatable, intent(out) :: compdatabase
type(compdataINTER), dimension(:), allocatable, intent(out) :: compdatabase_interac
integer :: nl, ndieze, ntot, ok, i, j, iostatus
character(len=1) :: line1
character(len=300) :: line
logical :: here
!
!-----------------------------------------------------------------------------------------------------------------------------------
! Readin of the FIRST data file:
! We check that the PC-SAFT parameters database is really present in the directory of work:
inquire(file=datafilename(1:len_trim(datafilename)),exist=here)
if ( .NOT. here ) then
write(6,*) ''
call wf(' > In module: ', 'bright red', mod_name(1:len_trim(mod_name)), 'bright yellow')
write(6,*) ' > In subroutine ''', subname(1:len_trim(subname)), ''' called by ''', sub(1:len_trim(sub))
write(6,*) ' the file ''', datafilename(1:len_trim(datafilename)), ''' is not present'
write(6,*) ''
stop
end if
! We count the number of species taken into account by the database:
nl=0; ndieze=0
open(unit=1,status='unknown',form='formatted',file=datafilename(1:len_trim(datafilename)))
do
read(1,'(A)',iostat=iostatus) line1
if (iostatus < 0) then
exit
end if
if ( line1 /= '#' ) then
nl=nl+1
else
ndieze=ndieze+1
end if
end do
close(unit=1)
! We allocate the memory for the database:
allocate(compdatabase(1:nl),stat=ok); call alloc_error(subname,'compdatabase',ok)
! We read the database file:
open(unit=1,status='unknown',form='formatted',file=datafilename(1:len_trim(datafilename)))
ntot=nl+ndieze
nl=0
do i= 1, ntot
read(1,'(A)') line
if ( line(1:1) /= '#' ) then
nl=nl+1
!# | (kg.mol^-1) | "m" | "sig" (Ang) | e/kB (K) |
!N2 28.01340E-03 1.2414 3.2992 89.2230 NIST Used by Tan et al. (2013)
read(line,'(A15,ES17.5,F12.4,F14.4,F15.4,A100)') compdatabase(nl)%name, compdatabase(nl)%molmass, &
compdatabase(nl)%numseg, compdatabase(nl)%segdiam, &
compdatabase(nl)%segenerg, compdatabase(nl)%ref
end if
end do
close(unit=1)
!-----------------------------------------------------------------------------------------------------------------------------------
! Readin of the SECOND data file:
! We check that the PC-SAFT parameters database is really present in the directory of work:
inquire(file=datafilename_inter(1:len_trim(datafilename_inter)),exist=here)
if ( .NOT. here ) then
write(6,*) ''
call wf(' > In module: ', 'bright red', mod_name(1:len_trim(mod_name)), 'bright yellow')
write(6,*) ' > In subroutine ''', subname(1:len_trim(subname)), ''' called by ''', sub(1:len_trim(sub))
write(6,*) ' the file ''', datafilename_inter(1:len_trim(datafilename_inter)), ''' is not present'
write(6,*) ''
stop
end if
! We count the number of species taken into account by the database:
nl=0; ndieze=0
open(unit=1,status='unknown',form='formatted',file=datafilename_inter(1:len_trim(datafilename_inter)))
do
read(1,'(A)',iostat=iostatus) line1
if(iostatus < 0) then
exit
end if
if ( line1 /= '#' ) then
nl=nl+1
else
ndieze=ndieze+1
end if
end do
close(unit=1)
! We allocate the memory for the database:
allocate(compdatabase_interac(1:nl),stat=ok); call alloc_error(subname,'compdatabase_interac',ok)
! We read the database file:
open(unit=1,status='unknown',form='formatted',file=datafilename_inter(1:len_trim(datafilename_inter)))
ntot=nl+ndieze
nl=0
do i= 1, ntot
read(1,'(A)') line
if ( line(1:1) /= '#' ) then
nl=nl+1
read(line,'(2A15,2F12.4,A100)') compdatabase_interac(nl)%name1, compdatabase_interac(nl)%name2, &
compdatabase_interac(nl)%a, compdatabase_interac(nl)%b, &
compdatabase_interac(nl)%ref
!write(6,'(2A15,2F12.4,A100)') compdatabase_interac(nl)%name1, compdatabase_interac(nl)%name2, &
! compdatabase_interac(nl)%a, compdatabase_interac(nl)%b, &
! compdatabase_interac(nl)%ref
end if
end do
close(unit=1)
! We search for duplicates:
do i= 1, nl
!write(6,*) '''', compdatabase_interac(i)%name1, '''', '|', '''', compdatabase_interac(i)%name2, ''''
do j= 1, nl
! Cas o deux lignes se rptent simplement :
if ( (compdatabase_interac(i)%name1 .eq. compdatabase_interac(j)%name1) .AND. (i /= j) ) then
if ( compdatabase_interac(i)%name2 .eq. compdatabase_interac(j)%name2 ) then
write(6,*) ''
call wf(' > In module: ', 'bright red', mod_name(1:len_trim(mod_name)), 'bright yellow')
write(6,*) ' > In subroutine ''', subname(1:len_trim(subname)), ''' called by ''', sub(1:len_trim(sub))
write(6,*) ' there is a duplicate:'
write(6,*) ' ', i, compdatabase_interac(i)%name1, compdatabase_interac(i)%name2
write(6,*) ' ', j, compdatabase_interac(j)%name1, compdatabase_interac(j)%name2
write(6,*) ''
stop
end if
end if
!
if ( (compdatabase_interac(i)%name1 .eq. compdatabase_interac(j)%name2) .AND. (i /= j) ) then
if ( compdatabase_interac(i)%name2 .eq. compdatabase_interac(j)%name1 ) then
write(6,*) ''
call wf(' > In module: ', 'bright red', mod_name(1:len_trim(mod_name)), 'bright yellow')
write(6,*) ' > In subroutine ''', subname(1:len_trim(subname)), ''' called by ''', sub(1:len_trim(sub))
write(6,*) ' there is a duplicate:'
write(6,*) ' ', i, compdatabase_interac(i)%name1, compdatabase_interac(i)%name2
write(6,*) ' ', j, compdatabase_interac(j)%name1, compdatabase_interac(j)%name2
write(6,*) ''
stop
end if
end if
end do
end do
return
end subroutine read_database_pcsaft
!===================================================================================================================================
!-----------------------------------------------------------------------------------------------------------------------------------
!
! -- Computation of the pressure P in the frame of the PC-SAFT theory --
!
! >>> AS A FUNCTION OF "rho" and "T" <<<
!
!-----------------------------------------------------------------------------------------------------------------------------------
! D. Cordier, Institut UTINAM, Besançon, France.
! June 2sd 2014 -- [email protected]
!
! References: * Sunil Kumar Maity, "Modeling and simulation of solid-liquid equilibrium by
! perturbed-chain statistical associating fluid theory", Master of Technology in Chemical Engineering,
! 2003,
! Indian Institute of Technology, Kharagpur, India.
! WARNING: a lot of mistakes/typos in this references!
! * Abdelkrim Belkadi, Thse de l'Universit de Toulouse, 2008.
! * Gross & Sadowski, Ind. Eng. Chem. Res. 2001, 40, 1244-1260.
!
! - 15 avril 2020 : - passage des constantes numriques de l'criture '1.d0' vers l'criture '1._dp'
! - mise dans le module 'NEW_MOD_PCSAFT'.
!
subroutine P_pcsaft_rho(sub, ncall, species_prop, k_ij, T, rho_MKSA, Psaft, Z, muskT_k, ares)
use utils_DC
implicit none
character(len=*), parameter :: subname = 'P_pcsaft_rho'
! Inputs/outputs of the program:
character(len=*), intent(in) :: sub
integer, intent(in) :: ncall
character(len=len_trim(sub)+len_trim(subname)+1) :: subtemp
type(compprop), dimension(:), intent(in) :: species_prop ! The properties of studied species.
real(dp), dimension(:,:), intent(in) :: k_ij ! interaction parameters at the given temperature T
real(dp), intent(in) :: T ! Temperature (in K).
real(dp), intent(in) :: rho_MKSA ! the density (kg.m^-3)
real(dp), intent(out) :: Psaft ! The pressure (in Pa) computed in the frame of the PC-SAFT
real(dp), intent(out) :: Z ! The compressibility factor (no unit)
real(dp), dimension(:), intent(out) :: muskT_k ! The residual chemical potential divided by kT (no unit)
real(dp), intent(out), optional :: ares ! Residual (compare to ideal gas value) Helmoltz energy (J.mol^-1).
real(dp) :: rho ! density in "number of molecules / angstrom^3 (NMA)".
real(dp) :: eta ! reduced density (no unit).
logical :: opt_outputs
integer :: i, j, k, n, ok
real(dp), save :: m_bar
real(dp), dimension(:), allocatable :: d
real(dp), dimension(0:3) :: zeta
real(dp), dimension(:,:), allocatable :: zeta_xk
real(dp), save :: sum_ximidi0, sum_ximidi1, sum_ximidi2, sum_ximidi3
real(dp), dimension(:), allocatable :: gii_hs
real(dp), dimension(:), allocatable :: rho_dgii_hs
real(dp), dimension(:,:), allocatable :: dghs_ii_dxk
real(dp) :: C_1
real(dp), dimension(:,:), allocatable :: sigma, epsilon
real(dp), dimension(0:2,0:6) :: a, b
real(dp), dimension(0:6) :: a_i, b_i
real(dp) :: I_1, I_2, detaI1sdeta, detaI2sdeta
real(dp) :: C_2
real(dp) :: m2_epsilon_sigma3_bar, m2_epsilon2_sigma3_bar
real(dp) :: Z_hs, Z_hc, Z_disp, Zassoc
real(dp) :: molar_mass_bar
real(dp) :: a_res, a_hc, a_hs, a_disp ! reduced Helmholtz energy (no unit)
real(dp) :: sum_xi_mim1_ln_gii, sum_xi_mim1_ghsiim1_dghsiisxk, sum_xj_dares_dxj
real(dp), dimension(:), allocatable :: da_hs_s_dxk, dahc_dxk, dares_dxk
real(dp), dimension(:), allocatable :: m2_epsilon_sigma3_bar_xk, m2_epsilon2_sigma3_bar_xk
real(dp), dimension(:), allocatable :: C_1_xk
real(dp), dimension(:,:), allocatable :: ai_xk, bi_xk
real(dp), dimension(:), allocatable :: I_1_xk, I_2_xk, dadisp_xk
real(dp), dimension(:), allocatable :: muskT_k_assoc
!---------------------------------------------------------------------------------------------------------------------------------
! Call path of this subroutine:
subtemp= sub(1:len_trim(sub))//'/'//subname(1:len_trim(subname))
!---------------------------------------------------------------------------------------------------------------------------------
if ( rho_MKSA <= 0._dp ) then
write(6,*) ''
call wf(' > In module: ', 'bright red', mod_name(1:len_trim(mod_name)), 'bright yellow')
write(6,*) ' > In subroutine ''', subname(1:len_trim(subname)), ''' called by ''', subtemp(1:len_trim(subtemp)), ''''
write(6,*) ' the considered density "rho" is <= 0.!'
write(6,*) ' We have: "rho" = ', rho_MKSA
write(6,*) ''
stop
end if
!---------------------------------------------------------------------------------------------------------------------------------
! We compute the mean number of segment "m_bar" and various quantities used in the following:
m_bar= 0._dp ! cf. Gross & Sadowski (2001)
do i= 1, size(species_prop)
m_bar= m_bar + species_prop(i)%x * species_prop(i)%numseg
end do
! Computation of the d_i's:
allocate(d(1:size(species_prop)),stat=ok); call alloc_error(subname,'d',ok)
do i= 1, size(species_prop)
! Eq. (4.2.9) de Maity (2003) et aussi Eq. [II-136] p. 50 de Belkadi (2008)
d(i)= species_prop(i)%segdiam * (1._dp - 0.12_dp * exp(-3._dp * species_prop(i)%segenerg / T)) ! En Angstrom
!write(6,*) ' d(',i,')= ', d(i)
end do
sum_ximidi0= 0._dp
do i= 1, size(species_prop)
! cf. Eq. (4.2.8) p. 26 Maity (2003) and Eq. [II-131] p. 49 Belkadi (2008)
sum_ximidi0= sum_ximidi0 + species_prop(i)%x * species_prop(i)%numseg
end do
sum_ximidi1= 0._dp
do i= 1, size(species_prop)
! cf. Eq. (4.2.8) p. 26 Maity (2003) and Eq. [II-131] p. 49 Belkadi (2008)
sum_ximidi1= sum_ximidi1 + species_prop(i)%x * species_prop(i)%numseg * d(i)
end do
sum_ximidi2= 0._dp
do i= 1, size(species_prop)
! cf. Eq. (4.2.8) p. 26 Maity (2003) and Eq. [II-131] p. 49 Belkadi (2008)
sum_ximidi2= sum_ximidi2 + species_prop(i)%x * species_prop(i)%numseg * d(i)**2
end do
sum_ximidi3= 0._dp
do i= 1, size(species_prop)
! cf. Eq. (4.2.20) Maity (2003) and Eq. [II-124] p. 48 Belkadi (2008)
sum_ximidi3= sum_ximidi3 + species_prop(i)%x * species_prop(i)%numseg * d(i)**3
end do
!---------------------------------------------------------------------------------------------------------------------------------
! Computation of the density "rho" from "MKSA" to "number of molecules / angstrom^3 (NMA)" :
molar_mass_bar= 0._dp
do i= 1, size(species_prop)
molar_mass_bar= molar_mass_bar + species_prop(i)%x * species_prop(i)%molmass ! Mean molar mass in kg.mol^-1
end do
rho= rho_MKSA / molar_mass_bar * Navo * angstrom_meter**3._dp ! density in "number of molecules / angstrom^3 (NMA)"
! Reduced density (no unit):
eta= rho/6._dp*pi * sum_ximidi3 ! cf. Eq. (4.2.20) Maity (2003) and Eq. [II-124] p. 48 Belkadi (2008)
! unit: number of molecules / angstrom^3 ! V
!---------------------------------------------------------------------------------------------------------------------------------
allocate(gii_hs(1:size(species_prop)),stat=ok); call alloc_error(subname,'gii_hs',ok)
allocate(rho_dgii_hs(1:size(species_prop)),stat=ok); call alloc_error(subname,'rho_dgii_hs',ok)
allocate(sigma(1:size(species_prop),1:size(species_prop)),stat=ok); call alloc_error(subname,'sigma',ok)
allocate(epsilon(1:size(species_prop),1:size(species_prop)),stat=ok); call alloc_error(subname,'epsilon',ok)
allocate(zeta_xk(0:3,1:size(species_prop)),stat=ok); call alloc_error(subname,'zeta_xk',ok)
allocate(da_hs_s_dxk(1:size(species_prop)),stat=ok); call alloc_error(subname,'da_hs_s_dxk',ok)
allocate(dahc_dxk(1:size(species_prop)),stat=ok); call alloc_error(subname,'dahc_dxk',ok)
allocate(dares_dxk(1:size(species_prop)),stat=ok); call alloc_error(subname,'dares_dxk',ok)
allocate(dghs_ii_dxk(1:size(species_prop),1:size(species_prop)),stat=ok); call alloc_error(subname,'dghs_ii_dxk',ok)
allocate(m2_epsilon_sigma3_bar_xk(1:size(species_prop)),stat=ok); call alloc_error(subname,'m2_epsilon_sigma3_bar_xk',ok)
allocate(m2_epsilon2_sigma3_bar_xk(1:size(species_prop)),stat=ok); call alloc_error(subname,'m2_epsilon2_sigma3_bar_xk',ok)
allocate(C_1_xk(1:size(species_prop)),stat=ok); call alloc_error(subname,'C_1_xk',ok)
allocate(ai_xk(0:6,1:size(species_prop)),stat=ok); call alloc_error(subname,'ai_xk',ok)
allocate(bi_xk(0:6,1:size(species_prop)),stat=ok); call alloc_error(subname,'bi_xk',ok)
allocate(I_1_xk(1:size(species_prop)),stat=ok); call alloc_error(subname,'I_1_xk',ok)
allocate(I_2_xk(1:size(species_prop)),stat=ok); call alloc_error(subname,'I_2_xk',ok)
allocate(dadisp_xk(1:size(species_prop)),stat=ok); call alloc_error(subname,'dadisp_xk',ok)
!---------------------------------------------------------------------------------------------------------------------------------
!!! Beginning of the computation of the compressibility factor "Z":
!!!rho= 6.d0/pi * eta / sum_ximidi3 ! cf. Eq. (4.2.20) Maity (2003) and Eq. [II-124] p. 48 Belkadi (2008)
!!! ! unit: number of molecules / angstrom^3 ! V
! Voir Eq. (9) de Gross & Sadowski (2001)
zeta(0)= pi/6._dp * rho * sum_ximidi0 ! Angstrom^-3 ! V
zeta(1)= pi/6._dp * rho * sum_ximidi1 ! Angstrom^-2 ! V
zeta(2)= pi/6._dp * rho * sum_ximidi2 ! Angstrom^-1 ! V
zeta(3)= pi/6._dp * rho * sum_ximidi3 ! no unit ! V
! Eq. (8) de Gross & Sadowski (2001) et ((4.2.26 de MAity (2003) :
Z_hs= zeta(3)/(1._dp-zeta(3)) + &
3._dp*zeta(1)*zeta(2)/zeta(0)/(1._dp-zeta(3))**2 + &
(3._dp*zeta(2)**3-zeta(3)*zeta(2)**3)/zeta(0)/(1._dp-zeta(3))**3 ! V
! We compute the term "g_ii^hs", cf. for instance (4.2.7) p. 26 of Maity (2003):
do i= 1, size(species_prop)
! ATTENTION: cette formule est (4.2.7) p. 26 de Maity (2003) en accord avec [II-129] p. 49 de Belkadi (2008) mais
! dans le papier original de Gross & Sadowski (2001) il n'y a pas de carr au dnominateur du deuxime terme.
! en fait il y a de carr dans la formule (A.7) de G&S.
gii_hs(i)= 1._dp/(1._dp-zeta(3)) + d(i)*d(i)/(d(i)+d(i)) * 3._dp*zeta(2)/(1._dp-zeta(3))**2 + (d(i)*d(i)/(d(i)+d(i)))**2 * &
2._dp*zeta(2)**2/(1._dp-zeta(3))**3 ! V
end do
do i= 1, size(species_prop)
! Voir (4.2.27) p. 29 de Maity (2003) et (A.27) de G&S:
j=i
rho_dgii_hs(i)= zeta(3)/(1._dp-zeta(3))**2 + d(i)*d(j)/(d(i)+d(j)) * &
(3._dp*zeta(2)/(1._dp-zeta(3))**2 + 6._dp*zeta(2)*zeta(3)/(1._dp-zeta(3))**3) + &
(d(i)*d(j)/(d(i)+d(j)))**2 * &
(4._dp*zeta(2)**2/(1._dp-zeta(3))**3 + 6._dp*zeta(2)**2*zeta(3)/(1._dp-zeta(3))**4) ! V
end do
! On calcule Z_hc:
! Eq. (4.2.25) p. 28 de Maity (2003) et aussi Eq. (5) de Gross & Sadowski (2001):
Z_hc= m_bar* Z_hs
do i= 1, size(species_prop)
Z_hc= Z_hc - species_prop(i)%x * (species_prop(i)%numseg-1._dp) /gii_hs(i) * rho_dgii_hs(i) ! V
end do
! We compute Z^disp:
! Eq. (4.2.11) p. 26 de Maity (2003), ATTENTION dans G&S il manque l'exposant <<-1>> !
C_1= 1._dp/(1._dp + m_bar * (8._dp*eta-2._dp*eta**2)/(1._dp-eta)**4 + (1._dp-m_bar) * &
(20._dp*eta-27._dp*eta**2+12._dp*eta**3-2._dp*eta**4)/(1._dp-eta)**2/(2._dp-eta)**2) ! V
do i= 1, size(species_prop)
do j= 1, size(species_prop)
sigma(i,j) = 1._dp/2._dp * (species_prop(i)%segdiam + species_prop(j)%segdiam)
! ATTENTION 1 : on met une valeur fixe pour k_ij --->>> il faudrait faire une table initiale o on a tous les k_ij
! ce que je mets l est destin aux premiers tests.
! ATTENTION 2 : - la formule utlise ici est celle de Maity (2003) (4.2.15) p. 27
! - Belkadi (2008) en propose une autre (cf. p. 38)
! - Gross & Sadowski (2001) Eq. (24) est encore autre chose !
! ATTENTION : les calculs des 'k_ij' la temprature de travail 'T' se fait dans la routine 'pcsaft'
epsilon(i,j)= sqrt(species_prop(i)%segenerg * species_prop(j)%segenerg) * (1._dp - k_ij(i,j))
!print*, ' k_ij(i,j)= ', k_ij(i,j)
end do
end do
call univ_cst(subname,1,a,b)
do i= 0, 6
! Eq. (4.2.18) p. 27 de Maity (2003), aussi Eq. (18) de Gross & Sadowski (2001) et galement [II-145] et [II-146] p. 51 de Belkadi
a_i(i)= a(0,i) + (m_bar-1._dp)/m_bar * a(1,i) + (m_bar-1._dp)/m_bar * (m_bar-2._dp)/m_bar * a(2,i) ! V
b_i(i)= b(0,i) + (m_bar-1._dp)/m_bar * b(1,i) + (m_bar-1._dp)/m_bar * (m_bar-2._dp)/m_bar * b(2,i) ! V
end do
! the I1 and I2:
I_1= 0._dp
I_2= 0._dp
do i= 0, 6
! Eq. (4.2.16) (qui est fausse cf. (4.1.20) p. 21 et (4.2.15) de Maity (2003); et aussi Eq. (16) & (17) de Gross & Sadowski (2001)
! [II-143]et [II-144] p. 51 de Belkadi
I_1= I_1 + a_i(i) * eta**i
I_2= I_2 + b_i(i) * eta**i
end do
detaI1sdeta= 0._dp
detaI2sdeta= 0._dp
do i= 0, 6
! Eq. (4.2.29) & (4.2.30) de Maity (2003)
detaI1sdeta= detaI1sdeta + a_i(i) * (i+1._dp) * eta**i ! V
detaI2sdeta= detaI2sdeta + b_i(i) * (i+1._dp) * eta**i ! V
end do
! Eq. (4.2.31) p. 29 de Maity (2003)
C_2= -C_1**2 * (m_bar * (-4._dp*eta**2+20._dp*eta+8._dp)/(1._dp-eta)**5 + (1._dp-m_bar) * &
(2._dp*eta**3+12._dp*eta**2-48._dp*eta+40._dp) &
/ (1._dp-eta)**3 / (2._dp-eta)**3 ) ! V
! Calcul de <m^2 Epsilon sigma^3> :
m2_epsilon_sigma3_bar= 0._dp
do i= 1, size(species_prop)
do j= 1, size(species_prop)
! Eq. (4.2.12) p. 26 de Maity (2003)
m2_epsilon_sigma3_bar= m2_epsilon_sigma3_bar + species_prop(i)%x * species_prop(j)%x * &
species_prop(i)%numseg * species_prop(j)%numseg * &
epsilon(i,j)/T * sigma(i,j)**3 ! V
end do
end do
! Calcul de <m^2 Epsilon^2 sigma^3> :
m2_epsilon2_sigma3_bar= 0._dp
do i= 1, size(species_prop)
do j= 1, size(species_prop)
! Eq. (4.2.13) p. 27 de Maity (2003)
m2_epsilon2_sigma3_bar= m2_epsilon2_sigma3_bar + species_prop(i)%x * species_prop(j)%x * &
species_prop(i)%numseg * species_prop(j)%numseg * &
epsilon(i,j)**2/T**2 * sigma(i,j)**3 ! V
end do
end do
! On finalise le calcul de Z^disp:
! Eq. (4.2.28) p. 29 de Maity (2003) et (A.28) de G&S :
!print*, ' Dans Zdisp :'
!print*, ' > rho = ', rho
!print*, ' > detaI1sdeta = ', detaI1sdeta
!print*, ' > m2_epsilon_sigma3_bar = ', m2_epsilon_sigma3_bar
!print*, ' > C_1 = ', C_1
!print*, ' > detaI2sdeta = ', detaI2sdeta
!print*, ' > C_2 = ', C_2
!print*, ' > eta = ', eta
!print*, ' > I_2 = ', I_2
!print*, ' > m2_epsilon2_sigma3_bar= ', m2_epsilon2_sigma3_bar
Z_disp= -2._dp*pi*rho*detaI1sdeta*m2_epsilon_sigma3_bar - pi*rho*m_bar*(C_1*detaI2sdeta+C_2*eta*I_2)*m2_epsilon2_sigma3_bar ! V
! Calcul final de Z :
! Eq. (4.2.24) p. 28 de Maity (2003) et (A.24) de G&S :
Z= 1._dp + Z_hc + Z_disp
! End of the computation of the compressibility factor "Z":
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
! Beginning of the computation of the pressure "Psaft":
Psaft= (1._dp + Z_hc + Z_disp) * kB * T * rho / angstrom_meter**3.0_dp ! with rho in "number of molecules / angstrom^3", Psaft in pascal.
! End of the computation of the pressure "Psaft".
!---------------------------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------------------------
! Beginning of the computation of the chemical potential "mu_res/kT" (Eq. 4.2.33 of Maity, 2003):
!
a_hs= 1._dp/zeta(0) * ( 3._dp*zeta(1)*zeta(2)/(1._dp-zeta(3)) + & ! Eq. (A.6) de Gross & Sadowski (2001)
zeta(2)**3/zeta(3)/(1._dp-zeta(3))**2 + & ! et aussi Eq. (4.2.6) de Maity (2003)
(zeta(2)**3/zeta(3)**2-zeta(0)) * log(1._dp-zeta(3)) & !
) ! V !
sum_xi_mim1_ln_gii= 0._dp
do i= 1, size(species_prop)
sum_xi_mim1_ln_gii= sum_xi_mim1_ln_gii + species_prop(i)%x * & !
(species_prop(i)%numseg - 1._dp) * & !
log(gii_hs(i)) ! V
end do
a_hc= m_bar * a_hs - sum_xi_mim1_ln_gii ! cf. Eq. (A.4) p. 1256 de Gross & Sadowski (2001) et aussi Eq. (4.2.4) de Maity (2003) ! V