Child pages
  • Heat Transfer - 2D
Skip to end of metadata
Go to start of metadata
#include "SCmathlib.h"
#include <iostream>
#include "SCchapter9.h"

int nx = 10;
int ny = 10;

int Ind(int i,int j)
{
    return i*(ny-1)+j;
};

int main()
{        
    int n = (nx-1)*(ny-1);    
        
    SCVector b(n);
    SCMatrix A(n,n);
    
    double Hx = 20;
    double Hy = 20;    
  
    double dx = Hx/nx;
    double dy = Hy/ny;

    double Tx0=10;
    double Tx1=10;
    double Ty0=10;
    double Ty1=10;

    double q=10.0;
    
    for(int i=0;i<nx-1;i++)
    {
        for(int j=0;j<ny-1;j++)
        {
            if((i==0)&&(j==0))
            {
                //corner
                // du^2/dx^2
                int row = Ind(i,j);
                A(row,Ind(i+1,j)) = 1*dy*dy; A(row,Ind(i,j)) += -2*dy*dy; 
                //A(row,Ind(i-1,j)) = 1*dy*dy;
                // du^2/dy^2
                A(row,Ind(i,j+1)) = 1*dx*dx; A(row,Ind(i,j)) += -2*dx*dx; 
                //A(row,Ind(i,j-1)) = 1*dx*dx;
                b(row) += (-1)*1*dx*dx*Tx0 + (-1)*1*dy*dy*Ty0;
            }
            else if((i==nx-2)&&(j==0))
            {
                //corner
                // du^2/dx^2
                int row = Ind(i,j);
                //A(row,Ind(i+1,j)) = 1*dy*dy; 
                A(row,Ind(i,j)) += -2*dy*dy; A(row,Ind(i-1,j)) = 1*dy*dy;
                // du^2/dy^2
                A(row,Ind(i,j+1)) = 1*dx*dx; A(row,Ind(i,j)) += -2*dx*dx; 
                //A(row,Ind(i,j-1)) = 1*dx*dx;            
                b(row) += (-1)*1*dx*dx*Tx1 + (-1)*1*dy*dy*Ty0;
            }
            else if((i==nx-2)&&(j==ny-2))
            {
                //corner
                // du^2/dx^2
                int row = Ind(i,j);
                A(row,Ind(i,j)) += -2*dy*dy; A(row,Ind(i-1,j)) = 1*dy*dy;
                //A(row,Ind(i+1,j)) = 1*dy*dy;
                // du^2/dy^2
                A(row,Ind(i,j)) += -2*dx*dx; A(row,Ind(i,j-1)) = 1*dx*dx;            
                //A(row,Ind(i,j+1)) = 1*dx*dx;
                b(row) += (-1)*1*dx*dx*Tx1 + (-1)*1*dy*dy*Ty1;
            }
            else if((i==0)&&(j==ny-2))
            {
                //corner
                // du^2/dx^2
                int row = Ind(i,j);
                A(row,Ind(i+1,j)) = 1*dy*dy; A(row,Ind(i,j)) += -2*dy*dy; 
                //A(row,Ind(i-1,j)) = 1*dy*dy;
                // du^2/dy^2
                A(row,Ind(i,j)) += -2*dx*dx; A(row,Ind(i,j-1)) = 1*dx*dx;            
                //A(row,Ind(i,j+1)) = 1*dx*dx;
                b(row) += (-1)*1*dx*dx*Tx0 + (-1)*1*dy*dy*Ty1;
            }
            else if(i==0)
            {
                //side 
                // du^2/dx^2
                int row = Ind(i,j);
                A(row,Ind(i+1,j)) = 1*dy*dy; A(row,Ind(i,j)) += -2*dy*dy;
                //A(row,Ind(i-1,j)) = 1*dy*dy;
                // du^2/dy^2
                A(row,Ind(i,j+1)) = 1*dx*dx; A(row,Ind(i,j)) += -2*dx*dx; A(row,Ind(i,j-1)) = 1*dx*dx;
                b(row) += (-1)*1*dy*dy*Tx0;             
            }
            else if( i == nx-2 )
            {
                //side             
                // du^2/dx^2
                int row = Ind(i,j);
                A(row,Ind(i,j)) += -2*dy*dy; A(row,Ind(i-1,j)) = 1*dy*dy;
                //A(row,Ind(i+1,j)) = 1*dy*dy; 
                // du^2/dy^2
                A(row,Ind(i,j+1)) = 1*dx*dx; A(row,Ind(i,j)) += -2*dx*dx; A(row,Ind(i,j-1)) = 1*dx*dx;            
                b(row) += (-1)*1*dy*dy*Tx1;             
            }
            else if(j==0)
            {
                //side             
                // du^2/dx^2
                int row = Ind(i,j);
                A(row,Ind(i+1,j)) = 1*dy*dy; A(row,Ind(i,j)) += -2*dy*dy; A(row,Ind(i-1,j)) = 1*dy*dy;
                // du^2/dy^2
                A(row,Ind(i,j+1)) = 1*dx*dx; A(row,Ind(i,j)) += -2*dx*dx; 
                //A(row,Ind(i,j-1)) = 1*dx*dx;
                b(row) += (-1)*1*dx*dx*Ty0;              
            }
            else if(j==ny-2)
            {
                //side             
                // du^2/dx^2
                int row = Ind(i,j);
                A(row,Ind(i+1,j)) = 1*dy*dy; A(row,Ind(i,j)) += -2*dy*dy; A(row,Ind(i-1,j)) = 1*dy*dy;
                // du^2/dy^2
                A(row,Ind(i,j)) += -2*dx*dx; A(row,Ind(i,j-1)) = 1*dx*dx;  
                //A(row,Ind(i,j+1)) = 1*dx*dx; 
                b(row)=(-1)*1*dx*dx*Ty1;
            }
            else 
            {
                //interior point             
                // du^2/dx^2
                int row = Ind(i,j);
                A(row,Ind(i+1,j)) = 1*dy*dy; A(row,Ind(i,j)) += -2*dy*dy; A(row,Ind(i-1,j)) = 1*dy*dy;
                // du^2/dy^2
                A(row,Ind(i,j+1)) = 1*dx*dx; A(row,Ind(i,j)) += -2*dx*dx; A(row,Ind(i,j-1)) = 1*dx*dx;
            }            
        };
    };
    /* heat source */
    int ix = nx/2-1;
    int iy = ny/2-1; 
    int row = Ind(ix,iy);
    b(row) += -q*dx*dx*dy*dy;   
          
    int info = 0; 
    GaussElimination(A,b,info);
    
    for(int j=0;j<=ny;j++)
    {        
        std::cout<<" "<<std::scientific<<Tx0;
    }
    std::cout<<"\n";
    for(int i=0;i<nx-1;i++)
    {
        std::cout<<" "<<std::scientific<<Ty0;    
        for(int j=0;j<ny-1;j++)
        {
            std::cout<<" "<<std::scientific<<b(Ind(i,j));
        }
        std::cout<<" "<<std::scientific<<Ty1<<"\n";
    }
    for(int j=0;j<=ny;j++)
    {
        std::cout<<" "<<std::scientific<<Tx1;
    };
        
    return 0;
};
  • No labels