Solving Steady State 2D Heat Conduction Equation using C++ OOP

Objective :

The challenge is to solve the 2D Steady Heat Conduction equation using C++, 

where alpha is the thermal diffusivity

About equation :

In physics and mathematics, the heat equation is a partial differential equation that describes how the distribution of some quantity (such as heat) evolves over time in a solid medium, as it spontaneously flows from places where it is higher toward places where it is lower. It is a special case of the diffusion equation.

Abstract: This challenge will focus on solving the equation with three solvers jacobian, gauss-seidel ,and SOR iteration for two cases steady and transient.

While we use iterative solvers to solve steady-state PDE and both explicit and implicit methods are used to solve transient PDE.

Given,

Generating a 2D  profile 

I have generated a 2D profile of temperature which gets analyzed. The temperature at 4 boundaries are given as below

Top=400K

Bottom=900k

Left=600K

Right=800K

These values are given as initial and boundary conditions for the domain .

Solving 2D heat conduction equation of steady state: Here the temperature doesn't change with respect to time  

Steps:

Discretization: The discretized equation of the steady-state is

om further simplification we get

where,

 As we don't know the value sat some nodes like T( i+1) and T(j+1) we take guess values and we get a set of coupled equations that are to be solved at once,as when there

are guess values we use iterative solvers and those values converge under solver and we get a converged solution. we will say that our solution is converged when new iteration values and old values are same.

Here we use two types of iterative solvers, 

1) Jacobian Iterative Solver :

In the Jacobi method, each iteration is done using a set of previous values, even if new values are available in the middle of an iteration. Like when you are solving for n+1 iteration and for i th node you will use values of previous(n) iteration even if you have the updated values of the nodes in the current iteration.so this method is called the method of simultaneous displacements.

The 2D steady state heat conduction equation for this method is 

2) Gauss-Seidel iterative solver :

 In Gauss-Seidel, as soon as you find a new value in the current iteration you will use it for finding other values in the current iteration.

Like when you are solving for n+1 iteration and for i th node you will use all the available values up to i-1 node of current iteration(n+1) and use i+1 node values and above values of previous iteration(n).So we get more converged values. so this method is called the method of spontaneous displacements.

The 2D steady state heat conduction equation for this method is 

Input Variable: 

1. Nodes: we are going to take no. of nodes in X-direction,nx, and in Y-direction,ny=10
2. Length of the domain: This is taken as 1 along the x-axis and y-axis.
3. Given, Tolerance is assumed to be 1e-7.
4. And we also assumed the error as 999. we can assume any value but it must be greater than the required tolerance.

By taking these as inputs we create a function for each solver and use those functions in the main program and solve the 2D steady-state equation on a square domain using C++.

Main Program :

// Solving the steady state 2D heat conduction problem using iterative solvers over a 
 //square domain.

 //main function

#include"tempheader.h"

