# CFD Code Snippets—1: 1D transient conduction, FTCS, Dirichlet BCs

This semester, I am teaching an elective UG course on CFD. For illustrating ideas concerning the title problem, I began by borrowing the Python code written by Prof. Lorena Barba (general URL: [^]; the step concerning diffusion: [^] )), and then altered it a bit. I am sharing my modified version here.

In particular, here is a summary of my changes: (i) The Dirichlet BCs are now applied (so that time-stepping can go on for an indefinitely long period), (ii) graph plotting is done right within the time-marching loop, (iii) three different initial conditions can be used, etc. Variable names have been changed, and in fact, even comments have been expanded/corrected. (For instance, I use the term “thermal diffusivity” instead of“viscosity,” and use the actual material property value for it.)

This code can be used for the laboratory component (aka “Term-Work”) for the final-year elective course on CFD, in the S. P. Pune University.

This is the same University that denies me approval as a Professor, now because they think that I don’t have the required amount of work-experience (i.e. 10 years of both industrial, research and teaching experience put together).

The thing is, I didn’t always work for companies with annual turnover of 15 crores or more.

According to the PR guy of my employers, this university is now asking for my previous employers (e.g. startups in India/abroad) to certify in writing on their original letterheads not only the fact that they had employed me, but also to identify whether they had had 15 crores or more of annual turnover during the times when they had employed me or not. The rule is new, but applies retrospectively, i.e., it is supposed to go back even to the late 1990s or early ’00s times, when I worked in the software startups, sometimes with barely 1 or 2 employees too in their initial stages.

Most of these companies are by now wound up (e.g. QueriSoft grew to 500+ people whereas I was may be 13th/14th employee in Pune; nirWANa Inc. and Imagine Tech. where I was the fourth employee, and the second employee etc.—with nirWANa growing to US $1 Million in VC funding (coming from Kanwal Rekhi) but only after I left (at which point it had grown to 60+ employees)). I cannot have their certificates either because they can no longer be contacted, or even if the VC funding was US$ 1 Million, the sales turnover was not, etc.

Thus, the quantum of my total experience, according to the certifiable madness of the “authorities” at the S. P. University of Pune, suddenly dwindles down.

This rule, thus, allows the S. P. University of Pune to happily debar me from becoming a Professor.

It also allows employers to employ people like me not following the UGC scale, but on a consolidated, lower-level, salary.

I do not know precisely who among the two parties wishes to emphasize the rule more, to short-sell me. However, here, I am inclined to believe that it is the former. After all, it has been the government-run COEP and the (now S.P.) University of Pune who have happily supported people like Dr. Ashok Ghatol and Dr. (!) G. K. Kharate.

If professors who have worked nowhere else but in the private engineering colleges affiliated to the university rules are to be applied the same rule, then many of them would disqualify. This contradiction apparently does not matter to the University “authorities.” The rule is retrospective in action, but only for the new applicants, as per what our PR gentleman told me.

Anyway, I am going to share my code, but let me do so with a forthright moral denunciation of those office-bearers of the S. P. Pune University (and/or Maharashtra government) who are (or have been) responsible for framing such idiotic, harassing, anti-life and overall morally putrid rules. (For all I know, there could even be a casteist angle to it all; after all, it’s S. P. Pune University!)

Anyway, here is the Python code. It’s been written in v. 2.7 of Python. You will also need numpy and matplotlib installed. Feel free to use any IDE; I use PyCharm Community Edition.

#################################################################
# 1D Transient Conduction, using FTCS, with Dirichlet BCs.
# Based on the code by Prof. Lorena Barba (12 Steps to NS).
# Shared with an explicit identification of how morally putrid
# the S. P. Pune University is.
# Version: 27-Jan-2016.

#################################################################
# Importing Routines
# Use numpy arrays for efficiency
import numpy
# Use a neat plotting library
from matplotlib import pyplot  as plt

#################################################################
#  Program Data. All units are in SI
#
# Domain size, in meters
domLen = 0.1
# Number of cells
nc = 20
# Number of nodes
nx = nc + 1
# Compute the width of each cell
dx = domLen / nc
# Thermal diffusivity. This value is for Alumninium, in m^2/s
alpha = 8.418e-5
# Physical time over which the simulation is to be run
totalTime = 10.0
# Temperatures in degree Celcius, to use in constructing
# the initial and boundary conditions
tempLevel_1 = 0.0
tempLevel_2 = 100.0

# Parrot out the main input data
print("Length of the domain, in m: %lf" % domLen)
print("Number of cells: %d. Number of nodes: %d" % (nc, nx))
print("Therefore, the cell-width, in m: %lf" % dx)
print("Thermal diffusivity, in m^2/s: %lf" % alpha)
print("Duration (in physical seconds) for the simulation: %lf" % totalTime)
print("Temperature levels (deg. C): Level 1: %lf, Level 2: %lf" % (tempLevel_1, tempLevel_2))
# Main program data over
#################################################################

