ТРЕХМЕРНОЕ МОДЕЛИРОВАНИЕ АКУСТИЧЕСКОГО ПОЛЯ МЕТОДОМ КОНЕЧНЫХ РАЗНОСТЕЙ ВО ВРЕМЕННОЙ ОБЛАСТИ
,
Тихоокеанский государственный университет, г. Хабаровск
E-mail: *****@***edu. ru
В статье рассматривается применение метода конечных разностей во временной области для моделирования распространения ультразвука в твердом теле. Расчет проводился в трехмерном пространстве. Алгоритм позволяет учитывать размер и границы реальных объектов. Программное обеспечение может быть использовано при решении задач в области неразрушающего контроля для прогнозирования и управления волновым фронтом.
Технологии ультразвукового сканирования приобрели огромную популярность в области неразрушающего контроля и быстро развиваются. Применение компьютерной техники и современного программного обеспечения позволяет значительно количественно и качественно расширить возможность существующих методов неразрушающего контроля и создание новых. Огромную популярность в неразрушающем контроле получили численные методы решения волновых уравнений. Для решения дифференциальных уравнений применяют метод конечных разностей во временной области (КРВО) [1,2]. В статьях [3,4] автор рассматривает моделирование распространения акустического сигнала в двумерном случае, используя метод КРВО. Целью настоящей работы является разработка эффективных алгоритмов моделирования распространения ультразвукового сигнала в трехмерном случае с учетом свойств и геометрии материала, внешних и внутренних дефектов и создание соответствующего программного обеспечения.
Авторами разработан алгоритм, и на основе его реализован программный комплекс для построения акустического поля с учетом геометрии тела, материала и дефектов в трехмерном пространстве. В основе алгоритма расчета используется модель Кельвина вязкоупругой среды [3]:
(1)
где: σij – тензор напряжений [Па], εij – тензор деформации, Сijkl – тензор упругости [Па], ηijkl – тензор вязкости [Па*с], t –время [с].
А также второй закон Ньютона для сплошной среды:
(2)
где: ρ – плотность материала, vi – компоненты вектора скорости, Fi – компоненты вектора сторонних сил, воздействующих на среду.
Все деформации в алгоритме считаются линейными, поэтому тензор деформации εij имеет следующую зависимость от вектора смещения среды Ui:
(3)
Для расчета с учетом геометрии тела, материала и внутренних дефектов авторами реализован программный комплекс в среде программирования Visual Studio на языке C#. Более подробно алгоритм представлен авторами в статьях [5, 6].
Для уменьшения времени расчета была использована технология CUDA для расчета на GPU, что позволило сократить время вычислений в 20 раз в сравнении с расчетом на центральном процессоре. Для реализации параллельных вычислений на ГПУ массив входных данных разбивается на ряд блоков. Блоки образуют сетку, которая имеет размер M × N. Для обеспечения корректности вычислений имеется пересечение между соседними блоками, что позволяет осуществить обмен граничными данными между блоками. Обмен данными между блоками обычно производится при помощи использования разделяемой памяти (использования разделяемой памяти удается избежать в том случае, когда одна нить обрабатывает данные в одной ячейке сетки).
Для проверки используемого алгоритма и созданного на его основе ПО было проведено несколько численных экспериментов расчета акустических полей в трехмерной среде с различным расположение излучателей и наличием дефектов:
1. Расчет распространения ультразвукового сигнала в параллелепипеде размером 15×15×15 миллиметров. Точечный импульсный источник возмущения S1 расположим в середине ребра параллелепипеда (рис. 1а). Для наглядности результат расчета показан в выбранном сечении d. Так как амплитуды разных типов волн сильно отличаются и имеются большие изменения амплитуды в пределах фронта, то при выводе проводится выравнивание амплитуд, как 
. На рисунке 1б показан результат расчета в плоскости d в момент времени 1 мкс. Частота дискретизации Δx1 = Δx2 = Δx3 = 0,050 мм, временной шаг Δt = 0,5 нс.
Полученную картину можно интерпретировать следующим образом: Источник силы возбуждает продольные волны 1 (скорость распространения из расчета: расстояние, пройденное волной относительно времени кадра v = 5970 м/с) и поперечные волны 2 (скорость распространения v = 3050 м/c). На поверхности образуются волны Рэлея 3 и головная волна 4. Головная волна генерирует фронт поперечных волн 5, который, согласно теории, распространяется под третьим критическим углом, равным для стали 33 градусам. Полученные результаты для трехмерного случая согласуются с результатами двумерного [7].
Рис 1. а) Стальной параллелепипед размером 15x15x15 мм. S1 – источник возмущения, d – секущая плоскость. б) Волны точечного источника, расположенного на поверхности, стали. (λ=109ГПа, μ= 83 ГПа, ρ = 7.8 КГ/дм3. Длительность импульса 0,01 мкс. Время τ = 1 мкс.)
2. Расчет волнового фронта при наличии диагональной границы раздела сред интерференционная картина от точечного излучателя S3 на рисунке 2. Среда 1 представляет собой боросиликатное стекло, а среда 2 – плексиглас. Плотность боросиликатного стекла – 2230 кг/м3, скорость продольной волны – 5640 м/c. У плексигласа плотность – 1180 кг/м3, скорость продольной волны 2600 м/c. Из рисунка видно, что численное моделирование хорошо согласуется с расчетами на основе законов геометрической оптики. Волновой фронт частично отражается от границы раздела сред, частично преломляется.
Рис 2. Волны источника в плоскости d (Длительность импульса 0,01 мкс.) Частота дискретизации Δx1 = Δx2 = Δx3 = 0,050 мм2, Δt = 0,5 нс. Время a) τ = 0,5 мкс b) τ = 1 мкс
3. Распространение акустического сигнала при наличии дефекта. Расположим импульсный точечный источник S, как показано на рисунке 3. В плоскости w располагается дефект n (рис. 4). Размер полости: 10x10x0,1 мм.
Рис 3. Стальной параллелепипед размером 15x15x15 мм. S – источник возмущения. d, w – секущие плоскости, n – дефект в виде плоскости, заполненной вакуумом.
На рис. 4 приведены результаты расчета для прямоугольного излучателя S размером 3x1 мм в плоскостях d (рис.4а и b) и плоскости, перпендикулярной w и d (рис.4с и d). Для двух моментов времени Δt = 0,5 мкс и 1 мкс. В любой точке среды в данный момент времени колебания представляют собой сумму волн, пришедших от различных элементов поверхности излучателя. Так как длина пробега от элемента поверхности до излучателя разная, то образуется сложный волновой фронт колебаний. Происходит отражение от дефекта, а также дифракция на краях дефекта.
Рис 4. Волны прямоугольного источника размером 3x1 мм в плоскости d (a, b), и в плоскости перпендикулярной d и w (c, d). (λ=109 ГПа, μ= 83 ГПа, ρ = 7.8 кг/дм3. Длительность импульса 0,01 мкс.) Частота дискретизации Δx1 = Δx2 = Δx3 = 0,05 мм, Δt = 0,5 нс. Время a) 0,5 мкс b) 1 мкс c) 0,5 мкс d) 1 мкс.
Сравнение рисунков 4а и 4b с рисунками 4с и 4d показывает влияние размеров излучателя по разным направлениям на акустическое поле и свидетельствует о достоверности полученных результатов для трехмерных расчетов.
В программном комплексе реализована возможность установки виртуальных приемников сигнала. В точке установки приемника записываются все значения вектора скорости и тензора деформации для каждого момента времени в файл формата. txt для последующей обработки (получение временной зависимости тензора деформации, расчет коэффициента затухания сигнала).
Заключение
Таким образом, на основе модели Кельвина для вязкоупругой среды и метода КРВО с использованием технологии CUDA для расчета на GPU, создан программный комплекс, позволяющий в реальном времени проводить расчеты акустического поля в трехмерном пространстве. Программное обеспечение позволяет в реальном времени проводить построение волнового фронта с учетом материала, геометрии тела и наличии различных дефектов. Проведена качественная верификация расчетов волнового фронта, показавшая достоверность получаемых результатов для различных моделей. Созданное ПО позволяет рассчитывать акустическое поле, создаваемое системой излучателей (фазированной решеткой). В дальнейшем планируется провести исследование особенностей распространения акустического сигнала в деталях сложной формы и композиционных материалах.
Л И Т Е Р А Т У Р А
1. R. Madariaga. Dynamic of an expanding circular fault. Bull. Seism. Soc. Am., vol. 66(3), pp. 639-666, 1976.
2. J. G. Maloney and K. E. Cummings, “Adaptation of FDTD techniques to acoustic modeling”, 11th Annual Review of Progress in Applied Computational Electromagnetics, Vol. 2, pp. 724-731, Monterey, CA, March 1995.
3. . Решение динамических задач акустики методом конечных разностей во временной области. Основные соотношения. Анализ погрешностей // Дефектоскопия. – 2005. – №3. – С. 12–26.
4. J.Virieux. SH-wave propagation in heterogeneous media: Velocity - stress difference method. Geophysics, vol.49 (11), pp.1933-1942, November 1984.
5. Rimlyand V. I., Avdeev D. A. Simulations of ultrasound propogation in solids of various shapes. // ISCV22. – Florence, Italy - 2015.
6. , Моделирование распространения ультразвука в твердых телах различной формы // Информатика и системы управления. – 2015. –№2. – С. 15 – 22
7. , . Моделирование распространения ультразвука в твердом теле // Физика: фундаментальные и прикладные исследования, образование: Материалы Всероссийской молодежной конференции. – Благовещенск: Амурский гос. Ун-т, 2014. – С. 3 – 6.
