УДК 593.3
ОПРЕДЕЛЕНИЕ ГЕОМЕТРИЧЕСКИХ ПАРАМЕТРОВ ПОЛОСТИ УПРУГОГО ЦИЛИНДРА ПО РАССЕЯННОМУ АКУСТИЧЕСКОМУ ПОЛЮ
С.А. Скобельцын, Н.Ю. Пешков
Разработан метод определения геометрических параметров (положение центра, длины стороны, угла ориентации) полости квадратного сечения в упругом эллиптическом цилиндре по известному рассеянному полю плоской гармонической звуковой волны. Проведена проверка схемы идентификации параметров полости в ходе численного эксперимента. Оценено влияние погрешности измерительных приборов на точность определения геометрических характеристик. Изложенный алгоритм отличается универсальностью и может быть применен для идентификации произвольного количества параметров упругого препятствия.
Ключевые слова: полость с квадратным сечением, рассеянное поле, потенциал смещений, датчики, алгоритм Хука-Дживса, погрешность измерения.
Определение параметров рассеянного звукового поля на основе известных параметров падающей волны, геометрических параметров упругого препятствия и свойств его материала, а также свойств материала содержащей среды не представляет практического интереса.
Однако, существенный интерес представляют задачи определения параметров упругого препятствия по полностью или частично известному рассеянному полю. Такие задачи относятся к классу обратных задач [1, 2, 3]. В теоретическом плане полное решение таких задач представляет существенные трудности. Влияние многих искомых параметров на результат решения задачи о рассеянии звука является нелинейным, что не позволяет использовать аппарат линейных интегральных уравнений. Для решения ряда практических задач достаточно эффективными оказываются вариационные методы, при использовании которых ищется приближенное решение, доставляющее экстремальное значение функционала, величина которого характеризует степень соответствия приближенного решения точному.
Целью представленного ниже исследования является определение геометрических параметров упругой полости с квадратным сечением (положения центра (х0 Уо), длины стороны d и угла ориентации у) по рассеянному полю плоской гармонической звуковой волны.
Предполагается, что бесконечный изотропный однородный упругий цилиндр с эллиптическим сечением I с полуосями а и Ь (а - большая, Ь -малая), материал которого имеет плотность р, модуль Юнга Е и коэффициент Пуассона V, содержащий полость с квадратным сечением II, запол-
148
ненную жидкой (газообразной) средой с плотностью pi и скоростью звука ci, находится в жидкой среде III с плотностью р 0 и скоростью звука c0 (рис. 1).
Рис. 1. Геометрия задачи
Из среды III на эллиптический цилиндр по направлению нормали к малой полуоси его сечения падает монохроматическая гармоническая плоская звуковая волна единичной амплитуды, в которой потенциал смещений частиц жидкой (газообразной) среды имеет вид
Yp = ei(ko r-wt), (1)
где r - радиус-вектор, w = 2p/0 - круговая частота, t - время, а k0 - волновой вектор, определяющий направление распространения падающей волны такой, что |k0| = &0 =w/ c0 - волновое число в содержащей среде
III. Предполагается, что направление оси Ox совпадает с направлением распространения падающей волны. Тогда (1) примет вид:
Y = ei (k о x-wt )
В результате взаимодействия с препятствием волна искажается и образуется рассеянная волна, в которой потенциал смещений частиц жидкой (газообразной) среды обозначим через Ys.
Подразумевается, что в рассеянном поле известно полное акустиче-2
ское давление p = pw Yq, где Yq = Yp + Y5. Величина p фиксируется
конечным числом датчиков. Датчики расположены в узловых точках кругового контура Г радиуса r. Крестовыми маркерами и подписями Vk
(k = 1, K) изображены позиции датчиков. Предполагается, что сами датчики не вносят искажений в процесс рассеяния.
Для определения геометрических параметров полости используется подход, описанный в работах [4, 5], где предлагается вариационная постановка обратной задачи рассеяния звука при определении параметров упругого препятствия.
Подразумевается, что выполняется измерение полного акустического давления (значения потенциала Yq) в процессе рассеяния плоской звуковой волны упругим эллипсом в соответствии со схемой, приведенной на рис. 1, для фиксированных значений параметров полости с квадратным сечением z, = [xQi,yQj,dj,у, J из интервалов [xbХ2], [ybУ2], \dl,d2] и [gi,У2]
. Требуется найти действительные значения параметров полости на основании измеренных данных.
Вектором То (zj) обозначена совокупность измеренных значений
потенциала Yq для параметров zj в точках Vk, т.е.
Т 0 (Zj )=[Yq ji Yq j 2 ^ Yq jK ]. (2)
Через Уообозначается фактическое значение потенциала Yq, в точке Vk, рассчитанное для параметров zj . Очевидно, что в случае отсутствия ошибок измерения
Yq jk = y q jk.
В общем случае измеренные значения Yq jk включают погрешность, которая может быть связана, например, с погрешностями используемых измерительных приборов. Кроме того, сами теоретические (расчетные) значения yqjk могут содержать ошибки из-за неточности задания
параметров модели. Предполагается, что все погрешности в расчетах отнесены к измеренным значениям Yq jk, представленные в виде суммы
y0 Jk =y 0 Jk ^
где e - некоторая случайная величина с нулевым математическим ожиданием и средним квадратическим отклонением e q (полагается, что e имеет нормальное распределение N (0; e q )).
Совокупность рассчитанных значений для заданных параметров 2/ представлена вектором у о (2/ )—У 0/1 У 0 /2 У 0/к ], построенным аналогично (2).
Близость выбранных для вычислений значений параметров 2j к действительным можно охарактеризовать нормой
где у 012
(2 * ) =
)=|* 0 (2, )- У о (2 * К I К к-У 0 к 1
\ к—1
У01 У02
* тс *
У0к \ . Очевидно, что при 2/ — г в случае
— 0.
отсутствия ошибок (е 0 — 0) получаем *)— Т о (2 *)- у о (2 *,
Таким образом, задачу определения параметров квадратной полости можно сформулировать следующим образом: найти такие значения *0,
>0, Л и 7 из интервалов [л^, Х2 ], [>1, >2 ], Ль Л2 ] и [71,7 2 ], при которых совокупное отклонение расчетных значений потенциала ^0 - у о от измеренных значений Т о будет минимальным, или формально
е(--)
® Ш1П.
х0
е[*1,*2] >0е[>1,>2] ] 7е[71,72]
(3)
Задача (3) является нелинейной задачей оптимизации. Ее решение
' *
2 — 2 является оценкой действительных параметров 2 .
Для оценки влияния параметров схемы измерений на точность идентификации характеристик квадратной полости проведены численные эксперименты для значений отклонения е0 — 0.005, 0.01 при действительных значениях параметров 6,2 — [0 м ,0 м ,0.2 м ,2я/ 9].
Предполагалось, что наблюдается рассеяние плоской звуковой волны для значений параметров 2 из интервалов [-0.2,0.2], [-0.2,0.2],
10-3,0.25] и [0, я/ 2].
Количество датчиков К предполагалось равным 11 и 22. Полярные углы датчиков ф^Г определялись формулой
ф\ —(к - 1)бф, где постоянный шаг Дф — 2я/К .
Поиск локального минимума функции (3) от параметров 2 осуществлялся с помощью метода Хука-Дживса [6]. Значения полного акустического давления ^0 на границе Г были получены численно в математическом пакете СОМБОЬ МиШрИуБ^Б [7].
151
При выполнении приближенного решения, выбирались следующие механические и геометрические характеристики:
- размеры полуосей эллипса: а = 0.8 м и Ь = 0.4;
- механические параметры материала эллипса - меди: плотность
р = 8920 кг/м3 ; модуль Юнга Е = 11 • 1010 Па ; коэффициент Пуассона у = 0.3.
- механические параметры жидкости среды II - метанола: плотность р1 = 792 кг/м3 ; скорость звука С1 = 1143 м/с.
Предполагалось, что в акустической среде, заполненной водой с плотностью р 0 = 1000 кг/м3, скоростью звука с0 = 1403 м/с, распространяется плоская гармоническая звуковая волна с частотой /0= 1350 Гц. Для решения методом конечных элементов в ней выделяется слой III, окружающий эллиптический цилиндр. Внешняя поверхность слоя представляет собой круговой цилиндр с радиусом = 8 м и осью, совпадающей с осью упругого эллиптического цилиндра.
Множество точек измерения V\ размещается в слое III на окружности с радиусом / = 1 м.
Результаты расчетов и их визуализация представлены на рис. 2 - 5 и в табл. 1 - 4. Рассмотрены четыре случая действительного расположения полости со стороной й = 0.21 м в цилиндре. Предполагалось, что центр полости расположен в начале координат, а угол у в первом случае равен 41°, во втором - 45°, в третьем - 19°, в четвертом - 13°. Результаты определения положения полости для первого случая представлены на рис. 2 и в табл. 1, для второго - на рис. 3 и в табл. 2, для третьего - на рис. 4 и в табл. 3, для четвертого - на рис. 5 и в табл. 4.
а
Рис. 2. Визуализация итераций решения и зависимость ) для К = 11 и е0 = 0.01: а - визуализация итераций алгоритма решения
(начало)
152
б
Рис. 2. Визуализация итераций решения и зависимость е(г) для К = 11 и го = 0.01: б - зависимость ) от номера итерации
(окончание)
Рис. 2, а - 5, а иллюстрируют решения, полученные на промежуточных итерациях алгоритма Хука-Дживса. Эллипс представляет собой внешний контур сечения цилиндра I. Квадратные контуры показывают сечение полости с параметрами, полученными на очередном шаге. Число на контуре сечения полости на рис. 2, а - 5, а обозначает номер шага у алгоритма поиска минимума (3).
На рис. 2, б - 5, б показаны зависимости нормы г( 2) шага алгоритма от порядкового номера итерации решения. Круговыми маркерами отмечены расчетные значения е^у). Полученные расчетные точки для наглядности соединены ломаной линией. Наличие локальных максимумов на графиках е(2у) показывает, что используемый алгоритм не обладает монотонностью. В частности, это может объясняться сложной зависимостью измеряемых величин акустического давления и совокупностью параметров полости.
а
Рис. 3. Визуализация итераций решения и зависимость е(2) для К = 22 и г о = 0.01: а - визуализация итераций алгоритма решения
(начало)
153
5 10 15 20 25 30 б
Рис. 3. Визуализация итераций решения и зависимость г(г) для К = 22 и г о = 0.01: б - зависимость г(г) от номера итерации
(окончание)
Табл. 1 - 4 отображают значения характеристик полостей и нормы г( г) на каждом шаге алгоритма решения.
Таблица 1
Значения характеристик полости с квадратным сечением и нормы г(г) для К = 11 и г о = 0.01
№ а 7 х 0 У 0 г № а 7 х 0 У 0 £
1 0,05 0 -0,2 -0,2 0,5282 20 0,225 62,5 0,025 0,025 0,0508
2 0,05 50 -0,2 -0,2 0,317 21 0,225 56,25 0,025 0,025 0,0477
3 0,15 50 -0,2 -0,2 0,1618 22 0,225 56,25 0,025 -0,0125 0,0364
4 0,15 50 -0,2 -0,05 0,0999 23 0,225 50 0,025 -0,0125 0,0337
5 0,15 50 -0,2 0,1 0,2046 24 0,225 37,5 0,025 -0,0125 0,0307
6 0,15 75 -0,2 0,1 0,1601 25 0,225 31,25 0,025 -0,0125 0,0315
7 0,15 75 -0,05 0,1 0,1109 26 0,225 31,25 0,025 0,025 0,0311
8 0,15 75 -0,05 -0,05 0,0899 27 0,2125 37,5 0,025 -0,0125 0,0304
9 0,25 75 -0,05 -0,05 0,0834 28 0,2125 40,625 0,025 -0,0125 0,0292
10 0,25 75 0,1 -0,05 0,0739 29 0,2125 40,625 0,025 0,0062 0,0284
11 0,25 50 0,1 -0,05 0,0634 30 0,2125 43,75 0,025 0,025 0,0341
12 0,25 50 0,1 -0,05 0,0634 31 0,2125 40,625 0,025 0,025 0,032
13 0,25 50 0,025 -0,05 0,0486 32 0,2125 40,625 0,0062 0,025 0,032
14 0,2 50 -0,05 -0,05 0,052 33 0,2125 40,625 0,0062 0,0062 0,0288
15 0,2 50 0,025 -0,05 0,0409 34 0,2125 39,0625 0,025 0,0062 0,0284
16 0,2 50 0,025 -0,05 0,0409 35 0,2125 39,0625 0,0156 0,0062 0,0274
17 0,225 50 0,025 -0,05 0,0381 36 0,2062 37,5 0,0062 0,0062 0,027
18 0,225 56,25 0,025 -0,05 0,0378 37 0,2062 39,0625 0,0062 0,0062 0,0269
19 0,225 56,25 0,025 -0,0125 0,0364 38 0,2 40,625 -0,0031 0,0062 0,0267
Таблица 2
Значения характеристик полости с квадратным сечением и нормы е(г) для К = 22 и е0 = 0.01
№ а 7 х 0 У 0 е № а 7 х 0 У 0 е
1 0,05 0 -0,2 -0,2 0,6085 18 0,2 56,25 -0,0125 -0,0125 0,0619
2 0,05 50 -0,2 -0,2 0,4605 19 0,175 62,5 -0,05 0,025 0,0769
3 0,15 50 -0,2 -0,2 0,2412 20 0,175 56,25 -0,05 0,025 0,0678
4 0,15 50 -0,05 -0,2 0,2403 21 0,175 56,25 -0,05 -0,0125 0,0645
5 0,15 50 -0,05 -0,05 0,0982 22 0,2 50 -0,0125 -0,0125 0,0589
6 0,25 25 0,1 0,1 0,1381 23 0,2 50 0,025 -0,0125 0,0569
7 0,25 25 0,1 -0,05 0,107 24 0,225 43,75 0,0625 -0,0125 0,0618
8 0,2 50 -0,05 -0,05 0,0841 25 0,225 43,75 0,025 -0,0125 0,0585
9 0,2 62,5 -0,05 -0,05 0,0829 26 0,2125 50 0,025 -0,0125 0,0561
10 0,2 62,5 0,025 -0,05 0,0666 27 0,2125 46,875 0,025 -0,0125 0,0553
11 0,25 62,5 0,1 -0,05 0,083 28 0,2125 43,75 0,025 -0,0125 0,0555
12 0,25 62,5 0,025 -0,05 0,0739 29 0,2125 46,875 0,025 -0,0125 0,0553
13 0,2 50 0,025 -0,05 0,0655 30 0,2125 45,3125 0,025 -0,0125 0,0552
14 0,25 37,5 0,025 -0,05 0,0801 31 0,2125 45,3125 0,025 -0,0125 0,0552
15 0,25 50 0,025 -0,05 0,0775 32 0,2125 46,0938 0,025 -0,0125 0,0552
16 0,2 56,25 0,025 -0,05 0,065 33 0,2125 46,0938 0,025 -0,0125 0,0552
17 0,2 56,25 -0,0125 -0,05 0,0642
■0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8
а
Рис. 4. Визуализация итераций решения и зависимость е^) для К = 11 и е0 = 0.005: а - визуализация итераций алгоритма решения
(начало)
155
5 10 15 20 25 30 35 40 45
б
Рис. 4. Визуализация итераций решения и зависимость г(г) для К = 11 и г0 = 0.005: б - зависимость г(г) от номера итерации
(окончание)
■0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8
а
5 10 15 20 25 30 35
б
Рис. 5. Визуализация итераций решения и зависимость г(г) для К = 22 и г0 = 0.005: а - визуализация итераций алгоритма решения; б - зависимость г(г) от номера итерации
Таблица 3
Значения характеристик полости с квадратным сечением и нормы е^) для К = 11 и е0 = 0.005
№ а 7 *0 У0 е № а 7 *0 У0 е
1 0,05 0 -0,2 -0,2 0,5221 24 0,2 9,375 -0,05 -0,0125 0,0236
2 0,05 50 -0,2 -0,2 0,3114 25 0,2 9,375 -0,0312 -0,0125 0,0216
3 0,15 50 -0,2 -0,2 0,1548 26 0,2 9,375 -0,0312 0,0062 0,0205
4 0,15 50 -0,2 -0,05 0,0976 27 0,2125 6,25 -0,0125 0,025 0,0252
5 0,15 50 -0,2 0,1 0,2127 28 0,2125 9,375 -0,0125 0,025 0,0231
6 0,15 75 -0,2 0,1 0,1665 29 0,2125 9,375 -0,0125 0,0062 0,0183
7 0,15 75 -0,05 0,1 0,1116 30 0,225 12,5 0,0063 0,0062 0,0179
8 0,15 75 -0,05 -0,05 0,086 31 0,225 12,5 0,0063 -0,0125 0,0176
9 0,25 75 -0,05 -0,05 0,0774 32 0,2375 18,75 0,025 -0,0312 0,0266
10 0,25 75 0,1 -0,05 0,0736 33 0,2375 18,75 0,0063 -0,0312 0,0262
11 0,25 50 0,1 -0,05 0,0619 34 0,2375 18,75 0,0063 -0,0125 0,0221
12 0,25 50 0,1 -0,05 0,0619 35 0,225 15,625 0,0063 -0,0125 0,017
13 0,25 50 0,025 -0,05 0,0385 36 0,225 15,625 0,0063 0,0062 0,0164
14 0,2 50 -0,05 -0,05 0,0451 37 0,2125 18,75 0,0063 0,025 0,02
15 0,2 37,5 -0,05 -0,05 0,0424 38 0,2125 21,875 0,0063 0,025 0,0183
16 0,2 37,5 -0,05 0,025 0,0375 39 0,2125 21,875 -0,0125 0,025 0,0179
17 0,2 25 -0,125 0,1 0,0969 40 0,2125 21,875 -0,0125 0,0062 0,015
18 0,2 12,5 -0,125 0,1 0,0842 41 0,2 25 -0,0312 0,0062 0,0201
19 0,2 12,5 -0,05 0,1 0,0497 42 0,2 25 -0,0125 0,0062 0,0178
20 0,2 12,5 -0,05 0,025 0,0252 43 0,2 25 -0,0125 0,025 0,0176
21 0,2 12,5 -0,05 -0,0125 0,0238 44 0,2125 18,75 -0,0125 0,0062 0,0144
22 0,2 6,25 -0,05 -0,05 0,0316 45 0,2125 18,75 -0,0125 0,0062 0,0144
23 0,2 6,25 -0,05 -0,0125 0,0241
Таблица 4
Значения характеристик полости с квадратным сечением и нормы е^) для К = 22 и е0 = 0.005
№ а 7 *0 У0 е № а 7 *0 У0 е
1 2 3 4 5 6 1 2 3 4 5 6
1 0,05 0 -0,2 -0,2 0,6166 20 0,2 6,25 -0,05 -0,05 0,048
2 0,05 50 -0,2 -0,2 0,4701 21 0,2 6,25 -0,05 -0,0125 0,0352
3 0,15 50 -0,2 -0,2 0,2273 22 0,2 6,25 -0,05 -0,0125 0,0352
4 0,15 50 -0,2 -0,05 0,1333 23 0,2 6,25 -0,0312 -0,0125 0,0339
5 0,25 75 -0,2 0,1 0,2504 24 0,2125 6,25 -0,0125 -0,0125 0,0301
6 0,25 75 -0,05 0,1 0,1446 25 0,2125 9,375 -0,0125 -0,0125 0,0293
7 0,25 75 -0,05 -0,05 0,1076 26 0,2188 9,375 -0,0125 -0,0125 0,0293
8 0,25 75 0,1 -0,05 0,0916 27 0,2188 10,9375 -0,0125 -0,0125 0,0292
9 0,25 50 0,1 -0,05 0,0778 28 0,2188 10,9375 -0,0125 -0,0031 0,0284
10 0,25 50 0,1 -0,05 0,0778 29 0,2188 12,5 -0,0125 0,0062 0,0293
11 0,25 50 0,025 -0,05 0,0584 30 0,2188 14,0625 -0,0125 0,0062 0,0292
12 0,2 50 -0,05 -0,05 0,0636 31 0,2188 14,0625 -0,0031 0,0062 0,0281
Окончание таблицы 4
1 2 3 4 5 6 1 2 3 4 5 6
13 0,2 37,5 -0,05 -0,05 0,0634 32 0,2188 14,0625 -0,0031 -0,0031 0,0279
14 0,2 37,5 -0,05 0,025 0,0528 33 0,2188 18,75 0,0062 -0,0031 0,0288
15 0,2 25 -0,125 0,1 0,1351 34 0,2188 18,75 -0,0031 -0,0031 0,0285
16 0,2 12,5 -0,125 0,1 0,1155 35 0,2188 18,75 -0,0031 0,0062 0,0282
17 0,2 12,5 -0,05 0,1 0,0721 36 0,2188 13,2812 -0,0031 -0,0031 0,0279
18 0,2 12,5 -0,05 0,025 0,0379 37 0,2188 13,2812 -0,0031 -0,0031 0,0279
19 0,2 12,5 -0,05 -0,0125 0,0369
Рис. 2, б - 5, б и значения e(z) из табл. 1 - 4 и наглядно показывают, как уменьшение ошибки e 0 измерения потенциала Y0 приводит к уменьшению погрешности определения геометрических параметров полости.
Анализ результатов показывает, что предложенная схема измерений может быть использована для идентификации параметров z полости квадратного сечения с удовлетворительной точностью. Точность определения характеристик полости зависит от количества датчиков K, в которых осуществляется измерение значений потенциала Yo, и ошибки измерения потенциала Yo - e o.
Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований и Правительства Тульской области (проект № 16-41-710083).
Список литературы
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. Иванов В.И., Скобельцын С.А. Моделирование задачи идентификации положения полости в упругом препятствии по рассеянному звуковому полю // Известия Тульского государственного университета. Естественные науки. 2011. Вып. 3. С. 74 - 86.
5. Толоконников Л.А., Ходюшина Е.В. Определение радиуса концентрической полости упругой сферы по известному рассеянному акустическому полю // Известия Тульского государственного университета. Естественные науки. 2015. Вып. 3. С. 211 - 218.
6. Хук Р., Дживс Т.А. Прямой поиск решения для числовых и статических проблем, 1961. С. 212—219.
7. LiveLink for MATLAB User's Guide. Stockholm: COMSOL AB, 2015. 318 с.
Скобельцын Сергей Алексеевич, канд. физ.-мат. наук, доцент, [email protected], Россия, Тула, Тульский государственный университет,
Пешков Никита Юрьевич, аспирант, инженер-программист, [email protected], Россия, Тула, Developer Express Inc.
DETERMINATION OF THE GEOMETRIC PARAMETERS OF AN ELASTIC CYLINDER. 'S CA VITY USING THE SCA TTERED ACOUSTIC FIELD
S.A. Skobel'tsyn, N.Y. Peshkov
The method of determining the geometric parameters (a center's position, a side length, an orientation angle) of an elastic elliptic cylinder's cavity, which have a square cross section, by using the known scattered field of the plane harmonic acoustic wave has been developed. Verification of the scheme for identification of the cavity's parameters in a numeric experiment has been carried out. The influence of a measuring instruments' error on the accuracy of the geometric characteristics' determination has been evaluated. Explained algorithm is universal and can be applied to identify an arbitrary count of an elastic body's parameters.
Key words: cavity with square cross section, scattered field, displacements potential, sensors, Hooke-Jeeves algorithm, measurements error.
Skobel 'tsyn Sergey Alekseevich, candidate of physico-mathematical sciences, docent, [email protected], Russia, Tula, Tula State University,
Peshkov Nikita Yurievich, postgraduate, software engineer, [email protected], Russia, Tula, Developer Express Inc.