Jet actuator arrays, turning microjets into MIDIjets

So I’m currently working on this research problem: Microjets in cross flow for disturbance-based flow control. Jets in crossflow have some promise to be a viable flow control technique in aerodynamic applications, but it’s still in its early-mid research stages, where the technology has good theoretical support (i.e. it should work) and some experimental successes (it does work given several lab constraints, in very simple problems). Part of my thesis work will be to further the experimental support side of things.

Structure of a generic jet in crossflow (Source: Coussement et. al, 2012)

But when working with complex curved shapes (like any realistic aerodynamic surface) it is not clear where we should place a jet in the surface. Where is separation going to happen? Where to place the jets to prevent it from happening / make it happen earlier? Maybe we want to excite boundary layer waves, like the Tollmien-Schlichting waves? From the computational/theoretical standpoint, there’s some heavy-duty stability analysis that could potentially give possible “sensitive” locations for the jets. I’m not fully a computational person myself, but my current opinion given what I’ve seen so far is that we have too many assumptions we need to trust are approximate enough (i.e., Navier-Stokes are linear equations, jets produce content in the unstable eigenmodes of the flow,  we did resolve the relevant flow structures with the simulation results, cows are spherical, etc.). Again, this is not my specialty, so that’s probably why I find it hard to believe in the effectiveness of that approach.

But from the experimental (wind tunnel) standpoint, we need to drill a physical hole in the aerodynamic surface and route a pipe from inside to blow the jet. That requires some work, but more importantly it takes precious testing time when you’re testing your jet configurations. If all you were able to come up with were ineffective or mildly effective actuator patterns, that’s what you’re stuck with. And you’ll never know how close you got because you can only afford a few data points in the experiment. Furthermore, the background fluid dynamics knowledge required to come up with effective patterns requires decades of study and experience – which I don’t have. So I suggested: Why don’t we manufacture a reconfigurable actuator array and let a computer run thousands of pattern configurations? We could potentially abstract the jet placement problem from the fluid mechanics realm into a (rather complex, I admit) optimization problem. More jet configurations will be explored, increasing the confidence on the solutions found. And with the beauty of advanced flow diagnostics, we can even learn new physics from these solutions.

But then you might ask: Why don’t you just do CFD on these jets? Well, turns out that in order to perform any simulation work with jets in cross flow we need an obscene amount of resolution, which increases the computational time to an extent it is just easier to do the experiment. When it involves multiple jet configurations, you really need to be able to discard multiple runs, which requires them to be cheap. It’s a similar thing with AI. AI is only possible now because now each iteration is cheap to run, even though the math and the theoretical foundations come from several decades ago.

So this is the road I’m going through now, basically making microjet actuator studies cheap to run so we can discard most of them and try random stuff until we hit the jackpot. But, even though the prospect of having to build a reconfigurable array of jets with hundreds of jets may sound like a rather daunting task, there’s some fun along the path. And this is the point of this post!

Solenoid array under construction (62 solenoids shown here)

I’m building a manifold with 100 solenoids that can be individually controlled by a reconfigurable signal generator I designed and built (pictures below!). The signal generator board is based on a PIC32MZ (Design here) and has effectively 108 channels. I was able to update all channels simultaneously at 24900 Samples/s (well, there’s a 200ns delay between physical uC ports, but that’s virtually instantaneous from the mechanical standpoint). I designed it such that the board appears as a USB serial COM port in your computer, which then can receive messages through either a serial terminal or a serial interface on Matlab or C++, for example. This gives me a lot of control over the jets.

While putting all of this together and seeing the results of the system I built, I figured: Hey, I can turn this into a musical instrument! Of course, a rather crude one, because my bandwidth is very low (like <200Hz). But I decided anyways to code up a MIDI driver for this jet array and then change the notes to fit the bandwidth by shifting a few octaves on the song. The result is rather crude but it was so much fun to play with! MIDI files, for the uninitiated, are like a digital version of a sheet music. It contains a table of notes and the timing when they should be played and how long for. My job was to simply convert the digital instructions into the protocol I came up with for my serial communication and stream the instructions to the USB serial port.

So here’s a few songs I was able to play to a level where I believe people can actually reconize them: See if you can! (Answers in the description of the video). If you want to have more info on how I did it, perhaps you might consider following my research on ResearchGate and maybe a few years from now an academic paper on this topic will come from that! =)




Finding Vortex Cores from PIV fields with Gamma 1

Vortex core tracking is a rather niche task in fluid mechanics that is somewhat daunting for the uninitiated in data analysis. The Matlab implementation by Sebastian Endrikat (thanks!), which can be found here, inspired me to dive a little deeper. His implementation is based on the paper “Combining PIV, POD and vortex identification algorithms for the study of unsteady turbulent swirling flows” by Laurent Graftieaux, which was probably one of the first to perform vortex tracking from realistic PIV fields. The challenge is that when PIV is used, noise is introduced in the velocity fields due to the uncertainties related to the cross-correlation algorithm that tracks the particles. This noise, added to the fine-scale turbulence inherent to any realistic flow field encountered in experiments, makes vortex tracking through derivative-based techniques (such as λ2, Q criterion and vorticity) pretty much impossible.


Computational results are less prone to this effect of the noise and usually are tamer in regards to vortex tracking, though fine-scale turbulence can also be a problem. The three-dimensionality of flow fields doesn’t help. But many relevant flow fields can be “deemed” vortex dominated, where an obvious vortex core is present in the mean. Wingtip vortices are a great example of these vortex-dominated flow fields, though there are many other examples in research from pretty much any lift-generating surface.


As part of my PhD research I’m performing high speed PIV (Particle Image Velocimetry) on the wake of a cylinder with a slanted back (maybe a post later about that?). This geometry has a flow field that shares similarities with military cargo aircraft, but is far enough from the application to be used in publicly-available academic research. The cool part is that it forms a vortex pair, which is known to “wander”. The beauty of having bleeding-edge research equipment is that we can visualize these vortices experimentally in a wind tunnel. But how to turn that into actual data and understanding?


