Heat transfer equation in transient case - computing the result with triple summation - for loop

Question

Michael50 2018-3-13

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/388047-heat-transfer-equation-in-transient-case-computing-the-result-with-triple-summation-for-loop

编辑： Michael50 2018-3-19

I try to compute transient temperature, T (t) from the below heat transfer equation (for unsteady case). There are three indices for the triple summation and for indices n & m, N0 = 120 (the minimum limit for better accuracy of the results) . And for index k used in gamma1 and gamma2, N0 = 10. Gamma_1(k) and Gamma_2(k) are the roots of some equation loaded from the text files with high decimal precision because the roots need more decimal precision for zero of the function, f(x) =0

Can somebody check the below code for triple summation? I have defined all the other variables in global declaration. I am not sure about the three nested for loops used for triple summation because I do not get the expected result for transient temperature, for different times like 100 us, 1 ms, 100 ms and so on... Below I give the complete heat transfer equation for transient analysis (the equation is easy just looks complicated because of the notations) and the code.

   clear; clc
global lambda
global alpha
global rho
global c
global Lx
global Ly
global Lz
global x
global y
global P
lambda = 399; % Thermal conductivity W/Km
alpha = 20; % Heat transfer coefficient W/Km^2
rho = 8940; % density in kg/m^3
c = 385; % specific heat in J/KgC
Lx = 50e-03; % length of heat sink in x direction in meters
Ly = 50e-03; % in meters
Lz = 0.5e-03; % thickness of the heat sink in meters
x = 25e-03; % middle point in meters
y = 25e-03; % middle point in meters
P = 968e-03; % Average power loss in Watt
Area = 0.56e-6; % area of the heat sources in meter^2
z = 269e-6; %t thickness of the silicon die in meters
x1 = pi.*1e-03/(2*Lx) % for n* pi* x1 /(2*Lx)
y2 = pi.*1e-03/(2*Ly) % m* pi* x2 /(2*Ly)
%x2a = pi.*50e-03/(2*Lx) % only included when needed
%y2a = pi.*50e-03/(2*Ly) % only included when needed
const = 8 / (Lx*Ly*Lz) % constant multiplied after the triple summation
load Roots_Y1.txt
format long % format long is used because the roots need more decimal precision for zero of the function, f(x) =0
load Roots_Y2.txt
format long 
Y1 =  Roots_Y1;
Y2 =  Roots_Y2;
sum_1 = 0;
for n = 1:128
    for m = 1:128
       for k = 1:10
        gamma = sqrt((n.*pi./Lx)^2 + (m.*pi./Ly)^2);
          Pnm = (1/(n*m*pi.^2)*(P/Area).*sin(n*x1).*sin(m*y2));
          A = [((2*Y1(k)/Lz*sin(Y1(k))*cosh(gamma*Lz/2)) + gamma*cos(Y1(k))*sinh(gamma*Lz/2))*cos(Y1(k)*(z-Lz)/Lz)] / [(1+((sin(2*Y1(k)))/(2*Y1(k)))) *((gamma^2)+(2*Y1(k)/Lz)^2) *((gamma*lambda*sinh(gamma*Lz/2)) + (alpha*cosh(gamma*Lz/2)))];     
          B = [((-gamma.*sin(Y2(k)) *cosh(gamma.*Lz/2)) + (2*Y2(k)/Lz) *cos(Y2(k)) *sinh(gamma.*Lz/2)) *sin(Y2(k)*(z-Lz)/Lz)] / [(1-((sin(2*Y2(k)))/(2*Y2(k)))) *((gamma.^2)+(2*Y2(k)/Lz)^2) *((gamma.*lambda.*cosh(gamma.*Lz/2)) + (alpha.*sinh(gamma.*Lz/2)))]; 
          p1 = (lambda/(rho*c))*[(n*pi./Lx)^2 + (m*pi/Ly)^2 + (2*Y1(k)/Lz)^2];
          p2 = (lambda/(rho*c))*[(n*pi./Lx)^2 + (m*pi/Ly)^2 + (2*Y2(k)/Lz)^2];
          t = 1e-3; % time at which the temperature is calculated at the centre of the substrate
          Temp = Pnm *((cos(n.*pi.*x/Lx))*(cos(m.*pi.*y/Ly)))*[(A*exp(-p1*t))+ (B*exp(-p2*t))]        
          sum_1 = sum_1 + Temp
         end 
      end
  end
T = sum_1
T_tr = (const * T) + 25 % ambient temperature 25 deg Celcius

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

请先登录，再进行评论。

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

Answer 1

Abraham Boayue 2018-3-15

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/388047-heat-transfer-equation-in-transient-case-computing-the-result-with-triple-summation-for-loop#answer_310130

Michael, You seemed to be on the right track, however, your temperature is a function of time. I see that you set t equal to a single digit number. Why did you do that? I think t should be defined as a vector. In the transient state, the temperature should be changing with time.

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

Torsten 2018-3-16

Maybe you could include a link to the book.

Michael50 2018-3-19

编辑：Michael50 2018-3-19

在 MATLAB Online 中打开

