Расчетная оценка кинетики разупрочнения слоистых композитов на основе стохастического микро-мезо-моделирования
Автор: Шаблей Александра Андреевна, Сапожников Сергей Борисович, Шипулин Леонид Викторович
Рубрика: Численные методы моделирования
Статья в выпуске: 4 т.17, 2017 года.
Бесплатный доступ
Композитные материалы широко применяются в современной промышленности (авиа- и автомобилестроение, строительство и др.). Следует отметить, что современные высокопрочные и легкие композитные материалы (стекло-, угле- или органопластики) существенно дороже стали и алюминиевых сплавов, проявляют свойства упругой и прочностной анизотропии, имеют нелинейность деформирования и сравнительно малые (по сравнению с металлами) деформации разрушения. Наиболее важной задачей при проектировании конструкций из композитных материалов представляется минимизация ее массы без потери прочностных свойств. В первой части работы был разработан программный комплекс на языке С# для создания конечно-элементных моделей однонаправленного волокнистого композитного материала со случайно расположенными прямолинейными волокнами. Программный комплекс позволяет построить композитный материал с укладкой [0°/90°/0°], задавая входные параметры: размеры слоев и волокон в 90 слое. Программный комплекс выдает текстовый файл с расширением *.lgw, который в дальнейшем используется в численных расчетах. Во второй части работы представлены методика моделирования однонаправленного волокнистого композитного материала со случайно расположенными прямолинейными волокнами на микро- и мезоуровне. Задачи разного уровня рассматривали на примере задачи растяжения композитной панели. Механические характеристики в мезомодели подбирались на основании данных из модели микроуровня. В мезомодели постепенное разрушение серединного слоя (90° слой) задавали с помощью критерия Stochastic Failure (случайный разброс прочностных свойств по Мотту). Расчетные кривые деформирования, полученные в микро- и мезомоделях, хорошо коррелируют друг с другом.
Стеклопластик, микротрещины, разупрочнение, метод конечных элементов, стохастическая микро- и мезомеханика, слоистый композит
Короткий адрес: https://sciup.org/147151765
IDR: 147151765 | DOI: 10.14529/engin170406
Текст научной статьи Расчетная оценка кинетики разупрочнения слоистых композитов на основе стохастического микро-мезо-моделирования
Около 20 лет назад в научной литературе по композитным материалам произошла смена направления исследований с прогноза нагрузок разрушения на анализ кинетики нелинейного деформирования вплоть до полного разрушения. Для слоистых композитов это проявилось в активизации работ по анализу появления и развития трансверсальных трещин в волокнистых композитах как макрооднородных анизотропных средах. На основе таких гомогенизированных моделей делались попытки расчета систем трещин (аналитические модели), а также вводились заданные вариации механических свойств и оценивалось влияние таких вариаций на кинетику процесса трещинообразования. Более адекватные численные микроструктурные модели, в которых рассматривали случайное распределение волокон в матрице, строились без введения дополнительных гипотез о характере разброса механических свойств в гомогенизированных моделях. На таких моделях были получены реалистичные картины трансверсальных микротрещин, однако трудоемкость вычислений была весьма велика. Необходимость корректного учёта микроповреждений в слоях, приводящих к снижению общей жёсткости конструкции и к перераспределению напряжений в слоях – необходимый этап в расчете кинетики разупрочнения и прогнозе полного разрушения.
На поиск решения данной проблемы, а именно учёта микроповреждений и оценки их кинетики, направлен международный проект WWFE-III (World-Wide Failure Exercise). Организаторы проекта по результатам двух предыдущих этапов выявили необходимость определения степени корректности методов моделирования различных аспектов микроповреждений в композитных материалах [1]. Участники проекта исследовали растрескивание ортогонального слоя в слоистых композитах, составленных из чередующихся слоев с ортогональными направлениями укладки волокон (структуры типа 0m/90n, где m и n – число соответствующих слоев) [2].
В рамках проекта было проведено сравнение двенадцати основных международно признанных методов моделирования повреждений в композитах. Результаты показали, что в некоторых моделях были учтены лишь частные аспекты повреждения композитов [3, 4]. Но при создании других моделей были предприняты попытки учесть практически все возможные аспекты повреждения композитов [5–9] По результатам работ участников организаторы заключили, что необходимо продолжать исследования накопления микроповреждений [9].
Алгоритмы для моделирования микроструктуры применяются исследователями уже достаточно давно [10–13], но вопрос не потерял актуальности и в настоящее время [14–16].
В данной работе геометрия микромодели создана при помощи специализированной программы (разработчики Л.В. Шипулин, А.В. Игнатова), позволяющей составить набор команд в виде текстового файла *.lgw для построения двумерной геометрии фиксированной толщины в пакете ANSYS. В этой программе можно задать необходимые геометрические параметры образца, диаметры волокон, минимальное расстояние между волокнами, объемную долю волокон в материале, коэффициент вариации – степень случайного разброса укладки волокон по толщине 90° слоя, а также тип упаковки волокон в образце для регулярной сетки волокон.
Для проведения исследований необходимо построить трехмерную модель в системе конечноэлементного анализа ANSYS. Построение такой трехмерной модели вручную сопряжено с рядом сложностей. Во-первых, это большое количество волокон (до нескольких тысяч), присутствующих во внутреннем слое слоистой системы, и, следовательно, большой объем работы по подготовке конечно-элементной модели. А во-вторых, построенная вручную конечно-элементная модель является лишь одним вариантом случайного распределения волокон и для построения дру- гих вариантов распределения волокон придется создавать модель заново.
Построение расчетной схемы для программного модуля
Для разработки программного модуля генерации геометрической модели среды предлагается расчетная схема, представленная на рис. 1. В соответствии с расчетной схемой модель среды представляет собой параллелепипед, состоящий из трех параллельных слоев. Эти слои представ- ляют пакет с направлениями укладки волокон [0/90/0]. Внешние слои имеют толщину b1 и b3, а внутренний слой толщиной b2, содержит неизвестное количество n параллельных волокон различного диаметра d1, d2, ..., dn. Между волокнами всегда должно быть некоторое минимальное расстояние lmin. Ширина a и длина c для каждого слоя одинаковая, волокна также проходят на всю глубину модели и имеют длину c. В программном пакете предусмотрена возможность реализации различных схем упаковки волокон: случайной по равномерному закону распределения, кубической и гексагональной в двух направлениях.

