Ультразвуковой анемометр на двух HC-SR04

Ранее был сделан прототип анемометра из одного ультразвукового дальномера HC-SR04. Он умел рассчитывать проекцию скорости ветра на линию между приемником и передатчиком. Для получения вектора скорости ветра на плоскости (2D) требуется вторая координата, которую мы получим, если добавим второй датчик перпендикулярно первому. В этом случае можно закрепить анемометр стационарно — отпадает необходимость использовать флюгер и как-то организовывать подвижные контакты.

Первая версия

Сказано — сделано, причем основательно.

Из обрезков полипропиленовых труб сварил крестовину. Все датчики отпаял и удлинил проводами, которые проложил внутри труб. Расстояние между датчиками получилось 70 см. 

Код программы такой.

#include <dht.h>
#include <OneWire.h>
#include <DallasTemperature.h>
#define Trig 4
#define Echo 2
#define Trig2 8
#define Echo2 12
#define ONE_WIRE_BUS 7
#define Steps

dht DHT;

#define DHT21_PIN 0
static const float defDist  = .6985; // m
static const float defDist2 = .713; // m

// Setup a oneWire instance to communicate with any OneWire devices (not just Maxim/Dallas temperature ICs)
OneWire oneWire(ONE_WIRE_BUS);

// Pass our oneWire reference to Dallas Temperature. 
DallasTemperature sensors(&oneWire);
void setup() 
{ 
pinMode(Trig, OUTPUT);
pinMode(Echo, INPUT);
pinMode(Trig2, OUTPUT);
pinMode(Echo2, INPUT);
Serial.begin(57600); 
  // Start up the library
  sensors.begin();
Serial.println("X Distance  Tds18820  Tcalc Tdht  Hum V");
} 
unsigned long impulseTime=0; 

void loop() 
{
      // READ DATA
    //Serial.print("DHT21, \t");
    int chk = DHT.read21(DHT21_PIN);
    float DHTtemp = 10; 
    float DHThum = 50;
    switch (chk)
    {
    case DHTLIB_OK:
//        Serial.print("OK,\t");
      DHTtemp =DHT.temperature; 
      DHThum = DHT.humidity;
      break;
    default:
        Serial.print("DHT Error,\t");
        break;
    }

    // DISPLAY DATA
//    Serial.print(DHThum, 1);
//    Serial.print(",\t");
    //Serial.println(DHTtemp, 1);

    sensors.requestTemperatures(); // Send the command to get temperatures DS18820

  float dist = 0;
  float dist2 = 0;
  float temp = sensors.getTempCByIndex(0); //DHTtemp;

unsigned long impulseTime=0; 
unsigned long impulseTime2=0; 
int N=250;
for (int i = 0; i <N; i++) { 
digitalWrite(Trig, HIGH); delayMicroseconds(10); digitalWrite(Trig, LOW); 
impulseTime +=pulseIn(Echo, HIGH); delay(50); 
digitalWrite(Trig2, HIGH); delayMicroseconds(10); digitalWrite(Trig2, LOW); 
impulseTime2 +=pulseIn(Echo2, HIGH); delay(50); 
} 
//float P = 101325; 
float P = 761 * 133.3; 
float M = (28.95-10.934*DHT.humidity/100*(133.3*4.579*exp(17.14*temp/(235.3+temp)))/P)/1000; 
float R= 8.31447; float X = 1.4 * R/M ; //X = 287; 
float c = sqrt( X *(temp+273.15)); 
dist = impulseTime * c / 1e6 /N; 
dist2 = impulseTime2 * c / 1e6 /N; 
float Speed_of_sound = defDist*N/impulseTime * 1e6; 
float Speed_of_sound2 = defDist2*N/impulseTime2 * 1e6; 
float Tcalc = sq(Speed_of_sound)/X - 273.15; 
float v = c- Speed_of_sound; 
float v2 = c- Speed_of_sound2; float v3 = sqrt(sq(v) + sq(v2)); 
int wd = int(atan(-v/v2)*180/3.1416); 
if (v>0) {wd+=90;} else {wd+=270;}

//Serial.println("X Distance  Tds18820  Tcalc Tdht  Hum V");
//Serial.println(String(impulseTime) + char(9) + String(impulseTime2));
Serial.println(String(impulseTime) + char(9) + String(impulseTime2) + char(9) + String(dist, 5) + char(9) + String(dist2, 5) + char(9) + String(temp) + char(9) + String(Tcalc) + char(9) + String(DHTtemp)+ char(9) + String(DHThum) + char(9) + String(v) + char(9) + String(v2) + char(9) + String(v3) + char(9) + String(wd));

}

