I'm a brazilian Mechanical Engineer and PhD, currently a Postdoctoral Scholar at Florida State University, working in various projects in aerodynamics, aeroacoustics and flow control. In the past I've worked as a refrigeration systems designer and later as an R&D specialist at a refrigeration contracting company, researching for new products to push the industrial refrigeration market technologies forward.
Currently I'm working with experimental aerodynamics at FSU, using state-of-the-art equipment and techniques to observe the most intricate flow physics!

Just wanted to put this out as this work in progress evolves! This is something really cool that I hope will make some interesting impacts in the active flow control community.

In 2020, I’ve built and tested a solenoid array driver to individually toggle individually 108 solenoids. This work was published here in the AIAA Scitech 2021 conference. You can access my personal version of the paper on ResearchGate, if you don’t have an AIAA subscription. In summary, we utilized this individually addressable solenoid array to tackle the active flow control actuator placement problem, which turns out is rather difficult because flows display very non-linear behavior. We found an unintuitive pattern of actuators that improved the cost function (circulation of the vortex in the wake of a cylinder, which is proportional to drag); something we would never have guessed.

So now we are going to attempt to apply this technique to an even more difficult flow control problem: Supersonic jet noise suppression. The principle is the same, we have several channels that we can address individually, and we will use an algorithm that measures the overall sound pressure level (OASPL) of a supersonic jet and tries to minimize it by changing the actuator locations and parameters. This is how the actuator array looks like:

This is a work-in-progress post. As things progress, I’ll update this. I’m excited to see the results!

I’m typing this out of a frustration that I’m sure many of my colleagues have already felt when assembling a shadowgraph or a Schlieren flow visualization system: Adjusting the LED focusing lenses! Let’s look at the optics of this problem and try to see the rather unintuitive trade-offs we are making when adjusting these lenses in practice.

The objective function

There are two objective functions that are usually desirable when focusing an LED system: Firstly, we want to collect as much light as possible. Second, we want to focus the LED to a point as small as possible, because we want to pass this light through a pinhole. In fact, the two objectives are the same: Focusing to a pinhole without cutting much light is the same as maximizing the amount of light (first objective). Thus, there’s really only one objective function. Let’s look at the standard collimator arrangement utilized in a Z-type shadowgraph:

Light losses

There are three locations where light can be lost that we need to consider in this type of setup:

Light that is not collected by the lens

Light that is blocked by the pinhole

Light that is not collected by the parabolic mirror

If you are not thinking of these issues while setting up your shadowgraph, then you will lose a lot of light – perhaps more than you can afford to. The result, then, is an alignment struggle!

The reason this problem is rather unintuitive is because we usually don’t think about (3) because we’re focused into getting as much light into the pinhole. But the light lost in the parabolic mirror (3) really is the key to designing the optimal shadowgraph LED lens!

the ideal shadowgraph collector – single lens design

Let’s look at this problem backwards: Start with a parallel beam coming into the parabolic mirror. This beam will focus at where the pinhole should be, which only depends on the parabolic mirror focal length. The convergence angle of the beams, however, depends on the f-number of the parabolic mirror. As we’re going to see, the f-number of the mirror is crucial here.

We seek to put a lens after the focal point such that it focuses at the pinhole location with a beam convergence exactly equal to the f-number of the mirror. If we do so, there will be no loss of light at the parabolic mirror, meaning we literally eliminated a whole source of loss of light!

Let’s do some basic optics math to understand what are the parameters of this problem. We start assuming the f-number of the mirror is fixed, since the mirror is usually the most expensive – and less interchangeable – part of the shadowgraph. Then, we just need to select a lens diameter and focal length to focus the LED positioned at to a point at from the lens. The LED will be considered a circular source of diameter :

Using the thin lens equation, we have the following relationship between the LED distance and the image distance :

I’ll isolate the LED distance for now, as it is an important variable when aligning the lens assembly:

We know that the beam downstream of the lens has to have a defined f-number :

Thus, we get an equation for the LED distance as a function of known variables:

Well, we technically only know . In fact, and are the lens parameters we’re after. Now let’s look at the illumination by examining the light cone half-angle that the lens collects:

The cone half-angle is related to the solid angle by the following relation:

Well, most LED emitters (chip-on-board) are plane emitters. I am not considering, for simplicity purposes, the ones that have a built-in lens. Those usually have poor optical quality for shadowgraph applications and should be avoided, anyways.

Since the solid angle of a plane emitter is , then we get that the intensity collected by a LED lens is:

This is what we want to maximize, which means we want to maximize . One might think the obvious way to increase is to increase the lens diameter. But it’s not that simple! Let’s rearrange the first boxed equation to see why:

Let’s multiply and divide the right-hand side by :

We’ll call , the f-number of the lens.

Divide both sides by again, now:

Which basically means that for a given lens f-number, the LED distance from the lens scales linearly by the constant in the right-hand side. Thus, increasing the lens diameter without changing its f-number will simply move the lens farther from the LED, cancelling out the effect of light collection.

Thus, the lens f-number is really the only relevant variable in determining the light collection efficiency – assuming the lens is reasonably larger than the LED. Let’s plot the trend!

The plot above is for the following parameters: mm, mm, . Note the efficiency is inversely proportional to the f-number of the lens . Also, the spot size at the focal point is also smallest at the lowest f-number. This is great, because both objectives converge. A low f-number lens will accomplish the highest collection efficiency and the smallest spot size, which is precisely what we want for a shadowgraph.

Below, I’m plotting the same trend for a mirror f-number of 10. As we can see, the collection efficiency is weakly dependent on the mirror f-number; which is very likely to be between 5 and 10 for a typical shadowgraph. The spot size is even smaller for the higher f-number mirror, which is even more beneficial.

Unfortunately, the collection efficiency for a single-lens system is really limited to low single-digit percent because lenses with f-numbers of 1 or less are pretty hard to manufacture. Especially achromatic lenses with low spherical aberration, which are required for shadowgraphs with white LED light sources. But the lesson to be learned here is very simple: Just get your lowest f-number lens and make sure the beam f-number is close to the mirror f-number. Higher collection efficiencies might be possible with several lenses, but I haven’t really looked carefully into the math. Coming up, perhaps! In any case, hope this helps a fellow researcher!

Have you ever wanted to make a nice timelapse where you have some camera motion? If you look in the market, there’s some (rather expensive) motorized sliders that you can buy. These usually go for $300-$600, but they don’t have much travel. For example, this model has a 31.5″ (0.8m) travel and costs $258. This one is $369 and has a travel of 47″ (1.2m). I’m sure they’ve got much better construction and much nicer features than what I’ll present here.

But truth of the matter is – we are makers and we like to make our own equipment! So I came up with this while discussing with my friend how to make a timelapse where the camera circles around a subject. Since the radius is usually large, we would need ludicrously large equipment, and obviously it wouldn’t fit in our car. So the obvious alternative is to “loosen up” the requirement of rigidity of the setup and let the camera shake as it moves. Because at the end of the day, we’re in the computing era and we can fix the footage later, right? Well, see the video below for yourself! (Also in the end there’s a little compilation of a few timelapses I was able to take with this rig)

