CC BY
16
DOI: 10.15514/ISPRAS-2020-32(5)-14
Моделирование инфразвукового пистонфона
Д.В. Головин, ORCID: 0000-0002-9488-6505 <golovin@vniiftri.ru>
Всероссийский научно-исследовательский институт физико-технических и радиотехнических измерений, Госстандарт, 141570, Московская область, Солнечногорский район, г.п. Менделеево
Аннотация. Представлены результаты численного моделирования прикладной акустической задачи -моделирования газовых процессов, протекающих в измерительной камере инфразвукового пистонфона 3202 при различных частотах колебаний поршня (0,1 - 1000 Гц) и характеризующихся крайне малыми значениями числа Маха (9,110"7^9,110"3). Моделирование проведено с использованием квазигазодинамических (КГД) и квазигидродинамических (КГиД) уравнений вязкого сжимаемого теплопроводного газа, определены границы устойчивости алгоритмов КГД и КГиД в данной задаче.
Ключевые слова: акустика; инфразвук; численное моделирование; квазигазодинамические уравнения; квазигидродинамические уравнения
Для цитирования: Головин Д.В. Моделирование инфразвукового пистонфона. Труды ИСП РАН, том 32, вып. 5, 2020 г., стр. 181-198. DOI: 10.15514/ISPRAS-2020-32(5)-14
Благодарности. Автор благодарен Елизаровой Т.Г. за поддержку, консультации и обсуждение результатов при проведении данной работы.
Simulation of Infrasound Pistonphone
D.V. Golovin, ORCID: 0000-0002-9488-6505 <golovin@vniiftri.ru>
All-Russian Scientific Research Institute of Physical Technical and Radiotechnical Measurements, Rosstandart, Mendeleevo, Solnechnogorsky district, Moscow region, 141570
Abstract. There are presented the results of numerical simulation of an applied acoustic problem - modeling of gas processes occurring in the measuring chamber of the infrasound pistonphone 3202 at different frequencies of piston oscillation (0.1 - 1000 Hz) and characterized by extremely small Mach numbers (9.110-7^9.110-3). The simulation was performed using quasi-gas-dynamic (QGD) and quasi-hydrodynamic (QHD) equations of a viscous compressible heat-conducting gas with the use of a time-explicit difference scheme, all spatial derivatives was approximated by central differences. It is shown that QGD and QHD models can be used for a simulation of applied acoustics and, in particular, to the simulation of infrasonic pistonphone: the stability limits of the QGD and QHD algorithms in this problem were determined, the dependence of sound pressure on the tuning parameter a is investigated and it is shown that this dependence is quite small. The spectra of sound pressure at the control point calculated by QGD and QHD are given, their dependence on the tuning parameter a is shown, both models equally predict the value of the sound pressure amplitude at the fundamental frequency oscillations. At the end of the article, the sound pressure at the control point at the fundamental frequency oscillations obtained by using QGD and QHD is compared with the values calculated by using semi-empirical formula of sound pressure at closed volume for a case of small oscillations using the polytropic index obtained by Henry Gerber instead of the adiabatic coefficient.
Keywords: acoustic; infrasound; numerical simulation; quasi gas dynamic equations; quasi hydro dynamic equations
For citation: Golovin D.V. Simulation of infrasound pistonphone. Trudy ISP RAN/Proc. ISP RAS, vol. 32, issue 5, 2020, pp. 181-198 (in Russian). DOI: 10.15514/ISPRAS-2020-32(5)-14
Acknowledgements. The author is grateful to T.G. Elizarova for support, advice and discussion of the results in carrying out this work.
1. Введение
Актуальность исследования обусловлена появлением за последние 10 лет нового поколения инфразвуковых микрофонов, применяемых на станциях мониторинга за соблюдением Договора о всеобъемлющем запрещении ядерных испытаний (см., например, [1-3]). В России аналогом являются дифференциальные конденсаторные микробарометры К-304-АМ1И и ISGM-03M (вариант К-304-АМ1) с диапазоном частот вплоть до 0,001 - 63 Гц (рабочий диапазон частот от 0,02 до 4,0 Гц) и максимальной амплитудой входного сигнала не менее 100 Па, разработанные и производимые НТЦ «Геофизические измерения». Таким образом, для поверки и калибровки данных средств измерений необходимо разработать лабораторную установку, реализующую один из первичных (абсолютных) методов калибровки с рабочим диапазоном частот от 1 мГц до 50-100 Гц.
Рис. 1. Общий вид пистонфона 3202 Fig. 1. General view of pistonphone 3202 На таких частотах представляется удобным использовать классический метод пистонфона [4-7].Во ВНИИФТРИ в 1960-ые гг. был разработан пистонфон 3202, обладающий рабочим диапазоном частот 0,1 - 100 Гц (рис. 1). Поршень в основании измерительной камеры посредством электродинамической системы с подвижной катушкой приводится в движение, изменяет величину рабочего объема и тем самым создает колебания давления воздуха, действующее на устанавливаемый на верхку камеры эталонный измерительный микрофон. Перед проведением измерений проводится проверка на герметичность камеры с помощью избыточного давления, нагнетаемого резиновой грушей, и водяного манометра, соединенного с внутренним объемом камеры через штуцер. По величине двойного хода поршня, фиксируемого с помощью микроскопа по нанесенной на шток поршня светящейся метке, диаметра поршня и объема используемой измерительной камеры рассчитывается амплитуда звукового давления P, воздействующего на калибруемое средство измерения:
oirnjx
V 2'
где С = kh/sm(kh) - коэффициент волнового распределения звукового давления вдоль оси камеры; k - волновое число; h - высота камеры; N - коэффициент политропы, зависящий от частоты колебаний; Po - атмосферное давление в момент герметизации камеры; S - площадь поперечного сечения поршня пистонфона; V - объем камеры; X - двойная амплитуда колебаний поршня.
Коэффициент политропы N рассчитывается согласно исследованию Г. Гербера (H. Gerber) [8] для случая жесткого поршня по следующим формулам:
Nn
Е " v
N =1 +(у- 1)Е, Е = D^
_MJLa_ 1 2nfl2>
где у - коэффициент адиабаты; Е - комплексная передаточная функция температуры; D,Nn,Mn,l - члены, зависящие от геометрии камеры; а - коэффициент температуропроводности при постоянном объеме; f - частота колебаний воздуха в камере: i - мнимая единица.
Для камеры цилиндрической формы D, Nn, Mn и l имеют следующий вид:
8 1 (m + 1/2)2П2 + rtR2 V
D _? , Nmn I 2 , Mmn /о n 7T? , ^
n (m + 1/2)x2 (2R + 1)2 $
где m=0,1,2,..; n=0,1,2,..; xn - n-ый корень функции Бесселя нулевого порядка; R = -отношение внутренней высоты h к внутреннему диаметру d камеры; V и S - внутренний объем и внутренняя площадь стенок камеры.
Подобный расчет автоматически несет в себе допущения, принятые в [8] при выведении формулы зависимости коэффициента политропы от частоты колебаний (например, рассматривается случай малых колебаний при осреднении плотности, температуры и давления газа по всему объему камеры, т.е. не рассматриваются их объемные распределения). За последние 15 лет появилось несколько работ [9-11], посвященных исследованию результатов [8] и показывающих, что они недостаточно хорошо соответствуют современному уровню точности акустических измерений и что необходимо пересмотреть международный стандарт о первичной калибровке измерительных лабораторных микрофонов типа LS1 и LS2 методом взаимности IEC 61094-2 [12] в части калибровки микрофонов на инфразвуковых частотах, в которой используются результаты работы Г. Гербера. Таким образом, целью данной работы является проверка совпадения результатов численного моделирования звукового давления для сжимаемого вязкого теплопроводного воздуха со значением, вычисляемым по формуле (1), и результаты которой могут послужить для оптимизации характеристик пистонфона 3202.
2. Постановка задачи и газодинамические параметры
Измерительная камера пистонфона представляет собой цилиндр диаметром d=69,97 мм и высотой h=66,06 мм, поршень в основании камеры имеет диаметр dp=20 мм. Рассматриваются частоты колебаний поршня 0,1 - 1000 Гц, при которых поршень колеблется с амплитудой X=0,5 мм (необходимо отметить, что при таком колебании поршня в рассматриваемом диапазоне частот задача находится в рамках линейной акустики, поэтому моделирование не должно показать каких-либо сильных нелинейных эффектов). Начальные параметры воздуха: температура 70=296,15 K, давление p0=101325 Па, плотность р0=1,186 кг/м3, коэффициент вязкости д=1,8310-5 Пас, газовая постоянная R=288,5 Дж/(кгК), коэффициент адиабаты у=1,4, коэффициент теплопроводности /=0,0254 Вт/(м-К), число Прандтля Pr=0,728, показатель межмолекулярного взаимодействия ю=0,74.
Для описания движения поршня требуется учет движения границы и применение подробных пространственных сеток. Желательно также провести сопоставление между собой результатов, полученных при помощи различных численных подходов. Все это сложно выполнимо в рамках индивидуального кода, но может быть сделано в рамках имеющихся открытых пакетов, например, пакета OpenFoam. В данной работе выполнен упрощенный вариант постановки задачи в плане ее геометрии и условий на границе.
Рис. 2. Расчетная область измерительной камеры пистонфона Fig. 2. Calculation area of the pistonphone measurement coupler
Область расчета (рис. 2) полагается цилиндрической, а колебательно течение воздуха в камере пистонфона - осесимметричным. Для численного решения используются полные уравнения Навье-Стокса для вязкого сжимаемого теплопроводного газа в отсутствие внешних сил и источников тепла с применением их регуляризации двух видов - КГД и КГиД моделей [13-14]:
dp 1d(rjmr) djn
dt + г дг
■ + ■
дг
= 0,
Пт
д(риг) 1д(г]тгиг) д(]т2иг) др 1д(гПгг) дП2Г ^
--1----1---1--—---1----,
дЬ г дг дг дг г дг дг г
д(рщ) 1 д(г]'тги2) д(]тгиг) др _ 1 д(гПГ2) дПгг
дЬ г дг дг дг г дг дг
дЕ 1 д(г]'тгН) д(]тгН) 1 д(гдг) ддг
дЬ г дг дг г дг дг
1 д д — г~дг [Г(Пггиг + ПГ2и2)\ + — (П2гиг + П^щ),
где иг и и - проекции вектора скорости и на оси г и г; Е - полная энергия единицы объема; Н - полная удельная энтальпия:
Е=р-
ul + U2
+
р
Н =
Е + р
2 ' у — 1' р — рЯТ.
Компоненты вектора плотности вычисляются по формулам:
]тг — Р(иг — шr), ]'тг — Р(иг —
Величины шг и , компоненты тензора вязкости П и проекции на оси г и г теплового потока в КГД и КГиД считаются по-разному.
Граница расчетной области состоит из оси симметрии, поршня и стенок камеры. Движение поршня в произвольный момент времени / описывается с помощью граничного условия для азимутальной скорости:
ир151оп — иР^оп — 2л^Х • соб(2л¡Ь).
Также на границе поршня задано условие прилипания для радиальной скорости, постоянство температуры и непроницаемости:
ur = 0, Т = Т0,
др дг
= 0.
На оси симметрии заданы условия симметрии:
д u
др_0 дг
= 0,
д р
ur = 0, -:т=0.
дг дг
На стенках камеры - условия прилипания для скорости, постоянства температуры и непроницаемости:
ur = 0, uz = 0,
д р
Т = Т0, j-=0, дп
где n - нормаль к поверхности.
В качестве колебания давления Ap(t) (и, впоследствии, звукового давления) в момент времени t рассматривается величина
hp(t) = p(t)-p0,
где p(t) - полное давление в момент времени t; p0 - атмосферное давление. 2.1 Уравнения КГД
Величины ш* и , компоненты тензора вязкости П и проекции на оси r и z теплового потока в КГД и КГиД считаются по-разному. В случае КГД формулы расчета выглядят следующим образом:
ш.
QGD _
1 д д др -— (гри+ — (rpUrUz) + —
г дг дг дг.
ш,
QG D
1 д д др' -— (rpUrUz) + — (гри+ — г дг дг дг.
п
QGD _ nNS
Пгг + Ur^r + R
п Q
QGD _ rjNS
п
QGD _ j-TNS
Пгт + UZK,
zr lizr
j-fQGD j-rNS I D* пфф пфф + R ,
nZQzGD = nz + UzV* + R*
дUr 2
nNS = 2p —---p divu,
дг 3
'rz 'zr "\дг + дг)'
N4 ur 2 Пфф = 2p—— — p divu,
дщ 2
nNS = 2p —---p divu.
zz ^ дг 3 й
При этом величины ш*, ш** и R* вычисляются так:
д ur д ur д р
Ш* = Т pUr~ + puz —— + -— r д z д д
ш =
R* =
д u
д uz д р
PUr^ + PUz^ + Tz_ д р д р
ur— + uz — +YP divu r д z д
Формулы для компонент теплового потока:
Ri = тр
qQGD = qNS — urRq,
qNS "
qr Xдг, ur д р uz д р
qQGD =qNzS-UzR*
-д(Р)+^д(Р)+pur д (1)+puz д (1)
1дг(р) y — 1дг (p) р дг (р) р z дг(р)
р
р
*
2.2 Уравнения КГиД
По своему виду КГиД уравнения проще уравнений КГД, и уравнения замыкания в ней выглядят как:
ohd ( диг диг l др\ 0HD ( duz duz l др\
=Ли^+и^+рд;) =T[Ur~fr+u*~dF+pw
Чг = 4rs 4z = qNzs
n?rHD = nNs + ur«0HD, n?zHD = nrNs + ur«°°HD,
n0HD =nNS +u ^0hd n°HD = nNs
zr — lLrz ~ Uzwr ,
n°HD _ nNs , ,,°HD Uzz = nzz + uz^z .
Следует отметить, что при т=0 уравнения КГиД и КГД переходят в уравнения Навье-Стокса.
2.3 Сеточная аппроксимация уравнений и обезразмеривание газодинамических переменных
Для сеточной аппроксимации уравнений используются явные аппроксимации по времени и центральные разностные производные по пространству на равномерной сетке с одинаковыми шагами по координатным направлениям [13, 5.5]. Регуляризующий параметр т был задан как
h
т = а—, с
где а - схемный параметр, не изменяющийся в процессе расчета; h - шаг сетки; c - скорость звука в невозмущенной среде.
Как будет показано далее, шаг интегрирования по времени необходимо выбирать в соответствии с критерием Куранта.
Для удобства работы с уравнениями основные величины приведены к безразмерному виду для выделения параметров подобия, от которых зависит решение задачи (числа Маха и Рейнольдса).
Связь между размерными и безразмерными величинами (обозначены тильдой) выглядит следующим образом:
р = рр0, р = pp0(2nfX)2, и = и ■ 2nfX,
_ dp _ (2nfX)2
r = fdp, z = zdp, t = t — , T = T-7^.
Такое обезразмеривание не изменяет вид уравнений КГД и КГиД, но изменит вид уравнения состояния воздуха:
„ = Р? р Y .
Безразмерные коэффициенты вязкости, теплопроводности и параметр т (знак тильда у безразмерных величин опущен) вычисляются следующим образом:
l V
р = —(м2Т)ш, х= „ , T = ahM,
р Rey J л Pr(Y-l)
где M и Re - числа Маха и Рейнольдса:
2nfX p0(2nfX)dv
M = —!—, Re = v.
с р0
Для рассматриваемых частот 0,1 - 1000 Гц значения чисел Маха и Рейнольдса меняются от 9,1-10-7 и 0,204 до 9,1-10-3 и 2040 соответственно.
Головин Д.В. Моделирование инфразвукового пистонфона. Труды ИСП РАН, том 32, вып. 5, 2020 г., стр. 181-198
3. Результаты моделирования
3.1 Устойчивость уравнений явной аппроксимации по времени КГД и КГиД и влияние параметра а
Для случая явной аппроксимации по времени используемых уравнений, согласно теоретическим оценкам [12-13], шаги по времени йЬ и пространству к подчиняются критерию Куранта:
М = (2)
иmax
где 0 < р < 1 - число Куранта; Ьт1П - минимальный шаг применяемой сетки; с„ максимальная локальная скорость звука в расчетной области на текущем шаге по времени. Для рассматриваемой задачи, как было сказано выше, характерны крайне малые значения числа Маха, поэтому, при наличии в уравнениях регуляризующих слагаемых, вопрос изучения границ применимости формулы (2) является естественным.
Исследование устойчивости проводилось следующим образом: для двух частот колебаний поршня (10 и 100 Гц) и двух равномерных сеток (53х100 и 79х150) фиксировалось число а в диапазоне от 0,1 до 1,0 и выбиралось такое число р, при котором расчет являлся устойчивым вне зависимости от наличия схемных осцилляций в распределениях плотности, давления, температуры и скоростей. Дополнительным критерием являлась устойчивость колебаний давления с течением времени. Границы устойчивости определялись относительно грубо: если при зафиксированном а расчет расходился, то первая значащая цифра числа Куранта р уменьшалась на единицу. Таким образом подбиралось максимально возможное значение числа Куранта р устойчивого счета для заданного значения а.
0.1 0.2 03 0.4 0.S 0.6 0.7 О.в 0.9
а
-о ОНО, 10 Hz QHD. 1С0 Hz -9-QGD. 10 Hz ""-QGD, 100 Hz
0.05 !
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.S 0.3 1
Рис. 3. Зависимость числа Куранта 1 от параметра а на частотах 10 Гц и 100 Гц
для сеток 53x100 и 79x150 Fig. 3. The dependence of Courant number 1 from parameter а at the frequencies of 10 Hz and 100 Hz,
the grids 53x100 and 79x150
На рис. 3 представлены результаты этого исследования. Как видно из рисунков, даже при таких малых значениях чисел Маха и Рейнольдса выполняется условие Куранта (2), а зависимость в(а) определяется размером сетки и в рамках заданной точности подбора в не зависит от частоты колебаний и одинакова для моделей КГД и КГиД.
На рис. 4-7 приведены колебания давления в контрольной точке на частотах 10 и 100 Гц для моделей КГД и КГиД.
Рис. 4. Колебания давления в контрольной точке для различных значений параметра а на частоте колебаний поршня 10 Гц для модели КГД, сетка 53x100 Fig. 4. Pressure deviation in control point at different value ofparameter a at the piston frequency oscillation of 10 Hz, QGD model, the grid 53x100
Рис. 5. Колебания давления в контрольной точке для различных значений параметра а на частоте колебаний поршня 100 Гц для модели КГД, сетка 53x100 Fig. 5. Pressure deviation in control point at different value ofparameter a at the piston frequency oscillation of 100 Hz, QGD model, the grid 53x100
Рис. 6. Колебания давления в контрольной точке для различных значений параметра а на частоте колебаний поршня 10 Гц для модели КГиД, сетка 53x100 Fig. 6. Pressure deviation in control point at different value ofparameter а at the piston frequency oscillation of 10 Hz, QHD model, the grid 53x100
Рис. 7. Колебания давления в контрольной точке для различных значений параметра а на частоте колебаний поршня 100 Гц для модели КГиД, сетка 53x100 Fig. 7. Pressure deviation in control point at different value ofparameter а at the piston frequency oscillation of 100 Hz, QHD model, the grid 53x100
По мере уменьшения числа а для модели КГиД при частоте колебаний поршня 100 Гц увеличиваются осцилляции на кривой колебания давления (рис. 7). Если еще больше уменьшить значение а (меньше 0,1), то устойчивых колебаний давления не получится, хотя для модели КГиД в области а<0,1 получаются гладкие распределения плотности, давления, температуры и скоростей (при а>0,1 на сетке 53х100 на исследованных частотах присутствуют незначительные схемные осцилляции, уменьшающиеся при увеличении размеров сетки).
Рис. 8. Фурье-спектр давления в контрольной точке для различных значений параметра а на частоте колебаний поршня 10 Гц для моделей КГД и КГиД, сетка 53x100 Fig. 8. Fourier pressure amplitude spectrum in control point at different value ofparameter а at the piston frequency oscillation of 10 Hz, QGD and QHD models, the grid 53x100
Рис. 9. Фурье-спектр давления в контрольной точке для различных значений параметра а на частоте колебаний поршня 100 Гц для моделей КГД и КГиД, сетка 53x100 Fig. 9. Fourier pressure amplitude spectrum in control point at different value ofparameter a at the piston frequency oscillation of 100 Hz, QGD and QHD models, the grid 53x100 На рис. 8-9 представлены спектры колебаний давления в контрольной точке, соответствующие рис. 4-7. Модели КГД и КГиД предсказывают достаточно одинаковые значения амплитуды давления на основной частоте колебаний поршня, при определенных значениях числа а даже совпадают значения вторых гармоник - на частоте 10 Гц при а=1 спектры давления для КГД и КГиД очень хорошо совпали друг с другом (рис. 8). Если выбрать контрольную точку в другой части камеры, то на частотах 40-100 Гц значение
амплитуды колебаний давления основной частоты будут отличаться до 1 Па при расчете по КГД и КГиД.
Также была исследована зависимость амплитуды колебания давления в контрольной точке от значения параметра а - определенное влияние а присутствует, однако его сложно назвать серьезным (рис. 10).
Рис. 10. Зависимость амплитуды колебаний давления на частоте колебания поршня в контрольной
точке от параметра а на частотах 10 Гц и 100 Гц для моделей КГД и КГиД, сетка 53x100 Fig. 10. The dependence ofpressure deviation amplitude at the piston frequency oscillation in control point from parameter а at the frequencies of 10 Hz and 100 Hz, the QGD and QHD models, the grid 53x100
В заключение остается добавить, что при а = т = 0 не удалось подобрать такого значения шага по времени, чтобы вычисления были устойчивыми.
3.2 Звуковое давление в камере пистонфона
Теперь, убедившись на основе результатов подраздела 3.1, что модели КГД и КГиД действительно предсказывают колебания давления, слабо зависящие от настроечного параметра а и которые с хорошей точностью можно считать синусоидальными, далее будем использовать термин «звуковое давление».
В контрольной точке (см. рис. 2) на основе вычислений от 5 (на частотах 0,1 - 10 Гц для КГД и 1 - 10 Гц для КГиД) до 20 (в остальном диапазоне частот) периодов колебаний на сетке 53x100 определялся Фурье-спектр звукового давления и значение амплитуды звукового давления на основной частоте колебаний, которое затем сравнивалось с рассчитанным по формуле (1). Расчеты проводились при а=0,5.
Примеры линий тока и распределения звукового давления согласно КГД показаны для двух моментов времени на рис. 11-14.
0.01 D.02 0.03 0.04 0.05 0.06
Height, tn
Рис. 11. Линии тока и изолинии звукового давления на частоте f=10 Гц в расчетной области на сетке
53x100, параметры а=0,5 и в=0,3, число периодов с начала колебаний - 2,25. Расчет согласно КГД Fig. 11. Streamlines and sound pressure isolines at the frequency of 10 Hz at the calculation area on the grid 53x100, parameters a=0.5 and Courant number в=0.3, the number of periods since the beginning of oscillations - 2.25. Simulation according to QGD
0.01 0.02 0.03 0.04 0.05 0.06
Height, m
Рис. 12. Линии тока и изолинии звукового давления на частоте f=10 Гц в расчетной области на сетке
53x100, параметры а=0,5 и в=0,3, число периодов с начала колебаний - 2,75. Расчет согласно КГД Fig. 12. Streamlines and sound pressure isolines at the frequency of 10 Hz at the calculation area on the grid 53x100, parameters a=0.5 and Courant number в=0.3, the number of periods since the beginning of oscillations - 2.75. Simulation according to QGD
0.01 D.02 0.03 0.04 D.05 0.06
Height., in
Рис. 13. Линии тока и изолинии звукового давления на частоте f=100 Гц в расчетной области на сетке 53x100, параметры а=0,5 и в=0,3, число периодов с начала колебаний - 2,25. Расчет согласно
КГД
Fig. 13. Streamlines and sound pressure isolines at the frequency f=100 Hz at the calculation area on the grid 53x100, parameters a=0.5 and Courant number в=0.3, the number of periods since the beginning of oscillations - 2.25. Simulation according to QGD
0.01 0.02 0.03 0.04 0.05 0.06
Height, m
Рис. 14. Линии тока и изолинии звукового давления на частоте f=100 Гц в расчетной области на сетке 53x100, параметры а=0,5 и в=0,3, число периодов с начала колебаний - 2,75. Расчет согласно
КГД
Fig. 14. Streamlines and sound pressure isolines at the frequency f=100 Hz at the calculation area on the grid 53x100, parameters a=0.5 and Courant number в=0.3, the number of periods since the beginning of oscillations - 2.75. Simulation according to QGD На рис. 15-18 приведены распределения звукового давления и линии тока, рассчитанные по КГиД.
Pressure deviation. Pa. f=10 Hz. NritNz=53x100. alpha-G.S, !m!';i С ."( period 2.25
-j-----!— -1- - МШ
ej.s
84.BM
84В9в
B4.B94
А4.Й&2
B4.es
04.3ЙВ
BJ.SBb
84.S&4
EW.BB2
0.01 0.02 0.03 0.04 0.05 0.06
Height, tn
Рис. 15. Линии тока и изолинии звукового давления на частоте f=10 Гц в расчетной области на сетке 53x100, параметры а=0,5 и в=0,3, число периодов с начала колебаний - 2,25. Расчет согласно КГиД Fig. 15. Streamlines and sound pressure isolines at the frequency f=10 Hz at the calculation area on the grid 53x100, parameters a=0.5 and Courant number в=0.3, the number of periods since the beginning of oscillations - 2.25. Simulation according to QHD
0.01 0.02 0.03 0.04 0.05 0.06
Height, m
Рис. 16. Линии тока и изолинии звукового давления на частоте f=10 Гц в расчетной области на сетке 53x100, параметры а=0,5 и в=0,3, число периодов с начала колебаний - 2,75. Расчет согласно КГиД Fig. 16. Streamlines and sound pressure isolines at the frequency f=10 Hz at the calculation area on the grid 53x100, parameters a=0.5 and Courant number в=0.3, the number of periods since the beginning of oscillations - 2.75. Simulation according to QHD
0.01 D.02 0.03 0.04 D.0S 0.06
Height., in
Рис. 17. Линии тока и изолинии звукового давления на частоте f=100 Гц в расчетной области на сетке 53x100, параметры а=0,5 и в=0,3, число периодов с начала колебаний - 2,25. Расчет согласно
КГиД
Fig. 17. Streamlines and sound pressure isolines at the frequency f=100 Hz at the calculation area on the grid 53x100, parameters a=0.5 and Courant number в=0.3, the number of periods since the beginning of oscillations - 2.25. Simulation according to OHD
0.01 D.02 0.03 0.04 0.05 0.06
Height, tn
Рис. 18. Линии тока и изолинии звукового давления на частоте f=100 Гц в расчетной области на сетке 53x100, параметры а=0,5 и в=0,3, число периодов с начала колебаний - 2,75. Расчет согласно
КГиД
Fig. 18. Streamlines and sound pressure isolines at the frequency f=100 Hz at the calculation area on the grid 53x100, parameters a=0.5 and Courant number в=0.3, the number of periods since the beginning of oscillations - 2.75. Simulation according to QHD Хотя на кривых звукового давления для расчета по КГиД при а=0,5 практически нет осцилляций (только в начале первого периода колебаний, как на рис. 6-7), линии тока и изолинии звукового давления говорят о наличии схемных осцилляций на используемой сетке, как упоминалось выше. Также из приведенных выше рисунков можно видеть различие в предсказаниях изолиний звукового давления и линий тока у КГД и КГиД, но это
практически не сказывается на значении звукового давления в контрольной точке на основной частоте колебаний, в чем можно убедиться по спектрам на рис. 8-9.
Рис. 19. Сравнение результатов численного моделирования звукового давления на сетке 53х100 с
вычислением по формуле (1) Fig. 19. Comparison of numerical simulation results of sound pressure on the grid 53x100 with calculation
by the formula (1)
Звуковое давление для модели КГД вычислено в диапазоне частот 0,1 - 1000 Гц, для КГиД -от 1 до 250 Гц (на рис. 19). Это различие вызвано прежде всего тем, что для моделирования частот ниже 1 Гц требуется существенный объем времени. При этом на частотах выше 100 Гц размах схемных осцилляций распределения звукового давления для расчета по КГиД становится значительным, хотя вычислительная устойчивость сохраняется вплоть до исследованной частоты 1000 Гц. Следует отметить, что в области частот 1 - 250 Гц результаты моделирования с помощью КГиД уравнений очень хорошо совпадают со значениями, полученными применением уравнений КГД, и, основываясь на этих данных, можно ожидать, что совпадут также и на частотах 0,1 - 1 Гц.
Если предположить, что зависимость звукового давления от частоты слабо зависит от размера расчетной сетки, то, возможно, имеет смысл сравнить относительное изменение кривых звукового давления для аналитической формулы (1), КГД и КГиД:
A=20\g(-^) [дБ/dB],
где Л - изменение амплитуды звукового давления Р в контрольной точке на частоте колебаний поршня относительно ее значения Pref на опорной частоте 100 Гц. Соответствующий график представлен на рис. 20. Как по нему видно, до частоты 20 Гц относительное изменение давления для численного моделирования и вычисленного по формуле (1) хорошо согласуются друг с другом, но на частоте 0,1 Гц различие составляет около 0,2 дБ (примерно 2 %). При этом результаты моделирования по КГД и КГиД в рассчитанном для КГиД диапазоне частот 1 - 250 Гц практически полностью совпали друг с другом.
Рис. 20. Сравнение изменения величины звукового давления относительно его значения на частоте
100 Гц для КГД, КГиД и расчета по формуле (1) Fig. 20. Comparison of changes of sound pressure value relative to its value at the frequency of 100 Hz for
QGD, QHD and calculation by the formula (1)
4. Заключение
Изначально КГД и КГиД модели были выбраны из-за простоты написания явной разностной схемы, ее распараллеливания и того факта, что эти уравнения хорошо описывают нестационарные течения. И оказалось, что модели КГД и КГиД можно применять к задачам акустики без использования специальных методов и лимитеров, что появляются при моделировании акустических колебаний посредством уравнений Навье-Стокса. При этом анализ устойчивости обеих моделей для чисел Маха, характерных в этой задаче, выполнен впервые.
Полученная разность между значением, рассчитанным по формуле (1), и значениями, рассчитанными для сжимаемого вязкого теплопроводного воздуха по КГД и КГиД, оказалась меньше 0,05 дБ в диапазоне частот 20 - 1000 Гц, но по мере снижения частоты разности между аналитическим и полученными с помощью КГД и КГиД моделей значениями увеличиваются и достигают 0,2 дБ (около 2 %) на частоте 0,1 Гц. При этом, как показано выше, результаты моделирования по КГД и КГиД моделям слабо зависят от параметра а в его исследованной области значений и достаточно точно совпадают друг с другом. При подтверждении результатов моделирования с помощью эксперимента это будет свидетельствовать как в пользу корректности применения КГД и КГиД для задач акустики и сверхмедленных течений (следует напомнить, что значения чисел Маха в рассмотренной области частот - 9,Г10_7^9,Г10~3), так и наличия систематической погрешности для калибруемых микрофонов при использовании акустической теории Г. Гербера на инфразвуковых частотах.
Дополнительной сложностью является оценка численного моделирования - прямое сравнение с результатами эксперимента на момент написания статьи невозможно из-за подготовки пистонфона 3202 к его модернизации и расширения рабочего частотного диапазона до 0,001 - 100 Гц. Схожей задачей - по своему исполнению - является численное моделирование камер связи установки первичной калибровки микрофонов Bruel&Kjaer Type 9699, в которой реализуется метод взаимности для абсолютной калибровки измерительных лабораторных микрофонов типа LS1 и LS2. Результаты такого моделирования уже можно сравнить с экспериментом, но подобное исследование, по сути, будет являться самостоятельным, и выходит за рамки данной работы.
В заключение следует отметить, что серьезным недостатком примененной явной разностной схемы является то, что для численного моделирования частот колебаний ниже 1 Гц тратится большое количество времени: если для расчета 20 периодов колебаний с частотой 100 Гц потребовалось около 25 минут на ПК с ЦПУ Intel i9-9900K (4,8 ГГц по всем ядрам), то расчет 5 периодов колебаний давления с частотой 0,1 Гц занял примерно 3,5 суток.
Список литературы / References
[1]. Merchant B. John. Hyperion 5113/A Infrasound Sensor Evaluation. Technical Report SAND2015-8097, Sandia National Lab., 2015, 43 p.
[2]. Slad George William, Merchant Bion J. Chaparral Model 60 Infrasound Sensor Evaluation. Technical Report SAND-2016-1902, Sandia National Lab., 2016, 42 p.
[3]. Merchant Bion J., McDowell Kyle D. MB3a Infrasound Sensor Evaluation. Technical Report SAND2014-20108, Sandia National Lab., 2014, 59 p.
[4]. Коньков А.В. О методе пистонфона. Сборник научных трудов ВНИИФТРИ, вып. 23 (33), 1975 г., стр. 5-11 / Konkov A.V. On the pistonphone method. Collection of scientific papers of VNIIFTRI, issue 23 (33), 1975, pp. 5-11 (in Russian).
[5]. Головин Д.В., Коньков А.В. Влияние внешних условий на звуковое давление в камере пистонфона в инфразвуковом диапазоне частот. Измерительная техника, №9, 2020 г., стр. 67-72 / Golovin D.V., Konkov A.V. The influence of external conditions on sound pressure in the pistonphone at the infrasound frequency range. Measuring technology, №9, 2020, pp. 67-72. (in Russian).
[6]. Zhang Fan, He Wen, He Longbiao, Rong Zuochao. Acoustic properties of pistonphones at low frequencies in the presence of pressure leakage and heat conduction. Journal of Sound and Vibration, vol. 358, 2015, pp. 324-333.
[7]. He Wen, He Longbiao, Zhang Fan, Rong Zuochao, and Jia Shushi. A dedicated pistonphone for absolute calibration of infrasound sensors at very low frequencies. Measurement Science and Technology, vol. 27, no.2, article no. 025018.
[8]. Gerber H. Acoustic Properties of Fluid-Filled Chambers at Infrasonic Frequencies in the Absence of Convection. Journal of Acoustical Society of America, vol. 36, issue 8, 1964, pp. 1427-1434.
[9]. Guianvarc'h C., Durocher J.-N., Bruneau M., Bruneau A.-M. Acoustic transfer admittance of cylindrical cavities. Journal of Sound and Vibration, vol. 292, issues 3-5, 2006, pp. 595-603.
[10]. Jackett R., Avison J. Realizing the primary standard for sound pressure: The trouble with IEC 61094-2. In Proc. of the 44rd International Congress on Noise Control Engineering, 2015, 9 p.
[11]. Vincent P., Rodrigues D., Larsonnier F., Guianvarc'h C., Durand S. Acoustic transfer admittance of cylindrical cavities in infrasonic frequency range. Metrologia, vol. 56, no. 1, 2019, article no. 015003.
[12]. IEC 61094-2:2009 Electroacoustics - Measurement microphones - Primary method for pressure calibration of laboratory standard microphones by the reciprocity technique. International Electrical Commission, Geneva, 2009
[13]. Елизарова Т.Г. Квазигазодинамические уравнения и методы расчета вязких течений. М., Научный Мир, 2007 г., 350 стр. / Elizarova T.G. Quasi-Gas Dynamic Equations. Springer Science & Business Media, 2009, 286 p.
[14]. Шеретов Ю.В. Динамика сплошных сред при пространственно-временном осреднении. М., Ижевск, НИЦ «Регулярная и хаотическая динамика», 2009 г., 400 стр. / Sheretov Y.V. Continuum dynamics under spatiotemporal averaging. M., Izhevsk, Scientific and Publishing Center «Regular and Chaotic Dynamics», 2009, 400 p. (in Russian).
Информация об авторе / Information about author
Дмитрий Витальевич ГОЛОВИН - младший научный сотрудник отдела акустики. Его научные интересы включают калибровку средств измерения звукового давления в воздухе на инфразвуковых частотах и численное моделирование задач акустики.
Dmitrii VitaTevich GOLOVIN - researcher of acoustic laboratory. His research interests include calibration of measuring instruments of sound pressure in the air at infrasound frequencies and simulation of applied acoustics.