Два последних числа дают искомую горизонтальную скорость и направление ветра. Направление рассчитывается в виде азимута к направлению на север и дается в градусах. Вращение по часовой стрелке.

Увы, результаты меня разочаровали.

При усреднении в 25 измерений, показания в спокойном воздухе прыгают в среднем до 1.5 м/с, при этом измерения выдаются примерно раз в сек. Если усреднить в 10 раз больше показаний ситуация улучшается, но кардинально проблему не решает. К тому же судя по графику скоростей в двух осях, одна пара датчиков фонит существенно больше другой. Скорее всего дело в проводах, которыми я удлинил датчики. Придется переделывать.

Вторая версия

Есть еще одна причина все переделать. Как отмечалось в первой теоретической части, скорость звука изменится на 1 м/с при изменении температуры примерно на 1.5 °С. Погрешности измерений по обоим осям складываются. Нужно понимать, что порывы теплого или холодного воздуха могут существенно исказить показания такого анемометра. Нет смысла в показаниях 4 м/с при легком дуновении теплого ветерка. Из диаграммы натурного эксперимента видно, что даже медленное изменение температуры вызывает дрейф измеренной скорости, а быстрое изменение температуры на 1 градус скачком поменяло измеренную скорость ветра на 1.5 м/с, в то время как датчик температуры медленно отрабатывает это изменение. Важно заметить, что эксперимент этот проходил прямо у меня на столе и изменение температуры было естественным — я ничего не трогал и искусственно ничего не нагревал.

И тут на помощь приходит тот же принцип, что и при измерении расстояния. Если помним, датчики у оригинального HC-SR04 расположены вместе, поэтому результаты не зависят от наличия ветра. Если измерить скорость звука на известном расстоянии сначала в одном направлении, а затем в другом, то разница этих двух показаний, деленная пополам и будет искомой скоростью ветра в проекции на эту ось. При этом, изменение температуры в диапазоне ±25°С дает погрешность ±4%, что абсолютно не критично и мы можем обойтись вообще без термометра. Да и зачем нам термометр? Если мы знаем время прохождения сигнала в обоих направлениях, то по формулам из прошлой статьи мы легко вычислим температуру, а значит сможем уточнить скорость ветра.
Есть лишь одна маленькая загвоздка — придется использовать два HC-SR04 на одной оси. В промышленных образцах датчики попеременно выполняют роль приемника и передатчика. В нашем случае для этого придется подключить пищалки напрямую к arduino и программно генерировать 8 импульсов 40 кГц на одной, после чего вычленять их из другой. Зная про определенные сложности на этом пути, мне представляется проще купить еще 2 датчика по 55 рублей и попытаться обойтись малой кровью. Этим я займусь в следующий раз. А пока на двух датчиках сделаю измерение скорости ветра по одной оси и измерение температуры в такой конфигурации. Главная проблема здесь убрать помехи, которые дают такой большой разброс показаний в спокойном воздухе.

Конструкция

Вооружившись паяльником конструкция была беспощадно распаяна на составляющие. Новую версию решил не делать так основательно, а зря. Никогда не угадаешь, где найдешь, где потеряешь. Получилось как-то так.

Во-первых, приемник расположил как можно ближе к плате, а передатчик удалил всего лишь на 20 см. Второй комплект перевернул на 180 градусов и пищалки скрепил попарно изолентой. Чем точнее соблюсти соосность обоих пар датчиков, тем лучше. В идеале мы должны получить абсолютно идентичные показания скорости прохождения сигнала в обоих направлениях в спокойном воздухе. Натурные испытания подтвердили нашу теорию. В такой конфигурации получается мало помех и весьма точные показания независимо от температуры, что подтверждается графиком ниже.

