Micro-level water-resources engineering—8: Measure that water evaporation! Right now!!

It’s past the middle of May—the hottest time of the year in India.

The day-time is still lengthening. And it will continue doing so well up to the summer solstice in the late June, though once monsoon arrives some time in the first half of June, the solar flux in this part of the world would get reduced due to the cloud cover, and so, any further lengthening of the day would not matter.

In the place where I these days live, the day-time temperature easily goes up to 42–44 deg. C. This high a temperature is, that way, not at all unusual for most parts of Maharashtra; sometimes Pune, which is supposed to be a city of a pretty temperate climate (mainly because of the nearby Sahyaadris), also registers the max. temperatures in the early 40s. But what makes the region where I currently live worse than Pune are these two factors: (i) the minimum temperature too stays as high as 30–32 deg. C here whereas in Pune it could easily be falling to 27–26 deg. C even during May, and (ii) the fall of the temperatures at night-time proceeds very gradually here. On a hot day, it can easily be as high as 38 deg C. even after the sunset, and even 36–37 deg. C right by the time it’s the mid-night; the drop below 35 deg. C occurs only for the 3–4 hours in the early morning, between 4 to 7 AM. In comparison, Pune is way cooler. The max. temperatures Pune registers may be similar, but the evening- and the night-time temperatures fall down much more rapidly there.

There is a lesson for the media here. Media obsesses over the max. temperature (and its record, etc.). That’s because the journos mostly are BAs. (LOL!) But anyone who has studied physics and calculus knows that it’s the integral of temperature with respect to time that really matters, because it is this quantity which scales with the total thermal energy transferred to a body. So, the usual experience common people report is correct. Despite similar max. temperatures, this place is hotter, much hotter than Pune.


And, speaking of my own personal constitution, I can handle a cold weather way better than I can handle—if at all I can handle—a hot weather. [Yes, in short, I’ve been in a bad shape for the past month or more. Lethargic. Lackadaisical. Enervated. You get the idea.]


But why is it that the temperature does not matter as much as the thermal energy does?

Consider a body, say a cube of metal. Think of some hypothetical apparatus that keeps this body at the same cool temperature at all times, say, at 20 deg. C.  Here, choose the target temperature to be lower than the minimum temperature in the day. Assume that the atmospheric temperature at two different places varies between the same limits, say, 42 to 30 deg. C. Since the target temperature is lower than the minimum ambient temperature, you would have to take heat out of the cube at all times.

The question is, at which of the two places the apparatus has to work harder. To answer that question, you have to calculate the total thermal energy that has be drained out of the cube over a single day. To answer this second question, you would need the data of not just the lower and upper limits of the temperature but also how it varies with time between two limits.


The humidity too is lower here as compared to in Pune (and, of course, in Mumbai). So, it feels comparatively much more drier. It only adds to the real feel of a real hot weather.

One does not realize it, but the existence of a prolonged high temperature makes the atmosphere here imperceptibly slowly but also absolutely insurmountably, dehydrating.

Unlike in Mumbai, one does not notice much perspiration here, and that’s because the air is so dry that any perspiration that does occur also dries up very fast. Shirts getting drenched by perspiration is not a very common sight here. Overall, desiccating would be the right word to describe this kind of an air.

So, yes, it’s bad, but you can always take precautions. Make sure to drink a couple of glasses of cool water (better still, fresh lemonade) before you step out—whether you are thirsty or not. And take an onion with you when you go out; if you begin to feel too much of heat, you can always crush the onion with hand and apply the juice onto the top of your head. [Addendum: A colleague just informed me that it’s even better to actually cut the onion and keep its cut portion touching to your body, say inside your shirt. He has spent summers in eastern Maharashtra, where temperatures can reach 47 deg. C. … Oh well!]

Also, eat a lot more onions than you normally do.

And, once you return home, make sure not to drink water immediately. Wait for 5–10 minutes. Otherwise, the body goes into a shock, and the ensuing transient spikes in your biological metabolism can, at times, even trigger the sun-stroke—which can even be fatal. A simple precaution helps avoid it.

