The first-order (i.e. one-way) wave equation in 1D is given by:

.

where is the transverse displacement now, not temperature, and is a constant.

The equation describes a pulse moving at a constant speed to the right. If , the pulse moves to the left. Gilbert Strang [^] describes this equation with terms like (momentary) “flash of light” and “wall of water.”

It is with this problem that Prof. Lorena Barba chooses to begin [^] her 12 Steps program [^]. However, I am not so sure whether it is better to begin the course with convection rather than with diffusion.

The differential terms in convection are all of the first-order whereas those in the diffusion are mixed: first-order in time and second-order in space. It would thus appear that convection is a simpler phenomena than diffusion. In a certain sense, it is.

However, it so happens that for a course on *CFD*, the linear convection is only a stepping stone to the nonlinear convection. Seen from this viewpoint, it’s better to first cover diffusion, and only then go over to convection. This way, it becomes easier for the student to compare and contrast the nonlinear convection from the linear convection.

Currently, one way or the other, I am constrained to teach only according to the SPPU syllabus. So, one way or the other, the first-order wave equation could have come only after the diffusion equation.

No, don’t give too much credit to the S. P. Pune University. The sequence wouldn’t be a conscious choice of the syllabus framers. The most probable scenario is that they must have just arranged and rearranged the material until they were satisfied that every committee-member’s suggestion was taken into consideration without any consideration of merit, and that all the six units carried more or less the same amount of material. Or, less likely, until the point that they themselves got tired of arguing with each other. In any case, *none* could possibly have thought of ordering the material according to a proper hierarchy.

Anyway, here we go with the Python code. (I hurriedly rewrote it afresh after the recent HDD crash, and so, comments and all could be refined a bit further. As to bugs: You be the tester.)

################################################################# ''' First-Order Wave Equation, using First-Order Upwind Scheme Code written by and copyright (c)Ajit R. Jadhav. Based on the code by Prof. Lorena Barba (12 Steps to NS) For theory, also see notes/resources on the 'net, such as: Gilbert Strang: http://ocw.mit.edu/courses/mathematics/18-086-mathematical-methods-for-engineers-ii-spring-2006/readings/am52.pdf Olver: http://www.math.umn.edu/~olver/num_/lnp.pdf Gardner: https://math.la.asu.edu/~gardner/wave.pdf and others. Version: 29-March-2016. ''' ################################################################# import numpy from matplotlib import pyplot as plt import random ################################################################# # Program Data. All units are in SI # # Physical time over which the simulation is to be run, in s totalTime = 2.0 # Domain size, in meters domLen = 1.0 # Signal velocity of the pulse, m/s c = 1.0 # Number of cells nc = 50 # Number of nodes nx = nc + 1 # Compute the width of each cell dx = domLen / nc # Vertical (i.e. transverse) wave displacements in m, to use # in constructing the initial and boundary conditions T_Level_1 = 0.0 T_Level_2 = 0.1 # Create the numpy arrays, initialized to zero, for holding # the primary variable values (here, transverse displacements) # 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) # Compute the duration of each time-step, and the number of time-steps required # Stability criterion: CFL \equiv c \dfrac{\Delta t}{\Delta x} # CFL <= 1.0 for stability, but the value of 1.0 here is special! # Try other values of CFL < 1.0, a few values > 1.0, and a value of exactly 1.0 # Explain the smoothening at CFL < 1.0 # Explain the instability at CFL > 1.0 # Explain the ``oddity'' of the result at CFL = 1.0 CFL = 0.25 print("Assumed stability criterion value: %lf" % CFL) # Given the above CFL value, back-calculate the time-step duration dt = CFL * dx / c 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) ################################################################# # Setting up of the initial condition: # # Uncomment the required IC and keep others commented ######################### # IC1: A single pulse in the middle # TL = T_Level_1 # TR = T_Level_1 # T[ nx/2 ] = T_Level_2 ######################### # IC2: A step function in the middle # TL = T_Level_1 # TR = T_Level_1 # T[ nx/3.0 : 2.0*nx/3.0 ] = T_Level_2 ######################### # IC3: Left half at tempLevel_1; right half at tempLevel_2 # TL = T_Level_1 # TR = T_Level_2 # T[0: nx / 2] = TL # T[nx / 2 + 1: nx] = TR ######################### # IC4: An arbitrary function random.seed(10) TL = T_Level_1 TR = T_Level_1 for i in range(nc/5,nc/5+5): T[i] = T_Level_2 * random.random() print("Initial condition is:") print(T) # # Setting up of the initial condition is over ################################################################# # Create the x-axis vector, for plotting purposes xAxis = numpy.linspace(0, domLen, nx) # Total time elapsed---physical time tt = 0.0 plt.title("1D First-Order Wave Equation using\nFirst-Order Upwind Scheme. CFL: %5.2lf. Time: %5.2lf s" % (CFL, tt)) plt.xlabel("Distance, in m") plt.ylabel("Transverse Displacement, in m") plt.xlim(0.0, domLen) plt.ylim(T_Level_1, T_Level_2 * 1.5) plt.plot(xAxis, T, 'r') ################################################################# # Time-stepping, using upwind scheme # for n in range(nt): Tn = T.copy() # The trick on the following line is specific to Python's negative indexing for i in range(0, nx): # The following commented line would be the usual way of writing the loop # for i in range(1,nx): # The following step is where the first-order upwind scheme computation is done # Modify this one line to implement other schemes like Lax-Wendroff etc. T[i] = Tn[i] - CFL * (Tn[i] - Tn[i - 1]) tt = dt * (n + 1) plt.clf() plt.title("1D First-Order Wave Equation using\nFirst-Order Upwind Scheme. CFL: %5.2lf. Time: %5.2lf s" % (CFL, tt)) plt.xlabel("Distance, in m") plt.ylabel("Transverse Displacement, in m") plt.xlim(0.0, domLen) plt.ylim(T_Level_1, T_Level_2 * 1.5) plt.plot(xAxis, T, 'k') plt.pause(0.1) plt.plot(xAxis, T); plt.show() print("done")

Some things to try:

*Straight-forward:*As mentioned in the code, change the CFL to various values all of which are strictly less than 1.0. Then, change CFL to a few values greater than 1.0. Draw conclusions about stability.*Not so straight-forward:*Then, also set the CFL exactly to 1.0, and see what happens. Did you expect the outcome? Why or why not?*Programming-related:*Implement the Lax-Wendroff/MacCormack scheme.- The smoothing that occurs at CFL values < 1.0 is due to a phenomenon called false diffusion. Explain the similarity and differences of the false diffusion from the true diffusion (say as modeled in the steps 1 and 2 of this series).

**A Song I Like:**

(Marathi) “dis chaar jhaale mann…”

Music: Ashok Patki

Lyrics: Soumitra

Singer: Sadhana Sargam

[Minor editing may happen…]

[E&OE]