Вначале я пробовал просто дуть по направлению от синей пары к черной. Моих легких явно недостаточно. Но любопытный факт — воздух в легких успел нагреться на 1°, что раньше вызвало бы скачок скорости на 1.5 м/с, т.к. DS18B20 просто ничего не заметил. Отметим, что мои легкие способны дать всего лишь 0.5 м/с. Дальше я включил большой напольный вентилятор и направил все также от синего к черному. Видно как пошел более прохладный воздух из глубины комнаты и даже DS18B20 начал отрабатывать это снижение, но теперь его значения не используются для расчета скорости. Сделал открытие, что мой вентилятор дует со скоростью около 2 м/с. Дальше в течение паузы видим постепенное увеличение температуры и отличную корреляцию между рассчитанной и измеренной температурой. В конце поставил вентилятор с другой стороны и получил 2 м/с в обратном направлении с падением температуры. Ура, товарищи, это работает!

Программа расчета скорости ветра

#include <dht.h>
#include <OneWire.h>
#include <DallasTemperature.h>
#define Trig 4         // HC-SR04 №1
#define Echo 2
#define Trig2 8        // HC-SR04 №2
#define Echo2 12
#define ONE_WIRE_BUS 7 // DS18B20
#define Steps

dht DHT;

#define DHT21_PIN 0    // DHT21
static const float defDist  = .2121; // m
static const float defDist2 = .2121; // m
float Tcalc = 0;

// Setup a oneWire instance to communicate with any OneWire devices (not just Maxim/Dallas temperature ICs)
OneWire oneWire(ONE_WIRE_BUS);

// Pass our oneWire reference to Dallas Temperature. 
DallasTemperature sensors(&oneWire);
void setup() 
{ 
pinMode(Trig, OUTPUT);
pinMode(Echo, INPUT);
pinMode(Trig2, OUTPUT);
pinMode(Echo2, INPUT);
Serial.begin(57600); 
  // Start up the library
  sensors.begin();
Serial.println("X Distance  Tds18820  Tcalc Tdht  Hum V");
} 
unsigned long impulseTime=0; 

void loop() 
{
  float temp = 0;
  float DHTtemp = 0; 
  float DHThum = 50;
  // READ DHT DATA
  int chk = DHT.read21(DHT21_PIN);
  if (chk == DHTLIB_OK) 
    {
      DHTtemp =DHT.temperature; 
      DHThum = DHT.humidity;
    }
  if (sensors.getDeviceCount() > 0)
    {
      sensors.requestTemperatures(); // Send the command to get temperatures DS18820
      temp = sensors.getTempCByIndex(0); //DHTtemp;
    }
  float dist = 0;
  float dist2 = 0;

  unsigned long impulseTime=0; 
  unsigned long impulseTime2=0; 
  int N=50;
  for (int i = 0; i <N; i++)
  {

  digitalWrite(Trig, HIGH); 
  delayMicroseconds(10);
  digitalWrite(Trig, LOW);
  impulseTime +=pulseIn(Echo, HIGH);
  delay(50);
  digitalWrite(Trig2, HIGH); 
  delayMicroseconds(10);
  digitalWrite(Trig2, LOW);
  impulseTime2 +=pulseIn(Echo2, HIGH);
  delay(50);
  
  }

  //float P = 101325;
  float P = 761 * 133.3;
  float M = (28.95-10.934*DHThum/100*(133.3*4.579*exp(17.14*Tcalc/(235.3+Tcalc)))/P)/1000;  //M = 0.02895;
  float R= 8.31447;
  float X = 1.4 * R/M ;     
  Tcalc = sq((defDist+defDist2)*N/(impulseTime+impulseTime2) * 1e6)/X - 273.15;
  float c = sqrt( X *(Tcalc+273.15));
  dist = impulseTime * c / 1e6 /N;
  dist2 = impulseTime2 * c / 1e6 /N;

  float Speed_of_sound = defDist*N/impulseTime * 1e6;
  float Speed_of_sound2 = defDist2*N/impulseTime2 * 1e6;

  float v = (Speed_of_sound-Speed_of_sound2)/2;
  //float v2 = c- Speed_of_sound2;
  //float v3 = sqrt(sq(v) + sq(v2));
  //int wd = int(atan(-v/v2)*180/3.1416);
  //if (v>0) {wd+=90;} else {wd+=270;}

  Serial.println(String(impulseTime) + char(9) + String(impulseTime2) + char(9) + String(dist, 5) + char(9) + String(dist2, 5) + char(9) 
  + String(temp) + char(9) + String(Tcalc) + char(9) + String(DHTtemp)+ char(9) + String(DHThum)+ char(9) + String(M,5) + char(9) + String(v));

}

