CC BY
31
УДК 519.6
М. П. Г а л а н и н, Е. Б. Савенков, Ю. М. Т е м и с, И. А. Щеглов, Д. А. Яковлев
ПРИМЕНЕНИЕ МЕТОДА КОНЕЧНЫХ СУПЕРЭЛЕМЕНТОВ ДЛЯ РАСЧЕТА ХАРАКТЕРИСТИК ДИСПЕРСНО-АРМИРОВАННЫХ КОМПОЗИЦИОННЫХ МАТЕРИАЛОВ
Представлен программный комплекс, предназначенный для решения задачи теории упругости композиционных материалов в трехмерном случае. В качестве расчетного метода использован метод конечных суперэлементов Р.П. Федоренко. Приведены результаты численного решения для различных упругих и геометрических параметров материалов, составляющих композиционный материал.
Метод конечных суперэлементов (МКСЭ) предложен в работах [13] для задач, решения которых имеют локальные особенности геометрического или физического типа. При этом особенности геометрического типа порождаются резкими изменениями границы области и наличием в ней пустот, а физического типа — присутствием подобластей с резко отличными параметрами материала. Предполагается, что размер особенностей существенно меньше как размера исходной расчетной области, так и приемлемого шага расчетной сетки. Ранее метод применялся для решения ряда сложных задач диффузии, теории упругости, кинетики ядерных реакторов и ряда других.
В работах [4-9] предложен вариант теоретического исследования МКСЭ, который позволил формально строить аппроксимации метода для достаточно широкого класса эллиптических задач, в частности для линейных задач теории упругости.
В настоящее время существует несколько различных способов описания и анализа свойств композитов. Такими подходами являются, например, метод асимптотического осреднения Бахвалова-Победри [10], различные микромеханические [11-15] и смесевые модели [16-19], широко применяемые в инженерной практике. Эти модели позволяют определять приведенные упругие свойства композита в зависимости от свойств его структурных элементов. Краткий обзор этих моделей приведен в [20]. Там же подробно описаны процедура и некоторые результаты применения МКСЭ для определения приведенных упругих свойств композитов.
^Работа выполнена при частичной финансовой поддержке Российского фонда фундаментальных исследований (проект РФФИ № 06-01-00421).
Актуальность МКСЭ для решения описанной задачи связана с тем, что он позволяет проводить прямое численное моделирование напряженно-деформированного состояния композиционного материала с непосредственным учетом большого числа включений. Решение такой задачи обычными конечно-разностными или конечно-элементными методами представляет большую трудность, связанную с большой размерностью конечномерной задачи. Так, даже в области с одним включением необходимо использовать очень подробные сетки, позволяющие хорошо разрешить решение в окрестности контакта матрицы и включения. Число узлов в такой сетке может доходить до 105 и более. Понятно, что решить задачу в области с несколькими включениями на сетке с таким шагом очень сложно.
Помимо этого в работе [20], продолжением которой является данная работа, приведен алгоритм оценки приведенных характеристик изучаемого композиционного материала. Он основан на сравнении численного решения задачи в области с неоднородными упругими свойствами, полученного с помощью МКСЭ, с аналитическим решением задачи в однородной изотропной области с неизвестными приведенными упругими свойствами, подлежащими определению. Этот алгоритм использовался и в данной работе.
Одной из таких задач, для которых можно ожидать эффективного использования МКСЭ, является задача расчета напряженно-деформированного состояния композиционных материалов. Отметим, что использование МКСЭ для подобных задач проводилось и раньше [1-3]. Схема МКСЭ для решения указанной задачи приведена в [20]. Там же представлены результаты первых расчетов, указывающих на эффективность метода, проведено сравнение численных результатов с экспериментальными и известными теоретическими оценками. Получено хорошее совпадение расчетных значений как с данными эксперимента, так и с теоретическими оценками. В частности, использовались двусторонние оценки Фойгта-Рейсснера и Хашина-Штрикмана, основанные на двойственных вариационных формулировках задач теории упругости, которые содержат функционалы с совпадающими стационарными значениями, но достигающими в своих стационарных точках альтернативных экстремумов.
В качестве объекта исследования выбраны дисперсно-армированные композиционные материалы, для которых рассмотрена задача расчета напряженно-деформированного состояния композитного тела с помощью МКСЭ и, как следствие, задача расчета упругих характеристик композитов. Для решения указанной задачи разработан программно-вычислительный комплекс, построенный на модульном
принципе, позволяющий в автоматизированном режиме подготавливать и запускать расчет, а также анализировать его результаты. Комплекс также позволяет проводить расчеты для композиционных материалов с включениями в виде коротких волокон (эллипсоидальной формы).
Постановка задачи и аппроксимации МКСЭ. Будем считать, что некоторый материал — матрица — армирован включениями, форма которых близка к сферической. Поведение материала считаем упругим.
В общем случае представляет интерес детальная информация о напряженно-деформированном состоянии композита и его приведенных (усредненных) свойствах в зависимости от упругих и геометрических параметров матрицы и включений, например: объемная доля включений в материале, их геометрия и ориентация и т.д.
Рассмотрим постановку задачи, подробное описание которой приведено в работе [21]. Требуется определить поле перемещений и = {иъ}, г = 1,3 сплошной среды в ограниченной трехмерной области П, на границе дП которой заданы кинематические и динамические граничные условия. Будем считать, что в области П задана ортогональная декартова система координат Ох\х2х3. В дальнейшем будем использовать правило суммирования по повторяющимся индексам. Производную некоторой величины А по переменной х ъ будем обозначать через А,ъ, т.е. дА/дхг = А,,,.
Система соотношений, описывающих напряженно-деформированное состояние упругого тела, состоит из уравнения равновесия
' + /г = 0
материальных соотношений, выражающих обобщенный закон Гука,
= Е ( + » А ^
аг} , I £г} + п £ккАг} I ,
' ^ + 1 V ' 1 - V
и кинематических соотношений
£,} = 1/2 (и' + Щ,) .
Здесь а — тензор упругих напряжений; £ — тензор деформаций; а,', £,' — их компоненты в выбранной системе координат; / — объемная плотность внешних сил; Е — модуль Юнга; ^ — коэффициент Пуассона, являющиеся заданными функциями точки пространства; Аг} — компоненты единичного тензора.
Будем считать, что граница дП области П состоит из двух частей, т.е. дП = Гд и Гп. На части Гд заданы кинематические граничные условия
и\Тд = ид, и \Гд = ид,г,
а на части Гп границы заданы естественные граничные условия
^|Г„ = }д, 0г, V |г„ = ,
где
0 = 0V = {аг, V,};
V = {V,} — внешняя нормаль к границе области.
МКСЭ может рассматриваться как обычный метод Галерки-на-Бубнова с использованием специальных базисных функций. Поэтому конечномерная задача может быть получена следующим образом. Необходимо: 1) задать разбиение расчетной области на суперэлементе (СЭ); 2) задать на границах СЭ граничные функции, аппроксимирующие следы решения на границах СЭ; 3) каждую такую граничную функцию продолжить во внутренность СЭ решением исходного уравнения с указанной функцией как граничным условием первого рода, для чего внутри каждого СЭ независимо ввести расчетную сетку (см. [22, 23]) и использовать обычный метод конечных элементов (см. [24, 25]); 4) систему линейных алгебраических уравнений для определения узловых значений СЭ получить из условий, формально совпадающих с условиями метода Галеркина-Бубнова с рассчитанными функциями в качестве базисных.
Отметим, что, так как для СЭ базисные функции рассчитываются численно, то, вообще говоря, суперэлементы могут иметь произвольную форму.
Естественным требованием при разбиении области на СЭ в случае конкретной рассматриваемой задачи является расположение каждого включения композита целиком внутри одного СЭ. Это означает, что на границах СЭ упругие свойства среды не меняются, и, значит, на границе СЭ решение является относительно гладким. Это позволяет использовать для аппроксимации перемещений на границах СЭ простейшие билинейные базисные граничные функции. Для расчета решения внутри СЭ используется достаточно мелкая сетка для того, чтобы описать особенности решения, в частности, на контактной границе между матрицей и включением.
В общем случае все СЭ могут быть разными, каждый СЭ и соответствующее ему включение могут иметь разную форму и различные значения упругих параметров. В представленных расчетах все СЭ имеют форму куба, включения могут иметь разную форму, упругие коэффициенты матрицы и включения не меняются от СЭ к СЭ.
Программный комплекс. Для решения поставленной задачи разработан программно-вычислительный комплекс.
Комплекс имеет модульную структуру, обмен данными между модулями организован с помощью файлов ввода-вывода. В состав
комплекса входят: 1) препроцессор, основанный на пользовательском графическом интерфейсе (GUI); 2) основная расчетная программа; 3) постпроцессор.
Модульное исполнение комплекса допускает независимую замену и доработку перечисленных выше программ при условии сохранения форматов обмена данными.
Препроцессор поддерживает следующие функции.
1. Организация расчета в виде проекта: для расчета автоматически создается папка, где содержатся файл проекта с параметрами расчета, файлы с конечно-элементными (КЭ) сетками СЭ, файлы данных расчетных программ.
2. Сохранение и загрузка всех параметров расчета в файле проекта.
3. Возможность генерации и сохранения файлов для расчетной программы, содержащих сведения о конечно-элементных сетках СЭ, СЭ сетке и граничных условиях в зависимости от введенных параметров, а также изменение параметров расчета и перегенерация соответствующих файлов.
4. Графическое представление и статистический анализ построенных КЭ и СЭ сеток.
5. Автоматизированная генерация параметров СЭ и распределение СЭ по расчетной области в соответствии с заданным вероятностным законом распределения отдельных параметров.
6. Автоматизированный подбор конфигурации расчетной области по параметрам СЭ и заданной объемной плотности включений.
7. Генерация СЭ с эллипсоидальными включениями. При этом углы наклона осей эллипсоида и величины полуосей могут варьироваться в автоматизированном или ручном режиме. После сохранения изменений происходит автоматическая перегенерация всех связанных объектов таких, как размеры области СЭ и КЭ сетки СЭ.
На первом этапе при генерации КЭ сеток СЭ пользователь может выбрать параметры разбиения: число разбиений ребра СЭ nR = nx = = ny = nz и шаг h разбиения КЭ сетки СЭ [22-23]. Одинаковый размер СЭ и единые параметры дискретизации позволяют состыковывать СЭ без пересчета граничных узлов. При формировании расчетной области существует возможность добавления в нее СЭ без включений.
Примерный вид области размера 2 х 2 х 2, сгенерированной в препроцессоре и содержащей СЭ четырех типов, представлен на рис. 1.
Отметим, что схема расчета базисных функций СЭ идеально пригодна для параллельных вычислений, так как в большинстве случаев это самая ресурсоемкая операция, а расчеты отдельных типов СЭ не зависят друг от друга. Поэтому поддержка распараллеливания этой процедуры позволит многократно сократить общее время расчета.
Рис. 1. Расчетная область композитного тела (копия экрана препроцессора)
Результаты численного моделирования. Результатом работы основной расчетной программы является файл, содержащий решения в узлах СЭ сетки. На его основе проводится расчет приведенных параметров композита. Минимальные, максимальные и средние (по узлам суперэлементов; на графиках — min, max, avg соответственно) значения модуля Юнга E = E3 и коэффициентов Пуассона ц1, ц2 сохраняются в файле.
Это набор упругих параметров, которые можно определить для композиционного материала, сравнивая численное решение с точным решением задачи об одноосном растяжении однородного изотропного упругого стержня в направлении х3.
Пусть на торцах заданы нормальные напряжения
Yv о = pv,
где v = (0, 0, 1) или v = (0, 0, -1) — внешняя нормаль к границе области. На боковой поверхности области задано условие
Yv о = 0.
Если известны модуль Юнга E и коэффициент Пуассона ц среды, перемещение точки среды, совпадающей с началом координат (примем его равным нулю), то решение этой задачи известно и в терминах перемещений точек среды имеет вид:
и1 =
E
Х1, u2 =
цр p
—х2, u3 = — Х3.
E 3 E
(1)
Таким образом, если известны перемещения какой-либо точки стержня и известна нагрузка р, то с помощью соотношений (1) можно определить упругие параметры:
Хз Ui Хз U2 Хз
E = —p, Ц =---=---.
U3 Х1 U3 Х2 u3
(2)
Подставив в выражения (2) численное решение этой задачи, полученное с помощью МКСЭ для композиционного материала, вычислим значения упругих параметров в узлах суперэлементной сетки и, следовательно, min, max, avg значения для всей области (подробнее см. работу [20]).
Отметим, что если при разбиении области используется несколько различных типов СЭ, отличающиеся размером включений, то для небольшого числа СЭ в области (порядка 103 ... 104) не выполняется соотношение ц = ц1 = ц2, где
Щ Хз 42 Хз
Ц1 =---, Ц2 =---. (3)
Х1 43 Х2 43
Поэтому для такой серии расчетов приводятся результаты для Е = Ез, Цх = ць Цу = ц (рис. 10-13).
Полученные в результате расчетов приведенные упругие параметры хорошо согласуются с двусторонними оценками Хаши-на-Штрикмана.
Варианты расчетов. Проведен ряд расчетов, результаты которых позволяют проанализировать поведение композитного тела в зависимости от составляющих его материалов, размеров включений и распределения включений по области, а также выяснить зависимость решения от параметров дискретизации СЭ и размера расчетной области П. Расчеты группируются следующим образом.
1. Анализ зависимости значений приведенных упругих параметров композиционного материала: а) от значений упругих параметров включений СЭ; б) от размеров включения при неизменных значениях упругих параметров включений СЭ; в) от значений параметров разбиения КЭ сетки СЭ при неизменных значениях размеров включения и упругих параметров включений СЭ; г) от размера расчетной области (N0 = N х N х N) при неизменных параметрах материалов и заданном СЭ (размеры, сетка).
2. Анализ изменения результатов расчетов в зависимости от варианта распределения нескольких типов СЭ по области при неизменных значениях размеров включения, упругих параметров и расчетных сетках СЭ.
Область П является прямоугольным параллелепипедом, содержащим N0 = N х N х N СЭ. При генерации КЭ сеток СЭ задаются ориентировочный шаг сетки к и число NR разбиений ребра ячейки СЭ. Чем больше NR, тем больше узлов у КЭ сетки СЭ. Например, для отдельного СЭ зависимость числа узлов и ребер тетраэдральной сетки от NR представлена в табл. 1.
Во всех расчетах используются ячейки композита со сферическим включением различного радиуса, центр которого совпадает с центром СЭ. Все СЭ имеют кубическую форму. Для удобства модуль Юнга включения Ef нормирован на модуль Юнга матрицы, т.е. во всех расчетах Ет = 1,0.
При проведении расчетов варьировались следующие параметры задачи: размер расчетной области; параметры материалов включения
Таблица 1
Базовые параметры КЭ сетки СЭ
Число разбиений ребра СЭ Число узлов Число ребер
4 1159 7194
6 3597 23096
8 8171 53470
10 15553 103020
и матрицы — Ef ,Ет, цf, цт; параметры разбиения КЭ сетки СЭ — к, NR; распределение фиксированного набора СЭ по области; число различных СЭ, составляющих область.
В табл. 2 приведены базовые параметры для расчетов, в которых
изменялись один или несколько из представленных параметров.
Таблица 2
Базовые параметры расчетов
Параметр Обозначение Значение
Модуль Юнга матрицы E 1,0
Модуль Юнга включения Ef 6,5
Коэффициент Пуассона матрицы Mm 0,25
Коэффициент Пуассона включения Mf 0,33
Длина ребра СЭ i 1,0
Радиус включения R 0,3
Число разбиений ребра СЭ Nr 10
Число типов СЭ в области Nse 1
Длина ребра расчетной области в СЭ N (Nn) 4 (64)
Рассмотри результаты различных расчетов, выполненных в соответствии с приведенной выше группировкой.
Анализ результатов расчетов. Разные значения модуля Юнга материала включения. Как указывалось, модуль Юнга включения Ef нормирован на модуль Юнга матрицы Ет = 1,0. Использован один тип СЭ с заданными параметрами КЭ разбиения (шаг и число разбиений ребра СЭ), коэффициенты Пуассона матрицы цт и включения цf остаются неизменными: цт = цf = 0,33. Значения модуля Юнга включения Ef варьировались в диапазоне [1,0; 15,0].
На рис. 2-11 представлены графики зависимости минимальных, максимальных и средних значений приведенных параметров.
Как следует из рис. 2, армирование включениями с большим Ef позволяет увеличить общий средний модуль Юнга композитного тела на 30%. Данная зависимость является нелинейной. Например, чтобы получить 15 % прирост, модуль Юнга включения должен быть в 4 раза
Рис. 2. Зависимость модуля Юнга композита от модуля Юнга включения:
♦ — тот; □ — шах; А — аьд
больше, чем модуль Юнга материала матрицы (Ef = 4,0), а для 28 % прироста Ef должен быть равным 15,0.
На рис.3 показана зависимость коэффициента Пуассона ц от Ef. Коэффициент Пуассона монотонно убывает с ростом параметра, т.е. чем больше модуль Юнга армирующих волокон, тем меньше изменение размера области в поперечном направлении, вдоль Ох\ и Ох2.
Зависимость упругих параметров композитов от размеров включений. В каждом расчете используется один тип СЭ со сферическим включением, радиус которого изменяется в диапазоне [0,1; 0,45], параметры материалов и разбиения сетки остаются неизменными: Ef = 6,5, цт = 0,33, ц^ = 0,25. Предельным значением радиуса включения является Я = 0,5.
Графики полученных зависимостей представлены на рис. 4-5. Изменение радиуса включения сильно влияет на изменение приведенных параметров композита для Я > 0,3. Например, если при
/
Рис. 3. Зависимость коэффициента Пуассона композита от модуля Юнга включения:
♦ — min; □ — max; А — avg
Е .2,1 -1,7 -1,3 -
0,9 0 0,1 0,2 0,3 R
Рис. 4. Зависимость модуля Юнга композита от радиуса включения R:
♦ — min; □ — max; А — avg
/ ~
0,3 -
o,zs -
0,2 -1-1-1-1-
О 0,1 0,1 0,3 0,4- R
Рис.5. Зависимость коэффициента Пуассона композита от радиуса включения R:
♦ — min; □ — max; А — avg
R = 0,1 за счет малой объемной плотности включений pf ~ 0,004 влияние армирующего материала на композит ничтожно мало, то при R = 0,45 модуль Юнга Е увеличивается более чем в два раза (pf ~ 0,57), а коэффициент Пуассона ц приближается к значению коэффициента Пуассона включения цf = 0,25 (см. рис. 5).
Разные параметры дискретизации суперэлемента. В данной группе расчетов изменяется число разбиений NR = 4 ... 14 ребра области, ограничивающей СЭ. От этого параметра зависит число узлов и тетраэдров сетки СЭ.
Для определения относительной величины отрезка [min; max] введем оценку "разброс параметра P":
) = max Pt - min Рг max Pi '
где max Pi, min Pi — максимальное (не равное нулю) и минимальное значения вычисляемого упругого параметра P Е {Е,ц1, ц2} для расчета с номером i.
С увеличением числа узлов сетки происходит уменьшение разброса упругих параметров 1(Е), 1(ц). Графики полученных зависимостей упругих параметров материала представлены на рис. 6-7. Разбросы 1(Е), 1(ц) с ростом параметра NR снизились с 0,036 до 0,002 и с 0,024 до 0,001 соответственно.
Рис. 6. Зависимость модуля Юнга композита от числа разбиений ребра СЭ E(Nr):
о — min; □ — max; А — avg
Рис. 7. Зависимость коэффициента Пуассона композита от числа разбиений ребра СЭ:
о — min; □ — max; А — avg
Разный размер расчетной области. Пусть расчетная область состоит из N х N х N СЭ, где N изменяется в диапазоне 2... 30, а общее число СЭ в области изменяется в диапазоне 8... 27000. Радиус включения R = 0,4. Как и в других расчетах, параметры материалов матрицы и включений имеют следующие значения: Em = 1,0, Ef = 6,5, ¡лm = 0,33, ¡f = 0,25; число разбиений ребра ячейки СЭ NR = 10. Чтобы избежать влияния распределения СЭ, все расчеты проводились для одного фиксированного типа СЭ.
На рис. 8-9 приведены графики зависимостей приведенных упругих параметров композитного тела от параметра N. Как видно на рисунках, с увеличением N происходит сужение интервала [min; max] (это следует из увеличения размеров области и числа СЭ, по которым происходит усреднение) и плавное незначительное изменение средних значений — увеличение E и уменьшение
Рис. 8. Зависимость модуля Юнга композита от параметра N:
о — min; □ — max; А — avg
F ~
0,29 -0,28 -
0,27
О 5 10 15 20 25 дО N
Рис. 9. Зависимость коэффициента Пуассона композита от параметра N:
о — min; □ — max; А — avg
Для анализа зависимости решения от размера области введем оценку:
где а — вычисленное значение, а г,] принимают те же значения, что и N.
Для графиков avg, изображенных на рис. 8-9, оценка е составила соответственно е(Е) = 3,5 • 10-5 и е(д) = 1,1 • 10-4, т.е. в целом поведение графиков средних значений не зависит от значения N, а для оценки упругих параметров допустимо использовать значения N ^ 4.
Несколько вариантов распределения различных СЭ по области. Расчеты проведены для области фиксированного размера N = 4, П = 64 и разных распределений заданных типов СЭ (^е = 3) по области.
Разброс СЭ при этом по области имеет псевдослучайный характер. Проведен расчет для восьми случаев распределения, при котором сохранялась объемная плотность включений, рассчитанная для всей области. Как следует из графиков, представленных на рис. 10, а, 11, а, распределение СЭ по области незначительно влияет на средние значения упругих параметров композита.
Аналогичные выводы можно сделать и для коэффициентов Пуассона (см. рис. 10, б, 11, б).
Сравнивая графики на рис. 10-11 с соответствующими графиками на рис. 2-9, где использовался единственный тип СЭ, можно с уверенностью сделать вывод о том, что большой разброс размеров СЭ, как в случае данной группы расчетов, приводит к большому разбросу минимальных и максимальных приведенных параметров СЭ относительно среднего значения.
Можно ожидать, что для больших областей и большего числа типов СЭ зависимость результатов расчетов от распределения СЭ по области будет еще меньше; приведем результаты аналогичных расчетов для варианта N = 10 (см. рис. 11). Графики зависимостей коэффициентов Пуассона от варианта распределения для случая N =10 изображены на рис. 11, б.
Рис.10. Зависимость модуля Юнга Е Рис.11. Зависимость модуля Юнга Е (я) и коэффициентов Пуассона м®, Му (я) и коэффициентов Пуассона М®, Му (б) от варианта распределения при (б) от варианта распределения при N = 4: N =10:
х — min; о — max; А — avg
х — min; о — max; А — avg
Для графиков ауд(Е), ауд(цх), ауд(^у), изображенных на рис. 10, получаем
е(Ен=4) = 0,046, е(^=4) = 0,07, е(^=4) = 0,07, а для соответствующих графиков в случае N = 10 имеем (см. рис. 11) е(Елт=1о)= 0,011, е(^=10 ) = 0,012, е(^=10)= 0,011.
Таким образом, для области, состоящей из 1000 СЭ, влияние различных вариантов распределений СЭ по области на среднее значение упругого параметра минимально.
Заключение. Разработанный алгоритм МКСЭ применен для решения задачи теории упругости композиционных материалов в трехмерном случае. Для решения задачи разработан программно-вычислительный комплекс, позволяющий в автоматизированном режиме проводить расчеты упругих параметров композиционных материалов.
Программы, входящие в его состав, могут быть легко модифицированы или дополнены новыми функциональными возможностями для поддержки дальнейших исследований в этой области.
С помощью комплекса проведены численные исследования свойств композиционных материалов, получены оценки упругих параметров в зависимости от условий расчета.
Для дальнейшего развития МКСЭ и программного комплекса требуются разработка методов разделения областей, эффективных алгоритмов решения краевых задач в областях простой формы и эффективная реализация этих алгоритмов на многопроцессорных и параллельных ЭВМ.
СПИСОК ЛИТЕРАТУРЫ
1. Страховская Л. Г., Федоренко Р. П. Об одном варианте метода конечных элементов // ЖВМиМФ. - 1979. - Т. 19, № 4. - С. 950-960.
2. Страховская Л. Г., Федоренко Р. П. Расчет диффузии в многосвязной области методом конечных суперэлементов // Препринт ИПМ им. М.В. Келдыша АН СССР - 1987. - № 171. - 26 с.
3. Страховская Л. Г., Федоренко Р. П. Расчет напряжений в композитном теле методом конечных суперэлементов // Препринт ИПМ им. М.В. Келдыша АН СССР - 1994. - № 97. - 26 с.
4. Г а л а н и н М. П., Савенков Е. Б. Совместное использование метода конечных элементов и метода конечных суперэлементов // Препринт ИПМ им. М.В. Келдыша РАН. - 2004. - № 13. - 34 с.
5. Галанин М. П., Савенков Е. Б. О связи метода конечных суперэлементов Федоренко и проекционно-сеточных методов // Препринт ИПМ им. М.В. Келдыша РАН. - 2001. - № 67. - 35 с.
6. G a l a n i n M., Savenkov E. Fedorenko finite superelement method as special Galerkin approximation // Mathematical Modelling and Analysis. - 2002. -V. 7, № 1.-P. 41-50.
7. Галанин М. П., Савенков Е. Б. К обоснованию метода конечных суперэлементов Федоренко // ЖВМиМФ. - 2003. - Т. 43, № 5. - С. 713-729.
8. Галанин М., Савенков Е. Метод конечных суперэлементов для задачи о скоростном скин-слое // Труды международной конференции по вычислительной математике МКВМ-2004. 21-25 июня, 2004 г., Новосибирск, Академгородок, Россия. - С. 455-460.
9. Г а л а н и н М. П., С а в е н к о в Е. Б. Метод конечных суперэлементов для задачи о скоростном скин-слое // Препринт ИПМ им. М.В. Келдыша РАН. -2004. - № 3. - 32 с.
10. Бахвалов Н. С., Панасенко Г. П. Осреднение процессов в периодических средах. - М.: Наука, 1984. - 352 с.
11. Кристенсен Р. Введение в механику композитов. - М.: Мир, 1982. -334 с.
12. Скворцов В. Р. Основы механики композитов. - СПб.: Изд. центр СПбГМТУ, 1995. - 111 с.
13. A b o u d i J. Micromechanical analysis of composites by the method of cells // Appl. Mech. Rev. - 1989. - Vol. 42, № 7. - P. 193-221.
14. A b o u d i J. A continium theory for fiber-reinforced elastic-viscoplastic composites // Int. J. Engng. Sci. - 1982. - V. 20, № 5. - P. 605-621.
15. Fleming W. J., D o w s o n A. L. Prediction of the fatigue life of an aluminium metal matrix composite using the theory of cells // Science and Engineering of Composite Materials. - 1999. -V. 8, № 4. - P. 181-189.
16. Ефиме н ко А. В., Зарубин В. С., Кувыркин Г. Н. О двусторонних оценках термомеханических и теплофизических свойств неоднородных материалов // Динамика, прочность и износостойкость машин. Международный журнал на электронных носителях. - 1996. № 2. - С. 3-7.
17. З а р у б и н В. С., Кувыркин Г. Н. Прогнозирование теплофизических и термоупругих характеристик композитов // Вестник МГТУ им. Н.Э. Баумана. Сер. Машиностроение. - 1994. - С. 78-83.
18. Efimenko A. V., Kuvyrkin G. N. New estimates of the effective elastic moduli of two-component composites // Izv. RAN. Mekhanika Tverdogo Tela. -1994. - Vol. 29, № 1. - P. 18-26.
19. Бураго Н. Г., Галанин М. П., Кувыркин Г. Н. Основные вариационные принципы механики сплошной среды. - М.: Изд-во МГТУ им. Н.Э. Баумана, 2004. - 79 с.
20. Галанин М. П., Савенков Е. Б., Т е м и с Ю. М. Метод конечных суперэлементов для задач теории упругости. // Препринт ИПМ им. М.В. Келдыша РАН. -2004. -№ 38. - 38 с.
21. Г а л а н и н М. П., Савенков Е. Б., Т е м и с Ю. М., Щеглов И. А., Яковлев Д. А. Применение метода конечных суперэлементов для расчета напряженно-деформированного состояния композиционных материалов // Препринт ИПМ им. М.В. Келдыша РАН. - 2006. - № 39. - 31 с.
22. Г а л а н и н М. П., Щеглов И. А. Разработка и реализация алгоритмов трехмерной триангуляции сложных пространственных областей: прямые методы // Препринт ИПМ им. М.В. Келдыша РАН. - 2006. № 10. - 32 с.
23. Г а л а н и н М. П., Щеглов И. А. Разработка и реализация алгоритмов трехмерной триангуляции сплошных пространственных областей: итерационные методы // Препринт ИПМ им. М.В. Келдыша РАН. - 2006. - № 9. - 32 c.
24. М а р ч у к Г. И., А г о ш к о в В. И. Введение в проекционно-сеточные методы. - М.: Наука, 1981. - 416 с.
25. С ь я р л е Ф. Метод конечных элементов для эллиптических задач. - М.: Мир, 1980.- 512 с.
Статья поступила в редакцию 9.02.2007
