Численное моделирование процесса работы виброакустического датчика вязкости камертонного типа Текст научной статьи по специальности «Физика»

УДК 536.25

Численное моделирование процесса работы виброакустического датчика вязкости камертонного типа

Андрей В.Минаков* Александр А.Дектерев^

Институт инженерной физики и радиоэлектроники, Сибирский федеральный университет, Свободный 79, Красноярск, 660041,

Россия

Андрей А.Гаврилов*

Институт теплофизики СО РАН, Лаврентьева 1, Новосибирск, 630090,

Россия

Валерий Я.Рудяк

Baker Hughes Inc., Novosibirsk Science Center, Кутателадзе 4а, Новосибирск, 630090,

Россия

Получена 18.08.2009, окончательный вариант 25.09.2009, принята к печати 10.10.2009 Проведено детальное изучение 'работы датчика вязкости камертонного типа при помощи численного моделирования. В результате установлена линейная зависимость сдвига резонансной частоты колебаний датчика от плотности среды, в которую он погружен. Обнаружена линейная зависимость величины, обратной к добротности колебаний датчика, от квадратного корня из динамической вязкости жидкости. Получено хорошее количественное совпадение с экспериментальными данными.

Ключевые слова: аэроакустика, виброакустическая вискозиметрия, численное моделирование, метод контрольного объема.

Введение

Виброакустические датчики вязкости (плотности, давления) благодаря своей высокой чувствительности, малому размеру и низкой стоимости нашли широкое применение в различных отраслях человеческой деятельности. На сегодняшний день в мире существует множество типов подобного рода устройств [1]. Наиболее распространенными из них являются так называемые кварцевые TSM (thickness-shear шс^е)-резонаторы. Принцип действия таких устройств весьма прост и основан на зависимости резонансных характеристик датчика, от физических свойств среды, в которую он помещен. Попытки теоретического исследования данного эффекта представлены во многих источниках [1—2]. Как правило, эти исследования

* e-mail: tov-andrey@yandex.ru t e-mail: dekterev@mail.ru ^e-mail: sigma-cfd@torins.ru © Siberian Federal University. All rights reserved

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

В данной работе в полной постановке было проведено численное моделирование работы датчика вязкости камертонного типа (turning fork sensor), который используется для измерения физических свойств нефти.

1. Математическая модель

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

Поведение датчика описывается уравнениями движения упругой среды в пространстве, которое в общем случае можно записать как [3]:

д2 u

Ps др = VT + Ps F, (1)

где u — вектор деформации, ps — плотность упругой среды, F — вектор массовых сил, T — тензор напряжений, компоненты которого связаны с компонентами тензора деформации e-законом Гука.

Компоненты тензора деформации в случае малых деформаций можно определить следующим образом:

— 1 (dUi + дUfe А i к 2 V дх к dxi ) '

В случае изотропной среды закон Гука имеет вид:

e E Л a -И

Ti k — —- £¿ k + ^-SllOi к ,

1 + a \ 1 — 2a J

где a — коэффициент Пуассона, E — модуль упругости (модуль Юнга).

Поведение окружающей датчик жидкости описывается уравнениями Навье-Стокса:

- уравнение неразрывности

| + V-(PV)— 0;

- уравнение закона сохранения импульса

^ + V (pV • V) — —Vp + Vt, (2)

где V — скорость жидкости, p — плотность жидкости, p — давление, ц — молекулярная вязкость, т — тензор вязких напряжений, составляющие которого имеют следующий вид:

т = / dVi + dv¿_ 2 д

ij ^ \ c)xj дхц 3 ij дхк J

В силу малости рассматриваемых нами колебаний (характерная амплитуда колебаний датчика порядка 10-8 м) пренебрегаем конвективным членом V (pV • V) в уравнении (2). Поскольку относительные изменения плотности и давления малы, уравнения гидродинамики

можно линеаризовать [4]. А так как для данной задачи ротор скорости равен нулю везде, кроме узкого поверхностного слоя толщиной 5 = \ — (порядка 10-6 м), то удобнее всего

V Рш

