Multi-region Finite Difference- Interface between materials

Click For Summary
The discussion focuses on developing a multi-region diffusion code in Mathematica, specifically addressing the challenge of modeling an interface between two materials with differing properties. The original numerical code produces smooth curves even when material properties differ, indicating a potential issue with handling the interface condition. Participants suggest that a jump condition is necessary at the interface to ensure current conservation, expressed as D_{1}J_{1}(x)=D_{2}J_{2}(x). Implementing this condition at the boundary and using different difference schemes for each region can create the desired discontinuity. The problem is contextualized within a neutron diffusion model in a slab with two materials.
timman_24
Messages
52
Reaction score
0
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:

attachment.php?attachmentid=35045&stc=1&d=1304197454.png


The numerical code I am running is (Mathematica):

Code:
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:

attachment.php?attachmentid=35046&stc=1&d=1304197454.png


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
 

Attachments

  • Analyticalpart1.png
    Analyticalpart1.png
    2.3 KB · Views: 549
  • Analyticalpart1num.png
    Analyticalpart1num.png
    2.3 KB · Views: 616
Physics news on Phys.org
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.
 
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.
 

Similar threads

  • · Replies 6 ·
Replies
6
Views
4K
  • · Replies 13 ·
Replies
13
Views
3K
  • · Replies 3 ·
Replies
3
Views
2K
  • · Replies 7 ·
Replies
7
Views
5K
  • · Replies 5 ·
Replies
5
Views
2K
  • · Replies 16 ·
Replies
16
Views
1K
  • · Replies 4 ·
Replies
4
Views
1K
  • · Replies 3 ·
Replies
3
Views
3K
  • · Replies 7 ·
Replies
7
Views
2K
  • · Replies 1 ·
Replies
1
Views
2K