Set and Drift; or, the Speed and Direction of the Tidal Current (NMEA-0183 VDR)

Saltstraumen, Norway
Tucked away in a small fjord in Finland, the tidal current can reach 22 knots. Following the theoretical maximum hull speed formula, you'd need a boat almost 270 feet long just to make way in that current. In the San Francisco Bay, the speeds routinely exceed 4 knots (now we're looking at a much more realistic 9 foot boat to reach that speed). But when you're beating into the wind with two reefs in and maybe a little too much headsail, that current is going to push you off. It's the same thing in an airplane; if you have a 30 knot direct cross-wind from the right, you are now moving 30 knots to the left, regardless of how fast you're moving forward. Kind of a problem when you're trying to land. It's also kind of a problem when you're trying to sail a course, and your nose is pointed off by 30 degrees. Even the simple knowledge of how strong and where the current is can make or break your sail plan.

The direction and speed of the current, formally titled Set and Drift respectively, are really easy to figure out if you are given the right information. Over the past few blog posts, I've gone over how to get that information, and more importantly--how to manipulate it.

Recommended Gear

NMEA-0183 Instruments and VDR

The assumption is that you've followed all my previous posts (MPU-9150 Setup, DST-800 Setup) and have everything working swimmingly. What about the GPS, though? I followed David's writeup and soldered a wire to the chip where the output comes through, then kept the original usb cable until the very end where I cut it and ended up with a VCC, Ground, and the NMEA wire. There are probably much easier ways to do this, so... have at it, and report back.

Those are the only three instruments you need: Compass, paddlewheel  and GPS. The most difficult part of this is the paddlewheel. If you're boat is big and heavy, you'll need to haul it out and install it while on the hard. Mine was over due for a new paint job, so the timing worked out pretty well. Good time to check everything else too. Anyway, the point is that we need four numbers: heading, waterspeed (which is not the speed of the water; it is the speed of your hull through the water--sometimes I refer to it as hullspeed, but it is technically called waterspeed in the industry), course made good, and speed over ground. In more simple terms, we need the angle and magnitude of two vectors: our boat through the water, and our boat over the ground.

Here are some handy sketches to help visualize what we're trying to figure out, where
  • hdt = true heading
  • hsp = hullspeed, or speed through the water
  • cmg = course made good
  • sog = speed over ground
