neumann boundary condition in diffusion on a disk and ring

I want to solve the diffusion equation on a disk. The disk has two regions. One region is a disk with a radius of 1 and the other region is a ring around the disk. You can see the geometry below:

I want to apply neumann condition on the edges E1 to E9 and set the flux zero at these edges. To set the flux zero I use the defualt neumann condition:

applyBoundaryCondition(model,'neumann','Edge',[1,2,3,4,5,6,7,8,9]);

The initial condition at time zero looks like this:

function uinit = u0fun(location)
M = length(location.x);
uinit = zeros(1,M);
uinit(1,:) = 100.*exp(-((location.x).^2+location.y.^2)/0.1);
end

here a plot for it:

The solution at time 100 looks below and has the same value everywhere (look at the color map):

I expected becuase of the boundary condition there will not be any flux to the ring around the circle. However, as the solution indicates there is indeed flux to the ring.

Where does the problem lie?

Below is the code I used to solve the diffusion equation:

clear all
close all
clc

n_eqs = 1;
model = createpde(n_eqs);
% First number should be 1 for a circle. The second and thrid one belongs
% to the center of the circle and the last one is the radius
c1 = [1;0;0;1];
c2 = [1;0;0;1.2];

gd = [c1,c2];
ns = char('c1','c2');
ns = ns';
sf = 'c1 + (c2-c1)';
[dl,bt] = decsg(gd,sf,ns);

pdegplot(dl,'EdgeLabels','on','FaceLabels','on')

geometryFromEdges(model,dl);

applyBoundaryCondition(model,'neumann','Edge',[1,2,3,4,5,6,7,8,9]);

specifyCoefficients(model,'m',0,...
                          'd',1,...
                          'c',1,...
                          'a',0,...
                          'f',0);
                      
mesh_default = generateMesh(model,'Hmax',0.1);    

% figure
% pdemesh(mesh_default)


%Initial condition
u0 = @u0fun;
setInitialConditions(model,u0);

tlist = linspace(0,10,100);
results = solvepde(model,tlist);
res = results.NodalSolution;


%Initial condition
figure
pdeplot(model,'XYData',results.NodalSolution(:,1))

%Solution at time 100
figure
pdeplot(model,'XYData',results.NodalSolution(:,100))
                   

[Torsten]

My guess is that the boundary condition is ignored at edges 1,2,3 and 4 because they are in the interior of the domain.

Answer 1

[Alan Weiss]

I think that you might be making a modeling error. If the flux across the boundaries is zero, then perhaps you should not have the inner disk be part of the problem. Instead, make the problem be just for the outer ring. If, at the end, you want to extend the solution to the inner disk, then use the solution at the inner edge of the ring be a boundary condition for the solution of the disk.

As Torsten said, PDE Toolbox allows no boundary conditions between subdomains.

Alan Weiss

MATLAB mathematical toolbox documentation

[Minxiang Ji]

I want to ask a question. If I have a boundary condition on the interior boundarys such as

n is the outer normal vector on interior boundary. So how can I use MATLAB to solve such a problem? Can I set a initial geuss of the boundary and alternately calcutate on the two subdomains until the convergence?