How to implement time dependent heat generation

Question

John McGrath 2025-5-6

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2176945-how-to-implement-time-dependent-heat-generation

评论： John McGrath 2025-5-7

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")

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Torsten 2025-5-6

编辑：Torsten 2025-5-6

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 2025-5-7

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2176945-how-to-implement-time-dependent-heat-generation#answer_1564784

编辑：Torsten 2025-5-7

在 MATLAB Online 中打开

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)

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

John McGrath 2025-5-7

This is terrific. Thanks for your help, Torsten

请先登录，再进行评论。

How to implement time dependent heat generation

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

采纳的回答

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

How to implement time dependent heat generation

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

采纳的回答

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论