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)1
#   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]

Conservation of angular momentum isn’t [very] fundamental!

What are the conservation principles (in physics)?

In the first course on engineering mechanics (i.e. the mechanics of rigid bodies) we are taught that there are these three conservation principles: Conservation of: (i) energy, (ii) momentum, and (iii) angular momentum. [I am talking about engineering programs. That means, we live entirely in a Euclidean, non-relativistic, world.]

Then we learn mechanics of fluids, and the conservation of (iv) mass too gets added. That makes it four.

Then we come to computational fluid dynamics (CFD), and we begin to deal with only three equations: conservation of (i) mass, (ii) momentum, and (iii) energy. What happens to the conservation of the angular momentum? Why does the course on CFD drop it? For simplicity of analysis?

Ask that question to postgraduate engineers, even those who have done a specialization in CFD, and chances are, a significant number of them won’t be able to answer that question in a very clear manner.

Some of them may attempt this line of reasoning: That’s because in deriving the fluids equations (whether for a Newtonian fluid or a non-Newtonian one), the stress tensor is already assumed to be symmetrical: the shear stresses acting on the adjacent faces are taken to be equal and opposite (e.g. \sigma_{xy} = \sigma_{yx}). The assumed equality can come about only after assuming conservation of the angular momentum, and thus, the principle is already embedded in the momentum equations, as they are stated in CFD.

If so, ask them: How about a finite rotating body—say a gyroscope? (Assume rigidity for convenience, if you wish.) Chances are, a great majority of them will immediately agree that in this case, however, we have to apply the angular momentum principle separately.

Why is there this difference between the fluids and the finite rotating bodies? After all, both are continua, as in contrast to point-particles.

Most of them would fall silent at this point. [If not, know that you are talking with someone who knows his mechanics well!]


Actually, it so turns out that in continua, the angular momentum is an emergent/derivative property—not the most fundamental one. In continua, it’s OK to assume conservation of just the linear momentum alone. If it is satisfied, the conservation of angular momentum will get satisfied automatically. Yes, even in case of a spinning wheel.

Don’t believe me?

Let me direct you to Chad Orzel; check out here [^]. Orzel writes:

[The spinning wheel] “is a classical system, so all of its dynamics need to be contained within Newton’s Laws. Which means it ought to be possible to look at how angular momentum comes out of the ordinary linear momentum and forces of the components making up the wheel. Of course, it’s kind of hard to see how this works, but that’s what we have computers for.” [Emphasis in italics is mine.]

He proceeds to put together a simple demo in Python. Then, he also expands on it further, here [^].


Cool. If you think you have understood Orzel’s argument well, answer this [admittedly deceptive] question: How about point particles? Do we need a separate conservation principle for the angular momentum, in addition to that for the linear momentum at least in their case? How about the earth and the moon system, granted that both can be idealized as point particles (the way Newton did)?

Think about it.


A Song I Like:

(Hindi) “baandhee re kaahe preet, piyaa ke sang”
Singer: Sulakshana Pandit
Music: Kalyanji-Anandji
Lyrics: M. G. Hashmat

 

[E&OE]

From the horses’ mouths

My first choice for the title was: “From the Nobel Laureate’s Mouth”; I had spotted only the opinion piece by Professor David Gross in yesterday’s Indian Express [^]. Doing the ‘net search today for the URI link to provide here, I found that there also were three other Nobel laureates, also joined by one Fields Medalists. And they all were saying more or less the same thing [^].

… That way, coming from a Marathi-medium schooling background, I had always had a bit of suspicion for the phrase “from the horse’s mouth.” It seemed OK to use in the news reports when, say, a wrong-doer admits his wrong. But purely going by the usage, I could see that the phrase would also be used in the sense: “from the top-gun himself,” or “from the otherwise silent doer himself.” This guess turns out to be right [^]. Further, since there were as many as five “horses” here, the word to be used would have to be in the plural, and if you say it aloud: “From the horses’ mouths” [go ahead, say it aloud, sort of like:“horseses” mouth) it really sounds perfect (for something to be posted on the ‘net).

So, that’s how comes the title.

As to the horses’ thoughts… Ummm…

[But please, please, give me just a moment to get back to the title again, and congratulate me for not having chosen a title like: “From Dave Himself.” You see, Professor David Gross had visited COEP in 2013, and I might have been, you know, within 50 meters of where he was sitting. I mean, of all places, in the COEP campus! Right in the COEP campus!! [^]. Obviously, you must compliment me for my sense of restraint, of making understatements.]

