ME6043: Thermal Engineering Fundamentals, Fall 2022
CFD project (Due: Dec 5, 2022)
In this project, you will numerically solve two problems by writing a computer code. Both
problems involve simple partial differential equations (PDEs) with specified boundary/initial
conditions. Every student must work on this project individually without any collaborations.
To execute the project, you can use any programming language you want, e.g. C/C++, Python,
Fortran, etc. You can also use Matlab. Regardless of the programming model, you have to
write the code yourself and cannot use inbuilt libraries/packages for PDE solvers. However,
in both problems you will have to solve a linear system of the form Ax = b, which requires a
direct or indirect matrix inversion. For this inversion, you can use standard routines/libraries
available for your choice of programming language. While the problems are overall quite
straightforward, code writing and debugging can sometimes take up lot of time. Hence, start
early and make steady progress. It might be essentially impossible to finish the project if you
start a week or few days before the deadline.
Problem 1 – Unsteady Couette Flow
In this problem, you will solve the unsteady Couette flow using the methods we learned in class.
For Couette flow, under the parallel and fully-developed flow assumption, the Navier-Stokes
equations simplify to the following PDE:
∂u
∂t
= ν
∂2u
∂y2
, (1)
where u(y, t) is the streamwise velocity, satisfying the boundary conditions
u(0, t) = 0 , u(h, t) = U . (2)
Since the problem is unsteady, we also need an initial condition, which we will take as:
u(y, 0) = U
y
h
+ U sin
(piy
h
)
(3)
Methodology
The preliminary step in solving a PDE numerically typically involves non-dimensionalizing it.
Restate the problem using the dimensionless variables u∗ = u/U , y∗ = y/h and t∗ = νt/h2.
Henceforth, we will drop the superscript ∗ for convenience and all variables are understood to
be non-dimensional (unless stated otherwise).
The next step is to discretize the domain to generate a grid (or mesh). For the current
problem, we will simply use a uniformly spaced grid in y-direction. We will denote the number
of grid points as N . Note, the grid spacing ∆y = 1/(N − 1), and the grid point location
yj = (j − 1)∆y, for j = 1, 2, ..N .
1
Next we will discretize the equation using the nite-dierence approximation. We will use
2nd order central dierence for the spatial derivative and 1st order backwarddierence for the
temporal derivative:
un+1j unj
t
=
un+1j +1 2u
n+1
j + u
n+1
j 1
( y)2 : (4)
Here, n is the the time step number, such that time tn = n t, for n = 0 ; 1; 2::: This is an
unconditionally stable implicit scheme where un+1j can be obtained by solving a system of
linear equations:
ru n+1j 1 (2r + 1) un+1j + ru n+1j +1 = unj ; with r = t=( y)2 : (5)
The system of equations can be solved by writing them in the vector formAx = b, where
x is a vector containing un+1j , b is a vector based onunj (and also the boundary conditions),
and A is a tridiagonal matrix containing the coecients. Note, the solution is only needed at
interior points (2 j N 1), so one must be mindful of the boundary conditions, which will
aect some entries in the vector b. In principle, x can be obtained by simply inverting A, i.e.,
x = A 1b. Since the A matrix is xed, it can be inverted once at the beginning and reused
for time marching. However, given A is tridiagonal, one can also use the Thomas algorithm,
which is more computationally ecient and stable compared to matrix inversion. For the
current assignment, you are free to use either of the two methods. You can utilize publicly
available subroutines/functions/libraries for implementing matrix inversion.
Exact solution and error analysis
It can be shown that the PDE we are trying to solve, with the accompanying boundary/initial
conditions, has the following exact solution:
ue(y; t) = y + exp( 2t) sin( y ) : (6)
It is evident that u(y; 1 ) = y, which corresponds to the steady Couette ow solution (the
linear velocity prole that we learnt about earlier in the course). Using this knowledge, we can
compute two kinds of errors as we numerically solve the equation.
The rst error is computed w.r.t. the exact solution, dened as:
E1(tn ) =
0
@ 1
N j
X
j
junj ue(yj ; tn )j2
1
A
1=2
: (7)
This is simply the root-mean-square (RMS) error which tells us how accurate our solution is.
Note, the boundary points are not used in calculating the error, since the solution there is
always imposed to be exact. Thus,N j = N 2. The second error is based on the steady-state
solution:
E2(tn ) =

