Program refuerzo

real(8) Tm,Tc,Tf,k,Cp,rhol,f,k0,deltaH,rhoc,lx1,lx2, ly1, ly2, dx, dy, dt, r,phi1, phi2,kc1,kc2,E,phimax,g,rechupe,contraccion
integer nx, ny, i, j, z, tiempofinal,ciclo,b,nodosl1,nodosl2
real(8),allocatable, dimension (:,:)::T, AUXT,C
real(8),allocatable,dimension(:)::VT,Ve

!---------------------variables Geometricas-----------------------------
lx1=0.195/2;lx2=0.005
ly1=0.005;ly2=0.005         !Dimensiones del refuerzo (en metros)
dx=0.001; dy=0.001          !Distancia entre nodos(en metros)
dt=1.
tiempofinal=1000

nx=floor((2*lx1+lx2)/dx+0.5);
ny=floor((ly1+ly2)/dy+0.5)  !Cantidad de nodos
nx1=floor(lx1/dx+0.5);nx2=floor(lx2/dx+0.5)
ny1=floor(ly1/dy+0.5);ny2=floor(ly2/dy+0.5)



!--------------------Variables de Proceso-------------------------------
Tm=50+273                   !Temp molde (K)
Tc=110+273                  !Temp cristalizacion(K)
Tf=230+273                  !Temp Fundido(K)

k=0.46                      !Conduct termica(W/m.K)
Cp=2475                     !(J/Kg.k)
rhol=867                    !Densidad del liquido(Kg/m³)


!--------------------Parametros de Abrami-------------------------------
f=2.2506                    !Es el 'n' de abrami
k0=3.014E9
deltaH=188370               !delta de cristalizacion(J/Kg)
rhoc=983                    !Densidad del solido(Kg/m³)
E=14575                     !Energia de cristalizacion (J/mol)²
phimax=0.31

r=k*dt/(rhol*Cp*dx**2)

allocate(T(ny,nx), AUXT(ny,nx),C(ny,nx))
Call Forma1(T,Tc,Tf,nx,ny,nx1,nx2,ny1,ny2)
Call GraficaForma(T,nx,ny)
z=tiempofinal
ciclo=0
Do i=1,ny
 do j=1,nx
  C(i,j)=0.
 end do
end do

Do while (z==tiempofinal)
AUXT(:,:)=T(:,:)
ciclo=ciclo+1
Do i=2, ny1-1
 Do j=2, nx-1
  T(i,j)=r*AUXT(i+1,j)+r*AUXT(i-1,j)+(-4*r+1)*AUXT(i,j)+r*AUXT(i,j+1)+r*AUXT(i,j-1)
  if(T(i,j).le.Tc) then
   if (C(i,j)==0.) then
    kc2=k0*exp(-E/(8.314*(Tf-T(i,j))))
    phi2=1-exp(-kc2*(dt*ciclo)**f)
    C(i,j)=1.
    end if
  end if
  if(phi2.le.phimax) then
    kc1=k0*exp(-E/(8.314*(Tf-AUXT(i,j))))
    phi1=1-exp(-kc1*(dt*(ciclo-1))**f)
    g=(rhoc/(rhol*Cp))*deltaH*(phi2-phi1)
    T(i,j)=T(i,j)+g
  end if
  end do
end do
Do i=ny1,ny-1
 Do j=nx1+1, (nx1+nx2)-1
  T(i,j)=r*AUXT(i+1,j)+r*AUXT(i-1,j)+(-4*r+1)*AUXT(i,j)+r*AUXT(i,j+1)+r*AUXT(i,j-1)
  if(T(i,j).le.Tc) then
   if (C(i,j)==0) then
    kc2=k0*exp(-E/(8.314*(Tf-T(i,j))))
    phi2=1-exp(-kc2*(dt*ciclo)**f)
    C(i,j)=1
    end if
  end if
  if(phi2.le.phimax) then
    kc1=k0*exp(-E/(8.314*(Tf-AUXT(i,j))))
    phi1=1-exp(-kc1*(dt*(ciclo-1))**f)
    g=(rhoc/(rhol*Cp))*deltaH*(phi2-phi1)
    T(i,j)=T(i,j)+g
  end if
 end do
end do

If ((T((floor((ny1/2)+0.5))+1,nx1))<= Tc) then
Call GraficaForma(T,nx,ny)
z=tiempofinal+1
end if
end do

nodosl1=0
nodosl2=0
Do i=2, ny1-1
 Do j=2, nx-1
    if(T(i,j)>Tc) then
    nodosl1=nodosl1+1

  end if
 end do