OK. As to their thoughts… Umm….

I think these guys are being way too optimistic. Also naive.

Without substantial economic reforms, I see no possibility of the Indian Science in general undergoing any significant transformation yet again. And substantial economic reforms aren’t happening here any more. In fact, no one is even talking about it, any more. [Check out Arnab’s hours, or Sardesai’s, or Dutt’s, if you want to find out what they are talking about. [I don’t, because I know.]]

It was the 1991 that could propel, say a Mashelkar into prominence several years later, and help transform the 70+ CSIR labs from something like less than 100 patents a year, to thousands of them per year—all within a matter of a few years [less than a decade, to be sure]. If the same momentum were to be kept, the figure should have gone up to at least tens of thousands of patents by the CSIR labs alone—and with a substantial increase in the share of the international patents among them. Ditto, for the high-quality international journal papers.

Why didn’t any of it happen? Plain and clear. The momentum created by the economic liberalization of the early 1990s has been all but lost. Come on, face it, 1991 was twenty-five years ago.

To an anthropologist, 25 years is like an entire generation! More than enough of a time to lose any half-hearted momentum (which, despite the hysterical Indian press, the liberalization in the early 1990s was).

It’s been years that we entered the staleness 2.0 of the mixed economy 1.0. Even today, the situation continues “as is,” despite a change of regime in New Delhi. Yes, even under “Modiji.” [I am quoting Professor Gross—I mean the word.]

But, yes, the five gentlemen were also being realistic: Each one of them emphasized decades.

Decades of sustained efforts would have to go in, before the fruits could begin to be had. [But you know that decades isn’t a very long period—just recall what was happening to India’s economy some two decades ago—in the mid ’90s.]

Talking of how realistic they actually were being, Haroche even pointed out the lack of freedom in China [obvious to any one outside of California], and its presence in Europe [I don’t know about that] and in India [yeah, right!].

But anyway, it’s nice to hear something like this being highlighted after an Indian Science Congress, rather than, say, “vimaanshaastra.”

Both happened during “Modiji”’s tenure. So what is it that really accounts for the difference? I have no idea. (It can’t be a “pravaasi” whatever, to be sure; they would be too busy booking the next Olympics-size stadium.)

Whoever within the organizers of the Congress was responsible for the difference, compliments are due to him. (Hindi) “der se kiyaa lekin kuchh achhaa hi to kiyaa.”

In the meanwhile, bring out your non-programmable desk calculators and do some exercises: 0.8 \times \dots, 2.7 \times \dots, 4.4 \times \dots and 2.1 \times \dots. Oh well, you will have to refer to the ‘net.

OK then. Find out also the R&D spending by, say, (i) Baba Ramdev’s pharmaceutical industries, (ii) the top or most well-established five industrial groups in India (Reliance, Tatas, Mittals, whoever…), and (iii) the top three (or five) Indian IT firms. Compare them to those in the advanced countries. Let your comparisons be comparable: pharma to pharma; oil, steel and engineering (and salt!) to oil, steel and engineering (and salt!); IT to IT [engineering IT to engineering IT]; overall (GDP) to overall (GDP).

And, never forget that bit about freedom. Don’t just count the beans “spent” on research. Think also about whether it is the government spending or the private spending, and where the expenditure occurs (in private universities, private labs, independently run government labs, public universities in a country with a past of a private control, etc., or in the in-service-pensioner’s-paradises with something like “laboratory” in their titles).


But why didn’t the “horses” cite any specific statistics about how many Indian students go abroad for their graduate studies, and choose to permanently settle there—their trends?

Obvious: Nobel and Fields laureates (and in fact any visiting dignitaries to any country (and in fact any visitors to a foreign country)) generally tend to be more polite, and so tend to make understatements when it comes to criticism (of that host country). That’s why.


A Song I Like:

(Hindi) “kahin naa jaa…”
Music: R. D. Burman
Lyrics: Majrooh Sultanpuri
Singers: Kishore Kumar, Lata Mangeshkar

[E&OE]

My answers to the Edge questions

I was supposed to come back and write a better version of the QM-related comment which I had made at Scott Aaronson’s blog (see my last post). However, I am not going to do that right away because of two reasons: (i) I think that leaving aside a typo or two, the comment isn’t in all that bad a condition, and (ii) I have begun thinking about some of the issues involved, and so would like to come back only after getting them straightened out (with a better order of presentation).

In the meanwhile, I decided to check out my blog-ideas folder (yes, I do “maintain” one), and found that I already had some material for a post.