Also, if you look at the BOM; I’ve spent about $50 (plus the tripod) for what is effectively infinite travel – as long as there’s no obstacles in the way and the terrain is smooth enough. A much better deal in my opinion! I’ve tested this design only on pavement and stone sidewalks, and it works “okay”. It is wonderful in cement pavement. But it will fail in gravel, grass or anything too rough – the wheels I designed are simply too small. But perhaps this can prompt you to make a similar design for your own application!

So you can watch the video below or just follow the instructable images. The system is pretty straightforward. Enjoy!

~0.2kg of ABS filament [$4] – DON’T USE PLA if you intend to leave this in your car!

[Optional] 3.5mm female jack (Link) [$4 ea] for shutter release

[Optional] Shutter release cable (Link) [$7 ea] – MAKE SURE YOU BUY THE CORRECT CABLE FOR YOUR CAMERA!

So total estimated cost, without random materials like PCB, wires, solder, etc is about $47 (without the optionals). Might not be the cheapest project, but considering there’s not really any options to do this on the cheap (to my knowledge), this is quite decent!

Print the Parts!

Alright – So the first step is to print the design. I published the design in Thingiverse, both STL files and STEP files if you perhaps want to modify the design for your own necessity. The design I came up with is a clamp mechanism that simply adds wheels to each leg of the tripod. The front “driver” wheel has a stepper motor and the board that controls everything.

he driven wheels are actually quite straightforward to assemble. Screw the clamps to the main part and hammer the 5mm pin for the wheel in place. The wheel is bolted together for increased contact surface. Put the bearings on it and stick it to the shaft! If you want, you can use an end cap to prevent the wheel from falling off. My bearings were well-stuck to the shaft, so I skipped that part, but the caps are in the Thingiverse design.

Make the Driver Assembly

The top part of the driver assembly is the same as the other wheels. The “foot” mounts on the swivel plate, which has two arc-shaped slots for the driver assembly. The first two pictures show the first iteration of the driver, but I ended up revising it (last picture) to improve stability. The second version is the one published on Thingiverse.

Tripod Is Fully Converted! Now to the Brain:

So I’m using an Arduino Nano for this project. The first iteration used the paintbrush schematic in the latter pictures, but I finally found some time to make a PCB design of it. Link to the EasyEDA project below:

As you can see, we just connect the components together. The arduino serves 2 purposes: Step the stepper and enable an interface with the Bluetooth. The interface will be done with an app, which saves us the trouble of building the buttons/etc that goes into making an interface. It also gives us a lot of freedom to change what the system does!

I added in a revision a feature for connecting a shutter release cable to your camera, which usually converts your camera connector into a microphone 3.5mm jack. With this cable, the Arduino will be able to trigger the camera. This can be used for night timelapses, so we are not driving the tripod while the camera is exposing. It’s an optional feature, though.

Download the App to Configure the Timelapse Parameters

So I built a little app with MIT app inventor (see file attached, I had to “hack” the upload due to file extension restrictions in the Instructables website. Change the extension from txt to zip). You can upload the *.aia file to MIT app inventor to see how the app is built now, and you can install the *.apk file on your phone if you want to use it [the interface might be a little jagged in your phone because I’m not a professional app developer!]. If you’re afraid of installing outside apps on your phone (which you should be), then you can compile the MIT app inventor project and get a QR code, which you can use to install the app on your phone instead.

The app basically builds a string of 4 bytes that is sent to the Bluetooth HC05 device through the Serial interface. This 4-byte string has encoded in it the 5 parameters of the system, as follows:

VVVVVVVV VVVVVVVV MDDDDDDD SWWWWWWW

Where VVVVVVVV VVVVVVVV is a 16 bit integer that encodes the 10*velocity inputted by the user in the app, M is a bit encoding whether the motor is on or off, DDDDDDD is a 7 bit integer with a time (in seconds) for driving the system in stop-and-go mode, S is the bit that turns on the stop-and-go mode, and WWWWWWW is the 7 bit integer with the wait time of the stop-and-go mode.

I encoded it like that because the MIT app inventor has a built-in function that sends this 4-byte string. This means I can send one command only, without having to handle a pipeline or anything like that. Then the Arduino just needs to decode this and set the parameters for the stepper, as we’ll see in the next step!

Program the Arduino

Well, there’s not much here. The firmware is attached. It’s basically the byte manipulation I discussed in the last step! I also implemented the shutter release. This should be able to trigger the camera for beautiful long exposure photos and eliminates the need for an intervalometer (as the Arduino itself can serve as an intervalometer).

Post-processing for Stabilization of the Footage

As discussed, we can’t really do much mechanically to prevent shaking. Perhaps larger wheels with a softer material will shake a little less, but they become unwieldy when storing the parts and I don’t think they would stabilize the footage to the point of no need for post-processing. If you do build a variant with larger wheels, I’d like to hear from you if you could skip the post-processing step!

I used Adobe Lightroom to convert the pictures and then Adobe Premiere to post-process the footage. Warp Stabilizer does the job incredibly well – I was stunned! If you use a large value for smoothness (like 100% or more), the stabilization is really decent. Since the warp stabilizer filter in Premiere is quite simple to use, I believe I can’t really add much here. Just apply it to your footage and enjoy the result!

Thanks for reading this – please let me know if you found this useful and if you build one for yourself! I would love to see footage you generated with this idea!!

(Optional) Adapt for a 3D Printer

Another thing that this timelapse system can do is a camera circling around a 3D printer. You need to use a connection to the 3D printer, in this case. Just connect the printer trigger signal to the D6 pin in the Arduino and use the code attached instead. Now, when the printer is done with the layer and parked somewhere, you can usually add a custom M42 G code for toggling a pin in the printer board that will acknowledge the Arduino to move the camera and take a picture.

For those acquainted with POD (Proper Orthogonal Decomposition), it’s very easy to throw POD at everything data analysis related. It’s such a powerful technique! For those who aren’t acquainted, I really recommend. I’ll introduce with a refresher. If you’re good, follow to the next section! We can decompose a set of N dimensional data points, which we will call here , into three matrices. Since I’m not too fond of the “fancy mathematical jargon”, let’s just start with a practical example. I have a movie of particle image velocimetry (PIV), which is a data set that contains the information of a flow in a plane where I was looking at. This movie is a vector field, meaning it is a 2D set of data (like a screen), but each “pixel” has three numbers to describe the information contained (Velocity measured, in , and direction). See the movie below and think this could very well be simulation data or any other data set!

Therefore, the movie is a 4D data set: Let’s count the dimensions: (1) Pixels in the X direction, (2) pixels in the Y direction, (3) vector component (, or ) and (4) time. Since POD can only handle normal matrices, we have to do some remodeling of the data set to fit matrix algebra, which can only handle 2D arrays of numbers. In order to do that, we need to flatten all dimensions but one. Which one should be left unflattened is up to you, but in “movies”, it makes a lot of sense to leave the time dimension unflattened, since we know that things are changing through time. Technically, you could do the POD in the spatial dimension instead, and perhaps get some insight from it, too, but for this example let’s stick to the more intuitive case where dimensions (1, 2, 3) are flattened into a column vector and dimension (4), the time dimension, is left unflattened.

