Transient Neumann boundary condition

Question

Mickael 2020-5-1

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/522315-transient-neumann-boundary-condition

评论： Mickael 2020-5-5

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Ravi Kumar 2020-5-2

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/522315-transient-neumann-boundary-condition#answer_429705

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

6 个评论
显示 4更早的评论隐藏 4更早的评论

Mickael 2020-5-4

在 MATLAB Online 中打开

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

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

在 MATLAB Online 中打开

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

请先登录，再进行评论。

Transient Neumann boundary condition

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

回答（1 个）

6 个评论
显示 4更早的评论隐藏 4更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Transient Neumann boundary condition

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

回答（1 个）

6 个评论 显示 4更早的评论隐藏 4更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

6 个评论
显示 4更早的评论隐藏 4更早的评论