Question

Develop Fortran code for double stage pulse tube cryocooler to achieve 20 k temperature using online configuration with regenerative matrix heat exchanger...

Develop Fortran code and plot the results in Sicilian for validation

Answer 1

C double stage pulse tube cryocooler to achieve 20 k temperature using online configuration with regenerative matrix heat exchanger.
   IMPLICIT NONE
   integer::t1,t2,Re,Re1,Re2
   Real::p1,p2,Vf1,Vf2,Vg1,Vg2,V1,V2,hf1,hf2,hg1,hg2,h1,
*h2,uf1,uf2,u1,u2,ug1,ug2,Area,Vel1,Vel2,x1,x2,x,d,a,b,c,fm,f1,f2,
*inc_length1,inc_length2,cum_length1,cum_length2,Vm,
*e=2.718281828,pie=3.141592654,dia=1.63E-3,w=0.010,tc=40,te=5,
* t,p,h,V,u,inc_length, cum_length
   Area=(pie*dia**2)/4
   print*, "Area=", Area
   d=(4*w)/(pie*dia**2)
   print*, "d=", d

   t1=40
   p1=1000*E**(15.06-(2418.4/(t1+273.15)))
   p=p1/1000
   print*, "Pressure=", p1/1000
   Vf1=(0.777+0.002062*t1+0.00001608*t1**2)/1000
   PRINT*, "specific volume of saturated liquid=", Vf1
   hf1= 200+1.172*t1+0.001854*t1**2
   print*, "enthalpy of saturated liquid=", hf1
   uf1=0.0002367-(1.715E-6*t1)+(8.869E-9*t1**2)
   print*, "viscosity of saturated liquid=", uf1
   V1=Vf1
   h1=hf1
   u1=uf1
   Vel1=d*V1
   print*, "velocity=", Vel1
   Re1=(Vel1*dia)/(V1*u1)
   print*, "Reynolds Number=", Re1
   f1=0.33/(Re1)**0.25
   print*, "Friction Factor=", f1

   inc_length1=0
   cum_length1=0
  
      
   do 10 t2=39,5,-1
   print*, "t2=",t2  
   p2=1000*E**(15.06-(2418.4/(t2+273.15)))
   p=p2/1000
C   print*, "Pressure=", p2/1000
   Vf2=(0.777+0.002062*t2+0.00001608*t2**2)/1000
c   PRINT*, "specific volume of saturated liquid=", Vf2
   Vg2=(-4.26+(94050*(t2+273.15))/p2)/1000
c   print*, "Specific Volume of saturated vapour=", Vg2
   hf2= 200+(1.172*t2)+(0.001854*t2**2)
c   print*, "enthalpy of saturated liquid=", hf2
   hg2=405.5+(0.3636*t2)-(0.002273*t2**2)
c   print*, "enthalpy of saturated vapour=", hg2
   uf2=(0.0002367)-(1.715E-6*t2)+(8.869E-9*t2**2)
c   print*, "viscosity of saturated liquid=", uf2
   ug2=(11.945E-6)+(50.06E-9*t2)+(0.2560E-9*t2**2)
c   print*, "viscosity of saturated vapour=", ug2
c
   a=0.5*(Vg2-Vf2)*(Vg2-Vf2)*d**2
c   print*, "a=", a
   b=1000.0*(hg2-hf2)+(Vf2*(Vg2-Vf2)*d**2)
c   print*, "b=", b
   c=1000.0*(hf2-h1)+((0.5*d**2)*(Vf2**2))-(0.5*Vel1**2)
