Projects‎ > ‎Arduino Yun‎ > ‎

K-Type thermocouple correction (NIST)

    1. // corrected temperature reading for a K-type thermocouple
    2. // allowing accurate readings over an extended range
    3. // http://forums.adafruit.com/viewtopic.php?f=19&t=32086&p=372992#p372992
    4. // assuming global: Adafruit_MAX31855 thermocouple(CLK, CS, DO);
    5. float correctedCelsius(){
    6.    
    7.    // MAX31855 thermocouple voltage reading in mV
    8.    float thermocoupleVoltage = (thermocouple.readCelsius() - thermocouple.readInternal()) * 0.041276;
    9.    
    10.    // MAX31855 cold junction voltage reading in mV
    11.    float coldJunctionTemperature = thermocouple.readInternal();
    12.    float coldJunctionVoltage = -0.176004136860E-01 +
    13.       0.389212049750E-01  * coldJunctionTemperature +
    14.       0.185587700320E-04  * pow(coldJunctionTemperature, 2.0) +
    15.       -0.994575928740E-07 * pow(coldJunctionTemperature, 3.0) +
    16.       0.318409457190E-09  * pow(coldJunctionTemperature, 4.0) +
    17.       -0.560728448890E-12 * pow(coldJunctionTemperature, 5.0) +
    18.       0.560750590590E-15  * pow(coldJunctionTemperature, 6.0) +
    19.       -0.320207200030E-18 * pow(coldJunctionTemperature, 7.0) +
    20.       0.971511471520E-22  * pow(coldJunctionTemperature, 8.0) +
    21.       -0.121047212750E-25 * pow(coldJunctionTemperature, 9.0) +
    22.       0.118597600000E+00  * exp(-0.118343200000E-03 * 
    23.                            pow((coldJunctionTemperature-0.126968600000E+03), 2.0) 
    24.                         );
    25.                         
    26.                         
    27.    // cold junction voltage + thermocouple voltage         
    28.    float voltageSum = thermocoupleVoltage + coldJunctionVoltage;
    29.    
    30.    // calculate corrected temperature reading based on coefficients for 3 different ranges   
    31.    float b0, b1, b2, b3, b4, b5, b6, b7, b8, b9, b10;
    32.    if(thermocoupleVoltage < 0){
    33.       b0 = 0.0000000E+00;
    34.       b1 = 2.5173462E+01;
    35.       b2 = -1.1662878E+00;
    36.       b3 = -1.0833638E+00;
    37.       b4 = -8.9773540E-01;
    38.       b5 = -3.7342377E-01;
    39.       b6 = -8.6632643E-02;
    40.       b7 = -1.0450598E-02;
    41.       b8 = -5.1920577E-04;
    42.       b9 = 0.0000000E+00;
    43.    }
    44.    
    45.    else if(thermocoupleVoltage < 20.644){
    46.       b0 = 0.000000E+00;
    47.       b1 = 2.508355E+01;
    48.       b2 = 7.860106E-02;
    49.       b3 = -2.503131E-01;
    50.       b4 = 8.315270E-02;
    51.       b5 = -1.228034E-02;
    52.       b6 = 9.804036E-04;
    53.       b7 = -4.413030E-05;
    54.       b8 = 1.057734E-06;
    55.       b9 = -1.052755E-08;
    56.    }
    57.    
    58.    else if(thermocoupleVoltage < 54.886){
    59.       b0 = -1.318058E+02;
    60.       b1 = 4.830222E+01;
    61.       b2 = -1.646031E+00;
    62.       b3 = 5.464731E-02;
    63.       b4 = -9.650715E-04;
    64.       b5 = 8.802193E-06;
    65.       b6 = -3.110810E-08;
    66.       b7 = 0.000000E+00;
    67.       b8 = 0.000000E+00;
    68.       b9 = 0.000000E+00;
    69.    }
    70.    
    71.    else {
    72.       // TODO: handle error - out of range
    73.       return 0;
    74.    }
    75.    
    76.    return b0 + 
    77.       b1 * voltageSum +
    78.       b2 * pow(voltageSum, 2.0) +
    79.       b3 * pow(voltageSum, 3.0) +
    80.       b4 * pow(voltageSum, 4.0) +
    81.       b5 * pow(voltageSum, 5.0) +
    82.       b6 * pow(voltageSum, 6.0) +
    83.       b7 * pow(voltageSum, 7.0) +
    84.       b8 * pow(voltageSum, 8.0) +
    85.       b9 * pow(voltageSum, 9.0);
    86. }
    Alternatively (check both):
    1.     #include <SPI.h>
    2.     #include "Adafruit_MAX31855.h"
    3.  
    4.     #define DO   12
    5.     #define CS   11
    6.     #define CLK  10
    7.     Adafruit_MAX31855 thermocouple(CLK, CS, DO);
    8.  
    9.     void setup() {
    10.        Serial.begin(9600);
    11.        Serial.println("MAX31855 test");
    12.        // wait for MAX chip to stabilize
    13.        delay(500);
    14.     }
    15.     void loop() {
    16.        // Initialize variables.
    17.        int i = 0; // Counter for arrays
    18.        double internalTemp = thermocouple.readInternal(); // Read the internal temperature of the MAX31855.
    19.        double rawTemp = thermocouple.readCelsius(); // Read the temperature of the thermocouple. This temp is compensated for cold junction temperature.
    20.        double thermocoupleVoltage= 0;
    21.        double internalVoltage = 0;
    22.        double correctedTemp = 0;
    23.  
    24.         // Check to make sure thermocouple is working correctly.
    25.        if (isnan(rawTemp)) {
    26.         Serial.println("Something wrong with thermocouple!");
    27.       }
    28.        else {
    29.           // Steps 1 & 2. Subtract cold junction temperature from the raw thermocouple temperature.
    30.           thermocoupleVoltage = (rawTemp - internalTemp)*0.041276;  // C * mv/C = mV
    31.  
    32.           // Step 3. Calculate the cold junction equivalent thermocouple voltage.
    33.  
    34.           if (internalTemp >= 0) { // For positive temperatures use appropriate NIST coefficients
    35.              // Coefficients and equations available from http://srdata.nist.gov/its90/download/type_k.tab
    36.  
    37.              double c[] = {-0.176004136860E-01,  0.389212049750E-01,  0.185587700320E-04, -0.994575928740E-07,  0.318409457190E-09, -0.560728448890E-12,  0.560750590590E-15, -0.320207200030E-18,  0.971511471520E-22, -0.121047212750E-25};
    38.  
    39.              // Count the the number of coefficients. There are 10 coefficients for positive temperatures (plus three exponential coefficients),
    40.              // but there are 11 coefficients for negative temperatures.
    41.              int cLength = sizeof(c) / sizeof(c[0]);
    42.  
    43.              // Exponential coefficients. Only used for positive temperatures.
    44.              double a0 =  0.118597600000E+00;
    45.              double a1 = -0.118343200000E-03;
    46.              double a2 =  0.126968600000E+03;
    47.  
    48.  
    49.              // From NIST: E = sum(i=0 to n) c_i t^i + a0 exp(a1 (t - a2)^2), where E is the thermocouple voltage in mV and t is the temperature in degrees C.
    50.              // In this case, E is the cold junction equivalent thermocouple voltage.
    51.              // Alternative form: C0 + C1*internalTemp + C2*internalTemp^2 + C3*internalTemp^3 + ... + C10*internaltemp^10 + A0*e^(A1*(internalTemp - A2)^2)
    52.              // This loop sums up the c_i t^i components.
    53.              for (i = 0; i < cLength; i++) {
    54.                 internalVoltage += c[i] * pow(internalTemp, i);
    55.              }
    56.                 // This section adds the a0 exp(a1 (t - a2)^2) components.
    57.                 internalVoltage += a0 * exp(a1 * pow((internalTemp - a2), 2));
    58.           }
    59.           else if (internalTemp < 0) { // for negative temperatures
    60.              double c[] = {0.000000000000E+00,  0.394501280250E-01,  0.236223735980E-04, - 0.328589067840E-06, -0.499048287770E-08, -0.675090591730E-10, -0.574103274280E-12, -0.310888728940E-14, -0.104516093650E-16, -0.198892668780E-19, -0.163226974860E-22};
    61.  
    62.              // Count the number of coefficients.
    63.              int cLength = sizeof(c) / sizeof(c[0]);
    64.  
    65.              // Below 0 degrees Celsius, the NIST formula is simpler and has no exponential components: E = sum(i=0 to n) c_i t^i
    66.              for (i = 0; i < cLength; i++) {
    67.                 internalVoltage += c[i] * pow(internalTemp, i) ;
    68.              }
    69.           }
    70.  
    71.           // Step 4. Add the cold junction equivalent thermocouple voltage calculated in step 3 to the thermocouple voltage calculated in step 2.
    72.           double totalVoltage = thermocoupleVoltage + internalVoltage;
    73.  
    74.           // Step 5. Use the result of step 4 and the NIST voltage-to-temperature (inverse) coefficients to calculate the cold junction compensated, linearized temperature value.
    75.           // The equation is in the form correctedTemp = d_0 + d_1*E + d_2*E^2 + ... + d_n*E^n, where E is the totalVoltage in mV and correctedTemp is in degrees C.
    76.           // NIST uses different coefficients for different temperature subranges: (-200 to 0C), (0 to 500C) and (500 to 1372C).
    77.           if (totalVoltage < 0) { // Temperature is between -200 and 0C.
    78.              double d[] = {0.0000000E+00, 2.5173462E+01, -1.1662878E+00, -1.0833638E+00, -8.9773540E-01, -3.7342377E-01, -8.6632643E-02, -1.0450598E-02, -5.1920577E-04, 0.0000000E+00};
    79.  
    80.              int dLength = sizeof(d) / sizeof(d[0]);
    81.              for (i = 0; i < dLength; i++) {
    82.                 correctedTemp += d[i] * pow(totalVoltage, i);
    83.              }
    84.           }
    85.           else if (totalVoltage < 20.644) { // Temperature is between 0C and 500C.
    86.              double d[] = {0.000000E+00, 2.508355E+01, 7.860106E-02, -2.503131E-01, 8.315270E-02, -1.228034E-02, 9.804036E-04, -4.413030E-05, 1.057734E-06, -1.052755E-08};
    87.              int dLength = sizeof(d) / sizeof(d[0]);
    88.              for (i = 0; i < dLength; i++) {
    89.                 correctedTemp += d[i] * pow(totalVoltage, i);
    90.              }
    91.           }
    92.           else if (totalVoltage < 54.886 ) { // Temperature is between 500C and 1372C.
    93.              double d[] = {-1.318058E+02, 4.830222E+01, -1.646031E+00, 5.464731E-02, -9.650715E-04, 8.802193E-06, -3.110810E-08, 0.000000E+00, 0.000000E+00, 0.000000E+00};
    94.              int dLength = sizeof(d) / sizeof(d[0]);
    95.              for (i = 0; i < dLength; i++) {
    96.                 correctedTemp += d[i] * pow(totalVoltage, i);
    97.              }
    98.           } else { // NIST only has data for K-type thermocouples from -200C to +1372C. If the temperature is not in that range, set temp to impossible value.
    99.              // Error handling should be improved.
    100.              Serial.print("Temperature is out of range. This should never happen.");
    101.              correctedTemp = NAN;
    102.           }
    103.  
    104.           Serial.print("Corrected Temp = ");
    105.           Serial.println(correctedTemp, 5);
    106.           Serial.println("");
    107.  
    108.        }
    109.  
    110.        delay(1000);
    111.  
    112.     }




I thought it may be helpful to point out that the first part of the code, inverting and linearizing the cold junction comp, is not strictly necessary. The fixed Seebeck coefficient used in the 31855 for CJC (40.73 uV/C) is chosen to afford less than 0.1C integral linearity error _over the service temperature range of the chip_ . By comparison, the ADC resolution corresponds to 0.25C, and the internal thermometer's absolute error can be ±2C to start with. 

We therefore suffer no loss of accuracy by substituting the following line:
CODE: SELECT ALL | TOGGLE FULL SIZE
    internalVoltage = internalTemp * 0.04073 ; // (mV) MAX38155 comp scale is 40.73 uV/C

for the entirety of Step 3. This saves considerable memory and execution time! 

(The pseudo-floating point exp() and pow() are huge time sucks for a mere UNO).