MoDeNa  1.0
Software framework facilitating sequential multi-scale modelling
gc_method.f90
1 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
2 ! MODULE GC_DATA
3 !
4 ! This module contains the experimental data as well as the calculated
5 ! properties
6 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
7 !
8 MODULE gc_data
9 
10  implicit none
11  save
12 
13  !-----------------------------------------------------------------------------
14  ! storing experimental data
15  !-----------------------------------------------------------------------------
16 
17  INTEGER, PARAMETER :: mmx = 3500
18 
19  INTEGER :: nlv,nliq
20  REAL, DIMENSION(MMX) :: plv, tlv, pliq, tliq, vliq, v_crit1, v_crit2
21  CHARACTER (LEN=30), DIMENSION(MMX) :: gc_comp1, gc_comp2
22 
23  !-----------------------------------------------------------------------------
24  ! storing/transferring calculated properties
25  !-----------------------------------------------------------------------------
26 
27  INTEGER :: n_liq(20), n_lv(20), k_nr
28  REAL, DIMENSION(MMX) :: plvcal, v_cal
29  CHARACTER (LEN=30) :: subst(20)
30 
31 END MODULE gc_data
32 
33 
34 
35 
36 
37 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
38 ! MODULE GC_GROUP
39 !
40 ! This module contains parameters detailing a GC-group
41 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
42 !
43 MODULE gc_group
44 
45  USE basic_variables, ONLY: nc
46  implicit none
47  save
48 
49  !-----------------------------------------------------------------------------
50  ! N_comp: maximum number of components i considered in the GC-fit
51  ! N_grs : maximum number of groups j for a component
52  ! Index j runs j = 1, gc_tot_nr_groups(i) for specific comp. i
53  !-----------------------------------------------------------------------------
54 
55  INTEGER, PARAMETER :: n_comp = nc
56  INTEGER, PARAMETER :: n_grs = 20
57 
58  REAL, DIMENSION(N_comp,N_grs) :: m_gc, msig3, meps, mu_gc, mm_gc
59 
60  INTEGER, DIMENSION(N_grs) :: n_groups
61  INTEGER, DIMENSION(N_comp,N_grs) :: gc_groups
62  INTEGER, DIMENSION(N_comp) :: gc_tot_nr_groups
63  CHARACTER (LEN=20), DIMENSION(N_comp,N_grs) :: gc_grouptype
64  ! CHARACTER (LEN=30) :: comp_name(N_comp)
65 
66 END MODULE gc_group
67 
68 
69 
70 
71 
72 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
73 !
74 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
75 
76 MODULE mod_para_transfer
77  implicit none
78  save
79 
80  REAL, ALLOCATABLE :: para_transfer(:)
81 
82 END MODULE mod_para_transfer
83 
84 
85 
86 
87 
88 
89 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
90 ! SUBROUTINE GC_PARAMETER
91 !
92 ! This subroutine contains the group contribution parameters of the
93 ! individual chemical groups. It also combines the group-parameters of
94 ! a molecule to the molecular PCP-SAFT parameters.
95 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
96 
97 SUBROUTINE gc_parameter ( n_species, comp_name )
98 
99  USE basic_variables
100  USE gc_group
101  USE mod_para_transfer
102  IMPLICIT NONE
103 
104  !-----------------------------------------------------------------------------
105  INTEGER, INTENT(IN) :: n_species
106  CHARACTER (LEN=30), INTENT(IN) :: comp_name(n_comp)
107  !-----------------------------------------------------------------------------
108  INTEGER :: i, j, k, jj, no
109  INTEGER, DIMENSION(nc) :: nhb_typ
110  INTEGER, DIMENSION(nc,nsite) :: nhb_no
111  REAL, DIMENSION(nc,nc,nsite,nsite) :: eps_hb
112  REAL, DIMENSION(nc,nc) :: kap_hb
113  REAL, DIMENSION(N_comp,N_grs) :: epsilon_hb, kappa_hb
114  !-----------------------------------------------------------------------------
115 
116  n_groups = 0
117  gc_grouptype = ' '
118  gc_tot_nr_groups = 0
119 
120  m_gc = 0.0
121  msig3 = 0.0
122  meps = 0.0
123  mu_gc = 0.0
124  mm_gc = 0.0
125  epsilon_hb = 0.0
126  kappa_hb = 0.0
127 
128  !-----------------------------------------------------------------------------
129  ! draw structural groups of molecules from file
130  !-----------------------------------------------------------------------------
131 
132  CALL get_gc_groups ( n_species, comp_name ) !, gc_tot_nr_groups, gc_groups, gc_grouptype)
133  !write (*,*) n_species, comp_name
134 
135  !-----------------------------------------------------------------------------
136  ! group-contribution parameters
137  !-----------------------------------------------------------------------------
138 
139  DO i = 1, n_species
140  DO j = 1, gc_tot_nr_groups(i)
141  IF (gc_grouptype(i,j) == '-CH3') THEN ! CH3-
142  m_gc(i,j) = 0.695765
143  msig3(i,j)= 34.0119
144  meps(i,j) = 147.516
145  mm_gc(i,j) = 15.035
146  ELSE IF (gc_grouptype(i,j) == '-CH2-') THEN ! -CH2-
147  m_gc(i,j) = 0.429404
148  msig3(i,j)= 25.08362
149  meps(i,j) = 108.1502
150  mm_gc(i,j) = 14.02675 ! ?
151  ELSE IF (gc_grouptype(i,j) == '-CH<(branch)') THEN
152  m_gc(i,j) = 0.143
153  msig3(i,j)= 16.41211
154  meps(i,j) = 49.71252
155  mm_gc(i,j) = 13.015 ! ?
156  ELSE IF (gc_grouptype(i,j) == '>CH<(2branch)') THEN
157  m_gc(i,j) = -0.6700
158  msig3(i,j)= 3.82853
159  meps(i,j) = -72.1456
160  mm_gc(i,j) = 12.003 ! ?
161  ELSE IF (gc_grouptype(i,j) == 'CH4') THEN ! methane
162  m_gc(i,j) = 1.000
163  msig3(i,j)= 3.70388767**3
164  meps(i,j) = 150.033987
165  mm_gc(i,j) = 16.043
166  ELSE IF (gc_grouptype(i,j) == '-Cl') THEN
167  m_gc(i,j) = 1.5514 / 2.0
168  msig3(i,j)= 1.5514 / 2.0 * 3.3672**3
169  meps(i,j) = 265.67 * 1.5514 / 2.0
170  mm_gc(i,j) = 35.453
171  ELSE IF (gc_grouptype(i,j) == 'N(N2)') THEN ! nitrogen
172  m_gc(i,j) = 1.205275098 / 2.0
173  msig3(i,j)= 1.205275098 / 2.0 * 3.3129702**3
174  meps(i,j) = 90.9606924 * 1.205275098 / 2.0
175  mm_gc(i,j) = 14.005
176  ELSE IF (gc_grouptype(i,j) == '-CH2-(cyclic6)') THEN ! -CH2- C6-cyclic
177  m_gc(i,j) = 0.4217
178  msig3(i,j)= m_gc(i,j) * 3.84990887**3
179  meps(i,j) = m_gc(i,j) * 278.108786
180  mm_gc(i,j) = 14.02675
181  ELSE IF (gc_grouptype(i,j) == '-CH2-(cyclic5)') THEN ! -CH2- C6-cyclic
182  m_gc(i,j) = 0.473
183  msig3(i,j)= m_gc(i,j) * 3.71139254**3
184  meps(i,j) = m_gc(i,j) * 265.828755
185  mm_gc(i,j) = 14.02675
186  ELSE IF (gc_grouptype(i,j) == '-CH-(aromatic)') THEN ! -CH- C6-aromatic
187  m_gc(i,j) = 0.4108
188  msig3(i,j)= 19.93845
189  meps(i,j) = 118.0598
190  mm_gc(i,j) = 13.015
191  ELSE IF (gc_grouptype(i,j) == '-C-R(aromatic)') THEN ! -C-R C6-aromatic
192  m_gc(i,j) = 0.026654
193  msig3(i,j)= 14.1613
194  meps(i,j) = 52.5285
195  mm_gc(i,j) = 12.015
196  ELSE IF (gc_grouptype(i,j) == '-O-ether') THEN
197  m_gc(i,j) = 0.81453
198  msig3(i,j)= 10.94506
199  meps(i,j) = 162.6272
200  mu_gc(i,j) = 1.7488
201  mm_gc(i,j) = 15.998
202  ELSE IF (gc_grouptype(i,j) == '-OH') THEN
203  m_gc(i,j) = 0.6409
204  msig3(i,j)= m_gc(i,j) * ( 2.7707 )**3
205  meps(i,j) = m_gc(i,j) * 376.99
206  mu_gc(i,j) = 1.7398
207  mm_gc(i,j) = 17.007
208  epsilon_hb(i,j) = 2374.0
209  kappa_hb(i,j) = 0.00698
210  m_gc(i,j) = 1.4633
211  msig3(i,j)= m_gc(i,j) * ( 2.5320 )**3
212  meps(i,j) = m_gc(i,j) * 198.58
213  mu_gc(i,j) = 1.3843
214  mm_gc(i,j) = 17.007
215  epsilon_hb(i,j) = para_transfer(1) !2431.3
216  kappa_hb(i,j) = 0.03256
217  ELSE IF (gc_grouptype(i,j) == 'CO_Acetone') THEN
218  m_gc(i,j) = 1.138
219  msig3(i,j)= 14.0
220  meps(i,j) = 200.0
221  mu_gc(i,j) = 2.88
222  mm_gc(i,j) = 28.010
223  ELSE IF (gc_grouptype(i,j) == '-CH2-COO-CH2-') THEN ! -CH2-COO-CH2- ester-group (from acid-side)
224  m_gc(i,j) = 2.0000
225  msig3(i,j)= 82.684
226  meps(i,j) = 491.00
227  mu_gc(i,j) = 3.674
228  mm_gc(i,j) = 72.06 ! geschaetzt
229  ! m_gc(i,j) = 2.5 + ( para_transfer(1) - 10.0 )
230  ! msig3(i,j)= 100.0 + ( para_transfer(2) - 10.0 )*20.0
231  ! meps(i,j) = 600.0 + ( para_transfer(3) - 10.0 )*100.0
232  ! mu_gc(i,j) = 2.0 + ( para_transfer(4) - 10.0 )
233  ! mm_gc(i,j) = 72.06 ! geschaetzt
234  ELSE IF (gc_grouptype(i,j) == '-CH2-COO-CH3') THEN ! -CH2-COO-CH3 methyl ester-group (from acid-side)
235  m_gc(i,j) = 2.5386
236  msig3(i,j)= 91.500
237  meps(i,j) = 582.58
238  mu_gc(i,j) = 3.415
239  mm_gc(i,j) = 73.07 ! geschaetzt
240  ! write (*,*) m_gc(i,j),msig3(i,j),meps(i,j),mu_gc(i,j)
241  ! stop
242  ELSE
243  WRITE (*,*) 'GC_PARAMETER: GC-parameter not found !', gc_grouptype(i,j)
244  stop
245  END IF
246  END DO
247  END DO
248 
249  !-----------------------------------------------------------------------------
250  ! calculate pure component parameters (i.e. molecular parameters)
251  !-----------------------------------------------------------------------------
252 
253  DO i = 1, n_species
254  parame(i,1) = 0.0
255  parame(i,2) = 0.0
256  parame(i,3) = 0.0
257  parame(i,6) = 0.0
258  mm(i) = 0.0
259  DO j = 1, gc_tot_nr_groups(i)
260  parame(i,1) = parame(i,1) + m_gc(i,j) * REAL(gc_groups(i,j))
261  parame(i,2) = parame(i,2) + msig3(i,j) * REAL(gc_groups(i,j))
262  parame(i,3) = parame(i,3) + meps(i,j) * REAL(gc_groups(i,j))
263  parame(i,6) = parame(i,6) + mu_gc(i,j) * REAL(gc_groups(i,j))
264  mm(i) = mm(i) + mm_gc(i,j) * REAL(gc_groups(i,j))
265  if ( epsilon_hb(i,j) /= 0.0 ) then
266  if ( n_species > 1 ) write (*,*) 'extend GC method'
267  if ( n_species > 1 ) stop
268  eps_hb = 0.0
269  eps_hb(1,1,1,2) = epsilon_hb(i,j)
270  eps_hb(1,1,2,1) = eps_hb(1,1,1,2)
271  kap_hb(i,i)= kappa_hb(i,j)
272  nhb_typ(1) = 2
273  nhb_no(1,1) = 1
274  nhb_no(1,2) = 1
275  IF (nhb_typ(i) /= 0) THEN
276  parame(i,12) = REAL(nhb_typ(i))
277  parame(i,13) = kap_hb(i,i)
278  no = 0
279  DO jj=1,nhb_typ(i)
280  DO k=1,nhb_typ(i)
281  parame(i,(14+no))=eps_hb(i,i,jj,k)
282  no=no+1
283  END DO
284  END DO
285  DO jj=1,nhb_typ(i)
286  parame(i,(14+no))=REAL(nhb_no(i,jj))
287  no=no+1
288  END DO
289  ELSE
290  DO k=12,25
291  parame(i,k)=0.0
292  END DO
293  END IF
294  end if
295  END DO
296  parame(i,2) = ( parame(i,2)/parame(i,1) )**(1.0/3.0)
297  parame(i,3) = parame(i,3)/parame(i,1)
298  END DO
299 
300  DO i = 1, n_species
301  IF (parame(i,1) == 0.0) THEN
302  WRITE (*,*) 'error in GC_PARAMETER ',comp_name(i)
303  stop
304  END IF
305  END DO
306 
307  ! write (*,*) parame(1,1),parame(1,2),parame(1,3),parame(1,6),mm(1)
308  ! pause
309 
310 END SUBROUTINE gc_parameter
311 
312 
313 
314 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
315 ! SUBROUTINE GET_GC_GROUPS
316 !
317 ! This subroutine assignes chemical groups to each species. The species
318 ! is identified by name. The output of the subroutine is:
319 !
320 ! * gc_grouptype : type of a group
321 ! * gc_groups : numer of groups of a certain type
322 ! * gc_tot_nr_groups : total number of groups of a substance i
323 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
324 
325 SUBROUTINE get_gc_groups ( n_species, comp_name )
326 
327  USE gc_group
328  USE basic_variables, ONLY: nc
329  use utilities, only: file_open
330  IMPLICIT NONE
331 
332  !-----------------------------------------------------------------------------
333  INTEGER, INTENT(IN) :: n_species
334  CHARACTER (LEN=30), INTENT(IN) :: comp_name(n_comp)
335 
336  !-----------------------------------------------------------------------------
337  INTEGER :: i, j, k
338  INTEGER :: n_comps, tot_n_groups !, n_groups(N_grs)
339  CHARACTER (LEN=30) :: c_name
340  CHARACTER (LEN=50) :: gc_file
341  CHARACTER (LEN=20) :: dum, grouptype(n_grs)
342  !-----------------------------------------------------------------------------
343 
344  grouptype = ' '
345  gc_groups = 0
346 
347  gc_file = './input_file/gc_method/gc_groups.inp'
348  CALL file_open(gc_file,85)
349  READ(85,*) dum, n_comps
350  READ(85,*) dum
351  READ(85,*) dum
352  DO k = 1, n_comps
353  READ(85,*) c_name, tot_n_groups, n_groups(1), dum, grouptype(1)
354  DO j = 2, tot_n_groups
355  READ(85,*) n_groups(j), dum, grouptype(j)
356  END DO
357  DO i = 1, n_species
358  IF (c_name == comp_name(i) ) THEN
359  gc_tot_nr_groups(i) = tot_n_groups
360  DO j = 1, gc_tot_nr_groups(i)
361  gc_groups(i,j) = n_groups(j)
362  gc_grouptype(i,j) = grouptype(j) ! trim( adjustl( grouptype(j) ))
363  END DO
364  END IF
365  END DO
366  END DO
367  CLOSE (85)
368 
369 END SUBROUTINE get_gc_groups
370 
371 
372 
373 
374 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
375 ! SUBROUTINE GC_RESIDUAL
376 !
377 ! Subroutine for fitting pure component parameters. The residual between
378 ! calculated values and experimental data is calculated and stored in
379 ! vector 'fvec'. The vector 'fvec' serves as the objective function for
380 ! the parameter regression.
381 ! GC_RESIDUAL is called by GC_FITTING via the Levenberg-Marquardt routine.
382 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
383 
384 SUBROUTINE gc_residual (lm_m_dat, lm_n_par, para, fvec, iflag)
385 
386  USE parameters, ONLY: rgas
387  USE basic_variables
388  USE gc_data
389  USE gc_group, ONLY: n_comp
390  USE mod_para_transfer
391  use utilities, only: dens_calc, si_dens
392  IMPLICIT NONE
393 
394  !-----------------------------------------------------------------------------
395  ! INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15, 307)
396  INTEGER, INTENT(IN) :: lm_m_dat
397  INTEGER, INTENT(IN) :: lm_n_par
398  REAL, INTENT(IN) :: para(:)
399  REAL, INTENT(IN OUT) :: fvec(:)
400  INTEGER, INTENT(IN OUT) :: iflag
401 
402  !-----------------------------------------------------------------------------
403  INTEGER :: i, j_dat
404  INTEGER :: converg, crit_dat
405  REAL :: weigh, tc
406  REAL :: toterr, plv_dum, v_dum
407  REAL :: density(np), w(np,nc)
408  REAL :: rho_phas(np)
409  LOGICAL :: scan, vapor, nocon = .false.
410  CHARACTER (LEN=30) :: comp_name(n_comp)
411  CHARACTER (LEN=4) :: char_len
412  !-----------------------------------------------------------------------------
413 
414  ALLOCATE( para_transfer(lm_n_par) )
415  para_transfer = para
416 
417  fvec = 0.0
418 
419  parame(1,5) = 0.12
420 
421  WRITE (char_len,'(I3)') lm_n_par
422  IF (iflag == 1) WRITE (*,'(a,'//char_len//'(1x,G21.14))') ' parameter ', (para(i),i=1,lm_n_par)
423 
424  !-----------------------------------------------------------------------------
425  ! liquid densit data
426  !-----------------------------------------------------------------------------
427  DO i = 1, nliq
428 
429  comp_name(1) = gc_comp1(i)
430  IF ( i == 1 ) THEN
431  crit_dat = 0
432  tc = 0.0
433  CALL gc_parameter ( 1, comp_name )
434  ELSE
435  IF (gc_comp1(i) /= gc_comp1(i-1)) THEN
436  crit_dat = 0
437  tc = 0.0
438  CALL gc_parameter ( 1, comp_name )
439  ENDIF
440  END IF
441 
442 
443  IF ( pliq(i) > 0.0 ) THEN
444 
445  !-----------------------------------------------------------------------
446  ! case 1: calculate rho to given (p,T)
447  !-----------------------------------------------------------------------
448 
449  nocon = .false.
450  scan = .false.
451 
452  j_dat = i
453  nphas = 1
454  t = tliq(i)
455  p = pliq(i)
456  xi(1,1) = 1.0
457  IF (vliq(i) <= v_crit1(i)) THEN
458  vapor = .false.
459  densta(1) = 0.5
460  ! IF (d_kond(1,i).NE.0.0) densta(1) = d_kond(1,i)
461  ELSE
462  vapor = .true.
463  densta(1) = 0.002
464  ! IF (d_kond(1,i).NE.0.0) densta(1) = d_kond(1,i)
465  END IF
466 
467 38 CONTINUE
468  IF (scan) densta(1) = densta(1) + 0.002
469 
470  CALL dens_calc ( rho_phas )
471  CALL si_dens ( density, w )
472  v_cal(i) = 1.0 / density(1)
473 
474  IF ((v_cal(i)-vliq(i))*1.d2/vliq(i) < 5.0 .AND. .NOT. scan) THEN
475  d_kond(1,i) = dense(1)
476  ELSE
477  d_kond(1,i) = 0.0
478  END IF
479 
480 
481  IF ( vapor .AND. (1.0/density(1)) <= (v_crit1(i)/1.1) ) THEN
482  IF ( densta(1) <= 0.228 ) THEN
483  scan = .true.
484  GO TO 38
485  ELSE
486  nocon = .true.
487  END IF
488  END IF
489  IF ( dense(1) > 0.55 ) WRITE (*,*) 'density ',i,dense(1)
490 
491  weigh = 0.5
492  IF (nocon) weigh = weigh / 4.0
493 
494  fvec(j_dat) = weigh * (v_cal(i)-vliq(i))/ &
495  ( sqrt(abs(vliq(i))) * sqrt(abs(v_cal(i))) )
496 
497 
498 
499  ELSE
500 
501  !-----------------------------------------------------------------------
502  ! case 2, when pliq(i) = 0, indicating rho_L at coexistence
503  !-----------------------------------------------------------------------
504 
505  nphas = 2
506 
507  j_dat = i
508 
509  val_init(1) = 0.5
510  val_init(2) = 1.e-8
511  val_init(3) = tliq(i)
512  val_init(4) = 10.0
513  IF (d_kond(1,j_dat) /= 0.0) val_init(1)=d_kond(1,j_dat)
514  IF (d_kond(2,j_dat) /= 0.0) val_init(2)=d_kond(2,j_dat)
515  IF (d_kond(1,j_dat) /= 0.0 .AND. d_kond(2,j_dat) /= 0.0) val_init(4) = plv_kon(j_dat)
516  val_init(5) = 0.0 ! mole fraction lnx(1,1)
517  val_init(6) = 0.0 ! mole fraction lnx(2,1)
518 
519  CALL pure_equilibrium_fit ( crit_dat, tc, tliq(i), pliq(i), converg, v_cal(i), plv_dum )
520 
521  IF (converg == 1) THEN
522  d_kond(1,j_dat) = dense(1)
523  d_kond(2,j_dat) = dense(2)
524  plv_kon(j_dat) = val_conv(4)
525  ELSE
526  d_kond(1,j_dat) = 0.0
527  d_kond(2,j_dat) = 0.0
528  plv_kon(j_dat) = 10.0
529  plv_dum = 1.e5
530  v_cal(i) = vliq(i) * 0.75
531  IF (crit_dat == 1 .AND. tliq(i) > tc ) v_cal(i) = vliq(i)*1.1
532  END IF
533 
534  weigh = 0.5
535  IF ( nocon ) weigh = weigh / 4.0
536 
537  fvec(j_dat) = weigh * (v_cal(i)-vliq(i))/ &
538  (sqrt(abs(vliq(i)))*sqrt(abs(v_cal(i))))
539 
540  END IF
541 
542  END DO
543  !-----------------------------------------------------------------------------
544  ! end: liquid density data
545  !-----------------------------------------------------------------------------
546 
547 
548 
549  !-----------------------------------------------------------------------------
550  ! vapor pressure data
551  !-----------------------------------------------------------------------------
552  DO i = 1, nlv
553 
554  comp_name(1) = gc_comp2(i)
555  IF ( i == 1 ) THEN
556  crit_dat = 0
557  tc = 0.0
558  CALL gc_parameter ( 1, comp_name )
559  ELSE
560  IF(gc_comp2(i) /= gc_comp2(i-1)) THEN
561  crit_dat = 0
562  tc = 0.0
563  CALL gc_parameter ( 1, comp_name )
564  ENDIF
565  END IF
566 
567  nphas = 2
568 
569  j_dat = i + nliq
570 
571  val_init(1) = 0.5
572  val_init(2) = 1.d-8
573  val_init(3) = tlv(i)
574  val_init(4) = plv(i)
575  IF (d_kond(1,j_dat) /= 0.0) val_init(1) = d_kond(1,j_dat)
576  IF (d_kond(2,j_dat) /= 0.0) val_init(2) = d_kond(2,j_dat)
577  IF (d_kond(1,j_dat) /= 0.0 .AND. d_kond(2,j_dat) /= 0.0) val_init(4) = plv_kon(j_dat)
578  val_init(5) = 0.0 ! mole fraction lnx(1,1)
579  val_init(6) = 0.0 ! mole fraction lnx(2,1)
580 
581  CALL pure_equilibrium_fit ( crit_dat, tc, tlv(i), plv(i), converg, v_dum, plvcal(i) )
582 
583  IF (converg == 1) THEN
584  d_kond(1,j_dat) = dense(1)
585  d_kond(2,j_dat) = dense(2)
586  plv_kon(j_dat) = val_conv(4)
587  ELSE
588  d_kond(1,j_dat) = 0.0
589  d_kond(2,j_dat) = 0.0
590  plv_kon(j_dat) = plv(i)
591  IF ( crit_dat == 1 .AND. tc > tlv(i) ) THEN
592  plvcal(i) = plv(i) * 0.75
593  ! IF (i > 2) plvcal(i) = EXP( LOG(plvcal(i-1)) + (LOG(plvcal(i-1))-LOG(plvcal(i-2)))
594  ! /(1.0/tlv(i-1)-1.0/tlv(i-2)) * (1.0/tlv(i)-1.0/tlv(i-1)) )
595  WRITE (*,'(a,i4,4G15.5)') ' loose ! ', i, tlv(i), t, tc, plv(i)
596  ELSE
597  plvcal(i) = plv(i)
598  IF (i > 2) THEN
599  IF (plvcal(i-2) > 0.0 .AND. plvcal(i-1) > 0.0 .AND. &
600  tlv(i-1) /= tlv(i-2)) THEN
601  plvcal(i) = exp( log(plvcal(i-1)) + (log(plvcal(i-1))-log(plvcal(i-2))) &
602  /(1.0/tlv(i-1)-1.0/tlv(i-2)) *(1.0/tlv(i)-1.0/tlv(i-1)) )
603  END IF
604  END IF
605  ! WRITE (*,*) ' supercritical at point ', i, tlv(i), tc
606  END IF
607  END IF
608 
609  fvec(j_dat) = 0.2 * abs( plvcal(i)-plv(i) )/ &
610  ((sqrt(abs(plv(i)))*sqrt(abs(plvcal(i)))))**0.85
611 
612  IF (tc < tlv(i).AND.crit_dat == 1) THEN
613  fvec(j_dat) = fvec(j_dat) + 5.0 * abs(tlv(i)-tc) + 2.0
614  WRITE(*,'(a,i3,1x,a22,2(f15.5))') ' skipped ! ',i,gc_comp2(i), tlv(i),tc
615  END IF
616 
617  END DO
618  !-----------------------------------------------------------------------------
619  ! end: vapor pressure data
620  !-----------------------------------------------------------------------------
621 
622 
623  toterr = 0.0
624  DO i = 1, nlv+nliq
625  toterr = toterr + fvec(i)**2
626  END DO
627  WRITE (*,*) '--- error ---',toterr
628 
629  IF ( j_dat /= lm_m_dat ) WRITE (*,*) 'ERROR in GC_RESIDUAL', j_dat, lm_m_dat
630 
631  DEALLOCATE( para_transfer )
632 
633 END SUBROUTINE gc_residual
634 
635 
636 
637 
638 
639 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
640 ! SUBROUTINE GC_FITTING
641 !
642 ! Driver for fitting group-contribution parameters of the PCP-SAFT model.
643 ! It uses a Levenberg-Marquardt algorithm to minimize deviations of
644 ! calculated properties to experimental data.
645 ! Under construction: two issues are so far implemented very temprarily
646 ! and need upgrading: (1) the components that are considered in the
647 ! regression are specified in the following subroutine (hard wired !!!)
648 ! and (2) the number of adjustable parameters (lm_n_par) is hard-wired in
649 ! this subroutine !!
650 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
651 
652 SUBROUTINE gc_fitting
653 
654  USE basic_variables
655  USE levenberg_marquardt
656  USE gc_data
657  IMPLICIT NONE
658 
659  !-----------------------------------------------------------------------------
660  INTERFACE
661  SUBROUTINE gc_residual(lm_m_dat, lm_n_par, para, fvec, iflag)
662  IMPLICIT NONE
663  INTEGER, INTENT(IN) :: lm_m_dat, lm_n_par
664  REAL, INTENT(IN) :: para(:)
665  REAL, INTENT(IN OUT) :: fvec(:)
666  INTEGER, INTENT(IN OUT) :: iflag
667  END SUBROUTINE gc_residual
668  END INTERFACE
669 
670  !-----------------------------------------------------------------------------
671  INTEGER :: info
672  INTEGER :: lm_n_par ! number of adjustable parameters
673  INTEGER :: lm_m_dat ! number of data points
674  INTEGER, ALLOCATABLE :: ipvt(:)
675  REAL, ALLOCATABLE :: para(:)
676  REAL, ALLOCATABLE :: fvec(:)
677  REAL :: tol = 1.0e-6
678  REAL :: epsfcn
679  REAL :: lm_factor
680 
681  INTEGER :: i, k
682  INTEGER :: n0
683  REAL :: av_dev, max_dev, rms
684  !-----------------------------------------------------------------------------
685 
686 
687  ncomp = 1
688 
689  ! numerical constants for Levenberg-Marquardt algorithm
690  ! ----------------------------------------------------------------------
691  scaling(1) = 1.e2
692  IF ( num == 0 ) epsfcn = 1.0e-6**2 ! sqrt of relat. step size (finite differences)
693  IF ( num == 1 ) epsfcn = 1.0e-5**2 ! sqrt of relat. step size (finite differences)
694  lm_factor = 0.01 ! maximum initial step size (parameter iteration)
695 
696  ! ======================================================================
697  lm_n_par = 4
698  ! ======================================================================
699 
700 
701  ! initialize parameter vectors
702  ! ----------------------------------------------------------------------
703  ALLOCATE( fvec(lm_m_dat) )
704  ALLOCATE( para(lm_n_par) )
705  ALLOCATE( ipvt(lm_n_par) )
706 
707  DO i = 1, lm_n_par
708  para(i) = 10.0
709  END DO
710 
711  ! ----------------------------------------------------------------------
712  ! read input-files
713  ! ----------------------------------------------------------------------
714 
715  CALL gc_datin ( lm_m_dat )
716 
717  ! ----------------------------------------------------------------------
718  ! adjust parameters (Levenberg-Marquardt scheme)
719  ! ----------------------------------------------------------------------
720 
721  CALL lmdif1(gc_residual,lm_m_dat,lm_n_par,para,fvec,tol,epsfcn,lm_factor,info,ipvt)
722 
723 
724 
725  ! ----------------------------------------------------------------------
726  ! writing results to files
727  ! ----------------------------------------------------------------------
728 
729  DO k = 1, k_nr
730 
731  OPEN (23,file='./output_file/gc_method/'//trim( adjustl( subst(k) ))//'.dat')
732  WRITE(23,*) subst(k)
733  WRITE(23,*) '----------------------------------------------------'
734  ! --------------------------------------------------------------------
735  ! writing experimental liquid density data
736  ! --------------------------------------------------------------------
737  WRITE(23,*) ' liquid density data'
738  WRITE(23,*) ' no T/K p/bar rho/(kg/m**3) rho_cal/(kg/m**3) err%'
739  av_dev = 0.0
740  max_dev = 0.0
741  rms = 0.0
742  n0 = sum( n_liq(1:k-1) )
743  DO i = n0+1 , n0+n_liq(k)
744  WRITE(23,'( I3, 2X, 5(2x,G14.6) )') i-n0, tliq(i), &
745  pliq(i)/1.d5,1.0/vliq(i),1.0/v_cal(i), (vliq(i)/v_cal(i)-1.0)*100.0
746  max_dev = max( max_dev, abs(vliq(i)/v_cal(i)-1.0)*100.0 )
747  rms = rms + (vliq(i)/v_cal(i)-1.0)**2
748  av_dev = av_dev + abs(vliq(i)/v_cal(i)-1.0)*100.0
749  END DO
750  av_dev = av_dev / REAL(n_liq(k))
751  rms = sqrt( rms / REAL(n_liq(k)) )*100.0
752  WRITE(23,*) ' '
753  WRITE(23,*) 'Deviation of calculated to exp. density data'
754  WRITE(23,*) ' '
755  WRITE(23, '(3(A, t35, F7.3, A/))') 'Max. Deviation:',max_dev, ' %', &
756  'Average Deviation: ',av_dev, ' %', 'RMS: ',rms, ' %'
757  WRITE(23,*) ' '
758 
759  ! --------------------------------------------------------------------
760  ! writing experimental vapor pressure data
761  ! --------------------------------------------------------------------
762  WRITE(23,*) ' vapor pressure data'
763  WRITE(23,*) ' no T/K psat_exp/bar psat_cal/bar err%'
764  av_dev = 0.0
765  max_dev = 0.0
766  rms = 0.0
767  n0 = sum( n_lv(1:k-1) )
768  DO i = n0+1, n0+n_lv(k)
769  WRITE(23,'( I3, 2X, 4(2x,G14.6) )') i-n0, tlv(i), plv(i)/1.e5, &
770  plvcal(i)/1.e5,(plvcal(i)-plv(i))/plv(i)*100.0
771  max_dev = max( max_dev, abs(plvcal(i)-plv(i))/plv(i)*100.0 )
772  rms = rms + ( (plvcal(i)-plv(i))/plv(i) )**2
773  av_dev = av_dev + abs(plvcal(i)-plv(i))/plv(i) * 100.0
774  END DO
775  av_dev = av_dev / REAL(n_lv(k))
776  rms = sqrt( rms / REAL(n_lv(k)) )*100.0
777  WRITE(23,*) ' '
778  WRITE(23,*) 'Deviation of calculated to exp. vapor pressure data'
779  WRITE(23,*) ' '
780  WRITE(23, '(3(A, t35, F7.3, A/))') 'Max. Deviation:',max_dev, ' %', &
781  'Average Deviation: ',av_dev, ' %', 'RMS: ',rms, ' %'
782  WRITE(23,*) ' '
783 
784 
785  CLOSE(23)
786  END DO
787 
788  DEALLOCATE( fvec, para, ipvt )
789 
790 END SUBROUTINE gc_fitting
791 
792 
793 
794 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
795 ! SUBROUTINE GC_DATIN ( lm_m_dat )
796 !
797 ! Reading experimental data from pure component input files. Temporarily,
798 ! the species to be considered in the regression are specified here !!
799 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
800 
801 SUBROUTINE gc_datin ( lm_m_dat )
802 
803  USE basic_variables
804  USE gc_data
805  IMPLICIT NONE
806 
807  !-----------------------------------------------------------------------------
808  INTEGER, INTENT(OUT) :: lm_m_dat
809  !-----------------------------------------------------------------------------
810  INTEGER :: datnr, i, k, n0
811  REAL :: v_crit, tlv_tmp, plv_tmp
812  CHARACTER (LEN=1) :: info*80
813  !-----------------------------------------------------------------------------
814 
815  eos = 1
816 
817  !$$$ k_nr = 7
818  !$$$ subst(1) = 'propane'
819  !$$$ subst(2) = 'butane'
820  !$$$ subst(3) = 'pentane'
821  !$$$ subst(4) = 'hexane'
822  !$$$ subst(5) = 'heptane'
823  !$$$ subst(6) = 'octane'
824  !$$$ subst(7) = 'decane'
825 
826  !$$$ k_nr = 4
827  !$$$ subst(1) = 'toluene'
828  !$$$ subst(2) = 'ethylbenzene'
829  !$$$ subst(3) = 'npropylbenzene'
830  !$$$ subst(4) = 'nbutylbenzene'
831 
832 !!$k_nr = 8
833 !!$subst(1) = 'methyl-propanoate'
834 !!$subst(2) = 'ethyl-propanoate'
835 !!$subst(3) = 'npropyl-propanoate'
836 !!$! subst(6) = 'nbutyl-propionate'
837 !!$subst(4) = 'methyl-butanoate'
838 !!$subst(5) = 'ethyl-butanoate'
839 !!$subst(6) = 'npropyl-butanoate'
840 !!$subst(7) = 'methyl-decanoate'
841 !!$subst(8) = 'methyl-dodecanoate'
842 
843  k_nr = 1
844  subst(1) = 'dimethyl-ether'
845 
846 
847  nliq = 0
848  nlv = 0
849 
850  DO k = 1, k_nr
851 
852  OPEN (22,file='./input_file/gc_method/pure_comp/' &
853  //trim( adjustl( subst(k) ))//'.dat')
854 
855  !--------------------------------------------------------------------------
856  ! reading experimental pure component data
857  !--------------------------------------------------------------------------
858  DO i=1,16
859  READ(22,'(A80)') info
860  END DO
861  READ(22,*) mm(1), v_crit ! note molecular mass mm(i) is ignored (overwritten later)
862 
863  DO i=1,15
864  READ(22,'(A80)') info
865  END DO
866  READ(22,*) n_liq(k) ,n_lv(k)
867  READ(22,'(A80)') info
868  READ(22,'(A80)') info
869 
870 
871  !--------------------------------------------------------------------------
872  ! reading experimental liquid density data
873  !--------------------------------------------------------------------------
874  n0 = nliq
875  DO i = n0+1 , n0+n_liq(k)
876  READ(22,*) datnr,pliq(i),tliq(i),vliq(i)
877  nliq=i
878  gc_comp1(i) = subst(k)
879  v_crit1(i) = v_crit
880  END DO
881  READ(22,'(A80)') info
882  READ(22,'(A80)') info
883 
884  !--------------------------------------------------------------------------
885  ! reading experimental vapor pressure data
886  !--------------------------------------------------------------------------
887  n0 = nlv
888  DO i = n0+1, n0+n_lv(k)
889  READ(22,*) datnr,plv_tmp,tlv_tmp
890  IF (plv_tmp > 20.0) THEN
891  nlv = nlv + 1
892  tlv(nlv) = tlv_tmp
893  plv(nlv) = plv_tmp
894  gc_comp2(nlv) = subst(k)
895  v_crit2(nlv) = v_crit
896  END IF
897  END DO
898  n_lv(k) = nlv-n0
899 
900  CLOSE(22)
901 
902  END DO
903 
904  lm_m_dat = nliq + nlv
905 
906 END SUBROUTINE gc_datin
907 
parameters
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains constant...
Definition: modules.f90:6
basic_variables
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains paramete...
Definition: modules.f90:29