neumann boundary condition in diffusion on a disk and ring

10 views (last 30 days)
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))
  1 Comment
Torsten
Torsten on 22 Mar 2019
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.

Sign in to comment.

Answers (1)

Alan Weiss
Alan Weiss on 22 Mar 2019
Edited: Alan Weiss on 22 Mar 2019
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
  1 Comment
Minxiang Ji
Minxiang Ji on 6 Jan 2021
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?

Sign in to comment.

Tags

Products


Release

R2018b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!