Нелинейные процессы в ультразвуковых пучках вследствие отсутствия физической дисперсии в большинстве звукопрозрачных сред представляют собой сложные пространственно-временные явления, описываемые квазилинейными уравнениями со степенным характером нелинейных членов. Оптимальной возможностью изучения нелинейных волновых эффектов является применение методов математического моделирования.
Основными физическими процессами, сопровождающими распространение волновых пучков большой амплитуды, являются дифракционные, диссипативные, нелинейные. Влиянием дисперсии скорости звука в среде в большинстве практически важных случаев можно пренебречь. Указанные процессы могут быть учтены в рамках квазилинейного уравнения Хохлова-Заболотской-Кузнецова [5, 6]. В последнее время интерес к изучению сложных пространственных нелинейных волновых процессов связан с широкими перспективами их практического использования [3, 4]. При этом, как правило, результаты по существующим моделям носят частный характер и встречаются в научной литературе по нелинейной акустике [1, 2].
Рассмотрим кратко основные особенности математических моделей для решения уравнения нелинейных звуковых пучков. Бахвалов Н.С. с соавторами разработал численный метод для решения уравнения ХЗК на основе схемы расщепления по физическим процессам для осесимметричных источников. В предложенном алгоритме дифракционные и диссипативные процессы вычислялись в частотной области, а нелинейный член − во временной области, при этом использовался метод Годунова со специальной разностной схемой для решения уравнений в частных производных [1]. В 1984 году Аанонсен с соавторами [13] разработал алгоритм решения ХЗК в частотной области. Алгоритм использует преобразование Фурье для конвертации ХЗК в набор сдвоенных дифференциальных уравнений в частных производных. Позже Бейкер с соавторами [14] обобщил частотный алгоритм Аанонсена для неосесимметричных источников и обобщил метод для импульсных волн (код Бергена). В 2005 г. Янг и соавторы [16] обобщили Техасский код на случай неосесимметричных источников, при этом использовалась схема разложения по поперечным пространственным направлениям. В работе [15] был предложен подход для аппроксимации дифракционного члена без разложения его по боковым осям х и y.
Важным недостатком разновидностей Техасского кода является то, что при решении дифракционной задачи направление счета совпадает с направлением роста временной переменной, при этом используются граничные условия первого рода для 
начального и конечного момента времени, такие же, как и при расчете диссипации. Вследствие этого в указанных моделях выполнение правого граничного условия по времени (для конечного момента времени) при расчете дифракционной задачи не гарантируется.
Постановка задачи
Для описания распространения звуковых пучков конечной амплитуды в нелинейно-диссипативной среде использовано уравнение Хохлова-Заболотской-Кузнецова:
(1)
с начальным и граничными условиями:
, (2)
условие периодичности сигнала:
, (3)
условие отсутствия источников энергии в бесконечно удаленной точке:
(4)
где 
– величина скорости частиц среды, Г – диссипативный параметр, θ – время в сопровождающей системе координат, z – нормированное расстояние, N – параметр уравнения, характеризующий соотношение нелинейной и дифракционной длин волны, 
поперечный лапласиан в декартовой системе координат
.
Рис. 1. Расчетная область
Для решения поставленной задачи в работе использованы методы расщепления по физическим процессам, интегро-интерполяционный метод, экономичные прямые методы решения сеточных уравнений, хорошо зарекомендовавшие себя при решении задач гидродинамики [10], метод решения комплексной матрицы специального вида [7].
Программная реализация
Разработка программного комплекса, реализующего построенную математическую модель, предполагала достижение следующих основных целей: создание мощного инструмента для проведения исследований в области теории нелинейных волновых процессов; разработка современных программных средств для выполнения инженерных расчетов при создании новых образцов техники, основанной на принципах нелинейной акустики. Программный комплекс реализует два основных подхода к описанию нелинейных волновых процессов – полевой и спектральный. Предложенные программно-алгоритмические решения позволяют выполнять расчеты поля скорости частиц среды в акустической волне и порядка нескольких тысяч гармонических составляющих нелинейно искажающегося временного профиля. В структуру программного комплекса входят следующие блоки: блок управления, в котором осуществляются следующие действия: выделение памяти, ввод начальных условий, вывод данных; блок расчета скорости частиц среды с учетом диссипации и нелинейности процесса распространения волновых пучков; блок построения сеточных уравнений для расчета скорости частиц среды; блок дискретного преобразования Фурье, в которых для прямого и обратного направлений преобразований для каждого сечения 
как для вектора применяется быстрое преобразование Фурье; блок расчета скорости частиц среды с учетом дисперсии; блок расчета комплексной СЛАУ специального вида; блок расчета СЛАУ методом циклической прогонки [9].
Схема алгоритма программы представлена на рис. 2. Для тестирования программного комплекса использовалась модельная задача распространения первоначально синусоидального звукового пучка конечной амплитуды с гауссовым поперечным распределением.
Рис. 2. Структурная схема алгоритма программы
Такая задача хорошо исследована и может считаться эталонной для проверки правильности работы построенной математической модели и ее программной реализации. При N>>1 приведенная задача имеет известное аналитическое решение [4].
Результаты численных экспериментов для гауссового пучка
Рассмотрим распространение звукового пучка с начальным значением скорости частиц среды:
.
Моделирование производилось с параметрами: M = 4096, a = 4,0, b = 4,0, Г = 0,001, N = 0,4, Px = Py = 101. На рис. 3 представлена зависимость скорости распространения волны от двух переменных 
при фиксированных значениях z: z = 0,0, 0,5, 1,0, 2,0. Вертикальная ось соответствует r, горизонтальная – θ. С ростом z происходит расширение пучка вдоль r (диссипация энергии), положительный и отрицательный фронты сближаются (эффекты нелинейности).
Рис. 3. Зависимость функции v при y = 0 на расстоянии z, равном: а) z = 0, б) z = 0,5, в) z = 1, г) z = 2
Рис. 4. Зависимость первой (синяя), второй (зеленая) и третьей (красная) гармоник от расстояния z на оси симметрии
С дальнейшим распространением вдоль оси z отрицательный волновой фронт «набегает» на положительный. С увеличением z энергия перераспределяется между гармониками из-за нелинейности процесса, степень нелинейности регулируется параметром Г. На рис. 4 представлены зависимости амплитуд первых трех гармоник от z.
Пример применения программного комплекса для исследования сложных нелинейных пространственных эффектов
Построенная математическая модель позволяет значительно расширить возможности моделирования различных физических ситуаций при изучении закономерностей распространения нелинейных волн при их фокусировке. В частности, например, можно выполнять моделирование сходящихся звуковых пучков не только в точку (такая возможность имелась в модели [5]), но и имеющих цилиндрический характер фокуса.
Моделируем распространение звукового пучка с начальным условием:
.
Вычисление выполнено при данных параметрах модели: Г = 0,001, N = 0,2, Nx = Ny = 81, a = 8,0, b = 8,0. Рис. 5 иллюстрирует зависимости скорости распространения частиц среды от θ на оси симметрии (x = 0, y = 0) при z: z = 0,0, 0,5, 1,0, 2,0, 4,0. В точке z ≈ 1,0 происходит фокусировка звукового пучка.
Рис. 5. Скорость движения частиц среды в зависимости от θ на оси симметрии (x = 0, y = 0) при различных z: z = 0,0, 0,5, 1,0, 2,0, 4,0
Рис. 6. Зависимость поля колебательной скорости от x и времени θ в сечении y = 0 при различных фиксированных z: z = 0,0, 0,5, 1,0, 2,0, 4,0
С расстоянием происходит значительное увеличение градиентов поля и концентрация энергии в фокальной области x = 0. Соответствующее распределение амплитуды положительной и отрицательной части временного профиля представлено на рис. 7.
Рис. 7. Распределение амплитуды положительной и отрицательной части временного профиля
Верификация результатов работы программного комплекса
Для проверки корректности работы разработанного программного комплекса выполнялись расчеты при стандартных начальных условиях в виде гауссова звукового пучка с гармонической временной зависимостью при различных значениях внешних параметров N и Г модельного уравнения ХЗК. Такая методика проверки позволяет выполнить ее для всех предельных случаев соотношения масштабов проявления нелинейных, дифракционных и диссипативных процессов. Проверка работы программного комплекса выполнялась для промежуточных значений внешних параметров уравнения ХЗК. Для этого использовались данные вычислительных экспериментов из [5]. Результаты наложения волновых профилей при соответствующих значениях параметров N, Z, Г представлены на рис. 8.
Рис. 8. Результаты наложения волновых профилей: 1 – результаты [5], 2 – результаты экспериментов, основанных на данной работе
Так как в каждом из предельных случаев разработанное программное обеспечение дает погрешность всего лишь в пределах 2 %, то можно предположить, что значительное расхождение с источником [5] не связаны с ошибочностью результатов работы разработанного программного комплекса.
Выводы
В работе предложен метод для решения систем линейных алгебраических уравнений с комплексной матрицей специального вида. Данный метод позволил модифицировать ранее разработанную конечно-разностную модель решения уравнения Хохлова-Заболотской-Кузнецова, значительно расширив возможности моделирования различных физических явлений при изучении закономерностей распространения нелинейных волн в квадратично-нелинейных средах, не обладающих физической дисперсией. Разработанный программный комплекс может быть использован для выполнения инженерных расчетов при подготовке новой ультразвуковой аппаратуры, работающей на принципах нелинейной акустики.
Работа выполнена при финансовой поддержке РФФИ по проектам № 16-37-00129, № 15-01-08619, № 15-07-08626.