Моделирование задачи идентификации положения полости в упругом препятствии по рассеянному звуковому полю Текст научной статьи по специальности «Физика»

Известия Тульского государственного университета Естественные науки. 2011. Вып. 3. С. 74-86 Механика

УДК 539.3

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

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

Аннотация. Предлагается вариационная модель решения обратной задачи рассеяния звука упругим шаром с неконцентрической полостью. При решении определяется угловое положение центра полости в плоскости у = 0. Алгоритм определения строится на основе решения прямой задачи о рассеянии плоской звуковой волны упругим препятствием с полостью, заполненной жидкостью. Показана устойчивость решения задачи при различном числе точек измерения давления и в условиях действия ошибок измерения. Оценена точность нахождения положения полости при различных величинах ошибок.

Ключевые слова: рассеяние звука, гармоническая плоская волна, упругое препятствие, геометрическая обратная задача рассеяния.

Одна из основных целей решения задач о рассеянии звука различными препятствиями — создание теоретической основы для разработка практических технологий ультразвуковой дефектоскопии, гидроакустики, медицинской томографии и др. В этих важнейших областях приложения акустических методов большая часть задач сводится к идентификации параметров рассеивателя (расположения, геометрической конфигурации, свойств материала) по отраженному звуковому полю. Такие задачи относятся к классу обратных задач [1-3]. В теоретическом плане полное решение таких задач представляет существенные трудности. Влияние многих искомых параметров на результат решения задачи о рассеянии звука является нелинейным, что не позволяет использовать аппарат линейных интегральных уравнений.

Для решения ряда практических задач достаточно эффективными оказываются вариационные методы, при использовании которых ищется приближенное решение, доставляющее экстремальное значение

* Работа выполнена при финансовой поддержке РФФИ (проект № 11-01-97509-р-центр-

а).

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

Общий принцип использования такого подхода к решению обратных задач рассеяния звука был декларирован в [4]. В данной работе рассматривается решение одной обратной задачи о рассеянии звука упругим шаром с полостью.

Пусть в однородной идеальной жидкости с плотностью ро и скоростью звука Со распространяется плоская звуковая волна. Предполагается, что эта волна является гармонической с зависимостью от времени вида ехр(—гиЬ), где и — круговая частота.

Волна рассеивается упругим шаром радиуса Е1, внутри которого находится сферическая полость радиуса К2 < К\, заполненная идеальной жидкостью. Обозначим через р, ц и Л плотность и модули упругости Ламе материала упругой среды, а через р1 и С1 плотность и скорость звука идеальной жидкости в полости. В общем случае центры шара О и полости О1 не совпадают, т.е. допускается, что поверхности препятствия и полости не являются концентрическими поверхностями.

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

О и О1 .

Для решения задачи введем декартову систему координат х, у, г так, чтобы падающая волна распространялась вдоль оси г, а центры шара О и полости О1 располагались в плоскости у = 0 (рис. 1).

Символом Оо на рисунке обозначена внешняя акустическая среда, О обозначает область шара, заполненную упругой средой. Через Г обозначена

сферическая поверхность препятствия, а Г1 — поверхность полости П1. Будем считать, что потенциал скоростей в падающей волне Фр имеет вид

Фр = ехр(г(ког — и£)), (1)

где к0 = и/с0 — волновое число в содержащей жидкости; £ — время; (амплитуда потенциала без ограничений общности полагается равной единице).

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

Vр = gradФp. (2)

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

Если считать, что расстояние между центрами полости и тела фиксировано и равно Ко = |ОО11, то искомым параметром, определяющим положени е полости П1 в плоскости у = 0, будет являться угол во между осью г и радиус-вектором точки О1 в системе координат х, у, г (рис. 2).

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

Вначале рассмотрим некоторые основные моменты решения прямой задачи о рассеянии плоской звуковой волны (1) указанным препятствием. Полное решение такой задачи получено в работе [5].

Для аналитического представления решения задачи о рассеянии вводится дополнительная декартова система координат Х1, у1, г1 с началом в точке

О1 так, чтобы ее оси были сонаправлены с осями системы х, у, г. В этом случае радиус-векторы г и Г1 некоторой точки пространства М в системах х, у, г и Х1, у1, г1 соответственно связаны так, что

т = 00і + Ті.

Далее с системами х, у, г и хі, уі, гі связываются сферические системы координат т, 9, ф и ті, 9і, фі. В этих системах поверхности Г и Гі задаются уравнениями т = Яі, Ті = Я і соответственно.

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

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

д,т' 1 ко т

ДФ5 + кОТ = 0, (3)

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

