Моделирование процесса измерения уровня жидкости акустическим методом, учитывающим факторы среды Текст научной статьи по специальности «Математика»

Электронный журнал «Техническая акустика» http://www .ejta.org 2007, 10

А. Н. Громов

ОАО «НИИТеплоприбор», Москва, Россия e-mail: an_gromov@rambler.ru

Моделирование процесса измерения уровня жидкости акустическим методом, учитывающим факторы среды

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

ВВЕДЕНИЕ

В промышленности распространенной задачей является измерение уровня веществ в емкостях. По методам измерения уровня имеется большое число публикаций [1-12], и в силу актуальности задачи их число продолжает расти.

Наиболее широко применяют приборы бесконтактного измерения уровня радиолокационного и ультразвукового типов.

Принцип действия радарного уровнемера основан на радиолокационном методе измерения дальности объекта наблюдения, т.е. на облучении контролируемой поверхности радиоволнами и приеме отраженного радиосигнала.

Наиболее естественным представляется импульсный метод, заключающийся в измерении времени запаздывания принятого импульса относительно излученного. Однако, в силу “огромной” скорости радиосигнала, очевидны технические трудности реализации этого метода. Импульсные методы применяются лишь там, где отражательная способность продукта сравнительно высока, диапазон измеряемых уровней не велик и не требуется высокая точность.

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

Получена 21.01.2007, опубликована 04.05.2007

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

В статье [12], описывается история появления радарных уровнемеров. Рассмотрен принцип их действия с оценкой вклада отдельных узлов уровнемеров в формирования технических и стоимостных параметров.

На российском рынке наиболее известны иностранные производители радарных уровнемеров. Лидирует отделение шведского концерна Saab - Saab Tank Control, которое поставляет уровнемеры для измерения уровня жидкостей в резервуарах хранилищ нефти, сжиженных газов и продуктов нефтепереработки. Продукция этой фирмы отличается высокими метрологическими характеристиками, сохраняющими стабильность в течение всего срока службы без поверок и перекалибровок. Стоимость приборов такого типа превышает 10 тыс. долл. По этой причине они используются преимущественно в экономически благополучной нефтехимической отрасли. Именно поэтому активно развиваются различные ультразвуковые методы измерения уровня и уровня жидкости в резервуарах в частности [4-7].

Одним из направлений является метод акустической локации, который использует для определения скорости звука c в воздушной среде над жидкостью эталонное плечо [10-11]. Так как скорость акустических волн на шесть порядков меньше скорости радиоволн, то в акустических локаторах, в отличие от радарных уровнемеров, измеряется время запаздывания принятого импульса относительно излученного, и уровень вычисляется по формуле

L =1 c ■ t,

2

где L — уровень, c = Const — скорость распространения волны, t — время прохождения сигнала.

Применение эталонного плеча, при условии c = Const вдоль оси емкости, позволяет свести погрешности измерения к погрешности определения длины эталонного плеча. Однако на практике условие c = Const зачастую не выполняется. Это происходит, например, при использовании данного метода для измерения уровня тяжелых нефтепродуктов. В работе [3] установлено, что погрешность измерения уровня акустическим эхо-методом, из-за влияния среды на скорость распространения акустической волны, может составить порядка 5%. Измерения проводились в технологическом резервуаре с мазутом вместимостью 1000 тонн и высотой 10.7 метра. Показано, что на 80% эта погрешность обусловлена наличием градиента температуры по высоте емкости. Полное изменение температуры может достигать 40-50° С.

В работе [1] предложен метод измерения уровня, учитывающий физические свойства среды на пути распространения акустического сигнала и таким образом свободный от основного недостатка традиционных локационных способов измерения

уровня. Особенность метода в том, что задача измерения уровня L поставлена как задача математическая. А именно, как задача определения длины отрезка [0, L], на котором порождена краевая задача для одномерного уравнения, описывающего распространение акустических колебаний в плоскослоистой среде, с заранее не известной зависимостью скорости звука от вертикальной координаты.

Таким образом, для вертикальной составляющей давления (радиальное волновое число полагается равным нулю) имеем одномерное уравнение Гельмгольца

