cylinder_model.f95

!#######################################################################
!Code to apply the cylinder model for computing the radiative absorption
!and net radiative loss.

!The code expects an input file in the .kg format (post-processed 
!flamelet computed by the FlameMaster solver), and outputs a file with 
!the radiative emission, radiative absorption and net radiative loss.

!Both gray and non-gray (using the WSGG model with correlations by
!Bordbar et al., 2014) calculations are supported.

!Developed by Guilherme Fraga, Bifen Wu and Xinyu Zhao at UCONN.
!Contact: guilherme.fraga@uconn.edu
!			 bifen.wu@uconn.edu
!			 xinyu.zhao@uconn.edu
!#######################################################################
program cylinder_model

implicit none
character(200) :: input_file,output_file
integer,parameter :: eb=kind(0.d0)
integer :: i,j,l,iinv,itheta,ipsi,ngg
integer :: ncells,ntheta,npsi,ntau
integer :: input_unit,output_unit
real(eb) :: asum,dr,tau0,Gint,dtheta,dpsi
real(eb) :: dummy_eb,cp,lambda,rho,Mmix
real(eb) :: MCO2,MCO,MH2O,YCO2,YCO,YH2O
real(eb) :: g1,g2
real(eb) :: lbl_t(126),lbl_h2o(126),lbl_co2(126),lbl_co(126)
real(eb),parameter :: sigma = 5.669e-8_eb
real(eb),parameter :: pi = 3.14159_eb
real(eb),parameter :: small = 1.e-10_eb
real(eb),allocatable,dimension(:) :: a,Sb,tau,theta,psi,Rad_abs,Rad_emi
real(eb),allocatable,dimension(:) :: r,Z,XCO2,XCO,XH2O,T,alpha,chi
logical :: wsgg_model = .false.

!-----------------
!Input parameters
!-----------------
input_file = 'N-C7H16_p01_0chi0.00902576tf0371to0300Tst1046.kg'
output_file = 'output.csv'
wsgg_model = .false.

if (wsgg_model) then
	ngg = 4
else
	ngg = 1
endif

!------------------------------------
! Discretization parameters (part 1)
!------------------------------------
ncells = 402
ntheta = 50
npsi = 100
ntau = 100

!-------------------
! Allocating arrays
!-------------------
allocate(a(1:ncells))
allocate(r(1:ncells))
allocate(Z(1:ncells))
allocate(XCO(1:ncells))
allocate(XCO2(1:ncells))
allocate(XH2O(1:ncells))
allocate(alpha(1:ncells))
allocate(chi(1:ncells))
allocate(T(1:ncells))
allocate(Sb(1:ncells))
allocate(tau(1:ncells))
allocate(theta(1:ntheta))
allocate(psi(1:npsi))
allocate(Rad_abs(1:ncells))
allocate(Rad_emi(1:ncells))

!------------------------------------
! Discretization parameters (part 2)
!------------------------------------
dtheta = pi/real(ntheta-1,eb)
theta(1) = 0._eb
do i=2,ntheta
	theta(i) = theta(i-1) + dtheta
enddo

psi(1) = -0.5_eb*pi
dpsi = 2*pi/real(npsi-1,eb)
do i=2,npsi
	psi(i) = psi(i-1) + dpsi
enddo

!------------------
! Reading the data
!------------------
!Getting and opening the unit
input_unit = 20
open(file=input_file,unit=input_unit)

!Skipping the first two lines of data
read(input_unit,*)
read(input_unit,*)

!Reading
do i=1,ncells
	iinv = 1+ncells-i
	read(input_unit,*) Z(iinv),T(iinv),(dummy_eb,l=1,7),&						!This may need to be modified depending on the fuel
		YH2O,YCO2,(dummy_eb,l=1,2),YCO,(dummy_eb,l=1,22),&
      Mmix,(dummy_eb,l=1,2),chi(iinv),rho,lambda,cp

	!Computing the local mole fractions
	MH2O = 18._eb
   MCO2 = 44._eb
   MCO = 28._eb
	XH2O(iinv) = YH2O*Mmix / (MH2O + small)
   XCO2(iinv) = YCO2*Mmix / (MCO2 + small)
   XCO(iinv) = YCO*Mmix / (MCO + small)
    
   !Computing the local thermal diffusivity
   alpha(iinv) = lambda/(rho*cp + small)

enddo

close(input_unit)

!--------------
! Precomputing 
!--------------
if (.not.wsgg_model) call load_kp_data()			!Load data for the Planck-mean absorption coefficient