V = grad(Фp + Ф8). (4)

Упругие колебания внутри препятствия (в области П) можно описывать потенциалами Ф, р продольных и поперечных волн соответственно, которые должны удовлетворять волновым уравнениям:

ДФ + к2Ф = 0, (5)

ДФ + Х2Ф = 0, (6)

! Р Р

где к = ------, х = и \------волновые числа продольных и поперечных

у Л + 2ц у ц

волн в упругом материале препятствия.

При этом смещение частиц в упругом материале П определяется

выражением _ _

и = gradФ + го1Ф. (7)

Возбужденные внутри полости П1 звуковые колебания будем описывать потенциалом скоростей частиц жидкости Фс, который должен удовлетворять уравнению Гельмгольца

ДФс + к2Фс = 0, (8)

1 и

где к1 =------волновое число для жидкости в полости.

с1

При этом скорость частиц акустической среды в полости будет

определяться выражением __________

V1 = gradФc. (9)

Решение уравнений (3), (5), (6), (8) проводится с учетом граничных условий.

На поверхности Г должны выполняться условия равенства нормальных скоростей и напряжений в акустической среде По и упругом материале П:

—шиг |г=я1 = Vr, Стт |Г=Я1 = —р, атв |Г=Я1 = 0, аг^\г=Я1 = 0, (10)

где величины слева записаны для П, а справа — для По; иг, Уг — радиальные компоненты векторов (7), (4) в системе г, в, ф; р = шр0(Фр + Ф*) — акустическое давление в жидкости По; игг, аг$, аг^ — компоненты тензора напряжений в упругом материале П в системе координат г, в, ф. Аналогичные условия на поверхности Г1 имеют вид:

—ши(1) = ^(1), а"1" = — р1, &Г1 =0, а(1) =0,

' Г1 = Я2 " " Г1=Я2 " Г1=Я2 Г1=Я2

(11)

где величины слева записаны для П, а справа — для П1; верхний индекс (1) показывает, что компоненты векторов и тензоров записываются в системе координат г1, в1, ф1; р1 = шр0Фс — акустическое давление в жидкости П1.

Кроме того, для потенциала скоростей в рассеянной волне Ф* должны выполняться условия излучения на бесконечности [6].

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

ГО П

Ф* = ]Т ]Т АптНп(ког)Рт(еавв)егт^, (12)

п=—го т=—п

где Нт(х) — сферическая функция Ханкеля первого рода порядка п; Р^(х) — присоединенная функция полинома Лежандра степени п порядка т; Апт — неизвестные коэффициенты, подлежащие определению из граничных условий (10), (11).

В аналогичной форме можно получить и решение уравнения (8) в сферической системе координат г1, в1, ф1

ОО

Фс = ]Т ]Т БптЗп(к1п)Рт(сазв1)егт^. (13)

= т=

Здесь в качестве решения уравнения Бесселя выбрана функция ]п(х), поскольку потенциал (13) должен быть ограниченной функцией при х = 0. Решение уравнения (5) ищется в виде

Ф = Ф1 + Ф2, (14)

где Фі — решение уравнения (5) в сферической системе координат т, 9, ф, а Ф2 — решение того же уравнения во второй сферической системе координат гі, 9і, фі.

В Фі в качестве функции Бесселя будем использовать функцию іп(х) — сферическую функцию Бесселя первого рода порядка п , которая является регулярной при х = 0. Это необходимо потому, что в том случае, когда точка О не попадает внутрь области Пі, то ,]п(х) надо будет вычислять для х = 0. Таким образом, аналогично (13) Фі представим в виде

оо п

Фі = Е Е СпшІп(кт)РПт(соз9)егт^. (15)

п=-о т=—п

Второе слагаемое в представлении (14) строится так, чтобы оно соответствовало бегущей волне, отходящей от границы Гі. Таким образом, потенциал Ф2 ищется в виде

о п

Ф2 = Е Е °птЬп(кті)Рт(еОв9і)егт^. (16)

п=-о т=—п

Аналогично (14) решение уравнения (6) будем искать в виде

Ф = Фі + Ф2, (17)

где Фі — решение уравнения (6) в сферической системе координат т, 9, ф, а Ф2 — решение того же уравнения во второй сферической системе координат Ті, 9і, фі. _

Согласно [8] каждую из векторных функций Ф- (і = 1,2) представим через 2 скалярные Ф-д, Ф-2 в виде

Ф- = тФ^ё. + го1(тФ-2ёг), (18)

где ёг — вектор физического базиса координаты т (при і = 1) или

координаты ті (при і = 2).