Рис. 1. Расчетная схема
Структура разрабатываемого программного модуля
В соответствии с поставленной задачей был разработан программный модуль, предназначенный для виртуального геометрического построения слоев и волокон и последующей записи этих данных в текстовый файл формата *.lgw. Структурно-программный модуль состоял из следующих основных блоков:
-
1. Блок ввода исходных данных. В этом блоке пользователь вводит в компоненты TextBox исходные данные, которые были описаны выше, а именно: a , b 1 , b 2 , b 3 , c , d и l min . Также в форме пользователь должен указать схему упаковки волокон, которую необходимо получить в геометрической модели: случайную по равномерному закону распределения волокон, кубическую или
- гексагональную (горизонтальную или вертикальную). Для выбора схемы упаковки применяется компонент RadioButton. По завершении ввода исходных данных и выбора схемы упаковки пользователь запускает расчет нажатием соответствующей кнопки.
-
2. Формирование геометрической модели среды. Поскольку в каждом поперечном к направлению волокон сечении расположение слоев и волокон не будет отличаться, то целесообразно рассматривать одно сечение, и в этом сечении производить все геометрические построения, а уже потом распространять это сечение на всю трехмерную модель. Таким образом, задача построения геометрической модели сводится к построению плоского сечения, перпендикулярного направлению волокон.
-
3. Создание файлов с результатами. После создания геометрической модели необходимо все сформированные данные сохранить в файле формата *.lgw для импорта в ANSYS. В соответствии с формальными правилами составления файла для командной строки файл должен начинаться со строки «/PREP7» и заканчиваться строками «vovlap, all» и «vglue, all». Между начальными и заключительными строками находятся команды создания конечно-элементной модели. Для создания в конечно-элементной модели прямоугольного параллелограмма (слоя) применяется команда «BLC4,x,y,a,b,c», где x и y – начальные координаты параллелограмма, a, b и c – его ширина, высота и глубина соответственно. Для создания цилиндра (волокна) применяется команда «CYL4,x 0 ,y 0 ,r,,,,c», где x 0 и y 0 – координаты центра образующей окружности цилиндра, r – радиус этой окружности, c – длина цилиндра. Полученный файл формата *.lgw содержит информацию обо всей конечно-элементной модели трехслойного волокнистого материала.
Для построения поперечного сечения в первую очередь построили сечения слоев, представленные прямоугольниками, ширина a которых одинакова, а толщина каждого слоя b 1 , b 2 и b 3 индивидуальна. Слои располагаются один над другим. Во внутреннем слое должна быть сформирована совокупность волокон, если рассматривать одно сечение, то волокна будут представлены окружностями. При этом следует помнить, что необходимо разработать порядок создания различных схем упаковок волокон: случайную по равномерному закону распределения, кубическую и гексагональную. Геометрически задача сводится к заполнению прямоугольника окружностями по заданной схеме упаковки. В информационном виде волокна представлены в виде массивов координат центров и значений диаметра волокна.
Для кубической и гексагональной (в двух направлениях) схем упаковки внутренний слой в виде прямоугольника заполняется окружностями по общеизвестным формулам, в результате чего формируются массивы координат центров и диаметров волокон. Более сложной является задача случайной упаковки волокон по равномерному закону распределения. Дело в том, что при создании таких волокнистых материалов задают объемную долю волокон, следовательно, примерно такой же объем волокон необходимо получить и в виртуальной модели. Решение этой задачи выполнено следующим образом. Пользователь задает дополнительный исходный параметр – объемную долю волокон во внутреннем слое, например, 0,5. В начальный момент времени в виртуальной модели текущий объем волокон равен нулю, по мере генерирования волокон этот объем увеличивается и приближается к заданному. Генерация положения каждого отдельного волокна заключается в следующем: создается окружность, координаты центра которой генерируются по равномерному закону распределения, а диаметр – по нормальному закону распределения с заданным коэффициентом вариации. После генерации этих параметров происходит проверка на пересечение с имеющимися волокнами, и, если такие пересечения имеются – параметры генерируются заново. Если сгенерированное волокно не пересекается с существующими волокнами – оно фиксируется в массиве волокон, а текущий объем волокон увеличивается на величину объема этого волокна. В том случае, если текущий объем волокон превысил требуемый пользователем объем, создание волокон останавливается и модель считается готовой.
Помимо вышеуказанного файла формируются и другие:
-
– файл с данными о значениях диаметров каждого волокна для статистической оценки выборки;
-
– файл со справочными данными о результатах построения геометрической модели (объемная доля волокон);
-
– графический файл с визуализацией геометрической модели (одного сечения).
Реализация модуля моделирования слоистого материала
По описанным выше структурным блокам разработана программа на языке C#. Интерфейс программного модуля представлен на рис. 2. При использовании программного модуля пользователь задает исходные данные – параметры геометрии и тип упаковки волокон (см. рис. 2). После ввода исходных данных пользователь запускает расчет, в ходе которого формируется файл *.lgw для импорта в ANSYS. Пример такого файла представлен на рис. 3.
Высота 1 слоя (Ь1). мм
Толщина (с), мм
Параметры геометрам:
Ширина (а), мм
Тип упаковки волокон:
А случайная по нормальному закону 0,05 коэффициент вариации
Высота 2 слоя (Ь2). мм © регулярная кубическая
Высота 3 слоя (ЬЗ). мм © регулярная гексагональная вдоль
регулярная гексагональная вдоль
Диаметр (d), мм Расстояние между W,WI волокнами (min), мм 0,5 Доля волокон в объеме
Запустить расчет
Рис. 2. Интерфейс модуля моделирования слоистого волокнистого материала
Рис. 3. Файл, созданный в программном модуле и импортируемый в ANSYS
Также в результате расчета получен файл с визуализацией результатов моделирования волокнистого материала (рис. 4а). На рис. 4а видно два внешних слоя и один внутренний, в котором по случайному закону распределения расположены волокна. Для иллюстрации возможностей моделирования других схем упаковки волокон при тех же самых геометрических параметрах модели проведено построение кубической (тетрагональной) упаковки волокон (рис. 4б), гексагональной продольной (рис. 4в) и гексагональной поперечной (рис. 4г).




