Численное моделирование упругопластического деформирования пространственно-армированных композитов Текст научной статьи по специальности «Физика»

УДК 539.3

Ю. И. Д и м и т р и е н к о, А. П. С о к о л о в, Ю. В. Ю р и н

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ УПРУГО-ПЛАСТИЧЕСКОГО ДЕФОРМИРОВАНИЯ ПРОСТРАНСТВЕННО-АРМИРОВАННЫХ КОМПОЗИТОВ

Предложена модель упругопластического деформирования композиционных материалов с пространственной структурой армирования, основанная на методе асимптотического усреднения. Приведено численное решение локальных задач деформационной теории пластичности на базе метода упругих решений с переменными модулями упругости с применением метода конечных элементов. Предложен метод расчета эффективных упругопластических характеристик композитов. Приведены результаты решения задачи о напряженно-деформированном состоянии композитов с ортогональной 3Б-структурой армирования в соответствии с теорией пластичности Ильюшина.

Е-шаП: dimit.bmstu@gmail.com

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

Композиционные материалы с пространственными структурами армирования в настоящее время широко применяют в технике. Особый интерес представляют композиты на базе металлической матрицы, которые можно эффективно использовать в высокотемпературных конструкциях двигательных установок. Существующие методы расчета эффективных упругопластических характеристик [1—3], как правило, являются приближенными и не обеспечивают необходимой точности расчетов. В этой связи перспективным является метод асимптотического усреднения [4—7], который позволяет найти точные (в математическом смысле) эффективные характеристики с помощью решения так называемой задачи на ячейке периодичности. Модификации метода конечных элементов (МКЭ) для решения этих задач предложены в работах [8—16]. Цель настоящей работы — развитие метода решения задачи упругопластического деформирования на ячейке периодичности, основанного на варианте метода упругих решений с переменными модулями упругости и численной реализации этого метода с помощью конечно-элементных (КЭ) технологий.

Метод асимптотического усреднения для упругопластических задач. Рассмотрим модель ячейки периодичности пространственно-

армированного композита (рис. 1), занимающей в пространстве Я3 область V с внешней поверхностью £ . Композит состоит из матрицы и N - 1 типов волокон, каждый из которых ориентирован по своему направлению. Полагаем, что волокна и матрица являются изотропными упру-гопластическими в соответствии с деформационной теорией пластичности А.А. Ильюшина (теория малых упруго-пластических деформаций) [17]. Обозначим Vа, а =1,..., N, области, занятые

а-м компонентом композита, а-м типом волокна или матрицей (для матрицы будем использовать индекс а = ]); обозначим также £а поверхность областей Vc

Рис. 1. Ячейка периодичности и КЭ-сетка для про-странственно-армированно-го композита с ортогональной ЭБ-схемой

£aN — поверхности

контакта матрицы и волокон, волокна считаем не контактирующими между собой, а £ае — часть поверхности £ композита, занятая

а-м компонентом (причем £а = £а] и £ае —

для волокон и

N-1

£N = ^ £а] и£]е — для матрицы). Тогда в каждой области

а=1

Vа, а-1,..., ], можно рассмотреть следующую задачу теории пла-

стичности:

< у = 0

ау =

sa =

iJ

ua = uN J "a

та (е® )

2 ( + иа) 0)

(о® -а]) = 0 на поверхности £а], иа = и е на поверхности £ае, а^п^ = Б е на поверхности £ае.

п а а а

Здесь ui , а^, е^ — компоненты вектора перемещения, тензоров

напряжения и деформации в а-м компоненте композита; последние два условия в системе (1) отражают условия идеального контакта матрицы и волокон, а Т а (е® ) — нелинейная тензорная функция

определяющих соотношений пластичности матрицы и волокон.

В случае малых упругопластических деформаций, согласно теории пластичности Ильюшина, функции Т* (е1 ) в задаче (1) для волокон и матрицы имеют вид

а _ /~<сс („ос ч,-,«

] _ С]к1 ( тп )5Ы ,

2 Оа(1 - <) - Ка] 55 + Оа( - <) +5а5]к),

са ,а ) _ ijkl\ mn '

где Оа — модуль сдвига; Ка — модуль объемного сжатия; а>"(е^ ) — функция пластичности Ильюшина, зависящая от второго инварианта — интенсивности тензора деформаций 8а=8"е" -—е*2

