Commit: 14bab3c9776bff028afead1fb45e75b0fd83ac94 - myprojects-hg-reborn (hg)

Commit MetaInfo

• Revision: 14bab3c9776bff028afead1fb45e75b0fd83ac94 (tree)
• Zeit: 2007-01-30 20:37:41
• Autor: iselllo
• Commiter: iselllo

Log Message

This time a major update of the code R-codes/plot-boundary.R. Now the
code is able to read a velocity profile and a concentration profile and
work out the axial flux considering separately the contributions of the
layer and of the core. Careful about the input data files: temperature
and velocities may be defined on different domains and have a different
number of non-zero elements.

Ändern Zusammenfassung

• modified: R-codes/plot-boundary.R

Diff

diff -r 9ae8f6008c65 -r 14bab3c9776b R-codes/plot-boundary.R
--- a/R-codes/plot-boundary.R	Tue Jan 30 09:40:27 2007 +0000
+++ b/R-codes/plot-boundary.R	Tue Jan 30 11:37:41 2007 +0000
@@ -2,17 +2,30 @@
 solution<-read.table("boundary.txt",header=TRUE,nrow=20000)
 solution<-as.matrix(solution)
 delta<-abs(solution[300,1]-solution[299,1])
-#zero<-which(solution[ ,2]==0)
-#zero<-as.vector(zero)
-#zero<-zero[1]
-#end_sol<-zero-1
+
+
+# routine for integration
+def_integral<-function(integrand,h)
+{
+n<-length(integrand)
+myseq<-seq(1,n-4,4)
+integral<-0
+for (i in myseq)
+{
+integral<- integral + h/45.0*(14.0*integrand[i]+64.0*integrand[i+1] +24.0*integrand[i+2]+64.0*integrand[i+3]+14.0*integrand[i+4])
+
+}
+return(integral)
+}
+
+
 vecdim<-dim(solution)
 end_sol<-vecdim[1]
 solution<-solution[1:end_sol, ] # numerical solution up to the wall
                                # where it is still non-zero and it has to be corrected
 r_wall<-0.15895
 layer_thickness<-2.5e-3
-viscous_length<-2.5e-3/30 # either I first define the layer thickness and then I found out how many viscous lengths it amounts to, or I start from the viscous length times the number of lengths amounting to the layer thickness.
+viscous_length<-2.5e-3/150 # either I first define the layer thickness and then I found out how many viscous lengths it amounts to, or I start from the viscous length times the number of lengths amounting to the layer thickness.
 
 
 
@@ -30,7 +43,10 @@
 n_viscous<-as.integer(layer_thickness/viscous_length)
 n_delta<-as.integer(layer_thickness/delta)
 
-end_sol<-end_sol-n_delta  
+end_sol<-end_sol-n_delta-1 # I added a -1 for better compatibility wrt what I am doing below
+
+print("end_sol is")
+print(end_sol)
 
 solution<-solution[1:end_sol, ] # numerical solution up to the point where it is joined to 
             # the analytical expression
@@ -46,20 +62,27 @@
 
 wall_arr<-wall_fun()/norm_factor*solution[end_sol,2]
 
+
+
 pdf("wall_profile.pdf")  
 plot(y_plus,wall_arr,"l",lwd=2,col="blue",ylab="velocity profile",xlab="distance from wall in wall units")
 dev.off()
 
+
+
 y_plus<-y_plus*viscous_length
 y_plus<-sort(y_plus,decreasing=FALSE)+solution[end_sol,1]
 
 
+
+
 pdf("wall_profile2.pdf")  
 plot(y_plus ,wall_arr,"l",lwd=2,col="blue",ylab="velocity profile",xlab="distance from centre [m]")
 
 dev.off()
 
 
+
 extension<-cbind(y_plus,wall_arr)
 
 new_sol<-rbind(solution[1:end_sol-1, ],extension)
@@ -75,4 +98,119 @@
 fn <- paste("layer_thickness",layer_thickness,".txt",sep="")
 write.table(new_sol,fn,row.names=FALSE)
 
