# # Cryo surgery # # Access ODE integrator library("deSolve"); # # Access functions for numerical solutions setwd("c:/R/bme_pde/chap8"); source("cryo_1.R"); source("kc.R"); source("Cp.R"); source("Qm.R"); source("dkc.R"); # # Model parameters rc=5; zc=10; Ts=37; Tp=-196; Qmt=0.0042; Cb=3.6; Cu=3.6; Cf=1.8; kf=0.02; ku=0.005; Ql=250; Tml=-8; Tmu=-1; # # Radial grid nr=6;dr=rc/(nr-1);drs=dr^2; r=seq(from=0,to=rc,by=dr); drs=dr^2; # # Axial grid nz=11;np=6; dz=zc/(nz-1);dzs=dz^2; z=seq(from=0,to=zc,by=dz); # # Grid in t nout=5;tf=120; tout=seq(from=0,to=tf,by=tf/(nout-1)); # # Display selected parameters cat(sprintf( "\n\n nr = %2d np = %2d nz = %2d\n",nr,np,nz)); # # Initial condition u0=rep(0,nr*nz);T0=Ts; for(i in 1:nr){ for(j in 1:nz){ u0[(i-1)*nz+j]=T0; } } ncall=0; # # ODE integration out=lsodes(y=u0,times=tout,func=cryo_1,parms=NULL); nrow(out) ncol(out) # # Calls to ODE routine cat(sprintf("\n\n ncall = %5d\n\n",ncall)); # # Plot z profiles T_zplot1=matrix(0,nrow=nz,ncol=nout); T_zplot2=matrix(0,nrow=nz,ncol=nout); T_zplot3=matrix(0,nrow=nz,ncol=nout); T_zplot4=matrix(0,nrow=nz,ncol=nout); T_zplot5=matrix(0,nrow=nz,ncol=nout); T_zplot6=matrix(0,nrow=nz,ncol=nout); for(it in 1:nout){ for(j in 1:nz){ T_zplot1[j,it]=out[it, j+1]; T_zplot2[j,it]=out[it, nz+j+1]; T_zplot3[j,it]=out[it,2*nz+j+1]; T_zplot4[j,it]=out[it,3*nz+j+1]; T_zplot5[j,it]=out[it,4*nz+j+1]; T_zplot6[j,it]=out[it,5*nz+j+1]; } T_zplot1[np,it]=Tp; T_zplot1[1,it] =Ts; T_zplot2[1,it] =Ts; T_zplot3[1,it] =Ts; T_zplot4[1,it] =Ts; T_zplot5[1,it] =Ts; T_zplot6[1,it] =Ts; } # # Automatic vertical scaling for each plot par(mfrow=c(3,2)); matplot(x=z,y=T_zplot1,type="l",xlab="z", ylab="T(r=0,z,t)",xlim=c(0,z[nz]),lty=1, main="T(r=0,z,t),t=0,30,...,120",lwd=2); matplot(x=z,y=T_zplot2,type="l",xlab="z", ylab="T(r=1,z,t)",xlim=c(0,z[nz]),lty=1, main="T(r=1,z,t),t=0,30,...,120",lwd=2); matplot(x=z,y=T_zplot3,type="l",xlab="z", ylab="T(r=2,z,t)",xlim=c(0,z[nz]),lty=1, main="T(r=2,z,t),t=0,30,...,120",lwd=2); matplot(x=z,y=T_zplot4,type="l",xlab="z", ylab="T(r=3,z,t)",xlim=c(0,z[nz]),lty=1, main="T(r=3,z,t),t=0,30,...,120",lwd=2); matplot(x=z,y=T_zplot5,type="l",xlab="z", ylab="T(r=4,z,t)",xlim=c(0,z[nz]),lty=1, main="T(r=4,z,t),t=0,30,...,120",lwd=2); matplot(x=z,y=T_zplot6,type="l",xlab="z", ylab="T(r=5,z,t)",xlim=c(0,z[nz]),lty=1, main="T(r=5,z,t),t=0,30,...,120",lwd=2); # # Vertical scale set as -200 <= T <= 40 par(mfrow=c(3,2)); matplot(x=z,y=T_zplot1,type="l",xlab="z", ylab="T(r=0,z,t)",xlim=c(0,z[nz]),lty=1, main="T(r=0,z,t),t=0,30,...,120",lwd=2, ylim=c(-200,40)); matplot(x=z,y=T_zplot2,type="l",xlab="z", ylab="T(r=1,z,t)",xlim=c(0,z[nz]),lty=1, main="T(r=1,z,t),t=0,30,...,120",lwd=2, ylim=c(-200,40)); matplot(x=z,y=T_zplot3,type="l",xlab="z", ylab="T(r=2,z,t)",xlim=c(0,z[nz]),lty=1, main="T(r=2,z,t),t=0,30,...,120",lwd=2, ylim=c(-200,40)); matplot(x=z,y=T_zplot4,type="l",xlab="z", ylab="T(r=3,z,t)",xlim=c(0,z[nz]),lty=1, main="T(r=3,z,t),t=0,30,...,120",lwd=2, ylim=c(-200,40)); matplot(x=z,y=T_zplot5,type="l",xlab="z", ylab="T(r=4,z,t)",xlim=c(0,z[nz]),lty=1, main="T(r=4,z,t),t=0,30,...,120",lwd=2, ylim=c(-200,40)); matplot(x=z,y=T_zplot6,type="l",xlab="z", ylab="T(r=5,z,t)",xlim=c(0,z[nz]),lty=1, main="T(r=5,z,t),t=0,30,...,120",lwd=2, ylim=c(-200,40));