Программа будет работать и без датчиков DHT-21 и DS18B20. DS18B20 для вычислений в этом коде нигде не задействован — только выводится в терминал как эталон. Без датчика влажности температура будет рассчитываться как для воздуха с 50% влажностью. На практике это вносит очень маленькую погрешность. На измерения скорости ветра эти датчики вообще не оказывают никакого влияния.

Собственно это все что можно выжать из двух HC-SR04. Для получения вектора скорости ветра на плоскости нужно добавить еще 2 датчика перпендикулярно первым и по формулам первой версии получить полную скорость и направление. Этим займусь как только приедут заказанные дополнительные датчики.

P.S.

Датчики давно приехали, конструкцию переделывал еще 2 раза и в конце концов он заработал как надо, но до крыши этот ультразвуковой анемометр так и не доехал, поэтому до сих пор продолжения и не написал, хотя идея рабочая.

P.P.S. 2018

По многочисленным просьбам выкладываю итоговый скетч, который не требует никаких библиотек (кроме стандартной EEPROM) и работает с 4 датчиками. Код со всякими вкусностями типа встроенной калибровки и сохранением калибровочных значений в энергонезависимую память. И самое главное. Описанная выше проблема с погрешностями по одной из осей была связана не с проводами, а с работающими в одной комнате с датчиками импульсными блоками питания компьютера, монитора и т.п. (их схема преобразования работает на близкой частоте 40 кГц). Я остановился на проблеме выноса датчика на улицу подальше от помех (с передачей данных по блютус). В остальном это работает. Это версия для распаянных датчиков, но есть способ не распаивать. Если вернусь к проекту — реализую.
Для этого кода неважно какое расстояние между датчиками. Нужно поместить устройство в безветренное пространство (и без импульсных помех) и через терминал несколько раз отдать 2 команды:

t21.5
u

Первая — текущая температура по эталонному термометру (любой уличный), вторая — говорит контроллеру что сейчас скорость ветра 0. Согласно этим данным он вычислит расстояние между датчиками и запишет их в EEPROM. Все дальнейшие измерения будут отталкиваться от этих значений.

// WindSpeed v.4 - Анемометр
// Copyright  Evgeny Istomin gena@regimov.net, blog.regimov.net
// Используются 4 датчика HC-SR04. Приемник и передатчик разнесены на противоположные концы крестовины
// На схеме изображено положение приемников для правильного расчета направления и силы ветра
//                    HC-SR04 №1
//                    Север (0 гр)
//                         o
//                         |
//                         |
//   HC-SR04 №2    o-------|---------o  HC-SR04 №4
//   Запад (270гр)         |           Восток (90 гр)
//                         |
//                         o
//                    HC-SR04 №3
//                    Юг (180 гр)
// при выборе материала крестовины руководствоваться http://temperatures.ru/pages/temperaturnyi_koefficient_lineinogo_rasshireniya
// лучший выбор материала - труба инвар 36H, но и обычное железо вполне годится :-)

#define DEFINE_DISTANCE 0.22  // примерное расстояние между датчиками, измеренное линейкой, в метрах.
#define MES_PAUSE       90     // Пауза между измерениями для затухания отражений. В реальности нужно не менее 1, мс.
#define MES_AVERAGE     8     // сколько измерений усреднять для температуры. 
#define PRINT_PERIOD    500  // период вывода измерений в терминал, ms


#include <EEPROM.h>
#define FALSE 0
#define TRUE  1
#define Echo1 2
#define Echo2 3
#define Echo3 8
#define Echo4 5
#define Trig1 6        // HC-SR04 №1
#define Trig2 7        // HC-SR04 №2
#define Trig3 4        // HC-SR04 №3
#define Trig4 9        // HC-SR04 №4
#define Pow1  10
#define Pow2  11
#define Pow3  12
#define Pow4  13

