How to implement time dependent heat generation

The code at the end of this message produces the transient temperature history for a square with constant temperature boundary conditions and constant heat generation. I would like help to modify the code to produce the transient temperature history when the heat generation is time dependent. A specific example would be to have the heat generation switched off and on from 0 to 4000. The logic would be something like that below:

tlist = 0:0.1:10;

% Initialize x with zeros
Qgen = zeros(size(tlist));

% Assign Qgen = 0 for 0 <= t < 2
Qgen(tlist >= 0 & tlist < 2) = 0;

% Assign Qgen = 4000 for 2 <= t <= 4
Qgen(tlist >= 2 & tlist <= 4) = 4000;

% Assign Qgen = 0 for 4 <= t <= 6
Qgen(tlist >= 4 & tlist <= 6) = 0;

% Assign Qgen = 4000 for 6 <= t <= 8
Qgen(tlist >= 6 & tlist <= 8) = 4000;

% Assign Qgen = 0 for 8 <= t <= 10
Qgen(tlist>= 8 & tlist <= 10) = 0;

Judging by the MatLab PDE Users Guide, I think that I should be using something like:

fcoeff= @(location, state) myfuncWithAdditonalArgs(location, state, arg1, arg2…..)

and then

SpecifyCoefficents(thermalmodel, f@fcoeff)

where the location would be F1 and the state and arg parts would be set up to deal with the off-on logic defined above. Is this the correct approach and if so, can you please help me with the correct way to implement the off-on logic example given above?

Thanks

%% Create thermal model PDE

thermalmodel = createpde(1);

%% Create Geometry

R2= [3,4,-1.5,1.5,1.5,-1.5,-1.5,-1.5,1.5,1.5]';
geom=[R2]
g=decsg(geom)
model= createpde
geometryFromEdges(thermalmodel,g);
pdegplot(thermalmodel,"EdgeLabels","on","FaceLabels","on")
xlim([-2 2])
axis equal

%% Generate and plot mesh

generateMesh(thermalmodel);

figure
pdemesh(thermalmodel)
title("Mesh with Quadratic Triangular Elements")

%% Apply BCs

applyBoundaryCondition(thermalmodel, "dirichlet",Edge=[1],u=0);
applyBoundaryCondition(thermalmodel, "dirichlet",Edge=[2],u=0);
applyBoundaryCondition(thermalmodel, "dirichlet",Edge=[3],u=0);
applyBoundaryCondition(thermalmodel, "dirichlet",Edge=[4],u=0);

%% Apply thermal properties

rho1= 1                                                         % 
cp1= 1                                                           % 
rhocp1= rho1*cp1                                   %
k1= 1                                                             % W/mK

%% Define uniform volumetric heat generation rate

Qgen= 4000                              % W/m3                   

%% Define coefficients for generic Governing Equation to be solved

f= [Qgen]

specifyCoefficients(thermalmodel, m=0, d=rhocp1, c=k1, a=0, f=f, Face=1);

coeffs= thermalmodel.EquationCoefficients;
findCoefficients(coeffs, Face=1)

% Define initial condition

setInitialConditions(thermalmodel, 20);

% Define time duration and intervals

tlist = 0:0.1:10;

% Solve PDE

thermalresults = solvepde(thermalmodel,tlist);

T = thermalresults.NodalSolution;

getClosestNode = @(p,x,y) min((p(1,:) - x).^2 + (p(2,:) - y).^2);

% Call this function to get a node at the center of the square
[~,nid] = getClosestNode(thermalresults.Mesh.Nodes,0,0);                   

% Plot temperature history at location x,y

figure
plot(tlist, T(nid,:)); 
% plot(tlist, sol(nid,:)); 

grid on
title("T(nid,:) Middle of square")
xlabel("Time, seconds")
ylabel("Temperature, degrees-Celsius")

[Torsten]

Create a loop in which you alternate between 0 and 4000 for f and in which you take the solution for the last time of the previous 2-seconds-step as initial temperature distribution for the next 2-seconds-step. This way, you don't need a coefficient function f and you circumvent the problem of discontinuous heat generation over time.

Answer 1

[Torsten]