That’s where the Gamma 1 tracking comes into play. Gamma 1 is great because it’s an integral quantity. It is also very simple to describe and understand: If I have a vector field and I’m at the vortex core, I can define a vector from me to any point in this vector field (this vector is called by Graftieaux, “PM“). The angle between this vector and the velocity vector at that arbitrary point would be exactly 90º if the vortex was ideal and I was at the vortex core. Otherwise, it would be another angle. So if I just look at many vectors around me I just need to find the mean of the sine of these two vectors. This quantity should peak at the vortex core. That’s Gamma 1, brilliant!

Sebastian Endrikat did a pretty good job at implementing Graftieaux’s results, and I used his code a lot. But since each run I have has at least 5000 velocity fields, his code was taking waaaay too long. Each field would take 4.5 seconds to parse in a pretty decent machine! So I decided to look back at the math. And I realized that the same task can be accomplished by two convolutions after some juggling. A write-up of that is below:

PDF File with the math


The result, though, is really impressive. Now each field takes 5 milliseconds (3 orders of magnitude better!) to parse in the same machine. So good I made a video of the vortex core. Here it is:



I’m really thankful amazing people like Graftieaux and Endrikat are in the academic community publishing this stuff. Standing over the shoulders of giants!

The floating light bulb: Theory vs Practice

Yes – You can go to today and buy one of these gimmicky toys that float a magnet in the air. Some of which will even float a circuit that can light an LED and become a floating light bulb. A floating light bulb that powers on with wireless energy? What a time to live!

A quote from Arthur C. Clarke, who wrote “The Sentinel” (which later on became the basis for the science fiction movie “2001: A Space Odyssey”), goes along the lines:

“Any sufficiently advanced technology is indistinguishable from magic.”

This is what led me to the Engineering path. Because, if the advanced technology is indistinguishable from magic, who creates the technology is a real-life wizard. Who creates the technology? The engineers and scientists all around this world. So let me complement his quote with my own thoughts:

“Any sufficiently advanced technology is indistinguishable from magic. Therefore, engineers and scientists are the true real-life wizards.

Of course, if I’m writing about it is because I went through the engineering exercise. And boy, I thought it was an “easy” project. You see these projects of floating stuff around the internet, but nobody speaks about what goes wrong. So here we’ll explore why people spend so much time tweaking their setup and what are the traps along the way.

But first, some results to motivate you to read further:


Prof. Christian Hubicki was kind enough to let me pursue this as a graduate course project in Advanced Control Systems class at FSU, so I ended up with a “project report” on it. It is in the link below:

Full 11-page report with all the data

But if you don’t want to read all of that, here’s a list of practical traps I learned during this project:

  1. DON’T try fancy control techniques if you don’t have fast and accurate hardware. This project WILL require you more than a 10-bit ADC and more than 3-5 kS/s. The dynamics are very fast because the solenoid is fast. And you want a fast solenoid to be able to control the levitating object! Unless you can have a large solenoid inductance and a rise time in the order of ~100ms, there’s no way an Arduino implementation can control this. I’m think a nice real-time DAQ controller (like the ones offered by NI) could work here. But an Arduino is just too strapped in specs to cut it! The effects of sampling and digitization are too restrictive. It MIGHT work in some specific configurations, but it is not a general solution (and certainly it didn’t work for me).
  2. Analog circuits are fast – why not use them? Everyone (in the hobby electronics world) thinks Arduino is a silver bullet for everything. Don’t forget an op-amp is 100’s of times faster than a digital circuit!
  3. Bang-Bang! You see many implementations in the web use a hysteresis (or bang-bang) controller. The bang-bang controller is ideal for cheap projects because it deals well with the non-linearities gracefully. But it is not bullet-proof either. It will become unstable even with high bandwidth in some cases if the non-linearity is strong enough.
  4. Temperature Effects: The dynamic characteristics of your solenoid will change as it heats up (you’re dissipating power to turn it on!). So it can get very confusing if you have, say, a PID controller, to tune the gains because the gains will be different depending on the temperature of the coil. Since this effect is very slow (order 10 minutes!) it can result in you chasing your own tail because you’re tuning a plant that is changing with time!
  5. The wireless TX introduces noise! This one is particular to this project: If you’re using a Hall effect sensor to sense the presence of the floating object (by its magnetic field), then your Hall sensor will also measure the solenoid strength! Apart from that, the TX is also generating a high-frequency magnetic field, which will also be measured in the Hall effect sensor signal. The effect of the TX is very small (~2mV) buy it appears in the scope. The problem is that Arduinos don’t have low-pass filtering in their ADC inputs. So anything higher than the sampling rate will appear as an “Aliased” signal, which is very nasty to deal with.
  6. Make sure your solenoid can lift your object and more. This is an obvious one but I think it is easy to overlook because you need to over-design it. I designed my solenoid to lift 100 grams of weight. But in the end, I could only work with 35 grams because the controller needed a lot of space to work. So overdesign is really crucial here. I ended up shaving a lot of mass from the floating object because I couldn’t lift the original design’s mass!

I’d like to put a more complete tutorial on making this, but since I already invested a lot of time in putting the report together, I think if you put some time on reading it and the conclusions from the measurements/simulations you will be able to reproduce this design or adapt the concepts to your design. Let me know if you think this was useful or maybe if you need any help!

Sub-microsecond Schlieren photography

(Edit: My entry on the Gallery of Fluid Motion using this technique is online!)

For the ones not introduced to the art of Schlieren photography, I can assure you it was incredibly eye-opening and fascinating to me when I learned that we can see thin air with just a few lenses (or even just one mirror as Josh The Engineer demonstrated here on a hobby setup).