При этом каждая из функций Ф-, (і = 1,2; д = 1,2) удовлетворяет скалярному уравнению вида (6).

Решая скалярные уравнения вида (6) в координатах т, 9, ф, представим Фіі и Фі2 в виде аналогичном (15):

о

Фі, = Е Е Ептіп(хт)Рт(соз9)егт*. (19)

п=-о т=—п

Решая уравнения вида (6) в координатах т, 9, ф, представим Ф2і и Ф22 в виде аналогичном (16):

о

ф2, = Е Е Р,пшК(хті)Рт(со89і)е^1. (20)

п=-о т=—п

Коэффициенты Апт1 Впт1 Спт1 Опт1 ЕдПт1 РдПт в разложенияХ (12), (13),

(15), (16), (19), (20) подлежат определению из граничных условий (10), (11) с учетом (1), (4), (7), закона Гука и выражений компонентов тензора деформаций через компоненты вектора смещения и в сферической системе координат [9].

При рассмотрении граничных условий (10) будем использовать разложение плоской падающей волны (1) по сферическим гармоникам [10]. Кроме того, при удовлетворении этих граничных условий разложим представления (16), (20) по базисным решениям уравнений Гельмгольца в смещенной относительно Т\, 9\, ф1 сферической системе координат г, 9, ф. Для этого воспользуемся соответствующими формами теорем сложения для базисных решений в сферической системе координат [11].

Тогда при подстановке соответствующих представлений в функции, участвующие в граничных условиях (10), с учетом того, что на Г г = = Я\ граничные условия сведутся к равенствам рядов по сферическим гармоникам. Ввиду ортогональности последних из равенств рядов следует необходимость совпадения коэффициентов в рядах. Условия совпадения коэффициентов дадут четыре ряда линейных уравнений, связывающих коэффициенты Апт, Спт, Опт, Едпт, Ечпт для различных пар индексов п и т.

Аналогично при удовлетворении граничных условий (11) разложим представления (15), (19) по базисным решениям уравнений Гельмгольца в смещенной относительно г, 9, ф сферической системе координат Г1, 91, ф1. Соответственно из ортогональности сферическим гармоник уже на Г1 (г! = Я2) получим еще четыре ряда линейных уравнений, связывающих коэффициенты Bnm, Cnm, Dnm, Eqnm, Рдпт,.

Задавшись предельным значением N индекса п в суммах (12), (13), (15),

(16), (19), (20) и объединив уравнения, полученные из граничных условий (10) и (11) получим конечную систему линейных алгебраических уравнений для определения коэффициентов, Апт, Спт, Опт, Eqnm, ЕЧпт (п = 1, 2,...

..., N; т = -п, —п + 1,..., п; д = 1, 2). Тем самым получаем приближенное решение задачи о рассеянии звука шаром с полостью.

Реализовав представленный алгоритм решения прямой задачи, можно численно рассчитать коэффициенты Апт, а, значит, и определить теоретические значения давления рассеянного звукового поля в любой точке X(г,9,ф) области По по формуле

N п

Ро = шро^8 = гиро ^ ^2 АптЬп(ког)Рпт(сов9)егт1р. (21)

п=—N т=-п

Предполагается, что X представляет собой радиус-вектор точки пространства с координатами г, 9, ф.

Заметим, как коэффициенты Апт, так, естественно, и само расчетное значение давления Ро(Х) зависят от значения 9о, определяющего положение

полости внутри шара-препятствия. Поэтому будем обозначать их Апт(9о) и Ро(Х; 9о).

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

Пусть в результате наблюдений процесса рассеяния звука в окрестности препятствия измеряется давление в рассеянной волне во множестве М точек наблюдения V = {X1, X2,..., Xм}. Обозначим через Р(X; 91) функцию давления в рассеянной волне в точке X при положении полости, задаваемом углом 9о = 91. Тогда совокупность измеренных значений можно обозначить вектором р = (Р1, Р2,..., Рм), где ри = Р(Xи; 91).

Здесь через 91 обозначено действительное положение полости, при котором произведены измерения. В решаемой же задаче идентификации угол 91 выступает в качестве искомой неизвестной величины.

Очевидно, что если вычисления в (21) и измерения в точках V проведены достаточно точно, то должно наблюдаться совпадение наблюдаемых и расчетных значений давления при 9о = 91, т.е. ри = P0(Xи; 91). Тогда задачу определения угла 91 действительного положения полости по измеренным значениям давления можно представить как задачу нахождения такого угла 9о в расчетных значениях (21), при котором наблюдается минимальное совокупное отклонение расчетных значений от измеренных.