г)
а) б) в)
Рис. 4. Визуализация моделирования слоистого волокнистого материала для различных схем упаковки волокон: а – случайная упаковка по равномерному закону распределения;
б – кубическая упаковка волокон; в, г – гексагональная упаковка, продольная и поперечная
Микромоделирование
В дальнейших расчетах использовали модель со случайной упаковкой волокон. Такой выбор объясняется тем, что данный тип упаковки соответствует распределению волокон в реальном однонаправленном волокнистом композитном материале (ОВКМ).

Рис. 5. Геометрическая структура микромодели ОВКМ со свойствами, закладываемыми на этапе проектирования

Рис. 6. Структура ОВКМ (стеклопластик)
Для получения механических характеристик 90 ° слоя ОВКМ была сгенерирована конечно-элементная микромодель с помощью разработанной программы, которая позволяет получить текстовый файл с командами на языке APDL и загрузить в пакет ANSYS Mechanical APDL (рис. 5). Заданные параметры для построения геометрии представлены на рис. 2.
Программа воспроизводит геометрию ОВКМ с достаточной степенью реализма. Для качественного сравнения на рис. 6 представлена реальная микроструктура ОВКМ [17].
Назначение внутренних параметров компонентов моделей. Расчеты деформирования и разрушения полученной геометрической модели производили в пакете ANSYS Workbench (explicit formulation). Свойства материалов представлены в табл. 1, которые были взяты из библиотеки конечно-элементного пакета. Верхний и нижний 0° слои – это высокомодульный однонаправленный композит (Epoxy E-Glass UD). Материал матрицы – эпоксидная смола (Resin Epoxy), материал волокон – стекло (Glass) (рис. 7).
Рис. 7. Микромодель ОВКМ
Таблица 1
Характеристика материалов
Материал |
Плотность, кг/м3 |
Е, МПа |
μ |
Волокно |
2500 |
73 000 |
0,25 |
Матрица |
1160 |
3780 |
0,35 |
Однонаправленный композит* |
2000 |
40 000 |
0,3 |
* свойства указаны в направлении оси x и плоскости yx.
На рис. 8 представлен фрагмент КЭ сетки микромодели (1/7 от полной модели). Идеальный контакт между волокном и матрицей задан таким образом, что «отклеивание» волокна от матрицы не рассматривается как отдельный механизм разрушения элементов матрицы. Между слоями 0 ° и 90 ° задан контакт «bonded».
На рис. 9 представлены граничные условия для полной расчетной микромодели ОВКМ (считаем деформирование плоским, перемещения узлов в направлении оси z равны нулю). Литерами A, B, C и D обозначены граничные условия, перемещения, заданные для данной модели. К правому краю (C) прикладывали линейно возрастающее во времени перемещение вдоль оси X от нуля до 0,5 мм.