3

(еа = 8* 5] — первый инвариант тензора деформаций, а С*ш (етп) —

тензор модулей упругости, зависящий от тензора деформаций).

Теперь пусть композиционный материал обладает периодической структурой (см. рис. 1), у которого ячейка периодичности (ЯП) V

состоит из N компонентов а = 1,..., N. Введем малый параметр к = I / Ь «1 как отношение характерного размера ЯП к характерному размеру всего композита, а также глобальные хк и локальные £к координаты. Будем полагать, что матрица является связной областью. Обозначим также = п V — поверхности раздела матрицы и

волокон в ЯП.

В этом случае для такой структуры можно применить метод асимптотического усреднения [8—10, 16], согласно которому решение задачи (1) строят в виде асимптотических по степеням малого параметра к разложений вида

< = и(0)( хк) + к< (1)( хк, £) +...,

в] =88(0)(хк, £) + 8(1)(хк, £) +..., (2)

а* = а] (0)( хк, £) + ка* (1)( хк, £) +...,

причем по аргументу £ эти функции полагаются периодическими. Деформации и напряжения так называемого нулевого уровня имеют следующий вид:

]>=*< + 1 (+). <3»

_ 1

= ■

r«(0) =

ff

2 (u (0 + u <0))

Ja (%?0)), если ? e VaS, a = 1,..., N,

(5)

где •,1 = д»/дх и •/1 = д»/— производные по двум типам координат. При выводе формул (3)—(5) и далее используется правило дифференцирования асимптотических разложений.

Подставляя разложения (4) в систему (1), применяя правило дифференцирования и собирая члены при одинаковых степенях к, получаем так называемую локальную задачу на ячейке периодичности:

C7a(0) = 0 / j = 0,

a(0) =

Ja /_a(0)\ Jj \Skl ),

ra(0) =- , 1 (ua(1) , ua(1) Sj =Sj + 2 \U,/ j + Uj/,

ua(1) = uN (1)

(6)

на поверхности {&"(1) -а^(1) ^ = 0 на поверхности Е

(<(1)>

= 0,

,a(1)

= 0.

Здесь оператор усреднения обозначен как

N N

a=1 V

a=1v.

В задаче (6) условие

u

= 0 представляет собой условие периодичности функций на границе ячейки периодичности, а условие (<) = 0 определяется требованием единственности решения локальной задачи [5]. В силу периодичности функций мга(1) имеет место следующее соотношение:

^ =Ua (0)

Метод упругих решений для локальных задач на ячейке периодичности. Предположим, что волокна и матрица изотропные, ЯП является симметричной при зеркальном отражении относительно

трех координатных плоскостей 2 ^ = = 0}, а также симметричной

относительно поворота на угол ж вокруг каждой оси координат и симметричной при преобразовании центральной симметрии с центром в точке О (подробнее об этих преобразованиях см. [18]). Это ограничение назовем основным допущением о симметрии композита. Согласно этому допущению, вместо решения локальной задачи (6) для всей области V^ ЯП можно перейти к решению задачи в области

V 1/8 ЯП.

Задача (5) является нелинейной, поэтому для ее решения применим итерационный метод, являющийся разновидностью метода упругих решений [17]. Согласно этому методу, определяющие соотношения в системе (6) для компонентов композита линеаризуются следующим образом:

а{т} _ у"<а {т-1} а{т} = С]к1 ек1 ,

Са {т-1} = с* (_а{т-1})

С,]к1 ~ С,]к1 \5тп ),

где с"{т} и е"{т} — значения напряжений сС(0) и деформаций (0)

С а {т-1} ^

]к1 — тензоры модулей упругости компонентов композита на (т-1)-м шаге итерации. Обозначим также иа {т} — значения перемещений и] (1) на т-м шаге итерационного цикла. Тогда на т-м шаге итерации вместо задачи (4) получаем следующую линеаризованную задачу:

гга{т} = 0 / ] = 0,

а{т} _ /~ча {т-1} а{т} ~ С]к1 к1 ,

еа{т} = в + 1 (иа{т} + иа{т}\ е,] =е] + 2 \и, / ] +и] /, ),

и*{т} = и,^т} на поверхности 2 ^,

({т} -а (т})=0.