1
N j
X
i
junj ue(yj ; 1 )j2
! 1=2
; (8)
2
which evidently tells us how close we are to the steady-state solution. Once the steady-state is
reached, we can stop the time marching, since the solution will not change (except for round-off
errors). For this reason, we will obtain the solution until the condition E2 < is satisfied, where
is the allowed tolerance in the round-off error. For the current assignment, use = 10−7.
Make sure all the variables in your code are double precision (64 bits), otherwise you might
have issues with convergence.
For the current assignment, you are required to write a code which implements the nu-
merical scheme to obtain the solution at the grid points and time steps until the steady-state
solution is reached, as prescribed by the condition E2 < . Use N = 11, 21, 41, 81, and for each
case use ∆t = 0.001, 0.003, 0.01, 0.03, 0.1. Perform a detailed assessment on the behavior of
the error E1 and E2 as a function of both spatial and temporal resolution. Also assess how the
numerical solution compares with the exact solution for these cases. For instance, you should
investigate things like how the error(s) behaves as a function of N (or ∆y) and ∆t; how it
changes as a function of total number of time steps (or time itself) for various cases, or in
general, anything else that can think of. A portion of the total points will be awarded on the
basis of novelty of the analysis and the conclusions drawn from it.
Things to submit
1. The source code which solves the problem, with an accompanying readme file with clear
instructions on how to compile and run the code. The entire source code, with all
functions, subroutines, etc. should be in one file (which can be compiled and executed).
Everyone should write their own code. All the relevant parameters of the problem, e.g.
N , ∆t, , etc. should be changeable at the beginning the code.
2. A text/data file containing the sample solution corresponding to the case N = 21 and
∆t = 0.003, and an accompanying text/data file containing the errors – please do not
submit excel or word files. The first ‘solution’ file should contain six columns, viz. time
step number, time, y-location, numerical solution, exact solution (the columns should be
clearly readable). Whereas, the second ‘error’ file should have four columns, viz. time
step number, time, E1 error, E2 error. Please make sure that these files are formatted
to be human-readable, i.e., the columns are structured and clearly separated. Also, do
not truncate the output in any way, i.e., all the significant digits should be written out.
3. A project report (as a PDF file) discussing all the various cases/results, with appropriate
figures/graphs. You are free to present the analysis in the manner you want, but the
report should be no longer than 8 pages (standard formatting). The first page should
briefly outline the problem and the last page should contain a conclusion section, where
you summarize your findings (you do not have to fill the entire page – more is not always
better). Essentially, you have 6 pages to present your findings in a concise manner, but
at the same time ensure all the important and relevant results are conveyed clearly.
3
Problem 2 – Steady heat equation
In this problem, you will solve the steady heat equation in two dimensions. For this case, the
governing PDE for temperature T (x, y) is given by the Laplace equation
r 2T =
∂2T
∂x2
+
∂2T
∂y2
= 0 . (9)
The solution domain is rectangular and defined by 0 x Lx , 0 y Ly . To make things
easy, let’s assume Lx = Ly = L. The boundary conditions for the problem are
T (0, y) = T0 , T (L, y) = 0 (10)
∂T
∂y
(x, 0) = 0 ,
∂T
∂y
(x, L) = 0 (11)
Methodology
The first step here is to also appropriately non-dimensionalize the equations. You can do so
by using T ∗ = T/T0, x∗ = x/L and y∗ = y/L. Henceforth, the superscript ∗ is dropped for
convenience. To discretize the domain, we will assume uniform spacing in each direction, but
we will use different number of grid points in each direction. In x-direction, Nx points are used,
such that grid spacing is ∆x = 1/(Nx 1), whereas in y-direction Ny points are used, such
that grid spacing is ∆y = 1/(Ny 1). The coordinate of a grid point is given by xi = (i 1)∆x
and yj = (j 1)∆y, for i = 1, 2, ...Nx and j = 1, 2, ...Ny .
Next, the equations are discretized using 2nd order central differences for both derivatives:
Ti+1;j 2Ti;j + Ti−1;j
(∆x)2
+
Ti;j +1 2Ti;j + Ti;j −1
(∆y)2
= 0 . (12)
This can be rearranged to be written as
Ti+1;j (2 + 2r)Ti;j + Ti−1;j + rTi;j +1 + rTi;j −1 = 0 , with r =
(∆x)2
(∆y)2
. (13)
The above system of equations can once again be solved by writing them in a vector form
Ax = b, where x is a vector containing Ti;j and b is a vector which will receive contributions
from the boundary conditions. The matrix A in this case will have five nonzero diagonal
bands, whose location will depend on how you define the vector x. Note, when implementing
boundary conditions the above equation will be somewhat modified and unlike in problem 1,
all the entries in the same band in the matrix A will not necessarily be the same.
Implementing Dirichlet boundary conditions is straightforward and similar as problem 1.
For Neumann boundary conditions, you will have to use a finite difference approximation. For
instance, the boundary condition ∂T/∂y(x, 0) = 0 can be implemented by using first order
forward scheme in the following manner
0 =
∂T
∂y
(x, 0) =
∂T
∂y
(xi , y1) =
Ti;2 Ti;1
∆y
(14)
4
which essentially gives Ti,1 = Ti,2, allowing you to write the value of T at this boundary in
terms of an unknown interior point. A similar result can be obtained for the other boundary
at (x, Ly), but be mindful that this will require a backward approximation.
After forming the matrix A and the vector b, you can once again simply invert the matrix
A to obtain the entire solution in one go. You can use direct inversion or you can use other
indirect algorithms for sparse matrices; once again, you can use standard libraries for this.
Note, Thomas algorithm cannot be used anymore.
Exact solution and error analysis
For this problem, the exact solution is
T e(x, y) = 1− x . (15)
Since there is no time marching, the RMS error in this case is simply defined as
E =
0
@ 1
NiNj
X
i,j
|Ti,j − T e(xi, yj)|2
1
A
1/2
, (16)
where Ni = Nx − 2 and Nj = Ny − 2, since once again the boundary points are not counted.
You are required to write a code which implements the methodology described above and
obtain the solution and the error. Try Nx = 26, 51, 101, 201, 401 and for each case use Ny =
26, 51, 101, 201, 401 (i.e., a total of 25 runs). Once again, assess how the numerical solution
compares to exact solution and the error changes with number of grid points/grid spacing.
Also consider how cases with Nx 6= Ny compare to cases with Nx = Ny.
Things to submit
1. The source code, with an accompanying readme file with clear instructions on how to
compile and run the code. The entire source code, with all functions, subroutines, etc.
should be in one file (which can be compiled and executed). All the relevant parameters
of the problem, e.g. Nx, Ny, etc. should be changeable at the beginning the code.
2. A text/data file containing the sample solution corresponding to the case Nx = Ny = 101
together with the error. The file should contain 6 columns, viz. x-location, y-location,
numerical solution, exact solution and the RMS error (the columns should be clearly
readable). Again, please make sure that these files are formatted to be human-readable
(and no excel/word files).
3. A project report (as a PDF file) discussing all the various cases/results, with appropriate
figures/graphs. Once again, you are free to present the analysis in the manner you want,
but this report should be no longer than 6 pages (standard formatting). The first page
should briefly outline the problem and the last page should contain a conclusion section,
where you summarize your findings – giving you 4 pages to present your findings.
5
Bonus
For bonus points, redo problem 2 using a a fourth order central scheme to approximate the
second derivatives. Note, this will change the stencil from 3 to 5 points in each direction.
However, for grid points right next to the boundary, using a 5 point stencil requires an ad-
ditional point which will be outside the boundary. To avoid this issue, still use second order
central approximation for points just inside the boundary. Instead, to compensate for reduced
accuracy near the boundary, utilize a second order approximation to impose the Neumann
boundary conditions. For instance, at the boundary (x, 0), you will use a second order forward
approximation, which will allow you to obtain Ti,1 in terms of Ti,2 and Ti,3 (as opposed to only
Ti,2 earlier in problem 2). The submission requirement for this bonus problem is the same as
problem 2.