If we do that, we can construct a data matrix , which has dimensionality . , the number of rows of , is the number of entries necessary to build the flattened snapshot of the movie (in our case, dimensions 1, 2, 3 are flattened in the direction, so is the number of elements in each of those dimensions multiplied together). is the number of time snapshots in the movie. This matrix can be decomposed into three matrices, which – without getting too much in the detail – are related to solving the eigenvalue problem of . Since is , we can’t technically solve an eigenvalue problem, but instead we can find the singular values (and corresponding vectors) of this matrix – hence POD is also called Singular Value Decomposition. These are as important as eigenvectors, in the sense that they form a set of basis vectors, othonormal in nature, that optimally span the space . In lay words: The singular vectors of are a set of nice pictures that represent the most variance in , and the singular values are the square root of this variance – meaning if has most of its content changing in a certain way, POD will capture that in its first few singular vectors/values.

POD is usually represented as this decomposition:

Where is the “mode matrix”, and each of its columns contains the flattened version of the mode – which can be unpacked by reversing the process of flattening used to generate the columns of . is a diagonal matrix with the singular values, and the modern SVD algorithms (like svd() in Matlab) return in descending order (i.e., the first diagonal entry is the most important singular value). What we’re actually after here – for the current discussion, at least – is . This is the “right singular vector matrix”, and each of its rows contains a series of coefficients that, when multiplied by the corresponding mode column, gives the “movie” of only that mode playing. It, in a sense, is a “time series matrix” when the unflattened dimension is time, and its coefficients can be analyzed to interpret what each mode is doing to the system under analysis. For example, you can try to take the Fourier Transform of a row of to understand if there’s any periodic behavior related to the corresponding mode in . It’s quite useful!

Some caveats and my two cents

Performing SVD in low-rank systems (i.e., when most of the energy in the eigenvalues of is contained within the first few entries) can be quite insightful, but sometimes we need long time series to get converged statistics (say, for example in the case of Fourier Analysis of , we need as many entries in as possible). This is particularly the case in experiments, where noise can contaminate the spectrum and spectral averaging can “clean up” our data by a good amount. In general, Fourier analysis needs long data sets to capture periodicity with a good amount of certainty. Therefore, we would like to perform SVD in the whole matrix so we can get the whole matrix, right? Well, yes, but it is not uncommon to run out of memory! For example, doing SVD on one of our data sets where was flooded the 128GB memory of our data processing computer! SVD is, unfortunately, an algorithm that needs quite a bit of memory – more specifically, we need memory, where is the shorter dimension. But if we assume (really, not a bad assumption in many cases) that the system is low-rank (say, only the first ~300 modes really matter), we could potentially perform SVD in the matrix, which is quite doable even in a laptop because the memory complexity scales with . The problem, is obviously, that now our matrix has only entries, which is not enough to do good Fourier analysis on. So what can we do to solve this?

Well, let’s remember that each column of is a basis vector to describe the data set. Let’s say we compute a cheaper mode matrix like the big matrix – let’s call it (where is the number of reduced dimensions used, in our example , whereas ). We can project every snapshot (column) of back into to get its entry in . To see the math better, let’s show how it works in matrix algebra:

Each column entry of is a flattened snapshot. Let’s look at how we generate one snapshot from the SVD:

So we multiply the full mode matrix and singular value matrix by a column vector that corresponds to the contribution of each mode, . Note the dimensions all cancel out, so we can “stick” our reduced matrix and its corresponding singular vectors in the mix, right? Let’s do that:

Where are the time coefficients of the particular snapshot we’re looking at (and in a sense, is what we’re after). Let’s do some math now – premultiply both sides by :

We get an identity matrix because is orthonormal. So let’s pre-multiply everything by the inverse of once more to get rid of it:

Or, more neatly:

This is actually quite neat. A column vector of the reduced rank right singular vector matrix is obtainable from ANY corresponding snapshot by projecting it onto the subspace of . This means we can obtain cheaply by simply randomly choosing a fraction of our data set to compose , and then perform this projection to the remainder of the data set in order to get all the columns of , for the first modes at least. This would be a snapshot-by-snapshot process, without the memory overhead of SVD. Of course, you can use some more “mature” algorithms, like randomized SVD, but for me the issue is that you need to load the data matrix twice: Once to produce the randomized range finder, and another time to compute the randomized SVD. When you can’t even fit the data set in memory and the data transfer time is really long, you don’t want to duplicate the loading process of the data matrix . This approach, although approximate, removes that need.

Of course, you don’t get the full benefit of SVD. In the end, is an approximation of , thus the projection will leave some stuff out. also is using less snapshots to be constructed, so you’ll have an error related to the modes themselves with respect to . But when lack of memory literally prevents you from doing SVD altogether and you’re most interested in analyzing (i.e., only the first modes of ), then you can still do the analysis to a pretty decent extent. Also, it might be a compelling way of storing time-resolved, 3D computational results without all the hard drive storage necessary (like a compression). So let’s look at some really nasty, noise contaminated, experimental data sets:

Results

Ok, so we’re looking at the wake of a NACA0020 airfoil that I used in a study I’m working on. The airfoil is at 0º angle of attack, but it still sheds coherent structures because of the shear layer instability. We want to pick up what are the frequencies of these coherent structures with Time-Resolved PIV, but this technique – even though it’s amazing – is a little tough to get to low noise levels. The individual vectors in the snapshots have uncertainties of the order of 10% of the mean flow! But as we will see, this is cake for POD – as long as there’s low-rank behavior happening.

So we construct the matrix by randomly sampling vector fields from the full data set. I used randperm() in Matlab for this part. Since the computer I’m using can handle up to about 3500 snapshots, we will look at 3500 snapshots (the full SVD) versus 350 snapshots (the reduced SVD), and then project the full matrix onto the subspace of to get the time coefficients of the first 350 modes.

Comparison of the time series of the first POD mode of the wing at ; “Full POD” corresponding to 3500 snapshots; “Random POD” corresponding to 350 snapshots and the projection method described here

As shown in the figure here, there’s not really much difference between the time series of the reduced, “random” POD, versus the full-fledged POD with 3500 snapshots. The difference in SVD computation, however, is quite a big deal:

Data Matrix size

Memory Consumption (Matlab SVD ‘econ’)

Computational Time

Full POD

5.7 GB

252 s

Reduced POD

too small to measure (<100MB)

7.1 s

SVD memory complexity for the case studied

As we can see, the difference is quite large in this “toy, but practical case”. Perhaps we could use more snapshots in the reduced POD case, but I also wanted to push it a little. The mode shapes look indistinguishable from each other, but of course we could get an euclidean distance metric (I didn’t really bother). Looking at the spectrum, we can see that the dynamics are pretty much indistinguishable in the first modes:

I’m not fully into the CFD business, but I think this simple technique could be used to compress CFD time series data, since we already know most flows do display low-rank behavior (and we are already doing these decompositions anyways to find this low-rank behavior!). In other words, I think one could: (1) run a DNS/LES solution until their statistics are converged (as is usual in the business), (2) run for some period of time while saving enough random snapshots to generate the snapshot matrix – you still need to run for a long-enough time, though – and then (3) perform POD on the snapshot matrix, and finally (4) run for as long as necessary, and project every snapshot onto this random POD basis to get the matrix (low memory consumption). This enables one to save the matrix, along with the first few modes of onto a file, and in the future reconstruct the entire simulation time series (to reasonable accuracy) for further analysis. Obviously the analysis will be limited to the compression accuracy and the time span of the capture process, so one has to be aware of that!