It had so happened that while going through Scott’s post on the Edge question, especially the beginning of that post, an idea had struck me:

why not write my own one-line answers to all the Edge questions thus far?

So, I had gone to the Edge home page, and jotted down my answers in a .txt file. That was on 4th January 2016, between 4:17 PM–5:43 PM, IST (which means, one day before I wrote my comment on Scott’s blog). The answers were jotted down very, very rapidly—average time between reading a question and starting typing answer could be barely 10–15 seconds. That’s because we were very busy with the accreditation related work, remember? Though I don’t remember that particular day, I must have finished the day’s work and only then done this thing. Here I am going to mostly copy-paste its contents.

A couple of notes before we go over to my answers:

  1. I just copy-pasted the Edge questions in the order I saw them on their home page. The order of answers was, thus, determined more or less by their Web page layout. So, the answers are not in the chronological sequence of the years over which they were raised.
  2. I am reproducing the contents of a very whimsically written file (which was written, as I said, near the end of a long work-day), and without effecting any editing at all. Today, right within a week of writing the answers, I find that I want to change the answers to some of the questions. In any case, I would certainly want to make the answers more streamlined. But, today’s editing is minimal; whatever whimsicality that had crept in, I have decided to keep its original flavor as is.

So, here we go, with a copy-paste from that file:


2004 : WHAT’S YOUR LAW?

The front velocity in the linear diffusion is a finite quantity (and in the first-order analysis, it is a constant).

2003 : WHAT ARE THE PRESSING SCIENTIFIC ISSUES FOR THE NATION AND THE WORLD, AND WHAT IS YOUR ADVICE ON HOW I CAN BEGIN TO DEAL WITH THEM? – GWB

The nation (i.e. USA): I don’t care a hoot about. The position was reciprocally derived.
The world: More or less the same. More or less the same.
My advise (to any President of the USA): Embrace 100% pure capitalism.

2002 : WHAT IS YOUR QUESTION? … WHY?

In his evolution, exactly for how long has the modern man been in the junglee-like state? (Rough estimates go like: lakhs of years of the jungle state vs. 5700 years of civilization, if you want to believe in the Western/modern science.)

Why it’s important? From many viewpoints: What led to the dominance of reason implicit in civilization? To what extent are the biological structural features really important in the mind-body connection? In what subtle ways must we still be junglee-like?

2001 : WHAT NOW?

Revising my viewpoint concerning non-relativistic quantum mechanics.
Reaching a conclusion concerning the issue of “stress or strain—which one is more fundamental?”
Writing some toy CFD code

2001 : WHAT QUESTIONS HAVE DISAPPEARED?

Whether the following phenomena are real or not: (i) rebirth, (ii) tele-many things such as telepathy, mind-reading, etc. (iii) psychic attacks, (iv) the actual American depravity in practice. Also the question of where Osama was hiding. Answer: In a cantonment town in a mostly army-ruled country that is a friend of USA.

2000 : WHAT IS TODAY’S MOST IMPORTANT UNREPORTED STORY?

How telecommunications actually changed the attitudes of rural people in India

1999 : WHAT IS THE MOST IMPORTANT INVENTION IN THE PAST TWO THOUSAND YEARS?

I don’t answer such silly questions. Idea? I would have attempted; generalization is possible with ideas. But not with inventions.

1998 : WHAT QUESTIONS ARE YOU ASKING YOURSELF?

Are the reported incidents in which some people could predict the future events true? If yes, what could be the mechanism in the pycho-somatic i.e. yogic terms? What implication does it hold for free will (properly defined, most fundamentally, only as the choice to focus or to evade), and the nature of soul and of consciousness.

2015: WHAT DO YOU THINK ABOUT MACHINES THAT THINK?

It’s a stupid idea. Thinking requires a conceptual (i.e. volitional) consciousness, which in turn requires life.

2014: WHAT SCIENTIFIC IDEA IS READY FOR RETIREMENT?

None. Proper ideas belonging to proper science never retire, they get subsumed away.

There are any number of ideas that should be retired because they were neither proper ideas nor scientific. For instance, the idea of an infinite physical universe. Or, the idea that ascribing anything quantitative to the universe (as a whole) is epistemologically valid. Etc. Also, string theory.

2013 : WHAT *SHOULD* WE BE WORRIED ABOUT?

The String Theory.

More seriously, government control of science and education.

2012: WHAT IS YOUR FAVORITE DEEP, ELEGANT, OR BEAUTIFUL EXPLANATION?