!Converting from Z to r
r(1) = 0.0_eb
r(2) = sqrt((alpha(2)+alpha(1))/(0.5_eb*(chi(2)+chi(1))+small))*dabs(Z(2)-Z(1))
do i=3,ncells-1
	r(i) = r(i-1) + (0.5*sqrt(2._eb*alpha(i)/(chi(i)+small)) + 0.5*sqrt(2._eb*alpha(i-1)/(chi(i-1)+small)))*&
		abs(Z(i)-Z(i-1))
enddo
r(ncells) = r(ncells-1) + sqrt((alpha(ncells)+alpha(ncells-1))/(0.5_eb*(chi(ncells)+chi(ncells-1))+small))*&
	dabs(Z(ncells)-Z(ncells-1))

Rad_abs = 0._eb
Rad_emi = 0._eb
do j=1,ngg

	!Defining radiative properties
	do i=1,ncells
		if (wsgg_model) then 
			!Source function
			Sb(i) = get_aj(T(i),XH2O(i),XCO2(i),j)*sigma*(T(i)**4._eb)/pi
			
			!Planck-mean absorption coefficient
			a(i) = get_kj(XH2O(i),XCO2(i),j)
		
		else
			!Source function
			Sb(i) = sigma*(T(i)**4._eb)/pi
			
			!Planck-mean absorption coefficient
			a(i) = get_kp(T(i),XH2O(i),XCO2(i),XCO(i))
		endif
	enddo

	!Computing predefined optical thickness
	asum = 0._eb
	do i=2,ncells
		dr = r(i)-r(i-1)
		tau(i) = tau(i-1) + 0.5_eb*(a(i) + a(i-1))*dr
	enddo
	tau0 = tau(ncells)                        !Optical thickness tau0

	!Computing all integrals needed for G
	do i=1,ncells
		Gint = 0._eb
		do itheta=1,ntheta
				do ipsi=1,npsi
					if (psi(ipsi).le.(pi/2._eb)) then
						g1 = g_integral(i,psi(ipsi),theta(itheta),&
								tau(i)*abs(sin(psi(ipsi))),tau0,.true.)
						g2 = g_integral(i,psi(ipsi),theta(itheta),&
								tau(i)*abs(sin(psi(ipsi))),tau(i),.false.)
					else
						g1 = g_integral(i,psi(ipsi),theta(itheta),&
								tau(i),tau0,.true.)
						g2 = 0._eb
					endif
					Gint = Gint + (g1+g2)*sin(theta(itheta))*dtheta*dpsi
			enddo
		enddo
		Rad_abs(i) = Rad_abs(i) + a(i)*Gint
		Rad_emi(i) = Rad_emi(i) + 4._eb*pi*a(i)*Sb(i)
	enddo
enddo

output_unit = 21
open(file=output_file,unit=output_unit)
write(output_unit,'(100(a,:,","))') 'Z','T','XH2O','XCO2','XCO',&
	'CHI','ALPHA','SB','KP','R','TAU','ABS','EMI','NET'
do i=1,ncells
	write(output_unit,'(100(e22.15,:,","))') Z(i),T(i),XH2O(i),&
    	XCO2(i),XCO(i),chi(i),alpha(i),Sb(i),a(i),r(i),tau(i),Rad_abs(i),Rad_emi(i),Rad_emi(i)-Rad_abs(i)
enddo

contains

subroutine load_kp_data()

	character(200) :: h2o_file,co2_file,co_file
	integer :: h2o_unit,co2_unit,co_unit
	integer :: ii,int_t
	
	!File names
	h2o_file = 'kpl_h2o.dat'
	co2_file = 'kpl_co2.dat'
	co_file = 'kpl_co.dat'
	
	!File units
	h2o_unit = 50
	co2_unit = h2o_unit+1
	co_unit = co2_unit+1

	!Open units
	open(file=h2o_file,unit=h2o_unit)
	open(file=co2_file,unit=co2_unit)
	open(file=co_file,unit=co_unit)
	
	!Reading files
	do ii=1,126
		read(h2o_unit,*) lbl_t(ii),lbl_h2o(ii)
		read(co2_unit,*) lbl_t(ii),lbl_co2(ii)
		read(co_unit,*) lbl_t(ii),lbl_co(ii)
	enddo
	
	!Close units
	close(h2o_unit)
	close(co2_unit)
	close(co_unit)

endsubroutine load_kp_data

real(eb) function get_kp(tmp,xxh2o,xxco2,xxco)

	real(eb),intent(in) :: tmp,xxh2o,xxco2,xxco
	real(eb) :: kp_h2o,kp_co2,kp_co
	
	kp_h2o = LinInterp(tmp,lbl_t,lbl_h2o)
	kp_co2 = LinInterp(tmp,lbl_t,lbl_co2)
	kp_co = LinInterp(tmp,lbl_t,lbl_co)
	
	get_kp = (kp_h2o*xxh2o + kp_co2*xxco2 + kp_co*xxco)*100._eb

