Multi-region Finite Difference- Interface between materials

[timman_24]

I am writing a multi-region diffusion code. The two regions have different material properties, so the analytical solution shows a discontinuity at the interface between the regions.

As can be seen here:

The numerical code I am running is (Mathematica):

While[converge > .00001,
 count = count + 1;
 list1[[1]] = list1[[2]];
 For[i = 2, i < Ai + 1, i++,
  list1[[i]] = Dw (list1[[i - 1]] + list1[[i + 1]])/(H^2 Ea1 + 2 Dw)];
 For[j = Ai + 1, j < (Ai + Bi), j++,
  list1[[j]] = .5*(-H^2 (list1[[j]] Ea2/Dc - S/Dc) + list1[[j - 1]] + 
      list1[[j + 1]])];
 converge = Max[Abs[list2 - list1]];
 list2 = list1;]

This works great if the material properties between the two regions are identical, but if I use differing material properties, I still get a smooth curve:

Is there a trick to getting this to work with finite difference method? How do I deal with this type of interface condition?

Thanks guys

 

[AiRAVATA]

Colud you please post the problem you're trying to solve? I'm not sure, but it looks like you're missing some sort of jump condition. You should check if your algorithm is correct when j = Ai, because that's where the coupling occurs.

 

[timman_24]

I ended up adding an interface condition to the solver.

To conserve current, we need this condition to be satisfied at the boundary:

$$D_{1}J_{1}(x)=D_{2}J_{2}(x)$$

To satisfy this condition, I stopped at the interface and applied this condition at that point. Then I continued. I used a backward difference scheme for the $$J_{1}$$ and a forward difference approx for the $$J_{2}$$ and rearranged for i.

This gave me the kink I was looking for.

BTW, this was a neutron diffusion model in a slab with two materials.