For the same reason, take care to sit down in the shade of a tree for a few minutes before you eat that slice of water-melon. Water-melon is nothing but more than 95% water, thrown with a little sugar, some fiber, and a good measure of minerals. All in all, good for your body because even if the perspiration is imperceptible in the hot and dry regions, it is still occurring, and with it, the body is being drained of the necessary electrolytes and minerals. … Lemonades and water-melons supply the electrolytes and the minerals. People do take care not to drink lemonade in the Sun, but they don’t always take the same precaution for water-melon. Yet, precisely because a water-melon has so much water, you should take care not to expose your body to a shock. [And, oh, BTW, just in case you didn’t know already, the doctor-recommended alternative to Electral powder is: your humble lemonade! Works exactly equivalently!!]


Also, the very low levels of humidity also imply that in places like this, the desert-cooler is effective, very effective. The city shops are full of them. Some of these air-coolers sport a very bare-bones design. Nothing fancy like the Symphony Diet cooler (which I did buy last year in Pune!). The air-coolers locally made here can be as simple as just an open tray at the bottom to hold the water, a cube made of a coarse wire-mesh which is padded with the khus/wood sheathings curtain, and a robust fan operating [[very] noisily]. But it works wonderfully. And these local-made air-coolers also are very inexpensive. You can get one for just Rs. 2,500 or 3,000. I mean the ones which have a capacity to keep at least 3–4 people cool.(Branded coolers like the one I bought in Pune—and it does work even in Pune—often go above Rs. 10,000. [I bought that cooler last year because I didn’t have a job, thanks to the Mechanical Engineering Professors in the Savitribai Phule Pune University.])


That way, I also try to think of the better things this kind of an air brings. How the table salt stays so smoothly flowing, how the instant coffee powder or Bournvita never turns into a glue, how an opened packet of potato chips stays so crisp for days, how washed clothes dry up in no time…

Which, incidentally, brings me to the topic of this post.


The middle—or the second half—of May also is the most ideal time to conduct evaporation experiments.

If you are looking for a summer project, here is one: to determine the evaporation rate in your locality.

Take a couple of transparent plastic jars of uniform cross section. The evaporation rate is not very highly sensitive to the cross-sectional area, but it does help to take a vessel or a jar of sizeable diameter.

Affix a mm scale on the outside of each jar, say using cello-tape. Fill the plastic jars to some level almost to the full.

Keep one jar out in the open (exposed to the Sun), and another one, inside your home, in the shade. For the jar kept outside, make sure that birds don’t come and drink the water, thereby messing up with your measurements. For this purpose, you may surround the jar with an enclosure having a coarse mesh. The mesh must be coarse; else it will reduce the solar flux. The “reduction in the solar flux” is just a fancy [mechanical [thermal] engineering] term for saying that the mesh, if too fine, might cast too significant a shadow.

Take measurements of the heights of the water daily at a fixed time of the day, say at 6:00 PM. Conduct the experiment for a week or 10 days.

Then, plot a graph of the daily water level vs. the time elapsed, for each jar.

Realize, the rate of evaporation is measured in terms of the fall in the height, and not in terms of the volume of water lost. That’s because once the exposed area is bigger than some limit, the evaporation rate (the loss in height) is more or less independent of the cross-sectional area.

Now figure out:

Does the evaporation rate stay the same every day? If there is any significant departure from a straight-line graph, how do you explain it? Was there a measurement error? Was there an unusually strong wind on a certain day? a cloud cover?

Repeat the experiment next winter (around the new year), and determine the rate of evaporation at that time.

Later on, also make some calculations. If you are building a check-dam or a farm-pond, how much would be the evaporation loss over the five months from January to May-end? Is the height of your water storage system enough to make it practically useful? economically viable?


A Song I Like:

(Hindi) “mausam aayegaa, jaayegaa, pyaar sadaa muskuraayegaa…”
Music: Manas Mukherjee
Singers: Manna Dey and Asha Bhosale
Lyrics: Vithalbhai Patel

CFD Code Snippets—2: Transients in the Couette flow, Crank-Nicolson

Here is some Python code for modeling transients in the pure shear-driven Couette flow between two infinite horizontal flat plates.

Initially, both the plates are stationary. Then, at time t = 0, the upper plate is suddenly set into motion with a constant speed u = ue. Due to the no-slip boundary condition between the fluid and the upper plate, the topmost layer of the fluid, too, is set into motion with the same speed, i.e., u = ue. Due to viscosity (i.e. friction internal to the fluid) this topmost layer tugs on to the next layer to move in the same direction, too. (That’s sort of like how when you pull a book from the middle of the stack, the books below it tend to come out, too.) Thus, the next layer also begins to move. Then, this second layer, in turn, nudges the third layer to move. …

