Finite Difference Method for non-square grid

[nazmulislam]

Hi,

I have written some codes for the finite difference solution of diffusion equation (\frac{\partial c}{\partial t}= D {\nabla^2 c}, where c is the species concentration and D is the diffusion coefficient) as follows:

DO k= 1, tsteps+1
DO i = 2, zsteps
DO j = 2, rsteps
cnew(i,j) = c(i,j) + D*dt*((c(i+1,j) - 2.D0*c(i,j) + c(i-1,j))/dz**2 + (c(i,j+1)-2.D0*c(i,j) + c(i,j-1))/dr**2 +
(c(i,j+1) - c(i,j-1))/(r(j)*2.D0*dr))
END DO
END DO
END DO

where 'dt' is the step size in time, 'dz' is the z-step size, 'dr' is the r-step size, 'tsteps' is the number of time steps, 'zsteps' is the number of steps in the z-direction and 'rsteps' is the number of steps in the r-direction.

The code is working only for the square grid (dz=dr), but not working for the non-square grid space.

Can anybody help me in this matter.

Cheers

 

[RUber]

Normally, finite difference codes require square grids to work. If your grids are spaced in non-square segments, you will need to adjust your approximations for the derivatives to adjust for the transformation of dy or dx to the new directions.

 

[Chestermiller]

Hi,

I have written some codes for the finite difference solution of diffusion equation (\frac{\partial c}{\partial t}= D {\nabla^2 c}, where c is the species concentration and D is the diffusion coefficient) as follows:

DO k= 1, tsteps+1
DO i = 2, zsteps
DO j = 2, rsteps
cnew(i,j) = c(i,j) + D*dt*((c(i+1,j) - 2.D0*c(i,j) + c(i-1,j))/dz**2 + (c(i,j+1)-2.D0*c(i,j) + c(i,j-1))/dr**2 +
(c(i,j+1) - c(i,j-1))/(r(j)*2.D0*dr))
END DO
END DO
END DO

where 'dt' is the step size in time, 'dz' is the z-step size, 'dr' is the r-step size, 'tsteps' is the number of time steps, 'zsteps' is the number of steps in the z-direction and 'rsteps' is the number of steps in the r-direction.

The code is working only for the square grid (dz=dr), but not working for the non-square grid space.

Can anybody help me in this matter.

Cheers

There could be two sources of your difficulty.

1. How are you handling the finite differencing for the boundary condition at r = 0.
2. You are using an explicit finite difference scheme in time, so there may be an issue with numerical stability.

Chet

 

[nazmulislam]

Normally, finite difference codes require square grids to work. If your grids are spaced in non-square segments, you will need to adjust your approximations for the derivatives to adjust for the transformation of dy or dx to the new directions.

Thanks for your suggestion.

 

[Chestermiller]

Normally, finite difference codes require square grids to work. If your grids are spaced in non-square segments, you will need to adjust your approximations for the derivatives to adjust for the transformation of dy or dx to the new directions.

From the computer code he provided, it seems like this has already been done.

Chet

 

[Chestermiller]

By the way, what is the meaning of the parameter D0?

Chet

 

[nazmulislam]

There could be two sources of your difficulty.

1. How are you handling the finite differencing for the boundary condition at r = 0.
2. You are using an explicit finite difference scheme in time, so there may be an issue with numerical stability.

Chet

Many thanks for your response. I have imposed symmetry condition, cnew(i,1)=cnew(i,2) at r=0 (since I am calculating concentration within a straight tube). Also, I have used no-flux condition, cnew(i, rsteps+1)=cnew(i, rsteps) at r=r_max (i.e., on the tube surface), and periodic condition cnew(1, j)=(cnew(2, j)+cnew(zsteps, j))/2.D0 and cnew(zsteps+1,j)=cnew(1, j).

I have set the step size in time as dt=(dz**2+dr**2)/(8.D0*d) which is the stability condition for two spatial dimensions, where 'd' is the diffusion coefficient.