# Compute the duration of each time-step, and the number of time-steps required
# Stability criterion: r = \alpha \dfrac{\Delta t}{(\Delta x)^2}
# r <= 0.5 for stability
r = .35
print("Assumed stability criterion value: %lf" % r)
# Time-interval between steps, in seconds
dt = r * dx ** 2 / alpha
print("Duration of each time-step in s, after applying the stability criterion: %lf" % dt)
# Compute the total no. of time-steps (an integer) to fit within the given time-interval
nt = int(totalTime / dt)
print("Number of time-steps required to step through: %d" % nt)
# No. of time-steps to skip for plotting
stepsToSkipInPlotting = 10

# Create the numpy arrays, initialized to zero, for holding
# the primary variable values (here, temperature)

# The main array, holds the values after time-marching through one step
T = numpy.zeros(nx)
# A temporary array, holds the values at the previous time-step
Tn = numpy.zeros(nx)
print(T)
#################################################################
#  Setting up of the initial condition:
#
#  Uncomment the required IC and keep others commented

#########################
# IC1: A single pulse in the middle
# TL = tempLevel_1
# TR = tempLevel_1
# T[ nx/2 ] = tempLevel_2

#########################
# IC2: A step function in the middle
# TL = tempLevel_2
# TR = tempLevel_2
# T[ nx/3.0 : 2.0*nx/3.0  ] = tempLevel_1

#########################
# IC3: Left half at tempLevel_1; right half at tempLevel_2
TL = tempLevel_1
TR = tempLevel_2
T[0: nx / 2] = TL
T[nx / 2 + 1: nx] = TR

#
# Setting up of the initial condition is over
#################################################################

# Create the x-axis vector, for plotting purposes
#   numpy.linspace(start, stop, num, endpoint=True, retstep=False, dtype=None)
#   Returns evenly spaced numbers over a specified interval.
xAxis = numpy.linspace(0, domLen, nx)
plt.title("1D Transient Conduction, using FTCS")
plt.xlabel("Distance, in m")
plt.ylabel("Temperature, in deg C")
plt.plot(xAxis, T, 'r')

print("Initial condition is:")
print(T)

#################################################################
#  Time-stepping, using FTCS
#
# Multiplying factor to be used at each time-step
MF = alpha * dt / dx ** 2

# Python range(nSize) returns an array of size nSize
# Thus, range(3) returns [0,1,2]
for n in range(nt):
# We save the values of T at time n, to Tn,
Tn = T.copy()
# and apply the Dirichlet boundary conditions
# to the solution at the time-step n
Tn[0] = TL
Tn[nx - 1] = TR

# Then, we advance time to (n+1) and set T with it
# A note on the range() function, and the index values.
# -- The index to the right-most node in the arrays
#    T and Tn is: nx-1. If we specify range( nx-1 ),
#    then i will at most become _equal_ to nx-2
#    i.e., one less than what we specified.
for i in range(1, nx - 1):
# This is the FTCS algorithm
T[i] = Tn[i] + MF * (Tn[i - 1] - 2.0 * Tn[i] + Tn[i + 1])

# We skip stepsToSkipInPlotting' number of steps, in between two plots
if n and not (n % stepsToSkipInPlotting):
# Uncomment the following line for erasing out any earlier plots
# plt.clf()
# Uncomment the following line if you want to print T to console
# print( T )
plt.plot(xAxis, T, 'k')
plt.pause(0.1)

# The solution after nt number of time-steps
print("The solution after %d steps i.e. %lf physical seconds is: " % (nt, nt * dt))
print(T)
plt.plot(xAxis, T, 'b')
plt.show()

print("done")
#################################################################
#################################################################


Some things to try:

1. Try different values of the parameter <code>r</code>. For instance, set <code>r = 0.55</code> and see the instability develop.
2. Try mesh refinement.
3. Try changing material properties.
4. Try different initial conditions. For instance, try triangular profile, or a sine-wave profile.
5. Try setting a higher temperature, and keep the r’ parameter very close to 0.5. See whether it is the initial condition of a triangle or of a sharp pulse which first develops instability. Try to understand the fact that the von Neumann stability analysis is only linear, and therefore provides only an initial estimate for the stability criterion.
6. Compare in detail with the analytical solutions. For instance, the analytical solution for the initial triangular profile is worked out in Kreyszig.
7. If inclined towards programming, try using fft routines to work out analytical solutions, have the program save both the FTCS and the FFT-based solutions to output files, and then write short scripts to plot the errors as functions of, e.g., the shape of the initial profile, space (location), time-step, mesh-refinement, `r’ parameter, etc.

[May be one more editing pass, sometime later this week or so.]

[E&OE]