CFD || Solving Transient State 2D Heat Conduction Equation

OBJECTIVE

To solve a Transient 2D Heat conduction equation using different Iterative techniques (Jacobi,Gauss Seidal,SOR) and also carry out stability analysis in the transient problem.

ASSUMPTION

• Assume that the domain is a unit square.
• Assume nx = ny [Number of points along the x direction is equal to the number of points along the y direction]
• Boundary conditions for steady case - Left:400K, Top:600K, Right:800K, Bottom:900K

OVERVIEW

We know 2D Transient State Heat Conduction Eqn as

`(partialT)/(partialt)-alpha[(partialT)/(partialx^2)+(partialT)/(partialy^2)]=0`

Let,

• `T_p = T(i,j);`
• `T_L = T(i-1,j);`
• `T_R = T(i+1,j);`
• `T_T = T(i,j+1);`
• `T_B = T(i,j-1);`

Discretizing the above equation

`(T_p^(n+1)-T_p^n)/(Deltat)- alpha[(T_R-2T_p+T_L)/(Deltax^2)+(T_T-2T_p+T_B)/(Deltay^2)]^n=0`

`impliesT_p^(n+1)= T_p^n+alphaDeltat[(T_R-2T_p+T_L)/(Deltax^2)+(T_T-2T_p+T_B)/(Deltay^2)]^n` ......(Explicit Eqn)

Similarly;

`T_p^(n+1)= T_p^n+alphaDeltat[(T_R-2T_p+T_L)/(Deltax^2)+(T_T-2T_p+T_B)/(Deltay^2)]^(n+1)`......(Implicit Eqn)

To above implicit equation apply the Jacobi,Gauss seidal, SOR algorithm as shown in Steady State Simulation previously

PROGRAM

