EML 5725 Computational Fluid Dynamics Spring 2003
Groups Project 1

Write a program that solves the following unidirectional flow problem for the unknown velocity v(x,t):

P.D.E.:

Boundary conditions:

Initial condition:

where the viscosity , (not to be confused with the velocity v), and V₀, V₁, and V₂ are all given constants.

You first need to create an input file ftcs.in for your program with the data for the problem to be solved. In the file, initially insert the following values of the variables followed by their descriptions:

0.3     coefficient of viscosity nu
2.0     length l
3.0     magnitude of the parabolic part of the initial velocity V_0
2.0     velocity at the left end V_1
1.0     velocity at the right end V_2
4       maximum value of spatial index j (i.e. number of mesh cells) j_max
0.5     nondimensional mesh viscosity number Nm = nu dt / dx^2
0.5     time interval between doing output t_out

As far as is concerned, your program will use this value to compute the time step. The scheme is stable if .

As far as is concerned, your program should output the velocity distribution v as a function of x at times and . Then use the software of your choice to plot these 5 velocity curves in a single v,x-plot.

Now write the program. Give the program file the name ftcs.c, ftcs.cpp, or ftcs.f, depending on what language you are using. The structure of the program should be as follows, in order of appearance:

1.

If you are using C or C++, first the usual C includes. Also a define of a constant j_dim (alternatively, make it a global constant), with a value equal to the maximum number of mesh points that you would ever want to use, say 1,001. And of course the void main(void) starting lines.

In Fortran, the first lines should be a program statement, as well as an implicit none statement, followed by a declaration of an integer j_dim followed by a parameter statement making j_dim a parameter with a value equal to the maximum number of mesh cells that you would ever want to use, say 1,000.

2.

Next, declare all variables, including the physical variables nu, l, V_0, V_1, and V_2, and the computational variables, such as the mesh indices and their maximum values j, n, j_max and n_max, the mesh sizes dx and dt, etcetera.

You also want to declare a variable Nm to contain the value of the nondimensional group .

You also need to declare storage for the velocities to be computed. In C or C++, declare the storage as float v[j_dim][2], while in Fortran, declare it as real v(0:j_dim,0:1). Note from this that only two time planes are stored. In particular, any velocity vⁿ_j must be stored as v[j][n%2] in C/C++ and as v(j,mod(n,2)) in Fortran during the computations. When computing new times, the new velocities overwrite the older values.

You may also want to declare an array x of dimension j_dim to hold the mesh x-values and a variable t to hold time.

You also need to declare a variable n_out which is the number of time steps to do before doing output and a corresponding time interval t_out

Each declaration must be preceded by a comment line clearly describing what the variable respresents to anyone who might read the program.

For the remaining parts of the program, precede every small thing you do by a comment line explaining what you are doing.

3.

After the declarations, you will want to open an output file for the results to be written to. You also want to open the input file ftcs.in and read the values of nu, l, V_0, V_1, V_2, j_max, Nm, and t_out from it.

Print each value out to the output file, with its description, as it is read in from the input file. Check each value read in for acceptability. Is nu greater then zero? Is j_max greater than 1? Is it less than j_dim, so that we have enough storage to store all values computed? Etcetera. If an unacceptable value is found, print out a message what precisely is wrong and terminate program execution.

You probably also want to open some plot file for the plot of the velocities v at the requested 5 times against x.

4.

Compute various derived program quantities, such as dx, dt, n_out, etcetera, using what you read in from the file. The maximum value of n to compute, n_max, will be 4 n_out since there is no point in keeping computing if you are not going to plot the results.

5.

Using the given initial condition and the two boundary conditions, store the j_max+1 initial velocity values at time t=0 in array locations v[j][0] or v(j,0).

Write these initial velocities also to your plot file.

To your output file, write a single line saying that the current time t=0, followed by the minimum and maximum initial temperatures.

6.

The final part of the program will be a big for or do loop with the main purpose to compute the line of vⁿ⁺¹_j velocities for all j-values from already known velocities vⁿ_j. The loop will start at n=0 and terminate at . In short, the first time through this loop we compute the v¹_j-values from the already stored v⁰_j-values; the second time through the loop the v²_j-values from the now stored v¹_j-values; etcetera, until the last time where we compute the from the .

The subparts inside this n-loop will be

(a)

``Compute'' the boundary value vⁿ⁺¹₀.

(b)

Inside a loop that varies the value of j from 1 to j_max-1, compute the vⁿ⁺¹_j from the vⁿ_j using the following formula:

Program in the solution of this equation without precomputing any quantities and only using the variables shown in the formula. (No N.)

(c)

``Compute'' the boundary value .

(d)

Search through the just computed vⁿ⁺¹_j-values to find the smallest and largest velocities in the flow at that time and print out a single line with the time tⁿ⁺¹ as well as these smallest and largest velocities at that time.

(e)

Check whether n+1 is a multiple of n_out and if so, write all the just computed velocity values vⁿ⁺¹_j to your plot file. Also print a message to the output file telling you which time you are plotting; it will not be exactly the time you want because of the finite time step.

Run your program and plot the results using your favorite plotting program. Plot the computed velocities as individual points, not as a line. Best is to use a different plot symbol for each of the five times.

Increase the value of j_max by a factor 2 and see what happens to the computed solution. Increase it by another factor 2 and see.

Now change back to the original j_max and make the value of Nm equal to 0.25. What happens to the solution?

Now make Nm equal to 1. What happens to the solution? Try increasing j_max to get more accuracy. Does it help?

Comment on all results.