УДК 539.3
Ю. И. Димитриенко, А. Н. Морозов, А. П. Соколов, Е. С. Ничеговский
МОДЕЛИРОВАНИЕ ЭФФЕКТИВНЫХ ПЬЕЗОЭЛЕКТРОУПРУГИХ КОМПОЗИЦИОННЫХ МАТЕРИАЛОВ
Предложен конечно-элементный метод решения локальных задач теории пьезоупругости на ячейке периодичности для композитов со сложными пространственными структурами армирования. Дано теоретическое обоснование метода. Разработан программный комплекс для решения локальных задач и расчета эффективных пьезоупругих характеристик композиционных материалов. Проведенное сравнение результатов расчетов с аналитическим решением для слоистого композита показало высокую точность разработанного метода.
E-mail: [email protected]
Ключевые слова: композиционные материалы, пьезоэлектроупругость,
метод асимптотического осреднения, метод конечного элемента.
В последнее время, особенно в связи с развитием нанотехнологий, резко возрос интерес к разработке методов моделирования новых материалов с необычными свойствами. Широкими возможностями для создания таких методов обладает метод асимптотического осреднения (МАО), называемый также методом гомогенизации [1-3]. В работах [4-7] описаны методы численного решения так называемых локальных задач, которые возникают в МАО, и, вообще говоря, являются интегродифференциальными. В настоящей работе предложено развитие этого метода для моделирования композиционных материалов с пьезоупругими свойствами.
Задача на ячейке периодичности. Рассмотрим композиционный материал, обладающий периодической структурой, ячейка периодичности (ЯП) V которого состоит из N фаз Va, а = 1 ,...,N. Введем малый параметр к = l/L ^ 1, где l — характерный размер ЯП, L — характерный линейный размер всего композита. Введем декартовы координаты xk, глобальные координаты Xk = xk/L и локальные координаты = xk/к. Рассмотрим для такой структуры задачу пьезо-упругости, предполагая выполненными условия идеального контакта на поверхности раздела между компонентами:
0;
dtj + Pe =0;
^j Cjjßki + vkij ek; £ij ui,j + j ;
(1) (2) (3)
(4)
di = — efe; efe = ;
[Ui] = 0, [^ij ] n, = 0 на Eij;
[<^] =0, [di] ni = 0 на Ei
ij
u | ^ Ue, Cij | ^ nj Si
dini I v dn 1
(5)
(6)
(7)
(8) (9)
(10)
Здесь (1) — уравнения равновесия, (2) — уравнение Гаусса, (3) — соотношения упругости с учетом пьезоупругого эффекта, (4) — соотношения Коши, (5) — электрические соотношения с учетом обратного пьезоэффекта, (6) — выражение для электрического потенциала, (7) — условия идеального механического контакта на границе раздела компонентов, (8) — соотношения идеального электрического контакта на границе раздела, (9) — граничные кинематические и статические условия на частях внешней поверхности £и и , (10) — электрические граничные условия, О^* — компоненты тензора модулей упругости, £¿7 — компоненты тензора малых деформаций, м^* — компоненты тензора пьезоэлектрических коэффициентов, а^ — компоненты тензора напряжений, э^- — компоненты тензора диэлектрических постоянных, е — компоненты вектора напряженности электрического поля, ^ — компоненты вектора электрической индукции; ^ — электрический потенциал.
Согласно МАО решение задачи пьезоупругости строится в виде асимптотических разложений:
и=и(0)(х*) + киг(1)(х*, <г) + к2 ...; = 40)(хк 1) + «^(х* 1) + к2 ...; = а(0)(х* ,<г) + «^(х*1) + к2 ...; р = ^(0)(х*) + к^(1)(х*,<) + к2 ... .
(11)
Подставляя (11) в (10), получаем в нулевом приближении локальную задачу пьезоупругости
j = 0;
-d(0) = 0-
di/i = °>
^(0) = C _(0) + v e(0).
-(0) _ -
ij (0)
ij
,(1) , „,(!)
= £,', +--U / + U . /.
ij 2 V i/ j/i
d(0) = „ ^(0) _ э e(0).
di = ^ifeefe .
(0)
e(0) = w(1) + ^.
(12)
k(1)] = 0, <U(1)> = 0;
(0)
n, = 0, [^(1)] = 0, [d(0)] ni = 0, <^(1)> =0; [[u(1)]] =0; [[^(1)]] =0.
1
Здесь Ец = (и(Ц + и) /2 — компоненты тензора осредненных деформаций композита и ек = еко) — компоненты осредненного вектора напряженности электрического поля (в локальной задаче они являются заданными величинами); [[ ]] — условие периодичности функций на границе ЯП; ( ) — операция осреднения на ЯП; , I = д/дхХ и / = д/д£1 — производные по двум типам координат; £^ав — поверхность контакта двух компонентов композита в пределах отдельно взятой ЯП, т.е.
££ав = ^ ^а.
Преобразование локальной задачи к задачам "классического"
типа. Задача содержит 9 входных функций Ец и , будем искать ее решение в виде суммы по этим функциям:
3 3 3 3
и(1) = ^ иг(рд) + ^ ^ = ^2 ^(г) + ^(Р^ (13)
р,д=1 г=1 г=1 р>9=1
где
иг(рд) = —Ерд (4р + $р) + ^(р^Ь ф(т) = — бг + (14)
а и^(рд) и Ф(г) — новые неизвестные функции. Подставляя (13) и (14)
ЛЛ, (0) (0) (0) 7(0)
в (12), находим выражения для е(0) , ек , ац и а\ :
3 3
р(0) = V + V •
fcj Z^ fcij(pq) + fcij(r);
p,q=1 r=1
j«) _ 1 + Uj(pq)/i)
p,q
ei°(r) = 1 + uj(r)i);
e(0) _ V e(0) + V e(0) • (15)
r=1 p,q=1
e(0) _ xT/ • e(0) _ .. •
33
.(0) _ V^ a0) + v^ J0) •
p,q=1 r=1
d(0) _ V d(0) + V d(0)
r=1 p,q=1
Тогда задачу (12) можно разделить на 9 отдельных связанных задач теории пьезоэлектроупругости:
L
„j у = 0;
„(0) = C _(0) + v e(0) ;
„ij(pq) Cijk1 fcfei(pq) + Vkij ek(pq);
(0)
pq
-d(0) )/. = 0
i(pq)/i
d(0) = v Д0) _ _ e
di(pq) Vik1fcfe1(pq) ek(pq)
Lr
„ij(r)/j = 0;
(16)
„(0) = C _(0) + v e(0)
„ij(r) = Cijk1 fcfe1(r) + Vkij efe(r)
d
(0) ; r);
(0)
i(r)/i )
i(r)
= 0;
d(0) = _(0) _ э e(0) di(r) = vifeifcfe1(r) efe(r).
С помощью представления решения в форме (13)-(15) уравнения (16) можно сформулировать не для всей ЯП, а для 1/8 ее части, граничные условия для этих систем уравнений на 1/8 ЯП имеют следующий вид.
Для задач : при р = д
Цфд) = (£рд/2)4Р, на Е£, иг(рд) = 0 на Ег;
£?(р«) = 0, £*(рд) = 0 на Ег и Е^ (г = = к = г); при р = д
иг(рд) = (ёгр/4)^гр на Е7, Ц^) = 0 на Е^;
(17)
jq) = 0, Ua(pq) = 0 на E, и Ej (i,j = {p,q});
j(pq)
а
i(pq)
k(pq)
та
■j(pq)
Sapq) = 0, S"(pq) = 0, Ua(pq) = 0 на Efe и E'fc (i = j = k = i);
при p = q ^(pp)|Sp = 0, ^(pp)|s^ = 0,5ep;
пРи p = q P(pq)|q |s =0, ^(pq)|q | =0.
(18) (19)
Для задач Lr :
Ui(r) = (err/2)^ir на Ei, Ui(r) = 0 на Ei; Sj(r) = 0, Sfe(r) = 0 на Ei и Ei (i = j = k = i);
^(r)!s = 0; ^(r)|s, = 0,5er.
(20) (21)
В формулах (17)-(21) S,(pq) = „,7(pq)n (i = j) — компоненты векторов усилий; Ei = = 0} — координатные поверхности; Ei = {6 = 1/2} — внешние поверхности 1/8 ЯП.
Вычисление эффективных характеристик композита. После решения локальных задач Lpq, Lr могут быть вычислены осредненные поля напряжений и электрической индукции
<„(0'> = £<„SU> + Вj>>; -<di0,> = £«',> + Е<0-
p,q r r p,q
(22)
Соотношения между этими осредненными полями и средними полями деформаций и напряженности электрического поля образуют искомые эффективные определяющие соотношения композита
3 3 3 3
<„(.)> ^ ^ Сijpqepq + ^ ^ Vrijer; <d( )> ^ ^ Vipqepq ^ ^ эirer. p,q=1 r=1 p,q=1 r=1
(23)
Из соотношений (23) находим компоненты Cijk1 эффективного тензора модулей упругости, э, — тензора диэлектрических постоянных и компоненты эффективного тензора пьезоэлектрических коэффициентов композита vijk:
Сijpq = <„SL) >->=0; эir = -у <d(2 =0
pq r
33
1 ' (0) \ , . V /„(0)
Vrij = F j) >-- =0 + E <„ST) >-r=0 ; (24)
1
/ . - t \ >-r =0 I / , - , , p,q r=1
33
Vipq = ~ | E H0) ^ =0 ^ <d((0)q) >-q =0
e
pq
\Г=1 Р,9=1
Здесь (а^))ер=0 — осредненные напряжения в задаче при одной ненулевой компоненте и всех нулевых компонентах ер, а (а(0) )ёт.=0 — осредненные напряжения в задаче при одной ненулевой компоненте £рд и одной ненулевой компоненте ег. Остальные обозначения аналогичны.
Применение метода конечных элементов для решения локальных задач. Вариационная формулировка локальных задач имеет вид
I ¿V = у Ш^Е, /¿еМ ¿V = 0, (25)
V, Ея К
где обозначены координатные строки из неизвестных функций: ит =
= [и^), и2(м), иЧт)] — перемещений; ат = [а^р^ а202}рд) ... а(0зЦ - напряжений; ет = [е(п)(рд),е2<2)(Р5)... ^(м) ] - деформаций; 8т = = — поверхностных сил; ет=
,(о) e(0) e(0)
— электрической напряженности; dт = ¿("Р^, ¿2^), ¿3^) — электрической индукции.
Аналогично формулируется вариационная постановка задачи Ьг. Запишем в матричном виде соотношения Коши в (15): е = ЮИ, где Ю — матрица линейных дифференциальных операторов, а также определяющие соотношения в (16): а = Се + Уе, d = Уе — Эе, где С
— матрица 6 х 6 из компонент тензора модулей упругости компонентов композита [8], У — матрица 6 х 3 коэффициентов пьезоупругости, Э — матрица 3 х 3 коэффициентов диэлектрической проницаемости, а е = где — дифференциальный оператор.
Аппроксимируем перемещения И и электрический потенциал р в КЭ линейными функциями координат: И = Фq, р = фту, где я и у — координатные столбцы перемещений и электрических потенциалов в узлах КЭ; Ф и фт — матрица и столбец функций формы. Тогда из вариационного уравнения (25) получаем итоговую разрешающую систему
линейных алгебраических уравнений: К я = I, где К = J В тС В дУ,
г V
Ке = В^ЭВедУ — матрицы жесткости, В = ЮФ, Ве = Юеф;
с = (С У);В= (В Ве);(Я):* = (/);(26)
/ — столбец нагрузок и /е — столбец электрической индукции для КЭ.
Составляя глобальную СЛАУ для всей 1/8 ЯП и решая ее, находим перемещения я и электрические потенциалы у в узлах. Для решения СЛАУ применялись методы сопряженного градиента. Для генерации конечно-элементных сеток, решения задач, проведения расчетов эффективных характеристик и визуализации результатов использовались программные технологии, разработанные на кафедре "Вычислительная математика и математическая физика" МГТУ им. Н.Э. Баумана.
Решение тестовой задачи для слоистого композита. В качестве тестовой задачи был рассчитан двухслойный композит, для которого известно точное аналитическое решение локальных задач (12) [9]. Объемная концентрация слоев равна 0,5. В качестве 1-го слоя взяты константы турмалина, обладающего тригональным (В-ромбоэдрическим) классом симметрии с группой анизотропии 3т
(0(8 по классификации [8]) (табл. 1), в качестве 2-го слоя взяты характеристики а-кварца (табл. 2), также обладающего тригональным классом симметрии, но с группой анизотропии 32 (0(8).
Таблица 1
Значения независимых компонент тензора модулей упругости Суиг, ГПа, тензора пьезоэлектрических коэффициентов, vijk, Кл/м2, и тензора диэлектрических постоянных э^, 10-10 Ф/м, для турмалина
С1111 С1122 С1133 С1123 С3333 С2323 V131 v222 v331 v333 Э11 Э33
270 69 8,8 -7,8 161 67 0,247 -0,018 0,103 0,32 0,598 0,489
Таблица 2
Значения независимых компонент тензора модулей упругости Суиг, ГПа, тензора пьезоэлектрических коэффициентов Vijk, Кл/м2, и тензора диэлектрических постоянных э^, 10-10 Ф/м, для а-кварца
С1111 С1122 С1133 С1123 С3333 С2323 V111 v123 Э11 Э33
86,8 6 11,9 -17,8 101 57,8 0,171 0,04 0,4 0,41
Значения эффективных пьезоэлектроупругих характеристик двухслойного композита композита, вычисленные по разработанному методу и полученные на основе аналитического решения [9], приведены в табл. 3. Имеет место достаточно хорошее совпадение результатов — относительная ошибка не превышает 4%.
Таблица 3
Значения независимых компонент эффективных тензоров модулей упругости Суи, ГПа, пьезоэлектрических коэффициентов Суи, Кл/м2, и диэлектрических постоянных Эу, 10-10 Ф/м, для двухслойного композита турмалин/а-кварц
Вид решения С1111 ©1133 ©111 V123 ©11 ©33
Аналитическое 178 10,8 0,087 0,022 0,501 0,448
Численное 171 10,5 0,09 0,027 0,496 0,444
Расчет для ортогонально армированного волокнистого композита. Решены локальные задачи для волокнистого композита, состоящего из матрицы, армированной волокнами, ориентированными по трем координатным направлениям. В качестве матрицы выбран турмалин, характеристики волокон соответствуют а-кварцу.
На рис. 1-4 приведены распределения напряжений ^О),^ в 1/8 ЯП, полученные в результате решения задач Ьт. На рис. 5 показано
Рис. 1. Распределение напряжений о^^, МПа, в 1/8 ЯП в задаче L11
Рис. 2. Распределение напряжений ^(2(12)' МПа, в 1/8 ЯП
Рис. 3. Распределение напряжений ^(3(13), МПа, в 1/8 ЯП
Рис. 5. Распределение компоненты ¿1(1) вектора электрической индукции в 1/8 ЯП композита
(0)
распределение компоненты ¿1(1) вектора электрической индукции в 1/8 ЯП композита в задаче Ь1. В табл.4 приведены значения компонент эффективных тензоров модулей упругости, пьезоэлектрических коэффициентов и диэлектрических постоянных композита.
Выводы. Разработан метод решения локальных задач пьезоупру-гости на ячейке периодичности для композитов со сложными структурами армирования, основанный на переходе от задачи для полной ЯП к существенно более простой задаче на 1/8 ЯП. Разработан конечно-элементный метод расчета эффективных пьезоэлектроупругих харак-
Таблица 4
Значения независимых компонент эффективных тензоров модулей упругости С^и, ГПа, пьезоэлектрических коэффициентов , Кл/м2, и диэлектрических постоянных Зу, 10-10 Ф/м, для трехмерного ортогонально-армированного композита турмалин/а-кварц
Сии С1122 С1133 С3333 С111 С123 С333 С131 Э11 С33
123 28,9 12,6 124,6 0,155 0,024 0,102 0,079 0,465 0,437
теристик композиционных материалов. Проведенное сравнение результатов расчетов с аналитическим решением для слоистого композита продемонстрировало высокую точность разработанного метода расчета эффективных пьезоупругих характеристик композитов. Показано, что разработанный метод может быть использован для решения локальных задач и расчета эффективных характеристик композитов со сложными пространственными структурами армирования.
СПИСОК ЛИТЕРАТУРЫ
1.Бахвалов Н. С., Панасенко Г. П. Осреднение процессов в периодических средах. - М.: Наука, 1984. - 352 с.
2. Санчес-Паленсия Э. Теория колебаний и неоднородные среды. - М.: Мир, 1984. - 472 с.
3. П о б е д р я Б. Е. Механика композиционных материалов. - М.: Изд-во МГУ - 1984.-356 с.
4. Димитриенко Ю. И. Механика композиционных материалов при высоких температурах. - М.: Машиностроение, 1997. - 367 с.
5. Димитриенко Ю. И., Кашкаров А. И. Конечно-элементный метод для вычисления эффективных характеристик пространственно-армированных композитов // Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки. -2002. - № 2. - С. 95-108.
6. Димитриенко Ю. И., Кашкаров А. И., Макашов А. А. Конечно-элементный расчет эффективных упругопластических характеристик композитов на основе метода асимптотического осреднения // Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки. - 2007. - № 1. - С. 102-116.
7. Димитриенко Ю. И., Соколов А. П. Разработка системы автоматизированного вычисления эффективных упругих характеристик композитов. // Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки. - 2008. - № 2. -С. 57-67.
8. Димитриенко Ю. И. Тензорное исчисление. - М.: Высш. шк. - 2001. -575 с.
9. Каралюнас Р. И. Эффективные термопьезоэлектрические свойства слоистых композитов // Механика композитных материалов. - 1990. - № 5. - С. 823830.
Статья поступила в редакцию 17.05.2010
Юрий Иванович Димитриенко родился в 1962 г., окончил в 1984 г. МГУ им. М.В.Ломоносова. Д-р физ.-мат. наук, профессор, заведующий кафедрой "Вычислительная математика и математическая физика" МГТУ им. Н.Э. Баумана, действительный член Академии инженерных наук. Автор более 160 научных работ в области вычислительной механики, нелинейного тензорного анализа, термомеханики композитов, математического моделирования в материаловедении.
Yu.I. Dimitrienko (b.1962) graduated from the Lomonosov Moscow State University in 1984. D. Sc. (Phys.-Math.), professor, head of "Computational Mathematics and Mathematical Physics" department of the Bauman Moscow State Technical University. Full member of the Russian Academy of Engineering Sciences. Author of more than 160 publications in the field of computational mechanics, nonlinear tensor analysis, thermomechanics of composite materials, mathematical simulation in science of materials.
Андрей Николаевич Морозов родился в 1959 г., окончил МВТУ им. Н.Э. Баумана в 1982 г., д-р физ.-мат. наук, профессор, заведующий кафедрой "Физика" МГТУ им. Н.Э. Баумана. Автор более 100 научных работ в области прецизионных измерений и физической кинетики.
A.N. Morozov (b. 1959) graduated from the Bauman Moscow Higher Technical School in 1982. D. Sc. (Phys.-Math.), professor, head of "Physics" department of the Bauman Moscow State Technical University. Author of more than 100 publications in the field of high precision measuring systems and physical kinetics theory.
Александр Павлович Соколов родился в 1983 г., окончил в 2005 г. МГТУ им. Н.Э. Баумана. Канд. физ.-мат. наук, доцент кафедры "Вычислительная математика и математическая физика" МГТУ им. Н.Э. Баумана. Автор 12 научных работ в области применения численных методов в механике композитов.
A.P. Sokolov (b. 1983) graduated from the Bauman Moscow State Technical University in 2005. Ph. D. (Phys.-Math.), assoc. professor of "Computational Mathematics and Mathematical Physics" department of the Bauman Moscow State Technical University. Author of 12 publications in the field of application of numerical methods in mechanics of composites.
Евгений Сергеевич Ничеговский родился в 1985 г., окончил в 2007 г. МГТУ им. Н.Э. Баумана. Аспирант кафедры "Вычислительная математика и математическая физика" МГТУ им. Н.Э. Баумана. Автор ряда научных работ в области численных методов в механике композитов.
E.S.Nichegovskiy (b. 1985) graduated from the Bauman Moscow State Technical University in 2007. Post-graduate of "Computational Mathematics and Mathematical Physics" department of the Bauman Moscow State Technical University. Author of some publications in the field of numerical methods in composite mechanics.