## Numerical solution of PDE:s, Part 3: 2D diffusion problem

In the earlier posts related to PDE numerical integration, I showed how to discretize 1-dimensional diffusion or heat conduction equations either by explicit or implicit methods. The 1d model can work well in some situations where the symmetry of a physical system makes the concentration or temperature field practically depends on only one cartesian coordinate and is independent of the position along two other orthogonal coordinate axes.

The 2-dimensional version of the diffusion/heat equation is

assuming that the diffusion is isotropic, i.e. its rate does not depend on direction.

In a discretized description, the function C(x,y,t) would be replaced by a three-index object

assuming that one of the corners of the domain is at the origin. Practically the same can also be done with two indices, as in

where $N_x$ is the number of discrete points in x-direction.

Using the three-index version of the 2D diffusion equation and doing an explicit discretization, we get this kind of a discrete difference equation:

which can be simplified a bit if we have $\Delta x = \Delta y$ .

Below, I have written a sample R-Code program that calculates the evolution of a Gaussian concentration distribution by diffusion, in a domain where the x-interval is 0 < x < 6 and y-interval is 0 < y < 6. The boundary condition is that nothing diffuses through the boundaries of the domain, i.e. the two-dimensional integral of C(x,y,t) over the area $[0,6]\times [0,6]$ does not depend on t.

library(graphics) #load the graphics library needed for plotting

lx <- 6.0 #length of the computational domain in x-direction
ly <- 6.0 #length of the computational domain in y-direction
lt <- 6 #length of the simulation time interval
nx <- 30 #number of discrete lattice points in x-direction
ny <- 30 #number of discrete lattice points in y-direction
nt <- 600 #number of timesteps
dx <- lx/nx #length of one discrete lattice cell in x-direction
dy <- ly/ny #length of one discrete lattice cell in y-direction
dt <- lt/nt #length of timestep

D <- 1.0 #diffusion constant (assumed isotropic)

Conc2d = matrix(nrow=ny,ncol=nx)
DConc2d <- matrix(nrow=ny, ncol=nx) #a vector for the changes in concentration during a timestep
xaxis <- c(0:(nx-1))*dx #the x values corresponding to the discrete lattice points
yaxis <- c(0:(ny-1))*dy #the y values corresponding to the discrete lattice points

kappax <- D*dt/(dx*dx) #a parameter needed in the discretization
kappay <- D*dt/(dy*dy) #a parameter needed in the discretization

for (i in c(1:ny)) {
for (j in c(1:nx)) {
Conc2d[i,j] = exp(-(i*dy-3)*(i*dy-3)-(j*dx-3)*(j*dx-3)) #2D Gaussian initial concentration distribution
DConc2d[i,j] <- 0 #all initial values in DConc vector zeroed
}
}

for (j in c(1:nt)) { #main time stepping loop

for(k in c(1:nx))
{
Conc2d[1,k] <- Conc2d[k,2] #fluxes through the boundaries of the domain are forced to stay zero
Conc2d[nx,k] <- Conc2d[nx-1,k]
}

for(k in c(1:ny)) {
Conc2d[k,1] <- Conc2d[k,2]
Conc2d[k,nx] <- Conc2d[k,nx-1]
}
for (k in c(2:(ny-1))) {
for (l in c(2:(nx-1))) {
DConc2d[k,l] <- kappax*(Conc2d[k,l-1]-2*Conc2d[k,l]+Conc2d[k,l+1]) + kappay*(Conc2d[k-1,l]-2*Conc2d[k,l]+Conc2d[k+1,l]) #time stepping
}
}

for (k in c(2:(ny-1))) {
for (l in c(2:(nx-1))) {
Conc2d[k,l] <- Conc2d[k,l]+DConc2d[k,l] #add the changes to the vector Conc
}
}
k <- 0
l <- 0

if(j %% 3 == 1) { #make plots of C(x,y) on every third timestep
jpeg(file = paste("plot_",j,".jpg",sep=""))
persp(yaxis,xaxis,Conc2d,zlim=c(0,1))
title(paste("C(x,y) at t =",j*dt))
dev.off()
}
}

To plot the distribution C(x,y) with a color map instead of a 3D surface, you can change the line

persp(yaxis,xaxis,Conc2d,zlim=c(0,1))

to

image(yaxis,xaxis,Conc2d,zlim=c(0,0.3))

The two kinds of graphs are shown below for three different values of t.

Figure 1. Surface plots of the time development of a Gaussian mass or temperature distribution spreading by diffusion.

Figure 2. Time development of a Gaussian mass or temperature distribution spreading by diffusion, plotted with a red-orange-yellow color map.

