Конечноэлементное моделирование конструкций из неоднородных материалов с усложнёнными свойствами Текст научной статьи по специальности «Математика»

УДК 539.3+534.1 й01: 10.12737/2015

Конечноэлементное моделирование конструкций из неоднородных

и *

материалов с усложненными свойствами

П. А. Оганесян, А. Н. Соловьёв

(Донской государственный технический университет)

Рассмотрены два алгоритма конечноэлементного моделирования неоднородных свойств, в том числе при связности механических и электрических полей. В первом случае неоднородные механические и пьезоэлектрические свойства тела задаются в виде аналитических зависимостей в некоторых классах функций. Во втором предполагается, что функции, описывающие неоднородные свойства, известны в дискретном наборе точек. Неоднородные свойства учитываются при вычислении элементных матриц. Для этого во втором случае интерполяция неоднородных свойств внутри тела осуществляется двумя способами: с помощью сплайнов, реализующих модель тонкой упругой пластины, и сплайнов Ж. Менге. Предложенные алгоритмы реализованы в виде модулей конечноэлементного комплекса ACELAN. Разработанное программное обеспечение может быть использовано при решении статических и динамических задач для неоднородных упругих и электроупругих тел. Приведены соответствующие примеры.

Ключевые слова: функционально неоднородные упругие и электроупругие материалы, сплайны, метод конечных элементов.

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

В современных исследованиях выделяются два класса функционально неоднородных материалов. В первом случае неоднородность возникает при поляризации керамики, когда электродное покрытие имеет сложную геометрию [1]. В работе [2] рассматривается неоднородность другой природы. В неоднородном материале имеется трещина, которая подвергается структурному анализу при помощи метода Галёркина — Петрова. Вблизи трещины все свойства изменяются по одному закону (предполагается увеличение в е0,5 раз). Исключение — тепловые коэффициенты, которые просто экспоненциально убывают по мере приближения к трещине.

В настоящей статье разработан модуль конечноэлементного комплекса АСБЦАЫ для моделирования функционально неоднородных материалов, в том числе пьезоэлектрических керамик с неоднородной поляризацией. При этом функция распределения неоднородных свойств материала, зависящая от положения точки в теле, предполагается известной в виде либо функциональной зависимости, либо дискретного набора значений величин соответствующих физических свойств в определённых узлах. Описана методика интерполяции этих свойств на используемую конечноэлементную сетку. Приведены примеры решения задач для тел с неоднородными свойствами.

* Работа выполнена при финансовой поддержке РФФИ (гранты 13-01-00196-а, 13-01-00943-а, 13-01-12006-а).

Постановка задач для тел с неоднородными свойствами. Математическая модель составного упругого, электроупругого и акустического тела с неоднородными свойствами твёрдых тел состоит из краевой задачи [3]: — для упругих и электроупругих тел

рми + а^рр - V-а = f], V-D = 0, (1)

а = cEj -(е + р,;е)-eT • E, D + ^D = e; -(е + ^£) + £** • E, (2)

£ = (vu + VuT)/2, E = -Уф; (3)

— для акустической среды

—p + V • v = 0; v = Уф, pV = V • а ; а = -pi + bVv. (4)

p c

j

Здесь а — тензор напряжений; pJ — плотность тела; e — тензор деформаций; u — вектор перемещений; D — вектор электрической индукции; E — вектор напряжённости электрического поля; fj — вектор массовых сил; ф — электрический потенциал; adj, pdj, ;d — коэффициенты демпфиро-

E T S «

вания; cE, eT, £aS — тензоры упругих констант, пьезомодулей и диэлектрических проницаемо-

стей; индекс j обозначает номер тела в модели. В уравнениях для акустической среды: c — скорость звука; v — вектор скоростей; ф — потенциал скоростей; p — звуковое давление, I — единичный тензор.

К уравнениям (1)—(4) добавляются механические, электрические и акустические граничные условия. Тензоры cE, eTj, eaSj и плотность тел будем считать функциями от положения точки в теле:

p* = Ppk (x), cE = cE (X), £aS = £aS (x), eT = eT (x). (5)

Для учёта степени поляризации предлагается использовать следующую зависимость механических, электрических свойств материалов:

— для тензоров cEj и eaS — g = gi + || (ga - gi);

T , (6)

— для тензора eT — g = | ga.

Здесь через g обозначены компоненты соответствующих тензоров. При этом индексом /обозначены компоненты тензоров для изотропного состояния, а индексом а — для анизотропного. Модуль вектора поляризации || влияет на то, какие свойства будут проявляться в большей

степени. Причём тензор пьезомодулей eTj будет нулевым для изотропных тел.

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

Для блока восстановления функций реализованы несколько алгоритмов построения сплайнов. Во-первых, это thin-plate spline [4], отвечающий модели бесконечно протяжённой пластины, деформируемой в некотором наборе точек изгибом. К решению задачи нас приведёт минимизация функционала, описывающего полную свободную энергию изогнутой упругой пластины. Вид этого функционала:

J2 ( Ф )=П

