module GIBBS_1p4 ! objetivo: módulo com as sub-rotinas funcionais do programa Gibbs_1p4 !------------------------------------------------------------------------------- ! Implementação da versão atual (1.4): ! Programador: Carlos Henrique Marchi ! Local: LENA/DEMEC/UFPR, Curitiba, PR ! Período: 24 de junho a 8 de novembro de 2005 ! Última alteração: 24 de fevereiro de 2006 ! Implementação da versão atual (1.3): ! Programador: Carlos Henrique Marchi ! Local: LENA/DEMEC/UFPR, Curitiba, PR ! Período: 14 de janeiro a 24 de fevereiro de 2005 ! Última alteração: 3 de maio de 2005 ! Implementação da versão atual (1.2): ! Programador: Carlos Henrique Marchi ! Local: LENA/DEMEC/UFPR, Curitiba, PR ! Período: 10 a 13 de janeiro de 2005 ! Última alteração: 13 de janeiro de 2005 ! Implementação da versão atual (1.1): ! Programador: Carlos Henrique Marchi ! Local: LENA/DEMEC/UFPR, Curitiba, PR ! Período: 7 a 7 de janeiro de 2005 ! Última alteração: 7 de janeiro de 2005 ! Implementação da versão inicial (1.0): ! Programador: Carlos Henrique Marchi ! Local: LENA/DEMEC/UFPR, Curitiba, PR ! Período: 4 a 6 de janeiro de 2005 !------------------------------------------------------------------------------- ! Arquivos necessários para usar este módulo: ! Gibbs_coeficientes_cp_h_g.dat ! Gibbs_coeficientes_ki_mi.dat !------------------------------------------------------------------------------- use PORTLIB use IFPORT implicit none !------------------------------------------------------------------------------- ! *** constantes *** real*8,parameter,private :: po = 101325.0d+0 ! pressão ao nível do mar (Pa) real*8,parameter,private :: Ru = 8.314510d+0 ! constante universal dos gases (J/mol.K) !------------------------------------------------------------------------------- ! *** variáveis do bloco TERMO *** !------------------------------------------------------------------------------- integer,private :: Ne ! número de espécies químicas integer,private :: i ! número da espécie química ! 1 = H2O 2 = O2 3 = H2 4 = OH ! 5 = O 6 = H 7 = HO2 8 = H2O2 integer,private :: base ! base das propriedades: 1=mássica; 2=molar real*8, private :: T ! temperatura (K) ! variáveis de cada espécie química i character*4,dimension(:),allocatable,private :: especie ! nome da espécie química real*8,dimension(:),allocatable,private :: cpi ! calor específico congelado a ! pressão constante (J/mol.K) ou (J/kg.K) real*8,dimension(:),allocatable,private :: hi ! entalpia (J/mol) ou (J/kg) real*8,dimension(:),allocatable,private :: Moli ! Mol (kg/kmol) real*8,dimension(:),allocatable,private :: mii ! viscosidade absoluta (Pa.s) real*8,dimension(:),allocatable,private :: ki ! condutividade térmica (W/m.K) real*8,dimension(:),allocatable,private :: gi ! energia livre de Gibbs na base ! molar (J/mol) ou (J/kg) ! coeficientes empíricos de NASA TM-4513 (1993) ! primeiro índice = número da espécie química i ! segundo índice = número do coeficiente l ! terceiro índice = faixa de temperatura m real*8,dimension(:,:,:),allocatable,private :: ca ! coeficientes para cpi, hi, gi real*8,dimension(:,:,:),allocatable,private :: cb ! coeficientes para mii real*8,dimension(:,:,:),allocatable,private :: cc ! coeficientes para ki integer,private :: l ! número do coeficiente integer,private :: m ! faixa de temperatura; 1 para T >= 1000 K; 2 para T < 1000 K !------------------------------------------------------------------------------- ! *** variáveis do bloco CONGELADO *** !------------------------------------------------------------------------------- ! variáveis da mistura de Ne espécies químicas real*8,private :: cp ! calor específico congelado a pressão constante (J/kg.K) real*8,private :: Mol ! massa molecular (kg/kmol) real*8,private :: n ! número total final de moles dos produtos (mol) real*8,private :: R ! constante da mistura de gases (J/kg.K) real*8,private :: gama ! razão entre os calores específicos congelados (adim.) real*8,private :: kc ! condutividade térmica (W/m.K) real*8,private :: mi ! viscosidade absoluta (Pa.s) ! variáveis de cada espécie química i real*8,dimension(:),allocatable,private :: ni ! número de moles (mol) real*8,dimension(:),allocatable,private :: Xi ! fração molar (adim.) real*8,dimension(:),allocatable,private :: Yi ! fração mássica (adim.) ! variável auxiliar para cálculo de kc e mi (Ne,Ne) real*8,dimension(:,:),allocatable,private :: fi ! fator molar (adim.) !------------------------------------------------------------------------------- ! *** variáveis do bloco EQUILIBRIO *** !------------------------------------------------------------------------------- integer,private :: modelo ! número do modelo de reações em equilíbrio químico ! modelo = 0: 3 espécies químicas (H2O, O2, H2) ! sem reações de dissociação ! modelo = 1: 3 espécies químicas (H2O, O2, H2) ! com 1 reação de dissociação ! Reação de dissociação: 2*H2 + O2 --> 2*H2O ! modelo = 2: 4 espécies químicas (H2O, O2, H2, OH) ! com 2 reações de dissociação ! Reação de dissociação 1: 2*H2 + O2 --> 2*H2O ! Reação de dissociação 2: H2 + O2 --> 2*OH ! modelo = 3: 6 espécies químicas (H2O, O2, H2, OH, O, H) ! com 4 reações de dissociação ! Reação de dissociação 1: H + OH --> H2O ! Reação de dissociação 2: H + H --> H2 ! Reação de dissociação 3: O + O --> O2 ! Reação de dissociação 4: O + H --> OH ! modelo = 4: 6 espécies químicas (H2O, O2, H2, OH, O, H) ! com 4 reações de dissociação ! Reação de dissociação 1: 2*H2 + O2 --> 2*H2O ! Reação de dissociação 2: H + H --> H2 ! Reação de dissociação 3: O + O --> O2 ! Reação de dissociação 4: O + H --> OH ! modelo = 5: 6 espécies químicas (H2O, O2, H2, OH, O, H) ! com 8 reações de dissociação ! Reação de dissociação 1: H + OH --> H2O ! Reação de dissociação 2: H + H --> H2 ! Reação de dissociação 3: O + O --> O2 ! Reação de dissociação 4: O + H --> OH ! Reação de dissociação 5: O + OH --> H + O2 ! Reação de dissociação 6: H + OH --> O + H2 ! Reação de dissociação 7: H2 + OH --> H + H2O ! Reação de dissociação 8: OH + OH --> O + H2O ! modelo = 6: 6 espécies químicas (H2O, O2, H2, OH, O, H) ! com 4 reações de dissociação ! Reação de dissociação 1: H2 + OH --> H2O + H ! Reação de dissociação 2: OH + OH --> H2O + O ! Reação de dissociação 3: H2 + O --> H + OH ! Reação de dissociação 4: O2 + H --> O + OH ! modelo = 7: 6 espécies químicas (H2O, O2, H2, OH, O, H) ! com 8 reações de dissociação ! Reação de dissociação 1: H2 + OH --> H2O + H ! Reação de dissociação 2: OH + OH --> H2O + O ! Reação de dissociação 3: H2 + O --> H + OH ! Reação de dissociação 4: O2 + H --> O + OH ! Reação de dissociação 5: H + OH --> H2O ! Reação de dissociação 6: O + H --> OH ! Reação de dissociação 7: O + O --> O2 ! Reação de dissociação 8: H + H --> H2 ! modelo = 9: 8 espécies químicas (H2O, O2, H2, OH, O, H, HO2, H2O2) ! com 18 reações de dissociação ! Reação de dissociação 1: H + OH --> H2O ! Reação de dissociação 2: H2 --> H + H ! Reação de dissociação 3: O2 --> O + O ! Reação de dissociação 4: H + O2 --> HO2 ! Reação de dissociação 5: H2O2 --> OH + OH ! Reação de dissociação 6: H2 + O2 --> OH + OH ! Reação de dissociação 7: OH + H2 --> H2O + H ! Reação de dissociação 8: H + O2 --> OH + O ! Reação de dissociação 9: 0 + H2 --> OH + H ! Reação de dissociação 10: H + 2*O2 --> HO2 + O2 ! Reação de dissociação 11: OH + HO2 --> H2O + O2 ! Reação de dissociação 12: H + HO2 --> OH + OH ! Reação de dissociação 13: O + HO2 --> O2 + OH ! Reação de dissociação 14: OH + OH --> O + H2O ! Reação de dissociação 15: H + HO2 --> H2 + O2 ! Reação de dissociação 16: HO2 + HO2 --> H2O2 + O2 ! Reação de dissociação 17: H2O2 + H --> HO2 + H2 ! Reação de dissociação 18: H2O2 + OH --> H2O + HO2 ! modelo = 10: 8 espécies químicas (H2O, O2, H2, OH, O, H, HO2, H2O2) ! com 6 reações de dissociação ! Reação de dissociação 1: H + OH --> H2O ! Reação de dissociação 2: H + H --> H2 ! Reação de dissociação 3: O + O --> O2 ! Reação de dissociação 4: O + H --> OH ! Reação de dissociação 5: H + O2 --> HO2 ! Reação de dissociação 6: H2O2 --> OH + OH ! modelo = 11: 6 espécies químicas (H2O, O2, H2, OH, O, H) ! com 3 reações de dissociação ! Reação de dissociação 1: H + OH --> H2O ! Reação de dissociação 2: H2 --> H + H ! Reação de dissociação 3: O2 --> O + O ! modelo = 12: 8 espécies químicas (H2O, O2, H2, OH, O, H, HO2, H2O2) ! com 5 reações de dissociação ! Reação de dissociação 1: H + OH --> H2O ! Reação de dissociação 2: H2 --> H + H ! Reação de dissociação 3: O2 --> O + O ! Reação de dissociação 4: H + O2 --> HO2 ! Reação de dissociação 5: H2O2 --> OH + OH real*8, private :: p ! pressão (atm) real*8, private :: p_si ! pressão (Pa) real*8, private :: ro ! massa específica (g/cm3) real*8, private :: ro_si ! massa específica (kg/m3) real*8, private :: C ! concentração (mol/cm3) real*8, private :: OF ! razão em massa oxidante/combustível (adim.) real*8, private :: equiv ! razão de equivalência (adim.) real*8, private :: a ! número de moles do reagente O2 (mol) ! número de moles do reagente H2 = 1 mol real*8, private :: b ! número de moles do produto H20 sem dissociação (mol) real*8, private :: d ! número de moles do produto O2 sem dissociação (mol) real*8, private :: f ! número de moles do produto H2 sem dissociação (mol) real*8, private :: no ! número total inicial de moles dos produtos (mol) real*8, private :: tol_e ! tolerância para resolver e(j) real*8, private :: tol_n ! tolerância para resolver n integer,private :: itimax ! número máximo de iterações para resolver e(j) integer,private :: itemax ! número máximo de iterações para resolver n integer,private :: iti ! número da iteração interna para obter e(j) integer,private :: ite ! número da iteração externa para obter n ! variáveis de cada espécie química i real*8,dimension(:),allocatable,private :: Ci ! concentração (mol/cm3) ! variáveis das Nr reações de dissociação integer,private :: Nr ! número de reações químicas de dissociação integer,private :: j ! número da reação química de dissociação real*8,dimension(:),allocatable,private :: dg ! variação dos gi da reação j real*8,dimension(:),allocatable,private :: k ! constante de equilíbrio da ! reação j baseada na pressão parcial real*8,dimension(:),allocatable,private :: e ! taxa de dissociação da reação j real*8,dimension(:),allocatable,private :: de ! variação de e entre duas iterações !------------------------------------------------------------------------------- ! *** variáveis do bloco CÂMARA *** !------------------------------------------------------------------------------- real*8,private :: Tc ! temperatura de combustão (K) integer,private :: itcmax ! número máximo de iterações para resolver Tc real*8,private :: HR ! entalpia total dos reagentes (J) real*8,private :: HP ! entalpia total dos produtos (J) real*8,private :: T_O ! temperatura de injeção do oxidante na câmara (K) real*8,private :: T_F ! temperatura de injeção do combustível na câmara (K) real*8,private :: h_O ! entalpia de injeção do oxidante na câmara (J/mol) real*8,private :: h_F ! entalpia de injeção do combustível na câmara (J/mol) real*8,private :: tol_Tc ! tolerância para obter Tc real*8,private :: dTc ! variação de Tc entre duas iterações integer,private :: itc ! número da iteração para obter Tc integer,private :: tipo_hr ! tipo de HR: 1=automático; 2=prescrito !------------------------------------------------------------------------------- ! *** variáveis do bloco TAXA FINITA *** !------------------------------------------------------------------------------- ! modelo = 5: 6 espécies químicas (H2O, O2, H2, OH, O, H) ! com 8 reações de dissociação ! Reação de dissociação 1: H + OH + M --> H2O + M ! Reação de dissociação 2: H + H + M --> H2 + M ! Reação de dissociação 3: O + O + M --> O2 + M ! Reação de dissociação 4: O + H + M --> OH + M ! Reação de dissociação 5: O + OH --> H + O2 ! Reação de dissociação 6: H + OH --> O + H2 ! Reação de dissociação 7: H2 + OH --> H + H2O ! Reação de dissociação 8: OH + OH --> O + H2O ! Barros et al. (1990) ! modelo = 7: 6 espécies químicas (H2O, O2, H2, OH, O, H) ! com 8 reações de dissociação ! Reação de dissociação 1: H2 + OH --> H2O + H ! Reação de dissociação 2: OH + OH --> H2O + O ! Reação de dissociação 3: H2 + O --> H + OH ! Reação de dissociação 4: O2 + H --> O + OH ! Reação de dissociação 5: H + OH + M --> H2O + M ! Reação de dissociação 6: O + H + M --> OH + M ! Reação de dissociação 7: O + O + M --> O2 + M ! Reação de dissociação 8: H + H + M --> H2 + M ! NASA TP-2725 (1987) ! modelo = 31: 6 espécies químicas (H2O, O2, H2, OH, O, H) ! com 4 reações de dissociação ! Reação de dissociação 1: H + OH + M --> H2O + M ! Reação de dissociação 2: H + H + M --> H2 + M ! Reação de dissociação 3: O + O + M --> O2 + M ! Reação de dissociação 4: O + H + M --> OH + M ! Barros et al. (1990) ! modelo = 32: 6 espécies químicas (H2O, O2, H2, OH, O, H) ! com 4 reações de dissociação ! Reação de dissociação 1: H + OH + M --> H2O + M ! Reação de dissociação 2: H + H + M --> H2 + M ! Reação de dissociação 3: O + O + M --> O2 + M ! Reação de dissociação 4: O + H + M --> OH + M ! NASA TP-2725 (1987) ! modelo = 10: 8 espécies químicas (H2O, O2, H2, OH, O, H, HO2, H2O2) ! com 6 reações de dissociação ! Reação de dissociação 1: H + OH + M --> H2O + M ! Reação de dissociação 2: H + H + M --> H2 + M ! Reação de dissociação 3: O + O + M --> O2 + M ! Reação de dissociação 4: O + H + M --> OH + M ! Reação de dissociação 5: H + O2 + M --> HO2 + M ! Reação de dissociação 6: H2O2 + M --> OH + OH + M ! Reações 1 a 4, de Barros et al. (1990) ! Reações 5 e 6, de CHEMKIN (1990) ! modelo = 9: 8 espécies químicas (H2O, O2, H2, OH, O, H, HO2, H2O2) ! com 18 reações de dissociação ! Reação de dissociação 1: H + OH + M --> H2O + M ! Reação de dissociação 2: H2 + M --> H + H + M ! Reação de dissociação 3: O2 + M --> O + O + M ! Reação de dissociação 4: H + O2 + M --> HO2 + M ! Reação de dissociação 5: H2O2 + M --> OH + OH + M ! Reação de dissociação 6: H2 + O2 --> OH + OH ! Reação de dissociação 7: OH + H2 --> H2O + H ! Reação de dissociação 8: H + O2 --> OH + O ! Reação de dissociação 9: 0 + H2 --> OH + H ! Reação de dissociação 10: H + 2*O2 --> HO2 + O2 ! Reação de dissociação 11: OH + HO2 --> H2O + O2 ! Reação de dissociação 12: H + HO2 --> OH + OH ! Reação de dissociação 13: O + HO2 --> O2 + OH ! Reação de dissociação 14: OH + OH --> O + H2O ! Reação de dissociação 15: H + HO2 --> H2 + O2 ! Reação de dissociação 16: HO2 + HO2 --> H2O2 + O2 ! Reação de dissociação 17: H2O2 + H --> HO2 + H2 ! Reação de dissociação 18: H2O2 + OH --> H2O + HO2 ! CHEMKIN (1990) ! variáveis de cada espécie química i real*8,dimension(:),allocatable,private :: Wi ! taxa de geração de massa por ! unidade de volume (kg/m3.s) ! variáveis das Nr reações de dissociação real*8,dimension(:),allocatable,private :: kf ! taxa de reação no sentido direto real*8,dimension(:),allocatable,private :: kb ! taxa de reação no sentido inverso real*8,dimension(:),allocatable,private :: gam ! fator gama de teta real*8,dimension(:),allocatable,private :: lam ! fator lambda de teta real*8,dimension(:),allocatable,private :: teta ! fator de Wi !------------------------------------------------------------------------------- ! *** variáveis do bloco MACH1D *** !------------------------------------------------------------------------------- ! variáveis da mistura de Ne espécies químicas real*8,private :: cpr ! calor específico reativo a pressão constante (J/kg.K) real*8,private :: gamae ! razão efetiva entre os calores específicos (adim.) ! variáveis de cada espécie química i real*8,dimension(:),allocatable,private :: Yi1 ! fração mássica a T=T+dT/2 (adim.) ! ou fração molar a p=p+dp/2 real*8,dimension(:),allocatable,private :: Yi2 ! fração mássica a T=T-dT/2 (adim.) ! ou fração molar a p=p-dp/2 real*8,dimension(:),allocatable,private :: si ! entropia (J/mol.K) !------------------------------------------------------------------------------- ! *** váriaveis auxiliares (transporte de informações) *** integer,private :: vol ! contador (volume) integer,private :: nvol ! número total de volumes integer,private :: itgeral ! iteração global geral integer,private :: max_it ! número máximo de iterações integer,parameter,private :: arq_ej = 88 ! número do arquivo para escrita dos ej's integer,parameter,private :: arq_ej_a = 89 ! número do arquivo para leitura dos ej's real*8,dimension(:,:),allocatable :: ej ! auxiliar - taxa de dissociação da reação j (guardado) !------------------------------------------------------------------------------- ! *** váriaveis auxiliares *** integer,parameter,private :: in = 7 ! número do arquivo de entrada (dados) integer,private :: out ! número do arquivo de saída principal integer,private :: sai ! =1, escreve variáveis secundárias durante as iterações integer,private :: dos ! acessa prompt dos contains !------------------------------------------------------------------------------- ! ****************************************************************************** ! *** Rotinas do bloco TERMO *** ! ****************************************************************************** ! Objetivo: calcular propriedades termodinâmicas e de transporte de cada espécie ! química, na base mássica ou molar !------------------------------------------------------------------------------- subroutine GIBBS_TERMO_dados_recebe ( outt, Net, baset, Tt ) integer,intent(in) :: outt ! número do arquivo de saída principal integer,intent(in) :: Net ! número de espécies químicas integer,intent(in) :: baset ! base das propriedades: 1=mássica; 2=molar real*8, intent(in) :: Tt ! temperatura out = outt Ne = Net base = baset T = Tt end subroutine GIBBS_TERMO_dados_recebe !------------------------------------------------------------------------------- subroutine GIBBS_TERMO_dados_escreve ( out ) integer,intent(in) :: out ! número do arquivo de saída write(out,1) Ne, T, Ru 1 format(/, t1,'Número de espécies químicas =', i4, & /, t1,'Temperatura (K) =', 1pe22.9, & /, t1,'Constante universal dos gases (J/mol.K) = ', 1pe22.9) end subroutine GIBBS_TERMO_dados_escreve !------------------------------------------------------------------------------- subroutine GIBBS_le_coeficientes_cp_h_g allocate ( ca(Ne,7,2), especie(Ne), Moli(Ne) ) open(in,file='Gibbs_coeficientes_cp_h_g.dat') do i = 1, Ne read(in,*) especie(i) read(in,*) read(in,*) Moli(i) read(in,*) do m = 1, 2 do l = 1, 7 read(in,*) ca(i,l,m) end do read(in,*) end do end do close(in) end subroutine GIBBS_le_coeficientes_cp_h_g !------------------------------------------------------------------------------- subroutine GIBBS_escreve_coeficientes_cp_h_g ( out ) integer,intent(in) :: out ! número do arquivo de saída write(out,1) 1 format(/, t1,'espécie', t16,'massa molecular (kg/kmol)') do i = 1, Ne write(out,2) especie(i), Moli(i) 2 format( t3,a4, t9,1pe22.9 ) end do write(out,3) 3 format(/, 'Coeficientes para calor específico, entalpia e energia livre', & /, 'Para temperatura menor do que 1000 K', & /, t1,'espécie', t10,'a1', t26,'a2', t42,'a3', t58,'a4', t74,'a5', & t90,'a6', t106,'a7' ) do i = 1, Ne write(out,4) especie(i), (ca(i,l,2),l=1,7) 4 format( t3,a4, 1x, 7(1pe16.8) ) end do write(out,5) 5 format('Para temperatura maior ou igual a 1000 K', & /, t1,'espécie', t10,'a1', t26,'a2', t42,'a3', t58,'a4', t74,'a5', & t90,'a6', t106,'a7' ) do i = 1, Ne write(out,4) especie(i), (ca(i,l,1),l=1,7) end do end subroutine GIBBS_escreve_coeficientes_cp_h_g !------------------------------------------------------------------------------- subroutine GIBBS_fecha_memoria_ca_especie_Moli deallocate ( ca, especie, Moli ) end subroutine GIBBS_fecha_memoria_ca_especie_Moli !------------------------------------------------------------------------------- subroutine GIBBS_le_coeficientes_ki_mii allocate ( cb(Ne,4,2), cc(Ne,4,2) ) open(in,file='Gibbs_coeficientes_ki_mi.dat') do i = 1, Ne read(in,*) read(in,*) do m = 1, 2 read(in,*) ( cb(i,l,m), l = 1, 4 ) end do do m = 1, 2 read(in,*) ( cc(i,l,m), l = 1, 4 ) end do read(in,*) end do close(in) end subroutine GIBBS_le_coeficientes_ki_mii !------------------------------------------------------------------------------- subroutine GIBBS_escreve_coeficientes_ki_mii ( out ) integer,intent(in) :: out ! número do arquivo de saída write(out,1) 1 format(/, 'Coeficientes para viscosidade', & /, 'Para temperatura menor do que 1000 K', & /, t1,'espécie', t10,'b1', t26,'b2', t42,'b3', t58,'b4' ) do i = 1, Ne write(out,2) especie(i), (cb(i,l,2),l=1,4) 2 format( t3,a4, 1x, 4(1pe16.8) ) end do write(out,3) 3 format('Para temperatura maior ou igual a 1000 K', & /, t1,'espécie', t10,'b1', t26,'b2', t42,'b3', t58,'b4' ) do i = 1, Ne write(out,2) especie(i), (cb(i,l,1),l=1,4) end do write(out,4) 4 format(/, 'Coeficientes para condutividade térmica', & /, 'Para temperatura menor do que 1000 K', & /, t1,'espécie', t10,'c1', t26,'c2', t42,'c3', t58,'c4' ) do i = 1, Ne write(out,2) especie(i), (cc(i,l,2),l=1,4) end do write(out,5) 5 format('Para temperatura maior ou igual a 1000 K', & /, t1,'espécie', t10,'c1', t26,'c2', t42,'c3', t58,'c4' ) do i = 1, Ne write(out,2) especie(i), (cc(i,l,1),l=1,4) end do end subroutine GIBBS_escreve_coeficientes_ki_mii !------------------------------------------------------------------------------- subroutine GIBBS_fecha_memoria_cb_cc deallocate ( cb, cc ) end subroutine GIBBS_fecha_memoria_cb_cc !------------------------------------------------------------------------------- subroutine GIBBS_cpi_aloca_memoria allocate ( cpi(Ne) ) end subroutine GIBBS_cpi_aloca_memoria !------------------------------------------------------------------------------- subroutine GIBBS_cpi_massa_calculo ! cpi (J/kg.K) call GIBBS_cpi_molar_calculo cpi = cpi * 1.0d+3 / Moli end subroutine GIBBS_cpi_massa_calculo !------------------------------------------------------------------------------- subroutine GIBBS_cpi_molar_calculo ! cpi (J/mol.K) if ( T >= 1.0d+3 ) then m = 1 else m = 2 end if do i = 1, Ne cpi(i) = ( ca(i,1,m) + ca(i,2,m)*T + ca(i,3,m)*(T**2) & + ca(i,4,m)*(T**3) + ca(i,5,m)*(T**4) ) * Ru end do end subroutine GIBBS_cpi_molar_calculo !------------------------------------------------------------------------------- subroutine GIBBS_cpi_escreve ( out ) integer,intent(in) :: out ! número do arquivo de saída if ( base == 1 ) then write(out,1) 1 format(/, t1,'espécie', t16,'calor específico à pressão constante (J/kg.K)') else write(out,2) 2 format(/, t1,'espécie', t16,'calor específico à pressão constante (J/mol.K)') end if do i = 1, Ne write(out,3) especie(i), cpi(i) 3 format( t3,a4, t10,1pe22.9 ) end do end subroutine GIBBS_cpi_escreve !------------------------------------------------------------------------------- subroutine GIBBS_cpi_envia ( cpit ) real*8,dimension(Ne),intent(out) :: cpit cpit = cpi end subroutine GIBBS_cpi_envia !------------------------------------------------------------------------------- subroutine GIBBS_cpi_fecha_memoria deallocate ( cpi ) end subroutine GIBBS_cpi_fecha_memoria !------------------------------------------------------------------------------- subroutine GIBBS_hi_aloca_memoria allocate ( hi(Ne) ) end subroutine GIBBS_hi_aloca_memoria !------------------------------------------------------------------------------- subroutine GIBBS_hi_massa_calculo ! hi (J/kg) call GIBBS_hi_molar_calculo hi = hi * 1.0d+3 / Moli end subroutine GIBBS_hi_massa_calculo !------------------------------------------------------------------------------- subroutine GIBBS_hi_molar_calculo ! hi (J/mol) if ( T >= 1.0d+3 ) then m = 1 else m = 2 end if do i = 1, Ne hi(i) = ( ca(i,1,m) + ca(i,2,m)*T/2.0d0 + ca(i,3,m)*(T**2)/3.0d0 & + ca(i,4,m)*(T**3)/4.0d0 + ca(i,5,m)*(T**4)/5.0d0 & + ca(i,6,m)/T ) * Ru * T end do end subroutine GIBBS_hi_molar_calculo !------------------------------------------------------------------------------- subroutine GIBBS_hi_escreve ( out ) integer,intent(in) :: out ! número do arquivo de saída if ( base == 1 ) then write(out,1) 1 format(/, t1,'espécie', t16,'entalpia (J/kg)') else write(out,2) 2 format(/, t1,'espécie', t16,'entalpia (J/mol)') end if do i = 1, Ne write(out,3) especie(i), hi(i) 3 format( t3,a4, t10,1pe22.9 ) end do end subroutine GIBBS_hi_escreve !------------------------------------------------------------------------------- subroutine GIBBS_hi_envia ( hit ) real*8,dimension(Ne),intent(out) :: hit hit = hi end subroutine GIBBS_hi_envia !------------------------------------------------------------------------------- subroutine GIBBS_hi_fecha_memoria deallocate ( hi ) end subroutine GIBBS_hi_fecha_memoria !------------------------------------------------------------------------------- subroutine GIBBS_gi_aloca_memoria allocate ( gi(Ne) ) end subroutine GIBBS_gi_aloca_memoria !------------------------------------------------------------------------------- subroutine GIBBS_gi_massa_calculo ! gi (J/kg) call GIBBS_gi_molar_calculo gi = gi * 1.0d+3 / Moli end subroutine GIBBS_gi_massa_calculo !------------------------------------------------------------------------------- subroutine GIBBS_gi_molar_calculo ! gi (J/mol) if ( T >= 1.0d+3 ) then m = 1 else m = 2 end if do i = 1, Ne gi(i) = ( ca(i,1,m)*(1.0d0-dlog(T)) - ca(i,2,m)*T/2.0d0 - ca(i,3,m)*(T**2)/6.0d0 & - ca(i,4,m)*(T**3)/12.0d0 - ca(i,5,m)*(T**4)/20.0d0 & + ca(i,6,m)/T - ca(i,7,m) ) * Ru * T end do end subroutine GIBBS_gi_molar_calculo !------------------------------------------------------------------------------- subroutine GIBBS_gi_escreve ( out ) integer,intent(in) :: out ! número do arquivo de saída if ( base == 1 ) then write(out,1) 1 format(/, t1,'espécie', t16,'energia livre de Gibbs (J/kg)') else write(out,2) 2 format(/, t1,'espécie', t16,'energia livre de Gibbs (J/mol)') end if do i = 1, Ne write(out,3) especie(i), gi(i) 3 format( t3,a4, t10,1pe22.9 ) end do end subroutine GIBBS_gi_escreve !------------------------------------------------------------------------------- subroutine GIBBS_gi_envia ( git ) real*8,dimension(Ne),intent(out) :: git git = gi end subroutine GIBBS_gi_envia !------------------------------------------------------------------------------- subroutine GIBBS_gi_fecha_memoria deallocate ( gi ) end subroutine GIBBS_gi_fecha_memoria !------------------------------------------------------------------------------- subroutine GIBBS_mii_aloca_memoria allocate ( mii(Ne) ) end subroutine GIBBS_mii_aloca_memoria !------------------------------------------------------------------------------- subroutine GIBBS_mii_calculo ! mi (Pa.s) if ( T >= 1.0d+3 ) then m = 1 else m = 2 end if do i = 1, Ne mii(i) = dexp( cb(i,1,m)*dlog(T) + cb(i,2,m)/T & + cb(i,3,m)/(T**2) + cb(i,4,m) ) * 1.0d-7 end do end subroutine GIBBS_mii_calculo !------------------------------------------------------------------------------- subroutine GIBBS_mii_escreve ( out ) integer,intent(in) :: out ! número do arquivo de saída write(out,1) 1 format(/, t1,'espécie', t16,'viscosidade absoluta (Pa.s)') do i = 1, Ne write(out,2) especie(i), mii(i) 2 format( t3,a4, t10,1pe22.9 ) end do end subroutine GIBBS_mii_escreve !------------------------------------------------------------------------------- subroutine GIBBS_mii_envia ( miit ) real*8,dimension(Ne),intent(out) :: miit miit = mii end subroutine GIBBS_mii_envia !------------------------------------------------------------------------------- subroutine GIBBS_mii_fecha_memoria deallocate ( mii ) end subroutine GIBBS_mii_fecha_memoria !------------------------------------------------------------------------------- subroutine GIBBS_ki_aloca_memoria allocate ( ki(Ne) ) end subroutine GIBBS_ki_aloca_memoria !------------------------------------------------------------------------------- subroutine GIBBS_ki_calculo ! ki (W/m.K) if ( T >= 1.0d+3 ) then m = 1 else m = 2 end if do i = 1, Ne ki(i) = dexp( cc(i,1,m)*dlog(T) + cc(i,2,m)/T & + cc(i,3,m)/(T**2) + cc(i,4,m) ) * 1.0d-4 end do end subroutine GIBBS_ki_calculo !------------------------------------------------------------------------------- subroutine GIBBS_ki_escreve ( out ) integer,intent(in) :: out ! número do arquivo de saída write(out,1) 1 format(/, t1,'espécie', t16,'condutividade térmica (W/m.K)' ) do i = 1, Ne write(out,2) especie(i), ki(i) 2 format( t3,a4, t10,1pe22.9 ) end do end subroutine GIBBS_ki_escreve !------------------------------------------------------------------------------- subroutine GIBBS_ki_envia ( kit ) real*8,dimension(Ne),intent(out) :: kit kit = ki end subroutine GIBBS_ki_envia !------------------------------------------------------------------------------- subroutine GIBBS_ki_fecha_memoria deallocate ( ki ) end subroutine GIBBS_ki_fecha_memoria !------------------------------------------------------------------------------- subroutine GIBBS_envia_Ru_Moli_especie ( Rut, Molit, especiet ) real*8,intent(out) :: Rut real*8,dimension(Ne),intent(out) :: Molit character*4,dimension(Ne),intent(out) :: especiet Rut = Ru Molit = Moli especiet = especie end subroutine GIBBS_envia_Ru_Moli_especie !------------------------------------------------------------------------------- ! ****************************************************************************** ! *** rotinas do bloco CONGELADO *** ! ****************************************************************************** ! Objetivo: calcular propriedades termodinâmicas e de transporte de uma mistura ! de espécies químicas !------------------------------------------------------------------------------- subroutine GIBBS_CONGELADO_dados_recebe ( outt, Net, Tt, nit ) integer,intent(in) :: outt ! número do arquivo de saída principal integer,intent(in) :: Net ! número de espécies químicas real*8, intent(in) :: Tt ! temperatura real*8,dimension(Net),intent(in) :: nit ! número de moles de cada espécie (mol) out = outt Ne = Net T = Tt call GIBBS_ni_Yi_aloca_memoria ni = nit call GIBBS_le_coeficientes_cp_h_g call GIBBS_cpi_aloca_memoria end subroutine GIBBS_CONGELADO_dados_recebe !------------------------------------------------------------------------------- subroutine GIBBS_ni_recebe ( nit ) real*8,dimension(Ne),intent(in) :: nit ! número de moles de cada espécie (mol) ni = nit end subroutine GIBBS_ni_recebe !------------------------------------------------------------------------------- subroutine GIBBS_ni_escreve ( out ) integer,intent(in) :: out ! número do arquivo de saída write(out,1) 1 format(/, t1,'espécie', t16,'número de moles (mol)' ) do i = 1, Ne write(out,2) especie(i), ni(i) 2 format( t3,a4, t10,1pe22.9 ) end do end subroutine GIBBS_ni_escreve !------------------------------------------------------------------------------- subroutine GIBBS_CONGELADO_dados_escreve ( out ) integer,intent(in) :: out ! número do arquivo de saída write(out,1) Ne, T, Ru 1 format(/, t1,'Número de espécies químicas =', i4, & /, t1,'Temperatura (K) =', 1pe22.9, & /, t1,'Constante universal dos gases (J/mol.K) = ', 1pe22.9 ) call GIBBS_ni_escreve ( out ) end subroutine GIBBS_CONGELADO_dados_escreve !------------------------------------------------------------------------------- subroutine GIBBS_ni_Yi_aloca_memoria allocate ( ni(Ne), Yi(Ne) ) end subroutine GIBBS_ni_Yi_aloca_memoria !------------------------------------------------------------------------------- subroutine GIBBS_ni_Yi_fecha_memoria deallocate ( ni, Yi ) end subroutine GIBBS_ni_Yi_fecha_memoria !------------------------------------------------------------------------------- subroutine GIBBS_n_calculo n = sum(ni) end subroutine GIBBS_n_calculo !------------------------------------------------------------------------------- subroutine GIBBS_Mol_calculo Mol = 0.0d0 do i = 1, Ne Mol = Mol + ni(i) * Moli(i) end do Mol = Mol / n end subroutine GIBBS_Mol_calculo !------------------------------------------------------------------------------- subroutine GIBBS_R_calculo R = 1000 * Ru / Mol end subroutine GIBBS_R_calculo !------------------------------------------------------------------------------- subroutine GIBBS_Yi_calculo do i = 1, Ne Yi(i) = ni(i) * Moli(i) / (n * Mol) end do end subroutine GIBBS_Yi_calculo !------------------------------------------------------------------------------- subroutine GIBBS_Yi_escreve ( out ) integer,intent(in) :: out ! número do arquivo de saída write(out,1) 1 format(/, t1,'espécie', t16,'fração mássica (adim.)' ) do i = 1, Ne write(out,2) especie(i), Yi(i) 2 format( t3,a4, t10,1pe22.9 ) end do end subroutine GIBBS_Yi_escreve !------------------------------------------------------------------------------- subroutine GIBBS_Yi_envia ( Yit ) real*8,dimension(Ne),intent(out) :: Yit Yit = Yi end subroutine GIBBS_Yi_envia !------------------------------------------------------------------------------- subroutine GIBBS_cp_calculo call GIBBS_cpi_massa_calculo cp = 0.0d0 do i = 1, Ne cp = cp + Yi(i) * cpi(i) end do end subroutine GIBBS_cp_calculo !------------------------------------------------------------------------------- subroutine GIBBS_gama_calculo gama = cp / (cp - R) end subroutine GIBBS_gama_calculo !------------------------------------------------------------------------------- subroutine GIBBS_CONGELADO_calculos1 call GIBBS_n_calculo call GIBBS_Mol_calculo call GIBBS_R_calculo call GIBBS_Yi_calculo call GIBBS_cp_calculo call GIBBS_gama_calculo end subroutine GIBBS_CONGELADO_calculos1 !------------------------------------------------------------------------------- subroutine GIBBS_CONGELADO_escreve1 ( out ) integer,intent(in) :: out ! número do arquivo de saída write(out,1) n, Mol, R, cp, gama 1 format(/, t1,'número total de moles (mol) =', 1pe22.9, & /, t1,'massa molecular (kg/kmol) =', 1pe22.9, & /, t1,'constante da mistura de gases (J/kg.K) = ', 1pe22.9, & /, t1,'calor específico congelado à pressão constante (J/kg.K) = ', 1pe22.9, & /, t1,'razão entre os calores específicos congelados (adim.) = ', 1pe22.9 ) call GIBBS_escreve_coeficientes_cp_h_g ( out ) end subroutine GIBBS_CONGELADO_escreve1 !------------------------------------------------------------------------------- subroutine GIBBS_CONGELADO_envia1 ( nt, Molt, Rt, cpt, gamat ) real*8,intent(out) :: nt, Molt, Rt, cpt, gamat nt = n Molt = Mol Rt = R cpt = cp gamat = gama end subroutine GIBBS_CONGELADO_envia1 !------------------------------------------------------------------------------- subroutine GIBBS_mi_kc_inicio call GIBBS_fi_Xi_aloca_memoria call GIBBS_le_coeficientes_ki_mii call GIBBS_mii_aloca_memoria call GIBBS_ki_aloca_memoria end subroutine GIBBS_mi_kc_inicio !------------------------------------------------------------------------------- subroutine GIBBS_fi_Xi_aloca_memoria allocate ( fi(Ne,Ne), Xi(Ne) ) end subroutine GIBBS_fi_Xi_aloca_memoria !------------------------------------------------------------------------------- subroutine GIBBS_fi_Xi_fecha_memoria deallocate ( fi, Xi ) end subroutine GIBBS_fi_Xi_fecha_memoria !------------------------------------------------------------------------------- subroutine GIBBS_Xi_calculo do i = 1, Ne Xi(i) = ni(i) / n end do end subroutine GIBBS_Xi_calculo !------------------------------------------------------------------------------- subroutine GIBBS_Xi_escreve ( out ) integer,intent(in) :: out ! número do arquivo de saída write(out,1) 1 format(/, t1,'espécie', t16,'fração molar (adim.)' ) do i = 1, Ne write(out,2) especie(i), Xi(i) 2 format( t3,a4, t10,1pe22.9 ) end do end subroutine GIBBS_Xi_escreve !------------------------------------------------------------------------------- subroutine GIBBS_Xi_envia ( Xit ) real*8,dimension(Ne),intent(out) :: Xit Xit = Xi end subroutine GIBBS_Xi_envia !------------------------------------------------------------------------------- subroutine GIBBS_CONGELADO_calculos2 integer j ! auxiliar real*8 fat ! auxiliar call GIBBS_mii_calculo call GIBBS_ki_calculo call GIBBS_Xi_calculo ! fi !write(out,1) !1 format(/, t1,'i', t11,'j', t26,'fi' ) do i = 1, Ne do j = 1, Ne fi(i,j) = (( 1 + ((mii(i)/mii(j))**(0.5d0)) * ((Moli(j)/Moli(i))**(0.25d0)) ) ** 2) & / ( ((1 + Moli(i)/Moli(j))**(0.5d0)) * dsqrt(8.0d0) ) !write(out,2) i, j, fi(i,j) !2 format( t3,i4, t13,i4, t20,1pe22.9 ) end do end do ! mi e kc mi = 0.0d0 kc = 0.0d0 do i = 1, Ne fat = 0.0d0 do j = 1, Ne fat = fat + Xi(j)*fi(i,j) end do mi = mi + Xi(i)*mii(i)/fat kc = kc + Xi(i)*ki(i)/fat end do end subroutine GIBBS_CONGELADO_calculos2 !------------------------------------------------------------------------------- subroutine GIBBS_CONGELADO_escreve2 ( out ) integer,intent(in) :: out ! número do arquivo de saída write(out,1) mi, kc 1 format(/, t1,'viscosidade absoluta (Pa.s) =', 1pe22.9, & /, t1,'condutividade térmica (W/m.K) = ', 1pe22.9 ) call GIBBS_escreve_coeficientes_ki_mii ( out ) end subroutine GIBBS_CONGELADO_escreve2 !------------------------------------------------------------------------------- subroutine GIBBS_CONGELADO_envia2 ( mit, kct ) real*8,intent(out) :: mit, kct mit = mi kct = kc end subroutine GIBBS_CONGELADO_envia2 !------------------------------------------------------------------------------- ! ****************************************************************************** ! *** rotinas do bloco EQUILIBRIO *** ! ****************************************************************************** ! Objetivo: calcular a composição de uma mistura de gases !------------------------------------------------------------------------------- subroutine GIBBS_EQUILIBRIO_dados_recebe ( outt, sait, modelot, itimaxt, & itemaxt, tol_et, tol_nt, Tt, & pt, OFt ) integer,intent(in) :: outt, sait, modelot, itimaxt, itemaxt real*8, intent(in) :: tol_et, tol_nt, Tt, pt, OFt out = outt sai = sait modelo = modelot itimax = itimaxt tol_e = tol_et itemax = itemaxt tol_n = tol_nt T = Tt p = pt OF = OFt end subroutine GIBBS_EQUILIBRIO_dados_recebe !------------------------------------------------------------------------------- subroutine GIBBS_EQUILIBRIO_inicio1 select case ( modelo ) case ( 0 ) Ne = 3 Nr = 0 case ( 1 ) Ne = 3 Nr = 1 case ( 2 ) Ne = 4 Nr = 2 case ( 3 ) Ne = 6 Nr = 4 case ( 4 ) Ne = 6 Nr = 4 case ( 5 ) Ne = 6 Nr = 8 case ( 6 ) Ne = 6 Nr = 4 case ( 7 ) Ne = 6 Nr = 8 case ( 9 ) Ne = 8 Nr =18 case (10 ) Ne = 8 Nr = 6 case (11 ) Ne = 6 Nr = 3 case (12 ) Ne = 8 Nr = 5 end select call GIBBS_le_coeficientes_cp_h_g call GIBBS_EQUILIBRIO_aloca_memoria equiv = Moli(2) / (2 * OF * Moli(3)) a = OF * Moli(3) / Moli(2) if ( a >= 0.5d0 ) then b = 1.0d0 d = a - 0.5d0 f = 0.0d0 else b = 2 * a d = 0.0d0 f = 1.0d0 - b end if no = b + d + f call GIBBS_ni_calculo call GIBBS_n_calculo e = 0.0d0 de = 0.0d0 end subroutine GIBBS_EQUILIBRIO_inicio1 !------------------------------------------------------------------------------- subroutine GIBBS_EQUILIBRIO_dados_escreve ( out ) integer,intent(in) :: out ! número do arquivo de saída write(out,1) modelo, Ne, Nr, T, p, OF, equiv, a, Ru, itimax, tol_e, itemax, & tol_n, a, b, d, f, no 1 format(/, t1,'modelo =', i4, & /, t1,'Número de espécies químicas =', i4, & /, t1,'Número de reações =', i4, & /, t1,'Temperatura (K) =', 1pe22.9, & /, t1,'pressão (atm) =', 1pe22.9, & /, t1,'razão O/F mássico (adim.) =', 1pe22.9, & /, t1,'razão de equivalência (adim.) =', 1pe22.9, & /, t1,'razão O/F molar (adim.) =', 1pe22.9, & /, t1,'Constante universal dos gases (J/mol.K) = ', 1pe22.9, & /, t1,'itimax =', i4, & /, t1,'tol_e =', 1pe22.9, & /, t1,'itemax =', i4, & /, t1,'tol_n =', 1pe22.9, & /, t1,'número de moles do reagente O2 (mol) =', 1pe22.9, & /, t1,'número de moles do reagente H2 = 1 mol', & /, t1,'número de moles do produto H2O sem dissociação (mol) =', 1pe22.9, & /, t1,'número de moles do produto O2 sem dissociação (mol) =', 1pe22.9, & /, t1,'número de moles do produto H2 sem dissociação (mol) =', 1pe22.9, & /, t1,'número total inicial de moles dos produtos (mol) =', 1pe22.9 ) end subroutine GIBBS_EQUILIBRIO_dados_escreve !------------------------------------------------------------------------------- subroutine GIBBS_EQUILIBRIO_aloca_memoria allocate ( ni(Ne), Yi(Ne), gi(Ne), cpi(Ne) ) allocate ( k(Nr), dg(Nr), e(Nr), de(Nr) ) end subroutine GIBBS_EQUILIBRIO_aloca_memoria !------------------------------------------------------------------------------- subroutine GIBBS_EQUILIBRIO_fecha_memoria deallocate ( ni, Yi, gi, cpi ) deallocate ( k, dg, e, de ) end subroutine GIBBS_EQUILIBRIO_fecha_memoria !------------------------------------------------------------------------------- subroutine GIBBS_EQUILIBRIO_calculos1 real*8 :: na ! número total de moles dos produtos da iteração anterior call GIBBS_gi_dg_k_calculo if ( sai == 1 ) write(out,1) 1 format(/, t1,'iteração', t18,'variação de n' ) na = no if ((itgeral > 1).and.(itgeral <= max_it)) then ! call GIBBS_MACH1D_abre_ej_a ! call GIBBS_MACH1D_leitura_ej_a ! call GIBBS_MACH1D_fecha_ej_a call GIBBS_MACH1D_compara_e end if do ite = 1, itemax ! cálculo dos e(j) de cada modelo select case ( modelo ) case ( 1 ) call GIBBS_EQUILIBRIO_Modelo_1 case ( 2 ) call GIBBS_EQUILIBRIO_Modelo_2 case ( 3 ) call GIBBS_EQUILIBRIO_Modelo_3 case ( 4 ) call GIBBS_EQUILIBRIO_Modelo_4 case ( 5 ) call GIBBS_EQUILIBRIO_Modelo_5 case ( 6 ) call GIBBS_EQUILIBRIO_Modelo_6 case ( 7 ) call GIBBS_EQUILIBRIO_Modelo_7 case ( 9 ) call GIBBS_EQUILIBRIO_Modelo_9 case (10 ) call GIBBS_EQUILIBRIO_Modelo_10 case (11 ) call GIBBS_EQUILIBRIO_Modelo_11 case (12 ) call GIBBS_EQUILIBRIO_Modelo_12 end select call GIBBS_ni_calculo call GIBBS_n_calculo ! teste de convergência if ( sai == 1 ) write(out,2) ite, dabs(n-na) 2 format( t1,i6, t11,1pe22.9 ) if ( dabs(n-na) <= tol_n .and. ite > 1 ) exit na = n end do ! modificação para escrita de ej ! call GIBBS_MACH1D_escreve_ej if ( (itgeral > 0).and.(itgeral <= max_it) ) call GIBBS_MACH1D_compara_ej call GIBBS_CONGELADO_calculos1 end subroutine GIBBS_EQUILIBRIO_calculos1 !------------------------------------------------------------------------------- subroutine GIBBS_EQUILIBRIO_escreve1 call GIBBS_ej_escreve ( out ) call GIBBS_ni_escreve ( out ) call GIBBS_Yi_escreve ( out ) call GIBBS_CONGELADO_escreve1 ( out ) end subroutine GIBBS_EQUILIBRIO_escreve1 !------------------------------------------------------------------------------- subroutine GIBBS_ej_escreve ( out ) integer,intent(in) :: out ! número do arquivo de saída write(out,1) 1 format(/, t1,'reação', t18,'taxa de dissociação') do j = 1, Nr write(out,2) j, e(j) 2 format( t3,i4, t10,1pe22.9 ) end do end subroutine GIBBS_ej_escreve !------------------------------------------------------------------------------- subroutine GIBBS_EQUILIBRIO_calculos2 call GIBBS_mi_kc_inicio call GIBBS_CONGELADO_calculos2 end subroutine GIBBS_EQUILIBRIO_calculos2 !------------------------------------------------------------------------------- subroutine GIBBS_EQUILIBRIO_escreve2 call GIBBS_CONGELADO_escreve2 ( out ) call GIBBS_Xi_escreve ( out ) end subroutine GIBBS_EQUILIBRIO_escreve2 !------------------------------------------------------------------------------- subroutine GIBBS_Ci_aloca_memoria allocate ( Ci(Ne) ) end subroutine GIBBS_Ci_aloca_memoria !------------------------------------------------------------------------------- subroutine GIBBS_Ci_fecha_memoria deallocate ( Ci ) end subroutine GIBBS_Ci_fecha_memoria !------------------------------------------------------------------------------- subroutine GIBBS_Ci_calculo do i = 1, Ne Ci(i) = ro * Yi(i) / Moli(i) end do end subroutine GIBBS_Ci_calculo !------------------------------------------------------------------------------- subroutine GIBBS_Ci_escreve ( out ) integer,intent(in) :: out ! número do arquivo de saída write(out,1) 1 format(/, t1,'espécie', t16,'Concentração (mol/cm3)' ) do i = 1, Ne write(out,2) especie(i), Ci(i) 2 format( t3,a4, t10,1pe22.9 ) end do end subroutine GIBBS_Ci_escreve !------------------------------------------------------------------------------- subroutine GIBBS_Ci_envia ( Cit ) real*8,dimension(Ne),intent(out) :: Cit Cit = Ci end subroutine GIBBS_Ci_envia !------------------------------------------------------------------------------- subroutine GIBBS_C_calculo C = sum(Ci) end subroutine GIBBS_C_calculo !------------------------------------------------------------------------------- subroutine GIBBS_EQUILIBRIO_calculos3 p_si = p * po ro_si = p_si / ( R * T ) ro = 1.0d-3 * ro_si call GIBBS_Ci_calculo call GIBBS_C_calculo end subroutine GIBBS_EQUILIBRIO_calculos3 !------------------------------------------------------------------------------- subroutine GIBBS_EQUILIBRIO_escreve3 write(out,1) po, p_si, ro, ro_si, C 1 format(/, t1,'pressão ao nível do mar (Pa) =', 1pe22.9, & /, t1,'pressão (Pa) =', 1pe22.9, & /, t1,'massa específica (g/cm3) =', 1pe22.9, & /, t1,'massa específica (kg/m3) =', 1pe22.9, & /, t1,'concentração (mol/cm3) =', 1pe22.9 ) call GIBBS_Ci_escreve ( out ) end subroutine GIBBS_EQUILIBRIO_escreve3 !------------------------------------------------------------------------------- subroutine GIBBS_gi_dg_k_calculo ! cálculo da energia livre de Gibbs de cada espécie (molar) call GIBBS_gi_molar_calculo ! cálculo dos dg(j) de cada modelo select case ( modelo ) case ( 1 ) dg(1) = 2 * gi(1) - gi(2) - 2 * gi(3) case ( 2 ) dg(1) = 2 * gi(1) - gi(2) - 2 * gi(3) dg(2) = 2 * gi(4) - gi(2) - gi(3) case ( 3 ) dg(1) = gi(1) - gi(4) - gi(6) dg(2) = gi(3) - 2 * gi(6) dg(3) = gi(2) - 2 * gi(5) dg(4) = gi(4) - gi(5) - gi(6) case ( 4 ) dg(1) = 2 * gi(1) - gi(2) - 2 * gi(3) dg(2) = gi(3) - 2 * gi(6) dg(3) = gi(2) - 2 * gi(5) dg(4) = gi(4) - gi(5) - gi(6) case ( 5 ) dg(1) = gi(1) - gi(4) - gi(6) dg(2) = gi(3) - 2 * gi(6) dg(3) = gi(2) - 2 * gi(5) dg(4) = gi(4) - gi(5) - gi(6) dg(5) = gi(2) + gi(6) - gi(4) - gi(5) dg(6) = gi(3) + gi(5) - gi(4) - gi(6) dg(7) = gi(1) + gi(6) - gi(3) - gi(4) dg(8) = gi(1) + gi(5) - 2 * gi(4) case ( 6 ) dg(1) = gi(1) + gi(6) - gi(3) - gi(4) dg(2) = gi(1) + gi(5) - 2 * gi(4) dg(3) = gi(4) + gi(6) - gi(3) - gi(5) dg(4) = gi(4) + gi(5) - gi(2) - gi(6) case ( 7 ) dg(1) = gi(1) + gi(6) - gi(3) - gi(4) dg(2) = gi(1) + gi(5) - 2 * gi(4) dg(3) = gi(4) + gi(6) - gi(3) - gi(5) dg(4) = gi(4) + gi(5) - gi(2) - gi(6) dg(5) = gi(1) - gi(4) - gi(6) dg(6) = gi(4) - gi(5) - gi(6) dg(7) = gi(2) - 2 * gi(5) dg(8) = gi(3) - 2 * gi(6) case ( 9 ) dg(1) = gi(1) - gi(4) - gi(6) dg(2) = 2 * gi(6) - gi(3) dg(3) = 2 * gi(5) - gi(2) dg(4) = gi(7) - gi(2) - gi(6) dg(5) = 2 * gi(4) - gi(8) dg(6) = 2 * gi(4) - gi(2) - gi(3) dg(7) = gi(1) + gi(6) - gi(3) - gi(4) dg(8) = gi(4) + gi(5) - gi(2) - gi(6) dg(9) = gi(4) + gi(6) - gi(3) - gi(5) dg(10)= gi(2) + gi(7) - 2 * gi(2) - gi(6) dg(11)= gi(1) + gi(2) - gi(4) - gi(7) dg(12)= 2 * gi(4) - gi(6) - gi(7) dg(13)= gi(2) + gi(4) - gi(5) - gi(7) dg(14)= gi(1) + gi(5) - 2 * gi(4) dg(15)= gi(2) + gi(3) - gi(6) - gi(7) dg(16)= gi(2) + gi(8) - 2 * gi(7) dg(17)= gi(3) + gi(7) - gi(6) - gi(8) dg(18)= gi(1) + gi(7) - gi(4) - gi(8) case (10 ) dg(1) = gi(1) - gi(4) - gi(6) dg(2) = gi(3) - 2 * gi(6) dg(3) = gi(2) - 2 * gi(5) dg(4) = gi(4) - gi(5) - gi(6) dg(5) = gi(7) - gi(2) - gi(6) dg(6) = 2 * gi(4) - gi(8) case (11 ) dg(1) = gi(1) - gi(4) - gi(6) dg(2) = 2 * gi(6) - gi(3) dg(3) = 2 * gi(5) - gi(2) case (12 ) dg(1) = gi(1) - gi(4) - gi(6) dg(2) = 2 * gi(6) - gi(3) dg(3) = 2 * gi(5) - gi(2) dg(4) = gi(7) - gi(2) - gi(6) dg(5) = 2 * gi(4) - gi(8) end select ! cálculo das constantes de equilíbrio baseadas nas pressões parciais if ( sai == 1 ) write(out,1) 1 format(/, t1,'reação', t18,'delta Gibbs', t40,'kp') do j = 1, Nr k(j) = dexp ( - dg(j) / ( Ru * T ) ) if ( sai == 1 ) write(out,2) j, dg(j), k(j) 2 format( t1,i6, t11,1pe22.9, t33,1pe22.9 ) end do end subroutine GIBBS_gi_dg_k_calculo !------------------------------------------------------------------------------- subroutine GIBBS_ni_calculo select case ( modelo ) case ( 0 ) ni(1) = b ni(2) = d ni(3) = f case ( 1 ) ni(1) = b + 2*e(1) ni(2) = d - e(1) ni(3) = f - 2*e(1) case ( 2 ) ni(1) = b + 2*e(1) ni(2) = d - e(1) - e(2) ni(3) = f - 2*e(1) - e(2) ni(4) = 2*e(2) case ( 3 ) ni(1) = b + e(1) ni(2) = d + e(3) ni(3) = f + e(2) ni(4) = - e(1) + e(4) ni(5) = - 2*e(3) - e(4) ni(6) = - e(1) - 2*e(2) - e(4) case ( 4 ) ni(1) = b + 2*e(1) ni(2) = d - e(1) + e(3) ni(3) = f - 2*e(1) + e(2) ni(4) = e(4) ni(5) = - 2*e(3) - e(4) ni(6) = - 2*e(2) - e(4) case ( 5 ) ni(1) = b + e(1) + e(7) + e(8) ni(2) = d + e(3) + e(5) ni(3) = f + e(2) + e(6) - e(7) ni(4) = - e(1) + e(4) - e(5) - e(6) - e(7) - 2*e(8) ni(5) = - 2*e(3) - e(4) - e(5) + e(6) + e(8) ni(6) = - e(1) - 2*e(2) - e(4) + e(5) - e(6) + e(7) case ( 6 ) ni(1) = b + e(1) + e(2) ni(2) = d - e(4) ni(3) = f - e(1) - e(3) ni(4) = - e(1) - 2*e(2) + e(3) + e(4) ni(5) = e(2) - e(3) + e(4) ni(6) = e(1) + e(3) - e(4) case ( 7 ) ni(1) = b + e(1) + e(2) + e(5) ni(2) = d - e(4) + e(7) ni(3) = f - e(1) - e(3) + e(8) ni(4) = - e(1) - 2*e(2) + e(3) + e(4) - e(5) + e(6) ni(5) = e(2) - e(3) + e(4) - e(6) - 2*e(7) ni(6) = e(1) + e(3) - e(4) - e(5) - e(6) - 2*e(8) case ( 9 ) ni(1) = b + e( 1) + e( 7) + e(11) + e(14) + e(18) ni(2) = d - e( 3) - e( 4) - e( 6) - e( 8) - e(10) & + e(11) + e(13) + e(15) + e(16) ni(3) = f - e( 2) - e( 6) - e( 7) - e( 9) + e(15) + e(17) ni(4) = - e( 1) + 2*e( 5) + 2*e( 6) - e( 7) + e( 8) + e( 9) & - e(11) + 2*e(12) + e(13) - 2*e(14) - e(18) ni(5) = 2*e( 3) + e( 8) - e( 9) - e(13) + e(14) ni(6) = - e( 1) + 2*e( 2) - e( 4) + e( 7) - e( 8) & + e( 9) - e(10) - e(12) - e(15) - e(17) ni(7) = e( 4) + e(10) - e(11) - e(12) - e(13) & - e(15) - 2*e(16) + e(17) + e(18) ni(8) = - e( 5) + e(16) - e(17) - e(18) case (10 ) ni(1) = b + e(1) ni(2) = d + e(3) - e(5) ni(3) = f + e(2) ni(4) = - e(1) + e(4) + 2*e(6) ni(5) = - 2*e(3) - e(4) ni(6) = - e(1) - 2*e(2) - e(4) - e(5) ni(7) = e(5) ni(8) = - e(6) case (11 ) ni(1) = b + e( 1) ni(2) = d - e( 3) ni(3) = f - e( 2) ni(4) = - e( 1) ni(5) = 2*e( 3) ni(6) = - e( 1) + 2*e( 2) case (12 ) ni(1) = b + e( 1) ni(2) = d - e( 3) - e( 4) ni(3) = f - e( 2) ni(4) = - e( 1) + 2*e( 5) ni(5) = 2*e( 3) ni(6) = - e( 1) + 2*e( 2) - e( 4) ni(7) = e( 4) ni(8) = - e( 5) end select end subroutine GIBBS_ni_calculo !------------------------------------------------------------------------------- subroutine GIBBS_EQUILIBRIO_Modelo_1 real*8 :: c1, c2, c3, c4, u, Fe, Fl ! solução de e(1) u = k(1)*p c1 = - 4*u c2 = 4*(u*(d+f)-n) c3 = - ( u*f*(4*d+f)+4*b*n ) c4 = (u*d*(f**2)) - ((b**2)*n) do iti = 1, itimax Fe = c1*(e(1)**3) + c2*(e(1)**2) + c3*e(1) + c4 Fl = 3*c1*(e(1)**2) + 2*c2*e(1) + c3 de(1) = - Fe / Fl e(1) = e(1) + de(1) if ( dabs(de(1)) <= tol_e ) exit end do if ( sai == 1 ) then write(out,1) iti 1 format('iti(j) = ', 1(i6)) write(out,2) de 2 format('de(j) = ', 1(1pe22.9)) write(out,3) e 3 format('e(j) = ', 1(1pe22.9)) end if end subroutine GIBBS_EQUILIBRIO_Modelo_1 !------------------------------------------------------------------------------- subroutine GIBBS_EQUILIBRIO_Modelo_2 real*8 :: c1, c2, c3, c4, q2, q3, u, Fe, Fl integer :: iti1, iti2 ! solução de e(1) u = k(1)*p q2 = d - e(2) q3 = f - e(2) c1 = - 4*u c2 = 4*(u*(q2+q3)-n) c3 = - ( u*q3*(4*q2+q3)+4*b*n ) c4 = (u*q2*((q3)**2)) - ((b**2)*n) do iti = 1, itimax Fe = c1*(e(1)**3) + c2*(e(1)**2) + c3*e(1) + c4 Fl = 3*c1*(e(1)**2) + 2*c2*e(1) + c3 de(1) = - Fe / Fl e(1) = e(1) + de(1) if ( dabs(de(1)) <= tol_e ) exit end do iti1 = iti ! solução de e(2) q2 = d - e(1) q3 = f - 2*e(1) c1 = k(2) - 4 c2 = -k(2)*(q2+q3) c3 = k(2)*q2*q3 do iti = 1, itimax Fe = c1*(e(2)**2) + c2*e(2) + c3 Fl = 2*c1*e(2) + c2 de(2) = - Fe / Fl e(2) = e(2) + de(2) if ( dabs(de(2)) <= tol_e ) exit end do iti2 = iti if ( sai == 1 ) then write(out,1) iti1, iti2 1 format('iti(j) = ', 2(i6)) write(out,2) de 2 format('de(j) = ', 2(1pe22.9)) write(out,3) e 3 format('e(j) = ', 2(1pe22.9)) end if end subroutine GIBBS_EQUILIBRIO_Modelo_2 !------------------------------------------------------------------------------- subroutine GIBBS_EQUILIBRIO_Modelo_3 real*8 :: c1, c2, c3, q4, q5, q6, u, Fe, Fl integer :: iti1, iti2, iti3, iti4 ! solução de e(1) u = k(1)*p q4 = e(4) q6 = - 2*e(2) - e(4) c1 = u c2 = - (u*(q4+q6)+n) c3 = u*q4*q6-b*n do iti = 1, itimax Fe = c1*(e(1)**2) + c2*e(1) + c3 Fl = 2*c1*e(1) + c2 de(1) = - Fe / Fl e(1) = e(1) + de(1) if ( dabs(de(1)) <= tol_e ) exit end do iti1 = iti ! solução de e(2) u = k(2)*p q6 = - e(1) - e(4) c1 = 4*u c2 = -(4*u*q6+n) c3 = u*((q6)**2) - f*n do iti = 1, itimax Fe = c1*(e(2)**2) + c2*e(2) + c3 Fl = 2*c1*e(2) + c2 de(2) = - Fe / Fl e(2) = e(2) + de(2) if ( dabs(de(2)) <= tol_e ) exit end do iti2 = iti ! solução de e(3) u = k(3)*p q5 = - e(4) c1 = 4*u c2 = - (4*u*q5+n) c3 = u*((q5)**2) - d*n do iti = 1, itimax Fe = c1*(e(3)**2) + c2*e(3) + c3 Fl = 2*c1*e(3) + c2 de(3) = - Fe / Fl e(3) = e(3) + de(3) if ( dabs(de(3)) <= tol_e ) exit end do iti3 = iti ! solução de e(4) u = k(4)*p q4 = - e(1) q5 = - 2*e(3) q6 = - e(1) - 2*e(2) c1 = u c2 = -(u*(q5+q6)+n) c3 = u*q5*q6 - q4*n do iti = 1, itimax Fe = c1*(e(4)**2) + c2*e(4) + c3 Fl = 2*c1*e(4) + c2 de(4) = - Fe / Fl e(4) = e(4) + de(4) if ( dabs(de(4)) <= tol_e ) exit end do iti4 = iti if ( sai == 1 ) then write(out,1) iti1, iti2, iti3, iti4 1 format('iti(j) = ', 4(i6)) write(out,2) de 2 format('de(j) = ', 4(1pe22.9)) write(out,3) e 3 format('e(j) = ', 4(1pe22.9)) end if end subroutine GIBBS_EQUILIBRIO_Modelo_3 !------------------------------------------------------------------------------- subroutine GIBBS_EQUILIBRIO_Modelo_4 real*8 :: c1, c2, c3, c4, q2,q3, q5, q6, u, Fe, Fl integer :: iti1, iti2, iti3, iti4 ! solução de e(1) u = k(1)*p q2 = d + e(3) q3 = f + e(2) c1 = - 4*u c2 = 4*(u*(q2+q3)-n) c3 = - ( u*q3*(4*q2+q3)+4*b*n ) c4 = (u*q2*((q3)**2)) - ((b**2)*n) do iti = 1, itimax Fe = c1*(e(1)**3) + c2*(e(1)**2) + c3*e(1) + c4 Fl = 3*c1*(e(1)**2) + 2*c2*e(1) + c3 de(1) = - Fe / Fl e(1) = e(1) + de(1) if ( dabs(de(1)) <= tol_e ) exit end do iti1 = iti ! solução de e(2) u = k(2)*p q3 = f - 2*e(1) q6 = - e(4) c1 = 4*u c2 = -(4*u*q6+n) c3 = u*((q6)**2) - q3*n do iti = 1, itimax Fe = c1*(e(2)**2) + c2*e(2) + c3 Fl = 2*c1*e(2) + c2 de(2) = - Fe / Fl e(2) = e(2) + de(2) if ( dabs(de(2)) <= tol_e ) exit end do iti2 = iti ! solução de e(3) u = k(3)*p q2 = d - e(1) q5 = - e(4) c1 = 4*u c2 = - (4*u*q5+n) c3 = u*((q5)**2) - q2*n do iti = 1, itimax Fe = c1*(e(3)**2) + c2*e(3) + c3 Fl = 2*c1*e(3) + c2 de(3) = - Fe / Fl e(3) = e(3) + de(3) if ( dabs(de(3)) <= tol_e ) exit end do iti3 = iti ! solução de e(4) u = k(4)*p q5 = - 2*e(3) q6 = - 2*e(2) c1 = u c2 = -(u*(q5+q6)+n) c3 = u*q5*q6 do iti = 1, itimax Fe = c1*(e(4)**2) + c2*e(4) + c3 Fl = 2*c1*e(4) + c2 de(4) = - Fe / Fl e(4) = e(4) + de(4) if ( dabs(de(4)) <= tol_e ) exit end do iti4 = iti if ( sai == 1 ) then write(out,1) iti1, iti2, iti3, iti4 1 format('iti(j) = ', 4(i6)) write(out,2) de 2 format('de(j) = ', 4(1pe22.9)) write(out,3) e 3 format('e(j) = ', 4(1pe22.9)) end if end subroutine GIBBS_EQUILIBRIO_Modelo_4 !------------------------------------------------------------------------------- subroutine GIBBS_EQUILIBRIO_Modelo_5 real*8 :: c1, c2, c3, q1, q2, q3, q4, q5, q6, u, Fe, Fl integer :: iti1, iti2, iti3, iti4, iti5, iti6, iti7, iti8 ! solução de e(1) u = k(1)*p q1 = b + e(7) + e(8) q4 = e(4) - e(5) - e(6) - e(7) - 2*e(8) q6 = - 2*e(2) - e(4) + e(5) - e(6) + e(7) c1 = u c2 = - (u*(q4+q6)+n) c3 = u*q4*q6-q1*n do iti = 1, itimax Fe = c1*(e(1)**2) + c2*e(1) + c3 Fl = 2*c1*e(1) + c2 de(1) = - Fe / Fl e(1) = e(1) + de(1) if ( dabs(de(1)) <= tol_e ) exit end do iti1 = iti ! solução de e(2) u = k(2)*p q3 = f + e(6) - e(7) q6 = - e(1) - e(4) + e(5) - e(6) + e(7) c1 = 4*u c2 = -(4*u*q6+n) c3 = u*((q6)**2) - q3*n do iti = 1, itimax Fe = c1*(e(2)**2) + c2*e(2) + c3 Fl = 2*c1*e(2) + c2 de(2) = - Fe / Fl e(2) = e(2) + de(2) if ( dabs(de(2)) <= tol_e ) exit end do iti2 = iti ! solução de e(3) u = k(3)*p q2 = d + e(5) q5 = - e(4) - e(5) + e(6) + e(8) c1 = 4*u c2 = - (4*u*q5+n) c3 = u*((q5)**2) - q2*n do iti = 1, itimax Fe = c1*(e(3)**2) + c2*e(3) + c3 Fl = 2*c1*e(3) + c2 de(3) = - Fe / Fl e(3) = e(3) + de(3) if ( dabs(de(3)) <= tol_e ) exit end do iti3 = iti ! solução de e(4) u = k(4)*p q4 = - e(1) - e(5) - e(6) - e(7) - 2*e(8) q5 = - 2*e(3) - e(5) + e(6) + e(8) q6 = - e(1) - 2*e(2) + e(5) - e(6) + e(7) c1 = u c2 = -(u*(q5+q6)+n) c3 = u*q5*q6 - q4*n do iti = 1, itimax Fe = c1*(e(4)**2) + c2*e(4) + c3 Fl = 2*c1*e(4) + c2 de(4) = - Fe / Fl e(4) = e(4) + de(4) if ( dabs(de(4)) <= tol_e ) exit end do iti4 = iti ! solução de e(5) u = k(5) q2 = d + e(3) q4 = - e(1) + e(4) - e(6) - e(7) - 2*e(8) q5 = - 2*e(3) - e(4) + e(6) + e(8) q6 = - e(1) - 2*e(2) - e(4) - e(6) + e(7) c1 = u - 1.0d0 c2 = - (u*(q4+q5)+q2+q6) c3 = u*q4*q5-q2*q6 do iti = 1, itimax Fe = c1*(e(5)**2) + c2*e(5) + c3 Fl = 2*c1*e(5) + c2 de(5) = - Fe / Fl e(5) = e(5) + de(5) if ( dabs(de(5)) <= tol_e ) exit end do iti5 = iti ! solução de e(6) u = k(6) q3 = f + e(2) - e(7) q4 = - e(1) + e(4) - e(5) - e(7) - 2*e(8) q5 = - 2*e(3) - e(4) - e(5) + e(8) q6 = - e(1) - 2*e(2) - e(4) + e(5) + e(7) c1 = u - 1.0d0 c2 = -(u*(q4+q6)+q3+q5) c3 = u*q4*q6 - q3*q5 do iti = 1, itimax Fe = c1*(e(6)**2) + c2*e(6) + c3 Fl = 2*c1*e(6) + c2 de(6) = - Fe / Fl e(6) = e(6) + de(6) if ( dabs(de(6)) <= tol_e ) exit end do iti6 = iti ! solução de e(7) u = k(7) q1 = b + e(1) + e(8) q3 = f + e(2) + e(6) q4 = - e(1) + e(4) - e(5) - e(6) - 2*e(8) q6 = - e(1) - 2*e(2) - e(4) + e(5) - e(6) c1 = u - 1.0d0 c2 = -(u*(q3+q4)+q1+q6) c3 = u*q3*q4 - q1*q6 do iti = 1, itimax Fe = c1*(e(7)**2) + c2*e(7) + c3 Fl = 2*c1*e(7) + c2 de(7) = - Fe / Fl e(7) = e(7) + de(7) if ( dabs(de(7)) <= tol_e ) exit end do iti7 = iti ! solução de e(8) u = k(8) q1 = b + e(1) + e(7) q4 = - e(1) + e(4) - e(5) - e(6) - e(7) q5 = - 2*e(3) - e(4) - e(5) + e(6) c1 = 4*u - 1.0d0 c2 = -(4*u*q4+q1+q5) c3 = u*((q4)**2) - q1*q5 do iti = 1, itimax Fe = c1*(e(8)**2) + c2*e(8) + c3 Fl = 2*c1*e(8) + c2 de(8) = - Fe / Fl e(8) = e(8) + de(8) if ( dabs(de(8)) <= tol_e ) exit end do iti8 = iti if ( sai == 1 ) then write(out,1) iti1, iti2, iti3, iti4, iti5, iti6, iti7, iti8 1 format('iti(j) = ', 8(i6)) write(out,2) de 2 format('de(j) = ', 8(1pe22.9)) write(out,3) e 3 format('e(j) = ', 8(1pe22.9)) end if end subroutine GIBBS_EQUILIBRIO_Modelo_5 !------------------------------------------------------------------------------- subroutine GIBBS_EQUILIBRIO_Modelo_6 real*8 :: c1, c2, c3, q1, q2,q3, q4, q5, q6, u, Fe, Fl ! solução de e(1) e(1) = 1.0d-3 u = k(1) q1 = b + e(2) q3 = f - e(3) q4 = - 2*e(2) + e(3) + e(4) q6 = e(3) - e(4) c1 = u - 1.0d0 c2 = -4*(u*(q3+q4)+q1+q6) c3 = u*q3*q4 - q1*q6 do iti = 1, itimax Fe = c1*(e(1)**2) + c2*e(1) + c3 Fl = 2*c1*e(1) + c2 de(1) = - Fe / Fl e(1) = e(1) + de(1) if ( dabs(de(1)) <= tol_e ) exit end do ! solução de e(2) u = k(2) q1 = b + e(1) q4 = - e(1) + e(3) + e(4) q5 = - e(3) + e(4) c1 = 4*u - 1.0d0 c2 = -(4*u*q4+q1+q5) c3 = u*((q4)**2) - q1*q5 do iti = 1, itimax Fe = c1*(e(2)**2) + c2*e(2) + c3 Fl = 2*c1*e(2) + c2 de(2) = - Fe / Fl e(2) = e(2) + de(2) if ( dabs(de(2)) <= tol_e ) exit end do ! solução de e(3) u = k(3) q3 = f - e(1) q4 = - e(1) - 2*e(2) + e(4) q5 = e(2) + e(4) q6 = e(1) - e(4) c1 = u - 1.0d0 c2 = - (u*(q3+q5)+q4+q6) c3 = u*q3*q5 - q4*q6 do iti = 1, itimax Fe = c1*(e(3)**2) + c2*e(3) + c3 Fl = 2*c1*e(3) + c2 de(3) = - Fe / Fl e(3) = e(3) + de(3) if ( dabs(de(3)) <= tol_e ) exit end do ! solução de e(4) u = k(4) q2 = d q4 = - e(1) - 2*e(2) + e(3) q5 = e(2) - e(3) q6 = e(1) + e(3) c1 = u - 1.0d0 c2 = -(u*(q2+q6)+q4+q5) c3 = u*q2*q6 - q4*q5 do iti = 1, itimax Fe = c1*(e(4)**2) + c2*e(4) + c3 Fl = 2*c1*e(4) + c2 de(4) = - Fe / Fl e(4) = e(4) + de(4) if ( dabs(de(4)) <= tol_e ) exit end do end subroutine GIBBS_EQUILIBRIO_Modelo_6 !------------------------------------------------------------------------------- subroutine GIBBS_EQUILIBRIO_Modelo_7 real*8 :: c1, c2, c3, q1, q2,q3, q4, q5, q6, u, Fe, Fl integer :: iti1, iti2, iti3, iti4, iti5, iti6, iti7, iti8 ! solução de e(5) u = k(5)*p q1 = b + e(1) + e(2) q4 = - e(1) - 2*e(2) + e(3) + e(4) + e(6) q6 = e(1) + e(3) - e(4) - e(6) - 2*e(8) c1 = u c2 = - (u*(q4+q6)+n) c3 = u*q4*q6-q1*n do iti = 1, itimax Fe = c1*(e(5)**2) + c2*e(5) + c3 Fl = 2*c1*e(5) + c2 de(5) = - Fe / Fl e(5) = e(5) + de(5) if ( dabs(de(5)) <= tol_e ) exit end do iti5 = iti ! solução de e(8) u = k(8)*p q3 = f - e(1) - e(3) q6 = e(1) + e(3) - e(4) - e(5) - e(6) c1 = 4*u c2 = -(4*u*q6+n) c3 = u*((q6)**2) - q3*n do iti = 1, itimax Fe = c1*(e(8)**2) + c2*e(8) + c3 Fl = 2*c1*e(8) + c2 de(8) = - Fe / Fl e(8) = e(8) + de(8) if ( dabs(de(8)) <= tol_e ) exit end do iti8 = iti ! solução de e(7) u = k(7)*p q2 = d - e(4) q5 = e(2) - e(3) + e(4) - e(6) c1 = 4*u c2 = - (4*u*q5+n) c3 = u*((q5)**2) - q2*n do iti = 1, itimax Fe = c1*(e(7)**2) + c2*e(7) + c3 Fl = 2*c1*e(7) + c2 de(7) = - Fe / Fl e(7) = e(7) + de(7) if ( dabs(de(7)) <= tol_e ) exit end do iti7 = iti ! solução de e(6) u = k(6)*p q4 = - e(1) - 2*e(2) + e(3) + e(4) - e(5) q5 = e(2) - e(3) + e(4) - 2*e(7) q6 = e(1) + e(3) - e(4) - e(5) - 2*e(8) c1 = u c2 = -(u*(q5+q6)+n) c3 = u*q5*q6 - q4*n do iti = 1, itimax Fe = c1*(e(6)**2) + c2*e(6) + c3 Fl = 2*c1*e(6) + c2 de(6) = - Fe / Fl e(6) = e(6) + de(6) if ( dabs(de(6)) <= tol_e ) exit end do iti6 = iti ! solução de e(1) u = k(1) q1 = b + e(2) + e(5) q3 = f - e(3) + e(8) q4 = - 2*e(2) + e(3) + e(4) - e(5) + e(6) q6 = + e(3) - e(4) - e(5) - e(6) - 2*e(8) c1 = u - 1.0d0 c2 = -(u*(q3+q4)+q1+q6) c3 = u*q3*q4 - q1*q6 do iti = 1, itimax Fe = c1*(e(1)**2) + c2*e(1) + c3 Fl = 2*c1*e(1) + c2 de(1) = - Fe / Fl e(1) = e(1) + de(1) if ( dabs(de(1)) <= tol_e ) exit end do iti1 = iti ! solução de e(2) u = k(2) q1 = b + e(1) + e(5) q4 = - e(1) + e(3) + e(4) - e(5) + e(6) q5 = - e(3) + e(4) - e(6) - 2*e(7) c1 = 4*u - 1.0d0 c2 = -(4*u*q4+q1+q5) c3 = u*((q4)**2) - q1*q5 do iti = 1, itimax Fe = c1*(e(2)**2) + c2*e(2) + c3 Fl = 2*c1*e(2) + c2 de(2) = - Fe / Fl e(2) = e(2) + de(2) if ( dabs(de(2)) <= tol_e ) exit end do iti2 = iti ! solução de e(3) u = k(3) q3 = f - e(1) + e(8) q4 = - e(1) - 2*e(2) + e(4) - e(5) + e(6) q5 = e(2) + e(4) - e(6) - 2*e(7) q6 = e(1) - e(4) - e(5) - e(6) - 2*e(8) c1 = u - 1.0d0 c2 = - (u*(q3+q5)+q4+q6) c3 = u*q3*q5 - q4*q6 do iti = 1, itimax Fe = c1*(e(3)**2) + c2*e(3) + c3 Fl = 2*c1*e(3) + c2 de(3) = - Fe / Fl e(3) = e(3) + de(3) if ( dabs(de(3)) <= tol_e ) exit end do iti3 = iti ! solução de e(4) u = k(4) q2 = d + e(7) q4 = - e(1) - 2*e(2) + e(3) - e(5) + e(6) q5 = e(2) - e(3) - e(6) - 2*e(7) q6 = e(1) + e(3) - e(5) - e(6) - 2*e(8) c1 = u - 1.0d0 c2 = -(u*(q2+q6)+q4+q5) c3 = u*q2*q6 - q4*q5 do iti = 1, itimax Fe = c1*(e(4)**2) + c2*e(4) + c3 Fl = 2*c1*e(4) + c2 de(4) = - Fe / Fl e(4) = e(4) + de(4) if ( dabs(de(4)) <= tol_e ) exit end do iti4 = iti if ( sai == 1 ) then write(out,1) iti1, iti2, iti3, iti4, iti5, iti6, iti7, iti8 1 format('iti(j) = ', 8(i6)) write(out,2) de 2 format('de(j) = ', 8(1pe22.9)) write(out,3) e 3 format('e(j) = ', 8(1pe22.9)) end if end subroutine GIBBS_EQUILIBRIO_Modelo_7 !------------------------------------------------------------------------------- subroutine GIBBS_EQUILIBRIO_Modelo_9 real*8 :: c1, c2, c3, q1, q2, q3, q4, q5, q6, q7, q8, u, Fe, Fl integer :: iti1, iti2, iti3, iti4, iti5, iti6, iti7, iti8, iti9, & iti10, iti11, iti12, iti13, iti14, iti15, iti16, iti17, iti18 ! solução de e(1) u = k(1)*p q1 = b + e( 7) + e(11) + e(14) + e(18) q4 = 2*e( 5) + 2*e( 6) - e( 7) + e( 8) + e( 9) & - e(11) + 2*e(12) + e(13) - 2*e(14) - e(18) q6 = 2*e( 2) - e( 4) + e( 7) - e( 8) & + e( 9) - e(10) - e(12) - e(15) - e(17) c1 = u c2 = - (u*(q4+q6)+n) c3 = u*q4*q6-q1*n do iti = 1, itimax Fe = c1*(e(1)**2) + c2*e(1) + c3 Fl = 2*c1*e(1) + c2 de(1) = - Fe / Fl e(1) = e(1) + de(1) if ( dabs(de(1)) <= tol_e ) exit end do iti1 = iti ! solução de e(2) u = k(2) q3 = f - e( 6) - e( 7) - e( 9) + e(15) + e(17) q6 = - e( 1) - e( 4) + e( 7) - e( 8) & + e( 9) - e(10) - e(12) - e(15) - e(17) c1 = 4*p c2 = 4*p*q6 + u*n c3 = p*((q6)**2) - u*q3*n do iti = 1, itimax Fe = c1*(e(2)**2) + c2*e(2) + c3 Fl = 2*c1*e(2) + c2 de(2) = - Fe / Fl e(2) = e(2) + de(2) if ( dabs(de(2)) <= tol_e ) exit end do iti2 = iti ! solução de e(3) u = k(3) q2 = d - e( 4) - e( 6) - e( 8) - e(10) & + e(11) + e(13) + e(15) + e(16) q5 = e( 8) - e( 9) - e(13) + e(14) c1 = 4*p c2 = 4*p*q5 + u*n c3 = p*((q5)**2) - u*q2*n do iti = 1, itimax Fe = c1*(e(3)**2) + c2*e(3) + c3 Fl = 2*c1*e(3) + c2 de(3) = - Fe / Fl e(3) = e(3) + de(3) if ( dabs(de(3)) <= tol_e ) exit end do iti3 = iti ! solução de e(4) u = k(4)*p q2 = d - e( 3) - e( 6) - e( 8) - e(10) & + e(11) + e(13) + e(15) + e(16) q6 = - e( 1) + 2*e( 2) + e( 7) - e( 8) & + e( 9) - e(10) - e(12) - e(15) - e(17) q7 = e(10) - e(11) - e(12) - e(13) & - e(15) - 2*e(16) + e(17) + e(18) c1 = u c2 = -(u*(q2+q6)+n) c3 = u*q2*q6 - q7*n do iti = 1, itimax Fe = c1*(e(4)**2) + c2*e(4) + c3 Fl = 2*c1*e(4) + c2 de(4) = - Fe / Fl e(4) = e(4) + de(4) if ( dabs(de(4)) <= tol_e ) exit end do iti4 = iti ! solução de e(5) u = k(5) q4 = - e( 1) + 2*e( 6) - e( 7) + e( 8) + e( 9) & - e(11) + 2*e(12) + e(13) - 2*e(14) - e(18) q8 = e(16) - e(17) - e(18) c1 = 4*p c2 = 4*p*q4 + u*n c3 = p*(q4**2) - u*q8*n do iti = 1, itimax Fe = c1*(e(5)**2) + c2*e(5) + c3 Fl = 2*c1*e(5) + c2 de(5) = - Fe / Fl e(5) = e(5) + de(5) if ( dabs(de(5)) <= tol_e ) exit end do iti5 = iti ! solução de e(6) u = k(6) q2 = d - e( 3) - e( 4) - e( 8) - e(10) & + e(11) + e(13) + e(15) + e(16) q3 = f - e( 2) - e( 7) - e( 9) + e(15) + e(17) q4 = - e( 1) + 2*e( 5) - e( 7) + e( 8) + e( 9) & - e(11) + 2*e(12) + e(13) - 2*e(14) - e(18) c1 = u - 4.0d0 c2 = -(u*(q2+q3)+4*q4) c3 = u*q2*q3 - (q4**2) do iti = 1, itimax Fe = c1*(e(6)**2) + c2*e(6) + c3 Fl = 2*c1*e(6) + c2 de(6) = - Fe / Fl e(6) = e(6) + de(6) if ( dabs(de(6)) <= tol_e ) exit end do iti6 = iti ! solução de e(7) u = k(7) q1 = b + e( 1) + e(11) + e(14) + e(18) q3 = f - e( 2) - e( 6) - e( 9) + e(15) + e(17) q4 = - e( 1) + 2*e( 5) + 2*e( 6) + e( 8) + e( 9) & - e(11) + 2*e(12) + e(13) - 2*e(14) - e(18) q6 = - e( 1) + 2*e( 2) - e( 4) - e( 8) & + e( 9) - e(10) - e(12) - e(15) - e(17) c1 = u - 1.0d0 c2 = -(u*(q3+q4)+q1+q6) c3 = u*q3*q4 - q1*q6 do iti = 1, itimax Fe = c1*(e(7)**2) + c2*e(7) + c3 Fl = 2*c1*e(7) + c2 de(7) = - Fe / Fl e(7) = e(7) + de(7) if ( dabs(de(7)) <= tol_e ) exit end do iti7 = iti ! solução de e(8) u = k(8) q2 = d - e( 3) - e( 4) - e( 6) - e(10) & + e(11) + e(13) + e(15) + e(16) q4 = - e( 1) + 2*e( 5) + 2*e( 6) - e( 7) + e( 9) & - e(11) + 2*e(12) + e(13) - 2*e(14) - e(18) q5 = 2*e( 3) - e( 9) - e(13) + e(14) q6 = - e( 1) + 2*e( 2) - e( 4) + e( 7) & + e( 9) - e(10) - e(12) - e(15) - e(17) c1 = u - 1.0d0 c2 = -(u*(q2+q6) + q4 + q5) c3 = u*q2*q6 - q4*q5 do iti = 1, itimax Fe = c1*(e(8)**2) + c2*e(8) + c3 Fl = 2*c1*e(8) + c2 de(8) = - Fe / Fl e(8) = e(8) + de(8) if ( dabs(de(8)) <= tol_e ) exit end do iti8 = iti ! solução de e(9) u = k(9) q3 = f - e( 2) - e( 6) - e( 7) + e(15) + e(17) q4 = - e( 1) + 2*e( 5) + 2*e( 6) - e( 7) + e( 8) & - e(11) + 2*e(12) + e(13) - 2*e(14) - e(18) q5 = 2*e( 3) + e( 8) - e(13) + e(14) q6 = - e( 1) + 2*e( 2) - e( 4) + e( 7) - e( 8) & - e(10) - e(12) - e(15) - e(17) c1 = u - 1.0d0 c2 = - (u*(q3+q5) + q4 + q6) c3 = u*q3*q5 - q4*q6 do iti = 1, itimax Fe = c1*(e(9)**2) + c2*e(9) + c3 Fl = 2*c1*e(9) + c2 de(9) = - Fe / Fl e(9) = e(9) + de(9) if ( dabs(de(9)) <= tol_e ) exit end do iti9 = iti ! solução de e(10) u = k(10)*p q2 = d - e( 3) - e( 4) - e( 6) - e( 8) & + e(11) + e(13) + e(15) + e(16) q6 = - e( 1) + 2*e( 2) - e( 4) + e( 7) - e( 8) & + e( 9) - e(12) - e(15) - e(17) q7 = e( 4) - e(11) - e(12) - e(13) & - e(15) - 2*e(16) + e(17) + e(18) c1 = u c2 = -(u*(q2+q6) + n) c3 = u*q2*q6 - q7*n do iti = 1, itimax Fe = c1*(e(10)**2) + c2*e(10) + c3 Fl = 2*c1*e(10) + c2 de(10) = - Fe / Fl e(10) = e(10) + de(10) if ( dabs(de(10)) <= tol_e ) exit end do iti10 = iti ! solução de e(11) u = k(11) q1 = b + e( 1) + e( 7) + e(14) + e(18) q2 = d - e( 3) - e( 4) - e( 6) - e( 8) - e(10) & + e(13) + e(15) + e(16) q4 = - e( 1) + 2*e( 5) + 2*e( 6) - e( 7) + e( 8) + e( 9) & + 2*e(12) + e(13) - 2*e(14) - e(18) q7 = e( 4) + e(10) - e(12) - e(13) & - e(15) - 2*e(16) + e(17) + e(18) c1 = u - 1.0d0 c2 = - (u*(q4+q7) + q1 + q2) c3 = u*q4*q7 - q1*q2 do iti = 1, itimax Fe = c1*(e(11)**2) + c2*e(11) + c3 Fl = 2*c1*e(11) + c2 de(11) = - Fe / Fl e(11) = e(11) + de(11) if ( dabs(de(11)) <= tol_e ) exit end do iti11 = iti ! solução de e(12) u = k(12) q4 = - e( 1) + 2*e( 5) + 2*e( 6) - e( 7) + e( 8) + e( 9) & - e(11) + e(13) - 2*e(14) - e(18) q6 = - e( 1) + 2*e( 2) - e( 4) + e( 7) - e( 8) & + e( 9) - e(10) - e(15) - e(17) q7 = e( 4) + e(10) - e(11) - e(13) & - e(15) - 2*e(16) + e(17) + e(18) c1 = u - 4.0d0 c2 = -(u*(q6+q7) + 4*q4) c3 = u*q6*q7 - q4*q4 do iti = 1, itimax Fe = c1*(e(12)**2) + c2*e(12) + c3 Fl = 2*c1*e(12) + c2 de(12) = - Fe / Fl e(12) = e(12) + de(12) if ( dabs(de(12)) <= tol_e ) exit end do iti12 = iti ! solução de e(13) u = k(13) q2 = d - e( 3) - e( 4) - e( 6) - e( 8) - e(10) & + e(11) + e(15) + e(16) q4 = - e( 1) + 2*e( 5) + 2*e( 6) - e( 7) + e( 8) + e( 9) & - e(11) + 2*e(12) - 2*e(14) - e(18) q5 = 2*e( 3) + e( 8) - e( 9) + e(14) q7 = e( 4) + e(10) - e(11) - e(12) & - e(15) - 2*e(16) + e(17) + e(18) c1 = u - 1.0d0 c2 = - (u*(q5+q7)+q2+q4) c3 = u*q5*q7-q2*q4 do iti = 1, itimax Fe = c1*(e(13)**2) + c2*e(13) + c3 Fl = 2*c1*e(13) + c2 de(13) = - Fe / Fl e(13) = e(13) + de(13) if ( dabs(de(13)) <= tol_e ) exit end do iti13 = iti ! solução de e(14) u = k(14) q1 = b + e( 1) + e( 7) + e(11) + e(18) q4 = - e( 1) + 2*e( 5) + 2*e( 6) - e( 7) + e( 8) + e( 9) & - e(11) + 2*e(12) + e(13) - e(18) q5 = 2*e( 3) + e( 8) - e( 9) - e(13) c1 = 4*u - 1.0d0 c2 = -(4*u*q4 + q1 + q5) c3 = u*q4*q4 - q1*q5 do iti = 1, itimax Fe = c1*(e(14)**2) + c2*e(14) + c3 Fl = 2*c1*e(14) + c2 de(14) = - Fe / Fl e(14) = e(14) + de(14) if ( dabs(de(14)) <= tol_e ) exit end do iti14 = iti ! solução de e(15) u = k(15) q2 = d - e( 3) - e( 4) - e( 6) - e( 8) - e(10) & + e(11) + e(13) + e(16) q3 = f - e( 2) - e( 6) - e( 7) - e( 9) + e(17) q6 = - e( 1) + 2*e( 2) - e( 4) + e( 7) - e( 8) & + e( 9) - e(10) - e(12) - e(17) q7 = e( 4) + e(10) - e(11) - e(12) - e(13) & - 2*e(16) + e(17) + e(18) c1 = u - 1.0d0 c2 = -(u*(q6+q7)+q2+q3) c3 = u*q6*q7 - q2*q3 do iti = 1, itimax Fe = c1*(e(15)**2) + c2*e(15) + c3 Fl = 2*c1*e(15) + c2 de(15) = - Fe / Fl e(15) = e(15) + de(15) if ( dabs(de(15)) <= tol_e ) exit end do iti15 = iti ! solução de e(16) u = k(16) q2 = d - e( 3) - e( 4) - e( 6) - e( 8) - e(10) & + e(11) + e(13) + e(15) q7 = e( 4) + e(10) - e(11) - e(12) - e(13) & - e(15) + e(17) + e(18) q8 = - e( 5) - e(17) - e(18) c1 = 4*u - 1.0d0 c2 = -(4*u*q7+q2+q8) c3 = u*((q7)**2) - q2*q8 do iti = 1, itimax Fe = c1*(e(16)**2) + c2*e(16) + c3 Fl = 2*c1*e(16) + c2 de(16) = - Fe / Fl e(16) = e(16) + de(16) if ( dabs(de(16)) <= tol_e ) exit end do iti16 = iti ! solução de e(17) u = k(17) q3 = f - e( 2) - e( 6) - e( 7) - e( 9) + e(15) q6 = - e( 1) + 2*e( 2) - e( 4) + e( 7) - e( 8) & + e( 9) - e(10) - e(12) - e(15) q7 = e( 4) + e(10) - e(11) - e(12) - e(13) & - e(15) - 2*e(16) + e(18) q8 = - e( 5) + e(16) - e(18) c1 = u - 1.0d0 c2 = -(u*(q6+q8)+q3+q7) c3 = u*q6*q8 - q3*q7 do iti = 1, itimax Fe = c1*(e(17)**2) + c2*e(17) + c3 Fl = 2*c1*e(17) + c2 de(17) = - Fe / Fl e(17) = e(17) + de(17) if ( dabs(de(17)) <= tol_e ) exit end do iti17 = iti ! solução de e(18) u = k(18) q1 = b + e( 1) + e( 7) + e(11) + e(14) q4 = - e( 1) + 2*e( 5) + 2*e( 6) - e( 7) + e( 8) + e( 9) & - e(11) + 2*e(12) + e(13) - 2*e(14) q7 = e( 4) + e(10) - e(11) - e(12) - e(13) & - e(15) - 2*e(16) + e(17) q8 = - e( 5) + e(16) - e(17) c1 = u - 1.0d0 c2 = -(u*(q4+q8) + q1 + q7) c3 = u*q4*q8 - q1*q7 do iti = 1, itimax Fe = c1*(e(18)**2) + c2*e(18) + c3 Fl = 2*c1*e(18) + c2 de(18) = - Fe / Fl e(18) = e(18) + de(18) if ( dabs(de(18)) <= tol_e ) exit end do iti18 = iti if ( sai == 1 ) then write(out,1) iti1, iti2, iti3, iti4, iti5, iti6, iti7, iti8, iti9, & iti10, iti11, iti12, iti13, iti14, iti15, iti16, iti17, iti18 1 format('iti(j) = ', 18(i6)) write(out,2) de 2 format('de(j) = ', 18(1pe22.9)) write(out,3) e 3 format('e(j) = ', 18(1pe22.9)) end if end subroutine GIBBS_EQUILIBRIO_Modelo_9 !------------------------------------------------------------------------------- subroutine GIBBS_EQUILIBRIO_Modelo_10 real*8 :: c1, c2, c3, q1, q2, q3, q4, q5, q6, q7, q8, u, Fe, Fl integer :: iti1, iti2, iti3, iti4, iti5, iti6 ! solução de e(1) u = k(1)*p q1 = b q4 = e(4) + 2*e(6) q6 = - 2*e(2) - e(4) - e(5) c1 = u c2 = - (u*(q4+q6)+n) c3 = u*q4*q6-q1*n do iti = 1, itimax Fe = c1*(e(1)**2) + c2*e(1) + c3 Fl = 2*c1*e(1) + c2 de(1) = - Fe / Fl e(1) = e(1) + de(1) if ( dabs(de(1)) <= tol_e ) exit end do iti1 = iti ! solução de e(2) u = k(2)*p q3 = f q6 = - e(1) - e(4) - e(5) c1 = 4*u c2 = -(4*u*q6+n) c3 = u*((q6)**2) - q3*n do iti = 1, itimax Fe = c1*(e(2)**2) + c2*e(2) + c3 Fl = 2*c1*e(2) + c2 de(2) = - Fe / Fl e(2) = e(2) + de(2) if ( dabs(de(2)) <= tol_e ) exit end do iti2 = iti ! solução de e(3) u = k(3)*p q2 = d - e(5) q5 = - e(4) c1 = 4*u c2 = - (4*u*q5+n) c3 = u*((q5)**2) - q2*n do iti = 1, itimax Fe = c1*(e(3)**2) + c2*e(3) + c3 Fl = 2*c1*e(3) + c2 de(3) = - Fe / Fl e(3) = e(3) + de(3) if ( dabs(de(3)) <= tol_e ) exit end do iti3 = iti ! solução de e(4) u = k(4)*p q4 = - e(1) + 2*e(6) q5 = - 2*e(3) q6 = - e(1) - 2*e(2) - e(5) c1 = u c2 = -(u*(q5+q6)+n) c3 = u*q5*q6 - q4*n do iti = 1, itimax Fe = c1*(e(4)**2) + c2*e(4) + c3 Fl = 2*c1*e(4) + c2 de(4) = - Fe / Fl e(4) = e(4) + de(4) if ( dabs(de(4)) <= tol_e ) exit end do iti4 = iti ! solução de e(5) u = k(5)*p q2 = d + e(3) q6 = - e(1) - 2*e(2) - e(4) q7 = 0.0d0 c1 = u c2 = - (u*(q2+q6)+n) c3 = u*q2*q6-q7*n do iti = 1, itimax Fe = c1*(e(5)**2) + c2*e(5) + c3 Fl = 2*c1*e(5) + c2 de(5) = - Fe / Fl e(5) = e(5) + de(5) if ( dabs(de(5)) <= tol_e ) exit end do iti5 = iti ! solução de e(6) u = k(6) q4 = - e(1) + e(4) q8 = 0.0d0 c1 = 4*p c2 = 4*q4*p + u*n c3 = (p*(q4**2)) - u*q8*n do iti = 1, itimax Fe = c1*(e(6)**2) + c2*e(6) + c3 Fl = 2*c1*e(6) + c2 de(6) = - Fe / Fl e(6) = e(6) + de(6) if ( dabs(de(6)) <= tol_e ) exit end do iti6 = iti if ( sai == 1 ) then write(out,1) iti1, iti2, iti3, iti4, iti5, iti6 1 format('iti(j) = ', 6(i6)) write(out,2) de 2 format('de(j) = ', 6(1pe22.9)) write(out,3) e 3 format('e(j) = ', 6(1pe22.9)) end if end subroutine GIBBS_EQUILIBRIO_Modelo_10 !------------------------------------------------------------------------------- subroutine GIBBS_EQUILIBRIO_Modelo_11 real*8 :: c1, c2, c3, q1, q2, q3, q4, q5, q6, u, Fe, Fl integer :: iti1, iti2, iti3 ! solução de e(1) u = k(1)*p q1 = b q4 = 0.0d0 q6 = 2*e( 2) c1 = u c2 = - (u*(q4+q6)+n) c3 = u*q4*q6-q1*n do iti = 1, itimax Fe = c1*(e(1)**2) + c2*e(1) + c3 Fl = 2*c1*e(1) + c2 de(1) = - Fe / Fl e(1) = e(1) + de(1) if ( dabs(de(1)) <= tol_e ) exit end do iti1 = iti ! solução de e(2) u = k(2) q3 = f q6 = - e( 1) c1 = 4*p c2 = 4*p*q6 + u*n c3 = p*((q6)**2) - u*q3*n do iti = 1, itimax Fe = c1*(e(2)**2) + c2*e(2) + c3 Fl = 2*c1*e(2) + c2 de(2) = - Fe / Fl e(2) = e(2) + de(2) if ( dabs(de(2)) <= tol_e ) exit end do iti2 = iti ! solução de e(3) u = k(3) q2 = d q5 = 0.0d0 c1 = 4*p c2 = 4*p*q5 + u*n c3 = p*((q5)**2) - u*q2*n do iti = 1, itimax Fe = c1*(e(3)**2) + c2*e(3) + c3 Fl = 2*c1*e(3) + c2 de(3) = - Fe / Fl e(3) = e(3) + de(3) if ( dabs(de(3)) <= tol_e ) exit end do iti3 = iti if ( sai == 1 ) then write(out,1) iti1, iti2, iti3 1 format('iti(j) = ', 3(i6)) write(out,2) de 2 format('de(j) = ', 3(1pe22.9)) write(out,3) e 3 format('e(j) = ', 3(1pe22.9)) end if end subroutine GIBBS_EQUILIBRIO_Modelo_11 !------------------------------------------------------------------------------- subroutine GIBBS_EQUILIBRIO_Modelo_12 real*8 :: c1, c2, c3, q1, q2, q3, q4, q5, q6, q7, q8, u, Fe, Fl integer :: iti1, iti2, iti3, iti4, iti5 ! solução de e(1) u = k(1)*p q1 = b q4 = 2*e( 5) q6 = 2*e( 2) - e( 4) c1 = u c2 = - (u*(q4+q6)+n) c3 = u*q4*q6-q1*n do iti = 1, itimax Fe = c1*(e(1)**2) + c2*e(1) + c3 Fl = 2*c1*e(1) + c2 de(1) = - Fe / Fl e(1) = e(1) + de(1) if ( dabs(de(1)) <= tol_e ) exit end do iti1 = iti ! solução de e(2) u = k(2) q3 = f q6 = - e( 1) - e( 4) c1 = 4*p c2 = 4*p*q6 + u*n c3 = p*((q6)**2) - u*q3*n do iti = 1, itimax Fe = c1*(e(2)**2) + c2*e(2) + c3 Fl = 2*c1*e(2) + c2 de(2) = - Fe / Fl e(2) = e(2) + de(2) if ( dabs(de(2)) <= tol_e ) exit end do iti2 = iti ! solução de e(3) u = k(3) q2 = d - e( 4) q5 = 0.0d0 c1 = 4*p c2 = 4*p*q5 + u*n c3 = p*((q5)**2) - u*q2*n do iti = 1, itimax Fe = c1*(e(3)**2) + c2*e(3) + c3 Fl = 2*c1*e(3) + c2 de(3) = - Fe / Fl e(3) = e(3) + de(3) if ( dabs(de(3)) <= tol_e ) exit end do iti3 = iti ! solução de e(4) u = k(4)*p q2 = d - e( 3) q6 = - e( 1) + 2*e( 2) q7 = 0.0d0 c1 = u c2 = -(u*(q2+q6)+n) c3 = u*q2*q6 - q7*n do iti = 1, itimax Fe = c1*(e(4)**2) + c2*e(4) + c3 Fl = 2*c1*e(4) + c2 de(4) = - Fe / Fl e(4) = e(4) + de(4) if ( dabs(de(4)) <= tol_e ) exit end do iti4 = iti ! solução de e(5) u = k(5) q4 = - e( 1) q8 = 0.0d0 c1 = 4*p c2 = 4*p*q4 + u*n c3 = p*(q4**2) - u*q8*n do iti = 1, itimax Fe = c1*(e(5)**2) + c2*e(5) + c3 Fl = 2*c1*e(5) + c2 de(5) = - Fe / Fl e(5) = e(5) + de(5) if ( dabs(de(5)) <= tol_e ) exit end do iti5 = iti if ( sai == 1 ) then write(out,1) iti1, iti2, iti3, iti4, iti5 1 format('iti(j) = ', 5(i6)) write(out,2) de 2 format('de(j) = ', 5(1pe22.9)) write(out,3) e 3 format('e(j) = ', 5(1pe22.9)) end if end subroutine GIBBS_EQUILIBRIO_Modelo_12 !------------------------------------------------------------------------------- ! ****************************************************************************** ! *** rotinas do bloco CÂMARA *** ! ****************************************************************************** ! Objetivo: calcular a temperatura de combustão e ! a composição de uma mistura gases !------------------------------------------------------------------------------- subroutine GIBBS_CAMARA_dados_recebe ( outt, sait, modelot, tipo_hrt, & itimaxt, itemaxt, itcmaxt, tol_et, & tol_nt, tol_Tct, T_Ot, T_Ft, & h_Ot, h_Ft, pt, OFt ) integer,intent(in) :: outt, sait, modelot, itimaxt, itemaxt, itcmaxt, & tipo_hrt real*8, intent(in) :: tol_et, tol_nt, tol_Tct, T_Ot, T_Ft, pt, & OFt, h_Ot, h_Ft out = outt sai = sait modelo = modelot tipo_hr = tipo_hrt itimax = itimaxt tol_e = tol_et itemax = itemaxt tol_n = tol_nt itcmax = itcmaxt tol_Tc = tol_Tct T_O = T_Ot T_F = T_Ft h_O = h_Ot h_F = h_Ft p = pt OF = OFt end subroutine GIBBS_CAMARA_dados_recebe !------------------------------------------------------------------------------- subroutine GIBBS_CAMARA_dados_escreve ( out ) integer,intent(in) :: out ! número do arquivo de saída write(out,1) modelo, Ne, Nr, p, OF, equiv, a, Ru, itimax, tol_e, itemax, & tol_n, itcmax, tol_Tc, h_F, h_O, a, b, d, f, no 1 format(/, t1,'modelo =', i4, & /, t1,'Número de espécies químicas =', i4, & /, t1,'Número de reações =', i4, & /, t1,'pressão (atm) =', 1pe22.9, & /, t1,'razão O/F mássico (adim.) =', 1pe22.9, & /, t1,'razão de equivalência (adim.) =', 1pe22.9, & /, t1,'razão O/F molar (adim.) =', 1pe22.9, & /, t1,'Constante universal dos gases (J/mol.K) = ', 1pe22.9, & /, t1,'itimax =', i4, & /, t1,'tol_e =', 1pe22.9, & /, t1,'itemax =', i4, & /, t1,'tol_n =', 1pe22.9, & /, t1,'itcmax =', i4, & /, t1,'tol_Tc =', 1pe22.9, & /, t1,'entalpia de injeção do combustível na câmara (J/mol) =', 1pe22.9, & /, t1,'entalpia de injeção do oxidante na câmara (J/mol) =', 1pe22.9, & /, t1,'número de moles do reagente O2 (mol) =', 1pe22.9, & /, t1,'número de moles do reagente H2 = 1 mol', & /, t1,'número de moles do produto H2O sem dissociação (mol) =', 1pe22.9, & /, t1,'número de moles do produto O2 sem dissociação (mol) =', 1pe22.9, & /, t1,'número de moles do produto H2 sem dissociação (mol) =', 1pe22.9, & /, t1,'número total inicial de moles dos produtos (mol) =', 1pe22.9 ) end subroutine GIBBS_CAMARA_dados_escreve !------------------------------------------------------------------------------- subroutine GIBBS_CAMARA_inicio call GIBBS_EQUILIBRIO_inicio1 allocate ( hi(Ne) ) if ( tipo_hr == 1 ) then T = T_F call GIBBS_hi_molar_calculo h_F = hi(3) T = T_O call GIBBS_hi_molar_calculo h_O = hi(2) end if HR = h_F + a * h_O end subroutine GIBBS_CAMARA_inicio !------------------------------------------------------------------------------- subroutine GIBBS_CAMARA_calculos1 real*8 :: T1, T2, Ta, dH if ( sai == 1 ) write(out,1) 1 format(/, t11,'itc', t21,'dTc', t43,'Tc') T1 = 3.0d+2 T2 = 6.0d+3 T = (T1 + T2) / 2 Ta = T1 do itc = 1, itcmax call GIBBS_EQUILIBRIO_calculos1 call GIBBS_hi_molar_calculo HP = 0.0d0 do i = 1, Ne HP = HP + ni(i)*hi(i) end do dH = HP - HR if ( dH < 0.0d0 ) then T1 = T else T2 = T end if dTc = T - Ta if ( sai == 1 ) write(out,2) itc, dTc, T 2 format( t11,i6, t21,1pe22.9, t43,1pe22.9, ' *** Tc *** ' ) if ( dabs(dTc) <= tol_Tc ) exit Ta = T T = (T1 + T2) / 2 end do Tc = T end subroutine GIBBS_CAMARA_calculos1 !------------------------------------------------------------------------------- subroutine GIBBS_CAMARA_escreve ( out ) integer,intent(in) :: out ! número do arquivo de saída write(out,1) Tc 1 format(/, t1,'Temperatura de combustão (K) =', 1pe22.9) end subroutine GIBBS_CAMARA_escreve !------------------------------------------------------------------------------- ! ****************************************************************************** ! *** rotinas do bloco TAXA FINITA *** ! ****************************************************************************** ! Objetivo: calcular a taxa de geração de massa por unidade de volume de ! cada espécie química !------------------------------------------------------------------------------- subroutine GIBBS_TAXA_FINITA_dados_recebe ( outt, sait, modelot, Tt, pt ) integer,intent(in) :: outt, sait, modelot real*8, intent(in) :: Tt, pt out = outt sai = sait modelo = modelot T = Tt p = pt end subroutine GIBBS_TAXA_FINITA_dados_recebe !------------------------------------------------------------------------------- subroutine GIBBS_TAXA_FINITA_inicio select case ( modelo ) case ( 5, 7 ) Ne = 6 Nr = 8 case ( 31, 32 ) Ne = 6 Nr = 4 case ( 10 ) Ne = 8 Nr = 6 case ( 9 ) Ne = 8 Nr =18 end select call GIBBS_le_coeficientes_cp_h_g call GIBBS_TAXA_FINITA_aloca_memoria dos = system('notepad Gibbs_1p3_Taxa_Finita2.in') ! lista dados open(in,file='Gibbs_1p3_Taxa_Finita2.in') do i = 1, Ne read(in,*) Yi(i) end do end subroutine GIBBS_TAXA_FINITA_inicio !------------------------------------------------------------------------------- subroutine GIBBS_TAXA_FINITA_dados_escreve ( out ) integer,intent(in) :: out ! número do arquivo de saída write(out,1) modelo, Ne, Nr, T, p, Ru 1 format(/, t1,'modelo =', i4, & /, t1,'Número de espécies químicas =', i4, & /, t1,'Número de reações =', i4, & /, t1,'Temperatura (K) =', 1pe22.9, & /, t1,'pressão (atm) =', 1pe22.9, & /, t1,'Constante universal dos gases (J/mol.K) = ', 1pe22.9 ) call GIBBS_Yi_escreve (out ) end subroutine GIBBS_TAXA_FINITA_dados_escreve !------------------------------------------------------------------------------- subroutine GIBBS_TAXA_FINITA_aloca_memoria allocate ( Yi(Ne), gi(Ne), Wi(Ne), Ci(Ne) ) allocate ( kf(Nr), kb(Nr), dg(Nr), k(Nr), gam(Nr), lam(Nr), teta(Nr) ) end subroutine GIBBS_TAXA_FINITA_aloca_memoria !------------------------------------------------------------------------------- subroutine GIBBS_TAXA_FINITA_fecha_memoria deallocate ( Yi, gi, Wi, Ci ) deallocate ( kf, kb, dg, k, gam, lam, teta ) end subroutine GIBBS_TAXA_FINITA_fecha_memoria !------------------------------------------------------------------------------- subroutine GIBBS_TAXA_FINITA_calculos1 p_si = p * po R = 0.0d0 do i = 1, Ne R = R + Yi(i) / Moli(i) end do R = 1.0d+3 * R * Ru ro_si = p_si / ( R * T ) ro = 1.0d-3 * ro_si call GIBBS_Ci_calculo call GIBBS_C_calculo end subroutine GIBBS_TAXA_FINITA_calculos1 !------------------------------------------------------------------------------- subroutine GIBBS_TAXA_FINITA_escreve1 write(out,1) po, p_si, ro, ro_si, R, C 1 format(/, t1,'pressão ao nível do mar (Pa) =', 1pe22.9, & /, t1,'pressão (Pa) =', 1pe22.9, & /, t1,'massa específica (g/cm3) =', 1pe22.9, & /, t1,'massa específica (kg/m3) =', 1pe22.9, & /, t1,'constante da mistura de gases (J/kg.K) = ', 1pe22.9, & /, t1,'concentração (mol/cm3) =', 1pe22.9 ) call GIBBS_Ci_escreve ( out ) end subroutine GIBBS_TAXA_FINITA_escreve1 !------------------------------------------------------------------------------- subroutine GIBBS_TAXA_FINITA_calculos2 ! cálculo da energia livre de Gibbs de cada espécie (molar) call GIBBS_gi_molar_calculo ! cálculo dos Wi(i) de cada modelo select case ( modelo ) case ( 5 ) call GIBBS_TAXA_FINITA_Modelo_5 case ( 7 ) call GIBBS_TAXA_FINITA_Modelo_7 case ( 31 ) call GIBBS_TAXA_FINITA_Modelo_31 case ( 32 ) call GIBBS_TAXA_FINITA_Modelo_32 case ( 10 ) call GIBBS_TAXA_FINITA_Modelo_10 case ( 9 ) call GIBBS_TAXA_FINITA_Modelo_9 end select end subroutine GIBBS_TAXA_FINITA_calculos2 !------------------------------------------------------------------------------- subroutine GIBBS_TAXA_FINITA_escreve2 write(out,1) 1 format(/, t1,'reação', t18,'delta Gibbs', t40,'kp', t62,'kf', t84,'kb') do j = 1, Nr write(out,2) j, dg(j), k(j), kf(j), kb(j) end do 2 format( t1,i6, t11, 4(1pe22.9) ) write(out,3) 3 format(/, t1,'reação', t18,'gama', t40,'lambda', t62,'teta') do j = 1, Nr write(out,2) j, gam(j), lam(j), teta(j) end do write(out,4) 4 format(/, t1,'espécie', t16,'Wi (kg/m3.s)' ) do i = 1, Ne write(out,5) especie(i), Wi(i) 5 format( t3,a4, t10,1pe22.9 ) end do end subroutine GIBBS_TAXA_FINITA_escreve2 !------------------------------------------------------------------------------- subroutine GIBBS_k_calculo ! cálculo das constantes de equilíbrio baseadas nas pressões parciais do j = 1, Nr k(j) = dexp ( - dg(j) / ( Ru * T ) ) end do end subroutine GIBBS_k_calculo !------------------------------------------------------------------------------- subroutine GIBBS_teta_calculo ! cálculo dos teta(j) do j = 1, Nr teta(j) = lam(j) * gam(j) end do end subroutine GIBBS_teta_calculo !------------------------------------------------------------------------------- subroutine GIBBS_TAXA_FINITA_Modelo_5 ! cálculo dos dg(j) de cada modelo dg(1) = gi(1) - gi(4) - gi(6) dg(2) = gi(3) - 2 * gi(6) dg(3) = gi(2) - 2 * gi(5) dg(4) = gi(4) - gi(5) - gi(6) dg(5) = gi(2) + gi(6) - gi(4) - gi(5) dg(6) = gi(3) + gi(5) - gi(4) - gi(6) dg(7) = gi(1) + gi(6) - gi(3) - gi(4) dg(8) = gi(1) + gi(5) - 2 * gi(4) ! cálculo das constantes de equilíbrio baseadas nas pressões parciais call GIBBS_k_calculo ! cálculo dos kf(j) kf(1) = 3.626D19 / T kf(2) = 7.500D17 / T kf(3) = 3.626D18 / T kf(4) = 3.626D18 / T kf(5) = 4.500D14 * dexp (- 30.d0/T) / ( T ** 0.5d0 ) kf(6) = 4.900D03 * dexp (-1950.d0/T) * ( T ** 2.8d0 ) kf(7) = 6.300D06 * dexp (-1490.d0/T) * ( T ** 2.0d0 ) kf(8) = 2.100D08 * dexp ( 200.d0/T) * ( T ** 1.4d0 ) ! cálculo dos kb(j) do j = 1, 4 kb(j) = kf(j) * ( (Ru*T)**(-1) ) / k(j) end do do j = 5, 8 kb(j) = kf(j) / k(j) end do ! cálculo dos gam(j) gam(1) = kf(1)*Ci(4)*Ci(6) - kb(1)*Ci(1) gam(2) = kf(2)*Ci(6)*Ci(6) - kb(2)*Ci(3) gam(3) = kf(3)*Ci(5)*Ci(5) - kb(3)*Ci(2) gam(4) = kf(4)*Ci(5)*Ci(6) - kb(4)*Ci(4) gam(5) = kf(5)*Ci(4)*Ci(5) - kb(5)*Ci(2)*Ci(6) gam(6) = kf(6)*Ci(4)*Ci(6) - kb(6)*Ci(3)*Ci(5) gam(7) = kf(7)*Ci(3)*Ci(4) - kb(7)*Ci(1)*Ci(6) gam(8) = kf(8)*Ci(4)*Ci(4) - kb(8)*Ci(1)*Ci(5) ! cálculo dos lam(j) do j = 1, 4 lam(j) = C end do do j = 5, 8 lam(j) = 1.0d0 end do ! cálculo dos teta(j) call GIBBS_teta_calculo ! cálculo dos Wi(i) Wi(1) = Moli(1) * (teta(1) + teta(7) + teta(8)) Wi(2) = Moli(2) * (teta(3) + teta(5)) Wi(3) = Moli(3) * (teta(2) + teta(6) - teta(7)) Wi(4) = - Moli(4) * (teta(1) - teta(4) + teta(5) + teta(6) + teta(7) + 2*teta(8)) Wi(5) = Moli(5) * (teta(6) + teta(8) - teta(4) - teta(5) - 2*teta(3)) Wi(6) = - Moli(6) * (teta(1) + teta(4) - teta(5) + teta(6) - teta(7) + 2*teta(2)) Wi = 1.0d+3 * Wi end subroutine GIBBS_TAXA_FINITA_Modelo_5 !------------------------------------------------------------------------------- subroutine GIBBS_TAXA_FINITA_Modelo_7 ! cálculo dos dg(j) de cada modelo dg(1) = gi(1) + gi(6) - gi(3) - gi(4) dg(2) = gi(1) + gi(5) - 2 * gi(4) dg(3) = gi(4) + gi(6) - gi(3) - gi(5) dg(4) = gi(4) + gi(5) - gi(2) - gi(6) dg(5) = gi(1) - gi(4) - gi(6) dg(6) = gi(4) - gi(5) - gi(6) dg(7) = gi(2) - 2 * gi(5) dg(8) = gi(3) - 2 * gi(6) ! cálculo das constantes de equilíbrio baseadas nas pressões parciais call GIBBS_k_calculo ! cálculo dos kf(j) kf(1) = 2.20D13 * dexp (- 5150.0d0/(1.986d0*T)) kf(2) = 6.30D12 * dexp (- 1000.0d0/(1.986d0*T)) kf(3) = 1.80D10 * dexp (- 8900.0d0/(1.986d0*T)) * T kf(4) = 2.20D14 * dexp (-16800.0d0/(1.986d0*T)) kf(5) = 8.40D21 / (T**2) kf(6) = 3.62D18 / T kf(7) = 1.90D13 * dexp ( 1790.0d0/(1.986d0*T)) kf(8) = 6.40D17 / T ! cálculo dos kb(j) do j = 1, 4 kb(j) = kf(j) / k(j) end do do j = 5, 8 kb(j) = kf(j) * ( (Ru*T)**(-1) ) / k(j) end do ! cálculo dos gam(j) gam(1) = kf(1)*Ci(3)*Ci(4) - kb(1)*Ci(1)*Ci(6) gam(2) = kf(2)*Ci(4)*Ci(4) - kb(2)*Ci(1)*Ci(5) gam(3) = kf(3)*Ci(3)*Ci(5) - kb(3)*Ci(4)*Ci(6) gam(4) = kf(4)*Ci(2)*Ci(6) - kb(4)*Ci(4)*Ci(5) gam(5) = kf(5)*Ci(4)*Ci(6) - kb(5)*Ci(1) gam(6) = kf(6)*Ci(5)*Ci(6) - kb(6)*Ci(4) gam(7) = kf(7)*Ci(5)*Ci(5) - kb(7)*Ci(2) gam(8) = kf(8)*Ci(6)*Ci(6) - kb(8)*Ci(3) ! cálculo dos lam(j) do j = 1, 4 lam(j) = 1.0d0 end do lam(5) = 17.d0*Ci(1) + 6.0d0*Ci(2) + 5.d0*Ci(3) + 12.5d0*(Ci(4) + Ci(5) + Ci(6)) lam(6) = 5.d0*Ci(1) + 5.0d0*Ci(2) + 5.d0*Ci(3) + 12.5d0*(Ci(4) + Ci(5) + Ci(6)) lam(7) = 5.d0*Ci(1) + 11.0d0*Ci(2) + 5.d0*Ci(3) + 12.5d0*(Ci(4) + Ci(5) + Ci(6)) lam(8) = 10.d0*Ci(1) + 1.5d0*Ci(2) + 4.d0*Ci(3) + 25.0d0*(Ci(4) + Ci(5) + Ci(6)) ! cálculo dos teta(j) call GIBBS_teta_calculo ! cálculo dos Wi(i) Wi(1) = Moli(1) * (teta(1) + teta(2) + teta(5)) Wi(2) = Moli(2) * (teta(7) - teta(4)) Wi(3) = Moli(3) * (teta(8) - teta(1) - teta(3)) Wi(4) = Moli(4) * (teta(3) + teta(4) + teta(6) - teta(1) - teta(5) - 2*teta(2)) Wi(5) = Moli(5) * (teta(2) + teta(4) - teta(3) - teta(6) - 2*teta(7)) Wi(6) = Moli(6) * (teta(1) + teta(3) - teta(4) - teta(5) - teta(6) - 2*teta(8)) Wi = 1.0d+3 * Wi end subroutine GIBBS_TAXA_FINITA_Modelo_7 !------------------------------------------------------------------------------- subroutine GIBBS_TAXA_FINITA_Modelo_31 ! cálculo dos dg(j) de cada modelo dg(1) = gi(1) - gi(4) - gi(6) dg(2) = gi(3) - 2 * gi(6) dg(3) = gi(2) - 2 * gi(5) dg(4) = gi(4) - gi(5) - gi(6) ! cálculo das constantes de equilíbrio baseadas nas pressões parciais call GIBBS_k_calculo ! cálculo dos kf(j) kf(1) = 3.626D19 / T kf(2) = 7.500D17 / T kf(3) = 3.626D18 / T kf(4) = 3.626D18 / T ! cálculo dos kb(j) do j = 1, Nr kb(j) = kf(j) * ( (Ru*T)**(-1) ) / k(j) end do ! cálculo dos gam(j) gam(1) = kf(1)*Ci(4)*Ci(6) - kb(1)*Ci(1) gam(2) = kf(2)*Ci(6)*Ci(6) - kb(2)*Ci(3) gam(3) = kf(3)*Ci(5)*Ci(5) - kb(3)*Ci(2) gam(4) = kf(4)*Ci(5)*Ci(6) - kb(4)*Ci(4) ! cálculo dos lam(j) lam = C ! cálculo dos teta(j) call GIBBS_teta_calculo ! cálculo dos Wi(i) Wi(1) = Moli(1) * teta(1) Wi(2) = Moli(2) * teta(3) Wi(3) = Moli(3) * teta(2) Wi(4) = Moli(4) * (teta(4) - teta(1)) Wi(5) = - Moli(5) * (2*teta(3) + teta(4)) Wi(6) = - Moli(6) * (teta(1) + 2*teta(2) + teta(4)) Wi = 1.0d+3 * Wi end subroutine GIBBS_TAXA_FINITA_Modelo_31 !------------------------------------------------------------------------------- subroutine GIBBS_TAXA_FINITA_Modelo_32 ! cálculo dos dg(j) de cada modelo dg(1) = gi(1) - gi(4) - gi(6) dg(2) = gi(3) - 2 * gi(6) dg(3) = gi(2) - 2 * gi(5) dg(4) = gi(4) - gi(5) - gi(6) ! cálculo das constantes de equilíbrio baseadas nas pressões parciais call GIBBS_k_calculo ! cálculo dos kf(j) kf(1) = 8.4D21 / (T**2) kf(2) = 6.4D17 / T kf(3) = 1.9D13 * dexp (1790.0d0/(1.986d0*T)) kf(4) = 3.62D18 / T ! cálculo dos kb(j) do j = 1, Nr kb(j) = kf(j) * ( (Ru*T)**(-1) ) / k(j) end do ! cálculo dos gam(j) gam(1) = kf(1)*Ci(4)*Ci(6) - kb(1)*Ci(1) gam(2) = kf(2)*Ci(6)*Ci(6) - kb(2)*Ci(3) gam(3) = kf(3)*Ci(5)*Ci(5) - kb(3)*Ci(2) gam(4) = kf(4)*Ci(5)*Ci(6) - kb(4)*Ci(4) ! cálculo dos lam(j) lam(1) = 17.d0*Ci(1) + 6.0d0*Ci(2) + 5.d0*Ci(3) + 12.5d0*(Ci(4) + Ci(5) + Ci(6)) lam(2) = 10.d0*Ci(1) + 1.5d0*Ci(2) + 4.d0*Ci(3) + 25.0d0*(Ci(4) + Ci(5) + Ci(6)) lam(3) = 5.d0*Ci(1) + 11.0d0*Ci(2) + 5.d0*Ci(3) + 12.5d0*(Ci(4) + Ci(5) + Ci(6)) lam(4) = 5.d0*Ci(1) + 5.0d0*Ci(2) + 5.d0*Ci(3) + 12.5d0*(Ci(4) + Ci(5) + Ci(6)) ! cálculo dos teta(j) call GIBBS_teta_calculo ! cálculo dos Wi(i) Wi(1) = Moli(1) * teta(1) Wi(2) = Moli(2) * teta(3) Wi(3) = Moli(3) * teta(2) Wi(4) = Moli(4) * (teta(4) - teta(1)) Wi(5) = - Moli(5) * (2*teta(3) + teta(4)) Wi(6) = - Moli(6) * (teta(1) + 2*teta(2) + teta(4)) Wi = 1.0d+3 * Wi end subroutine GIBBS_TAXA_FINITA_Modelo_32 !------------------------------------------------------------------------------- subroutine GIBBS_TAXA_FINITA_Modelo_9 ! cálculo dos dg(j) de cada modelo dg(1) = gi(1) - gi(4) - gi(6) dg(2) = 2 * gi(6) - gi(3) dg(3) = 2 * gi(5) - gi(2) dg(4) = gi(7) - gi(2) - gi(6) dg(5) = 2 * gi(4) - gi(8) dg(6) = 2 * gi(4) - gi(2) - gi(3) dg(7) = gi(1) + gi(6) - gi(3) - gi(4) dg(8) = gi(4) + gi(5) - gi(2) - gi(6) dg(9) = gi(4) + gi(6) - gi(3) - gi(5) dg(10)= gi(2) + gi(7) - 2 * gi(2) - gi(6) dg(11)= gi(1) + gi(2) - gi(4) - gi(7) dg(12)= 2 * gi(4) - gi(6) - gi(7) dg(13)= gi(2) + gi(4) - gi(5) - gi(7) dg(14)= gi(1) + gi(5) - 2 * gi(4) dg(15)= gi(2) + gi(3) - gi(6) - gi(7) dg(16)= gi(2) + gi(8) - 2 * gi(7) dg(17)= gi(3) + gi(7) - gi(6) - gi(8) dg(18)= gi(1) + gi(7) - gi(4) - gi(8) ! cálculo das constantes de equilíbrio baseadas nas pressões parciais call GIBBS_k_calculo ! cálculo dos kf(j) kf( 1) = 7.500D23 / ( T ** 2.6d0 ) kf( 2) = 2.230D12 * ( T ** 0.5d0 ) * dexp ( -92600.d0 / (1.986d0 * T) ) kf( 3) = 1.850D11 * ( T ** 0.5d0 ) * dexp ( -95560.d0 / (1.986d0 * T) ) kf( 4) = 2.100D18 / T kf( 5) = 1.300D17 * dexp ( -45500.d0 / (1.986d0 * T) ) kf( 6) = 1.700D13 * dexp ( -47780.d0 / (1.986d0 * T) ) kf( 7) = 1.170D09 * ( T ** 1.3d0 ) * dexp ( - 3626.d0 / (1.986d0 * T) ) kf( 8) = 5.130D16 * dexp ( -16507.d0 / (1.986d0 * T) ) / ( T ** 0.816d0 ) kf( 9) = 1.800D10 * T * dexp ( - 8826.d0 / (1.986d0 * T) ) kf(10) = 6.700D19 / ( T ** 1.42d0 ) kf(11) = 5.000D13 * dexp ( - 1000.d0 / (1.986d0 * T) ) kf(12) = 2.500D14 * dexp ( - 1900.d0 / (1.986d0 * T) ) kf(13) = 4.800D13 * dexp ( - 1000.d0 / (1.986d0 * T) ) kf(14) = 6.000D08 * ( T ** 1.3d0 ) kf(15) = 2.500D13 * dexp ( - 700.d0 / (1.986d0 * T) ) kf(16) = 2.000D12 kf(17) = 1.600D12 * dexp ( - 3800.d0 / (1.986d0 * T) ) kf(18) = 1.000D13 * dexp ( - 1800.d0 / (1.986d0 * T) ) ! cálculo dos kb(j) j = 1 kb(j) = kf(j) * ( (Ru*T)**(-1) ) / k(j) do j = 2, 3 kb(j) = kf(j) * Ru*T / k(j) end do j = 4 kb(j) = kf(j) * ( (Ru*T)**(-1) ) / k(j) j = 5 kb(j) = kf(j) * Ru*T / k(j) do j = 6, 9 kb(j) = kf(j) / k(j) end do j = 10 kb(j) = kf(j) * ( (Ru*T)**(-1) ) / k(j) do j = 11, 18 kb(j) = kf(j) / k(j) end do ! cálculo dos gam(j) gam( 1) = kf( 1)*Ci(4)*Ci(6) - kb( 1)*Ci(1) gam( 2) = kf( 2)*Ci(3) - kb( 2)*Ci(6)*Ci(6) gam( 3) = kf( 3)*Ci(2) - kb( 3)*Ci(5)*Ci(5) gam( 4) = kf( 4)*Ci(2)*Ci(6) - kb( 4)*Ci(7) gam( 5) = kf( 5)*Ci(8) - kb( 5)*Ci(4)*Ci(4) gam( 6) = kf( 6)*Ci(2)*Ci(3) - kb( 6)*Ci(4)*Ci(4) gam( 7) = kf( 7)*Ci(3)*Ci(4) - kb( 7)*Ci(1)*Ci(6) gam( 8) = kf( 8)*Ci(2)*Ci(6) - kb( 8)*Ci(4)*Ci(5) gam( 9) = kf( 9)*Ci(3)*Ci(5) - kb( 9)*Ci(4)*Ci(6) gam(10) = kf(10)*Ci(2)*Ci(6) - kb(10)*Ci(7) gam(11) = kf(11)*Ci(4)*Ci(7) - kb(11)*Ci(1)*Ci(2) gam(12) = kf(12)*Ci(6)*Ci(7) - kb(12)*Ci(4)*Ci(4) gam(13) = kf(13)*Ci(5)*Ci(7) - kb(13)*Ci(2)*Ci(4) gam(14) = kf(14)*Ci(4)*Ci(4) - kb(14)*Ci(1)*Ci(5) gam(15) = kf(15)*Ci(6)*Ci(7) - kb(15)*Ci(2)*Ci(3) gam(16) = kf(16)*Ci(7)*Ci(7) - kb(16)*Ci(2)*Ci(8) gam(17) = kf(17)*Ci(6)*Ci(8) - kb(17)*Ci(3)*Ci(7) gam(18) = kf(18)*Ci(4)*Ci(8) - kb(18)*Ci(1)*Ci(7) ! cálculo dos lam(j) lam(1) = 20.d0*Ci(1) + Ci(2) + Ci(3) + Ci(4) + Ci(5) + Ci(6) + Ci(7) + Ci(8) lam(2) = 6.d0*Ci(1) + 3.d0*Ci(3) + 2.d0*Ci(6) + Ci(2) + Ci(4) + Ci(5) + Ci(7) + Ci(8) lam(3) = C lam(4) = 21.d0*Ci(1) + 3.3d0*Ci(3) + Ci(4) + Ci(5) + Ci(6) + Ci(7) + Ci(8) lam(5) = C do j = 6, 18 lam(j) = 1.0d0 end do ! cálculo dos teta(j) call GIBBS_teta_calculo ! cálculo dos Wi(i) Wi(1) = Moli(1) * (teta( 1) + teta( 7) + teta(11) + teta(14) + teta(18)) Wi(2) = Moli(2) * (teta(11) + teta(13) + teta(15) + teta(16) - teta( 3) & - teta( 4) - teta( 6) - teta( 8) - teta(10) ) Wi(3) = Moli(3) * (teta(15) + teta(17) - teta( 2) - teta( 6) - teta( 7) & - teta( 9) ) Wi(4) = Moli(4) * (2*teta( 5) + 2*teta(6) + teta(8) + teta( 9) + 2*teta(12) & + teta(13) - teta(1) - teta(7) - teta(11) - 2*teta(14) & - teta(18) ) Wi(5) = Moli(5) * (2*teta( 3) + teta( 8) + teta(14) - teta( 9) - teta(13)) Wi(6) = Moli(6) * (2*teta( 2) + teta( 7) + teta( 9) - teta( 1) - teta( 4) & - teta( 8) - teta(10) - teta(12) - teta(15) - teta(17)) Wi(7) = Moli(7) * ( teta( 4) + teta(10) + teta(17) + teta(18) - teta(11) & - teta(12) - teta(13) - teta(15) - 2*teta(16)) Wi(8) = Moli(8) * ( teta(16) - teta( 5) - teta(17) - teta(18)) Wi = 1.0d+3 * Wi end subroutine GIBBS_TAXA_FINITA_Modelo_9 !------------------------------------------------------------------------------- subroutine GIBBS_TAXA_FINITA_Modelo_10 ! cálculo dos dg(j) de cada modelo dg(1) = gi(1) - gi(4) - gi(6) dg(2) = gi(3) - 2 * gi(6) dg(3) = gi(2) - 2 * gi(5) dg(4) = gi(4) - gi(5) - gi(6) dg(5) = gi(7) - gi(2) - gi(6) dg(6) = 2 * gi(4) - gi(8) ! cálculo das constantes de equilíbrio baseadas nas pressões parciais call GIBBS_k_calculo ! cálculo dos kf(j) kf(1) = 3.626D19 / T kf(2) = 7.500D17 / T kf(3) = 3.626D18 / T kf(4) = 3.626D18 / T kf(5) = 2.100D18 / T kf(6) = 1.300D17 * dexp ( -45500.d0 / (1.986d0 * T) ) ! cálculo dos kb(j) do j = 1, 5 kb(j) = kf(j) * ( (Ru*T)**(-1) ) / k(j) end do j = 6 kb(j) = kf(j) * ( Ru*T ) / k(j) ! cálculo dos gam(j) gam(1) = kf(1)*Ci(4)*Ci(6) - kb(1)*Ci(1) gam(2) = kf(2)*Ci(6)*Ci(6) - kb(2)*Ci(3) gam(3) = kf(3)*Ci(5)*Ci(5) - kb(3)*Ci(2) gam(4) = kf(4)*Ci(5)*Ci(6) - kb(4)*Ci(4) gam(5) = kf(5)*Ci(2)*Ci(6) - kb(5)*Ci(7) gam(6) = kf(6)*Ci(8) - kb(6)*Ci(4)*Ci(4) ! cálculo dos lam(j) lam = C ! cálculo dos teta(j) call GIBBS_teta_calculo ! cálculo dos Wi(i) Wi(1) = Moli(1) * teta(1) Wi(2) = Moli(2) * (teta(3) - teta(5)) Wi(3) = Moli(3) * teta(2) Wi(4) = Moli(4) * (teta(4) - teta(1) + 2*teta(6)) Wi(5) = - Moli(5) * (teta(4) + 2*teta(3)) Wi(6) = - Moli(6) * (teta(1) + teta(4) + teta(5) + 2*teta(2)) Wi(7) = Moli(7) * teta(5) Wi(8) = - Moli(8) * teta(6) Wi = 1.0d+3 * Wi end subroutine GIBBS_TAXA_FINITA_Modelo_10 !------------------------------------------------------------------------------- ! ****************************************************************************** ! *** rotinas do bloco MACH1D 5.0 *** ! ****************************************************************************** ! Objetivo: rotinas específicas para uso com o programa MACH1D 5.0 !------------------------------------------------------------------------------- subroutine GIBBS_MACH1D_envia1 ( cpt, gamat ) real*8,intent(out) :: cpt, gamat cpt = cp gamat = gama end subroutine GIBBS_MACH1D_envia1 !------------------------------------------------------------------------------- subroutine GIBBS_MACH1D_envia2 ( Net ) integer,intent(out) :: Net ! número de espécies químicas Net = Ne end subroutine GIBBS_MACH1D_envia2 !------------------------------------------------------------------------------- subroutine GIBBS_MACH1D_envia3 ( especiet ) character*4,dimension(Ne),intent(out) :: especiet ! nome da espécie química especiet = especie end subroutine GIBBS_MACH1D_envia3 !------------------------------------------------------------------------------- subroutine GIBBS_MACH1D_envia4 ( cpt, gamat, Rt, Yit, hit ) real*8, intent(out) :: cpt, gamat, Rt real*8,dimension(Ne),intent(out) :: Yit, hit cpt = cp gamat = gama Rt = R Yit = Yi hit = hi end subroutine GIBBS_MACH1D_envia4 !------------------------------------------------------------------------------- subroutine GIBBS_MACH1D_envia5 ( Rt ) real*8,intent(out) :: Rt Rt = R end subroutine GIBBS_MACH1D_envia5 !------------------------------------------------------------------------------- subroutine GIBBS_MACH1D_envia6 ( cpt, gamat, Rt, Yit ) real*8, intent(out) :: cpt, gamat, Rt real*8,dimension(Ne),intent(out) :: Yit cpt = cp gamat = gama Rt = R Yit = Yi end subroutine GIBBS_MACH1D_envia6 !------------------------------------------------------------------------------- subroutine GIBBS_MACH1D_envia7 (ej) real*8,dimension(Nr),intent(out) :: ej ej = e end subroutine GIBBS_MACH1D_envia7 !------------------------------------------------------------------------------- subroutine GIBBS_MACH1D_envia8 (Nrt) integer,intent(out) :: Nrt Nrt = Nr end subroutine GIBBS_MACH1D_envia8 !------------------------------------------------------------------------------- subroutine GIBBS_MACH1D_recebe1 ( Yit ) real*8,dimension(Ne),intent(in) :: Yit Yi = Yit end subroutine GIBBS_MACH1D_recebe1 !------------------------------------------------------------------------------- subroutine GIBBS_MACH1D_recebe2 ( Tt ) real*8,intent(in) :: Tt T = Tt end subroutine GIBBS_MACH1D_recebe2 !------------------------------------------------------------------------------- subroutine GIBBS_MACH1D_recebe3 (ivol) integer,intent(in) :: ivol vol = ivol end subroutine GIBBS_MACH1D_recebe3 !------------------------------------------------------------------------------- subroutine GIBBS_MACH1D_recebe4 (itt) integer,intent(in) :: itt itgeral = itt end subroutine GIBBS_MACH1D_recebe4 !------------------------------------------------------------------------------- subroutine GIBBS_MACH1D_recebe5 (n) integer,intent(in) :: n nvol = n end subroutine GIBBS_MACH1D_recebe5 !------------------------------------------------------------------------------- subroutine GIBBS_MACH1D_recebe6 (ej_a) real*8,dimension(:,:),intent(in) :: ej_a ej = ej_a end subroutine GIBBS_MACH1D_recebe6 !------------------------------------------------------------------------------- subroutine GIBBS_MACH1D_recebe7 (nitm) integer,intent(in) :: nitm max_it = nitm end subroutine GIBBS_MACH1D_recebe7 !------------------------------------------------------------------------------- subroutine GIBBS_MACH1D_alocacao_ej allocate ( ej(nvol,Nr) ) end subroutine GIBBS_MACH1D_alocacao_ej !------------------------------------------------------------------------------- subroutine GIBBS_MACH1D_desalocacao_ej deallocate ( ej ) end subroutine GIBBS_MACH1D_desalocacao_ej !------------------------------------------------------------------------------- subroutine GIBBS_MACH1D_compara_ej integer :: cont ! contador auxiliar (reação) do cont = 1, Nr ej(vol,cont) = e(cont) end do end subroutine GIBBS_MACH1D_compara_ej !------------------------------------------------------------------------------- subroutine GIBBS_MACH1D_compara_e integer :: cont ! contador auxiliar (reação) do cont = 1, Nr e(cont) = ej(vol,cont) end do end subroutine GIBBS_MACH1D_compara_e !------------------------------------------------------------------------------- subroutine GIBBS_MACH1D_abre_ej_a open(arq_ej_a, file ='ej_a.dat') end subroutine GIBBS_MACH1D_abre_ej_a !------------------------------------------------------------------------------- subroutine GIBBS_MACH1D_fecha_ej_a close(arq_ej_a) end subroutine GIBBS_MACH1D_fecha_ej_a !------------------------------------------------------------------------------- subroutine GIBBS_MACH1D_renomeia_ej integer :: cont cont = rename("ej.dat","ej_a.dat") end subroutine GIBBS_MACH1D_renomeia_ej !------------------------------------------------------------------------------- subroutine GIBBS_MACH1D_leitura_ej_a integer :: cont ! contador auxiliar (reação) integer :: cont1 ! contador auxiliar (volume) integer :: cont2 ! auxiliar do cont1 = 1, vol read(arq_ej_a,*) cont2, cont2, (e(cont), cont = 1, Nr) end do end subroutine GIBBS_MACH1D_leitura_ej_a !------------------------------------------------------------------------------- subroutine GIBBS_MACH1D_escreve_ej integer :: cont ! contador auxliar (reação) write(arq_ej,5) itgeral, vol, (e(cont), cont = 1, Nr) 5 format ( <2>(i8), (1pe25.15) ) end subroutine GIBBS_MACH1D_escreve_ej !------------------------------------------------------------------------------- subroutine GIBBS_MACH1D_abre_ej open(arq_ej, file ='ej.dat') end subroutine GIBBS_MACH1D_abre_ej !------------------------------------------------------------------------------- subroutine GIBBS_MACH1D_fecha_ej close(arq_ej) end subroutine GIBBS_MACH1D_fecha_ej !------------------------------------------------------------------------------- subroutine GIBBS_MACH1D_inicio1 allocate ( Yi1(Ne), Yi2(Ne), hi(Ne) ) end subroutine GIBBS_MACH1D_inicio1 !------------------------------------------------------------------------------- subroutine GIBBS_MACH1D_reativo1 Yi1 = Yi end subroutine GIBBS_MACH1D_reativo1 !------------------------------------------------------------------------------- subroutine GIBBS_MACH1D_reativo2 Yi2 = Yi end subroutine GIBBS_MACH1D_reativo2 !------------------------------------------------------------------------------- subroutine GIBBS_MACH1D_reativo3 ( Tt, dTY, cprt ) real*8,intent(in) :: Tt, dTY real*8,intent(out) :: cprt T = Tt call GIBBS_hi_massa_calculo cpr = 0.0d0 do i = 1, Ne cpr = cpr + hi(i) * ( Yi1(i) - Yi2(i) ) end do cpr = cpr / dTY cprt = cpr end subroutine GIBBS_MACH1D_reativo3 !------------------------------------------------------------------------------- subroutine GIBBS_MACH1D_reativo4 ( dpY, gamaf, Molt, pPa, gamaet ) real*8,intent(in) :: dpY, gamaf, Molt, pPa real*8,intent(out) :: gamaet real*8 :: Cg Cg = 0.0d0 do i = 1, Ne Cg = Cg + Moli(i) * ( Yi1(i) - Yi2(i) ) end do Cg = 1.0d0 + Cg * pPa / (dpY * Molt) gamae = gamaf / Cg gamaet = gamae end subroutine GIBBS_MACH1D_reativo4 !------------------------------------------------------------------------------- subroutine GIBBS_MACH1D_reativo5 call GIBBS_Xi_calculo Yi1 = Xi end subroutine GIBBS_MACH1D_reativo5 !------------------------------------------------------------------------------- subroutine GIBBS_MACH1D_reativo6 call GIBBS_Xi_calculo Yi2 = Xi end subroutine GIBBS_MACH1D_reativo6 !------------------------------------------------------------------------------- subroutine GIBBS_Xi_aloca_memoria allocate ( Xi(Ne) ) end subroutine GIBBS_Xi_aloca_memoria !------------------------------------------------------------------------------- subroutine GIBBS_si_aloca_memoria allocate ( si(Ne) ) end subroutine GIBBS_si_aloca_memoria !------------------------------------------------------------------------------- subroutine GIBBS_si_molar_calculo ! si (J/mol.K) if ( T >= 1.0d+3 ) then m = 1 else m = 2 end if do i = 1, Ne si(i) = ( ca(i,1,m)*dlog(T) + ca(i,2,m)*T + ca(i,3,m)*(T**2)/2.0d0 & + ca(i,4,m)*(T**3)/3.0d0 + ca(i,5,m)*(T**4)/4.0d0 + ca(i,7,m) ) * Ru end do end subroutine GIBBS_si_molar_calculo !------------------------------------------------------------------------------- subroutine GIBBS_MACH1D_camara_entropia ( st ) real*8,intent(out) :: st real*8 :: s1, s2, s3 call GIBBS_Xi_calculo call GIBBS_si_molar_calculo s1 = 0.0d0 s2 = 0.0d0 do i = 1, Ne s1 = s1 + Xi(i)*si(i) if ( Xi(i) > 0.0d0 ) s2 = s2 + Xi(i)*dlog(Xi(i)) end do st = 1.0d+3 * ( s1 - Ru*dlog(p) - Ru*s2 ) / Mol end subroutine GIBBS_MACH1D_camara_entropia !------------------------------------------------------------------------------- subroutine GIBBS_MACH1D_entropia ( Tt, pt, Rt, Yit, st ) real*8,intent(in) :: Tt, pt, Rt real*8,dimension(Ne),intent(in) :: Yit real*8,intent(out) :: st real*8 :: s1, s2 T = Tt p = pt R = Rt Yi = Yit Mol = 1.0d+3 * Ru / R Xi = Mol * Yi / Moli call GIBBS_si_molar_calculo s1 = 0.0d0 s2 = 0.0d0 do i = 1, Ne s1 = s1 + Xi(i)*si(i) if ( Xi(i) > 0.0d0 ) s2 = s2 + Xi(i)*dlog(Xi(i)) end do st = 1.0d+3 * ( s1 - Ru*dlog(p) - Ru*s2 ) / Mol end subroutine GIBBS_MACH1D_entropia !------------------------------------------------------------------------------- subroutine GIBBS_MACH1D_temperatura ( Tt, Yit, h_est, cp_est ) real*8,intent(in) :: Tt real*8,dimension(Ne),intent(in) :: Yit real*8,intent(out) :: h_est, cp_est T = Tt Yi = Yit call GIBBS_cp_calculo cp_est = cp call GIBBS_hi_massa_calculo h_est = 0.0d0 do i = 1, Ne h_est = h_est + Yi(i) * hi(i) end do end subroutine GIBBS_MACH1D_temperatura !------------------------------------------------------------------------------- subroutine GIBBS_MACH1D_finito_1 ( modelot, Tt, pt, Yit, Wit, cpt, gamat, Rt, hit ) integer,intent(in) :: modelot real*8, intent(in) :: Tt, pt real*8,dimension(Ne),intent(in) :: Yit real*8, intent(out) :: cpt, gamat, Rt real*8,dimension(Ne),intent(out) :: Wit, hit modelo = modelot T = Tt p = pt Yi = Yit call GIBBS_TAXA_FINITA_calculos1 call GIBBS_TAXA_FINITA_calculos2 Wit = Wi call GIBBS_cp_calculo call GIBBS_gama_calculo call GIBBS_hi_massa_calculo cpt = cp gamat = gama Rt = R hit = hi end subroutine GIBBS_MACH1D_finito_1 !------------------------------------------------------------------------------- subroutine GIBBS_MACH1D_finito_2 allocate ( Wi(Ne), Ci(Ne) ) allocate ( kf(Nr), kb(Nr), gam(Nr), lam(Nr), teta(Nr) ) end subroutine GIBBS_MACH1D_finito_2 !------------------------------------------------------------------------------- end module GIBBS_1p4