c   print*, "c=", c
   x1=(-b+ SQRT(b**2-4*a*c))/(2*a)
   x2=(-b- SQRT(b**2-4*a*c))/(2*a)
   x=max(x1,x2)
   print*, "dryness fraction=", x

   h2=hf2*(1-x)+x*hg2
   print*, "enthalpy=", h2
   V2=Vf2*(1-x)+x*Vg2
   print*, "specific volume=", V2
   u2=uf2*(1-x)+x*ug2
   print*, "viscosity=", u2
   Vel2=d*V2
   print*, "velocity=", Vel2
   Re2=(Vel2*dia)/(V2*u2)
   print*, "Reynolds Number=", Re2
   f2=0.33/(Re2)**0.25
   print*, "Friction Factor=", f2

   fm=(f1+f2)/2
   print*, "fm=", fm
   Vm=0.5*(Vel1+Vel2)
   print*, "Vm=", Vm
   inc_length2=(2*Area*dia*((p1-p2)-(w*(V2-V1))))/(fm*Vm*w)
   print*, "inc_length=", inc_length2
   cum_length2=cum_length1+inc_length2
   print*, "cum_length=", cum_length2

   p1=p2
   Vf1=Vf2
   Vg1=Vg2
   V1=V2
   hf1=hf2
   hg1=hg2
   h1=h2
   uf1=uf2
   ug1=ug2
   u1=u2
   Vel1=Vel2
   Re1=Re2
   f1=f2
   inc_length1=inc_length2
   cum_length1=cum_length2
   OPEN (UNIT=1, FILE='RESULTS1', STATUS='OLD')
   OPEN (UNIT=2, FILE='RESULTS2', STATUS='OLD')
   WRITE (2,*) t1,p,x,V1,h1,Vel1,inc_length1, cum_length1
   CLOSE(2)
   WRITE (1,*) t2,p,x,V2,h2,Vel2,inc_length2, cum_length2
10   continue
   close(1)
  
  
   stop
   end

C LENGTH OF CAPILLARY TUBE ID=1.63mm R-22 Tc=40 Te=5 w=0.010kg/s
C   AT THE INLET OF CAPILLARY TUBE
   IMPLICIT NONE
   integer::i,m
   Real :: Area, d, a(100), b(100), c(100),l,j(100), k(100),
* Vf(100),Vg(100),hf(100),hg(100),uf(100),ug(100),Re(100),
* Vel(100),f(100),V(100), h(100), u(100), x(100), Vm(100),
* fm(100), t(100), z(100), y(100), p(100),
* inc_length(100), cum_length(100)
   REAL:: E=2.718281828, pie=3.141592654, dia=1.63E-3, w=0.010
   t(1)=40
   print*, "t(1)=", t(1)
   z(1)=2418.4/(t(1)+273.15)
   y(1)=15.06-z(1)
   p(1)=1000*E**y(1)
   print*, "Pressure=", p(1)/1000
   x(1)=0

   Vf(1)=(0.777+0.002062*t(1)+0.00001608*t(1)**2)/1000
   V(1)=Vf(1)
   PRINT*, "specific volume of saturated liquid=", V(1)
   hf(1)= 200+1.172*t(1)+0.001854*t(1)**2
   h(1)=hf(1)
   print*, "enthalpy of saturated liquid=", h(1)
   uf(1)=0.0002367-(t(1)*1.715E-6)+(t(1)*t(1)*8.869E-9)
   u(1)=uf(1)
c   print*, "viscosity of saturated liquid=", u(1)
   Area=(pie*dia**2)/4
c   print*, "Area=", Area
   d= (w*4)/(pie*dia**2)
   Vel(1)=d*Vf(1)
   print*, "velocity=", Vel(1)
   Re(1)=(Vel(1)*dia)/(Vf(1)*uf(1))
c   print*, "Reynolds Number=", Re(1)
   f(1)=0.33/(Re(1))**0.25
c   print*, "Friction Factor=", f(1)
   inc_length(1)=0
   cum_length(1)=0

   do 10 i=1,35,1

   t(i+1)=t(i)-1
   print*, "t(i+1)=", t(i+1)
   m=i
   m=m+1

   z(m)=2418.4/(t(i+1)+273.15)
   y(m)=15.06-z(m)
   p(m)=1000*E**y(m)
   print*, "Pressure=", p(m)/1000.0
  
   Vf(m)=(0.777+0.002062*t(i+1)+0.00001608*t(i+1)*t(i+1))/1000.0
c   print*, "Vf=" ,Vf(m)
   Vg(m)= (-4.26+94050*(t(i+1)+273.15)/p(m))/1000.0
c   print*, "Vg=", Vg(m)
   hf(m)=200.0+1.172*t(i+1)+0.001854*t(i+1)**2
c   print*, "hf=", hf(m)
   hg(m)=405.5+0.3636*t(i+1)-0.002273*t(i+1)**2
c   print*, "hg=", hg(m)
   uf(m)=0.0002367-(t(i+1)*1.715E-6)+(t(i+1)*t(i+1)*8.869E-9)
c   print*, "uf=", uf(m)
   ug(m)= (11.945E-6)+(t(i+1)*50.06E-9)+(t(i+1)*t(i+1)*0.2560E-9)