For the initiated in the technique, its uses are obvious in the art and engineering of bleeding-edge aerodynamic technology. Supersonic flows are the favorites here, because the presence of shock waves that make for beautiful crisp images and help us understand and describe many kinds of fluid dynamics phenomena.

Schlieren image of a 2mm supersonic microjet taken at Florida State University FCAAP laboratory. Illumination time is 500 nanoseconds, taken with a Nikon D90 DSLR to demonstrate the potential for hobby applications. Note the crispiness of the image – the flow was effectively frozen.

What I’m going to describe in this article is a very simple circuit published by Christian Willert here but that most likely is paywalled and might have too much formalism for someone who is just looking for some answers. Since the circuit and the electrical engineering is pretty basic, I felt I (with my hobby-level electronics knowledge) could give it a go and I think you also should. I am also publishing my EasyEDA project if you want to make your boards (Yes, EasyEDA).

But first, let’s address the elephant in the room: Why should you care? Well, if you ever tinkered with a Schlieren/shadowgraph apparatus – for scientific, engineering or artistic purposes -you might be interested in taking sharper pictures. Obtaining sharper pictures of moving stuff works exactly like in regular photography. They can be achieved by reducing the aperture of the lens, by reducing the exposure time or by using a flash. The latter is when a pulsed light source really shines (pun intended!). The great part here is that the first two options involve reducing the amount of light – whereas the last option doesn’t (necessarily).

The not-so-great part is that camera sensors are “integrators”. This means they measure the amount of photons that happened to be absorbed given an amount of time. Therefore, what really matters is the total amount of photons you sent to the camera. Of course, if you sent an insanely large amount of photons in a very short instant, you would risk burning the camera sensor – but if you’re using an LED (as we are going to here), your LED will be long gone before that happens.

So the secret for high speed photography is to have insanely large amounts of light dispensed at once. That would guarantee everything will be as sharp as your optics allow. Since we don’t live in the world of mathematical idealizations, we cannot deliver anything “instantly”, and therefore we have to live with some finite amount of time. Brief enough is relative and depends on what you want to observe. For example, if you’re taking a selfie in a party, probably tens of milliseconds is brief enough to get sharp images. For taking a picture of a tennis player doing a high speed serve, you’re probably fine with tens or hundreds of microseconds. The technical challenges begin to appear when you’re taking pictures of really fast stuff (like supersonic planes) or at larger magnifications. The picture of the jet above is challenging in both ways: its magnification level is 0.7x (meaning the physical object is as projected in the sensor at 0.7x scale) and its speed is roughly 500 meters per second. In other words, the movement of the object (the Schlieren object) is happening at roughly 63.6 million px/second, which requires a really fast shutter to have any hopes to “freeze the flow”. If you’re fond to making simple multiplications in your calculator, the equation is very simple:


Where D is the object displacement in px/second, v is its velocity in physical units (i.e. m/s), M is the magnification achieved in the setup and s_{px} is the physical pixel size of your camera (i.e. s_{px}=5.5 \mu m for a Nikon D90).

I know, I know. These are very specialized applications. But who knows which kinds of high speed photography is happening right now in someone’s garage, right? The point is – getting a light source that is fast enough is very challenging. Some options, such as laser-pulsed plasma light sources, can get really expensive even if you make them yourself. But LEDs are a very well-established, reliable technology that has an incredibly fast rise time. And they can get very bright, too (well… kinda).

So what Willert and his coauthors did was very simple: Let’s overdrive a bright LED with 20 times its design current and hope they don’t explode. Spoiler alert: Some LEDs didn’t survive this intellectual journey. But they mapped the safe operational regions for overdriven LEDs of many different manufacturers. To name a few: Luminus Phlatlight CBT-120Luminus Phlatlight CBT-140Phillips LXHL-PM02, among others. These are raw LEDs, no driver included, rated for ~3.6-4V, and are incredibly expensive for an LED. The price ranges from $100 to $150, and they are usually employed in automotive applications. The powerful flash is, however, blinding. And if they do burn out, it can be harmful for the hobbyist’s pockets.

LED driver power section.

The driver circuit (which is available here) is very simple: An IRF3805 N-channel power MOSFET just connects the LED to a 24V power supply. Remembering the LED is rated for 4V – so it’s going to get a tiny bit brighter (sarcasm). Jokes apart, the LED (CBT-140) is rated for 28A continuous with very efficient heatsinking, which means we will definitely be overdriven. By how much we can measure with R2. Hooking a scope between Q1 and R2 is not harmful to the scope and allows to measure the current going through the LED (unless the currents exceed ~600A, then the voltage spike when the MOSFET turns off might be on the few tens of volts). We don’t want to operate at these currents anyways, because the LED will end up as in the figure below. There’s a trim pot (R3) that controls the MOSFET gate voltage, make sure pin 2 of U1 is giving a low voltage when tuning.

A sacrifice for science.

What is really happening is that C1 and C2 (C2 is optional) are being charged by the 24V power supply when the MOSFET is off. Then they discharge at the LED when the MOSFET is activated. No power supply will be able to push 200A continuously through an LED, so if the transistor turns on for too long, the power supply voltage will drop and the power supply will reset. Actually, this is one of the ways to tell if you melted the MOSFET (which happened to me once). The MOSFET needs to turn on in nanoseconds, which will require a decent amount of current (like 4-5 amps) just to charge the gate up. This means we need a driver IC – which in this case I’m using a UCC27424. Make sure to have as little resistance between the driver and the gate to minimize the time constant. The 1.5 Ohms is very close to giving 4A to the MOSFET. Since the gate capacitance is around 8nF, the MOSFET gate rise time is somewhat slow (12 ns).

Speaking about time constants, during the design I realized the time constants of the capacitor that discharges into the LED and the parasitic inductances in the path between the components will dictate the rise time of the circuit, at least for the most part. In my circuit, the time constant was measured to be 100ns, directly with a photodiode. This means we can do >1MHz photography, which is pretty amazing! Unfortunately the cameras that are capable of 1 million frames per second aren’t really accessible to mortals (except when said mortals work in a laboratory that happens to have them!).

