MoDeNa  1.0
Software framework facilitating sequential multi-scale modelling
main_prog.f90
Go to the documentation of this file.
1 
5 
23 
24 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
30 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
31 
32 PROGRAM pc_saft
33 
34  !-----------------------------------------------------------------------------
35  USE basic_variables
36  USE eos_constants
37  USE eos_variables, ONLY: dd_term, qq_term, dq_term
38  USE kij_fitting
39  USE bubblepoint
40  IMPLICIT NONE
41 
42  INTEGER :: option
43  !-----------------------------------------------------------------------------
44 
45  num = 0 ! (num=0: using analytical derivatives of EOS)
46  ! (num=1: using numerical derivatives of EOS)
47  ! (num=2: White's renormalization group theory)
48 
49  ensemble_flag = 'tp' ! this specifies, whether the eos-subroutines
50  ! are executed for constant 'tp' or for constant 'tv'
51 
52  rgt_variant = 'phase_cell' ! only for calc. including critical point renormalization:
53  ! choose variant as 'phase_cell' or 'isomorphic'
54 
55  dd_term = 'GV' ! dipole-dipole term of Gross & Vrabec (2006)
56  qq_term = 'JG' ! quadrupole-quadrupole term of Gross (2005)
57  dq_term = 'VG' ! dipole-quadrupole term of Vrabec & Gross (2008)
58  IF ( num /= 0 ) CALL set_default_eos_numerical
59 
60  !-----------------------------------------------------------------------------
61  ! EOS constants
62  !-----------------------------------------------------------------------------
63  CALL eos_const ( ap, bp )
64 
65 
66  !Perform bubble point calculation, i.e. determine
67  !the composition of a vapor phase which is in equilibrium
68  !with the liquid phase that was specified in the input file
69  CALL bubblepointcalculation
70 
71 END PROGRAM pc_saft
72 
73 
74 
75 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
76 !
77 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
78 
79 SUBROUTINE set_default_eos_numerical
80 
81  USE eos_numerical_derivatives
82  IMPLICIT NONE
83 
84  !-----------------------------------------------------------------------------
85 
86  ideal_gas = 'no' ! ( yes, no )
87  hard_sphere = 'CSBM' ! ( CSBM, no )
88  chain_term = 'TPT1' ! ( TPT1, HuLiu, no )
89  disp_term = 'PC-SAFT' ! ( PC-SAFT, CK, PT1, PT2, PT_MF, PT_MIX, no )
90  hb_term = 'TPT1_Chap' ! ( TPT1_Chap, no )
91  lc_term = 'no' ! ( MSaupe, OVL, no )
92  branch_term = 'no' ! ( TPT2, no )
93  ii_term = 'no'
94  id_term = 'no'
95 
96  subtract1 = 'no' ! (1PT, 2PT, no)
97  subtract2 = 'no' ! (ITTpolar, no)
98 
99 END SUBROUTINE set_default_eos_numerical
100 
101 
102 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
103 ! SUBROUTINE FLASH
104 !
105 ! This subroutine performes a 2-phase flash calculation for multi-
106 ! component mixtures. The feed-composition needs to be provided in
107 ! the file INPUT.inp. The subroutine START_VAR actually perfoms the
108 ! flash calculation.
109 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
110 
111 SUBROUTINE flash
112 
113  USE basic_variables
114  USE starting_values
115  IMPLICIT NONE
116 
117  !-----------------------------------------------------------------------------
118  INTEGER :: converg
119  !-----------------------------------------------------------------------------
120 
121  CALL read_input
122  nphas = 2
123 
124  CALL start_var (converg) ! gets starting values, sets "val_init"
125 
126  IF (converg == 1) CALL output
127  IF (converg /= 1) write (*,*) ' NO PHASE SPLIT DETECTED'
128 
129 END SUBROUTINE flash
130 
131 
132 
133 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
134 ! SUBROUTINE PURE_COMP
135 !
136 ! This subroutine performes VLE calculations (option 1) for pure
137 ! substances. It also allows the calculation of isotherms and
138 ! isobars of pure components (option 2 and 3, respectively)
139 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
140 
141 SUBROUTINE pure_comp
142 
143  USE basic_variables
144  use utilities
145  IMPLICIT NONE
146 
147  !-----------------------------------------------------------------------------
148  INTEGER :: i, converg, calcopt, state
149  REAL :: steps, end_x, end_p, end_t, start_p, start_t
150  REAL :: lnphi(np,nc), density(np), w(np,nc), t_sav
151  REAL :: tc, pc, rhoc
152  REAL :: pges, zges, b_cal
153  !-----------------------------------------------------------------------------
154  steps = 20.0
155 
156  CALL read_input
157  nphas = 2
158 
159  IF ( ncomp /= 1 ) THEN
160  write (*,*) 'SPECIFY ONLY ONE COMPONENT IN THE INPUT-FILE:'
161  write (*,*) ' ./input_file/INPUT.INP'
162  stop
163  ENDIF
164 
165  WRITE (*,*) ' CALCULATE VLE (1)'
166  WRITE (*,*) ' CALCULATE ISOTHERMS (2)'
167  WRITE (*,*) ' CALCULATE ISOBARS (3)'
168  WRITE (*,*) ' CALCULATE VIRIAL COEFFICIENTS (4)'
169  WRITE (*,*) ' CRITICAL POINT (5)'
170  READ (*,*) calcopt
171 
172  IF (calcopt == 1) THEN
173 
174  n_unkw = ncomp ! number of quantities to be iterated
175  it(1) = 'p' ! iteration of p
176  running = 't' ! Temp. is running variable in PHASE_EQUILIB
177 
178  t_sav=t
179 
180  val_init(1) = 0.5 ! starting density for a liquid phase
181  val_init(2) = 1.e-6 ! starting density for a vapor phase
182  val_init(3) = t_sav
183  val_init(4) = 1.e-1 ! default starting value for P: 0.1 Pa
184  IF (eos >= 4) val_init(4) = 1.e-3 ! default starting value for LJ-model
185  val_init(5) = 0.0 ! logar. mole fraction: lnx=0, x=1, phase 1
186  val_init(6) = 0.0 ! logar. mole fraction: lnx=0, x=1, phase 2
187 
188  write (*,*) ' CHOOSE END TEMPERATURE [units as in input-file]:'
189  READ (*,*) end_x
190  end_x = end_x + u_in_t
191 
192  outp = 1 ! output to terminal
193  CALL phase_equilib(end_x,steps,converg)
194  tc = t
195  pc = p
196  IF (eos < 4) CALL critical ( tc, pc, rhoc )
197  ! IF (eos >= 4 .AND. eos /= 7) CALL LJ_CRITICAL (tc,pc,rhoc)
198  ! IF (eos >= 4 .AND. eos /= 7) CALL eLJ_CRITICAL (tc,pc,rhoc)
199  ! IF (eos >= 4 .AND. eos /= 7) write(*,*) ' eLJ routine active !!!'
200  ! IF (eos >= 4 .AND. eos /= 7) write(*,*) ' eLJ routine active !!!'
201  ! IF (eos == 7) CALL SW_CRITICAL (tc,pc,rhoc)
202  WRITE(*,'(a,3(f16.5))') 'tc, pc, rhoc',tc-u_out_t, pc/u_out_p, rhoc
203  WRITE (40,'(5(2x,f18.8))') tc-u_out_t, pc/u_out_p, rhoc, rhoc
204  WRITE (40,*) ' '
205 
206  ELSEIF ( calcopt == 2 .OR. calcopt == 3 ) THEN
207 
208  nphas = 1
209  xi(1,1) = 1.0
210 
211  end_p = p
212  end_t = t
213  start_p = p
214  start_t = t
215  IF ( calcopt == 2 ) THEN
216  WRITE(*,*) ' SPECIFY END-PRESSURE [units as in input-file]'
217  READ (*,*) end_p
218  end_p = end_p*u_in_p
219  ELSE
220  WRITE(*,*) ' SPECIFY END-TEMPERATURE [units as input-file]'
221  READ (*,*) end_t
222  end_t = end_t+u_in_t
223  ENDIF
224  WRITE (*,*) ' SPECIFY LIQUID (1) or VAPOR (2) state'
225  READ (*,*) state
226  densta(1) = 0.5 ! starting density for a liquid phase
227  IF (state == 2) densta(1) = 1.e-5 ! density for a vapor phase
228 
229  DO i = 0, int(steps)
230  p = start_p + (end_p-start_p)*REAL(i)/steps
231  t = start_t + (end_t-start_t)*REAL(i)/steps
232  CALL fugacity (lnphi)
233  IF (num /= 2) CALL enthalpy_etc
234  CALL si_dens (density,w)
235  WRITE (*,'(8(2x,g16.9))') t-u_out_t,p/u_out_p,dense(1),density(1),cpres(1),enthal(1), &
236  gibbs(1), speed_of_sound(1)
237  WRITE(40,'(7(2x,g16.9))') t-u_out_t,p/u_out_p,density(1),cpres(1),enthal(1),gibbs(1), &
238  speed_of_sound(1)
239  ENDDO
240 
241  ELSEIF (calcopt == 4) THEN
242 
243  nphas = 1
244  xi(1,1) = 1.0
245  dense(1)= 1.e-7
246 
247  start_t = t
248  WRITE(*,*) ' SPECIFY END-TEMPERATURE [units as input-file]'
249  READ (*,*) end_t
250  end_t = end_t + u_in_t
251 
252  WRITE (*,*) ' T ,B / cm**3/mol'
253  WRITE(40,*) ' T ,B / cm**3/mol'
254 
255  DO i = 0, int(steps)
256  t = start_t + (end_t-start_t)*REAL(i)/steps
257  CALL p_calc (pges,zges)
258  CALL si_dens (density,w)
259 
260  b_cal = (zges-1.0)/density(1)*1000.0*mm(1)
261  WRITE (*,'(2(2x,f14.8))') t-u_out_t,b_cal
262  WRITE(40,'(2(2x,f14.8))') t-u_out_t,b_cal
263  ENDDO
264 
265  ELSE
266  tc = t
267  pc = p
268  IF ( eos < 4) CALL critical (tc,pc,rhoc)
269  WRITE (*,'(a,3(f16.5))') 'tc,pc,rhoc',tc-u_out_t, pc/u_out_p,rhoc
270  WRITE (40,'(5(2x,f18.8))') tc-u_out_t, pc/u_out_p, rhoc, rhoc
271  WRITE (40,*) ' '
272  ENDIF
273 
274 END SUBROUTINE pure_comp
275 
276 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
277 !
278 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
279 
280 SUBROUTINE temp
281 
282  USE basic_variables
283  USE eos_variables, ONLY: eta, rho, pges
284  IMPLICIT NONE
285 
286  !-----------------------------------------------------------------------------
287  INTEGER :: i
288  REAL :: lnphi(np,nc)
289  !-----------------------------------------------------------------------------
290  num = 1
291  CALL set_default_eos_numerical
292 
293  CALL read_input
294  nphas = 1
295  ensemble_flag = 'tv'
296  xi(1,1) = 1.0
297 
298  DO i = 1, 100
299  !densta(1) = 0.45
300  densta(1) = REAL(101-i) / 100.0 * 0.45
301  dense(1) = densta(1)
302  !p = 101.E5 - REAL(i) / 100.0 * 100.E5
303 
304  CALL fugacity (lnphi)
305  write (*,'(3G18.8)') lnphi(1,1)+log(rho), pges/1.e5, eta
306  !CALL ENTHALPY_ETC
307  END DO
308 
309 END SUBROUTINE temp
310 
311 
312 
313 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
314 ! SUBROUTINE RGT_FIT_PARAMETERS
315 !
316 ! ......
317 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
318 
319 SUBROUTINE rgt_fit_parameters
320 
321  USE basic_variables
322  USE fitting_rgt_parameters
323  IMPLICIT NONE
324 
325  INTERFACE
326  SUBROUTINE crit_dev ( fmin, optpara, n )
327  INTEGER, INTENT(IN) :: n
328  REAL, INTENT(IN) :: optpara(:)
329  REAL, INTENT(IN OUT) :: fmin
330  END SUBROUTINE crit_dev
331  END INTERFACE
332 
333  !-----------------------------------------------------------------------------
334  INTEGER :: nadjust, PRIN
335  REAL :: fmin, t0, h0, MACHEP
336  REAL, ALLOCATABLE :: optpara(:)
337 
338  ! INTEGER :: IERROR !,LW,PRIN
339  ! EXTERNAL CRIT_DEV,CRIT_DEV2 !,func ! deactivated func when making a transition to f90
340  ! REAL :: CRIT_DEV,CRIT_DEV2 !,WVEC(14*3),GVEC(3)
341  REAL :: pp(3) !,func ! ,pam(4,3),yam(4),xi(3,3)
342  !-----------------------------------------------------------------------------
343 
344  num = 2
345  CALL read_input
346 
347  IF (ncomp /= 1) THEN
348  write(*,*)'SPECIFY ONLY ONE COMPONENT IN THE INPUT-FILE'
349  stop
350  ENDIF
351 
352  nadjust = 3
353  ALLOCATE( optpara(nadjust) )
354 
355  critfit = 1
356 
357  tcf = 568.9 ! octane
358  pcf = 24.9d5
359  rcf = 232.0
360  tcf = 373.4 ! h2s
361  pcf = 89.70d5
362  rcf = 388.5
363  tcf = 304.21 ! co2
364  pcf = 73.825d5
365  rcf = 466.01
366  tcf = 416.25 ! chloromethane
367  pcf = 66.9d5
368  rcf = 363.0
369  tcf = 386.5 ! 11-difluoro-ethane (R152a)
370  pcf = 45.198d5
371  rcf = 367.9
372  tcf = 190.555 ! methane
373  pcf = 45.95d5
374  rcf = 162.2
375  ! tcf = 512.64 ! methanol
376  ! pcf = 81.35d5
377  ! rcf = 272.0
378  tcf = 425.1 ! butane
379  pcf = 38.e5
380  rcf = 3.92*58.123
381  tcf = 617.8 ! decane
382  pcf = 21.1d5
383  rcf = 232.
384  ! tcf = 369.8 ! propane
385  ! pcf = 42.5d5
386  ! rcf = 225.
387  tcf = 126.19 ! n2
388  pcf = 33.978d5
389  rcf = 312.0
390  tcf = 508.15 ! acetone
391  pcf = 47.62d5
392  rcf = 269.0
393 
394  ! optpara(1) = (2.53 -1.0 )
395  ! optpara(2) = (40.06 -20.0) /20.0
396  ! optpara(3) = (0.165 -0.12 ) *10.0
397 !!$ optpara(1) = (1.6 -1.0 ) *5.0
398 !!$ optpara(2) = (18.0 -10.0) /20.0
399 !!$ optpara(3) = (0.55 -0.4 ) *10.0
400 
401  optpara(1) = (2.1 -0.0 ) *1.0
402  optpara(2) = (35.0 -0.0 ) /20.0
403  optpara(3) = (0.23 -0.0 ) *5.0
404 
405 
406  ! t0 = 1.E-4
407  ! optpara(1) = 2.33 +0.2
408  ! optpara(2) = 30.06 /20.0 +0.2
409  ! optpara(3) = 0.155 *5.0
410  ! pam(1,1) = optpara(1)
411  ! pam(1,2) = optpara(2)
412  ! pam(1,3) = optpara(3)
413  ! yam(1) = CRIT_DEV (optpara,nadjust)
414  ! optpara(1) = 2.33 -0.2
415  ! optpara(2) = 30.06 /20.0 -0.2
416  ! optpara(3) = 0.155 *5.0
417  ! pam(2,1) = optpara(1)
418  ! pam(2,2) = optpara(2)
419  ! pam(2,3) = optpara(3)
420  ! yam(2) = CRIT_DEV (optpara,nadjust)
421  ! optpara(1) = 2.33 +0.0
422  ! optpara(2) = 30.06 /20.0 -0.0
423  ! optpara(3) = 0.155 *5.0 +0.1
424  ! pam(3,1) = optpara(1)
425  ! pam(3,2) = optpara(2)
426  ! pam(3,3) = optpara(3)
427  ! yam(3) = CRIT_DEV (optpara,nadjust)
428  ! optpara(1) = 2.33 -0.0
429  ! optpara(2) = 30.06 /20.0 -0.0
430  ! optpara(3) = 0.155 *5.0 -0.1
431  ! pam(4,1) = optpara(1)
432  ! pam(4,2) = optpara(2)
433  ! pam(4,3) = optpara(3)
434  ! yam(4) = CRIT_DEV (optpara,nadjust)
435  ! CALL amoeba(pam,yam,4,3,nadjust,t0,CRIT_DEV,IERROR)
436 
437 
438  ! pp(1) = optpara(1)
439  ! pp(2) = optpara(2)
440  ! pp(3) = optpara(3)
441  ! xi(1,1)=1.0
442  ! xi(2,1)=1.0
443  ! xi(1,2)=-1.0
444  ! xi(2,2)=1.0
445  ! xi(3,3)=1.0
446  ! t0 = 1.E-4
447  ! CALL powell(pp,xi,nadjust,nadjust,t0,IERROR,fmin)
448  ! write (*,*) 'done'
449  ! fmin=func (pp)
450 
451  pp(1) = optpara(1)
452  pp(2) = optpara(2)
453  pp(3) = optpara(3)
454  t0 = 1.e-4
455  ! CALL frprmn(pp,nadjust,t0,IERROR,fmin) ! deactivated this line when making a transition to f90
456  write (*,*) 'done'
457  ! fmin=func (pp) ! deactivated this line when making a transition to f90
458 
459 
460  t0 = 1.e-4
461  h0 = 0.4
462  prin = 0
463  machep = 1.e-12
464  !
465  call praxis( t0, machep, h0, nadjust, prin, optpara, crit_dev, fmin )
466  !call CRIT_DEV ( fmin, optpara, nadjust )
467  !
468  !! fmin = 1.E-3
469  !! LW = 14*nadjust
470  !! CALL TN (IERROR,nadjust,optpara,fmin,GVEC,WVEC,LW,CRIT_DEV2)
471 
472 
473  WRITE (*,'(a,1(f16.8))') 'deviation %',fmin*100.0
474  WRITE (*,*) optpara(1),optpara(2)*20.0,optpara(3)/5.0
475  WRITE (*,*) ' '
476  DEALLOCATE( optpara )
477 
478 END SUBROUTINE rgt_fit_parameters
479 
480 
481 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
482 !
483 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
484 
485 SUBROUTINE crit_dev ( fmin, optpara, nadjust )
486 
487  USE basic_variables
488  USE fitting_rgt_parameters
489  use utilities, only: critical
490  IMPLICIT NONE
491 
492  !-----------------------------------------------------------------------------
493  INTEGER, INTENT(IN) :: nadjust
494  REAL, INTENT(IN) :: optpara(nadjust)
495  REAL, INTENT(IN OUT) :: fmin
496  !-----------------------------------------------------------------------------
497  REAL :: tc,pc,rhoc
498  !-----------------------------------------------------------------------------
499 
500  llfit = optpara(1)
501  phifit = optpara(2) *20.0
502  chapfit = optpara(3) /5.0
503 
504  tc = tcf + 40.0
505  pc = pcf
506  rhoc = rcf
507 
508  WRITE (*,'(3(f18.12))') optpara(1),optpara(2)*20.0,optpara(3)/5.0
509  CALL critical ( tc, pc, rhoc )
510 
511  fmin = 100.0 * abs(1.0-tc/tcf )**2 &
512  + abs(1.0-pc/pcf )**2 &
513  + 1.e-1 * abs(1.0-rhoc/rcf)**2
514 
515  ! fmin = 10.0*ABS(1.0-optpara(1) )**2 + ABS(2.0-optpara(2) )**1.5 &
516  ! +1.E-1*ABS(3.0-optpara(3))**2.5
517 
518  WRITE (*,'(a,3(f16.5))') 'tc,pc,rhoc',tc-u_out_t, pc/u_out_p,rhoc
519  WRITE (*,'(a,1(f16.8))') 'deviation %',fmin*100.0
520 
521  WRITE (40,'(5(2x,f18.6))') tc-u_out_t, pc/u_out_p, rhoc, rhoc
522  WRITE (40,*) fmin,optpara(1),optpara(2)*20.0,optpara(3)/5.0
523 
524 END SUBROUTINE crit_dev
525 
526 
527 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
528 !
529 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
530 
531 SUBROUTINE lc_comp
532 
533  USE basic_variables
534  USE eos_variables, ONLY: alpha_nematic
535  USE eos_numerical_derivatives
536  IMPLICIT NONE
537 
538  !-----------------------------------------------------------------------------
539  INTEGER :: converg, ncompsave, phaseindex
540  REAL :: steps, end_x, startcon
541  REAL :: p_sav, t_sav
542  !-----------------------------------------------------------------------------
543  steps = 1.0
544  num = 1
545 
546  !-----------------------------------------------------------------------------
547  ! define the Helmholtz energy contributions
548  !-----------------------------------------------------------------------------
549  CALL set_default_eos_numerical
550  chain_term = 'HuLiu' ! ( TPT1, HuLiu, no )
551  disp_term = 'NO' ! ( PC-SAFT, CK, no )
552  hb_term = 'TPT1_Chap' ! ( TPT1_Chap, no )
553  lc_term = 'OVL' ! ( MSaupe, OVL, no )
554 
555  !-----------------------------------------------------------------------------
556  ! read input and start by re-setting the number of sustances =1 (only the pure LC)
557  !-----------------------------------------------------------------------------
558  CALL read_input
559  nphas = 2
560  ncompsave = ncomp
561  ncomp = 1
562 
563  !-----------------------------------------------------------------------------
564  ! for the Mayer-Saupe theory, this used to be it(1)='t' and running='p'
565  !-----------------------------------------------------------------------------
566  n_unkw = ncomp ! number of quantities to be iterated
567  it(1) = 'p' ! iteration of pressure
568  running = 't' ! Temp. is running variable in PHASE_EQUILIB
569 
570  t_sav = t
571  p_sav = p
572 
573  !-----------------------------------------------------------------------------
574  ! default starting values
575  !-----------------------------------------------------------------------------
576  val_init(1) = 0.35
577  val_init(2) = 0.35
578  val_init(3) = t_sav
579  val_init(4) = p_sav
580  val_init(5) = 0.0 ! logar. mole fraction: lnx=0, x=1, phase 1
581  val_init(6) = 0.0 ! logar. mole fraction: lnx=0, x=1, phase 2
582 
583  !-----------------------------------------------------------------------------
584  ! starting value for the OVL-order parameter
585  !-----------------------------------------------------------------------------
586  alpha_nematic = 10.0 ! starting value for alpha of OLV-Theory
587 
588  !-----------------------------------------------------------------------------
589  ! for diagrams: define end-point
590  !-----------------------------------------------------------------------------
591  end_x = t
592  IF (ncompsave == 1) THEN
593  !write(*,*)' CHOOSE END PRESSURE [units as in input-file]:'
594  !READ (*,*) end_x
595  !end_x = end_x*u_in_P
596  ENDIF
597 
598  !-----------------------------------------------------------------------------
599  ! perform phase equilibrium calculation
600  !-----------------------------------------------------------------------------
601  outp = 1 ! output to terminal
602  CALL phase_equilib(end_x,steps,converg)
603 
604  write (*,*) 'pressure=',p/1.e6,'MPa'
605  write (*,*) 'eta =',dense(1),dense(2)
606  write (*,*) 'alpha =',alpha_nematic
607 
608 
609  !-----------------------------------------------------------------------------
610  ! for mixtures: now extend the number of species
611  !-----------------------------------------------------------------------------
612  IF (ncompsave > 1) THEN
613 
614  ncomp = ncompsave
615  n_unkw = ncomp ! number of quantities to be iterated
616  startcon = log(0.000001)
617  write(*,*) ' specify desired end-composition of:'
618  write(*,*) ' component 2 (molefraction)'
619  write(*,*) ' component 2 is:',compna(2)
620  read (*,*) end_x
621  ! end_x = LOG(0.01)
622  ! end_x = 0.14861586
623  ! end_x = 0.05489185
624 
625  val_init(5) = log(1.0 - exp(startcon - 0.0))
626  val_init(6) = startcon - 0.0
627  val_init(7) = log(1.0 - exp(startcon))
628  val_init(8) = startcon
629 
630  write(*,*) ' Is this composition specified for the LC-phase (1)'
631  write(*,*) ' or for the isotropic phase (2)'
632  read (*,*) phaseindex
633 
634  IF (phaseindex == 1) THEN
635  it(2)='x22' ! iteration of mol fraction of comp.1 phase 2
636  running = 'x12' ! lnx(1,1) is running var. in PHASE_EQUILIB
637  ELSE
638  it(2)='x12' ! iteration of mol fraction of comp.1 phase 2
639  running = 'x22' ! lnx(1,1) is running var. in PHASE_EQUILIB
640  ENDIF
641  it(1)='t' ! iteration of temperature
642  sum_rel(1)='x11' ! summation relation: x12 = 1 - sum(x1j)
643  sum_rel(2)='x21' ! summation relation: x22 = 1 - sum(x2j)
644  outp = 1 ! output to terminal
645  steps=10.0
646 
647 
648  CALL phase_equilib(end_x,steps,converg)
649 
650  write(*,*)' '
651  ! write(*,*) ' SPECIFY TEMPERATURE (1)'
652  ! write(*,*) ' OR PRESSURE (2)'
653  ! READ (*,*) iso_x
654  ! IF (iso_x == 1) THEN
655  write(*,*) ' CHOOSE END TEMPERATURE [units as in input-file]:'
656  READ (*,*) end_x
657  end_x = end_x + u_in_t
658  it(1)='p' ! iteration of pressure
659  running='t' ! Temperature is running var. in PHASE_EQUILIB
660 
661  CALL phase_equilib(end_x,steps,converg)
662 
663  write (77,*) val_conv
664  ENDIF
665 
666 END SUBROUTINE lc_comp
667 
668 
669 
670 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
671 ! SUBROUTINE LC_3PHASE
672 !
673 ! This subroutine performes VLLE calculations for binary mixtures
674 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
675 
676 SUBROUTINE lc_3phase
677 
678  USE basic_variables
679  use utilities, only: si_dens
680  IMPLICIT NONE
681 
682  !-----------------------------------------------------------------------------
683  INTEGER :: i, converg
684  REAL :: density(np), w(np,nc)
685  !-----------------------------------------------------------------------------
686 
687  CALL read_input
688 
689  write(*,*) 'we assume the phase 1 to have a target composition !!'
690 
691  nphas = 3
692  n_unkw = 4 ! number of quantities to be iterated
693  it(1) = 'p ' ! iteration of pressure
694  it(2) = 'x11' ! iteration of temperature
695  it(3) = 'x21' ! iteration of mol fraction of comp.1 phase 2
696  it(4) = 'x31' ! iteration of mol fraction of comp.1 phase 3
697  sum_rel(1) = 'x12' ! summation relation: x22 = 1 - sum(x2j)
698  sum_rel(2) = 'x22' ! summation relation: x22 = 1 - sum(x2j)
699  sum_rel(3) = 'x32' ! summation relation: x32 = 1 - sum(x3j)
700 
701  val_init(0) = 1.e-8 ! set starting value for vapor density (3rd phase)
702  val_init(1) = 0.436817011049269
703  val_init(2) = 0.435497920638304
704  val_init(3) = 308.850000000000
705  val_init(4) = 1858758.08530477
706  val_init(5) = -0.160891853903475
707  val_init(6) = -1.90639042291856
708  val_init(7) = -0.173218595064226
709  val_init(8) = -1.83856034172500
710 
711  ! die folgende Definition ist eine voruebergehende Kruecke:
712  xi(3,1) = 0.0001 ! set starting value for 3rd phase - temporary
713  xi(3,2) = 1.0 - xi(3,1)
714  lnx(3,1) = log(xi(3,1))
715  lnx(3,2) = log(xi(3,2))
716  val_init(9) = lnx(3,1)
717  val_init(10)= lnx(3,2)
718 
719  loop_condition: DO
720 
721  CALL objective_ctrl (converg)
722 
723  IF (converg == 1) THEN
724  CALL si_dens (density,w)
725  write(*,*) ' '
726  write(*,*) ' 3-phase equilibrium calc. successful'
727  write(*,*) ' '
728  write(*,*) ' T, P'
729  write(*,*) val_conv(3)-u_out_t,val_conv(4)/u_out_p
730  write(*,*) ' x11 x21 x31'
731  write(*,*) exp(val_conv(5)),exp(val_conv(7)),exp(val_conv(9))
732  write(*,*) exp(val_conv(6)),exp(val_conv(8)),exp(val_conv(10))
733  write(*,*) ' eta1 , eta2 , eta3'
734  write(*,*) dense(1),dense(2),dense(3)
735  write(*,*) ' '
736  write(*,*) ' output written to file: ./output_file/3-phase.xlo'
737 
738  OPEN (67,file='./output_file/3-phase.xlo')
739  write(67,*) ' T, P'
740  write(67,*) val_conv(3)-u_out_t,val_conv(4)/u_out_p
741  write(67,*) ' x11 x21 x31'
742  write(67,*)exp(val_conv(5)),exp(val_conv(7)),exp(val_conv(9))
743  write(67,*)exp(val_conv(6)),exp(val_conv(8)),exp(val_conv(10))
744  write(67,*) ' rho1[kg/m3] rho2[kg/m3] rho3[kg/m3]'
745  write(67,*) density(1),density(2),density(3)
746  write(67,'(6e16.8)') val_conv(3), val_conv(4)/u_out_p, &
747  1.0-exp(val_conv(5)),exp(val_conv(7)),exp(val_conv(9))
748  ENDIF
749 
750  IF (t.lt.309.8) THEN
751  DO i=1,(4+ncomp*nphas)
752  val_init(i)=val_conv(i)
753  ENDDO
754  val_init(3) = val_init(3) + 0.2
755  t = val_init(3)
756  ELSE
757  EXIT loop_condition
758  ENDIF
759  END DO loop_condition
760 
761 END SUBROUTINE lc_3phase
762 
763 
764 
765 
766 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
767 ! SUBROUTINE JT_INVERSION
768 !
769 ! This sub. calculates Joule-Thomson curves for pure substances
770 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
771 
772 SUBROUTINE jt_inversion
773 
774  USE parameters, ONLY: rgas
775  USE basic_variables
776  use utilities, only: dens_calc, si_dens
777  IMPLICIT NONE
778 
779  !-----------------------------------------------------------------------------
780  INTEGER :: state
781  REAL :: zdT = 10.0
782  REAL :: zdT2, zges1, zges2, zges3, start_t, delta
783  REAL :: density(np), w(np,nc)
784  REAL :: rho_phas(np)
785  !-----------------------------------------------------------------------------
786 
787  CALL read_input
788  nphas = 1
789 
790  IF (ncomp /= 1) THEN
791  write(*,*)'SPECIFY ONLY ONE COMPONENT IN THE INPUT-FILE:'
792  write(*,*)' ./input_file/INPUT.INP'
793  stop
794  ENDIF
795 
796  xi(1,1) = 1.0
797  delta = 1.e-1
798  start_t = t
799  WRITE (*,*) ' SPECIFY LIQUID (1) or VAPOR (2) state'
800  READ (*,*) state
801  densta(1) = 0.5
802  IF (state == 2) densta(1) = 1.e-5
803 
804  DO WHILE ( abs(zdt) > 1.e-10 )
805 
806  t = start_t -delta
807  CALL dens_calc(rho_phas)
808  CALL si_dens (density,w)
809  zges1 = p/(density(1)*1000.0/mm(1)*rgas*t)
810  t = start_t +delta
811  CALL dens_calc(rho_phas)
812  CALL si_dens (density,w)
813  zges3 = p/(density(1)*1000.0/mm(1)*rgas*t)
814  t = start_t
815  CALL dens_calc(rho_phas)
816  CALL si_dens (density,w)
817  zges2 = p/(density(1)*1000.0/mm(1)*rgas*t)
818  zdt = (zges3-zges1)/(2.0*delta)
819  zdt2 = (zges3-2.0*zges2+zges1)/delta**2
820  start_t = t - zdt / zdt2
821  write (*,*) start_t ,p, zdt
822 
823  END DO
824  write (40,*) 'T/K p/kPa error=dz_dT'
825  write (40,*) start_t ,p/1000.0, zdt
826  write (*,*) 'result written to: output_file/output.xlo'
827 
828 END SUBROUTINE jt_inversion
829 
830 
831 
832 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
833 ! SUBROUTINE BINMIX
834 !
835 ! This subroutine performes VLE and LLE calculations for binary
836 ! mixtures. Allowed are regular or polymeric substances. The routine
837 ! is designed to derive phase equilibrium diagrams of binary systems.
838 ! Three calculation options are available:
839 ! option 1 (ISOTHERM): gives Pxy-diagram outputs. The pressure is
840 ! varied in 20 steps from the starting-pressure to an end-pressure
841 ! defined through by the user through the terminal. The step-size
842 ! in pressure is adjusted if the calculation approaches the
843 ! critical point (or an azeotropic point)
844 ! option 2 (ISOBAR): gives Txy-diagram outputs. Similar to option 1,
845 ! whereby the temperature is varied.
846 ! option 3 (MOLE FRACTION SCAN): Either, Txy- or Pxy-diagram outputs
847 ! can be selected this is a suitable option for azeotropic
848 ! mixtures.
849 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
850 
851 SUBROUTINE binmix
852 
853  USE basic_variables
854  USE starting_values, only: start_var
855  use utilities, only: select_sum_rel
856  IMPLICIT NONE
857 
858  !-----------------------------------------------------------------------------
859  INTEGER :: iso_x, converg
860  REAL :: steps, end_x, pinit
861  LOGICAL :: second
862  !-----------------------------------------------------------------------------
863 
864  steps = 20.0
865  second = .false.
866  bindiag = 1
867 
868  CALL read_input
869  nphas = 2
870  pinit = p
871 
872  CALL start_var(converg) ! gets starting values, sets "val_init"
873 
874  n_unkw = ncomp ! number of quantities to be iterated
875  it(2) = 'x21' ! iteration of mol fraction of comp.1 phase 2
876  sum_rel(1) = 'x12' ! summation relation: x12 = 1 - sum(x1j)
877  sum_rel(2) = 'x22' ! summation relation: x22 = 1 - sum(x2j)
878 
879 
880  IF (converg == 1) THEN
881  write (*,*) ' '
882  write(*,*)' CHOOSE YOUR CALCULATION OPTION:'
883  write(*,*)' ISOTHERM (1) or ISOBAR (2) or MOL-FRACTION-SCAN (3)'
884  READ (*,*) iso_x
885 
886  other_t: DO
887 
888  IF (iso_x == 1) THEN
889  write (*,*) 'Specify end pressure [units as in input-file]'
890  READ (*,*) end_x
891  end_x = end_x*u_in_p
892  it(1) = 'x11' ! iteration of mol fraction of comp.1 phase 1
893  running = 'p' ! Press. is running variable in PHASE_EQUILIB
894  CALL select_sum_rel (1,0,1)
895  CALL select_sum_rel (2,0,2)
896  ELSEIF (iso_x == 2) THEN
897  write (*,*) 'Specify end temperature [units as in input-file]'
898  READ (*,*) end_x
899  end_x = end_x + u_in_t
900  it(1) = 'x11' ! iteration of mol fraction of comp.1 phase 1
901  running = 't' ! Temp. is running variable in PHASE_EQUILIB
902  CALL select_sum_rel (1,0,1)
903  CALL select_sum_rel (2,0,2)
904  ELSEIF (iso_x == 3) THEN
905  write (*,*) ' YOU HAVE DECIDED TO VARY THE MOLEFRACTION'
906  write (*,*) ' CHOOSE AN END-CONCENTRAION, FOR EXAMPLE:'
907  write (*,*) ' CALCULATE TOWARDS xi(1) = 0 (0)'
908  write (*,*) ' CALCULATE TOWARDS xi(1) = 1 (1)'
909  write (*,*) ' COMPOUND #1 IS: ',compna(1)
910  READ (*,*) end_x
911  IF (end_x <= 0.0) THEN
912  end_x = 1.E-12
913  ELSE
914  end_x = 1.0 - 1.E-12
915  ENDIF
916  running = 'x11' ! xi(1,1) is running var. in PHASE_EQUILIB
917  steps = 80.0
918  write(*,*)' ISOTHERM (1) or ISOBAR (2)'
919  READ (*,*) iso_x
920  IF (iso_x == 1) it(1) = 'p' ! iteration of pressure
921  IF (iso_x == 2) it(1) = 't' ! iteration of temperature
922  ENDIF
923 
924  outp = 1 ! output to terminal
925  CALL phase_equilib(end_x,steps,converg)
926 
927  IF ( .NOT.second .AND. iso_x == 1 .AND. running /= 'x11') THEN
928  second = .true.
929  p = pinit
930  val_conv = 0.0
931  CALL start_var(converg) ! gets starting values, sets "val_init"
932  end_x = 10.0 ! End-Druck f¸r zweiten Durchlauf
933  IF (eos >= 4) end_x = 1.e-3
934  steps = 8.0
935  CALL phase_equilib(end_x,steps,converg)
936  ENDIF
937 
938  IF (iso_x == 1) THEN
939  write (*,*) ' '
940  write (*,*) ' CALCULATE ANOTHER TEMPERATURE ?'
941  write (*,*) ' CHOOSE (0 FOR NO, 1 FOR YES)'
942  READ (*,*) iso_x
943  IF (iso_x == 1) THEN
944  WRITE (40,*) ' '
945  xif = 0.0
946  p = pinit
947  second =.false.
948  write (*,*) 'Specify temperature [units as in input-file]'
949  READ (*,*) t
950  t = t + u_in_t
951  val_conv = 0.0
952  CALL start_var(converg) ! gets starting values, sets "val_init"
953  steps = 20.0
954  it(2) = 'x21' ! iteration of mol fraction of comp.1 phase 2
955  sum_rel(1) = 'x12' ! summation relation: x12 = 1 - sum(x1j)
956  sum_rel(2) = 'x22' ! summation relation: x22 = 1 - sum(x2j)
957  cycle other_t
958  ENDIF
959  ENDIF
960  EXIT other_t
961  END DO other_t
962 
963  ENDIF
964 
965  write (*,*) ' '
966  write (*,*) 'kij(1,2)',kij(1,2)
967 
968 END SUBROUTINE binmix
969 
970 
971 
972 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
973 !
974 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
975 
976 SUBROUTINE bincrit
977 
978  USE basic_variables
979  IMPLICIT NONE
980 
981  !-----------------------------------------------------------------------------
982  INTEGER :: i, steps_i
983  !-----------------------------------------------------------------------------
984 
985  CALL read_input
986  nphas = 2
987 
988  steps_i = 20
989 
990  OPEN (71,file='./output_file/PT-crit.xlo')
991 
992  DO i = 1, steps_i-1
993  xif(2) = REAL( i ) / REAL( steps_i )
994  xif(1) = 1.0 - xif(2)
995  dense(1) = 0.15
996  CALL heidemann_khalil
997  write (71,'(8(2x,f18.8))') dense(1), t-u_out_t,p/u_out_p, &
998  xif(1:ncomp)
999  ENDDO
1000 
1001  write (*,*) ' '
1002  write (*,*) 'kij(1,2)',kij(1,2)
1003 
1004 END SUBROUTINE bincrit
1005 
1006 
1007 
1008 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1009 ! SUBROUTINE HE_MIX
1010 !
1011 ! This subroutine calculates physical properties of binary mixtures.
1012 ! The concentration ranges from one pure compound in 40 steps to
1013 ! the other pure substance.
1014 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1015 
1016 SUBROUTINE he_mix
1017 
1018  USE basic_variables
1019  use utilities, only: enthalpy_etc, si_dens
1020  IMPLICIT NONE
1021 
1022  !-----------------------------------------------------------------------------
1023  INTEGER :: i, i_steps
1024  REAL :: density(np),w(np,nc)
1025  REAL :: lnphi(np,nc)
1026  !-----------------------------------------------------------------------------
1027 
1028  CALL read_input
1029 
1030  OPEN (68,file='./output_file/caloric.xlo')
1031  WRITE (68,*) 'phi1 phi2 cp_res[J/molK] h_res[J/mol] g_res[J/mol] rho t[K] p[Pa] xi(1)'
1032 
1033  IF (ncomp /= 2) THEN
1034  WRITE (*,*)'CALCULATION FOR BINARY MIXTURE ONLY!'
1035  stop
1036  ENDIF
1037  nphas = 1
1038 
1039  densta(1) = 0.5
1040 
1041  i_steps = 40
1042 
1043  DO i= 1,(i_steps+1)
1044  xi(1,1) = 1.0 - REAL(i-1)/REAL(i_steps)
1045  xi(1,2) = 1.0 - xi(1,1)
1046 
1047  CALL fugacity (lnphi)
1048  CALL enthalpy_etc
1049 
1050  CALL si_dens (density,w)
1051  WRITE (*,*) cpres(1),enthal(1),dense(1),density(1),t,p,xi(1,1)
1052  WRITE (*,*) gibbs(1),xi(1,1)
1053  WRITE (68,'(9G15.7)') exp(lnphi(1,1)),exp(lnphi(1,2)), &
1054  cpres(1),enthal(1),gibbs(1),density(1),t,p,xi(1,1)
1055 
1056  END DO
1057 
1058 END SUBROUTINE he_mix
1059 
1060 
1061 
1062 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1063 ! SUBROUTINE STATE_PROP
1064 !
1065 ! This subroutine calculates the physical properties of a one-phase
1066 ! (liquid) multicomp. mixture at given concentration.
1067 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1068 
1069 SUBROUTINE state_prop
1070 
1071  USE basic_variables
1072  use utilities, only: enthalpy_etc, si_dens
1073  IMPLICIT NONE
1074 
1075  !-----------------------------------------------------------------------------
1076  INTEGER :: i
1077  REAL :: density(np), w(np,nc)
1078  REAL :: lnphi(np,nc), aver_m
1079  !-----------------------------------------------------------------------------
1080 
1081  CALL read_input
1082 
1083  nphas = 1
1084  densta(1) = 0.5
1085 
1086  WRITE (*,*) ' This subroutine calculates the physical properties'
1087  WRITE (*,*) ' of a one-phase (liquid) multicomponent mixture at'
1088  WRITE (*,*) ' given concentration'
1089  WRITE (*,*) ' '
1090  WRITE (*,'(a,i3,a)') ' SPECIFY MOL FRACTION OF ALL',ncomp,' COMPONENTS'
1091  DO i=1,ncomp
1092  READ (*,*) xi(1,i)
1093  END DO
1094 
1095 
1096  OPEN (68,file='./output_file/state_point.xlo')
1097 
1098  CALL fugacity (lnphi)
1099  CALL enthalpy_etc
1100 
1101  CALL si_dens (density,w)
1102 
1103  aver_m = sum( xi(1,1:ncomp)*mm(1:ncomp) )
1104 
1105  WRITE (*,*) cpres(1),enthal(1),dense(1),density(1),t-u_out_t, p/u_out_p
1106  WRITE (*,*) ' '
1107  WRITE (*,*) (xi(1,i), i = 1,ncomp)
1108  WRITE (68,*) 'cp_res[J/molK] h_res[J/mol] g_res[J/mol] rho t p '
1109  WRITE (68,*) cpres(1),enthal(1),gibbs(1),density(1),t-u_out_t, p/u_out_p
1110  WRITE (68,*) ' '
1111  WRITE (68,*) 'x(i) '
1112  WRITE (68,*) (compna(i),i=1,ncomp)
1113  WRITE (68,*) (xi(1,i),i=1,ncomp)
1114  WRITE (68,*) ' '
1115  WRITE (68,*) 'w(i) '
1116  WRITE (68,*) ((xi(1,i)*mm(i)/aver_m),i=1,ncomp)
1117  WRITE (68,*) ' '
1118  WRITE (68,*) 'average molar mass'
1119  WRITE (68,*) aver_m
1120  write (69,*) dense(1),dense(2),p
1121 
1122 END SUBROUTINE state_prop
1123 
1124 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1125 ! SUBROUTINE P_LINE
1126 !
1127 ! This subroutine calculates the physical properties of a one-phase
1128 ! (liquid) multicomp. mixture at given concentration.
1129 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1130 
1131 SUBROUTINE p_line
1132 
1133  USE basic_variables
1134  use utilities, only: p_calc
1135  IMPLICIT NONE
1136 
1137  !-----------------------------------------------------------------------------
1138  REAL :: pges,zges
1139  !-----------------------------------------------------------------------------
1140 
1141  CALL read_input
1142 
1143  nphas=1
1144  dense(1)= 0.0
1145  xi(1,1) = 1.0
1146  write (*,*) t,p
1147 
1148 
1149  DO WHILE (dense(1) < 0.25)
1150 
1151  dense(1) = dense(1) + 0.01
1152  CALL p_calc (pges,zges)
1153  write (71,*) dense(1),pges
1154 
1155  END DO
1156 
1157 END SUBROUTINE p_line
1158 
1159 
1160 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1161 ! SUBROUTINE P_T_DIAGRAM
1162 !
1163 ! This subroutine calculates P-T-diagrams for binary and ternary
1164 ! mixtures.
1165 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1166 
1167 SUBROUTINE p_t_diagram
1168 
1169  USE basic_variables
1170  USE starting_values, only: start_var
1171  use utilities, only: molefrac, switch
1172  IMPLICIT NONE
1173 
1174  !-----------------------------------------------------------------------------
1175  INTEGER :: i, iso_x, converg, swtch, ncompsave
1176  REAL :: steps, end_x, end_x1, end_x2, startcon
1177  REAL :: w(np,nc)
1178  !-----------------------------------------------------------------------------
1179 
1180  CALL read_input
1181  nphas = 2
1182  ncompsave = ncomp
1183  ncomp = 2
1184 
1185  CALL start_var(converg) ! gets starting values, sets "val_init"
1186  ncomp = ncompsave
1187 
1188  n_unkw = ncomp ! number of quantities to be iterated
1189 
1190  write (*,*) ' '
1191  write (*,*) ' CHOOSE A MASS FRACTION FOR COMP. #1 :',converg
1192  write (*,*) ' COMPOUND #1 IS ',compna(1)
1193  READ (*,*) w(1,1)
1194  IF (ncomp == 2) THEN
1195  w(1,2) =1.0 - w(1,1)
1196  CALL molefrac (w,1,0)
1197  end_x1 = LOG(xi(1,1))
1198 
1199  write (*,*) ' '
1200  write(*,*)' FOR APPROACHING THE END CONCENTRATION FOR'
1201  write(*,*)' COMPONENT # 1 CHOOSE CALCULATION PATH: '
1202  write(*,*)' ISOTHERM (1) or ISOBAR (2)'
1203  READ (*,*) iso_x
1204  write (*,*) ' '
1205  write(*,*)' END CONC. LEFT OF CRIT. POINT (yes:1 / no:0)'
1206  READ (*,*) swtch
1207  IF (swtch == 1) THEN
1208  CALL switch
1209  DO i=1,8
1210  val_conv(i)=0.0
1211  ENDDO
1212  ENDIF
1213 
1214  IF (iso_x == 1) it(1)='p' ! iteration of pressure
1215  IF (iso_x == 2) it(1)='t' ! iteration of temperature
1216  it(2)='x21' ! iteration of mol fraction of comp.1 phase 2
1217  sum_rel(1)='x12' ! summation relation: x12 = 1 - sum(x1j)
1218  sum_rel(2)='x22' ! summation relation: x22 = 1 - sum(x2j)
1219  outp = 1 ! output to terminal
1220  running = 'l11' ! lnx(1,1) is running var. in PHASE_EQUILIB
1221  steps=10.0
1222  CALL phase_equilib(end_x1,steps,converg)
1223 
1224  ELSEIF (ncomp == 3) THEN
1225  write (*,*) ' '
1226  write (*,*) ' CHOOSE A MASS FRACTION FOR COMP. #3 :'
1227  write (*,*) ' COMPOUND #3 IS ',compna(3)
1228  READ (*,*) w(1,3)
1229  ! xi(1,2) =1.0 - xi(1,1) - xi(1,3)
1230  w(1,2) =1.0 - w(1,1) - w(1,3)
1231  CALL molefrac (w,1,0)
1232  end_x1 = LOG(xi(1,1))
1233  end_x2 = LOG(xi(1,3))
1234 
1235  ! if convergence problems are encountered, play around here:
1236  startcon = -12.0 ! starting value for lnx(1,3)
1237  steps=10.0
1238 
1239  val_init(5) = val_conv(5)
1240  val_init(6) = val_conv(6)
1241  val_init(7) = startcon
1242  val_init(8) = val_conv(7)
1243  val_init(9) = val_conv(8)
1244  ! just an attempt for an initial guess for val_init(10):
1245  val_init(10) = startcon + (end_x2-startcon)/steps
1246  write (*,*) ' '
1247  write(*,*)' END CONC. LEFT OF CRIT. POINT (yes:1 / no:0)'
1248  READ (*,*) swtch
1249  IF (swtch == 1) CALL switch
1250 
1251  write (*,*) ' '
1252  write(*,*)' FOR APPROACHING THE END CONCENTRATION FOR'
1253  write(*,*)' COMPONENT # 3 CHOOSE CALCULATION PATH: '
1254  write(*,*)' ISOTHERM (1) or ISOBAR (2)'
1255  READ (*,*) iso_x
1256  IF (iso_x == 1) it(1)='p' ! iteration of pressure
1257  IF (iso_x == 2) it(1)='t' ! iteration of temperature
1258  it(2)='x21' ! iteration of mol fraction of comp.1 phase 2
1259  it(3)='x23' ! iteration of mol fraction of comp.3 phase 2
1260  sum_rel(1)='x12' ! summation relation: x12 = 1 - sum(x1j)
1261  sum_rel(2)='x22' ! summation relation: x22 = 1 - sum(x2j)
1262  outp = 1 ! output to terminal
1263  running='l13' ! lnx(1,3) is running var. in PHASE_EQUILIB
1264 
1265  CALL phase_equilib(end_x2,steps,converg)
1266 
1267 
1268  write (*,*) ' '
1269  write(*,*)' FOR APPROACHING THE END CONCENTRATION FOR'
1270  write(*,*)' COMPONENT # 1 CHOOSE CALCULATION PATH: '
1271  write(*,*)' ISOTHERM (1) or ISOBAR (2)'
1272  READ (*,*) iso_x
1273  IF (iso_x == 1) it(1)='p' ! iteration of pressure
1274  IF (iso_x == 2) it(1)='t' ! iteration of temperature
1275  it(2)='x21' ! iteration of mol fraction of comp.1 phase 2
1276  it(3)='x23' ! iteration of mol fraction of comp.3 phase 2
1277  sum_rel(1)='x12' ! summation relation: x12 = 1 - sum(x1j)
1278  sum_rel(2)='x22' ! summation relation: x22 = 1 - sum(x2j)
1279  outp = 1 ! output to terminal
1280  running='l11' ! lnx(1,1) is running var. in PHASE_EQUILIB
1281  steps=10.0
1282  CALL phase_equilib(end_x1,steps,converg)
1283 
1284  ENDIF
1285 
1286 
1287  ! start calculation of P-T-diagram
1288  write(*,*)' '
1289  write(*,*) ' SPECIFY TEMPERATURE (1)'
1290  write(*,*) ' OR PRESSURE (2)'
1291  READ (*,*) iso_x
1292  IF (iso_x == 1) THEN
1293  write(*,*) ' CHOOSE END TEMPERATURE [units as in input-file]:'
1294  READ (*,*) end_x
1295  end_x = end_x + u_in_t
1296  it(1)='p' ! iteration of pressure
1297  running='t' ! Temperature is running var. in PHASE_EQUILIB
1298  ELSEIF (iso_x == 2) THEN
1299  write(*,*) ' CHOOSE END PRESSURE [units as in input-file]:'
1300  READ (*,*) end_x
1301  end_x = end_x*u_in_p
1302  it(1)='t' ! iteration of temperature
1303  running='p' ! Temperature is running var. in PHASE_EQUILIB
1304  ENDIF
1305  outp = 1 ! output to terminal
1306  steps=20.0
1307  CALL phase_equilib(end_x,steps,converg)
1308 
1309 END SUBROUTINE p_t_diagram
1310 
1311 
1312 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1313 ! SUBROUTINE TERN_DIAGRAM
1314 !
1315 ! This subroutine calculates ternary diagrams for given (T and P)
1316 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1317 
1318 SUBROUTINE tern_diagram
1319 
1320  USE basic_variables
1321  USE starting_values, only: start_var
1322  use utilities, only: molefrac, switch
1323  IMPLICIT NONE
1324 
1325  !-----------------------------------------------------------------------------
1326  INTEGER :: converg, spec_ph, spec_comp
1327  REAL :: steps, end_x
1328  !-----------------------------------------------------------------------------
1329 
1330  CALL read_input
1331  nphas = 2
1332 
1333  CALL start_var(converg) ! gets starting values, sets "val_init"
1334 
1335  n_unkw = ncomp ! number of quantities to be iterated
1336 
1337  write (*,*) ' '
1338  write (*,*) ' THE MOLE FRACTION OF OF PHASE 1 OR OF PHASE 2'
1339  write (*,*) ' CAN NOW BE CHANGED. SELECT:'
1340  write (*,*) ' PHASE 1 (1)'
1341  write (*,*) ' PHASE 2 (2)'
1342  ! write (*,*) ' FEED-PHASE (3)'
1343  READ (*,*) spec_ph
1344  write (*,*) ' '
1345  write (*,*) ' SELECT A COMPONENT TO BE SPECIFIED'
1346  write (*,*) ' COMPONENT #1 ',compna(1),' (1)'
1347  write (*,*) ' COMPONENT #2 ',compna(2),' (2)'
1348  write (*,*) ' COMPONENT #3 ',compna(3),' (3)'
1349  READ (*,*) spec_comp
1350 
1351  write (*,*) ' '
1352  CALL output
1353  write (*,*) ' '
1354  write (*,*) ' CHOOSE AN END MOLE FRACTION FOR COMP. #',spec_comp
1355  write (*,*) ' (CAUTION: THREE-PHASE EQUIL. ARE NOT DETECTED)'
1356  READ (*,*) end_x
1357  end_x = LOG(end_x)
1358 
1359  steps = 10.0
1360 
1361  it(1) = 'x11' ! iteration of mol fraction of comp.1 phase 1
1362  it(2) = 'x21' ! iteration of mol fraction of comp.1 phase 2
1363  it(3) = 'x23' ! iteration of mol fraction of comp.3 phase 2
1364  sum_rel(1) = 'x12' ! summation relation: x12 = 1 - sum(x1j)
1365  sum_rel(2) = 'x22' ! summation relation: x22 = 1 - sum(x2j)
1366 
1367  IF (spec_ph == 1) THEN
1368  IF (spec_comp == 1) running='l11' ! lnx(1,1) is running var. in PHASE_EQUILIB
1369  IF (spec_comp == 1) it(1)='x13' ! iteration of mol fraction of comp.3 phase 1
1370  IF (spec_comp == 2) running='l12' ! lnx(1,2) is running var. in PHASE_EQUILIB
1371  IF (spec_comp == 2) sum_rel(1)='x13' ! summation relation: x12 = 1 - sum(x1j)
1372  IF (spec_comp == 3) running='l13' ! lnx(1,3) is running var. in PHASE_EQUILIB
1373  ELSEIF (spec_ph == 2) THEN
1374  IF (spec_comp == 1) running='l21' ! lnx(2,1) is running var. in PHASE_EQUILIB
1375  IF (spec_comp == 1) it(2)='x13' ! iteration of mol fraction of comp.3 phase 1
1376  IF (spec_comp == 2) running='l22' ! lnx(2,2) is running var. in PHASE_EQUILIB
1377  IF (spec_comp == 2) it(2)='x13' ! iteration of mol fraction of comp.3 phase 1
1378  IF (spec_comp == 2) sum_rel(2)='x21' ! summation relation: x21 = 1 - sum(x2j)
1379  IF (spec_comp == 3) running='l23' ! lnx(2,3) is running var. in PHASE_EQUILIB
1380  IF (spec_comp == 3) it(2)='x13' ! iteration of mol fraction of comp.3 phase 1
1381  IF (spec_comp == 3) it(3)='x21' ! iteration of mol fraction of comp.1 phase 2
1382  ELSEIF (spec_ph == 3) THEN
1383  ! it(1)='fls' ! do a flash-calcul. for fixed feed-conc
1384  ! IF (spec_comp == 1) running='xF1'
1385  ! IF (spec_comp == 2) running='xF2'
1386  ! IF (spec_comp == 3) running='xF3'
1387  ENDIF
1388  outp = 1 ! output to terminal
1389 
1390  CALL phase_equilib(end_x,steps,converg)
1391 
1392 END SUBROUTINE tern_diagram
1393 
1394 
1395 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1396 ! SUBROUTINE BIN_3PHASE
1397 !
1398 ! This subroutine performes VLLE calculations for binary mixtures
1399 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1400 
1401 SUBROUTINE bin_3phase
1402 
1403  USE basic_variables
1404  USE starting_values, only: start_var
1405  use utilities, only: si_dens
1406  IMPLICIT NONE
1407 
1408  !-----------------------------------------------------------------------------
1409  INTEGER :: k, ph, iso_x, converg, i
1410  REAL :: density(np), w(np,nc)
1411  !-----------------------------------------------------------------------------
1412 
1413  CALL read_input
1414  nphas = 2
1415 
1416  !--- determine liquid-liquid equilibria first
1417  CALL start_var(converg) ! gets starting values, sets "val_init"
1418 
1419  IF (converg == 1) THEN
1420  write (*,*) ' '
1421  write (*,*) ' LLE detected for starting values'
1422  write (*,*) ' '
1423  write (*,*) ' starting 3-phase calculation'
1424  write (*,*) ' '
1425  ENDIF
1426 
1427  write(*,*)' CHOOSE YOUR CALCULATION OPTION:'
1428  write(*,*)' ISOTHERM (1) or ISOBAR (2)'
1429  READ (*,*) iso_x
1430  nphas = 3
1431  n_unkw = 4 ! number of quantities to be iterated
1432  it(1)='x11' ! iteration of mol fraction of comp.1 phase 1
1433  it(2)='x21' ! iteration of mol fraction of comp.1 phase 2
1434  it(3)='x31' ! iteration of mol fraction of comp.1 phase 3
1435  IF (iso_x == 1) it(4)='p ' ! iteration of pressure
1436  IF (iso_x == 2) it(4)='t ' ! iteration of temperature
1437  sum_rel(1)='x12' ! summation relation: x12 = 1 - sum(x1j)
1438  sum_rel(2)='x22' ! summation relation: x22 = 1 - sum(x2j)
1439  sum_rel(3)='x32' ! summation relation: x32 = 1 - sum(x3j)
1440 
1441  ! die folgende Definition ist eine voruebergehende Kruecke:
1442  xi(3,1) = 0.01 ! set starting value for 3rd phase - temporary
1443  xi(3,2) = 1.0 - xi(3,1)
1444  lnx(3,1) = log(xi(3,1))
1445  lnx(3,2) = log(xi(3,2))
1446 
1447  val_init(0) = 1.e-8 ! set starting value for vapor density (3rd phase)
1448 
1449  DO ph = 1,nphas
1450  DO k = 1,ncomp
1451  val_init(4+k+(ph-1)*ncomp) = lnx(ph,k)
1452  ENDDO
1453  ENDDO
1454 
1455  other_condition: DO
1456 
1457  CALL objective_ctrl (converg)
1458 
1459  IF (converg == 1) THEN
1460  CALL si_dens (density,w)
1461  write(*,*) ' '
1462  write(*,*) ' 3-phase equilibrium calc. successful'
1463  write(*,*) ' '
1464  write(*,*) ' T, P'
1465  write(*,*) val_conv(3)-u_out_t,val_conv(4)/u_out_p
1466  write(*,*) ' x11 x21 x31'
1467  write(*,*) exp(val_conv(5)),exp(val_conv(7)),exp(val_conv(9))
1468  write(*,*) exp(val_conv(6)),exp(val_conv(8)),exp(val_conv(10))
1469  write(*,*) ' eta1 , eta2 , eta3'
1470  write(*,*) dense(1),dense(2),dense(3)
1471  write(*,*) ' '
1472  write(*,*) ' output written to file: ./output_file/3-phase.xlo'
1473 
1474  OPEN (67,file='./output_file/3-phase.xlo')
1475  write(67,*) ' T, P'
1476  write(67,*) val_conv(3)-u_out_t,val_conv(4)/u_out_p
1477  write(67,*) ' x11 x21 x31'
1478  write(67,*)exp(val_conv(5)),exp(val_conv(7)),exp(val_conv(9))
1479  write(67,*)exp(val_conv(6)),exp(val_conv(8)),exp(val_conv(10))
1480  write(67,*) ' rho1[kg/m3] rho2[kg/m3] rho3[kg/m3]'
1481  write(67,*) density(1),density(2),density(3)
1482  write(67,'(6e16.8)') val_conv(3), val_conv(4)/u_out_p, &
1483  1.0-exp(val_conv(5)),exp(val_conv(7)),exp(val_conv(9))
1484  ENDIF
1485 
1486  IF (t.lt.210.15) THEN
1487  DO i=1,(4+ncomp*nphas)
1488  val_init(i)=val_conv(i)
1489  ENDDO
1490  val_init(3) = val_init(3) + 0.025
1491  t = val_init(3)
1492  ELSE
1493  EXIT other_condition
1494  ENDIF
1495 
1496  END DO other_condition
1497 
1498 END SUBROUTINE bin_3phase
1499 
1500 
1501 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1502 ! SUBROUTINE TERN_3PHASE
1503 !
1504 ! This subroutine performes VLLE calculations for binary mixtures
1505 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1506 
1507 SUBROUTINE tern_3phase
1508 
1509  USE basic_variables
1510  use utilities, only: si_dens
1511  IMPLICIT NONE
1512 
1513  !-----------------------------------------------------------------------------
1514  INTEGER :: k, ph, converg
1515  REAL :: density(np), w(np,nc)
1516  !-----------------------------------------------------------------------------
1517 
1518  CALL read_input
1519  nphas = 3
1520 
1521  lnx(1,1) = -4.375218662
1522  lnx(1,2) = -0.241363748
1523  lnx(1,3) = -1.600186935
1524  lnx(2,1) = -14.98669354
1525  lnx(2,2) = -0.395440676
1526  lnx(2,3) = -1.118968702
1527  lnx(3,1) = -30.0
1528  lnx(3,2) = -0.821288895
1529  lnx(3,3) = -0.579576292
1530 
1531  t = 449.0
1532  p = 16.e6
1533  val_init(3) = t
1534  val_init(4) = p
1535 
1536  val_init(1) = 0.4 ! set starting value for vapor density (1st phase)
1537  val_init(2) = 0.4 ! set starting value for vapor density (2nd phase)
1538  val_init(0) = 1.e-5 ! set starting value for vapor density (3rd phase)
1539 
1540  n_unkw = 6 ! number of quantities to be iterated
1541  it(1)='x21' ! iteration of mol fraction of comp.1 phase 2
1542  it(2)='x22' ! iteration of mol fraction of comp.2 phase 2
1543  it(3)='x31' ! iteration of mol fraction of comp.1 phase 3
1544  it(4)='x32' ! iteration of mol fraction of comp.2 phase 3
1545  it(5)='t ' ! temp
1546  it(6)='p ' ! pressure
1547  sum_rel(1)='x13' ! summation relation: x12 = 1 - sum(x1j)
1548  sum_rel(2)='x23' ! summation relation: x22 = 1 - sum(x2j)
1549  sum_rel(3)='x33' ! summation relation: x32 = 1 - sum(x3j)
1550 
1551 
1552  DO ph = 1,nphas
1553  DO k = 1,ncomp
1554  val_init(4+k+(ph-1)*ncomp) = lnx(ph,k)
1555  ENDDO
1556  ENDDO
1557 
1558  ! 10 continue
1559 
1560  CALL objective_ctrl (converg)
1561 
1562  IF (converg == 1) THEN
1563  CALL si_dens (density,w)
1564  write(*,*) ' '
1565  write(*,*) ' 3-phase equilibrium calc. successful'
1566  write(*,*) ' '
1567  write(*,*) ' T, P'
1568  write(*,*) val_conv(3)-u_out_t,val_conv(4)/u_out_p
1569  write(*,*) ' x_phase1 x_phase2 x_phase3'
1570  write(*,*)exp(val_conv(5)),exp(val_conv(8)),exp(val_conv(11))
1571  write(*,*)exp(val_conv(6)),exp(val_conv(9)),exp(val_conv(12))
1572  write(*,*)exp(val_conv(7)),exp(val_conv(10)),exp(val_conv(13))
1573  write(*,*) ' w_phase1 w_phase2 w_phase3'
1574  write(*,*)w(1,1),w(2,1),w(3,1)
1575  write(*,*)w(1,2),w(2,2),w(3,2)
1576  write(*,*)w(1,3),w(2,3),w(3,3)
1577  write(*,*) ' eta1 , eta2 , eta3'
1578  write(*,*) dense(1),dense(2),dense(3)
1579  write(*,*) ' '
1580  write(*,*) ' output written to file: ./output_file/3-phase.xlo'
1581 
1582  OPEN (67,file='./output_file/3-phase.xlo')
1583  write(67,*) ' T, P'
1584  write(67,*) val_conv(3)-u_out_t,val_conv(4)/u_out_p
1585  write(67,*) ' x_phase1 x_phase2 x_phase3'
1586  write(67,*)exp(val_conv(5)),exp(val_conv(8)),exp(val_conv(11))
1587  write(67,*)exp(val_conv(6)),exp(val_conv(9)),exp(val_conv(12))
1588  write(67,*)exp(val_conv(7)),exp(val_conv(10)),exp(val_conv(13))
1589  write(67,*) ' rho1[kg/m3] rho2[kg/m3] rho3[kg/m3]'
1590  write(67,*) density(1),density(2),density(3)
1591  write(67,'(6e16.8)') val_conv(3), val_conv(4)/u_out_p, &
1592  1.0-exp(val_conv(5)),exp(val_conv(7)),exp(val_conv(9))
1593  ENDIF
1594 
1595  write (*,*) ' '
1596  write (*,*) parame(1,13),parame(1,16)
1597  write (*,*) parame(1,3),kij(1,2),kij(1,3)
1598  write (*,*) 'modify parameter eps(1)/k'
1599  read(*,*) parame(1,3)
1600  write (*,*) 'modify parameter eps_assoc'
1601  read(*,*) parame(1,13)
1602  ! parame(1,16) = parame(1,15)
1603  write (*,*) 'modify parameter kij(1,2)'
1604  read(*,*) kij(1,2)
1605  kij(2,1) = kij(1,2)
1606  write (*,*) 'modify parameter kij(1,3)'
1607  read(*,*) kij(1,3)
1608  kij(3,1) = kij(1,3)
1609 
1610  ! IF (t.lt.473.15) THEN
1611  ! DO i=1,(4+ncomp*nphas)
1612  ! val_init(i)=val_conv(i)
1613  ! ENDDO
1614  ! val_init(3) = val_init(3) + 5.0
1615  ! t = val_init(3)
1616  ! GOTO 10
1617  ! ENDIF
1618 
1619 END SUBROUTINE tern_3phase
1620 
1621 
1622 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1623 ! SUBROUTINE PARITER_CRIT
1624 !
1625 ! This subroutine calculates the three pure component parameters of
1626 ! PC-SAFT for non-polar and non-associating components from the
1627 ! critical temp., the critical pressur and the acentric factor
1628 ! OMEGA.
1629 ! This procedure sacrifices the density in the liquid phase and
1630 ! sacrifices therefore also caloric properties, like cp and enthalpy
1631 ! (of vaporization).
1632 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1633 
1634 SUBROUTINE pariter_crit
1635 
1636  USE basic_variables
1637  USE utilities, only: critical, paus
1638  IMPLICIT NONE
1639 
1640  !-----------------------------------------------------------------------------
1641  INTEGER :: i, j, converg
1642  REAL :: sig_new, eps_new, fvek1, fvek2
1643  REAL :: sig_save, eps_save
1644  REAL :: df1dx, df2dx, df1dy, df2dy
1645  REAL :: betrag, attenu
1646  REAL :: tc, pc, rhoc, omega, seg_new, par1_sv, par2_sv
1647  REAL :: par3_sv, omeg_sv, omeg_nw, err_sv, err_nw, derrdr
1648  REAL :: steps, end_x
1649  !-----------------------------------------------------------------------------
1650 
1651 
1652  OPEN (10,file = './output_file/crit-par.out')
1653 
1654  write (*,*) ' '
1655  write (*,*) ' ENTER MOLECULAR MASS M / g/mol :'
1656  READ (*,*) mm(1)
1657  ! mm(1)=16.043
1658  write (*,*) ' ENTER CRITICAL TEMPERATURE Tc / K :'
1659  READ (*,*) tc
1660  ! tc = - 82.59
1661  write (*,*) ' ENTER CRITICAL PRESSURE Pc / bar :'
1662  READ (*,*) pc
1663  ! pc = 45.95
1664  write (*,*) ' ENTER ACENTRIC FACTOR omega :'
1665  READ (*,*) omega
1666  ! omega = 0.01
1667  ! tc = tc + 273.15
1668  pc = pc * 1.e5
1669 
1670  eos = 1
1671  ncomp = 1
1672  nphas = 2
1673 
1674  !----------------------------------------------------------
1675  sig_new = 4.2
1676  eps_new = 240.0
1677  seg_new = mm(1)*0.05
1678  !----------------------------------------------------------
1679  parame(1,4) = 0.0
1680  parame(1,5) = 0.12
1681  parame(1,6) = 0.0
1682  parame(1,7) = 0.0
1683  DO j=8,25
1684  parame(1,j) = 0.0
1685  ENDDO
1686  scaling(1) = 1.e6
1687 
1688  attenu = 1.0
1689 
1690  crit_loop: DO
1691 
1692  DO i = 1,2
1693 
1694  j = 1
1695  DO WHILE ( (abs(fvek1)+abs(fvek2)) > 2.6e-4 )
1696  j = j + 1
1697 
1698  parame(1,2) = sig_new
1699  parame(1,3) = eps_new
1700 
1701  !------------------------------
1702  xi(1,1) = 1.0
1703  parame(1,1) = seg_new + 0.01*REAL(i-1)
1704  t = tc
1705  p = pc
1706 
1707  CALL critical (tc,pc,rhoc)
1708 
1709  fvek1= (t-tc)/tc
1710  fvek2= (p-pc)/pc
1711 
1712  ! Bildung der numerischen Ableitung nach sigma:
1713  sig_save = parame(1,2)
1714  parame(1,2) = sig_save*1.0002
1715 
1716  CALL critical (tc,pc,rhoc)
1717 
1718  df1dx = (fvek1-(t-tc)/tc )/( sig_save-parame(1,2) )
1719  df2dx = (fvek2-(p-pc)/pc )/( sig_save-parame(1,2) )
1720  parame(1,2) = sig_save
1721 
1722 
1723 
1724  ! Bildung der numerischen Ableitung nach epsilon:
1725  eps_save = parame(1,3)
1726  parame(1,3) = eps_save*1.0002
1727 
1728  CALL critical (tc,pc,rhoc)
1729 
1730  df1dy = (fvek1-(t-tc)/tc )/( eps_save-parame(1,3) )
1731  df2dy = (fvek2-(p-pc)/pc )/( eps_save-parame(1,3) )
1732  parame(1,3) = eps_save
1733 
1734  betrag = df1dx*df2dy - df1dy*df2dx
1735  sig_new = parame(1,2) - attenu / betrag * (df2dy*fvek1+(-df1dy)*fvek2)
1736  eps_new = parame(1,3) - attenu / betrag * ((-df2dx)*fvek1+df1dx*fvek2)
1737 
1738  write (10,*) parame(1,2),parame(1,3),fvek1,fvek2,i
1739  write (*,*) parame(1,2),parame(1,3),fvek1,fvek2,i
1740  ! if (j == 15) stop
1741 
1742  ! IF (((ABS(fvek1) > 1.E-4).OR.(ABS(fvek2) > 1.E-4)).
1743  END DO
1744 
1745 
1746  ! --- calculate vapor pressure at 0.7*tc (for acentric factor)
1747  n_unkw = ncomp ! number of quantities to be iterated
1748  it(1) = 'p' ! iteration of mol fraction of comp.1 phase 1
1749  running = 't' ! Temp. is running variable in PHASE_EQUILIB
1750 
1751  val_init(1) = 0.5
1752  val_init(2) = 1.e-6
1753  val_init(3) = 0.7*tc
1754  val_init(4) = 1.e+2 ! default starting value for P: 100 Pa
1755  val_init(5) = 0.0 ! logar. mole fraction: lnx=0, x=1, phase 1
1756  val_init(6) = 0.0 ! logar. mole fraction: lnx=0, x=1, phase 2
1757 
1758  end_x = 0.7*tc
1759 
1760  steps=1.0
1761  outp = 0 ! output to terminal
1762  CALL phase_equilib(end_x,steps,converg)
1763  IF (converg /= 1) THEN
1764  call paus ('PARITER_CRIT: equilibrium calculation failed')
1765  ENDIF
1766 
1767  IF (i == 1) THEN
1768  par1_sv = parame(1,2)
1769  par2_sv = parame(1,3)
1770  par3_sv = parame(1,1)
1771  omeg_sv = -1.0 - log10(val_conv(4)/pc)
1772  err_sv = omeg_sv-omega
1773  write (10,*) 'error',i,err_sv
1774  write (*,*) 'error',i,err_sv
1775  ELSE
1776  omeg_nw = -1.0 - log10(val_conv(4)/pc)
1777  err_nw = omeg_nw-omega
1778  write (10,*) 'error',i,err_nw
1779  write (*,*) 'error',i,err_nw
1780  ENDIF
1781  write (10,*) 'omega',i,-1.0 - log10(val_conv(4)/pc)
1782  write (*,*) 'omega',i,-1.0 - log10(val_conv(4)/pc)
1783 
1784  END DO
1785 
1786  derrdr =(err_sv - err_nw)/(par3_sv-parame(1,1))
1787  seg_new = par3_sv - 0.9 * err_sv/derrdr
1788  write (10,*) 'derivative',derrdr
1789  write (10,*) 'new value for m',seg_new
1790  write (*,*) 'derivative',derrdr
1791  write (*,*) 'new value for m',seg_new
1792 
1793  IF (abs(err_sv) <= 3.e-5) EXIT crit_loop
1794 
1795  END DO crit_loop
1796 
1797  write (*,*) ' '
1798  write (*,*) ' Pure component parameters :'
1799  write (*,*) ' '
1800  write (*,*) ' sigma = ',par1_sv
1801  write (*,*) ' eps/k = ',par2_sv
1802  write (*,*) ' m = ',par3_sv
1803  write (10,*) ' '
1804  write (10,*) ' Pure component parameters :'
1805  write (10,*) ' '
1806  write (10,*) ' mm(i) = ',mm(1)
1807  write (10,*) ' parame(i,1) = ',par3_sv
1808  write (10,*) ' parame(i,2) = ',par1_sv
1809  write (10,*) ' parame(i,3) = ',par2_sv
1810 
1811  CLOSE (10)
1812 
1813 END SUBROUTINE pariter_crit
1814 
1815 
parameters
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains constant...
Definition: modules.f90:6
eos_constants
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains paramete...
Definition: modules.f90:200
eos_variables
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains paramete...
Definition: modules.f90:120
starting_values
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW Module STARTING_VALUES This m...
Definition: modules.f90:251
basic_variables
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains paramete...
Definition: modules.f90:29
eos_numerical_derivatives
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains paramete...
Definition: modules.f90:220
starting_values::start_var
subroutine, public start_var(converg)
subroutine start_var
Definition: starting_value.f90:65
pc_saft
program pc_saft
Starting Point of PC-SAFT program.
Definition: main_prog.f90:32