Thanks.

 

[nazmulislam]

By the way, what is the meaning of the parameter D0?

Chet

D0 means double precision. Forexample, 2.D0.

 

[Chestermiller]

I'm having trouble following what you did. Please first write out the spatial finite difference approximations you used (algebraically), including what you did at r = 0. Also, what is this about a periodic boundary condition?

Chet

 

[nazmulislam]

I'm having trouble following what you did. Please first write out the spatial finite difference approximations you used (algebraically), including what you did at r = 0. Also, what is this about a periodic boundary condition?

Chet

I have used the following central difference formula for the spatial derivatives: (written in Latex)

$$\frac{\partial^2 c}{\partial z^2} = \frac{c_{\rm i+1,j}^{k}-2 c_{\rm i,j}^{k}+c_{\rm i-1,j}^{k}}{(\Delta z)^2}+O[(\Delta z)^2]$$

$$\frac{\partial c}{\partial r} = \frac{c_{\rm i,j+1}^{k}-c_{\rm i,j-1}^{k}}{(\Delta z)^2}+O[(\Delta r)^2]$$

$$\frac{\partial^2 c}{\partial r^2}=\frac{c_{\rm i,j+1}^{k}-2 c_{\rm i,j}^{k}+c_{\rm i,j-1}^{k}}{(\Delta r)^2}+O[(\Delta r)^2]$$

and the forward time difference of the time derivative

$$\frac{\partial c}{\partial t}=\frac{c_{\rm i,j}^{k+1}-c_{\rm i,j}^{k}}{\Delta t}+O[(\Delta t)].$$

I have imposed the symmetry condition at r=0 as follows:

$$\frac{\partial c(z,r,t)}{\partial z}=0$$ at r=0.

Also, I have used the following periodic conditions:

$$c(z,r,t)|_{z=-L/2} = c(z,r,t)|_{z=L/2}$$,

$$\frac{\partial c(z,r,t)}{\partial z}|_{z=-L/2} = \frac{\partial c(z,r,t)}{\partial z}|_{z=L/2}$$

where 'L' is the length of the tube. Periodic condition indicate that the concentration is the same at two ends-left & right.

At first, i have calculated the concentrations at all interior points by using finite difference approximations that I wrote in the first message. Then, I update the concentrations at boundary by the boundary conditions.

Thanks.

 

[Chestermiller]

I have used the following central difference formula for the spatial derivatives: (written in Latex)

$$\frac{\partial^2 c}{\partial z^2} = \frac{c_{\rm i+1,j}^{k}-2 c_{\rm i,j}^{k}+c_{\rm i-1,j}^{k}}{(\Delta z)^2}+O[(\Delta z)^2]$$

$$\frac{\partial c}{\partial r} = \frac{c_{\rm i,j+1}^{k}-c_{\rm i,j-1}^{k}}{(\Delta z)^2}+O[(\Delta r)^2]\tag{1}$$

$$\frac{\partial^2 c}{\partial r^2}=\frac{c_{\rm i,j+1}^{k}-2 c_{\rm i,j}^{k}+c_{\rm i,j-1}^{k}}{(\Delta r)^2}+O[(\Delta r)^2]$$

and the forward time difference of the time derivative

$$\frac{\partial c}{\partial t}=\frac{c_{\rm i,j}^{k+1}-c_{\rm i,j}^{k}}{\Delta t}+O[(\Delta t)].$$

I have imposed the symmetry condition at r=0 as follows:

$$\frac{\partial c(z,r,t)}{\partial z}=0\tag{2}$$ at r=0.

Also, I have used the following periodic conditions:

c(z,r,t)|_{z=-L/2} = c(z,r,t)|_{z=L/2},

$$\frac{\partial c(z,r,t)}{\partial z}\right|_{z=-L/2} = \frac{\partial c(z,r,t)}{\partial z}\right|_{z=L/2}$$