The mathematical aspects of the non-relativistic quantum mechanics. (A proper theory for its physical aspects is not yet available—despite my own attempts to build one.)

2011: WHAT SCIENTIFIC CONCEPT WOULD IMPROVE EVERYBODY’S COGNITIVE TOOLKIT?

The idea that all physical theories should be local.

I wrote that down, and then realized that the question asks for a toolkit.

I think computational modeling as easily realizable using Python and its ecosystem. It should be a part of every scientist’s grooming up. (Python is already being used to introduce programming to the CBSE XI-XII students in India.)

2010: HOW IS THE INTERNET CHANGING THE WAY YOU THINK?

It hasn’t changed *my* way of thinking, much. I always was more conceptual, and always had relatively lower working memory (and therefore also lessser maths skills) as compared to the class-mates who I ended up competing or at least [being] compared with. Internet has made more people forgetful and overall more like me, simply because they can just google up a lot of the concrete information.

2009: WHAT WILL CHANGE EVERYTHING?

I have no idea. Also, I don’t think it really is important or even desirable to change everything.

2008: WHAT HAVE YOU CHANGED YOUR MIND ABOUT? WHY?

My view of quantum mechanics. I used to think that physicists were conceptually dumb. I now know that the original discoverers were not dumb—not even conceptually—but that their major flaw was either too much of careless[ness] about basic philosophical ideas, or too much of a diffidence towards bad philosophers.

Also, my view of Objectivists, Americans, and in fact, Western people in general. I no longer think that any of them could be fully relied on, practically speaking. They are not bad. [Here I meant morally bad.] In fact, they are pretty OK. They are just OK in the sense, they can be nice if you already are a practical success yourself. But if you are not a success, [they won’t bother to find that] hidden talent or spark in you. They are too crude and/or dumb and/or isolated, to be able to do that. Especially to an Indian-born. (They are rather like Brahmins in the Indian cities [and in the SF Bay Area]. The Brahmin-borns are not bad [by birth]. It’s just that a Maratha-born growing up rural areas cannot rely on them. He may receive acknowledgement of his talents from them, but not so reliably frequently to be a social norm. [Note today: It’s the elite vs non-elite phenomenon, really speaking, but the point is: I had never imagined as a young man that it would be as bad as it is, esp. in America.]

2007: WHAT ARE YOU OPTIMISTIC ABOUT?

That my new theory of QM will work out.

Actually, on most things, and speaking in terms of how I am habitually like, I am neither optimistic nor pessimistic. Habitually, if there is any one option towards which I do almost habitually get inclined, it is that: I like to know.

2006: WHAT IS YOUR DANGEROUS IDEA?

Dangerous in general? I have none. But if you mean those tiny cultural points of conversation that society people and people of culture (esp. writers of novels like Ayn Rand) use, you know, BS like: “he is too dangerous because he is too independent,” etc., then guess *my* most dangerous idea has been my sharing of the dimmed view of many people (Americans, Objectivists, Brahmins, army people, Indian University people, Westerners). Dangerous to me, that is.

2005: WHAT DO YOU BELIEVE IS TRUE EVEN THOUGH YOU CANNOT PROVE IT?

But I don’t always look for proofs, that’s the point. Validation is one thing, proofs is another. I do make sure to know the meaning and the roots of the concepts that I use. (And, sometimes, I have spent decades understanding their meaning and/or roots.)

But knowing the roots isn’t the same as proving, esp. to someone else, esp. to his satisfaction. You know how increasingly pointless and dangerous this enterprise gets.

Anyway, as the “cultural” kind of a fraud kind of an answer: My new theory of QM.

2016: WHAT DO YOU CONSIDER THE MOST INTERESTING RECENT [SCIENTIFIC] NEWS? WHAT MAKES IT IMPORTANT?

I can’t think of any. It has to be at least interesting, right? And, it has to be recent. Sorry. Can’t think of something that has both. At least not off the hand.


It would be fun knowing your answers.

It would also be fun writing my own answers afresh some time later on.


A Song I Like:

(Hindi) “dil lagaa kar hum yeh samajhe”

Singers (separate versions each): Mahendra Kapoor and Asha Bhosale
Music: C. Ramchandra
Lyrics: Shakeel Badayuni

[Another whimsicality of mine: Unlike most people whose opinions I usually find sensible, here I find that it is tough to decide which singer’s version is better. I mean to say, I don’t automatically find Asha’s version better. In fact, I like Mahendra Kapoor’s version [just] a shade better.]

[E&OE]