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

assembparabolicpdepdebound

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!

Best Answer

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

Related Solutions

MATLAB: Time Dependent Boundary Conditions in PDE Toolbox

Michael,

Before commenting on your f-function, let me make sure I understand the situation. This issue is independent of your original post, right? I assume you have resolved the original problem?

As far as your f_coeff function, I don't know of any reason why PDE Toolbox would necessarily have problems with coefficients that are discontinuous with respect to time. What problem are you seeing? One problem I see immediately is that your original geometry (your script above) has only one subdomain. Your f_coeff function is setting f(j) to be non-zero only for a non-existent subdomain 2.

Another issue concerns how PDE Toolbox detects whether a coefficient is time dependent or not. As a test, the code passes time=NaN and expects the function to return a result that contains at least one NaN. This is documented here: http://www.mathworks.com/help/pde/ug/pde-coefficients.html but is easy to overlook.

I suggest adding code like this to the top of your f_coeff function:

if(isnan(time))
f = NaN*ones(1,size(t,2));
end

Bill

MATLAB: Supporting documentation regarding PDE boundary files

Hi,

Your question is certainly a reasonable one. The boundary matrix format in those files is unquestionably difficult to decode.

I suggest an alternate approach instead of trying to decipher those example files. This example boundary file below will define a zero Neumann boundary condition on all edges of a model-- a square or something much more complicated.

    function [ q, g, h, r ] = boundaryFileZeroNeumann( p, e, u, time )
        %   Define a zero Neumann condition on all edges
        N = 1 % number of pde in the system
        ne = size(e,2); % number of boundary edges
        q = zeros(N^2, ne); % Neumann coefficient q is zero on all edges
        g = zeros(N, ne); % Neumann coefficient g is zero on all edges
        h = zeros(N^2, 2*ne); % Dirichlet h coefficient is zero at both ends of all edges
        r = zeros(N,2*ne); % Dirichlet r coefficient is zero at both ends of all edges
    end

The line N = 1 says that you have a single PDE (i.e. a scalar PDE) in your system. If you have more than one, you need to change this line. You can pass this function to assempde or the other high-level PDE Toolbox functions by defining

    b = @boundaryFileZeroNeumann;

This function is a very simple boundary file example. You can easily generalize this to an arbitrary combination of Dirichlet and Neumann conditions. Take a look at this documentation page if you are interested in doing that.

http://www.mathworks.com/help/pde/ug/boundary-conditions-for-scalar-pde.html

Bill

Best Answer

Related Solutions

MATLAB: Time Dependent Boundary Conditions in PDE Toolbox

MATLAB: Supporting documentation regarding PDE boundary files

Related Question