The chaotic Lorenzian waterwheel

  • Difficulty: 018/100
  • Author: jeremy theler
  • Keywords: PHASE_SPACE, CONST, NUMBER, VECTOR, SIZE, .=, IF, ENDIF, PRINT, TEXT, lorenz,

I first read about the Lorenzian waterwheel in James Gleick’s “Chaos: Making a New Science” back in 2004 when I was an undergraduate student. Some weeks after that, I built a physical version of the waterwheel at the University, measured its behavior, compared it to numerical results and prepared a poster about it (presented of course at the classical poster session of Física Experimental I at “Instituto Balseiro”: “http://www.ib.edu.ar”):

FIGURE: [018-waterwheel/jeremywheel] Me, myself and my experimental waterwheel in 2004 label{jeremywheel}

FIGURE: [018-waterwheel/poster] The resulting “poster”: “http://www.talador.com.ar/jeremy/papers/2004-rueda.pdf”

The chaotic waterwheel has somehow fascinated me ever since, and I certainly returned to its study in 2006/2007 during the writing of my “engineering thesis”: “http://www.talador.com.ar/jeremy/papers/2007-theler-grado-1sided.pdf” where in one of the chapters I apply a fuzzy logic controller to force the (numerical) wheel to rotate always in the same direction. Indeed, non-linear dynamics of nuclear reactors was the subject of my “masters thesis”: “http://www.talador.com.ar/jeremy/papers/2008-theler-maestria-1sided.pdf”, which was chosen after my delightful experience with chaotic systems.

This report shows how “wasora”: “http://www.talador.com.ar/jeremy/wasora” can help us building a numerical model of the Lorenzian chaotic waterwheel by just writing the differential equations of movement into an input file just as nearly we would write them in a sheet of paper, and let the DAE-solving libraries to take care of all the numerical details. In particular, this example focuses on a numerical waterwheel model that tries to resemble the experimental version I built almost ten years ago, which in turn is based on the description Gleick gives in his book. Other variations of waterwheels—both numerical and experimental—can be found online.

The Lorenzian waterwheel consists of an array of \(N\) cylindrical buckets of height \(h\) and diameter \(b\) equally spaced around a disc that is able to rotate around an axis, configuring a circle or radius \(R\). The particularity is that the buckets have a little hole at their bottom so water is drained at a certain rate. Water is poured from above at a volumetric rate \(Q\) by means of a faucet or a hose that is located exactly at the center of the wheel. When the buckets are empty, the wheel has a momentum of inertia \(I_0\) and the bearing around which it rotates provides a frictional damping.

FIGURE: [018-waterwheel/wheel3d] A 3D model of the experimental wheel

If the water flow \(Q\) is small, the buckets quickly empty and due to the damping in the bearing, the wheel eventually stops rotating. However, if the flow is high enough, due to the nonlinear nature of the system, the wheel may rotate in either direction and present what it is known as chaos, i.e. non-periodic long-term behavior and sensitivity to initial-conditions. Several studies show that this system can be related to the seminal Lorenz’ system (see the wasora example in section ref{005-lorenz}), whereas in this case the water flow \(Q\) acts as the control parameter \(r\) of the original system.

MOVIE: [realwheel.webm, width=320 height=240] The real chaotic waterwheel I built in action

The equations that describe the wheel’s movement are easily derived as follows. Of course, the angular acceleration \(\dot{\omega}\) is the ratio of the net instantaneous torque \(\tau\) and the instantaneous momentum of inertia \(I\):

!bt [ = ] !et

The angular position \(\theta\) of the wheel is defined as the location of bucket number one, defining \(\theta = 0\) as the vertical position above the wheel bearing and taking the positive direction as the clockwise direction. Let \(v_i\) be the instantaneous of water contained in bucket \(i\). Therefore

!bt \begin{align*} I &= I_0 + \rho R^2 \cdot \sum_{i=1}^{N} v_i \\ \tau &= -k \cdot \omega + \rho g R \cdot \sum_{i=1}^{N} v_i \cdot \sin \left( \theta + (i-1) \cdot \frac{2\pi}{N} \right) \end{align*}

!et

where of course \(\rho\) is the water density and \(g\) is the gravity acceleration. It should be noted that the equation above for the torque does not include the change of momentum between the water falling down the source and the bucket, which depending on the experimental configuration may not be negligible.

Now we have to give the rates at which the buckets are filled and emptied, which is in fact where fun part is located. On the one hand, the \(N\) buckets continually drain water at at rate that in principle is proportional to the square root of the water level through a constant \(p\). On the other hand, a certain bucket \(i\) is filled with water at the rate \(Q\) only when its diameter \(b\) is located beneath the flow up to the arc subtended by the head of the bucket. Using the angular coordinate \(\theta\), this condition can be written as \(\theta + (i-1)2\pi/N\) modulo \(2\pi\) is either less that \(\gamma = \arctan(\frac{b}{2R})\) or greater than \(2\pi - \gamma\) (make a drawing to get this approximation straight). However, as the buckets have a limited capacity, they cannot hold more than a certain maximum volume \(v_\text{max} = \pi b^2 h/4\). The rate of change of the volume contained in the \(N\) buckets is therefore

!bt [ i = -p v_i>0 + Q v_i < v ] !et

Summing up all the equations, we can arrive at a system of algebraic-differential equations of the form \(\mathbf{F}(\mathbf{\dot{x}}, \mathbf{x}, t) = 0\) of size \(4+N\) over the phase-space coordinates \(\theta\), \(\omega\), \(I\), \(\tau\) and \(v_i\) for \(i=1,\dots, N\) that can be easily handled by wasora:

!bt \begin{align*} \dot{\theta} &= \omega \\ \dot{\omega} &= \frac{\tau}{I} \\ I &= I_0 + \rho R^2 \cdot \sum_{i=1}^{N} v_i \\ \tau &= -k \cdot \omega + \rho g R \cdot \sum_{i=1}^{N} v_i \cdot \sin \left( \theta + (i-1) \cdot \frac{2\pi}{N} \right) \\ \dot{v}_i &= -p \sqrt{\frac{4 v_i}{\pi b^2}} \text{ if } v_i>0 + Q \text{ if } v_i < v_\text{max} \wedge \left[ \theta + (i-1) \frac{2\pi}{N} \bmod 2\pi < \gamma \vee \theta + (i-1) \frac{2\pi}{N} \bmod 2\pi > 2\pi-\gamma \right] \end{align*}

!et

wheel.was

To solve model the chaotic waterwheel with wasora the following input can be used, which is—at least from my humble point of view—pretty simple taking into account that the equations are not trivial. Indeed, as far as I know, this kind of input file—i.e. the one that wasora parses and understands—is the one that resembles the most what a human expert would write down in a sheet of paper when trying to model the wheel. That is to say, if \(\mbox{v}\) is known to be a vector of size \(N\), then writing v_dot(i) = something should be enough to state all the equations for \(i=1,\dots,N\). Any other DAE-solver I ever saw needs an explicit loop to sweep \(i\) between \(1\) and \(N\).

In this particular case, the volumetric flow \(Q\) is fixed to a constant value. As I said above, I once designed a linguistic fuzzy controller acting over \(Q\) to force the wheel to always rotate in the clockwise direction. Some time in the future I will re-write the fuzzy control code as a wasora plugin and illustrate its usage within this Real Book.

# lorenzian chaotic waterwheel
# DAE equations for wasora
#   with constant water flow (i.e. uncontrolled)
#

# main parameters (in SI units)
CONST p I0 b h rho R g k
N = 8      # number of buckets

end_time = 40   # [s]


# geometric parameters
R = 0.2         # distance from the buckets to the center [m]
b = 0.0435      # bucket diameter [m]
h = 0.086       # bucket height  [m]

# physical parameters
I0 = 0.04       # momentum of inertia of the bare wheel [kg m^3]
k = 0.025       # bearing damping [kg m^2 / s]
p = 55e-6       # bucket permeability [kg / (s * m^(1/2))]
rho = 1000      # water density [kg / m^3]
g = 9.8         # gravity [m / s^2]

# secondary parameters
CONST gamma v_max step

gamma = atan((b/2)/R)  # semiarc subtended by the bucket
v_max = pi*(b/2)^2*h   # total volume of each bucket
step = 2*pi/N          # angular step

# instantaneous volumes of the N buckets
VECTOR v SIZE N
PHASE_SPACE I tau theta omega v Q

# initial conditions
theta_0 = 1/10
I_0 = I0
tau_0 = -k*omega
Q_0 = 0.0002

# DAE equations
# angular position and acceleration
theta_dot .= omega
omega_dot .= tau/I

# instantaneous moment of inertia
I .= I0 + rho*R^2 * vecsum(v)

# instantaneous torque
tau .= -k*omega + rho*g*R * sum(v(i)*sin(theta + (i-1)*step), i, 1, N)
 
# instantaneous volumes of the buckets
v_dot(i) .= {  if(v(i)>0, -p*sqrt((4*v(i))/(pi*b^2))) 
             + if(v(i)<v_max &
                    ((mod(theta+(i-1)*step, 2*pi)<gamma) |
                     (mod(theta+(i-1)*step, 2*pi)>(2*pi-gamma))), Q) }

# volumetric flow
Q .= Q_0

# give some information on the output
IF in_static
 PRINT TEXT "\# R =" %g R
 PRINT TEXT "\# b =" %g b
 PRINT TEXT "\# h =" %g h
ENDIF

# write the output!
PRINT t omega theta v 
$ wasora wheel.was | qdp -o wheel-dat --xlabel "\$t\$ [s]" --ylabel "\$\\omega\$ [rad/s]" --ycol 2 --ti "\$\\omega\$"
$ 
wheel-dat.png
wheel-dat.png

The output this input generates is directly fed to the script qdp (Quick & Dirty Plot) which generates a plot of the angular speed \(\omega\) vs. time. Is resemblance with the Lorenz attractor discussed in section ref{005-lorenz} is remarkable.

A full analysis of the behavior of this system with respect to the control parameter \(Q\) can be done with wasora using the PARAMETRIC keyword and analyzing the output with different algorithms. However, this kind of study is beyond the scope of this “wasora Real Book”: “http://www.talador.com.ar/jeremy/wasora/realbook”.

online.was

So far with the equations and the plots. The main objective of the waterwheel is fun, and there cannot be fun without drawings. So here enters the wasora plugin besssugo, which provides means to draw graphical representations based on wasora’s objects. First of all, we have to make sure besssugo is correctly compiled and installed (i.e. the shared object file besssugo.so is located in a directory where the operating system can find it such as /usr/local/lib or $LD_LIBRARY_PATH):

!bc console $ wasora -p besssugo.so besssugo 0.2.1 trunk (2014-01-03 19:48:41 -0300 dirty) scientific video generator for wasora $ !ec

In this section we use besssugo at the same time we ask wasora to solve the equations, so we draw the wheel and the buckets at each step. First of all, we write a file called drawing.was (which we will reuse afterwads):

@@@CODE 018-waterwheel/drawing.was envir=wasora

Then we use the following file online.was as the main wasora input file:

# lorenzian chaotic waterwheel
# online solution and graphical representation

LOAD_PLUGIN besssugo
realtime_scale = 1

INCLUDE wheel.was
INCLUDE drawing.was

BESSSUGO_DRAW
$ 

The result cannot be shown here in the report, as it consists—besides the same text output as in the last example—on the appearance of a window containing a visual representation of the wheel and the buckets. Wait for a few more inputs until we generate and show a video of the numerical wheel. In the meantime, if you tried to run the case, you would notice that there is something weird about the time scale. In effect, each time on bucket is located below the flow, the time tends to slow down. This can be explained by noting that even that we asked wasora to try to run in realtime by setting realtime_scale to 1, the DAE solver adapts the time step according to the instantaneous stiffness of the system. Whenever a bucket is filled, the equations become stiffer, the time step decrements and wasora cannot advance in real time. Therefore the time scale shown in the graphical representation is distorted. Moreover, if we wanted to dump the frames into individual images in order to build a video, we would then have to tell the video encoder how the frames-per-second should be adapted as the time step of the DAE solution changed.

offline.was

This time, instead of having wasora to actually solve the equations whilst besssugo draws a graphical representation of the wheel, we first ask wasora to solve the system, dump the solution into a file and then load it with wasora+besssugo and perform a constant time-step graphical representation of the loaded data.

Instead of including the input file wheel.was we use the following input called readmap.was to read a text file provided in the commandline (which should be the ouput of wheel.was), interpolate it and leave it ready (i.e. map it) for drawing.was to understand it:

@@@CODE 018-waterwheel/readmap.was envir=wasora

Thefore, the new case use the same drawing keywords as the previous one but uses a constant time step (which is wasora’s default equal to \(1/16\)).

# lorenzian chaotic waterwheel
# offline solution and graphical representation
# reading the data from the file given in commandline

LOAD_PLUGIN besssugo
realtime_scale = 1

INCLUDE readmap.was
INCLUDE drawing.was

BESSSUGO_DRAW
$ wasora wheel.was > wheel.dat
$ #wasora offline.was wheel.dat
$ 

If you tried to run the case, you would now see a realtime representation of the chaotic waterwheel which should resemble the real model described above.

video.was

Besides generating a window with a graphical representation, besssugo can also dump the individual frames as bitmap files which can then be merged together to generate a standalone video, which can be used in presentations, posted online in webpages and used to illustrate scientific or engineering points with colleagues that do not want to install wasora and besssugo (who wouldn’t, by the way?).

In order to do that, we now use the following input file. This time, it is not needed to set the realtime_scale variable as the transient computation should be done as fast as possible, because the actual frame-per-seconds of the video is set by the encoding tool and not by wasora. It should be noted that besssugo relies on the SDL graphical library which provides routines to dump images in BMP format, so the temporary storage needed to run this case may be quite high.

# lorenzian chaotic waterwheel
# offline solution and graphical representation
# reading the data from the file given in commandline
INCLUDE readmap.was
INCLUDE drawing.was

FILE waterwheel waterwheel-%04g.png step_transient MODE w
BESSSUGO_DRAW
BESSSUGO_DUMP FILE waterwheel

IF done
  PRINT "generating video..."
# 16 FPS is because the default wasora dt = 1/16
  SHELL "./genvideo.sh waterwheel 16"
  SHELL "rm waterwheel-*.png"
ENDIF
$ cat genvideo.sh
#!/bin/sh
# generate a webm video out of png frames
if [ -z "$1" ] || [ -z "$2" ]; then
  echo usage: $0 basename fps
  echo finds files named basename-0000.png and builds a video named basename.avi
  echo with the prescribed number of frames per seconds
  exit
fi

avconv -y -f image2 -framerate $2 -i $1-%04d.png -vcodec libvpx $1.webm -v quiet
$ besssugo video.was wheel.dat
generating video... 
$ 

Finally, here is the video of the proposed transient and the again the experimental wheel video to compare them:

MOVIE: [waterwheel.webm, width=332 height=332] The numerical chaotic waterwheel solved by wasora and drawn by besssugo

MOVIE: [realwheel.webm, width=320 height=240] The real chaotic waterwheel I built in action