and the bold lines indicate the actual force vector acting upon the vessel, and each picture represents one of the three possible situations where a current can act upon your vessel (the port/starboard direction doesn't matter, because you can just mirror the image and the values will be identical).

Scenario 1: The current is pushing you back and to starboard

Scenario 2: the current is pushing your forward and to starboard

Scenario 3: same as Scenario 1, but the current is moving faster than your boat is moving 

How to Calculate the NMEA-0183 VDR Set and Drift Sentence

The math isn't that bad, plus I did it all for you. We're just doing simple trigonometry here, which is a little easier to understand by using vectors. If you really want to get down and dirty, I'll leave you to it--it'll take some time. But, given those four variables, the functions are as follows:




The astute reader may observe that you would need more than just those two equations to get the right answer under any conditions. But we'll take care of that in the actual Arduino code when we get there. For now, those two equations are the most basic representation. Set is measured in degrees relative from your heading; drift is measured in the same units as your hull speed and speed over ground (most likely knots). You would also have to operate in radians instead of degrees, in which case the 180 in the set formula becomes Pi.

What's important here is that these two equations are quite simple, and as long as you have the four data inputs, you can figure it out really quickly. And cheaply. This is the only equipment I could find that gives you similar information. And I can't even find a price; comparable items sell for $15,000 to $20,000, so that gives you an idea of what we're talking about here.

So, quite naturally, I am making the code available for free so you don't have to spend $20,000.

The Arduino Sketch to Calculate VDR

 ////////////////////////////////////////////////////////////////////////////  
 //  
 // This file is part of RTIMULib-Arduino  
 //  
 // Copyright (c) 2014-2015, richards-tech  
 //  
 // Permission is hereby granted, free of charge, to any person obtaining a copy of   
 // this software and associated documentation files (the "Software"), to deal in   
 // the Software without restriction, including without limitation the rights to use,   
 // copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the   
 // Software, and to permit persons to whom the Software is furnished to do so,   
 // subject to the following conditions:  
 //  
 // The above copyright notice and this permission notice shall be included in all   
 // copies or substantial portions of the Software.  
 //  
 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,   
 // INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A   
 // PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT   
 // HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION   
 // OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE   
 // SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.  
 //  
 // Further modified by King Tide Sailing (kingtidesailing.blogspot.com)  
 //  
 ////////////////////////////////////////////////////////////////////////////  
   
 #include <Wire.h>  
 #include "I2Cdev.h"  
 #include "RTIMUSettings.h"  
 #include "RTIMU.h"  
 #include "RTFusionRTQF.h"   
 #include "CalLib.h"  
 #include <EEPROM.h>  
 #include <SoftwareSerial.h>  
 #include <PString.h>  
 #include <nmea.h>  
   
 RTIMU *imu;               // the IMU object  
 RTFusionRTQF fusion;          // the fusion object  
 RTIMUSettings settings;         // the settings object  
 NMEA gps(GPRMC);            // GPS data connection to decode RMC sentence types (otherwise the serial output becomes too buggy if you decode all GPS sentences)  
 NMEA nmeaDecoder(ALL);         // NMEA data decoder for all sentence types (for the Transducer, and Anemometer)  
   
 // User defined variables (for your own boat and geographic location):  
 float magd = -14.0;   // magnetic deviation; if it's East, it's negative (14E in San Francisco 2015)   
   
 // Variables preset with recommended values (but can be changed):  
 #define DISPLAY_INTERVAL 500  // rate (in milliseconds) in which MPU9150 results are displayed  
   
 // Global variables:  
 unsigned long lastDisplay;  
 unsigned long lastRate;  
 int sampleCount;  
 float r = M_PI/180.0f; // degrees to radians  
 float d = 180.0f/M_PI; // radians to degrees  
   
 // Variables exchanged between instruments:  
 float ts;  // timestamp of GPS signal in UTC  
 float sog; // speed over ground  
 float cmg; // course made good  
 float hs;  // water speed, or hull speed through the water  
 float hdm; // magnetic heading  
 float hdt; // true heading  
   
 void setup()  
 {  
   int errcode;  
   
   Serial.begin(4800);   // output port to computer  
   Serial1.begin(4800);  // from GPS  
   Serial2.begin(4800);  // from DST-800 transducer  
   Wire.begin();      // from MPU-9150  
   imu = RTIMU::createIMU(&settings);  
   
   // Checking the IMU, if it fails, simply restarting the connection seems to work sometimes  
   if ((errcode = imu->IMUInit()) < 0) {  
     Serial.print("Failed to init IMU: "); Serial.println(errcode);  
   }  
   if (imu->getCalibrationValid())  
     Serial.println("Using compass calibration");  
   else  
     Serial.println("No valid compass calibration data");  
   
   lastDisplay = lastRate = millis();  
   sampleCount = 0;  
   
   // Slerp power controls the fusion and can be between 0 and 1  
   // 0 means that only gyros are used, 1 means that only accels/compass are used  
   // In-between gives the fusion mix. 0.02 recommended.  
     
   fusion.setSlerpPower(0.02);  
     
   // use of sensors in the fusion algorithm can be controlled here  
   // change any of these to false to disable that sensor  
     
   fusion.setGyroEnable(true);  
   fusion.setAccelEnable(true);  
   fusion.setCompassEnable(true);  
 }  
   
 // calculate checksum function (thanks to https://mechinations.wordpress.com)  
 byte checksum(char* str)   
 {  
   byte cs = 0;  
   for (unsigned int n = 1; n < strlen(str) - 1; n++)   
   {  
     cs ^= str[n];  
   }  
   return cs;  
 }  
   
 void loop()  
 {   
    
  // ***** BU-353 GPS ***** //  
   if (Serial1.available()) {  
    if (gps.decode(Serial1.read())) {  // if it's a valid NMEA sentence  
      Serial.println(gps.sentence()); // print the NMEA sentence  
      sog = atof(gps.term(7));    // speed over ground  
      cmg = atoi(gps.term(8));    // course made good  
      ts = atof(gps.term(1));     // timestamp  
     }  
   }  
  // ***** End BU-353 GPS ***** //  
   
  // ***** DST-800 Transducer ***** //  
   if (Serial2.available() > 0 ) {  
    if (nmeaDecoder.decode(Serial2.read())) {  
     char* title = nmeaDecoder.term(0);  
     if (strcmp(title,"VWVHW") == 0) {  
      hs = atof(nmeaDecoder.term(5));  
      }  
     Serial.println(nmeaDecoder.sentence());  
     }  
  }  
  // ***** End DST-800 Transducer ***** //  
   
  // ***** Integrated Navigation (MPU-9150 and derivatives) ***** //  
   unsigned long now = millis();  
   unsigned long delta;  
   int loopCount = 1;  
     
   while (imu->IMURead()) {  
     // this flushes remaining data in case we are falling behind  
     if (++loopCount >= 10)  
       continue;  
     fusion.newIMUData(imu->getGyro(), imu->getAccel(), imu->getCompass(), imu->getTimestamp());  
     sampleCount++;  
     if ((delta = now - lastRate) >= 1000) {      
       sampleCount = 0;  
       lastRate = now;  
     }  
     if ((now - lastDisplay) >= DISPLAY_INTERVAL) {  
       lastDisplay = now;   
       RTVector3 pose = fusion.getFusionPose();   
       RTVector3 gyro = imu->getGyro();   
       float roll = pose.y()*-1*d;   // negative is left roll   
       float pitch = pose.x()*d;    // negative is nose down   
       float yaw = pose.z()*d;     // negative is to the left of 270 magnetic   
       float rot = gyro.z()*d;     // negative is to left   
       hdm = yaw-90;       // 0 yaw = 270 magnetic; converts to mag degrees   
       if (yaw < 90 && yaw >= -179.99) {   
        hdm = yaw+270;   
       }   
       hdt = hdm-magd;      // calculate true heading   
       if (hdt > 360) {   
        hdt = hdt-360;   
       }   
       if (hdt < 0.0) {   
        hdt = hdt+360;   
       }   
       if (sog == 0.0) {        // this eliminates old artifacts in cmg when your speed over ground is zero  
        cmg = hdt;  
       }  
     
       // HDT True Heading //   
       char hdtSentence [23];    
       byte cst;   
       PString strt(hdtSentence, sizeof(hdtSentence));   
       strt.print("$HCHDT,");   
       strt.print(lround(hdt));   
       strt.print(",T*");    
       cst = checksum(hdtSentence);   
       if (cst < 0x10) strt.print('0');   
       strt.print(cst, HEX);   
       Serial.println(hdtSentence);   
         
       // VDR Set and Drift //  
        float set0 = (cmg-hdt)*r;     // gets the difference between boat's heading and cmg in radians  
        if (set0 < 0) {  
         set0 = set0+(2*M_PI);  
        }  
        float drift = sqrt(pow(hs,2)+pow(sog,2)-2*hs*sog*(cos(set0)));  
        float set1 = (M_PI-atan2((sog*sin(set0)),(hs-sog*cos(set0))))*d;  
        if (hs == 0 && sog == 0) {     // if you're not moving, then set = 0  
         set1 = set0;  
        }  
        float set = hdt+set1;  // converting to true instead of relative degrees  
        if (set > 360) {  
         set = set-360;  
        }  
        float setmag = set+magd;  
        if (setmag > 360) {   
         setmag = setmag-360;   
        }   
        if (setmag < 0.0) {   
         setmag = setmag+360;   
        }  
          
        char vdrSentence [43];   
        byte csv;  
        PString strv(vdrSentence, sizeof(vdrSentence));  
        strv.print("$INVDR,");  
        strv.print(lround(set));  
        strv.print(",T,");  
        strv.print(lround(setmag));  
        strv.print(",M,");  
        strv.print(drift);  
        strv.print(",N*");  
        csv = checksum(vdrSentence);  
        if (csv < 0x10) strv.print('0');  
        strv.print(csv, HEX);  
        Serial.println(vdrSentence);  
     }  
    }  
    // ***** End Integrated Navigation ***** //  
   
 }  

The big thing here is the VDR portion at the end, which spits out an NMEA sentence like this:

          1  2  3  4  5  6  7  
          |  |  |  |  |  |  |  
  $--VDR,x.x,T,x.x,M,x.x,N*hh  
 1. Course, True  
 2. T = True  
 3. Course, Magnetic  
 4. M = Magnetic  
 5. Knots  
 6. N = Knots  
 7. Checksum  
   

I messed around with this for a while, and realized that Magnetic degrees is kind of useless. So I changed it to relative degrees, and took my yacht out for some testing. It took a little while to get it totally correct, but...

The Results

Here's a screenshot of OpenCPN, where I found a forecasted current point and posted up to see how close I was.


As you can see, near this point in the San Francisco Bay, the tidally induced current is forecasted to be flowing towards 211 degrees true (its course), at a speed of 2.5 knots. This is a solid forecast using known techniques that produce good results.

So what did the King Tide Multiplexer say the current was at this point? Here's a screen shot shortly after:


It was saying the current was flowing towards 237 degrees at 1.89 knots, a difference of 25 degrees and half a knot. What gives? Well, a few things:

  • The point at which the current was forecasted is only at that particular point; as anyone who sails in the Bay knows, the current can change drastically and quite quickly over a small distance. So while I tried to get as accurate measurements as possible, the point at which I took the above screenshot isn't the same as the forecast point
  • Variations of water level can change how fast a current is moving based on the nearby bathymetry. I do not believe the current models account for the reduced inflow into the Bay due to the drought
  • I hadn't calibrated my heading sensor on my boat; it was only calibrated at my dining room table, which is a lot different than inside my nav desk
  • The heading sensor was pointing to port by 15 degrees (oops)
Here's another sample output, where I added the cross- and head-component of the current:


I edited the code to have the current direction reflect where it's coming from, as opposed to where it's flowing to. So right at that geographic point in the Bay, the current is coming from 041 (212 relative from my boat's heading, or from my port/aft quarter) at 2.57 knots. This also says that the current is pushing me 1.57 knots to starboard, and 2.17 knots forward. I just need a useful way of displaying this information now.

In any case, I will be taking it out for some more testing after I recalibrate the compass, and mount it more securely to the nav desk. And I'll be using a data logger in conjunction with OpenCPN and pass by the point precisely to see how close my math works out.

What to do next? Definitely find a plugin for OpenCPN that plots the current vector just like the CMG and HDT vector as shown. That way, I'll get a quick visual reference of what precisely the current is doing at a particular moment at a particular point.

EDIT: I have become aware that this is neglecting one important factor, called leeway, or the boat's drift due to being blown sideways by the wind. Unfortunately, the only way to figure out your leeway is through ops testing and deriving its value that way based on your heel angle, the relative wind direction, and speed. Maybe in a perfect world, I'd have all these things... but for now, this will do.

Comments

  1. What are the representation of negative speed current, please ?

    ReplyDelete
    Replies
    1. If you're getting a negative current value, then either my equation is wrong, or you may have accidentally made a typo. I have not received a negative value on my boat, but it is possible I typed it wrong above. But in any case, a negative value for the current would mean it is actually heading the opposite direction (for example, a -5 knot current flowing towards 090 degrees would actually be a 5 knot current flowing towards 270 degrees).

      Delete

Post a Comment

Popular Posts