temperature.f90

运载火箭用拉瓦尔喷管计算（此程序参数保密

subroutine temperature
use module1
real ::time,Bdiff(nz),beta(nz)
write(10,"(32(a10,','))")"tw-out(1)","tw-out(2)","tw-out(3)","tw-out(4)","tw-out(5)","tw-out(6)","tw-out(7)","tw-out(8)","tw-out(9)","tw-out(10)","tw-out(11)","tw-out(12)","tw-out(13)","tw-out(14)","tw-out(15)","tw-out(16)","tw-out(17)","tw-out(18)","tw-out(19)","tw-out(20)","tw-out(21)","tw-out(22)","tw-out(23)","tw-out(24)","tw-out(25)","tw-out(26)","tw-out(27)","tw-out(28)","tw-out(29)","tw-out(30)","tw-out(31)","time"
write(20,"(32(a10,','))")"tw-in(1)","tw-in(2)","tw-in(3)","tw-in(4)","tw-in(5)","tw-in(6)","tw-in(7)","tw-in(8)","tw-in(9)","tw-in(10)","tw-in(11)","tw-in(12)","tw-in(13)","tw-in(14)","tw-in(15)","tw-in(16)","tw-in(17)","tw-in(18)","tw-in(19)","tw-in(20)","tw-in(21)","tw-in(22)","tw-in(23)","tw-in(24)","tw-in(25)","tw-in(26)","tw-in(27)","tw-in(28)","tw-in(29)","tw-in(30)","tw-in(31)","time"
write(30,"(32(a10,','))")"r0(1)","r0(2)","r0(3)","r0(4)","r0(5)","r0(6)","r0(7)","r0(8)","r0(9)","r0(10)","r0(11)","r0(12)","r0(13)","r0(14)","r0(15)","r0(16)","r0(17)","r0(18)","r0(19)","r0(20)","r0(21)","r0(22)","r0(23)","r0(24)","r0(25)","r0(26)","r0(27)","r0(28)","r0(29)","r0(30)","r0(31)","time"
write(40,"(32(a10,','))")"alf(1)","alf(2)","alf(3)","alf(4)","alf(5)","alf(6)","alf(7)","alf(8)","alf(9)","alf(10)","alf(11)","alf(12)","alf(13)","alf(14)","alf(15)","alf(16)","alf(17)","alf(18)","alf(19)","alf(20)","alf(21)","alf(22)","alf(23)","alf(24)","alf(25)","alf(26)","alf(27)","alf(28)","alf(29)","alf(30)","alf(31)","time"
write(50,"(32(a10,','))")"s(1)","s(2)","s(3)","s(4)","s(5)","s(6)","s(7)","s(8)","s(9)","s(10)","s(11)","s(12)","s(13)","s(14)","s(15)","s(16)","s(17)","s(18)","s(19)","s(20)","s(21)","s(22)","s(23)","s(24)","s(25)","s(26)","s(27)","s(28)","s(29)","s(30)","s(31)","time"
write(60,"(32(a10,','))")"dsdt(1)","dsdt(2)","dsdt(3)","dsdt(4)","dsdt(5)","dsdt(6)","dsdt(7)","dsdt(8)","dsdt(9)","dsdt(10)","dsdt(11)","dsdt(12)","dsdt(13)","dsdt(14)","dsdt(15)","dsdt(16)","dsdt(17)","dsdt(18)","dsdt(19)","dsdt(20)","dsdt(21)","dsdt(22)","dsdt(23)","dsdt(24)","dsdt(25)","dsdt(26)","dsdt(27)","dsdt(28)","dsdt(29)","dsdt(30)","dsdt(31)","time"
write(70,"(32(a10,','))")"dsdt1(1)","dsdt1(2)","dsdt1(3)","dsdt1(4)","dsdt1(5)","dsdt1(6)","dsdt1(7)","dsdt1(8)","dsdt1(9)","dsdt1(10)","dsdt1(11)","dsdt1(12)","dsdt1(13)","dsdt1(14)","dsdt1(15)","dsdt1(16)","dsdt1(17)","dsdt1(18)","dsdt1(19)","dsdt1(20)","dsdt1(21)","dsdt1(22)","dsdt1(23)","dsdt1(24)","dsdt1(25)","dsdt1(26)","dsdt1(27)","dsdt1(28)","dsdt1(29)","dsdt1(30)","dsdt1(31)","time"
write(80,"(32(a10,','))")"pe(1)","pe(2)","pe(3)","pe(4)","pe(5)","pe(6)","pe(7)","pe(8)","pe(9)","pe(10)","pe(11)","pe(12)","pe(13)","pe(14)","pe(15)","pe(16)","pe(17)","pe(18)","pe(19)","pe(20)","pe(21)","pe(22)","pe(23)","pe(24)","pe(25)","pe(26)","pe(27)","pe(28)","pe(29)","pe(30)","pe(31)","time"
write(90,"(32(a10,','))")"te(1)","te(2)","te(3)","te(4)","te(5)","te(6)","te(7)","te(8)","te(9)","te(10)","te(11)","te(12)","te(13)","te(14)","te(15)","te(16)","te(17)","te(18)","te(19)","te(20)","te(21)","te(22)","te(23)","te(24)","te(25)","te(26)","te(27)","te(28)","te(29)","te(30)","te(31)","time"
write(100,"(32(a10,','))")"ma(1)","ma(2)","ma(3)","ma(4)","ma(5)","ma(6)","ma(7)","ma(8)","ma(9)","ma(10)","ma(11)","ma(12)","ma(13)","ma(14)","ma(15)","ma(16)","ma(17)","ma(18)","ma(19)","ma(20)","ma(21)","ma(22)","ma(23)","ma(24)","ma(25)","ma(26)","ma(27)","ma(28)","ma(29)","ma(30)","ma(31)","time"
write(110,"(32(a10,','))")"A(1)","A(2)","A(3)","A(4)","A(5)","A(6)","A(7)","A(8)","A(9)","A(10)","A(11)","A(12)","A(13)","A(14)","A(15)","A(16)","A(17)","A(18)","A(19)","A(20)","A(21)","A(22)","A(23)","A(24)","A(25)","A(26)","A(27)","A(28)","A(29)","A(30)","A(31)","time"