I leave the rest to you:

%% Create thermal model PDE
thermalmodel = createpde(1);
 
%% Create Geometry 
R2= [3,4,-1.5,1.5,1.5,-1.5,-1.5,-1.5,1.5,1.5]';
geom=[R2];
g=decsg(geom);
model= createpde;
geometryFromEdges(thermalmodel,g);
pdegplot(thermalmodel,"EdgeLabels","on","FaceLabels","on")
xlim([-2 2])
axis equal

%% Generate and plot mesh 
generateMesh(thermalmodel);
figure
pdemesh(thermalmodel)
title("Mesh with Quadratic Triangular Elements")

%% Apply BCs
applyBoundaryCondition(thermalmodel, "dirichlet",Edge=[1],u=0);
applyBoundaryCondition(thermalmodel, "dirichlet",Edge=[2],u=0);
applyBoundaryCondition(thermalmodel, "dirichlet",Edge=[3],u=0);
applyBoundaryCondition(thermalmodel, "dirichlet",Edge=[4],u=0);
 
%% Apply thermal properties
rho1= 1;                                                       % 
cp1= 1;                                                           % 
rhocp1= rho1*cp1;                                   %
k1= 1;                                                             % W/mK
 
%% Define uniform volumetric heat generation rate 
Qgen= 4000;                              % W/m3                   
 
%% Define coefficients for generic Governing Equation to be solved
tlist = 0:0.1:2;
specifyCoefficients(thermalmodel, m=0, d=rhocp1, c=k1, a=0, f=0, Face=1);
setInitialConditions(thermalmodel, 20);
thermalresults = solvepde(thermalmodel,tlist);
figure(1)
pdeplot(thermalresults.Mesh,"XYdata",thermalresults.NodalSolution(:,end))

●

getClosestNode = @(p,x,y) min((p(1,:) - x).^2 + (p(2,:) - y).^2);
 
% Call this function to get a node at the center of the square
[~,nid] = getClosestNode(thermalresults.Mesh.Nodes,0,0);    
T = thermalresults.NodalSolution(nid,:);
Time = tlist;

tlist = 2:0.1:4;
specifyCoefficients(thermalmodel, m=0, d=rhocp1, c=k1, a=0, f=Qgen, Face=1);
setInitialConditions(thermalmodel,thermalresults); 
thermalresults = solvepde(thermalmodel,tlist);
figure(2)
pdeplot(thermalresults.Mesh,"XYdata",thermalresults.NodalSolution(:,end))

T = [T,thermalresults.NodalSolution(nid,:)];
Time = [Time,tlist];

tlist = 4:0.1:6;
specifyCoefficients(thermalmodel, m=0, d=rhocp1, c=k1, a=0, f=0, Face=1);
setInitialConditions(thermalmodel,thermalresults); 
thermalresults = solvepde(thermalmodel,tlist);
figure(3)
pdeplot(thermalresults.Mesh,"XYdata",thermalresults.NodalSolution(:,end))

T = [T,thermalresults.NodalSolution(nid,:)];
Time = [Time,tlist];

tlist = 6:0.1:8;
specifyCoefficients(thermalmodel, m=0, d=rhocp1, c=k1, a=0, f=Qgen, Face=1);
setInitialConditions(thermalmodel,thermalresults); 
thermalresults = solvepde(thermalmodel,tlist);
figure(4)
pdeplot(thermalresults.Mesh,"XYdata",thermalresults.NodalSolution(:,end))

T = [T,thermalresults.NodalSolution(nid,:)];
Time = [Time,tlist];

tlist = 8:0.1:10;
specifyCoefficients(thermalmodel, m=0, d=rhocp1, c=k1, a=0, f=0, Face=1);
setInitialConditions(thermalmodel,thermalresults); 
thermalresults = solvepde(thermalmodel,tlist);
figure(5)
pdeplot(thermalresults.Mesh,"XYdata",thermalresults.NodalSolution(:,end))

T = [T,thermalresults.NodalSolution(nid,:)];
Time = [Time,tlist];
figure(6)
plot(Time,T)

[John McGrath]

This is terrific. Thanks for your help, Torsten