Well, the LED driver circuit is still in development – which means I’ll keep this post updated every now and then. But for now, it’s working well enough. The BOM cost is not too intimidating (~$60 at Digikey without the LED. Add the LED and we should be at ~$200), so a hobbyist can really justify this investment if it means an equivalent amount of hours of fun! Furthermore, this circuit implements a microcontroller that monitors and displays the LED and driver’s temperature. It features an auto shut-off, which disables the MOSFET driver if the temperature exceeds an operational threshold. The thermal limits are still to be evaluated, though.

Circuit board and a (crude) 3D printed case for the LED.

For now, I did my own independent tests,  and the results are very promising.  Below I’m showing a test rig to evaluate the illumination rise and fall times of the LED. The photodiode is a Thorlabs (forgot the model) that has a 1ns rise time if attached to a 50 ohm load. It’s internally biased, which is nice when you want to do a quick test.

Test rig for photodiode illumination response measurements

The results from the illumination standpoint are rather promising. Below a series of scope traces show that the LED lights up in a very short time and reaches a pretty much constant on state. The decay time, however, seems to be controlled by a phosphorescence mechanism that is probably because this is a white LED. Nevertheless, the pulses are remarkably brief.


Scope screen for LED illumination (blue curve) as seen by the photodiode. Yellow curve is current as measured by a 10mOhm resistor in the MOSFET source. Curves from 100ns, 300ns and 1000ns input pulse width, respectively.

The good thing about having high speed cameras is that now we’re ready to roll some experiments. By far, my favorite one is shown below. I was able to use the Schlieren setup to observe ultrasonic acoustic waves at 80kHz , produced by a micro-impinging jet (the jet is 2mm in diameter). The jet is supersonic, its velocity is estimated to be 400 m/s. Just to make sure you get what is in the video: The gray rectangle above is the nozzle. The shiny white line at the bottom is the impingement surface. The jet is impinging downwards, at the center of the image. The acoustic waves are the vertically traveling lines of bright and dark pixels. I was literally able to see sound! How cool is that?



Just as a final note. You might be discouraged to know that I am one of these mortals that happen to have access to a high-speed camera. But bear in mind, these pictures could have been taken with a regular DSLR. The only difference is that the frame sequence wouldn’t look continuous, because the DSLR frame rate is not synchronized with the phenomenon. Apart from that, everything else would be the same. You should give it a try! If you do, please let me know =)

My POV Display V1.0 – Wireless Powered

In a previous post I spoke about wireless power transfer and some engineering I was trying to do with it. This project proved to be very effective, and I’m quite happy with the results I got with it. I’m quite certain someone already had this idea, but the hobbyist endeavor is really useful anyways! And we all know the world is way too big, so sometimes it’s better not to bother to be innovative and just enjoy the build. Anyways, I thought maybe I could give my two cents to the DIY hackers community!

I wrote at Instructables (here!) a simplified guide to making this display. There’s no reason to cover it again in the blog, so I’ll just dive further here for those who want to play along at home (as Dave Jones from EEVBlog would say…). I’ll nevertheless talk about the build because that would be a shame not to be in this blog anyways:

1. The Design

So, as any other design, I begun with the question: What can I do with the process I have at my disposal and the budget I have? I knew a guy that had a laser cutter (he would charge me anyways, but at least I knew him). Also, I inevitably would have to lathe at least one part, because I don’t own a 3d printer and I need to connect my motor shaft to my display. So I set the budget to be ~100 USD (plus whatever I already had) and I begun with the overall characteristics:

“How many pixels I want?” – Well, I contacted a PCB supplier that told me he would cut 200x200mm PCBs no problem. I wanted to have the LEDs directly in the PCB, so the restriction was set: Radius of ~85mm -> Angle of ~150 degrees (no point in having the south pole as it would be covered). Run the numbers for a 5mm LED and we get ~40 LEDs max. Any electrical restriction? I had some PIC16F877As lying around. They have 4 full 8bit ports, which makes it really useful to easily toggle all LEDs at once in the routines (I’m not that great of a programmer, although I like to use the PIC!). I ended up using 4 ports, but only 6 bits each port. That was a total of 24 LEDs. It’s not great resolution, I know, but it’s a first prototype so bear with me ^^

“Update rate?” – Well, at least 20Hz to be convincing, right? I had an old fan motor lying around (haha!). I didn’t know its speed though. But I have a function generator, so I made a quick stroboscopic tachometer and measured its speed. The tachometer is just an LED and a resistor connected directly to the function gen’s output – Adjusting the frequency until you see the motor freeze. Then make sure that 1/2x the frequency doesn’t freeze the motor. It turned out to spin at 30Hz. Unfortunately, after loading it spun only at ~24.5Hz.

“Can I update these pixels fast enough?” – The individual signal rate would be easy to deal with – 30Hz * ~180px per revolution = 5.4kHz. How many processor cycles do I have to perform my calculations? (20MHz/4)/5.4kHz=926. That seems to be enough. But I still need to nail that update frequency (5.4kHz) quite well, to less than half a pixel, so the effect actually is convincing (otherwise there’ll be too much jitter in the image). That means, 5400±7.5Hz. Or, between 924.6 and 927.2 cycles. plus/minus 1.5 cycle. That’s more challenging. I mean, it can be done, but I think that’s possibly the limit for this PIC frequency. Unfortunately to have a faster clock I need to buy another microcontroller, which I don’t want to. So let’s give it a try =)

“How much power do I need if the display is fully on?” – That’s an easy one: 24 LEDs at 5V, 20mA – 2.4W. On the transmitter side, probably close to 10W given the inefficiency of a homemade device. The IRF630 should deal with this no problem.

