代写ME6043、代写C/C++, Python编程

日期：2022-11-27 10:27

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:

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

(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:

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:

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. (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)

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

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) = 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.