Математическое моделирование излучения акустической антенной на многопроцессорной системе*

Abstract

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

-

2. Постановка задачи ptt  c2p  f , (1) удовлетворяющего начальным условиям:

p(x, y, 0) 0 (x, y) , pt(x, y, 0) 1(x, y) (2) и граничным условиям:

p(x, y, t)  0 , при (x, y)  , (3) где f – функция-источник,  – граница расчетной области.

Наиболее эффективными методами для решения подобного вида задач являются сеточные методы. 3. Решение задачи t  {t n  nht , 0  n  Nt 1, lt  ht  Nt 1},  x  {xi  ihx , 0  i  Nx 1, lx  hx  Nx 1},  y  {y j  jhy , 0  j  N y 1, ly  hy  N y 1} , где n, i, j – индексы по временной координате и пространственным координатным направлениям Ox , Oy соответственно; ht , hx , hy – шаги по временной координате и пространственным координатным направлениям; Nt , Nx , N y – количество узлов по временной координате и пространственным координатным направлениям; lt , lx , ly – длина расчетной области по временной координате и пространственным координатным направлениям.

Для получения дискретной модели использован интегро-интерполяционный метод [ 4 ]. Разностная схема, аппроксимирующая уравнение (1): pin,j1  2 pin, j  pn1

ht2 i, j  c2 pi1, j  2hpx2i, j  pi1, j  c2 pi, j1  2hpy2i, j  pi, j1  fi,nj , (4) где pi, j   1 pin,j1  1 1  2  pin, j  2 pin, j1 ; 1, 2 – веса схемы [ 5 ].

Для улучшения «гладкости» решения сеточного решения будем предполагать, что ячейки заполнены не полностью [ 6 ]. Областью xy будем называть заполненную часть области Dxy x  xi1/2 , xi1/2 , y   y j1/2 , y j1/2  .

Также введем обозначения для следующих областей: D x  xi , xi1/2 , y   y j1/2 , y j1/2  , D2 x  xi1/2 , xi , y   y j1/2 , y j1/2  , 1 D3 x  xi1/2 , xi1/2 , y   y j , y j1/2  , D4 x  xi1/2 , xi1/2 , y   y j1/2 , y j  . Коэффициенты заполненности q0, q1, q2, q3, q4 для областей Dtxy , D1, D2 , D3, D4 вводятся следующим образом:

S S q0  Dxy , qi  Di , i  1, 4 ,

Sxy Si где S – площадь соответствующей части области; i – заполненная часть области Di .

Дискретный аналог уравнения (1) в случае граничных условий в форме Дирихле ( p  0 ) примет вид: pin,2j1  q0,i, j ht 2 pin, j  pn1

i, j  c2 min  q1,i, j , q2,i, j  pi1, j  2 pi, j  pi1, j  ht2 hx2 Дискретный аналог уравнения (1) в случае граничных условий в форме Неймона ( pn  0 ) (5) примет вид: q0,i, j

pin,j1  2hpt2in, j  pin,j1  c2q1,i, j pi1, hjx2 pi, j  c2q2,i, j pi, j hx2pi1, j  c2q3,i, j pi, j1hy2 pi, j  c2q4,i, j pi, j hy2pi, j1  q0,i, j fi,nj . 4. Решение модельной задачи

Схему итерационного двухслойного модифицированного попеременно-треугольного метода (МПТМ) запишем в виде [ 7–9 ]:

xn+1  xn  n1wn ,  D  R1  D1  D  R2  wn  r n , r n  Axn  f , где xn – вектор решения; w n – вектор поправки; А – оператор сеточного уравнения; D – диагональная часть оператора A;  – итерационный параметр; R1, R2 – верхняя и нижняя треугольные части оператора A; r n – вектор невязки; f – правая часть сеточного уравнения.

Для параллельной реализации адаптивного МПТМ были использованы методы декомпозиции области по одному направлению. Наиболее трудоемким расчетом с точки зрения построения параллельной программной реализации является расчет вектора поправки, который выполняется в два шага: 1)  D  R1  yn  r n , 2)  D  R2  wn  Dyn .

На первом шаге рассчитываются элементы вспомогательного вектора y n снизу вверх, а затем, зная его, на втором шаге находятся элементы вектора поправки w n сверху вниз. Схемы расчета вспомогательного вектора и вектора поправки представлены на рис. 1 (стрелками показаны направления счета и передачи, штриховкой отмечены расчетные узлы).

Рис. 1. Схемы расчета вспомогательного вектора yn (а) и вектора поправки wn (б) В схеме для расчета вектора y n только первый процессор не требует дополнительной информации и может независимо от других процессоров посчитать свою часть области, остальные процессоры ждут от предыдущего передачи элементов, стоящих в начале строки. Передача по одному элементу не оптимальна, так как появляются временные затраты, связанные с их организацией. Суммарное время на накладные расходы можно уменьшить путем увеличения объема передач. Данные рассуждения используются для расчета вектора поправки. 5. Параллельная реализация