Now that I had some idea of what I wanted to do and what would the limitations be, I dove into the PCB design. In my design the PCB spins and it looks like:


I confess routing the LED connections was quite difficult. I’m glad I chose not to use 52 LEDs now! For more LEDs sure I would need a 3D arrangement, or some sort of serial system. But for a beginner let’s keep it easy.

It was quite straightforward to design the case around this PCB. I chose to cut most pieces out of acrylic, but as I said – one of them had to be lathed. The connection between the motor and the spinning part. The design looks more or less like that:


It’s sad but at that time I didn’t harness the power of the Prusa. Thus, I couldn’t print the shaft connection. The entire project would be different had I had a 3d printer by then.

Design in hands, shot some e-mails to quote the parts. the costs were (as of 2017):


Yes, electronics are expensive in Brazil. And yes, hardware is cheap.


2. The Build

As I received the parts I began to build. Apart from hours of building, I didn’t actually have any trouble with it. Slowly building and checking the parts of the circuit were working as supposed, it wasn’t really eventful. The worst thing that happened was one open trail at one LED (which means the PCB manufacturer wasn’t really great), that I had to bridge. Below are some build pictures:



That was hours of fun, as you might be wondering. I know it’s not the best of its category, but it’s my child!


3. The programming

That’s the boring-to-talk about but fun-to-do part where you get to use the device you built and actually bring it to life. Yes, the issue of having a precise timing actually appeared during the programming. Also, it turned out the motor didn’t have a constant speed, but it oscillated between 24 and 25Hz. Fortunately I had a zero datum activated by an IR LED. That helped a lot but didn’t completely solve the jittering of the motor. Below are some “failed” tests while I was programming it:

I’m sure you want to see the finished product working. So there it is:


4. Conclusion

Well, as you can see this is more of a picture repository than anything else. If you want to take a closer look at the project, please check my Github page at The project files are all stored there.

I’d say the biggest lesson learned from this project is that processor speed is crucial for any display device. To update so many pixels, there’s just not enough time. I built an increased admiration for our high-resolution full HD screens. What an engineering feat!





Resonant Coupling Efficiency for Wireless Power Transfer

So I’ve been working on a hobby project that involves transferring power to a circuit wirelessly. The distance is quite short (<10mm), as I just need to have one part of the device spinning and the other stationary.

I decided to stick with LC resonant coupling circuits as a solution to transfer this power, as I can admit for some power loss and having low audible noise is a requirement. But while I was studying this circuit, I wanted to know what would be the best conversion efficiency attainable, and also I was interested in understanding why in the first place a resonant circuit is needed (why not just use a regular circuit like a transformer?).

To answer these questions, I was thinking of building the circuit model in a simulator (as in the picture shown below), and simply playing with the parameters should give me a good grasp on what was happening.

Circuit analyzed in this article

While it did help me understand some concepts, it still was not very clear what was the effect of the coupling coefficient (k) on the power dissipated in the load and the overall power transfer efficiency of the system. So I decided to step up a little bit and use Matlab. The simple code I wrote is in this link and you can also play with the parameters for your circuit.

Basically, I used the control systems toolbox to write the transfer function of the circuit. The transfer function I’m interested on, eta, is the power dissipated in the load ZL divided by the power generated in the source Vin. The remaining power should be dissipated in the source itself, through its own impedance ZS.

The plot of frequency response for many coupling factors (as shown below) is rather interesting: It shows there is indeed a dependence of peak power transfer efficiency as the coupling factor increases. But it shows also that you don’t need a very large coupling factor to get interesting efficiencies (for k=0.4, eta peaks at 77%, for example). This means we don’t need astonishingly good coupling (as in a transformer, for example, where k>0.98) to transfer power. I know, this conclusion has already been made by smarter researchers than me long ago, but it’s nice to arrive at it by myself anyways.

Frequency domain response of the power transfer efficiency for ZL=20Ω, ZS=1Ω, C1=C2=100nF, L1=5.5µH, L2-6.0µH.


The following picture shows the same information, but complies the peak efficiencies against the coupling factor k. It shows that below a certain coupling factor (~0.4) the power transfer efficiency plummets. Unfortunately, this means that for large distances the resonant circuit is not the solution (as everyone already knows!), since the coupling factors are absurdly small for these distances.

Even so, for my simple peasant application of transferring power through a 10mm air gap that’s more than enough enlightenment in the topic.

Peak efficiencies for the previous circuit parameters plotted versus k.

Just for the record, I measured the inductances of my homemade pair of coils (5.5µH and 6.0µH, as stated in the figure) and the coupling factor against distance (no capacitors, just the coils, using the open circuit voltage method, f=50 to 800 kHz) and the results are as shown below. According to the charts I generated, then, I should be fine, right? Well, it will depend whether I can make a MOSFET driver that will be able to attain a nice efficiency factor to multiply these maximum eta values I obtained!

Homemade coils and the coupling factor measurement results. The distance d is the gap between coils (I used a known thickness wood spacer)





A freezing tunnel (VRT) thermal model

As many refrigeration engineers around the world, I’ve always been curious to understand better how this machinery behaves transiently. For the ones who do not know the Variable Retention Tunnel, take a brief stop at Google before reading on, so you can appreciate better what I’m talking about!

In practice, we see that when the production rate increases, the heat load rises and the compressors might struggle to keep the evaporation temperature low. Therefore, the evaporator temperature rises, which makes the air temperature inside the freezing tunnel rise as well. The system is constantly looking for some equilibrium point, where the evaporator temperature fits for the current product mass flow rate.

But as each product cools down and eventually freezes, the heat transfer at its surface slows down because the surface temperature drops. This means some of the product inside the tunnel is not contributing to the heat load as much as the fresh, hot product that just entered the device. And all of this is happening simultaneously, which makes the appreciation and understanding of the process remarkably difficult.