c   print*, "ug=", ug(m)

   a(m)=0.5*(Vg(m)-Vf(m))*(Vg(m)-Vf(m))*d**2
c   print*, "a=", a(m)
   b(m)=1000.0*(hg(m)-hf(m))+Vf(m)*(Vg(m)-Vf(m))*d**2
c   print*, "b=", b(m)
   c(m)=1000.0*(hf(m)-h(m-1))+(0.5*d*d*Vf(m)**2)-(0.5*Vel(m-1)**2)
c   print*, "c=", c(m)
   j(m)=(-b(m)+ SQRT((b(m)*b(m))-(4*a(m)*c(m))))/(2*a(m))
   k(m)=(-b(m)- SQRT((b(m)*b(m))-(4*a(m)*c(m))))/(2*a(m))
   x(m)= max(j(m), k(m))

   print*, "dryness fraction=", x(m)
   h(m)=hf(m)*(1-x(m))+x(m)*hg(m)
   print*, "enthalpy=", h(m)
   V(m)=Vf(m)*(1-x(m))+x(m)*Vg(m)
   print*, "specific volume=", V(m)
   u(m)=uf(m)*(1-x(m))+x(m)*ug(m)
c   print*, "viscosity=", u(m)
   Vel(m)=d*V(m)
   print*, "velocity=", Vel(m)
   Re(m)=(Vel(m)*dia)/(V(m)*u(m))
c   print*, "Reynolds Number=", Re(m)
   f(m)=0.33/(Re(m))**0.25
c   print*, "friction factor=", f(m)
  
   fm(m)=(f(m-1)+f(m))/2
c   print*, "fm=", fm(m)
   Vm(m)=0.5*(Vel(m-1)+Vel(m))
c   print*, "Vm=", Vm(m)

   inc_length(m)=2*Area*dia*((p(m-1)-p(m)-(w*(V(m)-V(m-1))))
*   /(fm(m)*Vm(m)*w))
   print*, "inc_length=", inc_length(m)
   cum_length(m)=cum_length(m-1)+inc_length(m)
   print*, "cum_length=", cum_length(m)
10 continue
   stop
   end

C PROGRAM FOR CALCULATION OF LENGTH OF CAPILLARY TUBE(R-22)

   IMPLICIT NONE

   Real :: t,p,Vf,Vg,V,hf,hg,h,uf,ug,u,Area,Vel,x1,x2,x,d,a,b,c,Re,
*fm,f,inc_length. cum_length

   REAL::e=2.718281828,pie=3.141592654,dia=1.63E-3,w=0010,tc=40,te=5
  
   t=tc

10   p=1000*E**(15.06-(2418.4/(t+273.15)))
   print*, "Pressure=", p/1000

   Vf=(0.777+0.002062*t+0.00001608*t**2)/1000
   PRINT*, "specific volume of saturated liquid=", Vf

   Vg=(-4.26+(94050*(t+273.15))/p)/1000
   print*, "Specific Volume of saturated vapour=", Vg

   hf= 200+1.172*t+0.001854*t**2
   print*, "enthalpy of saturated liquid=", hf

   hg=405.5+0.3636*t-0.002273*t**2
   print*, "enthalpy of saturated vapour=", hg

   uf=0.0002367-(1.715E-6*t)+(8.869E-9*t**2)
   print*, "viscosity of saturated liquid=", uf

   ug=11.945E-6+50.06E-9*t+0.2560E-9*t**2
   print*, "viscosity of saturated vapour=", ug

   Area=(pie*dia**2)/4
   print*, "Area=", Area

   d=(4*w)/(pie*dia**2)
   print*, "d=", d
  
   t=t-1

   V=d*V1
   print*, "velocity=", V

   a=0.5*(Vg-Vf)*(Vg-Vf)*d**2
   print*, "a=", a
   b=1000.0*(hg-hf)+Vf*(Vg-Vf)*d**2
   print*, "b=", b
   c=1000.0*(hf-h)+(0.5*d*d*Vf**2)-(0.5*V**2)
   print*, "c=", c
   x1=(-b+ SQRT(b**2-4*a*c))/(2*a)
   x2=(-b- SQRT(b**2-4*a*c))/(2*a)
   x=max(x1,x2)
   print*, "dryness fraction=", x

   h=hf*(1-x)+x*hg
   print*, "enthalpy=", h
   V=Vf*(1-x)+x*Vg
   print*, "specific volume=", V
   u=uf*(1-x)+x*ug
   print*, "viscosity=", u

   Re=(V*dia)/(Vf*uf)
   print*, "Reynolds Number=", Re
   f=0.33/(Re)**0.25
   print*, "Friction Factor=", f
