Numerical modelling of acoustical processes by particles interaction.pdf Введение Численное моделирование акустических процессов широко используется для прогнозирования поведения ультразвуковых вибрационных систем, акустической томографии и ультразвуковой связи. Для моделирования акустических процессов широкое распространение получили классические методы численного моделирования, такие, как метод конечных разностей во временной области [1-3], метод конечных элементов [4-6], метод частиц в ячейках [7, 8]. Данные методы являются сеточными и требовательны к вычислительным ресурсам, так как решения приходится искать в узлах многомерной сетки. Также к недостаткам данных методов можно отнести накопление ошибок вычислений при выборе большого шага сетки, а выбор слишком маленького шага приведет к увеличению времени вычисления. Метод дискретных элементов применяется для расчета движения большого количества частиц, таких, как молекулы, песчинки, гравий, галька и прочих гранулированных сред. Преимущество метода заключается в возможности моделирования частиц с несферической поверхностью. В методе дискретных элементов частицы в исходном состоянии имеют начальную скорость. Метод конечных элементов хорошо подходит для линейных задач. В настоящее время метод широко используется как для решения электродинамических, так и акустических задач. Преимущества метода конечных разностей во временной области - это возможность моделирования нестационарных волновых процессов в вязкоупругих, неоднородных, анизотропных и нелинейных средах. В настоящее время для решения акустических задач в различных средах или в средах с деформируемыми границами все чаще находят применение бессеточные методы. Одним из таких методов является метод гидродинамики сглаженных частиц (Smoothed Particle Hydrodynamics (SPH)) [9-13]. Суть метода заключается в том, что моделируемая среда представляется как множество частиц, обладающих такими физическими параметрами, как скорость, плотность, давление. Метод SPH позволяет вычислять производные физических параметров частиц без вычислительной сетки и чаще всего находит применение для моделирования жидкостей. Существуют работы, в которых взаимодействие между частицами зависит от температуры для моделирования твердых тел, где взаимодействия между частицами описываются на основе молекулярной динамики [14, 15]. Разработана модель, которая описывает как жидкое, так и твердое состояние за счет изменения формы потенциальной кривой энергии как функции температуры. Представление среды в виде множества частиц также получило широкое распространение при описании акустических процессов в фононных кристаллах. В методе молекулярной динамики для описания движения и взаимодействия частиц используются законы классической механики. Временная эволюция системы взаимодействующих атомов или частиц отслеживается интегрированием их уравнений движения. Силы межатомного взаимодействия представляются в форме классических потенциальных сил (как градиент потенциальной энергии системы). Существует ряд работ, в которых применяется метод молекулярной динамики и его модифицированное представление для решения различных задач, таких, как расчет дисперсии фононов, для оценки нестабильности движения дислокаций, распространение звуковой волны в газообразном аргоне [16]. Чтобы рассмотреть все разнообразие акустических процессов, мы предлагаем моделировать твердые среды на фундаментальном уровне, представляя среду как трехмерную решетку частиц, с определенной массой динамика которых описывается уравнениями движения Ньютона, по аналогии с методом молекулярной динамики. Взаимодействие между частицами описывается функцией, являющейся зависимостью силы притяжения между частицами от расстояния. Свойства моделируемого материала задаются регулированием значения массы частиц и угла наклона кривой зависимости силы от расстояния вблизи точки равновесия. Силы взаимодействия между частицами характеризуются силовой функцией, которая определяется силой притяжения между частицами на расстоянии. В точке равновесия сила трения исчезает. Угол наклона силовой функции вблизи точки равновесия пропорционален модулю Юнга. В данной работе предлагается метод численного моделирования акустических процессов в твердых телах на основе массива частиц в кубической объемно-центрированной кристаллической решетке. Численное моделирование динамики частиц Предлагается моделировать акустические процессы, представляя твердые тела в виде набора частиц, расположенных в кристаллической решетке. Рассмотрим случай с кубической объемно-центрированной кристаллической решеткой. Частицы в рассматриваемой решетке располагаются периодически, но ближайшие частицы могут располагаться на расстоянии (при рассмотрении по осям рёбер куба) или (при рассмотрении по диагонали куба). Каждая частица имеет собственную массу и силу взаимодействия со своими соседями. Для каждой частицы мы рассчитываем скорость и координаты путем численного интегрирования ускорения частицы по времени. Ускорение частицы рассчитывается по второму закону Ньютона как отношение силы, действующей на частицу к ее массе. Масса частицы и сила взаимодействия определяют тип материала. Сила взаимодействия зависит от расстояния между частицами. Мы применяем разные силовые зависимости для соседних частиц на разных расстояниях: - для частиц на расстоянии d; - для частиц на расстоянии . Для на расстоянии сила равна нулю, для на расстоянии сила тоже равна нулю. Положение частиц стабильно, если нет внешних сил, действующих на них. Следующая ближайшая частица находится на расстоянии . Если расстояние равно или больше , сила взаимодействия считается равной нулю, так как мы учитываем только ближний порядок в расположении частиц. Силовая зависимость пропорционально масштабирована по оси x относительно . Такая зависимость характеризуется следующими формулами: (1) где и - расстояние между частицами; - сила взаимодействия между частицами на расстоянии - сила взаимодействия между частицами на расстоянии ; - коэффициент пропорциональности между силой и расстоянием для линейного приближения, измеренный в Н/м; и были найдены численно для обеспечения изотропии скорости звука. Все частицы ведут себя согласно закону Ньютона. Итерационный процесс вычисления координаты одной частицы выглядит следующим образом: , . (2) Здесь - ускорение частицы; - масса частицы; - сила, действующая на частицу и на соседние частицы; - скорость частицы в текущий момент времени ; - координата частицы; - шаг по времени в численном моделировании. Чтобы избежать самовозбуждения решетки, мы вводим силу «искусственного трения», которая пропорциональна скорости и направлена против нее. В реальных средах сила трения превращает упорядоченное движение частиц в хаотическое тепловое движение. Но для численного моделирования эти хаотические колебания обычно являются препятствием для правильного предсказания динамики частиц и распространения волн. Путем введения искусственного трения мы уменьшаем высокоскоростные хаотические колебания частиц решетки. Сила, действующая на конкретную частицу, рассчитывается как векторная сумма сил взаимодействия с соседними частицами, силы трения и силы внешнего возбуждения. Масса и силовая функция определяют плотность и скорость звука конкретного материала, что позволяет связать моделирование с процессами в реальных средах. Кристаллическая решетка и осевая зависимость силы взаимодействия могут привести к анизотропным эффектам. Для изотропных сред целесообразно использовать объемно-центрированную кубическую решетку. В такой решетке каждая частица имеет 14 ближайших соседей (6 на расстоянии и 8 на расстоянии ), что необходимо учитывать при расчете силы взаимодействия между частицами. Пока амплитуда возбуждения мала, поведение системы будет линейным, поскольку вблизи точки равновесия сила взаимодействия линейна. При увеличении силы возбуждения сдвиги между частицами будут увеличиваться и нелинейная часть функции силового взаимодействия будет оказывать влияние. Весь вычислительный процесс протекает во временной области, что позволяет учесть нелинейные эффекты. В линейном приближении мы можем рассчитать скорость звука в среде с учетом массы частиц и значения . Рассмотрим силы, действующие на одну частицу вдоль оси , и запишем уравнение движения частицы: , (3) где - смещение частицы в момент времени в плоскости ( - координата частицы в невозмущенном состоянии); - сдвиг соседней частицы на расстояния . Приближенно , следовательно, , которая является формой волнового уравнения, где скорость звука рассчитывается через выражение , следовательно, . Для задания плотности материала и скорости звука применялись следующие формулы: ; , где - плотность материала, - скорость звука в материале. Численное моделирование Для проверки возможности предложенного метода для моделирования волновых процессов предлагается рассмотреть фокусировку акустической волны линзой из меди, помещенной в алюминиевую фоновую среду. Если моделируемая плоская волна будет фокусироваться после прохождения через линзу в фокус линзы, это позволит судить о применимости модели к волновым процессам. Мы смоделировали алюминиевый куб с размерами 256×256×256 мм3 из 64×64×128 = = 524288 частиц. Внутрь алюминиевого куба была встроена медная линза с фокусным расстоянием 100 мм. Форма линзы была рассчитана аналитически для расстояния фокусировки линзы в 100 мм от ее вершины. На рис. 1 представлена смоделированная Z-составляющая скорости частиц при возбуждении вершины конструкции плоской волной на частоте 400 кГц. Хорошо виден эффект фокусировки. Но величина поля в точке фокусировки сравнима с полем, прошедшим мимо линзы, что объясняется отражением волн от границ линзы. Рис. 1. Численное моделирование акустической линзы, фокусирующей плоскую волну в точку Было проведено исследование устойчивости численного моделирования для разных временных шагов. Известно, что если шаг по времени слишком велик, ошибка в численном интегрировании дифференциальных уравнений возрастает. С другой стороны, если временной шаг слишком мал, накапливается ошибка из-за суммирования малых и больших чисел, что особенно важно при использовании чисел одинарной точности с плавающей запятой. Существует оптимальное значение шага во времени, когда ошибка численного моделирования наименьшая. Мы провели численное моделирование колебаний в алюминиевом параллелепипеде размерами 32×32×128 мм3, d = = 4 мм, время моделирования составило 20 мкс. Средний квадрат различий в координатах частиц между моделированием с меньшим и последующим большим шагом по времени был рассчитан по формуле , (4) где - координаты частиц в конце численного моделирования для меньшего временного шага; - координаты частиц в конце численного моделирования для следующего большего временного шага. Результаты моделирования для представлены на рис. 2. Согласно графику, представленному на рис. 2, лучшая точность достигается при временном шаге 3 нс для конкретного материала и пространственном шаге . Рис. 2. Исследование устойчивости для разных временных шагов Рассмотрим численную модель распространения акустических волн в твердом теле. В качестве фоновой среды возьмем алюминий, у которого скорость звука составляет 6400 м/с, а плотность - 2712 кг/м3. Моделирование проводилось в области 128×127×4 мм. В центре находился точечный источник и возбуждал поперечные колебания пластины. На центральном сечении расчетной области наблюдались значения поля. Начальное возмущение было направлено вдоль оси Z. Результат моделирования волны после 20 мкс показан на рис. 3. Рис. 3. Распространение сферической волны от точечного источника Волна, распространяющаяся в среде от точечного источника, почти цилиндрическая, что свидетельствует о незначительной анизотропии рассматриваемой кристаллической решетки. Мы провели численное моделирование распространения волн в цилиндрической структуре, состоящей из одной твердой среды, чтобы проверить способность предлагаемой модели прогнозировать резонансные частоты. Цилиндр (волновод) выполнен из титана (скорость 4141 м/с, плотность 4500 кг/м3), длина 80 мм, диаметр 8 мм. В качестве сигнала возбуждения мы применяем короткий импульс. Этот численный тест показал бы спектр резонансных частот моделируемой структуры, который должен соответствовать аналитической оценке полуволнового резонанса на частоте 25881 Гц. Во время моделирования контролировалась координата краевых точек. После численного моделирования мы рассчитали ее спектр . Спектр пространственных колебаний краевой точки волновода представлен на рис. 4. Рис. 4. Спектр импульсного сигнала на конце волновода На частоте 25781 Гц наблюдается резонанс. Эта частота близка к полуволновой резонансной частоте, равной 25881 Гц, согласно теоретической оценке , где = 80 мм - длина волновода. Для проверки возможности численной модели получить распределение амплитуды колебаний волновода более сложной формы было проведено моделирование составного волновода. При численном моделировании в качестве пьезоэлемента использовалась керамика ЦТБС-8 (шириной 20 мм), отражатель длиной 35 мм, болт длиной 78 мм и плавный переход в виде cos-функции длиной 64 мм, при этом общая длина волновода составила 180 мм (рис.5). Рис. 5. Конструкция ультразвукового волновода с cos(x)-переходом При численном моделировании задавалась скорость звука и плотность материала. Скорость звука в алюминии 5050 м/с, плотность 2712 кг/м3. В качестве пьезоэлементов использовалась керамика ЦТБС-8, скорость звука 3200 м/с, плотность 7600 кг/м3. В качестве возбуждающего сигнала применялся короткий импульс. Разными оттенками серого цвета представлены скорости, с которыми колеблются частицы. На конце волновода размещалась мониторинговая частица, которая регистрировала отклонения координаты и скорость. В результате проведенного численного моделирования был получен резонансный частотный спектр моделируемой структуры. Спектр краевых пространственных колебаний представлен на рис. 6. Рис. 6. Конструкция ультразвукового волновода с cos(x)-переходом По результатам численного моделирования можно сделать вывод, что предлагаемая модель предсказывает резонансные явления в волноводе. Заключение Предложена численная модель акустических процессов на основе динамики взаимодействующих частиц в кубической объемно-центрированной решетке. Мы показали способность модели описывать фокусировку акустической волны линзой, колебания волновода и распространение волны в твердом теле. Модель позволяет учитывать нелинейные эффекты в среде, задавая силовую функцию. На основе ультразвукового волновода показана возможность описания резонансных явлений. Также была сделана оценка оптимального временного шага для численного моделирования.

Скачать электронную версию публикации