При параллельной реализации использованы методы декомпозиции сеточных областей для вычислительно трудоемких задач, учитывающие архитектуру и параметры многопроцессорной вычислительной системы. Максимальная производительность МВС составляет 18,8 терафлопс. В качестве вычислительных узлов используются 128 однотипных 16-ядерных Blade-серверов HP ProLiant BL685c, каждый из которых оснащен четырьмя четырехъядерными процессорами AMD Opteron 8356 2.3GHz и оперативной памятью в объеме 32 ГБ. В таблице 1 приведены временные затраты выполнения одной итерации МПТМ на различных сетках, а также значения ускорения и эффективности для различного числа вычислительных ядер. Для расчетов использована технология распараллеливания MPI.

Таблица 1. Ускорение и эффективность работы параллельного варианта МПТМ 1 2 4 8 16 32 64 128 256 512 Время Ускорение Эффективность Время Ускорение Эффективность Время Ускорение Эффективность Время Ускорение Эффективность Время Ускорение Эффективность Время Ускорение Эффективность Время Ускорение Эффективность Время Ускорение Эффективность Время Ускорение Эффективность Время Ускорение Эффективность Из таблицы также видно, что для каждой из расчетных сеток ускорение принимает наибольшее значение при определенном значении вычислителей и при дальнейшем увеличении числа вычислительных ядер ускорение только уменьшается. Это связано с временными затратами на обмен данными между вычислителями. Следует отметить, что адаптивный МПТМ нашел свое применение при решении задач аэро-гидродинамики [ 10 ] и транспорта донных материалов [ 11–12 ]. 6. Результаты численных экспериментов

Программный компонент для расчета картин поля излучения акустических волн включает в себя следующие компоненты: блок ввода исходных данных; блок расчета геометрии; блок расчета коэффициентов сеточных уравнений; блок расчета функций правых частей сеточных уравнений; блок перехода на более грубую сетку; блок расчета сеточных уравнений на основе модифицированного попеременно-треугольного метода; блок расчета положения расчетного окна; блок учета граничных условий; блок вывода рассчитываемых функций давления; блок вывода спектра; блок расчета фазы; блок расчета градиента фазы; блок расчета направленности; блок вывода направленности.

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

Рис. 2. Распространение акустических волн от антенн с различными характеристиками направленности На рис. 3 представлена влияние наличия препядствия (неоднородностей) на картину распространения акустической волны.

Рис. 3. Рассеяние акустических волн на препятствии Для изучения отраженных сигналов удобно моделировать короткоимпульсные сигналы. На рис. 4 представлены картины распространения одиночного акустического сигнала и его отраженние от препятствия (неоднородности) в различные моменты времени. Рис. 4. Отраженный акустический сигнал 7. Результаты численных экспериментов

Требуется рассчитать акустические поля антенны с рабочими частотами: 1.5 кГц, 2.0 кГц, 2.5 кГц, 3.0 кГц. Приведенное давление: 10 кПа, 20 кПа, 50 кПа, 100кПа. Рабочие поверхности состоят из пьезокерамики ЦТС-19 и поляризованы радиально, электроды находятся на боковых поверхностях. На рис. 5 приведена геометрия излучающей акустической антенны.

Рис. 5. Геометрия излучающей антенны На рис. 6 приведены результаты расчета распространения колебательных процессов от антенны в случаях: непрерывного сигнала (а) и короткого импульса (б). Рабочая частота 3.0 кГц. Размер окна моделирования составлял 5×12 м. Временной интервал составлял 3 мс. Рис. 6. Результаты расчета распространения колебательных процессов от антенны: непрерывный сигнал (а) и короткий импульс (б) На рис. 7 приведены результаты расчета интенсивности и направления распространения колебательных процессов в ближнем поле антенны и вблизи элементов антенны в случае непрерывного сигнала.

Рис. 7. Результаты расчета направления распространения и интенсивности колебательных процессов вблизи элементов антенны (а) и в ближнем поле антенны (б) Для ускорения работы разработанного программного комплекса при расчете акустических полей в дальнем поле антенны переход на более грубую сетку. Рис. 8. Результаты расчета распространения колебательных процессов от антенны: временной интервал 3 мс (а) и 12 мс (б) На рис. 8(а) представлены результаты расчета акустического сигнала в ближнем поле антенны. Временной интервал составлял 3 мс. После чего осуществлен переход на сетку с 4 раза большими шагами по пространственным координатам. На рис. 8(б) приведены результаты расчета картины распространения звуковой волны через 12 мс после начала работы акустической антенны. Размер окна моделирования составлял 5×12 м на подробной сетке и 20×48 м на грубой сетке, что вполне достаточно для расчета интенсивности и направления распространения колебательных процессов в ближнем поле и на элементах конструкции антенны.

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

