Program empuxo1p1


! Diego
! versão 1.2.1:  8 Dez 14
! Corrigido erro no cálculo do volume da câmara de combustão
! Arquivos de saída são colocados na pasta output, via título que é uma variável de entrada
! do programa

! Diego
! versão 1.2:  5 Dez 14
! Adição do tipo de grão 4 (grão de ativação constante) com animação em gif 
! , melhor interação com o usuário (dados via arquivo de texto)
! , adicionado caso onde a pressão gerada na câmara é inferior a pressão ambiente
! assim, a pressão se torna a pressão ambiente, não gerando portanto força alguma

! Marchi
! versão 1.1:  4 Jun 05
! versão 1.0: 30 Mai 05
! última alteração: 1 Nov 2013

USE PORTLIB

implicit none

real*8  :: z1, z2, z3
real*8  :: Vol, Rg, Tc, Ab, At, a, n, rop, gama, ac, cte, Vola, Vol0
real*8  :: L, t, dt, x, r, tol
real*8  :: Ri, W, fmt, Ve, Fd, Cf, Fmed, Itot, Is, Ismed
real*8  :: pmed, pmax
real*8  :: p, pa, pao, paux, dp, Vtc, tq, z, g
real*8  :: Ai, Ae, Au, Voli, Vole, Volu, Re

real*8 :: Dg, Di, Da, Dgrao
real*8 :: r_ativ, t_inicio

real*8, dimension(:), allocatable :: Rativ

real*8, parameter :: pi = dacos(-1.d0)

integer :: passos, i, it, itmax, grao
integer :: ver, j

integer :: freq, passo
character*128 :: texto, arquivo
integer :: ativacao

t_inicio = 0.d0

i= system("Notepad entrada.txt")
! dados
open(1,file='entrada.txt')
read(1,*) Arquivo
read(1,*) Dg
read(1,*) Di
read(1,*) Da
read(1,*) L
read(1,*) grao
read(1,*) freq
read(1,*) r_ativ
read(1,*) gama
read(1,*) Rg
read(1,*) Tc
read(1,*) rop
read(1,*) a
read(1,*) n
read(1,*) pao
read(1,*) dt
read(1,*) passos
read(1,*) itmax
read(1,*) tol

! Volume inicial da câmara de combustão (m3)
Vol0 =    0.25d0*pi*Da**2*L

! Área de queima para o grão 1 (m2)
Ab   =    0.25d0*pi*Di**2  


! Área da garganta (m2)
At   =    0.25d0*pi*Dg**2  

! Volume total da câmara de combustão
Vtc  =    pi*0.25d0*Di**2*L

! L
! 1-Comprimento do grão (m)
! 2-Espessura do grão (m)
! 3-Metade da espessura do grão (m)
if(grao == 2) then
    W = L
    L = (Di-Da)*0.5d0
elseif(grao == 3) then
    W = L
    L = (Di-Da)*0.250d0
endif

! Raio interno do grão (m)
Ri   =    0.5d0*Da

! Raio externo do grão (m)
Re   =    0.5d0*Di

x      = 0.0d0
pa     = pao

ac  = dsqrt(gama*Rg*Tc)
cte = gama*((2.0d0/(gama+1.0d0))**((gama+1.0d0)/(2.0d0*(gama-1.0d0))))

! Aceleração gravitacional (m/s2)
g   = 9.80665d0

! para grão=4
if(grao==4) then

    ! Calculo do numero de discretizações para a ativação do grão-propelente
    ativacao = L/r_ativ/dt+1

    ! Discretização do raio do grão-propelente em 'ativacao' partes
    allocate(Rativ(ativacao))
    
    Rativ = 0.d0
    
    ! Dgrao = comprimento de cada pedaço de grão discretizado (m)
    Dgrao = L/ativacao
    ! Diferencas finitas!!
    ! i =1, x = 0
    ! i = ativacao+1, x= L
    write(*,*) Dgrao, ativacao
    
endif