endfunction get_kp

real(eb) function g_integral(i_in,psi_in,theta_in,mintau,maxtau,plus)

	integer,intent(in) :: i_in
   integer :: ii,nintpoints
   logical,intent(in) :: plus
   real(eb),intent(in) :: psi_in,theta_in,mintau,maxtau
   real(eb) :: gsum,dtau,taufrac,expterm,sqrtterm,taull,Sbll
   real(eb),allocatable,dimension(:) :: Sbl,taul
	logical :: isnumber,notinfty

	!Define integration parameters and arrays
   nintpoints = ntau
   allocate(Sbl(1:nintpoints))
   allocate(taul(1:nintpoints))

   !Define integration bounds
   dtau = (maxtau - mintau)/real(nintpoints-1,eb)

   !Define interpolated arrays
   taul(1) = mintau
   Sbl(1) = LinInterp(taul(1),tau,Sb)
	do ii=2,nintpoints
		taul(ii) = taul(ii-1) + dtau
		Sbl(ii) = LinInterp(taul(ii),tau,Sb)
   enddo
 
	!Compute the integral    
    gsum = 0._eb
    do ii=1,nintpoints-1
		taull = (taul(ii) + taul(ii+1))/2._eb
		Sbll = ( Sbl(ii) + Sbl(ii+1) )/2._eb
    
		sqrtterm = sqrt( taull**2._eb - (tau(i_in)*sin(psi_in))**2._eb )
		taufrac = taull / (sin(theta_in)* sqrtterm )
        	
		if (plus) then
			expterm = - (tau(i_in)*cos(psi_in) + sqrtterm )/sin(theta_in)
		else
			expterm = - (tau(i_in)*cos(psi_in) - sqrtterm )/sin(theta_in)
		endif
       	
      if (isnan(Sbll*taufrac*exp(expterm)*dtau)) then
			isnumber = .false.
		else
			isnumber = .true.
		endif
      
      if ((Sbll*taufrac*exp(expterm)*dtau).gt.huge(eb)) then
			notinfty = .false.
		else
			notinfty = .true.
      endif
      
      if (isnumber.and.notinfty) gsum = gsum + Sbll*taufrac*exp(expterm)*dtau
          
    enddo
    g_integral = gsum

    !Deallocate arrays
    deallocate(taul)
    deallocate(Sbl)

!stop
endfunction g_integral    

real(eb) function LinInterp(xtarget,xarray,yarray)

	integer :: ii,ntotal
	real(eb),intent(in) :: xtarget,xarray(:),yarray(:)
   real(eb) :: x1,y1,rxy
    
   ntotal = size(xarray)

   !Find bounds
   if (xtarget.le.xarray(1)) then
		x1 = 0._eb
      y1 = 0._eb
      rxy = yarray(1)/(xarray(1)+small)
    elseif(xtarget.ge.xarray(ntotal)) then
    	x1 = xarray(ntotal)
      y1 = yarray(ntotal)
      rxy = (yarray(ntotal) - yarray(ntotal-1))/(xarray(ntotal) - xarray(ntotal-1) + small)
	else
    	do ii=2,ntotal
			if ((xtarget.gt.xarray(ii-1)).and.(xtarget.le.xarray(ii))) then
				x1 = xarray(ii-1)
            y1 = yarray(ii-1)
            rxy = (yarray(ii) - yarray(ii-1)) / (xarray(ii) - xarray(ii-1) + small)
           endif
       enddo
    endif

    !Interpolate
    LinInterp = y1 + rxy*(xtarget - x1)

endfunction LinInterp

real(eb) function get_kj(xxw,xxc,jj)

	integer,intent(in) :: jj
	integer :: ii
	real(eb),intent(in) :: xxw, xxc
	real(eb) :: mol_rat,d(1:4,0:4),sum_kappa

	mol_rat = xxw/(xxc + small)

	!Trim mole ratio to the bounds of the correlations (0.01 - 4)
	if (mol_rat<0.01_eb) mol_rat = 0.01_eb
	if (mol_rat>4._eb)  mol_rat = 4._eb
		
	!Mounting the d's array
	d(1,0:4) = (/ 0.0340429_eb,  0.0652305_eb, -0.0463685_eb,  0.0138684_eb, -0.0014450_eb /)
	d(2,0:4) = (/ 0.3509457_eb,  0.7465138_eb, -0.5293090_eb,  0.1594423_eb, -0.0166326_eb /)
	d(3,0:4) = (/ 4.5707400_eb,  2.1680670_eb, -1.4989010_eb,  0.4917165_eb, -0.0542999_eb /)
	d(4,0:4) = (/ 109.81690_eb, -50.923590_eb,  23.432360_eb, -5.1638920_eb,  0.4393889_eb /)
		
	!Computing the absorption coefficient
	sum_kappa = 0.0
	do ii = 0,4
		sum_kappa = sum_kappa + d(jj,ii)*(mol_rat**real(ii,eb))
	enddo
	get_kj = sum_kappa*(xxw+xxc)			!Assuming atmospheric pressure	