So, recently while working on a project I had the idea: Why don’t I simulate the system so I could grasp better what is actually happening? I began, then, an endeavor towards this objective, and the concepts I could get in the path were quite enlightening.


1. The product comes in waves

Probably you are familiar with the product freezing curve shown below. It tells us the temperature evolution in time for some unfrozen product that is subject to a freezing process. The surface freezes first, then the low temperatures slowly (depending on the product’s thickness) reach the center of the product.

Typical food product freezing curve

Another thing this curve also tells us, although indirectly, is the power the product releases onto the environment air. As q'=h*A*(T_s-T_\infty), and assuming the convection coefficient and the product surface area don’t change, we can assure that the power released (q') by the product is proportional to the red curve in the figure.

The implications of this are quite interesting: As the curve shows, the temperature difference (and therefore, the power) drops by half as much in the first hour. So the product contributes the most to the heat load in the first hour, and then it “locks” the inner heat inside the product, making the heat transfer process much slower. As product is constantly coming in the tunnel during the loading hours, the effect looks like a wave of products, as shown in the video below:

The tunnel in this video has a 12h minimum retention time for any given product. This means during the working hours it is almost always full of “hot” product. When product stops coming, then it gradually stops “heating up” the tunnel and global heat load drops. But how does this translate in the total heat load transient profiles?


2. The daily oscillating heat load profile

Unfortunately, in practice we cannot measure the total heat load except with quite expensive equipment (and even so it would be quite inaccurate). But after simulating a whole tunnel, one can just integrate the data and get the total instantaneous heat load at each instant simulated.

Of course, the heat load profile will depend on the production schedule and on how long and when the compressors are turned off to save energy. In some countries, there is a strong variation on the energy fees depending on the hour of the day and on the season of the year. These variations are placed so the energetic system is not overloaded during the hours most people are home using electricity (normally early in the night), and companies are stimulated to turn off their machines when the fees are higher.