это сделать при помощи потенциала скорости у.

Используя все рассмотренные приближения, в итоге получаем систему уравнений акустики:

д2 У 2 л 4МлдУ /оч

- с Ау =3РА еъ; (3)

V = Уу;

р = -РЖ + ТАу'

где с — скорость звука в жидкой среде, р — относительное давление.

Таким образом, (1), (3) образуют полную систему уравнений рассматриваемой нами задачи. При этом нелинейная взаимосвязь между подсистемами (1) и (3) осуществляется через граничные условия.

1.1. Граничные условия для потенциала скорости жидкости

На поверхности жидкости, граничащей с датчиком, ставится условие равенства нормальной компоненты скорости жидкости и скорости стенки датчика:

ду < ■ \

тп = (и'п)'

где и — скорость деформации датчика, п — вектор нормали поверхности упругого тела. На твердой жесткой стенке ставится условие:

дУ = °-

дп

На открытой границе граничное условие имеет следующий вид:

+ 1 дУ \7 с дЪ

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

1.2. Граничные условия для вектора деформации твердого тела

На границе соприкосновения датчика и стенок камеры, в которой находится жидкость, задается фиксированное смещение (условие жесткого закрепления). На остальных поверхностях датчика задается сила, действующая на датчик со стороны жидкости Е$. Эта сила имеет следующий вид:

= Ер + ¥т

где Ер — сила давления, действующая по нормали к поверхности датчика и равная

Ер = -рп,

Ет — сила вязкости, направленная по касательной к поверхности датчика. Величину этой силы согласно [4] можно записать как

Гт = -^¿^т, (4)

0

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

Гйрр {■ 1 ¿ит4 Е = -у Iм т +

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

2. Численный алгоритм

Для компьютерной реализации математической модели использовался программный комплекс "стПс^". Возможности и численные алгоритмы, заложенные в нем, детально описаны в работе [5]. Здесь же мы лишь отметим некоторые основные моменты, связанные со спецификой рассматриваемой задачи.

Для решения описанных выше уравнений колебаний используется широко известный метод контрольного объема, суть которого заключается в разбиении расчетной области на контрольные объемы и интегрировании исходных уравнений сохранения по каждому контрольному объему для получения конечно разностных соотношений. Свойства метода достаточно подробно описаны в литературе [6-7].

Из всех представленных уравнений уравнение движения упругой среды (1) имеет наиболее общий вид, поэтому рассмотрим его дискретизацию. Остальные уравнения имеют более простую форму и дискретизируются аналогичным образом. Перепишем уравнение движения упругой среды в следующем виде:

д 2и

Рз= У ' + м(Уи)т + А • Е • 1т(Уи)) + р*Е, (5)

где р, А — коэффициенты Лямэ, значения которых связаны со значениями коэффициента Пуассона и модуля Юнга. Проинтегрируем уравнение (5) по объему расчетной ячейки и перейдем от объемных интегралов к поверхностным:

/ рзди¿V = I (рУи + р(Уи)т + А • Е • ¿г(Уи)) ¿Б + ^ рзЕсУ.

П дП П

При выполнении интегрирования будем использовать теорему о среднем, принимая в качестве среднего значения по объему — значение в центре ячейки, а в качестве среднего значения на грани — значение в центре грани:

(рзди) Зр = Е (руи + р(уи)т + А • Е • МУи)^ • Sf + (рзЕ)р Зр, (6)

где Зр — объем ячейки, — вектор площади грани. Здесь суммирование осуществляется по граням контрольного объема.

При аппроксимации члена диффузионного типа (^Уи + м(Уи)т + А • Е • Ьт(Уи)) . необходимо разбить его на явную и неявную части. Здесь под неявной частью оператора понимается та часть, которая заносится в коэффициенты разностного уравнения, а под явной — та часть, которая записывается в правую часть. Это является очень важным моментом, поскольку от способа этой разбивки будет существенным образом зависеть скорость сходимости разностной задачи. В данном случае наибольшей скорости сходимости удалось добиться при использовании следующей схемы:

(«Уи + м(Уи)т + А • Е • Ьт(Уи)) = ((2^ + А)Уи). +

1 --„->

неявно

+ (м(Уи)т + А • Е • Ьт(Уи) - (м + А)Уи) 1 • 4-^-'

явно

Дискретизация по времени осуществляется неявным способом, в программе реализованы следующие схемы аппроксимации второй временной производной: а) неявная схема второго порядка:

д2и дЬ2

{ \П+1 о / \П . { \П—1 ^ (и)р -2(и)р + (и)р .

б) неявная схема четвертого порядка:

2 (и)П+1 - 5 (и)Р + 4 (и)^-1 - 4 (и)п-2

д2и ^ ^У")р ' )р (7)

р '2

дЬ2

где т — шаг по времени, п — номер временного слоя.

Методические расчеты показали, что использование схемы четвертого порядка аппроксимации позволяет существенно повысить точность определения собственных частот колебаний упругого тела. Поскольку для данной задачи это обстоятельство очень важно, то все представленные ниже результаты были получены при помощи схемы (7).

Дальнейшая дискретизация уравнения (6) не представляет большого труда и подробно описана в литературе [7-8].

Верификация расчетного алгоритма проводилась на большом количестве различных задач упругости и аэроакустики. Здесь мы лишь перечислим рассмотренные тестовые задачи: задача о деформации нагруженной тонкой перфорированной пластины; расчет распространения волны деформации в металлической балке; определение собственных частот продольных и поперечных колебаний стержня; задача о звуковой волне, излучаемой точечным источником; задача об акустическом распаде разрыва.

Проведенное тестирование показало, что численные решения всех рассмотренных задач качественно и количественно хорошо согласуются с аналитическими решениями, экспериментальными данными и расчетами других авторов.

3. Моделирование работы датчика вязкости камертонного типа

Датчик, вырезанный из пьезокристалла Ы^ОЗ в форме камертона, применяется для измерения плотности и вязкости жидкой среды. Размеры датчика представлены на рис. 1.

Рис. 1. Модель датчика

Рис. 2. Схема колебаний датчика

Несмотря на то, что материал датчика является анизотропным, в силу того, что колебания в реальном датчике осуществляются практически в одном направлении (перпендикулярном оси датчика), то с хорошей степенью достоверности можно считать, что резонансные характеристики датчика (собственные частоты колебаний) будут определяться свойствами материала датчика именно в том направлении, в котором происходят колебания. Поэтому для моделирования работы реального датчика достаточно знать некий эффективный модуль упругости датчика Е (ГПа) в направлении, в котором происходят колебания. Выяснить, чему он равен, не составляет большого труда, поскольку из эксперимента нам известна собственная частота колебаний датчика в воздухе (25570Гц). Из теории мы знаем, что эта частота пропорциональна квадратному корню из модуля упругости. Сопоставив два этих факта, получаем, что модуль упругости в направлении преимущественных колебаний датчика должен быть равен 105,4 ГПа. В итоге для моделирования были выбраны следующие механические свойства материала датчика: модуль Юнга Е=105,4 ГПа, коэффициент Пуассона ^=0,3, плотность р=4700 кг/м3.

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

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

Рис. 3. Фрагменты расчетной сетки

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

Датчик помещен в закрытую емкость с жидкостью, расстояние от стенок датчика до стенок емкости — порядка нескольких миллиметров, что много меньше длины звуковых волн, генерируемых датчиком (несколько сантиметров), поэтому можно считать, что геометрия емкости не влияет на показания датчика. Методические расчеты подтверждают это предположение. В итоге размеры сосуда с жидкостью были взяты равными 10x5x18 мм.

Для расчета использовалось несколько разностных сеток различного пространственного разрешения. На рис. 3 показана сетка, содержащая 67x41x72 узла. Детализацию данной сетки можно считать оптимальной, поскольку полученное при помощи этой сетки решение хорошо согласуется с решением на более детальных сетках, и при этом на получение одного варианта расчета требуется приемлемое количество машинного времени (порядка 12 часов).

Для оценки влияния свойств жидкости на характеристики колебаний датчика была проведена серия расчетов колебаний датчика в вязкой жидкости. В расчетах варьировалась плотность (в диапазоне от 1 до 1113 кг/м3) и вязкость жидкой среды (в диапазоне от 0,0002 до 0,02 Пах с). Зависимость от каждого из этих параметров изучалась по отдельности.

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

Представленные на рисунках результаты приведены для варианта, в котором плотность

Рис. 4. Распределение вектора деформации в центральном сечении датчика в различные моменты времени (диапазон изменения от 0 до 8х10-8 м)

жидкости равна 660кг/м3, вязкость 0,0002 Пахс. Пространственное распределение и характер нестационарного поведения локальных величин (вектора деформации, вектора скорости жидкости, давления и т.д.) для других вариантов расчета в целом качественно аналогичные, естественно, с учетом того, что в более плотных и вязких жидкостях амплитуда и добротность колебаний датчика уменьшаются.

Различие рассмотренных вариантов расчета становится наиболее заметным, если рассматривать интегральные характеристики колебаний датчика. На рис. 6 приведено несколько спектров колебаний датчика для вариантов расчета с различной плотностью жидкой среды (по оси X отложен период колебаний, по оси У Фурье-преобразование полной механической энергии датчика). Вязкость жидкости при этом была фиксирована и равна 0,0002 Пахс. Здесь мы приводим только ту часть спектра, которая соответствует минимальной (основной) собственной частоте колебаний датчика. Другие частоты много больше и имеют меньшую интенсивность. Поэтому анализировать будем только поведение основной частоты. Из анализа графиков на рис. 6 следует, что увеличение плотности жидкости приводит к уменьшению добротности колебаний датчика, а также к уменьшению основной собственной частоты. Добротность колебаний датчика в данном случае определялась следующим образом: Q = //5, где / — собственная частота колебаний датчика, 5 — ширина резонансного пика в том месте, где его амплитуда падает в 1,414 раза.

Наблюдается так называемый сдвиг частот (разность основной частоты колебаний тела в воздухе и основной частоты колебаний тела в жидкой среде). На рис. 7 представлена зависимость величины сдвига основной частоты датчика от плотности жидкости, в которую он погружен. Как и следовало ожидать, зависимость линейная. Именно на этой зависимости основан принцип измерения виброакустическими датчиками плотности жидких сред.

Рис. 5. Изолинии относительного давления в центральных сечениях датчика в различные моменты времени (диапазон изменения от -19500 до 17700 Па)

1 1 1 1 кг.'мЗ 250 КГ/ЫЗ 1 г 1 1

500 КГ/«3 йбО кг/гаЗ \

1 (Н) 1 750 кг/йа * /К

ООП £>00 кг/ьв Щщ имД\л / ■.. - ■ - ^

Ц Ц К* мЗ X Т,с

1 .5-10

РРТ

ме

15 Л)5

У> 10

4 10

«¿10

гл

5 105 ЗЛО!"5

в-10'

Рис. 6. Зависимость спектра колебаний датчика от плотности жидкости

Рис. 7. Зависимость сдвига основной частоты колебаний датчика от плотности. Вязкость жидкости 0,0002 Пахс.

Влияние вязкости среды на характеристики колебаний датчика можно увидеть на рис. 8 (плотность жидкости 660кг/м3). Как видно, наличие вязкости приводит к затуханию колебаний датчика и, как следствие, к уширению резонансного спектра, причем величина этого уширения, как и в эксперименте, обратно пропорциональна квадратному корню из вязкости жидкости. Данный факт отражен на рис. 9, где Q — добротность колебаний датчика. Линейная зависимость служит хорошим поводом для использования подобного устройства в качестве прибора для измерения вязкости в широком диапазоне.

В заключение для сравнения сопоставим расчетную зависимость сдвига резонансной частоты колебаний датчика с экспериментальной зависимостью (рис. 10). Экспериментальные измерения резонансных характеристик датчика для пяти различных жидкостей (гексан, декан, додекан, этиленглюколь, раствор клейковины) проведены и предоставлены компанией Бакер Хьюз. Как видно из рис. 10, расчетные и экспериментальные данные достаточно хорошо согласуются между собой.

Выводы

Разработана и адаптирована математическая модель процесса акустического взаимодействия упругого колеблющегося тела с вязкой сжимаемой жидкостью. Предложенный алгоритм учитывает нестационарность и трехмерность процесса движения сенсора и жидкости, легко адаптируется к реальной геометрии любого сложного объекта. Тестирование разработанной математической модели и численного алгоритма на решении ряда задач теории упругости показало высокую его эффективность. Численные решения всех рассмотренных тестовых задач хорошо качественно и количественно согласуются с соответствующими аналитическими решениями, экспериментальными данными и расчетами других авторов. Данный численный алгоритм был использован для моделирования работы датчика вязкости камертонного типа. В результате обнаружена линейная зависимость сдвига основной соб-

б 1С

4 10 *

РРТ

а.ю

1>Ю"

л

1 1 1 1 \ 1

— 0 0002 Па ' с

---- 0 00034 Ш1 С |

" 0 0014 Па1 с р 1

эсга 0.003 Па1 с [у

0.02 Па*С Ё ¡/а 1 \\

пгттищниари^!!1""8^8^ т, с ! е -и-. -|а=нзг- ¡¡з ц ц—и-

3 ш

3-5 10

4 Е0

4.5 10

5 10

5.5 10

6 10

Рис. 8. Зависимость резонансных характеристик датчика от вязкости жидкости

Рис. 9. Зависимость обратной добротности колебаний датчика от вязкости. Плотность жидкости 660кг/м3.

^ 4000

Z

I-

о

о

{U :г

О х

0

1

(U X

о

о ф

о.

3000

2000

1000

ш

Г[

о

1 1 1 1 1 ООО расчет

##+ эксперимент

Ж""

-

О'

е- 1 1 1 1

200 400 600 S00 1000

плотность жидкости, кг/мЗ

1200

Рис. 10. Сопоставление расчетной и экспериментальной зависимостей относительного сдвига резонансной частоты

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

Данная работа частично финансировалась Бейкер Атлас 'российским научным центром.

Список литературы

[1] B.Robert, Quarts tuning fork resonators transducers, 8th Quarts devices conferences and exhibitions, 1986.

[2] L.F.Matsiev, Application of flexural mechanical resonators to high throughput liquid characterization, IEEE International Ultrasonics Symposium, San Juan, 2000, 8.

[3] Л.Д.Ландау, Е.М.Лившиц, Теория упругости, М., Наука, 1987.

[4] Л.Д.Ландау, Е.М.Лившиц, Гидродинамика, М., Наука, 1986.

[5] А.В.Минаков, А.А.Гаврилов, А.А.Дектерев, Численный алгоритм решения пространственных задач гидродинамики c подвижными твердыми телами и свободной поверхностью, Сибирский журнал индустриальной математики, 36(2008), №4, 95-105.

[6] С.Патанкар, Численные методы решения задач теплообмена и динамики жидкости, М., Энергоатомиздат, 1984.

[7] Ю.А.Быстров, С.А.Исаев, Н.А.Кудрявцев, А.И.Леонтьев, Численное моделирование вихревой интенсификации теплообмена в пакетах труб, СПб., Судостроение, 2005.

[8] J.H.Ferziger, M.Peric, Computational Methods for Fluid Dynamics, 3., rev. ed., Springer Verlag, Berlin, 2002.

Numerical Simulation of Working Process of Viscosity Turning Fork Sensor

Andrey V.Minakov Alexander A.Dekterev Andrey A.Gavrilov Valery Ya.Rudyak

A finite volume algorithm is developed and applied to compute working process tuning fork sensor. In this paper the results of numerical algorithm application are presented. The linear dependence of an oscillation frequency on density is found. The linear dependence of a quality factor of oscillations on the square root from dynamic viscosity of a fluid is found. The results show good agreement with experimental data and analytical solutions.

Keywords: aeroacoustics, TSM resonators, numerical simulation, finite volume method.