endfunction get_kj

real(eb) function get_aj(ttmp,xxw,xxc,jj)

	integer,intent(in) :: jj
	integer :: ii,kk
	real(eb),intent(in) :: ttmp,xxw,xxc
	real(eb) :: mol_rat,tref,c(1:4,0:4,0:4),sum_b,sum_c
	
	mol_rat = xxw/(xxc + small)
	tref = 1200._eb
	
	!Trim mole ratio to the bounds of the correlations (0.01 - 4)
	if (mol_rat<0.01_eb) mol_rat = 0.01_eb
	if (mol_rat>4._eb)  mol_rat = 4._eb
	
	!Mounting the c's array
	c(1,0,0:4) = (/  0.7412956_eb, -0.5244441_eb,  0.5822860_eb, -0.2096994_eb,  0.0242031_eb /)
	c(1,1,0:4) = (/ -0.9412652_eb,  0.2799577_eb, -0.7672319_eb,  0.3204027_eb, -0.0391017_eb /)
	c(1,2,0:4) = (/  0.8531866_eb,  0.0823075_eb,  0.5289430_eb, -0.2468463_eb,  0.0310940_eb /)
	c(1,3,0:4) = (/ -0.3342806_eb,  0.1474987_eb, -0.4160689_eb,  0.1697627_eb, -0.0204066_eb /)
	c(1,4,0:4) = (/  0.0431436_eb, -0.0688622_eb,  0.1109773_eb, -0.0420861_eb,  0.0049188_eb /)
	c(2,0,0:4) = (/  0.1552073_eb, -0.4862117_eb,  0.3668088_eb, -0.1055508_eb,  0.0105857_eb /)
	c(2,1,0:4) = (/  0.6755648_eb,  1.4092710_eb, -1.3834490_eb,  0.4575210_eb, -0.0501976_eb /)
	c(2,2,0:4) = (/ -1.1253940_eb, -0.5913199_eb,  0.9085441_eb, -0.3334201_eb,  0.0384236_eb /)
	c(2,3,0:4) = (/  0.6040543_eb, -0.0553385_eb, -0.1733014_eb,  0.0791608_eb, -0.0098934_eb /)
	c(2,4,0:4) = (/ -0.1105453_eb,  0.0464663_eb, -0.0016129_eb, -0.0035398_eb,  0.0006121_eb /)
	c(3,0,0:4) = (/  0.2550242_eb,  0.3805403_eb, -0.4249709_eb,  0.1429446_eb, -0.0157408_eb /)
	c(3,1,0:4) = (/ -0.6065428_eb,  0.3494024_eb,  0.1853509_eb, -0.1013694_eb,  0.0130244_eb /)
	c(3,2,0:4) = (/  0.8123855_eb, -1.1020090_eb,  0.4046178_eb, -0.0811822_eb,  0.0062981_eb /)
	c(3,3,0:4) = (/ -0.4532290_eb,  0.6784475_eb, -0.3432603_eb,  0.0883088_eb, -0.0084152_eb /)
	c(3,4,0:4) = (/  0.0869309_eb, -0.1306996_eb,  0.0741446_eb, -0.0202929_eb,  0.0020110_eb /)
	c(4,0,0:4) = (/ -0.0345199_eb,  0.2656726_eb, -0.1225365_eb,  0.0300151_eb, -0.0028205_eb /)
	c(4,1,0:4) = (/  0.4112046_eb, -0.5728350_eb,  0.2924490_eb, -0.0798076_eb,  0.0079966_eb /)
	c(4,2,0:4) = (/ -0.5055995_eb,  0.4579559_eb, -0.2616436_eb,  0.0764841_eb, -0.0079084_eb /)
	c(4,3,0:4) = (/  0.2317509_eb, -0.1656759_eb,  0.1052608_eb, -0.0321935_eb,  0.0033870_eb /)
	c(4,4,0:4) = (/ -0.0375491_eb,  0.0229520_eb, -0.0160047_eb,  0.0050463_eb, -0.0005364_eb /)

	!Computing the temperature coefficient of gray gas j
	sum_b = 0._eb
	do ii=0,4
		sum_c = 0._eb
		do kk=0,4
			sum_c = sum_c + c(jj,ii,kk)*mol_rat**(real(kk,eb))
		enddo
		sum_b = sum_b + sum_c*(ttmp/tref)**(real(ii,eb))
	enddo

endfunction get_aj

endprogram cylinder_model