CC BY
8
УДК 30-04
DOI: 10.14529/ mmp220205
СРАВНИТЕЛЬНЫЙ АНАЛИЗ ЭФФЕКТИВНОСТИ ПРОГРАММНОЙ РЕАЛИЗАЦИИ ПОЛУАНАЛИТИЧЕСКИХ МЕТОДОВ РАСЧЕТА ВОЛНОВЫХ ПОЛЕЙ В МНОГОСЛОЙНЫХ АНИЗОТРОПНЫХ КОМПОЗИТАХ
Е.В. Глушков1, Н.В. Глушкова1, М.В. Варелджан1
1 Кубанский государственный университет, г. Краснодар, Российская Федерация
Численно моделируется процесс возбуждения и распространения бегущих упругих волн в многослойных анизотропных пластинах. В основу алгоритмов, реализованных в программном комплексе, положены явные интегральные представления решения соответствующих краевых задач и выведенные из них асимптотические представления для дальней от источника зоны. В ближней зоне, в том числе и под источником, амплитудно-частотные характеристики волнового поля определяются с помощью численного интегрирования несобственных контурных интегралов. Рассматриваются три подхода к вычислению волновых полей; на основе численных примеров проводится сравнительный анализ их эффективности.
Ключевые слова: слоистый композит; нормальные моды; интегральные и асимптотические представления волновых полей; малозатратные алгоритмы.
Введение
Широкое использование композитных материалов делает актуальной задачу создания программных комплексов, предназначенных для анализа динамического поведения конструкций, изготовленных из таких материалов. В частности, проектирование и настройка SHM систем дистанционного ультразвукового мониторинга (SHM - Structural Health Monitoring, [1]) опирается на расчет волновых характеристик для инспектируемых тонкостенных анизотропных волноводных структур, таких как фюзеляж авиалайнера, стенки хранилища химреактивов, трубопровода, реактора и т.п. Для получения количественных (амплитудно-частотных) характеристик возбуждаемых в них бегущих волн, используются различные подходы, от классического модального анализа до конечно-элементной аппроксимации (МКЭ). Промежуточное положение в этом ряду занимает полуаналитический интегральный подход, базирующийся на явном интегральном представлении вектора смещений волнового поля u, возбуждаемого поверхностной нагрузкой q, приложенной в некоторой области П (рис. 1), через Фурье-символ K матрицы Грина рассматриваемой упругой слоистой структуры [2-4]. В настоящее время большое внимание уделяется созданию на этой основе программных комплексов с удобным пользовательским интерфейсом, ориентированных на развитие SHM технологий [5].
Явный вид матрицы K известен только для ограниченного круга однородных изотропных волноводов (слой, полупространство, двухслойные среды), поэтому, использование таких представлений для структур более сложного строения (многослой-ность, анизотропия, непрерывная стратификация и др.) требует также и разработки эффективных алгоритмов вычисления K. Основой для предлагаемых ниже методов использования интегрального представления являются уже имеющиеся алгоритмы
построения матрицы К для слоистых анизотропных волноводов [6,7]. Тестовые расчеты, проведенные на этой основе для волоконно-армированных композитных пластин, показали хорошее совпадение с результатами экспериментальных измерений возбуждаемых бегущих волн [8].
В настоящей статье на основе численных примеров для четырехслойного перекрестно-армированного композитного образца обсуждается сравнительная эффективность предлагаемых методов по вычислительным затратам и точности. Предварительно дается постановка задачи и краткое описание вывода основных соотношений.
нагрузки
1. Математическая модель
Рассматривается М-слойная упругая пластина с произвольной анизотропией слоев, к поверхности которой в области П приложена осциллирующая нагрузка ц(х,у)в-гш*, генерирующая поле установившихся гармонических колебаний и(х)е-гш4; ш = - круговая частота, f - частота (рис. 1). В частности, это может
быть композитная пластина толщины Н, изготовленная из волоконно-армированных трансверсально-изотропных слоев-препрегов, с приклеенным к поверхности круговым пьезоактуатором.
Вектор комплексной амплитуды смещений и удовлетворяет в каждом из слоев уравнениям линейной эластодинамики [9]. На внутренних границах раздела слоев выполняются условия непрерывности вектора смещений и и напряжений т. Внешние границы г = 0 и г = —Н свободны от напряжений за исключением области приложения нагрузки П, в которой т\х=0 = я(х,у). На бесконечности выполняются условия излучения, следующие из принципа предельного поглощения. Более подробное описание постановки задачи и вывода интегральных и асимптотических представлений решения можно найти в указанных выше работах [3-8] и цитируемой в них литературе.
Применение преобразования Фурье Т по горизонтальным координатам х и у позволяет получить явное представление решения в виде интеграла по контуру Г+, идущему в комплексной плоскости а вдоль положительной полуоси, отклоняясь от нее при обходе вещественных полюсов подынтегральной функции (п (рис. 2), с внутренним интегралом по угловой координате 7:
2п
и(х) = ——~т~г [ [ К (а, 7, Гу)е~1аг со^1~1р\1гу ас1а. (1)
(2п)2 } ] г+ 0
Рис. 2. Контур Г+ и развернутый контур Г = Г- + Г+, замкнутый в верхнюю полуплоскость
0 < 7,р < 2п;
Здесь x = (x, y, z) и использованы полярные координаты
х = г cos р, Q'i = a cos 7, г = \/х2 + у2, у = г sin р, Q'2 = a sin7 q¡ = v'n j • п;.
K = F[k] и Q = F[q] - Фурье-символы матрицы Грина k(x) и вектор-функции q(x, y).
Чтобы показать влияние анизотропии материала на зависимость характеристик возбуждаемых волн от направления распространения р, для расчетов выбран осесим-метричный источник: не завясящая от р нагрузка q = qc(r) = (qr, q^,qz), заданная в круговой области радиуса a (П: r < a). В изотропном случае такая нагрузка возбуждает осесимметричное волновое поле u = uc(r, z) = (ur,uip,uz), в то время как наличие анизотропии нарушает осевую симметрию, внося зависимость от р.
Здесь и далее нижний индекс 'с' указывает на запись в цилиндрических координатах. Связь между компонентами вектора в декартовых и цилиндрических координатах задается матрицей поворота C: u = Cuc, uc = C-1u,
(cos р — sin р 0 \
sin р cos р 0 , с-1(р) = C(—р), С-1(7)С(р) = С(р — y).
0 0 \j
В случае пленочного пьезоактуатора контактные напряжения концентрируются на границе области контакта, поэтому их действие обычно моделируется сосредоточенной нагрузкой qr = q0í(a — r), где q0 - амплитудный множитель, а $(a — r) - радиальная дельта-функция [10]. Фурье-символ такой нагрузки Qc(a) = (Qr, 0, 0), Qr = q0 iaJ1(aa), J1 - функция Бесселя.
В изотропном случае зависимость элементов матрицы K от 7 определяется сомножителями e±imY, m = 0,1, 2, что позволяет взять внутренние интегралы в явном виде, воспользовавшись представлением функции Бесселя [11]
2п
2nim J"(r) = / eÍ(r"^ (2)
0
Далее проводится разворот контура Г+ в полный контур Г = Г- + Г+ (рис. 2), его замыкание в верхнюю полуплоскость и замена интеграла суммой вычетов в полюсах Zп, попадающих в замкнутый контур. Вклад вычетов в вещественных и близких
к вещественным полюсах дает бегущие и слабозатухающие (псевдоповерхностные) волны, возбуждаемые нагрузкой q. Подробно техника их выделения из общего интегрального представления (1) как в изотропном, так и в анизотропном случае, описана в работах [3,4,6,7].
В анизотропном случае представление (1) для осесимметричного источника qc (г) принимает вид [4]:
ис(г, <р, £) = J Кс(а,(р,г)Цс(а)а(1ск, (3)
г+
в котором
2п
Кс = — ( С-1{^)К{а,1,г)С{1)егагсо^~1)(11. (4)
2п ] о
Произвольная зависимость К от 7 не позволяет проинтегрировать (4) аналитически и, воспользовавшись теорией вычетов, сразу получить асимптотическое представление в виде суммы бегущих волн. Тем не менее здесь также можно получить явное представление интеграла (3) в виде двойного ряда:
те
ис(г,^) = С-1Ы ^ ит(г,г)в-гт, (5)
1т
т=—те
^ 2п те
^^ г.
п=1 о
Функции Бесселя .]т(аг), обеспечивающие разворот контура и замену интеграла по а суммой вычетов, появляются здесь при использовании разложения экспоненты [11]
е—готс°8(7—^ = ^ .т(аг)егт(7—п/2). (6)
т=—те
Оба ряда по п и по т сходятся экспоненциально, первый - за счет убывания функций Ханкеля
а второй - ввиду убывания функций Бесселя как т—т в разложении экспоненты (6) при больших значениях порядка т:
1 еаг т
~ . - , т —> оо.
у } л/Ъгт V 2т) '
Однако представление (5) все-таки довольно громоздко и не дает непосредственного разложения по бегущим волнам, так как вычеты в каждом полюсе (п входят во все слагаемые ит. Более практичным и физически наглядным является асимптотическое разложения по бегущим волнам, полученное без перехода к функциям Бесселя и разворота контура [6,7]. Замена переменных
1) 7 = в + V + п/2 = в для п/2 <7 - 3п/2 и
2) а = -а, 7 = в + V — п/2 = в — п для |7 — < п/2
те
приводит интеграл (2) к виду
1
u
( / - ) K(a,e,z)Q(a,e)eiarsine(7)
(2п)2 , , , Г+ Г- 0
Во второй замене использовано также свойство K(—а, 7) = K(а, 7 + п) и Q(-а, 7) = Q(a,Y + п). В диапазоне интегрирования 0 < в < п синус в показателе экспоненты положителен, поэтому оба контура Г+ и Г- можно замкнуть вверх. При этом внутрь замкнутых контуров попадает не вся верхняя полуплоскость, а отдельно правая и левая четвертьплоскости (рис. 2). Поэтому после применения теоремы Коши в представлении для интеграла (7) кроме суммы вычетов остаются интегралы по мнимой полуоси [0,гто]. Однако их вклад большего порядка малости, чем у вычетов в вещественных полюсах, которые достаточно удержать для получения асимптотики бегущих волн в дальнем поле:
Nr
u(x) = Un (x) + Ud (x). (8)
n=1 n
Ura = ¿/ (9)
0
Rn = iZn res K (aAz)U=in Q(Zn,0)Zn (в), в = в + P + п/2, Sn (в) = Zn(0)sin e,
Nr - число вещественных полюсов, Ud - суммарный вклад вычетов в комплексных полюсах и интегралов по мнимой оси, причем un = O((Zr)-1/2), а ud = O((Zr)-1) при Zr ^ то (Z - характерное значение волнового числа).
Бегущие волны описываются асимптотикой интегралов (9), которая строится методом стационарной фазы [12]. Стационарные точки фазовой функции sn(e) определяются корнями уравнения
4(в) = Zn(0) cos в + zn(e) sin в = 0, (10)
которое эквивалентно уравнению ctg в = — Zn(в)/Zn(в). В изотропном случае Z^e) = 0, что приводит к уравнению cos в = 0, имеющему на отрезке [0,п] единственный корень в = п/2, не зависящий от направления p.
В анизотропном случае Zn(0) = 0, поэтому этот корень зависит от p, совпадая с п/2 только для главных направлений, когда Zn (в) = 0. Кроме того, кривые ctg в и —Zn/Zn могут пересекаться в нескольких точках, т.е. для определенных направлений p фазовое уравнение (10) может иметь несколько корней в'пт, m =1, 2,..., Mn. Вклад стационарных точек в'пт дает асимптотику
ura = J] a^p, z)eia™7\/CKl + O(Zr)"1), (r 00, (11)
m=1
p + n/2.
n
При вещественных (n слагаемые суммы (11) описывают квазицилиндрические бегущие волны, фазовые и групповые скорости которых cnm = u/snm и vnm = [dsnm/du]-1 определяются волновыми числами snm. В отличие от изотропного случая они зависят от направления распространения ip и совпадают с полюсами Zn только для главных направлений (вдоль главных осей материала). Число бегущих волн Mn, соответствующих одному полюсу Zn и их амплитудные множители anm также зависят от р.
При кусочно-постоянной зависимости упругих свойств от z общее решение уравнений (3) в области преобразования Фурье выписывается в каждом слое в явном виде с точностью до шести постоянных множителей tm = (^,tm , ■■■jtim), m = 1, 2,..., M. Вектор неизвестных коэффициентов t = (ti,t2,...,tM) длиной 6M определяется из системы линейных алгебраических уравнений
At = f, (12)
возникающей при удовлетворении граничных условий. Поэтому по построению определитель матрицы этой системы А = det(A) входит в знаменатель коэффициентов tm, а тем самым и в общий знаменатель элементов матрицы К: К = K(a1, а2, z)/A(a1, а2). Корни а = Zn характеристического уравнения
А(а1 ,а2,ш ) = А(а,7,ш) = 0 (13)
являются полюсами подынтегральной функции в представлениях (1) и (3) (числитель К, как и Фурье-символ нагрузки Q, полюсов не имеет). В анизотропном случае они зависят не только от частоты ш, но и от направления 7 в плоскости (а1,а2) : а =
2. Сравнительный анализ эффективности различных подходов 2.1. Особенности программной реализации
Предлагаются три основных способа использования интегрального представления (1):
1. Численное интегрирование по контуру Г+ в соответствии с представлением (3) с матрицей Грина (4), заданной интегралом по 7 (Int).
2. Использование представления в виде суммы вычетов (8), слагаемые которой задаются интегралами по угловой координате (9) (Res).
3. Явное аналитическое выражение (11) для бегущих волн (Asympt).
При этом наличие подпрограмм, реализующих алгоритмы построения матрицы K для произвольно анизотропных многослойных волноводов, позволяют работать с ее элементами как с известными функциями.
Каждый из этих подходов имеет свои плюсы и минусы, определяющие диапазон их применимости. Первый (Int) предполагает двукратное численное интегрирование по а и y и поэтому выглядит наиболее затратным. Для интегрирования по полубесконечному контуру Г+ используются хорошо зарекомендовавшие себя подпрограммы, реализующие численное интегрирование методом Симпсона с самонастраиваемым переменным шагом, который выбирается в зависимости от точности, достигнутой на предыдущем участке (Dinn4, Dinn5 и др.). Вычисление внутреннего интеграла по
Y менее затратно, так как отрезок интегрирования ограничен, а подынтегральная матрица-функция - гладкая. Он хорошо аппроксимируется квадратурными формулами метода Гаусса, причем требуемая точность достигается сравнительно небольшим числом гауссовых узлов, которое, однако, увеличивается по мере усиления осцилляции экспоненты е-гат cos(7-^) с ростом параметра ar.
В целом при умеренных значениях r (в ближней зоне) вычислительные затраты первого подхода вполне приемлемы, но с увеличением расстояния от источника они возрастают из-за усиливающейся осцилляции подынтегральной функции, оставаясь допустимыми при r меньше некоторого порогового значения r2. С другой стороны, асимптотическое представление (11) (Asympt) обеспечивает требуемую точность в дальней зоне r > ri, причем погрешность уменьшается с ростом r. На практике обычно диапазоны применимости подходов Int и Asympt перекрываются, и для r1 < r < r2 они дают одинаковый в пределах допустимой погрешности результат, взаимно контролируя корректность решения.
Второй подход (Res) занимает промежуточное положение. С одной стороны, аппроксимация интегралов по a суммой вычетов избавляет от численного интегрирования по полубесконечному контуру Г+. Остается только несколько интегралов (9), вычисление которых требует на порядок меньше затрат. С другой стороны, в рамках третьего подхода Asympt и эти затраты практически сходят на нет, так как интегралы (9) заменяются значениями подынтегральной функции в нескольких стационарных точках.
Следует, однако, учитывать, что взамен численного интегрирования появляется необходимость определения полюсов Zn и вычетов res K, то есть требуется организовать надежный поиск всех вещественных, а в некоторых реализациях, и ближайших комплексных корней характеристического уравнения (13) во всем диапазоне изменения y , а также разработать алгоритм вычисления вычетов для заданной численно матрицы K. Решение этих задач также требует больших затрат, причем не только вычислительных, но и аналитических и программистских.
Вычислительные затраты здесь можно сократить, учитывая, что Zn и res K не зависят от r, и для каждой частоты ш их значения можно найти и затабулировать один раз. Что касается вычетов, их вычисление не так сложно и затратно, как может показаться на первый взгляд, если воспользоваться конечно-разностным приближением
res K|a=Cn « (K|«=Cn+h - K|а=Сп_ь) h/2, (14)
дающим погрешность порядка O(h) при условии, что полюс Zn найден с точностью е на 2-3 порядка выше. Например, для нахождения вычетов с тремя верными значащими цифрами следует взять h = 10-4 при е = 10-6.
Слабым местом представления (14) является потеря точности вычисления K для значений a предельно близких к полюсу Zn, которая происходит из-за вырождения матрицы A системы (12) при обращении ее определителя в ноль при a = Zn. Этого недостатка нет у метода, предложенного в работе [4], который основан на разложении матрицы и векторов системы (12) в ряд Лорана в окрестности нуля определителя Zn, являющегося полюсом вектор-функции t(a):
A(a) = Ao + Ai(a - Zn) + •••, det Aq = 0, (15)
t(a) = t-1/(а - Zn) + to + t1 (а - Zn) + ■■■,
f(а) = fo + fl(а - Zn) + ■■■■
Искомый вычет res К |a=zn выражается через res t|a=zn = t-1, и задача сводится к определению этого вектора. Подстановка разложений (15) в систему (12) и приравнивание коэффициентов при одинаковых степенях (а — Zn) дает соотношения
Ao t-1 = 0, A11-1 + Aoto = fo
Первое из них означает, что t-1 - это собственный вектор матрицы Ao = A(Zn), который может быть определен стандартными методами с точностью до постоянного множителя c: t-1 = cm, где m - известный, нормированный определенным образом, собственный вектор. Неизвестный множитель c определяется из второго соотношения (16). Для этого достаточно скалярно домножить его на собственный вектор n сопряженной матрицы AQ . Соотношение принимает вид
c(A1m, n) + (Aoto, n) = (fo, n),
из которого в силу свойства сопряженной матрицы
(Aoto, n) = (to,AQn) = 0
однозначно фиксируется c = (fo, n)/(A1m, n).
При реализации данного алгоритма не возникает необходимости решения системы с плохо обусловленной матрицей, а конечно-разностная аппроксимация
A1 = A'(Zn) = [A(Zn + h) - A(Zn - h)]/(2h) численно устойчива при любых h.
(16)
2.2. Численные примеры
Для иллюстрации сравнительной эффективности описанных выше подходов ниже приводятся результаты расчетов для композитной пластины, составленной из волоконно-армированных слоев-препрегов с перекрестной схемой укладки [0°, 90°, 0°]: волокна двух внешних слоев ориентированы вдоль оси х, а двух внутренних - перпендикулярно к ней. Толщина препрегов к = 0, 56 мм, плотность р = 1522 кг/м3, общая толщина пластины Н = 4к = 2, 24 мм. Трансверсально-изотропные упругие свойства препрегов характеризуются упругими модулями (в обозначениях Фойгта [9], в ГПа):
С11 с12 с13 с22 с23 с33 с44 с55 с66
119,14 3,57 3,57 10,56 4,96 10,56 2,80 4,45 4,45
Пластина моделируется трехслойным пакетом слоев (М = 3) с толщинами к1 = к3 = к и к2 = 2к; источник колебаний (круговой пьезоактуатор) радиуса а = 8 мм -радиальной сосредоточенной нагрузкой дг = д0#(а — г).
На рис. 3 показаны зависимости вещественных полюсов (п от частоты f для данного образца, построенные в диапазоне до 0.8 МГц для двух значений угла 7: 7 = 0 и
3.5
3
2.5
5 2 S 2
1.5
О
1
0.5 0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 f, МГц f, МГц
Рис. 3. Частотные зависимости вещественных полюсов Zn для 7 = 0 и 7 = п/4
f=100 кГц
f=500 кГц
0 п
Зп/2
Зп/2
Рис. 4. Примеры угловых зависимостей полюсов при f = 0,1 и 0, 5 МГц
7 = п/4. В изотропном случае это были бы не зависящие от 7 дисперсионные кривые для волновых чисел бегущих волн Лэмба. Для рассматриваемой анизотропной пластины они, во-первых, формально не являются дисперсионными кривыми, так как в соответствии с представлениями (8) и (11) волновые числа впт хоть и выражаются через эти полюса, но не совпадают с ними. Во-вторых, вид кривых зависит от угла 7, что наглядно иллюстрируется круговыми диаграммами (п(7) для полюсов, дающих три первые (фундаментальные) моды: А0, 50 и 5Н0, на частоте f = 0,1 МГц и дополнительные высшие моды при f = 0, 5 МГц (рис. 4) (на рис. 3 эти частоты отмечены вертикальными пунктирными линиями).
Сравнительная эффективность разрабатываемых методов иллюстрируется графиками амплитуды радиальной и вертикальной компонент смещения |иг | и \иг | на нижней поверхности г = —Н в зависимости от расстояния г (рис. 5) и частоты f (рис. 6). Приводятся их абсолютные величины, отнесеннвге к нормирующей величине + \иг|2, вычисленной при г = 50 мм, f = 0,1 МГц и р = 0. Для контро-
uo
ur
ля достоверности на рис. 5 добавлены также кривые, полученные с помощью МКЭ пакета СОМБОЬ МиШрЬузюз (https://www.comsol.ru/) для области радиуса г0 = 100 мм, ограниченной с помощью поглощающих граничных условий РМЬ; максимальный размер ячейки сетки - 1, 5 мм, число элементов (степеней свободы) - 3 455 439. Для
0
r/u0
tp = 0 —COMSOL, 931 сек — Int, 387 сек д Res, 84 сек --Asympt, 22 сек
- v \
3 2.5 2 1.5 1
0.5 0
3 2.5 2 1.5 1
0.5 0
—COMSOL, 931 сек
(р = п/А — Int, 652 сек
a Res, 85 сек
■ - Asympt, 54 сек
I
—COMSOL, 931 сек| — Int, 371 сек Д Res, 82 сек Asympt, 24 сек
—COMSOL, 931 сек — Int, 206 сек a Res, 42 сек —Asympt, 11 сек
<р = 7г/4
—COMSOL, 931 се!
— Int, 341 сек д Res, 43 сек
— Asympt, 27 сек
—COMSOL, 931 сек ш = 7г/2 _Int, 197 сек VN г ' л Res, 41 сек
^ — Asympt, 12 сек
0 20 40 60 80 Г, ММ
100 о
20 40 60 80 100 Г, мм
u
u
0
z
Рис. 5. Иллюстрация сравнительной эффективности рассматриваемых методов: сопоставление зависимостей |иг| и \их| от расстояния г для трех направлений р = 0, п/4 и п/2 при / = 0,1 МГц
1.2
1
с „д 0.8
0.6
0.4
0.2
0
Я „ '/д V А <М г л ¡р = 0 —Int Д Res
--Asympt
д\
V ч
{ у/ | V
0.1
0.2
0.3 f, МГц
0.4
0.5
0.2 0.3 0.4 f, МГц
0
Рис. 6. То же, что и на рис. 5, но для зависимостей от частоты при г = 50 мм
сопоставления вычислительных затрат в легендах рис. 5 указано время счета каждой кривой на персональном компьютере с процессором AMD Ryzen 9 3900X, тактовая частота - 3,8 ГГц.
Здесь следует отметить, что для фиксированной частоты МКЭ (COMSOL) дает результат сразу для всех 600 точек rj, по которым строятся графики по r, в то
время как методы Int, Res и Asympt предполагают вычисление для каждой точки в отдельности (хотя здесь тоже используется определенная оптимизация, избавляющая от пересчета величин не зависящих от r, например, полюсов и вычетов в них). Тем не менее, время счета COMSOL все равно значительно превышает даже время самого затратного из трех подходов Int. На рис. 6 результаты COMSOL не приводятся, так как получение частотных зависимостей с помощью МКЭ требует многих часов или даже суток.
Анализ кривых на рис. 5 подтверждает сделанные выше предварительные выводы. В данных примерах точность вычисления радиальной компоненты (левый столбец графиков) выше, чем вертикальной. На этих графиках ясно видно наличие диапазона ri < r < r2, в котором совпадают результаты всех трех методов и COMSOL, в то время, как результаты Asymt отклоняются от остальных слева от этого диапазона, а Int - справа. Последнее - из-за того, что максимального числа узлов Гаусса (500), заданного для вычисления внутреннего интеграла по 7 в представлении (2), оказывается недостаточно для достижения требуемой точности при усиливающейся с ростом r осцилляции подынтегральной экспоненты. Точность асимптотических подходов Res и Asympt при вычислении uz в рассматриваемом диапазоне хуже, но и здесь наблюдается сходимость с ростом r, там где Int уже перестает работать.
Частотные зависимости на рис. 6 приводятся для не очень большого расстояния r = 50 мм, поэтому на низких и средних частотах заметно некоторое отклонение асимптотических результатов от результатов численного интегрирования, которое ожидаемо сходит на нет с ростом f.
Заключительные замечания и выводы
Дано сравнительное описание трех подходов к численному моделированию волновых процессов в слоистых анизотропных волноводах с заданным поверхностным источником колебаний, которые базируются на явных интегральных и асимптотических представлениях решения соответствующих краевых задач динамической теории упругости. Использование таких представлений позволяет не только существенно снизить вычислительные затраты по сравнению с широко используемым в инженерной практике методом конечных элементов (МКЭ), но и выделить физически наглядные представления для возбуждаемых бегущих волн. Анализ особенностей компьютерной реализации каждого из рассмотренных подходов (Int, Res и Asympt) и численных результатов, полученных как на их основе, так и с помощью МКЭ пакета COMSOL при моделировании ультразвуковых волн, возбуждаемых круговым пьезоактуатором в слоистой волоконно-армированной композитной пластине, позволяет сделать следующие выводы относительно их сравнительной эффективности и диапазонов применимости.
1. Int. Главное преимущество: применимость в ближней зоне, в том числе и в области приложения нагрузки (возможность рассчитывать динамическую реакцию, податливость, импеданс и т.п.); интегральное представление дает точное решение рассматриваемой задачи; точность численных результатов можно увеличить за счет увеличения времени счета.
Недостаток: сравнительно большие численные затраты, возникающие с ростом r или частоты ш; соответствующие ограничения диапазона применимости по r и ш.
2. Res. Значительное снижение затрат на численное интегрирование за счет появления дополнительной задачи поиска всех необходимых полюсов £n. Но неприменим в области приложения нагрузки r < a и некоторой ее окрестности (там, где вклад отброшенного слагаемого больше допустимой погрешности). Как и Int имеет ограничения сверху по r и w.
3. Asympt. Требует минимальных численных затрат (только на поиск полюсов Zn и вычисление res K). Дает явное представление для каждой моды бегущей волны, что позволяет не только определять их характеристики (скорость, амплитуду), но и исследовать распределение энергии источника между возбуждаемыми волнами и строить диаграммы направленности излучения энергии на бесконечность. Этим перекрываются возможности модального анализа, который также дает характеристики каждой бегущей волны, но только с точностью до постоянного множителя, т.е., без реальных значений амплитуды. Еще одно важное преимущество - явное выражение для групповой скорости vnm в заданном направлении р без необходимости предварительного определения угла отклонения вектора групповой скорости от волнового вектора, необходимого в рамках традиционной техники модального анализа.
Недостаток: неприменим в ближней зоне на расстояниях меньших нескольких длин волн от источника; требует больших аналитических и программистских затрат для надежной реализации.
Работа выполнена в рамках проекта FZEN-2020-0017 государственного задания Минобрнауки России.
Литература
1. Mitra, M. Guided Wave Based Structural Health Monitoring: a Review / M. Mitra, S. Gopalakrishnan // Smart Materials and Structures. - 2016. - V. 25, № 5. - P. 1-27.
2. Ворович, И.И. Динамические смешанные задачи теории упругости для неклассических областей / И.И. Ворович, В.А. Бабешко. - М.: Наука, 1979.
3. Бабешко, В.А. Динамика неоднородных линейно-упругих сред / В.А. Бабешко, Е.В. Глушков, Ж.Ф. Зинченко. - М.: Наука, 1989.
4. Глушкова, Н.В. Определение и учет сингулярных составляющих в задачах теории упругости: дисс. ... докт. физ.-мат. наук / Н.В. Глушкова. - Краснодар, 2000.
5. Глушков, Е.В. Возбуждение и распространение упругих волн в многослойных анизотропных композитах / Е.В. Глушков, Н.В. Глушкова, А.С. Кривонос // Прикладная математика и механика. - 2010. - Т. 74, № 3. - P. 419-432.
6. Glushkov, E. Forced Wave Propagation and Energy Distribution in Anisotropic Laminate Composites / E. Glushkov, N. Glushkova, A. Eremin // Journal of the Acoustical Society of America. - 2011. - V. 129, № 5. - P. 2923-2934.
7. Глушков, Е.В. Программный комплекс Waves-L для моделирования и визуализации волновых процессов в упругом слое / Е.В. Глушков, Н.В. Глушкова, С.И. Фоменко, А.А. Еремин, А.А. Евдокимов, О.И. Новиков // Вестник ЮУрГУ. Серия: Математическое моделирование и программирование. - 2019. - Т. 12, № 1. - С. 110-121.
8. Glushkov, E. Efficient Mathematical Representations for Computing the Forced Wave Dynamics of Anisotropic Laminated Composites / E. Glushkov, N. Glushkova, A. Eremin // CEAS Aeronautical Journal. - 2013. - V. 4, № 1. - P. 11-19.
9. Кристенсен Р. Введение в механику композитов / Р. Кристенсен. - М.: Мир, 1982.
10. Raghavan, A. Finite-Dimensional Piezoelectric Transducer Modeling for Guided Wave Based Structural Health Monitoring / A. Raghavan, C.E.S. Cesnik // Smart Materials and Structures. - 2005. - V. 14, № 6. - P. 1448-1461.
11. Бейтмен, Г. Высшие трансцендентные функции / Г. Бейтмен, А. Эрдейи. - М.: Наука, 1969.
12. Федорюк, М.В. Метод перевала/ М.В. Федорюк. - М.: Наука, 1977.
Евгений Викторович Глушков, доктор физико-математических наук, профессор, Институт математики, механики и информатики, Кубанский государственный университет (г. Краснодар, Российская Федерация), evg@math.kubsu.ru.
Наталья Вилениновна Глушкова, доктор физико-математических наук, профессор, Институт математики, механики и информатики, Кубанский государственный университет (г. Краснодар, Российская Федерация), nvg@math.kubsu.ru.
Михаил Владимирович Варелджан, аспирант, младший научный сотрудник, Институт математики, механики и информатики, Кубанский государственный университет (г. Краснодар, Российская Федерация), michael.vareldzhan.777@mail.ru.
Поступила в редакцию 27 июля 2021 г.
MSC 74J15 DOI: 10.14529/mmp220205
COMPARATIVE ANALYSIS OF SOFTWARE IMPLEMENTATION EFFICIENCY OF THE SEMI-ANALYTICAL METHODS FOR CALCULATING WAVE FIELDS IN MULTILAYER ANISOTROPIC COMPOSITES
E.V. Glushkov1, N.V. Glushkova1, M.V. Vareldzhan1
1 Kuban State University, Krasnodar, Russian Federation
E-mail: evg@math.kubsu.ru, nvg@math.kubsu.ru, michael.vareldzhan.777@mail.ru
We construct a numerical model of the excitation and propagation of traveling elastic waves in multilayer anisotropic plates. The algorithms implemented in the software package are based on explicit integral representations of the solution to the corresponding boundary value problems and the asymptotic representations derived from them for the zone farthest from the source. In the near-field zone, including under the source, the amplitude-frequency characteristics of a wave field are obtained by the numerical integrating of the improper path integrals. Three approaches to the calculation of wave fields are considered; a comparative analysis of their effectiveness is conducted based on numerical examples.
Keywords: layered composites; normal modes; integral and asymptotic representations of wave fields; low-cost algorythms.
References
1. Mitra M., Gopalakrishnan S. Guided Wave Based Structural Health Monitoring: a Review. Smart Materials and Structures, 2016, vol. 25, no. 5, pp. 1-27. DOI: 10.1088/09641726/25/5/053001
2. Vorovich I.I., Babeshko V.A. Dinamicheskie smeshannye zadachi teorii uprugosti dlya neklassicheskikh oblastej [Dynamic Mixed Problems of Elasticity for Nonclassical Domains]. Moscow, Nauka, 1979. (in Russian)
3. Babeshko V.A., Glushkov E.V., Zinchenko J.F. Dinamika neodnorodnyh linejno-uprugih sred [Dynamics of Inhomogeneous Linearly Elastic Media]. Moscow, Nauka, 1989. (in Russian)
4. Glushkova N.V. Determination and Accounting for Singular Terms in the Problems of Elasticity. PhD thesis, Krasnodar, 2010. (in Russian)
5. Glushkov E.V., Glushkova N.V., Fomenko S.I., Eremin A.A., Novikov O.I. Software Package Waves-L for Simulation and Visualization of Wave Processes in an Elastic Layer. Bulletin of the South Ural State University. Series: Mathematical Modelling, Programming and Computer Software, 2019, vol. 12, no. 1, p. 110-121. DOI: 10.14529/mmp190109
6. Glushkov Ye.V., Glushkova N.V., Krivonos A.S. The Excitation and Propagation of Elastic Waves in Multilayered Anisotropic Composites. Journal of Applied Mathematics and Mechanics, 2010, vol. 74, pp. 297-305. DOI: 10.1016/j.jappmathmech.2010.07.005
7. Glushkov E., Glushkova N., Eremin A. Forced Wave Propagation and Energy Distribution in Anisotropic Laminate Composites. Journal of the Acoustical Society of America, 2011, vol. 129, no. 5, pp. 2923-2934. DOI: 10.1121/1.3559699
8. Glushkov E., Glushkova N., Eremin A. Efficient Mathematical Representations for Computing the Forced Wave Dynamics of Anisotropic Laminated Composites. CEAS Aeronautical Journal, 2013, vol. 4, no. 1, p. 11-19. DOI: 10.1007/s13272-013-0064-1
9. Christensen R.M. Mechanics of Composite Materials. N.Y., Wiley-Interscience, 1979.
10. Raghavan A., Cesnik C.E.S. Finite-Dimensional Piezoelectric Transducer Modeling for Guided Wave Based Structural Health Monitoring. Smart Materials Structures, 2005, vol. 14, no 6. pp. 1448-1461. DOI: 10.1088/0964-1726/14/6/037
11. Bateman H., Erdelyi A. Higher Transcendental Functions. N.Y., Toronto London, McGraw-Hill Book Company, 1953.
12. Fedoryuk M. V. Metod perevala [Saddle-Point Method]. Moscow, Nauka, 1977. (in Russian)
Received July 27, 2021