Согласно предложенному в [16] варианту метода асимптотического усреднения, перемещения первого уровня и напряжения нулевого уровня при каждом значении т можно представить в виде следующих сумм:

3 3 3

..a{rn} _ ST iMm) „aim] V1 ^{m} V1 —a{m}

ui - Z ui (pq), ea } = Z etj(W), au } = Z ^yXpq)'

p,q=l p,q=1 p,q=1

если gk eF^, a-1,..., N, (7)

причем для функций u"^ для каждой комбинации (pq) выделяется линейная часть по локальным координатам:

<S = -spq fa? + Siqgp ) + U?™ ), (8)

где Sp — символ Кронекера; U^m} (g) — некоторые функции,

называемые псевдоперемещениями, для которых при каждом фиксированном наборе индексов (pq) получаем следующую задачу, называемую линеаризованной локальной задачей Lpq на ячейке периодичности:

_a{m} - о

°ч(pq)/j

a{m} _ sia {m-1} a{m} aij(pq) - Cijkl Sij(pq),

_a{m} - 1 Ta{m} л TTa{m} ^ (Q)

bij(pq)~ 2 \Ui(pq)/j j(pq)/i), (Q)

rra{m} TJN {m} „ v

Ui (pq) - Ui (pq) на поверхности 2gaN ,

(({m} {m} )и - 0 lCTij-(pq) °ij(pq))"j~

К системе (Q) добавим условия на координатных плоскостях

25 - {g' - 0} и на торцевых поверхностях ЯП 2' - g -1/2}, ' =

= 1, 2, 3, которые запишем следующим образом: при (p = q)

Uapq) - 1 ZpqSip , S*(pq) - 0, S"(pq) - ° на поверхности % , 1 * j * k* Ц (10)

при (p * q)

Ua - (1 / 4)^Sp, Sapq) - 0, Uapq) - 0 на поверхности 2j, U j - (P, q SKpq) - 0, Sapq) - 0, Uapq) - 0 на поверхности Ъ'к, i * j * k * i;

при (p = q)

UKpq) - 0, Sj(pq) - 0, Sk(pq) - 0 на поверхности 2i, J * j * k * (1 1)

при (p * q)

UKpq) = 0 Sj(pq) = Uk(pq) = 0 На ПОвеРХНОСТИ, 2 j i, j = {p, q},

а _/-ч oor _/-Ч TTa

i (pq) - 0, SJ ( pq)~ 0, Uk( pq)

S<( pq) = 0 Sj( pq) = 0 Uk( pq) = 0 На ПОвеРХНОСТИ 2k , i * j * k *

Вывод этих граничных условий для линеаризованных задач (7) такой же, как и для линейных задач на ЯП [8]. В настоящей работе введены обозначения для векторов сил (рд):

3

(рд)- ¿и ай (рд )пг-

г=1

Решение задачи (7)—(11) будем искать в области У^, представляющей 1/8 часть ЯП: У^= (^ > 0), где 2^ — поверхность контакта компонентов композита внутри области V^ : I^ =

Появление условий (10), (11) связано с требованием периодичности всех функций (1) на границе ЯП, а также наличием симметрии ЯП относительно трех координатных плоскостей для рассматриваемого типа композита.

Вариационная формулировка локальной задачи Ьрч. Для произвольного конечного объема V е V вариационная формулировка задачи Ьт (9) при фиксированных значенияхр и д имеет вид

\5е*о<Ш = \5и тЖ. (12)

V I

Здесь координатные столбцы псевдоперемещений и обозначены: напряжения с, деформация е и поверхностные силы 8 :

-|Т

и =.........................

тта{ш} тта{ш} тja{m} Ul(pq), U2(pq), U3(pq)

О

