CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Main CFD Forum

Lid Driven Unsteady using ADI method.

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Display Modes
Old   September 16, 2013, 06:37
Smile Lid Driven Unsteady using ADI method.
  #1
New Member
 
Rohan Sagar
Join Date: Sep 2013
Posts: 1
Rep Power: 0
unkn0wn is on a distinguished road
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";
*/


}
unkn0wn is offline   Reply With Quote

Reply

Tags
lid driven cavity, tmsa

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Lid Driven Cavity velocity profiles new_at_this Main CFD Forum 0 December 22, 2011 17:04
is there any parallel code for the famous Lid Driven Cavity flow? gholamghar Main CFD Forum 0 August 1, 2010 01:55
2D Lid Driven Cavity Flow simulation using MATLAB josephlm Main CFD Forum 3 June 24, 2010 01:58
Two sided Lid Driven Cavity Flow Perumal Main CFD Forum 1 January 22, 2007 14:27
pressure driven flow by pressure correction method justentered Main CFD Forum 0 December 30, 2003 00:52


All times are GMT -4. The time now is 16:53.