%solving for 2D Transient State condition
%dt/dt+(alpha*((dt^2/dx^2)+(dt^2/dy^2))
clear all
close all
clc

nx=10;
ny=10;
T=ones(nx,ny);

%BC
T(1,:)=900;    %Temp bottom
T(:,nx)=800;  %Temp right
T(ny,:)=600;   %Temp top
T(:,1)=400;   %Temp Left

%Intialisation
Told=T;
T_prev=T;
x=linspace(0,1,nx);
dx=x(2)-x(1);
y=linspace(0,1,ny);
dy=y(2)-y(1);

tol=1e-3;
error=10;
tol_stdy=0.01;
error_stdy=10;
alpha=1e-4;
dt=5;
k1=(alpha*dt)/(dx^2);
k2=(alpha*dt)/(dy^2);
cfl=alpha*(dt/dx^2);

solver=4;    %%%%(Choose the solver)
%_______________________________________________________________________



%_____Explicit method__________%
if (solver==1)    
   tic
    for k=1:1000        
        if(error_stdy>tol_stdy)        
            for i=2:nx-1
                for j=2:ny-1
                    T(i,j)=Told(i,j)+(k1*(Told(i+1,j)-(2*Told(i,j))+Told(i-1,j)))+(k2*(Told(i,j+1)-(2*Told(i,j))+Told(i,j-1)));
                end
            end
        Told=T;
        error_stdy=max(max(abs(T_prev-T)));
        T_prev=T;
        
        %plotting
       
        [x1,y1]=contourf(x,y,T);
        colormap(jet);
        clabel(x1,y1);
        xlabel('xaxis')
        ylabel('yaxis')
        text_title=sprintf('Temperature at time=%d Explicit',k); 
        title(text_title)
        pause(0.0001)
        
        end    
    end
    Time_Explicit=toc;
end

%________________________________________________________________________


%_________Jacobi method__________%
if (solver==2)    
    tic
    j_iter=1;
    term1=(1+2*k1+2*k2);
    term2=k1/term1;
    term3=k2/term1;
    for k=1:1000
        if(error_stdy>tol_stdy)     
            error=10;
            while(error>tol)
                for i=2:nx-1
                    for j=2:ny-1
                        H=(Told(i-1,j)+Told(i+1,j));
                        V=(Told(i,j-1)+Told(i,j+1));
                        T(i,j)=(T_prev(i,j)/term1)+(term2*H)+(term3*V);
                    end
                end
                    error=max(max(abs(T-Told)));
                    Told=T;
                    j_iter=j_iter+1;
            end
        error_stdy = max(max(abs(T_prev-T)));
        T_prev=T;
                
        
        %plotting
        
        [x1,y1]=contourf(x,y,T);
        colormap(jet);
        clabel(x1,y1);
        xlabel('xaxis')
        ylabel('yaxis')
        text_title=sprintf('Temperature at time step-%d Jacobi Iteration=%d',k,j_iter); 
        title(text_title)
        pause(0.0001)
       
        end
    end
     Total_Jacobi_Iter=j_iter;
     Time_Jacobi=toc;
end

%___________________________________________________________________________


%__________GS method___________%
if (solver==3)
    
    tic
    gs_iter=1;
    term1=(1+2*k1+2*k2);
    term2=k1/term1;
    term3=k2/term1;
    
    for k=1:1000
        if(error_stdy>tol_stdy)                
        error=10;        
            while(error>tol)
                for i=2:nx-1
                    for j=2:ny-1
                    H=(T(i-1,j)+Told(i+1,j));
                    V=(T(i,j-1)+Told(i,j+1));
                    T(i,j)=(T_prev(i,j)/term1)+(term2*H)+(term3*V);
                    end
                end
                    error=max(max(abs(T-Told)));
                    Told=T;
                    gs_iter=gs_iter+1;
            end
        error_stdy = max(max(abs(T_prev-T)));
        T_prev=T;        
        
        
        %plotting
        
        [x1,y1]=contourf(x,y,T);
        colormap(jet);
        clabel(x1,y1);
        xlabel('xaxis')
        ylabel('yaxis')
        text_title=sprintf('Temperature at time step-%d GS Iteration=%d',k,gs_iter); 
        title(text_title)
        pause(0.0001)
        end
    end
    Total_GS_Iter=gs_iter;
    Time_GS=toc;
end

%________________________________________________________________________    

%__________SOR METHOD__________%
if (solver==4)
    tic
    w=1.4;
    term1=(1+2*k1+2*k2);
    term2=k1/term1;
    term3=k2/term1;
    sor_iter=1;
   
    for k=1:1000
        if(error_stdy>tol_stdy)        
        error=10;        
            while(error>tol)
                for i=2:nx-1
                    for j=2:ny-1
                        H=(T(i-1,j)+Told(i+1,j));
                        V=(T(i,j-1)+Told(i,j+1));
                        T(i,j)=(w*((T_prev(i,j)/term1)+(term2*H)+(term3*V)))+((1-w)*Told(i,j));
                    end
                end
                    error=max(max(abs(T-Told)));
                    Told=T;
                    sor_iter=sor_iter+1;
            end
        error_stdy = max(max(abs(T_prev-T)));
        T_prev=T;        
        
        %plotting
        
        [x1,y1]=contourf(x,y,T);
        colormap(jet);
        clabel(x1,y1);
        xlabel('xaxis')
        ylabel('yaxis')
        text_title=sprintf('Temperature at time step-%d SOR Iteration=%d',k,sor_iter); 
        title(text_title)
        pause(0.0001)
        end
    end
    
     Total_SOR_Iter=sor_iter;
     Time_SOR=toc;
end


    

RESULTS

Time taken to reach steady state using Explicit Method = 101.62s

Time taken to reach steady state using Implicit Method(Jacobi) = 116.57s

Time taken to reach steady state using Implicit Method(Gauss-seidal) = 114.88s

Time taken to reach steady state using Implicit Method(SOR) = 96.78s

• Unstability

The above fig represent a unstability plot for transient problem solved using explicit method

CONCLUSION

1. The time taken to reach the steady state are in the order SOR<Explicit
2. The implicit function is unconditionally stable no matter what CFL number.
3. The explicit function is conditionally stable for CFL<0.28 (In this case).

Note: `CFL = alpha(dt)/(dx^2)`