-P" = k2 P, (1)

а

где k = —^- — волновое число, а — частота. c(z)

Уравнение (1) дополняется граничными условиями импедансного вида

P (0 ) + а0 P (0) = 0, P (L ) + aLP (L ) = 0 (2)

где a0,aL — действительные коэффициенты.

Если среда вдоль вертикальной оси OZ емкости однородна и скорость звука c (z) = Const = c0, то краевая задача (1), (2) — частный случай задачи Штурма-Лиувилля

с собственным значением Л = k2 = (a/c0)2. Известно [14], что промежуток [0, L], на

котором порождена краевая задача Штурма-Лиувилля, полностью определяется последовательностью ее собственных значений [Лп }, так как

L = п lim

n ^<Х)

(3)

■дп j

Для гармонической последовательности собственных частот {а n}, когда а0 = aL = 0

вместо предельного равенства (3) имеем соотношение n п c0

L = ■

(4)

где c0 — постоянная скорость звука; ап — гармонические собственные частоты.

Равенство (4) использовано в работе [2].

Для негармонического спектра {an}, когда c(z) = Const = c0, уровень L вычисляем

по формуле (3). Но в общем случае c(z)^ Const формула (3) уже не дает искомый геометрический параметр L , а определяет некоторую связанную с ним величину, которая, как показано в работе [1], позволяет найти уровень L , привлекая теорию решения обратных задач.

Для восстановления неизвестной зависимости c (z), то есть для решения обратной

задачи, недостаточно множества собственных частот {an} краевой задачи (1), (2). Пусть так же наблюдается решение (1) с граничными условиями

^ (0 ) + а0 P (0) = 0, P' (L ) + а^ (L ) = 0, (5)

где а0 — действительный коэффициент и а0 ^а0.

Обозначим через {сЪп} множество собственных частот краевой задачи (1), (5). В работе [1] показано, что восстановление неизвестной зависимости c (z )и определение уровня L сводится к решению обратной классической задачи Штурма-Лиувилля по двум спектрам {сп}, {(Ъп}, п = 1,2,... с последующим интегрированием обыкновенного

дифференциального уравнения. То есть собственные частоты {<эп }, {соп} , п = 1,2,...

составляют исходные данные для рассматриваемого метода измерения уровня.

Для реализации метода, математические основы которого изложены в [1], автором разработаны эффективные алгоритмы численного решения обратной задачи Штурма-Лиувилля, включая алгоритм вычисления нормировочных чисел по двум спектрам и алгоритм решения интегрального уравнения. Апробация метода путем физического моделирования проблематична, так как требует создания специальной установки. Численное моделирование на ПК лишено этого недостатка и позволяет оценить эффективность разработанных алгоритмов и возможность их использования во встраиваемых системах [17].

Известно [13], что успех решения обратной задачи Штурма-Лиувилля определяется возможностью вычислить по двум спектрам множество нормировочных чисел (квадраты норм собственных функций) {&п }, п = 1,2,.... Это объемная задача, которая требует отдельного рассмотрения. В данной работе нормировочные числа вычисляются по определению.

1. ИСХОДНЫЕ ДАННЫЕ ДЛЯ ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ

Для численного моделирования выбран экспоненциальный профиль скорости звука

а

^) = c0eL , а — безразмерный параметр, который относится к числу “решаемых” [15]. Решение краевой задачи (1), (2) выражается через функции Бесселя. Уравнения для собственных частот так же имеют вид комбинации этих функций.

Наличие аналитических зависимостей упрощает вычисление последовательности нормировочных чисел {осп } и собственных частот {сп }, которые составляют исходные

данные для решения обратной задачи и отыскания параметра L .

Краевая задача (1), (2) введением новой переменной и новой функции приводится к стандартной форме краевой задачи Штурма-Лиувилля. Приведение выполняется по формулам

'=$5 ■ ^=’ (6)

где £ — новая независимая переменная; Y (£) — новая функция.

После преобразования имеем

-7" + ~(А)У = Х У, £е[0,&],

(7)

гДе 7 (€) = -

а с0