end do
Do i=ny1,ny-1
 Do j=nx1+1, (nx1+nx2)-1
  if(T(i,j)>Tc) then

    nodosl2=nodosl2+1
  end if
 end do
end do
allocate(VT(nodosl1+nodosl2),Ve(nodosl1+nodosl2))
open(unit=4, file='Temperaturas_nodos_liquidos.dat')
b=1
Do i=2, ny1-1
 Do j=2, nx-1
    if(T(i,j)>Tc) then
    VT(b)=T(i,j)-273
    write(4,'(F15.3)')VT(b)
    b=b+1

  end if
 end do
end do


Do i=ny1,ny-1
 Do j=nx1+1, (nx1+nx2)-1
  if(T(i,j)>Tc) then
    VT(b)=T(i,j)-273
    write(4,'(F15.3)')VT(b)
    b=b+1

  end if
 end do
end do

close(4)



Call Colorear(T,nx,ny)
call system ("gnuplot -persist 'script.p'")
rechupe=0.
Ve=0.
Call Calculodey(VT,nodosl1,nodosl2,Ve)
Call Calcularechupe(nodosl1,nodosl2,rechupe,Ve)
contraccion=100*rechupe/2000

  write(*,*) ny1, nx1, floor(ny1/2 +0.5),nodosl1,nodosl2,rechupe
  write(*,*)'La contraccion volumetrica es',contraccion

Contains

Subroutine Forma1(T,Tc,Tf,nx,ny,nx1,nx2,ny1,ny2) !Esta subroutina arma la matriz de temperaturas con forma T
  real(8) Tc,Tf
  integer nx,ny,nx1,nx2,ny1,ny2
  real(8),dimension(ny,nx)::T

  T(:,:)=Tm
  T(2:ny/2,1)=Tc
  T(2:ny/2,nx)=Tc
  do i=2,ny/2
    T(i,2:nx-1)=Tf
  end do
  do i=ny1,ny1+ny2-1
    T(i,nx1+1:nx1+nx2-1)=Tf
  end do
 T(floor(ny/2+0.5),1:nx1)=Tm
 T(floor(ny/2+0.5),(nx1+nx2):nx)=Tm
end Subroutine
Subroutine GraficaForma(T,nx,ny) !escribe la matriz forma en un archivo
  integer nx,ny
  real(8),dimension(ny,nx)::T

  open(unit=2, file='forma1.dat')
  Do i=1,ny
    Do j=1,nx
      write(2,'(200F8.3)',advance='no')T(i,j)
    end do
    write(2,*)
  end do
  write(2,*)
  Do i=1,nx
    write(2,'(200I8)',advance='no') i
  end do

  Close(2)
end subroutine

Subroutine Colorear(T,nx,ny)
integer i,j
real(8)::T(ny,nx)

open(unit=3, file='grafico_de_temperaturas.dat')
! write(3,'(2A4,A15)')'x','y','Temperaturas'
  Do j=1,nx
   Do i=1,ny
    write(3,'(2I4,F115.3)')j,i,T(i,j)
    end do
    write(3,*)
   end do
   Close(3)
end subroutine

Subroutine Calcularechupe(nodosl1,nodosl2,rechupe,Ve)
integer nodosl1,nodosl2,i
real(8) rechupe,m,Vesolido
real(8):: Ve(nodosl1+nodosl2)


m=1                !gr
Vesolido=1.0513    !cm³/gr
rechupe=0.
Do i=1,(nodosl1+nodosl2)

rechupe=Ve(i)*m-Vesolido*m+rechupe

end do

end subroutine


Subroutine Calculodey (VT,nodosl1,nodosl2,Ve)
integer nodosl1,nodosl2,i
real(8):: VT(nodosl1+nodosl2),Ve(nodosl1+nodosl2)
real(8) suma,suma2,suma3

do i=1,nodosl1+nodosl2
suma=0.
suma2=0.
suma3=0.
suma=-0.0000000000000000000407514*(VT(i)**10)+0.0000000000000000000000222233*(VT(i)**11)
suma2=-1.75369+0.35938*VT(i)-0.0190579*(VT(i)**2)+0.00055385*(VT(i)**3)-0.00000984786*(VT(i)**4)
suma3=0.000000113272*(VT(i)**5)+0.00000000000445019*(VT(i)**7)-0.0000000000000151411*(VT(i)**8)+0.0000000000000000327439*(VT(i)**9)
Ve(i)=-0.00000000086681*(VT(i)**6)+suma+suma2+suma3

end do
Write(*,*)'Los valor de volumen especifico es=',Ve
end subroutine


End