Monday, December 21, 2015

How to Calibrate Water Speed for Vessel Heel/Sensor Angle (Hull Speed Through Water, NMEA 0183 VHW)

After a few rounds of ops testing my multiplexer, I've discovered that my readings are routinely off for the current. To test this, I sailed to known (predicted) current points and compared it to the calculated current from the King Tide Sailing (KTS) system. Here are two such sample points (070@2.40 means 70 degrees true at 2.40 knots):

Predicted Current: 070@2.40 knots. Calculated Current: 070.53@4.65 knots. 

Predicted Current: 060@2.28 knots. Calculated: 073.84@2.77 knots. Also I am .6 feet above the water?
The direction is pretty close (the second measurement is only 10 degrees off), but the first velocity measurement is off by a factor of almost 2. That's unacceptable, and needs to be corrected to satisfy my desire for good instruments.

Paddle Wheel Water/Hull Speed Sensor Calibration

My paddlewheel is part of a triducer that measures depth, temp, and speed (as well as a log, but that doesn't seem to work). However, because I screwed up and didn't fully investigate where to install my new sensor, I could only install it at a 23 degree angle from the vertical underneath a settee. That's really unfortunate, because the depth reading is actually saying the water is deeper than it really is due to the slant distance for the beam, but I actually correct for this offset using the KTS system (keep in mind this is only correct if the water floor is flat). But this also means my paddle wheel is off its desired axis by 23 degrees. Perhaps that means there needs to be another correction for offset angle to derive the actual speed. Which is true, because as the vessel speeds up, the boundary layer against the hull (where the sensor is) actually flows faster and faster than the actual speed of the vessel through the water. Vessel roll/heel has ta similar effect.

I found this handy website that does just that for their own sensor, with the table reproduced below.

Heelº Boat Speed (Knots)
5 10 15 20 25 30
0.0 -2.0 -3.9 -6.0 -7.8 -9.3
10º -0.2 -2.3 -4.0 -6.5 -9.6 -11.0
20º -0.4 -3.9 -6.1 -8.5 -11.5 -13.3

But how do you incorporate this table into the Arduino? It's more than just a linear relationship; at 0 degrees of heel (roll), the sensor error doesn't change at a constant rate. As with anything involving fluid dynamics, most relationships are either trigonometric (cosines and sines) or exponential. Fortunately, after dusting off a few college textbooks, I figured out how to derive a function for use in the Arduino code to correct for boat speed and roll.

Equation to Correct Water Speed for Vessel Speed & Heel/Roll

The key for interpolating between those points in the data table above lies in a Lagrangian polynomial. For those who didn't take any advanced mathematics classes, a Lagrangian polynomial is an equation that accurately describes a curve based on a few given data points. For us, those data points are given in the table (at 10 degrees of heel and at 15 knots, the correction is -4.0 knots). So we use the Lagrangian polynomial to calculate the correction for anywhere other than those data points (say, for 18.56 degrees of roll and 6.23 knots). Deriving the polynomial is left as an exercise to the reader, if you so choose to burn a weekend trying to figure it out. Or you could just keep reading.

For my vessel, I chose to only go up to 10 knots. Why? The theoretical hull speed of my 35 foot sailboat is 6.86 knots (~25' waterline). When would I ever be going faster than 10 knots?

However, the mathematically educated reader would point out that if you came up with a single equation to describe the table, you would actually get a positive correction between 0 and 5 knots, instead of a negative correction. This has to do with the polynomials, and blah blah blah, long story short is that you only use the Lagrangian equation to derive the corrections between 5 and 10 knots. For anything less than 5, it's actually a linear relationship. I know that's in direct contradiction with what I wrote earlier, but looking at the table, the corrections at 5 knots go from 0 to -0.2 to -0.4, and that's it--that's -0.2 knots for every 10 degrees of roll at 5 or less knots. For speeds under 5 knots, it really is a linear relationship.

So, on to the actual equations.

The first step is to figure out the corrections per degree of roll. This is done by doing a quick polynomial for the 5 and 10 knots columns (remember, I'm only going up to 10 knots because I'll never go faster than that through the water). Basically, here's the table I'm trying to fill in:

Obviously, at 0 knots, the correction is 0. But the first thing I need to do is fill in the 10 knot column and the 5 knot column. Then I can fill in horizontally from that. The 5 knot column is easy: it's a linear relationship, so that's solved by the simple equation "y = -.02*roll" and that's it. For the 10 knot column, we actually use a linear relationship between 0 degrees roll and 10 degrees roll because if we used the polynomial, the correction factor goes in the wrong direction. So that linear relationship is "y = -.03*roll - 2" all the way down to 10 degrees roll. So to figure out the polynomial between 10 and 20 degrees roll (and, in fact, beyond that), we come up with the following:

where x is degrees of roll. This produces the following graph:

where you can see that for a little bit, as the x-axis increases (roll), the correction factor goes the wrong way. As roll increases, the correction should always decrease, but it doesn't until about ~3 degrees or so. That's why I use a linear relationship until 10 degrees of roll, and then I pick up that equation, which, thanks to Wolfram Alpha, can be reduced down to this:


where x is roll in degrees. So now, given all this information, we can continue to fill in our table.


Now we need to fill in the horizontal rows, or the actual correction based on vessel speed. We use a similar method as before (linear for 5 and below knots (y=fiveknotcorrection*waterspeed/5), polynomial for more than 5 knots). Using the same equation format, we get something that looks like this (see the EDIT at the bottom for a link to the Wolfram Alpha post):


where x is the speed in knots, a is the 0 knot correction, b is the 5 knot correction, and c is the 10 knot correction. Since the 0 knot correction is always 0, we can eliminate that first term (multiplying anything by zero equals zero). This gives us a simplified equation as follows:


Using our new equation, we can fill in the rest of the rows based on our 5 and 10 knot corrections already calculated. This gives us the following table:


Of course, now that we have actual equations, we can calculate the correction values for up to 90 degrees of roll (which, at 10 knots, gives you a correction factor of -51.50 knots). So now how do we incorporate this into the Arduino code?

The Arduino Code for Vessel Water/Hull Speed Calibration/Correction

The first step is to figure out the 5 and 10 knot correction values for the vessel's roll. Remember, we use a linear equation for 5 knots and 10 knots (below 10 degrees roll), and then a polynomial for 10 knots (above 10 degrees roll). This is a function of roll, mind you. Once we know what our 5 knot and 10 knot corrections are for whatever roll/heel our vessel is presently at, we can plug those into the second formula, which is a function of speed. Remember, this is a linear equation up to 5 knots, and then the polynomial for anything above 5 knots.

 // VHW Vessel Water Speed Sentence //  
    if (strcmp(title,"VWVHW") == 0) {    
     wspworking = 0; // resets the "waterspeed is bad" counter back to 0 since it's good  
     float wsp0 = nmeaDecoder.term_decimal(5);  
     float wsp1 = wsp0;  
     if (wsp0 > 10){ // just in case the speed is above 10, we reset it back to 10  
      wsp1 = 10;     // so as not to get any weird correction values  
     }  
     float pwa = abs(tda - rll); // corrects for the triducer angle offset    
     float fkc = -.02*pwa;       // 5-knot correction for vessel roll  
     float tkc = -.03*pwa-2;     // 10-knot correction for vessel roll  
     if (pwa > 10){  
      tkc = pwa*(.035-.0065*pwa)-2;  
     }  
     float wspcor = fkc*wsp1/5;  
     if (wsp1>5){  
      wspcor = (wsp1*(tkc*(wsp1-5)-2*fkc*(wsp1-10)))/50;  
     }  
     wsp = wsp0+wspcor;  
     if (outputvhw == true) {    
      char vhwSentence [49];    
      byte csw;    
      PString strw(vhwSentence, sizeof(vhwSentence));    
      strw.print("$VWVHW,");    
      strw.print(hdt);    
      strw.print(",T,");    
      strw.print(hdm);    
      strw.print(",M,");    
      strw.print(wsp);    
      strw.print(",N,");    
      strw.print(wsp*1.852);    
      strw.print(",K*");    
      csw = checksum(vhwSentence);    
      if (csw < 0x10) strw.print('0');    
      strw.print(csw, HEX);    
      Serial.println(vhwSentence);    
      }    
    }  

where wsp0 is the raw waterspeed data, tda is triducer angle offset from the verticle (mine is slanted 23 degrees to starboard--this is necessary because at 23 degrees roll to starboard, the actual roll experience on the paddle wheel is 0), rll is vessel roll, and wspcor is the correction factor (always negative). If, for some reason, my boat is surfing along at faster than 10 knots, I included a quick little code (for wsp1) to only derive corrections as if the speed were only 10 knots. This is because we didn't include any data points from

Now all that's left is to test it out and see how it works. Hopefully, the new code will also include a 10 hz refresh rate for headings/attitude information, but I have to wait until the New Year to try it out. Expect a full write up and a new KTS Source Code then.

EDIT: I was curious to see the differences if I did the polynomial out to 15 knots instead of ending at 10. Here's the new equation for step one (calculating the 15 knot corrections where x is the roll):


which further simplifies to:



Now that we know what our 0, 5, 10, and 15 knot corrections are based on our roll, we can input that info to the following polynomial, which is a function based on speed (here is the Wolfram Alpha link if you want to play around with it):


which further simplifies to:


where b is the 5 knot correction, c is the 10 knot, and d is the 15 knot correction. As you can see, this can get messy really quick. When I plugged it into my table, the results got... weird, to say the least. I checked a few points, and basically sought to minimize the jump from 5 to 6 knots, since that's where it goes from linear to polynomial. If we include the 15 knot corrections, the jump increases by quite a bit more than just the 5/10 knot corrections--I don't like this, so I'm going to keep the method described originally (without a 15 knot correction factor).