c   inc_length=
  
  
  
   Vel=d*V
   print*, "velocity=", Vel
   Re=(Vel*dia)/(V*u)
   print*, "Reynolds Number=", Re
   f=0.33/(Re)**0.25
   print*, "friction factor=", f
  
   fm=(f+f)/2
   print*, "fm=", fm
   Vm=0.5*(Vel+Vel)
   print*, "Vm=", Vm

   inc_length=2*Area*dia*(p-p)-w*(V-V)/(fm*Vm*w)
   print*, "inc_length=", inc_length
   cum_length=inc_length+inc_length
   print*, "cum_length=", cum_length
10 continue
   stop
   end

C LENGTH OF CAPILLARY TUBE ID=1.63mm R-22 Tc=40 Te=5 w=0.010kg/s
C   AT THE INLET OF CAPILLARY TUBE
   IMPLICIT NONE
   Double Precision :: t1, z,y,p1, Vf1,V1, Vg1, hf1, hg1, uf1, ug1, d,
*Vel1, Re1,f1, t, p, Vf, Vg, hf, hg, uf, ug, Vel, Re, f, g, x, a,
* b,c, u, v, h,fm
   REAL:: E=2.718281828, pie=3.141592654, dia=1.63E-3, w=0.010
c   print*, "enter the value of temperature"
   t1=40
   z=2418.4/(t1+273.15)
   y=15.06-z
   p1=1000*E**y
   print*, "Pressure=", p1/1000
   Vf1=(0.777+0.002062*t1+0.00001608*t1**2)/1000
   V1=Vf1
   PRINT*, "specific volume of saturated liquid=", Vf1
   hf1= 200+1.172*t1+0.001854*t1**2
   print*, "enthalpy of saturated liquid=", hf1
   uf1=0.0002367-(t1*1.715E-6)+(t1*t1*8.869E-9)
   print*, "viscosity of saturated liquid=", uf1
   d= (w*4)/(pie*dia**2)
   Vel1=d*Vf1
   print*, "velocity=", Vel1
   Re1=(Vel1*dia)/(Vf1*uf1)
   print*, "Reynolds Number=", Re1
   f1=0.33/(Re1)**0.25
   print*, "Friction Factor=", f1

   t=39
   z=2418.4/(t+273.15)
   y=15.06-z
   p=1000*E**y
   print*, "Pressure=", p/1000
   Vf=(0.777+0.002062*t+0.00001608*t*t)/1000
   print*, "Vf=", Vf
   Vg= (-4.26+94050*(t+273.15)/p)/1000
   print*, "Vg=", Vg
   hf= 200+1.172*t+0.001854*t**2
   print*, "hf=", hf
   hg=405.5+0.3636*t-0.002273*t**2
   print*, "hg=", hg
   uf=0.0002367-(t*1.715E-6)+(t*t*8.869E-9)
   print*, "uf=", uf
   ug= (11.945E-6)+(t*50.06E-9)+(t*t*0.2560E-9)
   print*, "ug=", ug
   d= (w*4)/(pie*dia**2)
   a=0.5*(Vg-Vf)*(Vg-Vf)*d**2
   print*, "a=", a
   b=1000*(hg-hf)+Vf*(Vg-Vf)*d**2
   print*, "b=", b
   c=1000*(hf-hf1)+(0.5*d*d*Vf**2)-(0.5*Vel1**2)
   print*, "c=", c
   x=(-b+ SQRT((b*b)-(4*a*c)))/(2*a)
   print*, "dryness fraction=", x
   h=hf*(1-x)+x*hg
   print*, "enthalpy=", h
   v=vf*(1-x)+vg*x
   print*, "specific volume=", v
   u=uf*(1-x)+x*ug
   print*, "viscosity=", u
   Vel=d*v
   print*, "velocity=", Vel
   Re=(Vel*d)/(V*u)
   print*, "Reynolds Number=", Re
   f=0.33/(Re)**0.25
   print*, "friction factor=", f
  
  
   stop
   end