Well, perhaps I’ll try that in the future. But for now, I’ll leave the idea floating so someone can perhaps find some use to it.

I’m writing this because even though it’s already 2020, this kind of stuff still cannot be found anywhere in the internet!

Stacking layers of materials is something we all do in so many engineering applications! Electronic components, batteries, constructions panels, ovens, refrigerators and so many other layered materials that are inevitably under heat transfer. We want to predict the temperatures on the stack, at all layers, to ensure the materials are going to withstand the transfer of heat. In order to do that, any engineer worth their salt will already have used the Thermal/Electrical Circuit Analogy. And it works wonders!

But what if we want to do transient analyses? What if the materials are going to be stressed by some quick burst of power; or some sequence of bursts at a certain frequency? Of course, we can always just do some FEM and solve our issue – but one might find that this can become a multi-scale problem if there’s low frequencies AND high frequencies involved.

This is commonly the case in semiconductors. Applications I’m familiar with are related to LED pulsing, but I’m sure one will find many others. Another more obscure application I found in my research has to do with unsteady measurements of surface temperature when fluid is flowing on a surface. Some frequencies will be too fast and the stack of materials will not respond, so how do we proceed in that scenario? How to predict the temperature rises and – more importantly – understand what we need to change in the design in order to improve its frequency response?

Well, it turns out we can extend the thermal circuit analogy to transients by including the component’s thermal capacitance. Most layers in a stack would then be easily lumped into RC low-pass filters. This is the “basic approach” and it works – to a point.

So let’s examine it, to understand its pros and also its limitations:

The thermal capacitor

Capacitance in electrical circuits is the ability of a given electrical component to store charge. In the thermal analogy, this “charge” is heat – and the thermal capacitance is just the mass of the component times its heat capacity (). If we look at the energy conservation equation for a lumped element:

Which means that, for a given heat flow , a small capacitance leads to a large rate of temperature rise (large ). I will make this problem more 1-dimensional to consider the case where the stack of layers is much thinner than wide, because that’ll lead to interesting observations when we drop the “lumped element” part of this model. For now, let’s simply divide both left-hand and right-hand sides by the sheet area , to get everything in terms of heat fluxes:

The terms then are what we will call the specific component’s capacitance (), in the SI units of J/m^{2}.K. One can then represent a layer of a stack like the following RC circuit:

Obviously, this single-layer model is rather simple and does not constitute of much of an issue. The RC circuit has a time constant, , where . You can check the units: has units of m^{2}.K.s/J, which times the units of gives us a time constant in seconds. The cut-off frequency is , and we get a -20dB/decade low-pass filtering behavior.

It is really worth observing, though, that the capacitor is not connected to , but to an elusive “ground” temperature. This is because, as long as you have a non-zero , you will have heat transfer to the thermal capacitance. Therefore, the ground level can be considered absolute zero temperature – for didactic purposes. The temperature of this ground level is not important, though: since capacitors always block DC, its absolute value will never play a role.

Stacking them up (plus all the ways you will mess this up)

For thermal stacks, it is not straightforward how to assemble the circuits together. Do we put the RC circuits in series with each other? As it turns out – not really. To understand why, let’s think about what happens in a two-layer stack:

Looking at the picture above, we have a stack of two layers and 3 different circuits that could (potentially) represent it. The circuit (a) represents how we would construct the thermal RC network using the same strategy as for the single layer explained before. (a.1) is simply another representation of the same circuit, perhaps less “bulky”. The circuits (b) and (c) are wrong because they do not represent the thermal equations that describe the stack. Let’s see why:

The first layer of the stack has an input heat flux of , which according to the energy conservation equation applied to this lumped element, can be split into the heat flux that stays in the element () and the heat flux that goes through the element (). This heat flux then goes to the second layer and splits again into and . The energy conservation equations here are equivalent to the Kirchoff current laws applied to the nodes at and in the circuit representation (a):

This sounds easy because you already have the answer spoon-fed to you. But if you were left to your devices, I would bet many of you (including past me) would come up with the circuit representations (b) and (c). (b) doesn’t make any sense because capacitors block DC, therefore there would be no DC heat flux through the stack – total non-sense!

The circuit in (c) also does not make sense, but it is not clear why at a first glance. If you follow Kirchoff’s laws, however, you will see that (c) cannot be right, because in the node at the heat flux equation would have to be – which is simply not the correct energy balance.

It is not really very intuitive that the capacitor of the first layer of the stack should be connected to ground, instead of . In fact, I made quite a few times this mistake. But looking at the energy balance, we can’t really deny there is a correct equivalent circuit.

We can also solve the full heat pde with circuits!

Ok, so we tackled the single wall and the stacked double-wall problem. I think having gone through how this can go wrong, you can now build this to multi-walled systems easily. As you probably already realized, each new layer will behave as a new first-order system, and some very interesting behavior can appear as extra poles are introduced by the new layers. One thing, however, that is worth noting, is that this approach completely fails at high enough frequencies.

To understand why, let’s think about it a little further. If the temperature oscillations in are high frequency (i.e., their period is below the time constant), the lumped element model would predict a decay rate of -20dB/decade, which is what is expected from a standard first-order system.

But the problem in the (zero-dimensional) lumped element model is that the whole element’s temperature is assumed to be a constant value. At high enough frequencies, this won’t be the case anymore. We can see why the lumped element model becomes inaccurate by modeling the standard 1D heat equation over a slab:

The negative sign in the left-hand side is related to the direction of the coordinate system. One can recover the same equations by distributing the lumped element model related to each layer into the following circuit:

Where (thermal resistance per unit length) and (thermal capacitance per unit length). This circuit immediately follows from the circuit presented in (a.1) in the previous figure; since one can consider a single layer can be decomposed into an infinite stack of infinitely thin layers. Not all mechanical engineers will recognize, however, that this corresponds to a (simple) version of the 1-dimensional transmission line model. Robertson and Gross (1958, “An Electrical-Analog Method for Transient Heat-Flow Analysis”, Journal of Research of the National Bureau of Standards) did recognize this in their paper, for example. Funnily enough, their figure shows they also got the RC network wrong, since the first resistor is shown between the first capacitor and T1 – an indicator of how tricky this business can be!

A thermal layer is actually a transmission line

I invite you to read the article on Wikipedia about transmission lines if you’re not familiar with this concept, or read the book “Time-Harmonic Electromagnetic Fields”, by Harrington, 1961. Transmission line is an electrical model for pairs of conductors at high frequencies. At high enough frequencies, the pair of conductor’s distributed capacitance and inductance creates waves inside the conductor, which can be modeled by the wave equation. The element of the electrical transmission line is displayed below (figure from the Wikipedia article):

The electrical transmission line follows the 1-dimensional damped wave equation, which is derived from the Telegrapher’s equations:

We can easily see this is the wave equation by taking the Fourier Transform of the two equations:

And then replacing, say, the derivative of the first equation in the second:

Taking the IFT:

Which, obviously, is the wave equation. It’s damped by the distributed resistance and conductance . If we come back to the thermal stack and consider there’s no conductance and no inductance , we recover the heat equation:

Where the thermal diffusivity is ; as expected.

how this interpretation can be useful

Ok, by now you might be finding this all quite trivial and, perhaps, an over-complication of the heat transfer problem. Why would one need to invoke the far more complex “transmission line” to describe the behavior of such a simple system? Well, it basically allows us to easily determine the temperature transfer functions at different points in the stack (perhaps, easily is not the best word, but with a couple easily-programmable equations), without having to solve the differential equations numerically. Plus, it allows us a very cool insight about the behavior of the themal stack at higher frequencies. As it turns out, at high enough frequencies the first layer’s characteristic impedance – a material property – dominates the transfer function.

Thermal transmission lines in practice: the pulsed led

Let’s first analyze the following case that I find is the most practical for high frequency thermal stacks: An electronic power dissipation device that sources power in short bursts. I mentioned in a previous post that you can massively overdrive LEDs to produce a lot more light than “designed” for short bursts, which enables one to expose high speed cameras – like a camera flash. We’ll look at the Luminus Devices Phlatlight CBT-140 for this analysis. The reason why AC thermal analysis is important in pulsed LEDs is because, even though the mean temperatures of each component might be within spec due to the high thermal capacitance of the heat sink, the power burst will temporarily increase the temperature of the LED die, and we obviously also don’t want that to exceed the specs of the LED. Concerns about thermal expansion will not be addressed in this discussion – but you should also do your research on that!

The CBT-140 die looks like shown below. The named squares in the electrical analogy (bottom right) represent the “transmission lines” for each component.

We will consider the die is made of sapphire (many high power LED dies are made of it). This might be incorrect for this model – but I could not find any specs on that. The die has 0.31mm thickness; whereas the copper has 3.44mm thickness. The aluminium heat sink is pretty large. For the sake of this analysis we’ll consider 10mm. As we’ll detail later – it won’t matter much for the AC analysis.

So what we want to find is the “input impedance” as seen by the current source just above the die in the figure. The current source is a heat flux source, which in this case would be a square wave with some (low) duty cycle, since power is generated by the (micrometer thin) LED just above the sapphire die. Because there’s a cascade of 3 “coupled transmission lines”, we have to use the transmission line input impedance formula three times. From Wikipedia:

Where is the input impedance as seen by the node at and is the reflection coefficient due to the impedance mismatch of the aluminium slab and the convective boundary (also the load impedance) . in this case is the aluminium’s characteristic impedance, which in general has the form:

Note the characteristic impedance is inversely proportional to the square root of the frequency. This will explain the “half-order” decay (-10dB/decade) behavior in the frequency domain we’ll see later (in contrast with a first-order expected from these RC networks). The square root of predicts a phase lag of 45º at high frequencies (not 90º!). Also referring to the Wikipedia article, the propagation constant is a material property:

As we discussed, the math is not the “prettiest”, but this can all easily be hidden behind an Excel or Matlab formula. In this case, using the properties of the three materials shown in the table below and a convection coefficient of , we get:

Aluminium 6061

Copper

Sapphire

180 W/m K

287 W/m K

34.6 W/m K

1256 J/Kg K

376 J/Kg K

648 J/Kg K

2710 Kg/m^{3}

8800 Kg/m^{3}

3930 Kg/m^{3}

In order to find the “input impedance” at the node , we use the as the input load (), and then iterate. Since we’re looking for the maximum temperature in the stack, , let’s calculate the overall input impedance of the stack :

In this case we are driving the LED with a 500 nanosecond square pulse that repeats at 100 kHz. The drive current of the LED is 200A, and given the die area of 14 mm^{2} and a voltage drop of 20V in the LED, we have an instantaneous power of 285 MW/m^{2}. That’s a lot of power!! Well, as it turns out, the mean flux of power (in this case, at a duty cycle of 5%, 14.2 MW/m^{2}) appears in the spectrum at the DC component, but the sharp bursts of power also generate incredibly high amounts of high frequency power flux that really need to be considered.

Looking at the figure above, the solid, thick black line is the actual frequency transfer function of the whole thermal stack. The dashed blue line represents the lumped element assumption and it is a good model until about a frequency of 0.05 Hz. It then starts to fail when the aluminium starts to kill the waves very close to its surface due to its “large characteristc impedance”, . You can think about it like; its large capacitance combined with its (relatively) low conductivity prevents the thermal oscillations on one side of the material to immediately propagate to the opposite side – invalidating the lumped element assumption.

This issue happens once again for the copper (not clear exactly where) and then finally for the die itself (the sapphire) above the frequency of 40 Hz. Then the characteristic impedance of the sapphire dominates the transfer function – i.e., the sapphire is thick enough to be effectively “infinite”, as far as high frequency heat transfer is concerned (the sapphire is 0.31mm thick!!). The temperature waves then decay inside the sapphire, and very little oscillations can be detected at its opposite side. At those frequencies, it doesn’t matter anymore which material you put at the interface. It will not see any temperature oscillations.

All of this insight we get without needing to numerically solve the heat equation is quite fascinating. We now know that at high frequencies, the dominant part of the frequency response is a material property – meaning that a material choice here is crucial for high frequency heat transfer!

How important? Well, let’s run the numbers and see what the wave looks like. I removed the DC part of the wave, because this LED really needs some thermal mass attached to it (the aluminium block) as no amount of convection or heat-sinking will get rid of this insane amount of power. Therefore, let’s only focus on the oscillatory part of the wave:

We can see that we predict an overtemperature of almost 25 K for the die surface for this case. This is on top of the die heating up from the constant heat flux at DC. It is quite a lot and could (in fact, DOES) easily lead to die failure if one only reads the copper temperature. Therefore, one needs to account for this overtemperature when designing an emergency shutoff system – otherwise, it’ll only think about shutting off 25 degrees after the LED failed!

Temperature sensitive paint – high frequency measurements

So this one is quite an interesting application from the research standpoint! We use temperature sensitive paints to perform measurements in aerodynamic models. We paint an aerodynamic surface with this special paint, and it fluoresces under UV illumination, emitting a strong (usually red) light signal at colder temperatures and a weaker signal as the paint surface becomes hotter. An example of supplier of these paints is given here. Though this is quite a specialized application used mostly in high end aerodynamics research, there’s still a lot of interest into understanding its behavior.

So let’s say we want to perform unsteady measurements of heat transfer coefficient at the surface. We put the surface under a flow of nominally known temperature (for this example, a jet impinging on a surface) and then we use the TSP technique to image the temperature map over the impinging surface; getting a full temperature map. As the flow oscillates, we have fluctuations of this convection coefficient – which translate to measured fluctuations in temperature. Obviously, we might not be able to see them if the temperature does not fluctuate much because the substrate of the surface has too much thermal capacitance (as discussed throughout this article). In fact, now we know we have waves of temperature going through the material, which will be rapidly decaying. It might perhaps be the case that the waves are so weak that we cannot detect them. Therefore, we need to assess whether we can measure these waves by comparing the expected signal generated by the oscillations with the dynamic range of our camera.

Let’s start with a rather simplistic setup – to illustrate the issue. Say we have an impinging jet as in the figure below. The impinging jet is the flowing arrows at .

