Demo entry 6760695

fortran

   

Submitted by anonymous on Sep 18, 2018 at 04:11
Language: Fortran. Code size: 4.5 kB.

subroutine uvn_smax(u_spec,v_spec,ni_spec,tg,pg,nij_in,&
            switch,t,p,np,lamda,zj,mfrac)
use uvn_prop_mod, only: nc_uvn,iout_UVn_main,stab_1p_sol
use uvn_prop_mod, only: t_1p_sol

implicit none

!+++++++
! Purpose:
!   solve nUV flash with entropy maximization
!
! Epilog:
!   u_spec: specified internal energy [j]
!   v_spec: specified volume [m^3]
!   ni_spec: feeding composition of each species [mole]
!
!   np: number of phases
!
! Notes:
!   (1) start another phase only when instability occurs
!
!-------
!integer, intent(in) :: ind
real*8, intent(in) :: u_spec,v_spec,tg,pg
real*8, dimension(nc_uvn), intent(in) :: ni_spec
real*8, dimension(nc_uvn,2), intent(in) :: nij_in

logical, intent(out) :: switch
real*8, intent(out) :: t,p
integer, intent(out) :: np
real*8, dimension(nc_uvn,3), intent(out) :: mfrac
real*8, dimension(3), intent(out) :: lamda,zj

integer :: i,j,nc
integer, dimension(nc_uvn+2) :: dep_var
real*8 :: dum1,dum2
real*8, dimension(nc_uvn) :: z_test
real*8, dimension(nc_uvn) :: k_2p,k_triv,k_2p_stab,k_3p_stab
real*8, dimension(nc_uvn,3) :: nij,xij
real*8, dimension(3) :: uj,vj,mwmj
logical :: stable_2p,stable_3p,merge,conv

switch = .false.
t = 0.0d0
p = 0.0d0
np = 0
nc = nc_uvn

!initializations
mfrac = 0.0d0
lamda = 0.0d0
zj = 0.0d0
uj = 0.0d0
vj = 0.0d0
nij = 0.0d0
xij = 0.0d0
mwmj = 0.0d0

!initialize reference state values
!needed in calculations with chemical potential and its derivative
!call nUV_ref_reid

!**********step 1**********
!1. s maximization for 2 phase
np = 2

100 continue
!1.1 initial estimate of independent variables, uj,vj and nij
call uvn_init_2p(u_spec,v_spec,ni_spec,tg,pg,&
         uj(1:np),vj(1:np),nij(:,1:np))

!1.2 solve with s. maximization
call uvn_smax_sub(np,u_spec,v_spec,ni_spec,tg,pg,conv,&
       switch,uj(1:np),vj(1:np),zj(1:np),mwmj(1:np),&
	   nij(:,1:np),xij(:,1:np),t,p)


!First check if switch to nested TP method
if(switch)  then
!  write(iout_UVn_main,10) 
  return 
end if

!outputs
!write(iout_UVn_main,15) 
!write(*,15) 
!stop

!**********step 1**********

if(conv) then
!**********step 2**********
!2. check if 2 phase solution is stable
!   yes -> end of calculation
!    no -> continue to 3p calculation
!
! Only 1 phase is needed for stability test.

!merge 2 phases if possible
call uvn_merge_2p(uj,vj,nij,xij,merge)

!+++
!Note
!  (1) if merge occurs, trivial solution is found
!  (2) On the other hand, the stability already shows instability
!      further iterations are needed as the true solution is found.
!---
if(merge) then
  !if merge occurs but the stability shows unstable, need to start another iter.
  if(.not.stab_1p_sol) then
    conv = .false.
    goto 100
  end if
  write(iout_UVn_main,16) 
  np = np - 1
  call uvn_output4(np,t,p,uj(1:np),vj(1:np),zj(1:np),mwmj(1:np),&
                   nij(:,1:np),xij(:,1:np))
  return
else
  goto 21
  !goto 20
end if
!**********step 2**********
else
  goto 100
end if  !start with another iteration

!stop
21 continue

mfrac = xij

do i=1,nc
  k_2p(i) = mfrac(i,1) / mfrac(i,2) 
end do

!check stability of the reference phase (the 2nd phase)
z_test = mfrac(:,2)
k_triv = k_2p

call stability(T,P,z_test,k_triv,stable_2p,dum1,k_2p_stab,dum2)

!write(iout_UVn_main,25) 2,stable_2p

!stable 2 phase equilibrium state
if(stable_2p) then
  write(iout_UVn_main,26)   
  goto 20
end if

write(iout_UVn_main,30) 
!++++++++
!Normal output
20 continue

!convert to phase mole fraction
call uvn_convo(np,nij,lamda(1:np),mfrac(:,1:np))

!end of calculation, output results
!call uvn_output2(np,t,p,zj(1:np),lamda(1:np),mfrac(:,1:np))
!--------


10 format(/2x,'*******program switches to nested TP method')
15 format(/2x,'2 phase s max. successfully done...',&
       /2x,'*******test possible merge')
16 format(/2x,'merge is successfully done...',&
       /2x,'*******merge is found')
17 format(/2x,'merge is successfully done...',&
       /2x,'*******merge is not found',&
       /2x,'*******program goes to its stability test')
25 format(/2x,i4,2x,'phase stability test started...',&
      /6x,'global phase stability result:',l3)
26 format(/2x,'*******stable 2 phase found...')
30 format(/2x,'*******unstable 2 phase !!!...',&
      /2x,'*******program continues to 3 phase calculation...')
45 format(/2x,'instability is found at 3 phase mixture')
!!!!!!!!!
return
end subroutine uvn_smax

This snippet took 0.01 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).