У

2(ь-а со£)

Граничные условия (2) примут вид У'(0)-Л,У(0) = о, У'(А) + НУ (а ) = о,

( 0.5а ^ т~ ( 0.5а ^

)

— потенциал; X = т2 — собственное значение.

(8)

где к =-

а0 +•

ь

С0, Н =

аь +•

ь

С0е

Длина промежутка £,ь ^ L, на котором порождена краевая задача (7), (8), дается формулой, аналогичной (3)

А = Піт

( \ п

К'ЧХп )

(9)

где Хп — собственные значения краевой задачи (7), (8).

Для модельного профиля величину (9), согласно (6), вычисляем точно

л Ь а

А = — Iе Ь‘ж = Р(1-^),

п *

(10)

"0 0

ь

где р =-----------, в = е

-а

а с0

Переходя к переменной х є [0,п], приводим задачу (7), (8) к требуемому виду

- У" + д(х)У = ХУ , х є[0, ж\,

(11)

где 7 (х ) =

п

А

1 -в

--х

— потенциал; X =

/I

(п^2'

А

Функция д(х) восстанавливается в ходе решения обратной задачи Штурма-Лиувилля. Для граничных условий (8) имеем

Y'(0)- h1Y(0) = 0, Y'(п) + Ш(п) = 0, (12)

~ р(1 -в) ~ р(1 -в)

где Иг = Ъх----------------------------------------------------------------------, Н = Н- — действительные коэффициенты; р ив —

п

п

определены (10).

Если а0 = аъ = 0 , то из (8), (10), (12) следует

* =-Ьв Н = 1 -в

1 2п ’ 2пв

(13)

С другой стороны, для экспоненциального профиля скорости звука уравнение (1) сводится к уравнению для цилиндрических функций с индексом ноль

ё г.1 ё 2пл

-Р +------Р + т ■ Р = 0,

(14)

где п = ре — новая переменная.

Общее решение уравнения (14) записывается так

р(п) = А• Jo (°п) + В• 70 (сn), (15)

где J0(юг/) , У0(соц) — функции Бесселя порядка нуль первого и второго рода соответственно.

Решения уравнений (11) и (14) связаны формулой (6). Из формулы (15), получаем общее решение уравнения (11)

У (х) =

(

7X1

п

\\

1 -в

— х\

(

+ ВУ0

4Х (

п

Л

1 -в

— х I

1

п Л 2 ■-х I

.1 -в

(16)

где

4Х = А т.

п

Пусть 7(х,Л) те из решений (16), которые нормированы условиями У (0, Л) = 1 , У'(0,Л) = (17)

где \ — действительный коэффициент (13).

Функции 7(х,Л) при любом Л удовлетворяют первому из граничных условий (12). Поэтому собственные значения краевой задачи (11), (12) являются корнями функции Я(Л)

Я(Л) = Г(ж, Л) + НУ (п,Л), (18)

где Н — действительный коэффициент, определенный формулой (13).

Соответствующая собственная функция краевой задачи равна У(х,Лп), п = 0,1,2.... Из (16), (17), определив коэффициенты А , В, получаем

У (х,Х) = -->// У (4Хъ)■ З0 (Х(Ъ - х))-

З1 (Хъ)■ У0 (Х(Ъ - х)) ■ТЪ

(19)

где J1 (Л , У (Л — функции Бесселя порядка единица первого и второго рода

1 -в 1 -в Из формул (18),(19) следует

Я(Л) = -Пу/аЬЛ • [у (Ль) J1 (Ла)- J1 (VЛь)• у (¡Ла)]. (20)

Формула (20) дает уравнение для отыскания последовательности собственных значений {Лп }П=0, которое очевидно будет таким

Л = 0, у ((ЛЬ)) (Ла) - Jl ((ЛЬ) (4Ла) = 0, (21)

где а, Ь — коэффициенты, определенные в формуле (19).

Решение (19) остается ограниченным при Л = 0 и имеет вид

У (х,0)=[ЬЬХ] 2, <22>

то есть Л = Л = 0 является собственным значением краевой задачи (11),(12).

Согласно [13], нормировочные числа определяются равенством

П

а = |У2 (х,Лп), (23)

0

где У (х, Лп) — собственная функция краевой задачи (11),(12).

Формулы (19), (21), (23), приводят к следующему выражению для нормировочных

чисел:

b

an =~ n 2

1 -I п a%(((>-((>)-J,(О)(а-

n = 1,2,3... (24)

Нормировочное число а0 вычисляем непосредственно, подставив собственную функцию (22) в определение (23)

.1 + ß

(и — x \ux = n

о

О) = b11(b - x)dx = п1 + ß . (25)

2. ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ПРОЦЕССА ИЗМЕРЕНИЯ УРОВНЯ

Здесь приводятся результаты численного моделирования для среды с параметром а = 0.085 и Ь = 10 метрам. Значение параметра а выбрано так, что приращение скорости звука на длине емкости составляет 29.4 м/ с и соответствует приращению

( /-')400

температуры примерно на 50°С. Последовательности собственных частот ^]Яп} и

( ч400

нормировочных чисел \ап} которые при натурном моделировании являются

результатом измерения, здесь найдены по формулам (21), (24), (25). Первые и последние пять вычисленных значений представлены в таблице 1. Для собственных

частот приведены их дробные части у[лП — п .

Таблица 1. Собственные частоты и нормировочные числа

Номер Собственные Нормировочные Номер Собственные Нормировочные

п частоты у]— числа ап п частоты числа ап

0 0 3,01359204925 396 6,93642535е-07 1,57079633207

1 0,27434330е-03 1,57162227253 397 6,91895073е-07 1,57079633205

2 0,13729871е-03 1,57100315690 398 6,90156582е-07 1,57079633202

3 9,15482070е-05 1,57088827968 399 6,88426474е-07 1,57079633199

4 6,86652869е-05 1,57084805589 400 6,86705333е-07 1,57079633196

Сравним множество собственных частот экспоненциального профиля Еп с

гармонической последовательностью собственных частот /п = —- Гц, определяемых

2Ь

формулой (4). Круговые частоты ^1—, приведенные в таблице 1, пересчитываются согласно (11) и переводятся в герцы. Таким образом, Гп = ^—— Гц. Пусть 8п = Гп - /п.

2А

Для номера п = 1 находим 51 = ^ - /1 « 0.7 Гц. На рисунке 1 представлен график этой функции.

Рис. 1. Разность собственных частот экспоненциального профиля и собственных частот трубчатого резонатора с абсолютно жесткими крышками

График показывает, что сравниваемые частоты хорошо различимы, начиная с номеров первой десятки.

Получив исходные данные } , {ап , начинаем моделирование процесса

измерения уровня Ь . Сначала по исходным данным восстанавливаем потенциал д(х), т.е. решаем обратную задачу Штурма-Лиувилля. Затем, согласно [1], интегрируем дифференциальное уравнение

:гг2

=

(26)

п

2£’”£’ — 3^'

4£' 4

где д (£) = ' ± '

\Ч>ь У ^ь У

Интегрирование выполняется до значения £ () = %Ь . Тогда Ь = .

Для восстановления потенциала д(х), следуя [13], решаем интегральное уравнение

к(х,у) + р(х,у) +1 к(х, г)р(г,у) = 0 (о < у < х < п),

о

где Р (х, у) — ядро интегрального уравнения, которое имеет вид

1 1 ад ад

р (х у) = — со^у[я0х сы^у-----------------+ X 0;(и) + ^ И+ (у),

(27)

а0

п

(28)

п=1

п=1

где

В-(и) = - -Псов пи, Б+п {у) = - ~

2ап * п 2ап ^

----собпу .

п

и = х - у, 0 < и < п, V = х + у, 0 < у < 2п.

п

Уравнение (27) при каждом фиксированном х является интегральным уравнением типа Фредгольма (роль неизвестной функции играет функция к(х,у), а х играет роль параметра) и имеет единственное решение. Потенциал q(х) находим по формуле

Ч(х) =2-^-К (х,х), (29)

где К(х, у) — решение интегрального уравнения (27).

Ядро (28), интегрального уравнения (27), является медленно сходящимся рядом. В работе использованы асимптотические формулы для собственных частот и нормировочных чисел, чтобы преобразовать ядро ¥(х, у). Эти асимптотические формулы играют и роль регуляризирующего фактора, сужая множество функций среди которых ищется решение и обеспечивая тем самым устойчивость решения по

отношению к малым изменениям начальных условий.

Известно [13], что для бесконечно дифференцируемых функций q(x) имеют место бесконечные классические асимптотические разложения:

гт к Аз тг $2 $4 /о ач

^П = п + ^ + -3 +..., О, = - + -^ + -4- +..., (30)

п п 2 п п

где л\Хп , ап — собственные частоты и нормировочные числа краевой задачи (11), (12).

Численный анализ показал, что, сохраняя приемлемую точность, формулы (30) можно использовать, начиная с номера п = 1 и ограничиться четырьмя коэффициентами, что обеспечивает достаточную гибкость анализа точности решения (27). Коэффициенты А1,к3,к5,к7 и ё2,$4,ё6,ё8 найдены по методу обобщенной регрессии и их значения приведены в таблице 2.

Таблица 2. Коэффициенты асимптотических разложений

Коэффициенты Значения коэффициентов разложения ^Яу, Коэффициенты Значения коэффициентов разложения 1 / 2ап

к1 —,7468—402—8455— е-04 й— -1.67743442374836е-04

к3 -3,4013—335—54744е-07 й4 4.61927618275191е-007

к5 9,44690448318419е-10 й6 -2.1410365069608е-009

к7 8,58707037946961е-11 й8 1.65673220678243е-011

Известно [13] теоретическое значение коэффициента к1, которое обозначим кт .

кт =— П

(31)

где Ъ.х, Н — коэффициенты граничных условий (13).

Формула (31) позволяет оценить относительную погрешность вычисления коэффициента k1. Расчет дает

h-1.100% * 1.8 ,io-7% . kT

Такая высокая точность служит косвенным подтверждением правильности всех выполненных преобразований. Показано, что величина k1 устойчива по отношению к изменению объема выборки, используемой в расчете.

Далее приводится схема преобразования функций D-(u) (см. формулу 28). Все

результаты переносятся на функции D+ (v) заменой и ^ v. Пусть

k k3 k5 „ ч

А„ = -1+-3+-^..., (32)

n n n

где An — та составляющая собственной частоты (30), которая стремится к нулю, когда n .

Тогда функции cos^J^u, комбинации которых составляют ядро (28), можно записать в виде

cos^Á^u = cos A пи cos nu - sin A пи sin nu. (33)

В силу соотношения (32), функции cos Anu , sin Anu разлагаются по отрицательным степеням номера n

cos A nu = 1 + ¿ ^, sin A nu = ¿ ^n--^, (34)

i=1 n i=1 n

где c2i, s2i-1 — коэффициенты, которые подлежат определению.

Для первых коэффициентов разложений (34) найдено

С2 (и )=- и 2 , с4 (и ) = -k1k3u 2 + -214 k14u 4 ,

(u) = -[ к,к5 + — к32 їм2 +—к,3к3и4----— k,6u6,

i 1 5 2 3 ) 6 1 3 720

s1 (и) = k1u , s3(и) = —3u- 1k3u3, s5(и) = k5и -1 —ik3u3 +—1— —5u5,

6 2 120

s7 (и) = k7и -\ — k12k5 + — k1k32 |u 3 +—k14k3и5-1— k17и7,

7w 7 ^2 1 5 2 1 3) 24 1 3 5040 1

где k1, k3 , k5 , k7 — коэффициенты асимптотического разложения (32)

(35)

Собирая вместе формулы (33), (34) и асимптотику для величины 1 /2ап, получаем

n

n=1 k=1

Z7-.-/ \ г í sin nu , , / cos nu

Dn (u) = 2 M2k-1 (uE-lk-r + M2k(uE

_ x _______________________________ 2k

n=1 n n=1 n

+ rN , (36)

где rN — слагаемые с неучтенными степенями —-; (и) —

многочлены от

n' ' ‘ ' '

переменной и .

Многочлены Mt (и) имеют следующий вид:

M2k-1 (u )=-[ S2k-1 + d 2 S2k-3 + ... + d2k-2 S1 i , M2k (u )= C2k + d2 C2k-2 + ". + d2k , (37)

где d2, d4,... — коэффициенты асимптотического разложения 1/2an.

Г1,, ^ sin nu “ cos nu

Известно [16J, что суммы рядов L ~2k—г, L —2— выражаются через многочлены

_2k+1 _^2k

n=1n n=1 n

Бернулли. Формула (36) перепишется следующим образом

fx-w='^{UMkiuhr-, <38)

где Bk (и) — многочлены Бернулли.

Из формулы (38), ограничиваясь значением k = 7, находим

ад 14

L D-(u )= L eiu + rN, (39)

n=1 i=0

где ei, i = 0,...,14 — линейные комбинации коэффициентов многочленов (37).

В силу громоздкости выражений функции ei не выписаны в явном виде. В таблице 3 даны их числовые значения для рассматриваемой задачи (величины e10,...,e14 опущены в силу их малости).

Таблица 3. Значения функций e¡

Функция Значение функции Функция Значение функции

е0 -0.000378324724177147 Є5 5.29329003909441e-013

e1 0.00017204179135448 Є6 1.42645001113616e-012

Є2 2.60571505254777e-006 Є7 -1.54491479236586e-013

Є3 3.4514805373702e-008 Є8 6.98866109936978e-015

Є4 4.58803525031662e-010 Є9 -2.26764572914655e-019

Объединяя формулу (39) с аналогичной формулой для функции D+(v) и возвращаясь к переменным х, y, после несложных преобразований получим

8

F(y)=ZФ\(х)/i(y) + r, (40)

i= 1

где r — остаток ряда, содержащий неучтенные члены асимптотических разложений (30).

Функции <Pi (х) , y/i (х) имеют следующий вид

«+1 (х)= X2C,2ke¡х1 -2k , /к+i(x) = y2k, к = 0,1,...,7, (41)

i=2k

где СП — биноминальные коэффициенты.

Отбросив в формуле (40) остаток r и используя вырожденность, полученного ядра, сводим решение интегрального уравнения (27) к решению системы линейных алгебраических уравнений 8 порядка, что не требует больших вычислительных ресурсов и может быть реализовано встраиваемой системой. Решение интегрального уравнения находим в виде

K (Х У )=Z fi(х)/ (У X (42)

i=1

где fi (х) — функции, которые являются решениями системы линейных

алгебраических уравнений.

Система уравнений, определяющая функции fi (х), имеет следующий вид:

8

Z (ki(х) + Ski )fi(х) = -Рк (), (43)

i=1

х

где — символ Кронекера; 0кг(х) = ¡Pk (У) /i (y)dy •

0

В силу (41), интегралы от произведений срк (y)/ (y) вычисляются аналитически и дают функции 0ki (х) в виде многочленов от переменной х. Решение системы (43) в совокупности с (41),(42) позволяет вычислить функцию K(х, х), а затем по формуле (29) восстановить потенциал д(х).

Аналитическое выражение KT (х, х) функции K(х, х) находим, интегрируя (29) с потенциалом (11)

С п Y1

1

Кт (х,х) = ~ тП-в~х > (44)

1 -в

3 1 -в

где в — параметр, который определен формулой (10).

Формула (44) дает возможность оценить точность численного решения интегрального уравнения. На рис. 2 изображен график функции

= К(х,х)-Кт (х,х) _

Кт (х, х)

относительной погрешности вычисления функции К (х, х).

Рис. 2. Относительная погрешность решения интегрального уравнения (27)

Вычисление дает А(о)^-2.6 х 10-7%. Погрешность монотонно растет и

тах А(х) = А(^)« 0.1% . Это позволяет заключить, что полученное здесь представление (40) для ядра интегрального уравнения дает решение с «хорошей" точностью.

Восстановив потенциал д(х) по формуле (29), интегрируем уравнение (26). В таблице 4 приведены значения искомой функции в последних 10 точках.

Таблица 4. Результат интегрирования уравнения (26)

Переменная интегрирования гк Значения ^{гк ) Переменная интегрирования гк Значения ^{гк )

9.894 0.028668361 10.064 0.029140157

9.928 0.028762775 10.098 0.029234434

9.962 0.028857162 10.132 0.029328684

9.996 0.028951521 10.166 0.029422907

10.03 0.029045853 10.2 0.029517102

Из таблицы 4 следует, что £(9.996) = 0.028951521 является ближайшим к величине

А = p (1 — 0) = 2.8963112-10—2 значением и, следовательно, мы должны принять

измеряемый уровень L равным 9.996 метрам. Таким образом, искомая величина L найдена с абсолютной погрешностью 4 миллиметра (относительная погрешность

0.04.%).

ЗАКЛЮЧЕНИЕ

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

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

Полученные результаты позволяют сделать вывод, что сведение решения обратной задачи Штурма-Лиувилля к решению алгебраической системы линейных уравнений “невысокой” размерности обеспечивает достаточную точность предложенного метода измерения уровня и возможность его реализации в новом поколении интеллектуальных приборов [17].

ЛИТЕРАТУРА

1. Громов А. Н. Метод измерения уровня, учитывающий физические свойства среды на пути распространения акустического сигнала. Сборник трудов XI сессии Российского Акустического Общества, 2001, т. 2, с. 262-265.

2. Бардышев В. И., Громов Ю. И., Громов А. Н., Овчаренко А.Т., Римский-Корсаков А. В. Резонансный акустический уровнемер. Акуст. журн., 2000, т. 46, №3,

с. 320-324.

3. Римлянд В. И., Кондратьев А. И., Калинов Г. А. О точности измерения уровня жидкости в резервуарах акустическим эхо-методом. Акуст. журн., 2001, т. 47, №4, с. 564-566.

4. Клюев М. С., Клюев С. П., Краснобородько В. В. О погрешностях акустического измерения уровня жидкости и методах их снижения. Акуст. журн., 1999, т. 45, №6, с. 825-831.

5. Кабатчиков В. А. Ультразвуковой уровнемер. Патент РФ на изобретение №2064666, 1996.

6. Nyce D. S., Togneri M. G., Bulkowski R. S. Magnetostrictive position sensing probe with waveguide referenced to tip for determining fluid level in a container. US Patent №5848549, 1998.

7. Римлянд В. И., Калинов Г. А. Акустический тракт автоматизированной системы измерения уровня жидкости в резервуарах. Сборник трудов XI сессии российского акустического общества, 2001, т. 2, с. 265-274.

8. Викторов В. А. Резонансный метод измерения уровня, М.: Энергия, 1969.

9. Бобровников Г. Н., Катков А. Г. Методы измерения уровня, М.: Машиностроение, 1977.

10. Хамидуллин В. К. Ультразвуковые контрольно-измерительные устройства и системы. Л.: Изд-во Ленинградского ун-та, 1989, 249 с.

11. Stapleton et al. Ultrasonic liquid measuring device for use in storage tanks containing liquids having a non-uniform vapor density. US Patent №5085077, 1992.

12. Либерман В. В., Личков Г. Г. Радарные уровнемеры. Прошлое, настоящее, будущее. Промышленные АСУ и контроллеры, 2006, №8.

13. Левитан Б. М. Обратные задачи Штурма-Лиувилля. М.: Наука, 1984.

14. Романов В. Г. Обратные задачи математической физики. М.: Наука, 1984.

15. Бреховских Л. М., Годин О. А. Акустика слоистых сред. М.: Наука, 1989.

16. Прудников А. П., Брычков Ю. А., Маричев О. И. Интегралы и ряды. М.: Наука,

1981.

17. Аристов П. А., Громов А. Н., Курносов Н. М., Скрыпников С. Н., Хасиков В. В. Аппаратно-программное обеспечение нового поколения интеллектуальных полевых приборов технологии Fieldbus Foundation. Датчики и системы, 2001, №11, с. 10-14.