Рис. 8. Конечно-элементная микромодель

Рис. 9. Граничные условия для полной расчетной микромодели ОВКМ: A, B – ограничение перемещения по оси z; C – ограничение перемещения по оси x;
D – перемещение по оси x
Время расчёта (перемещение торца модели от 0 до 0,5 мм за 0,5 мc) задавали из условия малости доли кинетической энергии модели в полной энергии деформирования (квазистатическая постановка задачи в явной схеме интегрирования в пакете ANSYS). На каждом шаге решалась упругая задача, при достижении первым главным напряжением в матрице заданного предела прочности (50 МПа) конечный элемент из системы удалялся (технология death of finite element, принятая в ANSYS).
В результате расчета получены картины разрушения слоя ОВКМ. Распределение деформаций по оси X показано на рис. 10 для времени 0,5 мс.

Рис. 10. Картина растрескивания микромодели ОВКМ и распределение перемещений по оси X, перемещения указаны в мм
Мезомоделирование
Микромодели требуют больших вычислительных ресурсов, что не позволяет такой подход использовать для оценки поведения реальных объектов. Поэтому ниже рассмотрен способ перехода от свойств микромодели к мезомоделированию. В программном пакете 3D моделирования (SolidWorks) был создан образец ОВКМ, размеры которого соответствовали размерам образца микромодели ОВКМ (рис. 11). Основным отличием являлось лишь то, что волокна в 90° слое не моделировали явным образом. Механические свойства ОВКМ приведены в табл. 2 (использована база данных пакета ANSYS Workbench для стеклопластика Epoxy E-Glass UD с объемной долей волокон 50 %).
Аналогичным образом геометрическая мезомодель ОВКМ загружается в пакет ANSYS Workbench.
Таблица 2
Механические характеристики материала
Ортотропная упругость |
Ортотропные преде л ы прочности |
||
Модуль Юнга в нап равлении X, MPa |
40000 |
Растяжение в направлении X, MPa |
1100 |
Модуль Юнга в нап равлении Y, MPa |
10000 |
Растяжение в направлении Y, MPa |
33 |
Модуль Юнга в направлении Z, MPa |
10000 |
Растяжение в направлении Z, MPa |
33 |
Коэффициент Пуассона XY |
0,30 |
Сжатие в направлении X, M Pa |
–1000 |
Коэффициент Пуассона YZ |
0,40 |
Сжатие в направлении Y, M Pa |
–100 |
Коэффициент Пуассона XZ |
0,30 |
Сжатие в направлении Z, MPa |
–100 |
Модуль сдвига XY, MPa |
5000 |
Сдвиг XY, MPa |
50 |
Модуль сдвига YZ, MPa |
3846,2 |
Сдвиг YZ, MPa |
30 |
Модуль сдвига XZ, MPa |
5000 |
Сдвиг XZ, MPa |
50 |
Для материала серединного слоя задали случайный разброс прочностны х свойств п р и помо-

