CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Wiki > Source code archive - educational

Source code archive - educational

From CFD-Wiki

(Difference between revisions)
Jump to: navigation, search
m (One dimensional stedy state conduction with Heat generation code using C)
(Removed example code; putting it in separate page.)
Line 1: Line 1:
-
#include <stdio.h>
+
* [[One dimensional steady state conduction with heat generation | One dimensional steady state conduction with heat generation code in C]]
-
#include <stdlib.h>
+
-
 
+
-
#define MAXROWS 10000
+
-
 
+
-
double **GetMatrix(double, double, double, double, double, int);
+
-
void squareoutput(double **, int);
+
-
void nonsquareoutput(double **, int, int);
+
-
double **triangularize(double *[MAXROWS], int, int);
+
-
double *returnsolvector(double **c, int nrows);
+
-
void writeoutputvector(double *p, int nrows);
+
-
 
+
-
 
+
-
int main(int argc, char **argv)
+
-
{
+
-
+
-
/*
+
-
Variables chosen:
+
-
GAMMA:  Thermal Conductivity
+
-
L : Length of the domain
+
-
Dx : Node spacing
+
-
*/
+
-
 
+
-
double GAMMA, L, Dx, Ta , Tb, Area;
+
-
int N, row;
+
-
double **temp;
+
-
double **aug;
+
-
double *solvector;
+
-
 
+
-
+
-
printf("\nThis program works for only simple 1-D conduction");
+
-
printf("Is thermal Conductivity constant?");
+
-
+
-
+
-
printf("\nEnter the Value for Heat conductivity");
+
-
scanf("%lf", &GAMMA);
+
-
printf("\nEnter the Length of the computational domain");
+
-
scanf("%lf", &L);
+
-
printf("\nEnter the cross-sectional area of the computational domain");
+
-
scanf("%lf", &Area);
+
-
printf("\nEnter the temperature at the ends of the domain");
+
-
scanf("%lf %lf", &Ta, &Tb);
+
-
printf("\nEnter the number of Nodes u want to be chosen for analysis");
+
-
scanf("%d", &N);
+
-
 
+
-
temp=(double **) malloc(N*sizeof(double *));
+
-
    for(row=0;row<N;row++)
+
-
    temp[row]=(double *) malloc(N*sizeof(double));
+
-
 
+
-
aug= (double **) malloc ((N+1)*sizeof(double));
+
-
    for(row=0;row<N;row++)
+
-
    aug[row]=(double *) malloc((N+1)*sizeof(double));
+
-
 
+
-
solvector= (double *) malloc (N*sizeof(double));
+
-
+
-
Dx=L/(N);
+
-
 
+
-
printf("\nThe Node spacing is: %lf", Dx);
+
-
 
+
-
aug=GetMatrix(GAMMA, Dx, Ta, Tb, Area, N);
+
-
aug=triangularize(aug, N, N+1);
+
-
nonsquareoutput(aug, N, N+1);
+
-
printf("\n\n");
+
-
solvector=returnsolvector(aug,N);
+
-
writeoutputvector(solvector, N);
+
-
printf("\n\n\n");
+
-
return 0;
+
-
}
+
-
 
+
-
 
+
-
double **GetMatrix(double GAMMA, double Dx, double Ta, double Tb, double Area, int N)
+
-
{
+
-
int row, col;
+
-
double DxW, DxE;
+
-
double Aw, Ae, Ap;
+
-
+
-
/* Memory Allocation for the matrix */
+
-
 
+
-
double **matrix; double *solvector;
+
-
double **aug;
+
-
 
+
-
solvector= (double *) malloc (N*sizeof(double));
+
-
 
+
-
matrix=(double **) malloc(N*sizeof(double));
+
-
    for(row=0;row<N;row++)
+
-
    matrix[row]=(double *) malloc(N*sizeof(double));
+
-
 
+
-
aug= (double **) malloc ((N+1)*sizeof(double));
+
-
    for(row=0;row<N;row++)
+
-
    aug[row]=(double *) malloc((N+1)*sizeof(double));
+
-
 
+
-
 
+
-
/* Uniform Mesh has been chosen */
+
-
 
+
-
DxW=DxE=Dx;
+
-
+
-
Aw=(GAMMA/DxW)*Area;
+
-
Ae=(GAMMA/DxE)*Area;
+
-
 
+
-
/* Generation of the Matrix */
+
-
+
-
row=0;
+
-
for(col=0; col<N; col++)
+
-
{
+
-
if(col==row)
+
-
*(*(matrix+row)+col)= (Ae+2*Aw);
+
-
else if(col==row+1)
+
-
*(*(matrix+row)+col)= -Ae ;
+
-
else
+
-
*(*(matrix+row)+col)=0;
+
-
}
+
-
 
+
-
 
+
-
for (row=1; row<N-1; row++)
+
-
{
+
-
for(col=0; col<N; col++)
+
-
{
+
-
if(col==row+1)
+
-
*(*(matrix+row)+col)= -Aw;
+
-
else if (col==row-1)
+
-
*(*(matrix+row)+col)= -Ae;
+
-
else if (col==row)
+
-
*(*(matrix+row)+col)= (Aw+Ae);
+
-
else
+
-
*(*(matrix+row)+col)=0;
+
-
}
+
-
}
+
-
+
-
+
-
row=N-1;
+
-
for(col=0; col<N; col++)
+
-
{
+
-
if(col==row)
+
-
*(*(matrix+row)+col)= (2*Ae+Aw);
+
-
else if(col==row-1)
+
-
*(*(matrix+row)+col)=-Aw;
+
-
else
+
-
*(*(matrix+row)+col)=0;
+
-
}
+
-
 
+
-
/* Getting the solution vector */
+
-
 
+
-
*(solvector+0)= 2*Aw*Ta;
+
-
 
+
-
for (row=1; row<=N-2; row++)
+
-
*(solvector+row)=0;
+
-
 
+
-
*(solvector+N-1)= 2*Ae*Tb;
+
-
 
+
-
 
+
-
/* Getting the augmented matrix */
+
-
 
+
-
for (row=0; row<N; row++)
+
-
{
+
-
for(col=0; col<N; col++)
+
-
{
+
-
*(*(aug+row)+col) = *(*(matrix+row)+col);
+
-
}
+
-
}
+
-
+
-
for (row=0; row<N; row++)
+
-
*(*(aug+row)+N)= *(solvector+row);
+
-
 
+
-
 
+
-
return(aug);
+
-
}
+
-
+
-
 