#define T_ABS       273.15    // температура абсолютного нуля https://ru.wikipedia.org/wiki/%D0%90%D0%B1%D1%81%D0%BE%D0%BB%D1%8E%D1%82%D0%BD%D1%8B%D0%B9_%D0%BD%D1%83%D0%BB%D1%8C_%D1%82%D0%B5%D0%BC%D0%BF%D0%B5%D1%80%D0%B0%D1%82%D1%83%D1%80%D1%8B
#define PRINT_LOOP  PRINT_PERIOD/(4*(MES_PAUSE)) // сколько полных циклов пропустить перед выводом измерений в терминал

float defDist1 = DEFINE_DISTANCE;
float defDist3 = DEFINE_DISTANCE;
float defDist2 = DEFINE_DISTANCE;
float defDist4 = DEFINE_DISTANCE;
float Tcalc = 0;            // температура воздуха (расчетная)
const float DHThum = 50;    // % влажности
float M = 0.02895;          // молярная масса кг/моль https://ru.wikipedia.org/wiki/%D0%9C%D0%BE%D0%BB%D1%8F%D1%80%D0%BD%D0%B0%D1%8F_%D0%BC%D0%B0%D1%81%D1%81%D0%B0
const float R = 8.31447;    // Универса́льная га́зовая постоя́нная Дж/(моль*К) https://ru.wikipedia.org/wiki/%D0%A3%D0%BD%D0%B8%D0%B2%D0%B5%D1%80%D1%81%D0%B0%D0%BB%D1%8C%D0%BD%D0%B0%D1%8F_%D0%B3%D0%B0%D0%B7%D0%BE%D0%B2%D0%B0%D1%8F_%D0%BF%D0%BE%D1%81%D1%82%D0%BE%D1%8F%D0%BD%D0%BD%D0%B0%D1%8F
const float P = 761 * 133.3;// давление в Па. 101325 на уровне моря
float X = 1.4 * R / M ;
float c = sqrt( X * (Tcalc + T_ABS)); // скорость звука м/с https://ru.wikipedia.org/wiki/%D0%A1%D0%BA%D0%BE%D1%80%D0%BE%D1%81%D1%82%D1%8C_%D0%B7%D0%B2%D1%83%D0%BA%D0%B0
float impulseTime1 = defDist1 / c;
float impulseTime2 = defDist2 / c;
float impulseTime3 = defDist3 / c;
float impulseTime4 = defDist4 / c;
unsigned char count = 0;    // счетчик циклов