Hi Torsten, The book what I am referring to is for Hybrid circuits but as we know, the equations are valid for any circuits to find temperature. The book - Thermal modelling and optimization of power circuits by Andrzej Kos and Gilbert De Mey. My supervisor (the author of the book) suggested me to use the equations from this book. But now I am moving to other method since I am not getting the expected result using the transient equation for my problem. Now I am using PDE toolbox to model my transient thermal model as shown in the example https://www.mathworks.com/help/pde/examples/solving-a-heat-transfer-problem-with-temperature-dependent-properties.html (For my work, I import 3D model - .stl file) but now I am struggling to find out the node near the middle of a 3D model for transient heat transfer analysis to plot the temperature vs time. Since I import the 3D geometry, I do not know how matlab assigns xyz axis to my geometry. Please check my code and I have used comments to show it clearly. I attach my code and 3D geometry. Here is my code for quick reference

   clear;clc
% Computing the transient temeperature in the middle of the silicon die % (chip) 
% attached with the copper heat sink - 
thermalmodelT = createpde('thermal','transient'); % create PDE for transient thermal anaysis
gd =importGeometry(thermalmodelT,'ThermomodelT_Assembly.stl'); % importing geometry .stl file - dimensions in mm 
pdegplot(thermalmodelT,'FaceLabels','on','Celllabel','on'); % Viewing face and cell labels. So that boundary conditions can be implemented
xlim([-3,3])
ylim([-3,3])
thermalProperties(thermalmodelT,'ThermalConductivity',399,... % thermal - material properties of copper
                                'MassDensity',8940,...
                                'SpecificHeat',385,...
                                'Cell',1);
thermalProperties(thermalmodelT,'ThermalConductivity',120,... % thermal - material properties of silicon 
                                'MassDensity',2330,...
                                'SpecificHeat',800,...
                                'Cell',2);
thermalBC(thermalmodelT,'Face',[1,2,3,4,6], ...   % boundary conditions - covection coeff. - 20 W/(m^2•C) and T_amb = 25 degree Celcius
                       'ConvectionCoefficient',20,...
                       'AmbientTemperature',25)
%thermalBC(thermalmodelT,'Face',11,'HeatFlux',1728600)                                         
internalHeatSource(thermalmodelT,968000,'Cell',2) % power dissipation/ area of the chip = power density (Here it is 0.968 W/ 1e-6 m^2 = 968000 W/m^2) 
msh = generateMesh(thermalmodelT,'Hmax',0.2); % generate mesh size 
figure % plotting mesh modelling to check the mesh size suitable for the model
pdeplot3D(thermalmodelT);
axis equal
title 'Thermal model with Finite Element Mesh'
tlist = 0:100e-6:1; % time setup - let us say from 0 to 1 s in steps of 100us
thermalIC(thermalmodelT,25); % initial condition - at 25 degree C
R = solve(thermalmodelT,tlist); % solving the problem for the model
T = R.Temperature;
getClosestNode = @(p,x,y) min((p(1,:) - x).^2 + (p(2,:) - y).^2); %finding nearest nodes
[~,nid] = getClosestNode( msh.Nodes, 0.025, 0.025); % locating the node near the middle point ???????? (in my case, the middle point is 0.025 m and 0.025 m)
%(I do not know how to find that in my problem because I import my 3D model -I NEED HELP HERE) 
h = figure;
h.Position = [1 1 2 1].*h.Position;
subplot(1,2,1);
axis equal
pdeplot3D(thermalmodelT,'ColorMapData',T(:,end)) 
title 'Temperature, Final Time, Transient Solution'
subplot(1,2,2);
axis equal
plot(tlist, T(nid,:)); % transient temperature data at the node near the middle point ???????? (I do not know how to find that in my problem. I NEED HELP HERE)
grid on
title 'Temperature at the centre as a Function of Time';
xlabel 'Time, seconds'
ylabel 'Temperature, degrees-Celsius'

请先登录，再进行评论。

Answer 2

Abraham Boayue 2018-3-19

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/388047-heat-transfer-equation-in-transient-case-computing-the-result-with-triple-summation-for-loop#answer_310609

Hi Michael, I hope you have been able to come up with a solution to your heat transfer equation. If not, I have a simple question for you. Even though your triple for loops seemed to be well implemented from a programming perspective, are they really doing triple summation? If not, take a look at this image that I got from the internet. Can you verify by hand that the answer here is actually 45 and then test it on your function? You are going to be surprised of the result you get. I am certain that you will be able to fix your code after this exercise.

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

Michael50 2018-3-19

Hi Abraham, Before I have implemented triple summation in three for loops, I checked it with the simple calculations like this. Actually I have already seen this 'Performing a triple sum' video on youtube to verify my triple for loop implementation. In the above example, the answer is 396 (matlab and hand calculation also). In the video, he made a mistake and he has then accepted it in his comments below :) And till now I have not yet got the expected result. And now I have started with the implementation in partial differential equation toolbox for transient conduction heat transfer as shown in this example. https://www.mathworks.com/help/pde/examples/solving-a-heat-transfer-problem-with-temperature-dependent-properties.html. Please check above with my new code with using PDE toolbox in matlab

请先登录，再进行评论。