How can solve the 2d transient heat equation with nonlinear source term?

I am trying to solve the below problem of 2d transient heat equation. I have created the code for the simple case of 2d with constant temperatures at the boundary and it did work. Then, I included a convective boundary condition at the top edge, and symmetric boundary condition (dT/dn = 0) at the other three edges. At this time the problem appeared.

- When the outside temperature at the top edge is equal to the initial temperature of the plate (T_F=T0=300), it gave me totally unstable solution - When the outside temperature at the top edge is higher than the initial temperature of the plate (T_F>T0=300), a negative slope appears at the top edge, while it should be positive. Also the bottom edge appears to have a positive slope, which should be zero (not positive).

I am using the succesive over-relaxation iterative method, I have checked the finite difference formulation so many times and still do not know where the problem is. I have really been struggling with this problem for several days. After solving this issue, I will have to include the source term as well.

The code I have written is as follows:

%TRANSIENTHEAT EQUATION%

%%

clc clear all

i_max = 15; j_max = 20;

N = (j_max-2)*(i_max-2); %Number of unknowns

fprintf('\n The number of unknowns is %g \n ',N)

Lx = 0.3; % System size in x direction(length)

Ly = 0.2; % System size in y direction(length)

dx = Lx/(i_max-1); % Grid spacing in x %h%

dy = Ly/(j_max-1); % Grid spacing in y

dt=1; %Time increment

Lt=20; %Final time

t_max = Lt/dt+1

T0=300;

beta = dx/dy; %Grid aspect ratio

fprintf('\n The grid aspect ratio is %g \n ',beta)

change_want = 1e-4; % Stop when the change is given fraction

    x=linspace(0,Lx,i_max);

    y=linspace(0,Ly,j_max);

%% Properties

alpha = 97.1e-6;

rho = 2702;

cp = 903;

%% ------- Initialization and Initial Condition------ %%.

sx=alpha*dt/(dx^2);

sy=alpha*dt/(dy^2);

        T_old=zeros(j_max,i_max,t_max);

for j = 1:(j_max) for i = 1:(i_max) T_old(j,i,1)=T0; end end

%% ------------ Source Term --------------%%

Q = zeros(j_max,i_max,t_max); % initialization

%% Other Properties %%

k = 200;

h = 200;

T_F = 310; %%%the solution is unstable(fluctuating) when T_F = T0, and the solution does not have zero gradient at the bottom when T_F~=T0

c1 = (-h/k);

c2 = h*T_F/k;

%% ---------------- MAIN LOOP --------------%%

T_new = T_old;              

zeta = ((cos(pi/(i_max-1))+(beta^2*... 

 cos(pi/(j_max-1))))/(1+beta^2))^2;

omega = 2/zeta*(1-sqrt(1-zeta));

 fprintf('\n the optimum omega  is %g \n ',omega) %%!!%%SOR only

max_iter = i_max*j_max; % Set max to avoid excessively long runs

for t=1:((t_max-1))

for iter=1:max_iter*20 %iterations

temp = 0;

 for j=1:j_max

 for i=1:i_max