open(8,file="./output/" // trim(adjustl(arquivo)) // ".txt")

write(8,10)
10 format("#", t8,"i", t15,"tempo (s)", t30,"Fd (N)",   t45,"r (mm/s)", t60,"p (atm)", &
                  t75,"K (adim.)", t90,"fmt (kg/s)", t105,"Ve (m/s)",              &
                 t120,"x (mm)",   t135,"Cf (adim.)", t150,"Is (s)",   &
                 t165,"dp (Pa)",   t175,"iterações")

! instante inicial

p   = pao
Vol = Vol0
r   = a*(p**n)
t   = 0.0d0
fmt = 0.0d0
Ve  = 0.0d0
Fd  = 0.0d0
Cf  = 0.0d0
Is  = 0.0d0
dp  = 0.0d0
Itot= 0.0d0
Fmed= 0.0d0
pmed= 0.0d0 
pmax= 0.0d0
it  = 0
i   = 0
if ( grao == 2 ) Ab = 2*pi*W*(Ri + x)

write(8,20) i, t, Fd, r*1000, p/101325.0d0, Ab/At, fmt, Ve, x*1000, Cf, Is, dp, it
20 format (i8, 11(1pe15.4), i8)

! combustão

select case (grao)

  case ( 1, 2, 3, 4 )
    z1 = 1.0d0/(Rg*Tc)
    z2 = rop*a
    z3 = cte*At/ac

end select

do i = 1, passos

   Vola = Vol

   p  = pa

   r  = a*(p**n)

   x  = x + r*dt

   t  = i * dt
   
   select case (grao)
         
    case(1)
      Vol = Vol0 + Ab*x

    case(2)
      Ab  = 2*pi*W*(Ri + x)
      Vol = Vol0 + pi*W*( ((Ri+x)**2) - (Ri**2) )
      
    case(3)
      Ai  = 2 * pi * (Ri+x) * (W-2*x)
      Ae  = 2 * pi * (Re-x) * (W-2*x)
      Au  = 2 * pi * ( ((Re-x)**2) - ((Ri+x)**2) )
      Ab  = Ai + Ae + Au
      Voli= pi * W * ( (Ri+x)**2 )
      Vole= pi * W * ( (Re**2) - ((Re-x)**2) )
      Volu= pi * 2 * x * ( ((Re-x)**2) - ((Ri+x)**2) )
      Vol = Vol0 + Voli + Vole + Volu

    case (4)
        ! Atualiza a área atual de cada pedaço de grão
        do j=1, ativacao
            ! Se for além do máximo permitido (Diâmetro externo do grão) zera o valor de Rativ
            if( Rativ(j) > Re) Rativ(j) = 0.d0
            
            ! Se for diferente de zero, aumenta o raio do pedaço do grão
            if( Rativ(j) /= 0.d0 ) Rativ(j) = Rativ(j) + r*dt
            
            ! Caso a ativação de todo o grão ainda não aconteceu:
            if(t <= L/r_ativ) then        
                ! Percebe qual pedaço de grão ainda não foi ativado
                if( Rativ(j) == 0.d0 .and. j== i ) Rativ(j) = Ri       
            endif
        enddo
        
        ! Calculo das áreas de queima e do volume atual
        Ab = 0.d0
        do j=2, ativacao
            Ab = Ab + 0.5d0*(Rativ(j-1)+Rativ(j))
        enddo
        
        Vol= 0.d0
        do j=2, ativacao 
            if(Rativ(j-1) /= 0.d0 .and. Rativ(j) /= 0.d0) then
                Vol= Vol+ ( 0.5d0*(Rativ(j-1)+Rativ(j)) )**2 - Ri**2
            elseif(i>j) then
                Vol = Vol + Re**2 - Ri**2
            endif
        enddo
        Ab  = Ab*2*pi*dgrao 
        Vol = Vol0 + Vol*pi*dgrao
        
!        write(*,'(I10,3(1PE15.3))') i, Ab, vol, Vtc


        
    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
            if(Rativ(j) /= 0) then
                WRITE(2,*) dgrao*(j-1), Rativ(J)
            endif
        ENDDO
        if(Rativ(j) /= 0.d0 ) WRITE(2,*) dgrao*(j-1), Rativ(J-1)
        close(2)
        endif
    endif
   end select

   do it = 1, itmax

      paux = p
      
      p   = (Vola/Vol)*pa + ( z2*Ab*(p**n) - z3*p ) * dt / (Vol*z1) 

      dp = dabs(p-paux)

      if ( dp < tol ) exit
      
      ! Pressurizando, ainda não gerou pressão suficiente para aumentar a pressão na câmara
      if(p<=pao) then
        p = pao
        exit
      endif
      

   end do
   
   if(p>pao) pmed = pmed + p*dt
   
   if(p>=pmax) pmax = p
   
   ! Caso a pressão seja a pressão atmosférica e a força maior que zero (despressurizou)   
   if(p ==pao .and. Fd>0.d0) exit
   
   ! Resgata o tempo no qual começou a gerar força de empuxo
   if(p >=pao .and. Fd==0.d0) t_inicio = t

   fmt = cte*At*p/ac
   Ve  = ( 2*gama*Rg*Tc * (1.0d0 - ((pao/p)**((gama-1.0d0)/gama))) &
       / (gama-1.0d0) ) ** 0.5d0
   Fd  = fmt*Ve
   Itot= Itot + Fd*dt
   Is  = Fd/(fmt*g)
   Cf  = Fd/(p*At)
   Ismed = Ismed + Is*dt

   write(8,20) i, t, Fd, r*1000, p/101325.0d0, Ab/At, fmt, Ve, x*1000, Cf, Is, dp, it

    select case (grao)
    case(1,2,3)
        if ( x > L ) exit
   end select

   pa = p

end do
! despressurização

z  = cte*At*Rg*Tc/(ac*Vtc)
pa = p
Ab = 0.0d0
r  = 0.0d0
dp = 0.0d0
it = 0
tq = t

do j = i+1, passos

   t  = j * dt

   p  = pa * dexp(-z*(t-tq))

   if ( p < pao ) exit

   fmt = cte*At*p/ac

   Ve  = ( 2*gama*Rg*Tc * (1.0d0 - ((pao/p)**((gama-1.0d0)/gama))) &
       / (gama-1.0d0) ) ** 0.5d0

   Fd  = fmt*Ve
   Itot= Itot + Fd*dt
   Is  = Fd/(fmt*g)
   Cf  = Fd/(p*At)
   Ismed = Ismed + Is*dt

   write(8,20) j, t, Fd, r*1000, p/101325.0d0, Ab/At, fmt, Ve, x*1000, Cf, Is, dp, it

end do

close(8)

pmed = pmed / (t-t_inicio)

Fmed  = Itot / (t-t_inicio)

Ismed = Ismed / (t-t_inicio)

open(8,file="./output/" // trim(adjustl(arquivo)) // "_resultados.txt")
write(8,30) tq, t-tq, t-t_inicio, Itot, Fmed, Pmed/101325.0d0, Pmax/101325.0d0, Ismed, Ismed*g
30 format(/, 5x, "Tempo de queima total (s)     = ", F15.4, &
         //, 5x, "Tempo de despressurização (s) = ", F15.4, &
         //, 5x, "Tempo total com empuxo (s)    = ", F15.4, &
         //, 5x, "Empuxo médio (N)              = ", F15.4, &
         //, 5x, "Impulso total (N.s)           = ", F15.4, &
         //, 5x, "Pressão média (atm)           = ", F15.4, &
         //, 5x, "Pressão máxima (atm)          = ", F15.4, &
         //, 5x, "Impulso específico (s)        = ", F15.4, &
         //, 5x, "Velocidade de ejeção efetiva média (m/s)  = ", F15.4  )
close(8)

!ver = system("Notepad saida.txt")

! GRAFICO DA FORÇA
open(2,file="./output/" // trim(adjustl(arquivo)) // '_forca.gnu')

write(2,*) "set term png nocrop enhanced size 1280, 720"
write(2,*) "set out './output/" // trim(adjustl(arquivo)) // "_forca.png'"
write(2,*) "set style data lines"
write(2,*) "set grid"
write(2,*) "set xlabel 't (s)'"
write(2,*) "set ylabel 'F (N)'"
write(2,*) "set time"
write(2,'(A)') "plot './output/" // trim(adjustl(arquivo)) // ".txt' using 2:3 notitle with lines linetype -1 linewidth 2"

close(2)

i = system(".\gnuplot\gnuplot ./output/" // trim(adjustl(arquivo)) // "_forca.gnu")

! GRAFICO DA PRESSAO
open(2,file="./output/" // trim(adjustl(arquivo)) // '_pressao.gnu')

write(2,*) "set term png nocrop enhanced size 1280, 720"
write(2,*) "set out './output/" // trim(adjustl(arquivo)) // "_pressao.png'"
write(2,*) "set style data lines"
write(2,*) "set grid"
write(2,*) "set xlabel 't (s)'"
write(2,*) "set ylabel 'pressão (Atm)'"
write(2,*) "set time"
write(2,'(A)') "plot './output/" // trim(adjustl(arquivo)) // ".txt' using 2:5 notitle with lines linetype -1 linewidth 2"

close(2)

i = system(".\gnuplot\gnuplot ./output/" // trim(adjustl(arquivo)) // "_pressao.gnu")

! GRAFICO DA RAZAO DE QUEIMA (Kn)
open(2,file="./output/" // trim(adjustl(arquivo)) // '_Kn.gnu')

write(2,*) "set term png nocrop enhanced size 1280, 720"
write(2,*) "set out './output/" // trim(adjustl(arquivo)) // "_Kn.png'"
write(2,*) "set style data lines"
write(2,*) "set grid"
write(2,*) "set xlabel 't (s)'"
write(2,*) "set ylabel 'Kn (adim)'"
write(2,*) "set time"
write(2,'(A)') "plot './output/" // trim(adjustl(arquivo)) // ".txt' using 2:6 notitle with lines linetype -1 linewidth 2"

close(2)

i = system(".\gnuplot\gnuplot ./output/" // trim(adjustl(arquivo)) // '_Kn.gnu')

! GRAFICO Do fluxo de massa
open(2,file="./output/" // trim(adjustl(arquivo)) // '_fm.gnu')

write(2,*) "set term png nocrop enhanced size 1280, 720"
write(2,*) "set out './output/" // trim(adjustl(arquivo)) // "fm.png'"
write(2,*) "set style data lines"
write(2,*) "set grid"
write(2,*) "set xlabel 't (s)'"
write(2,*) "set ylabel 'Fluxo de Massa (kg/s)'"
write(2,*) "set time"
write(2,'(A)') "plot './output/" // trim(adjustl(arquivo)) // ".txt' using 2:7 notitle with lines linetype -1 linewidth 2"

close(2)

i = system(".\gnuplot\gnuplot ./output/" // trim(adjustl(arquivo)) // '_fm.gnu')

! GRAFICO Da velocidade de ejeção
open(2,file="./output/" // trim(adjustl(arquivo)) // '_ve.gnu')

write(2,*) "set term png nocrop enhanced size 1280, 720"
write(2,*) "set out './output/" // trim(adjustl(arquivo)) // "ve.png'"
write(2,*) "set style data lines"
write(2,*) "set grid"
write(2,*) "set xlabel 't (s)'"
write(2,*) "set ylabel 'Velocidade de Ejeção (m/s)'"
write(2,*) "set time"
write(2,'(A)') "plot './output/" // trim(adjustl(arquivo)) // ".txt' using 2:8 notitle with lines linetype -1 linewidth 2"

close(2)

i = system(".\gnuplot\gnuplot ./output/" // trim(adjustl(arquivo)) // "_ve.gnu")


if(freq>0 .and. grao==4) then
    OPEN(2,FILE='./temp/grao.txt')
    write(2,*) 0.d0, Ri
    write(2,*) L, Ri
    write(2,*) L, Re
    write(2,*) 0.d0, Re
    write(2,*) 0.d0, Ri
    close(2)

    open(2,file='./output/' // trim(adjustl(arquivo)) // 'plot_gif.gnu')

    write(2,*) "set terminal gif animate delay 10 size 1280, 720"
    write(2,*) "set out './output/" // trim(adjustl(arquivo)) // "progresso.gif'"
    write(2,*) "set style data lines"
    write(2,*) "set grid"
    write(2,*) "set xlabel 'x (m)'"
    write(2,*) "set ylabel 'r (m)'"
    write(2,'(A,F13.4,A,F13.4,A)') "set xrange [", -0.001d0, ":", L+.001d0, "]"
    write(2,'(A,F13.4,A,F13.4,A)') "set yrange [", -0.001d0, ":", Re+.001d0, "]"


    do i=1, passo
        if(i<10) then
            write(texto,'("0000",I1)') i
        elseif(i<100) then
            write(texto,'("000",I2)') i
        elseif(i<1000) then
            write(texto,'("00",I3)') i
        elseif(i<10000) then
            write(texto,'("0",I4)') i
        elseif(i<100000) then
            write(texto,'(I5)') i
        endif
        
        write(2,*) "plot './temp/"// trim(adjustl(texto)) //".txt' using 1:2 title 'Frente de chama' \"
        write(2,*) ", './temp/grao.txt' using 1:2 title 'Grão-propelente' lw 2 "

    enddo
    close(2)
    i = system(".\gnuplot\gnuplot ./output/" // trim(adjustl(arquivo)) // "plot_gif.gnu")
endif

!ver = system("Notepad saida2.txt")

ver = system("del .\temp\*.txt")


end Program empuxo1p1