Рис. 11. Геометрическая структура мезомодели ОВКМ со свойствами, закладываемыми на этапе расчетов
щи критерия Stochastic Failure (распределение Мотта, заложенное в пакет ANSYS). Данный кри- терий позволяет моделировать разброс прочностных свойств, близкий к разбросу прочностных свойств в материале с дефектами. Значение коэффициента стохастической вариации ɣ подбирали

Рис. 12. Расчет усилия, воспринимаемого слоем 90° микромодели
таким образом, чтобы максимальная сила, действующая в 90° сло е , была равна 12,5 Н (как в микромодели, рис. 12). В результате серии расчетов был получен коэффициент стохастической вариац и и ɣ, равны й 11.
На рис. 13 изображен фрагмент КЭ сетки мезомодели ОВКМ. Размеры элементов в 90° слое выбраны таким образом, чтобы количество элементов по толщине 90° слоя совпадало с количеством волокон по толщине слоя в микромодели. Здесь принято, что КЭ в мезомодели должен быть «представительным объектом» микромодели. В связи с этим упругие свойства КЭ в мезомодели совпадают с упругими свойствами ОВКМ, а прочностные имеют разброс. Все граничные условия и схема нагружения для мезомодели ОВКМ аналогичны граничным условиям и схеме нагружения микромодели ОВКМ (рис. 7).
На рис. 14 представлены результаты расчета пер е мещений для мезомо д ели (время 5 - 10-4с).
Разрушение (картина трещин) имеет тот же характер, что и у микромодели.
На рис. 15 показаны диаграммы деформирования (среднее напряжение – средняя деформация) микро- и мезомоделей (по 5 реализаций). На основании представленных диаграмм можно отметить качественное и количественное совпадение кривых с учетом, естественно, стохастики процесса разрушения. Из этого следует, что замена микромодели, расчет которой требует боль- ших вычислительных мощностей, мезомоделью возможна и позволяет переходить к расчетам более объемных моделей конструктивных элементов, получив практически в шесть раз меньшее время расчета одного варианта (при одной и той же конфигурации вычислительного кластера, состоящего из 48 ядер). Полученные данные по времени приведены в табл. 3. Исследования выполнены с использованием суперкомпьютерных ресурсов ЮУрГУ [18].
Рис. 13. Конечно-элементная мезомодель (верхний и нижний слои – высокомодульный однонаправленный композит с ориентацией волокон 0 ° ; серединный слой - высокомодульный однонаправленный композит, с волокнами, расположенными под углом 90° относительно оси x)