( Л2

Уф

dx2

+ 2

52ф dxdy

'aV2

5y2

dxdy ^ min. (7)

Решение такой задачи ф(xi, y i) = fi, i = 1,N имеет общий вид

N

ф(x,y) = YjCirf Inr2 + CN+1 + CN+2X + CN+3y ■

i=1

Здесь

r =(x - xi )2 +(y - yi )2

(8)

(9)

Чтобы определить коэффициенты с нужно решить следующую систему линейных алгебраических уравнений:

N

ф (xj У ) = Z cirj ln 4 + CN+1 + CN+2XJ + CN+зУj = fj , i=1

= (xj -xi)2+(yj - yi)2, j=ъй;

rj =ix

SC- = 0; Zxci = 0; ZyCi = 0.

(10)

Эта система не является положительно определённой (на главной диагонали — нули). Для её решения используется метод Гаусса с выбором ведущего элемента. Плюс этого подхода — простота реализации. Однако существуют и более эффективные способы построения сплайнов. Рассмотрим теорию интерполяционных С-сплайнов, которые позволяют получить функции с заданной степенью гладкости. Возникает следующая задача на минимизацию функционала гладкости:

Jm (ф)^ min, 0 = {ф eW2m (Rn) |ф(t ) = f, i = 1n) , (11)

( - V/2

где Jm (ф) = V>

1*2(0)

JZ m (°"ф )2 dO

V O PI=m

(12)

Здесь V, / = 1,^ — узлы интерполяции; fi, / = — заданные числа (значения восстанавливаемой функции в узлах интерполяции).

Вывод этого уравнения приводится в книге О. В. Ашкеназы [4]. При условии, что количество точек интерполирования хотя бы вдвое превосходит заданную степень гладкости т (в рассматриваемом случае это выполняется практически всегда), существует единственное решение задачи (12), общий вид которого:

N

ф (t) = ZCiEm (t, t ) + p(t) , t 6 R" .

С точностью до положительного множителя:

(. ч"/2-1i \2rn-" ,||

-1) u -V In u -V , при чётном ", .(-1)'

Em (u,v) =

("-1)/2 i |2m-"

■ u -v

при нечётном

Второе слагаемое является многочленом степени m -1:

P(t) = Z^v (t)6 PPm-1 .

(13)

(14)

(15)

Причём фу (V), V = 1, М, М = С"п+т_1 — это различные одночлены степени, меньшей или равной т. Чтобы определить коэффициенты сплайна, нужно решить следующую систему уравнений:

N M _

ZCjEm (t,tj ) + Z^v (t ) = f, i = 1,N;

" ' v=1

N _

Zc^v (tj ) = 0, v = 1,M.

j=1

j=1

Эта система решается с помощью алгоритма Ж. Менге. Некоторое подмножество точек выделяется для построения системы фундаментальных многочленов. Её порядок зависит от выбранной степени гладкости сплайна. Таким образом исходная система разделяется на две: первая (меньшего размера) — для нахождения системы многочленов и основная — для нахождения сплайна. Основная матрица в этом методе является положительно определённой и симметричной, что позволяет использовать быстрый и устойчивый метод квадратного корня для её решения. Программная реализация в АСЕАN. Модуль построения сплайнов — отдельный исполняемый файл с собственным интерфейсом, позволяющим использовать его для различных целей. Связь между частями комплекса осуществляется при помощи набора утилит. Это позволяет заменять отдельные части пакета независимо друг от друга. Для генерации входных данных, описывающих неоднородные свойства, была создана специальная графическая оболочка. В ней можно задать неоднородности для конкретной модели, построенной в ACELAN, и визуализировать их. Кроме этого проводится проверка введённых пользователем данных на непротиворечивость физической модели материала.

Чтобы пользователю было удобно оценивать поведение восстановленной с помощью сплайнов функции, разработан визуализатор на основе тепловых карт. К тепловой карте прилагаются отметки о максимальном и минимальном значении функции — им соответствуют определённые цвета. Карта строится на основе конечноэлементной сетки, поэтому на границах элементов могут возникать артефакты изображения. Это происходит при относительно большой разнице значений функции между двумя элементами. На этот случай предусмотрена возможность сглаживания изображения на основе классического фильтра Гаусса. Все представленные утилиты реализуют программные интерфейсы для взаимодействия с данными в форматах ACELAN и снабжены инструментами оценки эффективности модулей — записывается время выполнения всех операций, сложных с точки зрения вычислений.

Пример визуализации неоднородности представлен на рис. 1.

Рис. 1. Визуализация неоднородных свойств

Прямые задачи для неоднородных тел. Рассмотрим работу описанного ПО на примере вдавливания штампа (на верхнюю сторону которого приложено равномерное давление) в упругую неоднородную полосу. Размеры штампа — 0,4 х 0,1 см. К верхней границе приложены равномерно распределённые нормальные напряжения (интенсивность — 103 Н/м2). Материал штампа — изотропный модуль Юнга 2,0 х 1014 Н/м2. Коэффициент Пуассона — 0,29. Размеры полосы — 3 х 0,5 см. Нижний край закреплён. Материал — цирконат-титанат свинца (Р27"-4). Упругие свойства полосы изменяются по закону д = д0 ехр(29х0, где д — компоненты тензора упругих постоянных неоднородного тела, д0 — соответствующие компоненты для однородно поляризованной керамики.

На рис. 2 представлена геометрия задачи. Цветом выделяется неоднородность упругих свойств и распределение вертикального смещения на деформированном состоянии области.

А

хз

б)

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

На рис. 3, 4 и 5 представлены расчёты по излучению акустических волн в воду (область серого цвета на рис. 3, изображена только часть) с помощью пьезоэлемента в виде полупассивного биморфа. Верхняя часть биморфа изготовлена из меди и имеет длину 0,125 м при толщине 0,011 м. Нижняя часть состоит из керамики ЦТБС-3. Её длина 0,1 м, толщина — 0,01 м. Левый край биморфа закреплён. На верхней и нижней грани керамической пластинки заданы электрические потенциалы, образующие разность в 100 В.

Рис. 3. Визуализация неоднородности по электроупругим свойствам

Рассматривался случай неоднородной поляризации пьезокерамической пластины, при котором левый край поляризован полностью, а правый не поляризован (см. рис. 3). Модуль вектора поляризации || в соотношениях (6) изменяется по закону || (х) = -9х + 0,1.

На рис. 4 представлено распределение потенциала скоростей в жидкой среде при частоте 10 кГц.

б)

Рис. 4. Поле потенциала скорости в жидкости для неоднородной среды: а — однородный случай, б — неоднородный

На рис. 5 изображена амплитудно-частотная характеристика (АЧХ) вертикального смещения крайнего левого верхнего узла пластины.

Частота, кГц б)

Рис. 5. Результаты решения задачи для излучателя в акустической среде. АЧХ крайнего правого верхнего узла излучателя: а — для однородной среды, б— для неоднородной среды

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

1. Белоконь, А. В. Математическое моделирование необратимых процессов поляризации / А. В. Белоконь, А. С. Скалиух. — Москва : ФИЗМАТЛИТ. — 2010. — 328 с.

2. Fracture analysis in continuously nonhomogeneous magneto-electro-elastic solids under a thermal load by the MLPG / J. Sladek [et al.] // International Journal of Solids and Structures. — 2010. — Vol. 47. — Pp. 1381-1391.

3. Белоконь, А. В. Новые схемы конечноэлементного динамического анализа пьезоэлектрических устройств / А. В. Белоконь, А. В. Наседкин, А. Н. Соловьёв // Прикладная математика и механика. — 2002. — Т. 66, № 3. — С. 491-501.

4. Ашкеназы, В. О. Сплайн-поверхности. Основы теории и вычислительные алгоритмы / В. О. Ашкеназы. — Тверь : Тверской гос. ун-т, 2003. — 82 с.

Материал поступил в редакцию 20.05.2013.

References

1. Belokon, A. V., Skaliukh, A. S. Matematicheskoye modelirovaniye neobratimykh protsessov polyarizatsii. [Mathematical simulation of polarization irreversible processes.] Moscow : FIZMATLIT, 2010, 328 p. (in Russian).

2. Sladek, J., et al. Fracture analysis in continuously nonhomogeneous magneto-electro-elastic solids under a thermal load by the MLPG. International Journal of Solids and Structures, 2010, vol. 47, pp. 1381-1391.

3. Belokon, A. V., Nasedkin, A. V., Solovyev, A. N. Novyye skhemy konechnoelementnogo dinamicheskogo analiza pyezoelektricheskikh ustroystv. [New diagrams of finite-element dynamic analysis of piezoelectric devices.] Prikladnaya matematika i mekhanika, 2002, vol. 66, no. 3, pp. 491-501 (in Russian).

4. Ashkenazy, V. O. Splayn-poverkhnosti. Osnovy teorii i vychislitelnyye algoritmy. [Spline surfaces. Theoretical framework and computational algorithms.] Tver : Tverskoy gosudarstvennyy universitet, 2003, 82 p. (in Russian).

FE SIMULATION OF STRUCTURES FROM NONHOMOGENEOUS MATERIALS WITH COMPLICATED PROPERTIES*

P. A. Oganesyan, A. N. Solovyev

(Don State Technical University)

Two finite element simulation algorithms for heterogeneous properties including the connectivity of mechanical and electrical fields are considered. In the first instance, heterogeneous mechanical and piezoelectric properties of the body are given in the form of analytical relationships in certain classes of functions. The second method assumes that the functions describing the heterogeneous properties are known in a discrete set of points. The changes caused by inhomogeneous properties are taken into account under the calculation of the element matrices. For this, in the second case, the interpolation of heterogeneous properties within the body is implemented in two ways: using thin-plate splines, and J. Meinguet's splines. The proposed algorithms are implemented as modules of the ACELAN finite element complex. Examples of using the developed software for the calculation of the stress strain state of inhomogeneous elastic and electroelastic bodies are presented.

Keywords: functionally inhomogeneous elastic and electroelastic materials, splines, finite-element method.

* The research is done with the financial support from RFFI (grants 13-01-00196-a, 13-01-00943-a, 13-01-12006-a).

21