CC BY
51
Семенов Илья Сергеевич - Федеральное государственное автономное образовательное учреждение высшего профессионального образования «Южный федеральный университет»; e-mail: flanker555@yandex.ru; 347928, г. Таганрог, пер. Некрасовский, 44; тел.: 89085029807; кафедра высшей математики; аспирант.
Semenov Ilya Sergeevich - Federal State-Owned Autonomy Educational Establishment of Higher Vocational Education "Southern Federal University"; e-mail: flanker555@yandex.ru; 44, Nekrasovsky, Taganrog, 347928, Russia; phone: +79085029807; the department of higher mathematics; postgraduate student.
УДК 534.222
Д.В. Лапин, И.Н. Тетеревлев, О.А. Савицкий, Т.А. Чистякова
ПРОГРАММНЫЙ КОМПЛЕКС МОДЕЛИРОВАНИЯ ЗВУКОВЫХ ПУЧКОВ КОНЕЧНОЙ АМПЛИТУДЫ В КВАДРАТИЧНО-НЕЛИНЕЙНЫХ СРЕДАХ БЕЗ ФИЗИЧЕСКОЙ ДИСПЕРСИИ*
Разработан новый программный комплекс для моделирования звуковых пучков конечной амплитуды в квадратично-нелинейных средах без физической дисперсии (MSB). Программный комплекс позволяет выполнять научные исследования в области нелинейных волновых процессов, является мощным инструментом для решения инженерных задач разработки приборов и комплексов, основанных на принципах нелинейной акустики. Описаны принципы построения, математической модели, ее программно-алгоритмической реализации. Продемонстрированы возможности программного комплекса, определена область применения. Проведены тестовые расчеты и выполнена верификация результатов математического моделирования. Выполнено сравнение разработанного программного обеспечения с аналогичными программными комплексами.
Программный комплекс; математическое моделирование; звуковой пучок; волна конечной амплитуды.
D.V. Lapin, I.N. Teterevlev, O.A. Savitskiy, T.A. Chistyakova
SOFTWARE MODELING SOUND BEAM OF FINITE AMPLITUDE IN A QUADRATIC NONLINEAR MEDIUM WITHOUT PHYSICAL
DISPERSION
A new software package for the simulation of finite-amplitude sound beams in quadratic nonlinear media without physical dispersion. The software package allows you to perform scientific research in the field of nonlinear wave processes is a powerful tool for solving engineering problems of development of devices and systems based on the principles of non-linear acoustics. The principles of mathematical model, software and algorithmic implementation are described. The possibilities of the software system, sets out sphere of application. The test calculations are performed and verification of the results of mathematical modeling. The comparison of the developed software with the same program set.
Program complex; mathematical modeling; sound beam; finite amplitude wave.
Введение. Одним из перспективных направлений развития акустических и ультразвуковых технологий является применение нелинейных акустических эффектов. Акустическая аппаратура, в основе действия которой лежат принципы нелинейной акустики, обладает рядом уникальных возможностей, не свойственных традиционным акустическим устройствам. Примером удачного применения
* Работа выполнена при финансовой поддержке ФЦП «Научные и научно-педагогические кадры инновационной России на 2009-2013 гг.», проект №14.А18.21.0680.
нелинейных эффектов в акустике служит параметрическая антенна. На базе параметрических антенн в настоящее время разработаны и успешно эксплуатируются нелинейные гидролокаторы для работы на мелководье в узких каналах [1], параметрические профилографы [2], гидроакустические комплексы для охраны акваторий, параметрические эхоледомеры для подводного флота и др.
Большое разнообразие изученных нелинейных акустических эффектов в средах без дисперсии скорости звука [3, 4] открывает широкие возможности их использования в инновационных разработках. Например, при создании фазовых эхолокаторов для дистанционного зондирования и характеризации сред и объектов [5, 6, 7], измерителей нелинейных параметров материалов [8], в системах дальней локации с подавлением нелинейного поглощения [9], устройствах активного зву-когашения [10, 11] и др.
Вместе с тем, широкое внедрение нелинейных эффектов в прикладные разработки сдерживается сложностью описания нелинейных волновых процессов, отсутствием адекватных инженерных методик для расчета соответствующих технических устройств. В большинстве практически важных случаев отсутствуют точные аналитические решения уравнений нелинейной акустики, которые можно было бы положить в основу создания инженерных методик расчета.
Наиболее эффективным и экономически обоснованным подходом к разрешению обозначенной проблемы является математическое моделирование с применением современных вычислительных методов, алгоритмов и быстродействующих вычислительных средств, позволяющих выполнять прямое решение уравнений нелинейной акустики. Вместе с тем, в России в настоящее время отсутствуют программные комплексы для математического моделирования нелинейных волновых процессов, доступные для разработчиков акустической и ультразвуковой аппаратуры.
В научно-образовательном центре комплексных исследований и математического моделирования Южного федерального университета разработан программный комплекс математического моделирования нелинейных процессов в звуковых пучках волн конечной амплитуды (MSB) [12].
Разработанное программное обеспечение, в основе которого лежит численное решение уравнения акустики нелинейных звуковых пучков, позволяет выполнять вычислительный эксперимент при исследованиях нелинейных волновых процессов, а также производить инженерные расчеты основных параметров звуковых полей, их спектров, в широком диапазоне изменения начальных условий и внешних параметров исходного уравнения. В программном пакете MSB реализованы возможности трехмерной визуализации и вторичной обработки результатов моделирования.
Математическая модель. Рассмотрим принципы построения дискретной модели [13] и программно-алгоритмической реализации разработанного программного комплекса.
Процесс распространения звуковых пучков конечной амплитуды в квадратично-нелинейной диссипативной среде описывается уравнением Хохлова-Заболотской-Кузнецова, которое после применения процедуры обезразмеривания имеет вид:
где V = у( г, О, г) - величина скорости частиц среды; Г - диссипативный параметр;
О - время в сопровождающей системе координат; г - нормированное расстояние; N - параметр уравнения, характеризующий соотношение нелинейной и дифракционной длин волны; Д± — поперечный лапласиан
д dv dv ^ д 2v ^
М dZ " Vdd~ дМ J
(1)
Задача Коши для уравнения (1) решается при следующем начальном условии:
у(О,0, г) = V (в, г) (2)
и граничных условиях:
- условие периодичности сигнала
у(г,0,г) = у(г,2л,г), \в(г,0,г) = у'в(г,2л,г) (3)
- условие симметричности
уг(г,в,0 ) = 0, (4)
- условие отсутствия поля в бесконечно удаленной точке
V (г, в, ~) = 0. (5)
Решение дифференциального уравнения (1) ищется в области:
въ {0 <в< 2л,0 < Я < Я0 + Яьг,0 < г < г}.
Для разработки численной модели рассматриваемых волновых явлений отдается предпочтение методу расщепления по физическим процессам, так как двухслойная конечно-разностная схема, полученная путем замены частных производных их дискретными аналогами, не является устойчивой, а метод гармоник требует больших вычислительных затрат [14].
Расчетная область по пространственным направлениям X, у, г представляет
собой цилиндр (рис. 1).
Рис. 1. Расчетная область
Для построения разностной схемы уравнения (1), с учетом условий (2-5) будем использовать равномерную сетку:
ц?Н ={гг = пНг в = Л, гк = Я; п = О^, ] = 0М, к = 0Р; = I, МНв = 2л, РНГ = я}, где п, у , к - индексы по направлениям г, в, г соответственно; Н,,, кв, кг - шаги по направлениям г, в, г соответственно; N, М , Р - количество узлов сетки по направлениям г, в, г соответственно; I , Я - высота и радиус цилиндра соответственно.
Задача (1)-(5) решается в 2 этапа. На первом - учитываются нелинейные и диссипативные процессы, что соответствует решению уравнения Бюргерса, дифференциально-разностная аппроксимация которого имеет вид (6):
и — ип п дип+* ГЪ2ип+м
--и —--Г ^ , = 0. (6)
к7 дО до2
Уравнение (7) решается при следующих граничных условиях (с учетом периодичности волнового процесса):
и( г,0, г) = и (г,2л, г), (7)
и'о (г,0, г) = и'о (г,2я, г). (8)
На данном этапе осуществляется переход от пространственного слоя п по координате г к некоторому промежуточному слою п + о .
На втором этапе определяется величина скорости смещения частиц в волне с учетом дифракции, которая в рамках выбранного квазиоптического приближения заменяется диффузией энергии с мнимым коэффициентом в направлении, перпендикулярном направлению распространения волны:
д Л
дв
w - w h
N з
—A, w"+i, (9)
4 1
7 ^
где пг - шаг по пространству вдоль оси 7; М - величина скорости частиц среды в
поле волны на промежуточном временном слое; М - величина скорости волны на следующем временном слое.
Уравнение (9) решается при следующих граничных условиях:
м>'(г,д,0) = 0, м(г,О,^) = 0. (10)
На данном этапе осуществляется переход от некоторого промежуточного пространственного слоя к слою п +1.
Дискретная модель для исходной задачи представлена уравнениями (11) для задачи (6):
п+1 п п п+1 п п+1 п+1/2 Г) п+1/2 . п+1/2
и0,к — и0,к и1,ки1, к — им —1,киы —1, к т^и1,к — 2и0,к_ + им —1,к
2
р -~\,k----Q,k "M -1,k _ q j _ Q ;
h 4he К
n+1 n n n+1 n n+1 n+1/2 ry n+1/2 ■ n+1/2
Uj ,k - Uj ,k - Uj +1, kUj +1, k - Uj -1,kUj -1,k р Uj + 1,k - 2Uj, k + Uj -1,k _ Q
К 4h в '
1 < j <M-2; (11)
n+1 n n n+1 n n+1 n+1/2 ry n+1/2 . n+1/2
UM-1,k — UM-1,k UM-1,kUM-1,k — U1,kU1,k р U0,k — 2UM-1,k + UM-2,k _ q j _ m — 1
К 4Кв К
и уравнениями (12) для задачи (9):
n +1 n , т n+1/2 n+1/2
•V- Ci,k CJ,k д - C},Q
z«/ —-— _----—, k _ 0
h 2 К
ivjk
n+1 n ,T /' / -,4 n+1/2 n+1/2 / \ n+1/2 n+1/2 Л
k +11 Cn,k+1 -C ,k - f k -j - ^ '
cj ,k - cj ,k N
h 4
2J К2 l 21 h2
j
0 < k < P , (12)
j _Q. k _ P.
Коэффициенты c"j k на предыдущем временном слое находятся при помощи
алгоритма быстрого преобразования Фурье из поля wj k.
Разностная схема (12) абсолютно устойчива при: Л> 1/2и имеет второй порядок погрешности аппроксимации при Л = 1/2. Погрешность аппроксимации математической модели (11)-(12) равна O{hz + h2g j.
Особенности программно-алгоритмической реализации модели
Параллелизм. Для ускорения расчета модели используется факт наличия внутреннего параллелизма при решении уравнения Бюргерса и при вычислении серий прямых и обратных БПФ. Ввиду очевидной ориентации программы на рабочие станции, и того факта, что большинство современных рабочих станций являются системами с общей памятью (SMP), в качестве технологии распараллеливания была выбрана технология OpenMP.
Расчет начальных (граничных) условий. Одним из ключевых требований при разработке программы являлась возможность задавать граничные условия произвольным выражением (суперпозицией заранее определенных операторов и функций). При этом одно и то же выражение (введенное пользователем) должно быть вычислено для каждой точки границы (характерное число ~200 000 точек). Классический расчет с преобразованием в постфиксную нотацию и дальнейшим расчетом (алгоритм Э. Дейкстры) для такого количества точек занимает слишком много времени (~ 40 с), поэтому после преобразования в постфиксную нотацию, выражение «компилируется» в односвязный список, в который подставляются адреса используемых в выражении переменных и функций (операторы также реализованы в виде функций). Откомпилированное выражение может быть вычислено за один проход по списку и для обращения не требует разрешения имен переменных и функций. В результате характерное время вычисления граничного условия уменьшается на три порядка (до 0,04 с).
Программный интерфейс. Интерфейс программы позволяет задавать параметры моделирования, в том числе размеры расчетных сеток, шаги по пространственным переменным и времени, параметры Г и N.
Для задания начального условия пользователю предлагается ввести произвольное выражение с использованием предопределенных операторов, функций и скобок для задания порядка выполнения операций.
Кроме того, пользователь может задавать количество гармоник, которые могут быть интересны для анализа, имеет возможность сохранения пространственных распределений гармоник в файлы для дальнейшего анализа, а также указывать «контрольные точки», в которых нужно сохранять промежуточное решение.
Для удобства пользователя в программе предусмотрен журнал сообщений, в котором отображается информация о параметрах и ходе вычислительного процесса, а также индикатор выполнения моделирования.
Визуализация и средства анализа решения. В программном комплексе MSB обеспечена возможность графического представления информации о динамике развития нелинейных процессов во времени и пространстве.
Для визуализации начальных условий и результата решения задачи в программе MSB используется библиотека с открытым исходным кодом MathGL (mathgl.sourceforge.net). Большой набор инструментов и возможностей MathGL позволяет отображать двумерные и трехмерные поверхности с использованием освещения, строить изолинии, линии тока и т.д.
Интерфейсы визуализации программы позволяют строить сечения, отображать изолинии и управлять их количеством, осуществлять трехмерное вращение изображений относительно осей. Построенное изображение может быть сохранено в один из поддерживаемых форматов графических файлов (PNG, GIF, JPEG).
Для больших размеров сеток визуализация может занимать достаточно большое время, поэтому все окна приложения, в которых выполняется визуализация, вынесены в отдельные потоки.
Верификация результатов работы программного комплекса. Для проверки корректности работы программного комплекса (MSB) выполнялись расчеты при стандартных начальных условиях в виде гауссова звукового пучка с гармонической временной зависимостью при различных значениях внешних параметров N и Г модельного уравнения ХЗК. Такая методика проверки позволяет выполнить ее для всех предельных случаев соотношения масштабов проявления нелинейных (Lp), дифракционных (Ld) и диссипативных (¿з) процессов, сравнивая результаты численных расчетов с известными точными решениями и аналитическими результатами, полученными асимптотическими методами, а также с результатами ранее выполненных численных экспериментов других авторов [15].
В работе выполнялся качественный и количественный анализ результатов расчетов из различных источников. Для количественной оценки различий в поведении волнового профиля использовалась величина среднего квадратического отклонения с, выраженного в процентах по отношению к амплитуде исходной волны.
Для одномерного случая существуют решения Римана (уравнения простых волн) и точное решение уравнения Бюргерса, записанное в интегральной форме [16]. Как известно, уравнение ХЗК переходит в уравнения Римана и Бюргерса при следующих значениях параметров N и Г. Для уравнения простых волн (N —> 0,r —> 0); для уравнения Бюргерса (N — 0).
В работе был выполнен анализ формы волнового профиля на оси звукового пучка при N<<1, Г<<1, что соответствует преобладающему влиянию нелинейных процессов над дифракционными и диссипативными. Сравним ее с волновым профилем, являющимся решением уравнения простых волн для первоначально гармонической волны. Сравнение результатов работы программы MSB с точным аналитическим решением уравнения простых волн выполнялось при следующих значениях параметров модели Y = 0,001 и N = 0,001 при R = 0. Результаты сравнения представлены на рис. 2.
Видно, что имеет место не только качественное, но и хорошее количественное совпадение вплоть до области формирования ударного фронта.
Другим предельным случаем в модели распространения звуковых пучков конечной амплитуды является преобладание дифракционных процессов над нелинейными и диссипативными (N>>1). В этом случае для оценки погрешности работы программного комплекса MSB использовалось точное решение параболического уравнения теории дифракции для Гауссова пучка (нелинейность отсутствует):
ika2 A , ikr2
A(r, x) =--exp(---—). (13)
2 x + ika 2 x + ika
Стандартная запись решения (1) имеет вид
ka2Ai ^ k2a4r2 _ 2xkr2 2,, ,4
A(r, x) = —-2—?2 4 1/2 exp( —;—-TT - i (——-rT - arctg (—j-))), (14) (4x + k a ) 4x + k a 4x + k a ka
L, x t ka2 r Учитывая, что N = ——, z =_, Là =-, R = — получим:
lp lp 2 a
1
R2
R2 Nz
A(r > z) = t2 2 1Л1/2 exp(- " -i ( "T , - «rcig (Nz))). (15) (N z +1)1/2 N z +1 (Nz)2 +1
Рис. 2. Сравнение плосковолновых профилей
На рис. 3, 4 представлены осевые (R=0) распределения амплитуды и фазы основной гармоники в Гауссовом пучке соответственно. Сравниваемые зависимости получены путем расчета по формуле (15) и при помощи программного комплекса MSB (N=10).
Как видно из рис. 3 и 4 сравниваемые графики практически совпадают. Количественное значение погрешности для амплитуды и фазы составляет менее 0,5 %.
На следующем этапе верификации программного комплекса выполнялась проверка работы программы для промежуточных значений внешних параметров уравнения ХЗК. Для этого использовались данные вычислительных экспериментов из монографии [15].
Рис. 3. Сравнение пространственных распределений амплитуды основной гармоники при R=0 (N=10)
Рис. 4. Сравнение пространственных распределений фазы основной гармоники при R=0 (N=10)
Результаты сравнения волновых профилей при соответствующих значениях параметров N, Z, Г и данных [15] представлены на рис. 5.
Так как в каждом из предельных случаев программа MSB дает погрешность всего лишь в пределах 2 %, значительное расхождение с источником [15] связано с сравнительно небольшим количеством расчетных узлов, что видимо и является причиной низкой точности результатов [15].
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Максимов И.И., Савицкий О.А., Сахаров В.Л. Параметрический гидролокатор ближнего действия для охраны мелководных акваторий и узких каналов // Труды IX Всероссийской конференции «Прикладные технологии гидроакустики и гидрофизики». - СПб.: Наука, 2008. - C. 76-79.
2. Нагучев Д.Ш., Савицкий О.А., Сахаров В.Л. Предпосылки и концепция создания современных параметрических профилографов в ОКБ «РИТМ» ЮФУ // Известия ЮФУ. Технические науки.- 2008. - №. 12 (89). - С. 89-94.
3. Савицкий О.А. Взаимодействие волн с различными временными масштабами // Известия ЮФУ. Технические науки. - 2009. - №. 8 (97). - С. 89-93.
4. Савицкий О.А., Чистякова Т.А. Сжатие и декомпрессия импульсов при взаимодействии с низкочастотными волнами конечной амплитуды в звуковых пучках // Известия ЮФУ. Технические науки. - 2011. - № 8 (121). С. 122-128.
5. Гаврилов А.М., Савицкий О.А. Фазозависимые взаимодействия акустических волн.
- Таганрог: Изд-во ТТИ ЮФУ, 2010. - 362 с.
6. Гаврилов А.М., Савицкий О.А. Акустический импульсный локатор // Патент на изобретение RUS 2050558.
7. Гаврилов А.М., Савицкий О.А. Устройство для обнаружения и классификации объектов по акустической жесткости // Патент на изобретение RUS 2006877.
8. Gavrilov A.M., Germanenko O.N., Savitskij O.A. About one possibility of using the second harmonic for measurement of nonlinear media parameters // Акустический журнал. - 1995.
- Т. 41, № 3. - С. 500-501.
9. Гаврилов А.М., Савицкий О.А. Активное подавление нелинейного поглощения звука в квадратично-нелинейных средах без дисперсии // Акустический журнал. - 1997. - Т. 43, № 1. - С. 42.
10. Gavrilov A.M., Savitskij O.A. Absorption of sound by sound at degenerative interaction of acoustical waves // 16 th International Symposium on Nonlinear Acoustics. Nonlinear Acoustics at the Beginning of the 21 st Century. - 2002. - Т. 2. - С. 1043-1046.
11. Гаврилов А.М., Германенко О.Н., Савицкий О.А. Способ активного звукогашения // Патент на изобретение RUS 2185666 20.10.1999.
12. Савицкий О.А. Программа математического моделирования распространения звуковых пучков в квадратично-нелинейной среде без дисперсии (MSB) // Свидетельство о государственной регистрации программы для ЭВМ №2012614679. - 25.05.2012.
13. Савицкий О.А., Чистякова Т.А. Математическая модель распространения ультразвуковых пучков высокой интенсивности // Известия ЮФУ. Технические науки.. - 2010.
- №. 6 (107). - С. 168-174.
14. Чистякова Т.А., Савицкий О.А. Схема расщепления по физическим процессам для уравнения Хохлова-Заболоцкой-Кузнецова // Альманах современной науки и образования.
- 2008. - № 1.
15. Бахвалов Н.С., Я.М. Жилейкин, Е.А. Заболотская. Нелинейная теория звуковых пучков.
- М.: Наука, 1982. - 176 с.
16. Виноградова М.Б., Руденко О.В., Сухоруков А.П. Теория волн: Учебное пособие.
- Наука, 1990.
Статью рекомендовал к опубликованию д.т.н., профессор С.П. Тарасов.
Лапин Дмитрий Вадимович - Федеральное государственное автономное образовательное учреждение высшего профессионального образования «Южный федеральный университет»; e-mail: dmitri.lapin@gmail.com; 347928, г. Таганрог, пер. Некрасовский, 44; тел.: 88634371606; отдел математического моделирования и вычислительного эксперимента; начальник отдела.
Тетеревлёв Игорь Николаевич - e-mail: tin146@mail.ru; г. Таганрог, ул. Зои Космодемьянской, 9 Е; тел.: +79185199182; аспирант.
Савицкий Олег Анатольевич - e-mail: osav66@mail.ru; 347932, г. Таганрог, ул. Ломоносова, 57/1, кв. 57; тел.: 88634371741; кафедра высшей математики; доцент.
Чистякова Татьяна Алексеевна - e-mail: a_tanya84@mail.ru; 347928, г. Таганрог, пер. Некрасовский, 44; тел.: 88634371606; кафедра высшей математики; доцент.
Lapin Dmitriy Vadimovich - Federal State-Owned Autonomy Educational Establishment of Higher Vocational Education "Southern Federal University"; e-mail: dmitri.lapin@gmail.com; 44, Nekrasovskiy, Taganrog, 347928, Russia; phone: +78634371606; the department of mathematical modeling and computer simulation; head of the department.
Teterevlev Igor Nikolaevich - e-mail: tin146@mail.ru; 9 E, Zoya Kosmodemyanskaya street, Taganrog, Russia; phone: +79185199182; postgraduate student.
Savitskiy Oleg Anatoljevich - e-mail: osav66@mail.ru; 57/1, Lomonosov street, apr. 57, Taganrog, 347932, Russia; phone: +78634371741; the department of higher mathematics; associate professor.
Chistyakova Tatyana Alexeevna - e-mail: a_tanya84@mail.ru; 44, Nekrasovskiy, Taganrog, 347928, Russia; phone: +78634371606; the department of higher mathematics; associate professor.