A “thin” layer of temperature sensitive paint is painted on top of a 10mm thick aluminium plate. The layer of paint, although thin, will have an effect on the transfer function. Here’s what happens as we change the paint thickness. The solid black line represents a reasonable paint thickness (5 microns).

As we can see, the transfer function is highly dependent on the thickness of the paint: If the paint is too thin, we get dominance of the aluminium block on the transfer function, which reduces the signal gain very significantly at higher frequencies. As the paint layer thickens, it starts dominating the transfer function – but knowing the paint thickness might be a challenge.

A solution to this issue is to use an insulating film of similar properties to the paint binder (in this case, the binder used for the calculations was a Polyacrylic Acid polymer). Most importantly, if the two materials have similar characteristic impedances , then the transfer function obtained will not be as dependent on the paint thickness.

Polyacrylic acid has the following thermal properties (sources here and here). The film considered here is a PET film (this one).

Polyacrylic acid

PET Film

0.543 W/m K

0.21 W/m K

1248 J/Kg K

1465 J/Kg K

1440 Kg/m^{3}

1170 Kg/m^{3}

As we can see, the thermal impedance of the two materials is somewhat similar. This means less reflection will occur at the interface; which will help with a reduced dependence of the transfer function on the paint coat thickness. The circuit analyzed becomes the following:

This graph shows a great improvement; because the transfer function at higher frequencies is less dependent in the paint thickness. Not perfect – this would only be possible with a perfect match of thermal impedances. Nevertheless, this can definitely be used to make some measurements.

As a remark – the frequencies of interest in this work are around 20 kHz. The gain is then about -38dB (5 microns paint thickness); which is defined as . This corresponds to a signal gain of 0.011 – or; for every 1 degree of temperature fluctuation in , we get one hundreth of a degree of fluctuation in the paint surface temperature; which is quite a small signal. Since the TSP gives a signal of about 1.4%/ºC, we are talking about an intensity signal of , or 0.6 counts/ºC of a 12-bit camera where exactly saturates the camera (most likely not). I expect this paint measurement system to (hardly) detect temperature fluctuations of 10ºC pk-pk; considering good post-processing and noise reduction to reject camera noise and a very tame signal (like a pure frequency).

I’ll post updates as we perform experiments!

Tying it all back to the fourier number

So we see this interesting transfer function that switches between the lumped RC circuit assumption (with a -20dB/decade decay) to the continuous transmission line (with a -10dB/decade decay). Why does that happen and where is the “kink” in the graph?

Well, we’ll see this ties back to the thermal theory; once we find the kink frequency: The kink is where the -20dB/decade and the -10dB/decade cross. At the -20dB/decade regime, a slab of material has behaves as a capacitor and its thermal impedance is:

We have the definition of in the previous sections – so let’s focus on the math. The thermal impedance of the slab in the -10dB/decade regime is simply its characteristic impedance (another interpretation is that in the -10dB/decade regime the slab is effectively infinite).

The kink is located about where :

This gives us the radian frequency of the kink:

Now you should already be recognizing that we have an interesting non-dimensional coming up, as :

The non-dimensional is; obviously; the Fourier number:

We usually see the Fourier number with a time in the numerator , but this is the same kind of transient phenomenon. If the Fourier number is larger than 1, the “time scales” are large (low frequencies) and the heat transfer of the slab can be modeled by a simple capacitor. If the Fourier number is smaller than 1, then the “time scales” are small (high frequencies) and the heat transfer needs the heat equation to be modeled, which we are simply doing by using the transmission line model. The analogy is quite powerful and enables simple, quick analysis of these systems.

I boxed the equation because it is quite useful to determine whether our AC analysis needs to consider the transmission line model; knowing the frequencies of interest. Hope this provides the insight you needed for your analyses!

So in this (last) episode of our quest to visualize sound we have to do some sound-vis experiment; right? Well, I did the experiment with an ultrasonic manipulator I was working on a couple months ago. I built a Z-Type shadowgraph (turns out I didn’t have enough space on the optical table for a knife edge because the camera I used is so enormous) with the following characteristics (see part II for nomenclature):

For this first experiment I used a phased array device with 64 ultrasonic transducers:

When I give all transducers the same phase, they are supposed to produce a plane-like wave. This is how the wave looks like in this shadowgraph setup:

Where the oscillation levels are plus and minus 50 levels from the background level. Thus, we have and the expected pressure fluctuation at 40kHz is 257 Pa, corresponding to a SPL of 142.2dB. Quite loud!

A measurement with a very good microphone (B&K 4939) gave me a pressure fluctuation of 230 Pa – which is quite close to the expected value from the shadowgraph!

I made a little video explaining the details of the setup and there’s some really nice footage of acoustic manipulators I built on it. Perhaps you can find some inspiration on it:

So this post was rather short. Most of the content is in the video, so I hope you have a look!

As discussed in Parts I and II, we established that we can use a Schlieren or a shadowgraph apparatus to visualize sound waves. The shadowgraph is not as interesting an instrument, due to its stronger high-pass behavior. Nevertheless, both instruments are viable as long as one makes it sensitive enough to see the waves.

When reading the previous posts, the keen reader might be weary of my isothermal assumption when applying the ideal gas law. Shouldn’t temperature also be a spatial wave? The consequence, obviously, would be that the waves become sharper with increasing loudness and easier to visualize as the non-linearity kicks in. The reader would be rather right, the non-linearity effect can be significant in some cases. But for reasonable sound visualization, the pressure fluctuation amplitude has to be above about 1 order of magnitude of the ambient pressure in order for this effect to change the wave enough. Let’s look into the math:

For a compressible flow, assuming isentropic sound waves, we have:

Where is the specific heat ratio for the gas under analysis. A pressure wave of the form then becomes a density wave of the form:

Taking the derivative:

Staring at this expression for long enough reveals we can separate the non-linearity term from the remaining terms:

The non-linear function represents the wave steepening as the amplitude becomes sufficiently high. Looking at the phasor in the denominator () we see the derivative operator split the time fluctuations from the spatial component of the fluctuation. For our detectability purposes, what really matters is the minimum and maximum value of to get a sense of a peak-to-peak amplitude. In this case, the phasor minimum amplitude happens at:

Similarly, the maximum happens when the phasor at the denominator minimum point:

Note this value can even be zero, at about 194dB SPL, in which case these equations blow up to infinity. Nevertheless, the difference gives us the peak-to-peak amplitude of the non-linear derivative of density:

Finding the minimum SPL as a function of wavelength is not as simple, since the non-linear function makes it difficult to isolate . However, we can do the opposite – isolate the frequency and plot the points with the axes flipped. To do so, let’s isolate the wavenumber:

The factor of 2 comes from the fact that the calculations performed in the previous section considered the amplitude of the sine wave, not its peak-to-peak value. Now we can replace the equations for the Schlieren sensitivity to obtain:

(Schlieren)

For the shadowgraph, a second derivative needs to be taken. But as shown in the chart below, the non-linearity only becomes relevant above about 180dB. In the case you get to these SPL values, you probably have a lot more to worry on the compressibility side of things. But it comes as good news that we can use a simple relation (like the one derived in Part I) to draw insight about sound vis.

