Содержание к диссертации
Введение
1. Обзор методов и методик поиска мест сброса и контроля состава сточных вод . 11
1.1 Использование тепловой инфракрасной аэросъемки. 11
1.2 Использование томографического метода согласованной невзаимности (МСН). 13
1.3 Методы определения концентраций веществ в сточных водах. 17
1.4 Выводы. 19
2. Разработка и исследование математических моделей для оценки объемов коммунальных, промышленных и сельскохозяйственных сбросов . 20
2.1 Разработка математической модели для расчета объема выброса сточных вод. 20
2.2 Разработка модели распространения примесей в водоемах . 55
2.3 Выводы. 73
3. Разработка аппаратного обеспечения . 74
3.1 Разработка структуры ИИС. 74
3.2 Разработка модулей ИИС.
3.2.1 Разработка модуля измерения температур с использованием терморезисторов. 77
3.2.2 Разработка модуля измерения температур с использованием полупроводниковых датчиков . 82
3.2.3 Разработка модуля стабилизации положения измерительного шлейфа. 84
3.2.4 Использование сети MicroLAN для построения ИИС. 89
3.2.5 Разработка модуля позиционирования. 107
3.3 Описание результатов эксперимента. 110
3.4 Выводы. Ill
4 Анализ метрологических характеристик разрабатываемой ИИС. 112
4.1 Причины возникновения погрешностей. 112
4.2 Термоанемометрическая погрешность определения температуры воды . 113
4.3 Погрешности, вызванные электронными компонентами схемы. 122
4.4 Выводы к четвертой главе.
5 Заключение. 133
6 Список литературы. 135
7 Приложение А. Перечень методик определения компонентов в природных и сточных водах. 141
8 Приложение Б. Пример измерений температуры. 150
9 Приложение В. Фрагмент программы для работы с датчиками температуры. 154
10 Акт внедрения информационно- измерительной системы контроля загрязнений сточных вод.
- Использование томографического метода согласованной невзаимности (МСН).
- Разработка модели распространения примесей в водоемах
- Разработка модуля измерения температур с использованием полупроводниковых датчиков
- Термоанемометрическая погрешность определения температуры воды
Использование томографического метода согласованной невзаимности (МСН).
Акустическая томография широко признана в качестве перспективного метода мониторинга течений и поля температуры в источнике. Томография дает возможность зондировать районы, недоступные для контактных измерений. Объем получаемой информации возрастает квадратично с ростом числа томографических инструментов в отличие от линейного роста в случае локальных измерений.
Идею метода согласованного поля можно сформулировать как решение обратной задачи путем сведения ее к набору прямых задач. Определение неизвестных физических параметров среды распространения по методу согласованного поля осуществляется в два этапа. На первом этапе среда распространения описывается с помощью некоторого набора параметров и на ЭВМ моделируется распространение звука для выборки возможных значений этих параметров. На втором этапе результаты численного моделирования сопоставляются с экспериментальными данными. Метод согласованного поля построен на допущении, что чем ближе расчетные и экспериментальные звуковые поля, тем ближе модель среды, использованная в соответствующем расчете, к реальной среде распространения.
Для сопоставления полей вводится ценовая функция, которая служит количественной мерой их подобия. Ценовая функция зависит от экспериментальных данных и (неявно) от параметров, задающих пробную среду распространения. Варьируя эти параметры, отыскивается модель водохранилища, для которой ценовая функция достигает глобального минимума. Такая модель и представляет собой решение обратной задачи методом согласованного поля. Примером ценовой функции может служить среднеквадратичное отклонение рассчитанной акустической величины от ее значения в эксперименте. В подводной акустике поле часто измеряется как функция глубины, и, соответственно, при вычислении ценовой функции осреднение ведется по глубине. Сопоставление полей на одной глубине не дает информации, достаточной для отбора одной из множества моделей водохранилища. Инверсия также невозможна, если глубинная зависимость используемой акустической величины вырождена, например, она всегда постоянна или меняется линейно. Для успеха инверсии необходимо, чтобы в общем случае зависимость акустической величины от глубины была достаточно сложной. Сложная пространственно-частотная структура звукового поля в мелком море является недостатком с точки зрения традиционной томографии, поскольку она препятствует разрешению отдельных лучевых приходов. Она же превращается в достоинство при обработке методом согласованного поля и позволяет сделать выбор между близкими моделями водохранилища.
В качестве примера рассмотрим численный эксперимент для модели водохранилища, близкой к летним условиям к Саратовскому водохранилищу. Моделировалось распространение монохроматических акустических волн частоты 100 Гц между единичным трансивером и находящейся от него на расстоянии 2 км приблизительно вертикальной антенной трансиверов. В этом и последующих примерах нереалистично низкая частота звука выбрана для ускорения численных расчетов, т.е. является следствием ограниченных вычислительных возможностей, а не условием применимости самого метода. Глубина водохранилища была принята равной 40 м вдоль всей трассы распространения. Единичный трансивер был расположен в 4 м над дном. Рассматривались антенны с различным числом элементов, глубиной их положения и шагом между ними. На рис. 1.2 показаны профили скорости звука и проекции скорости течения на вертикальную плоскость, в которой расположены трансиверы. Скорости звука и течения предполагались независящими от горизонтальных координат. Невзаимность фазы, вычисленная для изображенного на рис. 16 профиля течения, использовалась в качестве "экспериментальных" данных. Пробные поля вычислялись для двухпараметрической модели течений, в которой переменными являлись амплитуды баротропной и бароклинной компонент скорости потока. Предполагалось линейное изменение бароклинной компоненты скорости с глубиной. Профиль скорости течения, который соответствует параметрам модели, минимизирующим ценовую функцию, дает оценку "истинного" течения по методу МСН.
На рис. 1.3 величина, обратная к ценовой функции, изображена для идеального случая, когда шумы и невязки отсутствуют, а антенна состоит из 161 трансивера, расположенного с постоянным шагом между поверхностью и дном водохранилища. В двумерном параметрическом пространстве эта величина имеет ярко выраженный максимум (т.е. ценовая функция достигает минимума), положение которого в точности соответствует истинному значению скорости течения. На рис. 1.36 показан суммарный эффект, оказываемый: 1) шумами, приводящими к отношению сигнал-шум в 10 дБ; 2) невязкой в расстоянии в 5 м; 3) отклонением антенны от вертикали на 10; 4) систематическими невязками и квазислучайными флуктуациями вертикального профиля скорости звука с амплитудами 0,25 и 0,5 м/с, соответственно; 5) ограничением длины антенны до 150 м и числа ее элементов до 16 (антенна охватывала нижнюю часть волновода от дна до глубины 250 м). Хотя минимум ценовой функции не столь острый, как в идеальном случае, он выражен достаточно четко и позволяет восстановить истинный профиль течения в пределах величины шага сетки в пространстве параметров, на которой вычислялись пробные значения невзаимности (в данном примере 0,1 м/с или 5% от скорости течения на поверхности).
Акустическая система, основанная на МСН, может быть использована для мониторинга переноса тепла. Дистанционное зондирование методом согласованной невзаимности полезно также для количественной оценки потоков радиоактивных и химических загрязнений от находящихся на дне источников.
Разработка модели распространения примесей в водоемах
Для получения численного решения системы (2.43) воспользуемся разностной схемой, приведённой в [3]. В этой разностной схеме использованы наиболее простые аппроксимации производных, а для решения системы алгебраических уравнений, получающихся в результате замены производных их разностными аппроксимациями, применяется одна из разновидностей метода итераций, называемый метод последовательной верхней релаксации (ПВР).
Запись алгоритма для всех уравнений системы (2.43) при этом довольно громоздка, поэтому ограничимся пояснениями только для наименее громоздкого третьего уравнения этой системы . Заменяя в этом уравнении производные их разностными аппроксимациями во всех внутренних узлах введённой выше сетки, получим + — + o)i,j = 0 (2.46) hx Ну
Заметим, что численное решение нам пока неизвестно, поэтому выражение (2.46) надо рассматривать как систему алгебраических уравнений для вычисления функции тока во всех внутренних узлах сетки (всего имеется (NX-2)(NY-2) уравнений для всех /=7,2... NX-2 nj=l,2... NY-2).
Запишем эту систему так, как требуется для её решения итерационным методом Зейделя, т.е. каждое уравнение системы должно быть разрешено относительно одного неизвестного hxlxr J r J/ h2 yui ly — [y/i + ij + y/i- l,j) + —у(у/,У + \ + y/i,j- l)+ C0i,j Если при решении этой системы методом Зейделя задать последовательность решения уравнений в порядке возрастания индексов і и у, то алгоритм метода будет выглядеть так s _ hx hy у/ /, J - 2 2 (2.47) + hx2 hy
Верхний индекс s здесь обозначает номер итерации. При этом необходимо обращать внимание на то, как расставлены эти индексы: если, как принято выше, уравнения системы (2.47) решаются в порядке возрастания индексов і и у, то при вычислении \j/jj значения функции тока при меньших индексах (i-1) и (/-1) уже известны на этой s -ой итерации, поэтому они уже используются в выражении (2.47), а значения функции тока при больших индексах і+l иу+1 известны пока только на предыдущей (s - 7)-ой итерации. Значение вихря в (2.47) использовано на s-ой итерации, так как предполагается, что значение вихря на 5-ой итерации уже известно.
Здесь qT, qro и q - параметры релаксации, а многоточиями заменены громоздкие выражения, в которые входят: число Прандтля, число Грасгофа, параметры сетки hx и hy, значения температур, вихря и функции тока в узлах сетки. Необходимо заметить, что при qr=l, qro=l и q l метод ПВР превращается в метод Зейделя.
Система алгебраических уравнений (2.49) представляет собой одну из возможных реализаций метода конечных разностей для решения системы уравнений конвекции (2.43). Всего в системе 3(NX-2)(NY-2) алгебраических уравнений, которые вместе с разностными аппроксимациями ГУ для температуры на верхней и нижней стенках полости (2.28), с ГУ для температуры на левой и правой стенках, с ГУ для функции тока на всех стенках (2.31) и с ГУ для вихря, пример которых для левой стенки (формула (2.37)), позволяют получить численное решение. Разностная аппроксимация ГУ для вихря очевидна на примере формулы (2.37): ( і ъ Л u)soj = (l-q »)a)s-lo,j + qc» --й/"1].,;-— s u (2.50) V 2 hx J Из (2.50) видно, что метод ПВР распространяется и на вычисление ГУ для вихря. Общее количество итераций, необходимое для решения задачи, зависит от ряда факторов: числа узлов сетки NX NY, физических параметров Рг и Gr, геометрического параметра H/L, выбора параметров релаксации qT, qCOj q и заданной итерационной погрешности EPS. Можно заметить, что зависимость от числа узлов сетки очевидна.
Задание итерационной погрешности EPS определяет окончание процесса итераций. В данной программе окончание процесса решения происходит тогда, когда значения неизвестных, полученных на s-ои итерации, отличаются от значений на предыдущей итерации не более чем на заданную величину. Конкретно итерационный процесс заканчивается тогда, когда будет выполнено неравенство EPSmax EPS (это неравенство можно назвать критерием сходимости итераций), где EPSmax определяется как максимальное из трёх чисел
Разности в числителях и средние величины в знаменателях вычисляются по значениям температуры, функции тока и вихря во всех узлах сетки. Для решения конкретных практических задач использовалась величина EPS =0,001.
Если расчёт выполняется для чисел Рэлея Ra 104, то все три параметра релаксации можно принять равными единице (при таком выборе получается приемлемое число итераций).
Если число Рэлея больше 10 , то необходимо подобрать оптимальные значения параметров релаксации по графикам, показанным на рис. (2.3, так как в противном случае может отсутствовать даже сходимость метода итераций. (Данные, представленные на рис. (2.3, следует рассматривать как ориентировочные, так как они получены экспериментально на сетке с NX=NY=20). .
Быстрота сходимости метода итераций зависит и от задания начальных значений искомых неизвестных (от "нулевого приближения"). В качестве начального приближения для температур в программе используется линейное по координате х распределение, а для функции тока и вихря
Необходимо предусмотреть расчёт интегральных чисел Нуссельта на левой (х=0) и правой (х=1) границах полости. Знание этих чисел позволяет оценить влияние конвекции на перенос тепла. Для такой оценки надо выполнить дополнительный расчёт для тех же условий, но с числом Грасгофа, соответствующим отсутствию конвекции (например, Gr=l). Отношение выданных программой чисел Нуссельта при расчёте с конвекцией и при расчёте без неё равно отношению потоков тепла через границы с "включённой" конвекцией и без конвекции (в так называемом случае "чистой теплопроводности").
Разработка модуля измерения температур с использованием полупроводниковых датчиков
Как видно из (3.4), результат измерения величины термосопротивления не зависит от параметров измерительного канала, напряжения питания моста и величины сопротивления R 5. Исключается также влияние нелинейности функции преобразования моста. Точность измерения Rt определяется главным образом точностью образцовых сопротивлений R1-R5.
Для уменьшения погрешностей, связанных с температурным дрейфом элементов описанной схемы, наиболее критичные из них (к изменению температуры) необходимо термостатировать. Это в первую очередь относится к образцовым сопротивлениям. В настоящее время наибольшее распространение получили термостаты, работающие с перегревом, т. е. поддерживающие температуру заведомо выше максимально возможной. На рисунке 3.6 приведен рисунок термостата в разрезе, разработанного для использования в измерительном модуле.
Корпус термостата выполнен из меди и отшлифован, в качестве нагревателя используется нихромовая обмотка, причем для уменьшения теплового сопротивления обмотка выполнена в три провода. Внутри термостата на уголках закреплены две платы с образцовыми мерами и схемой управления.
Наиболее просто можно реализовать старт-стопный режим работы термостата, более сложным является случай использования ШИМ управления термостатом [61], однако при этом можно обеспечить более высокую точность поддержания температуры. Корпус термостата помещен в кожух, выполненный из твердого пенопласта. Потребляемая мощность данной конструкции при нагреве до 50С не превышает 2-3 ватт. Время разогрева составляет порядка 5 минут. Для управления данным термостатом была разработана схема на базе интегральной микросхемы ТМП - 01 фирмы Analog Devices, (приведена на рис. 3.7.)
В этой схеме реализовано старт стоповое управление термостатом [61]. Информация о включении и выключении термостата передается на порт данных микропроцессора и далее, в случае надобности на центральный блок управления ИИС, в качестве коммутатора в нагревателе термостата используется мощный полевой транзистор.
В более простой схеме используются цифровые датчики температуры фирмы Dallas Semiconductor или National Semiconductor на базе однопроводной сети MicroLAN. Диапазон измеряемых температур в такой схеме составляет -40С +125 ОС с шагом ±,0625. Он состоит из ПЗУ с уникальным идентификационным номером, контроллера MicroLAN, температурного датчика и двух регистров для хранения верхнего и нижнего порогов температуры. Недостатком данного решения является большая инерционность системы. Таблица 3.
Датчики температуры и цифровые термометры фирмы National Semiconductor, Dallas Semiconductor, Analog Device Модель Назначение Темпер. диапазон и Погрешность LM35A Precision С Temp. Sensor -55до+150С 4 до 30 -1,-1.5/+1,+1.5С LM35C Precision С Temp. Sensor -40до+П0С 4 до 30 -1.5/+1.5С DS1820 Precision С Temp. Sensor -40до+110С 4 до 12 -0.5/+0.5С LM74 Precision С Temp. Sensor -40 до +125С 2 до 6 ±0,0625 AD420 Precision С Temp. Sensor -40до+125С 3 до 8 -0.5/+0.5С
Как и некоторые другие приборы iButton фирмы Даллас, термометр не содержит внутреннего источника, а использует "паразитное" питание от одно-опроводной шины. Однако при измерении температуры и записи данных в ЭППЗУ ток потребления микросхемы превышает 1 мА, а максимальный ток, который может обеспечить ведущий шины с помощью нагрузочного резистора 1.5... 5 кОм составляет 3,3... 1 мА. Противоречие разрешается или подключением внешнего источника питания или питанием по сети путем замены нагрузочного резистора на низкоомныи открытый ключ, который подает на 1-проводную шину напряжение питания +5 В (рис 3.8) на время преобразования температуры и записи данных в ЭППЗУ. В микросхему встроен детектор используемого режима питания. +5B Ведущий шины К AVdd U L ) +5Вtт DS1820 Однопроводная шина Рис. 3.8 "Паразитное" питание сенсора температуры с использованием ключа на основе полевого транзистора
Для измерения уровня погружения в заглубителе использован датчик с кварцевыми барочувствительными элементами фирмы Hewlett Packard (рисунок 3.9) с основной частотой 5 МГц, состоящий из разделительной мембраны 1, термостатирующего элемента 3, корпуса 4, кварцевого преобразователя 5, заполненного не сжимающей жидкостью 6. Подвод внешнего давления осуществляется через маслонаполненный капилляр, выполненный из нержавеющей стали и отверстие 2.[63]
Крутизна преобразователя давления в частоту составляет 1,5 Гц/пси (1пси = 6,89 10 Па). Порог чувствительности датчика составляет около 6,89 Па, что составляет 10" верхнего предела измеряемых давлений. При разработке электронной схемы датчика за основу была взята дифференциальная схема включения датчика [63], что позволяет значительно снизить погрешности нуля измерительного преобразователя.
Термоанемометрическая погрешность определения температуры воды
В соответствии с этим представляется возможным на основе моделей (4.2) более точно выполнить оценку искомых термометрических погрешностей. Ввиду того, что модель. (4.2) является нелинейной, а водяной поток оказывается в речном бассейне практически всегда сильно флуктуирующим , то для исследования динамических свойств термометров аналитические и приближенные методы вероятностного анализа оказываются малопригодными. Поэтому для получения адекватных оценок вероятностных характеристик рассматриваемых погрешностей следует использовать статистическое моделирование, суть которого заключается в следующем: - с помощью, численного генератора строятся реализации скорости V(t), принадлежащие случайному процессу с заданными вероятностными характеристиками; - численными методами определяются реализации выходного сигнала термометра u (t) при выходном сигнале V (t), конкретных параметрах модели о, So и постоянной температуре Т; - проводится статистическая обработка зависимости u(t) , в результате которой и получают оценки требуемых вероятностных характеристик термометрической погрешности: e=u(t). (4.5) Указанные расчеты выполняются для репрезентативных областей заданного набора вероятностных характеристик скорости потока Rv, параметров модели &о, So, по которым и получают табличные зависимости вида Rs=f(Rv,0o, s0), (4.6) где Rs- заданный набор вероятностных характеристик погрешности є. Последующая аппроксимация выражения (4.5) позволяет перейти от табличных зависимостей к аналитическим выражениям или графикам. При проведении статистического моделирования необходимо определить вероятностные характеристики скорости набегающего потока V(t), которую можно представить в виде V(t)= I W+co (t) I, (4.7) где W- постоянная средняя скорость движения воды относительно термодатчика; ю - случайные флуктуации вдоль направления средней скорости, которые могут принимать как положительные, так и отрицательные значения.
В результате зависимость V(t) также может оказаться отрицательной. Однако термоанемометрический эффект не зависит от знака, поэтому в выражении (4.7) использован знак абсолютной величины. Предположим, что флук 116 туации со являются нормальным случайным процессом с нулевым математическим ожиданием. Тогда их исчерпывающей вероятностной характеристикой является спектральная плотность Sw(co). Следует ожидать, что в наибольшей степени на термоанемометрическую погрешность будет влиять вид спектральной плотности и коэффициент вариации Cw Л, (4.8) где Dw = 2)SW (u))dco - дисперсия флуктуации; о со - круговая частота.
При задании спектральной плотности Sw(co) в данном случае использована двухвершинная аппроксимация, предложенная в работе [71] . Ее низкочастотный максимум занимает область от 0 до і = 2тг(10"2-5-102)град-с") и обусловлен крупномасштабными колебаниями скорости течений (синоптическая изменчивость, спуск воды из водохранилищ, инерционные и внутренние волны). Высокочастотный максимум обусловлен поверхностными волнами. В высокочастотной области функция Sw(a ), при со і со CON описывается выражением:
Таким образом, для определения функции Sw(p) и дисперсии Dw достаточно задать параметры V,Y,Z и с . При помощи параметра с можно изменять вид спектральной плотности Sw(co) при сохранении общей дисперсии Dw, что удобно при исследовании влияния формы функции плотности Sw(o) на характеристики термоанемометрической погрешности.
Выполнение статистического моделирования потребовало задание репрезентативного поля значений параметров модели (1): 0о,Єо, V0 и п. Показатель п задавался с использованием аппроксимации. В качестве значения V выбрано ОДм-с"1 , характерного для равнинных рек. Диапазон номинальных значений о для термометров весьма широк: от сотых долей секунды до десятков секунд. В качестве значений 0О взяты следующие: 0,1 с, 1 с и 5 с, т.е. достаточно малоинерционные термометры.
Основная мощность флуктуации скорости набегающего потока co(t) находится в полосе частот менее 0,5 Гц, которая сравнима с номинальной полосой пропускания термометра при 0О = 1 с. Если же o lc, что характерно для полупроводниковых термометров, то номинальная полоса пропускания термометра будет гораздо больше 0,5 Гц и его динамические свойства не повлияют на реальную термоанемометрическую погрешность. К сожалению, такие термометры обладают худшими статическими характеристиками.
Для определения характерных значений номинального перегрева so используем выражение (4.4). Учитывая, что коэффициент oto может быть выражен через число Нуссельта Nu, а площадь стержневого термометра S=7i-d-l, где 1-длина термометра, после некоторых преобразований, получим из формулы (4.4) выражение для So:
Для реальных стержневых термометров число Reo составляет 200 - 500, в этом случае, число Nuo лежит в интервале от 10 до 15. Поскольку для тех же термометров L = 0,01 до 0,2 м, то коэффициент М в выражении (4.14) расположен в пределах 0,14 - 0,58 Вт-К"1. Для термометра сопротивления в составе разрабатываемой станции коэффициент М=0,15, а для термометра ТСМ-50 он равен М = 0,55.
Рассеиваемая на термометре мощность Р зависит от параметров конкретной измерительной схемы. Как показано в предыдущей главе, необходимая точность достигается только при включении термометра сопротивления в мостовые измерительные схемы. В нашем случае сопротивление датчика лежит в пределах 400-500 Ом, а амплитуда напряжения питающего мост равна 15 вольтам. Изменяя питание моста в границах 5 -15 В, получим значение выделяющейся мощности в пределах от 0,012 до 0,112 Вт. Подставляя найденное значение Р и М в выражение (4.16), определим диапазон характерных значений для металлических термометров сопротивления в интервале от 2 10" до 7 10"2 К. Наиболее часто интервал значений So термометров более узок: от 0,01 до 0,05 К. Поэтому при моделировании задавались три значения Єо = 0,01; 0,05 и 0,1 К (последнее можно рассматривать как предельное значение для реальных термометров).
В условиях эксплуатации один и тот же датчик может функционировать при различных внешних условиях, т.е. при различных значениях вероятностных характеристик набегающего потока.
Чтобы исследовать влияние изменчивости вероятностных характеристик, определим, что С w изменялся в пределах от 0,2 до 4,0. Для каждого значения Sw с помощью параметра с в выражениях (4.10) и (4.11) задавались две спектральные плотности Sw(co), резко отличающиеся друг от друга по внешнему виду. При параметре С = 0, функция Sw(co) имеет один максимум в области волновых колебаний. При С = 1, функция Sw(co), имеет два максимума. Один из вариантов задания функции Sw(w) при С= 0 и С= 1 изображен на рис. 4.1.