Работа посвящена изучению волновых колебаний, а также построению комплекса программ, предназначенного для описания волновых процессов излучения в линейных антенных решетках с изменяемой геометрией. В основе предложенной математической модели лежит неоднородное волновое уравнение с соответствующими начальными и граничными условиями. Для решения поставленной задачи был выбран метод сеток. Дискретная модель была построена при помощи интегро-интерполяционного метода, при этом осуществлялся учет заполненности расчетных ячеек, что гарантировало выполнение основных законов сохранения на дискретном уровне. Полученные сеточные уравнения решены адаптивным модифицированным попеременно-треугольным методом вариационного типа, который имеет наилучшую скорость сходимости в классе двухслойных итерационных методов. Литература 12. Сухинов А.И., Чистяков А.Е., Проценко Е.А. Математическое моделирование транспорта наносов в прибрежных водных системах на многопроцессорной вычислительной системе // Вычислительные методы и программирование: новые вычислительные технологии. 2014. Т. 15. № 4. С. 610–620. Mathematical modeling of acoustic antenna radiation on multiprocessor system Don state technical University1, Kalyaev Scientific Research Institute of Multiprocessor Computer Systems at Southern Federal University2, N.N. Andreyev Acoustic Institute of the

Russian Academy of Sciences3 At present the development of optimal constructions of acoustic antennas is a relevant problem. For this it is necessary to calculate characteristics of antenna arrays. A discrete mathematical model of distribution of acoustic waves by the antennas array was designed in this paper. Grids that takes into account the occupancy of computational cells were used to improve the actual accuracy of solutions. A discrete model is based on the balance method. Results of numerical experiments for different forms of radiated antennas that operated in a given frequency range were given. Grids with dimensions from 105 nodes for twodimensional models to 109 nodes in three-dimensional models are required for numerical solution of these problems with the accuracy, required in designing tasks of acoustic antennas in view of edge effects. The need of multiple calculations and the limit of computing resources and, in particular, the memory of serial computers lead to design effective parallel computing methods for solving these problems on available multiprocessor systems, contained hundreds or thousands cores and available to design organizations.

References

Markov G.T. , Chaplin A.F. Vozbuzhdeniye elektromagnitnykh voln [Excitation of electromagnetic waves] . Moscow, Radio and Communications, 1983 .

universities. 4th ed. and ext.]. Moscow, Science Publishing, 1981 .

3. Sukhinov A.I. , Zuev V.N. , Semenistyy V.V. Uravneniya matematicheskoy fiziki [ Equations of mathematical physics].Taganrog, Publishing of the TRTU , 2005 .

4. Samarskiy A.A. Teoriya raznostnykh skhem [The theory of difference schemes] . Moscow, Science Publishing, 1989 .

5. Sukhinov , A.I. , Chistyakov , A.E. , Shishenya , A.V. Error estimate for diffusion equations solved by schemes with weights . Mathematical Models and Computer Simulations . 1 May 2014 . Vol. 6 , No. 3. P. 324 - 331 .

6. Sukhinov A.I. , Chistyakov A.E. , Fomenko N.A. Metodika postroyeniya raznostnykh skhem dlya zadachi diffuzii-konvektsii-reaktsii, uchityvayushchikh stepen' zapolnennosti kontrol'nykh yacheyek [The method of constructing difference schemes for problems of diffusion-convectionreaction, taking into account the degree of filling of the control cells]. Izvestiya Yuzhnogo federal'nogo universiteta . Tekhnicheskiye nauki [Izvestiya of SFEDU. Engineering sciences.] . 2013 . No. 4. P. 87 - 98 .

7. Samarski A.A. , Nikolayev E.S. Metody resheniya setochnykh uravneniy [Methods for solving grid equations] . Moscow, Science Publishing. 1978 . 588 P .

Differentsial'nyye uravneniya [Differential Equations] . 2004 . Vol. 40 , No. 7 . 953 P.

9. Sukhinov , A.I. , Chistyakov , A.E. Adaptive modified alternating triangular iterative method for solving grid equations with a non-self-adjoint operator . Mathematical Models and Computer Simulations . 1 July 2012 .Vol. 4 , No. 4 P. 398 - 409 .

10. Sukhinov A.I. , Chistyakov A.E. , Timofeeva E.F. , Shishenya A .V. Mathematical model for calculating coastal wave processes . Mathematical Models and Computer Simulations . 1 April 2013 . Vol. 5 . No. 2, P. 122- 129 .

11. Sukhinov A.I. , Chistyakov A.E. , Protsenko E.A. Mathematical modeling of sediment transport in the coastal zone of shallow reservoirs . Mathematical Models and Computer Simulations . 1 July 2014 . Vol. 6 , No. 4. P. 351 - 363 .

12. Sukhinov A.I. , Chistyakov A.E. , Protsenko E.A. Matematicheskoye modelirovaniye transporta nanosov v pribrezhnykh vodnykh sistemakh na mnogoprotsessornoy vychislitel'noy sisteme [Mathematical modeling of sediment transport in coastal water systems on multiprocessor computer system] . Vychislitel'nyye metody i programmirovaniye [Numerical Methods and Programming] . 2014 . Vol. 15 . P. 610 - 620 .