I started watching the CFD lectures provided by Boston University via iTunes, and thought I could share some bad code with the world. The lectures apparently consist of 12 steps to solving the Navier-Stokes equations, where the first one is to solve a 1D convection equation. The class used Python, but I will be using C++.

subject to some initial condition .The analytical solution to this problem is

We will approximate the left-hand side by a forward difference, and the right hand side by a backward difference. Read the post on finite-differences if the topic is unfamiliar to you.

The finite-differences we will use, together with their truncation errors are:

and

To simplify the final expression, we introduce the notation . Using a different letter for the computed solution indicates that the computed solution is different from the real one, . We are now able to rewrite equation (1) in its finite-difference form:

The initial condition makes sure that is known for all . We can therefore solve for to get,

I won’t go into details about the implementation, because it is rather simple. The initial condition is chosen to be a step-function. I think the reason Prof. Barba has chosen this function is that it is not continuous, and solving using finite-differences will give pretty bad results. The output is saved to a *.dat* file. I import the .dat file into MATLAB to generate plots. As expected, the results are not great. It’s quite instructive to play around with the step sizes, number of meshpoints and the initial condition. I apologize for not naming my variables as they are named in the above derivation.

/* * Solves 1D convection equation, also known as the one-way wave equation. */ #include <iostream> #include <fstream> #include <vector> using namespace std; void saveArray(vector< vector<double> > &u, int nx, int nt, double dx, double dt, int c); int main() { int nx = 20; int nt = 50; double dx = double(2)/(nx-1); double dt = 0.01; int c = 5; vector< vector<double> > u(nx,vector<double> (nt)); vector< vector<double> > un(nx,vector<double> (nt)); // Initial condition: for (int i=0; i <= nx-1; i++) { if (i*dx >= 0.5 && i*dx <= 1) { u[i][0] = 2; } else { u[i][0] = 1; } } // Finite-difference loop: for (int it=1; it<=nt-1; it++) { for (int k=0; k<=nx-1; k++) { un[k][it-1] = u[k][it-1]; } for (int i=1; i<=nx-1; i++) { u[0][it] = un[1][it-1]; u[i][it] = un[i][it-1] - c*dt/dx*(un[i][it-1]-un[i-1][it-1]); } } saveArray(u, nx, nt, dx, dt, c); return 0; } // Save array to file: void saveArray(vector< vector<double> > &u, int nx, int nt, double dx, double dt, int c) { ofstream myfile; myfile.open("1d_convection02.dat"); myfile << "%\t" << "nx = " << nx << "\tnt = " << nt << "\tdt = " << dt << "\tc = " << c << endl; for (int i=0; i<=nx-1; i++) { myfile << i*dx << "\t\t"; for (int j=0; j<=nt-1; j++) { myfile << u[i][j] << "\t\t"; } myfile << endl; } }

## 3 comments

Comments feed for this article

June 8, 2011 at 06:45

1D nonlinear convection equation (Inviscid Burgers) « DaFeda's Blog[...] In the last post, I posted code that solves the 1D linear convection equation. Here, I’ll just post the results of solving a 1D nonlinear convection equation, also known as the inviscid (i.e., zero viscosity) Burgers equation. It’s all part of the 12 step program to solving the Navier-Stokes equations, as taught by Prof. Barba at Boston University. The 1D nonlinear convection equation is of the form, [...]

June 9, 2011 at 06:15

1D Diffusion equation (Heat equation) « DaFeda's Blog[...] is the viscosity. The process is very similar to solving a 1D convection equation, but this time we need to approximate a second derivative. An easy way to get an expression that [...]

July 23, 2013 at 00:28

CFD Python: 12 steps to Navier-Stokes :: Lorena A. Barba Group[…] iTunesU and YouTube since 2010. In fact, at least two students have blogged about their experience: DaFeda in Norway (June 2011) and Sidafa Conde at UMass-Dartmouth (January […]