The goal of this project is to make a numerical model of the Couette flow - An incompressible flow between two parallell plates where the upper plate is travelling with a constant velocity U relative to the lower plate. The code will be written in C++ and will use ROOT for visualization

For an incompressible flow in the x-direction the governing equations (Navier-Stokes eqs.) reduces to$$\begin{array}{cccc}\rho \frac{\partial u}{\partial t}& =& \mu \frac{{\partial}^{2}u}{\partial t{}^{2}}& (1)\end{array}$$

Where u = u(y) is the velocity of the flow, ρ is density, μ is viscosity and t is time. U is the velocity of the upper plate, y is the distance from the lower plate, and D is the distance between the plates. We've also used the no-slip boundary conditions so that the flow velocity at the lower plate is zero and at the upper plate U.

$$\begin{array}{cccc}\rho \frac{\partial u}{\partial t}& =& \frac{1}{Re}\frac{{\partial}^{2}u}{\partial t{}^{2}}& (2)\end{array}$$

To solve this equation we follow the so-called “Crank-Nicolson” approach, a finite-difference technique which leads to the following equation:

$$\begin{array}{cc}{K}_{j}=\left[1-\frac{\u2206t}{(\u2206y{)}^{2}Re}\right]{u}_{j}^{n}+\frac{\u2206t}{2(\u2206y{)}^{2}Re}\left({u}_{j+1}^{n}+{u}_{j-1}^{n}\right)& (6)\end{array}$$

Re is the Reynolds number (=ρDu/μ), j is grid point index, n is iteration step so that the velocities on the left side of eq. (3) are unknown values (next, n+1, iteration step) as function of the known K values. If we chose N+1 grid points (in the y-direction) the separation between the points will be ∆y = D/N.

$\left[\begin{array}{ccccc}B& A& 0& 0& 0\\ A& B& A& 0& 0\\ 0& A& B& \mathrm{.}& 0\\ 0& 0& \mathrm{.}& \mathrm{.}& A\\ 0& 0& 0& A& B\end{array}\right]\left[\begin{array}{c}{u}_{2}^{n+1}\\ {u}_{3}^{n+1}\\ \mathrm{.}\\ \mathrm{.}\\ {u}_{N}^{n+1}\end{array}\right]=\left[\begin{array}{c}{K}_{2}\\ {K}_{3}\\ \mathrm{.}\\ \mathrm{.}\\ {K}_{N}-AU\end{array}\right]$

The program asks the user for values of Reynolds number, number of grid points N and time steps ∆t. Then it calculates (using Thomas' Algorithm) flow velocities in grid points at next time step. For each iteration step to the fourth power (1, 16, 81, 256..) a graph showing the velocity profile is drawn until the final, steady state is achieved. This happens when the velocities at all grid points doesn't deviate from the linear velocity distribution more than some accuracy epsilon (set to 0.001). At the end, the steady state solution is drawn.

The parameter $E=\frac{\u2206t}{Re(\u2206y{)}^{2}}$ is important for the numerical behavior of the code. The code will be stable as long as E is less than about 4000. although when E is large, we will get some nonphysical transient results but still obtain the correct steady state solution. Even for E equal to 4 or 5 we see some strange velocities near the upper plate after the first step.

For low values of E (<5) we have timewise accurate solutions, meaning that at same total relative time t, we have the same velocity distribution even though we use different sizes of the time steps. How many iteration steps we need is also dependent on E, and it turns out that E approximately equal to 10 is when we need the fewest amount of iteration steps to reach the steady state and therefore use the least computer power.

Underneath one can see screenshots of the velocity profile with E = 1 (figure 1), E = 4 (figure 2), E = 10 (figure 3) and the extreme case when E = 4000 (figure 4)

It has been good practice for me to get a feeling for simple CFD through working with the Couette Flow, and also to get familiar with some of the possibilities that comes with ROOT. In addition to what I did I also wanted to do a vector plot. Unfortunately I wasn't able to do this in ROOT.