Nonlinear conduction,where thermal conductivity depends on temparature

[annie.1991]

Hi everyone,
i am trying to solve and program a non linear differential equation in MATLAB where thermal conductivity depends on temparature.I am trying it to solve by explicit finite difference method.
the given equation is ∂²t/∂x²+∂²t/∂y² *k(t)= -q (x,y)

i have solved the equation taking k(t)= a-b*t,and when i further solved the equation it gave an quadratic equation. where temparature is the root of the equation.but unfortunately my program is running slow and not giving me precize results can anybody help me out in this.

n = 5;
x = linspace(0,1,n);
dx = x(2)-x(1);
y = x;
dy = dx;

Q =80;


M= zeros(n,n);
if (mod(n,2)==0)
M(n/2,n/2)=Q;
else
p = (n+1)/2 ;
M(p,p)= Q;
end


t = ones(n);

t(1,1:n) = 0;
t(n,1:n) = 0;
t(1:n,1) = 0;
t(1:n,n) = 0;
a = 1;
b = 0;
tic
error = 1; z = 0;

while error > 0.0000001
z = z+1;

tnew = t;

for i = 2:n-1
for j = 2:n-1

d = -((b*(tnew(i+1,j)+tnew(i-1,j)+ tnew(i,j+1)+tnew(i,j-1))) + (4*a));

e = ( a*(tnew(i+1,j)+tnew(i-1,j)+ tnew(i,j+1)+tnew(i,j-1))) + (M(i,j)* dx^2);

f = b*4;

answer=roots([f,d,e]);
t(i,j)= max(max(abs(answer)));
end
end
error = max(max(abs(tnew-t)));

end
toc

 

[jlefevre76]

Ummmm. I'm not sure, but since you're using a while loop, the tic and toc lines are probably unnecessary? Here's how I would handle it:

T i,j = [(Q i,j/k i,j)*dx^2 + T i+1,j + T i-1,j + T i,j+1 + T i,j-1]/4

Q i,j is the Q matrix you set up.

k i,j would be a matrix with values of k estimated from the last iteration, it would be convenient to store this in a matrix as well.

T i,j is the temperature at point i,j, and the others are offset from i,j by one in the + or - x or y directions (similar to what you already have).

As for the root finding lines that you threw in there, they are completely and totally unnecessary. As long as you update the values of k for each index of the matrix in each loop, the values will converge on their own over time. Your calculation of the error should work okay, but I don't think the number you chose for the error limit needs to be that small. 0.0000001 is extremely small, and it could take a long time to get it that low. Try something larger first, like 0.001. Then solve it down to 0.0001 and see if that changes anything, likely it won't make much of a difference to the final answer you get.

 

[annie.1991]

Thanks for your reply and your sujjestion,
Probably What i can undestand from your sujjestion is the linear heat conduction,but i am trying to do a non linear heat conduction where thermal conductivity (K) depends on Temparature (T).If you would have any better idea of solving this kind of equation please share with me

I didnt undestand what you are saying about K(i,j),,what exactly is this K(i,j) ?

 

[jlefevre76]

So, K (conductivity) has it's own matrix, where i and j are the indices (for any of the mentioned variables). Even if K is non-linear as a function of temperature, typically you can still use the same method and get convergence of your array.

I'm a little bit confused, above you have written k = a-b*t, which is linear. What is the equation that describes k as a function of temperature? True, the differential equation becomes non-linear after using a temperature dependent conductivity, but the conductivity itself appears to be linear.

Even if it isn't, as long as you have a function for conductivity versus temperature, it should converge (unless it's a really weird relationship in terms of conductivity going to zero or negative numbers).

 

[annie.1991]

i understood about k(i,j).if i take initial values of k as ones matrix and run the program with your given sujjestion then there will be no influence or no change in values of martix K.

in my code I have considerd K(T) = a - b*T where a>>b a and b are variables. i just taught i can see the variations of temperature depending on the values of a and b.I just wanted to write down in the equation where tempearture is changing with some other variables.

So If you feel like my idea doesn't works (or) my idea is wrong with the physics,then please share any other ideas which you have to solve this equation where i can see the conductivity changing with temperature.

 

[Chestermiller]

You formulation for the case of temperature-dependent thermal conductivity is incorrect. You should not have factored out the thermal conductivity from the derivatives. It should be:

\frac{∂}{∂x}\left(k\frac{∂T}{∂x}\right)+\frac{∂}{∂y}\left(k\frac{∂T}{∂y}\right)=-q
or equivalently
\frac{1}{k}\frac{∂}{∂x}\left(k\frac{∂T}{∂x}\right)+\frac{1}{k}\frac{∂}{∂y}\left(k\frac{∂T}{∂y}\right)=-\frac{q}{k}
or equivalently
\frac{∂^2T}{∂x^2}+\frac{∂^2T}{∂y^2}+\frac{∂lnk}{∂x}\frac{∂T}{∂x}+\frac{∂lnk}{∂y}\frac{∂T}{∂y}=-\frac{q}{k}

Chet

 

[jlefevre76]

You formulation for the case of temperature-dependent thermal conductivity is incorrect. You should not have factored out the thermal conductivity from the derivatives.

Ummmmm, I could be wrong, but Chet's derivation assumes k is a function of x and/or y? This is not the case you've laid out here, where k is only a function of temperature. True, temperature is also a function of the location, but that's why you calculate a new value for k at each iteration and location. (The derivation is not wrong, just unneccesary from what I've always seen done for a numerical method.) Given the conditions described, a>>b, it should have no trouble converging. Because you're solving the equation one element at a time, you can pull k out of it so long as k at each element gets calculated and is updated each iteration.

 

[Chestermiller]

Ummmmm, I could be wrong, but Chet's derivation assumes k is a function of x and/or y?

Not exactly. k is a function of temperature, and temperature is a function of x and y, so k is a function of x and y. In evaluating the derivatives involving lnk in the equation, you get k(T) at each of the grid points, and then use a central difference approximation to get the derivatives at the grid point under consideration. If you omit those derivatives of lnk, you will get the wrong answer to the problem. I have a lot of industrial experience in solving problems like this in the real world.

Another form of the equation that works is:

\frac{∂^2T}{∂x^2}+\frac{∂^2T}{∂y^2}+\frac{∂lnk}{∂T}\left[\left(\frac{∂T}{∂x}\right)^2+\left(\frac{∂T}{∂y}\right)^2\right]=-\frac{q}{k}

where ∂lnk/∂T is evaluated at the grid point under consideration.

Chet