Yes, it is turtles—but not all the way down.

The bottom-most fluid layer feels the tug not only of the layer immediately above it, but also the cent-per-cent tug of the no-slip boundary condition of the solid plate just below it. And the lower plate remains stationary throughout.

To cut a long story short, if the situation were steady-state, i.e., if the upper plate were always moving (at a constant speed ue) and the bottom plate stationary, then the velocity profile developed would be in the shape of a straight-line.

However, if both plates are initially stationary, and then if only the upper plate is suddenly moved, then the velocity profile would be initially like the Greek letter `\Gamma‘. With the passage of time, the profile would first lose the sharp 90^{\circ} angle, and so become a bit soft at that corner. Then, the curvature at the corner would become more and more diffuse. Then, the curved profile would become less and less curvy, and finally, it would reach the straight-line of the steady-state.

The code below shows how.

The Couette flow is important because it is one of the very, very (, very, very, very) rare cases in which the Navier-Stokes equations can (at all) be solved “analytically.”

The reason is, for this kind of a Couette flow—the one driven purely by the shear force of the upper plate (i.e., when no pressure gradients are applied in the horizontal direction)—the Navier-Stokes equations reduce to nothing but (our dear old friend): diffusion equation!

We did model the transient diffusion equation the last time, here [^]. However, the discretization approach used back then was FTCS—which is an explicit approach. The one now being used is: the Crank-Nicolson method—which is an implicit approach. FTCS is only conditionally stable, whereas CN is (at least for the diffusion equation, and at least within all reasonable limits) unconditionally stable. (Check out the Flow Science’s take on the explicit vs. implicit issue, here [^].) That makes it interesting.

The code once again comes with more comments on more issues than you want to read and/or know about.

The particular test case currently hard-coded is the one from John D. Anderson, Jr.’s celebrated text [^]. You can compare the output values produced by this program with the listing given in Anderson’s text (Table 9.1, p. 430.)

The Python code here even displays a MatPlotLib plot which you can directly compare with the fig. 9.4 in Anderson’s text (p. 431).

Anyway, here is the code, without further ado:

#################################################################
# Transients in the pure shear-driven Couette flow between two
# infinite flat plates, modeled using the Crank-Nicolson method.
# Here we generally follow the treatment and terminology as
# used in the text by John D. Anderson, Jr. However, please
# note, unlike in the book, here we use the _zero_-based
# index, _not_ 1-based. 
# Code written by and copyright (c)Ajit R. Jadhav.
# Version: 03-Feb-2016.
#################################################################

# Importing of modules and routines
# We use numpy arrays for efficiency
import numpy
# We use a neat plotting library
from matplotlib import pyplot as plt
from TDMASolver import TDMASolver

#################################################################
# Program Data. All units are in SI.
#
# Gap between the parallel plates, in m
D = 0.1
# Number of cells along the y-axis
nc = 20  # 4
# Number of nodes along the y-axis
N = nc + 1
# Compute the height of each cell
dy = D / nc
# Material properties
# The values below are for water
# Density, in kg / m^3
rho = 998.2
# viscosity, in Pa-s (or kg/(m.s))
mu = 8.90e-4
# speed of the moving (upper) plate, in m/s
ue = 0.05
# Compute the effective Reynolds' number
Re = rho * ue * D / mu
print("Calculated Re: %lf" % Re)

#################################################################
# Validation
# 1. The program output has been validated against manually
# worked out solutions for the first two time-steps, for these
# values of parameters: N=5, Re = 100.0, E = 1.0. For this set
# of values, and given the default truncation that happens in the
# Python print function, the convergence to the steady-state
# solution occurs at the 30th time-step.
#
# 2. Also validated with the case given in John D. Anderson's
# text. The values of parameters are: N=21, Re = 5000, nt =1.
# The program output matches the solution listing given on
# p. 430 in the book.
#################################################################

# For validation purposes, we will ignore the above-calculated
# Re, and instead set the variable with a simple value.
# Anderson's book uses 5000---quite close to the program data
# assumed above.
Re = 5000.0

# CN being implicit, has unconditional stability. However, for
# ease of comparisons, we define a parameter E anyway.
# $E \equiv \dfrac{(\Delta t)}{(Re)(\Delta y)^2})$
# Anderson's book uses E = 1.0
E = 1.0

# Number of time-steps to take
# Anderson's book shows graph up to 240 time-steps
nt = 241  # 36

print("Assumed Re: %lf, E: %lf, nt: %d" % (Re, E, nt))

# Calculate the time-step duration
dt = E * Re * dy * dy
print("Physical time elapsed in marching one computational time-step: %lf seconds." % dt)

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

# The main array, holds the values after time-marching through one step
u = numpy.zeros(N)
# A temporary array, holds the values at the previous time-step
un = numpy.zeros(N)

#################################################################
#  Setting up of the initial condition:
#
u[N - 1] = 1.0
print("Initial u: ", u)
#################################################################

# Create the y-axis vector, for plotting purposes
yAxis = numpy.linspace(0.0, 1.0, N)
plt.title("Transients in the pure shear-driven Couette flow,\nusing the Crank-Nicolson method")
plt.xlabel("Dimensionless horizontal velocity (u/ue)")
plt.ylabel("Dimensionless distance (y/D)")
plt.grid(True)
plt.plot(u, yAxis, 'r')

# Array containting the indices of time-steps at which the
# solution should be plotted:
ShowPlotInstants = [0, 2, 12, 36, 60, 240]

#################################################################
# Using the CN (Cranck-Nicolson) scheme, for the Couette Flow
#
# Compute A and B; they form a tri-diagonal matrix as below:
# B A
# A B A
#   A B A
A = -E / 2.0
B = 1.0 + E

print("A: %lf, B: %lf" % (A, B))

# a, b, c, d are the 1D arrays that define a tri-diagonal
# matrix. The tri-diagonal matrix in general looks like:
# a0|b0 c0 0 |  |x0| |d0|
#   |a1 b1 c1|  |x1|=|d1|
#   |0  a2 b2|c2|x2| |d2|
# a0 and c2 are not included in the tri-diagonal system, but
# the TDMASolver function assumes that these memory
# locations exist in the a and c arrays. Further, the
# algorithm employs the array indices in such a way that
# the first used element of the array a is a[1], not a[0],
# and similarly, the last used element of the array c is
# c[2], not c[1].
#
# For N-grid points, for the transient Couette flow problem,
# we get a system of N-2 X N-2 size, with the _global_ and
# _zero_-based _index_ going over [1, N-2]. Thus, if the
# number of grid points N = 5, we get [1,2,3].
a = numpy.zeros(N - 2)
b = numpy.zeros(N - 2)
c = numpy.zeros(N - 2)
d = numpy.zeros(N - 2)

# Initialize the a, b and c arrays
a[0:N - 2] = A
b[0:N - 2] = B
c[0:N - 2] = A

#################################################################
# Here begins the actual time-stepping...
nShowPlotIndex = 0
for n in range(nt):
    # We first save the known values of u at time n, to un,
    un = u.copy()
    # and apply the Dirichlet boundary conditions
    # to the solution at the time-step n
    un[0] = 0.0
    un[N - 1] = 1.0

    # The right hand-side vector on the CN is calculated here
    # Anderson calls the variable K
    # Note carefully the indices here. The indices used for Kj
    # are straight-forward from Anderson's book, except that here
    # they are redone as 0-based.
    # But the real change is on the index used on d. It is one
    # less, because here we don't at all store the first and the
    # last rows of the total system in it.
    for j in range(1, N - 1):
        Kj = (1.0 - E) * un[j] + E / 2.0 * (un[j + 1] + un[j - 1])
        d[j - 1] = Kj
    # For the same reason spelt above, note the indices: They are
    # respectively 1 less and 2 less.
    d[0] = d[0] - a[0] * un[0]
    d[N - 3] = d[N - 3] - c[N - 3] * un[N - 1]
    # print( d )

    # Realize, the solution array returned here has values only for
    # the interior nodes. Its elements do not include the two
    # boundary nodes.
    soln = TDMASolver(a, b, c, d)
    for j in range(1, N - 1):
        u[j] = soln[j - 1]

    # We skip plotting for all the time-steps except for those given
    # in the ShowPlotInstants array. The values are for the graph in
    # Anderson's text on p. 431
    if n == ShowPlotInstants[nShowPlotIndex]:
        nShowPlotIndex = nShowPlotIndex + 1
        # 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(u, yAxis, 'k')
        plt.pause(0.1)
        print("\nVelocity profile after time-step no. %d (or physical seconds: %lf)" % (n, float(n) * dt))
        print(u)

plt.plot(u, yAxis, 'b')
plt.show()

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

 

The above code calls on a TDMASolver() function written in a Python module of the same name. Here is the listing for that.

Initially, I tried writing my own TDMA (i.e. Thomas algorithm) code, but found that it was getting too time-consuming. So, I did the smart thing. I searched on the ‘net, found Ofan Tau’s code, here [^], and copy-pasted it!

One funny thing here was that Ofan in turn refers to some code on a Wikibooks page (here [^]), but while Ofan’s code is correct, the code on the Wikibooks page (at least the Python version of it) is not.

Incidentally, the best resource on TDMA aka Thomas’ algorithm that I found was one .PDF document by Prof. William T. Lee of the University of Limerick, Ireland, here [(.PDF) ^]. Keep it handy while going through this code.

#################################################################
# TDMASolver.py
#
# Solves a Tri-Diagonal Matrix system.
#
# a, b, c, d are the 1D arrays that define a tri-diagonal
# matrix. The tri-diagonal matrix in general looks like:
# a0|b0 c0 0 |  |x0| |d0|
#   |a1 b1 c1|  |x1|=|d1|
#   |0  a2 b2|c2|x2| |d2|
# a0 and c2 are not included in the tri-diagonal system, but
# the TDMASolver function assumes that these memory
# locations exist in the a and c arrays. Further, the
# algorithm employs the array indices in such a way that
# the first used element of the array a is a[1], not a[0],
# and similarly, the last used element of the array c is
# c[2], not c[1].
#
# The code here was taken from Ofan Tau's blog post:
# http://ofan666.blogspot.in/2012/02/tridiagonal-matrix-algorithm-solver-in.html
# Ofan's blog post makes reference to:
# https://en.wikibooks.org/wiki/Algorithm_Implementation/Linear_Algebra/Tridiagonal_matrix_algorithm
# It so happens that the _Python_ code at the wikibooks URL is
# _wrong_, but Ofan's code is correct (though it does not check
# for stability or even just diagonal dominance.)
# I have mostly kept Ofan's code as is; just added one or two
# comments.
#
# BTW, a wonderfully illustrated (and as far as I can tell,
# correct) reference for understanding TDMA is a .PDF document
# by Prof. William T. Lee (of Industrial Mathematics unit,
# University of Limerick, Ireland):
# http://www3.ul.ie/wlee/ms6021_thomas.pdf
#
#  -- Ajit R. Jadhav
# Version 03-February-2016
#################################################################

import numpy


#################################################################
# This helper function was written by Ajit R. Jadhav
# It uses the a, b and c vectors to fill a full-size matrix.
# It allows for matrix visualization. To be used mainly for
# debugging purposes.
def FillM(DM, a, b, c):
    N = len(a)
    DM[0][0] = b[0]
    DM[0][1] = c[0]
    for i in range(N - 1):
        DM[i][i - 1] = a[i]
        DM[i][i] = b[i]
        DM[i][i + 1] = c[i]
    DM[N - 1][N - 2] = a[N - 1]
    DM[N - 1][N - 1] = b[N - 1]
    return DM


#################################################################
# This function was shamelessly lifted from Ofan's blog-post:
# http://ofan666.blogspot.in/2012/02/tridiagonal-matrix-algorithm-solver-in.html,
# accessed 02-February-2016
#################################################################
def TDMASolver(a, b, c, d):
    '''
    TDMA solver, a b c d can be NumPy array type or Python list type.
    refer to http://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm

    # a, b, c, d are the 1D arrays that define a tri-diagonal
    # matrix. The triadiagonal matrix in general looks like:
    # a0|b0 c0 0 |  |x0| |d0|
    #   |a1 b1 c1|  |x1|=|d1|
    #   |0  a2 b2|c2|x2| |d2|
    # a0 and c2 are not included in the triadiagonal system, but
    # the TDMASolver function assumes that these memory
    # locations exist in the a and c arrays. Further, the
    # algorithm employs the array indices in such a way that
    # the first used element of the array a is a[1], not a[0],
    # and similarly, the last used element of the array c is
    # c[2], not c[1].

    # One test case is this:
    # a = [0, -1.0, -1.0, -1.0]
    # b = [4.0,  4.0,  4.0,  4.0]
    # c = [-1.0, -1.0, -1.0,  0.0]
    # d = [5.0,  5.0, 10.0, 23.0]
    # result: [2, 3, 5, 7]
    '''
    # number of equations
    nf = len(a)
    # print( "In TDMASolver.py: Size of the system: %d" % nf)

    # Make a copy the arrays
    ac, bc, cc, dc = map(numpy.array, (a, b, c, d))

    ##############################################
    # Make the matrix strictly upper triangular
    for i in range(1, nf):
        mc = ac[i] / bc[i - 1]
        bc[i] = bc[i] - mc * cc[i - 1]
        dc[i] = dc[i] - mc * dc[i - 1]

    ##############################################
    # The back-substitution pass:
    # Note the meaning of the -1 index in Python.
    xc = ac
    xc[-1] = dc[-1] / bc[-1]

    for j in range(nf - 2, -1, -1):
        xc[j] = (dc[j] - cc[j] * xc[j + 1]) / bc[j]

    # DM = numpy.zeros((nf,nf))
    # DM = FillM(DM, a, b, c)
    # print(DM)
    # del DM

    # delete variables from memory
    del bc, cc, dc

    return xc

Some things to try:

  1. Change E to 5.0, 10.0, and then even to 4000.0, and see if you get a solution like the one given in Anderson’s text (fig. 9.5, p. 934)
  2. Remove the hard-coded value of Reynolds’ number, and start using the calculated values for the same. Then, play with the physical parameters: use material properties for honey, glycerine, different SAE grades of motor oils, etc., and see if you get reasonable-looking “settling-in” time periods to approach steady-state or not.
  3. Also change the separation between the flat plates.
  4. What would happen if pressure gradients are applied to the fluid? In what way will the governing differential equation change? The velocity profile? Could you handle it using FDM?
  5. And, oh, yes! Do continue thinking about the SPPU (Savitribai Phule Pune University), and its past, present and future “authorities.”

[As usual, this post could change—improve or deprove—without any prior [or posterior] notice.]

[E&OE]

 

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.
# Code written by and copyright (c)Ajit R. Jadhav.
# 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
# About the linspace function
#   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]

The Mechanical-vs-Metallurgy “Branch-Jumping” Issue—Part III: A Couple of Stunning Papers

0. The earlier posts in this series may be found here [^] and here [^].

In this post, I am going to talk about two stunning papers.

1. The First Paper:

Before coming to the paper itself, let me discuss its background a bit.

The paper has a topical relevance to the current news; its subject matter in a way refers to those most unfortunate accidents on the Mumbai-Pune Expressway.

In the recent times, there has been a spate of one particular kind of accidents on the Expressway: a tire of a car (or of a tempo or even a container truck) suddenly bursts open, the car gets out of control, crosses over the median and lands in a lane on the opposite side of the traffic, and ends up having a head-on collision.

People have been discussing the matters, building awareness campaigns, trying to analyze causes. For instance, things like: over-speeding as the cause of the tires bursting open; the absence of a tall enough separating wall on either sides; and human factors such as loss of concentration, drunk-driving, etc.

It’s a fact that on the Expressway most every car exceeds the speed limit (of 80 km/hr). You can see many amateur YouTube videos of passenger buses and ordinary passenger cars registering 120–140 km/hr as a routine. BTW, while watching these videos, notice that in India we drive on the left-hand side, and so, the fastest traffic is (usually!) to be found in the rightmost lane.

When a car unfortunately develops a burst tire and lands up on the other side of the road, not only is it the fastest moving car in the forward direction but also the cars in whose path it lands up after crossing over the median, also happen to be the fastest moving ones on the opposite side. Thus, collisions can possibly occur at, say, 100 + 120 = 220 km/hr (even after discounting the reduction in speed during the crossover, say, from 140 to 100 km/hr), making the situation far more grave. Contrast this situation with that of a car moving at 80 km/hr in the left-most lane. Even if the driver loses control and the car veers on the left-hand side, it still collides on to the protective railing (having quite a bit of capacity to absorb the impact energy) at only, say, 60–70 km/hr. The difference of speeds in a head-on collision is, thus, three times over. Little wonder the percentage of fatalities in this kind of accidents has been so high.

Then, there are issues like the lights of the cars on the opposite side shining directly in your eyes, because there are no shades, no dividing walls that are high enough.

Though I am not at all an expert here, another issue I am toying with is the development of a peculiarly wavy kind of surface due to the combined effects of rolling and vibrations,  especially as caused by heavy weight container trucks travelling on a road. A similar effect can be seen also on the railway tracks, where a wavy profile gets developed. However, in the case of the rails, due to the homogeneity of the material, it becomes a somewhat simpler problem to analyse within a mechanics framework. Even then, let me tell you, it’s a very complex problem to analyse. IIT Kanpur, for instance, has been conducing research sponsored by railways on this problem for many years, and one could say that they are still in a somewhat early modelling stage. Now, when it comes to a steel-reinforced cement-concrete road, it becomes, from the applied mechanics viewpoint, a much more complex problem to model. When the road is made of concrete, I am guessing, the peaks of the wavy profile could perhaps be a bit a too sharp, thereby introducing enhanced high frequency components in the car-road interactions, vibrations and what effectively are small-scale but high-speed impacts. Now, if the tires are worn out and are conventional (i.e. if they carry tubes and are not of the tubeless variety), and if a driver has been over-speeding for a long time thereby overheating them, especially on this kind of a more sharply wavy surface, there could be increased chances of the tube of the tire bursting open. (The tubeless tires have a tendency to more gracefully handle an event of a similar kind; the running out of the air or nitrogen is less intense and not so catastrophic.)

So, people have been talking of things like that, and, as you can see, I too have joined them, wondering aloud about what possibly could be the causes of such accidents.

Against this backdrop, here is a paper I found in a recent Web search:

N. P. Dharmadhikari, D. C. Meshram, S. D. Kulkarni, S. M. Hambarde, A. P. Rao, S. S. Pimplikar, A. G. Kharat, and P. T. Patil (2010) “Geopathic stress: a study to understand its nature using Light Interference Technique,” Current Science, vol. 98, no. 5, pp. 695–697

I must say that the paper is simply stunning.

Check it for yourself, and if you think otherwise, do drop me a line. I am sure, though, that you, too, would find it simply stunning. Especially noticeable is the level of rigor in both experimentation and the kind of data that have been gathered therein. Do have a look.

I gather that “Current Science” is “the [sic] leading interdisciplinary science journal from India […] started in 1932 by the then stalwarts of Indian science such as CV Raman, Birbal Sahni, Meghnad Saha, Martin Foster and S.S. Bhatnagar.” I also gather that it is indexed in “Web of Science, Current Contents, Geobase, Chemical Abstracts, IndMed and Scopus” [^].

BTW, the fact that the journal is indexed by so many services, esp. Scopus, goes to satisfy the requirements clarified by Professor Dr. T. P. Singh (Director, Symbiosis Institute of Technology). Recently, in an informal chat, he had mentioned to me that he cannot hire me because I do not have two journal papers already published, and so, I cannot become a PhD guide at Symbiosis from day 1. He then added that it is now his requirement to hire only those people as Full Professors who have two papers published in Scopus/similar indexed journals already, so that they can become PhD guides right from day 1 (his emphasis) and then, he continued in the same conversation, that during the last semester (when I too had applied at Symbiosis) he had hired two people as Full Professors but without PhD degrees so they cannot act as PhD guides from day 1 (here, the emphasis is mine). BTW, at the time that he hired those two people, I, too, had applied for the same post (having PhD already in hand), but was not called for interview at all. … Yes, sometimes Prof. Singh somehow has a charming way of reminding me of the late Jaspal Bhatti, an engineer.

Anyway, coming back to this paper, please note, engineering college professors from JSPM, MIT and Sinhagad have participated in it. There are as many as 7 authors to this paper. (Isn’t this fact stunning by itself?) They all, I suppose, have been approved by the UGC committees for recruitment of professors. In fact, this research began with a PhD granted by the University of Pune, in the first place.

BTW, if you are interested, here is another article which carries more or less the same level of rigor but it is definitely much more readable: [^].

2. The Second Paper:

Another paper, too, has a bit of a topical importance, at least as far as this blog is concerned, in the light of the recent discussion on the Fourier theory and all.

The paper in question is this:

G. K. Kharate, V. H. Patil and N. L. Bhale (2007) “Selection of mother wavelet for image compression on basis of nature of image,” Journal of Multimedia, vol. 2, no. 6, pp. 44–51

This paper, too, is stunning.

Why do I find it stunning?

In this paper, the authors experimentally study the extent to which different types of images can be compressed, using wavelets transforms. In particular, the extent of compression that is possible for a natural image (a woman’s portrait) is compared and contrasted with certain artificial images consisting of simple regular patterns such as horizontal/vertical/oblique lines. The conclusions reached by the authors consist of the following:

“The results demonstrate that for line-based image percentage of zeros is more for db1 as compared to other wavelets and more energy is retained. It shows that the loss of information is less hence the quality is better. Also for natural image percentage of zeros is more for db4 as compared to other wavelets that indicates more compression and the loss of information is less.

If the import of the conclusion is not immediately clear, here are a few passages from the main text of the paper that may help you in ascertaining the nature of the paper:

“It is observed that for periodic functions, Fourier analysis is ideal. However, wavelet transforms are not restricted to only the periodic function, but for any function, provided it is admissible.”

I was left wondering: How about the computational complexity of FFT vis-a-vis wavelets—particularly because the paper deals with compression efficiency? Reading further on, I then saw this:

“We suggest the minimum two steps for the wavelet decomposition. In first step, decompose the image by using any mother wavelet. On the basis of energy content of approximation and detail coefficients, the image is [first] classified as line image or natural image. Then select the proper mother wavelet for further decompositions as mentioned above.” [Emphasis in bold is mine.]

How about entropy? Well, this aspect is not as clearly indicated.

The Journal of Multimedia is, of course, an international journal published from the USA; it is indexed in many services, including Scopus [^].

Anyway, the first author is the current Dean of the Faculty of Engineering at the University of Pune. I mean to say: he was first awarded a doctorate in engineering by the University of Pune (under the guidance of A. A. Ghatol); he was selected by a UGC committee for a Full Professor’s position; he was also democratically elected by the professors coming from all the engineering colleges affiliated to the university of Pune; he now sits on the UGC interview committees that decide whether someone (e.g. me) is to be at all interviewed or not.

Do note that it’s not his the then PhD student (i.e. V. H. Patil), but Kharate himself who is the first author of this paper. Not so stunning. Not any more.

But, yes, the paper, overall, is stunning. Its conclusion simply has a way to force you to do a double take. Its very appearance is against what you expect intuitively—you never expect a paper of this kind to appear in a typical journal covering topics such as multimedia, image or signal processing, or wavelets analysis.

As you might know, the wavelets transform is similar to the Fourier transform though it has a categorically less computational cost involved with it as compared to FFT. Also, take a moment to think about what kind of connections could at all be possible between computational complexity, information entropy and the compression efficiency.

Now, reconsider the cited text and the conclusion of this paper.

Hopefully, now you can see why I say that you must read this journal paper, too.

* * * * *    * * * * *   * * * * *   

BTW, I have finally been rejected in the interviewing process at the JPBTI Phaneesh Murthy’s company, iGATE. The reason communicated to me (over the phone) by Asad Kadri this late evening had to do with the fact that I don’t have documentary evidence for having actually done jobs at three places: (i) Mukand (in 1983), and (ii) Thermax (in 1983–84); at both these places I was a trainee engineer; and then, (iii) QueriSoft (in 1995–96 i.e. about 15+ years ago), a company which is closed by now. BTW, I had arranged for a former colleague to informally verify the fact of my having worked with Querisoft, and anyway did have a xerox copy of the appointment letter for it. However, iGATE apparently thinks that their end-client, Rio Tinto, would find it unacceptable that all of my employment cannot be verified. In any case, they didn’t reply back to that email written by that past colleague of mine (who later on worked as a director in a large American MNC, and has now taken voluntary retirement from the IT field).

On my part, I find both Rio Tinto and iGATE to be very funny companies—tinted, as if they were, with a very dogged and out-of-the-world sort of a sense of humour. (They took three months, two interviews with the VP, three interviews with HR, and many telephonic calls and emails exchanges, in order to reach the conclusion that they cannot appoint me because I do not carry experience certificates for companies where I worked two decades ago and more.)

So, my job search continues, and something on the academia side might materialize—something that allows me at least, say, survival. I will keep you posted.

… But, don’t forget only to see those two papers, na!

* * * * *    * * * * *   * * * * *  

A Song I Like:
[For the above-mentioned reasons, the inclusion this time round of this section is, once again, purely on the whim.]

(Hindi) “machalti hui hawaa main chham chham…” (“hamare sang sang chalein Ganga ki leheren”)
Singers: Kishore Kumar [an unlikely choice, esp. given early ’60s], and Lata Mangeshkar
Music: Chitragupt
Lyrics: Majrooh Sultanpuri [?]

[E&OE]