Comparison of minimum sensitivity for the linear (black) relationship derived in Part I and the compressible, non-linear relation derived here. Parameters are the same: m, mm, m and

This continues our saga we started at Part I (spoiler alert: you’re probably better off with a Schlieren). Thanks to a good discussion with my graduate school friend Serdar Seckin, I got curious about applying the same sensitivity criterion to a shadowgraph system. Turns out, Settles’ book also has the equations for contrast in the case of a shadowgraph (equation 6.10). It is even simpler, since the contrast is a purely geometrical relation:

Where is the distance between the “screen” (or camera focal plane, in case of a focused shadowgraph) and the Schlieren object. Since the ray deflection occurs in the two directions perpedicular to the optical axis ( and ), we get a contrast in the direction () and another in the direction (). The actual contrast seen in the image ends up being the sum of the two:

Since the ray deflection is dependent on the Schlieren object (as discussed in Part I), we can plug its definition to get:

Which is the “Laplacian” form of the Shadowgraph sensitivity. Laplacian between quotes, since it does not depend on the derivative. We can plug in the Gladstone-Dale relation and the ideal gas law to isolate the pressure Laplacian:

Replacing the traveling wave form , now obviously considering the two-dimensional wavenumber in the directions perpendicular to the optical axis:

Where is the squared wavenumber vector magnitude, projected in the plane. Applying the same criterion of levels for minimum contrast, we get the following relation for the minimum detectable pressure wave:

Which we can easily convert to minimum based on the wave frequency :

design consequences

The boxed formula for shadowgraph is very similar to the one for Schlieren found in Part 1. By far, the most important point is the dependence on [Shadowgraph] instead of [Schlieren]. This is obvious, since Schlieren detects first derivatives of pressure, whereas Shadowgraph detects second derivatives, popping out an extra factor.

From a detector apparatus perspective, however, this is not of the most desirable features: The shadowgraph acts as a second-order high-pass spatial filter, highlighting high frequencies and killing low frequencies. The Schlieren has a weaker, first-order high pass filtering behavior. But this means they have a crossing point, a frequency where the sensitivity is the same for both. The apparatus parameters are different, however we can do the exercise for the rather realistic apparatus given in Part 1. If we make m (i.e., the focal plane is 0.5 meters from the Schlieren object), we get the following curve for a 12-bit camera:

Minimum detectable sound levels – Schlieren versus Shadowgraph. Schlieren parameters: m, mm. Shadowgraph has m. Both consider and m.

If the news weren’t all that great for the Schlieren apparatus, the Shadowgraph is actually worse: The crossing point for this system occurs at kHz. There might be applications where these ultrasound frequencies are actually of interest, but I believe in most fluid dynamics problems the Schlieren apparatus is the to-go tool. Remembering that, even though I show SPL values above 170dB in this chart, the calculations considered a constant-temperature ideal gas, which is not the case anymore at those non-linear sound pressure levels.

With this, I’m declaring Schlieren a winner for sound-vis.

Tell me: Would you be excited to see sound? Well, I sure damn was when I first attempted to see the sound produced by an ultrasonic transducer! And still am, to be honest! So let’s learn how to do it with a Schlieren apparatus and the sensitivity considerations necessary in the design of a Schlieren optical system. For this, I’ll use my handy Schlieren and Shadowgraph techniques by no one else than the great Prof. Gary Settles. My PhD advisor, Prof. Alvi, was his PhD student – so I should get this right!

I’m a practical engineer, but I’ll derive the formulas in this section. If you want, you can skip ahead to the next section where I apply them and produce nice design charts.

the math

If you would like more details on the physics, please refer to Settles. This is more of a derivation section. We’ll start with the Gladstone-Dale law, which relates the index of refraction of a gas to its density:

Where is the index of refraction, (for air) is the wavelength-dependent Gladstone-Dale constant and is the density of the gas.

Assuming a parallel light beam, the beam deflection of a uniform Schlieren object of length is given by:

Where is the direction of deflection. This deflection is magnified in the focal point of the optical apparatus mirror/lens close to the camera/screen; where the knife-edge is placed, according to:

Where is the change in the focal point location in the direction perpendicular to the knife edge and is the focal length of the second optical element (lens/mirror). Then, the contrast change observed in the camera/screen is given by Equation 3.6. in Settles’ book:

We can use the ideal gas law to combine the expressions as a function of :

We can then assume we are looking for sinusoidal sound waves of the form :

OK. Now let’s stop here for a while. I think it goes without saying you’ll probably need a high speed camera to see the sound waves, as they aren’t among the mundane things you can use your regular DSLR to observe. Given that, we can consider whichever camera you’ll use has some bit depth . When setting up your Schlieren apparatus, you’ll probably use half of the dynamic range of the camera for the “background” frame. This means, your intensity is about half of the dynamic range of the camera. In bit count, you’ll have as the background, undisturbed luminance level. Now, let’s assume you need at least levels to distinguish the signal from the noise. Maybe you have some de-noising method, like POD, that improves your SNR and your sensitivity. How loud do the sound waves need to be in order to produce a signal given you have a background luminance?

Let’s isolate :

Using our bit depth, we get the minimum pressure disturbance to be:

Converting to sound pressure level ():

Evaluating the and replacing :

A capital is used for frequency to eliminate confusion with , which in this context is the mirror focal length. If you know the background level, , then you might want to expose it in the formula instead of the bit depth of the camera pixels:

usage in design

The boxed formula is extremely useful for sound visualization design, as it allows us to define whether we can even see the sound waves we want to. First, a few remarks: I separated the bit depth because the appears in so many places. In this case, we have it with a negative sign, which is awesome, since in the end we want to minimize to maximize sensitivity. The first fraction inside the is a constant, dependent on the gas and the reference pressure. For air at 25ºC, .

Note that the minimum SPL is inversely proportional to the frequency . This means, higher frequencies are easier to see. Ultrasonic frequencies, in particular, are somewhat straightforward to see, as we will discuss. Let’s plot as a function of frequency for a typical large-scale Schlieren apparatus:

Sensitivity of an air-based Schlieren aparatus with m, mm, m and

The black line represents the bit depth of most contemporary scientific cameras. As you can see, the news aren’t that great: For low frequencies, we have to scream really loud to see anything – and this is after blowing up the image such that we can only see 5 levels! For example, for a frequency of 5kHz we have about 128dB of minimum sound pressure level to be perceived by our camera. It’s not impossible, and with the advent of new data processing techniques like POD it is well feasible. But wouldn’t it be great to have 4 more bits per pixel (red dashed line)? That would bring the minimum SPL to 104dB, which would enable the visualization of a lot of phenomena (and even a rock concert!).

Well, I hope this helps you design your own Schlieren apparatus. Or maybe you lose hope altogether and quit – but at least this saved you time! If you want to design your own Schlieren setup to visualize sound waves, you can download the code that generated the chart above here. Anyhow, here’s a little video of how sound waves look like for a really loud (~140dB, ~5kHz) supersonic jet:

So you’ve probably already seen demos on Youtube showing this really weird “camera effect” where they stick a hose to a subwoofer and get the water to look like it’s being sucked back to the hose, seemingly against gravity.