time=0.0
do j=1,nz
   s(j)=0.0
   crev(j)=0.0
   do i=1,nr
      mtrt(i,j)=tini
      mtrtd(i,j)=tini
   end do
   mtrt(nr1,j)=tini
end do

do while(time<=1)
   call para1()
                                                                 write(30,"(32(f10.5,','))") (r0(j),j=1,nz),time
!-------------------划分新网格------------------------------------------
   do kk=1,nz
      dr(kk)=(ra-r0(kk))/float(Nr)
      do i=1,nr1
         ru(nr1+1-i,kk)=(i-1)*dr(kk)+r0(kk)
      end do
      do i=1,nr
         rp(i,kk)=(ru(i,kk)+ru(i+1,kk))/2.
      end do
      mtrt(nr1,kk)=(mtrt(nr,kk)+te(kk)*alf(kk)*dr(kk)/lanbt+crev(kk)*dr(kk)/lanbt)/(1+alf(kk)*dr(kk)/lanbt) 
      do i=2,nr
         mtrt(i,kk)=(mtrtd(i,kk)+mtrtd(i-1,kk))/2   
      end do
      mtrt(1,kk)=mtrtd(1,kk)
      do i=2,nr1
         dv(i-1,kk)=0.5*(ru(i,kk)+ru(i-1,kk))*dr(kk)
         an(i-1,kk)=ru(i-1,kk)*lanbt/dr(kk)
         as(i-1,kk)=ru(i,kk)*lanbt/dr(kk)
         ap0(i-1,kk)=rou*cmat*dv(i-1,kk)/dt
         ap(i-1,kk)=an(i-1,kk)+as(i-1,kk)+ap0(i-1,kk)
         b(i-1,kk)=ap0(i-1,kk)*(mtrt(i-1,kk)+mtrt(i,kk))/2
      end do
      call solve()   
   end do
                                                                  write(10,"(32(f10.5,','))") (mtrtd(1,j),j=1,nz),time
!--------------------网格点温度赋新值------------------------------------


   !------------更新温度参数---------------
   call para2()
  write(40,"(32(f10.5,','))") (alf(j),j=1,nz),time
  write(80,"(32(g10.5,','))") (pe(j),j=1,nz),time
  write(90,"(32(f10.5,','))") (te(j),j=1,nz),time
  write(100,"(32(f10.5,','))") (ma(j),j=1,nz),time
  write(110,"(32(f10.5,','))") (A(j),j=1,nz),time
   do j=1,nz
      tw(j)=mtrtd(nr,j)
   end do
                                                                    write(20,"(32(f10.5,','))") (tw(j),j=1,nz),time
   call solve6()  

   do j=1,nz
      crev(j)=-dmh2o(j)*dh1-dmco2(j)*dh2+qcor(j) 
   end do  
   time=time+dt
   nctr=nctr+1
   do j=1,nz
      s(j)=s(j)+dt*dmdt(j)/rou
   end do
															        write(50,"(32(f10.5,','))") (s(j),j=1,nz),time
                                                                    write(60,"(32(f10.5,','))") (dmdt(j)*1e3/rou,j=1,nz),time
																	
!==========================================
   do j=1,nz
      beta(j)=alf(j)/Cp
      bdiff(j)=Mc*rh2oe/Mh2o+Mc*rco2e/Mco2
      dmdt1(j)=bdiff(j)*beta(j)
   end do
                                                                    write(70,"(32(f10.5,','))") (dmdt1(j)*1e3/rou,j=1,nz),time

!==========================================
   do j=1,nz
      if((r0(j)+s(j))>ra)then 
        exit
      end if           
   end do
   write(*,*)time
end do     
end subroutine