Solving Steady State Heat Conduction Eqn w/ MATLAB

[range.rover]

HI guys,this is my first programming experience , i have developed an MATLAB code for steady state heat conduction equation , on governing equation
dt² /dx² + dt²/dy² = -Q(x,y)

i have solved this equation with finite difference method, As far as i know if we increase the mesh size it leads to decrease in the distance between the nodes and leads us to the true solution.i guess the final solution will be constant.

but my code is not converging or moving towerds true solution can anybody help me out on this. Here the solution is just varying with the mesh size.
what can be the reason behind this.

this is my code.

%%%%%% STRICTLY ONLY FOR ODD NUMBERS %%%%%%%%%


clear all
close all
clc
n = 101;% n grids and has n interior points per dimension
m = n+2;
x = linspace(0,1,m); % grid points x including boundaries
y = linspace(0,1,m); % grid points y including boundaries
dx = x(2)-x(1);
h = dx;
dy = dx;
k=1; %thermal conductivity
Q=8000; % heat source
Iint = 2:n+1;
Jint = 2:n+1;
%Right hand side matrix
B=zeros(n,n);
if (mod(n,2)==0)
B(n/2,n/2)=Q;
else
B((n+1)/2,(n+1)/2)=Q;
end

%Main Matrix
Told = zeros(m,m);

% boundary values
Told(1:m,1) = 0;
Told(1:m,m) = 0;
Told(1,1:m) = 0;
Told(m,1:m) = 0;
Tsoln = Told;
B(:,1) = B(:,1) + Told(Iint,1)*k/h^2; % boundary valuues are added from the left hand side matrix including k and h
B(:,n) = B(:,n) + Told(Iint,m)*k/h^2;
B(1,:) = B(1,:) + Told(1,Jint)*k/h^2;
B(n,:) = B(n,:) + Told(m,Jint)*k/h^2;

% Matrix formations
F = reshape(B,n*n,1);
I = speye(n);
e = ones(n,1);
T = spdiags([e -4*e e],[-1 0 1],n,n);
% v = full(T);
S = spdiags([e e],[-1 1],n,n);
% w = full(S);
% v = kron(I,T);
% w = kron(S,I);
% z = (kron(I,T) + kron(S,I));
% nny = full(z);
A = (kron(I,T) + kron(S,I))* k / h^2; % A is a tridiagonal spurse matix multiplied and dvided by h^2
% z = full(A);
% Tvec = abs(A\F); Solving AX = B ;

Tvec = abs(mldivide(A,F));
Tsoln(Iint,Jint) = reshape(Tvec,n,n);
sol = max(max(Tvec)) ;

% Ploting
figure,surf(x,y,Tsoln),
figure, contour(x,y,Tsoln),
title('Temperature (Steady State)'),xlabel('x'),ylabel('y'),colorbar
figure,pcolor(x,y,Tsoln), shading interp
title('Temperature (Steady State)'),xlabel('x'),ylabel('y'),colorbar

 

[Dr Transport]

What are your boundary conditions, more fully explain the problem and we can be more helpful

 

[range.rover]

thanks for your reply
i am solving in a two dimensional staedy state heat conduction equtation with a heat source at the center in square plate with dirichlet boundary conditions, i mean to say values on boundries of four sides are same and constant.
dt² /dx² + dt²/dy² = -Q(x,y)
i solved it by finite difference method
B(i,j) is a right hand side matrix with heat source at the center.
i am solving it by direct method, creating right hand side matrix and main spurse matrix. later i am dividing them for solution.
My discritised equation is
T(i,j) = ((T(i+1,j) +T(i-1,j) +T(i,j+1)+T(i,j-1) )*k)+(B(i,j)*dx^2)))/(4*k)

 

[gsal]

Why are you multiplying T by k/h^2 ?

If this is a 2D problem you are trying to solve, the cross sectional area shared by two adjacent cells is the length of the cell side times the depth...in 2D problems, one usually assumes 1 unit of length deep...

...and, so, I an thinking, you should be multiplying by k/(h*1)

...otherwise, you are changing the depth of your model and that is why you are getting different temperature every time.

 

[gsal]

Or...if I am way off, never mind my point...I don't think I understand how you are setting up the problem, where your cells are, where your nodes are and whether the boundary nodes are numbered and part of your program (they shouldn't be). I certainly don't understand your choice of m (n = n+2), or why you presume your stiffness matrix A to be tridiagonal (it is only triadiagonal for a 1D mesh where every node had only 2 neighbors).

 

[range.rover]

Hey gsal, I appreciate your reply