These factors make the heat load profile very rich and detailed. Below there is a picture of the same tunnel’s simulated profile during five days. The profile looks very repetitive, with the same peaks and troughs daily. During region A in the graph, the tunnel is empty and heat load builds up as product enters the tunnel. The heat load drops to zero in region E because the fans and the machine room are turned off between 18:00 and 21:00. Since the production is still ongoing during these times, the heat load peaks to B briefly during the machine room retake. This far exceeds the capacity of the compressors at the setpoint, and the machine room will struggle during this time. As product stops being fed into the system (at 0:00), the heat load stabilizes and, during the off-production hours (from 0:00 to 8:00, it drops until point D in the graph. This cycle runs around every day, and during the weekend break the power drops to zero (remaining only the conduction through insulation and the fan power). The weekly cycle begins again.

Simulated heat load profile of a typical VRT installation in Brazil.


3. The machine room struggles… But sometimes that’s okay

I thought that instead of assuming the machine room was “sized” to properly maintain the evaporator pressure at all times, I’d rather have it to “struggle” maintaining it at the setpoint during these high heat load periods of the day. “How would this affect the product?”, I asked myself. So let’s examine the chart below:

The tunnel starts the week, again, with virtually no heat load as the product’s temperature is already homogenized with the air temperature for staying so long during the weekend (region A). As product slowly loads the tunnel, the air temperatures (region B) rise but the machine room is still partially loaded and still keeping the evaporator temperature. In region C, the heat load of the product exceeded the machine room heat load, increasing the pressure of the evaporating fluid (and therefore, T_{evap}) even though the compressors are working at full load.

Simulated Evaporator and Air Temperatures for the same VRT installation.

In region D, the compressors are turned off for 3 hours, and the temperature readings are virtually meaningless. After turning on the compressors (region E), the huge heat load hits the machine room hard and makes the evaporator’s pressure rise all the way to a saturation temperature of -26ºC. Then as the heat load slowly settles down to a smaller value, the compressors begin to catch up.

This is a beautiful picture of the freezing tunnel and one that corresponds quite closely to what I saw in my experience. But if the compressors cannot keep the evaporator setpoint temperature, isn’t it a compressor sizing problem?

Well… It’s complicated. On one hand, the compressors normally are sized for a heat load that meets some average or “guessed” instantaneous criterion. Let’s say, for example, that one sizes the machine room capacity to meet the instantaneous heat load of an hour. The data used for this simple calculation is:

\dot{m}=24000 kg/h

\Delta H=275 kJ/kg

ExtraPower (OtherDevices)=257 kW

Safety Factor=f_s=1.15

The total power required for these compressors would, then, be:

P=f_s * (\dot{m} * \Delta H + Extra Power) = 2395 kW

Which corresponds to the magenta dotted line in the heat load chart. As the chart shows, the magenta line is below the peaks each and every day. But does this mean the compressors are undersized? This would mean they should actually be sized to a 9MW capacity, that occurs for a short half-an-hour. Is it worth it?

Well, it turns out, in practice, it’s not. The whole system has quite a lot of inertia. This tunnel stores 300 tons of product, which slows down the reaction time. Therefore, some transients can be accepted. But how far can we go?

Let’s say that another designer “wants” to size the same system for a daily average heat load. The system design heat load now will be reduced to 2395 kW*(21h/24h)=2095 kW. Let’s say we size our compressors for this heat load. This would mean we now are 300 kW (12.5%) shorter in compressor juice than we were before. So, obviously, we’re going to stay for longer in the “undersized” region, like the chart below displays. But what effect does this have in the product temperatures?

New evaporator temperature chart after reducing the compressor capacity by 12.5%

Well, let’s take a look and see what the product does. As different products stays different times in the tunnel (I’ll talk about this), then some products will exit at a lower temperature than others. This means there is some “statistical distribution” of the products, which is displayed in the next chart.

The chart actually shows the integral of statistical distribution, shows as a percentile. Let’s take the point (-30.2 ºC, 40%) in the upper chart as an example of how we read this chart. This point says that if I take a sample of the lowest 40% of the products outlet temperatures, the largest temperature would be -30.2ºC. In other words, 40% of the product is at -30.2ºC or lower.

Effect of a 12.5% reduction in compressor capacity on the product outlet temperatures (center)

The two charts above show that when we reduce the machine room compressor’s capacity, the average outlet temperature at the center of the products will be larger. But we need to compare it with a threshold, which is shown as a black dashed line in the chart (-18ºC). In the upper case (larger compressor), 99.6% of the product meets or exceeds the -18ºC threshold line. In the lower scenario (smaller compressor), a lower portion (97.1%) of the product meets the threshold line. To define whether these numbers are satisfactory or not is a decision the plant operator must make.


4. But I, the designer, shall dictate the capacity of the system!

The takeaway from the previous section is that the heat load of a given system (and this applies to other thermal systems as well) is not something we “decide”, but a product of the operational variables. Of course, we can control most of the operational variables, but we shouldn’t assume that’s always the case.

Another designer might say, for example, he “wants” to size the compressor to the weekly averaged heat load. He argues that the weekend time is idle, compressors are shut off unnecessarily during this time and he ought make the initial investment lower by putting more machines that will operate the whole week, therefore saving millions of dollars! Is that possible? Let’s see…

So the figure below shows how the temperatures behave for this compressor sizing. We can see that the oscillations are wilder, the peaks are roughly 5ºC above the previous simulations, and the machine room saturation temperature is constantly oscillating according to the tunnel heat load. But still at the weekend, the heat load drops and the compressors need to be shut off.

Temperature profile for a 1496 kW compressor (weekly sized).


So no improvement happened here, the compressors are still struggling during the week and being shut off during the weekend. Why? Well, that’s because we, as I previously mentioned, cannot dictate the heat load. The system does that. We just size for it.

So you’re probably curious to see how did the product temperatures perform under these “undersized” conditions. As shown in the figure below, we have now 75.4% of the product meeting the -18ºC center temperature criterion. It’s not remarkably bad, actually – but some clients will not be remarkably happy with some -10ºC measurements in their quality control spreadsheets. That’s the reason why we oversize!

Effect of a 37.5% reduction of compressor capacity in the product temperature statistical distribution.


5. Conclusion

This is sort of a results article. I believe I’ll write some other articles on the mathematical devices I used to put the simulator together. I hope it was enlightening for you, the reader, as much as it was for me. If you think that kind of simulation would be useful for you, feel free to contact me.

Thanks for reading so far!

Mass file renaming and relinking in Inventor Part I: The results

*Edit: New, translated, version posted on GitHub with install instructions. Check it out at


1. Introduction

If you’re a mechanical designer for a living, chances are that you already had to copy an existing design. Even higher are the chances that you already had to rename all the files in a given assembly to match an ERP needs or a coding system needs.

Well, there are many ways of doing this, but I found that all the built-in ways in Autodesk Inventor have some sort of caveat:

  • Manual renaming in explorer: This is a pain in the ***. Apart from renaming the files one by one, you have to keep track of the old-new name correlations so when you open the assemblies in Inventor, you can redo the links manually.
  • Excel VBA to rename files: This is a little bit less cumbersome. You relieve the pain of renaming the files one by one, but the relinking still needs to be done manually in Inventor. At least you can print the correlation table from Excel.
  • Using the Design Assistant: This tool is nice, but there is a missing functionality that keeps it from being something usable for large assemblies: One cannot paste a table in the treeview. Imagine you as a designer already having the old-new naming correlations in an Excel table, but unable to paste them in Design Assistant. You still have to do the renaming manually. Sometimes it hurts productivity even more as the files are organized according to the tree, so finding them is another issue.
  • Using the CopyDesign add-in sample from Inventor SDK: This is an example piece of software that comes with the Inventor API SDK. It does what it’s supposed to do: Copies a design and renames the files. The problem is, you need to have the new names beforehand. And new files (which are probably going to happen in a copied design) will have to be named correctly during design, which is not always possible under an ERP system.

2. My solution

So, now that I showed the available options and how they did not meet my needs, let me state what I want to do and what I came up with using the Inventor API.

I want to make a file list using Excel VBA (which will list the file names for the current assembly) and, using some made-up rule, assign new names to the existing files. Here I think Excel will always be way more flexible than any piece of software I would ever write, so I decided to leave the rule creation to the freedom of the user. A table like the one shown below will be the output from Excel as a CSV file and read by this program:

Input CSV file format for the renaming system

I open the main assembly file and import the correlation table on my API-brew software. It figures out what files need to be renamed if they match the name in the “from” column of the correlation table. It finds out whether there are any equally-named *.idw files in the path. It renames all the files and then redoes the links in all levels of the assembly. Finally, it saves the assembly.

That’s the procedure in a nutshell. The heavy work is done by the software. So how does it look in Inventor? Follow the step-by-step screenshots shown below.

2.1. Main assembly loaded in my API interface

***Edit: I changed the button name to “Batch Renaming“.

Screen of the “Batch Renaming” add-in showing the list of components referenced by the assembly

2.2. Loading the CSV file

CSV file being selected in the file browser

2.3. Correspondences loaded in the file being displayed on-screen

In this case, not all files were selected to be renamed as some parts belong to the library.

Files to be renamed are shown in the table in bold

2.4. Renamed files in Explorer

Files after the renaming in the explorer view…
… and in the Inventor treeview!

3. Conclusion

So the use of this tool saved us quite a lot of time on the renaming of files. For larger assemblies, of >700 parts, a reliable manual renaming took up to 3-4 hours and could be reduced to something like 10 minutes. If you think this tool might be useful for you too, please let me know in the comments!

I published it on GitHub in this address: Please check it out and let me know if it was useful for your job!

If you guys understand Portuguese, you might find this video I made in my language for my comrades useful. I’ll soon produce an English version with the updated code so more audience can understand it. **EDIT: there it is =)

