Transient Neumann boundary condition

Good evening,

I would like to simulate a heat transfer problem with the PDE toolbox and I am trying to apply a transient heat flux on one edge of a rectangle. The other edges have either adiabatic or constant boundary conditions.

Is it possible to put a transient heat flux function (Power(t) = P0 + P0 cos (2wt) with P0 the amplitude in W, w the pulsation in rad/s and t the time in s) directly with the applyBoundaryCondition function, do I have to use the setInitialConditions function to use the data coming from the last solvepde or am I totally wrong and I should use another method ?

Thank you for your help!

% Heat flux function of time which has ot be applied on edge 3
P = 0.0035;                                     % Power (W)
length = 0.0003;                                % Length (m)
f=2;                                            % Frequency (Hz)    
W = 2*pi*f;                                     % Pulsation (rad/s)
t_inc = 0.025;                                  % time increment
t_end = 5;                                      % time end
i=1;                                            % matrix increment
for t = 0:t_inc:t_end
    Power(i) = (P+P*cos(2*W*t))/length;         % Power matrix (W/m)
    time(i)=t;                                  % time matrix
    i=i+1;
end

% Geometry
numberOfPDE = 1;                                % Number of equation
model = createpde(numberOfPDE);                 % Create a PDE with numberOfPDE equation

% Define the dimension of the substrate
width = 0.0006; 
height = 0.0003;
gdm = [3 4 0 width width 0 0 0 height height]'; % Define the geometry
g = decsg(gdm, 'S1', ('S1')');                  % Decompose constructive solid
geometryFromEdges(model,g);

%% Set boundary conditions. 
setInitialConditions(model,Ta);                 % Specify the initial temperature on all nodes
applyBoundaryCondition(model,'dirichlet','edge',1,'u',Ta);
applyBoundaryCondition(model,'neumann','edge',3,'g',Power(1));

Answer 1

[Ravi Kumar]

You should be able to apply the Neumann BC using the power calculation that you have done in the beginning. But be sure to take care of units, you need heat flux, which is W/m^2, in 2-D case it would be W/m.

Regards,

Ravi

[Mickael]

Thank you for your reply and your remark. I am going to have a look on it.

[Mickael]

Following your suggestion I build a power function

P = 0.035;                                      % Power (W)
heater_width = 0.00001;                         % Heater width (m)
f = 1;                                          % Frequency (Hz)
Pfunction = @(location,state) (P + P*cos(4*pi*f*state.time))/heater_width;     % Heat flux vs time

then added it as a boundary condition in thermalBC function such as

thermalIC(model,Tini);
thermalBC(model,'Edge',1,'Temperature',Tini);
thermalBC(model,'Edge',4,'HeatFlux',Pfunction);

and solve it with time such as

tfinal = 5;
tstep = 0.025;
tlist = 0:tstep:tfinal;
results = solve(model,tlist);

however, the temperature along this edge 4 is not what is expected. One should observe an oscillation of the temperature instead of this :

So I wonder if the heat flux is taken correctly in consideration. What do you believe?

Thank you in advance,

Regards,

Mickaël

[Ravi Kumar]

Hi Mickael,

Could you post the compete reproduction code? I will take look.

Regards,

Ravi

[Mickael]

Hi Ravi,

Please find the code below. I have changed the geometry compared to the code above in order to be closer to the final sample. It doesn't look better and so if the code is ok, I wonder if the problem could come from the meshing which is not much refined close to the heater as sligth heater dimension variation leads to different T vs t figures...

Thanks again,

Regards,

Mickaël

clc
clear all
Tini = 294;        % Temperature of the system
model = createpde('thermal','transient');
% Geometry
width = 0.0007; 
height = 0.0003;
heater_width = 0.00001;
heater_height = 0.0000001;
r1 = [3 4 0 width width 0 0 0 height height];
r2 = [3 4 width/2-heater_width/2 width/2-heater_width/2 width/2+heater_width/2 width/2+heater_width/2 ...
    height height+heater_height height+heater_height height];
gdm = [r1; r2]';
g = decsg(gdm,'S1 + S2',['S1'; 'S2']');
geometryFromEdges(model,g);
% Thermal properties
thermalProperties(model,'Face',1,'ThermalConductivity',317,...
                        'MassDensity',1930,...
                        'SpecificHeat',128);
thermalProperties(model,'Face',2,'ThermalConductivity',1.5,...
                        'MassDensity',2650,...
                        'SpecificHeat',730);
% Determine the power function
P = 0.035;                                          % Power (W)
f = 2000;                                           % Frequency (Hz)
Pfunction = @(location,state) (P + P*cos(4*pi*f*state.time))/heater_width;     % Power vs time
% Boundary conditions
thermalIC(model,Tini);
thermalBC(model,'Edge',1,'Temperature',Tini);
thermalBC(model,'Edge',[4,5,6,8],'HeatFlux',Pfunction);
% Meshing
hmax = .00001;                      % maximum element size
msh = generateMesh(model,"Hgrad",2,'Hmax',hmax,'GeometricOrder',"quadratic");
figure 
pdeplot(model);
title 'substrate With Triangular Element Mesh'
xlabel 'X-coordinate, meters'
ylabel 'Y-coordinate, meters'
%Transient solution
tfinal = 5;
tstep = 0.025;
tlist = 0:tstep:tfinal;
model.SolverOptions.RelativeTolerance = 1.0e-3; 
model.SolverOptions.AbsoluteTolerance = 1.0e-4;
results = solve(model,tlist);
T = results.Temperature;
figure 
pdeplot(model,'XYData',T(:,end),'Contour','on','ColorMap','jet')
title(sprintf('Temperature In The substrate, Transient Solution( %d seconds)\n',tlist(1,end)))
% Find the temperature in the middle of the heater
Tmiddle = interpolateTemperature(results,[width/2;height+heater_height/2],1:numel(tlist));
figure
plot(tlist,Tmiddle)
title 'Temperature at the middle of the heater as a function of Time'
xlabel 'Time, s'
ylabel 'Temperature, K'
grid on
% Find temperature along sample vs time
y = 0:height/100:height;
x = width/2*ones(1,length(y));
Tmiddle_final = interpolateTemperature(results,x,y,tfinal);
figure
plot(Tmiddle_final)
title 'Temperature at the middle of the heater as a function of Time'
xlabel 'Y-coordinate, m'
ylabel 'Temperature, K'
grid on

[Ravi Kumar]

Hi Mickael,

I think your setup looks fine. You seem to have some unit inconsistency. In the power variation function, you have time constant of 1/(4*pi*2000) ~ 4E-5. So the dynamics is on the time-scale of 1E-5, but you are solving from 0 to 5, which also explains the initial fluctuation.

Regards,

Ravi

[Mickael]

Hi Ravi,

Thank you very much to point this. However, if you allow me, I still have a question regarding this problem because the temperature variation doesn't look like they should, the temperature variation are not a nice cosine function and above 5 periods there is no more periods... I would expect to observe cosine function whatever the number of period. I do not understand why increasing the final time would modify what has been calculated by the model at a lower time, is it expected ? is there still something that I am missing ?

Please find below the try I have performed with an arbitrary frequency of 100Hz. The time step is adjusted in order to always have the same number of point (100 here) to describe a period.

Nperiod = 2;
NpointPerPeriod = 100;
tfinal = Nperiod/(2*f);
tstep = tfinal/(Nperiod*NpointPerPeriod);
tlist = 0:tstep:tfinal;

For 2 period I obtain this figure :

For 5 period :

For 7 period:

For 10 period :

Thank you again for your very helpful help.

Best regards,

Mickaël