CC BY
28
УДК 539.3, 534.1
КОНЕЧНО-ЭЛЕМЕНТНОЕ МОДЕЛИРОВАНИЕ И ИДЕНТИФИКАЦИЯ НЕОДНОРОДНЫХ МАТЕРИАЛОВ В ACELAN
© 2013 г. П.А. Оганесян
Оганесян Павел Артурович - аспирант, кафедра теоретической и прикладной механики, факультет машин и оборудования агропромышленного комплекса, Донской государственный технический университет, пл. Гагарина, 1, г. Ростов н/Д, 344000, e-mail: wolwerine@yandex.ru.
Oganesyan Pavel Arturovich - Post-Graduate Student, Department of Theoretical and Applied Mechanics, Faculty of Machinery and Equipment of Agro-Industrial Complex, Don State Technical University, Gagarin Sq., 1, Rostov-on-Don, 344000, Russia, e-mail: wolwerine@yandex.ru.
Рассмотрено создание программного комплекса для определения неоднородных свойств конструкций из пъезоматериалов. Исследуется модель электроупругого тела с произвольной неоднородностью в континуальной и конечно-элементной постановке. Для решения коэффициентных задач предлагается применить генетический алгоритм. Изучены варианты использования некоторых стохастических алгоритмов. Для решения прямых задач применяется пакет АСЕЬАК. Реализованы программные интерфейсы между всеми частями комплекса. Приведены некоторые тестовые задачи. Получена оценка точности и эффективности при некоторых наборах параметров генетических алгоритмов.
Ключевые слова: неоднородные пьезоматериалы, обратные задачи, конечно-элементное моделирование.
The development of software for modeling of inhomogeneous properties of peizomaterials is observed in this paper. The model of electroelastic inhomogeneous body in analytical and finite-element formulation is discussed. We propose to use universal genetic algorithm to find coefficients. Possibility of using selected stochastic algorithms is studied. Modules ofpackage ACELAN are applied as direct problem solvers. Implemented software interfaces between finite element solvers and genetic algorithm. Input information can be presented by the set of eigen values or by the values of displacement at several points. The examples ofproblems solved with the help of the described software are presented. The accuracy and effectiveness of certain combinations of the parameters of genetic algorithms are researched.
Keywords: inhomogeneous piezomaterials, inverse problems, finite element modeling.
Существует множество примеров использования пьезоэлектрических преобразователей с неоднородными свойствами. Природа неоднородности может быть различной: неоднородные механические и электрические свойства материала или неоднородная поляризация, вызванная внешними факторами или конструкцией прибора. В монографии [1] разработаны методы определения неоднородных свойств электроупругих стержней, в работах А.О. Ва-тульяна и его учеников [2-5] - методы идентификации неоднородных свойств упругих тел. В данной работе представлен программный комплекс, позволяющий определять неоднородности в конструкциях из пьезоматериалов. Результаты, полученные с помощью комплекса, могут быть использованы для оптимизации приборов и конструкций из пьезома-териалов, оценки свойств путем неразрушающего анализа, поиска дефектов.
Постановка прямых и обратных задач
Комплекс состоит из трех основных частей. Для решения прямой задачи применяются разработанные модули программного комплекса ACELAN, позволяющие моделировать неоднородные свойства электроупругих материалов. Используемые математические модели базируются на краевых задачах линейной теории электроупругости и уравнений движения жидких и газообразных сред в акустическом приближении. Модель может состоять из нескольких тел. Введем индекс ] для нумерации тел и запишем уравнения и определяющие соотношения, которые описывают поведение исследуемого тела [5].
Для упругих и электроупругих тел Ррк» + - V-С = I] , V-Б = 0 , (1)
Е
о = c
J
(г + ßj) -ej • E , D + gdD = ej ■■ (г + gdг) + э> ■ E ,
£ = ^и + Vuг )/2, Е = где с - тензор напряжений; р^ - плотность тела; £ -тензор деформаций; и - вектор перемещений; Б -вектор электрической индукции; Е - вектор напряженности электрического поля; ^ - вектор массовых
сил; ф - электрический потенциал; а^, Р], дл -
коэффициенты демпфирования; сЕ, е ], э ] - тензоры упругих констант, пьезомодулей и диэлектрических проницаемостей; индекс ] отвечает номеру тела в модели. К уравнениям (1) добавляются механические, электрические и акустические граничные условия.
Тензоры сЕ, еТ , э^ и плотность тел будем считать
функциями от положения точки в теле:
рк = Ррк (х), сЕ = сЕ (х), э] = э ] (х), е ] = е] (х). (2)
На границах твердых тел предполагается задание механических и электрических граничных условий. На внутренних границах ставятся условия непрерывности соответствующих механических и электрических величин. На границах акустических сред предполагаются условия, соответствующие свободной границе, жесткой стенке, или импедансные условия, моделирующие бесконечную среду. Зависимость механических и электрических свойств материалов от координат (2) может возникнуть за счет неоднородной поляризации. Этот случай учитывается следующим образом:
для тензоров cE и э J
- g=g' + И (- g'),
для тензора е] - g = . (3)
Через g обозначены компоненты соответствующих тензоров, при этом индексом i обозначены компоненты тензоров для изотропного состояния, а индексом a - для анизотропного. Модуль вектора поляризации р| влияет на то, какие свойства будут проявляться в большей степени, причем тензор пьезомоду-лей еТ будет нулевым для изотропных тел.
В случае обратной задачи одна или несколько зависимостей, описанных в (2), считаются неизвестными. Для решения такой задачи требуется дополнительная информация. Ею может служить набор собственных частот модели или измеренные перемещения на некоторых участках поверхности модели. В случае использования собственных частот задача сводится к минимизации функционала:
S (Ф) =
1 --
Ф
Ф*
(4)
где Ф = (ф,Ф2,...,Фд,) - рассчитываемые частоты;
л*
Ф - эталонные частоты, измеренные в эксперименте; N - число рассматриваемых частот. В работе натурный эксперимент заменен численным и для получе-
*
ния Ф нужно решить прямую задачу с заданными неоднородными свойствами материала.
Разработка инструментов для решения задач
Предлагаемый алгоритм решения обратной задачи состоит из нескольких этапов. Первый из них - это построение модели прямой задачи в конечно -элементной постановке. Далее нужно построить функционал на основе информации о задаче, позволяющий оценивать решение прямой задачи на каждом шаге процедуры минимизации. Затем выбираются класс функций, в котором будет осуществлен поиск, и начальное приближение для каждого коэффициента искомой функции. В работе для минимизации функционала (4) используется генетический алгоритм (ГА)
2
[7], поэтому последним этапом подготовительных работ является выбор для него стратегии поиска.
Для решения прямой задачи используется конечно-элементная модель, реализованная в комплексе ACELAN. На основе первых двух уравнений из (1) запишем конечно-элементную модель задачи в векторной форме [6]: u(x, t) = NTU (x) • U(t) , tfx, t) = N^(x) • Ф^),
W(x, t) = N^ (x) • T(t) , где N u - матрица функций формы для поля перемещений; N^ - вектор функций
формы для электрического потенциала; N^ - вектор функций формы для потенциала скоростей в акустической среде. U(t), Ф(Х), ¥(t) - глобальные векторы соответствующих узловых степеней свободы.
Проведем замену а=[и,Ф,¥]т, теперь систему линейных алгебраических уравнений (СЛАУ) можно записать в виде M • a + C • a + K • a = F, где
f
M =
K =
мии 0 R t
0 R ищ 00 0 - M.
K - K
up
PP
0
0
0 0
- K
C =
Cuu
^d K^p R T
R ищ
R
0 0
0 - C
ищ 0
щщ)
w )
Fu
Fp+?d Fp| 0
Матрицы масс (Muu), демпфирования (Cuu) и жесткости (Kuu) симметричны и неотрицательны. Матрица Киф определяет пьезоэлектрические свойства, Кфф -диэлектрические, она симметрична и неотрицательно определена. Именно на этапе построения этих матриц происходит учет неоднородных свойств.
Комплекс ACELAN снабжён модулем задания неоднородностей, в котором пользователь может описать вид функции (2) с помощью набора коэффициентов для одного из классов функции (экспонента, полином, функция Хевисайда) либо задать набор значений в наборе точек тела, по которому будет восстановлена функция неоднородности. Имеется возможность проанализировать полученные неоднородности с помощью тепловых карт. Пользовательские данные проходят верификацию на непротиворечивость физическим законам модели. Для восстановления функций используются thin-plate сплайны и гладкие сплайны, получаемые по методу Ж. Менге [8]. Основным отличием между этими двумя методами является возможность задать степень гладкости восстанавливаемой функции при использовании метода Менге. Перейдем к дискретному описанию функции неоднородности для возможности применять стохастические алгоритмы. Однако при решении реальных задач следует учитывать количество параметров, используемых при такой дискретизации. Количество вычисляемых коэффициентов значительно влияет на эффективность ГА. Это налагает ограничения на поиск решения в виде полинома или сплайна. В случае, когда тело состоит из нескольких неоднородных материалов, возможны два способа решения. В первом подходе можно составить из нескольких наборов коэффициентов один вектор для создания популяции в ГА. Во втором - использовать итерационный алгоритм, в котором можно зафикси-
ровать значения всех неоднородных функций, кроме одной, и решать задачи последовательно, находя более точные значения по каждой из неоднородных функций. Выбор способа зависит от степени осведомленности исследователя о характере неоднород-ностей и связях между ними. Например, если известно, что неоднородные свойства пьезоматериала связаны не с его структурой, а с неоднородной поляризацией, то, используя соотношения (3), можно искать только одну векторную функцию P(x), описывающую неоднородность.
Для решения обратной задачи выбран модуль, реализующий ГА, описанный в [7]. Для применения его к конкретной задаче предусмотрена возможность подключения исполняемого файла, реализующего решение прямой задачи. В работе объектом популяции ГА выступает набор коэффициентов функции в заданном классе или множество значений этих функций в наборе точек в теле. При этом во втором случае размер популяции является переменным: вначале минимизация производится на грубой сетке, далее происходит уточнение. Необходимые протоколы передачи данных реализованы в соответствующей утилите на языке C# с применением фреймворка .Net 4.0. Это накладывает ограничения на использование комплекса под некоторыми операционными системами, но позволяет использовать широкий арсенал инструментов оптимизации и диагностики программы. Данная утилита производит запуск «решателя» прямой задачи, предварительно вычислив значения неоднородных функций для каждого элемента либо узла тела. Получив решение в виде набора собственных частот или таблицы смещений, утилита строит оценку предложенного решения на основе заранее выбранной нормы. После получения оценок для всех элементов текущей популяции формируется новая. Характер обновления зависит от параметров ГА, которые можно задавать в интерфейсе программы либо с помощью специального файла. Наиболее эффективные комбинации параметров алгоритма подбираются эмпирически.
Так как данные о задаче подаются в «решатель» в виде набора файлов, многочисленные вызовы «программы-решателя» для одной задачи, которые производит ГА, могут приводить к ошибкам доступа. Поэтому были разработаны дополнительные инструменты, обеспечивающие свободное чтение и запись всех необходимых файлов вне зависимости от того, какой из процессов комплекса в них нуждается. Для устранения ошибок доступа используется механизм поиска выполняемого процесса с заданным именем из библиотеки System.Diagnostics. Также добавлено измерение времени работы основных процессов, что позволило оптимизировать работу комплекса на этапе исследования тестовых задач.
В качестве альтернативного способа, не требующего аналитической формулировки процесса получения решения, рассматривался метод дифференциальной эволюции, описанный в [9]. В его основе стоит формирование новой популяции с использованием стохастического подхода - исследователь выбирает лишь один параметр метода (F), лежащий в диапазоне [0; 2] и отвечающий за скорость изменения популя-
F
ции. На каждом шаге метода выбираются три случайных элемента популяции, которые сформируют новый элемент. Метод хорошо себя зарекомендовал на многих классах задач [10]. Всего один входной параметр значительно облегчает работу, однако на практике количество решений прямой задачи возросло по сравнению с ГА, описанным выше, что негативно отразилось на времени решения задачи. Так как время, необходимое для работы конечно-элементной части комплекса, значительно превосходит время работы всех остальных модулей, количество решений прямой задачи является наиболее важным фактором, влияющим на общее время исследования. При использовании метода дифференциальной эволюции была реализована возможность сокращения числа решений прямой задачи при помощи сохранения значений целевой функции для особей, которые не изменялись при переходе к следующей популяции.
Примеры использования комплекса
Рассматривается задача о собственных колебаниях цилиндрического пьезокерамического (Р2Т-4) элемента, нижняя и верхняя поверхности которого полностью электродированы (рис. 1).
Значения 1,7653Е+000
8,7953Е-001 Ц
Рис. 2. Визуализация неоднородности в первом примере
Предположим, что, имея представление о характере неоднородности, необходимо вычислить коэффициенты с из формулы (5). Определим собственные частоты такого преобразователя с учетом (6). С помощью модуля ACELAN получаем набор собственных частот и собственных форм. На этом этапе исследования важно определить те частоты, которые наиболее заметно изменяются при данной неоднородности. В представленных численных экспериментах использовались 5 частот.
Передача данных в модуль ГА осуществляется с помощью утилиты. Настройка ГА производится при помощи графического интерфейса.
Поскольку все параметры подбираются эмпирически, нужно проводить несколько этапов численного эксперимента. В качестве оценки точности решения используется величина шах|-
S = -
max
N
(7)
Приведем результаты двух показательных численных экспериментов (табл. 1)
Таблица 1
Сравнение суперэлитной и классической стратегий для модели (рис. 1)
Показатель Суперэлитная Классическая
Область поиска [0, 20]х[-10, 0]х [0, 20]х[-10, 0]х
х[0,5, 5] х[0,5, 5]
Количество эпох 10 300
Размер популяции 10 20
Решение 9,5725; -0,2810; 7,9687; -1,0133;
0,849 0,94189
Относительная 8,5 3
ошибка (5), %
Рис. 1. Схема нагрузки
На электродах заданы электрические потенциалы, радиус пластины - 10-1 м, толщина - 10-2 м. Рассматривается цилиндр, неоднородный по плотности (например, при неоднородной пористости).
В первом примере (рис. 2) в качестве искомой функции неоднородности выберем р(х, у) = Ро/{х, у), Дх, у) = с X + С2У + Сз , (5)
где х, у - радиальная и осевая координаты.
Пусть эта функция описывает множитель, который будет изменять значение плотности относительно табличной для данного типа керамики. В тестовом примере выберем функцию f в виде
I (х, у) = 8х - 7 у +1. (6)
Стратегии отличаются параметрами (табл. 2)
Таблица 2
Характеристики стратегий
Стратегия Pc Pm Pi Pd Pt
Суперэлитная 0,85 0,8 0,1 0,5 0,01
Классическая 0,85 0,5 0,1 0,1 0,01
В табл. 2 Pc - вероятность кроссовера; Pm - мутации; Pi - инверсии; Pd - деструкции; Pt - транслокации.
Как видно из табл. 1, применялись две стратегии, характеристики которых представлены в табл. 2.
Высокие вероятности мутации и деструкции обеспечивают быстрое обновление наборов особей при суперэлитной стратегии. Это дает значительный выигрыш во времени. Важно также учитывать недетерминированность ГА, в силу которой два численных эксперимента с одинаковыми входными данными могут давать незначительно отличающиеся решения и использовать различное число эпох. Наименее точно рассчитан второй коэффициент, отвечающий за изменение плотности по оси у. В данной задаче толщина пластинки существенно меньше длины, поэтому вклад этого коэффициент в решение весьма мал. Такая ситуация может часто возникать при решении реальных задач, так как точно определить класс функции возможно далеко не всегда.
Далее рассматривается задача, в которой все искомые коэффициенты являются значащими, т.е. при полном соответствии искомой функции заданному шаблону. Во втором примере выберем искомую функцию в другом классе, оставив прежней геометрию задачи:
={2;х; ?2, (8)
с1 = 0,6, с2 = 0,05, с3 = 2,4 .
В качестве неоднородного свойства снова выбрана плотность.
Для оценки ошибки вычисляется расстояние между точками трехмерного пространства, образованного набором коэффициентов. Проведя аналогичные численные эксперименты, получаем следующие результаты (табл. 3).
Таблица 3
Сравнение суперэлитной и классической стратегий для модели (рис. 1)
Эксперимент, описанный в первой колонке табл. 3, предполагал, что исследователь знает, какая из частей модели обладает меньшей плотностью, что было отражено при выборе областей поиска. В случае применения классической стратегии поиска этого допущения не потребовалось - при одинаковых областях поиска для сI и с3 решение найдено с приемлемой точностью и правильно определено положение границы, разделяющей области с различными свойствами.
Выводы
Полученный комплекс можно использовать для решения обратных коэффициентных задач для электроупругих тел. Предлагается следующий алгоритм решения: определив тип неоднородности и возможные границы допустимых значений коэффициентов, провести решение по суперэлитной стратегии. Полученные данные следует использовать для формирова-
ния новых границ поиска для запуска классической стратегии с большим числом популяций и особей в популяции. Так можно добиться достаточно высокой точности. К недостаткам комплекса можно отнести его сегментированность, так как разные части реализованы на разных языках программирования, время их взаимодействия достаточно велико. Предполагается ускорение комплекса путем формирования единого исполняемого файла без внешних зависимостей и последующего использования параллельных вычислений [11]. Возможно также исследование эффективности восстановления искомой функции в виде ряда.
Литература
1. Ватульян А.О., Соловьев А.Н. Прямые и обратные задачи для однородных и неоднородных упругих и электроупругих тел. Ростов н/Д, 2008. 176 с.
2. Vatul'yan A.O. The theory of inverse problems in the linear mechanics of a deformable solid // J. of Applied Mathematics and Mechanics. 2010. Vol. 74, № 6. P. 648 - 653.
3. Ватульян А.О., Дударев В.В. О реконструкции неоднородных свойств пьезоэлектрических тел // Вычислительная механика сплошных сред. 2012. Т. 5, № 3. С. 259 - 264.
4. Бочарова О.В., Ватульян А.О. Обратные задачи для упругого неоднородного стержня // Изв. вузов. Сев.-Кавк. регион. Естеств. науки. 2008. № 3. С. 33 - 37.
5. Ватульян А.О., Половодова А.А. Определение неоднородной поляризации пьезокерамического стержня по данным акустического зондирования // Изв. вузов. Сев.-Кавк. регион. Естеств. науки. 2013. № 4. С. 21 - 25.
6. Белоконь А.В., Наседкин А.В., Соловьев А.Н. Новые схемы конечно-элементного динамического анализа пьезоэлектрических устройств // ПММ 2002. Т. 66, № 3. С. 491 - 501.
7. Баранов И.В., Ватульян А.О., Соловьев А.Н. Об одном генетическом алгоритме и его применении в обратных задачах идентификации упругих тел // Вычислительные технологии. 2006. Т. 11, № 3. С. 14 - 26.
8. Ашкеназы В.О. Сплайн-поверхности: Основы теории и вычислительные алгоритмы. Тверь, 2003. 82 с.
9. Storn R., Price K. Differential Evolution - A Simple and Efficient Heuristic for Global Optimization over Continuous Spaces // J. of Global Optimization. 1997. № 11. P. 341 - 359.
10. Storn R. System Design by Constraint Adaptation and Differential Evolution. Technical Report TR-96-039. ICSI, November 1996. URL: ftp.icsi.berkeley.edu (дата обращения: 14.09.2013).
11. Соловьев А.Н., Гервич Л.Р., Штейнберг Б.Я., Скали-ух А.С., Наседкин А.В., Сумбатян М.А. Параллельные решатели для пакета ACELAN : свидетельство о государственной регистрации программы для ЭВМ № 2011617952 от 11 октября 2011 г.
Показатель Суперэлитная Классическая
Область поиска [0,1, 1]х[0, 0,1] *[1, 5] [0,1, 10]х[0, 0,1] х[0,1, 1]
Количество эпох 10 300
Размер популяции 10 20
Решение 0,5756; 0,0531; 2,491 0,6344; 0,0529; 2,575
Относительная ошибка вычисления коэффициентов, % 8 10
Поступила в редакцию 7 октября 2013 г.