[oa{m} oa{m} oa{m} oa{m} / Г2 oa{m} / /2 0«{m} / /2T °11( pq )>°22( pq )> °33( pq)>°13( pq) ' V Z' °23( pq)/VZ>°12( pq) ' ^ Z J '

_ Г a{m} a{m} a{m} «{m} ,f2 ca{m} , Г2 ca{m} , Г2T _|_fc11(pq)' 22(pq)' 33(pq)> fc13(pq)'VZ> ¿23(pq)/VZ'fc12(pq) ' V ZJ '

S =

oCf oCf oCf

S1(pq), S2(pq), S3(pq)

Линеаризованные определяющие соотношения, следующие из второй группы уравнений (7), с использованием координатных столб-

цов записывают следующим образом: а - Се, где С — матрица упругости размером 6 х 6, составленная из компонентов тензора стандартным образом [18]. Соотношения Коши (третья группа уравнений в (9)) в матричном виде записывают следующим образом: е = DU, где D — матрица линейных дифференциальных опера-

торов дифференцирования (д1 = 5 / ):

Г 5! 0 0 1

0 5 2 0

0 0 5з

Б - 5 3/л/8 0 51 / л/8

0 53/Т8 5 2/>/8

ч5 2/^8 5,/^ 0 ,

С учетом указанных матричных соотношений вариационное уравнение (12) можно представить в виде

|CDUdV - \5итSdЕ.

V Е

Метод конечного элемента (МКЭ) для задач Ьрч, Аппроксимируя пседоперемещения U в КЭ линейными функциями: и - Фq, где q — координатный столбец псевдоперемещений в узлах КЭ; Ф(%к) — матрица функции формы, зависящая от типа КЭ, получаем разрешающую систему линейных алгебраических уравнений:

Кд - /,

где

К -1ВтCBdV, / -1 ФтSdЕ —

V Е

локальная матрица жесткости и столбец нагрузок; В - БФ.

Глобальную матрицу жесткости задачи составляют из локальной матрицы жесткости стандартным образом [19], после ее формирования к системе линейных алгебраических уравнений (СЛАУ) применяют граничные условия (10) и (11). Граничные условия идеального контакта (последняя группа соотношений (7)) не требуют специального учета, так как при данном варианте МКЭ они удовлетворяются автоматически.

Решая СЛАУ, находим перемещения q в узлах, по которым вычисляем: псевдоперемещения и - Фq, деформации е - Bq и напря-

жения а - CBq в КЭ. Для решения СЛАУ применяли метод сопряженных градиентов.

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

цикл решения соответствующей задачи до достижения условия сходимости решения:

ß,r=\

ZCT(«){m> -1} у

ßr( pq ) ßr ( pq ) ¿—i

f з Y1

\ßr=1

<J(a){m _1} ßr pq)

<8,

где 5-10 3,..., 10 5. Число итераций т -М, обеспечивающее выполнение данного условия, различается для различных задач Ьт и для различных значений ерд, но не превышает 10-15. Напряжения

при максимальном значении номера итерации т - М обозначали как а^рАрМ). Для обеспечения устойчивости решения задачи Ь11, Ь22 и Ь33

решали совместно как задачу с входными данными е11 е22 и ё"33.

Также совместно решали задачи Ь13 и Ь31, Ь12 и Ь21, Ь23 и Ь32. Далее проводили еще один цикл решения задач Ьрд, используя значения входных данных ё„п.

РЧ

Расчет эффективных упругопластических характеристик композиционного материала. После решения серии задач Ьт (7), (10), (12) указанным методом для всех значений pq проинтегрируем

напряжения а^М} по областям, занятым волокнами и матрицей:

а, }) -IЩМ }>- <>3)

р Л-1

в результате получим усредненные напряжения. В этом случае эффективные упругопластические соотношения, связывающие усредненные напряжения и деформации а, и ерд, можно записать в

следующем символическом операторном виде:

a<j =

Ъ [tpq) • (14)

Поскольку выше указан алгоритм вычисления средних напряжений а, (формулы (13), в которых арА){рМ)} вычислены методом, описан-

ным в предыдущем пункте) по заданным значениям усредненных деформаций ерд, фактически указан алгоритм нахождения значений

символического оператора (14).