//#include"temperature.h"
 int main() 
 { 
  temperature tempm  ;
tempm.tempinit() ;
tempm.L =1 ;
tempm.tolerance=1e-7 ;
tempm.j_iteration_number = 1 ;
tempm.g_iteration_number = 1 ;
tempm.error = 999 ; 
int n ;

// vector declared to recieve final temperature values from the function steady state jacobi
vector temperature_j  ;

//vector declared to receive the final temperature values from the function  steady state gauss - seidel
vector temperature_g  ;


//Asking the user for input for selecting the solver required
cout<<"Enter the solver number 1 : jacobi , 2 :gauss-seidel :  " ;
 cin>>n  ;
    if (n == 1)
    {

      // calling the function steady state jacobi
   temperature_j = tempm.steady_jacobian_method() ;

  //calling the function for printing the values to the console
     tempm.print( temperature_j ) ;

      cout<<"jacobi iteration count : "<

Program for temperature class header file :

//Header file for class declaration //tempheader.h

#pragma once
#include 
#include 
#include
using namespace std;
 #include 
 #include <bits/stdc++.h> 
 
class temperature   
{
public :

double L  ;
double nx   ; // rows
double ny   ;  // columns

double tolerance ;
int j_iteration_number ;
int g_iteration_number ;
double error ; 


//function for initialzing __ _variables
void tempinit() ;
 
vector steady_jacobian_method() ;

vector steady_gauss_method() ;

vector tempinitialization( vector &temp) ;

void print (vector temp) ;
 


} ;

 

 Program for definition of class - temperature , member functions :

  //program for the definition of class temperature member functions

#include"tempheader.h"  //including header file for declaration of class temperature 
  
// function for initializing variables
void temperature:: tempinit()
{
 cout<<"Enter no of grid points  in x direction:  " ;
 cin>>nx ;
cout<<"Enter no of  grid points  in y direction:  " ;
 cin>>ny ;

}

//function definition for printing values of temperature outside class

void temperature::print (vector temp)
{

for(int i = 0; i < nx; i++)
    {
        for(int j = 0; j < ny; j++)
        {
            cout <<std::setw(8)<<std::setiosflags(ios::left)<<temp[i][j] << "  ";
        }
        cout<< endl;
    }
} ;

//function definition outside of the class for jacobi method

vector temperature::steady_jacobian_method()

{    
 
    double dx=L/(nx-1);
double dy=L/(ny-1);
// TEMPERATURE matrix

/*
    Create a vector containing "nxf elements" (rows )
     and each vectors each of size "nyf".( columns) and initializing to 0 ,
    */ 
 
  vector temp( nx  , vector (ny , 0)); 
  
    
temp = tempinitialization(temp)  ; // initialized values
     
   
double k=(2*(pow(dx,2)+pow(dy,2))/(pow(dx,2)*pow(dy,2)));

  vector  Told=temp;

vector resultsub( nx  , vector (ny , 0)); 
   

j_iteration_number = 1 ;
// convergence loop
 while (error>tolerance)
 { 
 //nodal loop
  for(int i = 1; i < nx-1; i++)
 { 
 for(int j = 1; j < ny-1; j++)
{ 
    

temp[i][j] =(1/k)*((((Told[i-1][j]+Told[i+1][j])/pow(dx,2))+ ((Told[i][j-1]+Told[i][j+1]))/pow(dy,2))) ;


}
  
 
 }
  
  // tnew -told 
  
  for(int i = 1; i < nx-1; i++)
    {
        for(int j = 1; j < ny-1; j++)  
    {   
            resultsub[i][j] = temp[i][j] - Told[i][j];
    }


    }

//abs
 
  for(int i = 1; i < nx-1; i++)
    {
        for(int j = 1; j < ny-1; j++)  
    {   
            resultsub[i][j] = abs(resultsub[i][j]);
    }
    }
  
   double  tem  = 0  ;
   for(int i = 1; i < nx-1; i++)
    {
        for(int j = 1; j < ny-1; j++)  

            if(tem < resultsub[i][j])
            {
                tem  = resultsub[i][j];
            }
    }

error =   tem ; 

//%updating old values
Told =temp ;
 j_iteration_number=  j_iteration_number+1;


}
 
 
return temp ;
}

//function definition outside of the class for gauss seidel method

vector temperature::steady_gauss_method()

{    
 
    double dx=L/(nx-1);
double dy=L/(ny-1);
// TEMPERATURE matrix

/*
    Create a vector containing "nxf elements" (rows )
     and each vectors each of size "nyf".( columns) and initializing to 0 ,
    */ 
 
  vector temp( nx  , vector (ny , 0)); 
    
temp = tempinitialization(temp)  ; // initialized values
     
   
double k=(2*(pow(dx,2)+pow(dy,2))/(pow(dx,2)*pow(dy,2)));

  vector  Told=temp;

vector resultsub( nx  , vector (ny , 0)); 
   

g_iteration_number = 1 ;
// convergence loop
 while (error>tolerance)
 { 
 //nodal loop
  for(int i = 1; i < nx-1; i++)
 { 
 for(int j = 1; j < ny-1; j++)
{ 
    

temp[i][j] =(1/k)*((((temp[i-1][j]+Told[i+1][j])/pow(dx,2))+ ((temp[i][j-1]+Told[i][j+1]))/pow(dy,2))) ;


}
  
 
 }
  
  // tnew -told 
  
  for(int i = 1; i < nx-1; i++)
    {
        for(int j = 1; j < ny-1; j++)  
    {   
            resultsub[i][j] = temp[i][j] - Told[i][j];
    }


    }

//abs
 
  for(int i = 1; i < nx-1; i++)
    {
        for(int j = 1; j < ny-1; j++)  
    {   
            resultsub[i][j] = abs(resultsub[i][j]);
    }
    }
       //finding max
       
   double tem  = 0  ;
   for(int i = 1; i < nx-1; i++)
    {
        for(int j = 1; j < ny-1; j++)  

            if(tem < resultsub[i][j])
            {
                tem  = resultsub[i][j];
            }
    }

error =   tem ; 

//%updating old values
Told =temp ;
 g_iteration_number=  g_iteration_number+1;


}
 
 
return temp ;
}

//function definition outside of the class for temperature initialization over domain //boundary conditions and initial values

vector temperature::tempinitialization(vector &temp) 
{
 
      for(int i = 0; i < 1; i++)
    {
        for(int j = 0; j < ny; j++)
        {
            temp[i][j] =  400;
        }
    }
      

      for(int i = 0; i < nx; i++)
    {
        for(int j = 0; j < 1; j++)
        {
            temp[i][j] =  600;
        }
    }

    for(int i = 0; i < nx; i++)
    {
        for(int j = ny-1; j <ny; j++)
        {
            temp[i][j] =  200;
        }
    }

    for(int i = nx-1; i < nx; i++)
    {
        for(int j = 0; j < ny; j++)
        {
            temp[i][j] =  800;
        }
    }

 
      return temp ;
 }

   

Compiling the code  :

Output from the solvers  :

1) Jacobi method :

2) Gauss seidel method  :

In the above two results, you can observe that the number of iterations for Jacobi iteration method is high i.e ,323 iterations and then for the gauss seidel method with the number of iterations are 169 .