where 'L' is the length of the tube. Periodic condition indicate that the concentration is the same at two ends-left & right.

At first, i have calculated the concentrations at all interior points by using finite difference approximations that I wrote in the first message. Then, I update the concentrations at boundary by the boundary conditions.

Thanks.

You left out the delimiters on the LaTex expressions so I tried editing your post. I was able to edit most of the equations, but not all of them.

I have big issues with how you worked this problem numerically. I don't even know where to start.

First let me correct some of your equations above. Eqns. labeled 1 and 2 above are both incorrect. Eqn. 1 should have a $2Δr$ in the denominator, and Eqn. 2 should have a partial with respect to r in the denominator.

Now for your finite difference equations. You are trying to represent $\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial c}{\partial r}\right)$ by a second order accurate finite difference approximation. The relationship you used is fine for all the points except at the grid point at r = 0. However, a better approximation is:

$$\frac{r_{(j+1/2)}(c_{j+1}-c_j)-r_{(j-1/2)}(c_j-c_{j-1})}{r_j(Δr)^2}$$
EDIT: Actually, this gives you exactly the same approximation you already haveRather than solving for the concentrations in the interior region only and just applying a first-order extension at r = 0 and r = outer radius, you would do much better if you carried out the time-dependent integration at these locations as well.

At r = 0, the first derivative of concentration with respect to r is zero, so in the vicinity at r = 0, the concentration must be approximately quadratic with r, according to $c = c_0+(c_1-c_0)\frac{r^2}{(Δr)^2}$ to second order accuracy. This leads to the result that:
$$\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial c}{\partial r}\right)=\frac{4(c_1-c_0)}{(Δr)^2}$$
at r = 0, to second order accuracy. See if you can verify this result.

With this result, you can apply the time dependent integration not only at the interior points, but also at r = 0.

At the outer radius, you know that $c_{n+1}=c_{n-1}$. You can substitute this result into your finite difference equation and thereby include the point at the outer radius in the time dependent integration.

It looks like you have set this problem up to solve for the concentration distribution as a function of time, based on some pre-specified concentration distribution at time = 0. Is this correct?

Also, the necessity for periodic boundary conditions on z is a little confusing. Please tell me what the purpose of this is, and please tell me what the initial concentration distribution looks like.

Chet

 

[nazmulislam]

You left out the delimiters on the LaTex expressions so I tried editing your post. I was able to edit most of the equations, but not all of them.

I have big issues with how you worked this problem numerically. I don't even know where to start.

First let me correct some of your equations above. Eqns. labeled 1 and 2 above are both incorrect. Eqn. 1 should have a $2Δr$ in the denominator, and Eqn. 2 should have a partial with respect to r in the denominator.

Now for your finite difference equations. You are trying to represent $\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial c}{\partial r}\right)$ by a second order accurate finite difference approximation. The relationship you used is fine for all the points except at the grid point at r = 0. However, a better approximation is:

$$\frac{r_{(j+1/2)}(c_{j+1}-c_j)-r_{(j-1/2)}(c_j-c_{j-1})}{r_j(Δr)^2}$$
EDIT: Actually, this gives you exactly the same approximation you already haveRather than solving for the concentrations in the interior region only and just applying a first-order extension at r = 0 and r = outer radius, you would do much better if you carried out the time-dependent integration at these locations as well.

At r = 0, the first derivative of concentration with respect to r is zero, so in the vicinity at r = 0, the concentration must be approximately quadratic with r, according to $c = c_0+(c_1-c_0)\frac{r^2}{(Δr)^2}$ to second order accuracy. This leads to the result that:
$$\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial c}{\partial r}\right)=\frac{4(c_1-c_0)}{(Δr)^2}$$
at r = 0, to second order accuracy. See if you can verify this result.

With this result, you can apply the time dependent integration not only at the interior points, but also at r = 0.

