module graos use variaveis implicit none !-------------------------------------------------------------------------------------------------------- contains !-------------------------------------------------------------------------------------------------------- subroutine cilindrico ! Área de queima Ab = pi*ri**2 ! Cálculo do volume da câmara de combustão Vol = Ab*x !+Vol0 if(x>=L) Ab = 0 end subroutine !-------------------------------------------------------------------------------------------------------- subroutine tubular ! Área de queima Ab = 2*pi*L*(Ri + x) ! Cálculo do volume da câmara de combustão Vol0 = pi*ri**2*L ! Cálculo do volume da câmara de combustão Vol = Vol0 + pi*L*( ((Ri+x)**2) - (Ri**2) ) if(x>=Re-Ri) Ab = 0 end subroutine !-------------------------------------------------------------------------------------------------------- subroutine desinibido ! Área de queima ai = 2 * pi * (ri+x) * (l-2*x) ae = 2 * pi * (re-x) * (l-2*x) au = 2 * pi * ( ((re-x)**2) - ((ri+x)**2) ) ab = ai + ae + au ! Cálculo do volume da câmara de combustão Vol0 = pi*ri**2*L voli= pi * l * ( (ri+x)**2 ) vole= pi * l * ( (re**2) - ((re-x)**2) ) volu= pi * 2 * x * ( ((re-x)**2) - ((ri+x)**2) ) vol = vol0 + Voli + Vole + Volu if(x>=(Re-Ri)/2) Ab = 0 end subroutine !-------------------------------------------------------------------------------------------------------- subroutine cil_desinib ! Ãrea das coroas Au = 2 * pi * (Re-x) ** 2 ! Ãrea externa Ae = 2 * pi * (Re-x) * (L-2*x) Ab = Ae + Au ! Volume das coroas Volu= pi* 2 * x * ((Re-x)**2) ! Volume externo Vole= pi * ( Re**2 - (Re-x)**2)*L Vol = Vole + Volu !+Vol0 if(x>=Re) Ab = 0 if(x>=L) Ab = 0 end subroutine !-------------------------------------------------------------------------------------------------------- subroutine tub_2cdesinib ! Ãrea interna Ai = 2 * pi * (Ri+x) * (L-2*x) ! Ãrea das coroas Au = 2 * pi * ( ((Re)**2) - ((Ri+x)**2) ) ! Ãrea de queima total Ab = Ai + Au ! Volume interno Voli= pi * (L-2*x) * ( (Ri+x)**2 - Ri**2 ) ! Volume das coroas Volu= pi * 2 * x * ( ((Ri+x)**2) - (Ri**2)) Vol = Voli + Volu if(x>=Re-Ri) Ab = 0 if(x>=L) Ab = 0 end subroutine !-------------------------------------------------------------------------------------------------------- subroutine tub_1cdesinib ! Ãrea interna Ai = 2 * pi * (Ri+x) * (L-x) ! Ãrea das coroas Au = pi * ( ((Re)**2) - ((Ri+x)**2) ) ! Ãrea de queima total Ab = Ai + Au ! Volume interno Voli= pi * (L-x) * ( (Ri+x)**2 - Ri**2 ) ! Volume das coroas Volu= pi * x * ( ((Ri+x)**2) - (Ri**2)) Vol = Voli + Volu if(x>=Re-Ri) Ab = 0 if(x>=L) Ab = 0 end subroutine !-------------------------------------------------------------------------------------------------------- subroutine tc_inib Rme = Rme + r*dt Rma = Rma + r*dt if(Rma <= Re) then ! Ãrea interna Ab = Pi*(Rma+Rme)*sqrt(L**2+(Rme-Rma)**2) ! Volume interno Vol= pi * L/3 * ( Rma**2 + Rma*Rme + Rme**2 ) elseif(Rme <= Re) then Lx = 2*L/(Dma-Dme)*(Re-Rme) ! Ãrea interna Ab = Pi*(Re+Rme)*sqrt(Lx**2+(Rme-Re)**2) ! Parte constante Vol = pi*Re**2*(L-Lx) ! Volume interno Vol= Vol + pi * Lx/3 * ( Re**2 + Re*Rme + Rme**2 ) else Ab = 0 Vol = PI*Re**2*L endif end subroutine !-------------------------------------------------------------------------------------------------------- subroutine tc_desinib if(Rma <= Re) then ! Parcela, queima coroas Rme = Rme + r*dt *(Dma-Dme)/2/L Rma = Rma - r*dt *(Dma-Dme)/2/L ! Ãrea de queima das coroas Ab = Pi*(Re**2 - Rme ** 2) Ab = Ab + Pi*(Re**2 - Rma ** 2) ! Volume de queima das coroas Vol = Pi*(Re**2-(Rme-r*dt*(Dma-Dme)/2/L)**2) * r * dt - pi * r*dt/3 * & ( Rme**2 + Rme*(Rme-r*dt*(Dma-Dme)/2/L) + (Rme-r*dt*(Dma-Dme)/2/L)**2 ) Vol = Vol + Pi*(Re**2-(Rma+r*dt*(Dma-Dme)/2/L)**2) * r * dt - pi * r*dt/3 * & ( Rma**2 + Rma*(Rma+r*dt*(Dma-Dme)/2/L) + (Rma+r*dt*(Dma-Dme)/2/L)**2 ) ! Queima interna do tronco de cone Rme = Rme + r*dt Rma = Rma + r*dt ! Ãrea interna Ab = Ab + Pi*(Rma+Rme)*sqrt((L-2*x)**2+(Rme-Rma)**2) ! Volume interno Vol= Vol + pi * (L-2*x)/3 * ( Rma**2 + Rma*Rme + Rme**2 ) ! Parte constante Vol= Vol + pi*Re**2*(2*x) elseif(Rme <= Re) then if(Lnovo == -1) then Lnovo = 0 Lx = 2*L/(Dma-Dme)*(Re-Rme) endif ! Parcela, queima coroa Rme = Rme + r*dt *(Dma-Dme)/2/L ! Ãrea de queima da coroa Ab = Pi*(Re**2 - Rme ** 2) ! Volume de queima da coroa Vol = Pi*(Re**2-(Rme-r*dt*(Dma-Dme)/2/L)**2) * r * dt - pi * r*dt/3 * & ( Rme**2 + Rme*(Rme-r*dt*(Dma-Dme)/2/L) + (Rme-r*dt*(Dma-Dme)/2/L)**2 ) ! Queima interna do tronco de cone Rme = Rme + r*dt Lx = Lx - r*dt ! Ãrea interna Ab = Ab + Pi*(Re+Rme)*sqrt(Lx**2+(Rme-Re)**2) ! Volume interno Vol= Vol + pi * Lx/3 * ( Re**2 + Re*Rme + Rme**2 ) Lx = Lx - r*dt*Lx/Rme ! Parte constante Vol= Vol + pi*Re**2*(L-Lx) else Ab = 0 Vol = PI*Re**2*L endif end subroutine !-------------------------------------------------------------------------------------------------------- subroutine tc_1cinib_cig if(Dma > 0 .and. Dme > 0) then if(Rma <= Re) then ! Fase 1: Tronco de cone com coroa do lado de diâmetro maior desinibida + parte cigarro ! QUEIMA DA COROA ! Parcela, queima coroa Rma = Rma - r*dt *(Dma-Dme)/2/L ! Ãrea de queima da coroa Ab = Pi*(Re**2 - Rma ** 2) ! Volume de queima das coroas Vol = Pi*(Re**2-(Rma+r*dt*(Dma-Dme)/2/L)**2) * r * dt - & pi * r*dt/3 * ( Rma**2 + Rma*(Rma+r*dt*(Dma-Dme)/2/L) + (Rma+r*dt*(Dma-Dme)/2/L)**2 ) ! QUEIMA DO TRONCO DE CONE ! Queima interna do tronco de cone Rme = Rme + r*dt Rma = Rma + r*dt ! Ãrea interna Ab = Ab + Pi*(Rma+Rme)*sqrt((L-x)**2+(Rme-Rma)**2) ! Volume interno Vol= Vol + pi * (L-x)/3 * ( Rma**2 + Rma*Rme + Rme**2 ) ! Parte constante Vol= Vol + pi*Re**2*(x) ! QUEIMA DA PARTE CIGARRO if(Lcig >= x) then ! Queima da parcela cigarro Ab = Ab + Pi*(Rme)**2 ! Queima da parcela tubular Ab = Ab + 2*Pi*(Rme) * x ! Volume interno Vol= Vol + pi*(Rme)**2*x elseif(Lcig < x) then ! Queima da parcela tubular Ab = Ab + 2*Pi*(Rme) * Lcig ! Volume interno Vol= Vol + pi*(Rme)**2* Lcig endif ! write(*,*) Vol, Ab ; pause elseif(Rme <= Re) then ! Fase 2: R Chegou na parede, sem queima da coroa do lado de maior diâmetro ! QUEIMA DO TRONCO DE CONE ! Queima interna do tronco de cone Rme = Rme + r*dt Lx = 2*L/(Dma-Dme)*(Re-Rme) ! Ãrea interna Ab = Pi*(Re+Rme)*sqrt(Lx**2+(Rme-Re)**2) ! Parte constante Vol = pi*Re**2*(L-Lx) ! Volume interno Vol= Vol + pi * (Lx)/3 * ( Re**2 + Re*Rme + Rme**2 ) ! QUEIMA DA PARTE CIGARRO if(Lcig >= x) then ! Queima da parcela cigarro Ab = Ab + Pi*(Rme)**2 ! Queima da parcela tubular Ab = Ab + 2*Pi*(Rme) * x ! Volume interno Vol= Vol + pi*(Rme)**2*x elseif(Lcig < x) then ! Queima da parcela tubular Ab = Ab + 2*Pi*(Rme) * Lcig ! Volume interno Vol= Vol + pi*(Rme)**2* Lcig endif ! write(*,'(I3, 20(1PE12.3))') i, Rme, Rma, Re, L, Lx, Ab, Vol elseif(Lcig>=x) then ! write(*,'(I3, 20(1PE12.3))') i, Rme, Rma, Re, L, Lx, Ab, Vol ; pause ! Fase 3: r chegou na parede, queima apenas cigarro ! Queima da parcela cigarro Ab = Pi*(Re)**2 ! Parte constante Vol = pi*Re**2*(L+x) else Ab = 0 Vol = pi*Re**2*(L+Lcig) endif else ! Queima Apenas CIGARRO Ab = Pi*(Re)**2 ! Volume Desocupado Vol = pi*Re**2*(x) if(x>=Lcig) Ab = 0 endif end subroutine !-------------------------------------------------------------------------------------------------------- subroutine taxa_inibido !a = a * (1+vmax*qeros) !r = a * ( p**n ) ! variavel para analisar quandos dgrao já chegaram no dext do motor indice_grao1 = indice_grao do j=1, ativacao ! Se for além do máximo permitido (Diâmetro externo do grão) colocar o valor de Re no Rativ if( Rativ(j) > Re) then indice_grao = j Rativ(j) = Re endif enddo ! Calculo das áreas de queima e do volume atual Ab = 0.d0 ! Toda a superfície foi ativada if(i > ativacao) then do j=ativacao, indice_grao+2, -1 Ab = Ab + 0.5d0*(Rativ(j-1)+Rativ(j)) enddo ! Parte da superfície foi ativada else ! Se estiver no inicio da queima j=2 if(i==1) then Ab = Ab + 0.5d0*(Rativ(j-1)+Rativ(j)) endif do j=i, indice_grao+2, -1 Ab = Ab + 0.5d0*(Rativ(j-1)+Rativ(j)) enddo endif if(indice_grao - indice_grao1 > 0) Ab = Ab + (indice_grao - indice_grao1)*Re Vol= 0.d0 do j=2, ativacao Vol= Vol+ ( 0.5d0*(Rativ(j-1)+Rativ(j)) )**2 enddo Ab = Ab*2*pi*dgrao Vol = Vol*pi*dgrao ! Taxa de queima consome o propelente if(i <= ativacao) then j=1 Rativ(j) = Rativ(j) + dt*r do j=2, i Rativ(j) = Rativ(j) + dt*r enddo else do j=indice_grao+1, ativacao Rativ(j) = Rativ(j) + dt*r enddo endif if(freq > 0) then if(mod(i,freq) == 0) then passo = passo + 1 if(passo<10) then write(texto,'("0000",I1)') passo elseif(passo<100) then write(texto,'("000",I2)') passo elseif(passo<1000) then write(texto,'("00",I3)') passo elseif(passo<10000) then write(texto,'("0",I4)') passo elseif(passo<100000) then write(texto,'(I5)') passo endif OPEN(2,FILE='./temp/'// trim(adjustl(texto)) //'.txt') DO J=1, ativacao WRITE(2,*) dgrao*(j-1), Rativ(J) ENDDO WRITE(2,*) dgrao*(j-1), Rativ(J-1) close(2) endif endif end subroutine !-------------------------------------------------------------------------------------------------------- subroutine taxa_desinibido2 !a = a * (1+vmax*qeros) !r = a * ( p**n ) ! variavel para analisar quandos dgrao já chegaram no dext do motor indice_grao1 = indice_grao !write(*,*) indice_grao1,indice_grao2 do j=1, indice_grao2 !ativacao ! Se for além do máximo permitido (Diâmetro externo do grão) colocar o valor de Re no Rativ if( Rativ(j) >= Re) then indice_grao = j Rativ(j) = Re endif enddo ! Calculo das áreas de queima e do volume atual Ab = 0.d0 ! Toda a superfície foi ativada if(i > indice_grao2) then do j=indice_grao+1, indice_grao2 Ab = Ab + 0.5d0*(Rativ(j-1)+Rativ(j)) enddo ! Parte da superfície ainda não foi ativada else ! Se estiver no inicio da queima if(i<=2) Ab = Ab + 0.25d0*Rativ(1) do j=indice_grao+2, i Ab = Ab + 0.5d0*(Rativ(j-1)+Rativ(j)) enddo endif !write(*,*) indice_grao,indice_grao2,Ab !if(indice_grao - indice_grao1 > 0) Ab = Ab + (indice_grao - indice_grao1)*Re Vol= 0.d0 do j=2, ativacao Vol= Vol+ ( ( Rativ(j-1)+Rativ(j) )/2 )**2 enddo Ab = Ab*2*pi*dgrao Vol = Vol*pi*dgrao ! Taxa de queima consome o propelente no raio if(i <= ativacao) then ! Caso onde o grão não foi totalmente ativado do j=1, i Rativ(j) = Rativ(j) + dt*r enddo else ! Caso onde o grão já foi totalmente ativado, se indice_grão>0 chegou no diâmetro externo do j=indice_grao+1, indice_grao2 Rativ(j) = Rativ(j) + dt*r enddo endif !write(*,*) Ab, indice_grao;pause ! Taxa de queima consome o grão no comprimento if(i <= ativacao) then ! Caso onde o grão não foi totalmente ativado j=1 Rativ(j) = Rativ(j) + dt*r do j=2, i Rativ(j) = Rativ(j) + dt*r enddo ! Consome grão do lado da ativação do j=indice_grao+1, i if(x>(j-1)*dgrao)then ! Parte já foi consumida if(x >= j*dgrao) then Rativ(j) = Re else ! Parte está sendo consumida Ab = Ab + pi * ( Re**2 - Rativ(j)**2 ) endif endif enddo else if(xgrao == 0) xgrao = x ! Caso onde o grão já foi totalmente ativado, se indice_grão>0 chegou no diâmetro externo do j=indice_grao+1, indice_grao2 Rativ(j) = Rativ(j) + dt*r enddo do j=indice_grao+1, indice_grao2 ! Consome grão do lado da ativação if(x>(j-1)*dgrao)then ! Parte já foi consumida if(x >= j*dgrao) then Rativ(j) = Re else ! Parte está sendo consumida Ab = Ab + pi * ( Re**2 - Rativ(j)**2 ) endif elseif(x-xgrao>=L-(j-1)*dgrao)then ! Consome grão do lado oposto da ativação ! Parte já foi consumida if(x-xgrao>L-(j)*dgrao) then Rativ(j) = Re indice_grao2= j-1 else ! Parte está sendo consumida Ab = Ab + pi * ( Re**2 - Rativ(j)**2 ) endif endif enddo endif ! Criação do GIF if(freq > 0) then if(mod(i,freq) == 0) then passo = passo + 1 if(passo<10) then write(texto,'("0000",I1)') passo elseif(passo<100) then write(texto,'("000",I2)') passo elseif(passo<1000) then write(texto,'("00",I3)') passo elseif(passo<10000) then write(texto,'("0",I4)') passo elseif(passo<100000) then write(texto,'(I5)') passo endif OPEN(2,FILE='./temp/'// trim(adjustl(texto)) //'.txt') DO J=1, ativacao WRITE(2,*) dgrao*(j-1), Rativ(J) ENDDO WRITE(2,*) dgrao*(j-1), Rativ(J-1) close(2) endif endif end subroutine !-------------------------------------------------------------------------------------------------------- subroutine taxa_desinibido1 !a = a * (1+vmax*qeros) !r = a * ( p**n ) ! variavel para analisar quandos dgrao já chegaram no dext do motor indice_grao1 = indice_grao !write(*,*) indice_grao1,indice_grao2 do j=1, indice_grao2 !ativacao ! Se for além do máximo permitido (Diâmetro externo do grão) colocar o valor de Re no Rativ if( Rativ(j) >= Re) then indice_grao = j Rativ(j) = Re endif enddo ! Calculo das áreas de queima e do volume atual Ab = 0.d0 ! Toda a superfície foi ativada if(i > indice_grao2) then do j=indice_grao+1, indice_grao2 Ab = Ab + 0.5d0*(Rativ(j-1)+Rativ(j)) enddo ! Parte da superfície ainda não foi ativada else ! Se estiver no inicio da queima if(i<=2) Ab = Ab + 0.25d0*Rativ(1) do j=indice_grao+2, i Ab = Ab + 0.5d0*(Rativ(j-1)+Rativ(j)) enddo endif !write(*,*) indice_grao,indice_grao2,Ab !if(indice_grao - indice_grao1 > 0) Ab = Ab + (indice_grao - indice_grao1)*Re Vol= 0.d0 do j=2, ativacao Vol= Vol+ ( ( Rativ(j-1)+Rativ(j) )/2 )**2 enddo Ab = Ab*2*pi*dgrao Vol = Vol*pi*dgrao ! Taxa de queima consome o propelente no raio if(i <= ativacao) then ! Caso onde o grão não foi totalmente ativado do j=1, i Rativ(j) = Rativ(j) + dt*r enddo else ! Caso onde o grão já foi totalmente ativado, se indice_grão>0 chegou no diâmetro externo do j=indice_grao+1, indice_grao2 Rativ(j) = Rativ(j) + dt*r enddo endif !write(*,*) Ab, indice_grao;pause ! Taxa de queima consome o grão no comprimento if(i <= ativacao) then ! Caso onde o grão não foi totalmente ativado j=1 Rativ(j) = Rativ(j) + dt*r do j=2, i Rativ(j) = Rativ(j) + dt*r enddo ! Consome grão do lado da ativação do j=indice_grao+1, i if(x>(j-1)*dgrao)then ! Parte já foi consumida if(x >= j*dgrao) then Rativ(j) = Re else ! Parte está sendo consumida Ab = Ab + pi * ( Re**2 - Rativ(j)**2 ) endif endif enddo else if(xgrao == 0) xgrao = x ! Caso onde o grão já foi totalmente ativado, se indice_grão>0 chegou no diâmetro externo do j=indice_grao+1, indice_grao2 Rativ(j) = Rativ(j) + dt*r enddo do j=indice_grao+1, indice_grao2 ! Consome grão do lado da ativação if(x>(j-1)*dgrao)then ! Parte já foi consumida if(x >= j*dgrao) then Rativ(j) = Re else ! Parte está sendo consumida Ab = Ab + pi * ( Re**2 - Rativ(j)**2 ) endif !elseif(x-xgrao>=L-(j-1)*dgrao)then !! Consome grão do lado oposto da ativação ! ! Parte já foi consumida ! if(x-xgrao>L-(j)*dgrao) then ! Rativ(j) = Re ! indice_grao2= j-1 ! else ! ! Parte está sendo consumida ! Ab = Ab + pi * ( Re**2 - Rativ(j)**2 ) ! endif ! endif enddo endif ! Criação do GIF if(freq > 0) then if(mod(i,freq) == 0) then passo = passo + 1 if(passo<10) then write(texto,'("0000",I1)') passo elseif(passo<100) then write(texto,'("000",I2)') passo elseif(passo<1000) then write(texto,'("00",I3)') passo elseif(passo<10000) then write(texto,'("0",I4)') passo elseif(passo<100000) then write(texto,'(I5)') passo endif OPEN(2,FILE='./temp/'// trim(adjustl(texto)) //'.txt') DO J=1, ativacao WRITE(2,*) dgrao*(j-1), Rativ(J) ENDDO WRITE(2,*) dgrao*(j-1), Rativ(J-1) close(2) endif endif end subroutine !-------------------------------------------------------------------------------------------------------- end module graos