Моделирование решений задач акустики с использованием МКЭ Текст научной статьи по специальности «Физика»

Известия Тульского государственного университета

Естественные науки 2008. Выпуск 2. С. 132-145

= ИНФОРМАТИКА =

УДК 539.3

В.И. Иванов, С.А. Скобельцын

Тульский государственный университет

МОДЕЛИРОВАНИЕ РЕШЕНИЙ ЗАДАЧ АКУСТИКИ С ИСПОЛЬЗОВАНИЕМ МКЭ

Аннотация. Предлагается модель использования метода конечных элементов (МКЭ) для решения задачи о рассеянии плоской звуковой волны некруговым неоднородным упругим цилиндром. В области жидкости, прилегающей к цилиндру, выделяется часть с внешней поверхностью, имеющей форму кругового цилиндра. В этой части жидкости, а также внутри упругого цилиндра решение ищется численно с помощью МКЭ. Во внешней части жидкости решение представляется в аналитической форме в виде рядов по функциям Бесселя.

Метод конечных элементов является широко используемым инструментом решения практических задач гидродинамики и теории упругости [1, 2]. При решении задач, связанных с изучением звуковых колебаний в основном рассматриваются случаи ограниченных областей: акустические свойства помещений, собственные частоты сложных ограниченных объемов. Однако мало исследований такого рода посвящено изучению рассеяния звуковых волн сложными объектами в неограниченной среде. Применение МКЭ для решения таких задач позволит оценить влияние, как сложной формы рассеивателя, так и особенностей реальных сред: неоднородности, анизотропии, вязкости. Такие исследования могут служить средством анализа аналитических решений и инструментом разработки практических технологий в ультразвуковой дефектоскопии, гидролокации. В данной работе рассматривается решение с использованием МКЭ задачи о рассеянии плоской звуковой волны некруговым неоднородным анизотропным упругим цилиндром.

Пусть В однородной идеальной ЖИДКОСТИ С ПЛОТНОСТЬЮ Ро и скоростью звука со распространяется плоская звуковая волна. Предполагается, что падающая волна является гармонической с зависимостью от времени вида exp (—icot), где ш - круговая частота.

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

С учетом характера падения волны и типа неоднородности цилиндра движение частиц жидкости и препятствия будет являться плоским. Поэтому геометрия задачи может быть представлена рис. 1. На рисунке Фр - потенциал смещений в падающей волне; По - бесконечная область, занимаемая содержащей жидкостью; р, \ijki — плотность и модули упругости материала цилиндра.

Q0

Рис. 1. Геометрия задачи

Требуется найти полное поле установившихся упругих колебаний внутри цилиндра, звуковое поле рассеянных волн в ближней и дальней зонах. Решение будем искать в рамках линейных теорий гидродинамики и термоупругости.

Решение задачи будем проводить в декартовой системе координат, которую выберем так, чтобы ось Oz была направлена вдоль образующей цилиндра, а начало координат располагалось внутри области Г1. Таким образом, оси Ох и Оу будут находиться в плоскости нормального сечения цилиндра, а модули упругости и плотность материала неоднородного цилиндра будут являться функциями координат х и у

^ijkl = Р = р{х,у).

Направление оси Ох выберем совпадающим с направлением распространения падающей волны. Тогда потенциал смещений в падающей волне можно представить в виде

Фр = exp (i(kox — uj t)), (1)

где амплитуда потенциала без ограничений общности полагается равной единице, а ко = со/со - волновое число в содержащей жидкости.

Звуковые волны во внешней среде описываются полной системой линеаризованных дифференциальных уравнений для малых возмущений идеальной жидкости [3], из которых в случае гармонических колебаний для потенциала смещений получают уравнение Гельмгольца

ДФо + ^оФо = 0, (2)

где Фо = Фр + Фв - потенциал смещений в суммарном акустическом поле, Ф.ч потенциал смещений в рассеянном поле.

Потенциал Ф8 также удовлетворяет уравнению (2). При этом смещение частиц жидкости щ и давление ро в ней определяется соотношениями

no = grad Ф0 = grad (Фр + Фв), р0 = /?о^2Фо = /?о^2(ФР + Ф8). (3)

Под действием звуковой волны в упругом цилиндре возникают гармонические колебания. Колебания частиц материала цилиндра подчиняется общим линейным уравнениям движения сплошной среды, которые в случае плоских деформаций в декартовой системе координат имеют вид [4]