%% SOR method %%

      if  j == 1  

            if i==1

                T_new(j,i,t+1) = 1/(1+2*sx+2*sy+2*sy*dy*c1)*omega*(2*sy*T_old(j+1,i,t+1)-2*sy*dy*c2+ ...
                                +2*sx*T_old(j,i+1,t+1)+T_old(j,i,t)+dt*Q(j,i,t+1)/rho/cp) +  (1-omega)*T_old(j,i,t+1);

            elseif i==i_max

                T_new(j,i,t+1) = 1/(1+2*sx+2*sy+2*sy*dy*c1)*omega*(2*sy*T_old(j+1,i,t+1)-2*sy*dy*c2+ ...
                                +2*sx*T_old(j,i-1,t+1)+T_old(j,i,t)+dt*Q(j,i,t+1)/rho/cp) +  (1-omega)*T_old(j,i,t+1);
            else

              T_new(j,i,t+1) = 1/(1+2*sx+2*sy+2*sy*dy*c1)*omega*(2*sy*T_old(j+1,i,t+1)-2*sy*dy*c2+ ...
                                sx*T_old(j,i-1,t+1)+sx*T_old(j,i+1,t+1)+T_old(j,i,t)+dt*Q(j,i,t+1)/rho/cp) +  (1-omega)*T_old(j,i,t+1);
            end

      elseif j == j_max

          if i==1

                T_new(j,i,t+1) = 1/(1+2*sx+2*sy-2*sy*dy*c1)*omega*(2*sy*T_old(j-1,i,t+1)+2*sy*dy*c2+ ...
                                +2*sx*T_old(j,i+1,t+1)+T_old(j,i,t)+dt*Q(j,i,t+1)/rho/cp) +  (1-omega)*T_old(j,i,t+1);

            elseif i==i_max

                T_new(j,i,t+1) = 1/(1+2*sx+2*sy-2*sy*dy*c1)*omega*(2*sy*T_old(j-1,i,t+1)+2*sy*dy*c2+ ...
                                +2*sx*T_old(j,i-1,t+1)+T_old(j,i,t)+dt*Q(j,i,t+1)/rho/cp) +  (1-omega)*T_old(j,i,t+1);

            else

              T_new(j,i,t+1) = 1/(1+2*sx+2*sy-2*sy*dy*c1)*omega*(2*sy*T_old(j-1,i,t+1)+2*sy*dy*c2+ ...
                                sx*T_old(j,i-1,t+1)+sx*T_old(j,i+1,t+1)+T_old(j,i,t)+dt*Q(j,i,t+1)/rho/cp) +  (1-omega)*T_old(j,i,t+1);

          end

      else

        if (i==1) 

          T_new(j,i,t+1) = 1/(1+2*sx+2*sy)*omega*(sy*T_old(j+1,i,t+1)+sy*T_old(j-1,i,t+1)+ ...
                                 +2*sx*T_old(j,i+1,t+1)+T_old(j,i,t)+dt*Q(j,i,t+1)/rho/cp) +  (1-omega)*T_old(j,i,t+1);

      elseif (i==i_max) 

          T_new(j,i,t+1) = 1/(1+2*sx+2*sy)*omega*(sy*T_old(j+1,i,t+1)+sy*T_old(j-1,i,t+1)+ ...
                                 +2*sx*T_old(j,i-1,t+1)+T_old(j,i,t)+dt*Q(j,i,t+1)/rho/cp) +  (1-omega)*T_old(j,i,t+1);   
      else

      T_new(j,i,t+1) = 1/(1+2*sx+2*sy)*omega*(sy*T_old(j+1,i,t+1)+sy*T_old(j-1,i,t+1)+ ...
                                 +sx*T_old(j,i-1,t+1)+sx*T_old(j,i+1,t+1)+T_old(j,i,t)+dt*Q(j,i,t+1)/rho/cp) +  (1-omega)*T_old(j,i,t+1);

      end   

      end

     %--------------------------------------------------------------------------------%

      temp = temp + abs((T_new(j,i,t+1)-T_old(j,i,t+1))/T_new(j,i,t+1)); % Relative Error

      temp = temp + abs((T_new(j,i,t+1)-T_old(j,i,t+1))); % Relative Error

      T_old(j,i,t+1) = T_new(j,i,t+1);

   end

end

   T_old = T_new;   %%!!%%Reset T_old 

change(iter) = temp/(i_max-2)/(j_max-2);      % the average relative error of all interior points

if( change(iter) < change_want ) 

      fprintf('\n The total number of iterations excuted is %g \n',iter)

    disp('Desired accuracy achieved; breaking out of main loop');

    break;

end

end

    x=linspace(0,Lx,i_max);

    y=linspace(0,Ly,j_max);

     Z = flipud(T_new(:,:,t+1));

    surf(x,y,Z)   

     eval(['print -djpeg heat2d_' num2str(t) '.jpeg']);

end

[Mohamed Hussein]

Thanks for your answer.

The first link is quite good for solving similar problems. However, I have to apply the finite difference formulation in my code. Also, I need to subject the plate to derivative boundary conditions at all four edges ( convective at the top and the remaining are zero derivative).

My code was working for constant temperature boundary conditions, but when I applied the derivative boundary conditions it started to give me strange results.

[Sterling Baird]

The SOR algorithm looks really messy. Why all the if statements?

Answer 1

[John BG]

have you seen these links?

1.

https://uk.mathworks.com/help/pde/examples/inhomogeneous-heat-equation-on-a-square-domain.html?s_tid=srchtitle

2.

https://uk.mathworks.com/matlabcentral/fileexchange/40096-numerical-solutions-of-heat-equation?s_tid=srchtitle

3.

https://uk.mathworks.com/matlabcentral/fileexchange/59916-simple-heat-equation-solver?s_tid=srchtitle

there are more functions available in the file exchange.

If you find these lines useful would you please mark my answer as Accepted Answer?

To any other reader, if you find this answer of any help please click on the thumbs-up vote link,

thanks in advance for time and attention

John BG

<mailto:[e-mail] [e-mail]>

Answer 2

[AKASH BHOYAR]

@John BG @Mohamed Hussein How to solve when source term is having fixed intensity and conductivity varies as function of temp?