+
+#########################################################################
+#########################################################################
+#########################################################################
+#########################################################################
+#########################################################################
+#########################################################################
+
+
+
+
+# Now I try doing the same for the concentration. I have to use the same kind of correction 
+# close to the wall.
+
+solution2<-read.table("bound-conc.txt",header=TRUE,nrow=20000)
+solution2<-as.matrix(solution2)
+
+# now I re-define the solution
+
+
+#y_plus<-n_viscous
+
+#norm_factor<-wall_fun()
+
+#zero<-which(solution2[ ,2]==0)
+#zero<-as.vector(zero)
+#zero<-zero[1]                      #probably I can skip all this and keep 
+                      # the previous definition of end_sol (which I decreased by 1 anyway ) these complications arise since the velocity is defined up to the wall (though it it wrong in the layer) whereas the concentration is not defined (so it is zero) in the layer; therefore I have to pick up the right value [last non-zero value] for the concentration othrwise the union with the analytical prolongment gets screwed up
+#end_sol<-zero-1
+
+
+solution2<-solution2[1:end_sol, ] # numerical solution up to the wall
+                               # where it is still non-zero and it has to be corrected
+
+
+
+rm(y_plus)
+
+y_plus<-seq(n_viscous,0,length=15000)
+
+wall_arr<-wall_fun()/norm_factor*solution2[end_sol,2]
+
+
+
+y_plus<-y_plus*viscous_length
+y_plus<-sort(y_plus,decreasing=FALSE)+solution2[end_sol,1]
+
+
+extension2<-cbind(y_plus,wall_arr)
+
+# the following is the new solution for the concentration up to the wall
+new_sol2<-rbind(solution2[1:end_sol-1, ],extension2)
+
+pdf("new_solution_concentration.pdf")  
+plot(new_sol2[ ,1] ,new_sol2[ ,2],"l",lwd=2,col="blue",ylab="concentration profile",xlab="distance from centre [m]")
+dev.off()
+
+# now I want to start working out the integrals providing the incoming/outcoming flux along z. I divide the job: part outside the layer and inside the layer
+
+rho_v_no_layer<-solution[ ,2]*solution2[ ,2] 
+
+pdf("rho_v_outside_layer.pdf")  
+plot(solution[ ,1] ,rho_v_no_layer,"l",lwd=2,col="blue",ylab="rho*v",xlab="distance from centre [m]")
+dev.off()
+
+ # Now I have to multiply times 2*pi*r to take the radial integral
+
+rho_v_no_layer<-rho_v_no_layer*2*pi*solution[ ,1]
+
+pdf("rho_v_outside_layer*2pi*r.pdf")  
+plot(solution[ ,1] ,rho_v_no_layer,"l",lwd=2,col="blue",ylab="rho*v times 2pi*r",xlab="distance from centre [m]")
+dev.off()
+
+
+flux_no_layer<-def_integral(rho_v_no_layer,delta)
+print("the concentration flux outside the layer  is")
+print(flux_no_layer)
+
+
+# Now I need to work out the same thing with the layer.
+######################################
+######################################
+######################################
+######################################
+
+
+rho_v_layer<-extension[ ,2]*extension2[ ,2] 
+
+
+
+pdf("rho_v_layer.pdf")  
+plot(extension[ ,1] ,rho_v_layer,"l",lwd=2,col="blue",ylab="rho*v_layer",xlab="distance from centre [m]")
+dev.off()
+
+print("ok here1")
+
+ # Now I have to multiply times 2*pi*r to take the radial integral
+
+rho_v_layer<-rho_v_layer*2*pi*extension[ ,1]
+
+pdf("rho_v_outside_layer*2pi*r.pdf")  
+plot(extension[ ,1] ,rho_v_layer,"l",lwd=2,col="blue",ylab="rho*v_layer times 2pi*r",xlab="distance from centre [m]")
+dev.off()
+
+# Now I redefine the parameter delta I need for the integration
+
+delta<-abs(extension[200,1]-extension[201,1])
+
+flux_no_layer<-def_integral(rho_v_layer,delta)
+print("the concentration flux outside the layer  is")
+print(flux_no_layer)
+
+
+
+
 print("So far so good")
\ No newline at end of file