Рис. 14. Картина растрескивания мезомодели ОВКМ и распределение перемещений по оси X, перемещения указаны в мм
Таблица 3
Время |
Микромодель |
Мезомодель |
Процессорное время, мин |
115 |
20 |
Общее время, мин |
214 |
22 |
Заключение
В работе исследовано напряженное состояние и разрушение полимерной матрицы слоя с ориентацией 90 ° со случайным расположением волокон (микромодель) при квазистатическом нагружении ортогонально армированного стеклопластика типа [0 ° /90 ° /0 ° ].
При помощи метода конечных элементов ( explicit formulation , технология death of finite element ) рассмотрено накопление трансверсальных микроповреждений в структуре стеклопластика с объемной долей волокон около 50 % в плоской постановке в составе ортогонально армированного композита со случайным расположением волокон в слое 90 ° . В качестве критерия прочности при удалении конечных элементов использована величина первого главного напряжения в полимерной матрице.
Показано, что рассеянные микроповреждения в структуре однонаправленного стеклопластика и возникающие за счет их слияния крупные трещиноподобные дефекты, существенно снижают жесткость трансверсального слоя, что отражается на диаграмме деформирования композита в виде спадающей кривой.
Предложена оригинальная схема замещения микромодели эквивалентной по жесткости стохастической мезомоделью, в которой конечные элементы имеют различную прочность в соответствии с известным стохастическим распределением Мотта и разрушаются, снижая упругие свойства композита. При этом время расчета диаграммы деформирования существенно сокращается.
Исследование выполнено при финансовой поддержке РФФИ в рамках научного проекта № 17-08-01024 А.
Список литературы Расчетная оценка кинетики разупрочнения слоистых композитов на основе стохастического микро-мезо-моделирования
- Kaddour A.S., Hinton M.J., Smith P.A., Li S. The Background to the Third World-Wide Failure Exercise//Journal of Composite Materials, 2013, vol. 47, no. 20-21, pp. 2417-2426.
- Pinho S.T., Vyas G.M., Robinson P. Response and Damage Propagation of Polymer-Matrix Fibre-Reinforced Composites: Predictions for WWFE-III (Part A)//Journal of Composite Materials, 2013, vol. 47, no. 20-21, pp. 2595-2612.
- Kashtalyan M., Soutis C. Predicting Residual Stiffness of Cracked Composite Laminates Subjected to Multi-Axial Inplane Loading//Journal of Composite Materials, 2013, vol. 47, no. 20-21, pp. 2513-2524.
- Forghani A., Zobeiry N., Poursartip A., et al. A Structural Modelling Framework for Prediction of Damage Development and Failure of Composite Laminates//Journal of Composite Materials, 2013, vol. 47, no. 20-21, pp. 2553-2573.
- Laurin F., Carrere N., Huchette C., et al. A Multiscale Hybrid Damage and Failure Approach for Strength Predictions of Composite Structures//Journal of Composite Materials, 2013, vol. 47, iss. 20-21, pp. 2713-2747.
- Daghia F., Ladev‘eze P. Identification and Validation of an Enhanced Mesomodel for Laminated Composites within the WWFE-III//Journal of Composite Materials, 2013, vol. 47, no. 20-21, pp. 2675-2693.
- Flatscher Th., Schuecker C., Pettermann H.E. A Constitutive Ply Model Predicting Stiffness Degradation as Well as Inelastic Strain Accumulation and Its Application to WWFE-III (Part A)//Journal of Composite Materials, 2013, vol. 47, no. 20-21, pp. 2575-2593.
- Sapozhnikov S.B., Cheremnykh S.I. The Strength of FRP under a Complex Loading//Journal of Composite Materials, 2013, vol. 47, no. 20-21, pp. 2525-2552.
- Pinho S.T., Vyas G.M., Robinson P. Response and Damage Propagation of Polymer-Matrix Fibre-Reinforced Composites: Predictions for WWFE-III (Part A)//Journal of Composite Materials, 2013, vol. 47, no. 20-21, pp. 2595-2612.
- Kaddour A.S., Hinton M.J., Smith P.A., Li S. A Comparison between the Predictive Capability of Matrix Cracking, Damage and Failure Criteria for Fibre Reinforced Composite Laminates: Part A of the Third World-Wide Failure Exercise//Journal of Composite Materials, 2013, vol. 47, iss. 20-21, pp. 2749-2779.
- Berryman J.G. Random Close Packing of Hard Spheres and Disks//Physical Review A, 1983, vol. 27, no. 2, pp. 1053-1061.
- Clarke A.S., Wiley J.D. Numerical Simulation of the Dense Random Packing of a Binary Mixture of Hard Spheres: Amorphous Metals. Physical Review B, 1987, vol. 35, no. 14, pp. 7350-7357.
- He D., Ekere N.N., Cai L. Computer Simulation of Random Packing of Unequal Particles//Physical Review E, 1999, vol. 60, no. 6, pp. 7098-7104.
- Sanei S.H.R., Barsotti E.J., Leonhardt D., Fertig R.S. Characterization, Synthetic Generation, and Statistical Equivalence of Composite Microstructures//Journal of Composite Materials, 2017, vol. 51, no. 13, pp. 1817-1829.
- Wang W., Dai Y., Zhang C., Gao X., Zhao M. Micromechanical Modeling of Fiber-Reinforced Composites with Statistically Equivalent Random Fiber Distribution//Materials, 2016, vol. 9, no. 8, article number 624.
- Аношкин А.Н. Микромеханический анализ неупругого деформирования однонаправленных волокнистых композитов при многоосном нагружении и сдвиге//Механика композит. материалов. 2003 Т. 39, № 5. С. 575-586.
- Yu H., Longana M.L., Jalalvand M., Wisnom M.R., Potter K.D. Pseudo-Ductility in Intermingled Carbon/Glass Hybrid Composites with Highly Aligned Discontinuous Fibres//Composites: Part A, 2015, vol. 73, pp. 35-44.
- Костенецкий П.С., Сафонов А.Ю. Суперкомпьютерный комплекс ЮУрГУ//Параллельные вычислительные технологии (ПаВТ’2016): сб. ст. 10-й Междунар. науч. конф. (28 марта -1 апреля 2016 г.). Архангельск, 2016. С. 561-573.