Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Finite Difference Method for non-square grid

  1. Sep 14, 2015 #1

    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.

  2. jcsd
  3. Sep 15, 2015 #2


    User Avatar
    Homework Helper

    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.
  4. Sep 15, 2015 #3
    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.

  5. Sep 15, 2015 #4
    Thanks for your suggestion.
  6. Sep 15, 2015 #5
    From the computer code he provided, it seems like this has already been done.

  7. Sep 15, 2015 #6
    By the way, what is the meaning of the parameter D0?

  8. Sep 15, 2015 #7
    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.

  9. Sep 15, 2015 #8
    D0 means double precision. Forexample, 2.D0.
  10. Sep 15, 2015 #9
    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?

  11. Sep 15, 2015 #10
    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.

    Last edited: Sep 15, 2015
  12. Sep 15, 2015 #11
    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:

    EDIT: Actually, this gives you exactly the same approximation you already have

    Rather 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.

    Last edited: Sep 15, 2015
  13. Sep 15, 2015 #12
    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.

  14. Sep 16, 2015 #13
    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.

  15. Sep 16, 2015 #14
    Thanks for your suggestions, specially for the curve boundary.

Share this great discussion with others via Reddit, Google+, Twitter, or Facebook