To obtain an implicit differencing scheme, we need to write the discretized equation as a matrix-vector equation where $C_{i;j;k}$ is obtained from $C_{i;j;k+1}$ by backward time stepping. In the square matrix, there is one row (column) for every lattice point of the discrete coordinate system. Therefore, if there are N points in x-direction and N points in y-direction, then the matrix is an $N^2 \times N^2$ – matrix. From this it’s quite obvious that the computation time increases very quickly when the spatial resolution is increased.

Initially, one may think that the equation corresponding to diffusion in a really coarse $3 \times 3$ grid is the following one:

Where I’ve used the denotation $k_x = \frac{D\Delta t}{(\Delta x)^2}$ and $k_y = \frac{D\Delta t}{(\Delta x)^2}$ The problem with this is that there are unnecessary terms $-k_x$ in here, creating a cyclic boundary condition (mass that is diffusing through the boundary described by line x = L reappears from the boundary on the other side, x = 0. The correct algorithm for assigning the nonzero elements $A_{ij}$ of the matrix A is

1. $A_{ij} = 1 + 2k_x + 2k_y$ , when $i = j$
2. $A_{ij} = -k_x$ , when $j=i-1$ AND $i\neq 1$ (modulo $N_x$)
3. $A_{ij} = -k_x$ , when $j=i+1$ AND $i\neq 0$ (modulo $N_x$)
4. $A_{ij} = -k_y$ , when $j=i+N_x$ OR $j=i-N_x$

when you want to get a boundary condition which ensures that anything that diffuses through the boundaries is lost forever. Then the matrix-vector equation in the case of $3 \times 3$ lattice is

Unlike the linear system in the 1D diffusion time stepping, this is not a tridiagonal problem and consequently is slower to solve. An R-Code that produces a series of images of the diffusion process for a Gaussian concentration distribution in a $15 \times 15$ discrete lattice is given below.

library(graphics) #load the graphics library needed for plotting

lx <- 6.0 #length of the computational domain in x-direction
ly <- 6.0 #length of the computational domain in y-direction
lt <- 6 #length of the simulation time interval
nx <- 15 #number of discrete lattice points in x-direction
ny <- 15 #number of discrete lattice points in y-direction
nt <- 180 #number of timesteps
dx <- lx/nx #length of one discrete lattice cell in x-direction
dy <- ly/ny #length of one discrete lattice cell in y-direction
dt <- lt/nt #length of timestep

D <- 1.0 #diffusion constant

C = c(1:(nx*ny))
Cu = c(1:(nx*ny))
Conc2d = matrix(nrow=ny,ncol=nx)
xaxis <- c(0:(nx-1))*dx #the x values corresponding to the discrete lattice points
yaxis <- c(0:(ny-1))*dy #the y values corresponding to the discrete lattice points

kappax <- D*dt/(dx*dx) #a parameter needed in the discretization
kappay <- D*dt/(dy*dy) #a parameter needed in the discretization

A = matrix(nrow=(nx*ny),ncol=(nx*ny))

for(i in c(1:(nx*ny))) {
for(j in c(1:(nx*ny))) {
A[i,j] <- 0
if(i==j) A[i,j] <- 1+2*kappax+2*kappay
if(j==i+1 && (i%%nx != 0)) A[i,j] <- -kappax
if(j==i-1 && (i%%nx != 1)) A[i,j] <- -kappax
if(j==i+nx) A[i,j] <- -kappay
if(j==i-nx) A[i,j] <- -kappay
}
}

for (i in c(1:ny)) {
for (j in c(1:nx)) {
Conc2d[i,j] <- exp(-2*(i*dy-3)*(i*dy-3)-2*(j*dx-3)*(j*dx-3)) #2D Gaussian initial concentration distribution
}
}

for (j in c(1:nt)) { #main time stepping loop

for(k in c(1:ny)) {
for(l in c(1:nx)) {
C[(k-1)*nx+l] <- Conc2d[k,l]
}
}

Cu <- solve(A,C)

for(k in c(1:ny)) {
for(l in c(1:nx)) {
Conc2d[k,l] <- Cu[(k-1)*nx+l]
}
}

k <- 0
l <- 0

if(j %% 3 == 1) { #make plots of C(x,y) on every third timestep
jpeg(file = paste("plot_",j,".jpg",sep=""))
persp(y=xaxis,x=yaxis,z=Conc2d,zlim=c(0,1))
title(paste("C(x) at t =",j*dt))
dev.off()
}
}

Note that the function C(x,y) approaches a constant zero function as t increases, try modifying the code yourself to make a boundary condition that doesn’t let anything diffuse through the boundaries!