#include <iostream>
#include <fstream>
#include <iomanip>
#include <math.h>
#include <vector>
#include <float.h>
using namespace std;
void T_top();
void T_middle(int a);//a refers to the number of row
void T_bottom();
int row = 80, column = 80;
double density = 900.0, Cp = 500.0, Fo = 10; //For transient process
double L=0.6,W=0.6,lambda=4.5,Q=25000.0,q_right=1400.0,alpha=15.0,T_inf=0.0,Tw=20.0,omega=1.6;
double h=L/row;
double R=h/(2.0*lambda)+1.0/alpha;
const double errorLimit = 0.0000000001;
const double realTime = 100000; //For transient process
double maxError , error;
vector<double> A(column), B(column), C(column), D(column), X(column);
vector<vector<double>> T_old(row,vector<double>(column,25.0));
vector<vector<double>> T_new(row,vector<double>(column,25.0));
int main(){
double delta_t = Fo*density*Cp*h*h/lambda;
double T_center;
ofstream f;
f.open("result.txt");
if(!f.is_open()){
cout << "Failed opening result.txt" << endl;
return 1;
}
f << "time T_center difference" << endl;
for(int i =1; i*delta_t <= realTime; i++){
do{
maxError = DBL_MIN;
T_top();
for(int i=1; i<row-1; i++){
T_middle(i);
}
T_bottom();
}while(maxError >= errorLimit);
T_old.assign(T_new.begin(),T_new.end()); //T_old = T_new
T_center = (T_new[row/2-1][row/2-1]+T_new[row/2-1][row/2]+T_new[row/2][row/2-1]+T_new[row/2][row/2])/4;
f << fixed << left << setw(10) << setprecision(3) << i*delta_t << setw(10) << T_center << (T_center-25)/(461.62-25) << endl;
}
cout<<"write to result.txt!"<<endl;
f.close();
return 0;
}
void TDMA(){
B[0]=B[0]/A[0];
D[0]=D[0]/A[0];
for(int i=1;i<column;i++){
B[i]=B[i]/(A[i]-C[i]*B[i-1]);
D[i]=(D[i]+C[i]*D[i-1])/(A[i]-C[i]*B[i-1]);
}
X[column-1]=D[column-1];
for(int i=column-2;i>=0;i--){
X[i]=D[i]+B[i]*X[i+1];
}
}
void T_top(){
A[0]=2+h/lambda/R+1/Fo;
A[column-1]=2+h/lambda/R+1/Fo;
B[0]=1;
B[column-1]=0;
C[0]=0;
C[column-1]=1;
for(int i=1;i<column-1;i++){
A[i]=3+h/lambda/R+1/Fo;
B[i]=1;
C[i]=1;
}
D[0]=T_new[1][0]+h/lambda/R*T_inf+Q/lambda*h*h+1/Fo*T_old[0][0];
for(int i=1;i<column-1;i++){
D[i]=T_new[1][i]+h/lambda/R*T_inf+Q/lambda*h*h+1/Fo*T_old[0][i];
}
D[column-1]=T_new[1][column-1]+h/lambda/R*T_inf+Q/lambda*h*h+q_right*h/lambda+1/Fo*T_old[0][column-1];
TDMA();
for(int i=0;i<column;i++){ //SOR
T_new[0][i]=omega*X[i]+(1-omega)*T_new[0][i];
error=fabs((T_new[0][i]-X[i])/T_new[0][i]);
if(maxError<error){
maxError=error;
}
}
}
void T_middle(int a){ //a refers to the number of row
A[0]=3+1/Fo;
A[column-1]=3+1/Fo;
B[0]=1;
B[column-1]=0;
C[0]=0;
C[column-1]=1;
for(int i=1;i<column-1;i++){
A[i]=4+1/Fo;
B[i]=1;
C[i]=1;
}
D[0]=T_new[a-1][0]+T_new[a+1][0]+Q/lambda*h*h+1/Fo*T_old[a][0];
for(int i=1;i<column-1;i++){
D[i]=T_new[a-1][i]+T_new[a+1][i]+Q/lambda*h*h+1/Fo*T_old[a][i];
}
D[column-1]=T_new[a-1][column-1]+T_new[a+1][column-1]+Q/lambda*h*h+q_right*h/lambda+1/Fo*T_old[a][column-1];
TDMA();
for(int i=0;i<column;i++){
T_new[a][i]=omega*X[i]+(1-omega)*T_new[a][i];
error=fabs((T_new[a][i]-X[i])/T_new[a][i]);
if(maxError<error){
maxError=error;
}
}
}
void T_bottom(){
A[0]=4+1/Fo;
A[column-1]=4+1/Fo;
B[0]=1;
B[column-1]=0;
C[0]=0;
C[column-1]=1;
for(int i=1;i<column-1;i++){
A[i]=5+1/Fo;
B[i]=1;
C[i]=1;
}
D[0]=T_new[row-2][0]+2*Tw+Q/lambda*h*h+1/Fo*T_old[row-1][0];
for(int i=1;i<column-1;i++){
D[i]=T_new[row-2][i]+2*Tw+Q/lambda*h*h+1/Fo*T_old[row-1][i];
}
D[column-1]=T_new[row-2][column-1]+2*Tw+Q/lambda*h*h+q_right*h/lambda+1/Fo*T_old[row-1][column-1];
TDMA();
for(int i=0;i<column;i++){
T_new[row-1][i]=omega*X[i]+(1-omega)*T_new[row-1][i];
error=fabs((T_new[row-1][i]-X[i])/T_new[row-1][i]);
if(maxError<error){
maxError=error;
}
}
}