PDE Toolbox. Parabolic heat transfer: How to specify a forcing term to act only on the boundary?

I am interested in modeling transient heat transfer in a woody plant stem, and I have written the following code for the PDE Toolbox:

%Global constants
emiss=0.95; %Emissivity [Bark surface (%)]
StefanBoltz=5.670367e-8; %Stefan-Boltzmann constant (W/(m^2-K^4))

%Model constants
MC=0.30; %Moisture content, MC [Wood (%)]
Gb=0.418; %Specific gravity, [Wood (Dimensionless)]
Text=299.6967; %External temperature, T [Air (K)]
[wrho,wk,wcp]=thermprop(MC,Gb,Text);
th=60; %Heat up time, t (sec)
tc=1200; %Cool down time, t (sec)
Q=17133.1734201; %Heat source (W/m^2)

%Create model
model=createpde();

%Create geometry
C1=[1; 0; 0; 0.025];
[a]=decsg(C1);
geometryFromEdges(model,a);

%Specify boundary conditions
radfun = @(~,state) -emiss.*StefanBoltz.*(state.u.^4-(Text.^4));
applyBoundaryCondition(model,'edge',1:4,'g',radfun);

%Specify PDE coefficients
m = 0;
d = wrho.*wcp;
c = wk;
a = 0;
f = @(~,state) Q.*(heaviside(state.time)-heaviside(state.time-th));
specifyCoefficients(model,'m',m,'d',d,'c',c,'a',a,'f',f);

%Create mesh
generateMesh(model,'Hmax',0.002);

%Specify initial conditions
setInitialConditions(model,Text);

%Solve
dt=0.01;
t=0:dt:th+tc;
result=solvepde(model,t);

Currently, the forcing term (f coefficient) acts over the entire geometry for a 60 sec period defined by a linear combination of Heaviside functions, but I'd like to restrict it to act only at the boundary defined by x.^2+y.^2-0.025^2=0. Is it possible to do this using the PDE Toolbox?

Answer 1

[Daniel Burcham]

Thanks for your excellent suggestion. Ultimately, I found it easier to use the boundary conditions to define heat flux!

Answer 2

[michio]

It may not an elegant way. Since your boundary is defined in clear fashion, I think it's worth a try.

Instead of

 f = @(~,state) Q.*(heaviside(state.time)-heaviside(state.time-th));

define

 f = @(region,state) fcoeffunction(region,state,Q,th);

where

function f = fcoeffunction(region,state,Q,th)

    N = 1; % Number of equations
    nr = length(region.x); % Number of columns
    f = zeros(N,nr); % Allocate f

    Qt = Q.*(heaviside(state.time)-heaviside(state.time-th));
    dist2 = (region.x.^2) + (region.y).^2;
    index =  dist < 0.025^2 + eps & dist > 0.025^2 - eps;
    f(1,index) = Qt*ones(size(index));

The intention is for f to have non-zero values only where (region.x.^2) + (region.y).^2 is close "enough" to 0.025^2.

Since your script above is missing some of the function, I did not test the above function in your context, so please use this as a reference. Let me know how it goes.

The above follows an example shown in https://www.mathworks.com/help/pde/ug/f-coefficient-for-specifycoefficients.html