Is it possible to use non-constant Neumann boundary conditions with the parabolic pde solver?

1 次查看（过去 30 天）
Link
Link 2012-8-30
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!

请先登录，再回答此问题。

采纳的回答

Bill Greene
Bill Greene 2012-8-31
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 个评论
Link
Link 2012-8-31
Thank you very much! That helps immensely in my understanding of how to use functions like pdebound.
If I had a second Neumann boundary condition defined on the same body, would I just concatenate it onto g? And similarly, a zero flux boundary condition at the center would be of the form g = 0 for all xyMidEdge points?
Bill Greene
Bill Greene 2012-9-2
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

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Eigenvalue Problems 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by