when i discritised this equation "d²t/dx<sup>2 + d2t/dx2 = -Q " i wrote the linear eauation
((T(i+1,j) +T(i-1,j) +T(i,j+1)+T(i,j-1) -4*T(i,j))*k)/ h2 = Q,
above in my program h is the distance between the nodes
this has lead me to the tridiagonal matrix or a sparse matrix ,framing a system of linear equations ax = b

after seeing your message i tried it by using h^1, but it didnt work.i think we can not change that, as far as your second doubt, i am trying to solve n nodes of mtrix dimension of n+2,to set the boundary conditions to it.

if you still feel that there is some thing wrongwith my code. feel free to edit my program and do it in your way. i want to see that too.</sup>

 

[gsal]

Forgive me if I choose not to write the program for you. First, I don't use matlab; although it does not look like it would be much different from Python.

Anyway, setting up this problem with a few loops is not that difficult, but you really need to know exactly what you need to do...if I were you (I already was and did it), I would really do my homework and deduce the equations by hand, first.

They may sound like too many but I think you need at least a 3x3 grid (9 equations) to realized the pattern of the equations; otherwise, you should go for 4x4.

If you welcome my advice, forget about your dt^2/dx^2...equation for now.

I would take the 3x3 set of nodes, number them from 1 to 9 and setup Kirkchhoff's Equations of heat flowing away from each node toward its neighbors and BC "nodes" (these are not numbered, they are actual numeric values in the equations and the conductance connection to this temperature should half the one between actual nodes...you'll see).

...THEN, you will see a nice pattern on how each diagonal and off-diagonal entry looks like...I am thinking you really 4x4 cells, but you don't really need to set up all 16 equations, just the unique ones...10 ;-) ...4 corners, 4 sides, 1 diagonal (not corner), 1 off-diagonal (not side).

To me, the temperature difference (dT) between two nodes equal the heat flow (Q) times the thermal resistance (R) in between, and so: dT = QR, where R=ρL/A = L/(κA) = dx/(κ*dy*1) = 1/κ, if dx=dy...so, I don't see how h^2 came to be in your equations.

So, do Kirkchhoff's equations and organize them in the form Ax=b

• b is the known 1D vector of heat injections into the system (from center and boundary conditions)
• x is the 1D vector of temperatures to solve for
• A is the 2D stiffness (conductance) matrix

 

[range.rover]

i appreciate for your conversation,
i have done my homework and solved it numerically as you said taking 3*3 matrix and i happen to program my code in the same way i have calculated, now i get to know why you are facing the the problem,if you just run my program line to line you can check in MATLAB what am i doing, i have just done the same thing what you have written up there.
i have checked my program many times by calculatng with hands numerically by taking 3*3 and they turn out to be same.so i had no clue why this program is not converging if i inrcrease the mesh size.
let me clarify you again h^2 is distance between the nodes and it came in my equation while discritizing the main equation through Finite difference method.
i really didnt understand what are you talking about (dT).

 

[gsal]

I am a steady-state, finite-difference kind of guy and, so, I stay away from differential equations.

In my equation, dT is nothing else than (T₁ - T₂), for example:
T₁ - T₂ = Q₁₂R₁₂
and so
Qij = (Ti - Tj) / Rij
and so, the summation of the various (Ti - Tj) / Rij for all j connected to the same i is equal to injection Qi at node i.

Here, the denominator is thermal resistance connecting node i to node j, not just simply the distance between them, nor squared; like I illustrated before when I expanded/simplified dT=QR, in previous post, dx and dy cancel out leaving k/1, yielding the correct units on both sides of the equation...for a 2D problem, 1 unit deep, the units on both sides of Q = (T1-T2)/R are watts per meter squared (W/m²)...as it should be.

If you do a dimensional analysis in your linearized equation, the left hand side has temperature, times thermal conductivity divided by linear squared: K*W/Km/m² = W/m³...that's a volumetric heat density...how can that be the heat flow across the wall shared between two cells? The units need to be W/m².

So, I think I going back to what I said before...because the units of the left hand side are W/m³ and you have h² in there, the meaning of Q on the other side is not an absolute, fix amount of watts every time, instead it is a density...so, every time you change the size of the one element with heat generation, you change the amount of total heat being injected into the system and when you solve whatever problem you are solving, you of course get a different answer...that's all I can speculate.

 

[gsal]

Actually, let me take the part about W/m² half-way back...sure, steady-state heat rate through a shared wall can only be in W/m², but the heat from one node to the other (after multiplying by the area of the wall) should truly just be in plain W (watts), that's it.