I personally love this effect. In the case of the subwoofer, the effect is due to what is technically called “aliasing”. Aliasing is an effect important to all sorts of fields, from data analysis to telecommunications to image processing. In technical jargon, you get aliasing when you don’t satisfy the Nyquist criterion when sampling your signal. This might not be accessible to everyone, so I’ll explain it differently.

In the case of the hose stuck to the subwoofer, the speaker shakes the hose back and forth with a single frequency (a single tone) and generates a snaking/spiraling effect on the water stream. If the water stream is slow enough (that is, if its Reynolds number is low enough for the flow to be laminar), then no weird stuff (non-linear effects) occur, and we get a simple, single-toned spatial wave in the water jet. That being the case, the cycles repeat very nicely, becoming indistinguishable from each other. If you can fulfill this criterion, then aliasing also occurs in a nice manner, that is, if you happen to fail to satisfy the Nyquist criterion, you don’t get a jumbled mess but a nicely backward or forward motion that looks like it’s in slow motion.

It is a simple thing to do, but there’s some beautiful fluid dynamics on it. Generating repeatable patterns and laminar flows is not that simple, especially when you are engineering a device. If you attempt to reproduce the video linked, you’ll find yourself suffering through a parametric search of flow rate/shake amplitude until you get the right combination that displays a nice effect.

Here, I’ll discuss a different device, though – I’ll talk about the piddlers of Dr. Edgerton, that inspired awe in many people around the world – including myself. I have never personally seen one. But I understood what it was and that the effect was not just a camera artifact, but something that could be seen with the naked eye because they use a stroboscopic light to show the effect to the viewer. I have not – as of 2018 – found any instructions on how to make these. And since it turns out it’s quite simple, I think it should be popularized. Here’s my take:

the piddler design considerations no one talks about in practice

The water piddler is a device that generates a jet of water droplets. Water jets are naturally prone to breaking down into droplets – if you’re a man you know it! But on a more scientific tone: Water jets are subject to a fluid dynamic instability called the Rayleigh-Plateau instability. This document here is an incredible source that enables the prediction of what are the conditions for this unstable behavior without hassling you with all the complicated math behind. The Rayleigh-Plateau instability looks like shown below:

Naturally occurring Rayleigh-Plateau instability in my kitchen

It’s beautiful – but it is also not really a single frequency. There seems to be some level of repeatability to it, but not enough to make the strobe light trick work. The reason for this non-repeatability is the following curve:

Dispersion relationship for the Rayleigh-Plateau instability

This is the dispersion relationship, extracted from equation (23) of Breslouer’s work. It corresponds to an axisymmetric disturbance in the radius of the jet – , where is a small disturbance in the radius and is the streamwise coordinate. gives us a wavenumber of the ripple in the jet, and gives us a frequency of this disturbance. The dispersion relation normalizes the wavenumber by the original radius of the stream, . When , we get an exponentially growing disturbance in the radius, which eventually makes the jet break down into droplets. So the black curve in the chart shows that any disturbance between will grow, but higher frequency disturbances will not – they simply oscillate. Disturbances closer to the peak at will grow faster, which is an important design guideline when we want to break down the jet into a stream of droplets.

The problem, though, is that any other frequencies around the peak also grow. The peak is somewhat smooth, so there will be a lot of space for non-uniformity, especially when the disturbances themselves start at different amplitudes.

So what would be a good design procedure? Well, first, we need to make sure the jet will be laminar. One way to guarantee that is to make the Reynolds number of the nozzle that makes the stream lower than 2000. That guarantees the pipe flow is laminar, which in turn makes the stream laminar. Of course, this is a little limiting to us because we can only work with small jet diameters. You can try to push this harder, since the flow inside the stream tends to relaminarize as the stream exits the nozzle due to the removal of the no-slip condition generated by the nozzle wall.

The other constraint has to do with reasonable frequencies for strobing. You don’t want to use too low of a strobe frequency, because that is rather unbearable to watch. Strobe frequencies must be above 30Hz to be reasonably acceptable, but they only become imperceptible to the human eye about 60Hz. We get a design chart, then:

Design chart for f=60Hz, water/air, 25ºC

The chart shows the growth rates (real part of ) for combinations of realistic jet diameters and velocities, which are the actual design variables. The line of constant Reynolds number looks like a 1/x curve in this space. The white line shows the upper limit for laminar pipe flow. You want to be under the white line, as well as in the growing region, which is about the region enclosed by the dashed line. For higher frequencies, the slope of the black boundary decreases, meaning you need smaller diameters to make the strobe light work. For lower frequencies, the slope increases, improving the available parameter space, but too low frequencies will be uncomfortable to watch. In case you want to develop your own piddler, a Matlab implementation to generate the colorful chart above is here.

It is actually rather remarkable that the parameter space looks like this, because feasible diameter/frequency combinations actually will break down into droplets if excited with 60Hz – the line frequency in the US. Say, for example, for 3mm jet diameter and 1m/s speed, we have a high growth rate and the piddler will produce a nice effect. At 6mm, 0.5m/s, we still have laminar flow but the instabilities won’t grow at 60Hz (lower frequency instabilities grow instead). Thus, you’ll not get a good piddler out of that combination. You might be able, for example, to push the bounds a bit (which I did) and make the jet diameter 4.75mm and the jet speed about 1.2m/s. In that case, the Reynolds number is about 3200, which still makes a reasonably repeatable piddler pattern.

Another thing you can attempt (I did) is to try to use a more viscous fluid. More viscous fluids will increase the viable diameter/velocity combinations where the Reynolds number is still low by pushing the white line up and to the right. For example, propyleneglycol allows us to approximately double the diameter of the pipe. The problem, obviously, is that it’s incredibly messy!

now to the real world

The design map is a good guideline to start this, but there are a few tricks that no one really describes in the internet. I’ll save you weeks of suffering: The best way to generate the disturbance is to use a coffee machine pump. Yes, a coffee machine pump! It accomplishes the two tasks for this device: Recirculating the fluid and generating a strong enough, single-frequency disturbance such that you don’t really need to trust the Rayleigh-Plateau instability alone to generate the droplets.

This is the basic implementation of the piddler I built. The coffee machine pump is a ULKA-branded one (like this one). I believe it doesn’t really matter which model you use, since mine had too much flow rate and I had to reduce the flow with the flow control valve indicated in the schematic. These pumps are vibratory pumps. They function by vibrating a magnetic piston back and forth with a big electromagnet. The piston pumps a small slug of fluid, that goes through a check valve inside the pump. When the piston returns, the suction shuts off the check valve, preventing back-flow. Suction fills the piston cavity of new fluid, and the cycle repeats.

Since the piston is pumping the fluid in tiny slugs at the correct frequency (assuming you have a viable Rayleigh-Plateau instability on your design), an acoustic wave will go through the water until the nozzle, generating the intended velocity disturbance in the mean flow. It will not be a choppy flow, but an oscillating one, due to the really strong surface tension forces in water. The figure in the left shows that there’s a full wave of instability before the stream breaks down into droplets.

Now that we discussed the design, let’s go for a little demo. In this video, I’ll also go through the z-Transform, which is a cool mathematical tool for modeling discrete-time control systems. I used this piddler as a prop for the lecture!