до'хх д<7Ху 2 д<7Ху д(7уу 2 / л\

—— + —^ = -рш*их , = -рш*иу , 4)

Ох оу Ох оу

где иу - компоненты вектора смещений и в упругой среде цилиндра;

<7. (Тху, Оуу - КОМПОНеНТЫ 'ГеПЗОра ПИПрЯЖеПИЙ (Т.

В соответствии с обобщенным законом Гука [4] компоненты тензора напряжений выражаются через компоненты вектора смещений связью

= '^iзpq^pq1 (5)

где индексы /.р. ([ принимают значение 1, которое соответствует координате х, или 2, соответствующее координате у; еря - компоненты тензора малых деформаций:

_ ди* _9иу

“ “ дх ’ ху - 2 [дх + ду ) ’ уу ~ ду ■ ( 1

На границе сопряжения сред Г должны выполняться следующие граничные условия:

- должны совпадать нормальные компоненты смещения частиц взаимодействующих сред

ип |р = и0п 5 (7)

где ип = и ■ п; иоп = Щ ■ Щ п - единичный вектор нормали к границе Г, направленный в

нормальное напряжение в упругом цилиндре определяется давлением в окружающей среде

где апп = (п • сг) • п;

касательные напряжения в граничном слое цилиндра отсутствуют

где апт = (п ■ (т) • т. т единичный касательный вектор к границе Г (г • п = 0).

Кроме того, при г —оо (г = л/х2 + у2) для потенциала смещений в отраженной акустической волне Ф.ч должны выполняться условия излучения на бесконечности [5]

Таким образом, в математической постановке задача состоит в нахождении потенциала Ф.ч(.7\ у) в области Q{) и компонент вектора смещений их(х,у), ity(x. у) в области О, удовлетворяющих соответственно уравнениям (2) и (4) с учетом (5), (6), а также граничным условиям (7), (8), (9) на Г и условиям (10). Здесь и далее для искомых и заданных функций, описывающих движение, временной множитель exp (—iut) (как в падающей волне) для простоты не указываются.

Заметим, что, если уравнение (2) и может быть решено аналитически, то для произвольного типа анизотопии и неоднородности материала упругого цилиндра уравнения (4) не могут быть решены аналитически. Эффективным численным методом решения уравнений теории упругости для ограниченных областей (таких, как О в рассматриваемом случае) является метод конечных элементов. Однако традиционные технологии его использования предполагают классическое представление граничных условий - задание напряжений, смещений или соотношений между ними. В данной задаче в граничных условиях для уравнений (4) присутствует неизвестная функция Ф8(ж, у), которую нельзя искать традиционным МКЭ, т.к. область Q{). в которой нужно определить Ф.ч является неограниченной.

Далее предлагается подход, который позволяет на первом этапе решения исключить неизвестную функцию Ф.ч из граничных условий (7), (8), (9), затем решить задачу нахождения поля деформаций (и напряжений) в упругом цилиндре с использованием МКЭ. На последнем этапе определяется потенциал рассеянной волны Ф.ч. Общая идея такого подхода изложена в сообщениях [6, 7] и работе [8].

Используя известные решения задач о рассеянии звука на круговых жестких и упругих цилиндрах [5, 9], получим аналитическое решение уравнения

®пп 1г — POi

(8)

(9)

вида (2) для потенциала Ф.ч в цилиндрической системе координат, связанной с декартовой системой координат, введенной выше. Общее решение уравнения (2) с учетом условия (10) имеет вид

ОО

Ф»= $3 АтНт(к0гут'<’, (И)

т= — оо

где Нт(х) - функция Ханкеля первого рода порядка га; г, <р - цилиндрические координаты; Ат - неизвестные коэффициенты, подлежащие определению из граничных условий (7), (8).

В аналогичной форме можно представить потенциал (1) смещений в падающей волне [5]

ОО

фр= Е 1тЫк0г)е^, (12)

т=— оо

где 7т = /т; '1,п (.7-) - функция Бесселя порядка га.

Согласно МКЭ представим приближенно смещение в упругой среде (в области О) в виде

N

и = ^икфк(х,у), (13)

к=1

где (/д. - неизвестные узловые значения вектора и. фк(х,у) - известные координатные функции, образующих близкую к полной систему «почти» ортогональных функций; N - число узлов в области О после разбиения на

конечные элементы.

Далее, ограничившись в представлениях (11), (12) потенциалов Ф.ч. Фр только конечным числом слагаемых, полагая, например, вместо (11)

м

*«= £ АтНт{к0г)е™*, (14)

т= — М

и воспользовавшись одним из граничных условий (7)- (9), можно получить выражение коэффициентов Ат (га = — М, М) через узловые значения (/д..

Однако, при произвольной форме границы Г придется раскладывать в ряд Фурье по координате <р не только координатные функции о^(х. у), учитывая, что на границе Г х = х(<р), у = у(<р), но и функции из граничных

условий, зависящие от функций Бесселя из представлений (14) и аналогичного для Фр, поскольку в этом случае на границе и г = г(<р).

Для того, чтобы избежать этой процедуры «расширим» препятствие О до виртуального кругового препятствия О и О1, включив в его область и слой жидкости Ох, прилегающий к упругому цилиндру. Геометрическое изменив постановки задачи иллюстрируется на рис. 2

Рис. 2. К изменению постановки задачи

На рисунке область Пх представляет слой окружающей жидкости, прилегающий к упругому цилиндру. Его толщина имеет величину, сравнимую с размерами упругого цилиндра. Внутри границей слоя выступает поверхность упругого препятствия Г, а снаружи - поверхность кругового цилиндра Го с осью, совпадающей с осью О г выбранной ранее системы координат. Пусть внешний радиус слоя Пх равен г о.

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

Обозначим потенциал смещений частиц жидкости в области Пх через Ф. Тогда он, как и Ф5, должен удовлетворять уравнению Гельмгольца

ДФ + ^Ф = 0. (15)

Теперь на границе Г должны выполняться граничные условия, аналогичные условиям (7)-(9), с учетом замены щ и ро на й\ и р\ такие, что

{¡х = gradФ, рх = ро^2^. (16)

На внешней границе области Пх должны выполняться следующие условия:

- должны совпадать нормальные компоненты смещения частиц взаимодействующих сред

^1г|г0 = и0г] (17)

- должны совпадать давления во внешней и внутренней жидкостях

Р11г0 = Р®’ (18)

Таким образом, измененная математическая постановка сводится к необходимости решения уравнений (2), (4) и (15) при выполнении модифицированных граничных условий (7)-(9) на Г и условий (17), (18) на Го, а также условия излучения (10).

Для решения уравнения (15) в области жидкости Пх также, как и уравнений движения в П, будем использовать метод конечных элементов. Будем искать Ф в виде, аналогичном (13)

где N1 - число узлов в области Пх после разбиения на конечные элементы; Ф/с - узловые значения потенциала Ф; фь(х,у) - координатные функции (они, конечно, не совпадают с фь(х,у) из (13), но поскольку играют такую же роль, используем для них то же обозначение).

Пример разбиения области Пх и П на треугольные конечные элементы представлен на рис. 3.

Чтобы непосредственно приступить к поиску узловых значений [/^иФ^ в представлениях (13), (19) в соответствии с одной из обычных технологий МКЭ [1] необходимо исключить из граничных условий неизвестную функцию Ф5, точнее, коэффициенты Ат в ее разложении (11).

Для этого воспользуемся условием на давления на границе Го- Подставляя в (18) выражения (19), (11) и (12), получим

(19)

* У

->

X

Рис. 3. Разбиение на конечные элементы

где X = Го сов (р; у = Го ЭШ (р; го = кдГд.

С учетом выражений для ж и у на Го обозначим фк(х,у) в (20) так Фк(х,у) = фк(1р)- Используя очевидную периодичность функций фк{(р) с периодом 2тт. разложим координатные функции на Го в ряды Фурье оо 1 р2тт

Фк(ф) = ^2 4тегт^, Ще 4т = — фк^)е~гШ(И. (21)

2тг JQ

т= —оо и

Подставляя (21) в граничное условие (20) и используя ортогональность функций о1,п'р на [0. 2тг]. получим

V тп : Е ^кт^к — Тпг'Лп (-^о) АтНт (.20 ) • (^2)

к=1

Из последнего соотношения легко получить выражение для Ат

ЛГ1

-А-т — ®гв + ^ ^ Ьтк^к1 (23)

к=1

Где От — /Нт (-2о) 1 Ь-тк — -¡-кт / Нт (-2о) •

Рассмотрим далее правую часть граничного условия (17)

, з(фр + ф8)

Щг -

I Г —Г о

сю

= Е Ы7пЛ(2о) + АтН’гп(г0))ёт^: (24)

г=го т= —оо

ще штрих обозначает производную функции по аргументу.

Подставляя выражение (23) в (24), получим

оо N1 / оо \

“ог|г=г„= Е атИе^ + Ё ф» Е , (25)

гв=— оо &=1 \ гв=— оо /

где 5т = -2г'ут/(тгг0Нт(го)); Укт = коН'т(го)Ътк.

Введем обозначения

СЮ СЮ

яМ= Е Зт(¥')е™*’, = - (26)

т= — ос

Тогда (25) примет вид

гв= — оо гв= — оо

ЛГІ

11ог

</М-Еуг*мф*- (27)

г=г0

А-=1

Подставляя (27) в граничное условие (17) на Го, приведем его к виду

¿>ф Мі

* *=!

= гЫ- (28)

г=г0

Таким образом, для решения уравнений (4), (15) с использованием представлений неизвестных функций в виде (13), (19) требуется выполнение модифицированных граничных условий (7)-(9) и (28), которые в качестве неизвестных содержат только узловые значения IIк и Далее можно исполь-

зовать стандартные приемы МКЭ для решения краевой задачи (4), (15), (7)-(9), (28) в ограниченной области Пх и П. После нахождения узловых значений Vк и по формулам (23) определяются коэффициенты в реше-

нии (14) для рассеянной волны в неограниченной внешней области По-Предложенная технология была применена для решения задачи о рассеянии звука упругим цилиндром с эллиптическим сечением (рис. 4).

Рис. 4. Конечно-элементное разбиение Пх и П для случая эллиптического цилиндра

В качестве инструментального средства реализации МКЭ использовался математический пакет ЕЕМЬАВ 2.0 на базе Ма^АВ [10].

Рассматривалось рассеяние звука в воде (ро = 1000 кг/м3, со = 1485 м/с) упругим цилиндром, материал которого имеет плотность р = 2700 кг/м3, модуль упругости Е = 6.944 • Ю10 Н/м2, коэффициент Пуассона г/ = 0.335.

Анализ рассеянного звукового поля в дальней зоне проводится по нормированному значению потенциала смещения отраженных волн в соответствии с выражением

Пч>)

у/ко

оо

]Г (-г)тАте

гтср

т-

-оо

(29)

Зависимость Р{ср) при изменении ср от 0 до 2тт показывает распределение амплитуды отраженной волны в зависимости от угла наблюдения ср.

На рис. 5-10 представлены результаты расчета Р{ср) для различных соотношений между полуосями сечения цилиндра - а, направленной по оси х, и Ь, направленной по оси у. Расчеты проводились для случаев, когда а + Ь = 2, а частота и падающей волны такова, что ко (а + Ь)/2 = 2 и

волна распространяется в направлении оси Ох. На всех рисунках тонкой линией для сравнения представлена диаграмма направленности рассеянного поля -^о((/?) Для кругового упругого цилиндра, радиус которого ао равен (о + 6)/2 = 1. Кроме того в начале координат осей диаграммы также тонкой линией схематически изображены сечения эллиптического цилиндра О и кругового цилиндра, используемого для сравнения.

Ь/а = 3

Рис. 5. Диаграмма направленности Р(ср) при а = 0.5, Ь = 1.5

Ь/а = 1.857

Рис. 6. Диаграмма направленности Р(ср) при а = 0.7, Ъ = 1.3

(р=п/2

Ь/а = 1.222

р X \х

—> / : \

! | ^=0

—> 1 ■ 1 :

—» \

Рис. 7. Диаграмма направленности Р(ср) при а = 0.9, Ь = 1.1

^ 71/2 Ь/а = 0.818

У \

Р

-> / X

-» 1 V. \ ■ ^=0

-> 1 1 /

-» \

Рис. 8. Диаграмма направленности -Р(^) при а = 1.1, Ь = 0.9

^_7Г'2 Ь/а = 0.538

V / / '

Р —-

' / \| V. \

-» : \ \ (р=0

—» > I 1 / ;

-» ч \ 0

Рис. 9. Диаграмма направленности -Р(^) при а = 1.3, Ь = 0.7

Рис. 10. Диаграмма направленности -Р(у?) при а = 1.5, Ь = 0.5

Диаграммы показывают, что падение волны в направлении меньшей оси цилиндра с большим эксцентриситетом сечения приводит к существенному увеличению коэффициента обратного отражения, и уменьшению выраженности боковых лепестков. Это иллюстрируют диаграммы на рис. 5, 6. В случае распространения волны вдоль большей оси цилиндра с большим эксцентриситетом сечения амплитуда отраженной волны заметно уменьшается (рис. 9, 10), однако лепестки диаграммы становятся более выраженными.

На рис. 9, 11-13 представлены результаты расчета Р(<р) для различных направлений падения волны в случае, когда сечение эллиптического цилиндра вытянуто вдоль оси Ох и а = 1.3, 6 = 0.7. Направление распространения волны задается углом </?о между волновым вектором падающей волны и осью Ох. Подчеркнем, что на диаграммах рис. 5-10 сро = 0.

Ь/а = 0.333

ср=тг/2

Ь/а = 0.538

Рис. 11. Диаграмма направленности Р{}р)

при а = 1.3, Ь = 0.7, (ро = 30°

Рис. 12. Диаграмма направленности Р(Ф) при а = 1.3, Ь = 0.7, (ро = 60°

Рис. 13. Диаграмма направленности .Р(<р)

при а = 1.3, Ь = 0.7, (ро = 90°

Отклонение направления падения волны от оси симметрии цилиндра приводит к существенному нарушению симметрии в диаграмме Р{ср). Например, при (ро = 30° (рис. 11) в диаграмме появляется явно выраженный лепесток в направлении <р = —50° и существенный провал в противоположном направлении. При увеличении угла падения (рис. 12) этот лепесток расширяется и смещается в направление ближе к направлению противоположному волновому вектору падающей волны.

В предельном случае (ро = 90° (рис. 13), естественно, диаграмма f’(v) принимает вид диаграммы для случая а = 0.7, Ъ = 1.3 при vo = 0 (см. рис. 6) повернутой на 90°.

Таким образом, анализ полученных диаграмм позволяет оценить качественное и количественное влияние формы сечения цилиндра на перераспределение амплитуды отраженных волн по углу наблюдения.

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

Библиографический список

1. Галлагер Р. Метод конечных элементов. Основы / Р. Галлагер. -М.: Мир, 1984.

2. Полежаев В.И. Метод конечных элементов в механике вязкой жидкости / В.И. Полежаев, А.И. Простомолотов, А.И. Федосеев // Итоги науки и техн. ВИНИТИ. Сер. Механика жидкости и газа. -1987. -Т. 21. -С. 3-92.

3. Ландау Л.Д. Теоретическая физика, Т.б. Гидродинамика / Л.Д. Ландау, Е.М. Лиф-шиц -М.: Наука, 1988.

4. Новацкий В. Теория упругости / В. Новацкий -М.: Мир, 1975.

5. Скучик Е. Основы акустики, Т.2 / Е. Скучик -М.: Мир, 1976.

6. Скобельцын С.А. Подход к решению задач о рассеянии упругих волн с использованием МКЭ / С.А. Скобельцын // Тез. докл. Международной научной конф. «Современные проблемы математики механики, информатики». -Тула: ТулГУ, 2004. -С. 135-136.

7. Скобельцын С.А. Особенности конечно-элементной формулировки задач о рассеянии звука / С.А. Скобельцын, А.Н. Королев // Материалы Международной научной конференции «Современные проблемы математики, механики, информатики» -Тула: ТулГУ, 2007. С. 205-207.

8. Скобельцын С.А. Использование МКЭ для решения задачи о рассеянии звука ограниченной неоднородной анизотропной термоупругой пластиной / С.А. Скобельцын, А.Н. Королев // Вестник ТулГУ. Сер. Математика. Механика. Информатика. -2007. -Т.13. -Вып. 2. -С. 172-182.

9. Шендеров Е.Л. Волновые задачи гидроакустики / Е.Л. Шендеров. -Л.: Судостроение, 1972.

10. FEMLAB Reference Manual. - Stockholm: COMSOL AB, 2002.

Поступило 20.06.2008