Is it possible to use non-constant Neumann boundary conditions with the parabolic pde solver?
4 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
I am looking to solve the 2D heat equation T = T(r, theta, t) on a circle. I have a non-constant Neumann BC on the outside of the circle (r = a), where I have a heat flux as a function of theta (-k*dT/dx = q(theta) at r = a).
I know that I can decompose theta in to atan(y/x) -- in general though, I am unclear on how to use either the pdebound or assemb functions to prescribe these boundary conditions.
Thank you for your help!
0 comentarios
Respuesta aceptada
Bill Greene
el 31 de Ag. de 2012
Hi,
You are definitely on the right track in assuming that creating a pdebound function is a good way to define this BC. I created a simple example below that assumes you have an inward heat flux of 5*sin(theta), defines a pdebound function for this BC (I called it boundaryConditions), and then uses that to solve the heat equation on a circle. Hopefully this example will get you past some of the sticky issues.
Regards,
Bill
function transientHeatCircle( )
radius = 2;
gdm=[1 0 0 radius]';
g = decsg(gdm, 'C1', ('C1')');
[p,e,t]=initmesh(g);
c = 1; d = 2; a = 0; f = 0;
b = @boundaryConditions;
u0=0; tlist=0:.001:1;
u=parabolic(u0, tlist, b,p,e,t,c,a,f,d);
figure; pdeplot(p, e, t, 'xydata', u(:,end), 'mesh', 'on'); axis equal;
title 'Final Temperature Distribution'
n=4; %grid point at theta=pi/2
figure; plot(tlist, u(n,:)); grid on;
title(sprintf('Temperature at (%3.1f,%3.1f) as a Function of Time', ...
p(1,n), p(2,n)));
end
function [ q, g, h, r ] = boundaryConditions( p, e, u, time )
N = 1;
ne = size(e,2);
q = zeros(N^2, ne);
% calculate coordinates of edge mid-points
xy1 = p(:,e(1,:));
xy2 = p(:,e(2,:));
xyMidEdge = (xy1+xy2)/2;
g = 5*sin(atan2(xyMidEdge(2,:),xyMidEdge(1,:)));
h = zeros(N^2, 2*ne);
r = zeros(N, 2*ne);
end
2 comentarios
Bill Greene
el 2 de Sept. de 2012
My code handled the simple case where you have a single BC expression that applies to all edges on the boundary. A more general version of a boundary function might look something like this.
function [ q, g, h, r ] = boundFunc( p, e, u, time )
N = 1; % number of PDEs in the system
ne = size(e,2);
q = zeros(N^2, ne);
g = zeros(N, ne);
h = zeros(N^2, 2*ne);
r = zeros(N, 2*ne);
for i=1:ne
ei = e(5,i);
if(ei == 1)
% apply appropriate BCs on edge 1
else if(ei == 3)
% apply appropriate BCs on edge 3
else if(ei == ...)
% apply appropriate BCs on this edge
end
end
end
Más respuestas (0)
Ver también
Categorías
Más información sobre Eigenvalue Problems en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!