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

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!

Answer 1

[Bill Greene]

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

[Link]

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]

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