+
-
 
+
-
 
+
-
 
+
-
 
+
-
 
+
-
 
+
-
 
+
-
 
+
-
 
+
-
/* Displaying the matrix */
+
-
 
+
-
void squareoutput(double **c, int nrows)
+
-
{
+
-
int row, col;
+
-
for(row=0; row<nrows; row++)
+
-
{
+
-
printf("\n");printf("\n");
+
-
for (col=0; col<nrows; col++)
+
-
{
+
-
printf("%8.2lf", *(c[row]+col));
+
-
}
+
-
}
+
-
return;
+
-
}
+
-
 
+
-
void nonsquareoutput(double **c, int nrows, int ncols)
+
-
 
+
-
{
+
-
int row, col;
+
-
for(row=0; row<nrows; row++)
+
-
{
+
-
printf("\n");printf("\n");
+
-
for (col=0; col<ncols; col++)
+
-
{
+
-
printf("%12.2lf", *(c[row]+col));
+
-
}
+
-
}
+
-
return;
+
-
}
+
-
 
+
-
double **triangularize(double *a[MAXROWS], int nrows, int ncols)
+
-
{
+
-
int row, col, k, temp=0;
+
-
double coef;
+
-
double **tri;
+
-
tri=(double **)malloc(ncols*sizeof(double *));
+
-
    for(row=0;row<nrows;row++)
+
-
    tri[row]=(double*)malloc(ncols*sizeof(double));
+
-
 
+
-
 
+
-
{
+
-
int row, col;
+
-
for (row=0; row<nrows; row++)
+
-
for(col=0; col<ncols; col++)
+
-
*(tri[row]+col)=*(a[row]+col);
+
-
}
+
-
 
+
-
+
-
+
-
+
-
for(temp=0;temp<nrows-1; temp++)
+
-
{
+
-
if ((tri[temp]+temp)==0)continue;
+
-
else
+
-
{
+
-
for(row=temp; row<nrows; row++)
+
-
{
+
-
for(k=row+1; k< nrows; k++)
+
-
{
+
-
if( *(tri[k]+temp)==0)
+
-
continue;
+
-
else
+
-
coef= (*(tri[k]+temp))/(*(tri[temp]+temp));
+
-
for(col=0; col<ncols; col++)
+
-
*(tri[k]+col)=*(tri[k]+col)-(coef* *(tri[row]+col));
+
-
}
+
-
}
+
-
}
+
-
}
+
-
+
-
for(temp=nrows-1;temp>0; temp--)
+
-
{
+
-
if ((tri[temp]+temp)==0)continue;
+
-
else
+
-
{
+
-
for(row=temp; row>0; row--)
+
-
{
+
-
for(k=row-1; k>=0; k--)
+
-
{
+
-
if( *(tri[k]+temp)==0)
+
-
continue;
+
-
else
+
-
coef= (*(tri[k]+temp))/(*(tri[temp]+temp));
+
-
for(col=0; col<ncols; col++)
+
-
*(tri[k]+col)=*(tri[k]+col)-(coef* *(tri[row]+col));
+
-
}
+
-
}
+
-
}
+
-
}
+
-
 
+
-
return(tri);
+
-
}
+
-
 
+
-
double *returnsolvector(double **c, int nrows)
+
-
{
+
-
int temp;
+
-
double *solvector;
+
-
solvector= (double *) malloc (nrows*sizeof(double));
+
-
 
+
-
for(temp=0; temp<nrows; temp++)
+
-
{
+
-
*(solvector+temp)= *(*(c+temp)+nrows)/ *(*(c+temp)+temp);
+
-
}
+
-
return(solvector);
+
-
}
+
-
 
+
-
void writeoutputvector(double *p, int nrows)
+
-
{
+
-
int temp;
+
-
for(temp=0; temp<nrows; temp++)
+
-
printf("\nTemperature at Node %d = %lf", temp+1, *(p+temp));
+
-
}
+

Revision as of 08:01, 9 June 2006

My wiki