**Edit 2: Watch the installation process if you’re having trouble:

The freezing heat exchanger: Part III – Some insight

So in Part I and Part II we spoke about how to solve numerically the freezing heat exchanger. After some tinkering with the results I got, I realized something very simple: The maximum heat transfer rate in a freezing heat exchanger is not dependent on:

  • The ice thickness;
  • The non-freezing fluid heat transfer coefficient;
  • The material/thickness of the tube wall.

Which is quite unintuitive! Basically, the maximum heat transfer rate of a freezing HX happens when it is completely frozen. And then that happens, we know what’s the temperature of the liquid boundary: T_{wall}=T_{freeze} (Of course, assuming no supercooling).

This allows us to analyze a freezing heat exchanger based on the maximum heat transfer rate possible, in a way similar to the ε-NTU method. Let’s call it the ice-efficiency method, and name an ice-efficiency factor \varepsilon_{ice}, defined as:


Where q_{HX} is the effective heat transfer rate of the heat exchanger, while  q_{fully-frozen} is the heat transfer rate of the same heat exchanger, should it be fully frozen. Of course, q_{HX}<q_{fully-frozen}.

Actually, a similar approach has already been suggested as a “dimensionless heat transfer rate variable” q^* by Zerkle (1964) in his PhD thesis, so it’s not something really that innovative. The difference here is that I’m using the maximum heat exchangeable by the HX the size it is, instead of the maximum heat exchangeable by an “infinitely long” HX.

q_{fully-frozen} can be calculated using the ε-NTU method, assuming T_{wall}=T_{freeze}. First, let’s define q_{max}=q_{HX}/\varepsilon:


And the number of transfer units:

NTU=\frac{h_i A}{C_{min}}

An important thing to notice in the NTU equation is the fact that U, the global heat exchanger coefficient, has been replaced by h_i, the internal convection coefficient. This is because given a wall temperature (T_{freeze}), the heat transfer rate is now independent on the other thermal resistances. What this means is that, once fully frozen (which is the best case), the maximum heat transfer possible is only dependent on the liquid convection intensity.

Two takeaways come from this interpretation: The first is that (mild) solidification is desirable in a freezing heat exchanger. The second is that this sets up the operational limit. Of course, different from the ε-NTU method, the way to find \varepsilon_{ice} is not so clear. Maybe we can come up with a correlation for that?

The freezing heat exchanger: Part II – Solution

1. Introduction

So in Part I (link here) we devised a method to possibly solve the freezing heat exchanger problem. Here we’re going to look into the numerical method itself and what are the results we can get with this tool.

2. Solving it numerically

We began by stating the fundamental equations of this 2D flow, which need to be integrated in three dimensions (two space dimensions and one time dimension). By eliminating the need to solve for the flow field, the integrator will only need to solve for one space and one time dimension, greatly simplifying both the integration method needed and the number of computations to be made.

The solution method I devised is explained graphically in the diagram below:

Integrator diagram to solve for the freezing HX

This is the basic skeleton of the integrator. We can add a lot of stuff there to adapt the solution to our needs. For example, in my case I added the dynamic behavior of the shell side, so I could implement a controller for a motorized flow valve. This allowed me to understand the control behavior and how the HX responded to different conditions.

I implemented this method in VBA for MS Excel, and although I know this isn’t the best tool for this job it’s the one we can afford here. It doesn’t actually matter which programming language you use, as the calculations aren’t that intensive.

3. The results

Below I’m showing some sample results of what we can get with this simple methodology. The graphs are shown for a heat exchanger comprised of stainless steel piping, d_e=38.1mm and d_i=35.1mmL=5.0m and 16 tubes in total. It’s arranged in 4 circuits of 4 passes each, and the shell-side fluid is ammonia evaporating at -5ºC. The tube side fluid is water, which is pumped at a rate of 70 m³/h, entering at +2ºC.

From the graphs below, we can see that the ice builds up to a maximum thickness of 0.66mm at the outlet, which stabilizes after ~200 s. This means that in this condition, the heat exchanger will not clog with ice.


Graphs generated by solving the equations presented

So what happens if we reduce the water flow rate to half as much? The results for this condition are presented below. Look at how the ice layer builds up mich thicker now, going towards a thickness of 2.15mm and growing. Also, note in the instantaneous power graph (first graph) how the shell side power (power delivered to the fluid in the shell) is now much different to the tube side power (power delivered to the fluid in the tube). This is exactly due to the ice buildup, which requires heat to be released.

Even in this condition, it is not very likely that the HX will shut off with an ice plug. We can see that the process conditions must be very harsh in order for such a failure to occur.

Graphs for the new condition, with half as much flow rate in the heat exchanger.

So let’s see what happens if, due to a hypothetic process heat load drop, the temperature of the water drops to 1.0ºC instead. The simulation time now was increased to 1000s (~17min) to show the long-term effects. The ice now builds up to 4.5mm thick in the outlet. The entire heat exchanger is now frozen to some thickness. Even so, it seems to be in a stable condition, and if the pump can still move those 35m³/h through the HX, it’ll probably not clog up.

Well, we cannot assure anything about the integrity of the pipe, as the stresses in the pipe wall increase as ice builds up. For this, we need to perform a different calculation.

Yet another graph showing the ice thickness evolution when the water inlet temperature is 1ºC.

As a last graph I wanted to show, take a look at what happens if we simulate this HX, with 70m³/h of water coming in at +2ºC, if we begin with a constant ice thickness of 5mm:

Graphs for the same HX “thawing” with water to the steady-state conditions.

The ice thickness, as one might have expected, decrease to the operational condition. This means that we can simulate both HX freezing and thawing with this method, possibly opening up a way to simulate the behavior of ice banks.

I hope these plots displayed the capabilities of the method and inspired you to make your own version. If you have any doubts, feel free to send me an e-mail or a comment!


Continue to Part III >>