#################################################### ## Matematično modeliranje okoljskih procesov ## ## 3. vaja: 2D resevanje enacbe transporta toka ## ## ## ## Izračun in izris grafov ## ## ## ## Izdelal: Lovrenc Pavlin ## ## Datum: 7.4.2015 ## #################################################### ## Podatki ## b <- 59 #sirina kanala Q <- 29 #pretok h <- 1 #globina q0 <- 0.29 #dotok polutanta c0 <- 100 #zacetna koncentracija K <- 0.001 #mesalni koeficient ## Izracun hitrosti toka ## u <- Q/(b*h) ## Izracun koncentracij #x <- c(50,100,500,1000,5000,10000,15000,20000) x <- seq(50,20000,10) #y <- c(0,1,2,3,4,5,10,15,20) y <- seq(-20,20,0.1) z <- matrix(0,nrow=length(x),ncol=length(y)) for (j in 1:length(y)){ for (i in 1:length(x)){ z[i,j] <- q0*c0/(2*h*u*sqrt(pi*x[i]/u*K))*exp(-u*y[j]**2/(4*K*x[i])) } } ############# ## Grafi ## ############# # Graf celotnega opazovanega obmocja png(filename="konc_celotno.png",width=1000,height=740,units="px") filled.contour(x,y,z, nlevels=40, color.palette=rainbow, plot.title = title(main = "Koncentracija polutanta v kanalu", xlab = "stacionaža [m]", ylab = "širina [m]"), key.title = title(main = "konc.\n[%]")) dev.off() # Graf omejenega območja x[50,200],y=[-5,5] png(filename="konc_delno.png",width=1000,height=740,units="px") filled.contour(x,y,z, nlevels=40, color.palette=rainbow, xlim=range(50:200,finite=TRUE), ylim=range(-5:5,finite=TRUE), plot.title = title(main = "Koncentracija polutanta v kanalu", xlab = "stacionaža [m]", ylab = "širina [m]"), key.title = title(main = "konc.\n[%]")) dev.off() # Graf po prečnih prerezih png(filename="konc_precni.png",width=1000,height=740,units="px") cl=rainbow(8) #paleta barv plot(y,z[1,],type="l", col=cl[1],lwd=2, main="Koncentracija po prečnih presekih", xlab="y [m]", ylab="c [%]") lines(y,z[4,],type="l", col=cl[2],lwd=2) lines(y,z[30,],type="l", col=cl[3],lwd=2) lines(y,z[100,],type="l", col=cl[4],lwd=2) lines(y,z[300,],type="l", col=cl[5],lwd=2) lines(y,z[500,],type="l", col=cl[6],lwd=2) lines(y,z[700,],type="l", col=cl[7],lwd=2) lines(y,z[900,],type="l", col=cl[8],lwd=2) legend("topleft", lwd=2, legend = c(paste("X =",toString(x[1])),paste("X=",toString(x[4])),paste("X=",toString(x[30])),paste("X=",toString(x[100])),paste("X=",toString(x[300])),paste("X=",toString(x[500])),paste("X=",toString(x[700])),paste("X=",toString(x[900]))), col=cl) dev.off() # Graf po vzdolžnih prerezih png(filename="konc_vzdolzi.png",width=1000,height=740,units="px") cl2=rainbow(6) #paleta barv Y=c(201,211,221,231,261,301) plot(x,z[,Y[1]],type="l", col=cl2[1],lwd=2, main="Koncentracija po vzdolžnih presekih", xlab="x [m]", ylab="c [%]") lines(x,z[,Y[2]],type="l", col=cl2[2],lwd=2) lines(x,z[,Y[3]],type="l", col=cl2[3],lwd=2) lines(x,z[,Y[4]],type="l", col=cl2[4],lwd=2) lines(x,z[,Y[5]],type="l", col=cl2[5],lwd=2) lines(x,z[,Y[6]],type="l", col=cl2[6],lwd=2) legend("topright", lwd=2, legend = c(paste("y =",toString(y[Y[1]])), paste("y =",toString(y[Y[2]])), paste("y =",toString(y[Y[3]])), paste("y =",toString(y[Y[4]])), paste("y =",toString(y[Y[5]])), paste("y =",toString(y[Y[6]]))), col=cl2) dev.off()