Будем характеризовать совокупное отклонение расчетных значений Pо(Xи; 9о) от измеренных ри по всему множеству точек наблюдения V евклидовой нормой О(9о), определяемой следующим образом

О(9о) = \\р — Ро\\у, (22)

где Ро — вектор расчетных значений такой, что его компоненты Рои = = Ро^и; 9о); значок V показывает, что значение нормы существенно зависит от используемого множества точек наблюдения.

Следовательно задачу определения положения полости внутри упругого шара можно сформулировать как задачу нахождения такого значения 9о, которое доставляет минимальное значение нормы О(9о) при фиксированном векторе измерений р. Будем обозначать соответствующее значение 9*. Следовательно

шш О(9о) = Б(в*). (23)

во

Таким образом, в математической постановке задача идентификации положения полости в упругом шаре может быть сформулирована так: найти значение 9о = 9* такое, что

О(9о) ^ шш, (24)

где 0 ^ 9о ^ п; О(9о) определяется соотношением (22) при заданном

множестве точек наблюдения V; компоненты вектора Ро вычисляются по

(21); коэффициенты Апт = Апт(9о) (21) находятся из решения задачи о

рассеянии плоской звуковой волны (3), (5), (6), (8), решаемой с учетом граничных условий (10), (11).

Данная задача является нелинейной задачей оптимизации с одной переменной 9о. Ее решение позволяет получить значение 9*, которое является оценкой действительного значения угла положения полости 91, для которого получены измеренные значения давления р1, р2, ..., рм.

Точность решения будет зависеть от выбора частоты и падающей волны, от выбора совокупности точек наблюдения V и от точности измерения ри. Очевидно, что точность может быть повышена за счет увеличения числа точек наблюдения V и путем проведения измерений при различных частотах и>1, и>2, ..., ик. Заметим, что это не требует изменения общей схемы решения задачи (24), однако вызывает увеличение размерности векторов р, Ро и, соответственно, объема вычислений.

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

Рассматривалось препятствие, упругий материал которого имеет следующие характеристики: р = 2700 кг/м3, Л = 5.3 • 101С1 Н/м2, л = 2.6 • 101о Н/м2. Внешняя среда и жидкость в полости представляют собой идеальную жидкость с плотностью р1 = ро = 1000 кг/м3 и скоростью звука с1 = со = 1485 м/с. Соотношение радиусов полости и шара Е1 /К2 = 2. Расстояние от центра полости до центра шара определяется соотношением К1/Ко = 0.25. Частота падающей волны полагалась такой, что коК1 = 3.

Поскольку рассматривался вариант введения систем координат, связанных с центрами шара и полости так, как указано в постановке задачи, то центр полости О1 всегда находится в плоскости у = 0, имея значение азимутального угла ф = 0 или ф = п. Для того, чтобы в представлении результатов не оперировать значениями двух координат центра полости 9 и ф, введем условное обозначение угла 9о для точки О1 так, чтобы ее положение в полуплоскости ф = п задавалось значениями 9о > п на основе связи 9о = 2п — 9, где 9 — обычный полярный угол точки О1 в сферической системе координат г, 9, ф.

При численных исследованиях полагалось, что экспериментальное измерение давления производится на поверхности сферы с радиусом К таким, что К/К1 = 10. При этом система точек наблюдения X1, X2,..., Xм размещалась равномерно на окружности радиуса К в плоскости у = 0, начиная с точки со значением 9 = 0.

В качестве измеренных значений давления в рассеянной волне брались расчетные значения давления Ро^и; 91) с учетом возможных ошибок измерения. Так что величина измеренного давления определялась по формуле

ри = pо(X и; 91) + ^ (25)

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

л/2 —

а = а— тах\Ро(Хк; вг)\ (26)

2 к

а — параметр, характеризующий уровень шума измерений относительно максимального расчетного значения давления.

Определение значения 9* при решении задачи (24) выполнялось поиском минимального значения В(9о) в дискретном ряду значений В для 0 ^ 9о < 2^ с шагом Д9 = ^/180. Таким образом, потенциальная точность определения 9* составляла 1°.

На рис. 3 представлен результат определения 9* для случая, когда полость расположена под углом 91 = 40°, точки наблюдения размещены через 30° (т.е. число точек наблюдения М = 12), уровень погрешности измерений задан значением а = 0.01 в (26).

Рис. 3. Зависимости В(90) при 9\ = 40°, М =12, а = 0.01

На рисунке сплошной линией представлен график зависимости В(9о) с учетом псевдослучайной реализации погрешностей измерения е в (25). Для сравнения пунктирной линией представлена зависимость В(9о) без учета погрешностей измерений. Значения В на оси ординат приведены в относительных единицах после деления на ||Ро||у, а угол 9о на оси абсцисс

