On this page you’ll find some C++ source code I have written for the university. I am almost an mechanical engineer but I like programming, so don’t wonder if the code is not perfect. Its working without an error and it does its job :)
This code I had to write to calculate the temperature through a wall consisting of a two layer coating. Each coating has different values for the thermal conductivity and thickness. The inner wall was heated to a high temperature at t>0. For t=0 the wall had a constant and equal temperature everywhere. This was a unsteady heat transfer problem. One problem was the thickness of the wall and the required accuracy for the numerical solution. The second problem was that the complete heating of the wall to a steady state took rather long because of the thickness. This means the numerical calculation needs many POSITION and TIME nodes to be accurate and stable. To reduce the memory consumption for this code I designed the program in a way that only one matrix is needed to store the values at each POSITON node (idea of my supervising tutor). This made the code rather long by using many IF-statements.
Here is the code:
# include <iostream>
# include <math.h>
# include <iomanip>
# include <fstream>
# include <string>
# include <sstream>
# include <stdio.h>
# include <stdlib.h>
# include <time.h>
# include <algorithm>
using namespace std;
using std::cin;
using std::ofstream;
using std::endl;
/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
segmentation fault n. [Unix] 1. [techspeak] An error in which a running
program attempts to access memory not allocated to it and core dumps
with a segmentation violation error. This is often caused by improper
usage of pointers in the source code, dereferencing a null pointer, or
(in C) inadvertently using a non-pointer variable as a pointer.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/
int main()
{
clock_t start_nc = 0.0, finish_nc = 0.0, start_c = 0.0, finish_c = 0.0;
int i, m, nodes, e1, e2, e3, e4, e5, counter_beg, counter_end, delta_c;
double temp_nc_old[1201], temp_c_old[1201], temp_nc_new[1201], temp_c_new[1201], qstr_beg[1201], qstr_end[1201], rho[1201], cp[1201], lambda[1201];
double length, dx, t, t_h, dt, x, a1, a2, t1, time_nc, time_c, rho_aa, cp_aa;
// Deleting old files...
if (remove("Output/plot_data_non_conservative.dat")==0)
{
e1=0;
}
if (remove("Output/plot_data_conservative.dat")==0)
{
e2=0;
}
if (remove("Output/qstr_beg.dat")==0)
{
e3=0;
}
if (remove("Output/qstr_end.dat")==0)
{
e4=0;
}
if (remove("Output/analytical_tmax.dat")==0)
{
e5=0;
}
if ((e1==0) and (e2==0) and (e3==0) and (e4==0) and (e5==0))
{
std::cout << "\nAll files successfully deleted!" << endl;
std::cout << endl;
}
else
{
std::cout << "\nOne or more files could not be deleted!" << endl;
std::cout << endl;
}
// input & output
std::cout << "Number of position-nodes (multiple of 3 e.g. 301, 601, 901... but +1): ";
std::cin >> nodes;
std::cout << "Duration (in hours): ";
std::cin >> t_h;
// calculating the index of the node before wall 2 and converting it from ***double*** to ***int***
double node_before_wall2_double = floor ((nodes-1)/3.0);
int node_before_wall2_int;
node_before_wall2_int = static_cast <int> (node_before_wall2_double);
// calculating the usual stuff...
t = t_h*3600;
length = 0.75;
dx = length / (double) (nodes-1);
a1 = (5.2) / (4200*1300);
a2 = (0.8) / (2200*880);
// calculating tnodes_d and converting it into ***int***
dt = 0.4*((pow(dx, 2))/a1);
double tnodes_d = ceil(t / dt);
int tnodes;
tnodes = static_cast <int> (tnodes_d);
dt = (t/tnodes);
// setting average values
rho_aa=3200.0;
cp_aa=1090.0;
double lambda1 = 5.2;
double lambda2 = 0.8;
// calculating the analytical temperature and the heat flow density at the node before and after wall 2
double qstr_ana = ((1573.15 - 293.15)/((0.25/5.2)+(0.50/0.8)));
double T_12 = 293.15 + ((0.50/0.80)*qstr_ana);
double T_05 = floor (887.435);
double temp_c_new_fl;
double temp_nc_new_fl;
int c=0;
int d=0;
if (tnodes <= 501) // 50 till 500
{
delta_c = 1;
}
if ((tnodes > 501) and (tnodes <= 5001)) // 50 till 500
{
delta_c = 10;
}
if ((tnodes > 5001) and (tnodes <= 50001)) // 50 till 500
{
delta_c = 100;
}
if ((tnodes > 50001) and (tnodes <= 500001)) // 50 till 500
{
delta_c = 1000;
}
if ((tnodes > 500001) and (tnodes <= 5000001)) // 50 till 500
{
delta_c = 10000;
}
if ((tnodes > 5000001) and (tnodes <= 50000001)) // 50 till 500
{
delta_c = 100000;
}
// [i] ist for the position index, ***m*** for the time position in the loop
// file output
/*
ios::app ...All output operations are performed at the end of the file, appending the content to the current content of the file.
This flag can only be used in streams open for output-only operations.
*/
ofstream fobj_t_c("Output/plot_data_conservative.dat", ios::app);
ofstream fobj_t_nc("Output/plot_data_non_conservative.dat", ios::app);
ofstream fobj_qb("Output/qstr_beg.dat", ios::app);
ofstream fobj_qe("Output/qstr_end.dat", ios::app);
ofstream fobj_ana_tmax("Output/analytical_tmax.dat", ios::app);
ofstream fobj_time("Output/program_elapsed_time.dat", ios::app);
/*
*
* non-conservative and conservative block for *** m = 0 ***
*
*/
t1=0.0;
x=0.0;
for (i=0; i<=(nodes-1); i++)
{
if (i<((nodes-1)/3))
{
rho[i]=4200.0;
cp[i]=1300.0;
lambda[i]=5.2;
}
if (i>((nodes-1)/3))
{
rho[i]=2200.0;
cp[i]=880.0;
lambda[i]=0.8;
}
temp_nc_old[i]=293.15;
temp_c_old[i]=293.15;
// writing temp. values to file
fobj_t_nc << fixed << setprecision(4) << x << " " << (t1/3600) << " " << fixed << setprecision(2) << temp_nc_old[i] << endl;
fobj_t_c << fixed << setprecision(4) << x << " " << (t1/3600) << " " << fixed << setprecision(2) << temp_c_old[i] << endl;
if (i==(nodes-1))
{
fobj_t_nc << endl;
fobj_t_nc << endl;
fobj_t_c << endl;
fobj_t_c << endl;
}
x = x + dx;
}
// first values for qstr_beg and qstr_end for t1=0
x=0.0;
t1=0.0;
qstr_beg[0]=0.00;
fobj_qb << fixed << setprecision(4) << x << " " << (t1/3600) << " " << fixed << setprecision(2) << qstr_beg[i] << endl;
x=0.75;
qstr_end[nodes-1]=0.00;
fobj_qe << fixed << setprecision(4) << x << " " << (t1/3600) << " " << fixed << setprecision(2) << qstr_beg[i] << endl;
/*
*
* non-conservative block for *** m > 0 ***
*
*/
start_nc = clock(); // setting first value
t1=dt;
counter_beg=2;
counter_end=2;
for (m=1; (m<=tnodes); m++)
{
x=0.0;
for (i=0; i<=(nodes-1); i++)
{
if (m==1)
{
if (i==0)
{
temp_nc_old[i]=1573.15;
temp_c_old[i]=1573.15;
}
else
{
temp_nc_old[i]=293.15;
temp_c_old[i]=293.15;
}
}
else
{
if (m==2)
{
if ((i==0) or (i==(nodes-1)))
{
temp_nc_new[i]=temp_nc_old[i];
if (i==(nodes-1)) // calculate only AFTER all the other values in the walls have been calculated
{
temp_nc_new[i]=temp_nc_old[i];
temp_nc_new[node_before_wall2_int]=(lambda1*temp_nc_new[node_before_wall2_int-1] + lambda2*temp_nc_new[node_before_wall2_int+1])/(lambda1 + lambda2);
}
}
else
{
if (i<((nodes-1)/3))
{
temp_nc_new[i]=((a1*dt)/(pow(dx, 2.0)))*(temp_nc_old[i+1]-2.0*temp_nc_old[i]+temp_nc_old[i-1]) + temp_nc_old[i];
}
if (i>((nodes-1)/3))
{
temp_nc_new[i]=((a2*dt)/(pow(dx, 2.0)))*(temp_nc_old[i+1]-2.0*temp_nc_old[i]+temp_nc_old[i-1]) + temp_nc_old[i];
}
}
}
else
{
if ((i==0) or (i==(nodes-1)))
{
temp_nc_new[i]=temp_nc_old[i];
if (i==(nodes-1)) // calculate only AFTER all the other values in the walls have been calculated
{
temp_nc_new[i]=temp_nc_old[i];
temp_nc_new[node_before_wall2_int]=(lambda1*temp_nc_new[node_before_wall2_int-1] + lambda2*temp_nc_new[node_before_wall2_int+1])/(lambda1 + lambda2);
}
}
else
{
if (i<((nodes-1)/3))
{
temp_nc_new[i]=((a1*dt)/(pow(dx, 2.0)))*(temp_nc_new[i+1]-2.0*temp_nc_new[i]+temp_nc_new[i-1]) + temp_nc_new[i];
}
if (i>((nodes-1)/3))
{
temp_nc_new[i]=((a2*dt)/(pow(dx, 2.0)))*(temp_nc_new[i+1]-2.0*temp_nc_new[i]+temp_nc_new[i-1]) + temp_nc_new[i];
}
}
}
// heat flow density
if (m==counter_beg)
{
if (i==0)
{
qstr_beg[0] = - 5.2*((temp_nc_new[1] - temp_nc_new[0])/(dx));
fobj_qb << fixed << setprecision(4) << x << " " << (t1/3600) << " " << fixed << setprecision(2) << qstr_beg[i] << endl;
counter_beg = counter_beg + delta_c;
}
}
if (m==counter_end)
{
if (i==(nodes-1))
{
qstr_end[nodes-1] = - 0.8*((temp_nc_new[nodes-1] - temp_nc_new[nodes-2])/(dx));
fobj_qe << fixed << setprecision(4) << x << " " << (t1/3600) << " " << fixed << setprecision(2) << qstr_end[i] << endl;
counter_end = counter_end + delta_c;
}
}
if (m==(tnodes)) // last value for dot_q
{
if (i==0)
{
qstr_beg[0] = - 5.2*((temp_nc_new[1] - temp_nc_new[0])/(dx));
fobj_qb << fixed << setprecision(4) << x << " " << (t1/3600) << " " << fixed << setprecision(2) << qstr_beg[i] << endl;
}
if (i==(nodes-1))
{
qstr_end[nodes-1] = - 0.8*((temp_nc_new[nodes-1] - temp_nc_new[nodes-2])/(dx));
fobj_qe << fixed << setprecision(4) << x << " " << (t1/3600) << " " << fixed << setprecision(2) << qstr_end[i] << endl;
}
}
}
x = x + dx;
}
x=0.0;
for (i=0; i<=(nodes-1); i++)
{
// writing temp. values to file
if (m==tnodes) // for really long times, all curves would the have same values
{
fobj_t_nc << fixed << setprecision(4) << x << " " << (t1/3600) << " " << fixed << setprecision(2) << temp_nc_new[i] << endl;
if (((i==0) or (i==node_before_wall2_int) or (i==(nodes-1))) and (m==tnodes))
{
if (i==(node_before_wall2_int))
{
fobj_ana_tmax << fixed << setprecision(4) << x << " " << (t1/3600) << " " << fixed << setprecision(2) << T_12 << endl;
}
if ((i==0) or (i==(nodes-1)))
{
fobj_ana_tmax << fixed << setprecision(4) << x << " " << (t1/3600) << " " << fixed << setprecision(2) << temp_nc_new[i] << endl;
}
}
if (i==(nodes-1)) // print additional lines into the file, for gnuplot recognition of the new index
{
fobj_t_nc << endl;
fobj_t_nc << endl;
}
}
// test if T_ana = T_num
if (i==(2*(nodes-1)/3)) // for every ***m***
{
temp_nc_new_fl = floor(temp_nc_new[i]);
if ((temp_nc_new_fl>=T_05) and (d==0))
{
std::cout << endl;
std::cout << "Non-Conservative: Through heating after: " << fixed << setprecision(3) << (t1/3600) << " h\n" << endl;
d=d+1;
}
if ((temp_nc_new_fl<T_05) and (m==(tnodes))) // print only once!!!
{
std::cout << endl;
std::cout << "Non-Conservative: Through heating could not be accomblished yet!\n" << endl;
}
}
x = x + dx;
}
t1 = t1 + dt;
}
finish_nc = clock(); // setting second value
time_nc = (double)(finish_nc - start_nc) / CLOCKS_PER_SEC;
fobj_time << "elapsed time (non conservative): " << fixed << setprecision(4) << time_nc << " sec" << " " << "real time: " << fixed << setprecision(0) << t_h << " h" << " " << "position-nodes: " << nodes << " " << "time-nodes: " << tnodes << " " << "duration: " << t << " sec" << endl;
fobj_time << endl;
fobj_time << endl;
/*
*
* conservative block for *** m > 0 ***
*
*/
start_c = clock(); // setting third value
t1=dt;
for (m=1; (m<=tnodes); m++)
{
x=0.0;
for (i=0; i<=(nodes-1); i++)
{
if (m==1)
{
if ((i==0) or (i==(nodes-1)))
{
temp_c_new[i]=temp_c_old[i]; // set the old temp. value at the wall
}
else
{
if (i<((nodes-1)/3)) // means you're at x<0.25 which is also correct for the last node!
{
if (i==(node_before_wall2_int-1))
{
temp_c_new[i]=((a1*dt)/(pow(dx, 2.0)))*(temp_c_old[node_before_wall2_int] - 2.0*temp_c_old[node_before_wall2_int-1] + temp_c_old[node_before_wall2_int-2]) + temp_c_old[node_before_wall2_int-1];
}
else
{
temp_c_new[i]=((dt)/(rho[i]*cp[i]))*((lambda[i+1]*(temp_c_old[i+1]-temp_c_old[i])-lambda[i-1]*(temp_c_old[i]-temp_c_old[i-1]))/((pow(dx, 2.0)))) + temp_c_old[i];
}
}
if (i==((nodes-1)/3))
{
// --- WHAT ABOUT rho & cp --- ???
temp_c_new[i]=((dt)/(rho_aa*cp_aa))*((lambda[i+1]*(temp_c_old[i+1]-temp_c_old[i])-lambda[i-1]*(temp_c_old[i]-temp_c_old[i-1]))/((pow(dx, 2.0)))) + temp_c_old[i];
}
if (i>((nodes-1)/3))
{
if (i==(node_before_wall2_int+1))
{
temp_c_new[i]=((a2*dt)/(pow(dx, 2.0)))*(temp_c_old[node_before_wall2_int+2] - 2.0*temp_c_old[node_before_wall2_int+1] + temp_c_old[node_before_wall2_int]) + temp_c_old[node_before_wall2_int+1];
}
else
{
temp_c_new[i]=((dt)/(rho[i]*cp[i]))*((lambda[i+1]*(temp_c_old[i+1]-temp_c_old[i])-lambda[i-1]*(temp_c_old[i]-temp_c_old[i-1]))/((pow(dx, 2.0)))) + temp_c_old[i];
}
}
}
}
else
{
if ((i==0) or (i==(nodes-1)))
{
temp_c_new[i]=temp_c_old[i];
}
else
{
if (i<((nodes-1)/3))
{
temp_c_new[i]=((dt)/(rho[i]*cp[i]))*((lambda[i+1]*(temp_c_new[i+1]-temp_c_new[i])-lambda[i-1]*(temp_c_new[i]-temp_c_new[i-1]))/((pow(dx, 2.0)))) + temp_c_new[i];
}
if (i==(node_before_wall2_int-1))
{
temp_c_new[i]=((a1*dt)/(pow(dx, 2.0)))*(temp_c_new[node_before_wall2_int] - 2.0*temp_c_new[node_before_wall2_int-1] + temp_c_new[node_before_wall2_int-2]) + temp_c_new[node_before_wall2_int-1];
}
if (i==((nodes-1)/3))
{
// --- WHAT ABOUT rho & cp --- ???
temp_c_new[i]=((dt)/(rho_aa*cp_aa))*((lambda[i+1]*(temp_c_new[i+1]-temp_c_new[i])-lambda[i-1]*(temp_c_new[i]-temp_c_new[i-1]))/((pow(dx, 2.0)))) + temp_c_new[i];
}
if (i>((nodes-1)/3))
{
if (i==(node_before_wall2_int+1))
{
temp_c_new[i]=((a2*dt)/(pow(dx, 2.0)))*(temp_c_new[node_before_wall2_int+2] - 2.0*temp_c_new[node_before_wall2_int+1] + temp_c_new[node_before_wall2_int]) + temp_c_new[node_before_wall2_int+1];
}
else
{
temp_c_new[i]=((dt)/(rho[i]*cp[i]))*((lambda[i+1]*(temp_c_new[i+1]-temp_c_new[i])-lambda[i-1]*(temp_c_new[i]-temp_c_new[i-1]))/((pow(dx, 2.0)))) + temp_c_new[i];
}
}
}
}
// test if T_ana = T_num
if (i==(2*(nodes-1)/3)) // for every ***m***
{
temp_c_new_fl = floor(temp_c_new[i]);
if ((temp_c_new_fl>=T_05) and (c==0))
{
std::cout << "Conservative: Through heating after: " << fixed << setprecision(3) << (t1/3600) << " h\n" << endl;
c=c+1;
}
if ((temp_c_new_fl<T_05) and (m==(tnodes))) // print only once!!!
{
std::cout << "Conservative: Through heating could not be accomblished yet!\n" << endl;
}
}
// writing temp. values to file
if (m==tnodes)
{
fobj_t_c << fixed << setprecision(4) << x << " " << (t1/3600) << " " << fixed << setprecision(2) << temp_c_new[i] << endl;
if (i==(nodes-1))
{
fobj_t_c << endl;
fobj_t_c << endl;
}
}
x = x + dx;
}
t1 = t1 + dt;
}
finish_c = clock(); // setting fourth value
time_c = (double)(finish_c - start_c) / CLOCKS_PER_SEC;
fobj_time << "elapsed time (conservative): " << fixed << setprecision(4) << time_c << " sec" << " " << "real time: " << fixed << setprecision(0) << t_h << " h" << " " << "position-nodes: " << nodes << " " << "time-nodes: " << tnodes << " " << "duration: " << t << " sec" << endl;
fobj_time << endl;
fobj_time << endl;
}C++ Code written on openSUSE Linux.
$ whoami _ 