Оператор (14) можно конкретизировать, если имеется информация о типе геометрической симметрии ЯП композита и типе анизотропии волокон и матрицы. В силу основного допущения симметрии решение всех задач Ьрд будет обладать рассмотренными типами симметрии (т. е. не изменится во всей ЯП при указанных выше преобразованиях). Будем полагать, что диаметры волокон, ориентированных по разным координатным направлениям, одинаковы, тогда ЯП обладает симметрией относительно поворотов на угол п /2 . Это означает, что и оператор (14) должен обладать данным типом симметрии. Но перечисленные выше преобразования вместе с тождественным преобразованием образуют группу кубической симметрии [18], следовательно, по терминологии, введенной в [18], оператор (14) будет тензорной функцией, индифферентной относительно группы квазиизотропии (кубической симметрии). Кроме того, тензорные функции (2) матрицы и волокон являются квазилинейными, т. е. зависят только от первого и второго инвариантов тензора деформации е" е13, тогда и решение локальных задач (9), (10) также зависит только от линейных и квадратичных инвариантов тензора усредненных деформаций Щ. В группе кубической симметрии таких инвариантов только три [18], которые для тензора деформаций имеют вид

А = е11 + е22 + е33, 1е = (е11 _ е22) + (е11 ~е33) + (е33 _ е22) , А3 =

- ¿>23 &13 12'

Тогда функцию (14) можно представить в тензорном базисе группы квазиизотропии [18]:

3

= ( ерд) = ^5ч + Ё (2Р " Р3 )8>78]7ёгг + %Щ,

Г=1

где

РГ=РГ(-1е, 12 , !еъ ), У = 1,. ., 3, (15)

скалярные функции от трех инвариантов тензора усредненных деформаций.

Функции рг фактически представляют собой искомые эффективные упругопластические характеристики композита. Зная значе-

ния << и зрд, эти функции <у можно вычислить следующим образом. Запишем формулы (15) явным образом:

<11 = <<11\ + <22 = <<11\ + 2<2з22^ <13 = <3^13,

<33 = <11\ + 2<2833, <12 = <3 312, <23 = <3 З23,

а затем образуем из этих соотношений три инварианта тензора средних напряжений 1<, 1<, 1<, аналогичных инвариантам 1Еу, у = 1,..., 3,

в результате получим систему трех уравнений относительно трех функций <ру :

93 = 41< / 1з, 3<1 + 2< = 1< /1', 2< = 4 1</А,

решая которую находим искомые эффективные упругопластические функции композита (15). Численные значения этих функций вычисляют на основе решения локальных задач.

В области упругости матрицы и волокон композит также проявляет только упругие свойства, в этом случае функции <ру являются

константами:

<3 = 20, <1 = Я, <2 = Д

Отметим, что для квазиизотропного материала О фд и имеется три независимые упругие константы.

Результаты численного моделирования упругопластического деформирования композита. В вычислительном эксперименте рассматривался пространственно-армированный углеалюминиевый композит с алюминиевой матрицей и углеродными волокнами, которые для простоты предполагались изотропными и упругими. Константы упругости волокна были выбраны следующими: модуль Юнга Е = 250 ГПа, коэффициент Пуассона V = 0,25. Модель материала матрицы считали упругопластической, а функция пластичности Ильюшина имела вид

Г0 8м < 8т

< (зы) =! " '

" " "Ь-(зт+з-зту/но)/зм, зм >зт,

где s, Н0 — константы пластичности; зЦт — предельная деформация текучести. Характеристики матрицы были следующими: Е = 70 ГПа, V = 0,35, предел текучести матрицы зт = 0,4 -10-3, ^ = 0,4, Н0 = = 30 МПа.

Результаты моделирования приведены на рис. 2—5. На рис. 2—4 представлены поля микронапряжений ст"(11) в матрице композита,

полученные в результате решения задачи Ь11. Максимальные значения микронапряжений ст"1(11), ОщП), ОщП) достигаются на границе

раздела матрицы и волокон. Наибольшие значения имеют нормальные напряжения ОщП) (0,81 ГПа) и сдвиговые напряжения сгщП)

(0,11 ГПа) в матрице. На рис. 4 показано распределение остаточной пластической деформации еЦ(11) в матрице композита с 3Б-схемой

армирования композита после предварительного нагружения в пластической области. Максимум остаточной деформации достигается на границе раздела матрицы и волокон, а также в срединной части между волокнами.

Рис. 2. Распределение напряжения Сщп) и °12(и) в матрице ЯП композита с ортогональной 3.0-схемой армирования

Рис. 4. Распределение остаточной пластической деформации £щИ) в матрице композита ЯП

На рис. 5 показаны численно полученные диаграммы деформирования композита, определяющие зависимость <11 = Т1Х ^ при двух