At the outer radius, you know that $c_{n+1}=c_{n-1}$. You can substitute this result into your finite difference equation and thereby include the point at the outer radius in the time dependent integration.

It looks like you have set this problem up to solve for the concentration distribution as a function of time, based on some pre-specified concentration distribution at time = 0. Is this correct?

Also, the necessity for periodic boundary conditions on z is a little confusing. Please tell me what the purpose of this is, and please tell me what the initial concentration distribution looks like.

Chet

Sorry for my typing mistake in writing equations (1) & (2). Thanks for the correction. Also, I tried to write equation in Latex with single $..$ instead of writing with double $. For that reason it didn't work. Thanks for adding the delimiters. I have deleted '\right' from the last equation that I wrote previously and now it works and you can see it.

I have set uniform concentration, c=1, along centerline (at r=0) with some radius say, $r_0$ at t=0. This is my initial condition.

Although, I am writing the FDM codes for straight tube now, I will need to change the codes for some periodic tube. For a periodic tube, I assume the flux and the concentration will be the same at two ends of one period. For that reason, I have used the above mentioned periodic conditions.

Thanks.

 

[Chestermiller]

Sorry for my typing mistake in writing equations (1) & (2). Thanks for the correction. Also, I tried to write equation in Latex with single $..$ instead of writing with double $. For that reason it didn't work. Thanks for adding the delimiters. I have deleted '\right' from the last equation that I wrote previously and now it works and you can see it.

I have set uniform concentration, c=1, along centerline (at r=0) with some radius say, $r_0$ at t=0. This is my initial condition.

Although, I am writing the FDM codes for straight tube now, I will need to change the codes for some periodic tube. For a periodic tube, I assume the flux and the concentration will be the same at two ends of one period. For that reason, I have used the above mentioned periodic conditions.

Thanks.

For a periodic tube, there must be some symmetry where you can shift the z origin so that the dc/dz is zero at the two ends. This would be much easier to handle than actually applying a periodic boundary condition.

When you switch to periodic variation of the diameter with respect to z, I strongly recommend that you transform the radial coordinate to r*=r/R(z), where R is the tube radius and r* is the new radial coordinate. This will remap the problem onto a rectangular region again, and will put all the complexity into the differential equation where it can be handled much more easily (than a curved boundary). Of course, you will need to do the coordinate transformation carefully in the differential equation.

Your current stability problem may be related to the use of a periodic boundary condition. This may change the stability characteristics of the finite difference scheme. If you continue requiring a periodic boundary condition, I strongly recommend switching to a Backward Euler integration scheme in time, such that all the spatial differences on the right hand side of the equation are evaluated at k+1, rather than k. This will get rid of any instability problem, but will require you to solve a set of linear algebraic equations (for all the c's at k=1) at each time step.

Chet

 

[nazmulislam]

For a periodic tube, there must be some symmetry where you can shift the z origin so that the dc/dz is zero at the two ends. This would be much easier to handle than actually applying a periodic boundary condition.

When you switch to periodic variation of the diameter with respect to z, I strongly recommend that you transform the radial coordinate to r*=r/R(z), where R is the tube radius and r* is the new radial coordinate. This will remap the problem onto a rectangular region again, and will put all the complexity into the differential equation where it can be handled much more easily (than a curved boundary). Of course, you will need to do the coordinate transformation carefully in the differential equation.

Your current stability problem may be related to the use of a periodic boundary condition. This may change the stability characteristics of the finite difference scheme. If you continue requiring a periodic boundary condition, I strongly recommend switching to a Backward Euler integration scheme in time, such that all the spatial differences on the right hand side of the equation are evaluated at k+1, rather than k. This will get rid of any instability problem, but will require you to solve a set of linear algebraic equations (for all the c's at k=1) at each time step.

Chet

Thanks for your suggestions, specially for the curve boundary.

Cheers