///////////////////////////////////////////////////////////////////////////////////////////////////////////////
// простой фильтр-усреднитель
float filterA(float y1, float y)
{
  return  ((MES_AVERAGE - 1) * y1 + y) / MES_AVERAGE;
}
///////////////////////////////////////////////////////////////////////////////////////////////////////////////
// измеряем задержку прохождения звука между датчиками, сек
float measument(unsigned char Trig, unsigned char Echo, unsigned char Pow) {
  float y;
  digitalWrite(Pow, HIGH);
  delay(MES_PAUSE);
  digitalWrite(Trig, HIGH);
  delayMicroseconds(10);
  digitalWrite(Trig, LOW);
  y =  pulseIn(Echo, HIGH);
  if (count > PRINT_LOOP) Serial.print(String(y, 0) + char(9));
  digitalWrite(Pow, LOW);
  return y * 1e-6;
};
///////////////////////////////////////////////////////////////////////////////////////////////////////////////
// запоминаем в flash-памяти расстояния между датчиками
void StoreDefDist()
{
  EEPROM.put(0, defDist1);
  EEPROM.put(1 * sizeof(float), defDist2);
  EEPROM.put(2 * sizeof(float), defDist3);
  EEPROM.put(3 * sizeof(float), defDist4);
}
///////////////////////////////////////////////////////////////////////////////////////////////////////////////
// читаем из flash-памяти расстояния между датчиками
float GetDefDist(int adress)
{
  float dd;
  EEPROM.get(adress, dd);
  if (dd <= 0) dd = DEFINE_DISTANCE;
  return dd; 
} /////////////////////////////////////////////////////////////////////////////////////////////////////////////// 
// расчет скорости звука в зависимости от температуры, давления и влажности 
void GetC(float t) { 
  M = (28.95 - 10.934 * DHThum * 0.01 * (133.3 * 4.579 * exp(17.14 * t / (235.3 + t))) / P) / 1000; 
  X = 1.4 * R / M ; 
  c = sqrt( X * (t + T_ABS)); 
} /////////////////////////////////////////////////////////////////////////////////////////////////////////////// 
void setup() { 
  pinMode(Pow1, OUTPUT); 
  pinMode(Pow2, OUTPUT); 
  pinMode(Pow3, OUTPUT); 
  pinMode(Pow4, OUTPUT); 
  pinMode(Trig1, OUTPUT); 
  pinMode(Trig2, OUTPUT); 
  pinMode(Trig3, OUTPUT); 
  pinMode(Trig4, OUTPUT); 
  pinMode(Echo1, INPUT); 
  pinMode(Echo2, INPUT); 
  pinMode(Echo3, INPUT); 
  pinMode(Echo4, INPUT); 
  digitalWrite(Pow1, HIGH); 
  digitalWrite(Pow4, HIGH); 
  digitalWrite(Pow3, HIGH); 
  digitalWrite(Pow2, HIGH); 
  defDist1 = GetDefDist(0); // читаем из flash-памяти расстояния между датчиками 
  defDist2 = GetDefDist(1 * sizeof(float)); 
  defDist3 = GetDefDist(2 * sizeof(float)); 
  defDist4 = GetDefDist(3 * sizeof(float)); 
  Serial.begin(57600); 
  while (!Serial) { 
  ; // wait for serial port to connect. Needed for native USB port only 
  } 
  Serial.println("impT1\timpT3\timpT2\timpT4\tdist1\tdist3\tdist2\tdist4\tTcalc\tv1\tv2\tWD\tv3 " + String(PRINT_LOOP)); 
} 
/////////////////////////////////////////////////////////////////////////////////////////////////////////////// 
void loop() { 
  // период измерений = 1 / ( 4e-3 * (MES_PAUSE + 1)) 
  impulseTime1 = measument(Trig1, Echo1, Pow1); 
  impulseTime3 = measument(Trig3, Echo3, Pow3); 
  impulseTime2 = measument(Trig2, Echo2, Pow2); 
  impulseTime4 = measument(Trig4, Echo4, Pow4); 
  //if (count > MES_AVERAGE) Serial.print(String(impulseTime1*1e6) + char(9));

  Tcalc = filterA(Tcalc, sq((defDist2 + defDist4 + defDist1 + defDist3) /
                            (impulseTime1 + impulseTime3 + impulseTime2 + impulseTime4) ) / X - T_ABS);
  if ((Tcalc > 70) | (Tcalc < -50)) Tcalc = 0;

  GetC(Tcalc);
  float Speed_of_sound1 = defDist1 / impulseTime1 ;
  float Speed_of_sound2 = defDist2 / impulseTime2 ;
  float Speed_of_sound3 = defDist3 / impulseTime3 ;
  float Speed_of_sound4 = defDist4 / impulseTime4 ;

  float v1 = ((Speed_of_sound3 - Speed_of_sound1) / 2); 
  float v2 = ((Speed_of_sound2 - Speed_of_sound4) / 2); 
  float v3 = sqrt(sq(v1) + sq(v2));
  int wd = int(atan(v2 / v1) * 180 / 3.1416);
  if (v1 < 0) {
    wd += 180;
  } else if (v2 < 0) { 
    wd += 360; 
  } 
  if (count > PRINT_LOOP) {
    Serial.println(
      String(c, 5) + char(9) + String(Tcalc) + char(9) + String(v1) + char(9)
      + String(v2) + char(9) + String(wd)   + char(9) + String(v3)
    );
    count = 0;
  }

  while (Serial.available() > 0) {
    char inCh = Serial.read();
    // установка температуры. формат команды: t21.5
    if (inCh == 't') {
      String a = Serial.readString();
      Tcalc = a.toFloat();
      GetC(Tcalc);
    }
    // u - юстировка (установка на 0). формат команды: u
    if ((inCh == 't') | (inCh == 'u')) {
      defDist1 = impulseTime1 * c ;
      defDist2 = impulseTime2 * c ;
      defDist3 = impulseTime3 * c ;
      defDist4 = impulseTime4 * c ;
      StoreDefDist();
    }
  }
  count++;
}
///////////////////////////////////////////////////////////////////////////////////////////////////////////////

Последнее обновление 31.03.2022 в 12:14