различных концентрациях волокон, а также диаграммы деформирования алюминиевой матрицы и армирующих углеродных волокон. Композит с меньшей объемной долей волокон (46 %) проявляет большие пластические свойства, его диаграмма деформирования близка к диаграмме деформирования матрицы.

<тп Г Па

УВ / УВ/А1-54Я'*

УВ/Ад-46

AI

О 0,001 0,002 0,003 0,004 0,005 0,006 0,007 0,003

Рис. 5. Диаграммы деформирования углеродного волокна (УВ), алюминиевой матрицы (А1) и двух моделей ортогонально-армированных композитов с объемной долей волокон 46 (УВ/Аg-46) и 54 % (УВ/А1-54)

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

СПИСОК ЛИТЕРАТУРЫ

1. А д а м с Д. Ф. Упругопластическое поведение композитов / Композиционные материалы // Механика композиционных материалов. Т. 2: пер. с англ. М.: Мир. 1978.

2. В и л ь д е м а н В. Э., С о к о л к и н Ю. В., Т а ш к и н о в А. А. Механика неупругого деформирования и разрушения композиционных материалов. - М.: Наука. 1997.

3. К р и в о ш е и н а М. Н., К о б е н к о С. В., Т у ч Е. В. Усреднение свойств композиционных материалов при численном моделировании их разрушения // Физическая мезомеханика. - Т. 13. - № 2. - 2010. - С. 55-60.

4. П о б е д р я Б. Е. Механика композиционных материалов. - М.: Изд-во МГУ, 1984.

5. Б а х в а л о в Н. С., П а н а с е н к о Г. П. Осреднение процессов в периодических средах. М.: Наука, 1984.

6. С а н ч е с - П а л е н с и я Э. Неоднородные среды и теория колебаний. М.: Мир, 1984.

7. Д и м и т р и е н к о Ю. И. Механика композиционных материалов при высоких температурах. М.: Машиностроение, 1997.

8. Д и м и т р и е н к о Ю. И., К а ш к а р о в А. В. Конечно-элементный метод для вычисления эффективных характеристик пространственно-армированных композитов // Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки. 2002. № 2.

9. Д и м и т р и е н к о Ю. И. Метод многоуровневой гомогенизации иерархических периодических структур // Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки. 2002. № 2.

10. Д и м и т р и е н к о Ю. И., С о к о л о в А. П. Многомасштабное моделирование упругих композиционных материалов // Математическое моделирование. 2012. Т. 24. № 5. С. 3-20.

11. Д и м и т р и е н к о Ю. И., Д у б р о в и н а А. Ю., С о к о л о в А. П. Конечно-элементное моделирование усталостных характеристик композиционных материалов // Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки. Спец. Выпуск. Математическое моделирование. 2011. С. 34-50.

12. Д и м и т р и е н к о Ю. И., С о к о л о в А. П. Современный численный анализ механических свойств композиционных материалов // Изв. РАН. Физическая серия. Т. 75. 2011. № 11. С. 1551-1556.

13. Д и м и т р и е н к о Ю. И., М о р о з о в А. Н., С о к о л о в А. П., Н и ч е -г о в с к и й Е. С. Моделирование эффективных пьезоэлектроупругих свойств композиционных материалов // Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки. 2010. № 3. С. 86-97.

14. Д и м и т р и е н к о Ю. И., С о к о л о в А. П. Об упругих свойствах композиционных материалов // Математическое моделирование. 2009. Т. 21. № 4. С. 96-110.

15. Д и м и т р и е н к о Ю. И., С о к о л о в А. П. Автоматизация прогнозирования свойств композиционных материалов на основе метода асимптотического усреднения // Информационные технологии. 2008. № 8. С. 31-38.

16. Д и м и т р и е н к о Ю. И., К а ш к а р о в А. И., М а к а ш о в А. А. Конечно-элементный расчет эффективных упругопластических характеристик композитов на основе метода асимптотического осреднения // Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки. № 1. 2007. С. 102-116.

17. И л ь ю ш и н А. А. Пластичность. - М.: Физматлит. 2004.

18. Д и м и т р и е н к о Ю. И. Тензорное исчисление. М.: Высш. шк., 2001.

19. С е г е р л и н д Л. Применение метода конечного элемента: пер. с англ. М.: Мир. 1979.

Статья поступила в редакцию 03.07.2012.