| unkn0wn |
September 16, 2013 06:37 |
Lid Driven Unsteady using ADI method.
Hello all,
We were trying to solve unsteady code for Lid Driven Cavity. Using Alternative Direction Implicit method with TDMA Algorithm. We are getting solution upto few time steps but after that omega is blowing up.
Here the code we are trying.
Any ideas whats going wrong in our code.
Code:
#include<iostream>
#include<cstdlib>
#include<cmath>
#include<vector>
#include<cmath>
#include<fstream>
#include<string>
#define C_ERROR 0.000000001
#define grid 101
using namespace std;
vector< vector<double> > psi;
vector< vector<double> > omega;
vector< vector<double> > u;
vector< vector<double> > v;
double dx=0, dy=0, beta=0, U=0;
double dt=0, h=0;
void Initialise(){
psi.resize(grid, vector<double>(grid,0));
omega.resize(grid, vector<double>(grid,0));
u.resize(grid, vector<double>(grid,0));
v.resize(grid, vector<double>(grid,0));
for(uint i=0;i<u[0].size();i++)
u[0][i]=1;
dx=dy=(double)1/grid;
beta=dx/dy;
U=1;
h=(double)1/grid;
dt=0.001;
}
void Print(const vector< vector<double> > &v){
for(uint i=0;i<v.size();i++){
for(uint j=0;j<v[i].size();j++)
cout<<"["<<i<<"]"<<"["<<j<<"]"<<v[i][j]<<"\t";
cout<<endl;
}
cout<<endl;
}
void PrintSingle(const vector< double > &v){
for(uint i=0;i<v.size();i++){
cout<<"["<<i<<"]"<<v[i]<<"\t";
}
cout<<endl<<endl;
}
void SetBoundaryConditions(){
uint i=0;
for(i=1;i<grid-1;i++){
omega[0][i]=-(2*psi[1][i])/(dx*dx);
omega[grid-1][i]=-(2*psi[grid-2][i])/(dx*dx);
omega[i][0]=-(2*psi[i][1])/(dy*dy);
omega[i][grid-1]=-(2*psi[i][grid-2]+2*dy*U)/(dy*dy);
}
}
double Error(vector< vector<double> > a, vector< vector<double> > b){
double diffSum=0, nSum=0;
for(uint i=0;i<a.size();i++)
for(uint j=0;j<b.size();j++) {
nSum+=fabs(a[i][j]);
diffSum+=fabs(a[i][j]-b[i][j]);
}
return diffSum/nSum;
}
void thomas(vector<double>& a,
vector<double>& b,
vector<double>& c,
vector<double>& d,
vector<double>& f) {
size_t N = d.size();
for (int i=1; i<N; i++) {
b[i]=b[i]-c[i-1]*a[i]/b[i-1];
d[i]=d[i]-d[i-1]*a[i]/b[i-1];
}
f[N-1]=d[N-1]/b[N-1];
for (int i=N-1; i-- > 0; ) {
f[i] = (d[i]-c[i]*f[i+1])/b[i];
}
}
void WriteFile(double a,uint Re){
ofstream File;
File.open("unsteady.dat", ios::out | ios::app);
for(uint i=0;i<grid;i++){
for(uint j=0;j<grid;j++)
File<<a<<" "<<(i+1)*dx<<" "<<(j+1)*dy<<" "<<psi[i][j]<<" "<<omega[i][j]<<" "<<u[i][j]<<" "<<v[i][j]<<endl;
}
}
void Solve(double Re){
uint i=0, j=0, iter=0, mainiter=0;
double error=0;
vector< vector<double> > tempomega;
vector< vector<double> > temppsi;
do{
SetBoundaryConditions();
tempomega=omega;
mainiter++;
do{
//Computing PSI
temppsi=psi;
for(i=1;i<psi.size()-1;i++)
for(j=1;j<psi[i].size()-1;j++)
psi[i][j]=((psi[i+1][j]+psi[i-1][j]+beta*beta*(psi[i][j+1]+psi[i][j-1])+(dx*dx)*omega[i][j]))/(2*(1+beta*beta));
}while(Error(psi,temppsi)>C_ERROR);
//Print(psi);
//Solving u and v
for(i=1;i<grid-1;i++){
for(j=1;j<grid-1;j++){
u[i][j]=-(psi[i][j+1]-psi[i][j-1])/(2*dy);
v[i][j]=(psi[i+1][j]-psi[i-1][j])/(2*dy);
}
}
vector<double> a(grid-2,0);
vector<double> b(grid-2,0);
vector<double> c(grid-2,0);
vector<double> d(grid-2,0);
vector<double> x(grid-2,0);
for(j=1;j<grid-1;j++){
for(i=1;i<grid-1;i++){
a[i-1]=(-u[i][j]/h - 1/(Re*h*h));
b[i-1]=((1/dt)+(2/(Re*h*h)));
c[i-1]=((u[i][j]/h)-(1/(Re*h*h)));
d[i-1]=(-v[i][j]*((omega[i][j+1]-omega[i][j-1])/h) + (omega[i][j+1]+omega[i][j-1]-2*omega[i][j])/(Re*h*h) + (omega[i][j]*2)/dt);
}
omega[0][j]=(-v[0][j]*((omega[0][j+1]-omega[0][j-1])/h) + (omega[0][j+1]+omega[0][j-1]-2*omega[0][j])/(Re*h*h) + (omega[0][j]*2)/dt);;
omega[grid-1][j]=(-v[grid-1][j]*((omega[grid-1][j+1]-omega[grid-1][j-1])/h) + (omega[grid-1][j+1]+omega[grid-1][j-1]-2*omega[grid-1][j])/(Re*h*h) + (omega[grid-1][j]*2)/dt);
d[0]=d[0]-a[0]*(-v[0][j]*((omega[0][j+1]-omega[0][j-1])/h) + (omega[0][j+1]+omega[0][j-1]-2*omega[0][j])/(Re*h*h) + (omega[0][j]*2)/dt);
d[grid-3]=d[grid-3]-c[grid-3]*(-v[grid-1][j]*((omega[grid-1][j+1]-omega[grid-1][j-1])/h) + (omega[grid-1][j+1]+omega[grid-1][j-1]-2*omega[grid-1][j])/(Re*h*h) + (omega[grid-1][j]*2)/dt);
a[0]=0;
c[grid-3]=0;
thomas(a,b,c,d,x);
for(i=1;i<grid-1;i++){
omega[i][j]=x[i-1];
}
}
fill(a.begin(), a.end(), 0);
fill(b.begin(), b.end(), 0);
fill(c.begin(), c.end(), 0);
fill(d.begin(), d.end(), 0);
fill(x.begin(), x.end(), 0);
for(i=1;i<grid-1;i++){
for(j=1;j<grid-1;j++){
a[j-1]=(-v[i][j]/h)-(1/(Re*h*h));
b[j-1]=((1/dt)+(2/(Re*h*h)));
c[j-1]=(v[i][j]/h) -(1/(Re*h*h));
d[j-1]=(-u[i][j]*(omega[i+1][j]-omega[i-1][j]))/h + (omega[i+1][j]+omega[i-1][j]-2*omega[i][j])/(Re*h*h) + (2*omega[i][j])/dt;
}
omega[i][0]=(-u[i][0]*(omega[i+1][0]-omega[i-1][0]))/h + (omega[i+1][0]+omega[i-1][0]-2*omega[i][0])/(Re*h*h) + (2*omega[i][0])/dt;
omega[i][grid-1]=(-u[i][grid-1]*(omega[i+1][grid-1]-omega[i-1][grid-1]))/h + (omega[i+1][grid-1]+omega[i-1][grid-1]-2*omega[i][grid-1])/(Re*h*h) + (2*omega[i][grid-1])/dt;
d[0]=d[0]-a[0]*(-u[i][0]*(omega[i+1][0]-omega[i-1][0]))/h + (omega[i+1][0]+omega[i-1][0]-2*omega[i][0])/(Re*h*h) + (2*omega[i][0])/dt;
d[grid-3]=d[grid-3]-c[grid-3]*(-u[i][grid-1]*(omega[i+1][grid-1]-omega[i-1][grid-1]))/h + (omega[i+1][grid-1]+omega[i-1][grid-1]-2*omega[i][grid-1])/(Re*h*h) + (2*omega[i][grid-1])/dt;
a[0]=0;
c[grid-3]=0;
thomas(a,b,c,d,x);
for(j=1;j<grid-1;j++){
omega[i][j]=x[j-1];
}
}
Print(omega);
WriteFile(mainiter*dt,Re);
cout<<"Time Step : "<<mainiter*dt<<endl;
}while(Error(omega,tempomega)>C_ERROR);
}
int main(){
Initialise();
cout<<"Switching to Re=400\n";
Solve(100);
/*
//Test Case to cCheck Thomas Algorithm
vector<double> a(4,0);
vector<double> b(4,0);
vector<double> c(4,0);
vector<double> d(4,0);
vector<double> x(4,0);
d[0]=8;
d[1]=116;
d[2]=47;
d[3]=40;
a[0]=0;
a[1]=4;
a[2]=2;
a[3]=4;
b[0]=2;
b[1]=5;
b[2]=5;
b[3]=7;
c[0]=3;
c[1]=34;
c[2]=7;
c[3]=0;
thomas(a,b,c,d,x);
cout<<"A: "<<a[0]<<" "<<a[1]<<" "<<a[2]<<" "<<a[3]<<"\n";
cout<<"B: "<<b[0]<<" "<<b[1]<<" "<<b[2]<<" "<<b[3]<<"\n";
cout<<"C: "<<c[0]<<" "<<c[1]<<" "<<c[2]<<" "<<c[3]<<"\n";
cout<<"D: "<<d[0]<<" "<<d[1]<<" "<<d[2]<<" "<<d[3]<<"\n";
cout<<"X: "<<x[0]<<" "<<x[1]<<" "<<x[2]<<" "<<x[3]<<"\n";
*/
}
|