задан в градусах. Вертикальная штриховая линия на графике показывает значение угла в\ положения полости.

Из рисунка видно, что график В(во) без учета погрешностей измерений имеет строгий глобальный минимум (со значением 0) в точке во = 0\. Функция В(во) с учетом погрешностей имеет также глобальный минимум — в* — в окрестности в1, однако не совпадает с в\. Значение в* = 41°. Ошибка определения положения составляет 1°.

На рис. 4 представлены аналогичные зависимости для случая, когда погрешность измерений задана значением а = 0.1. На графике видно, что такой уровень погрешностей при М = 12 не позволяет хотя бы приближенно определить положение полости. Здесь ошибка определения составляет 99° (в* = 139°)!

Рис. 4. Зависимости В(в0) при в\ = 40°, М =12, а = 0.10

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

Из таблицы видно, что малое число точек наблюдения при погрешности измерений а = 0.1 в ряде случаев не дает возможности точно определить положение полости. При малых погрешностях измерения (а = 0.01) ошибка определения угла в1 не превышает 2°. Количество точек наблюдения М = 36 при рассмотренных уровнях погрешности измерений позволяет определять положение полости с точностью не хуже чем 3°.

Таблица 1

Результаты определения положения полости

0ъ° M а 0*0 ошибка

40 12 0.01 41 1

40 12 0.10 139 99

40 36 0.01 41 1

40 36 0.10 38 2

175 12 0.01 175 0

175 12 0.10 177 2

175 36 0.01 175 0

175 36 0.10 173 2

300 12 0.01 299 1

300 12 0.10 307 7

300 36 0.01 301 1

300 36 0.10 297 3

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

Список литературы

1. Агранович З.С., Марченко В.А. Обратная задача теории рассеяния. Харьков: Изд-во Харьковского ун-та, 1960. 268 с.

2. Ватульян А.О. Обратные задачи в механике деформируемого твердого тела. М.: Физматлит, 2007. 224 с.

3. Bui H.D. Inverse problems in the Mechanic of Materials: An Introduction. CRC Press, Boca Ration, FL, 1994. 318 p.

4. Скобельцын С.А. О вариационной постановке задачи идентификации материала упругого слоя по характеристикам прохождения звуковых волн // Вестник ТулГУ. Сер. Дифференциальные уравнения и прикладные задачи. 2009. Вып.1. С.82-91.

5. Толоконников Л.А., Филатова Ю.М. Дифракция плоской звуковой волны на упругом шаре с произвольно расположенной сферической полостью // Изв. ТулГУ. Естественные науки. 2010. Вып.1. С.114-122.

6. Исакович М.А. Общая акустика. М.: Наука, 1973. 496 с.

7. Кошляков Н.С., Глинер Э.Б., Смирнов М.М. Уравнения в частных производных математической физики. М.: Высшая школа, 1970. 712 с.

8. Морс Ф.М., Фешбах Г. Методы теоретической физики. Т.2. М.: Изд. иностр. лит., 1960. 886 с.

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

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

11. Ерофеенко В.Т. Теоремы сложения: Справочник. М.: Наука и техника, 1989. 255 с.

Иванов Валерий Иванович ([email protected]), д.ф.м.-н., профессор, декан, механико-математический факультет, Тульский государственный университет.

Скобельцын Сергей Алексеевич ([email protected]), к.ф.м.-н., доцент, кафедра прикладной математики и информатики, Тульский государственный университет.

Modelling of a identification problem of a cavity location in an elastic body by an scattered acoustic field

V. I. Ivanov, S. A. Skobeltsyn

Abstract. Variational model is proposed solving the inverse scattering problem of sound elastic sphere with non-concentric cavity. In the solution is determined by the angular position of the center of the cavity in the plane y = 0. The algorithm is based on the determination of the direct problem of scattering of plane sound waves by elastic obstacle with a cavity filled with fluid. The stability of solutions of the problem for different numbers of points of pressure measurement and according to the measurement errors. Evaluated the accuracy of location of the cavity for different values of errors.

Keywords: scattering of sound, harmonic plane wave, elastic obstacle, geometric inverse scattering problem.

Ivanov Valeriy ([email protected]), doctor of physical and mathematical sciences, professor, dean, mathematical and mechanical faculty, Tula State University.

Skobeltsyn Sergey ([email protected]), candidate of physical and mathematical sciences, associate professor, department of applied mathematics and computer science, Tula State University.

Поступила 12.10.2011