CC BY
216
УДК 539.422.24
Вестник СибГАУ 2014. № 4(56). С. 249-255
АНАЛИЗ РОСТА РАССЛОЕНИЙ В КОМПОЗИТНЫХ КОНСТРУКЦИЯХ
С. А. Чернякин, Ю. В. Скворцов
Самарский государственный аэрокосмический университет имени академика С. П. Королева (Национальный исследовательский университет) Российская Федерация, 443086, г. Самара, Московское шоссе, 34 E-mail: [email protected]
Проводится численное исследование разрушения образцов из композиционного материала с расслоением в квазистатической и усталостной постановках. Цель работы заключается в изучении закономерностей развития расслоений в многослойных композиционных материалах при статической и циклической нагрузке. Для моделирования роста расслоения в реальной конструкции ввиду сложности геометрии и условий нагружения наиболее подходящим инструментом является метод конечных элементов. При этом для исследования рассматриваемой задачи удобно использовать модель когезионной зоны. Основываясь на контактном взаимодействии слоёв-оболочек, предлагается подход к моделированию конструкций из композиционных материалов с расслоением. Верификация данного подхода проводится на образцах-балках DCB (Double Cantilever Beam) и ENF (End Notched Flexure), конфигурации которых описаны в стандартах ASTM. Сравнивая полученные результаты после численного моделирования задачи о квазистатическом разрушении образца из композиционного материала с расслоением с данными зарубежных авторов, отмечается хорошее согласование.
Поскольку МКЭ-пакет ANSYS® позволяет моделировать лишь квазистатический рост расслоений, авторами был разработан алгоритм усталостного накопления повреждений на фронте расслоения, использующий стандартные возможности программы. Алгоритм основан на внесении повреждений в контактные элементы с помощью приложения дополнительных перемещений к узлам элементов в области предполагаемого развития расслоения. Проверка работоспособности предложенного алгоритма рассматривается на примере усталостного разрушения DCB-образца. Погружение здесь осуществляется двумя противоположно направленными силами, приложенными к верхнему и нижнему слоям балки. Данные слои изготовлены из однонаправленного материала, армированного волокнами, которые ориентированы вдоль оси балки. Начальная длина расслоения равна расстоянию от точки приложения нагрузки до фронта дефекта. Полученные результаты хорошо согласуются с данными других авторов. Па основании полученных данных сделан вывод о работоспособности и достаточной точности предложенного подхода к моделированию задачи о квазистатическом и усталостном разрушении образцов из композиционных материалов с расслоением.
Ключевые слова: расслоение, усталость, модель когезионной зоны, конечно-элементный анализ, многослойные композиты.
Vestnik SibGAU 2014, No. 4(56), P. 249-255
ANALYSIS OF DELAMINATION PROPAGATION IN COMPOSITE STRUCTURES
S. A. Chernyakin, Yu. V. Skvortsov
Samara State Aerospace University 34, Moskovskoe shosse, Samara, 443086, Russian Federation E-mail: [email protected]
Numerical investigation of specimens fracture process of composite materials with delamination under quasistatic and fatigue case is performed. Objective of this work is to study regularities of delamination propagation in laminated composites under static and cyclic loading. Finite element method is the most appropriate apparatus for modeling of delamination propagation process at actual structure considering geometry complexity and loading conditions. Wherein CZM (cohesive zone model) is useful to investigate considered problem. Basing on the contact interaction of shell-layers, approach to modeling of structures of composite materials with delamination is proposed. Verification of this approach is carried out on the beam-samples DCB (Double Cantilever Beam) and ENF (End Notched Flexure), described in international standards ASTM. Good agreement is noted by comparing obtained results with data of foreign authors after numerical modeling of problem of quasistatic sample fracture of composite materials with delamination.
Fatigue damages implementation algorithm at delamination front, using standard capabilities of program, was developed, because CAE-system ANSYS® allows simulating quasistatic growth of delamination. Algorithm bases on damage introducing at contact elements through applying of additional displacements at assumed delamination propagation area. Efficiency examination of proposed algorithm was conducted by the example fatigue fracture of double cantilever beam. Loading is carried out by two opposite forces at the end, which are applied to upper and lower layers of beam. These layers were produced of unidirectional material, reinforced by fibers, which was oriented along the beam axis. Initial delamination length equals to distance from the loading point to the defect front. Obtained results are in good agreement with the data of other authors. Basing on received data we can conclude about efficiencyand sufficient accuracy of proposed algorithm to modeling quasistatic and fatigued fracture of samples of composite materials with delamination.
Keywords: delamination, fatigue, cohesive zone model, finite element analysis, laminated composite.
Введение. В композитных многослойных элементах конструкций часто встречаются дефекты типа расслоения, при которых нарушается связь между слоями, что нередко приводит к снижению их жёсткости и прочности. В условиях нормального отрыва (mode I) или сдвига (mode II) дефекты данного типа ведут себя примерно так же, как и трещины. Поэтому для задач расслоения обычно применяют методы линейной механики разрушения без каких-либо значительных изменений. Для характеристики напряжённо-деформированного состояния вблизи сингулярности можно воспользоваться конечным параметром механики разрушения - интенсивностью освобождения энергии G. Для нахождения данной величины часто используют метод виртуального закрытия трещины VCCT (Virtual Crack Closure Technique) [1-3].
В настоящее время одним из основных инструментов численного решения разнообразных задач математической физики, в том числе и задач моделирования роста расслоений, является метод конечных элементов [4; 5]. При этом использование интерфейсных элементов, основанных на модели когезионной зоны CZM (Cohesive Zone Model) [6-11], позволяет исследовать зарождение и развитие расслоения без задания начального дефекта и перестройки конечно-элементной сетки в процессе его распространения.
Одним из наиболее распространённых видов разрушения многих конструкций является усталостное разрушение, возникающее за счёт постепенного накопления повреждений при циклическом нагружении. Для задач многоцикловой усталости эффективнее применять подход к моделированию повреждений, основанный на использовании стратегии огибающей профиля нагружения [11; 12]. В данном случае в расчётах вместо реальной циклической нагрузки прикладываются лишь её максимальные значения.
В настоящей статье предлагается эффективный способ моделирования роста расслоений в слоистых композитах. Он основывается на CZM-подходе и реализован в CAE-системе ANSYS®. Его эффективность заключается в использовании элементов оболочки в сочетании с прикреплённым контактом, а также во внесении усталостных повреждений путём задания дополнительных перемещений, что позволяет рассматривать расслоения в реальных конструкциях, имеющих сложную геометрию.
Реализация усталостной модели когезионной зоны в МКЭ-программе. Численная схема CZM предназначена для моделирования эффектов так
называемой зоны процесса разрушения FPZ (Fracture Process Zone) во время роста расслоения. Данная зона располагается непосредственно перед фронтом дефекта и имеет конечный размер в направлении распространения расслоения. Основу CZM составляет физическое соотношение, связывающее напряжение сцепления с с раскрытием u взаимодействующих поверхностей расслоения и описывающее, по сути дела, диаграмму деформирования интерфейсного (когези-онного) материала. Наиболее простой и часто используемой на практике является простейшая билинейная диаграмма, состоящая из двух участков (рис. 1).
Рис. 1. Билинейная диаграмма деформирования
Вначале, когда относительное перемещение и не превышает значения и0, материал ведёт себя линейно упруго, и его жёсткость определяется коэффициентом К = с0 / и0. При и = и0 напряжение сцепления достигает локального предела поверхностной прочности с0. Затем идёт участок разупрочнения материала, характеризуемый накоплением необратимых повреждений. В дальнейшем (после разрушения материала) с не изменяется и остаётся равным нулю. Следует отметить, что площадь под диаграммой деформирования численно равна вязкости разрушения Ос.
При квазистатическом нагружении параметр повреждения d можно однозначно связать с максимальным относительным перемещением итах, достигнутым за всю предыдущую историю нагружения:
d =-
u - u„
(1)
Чтобы эта формула была справедливой не только для участка АС, но и для участка ОА, минимальное значение итах следует принять равным и0.
При усталостном нагружении формулировка CZM должна учитывать деградацию свойств когезионного материала как функцию числа циклов N. Как отмечалось ранее, для задач многоцикловой усталости следует использовать стратегию огибающей профиля нагружения. Такая огибающая для отнулевого цикла с постоянной амплитудой показана на рис. 2. Здесь сначала прикладывается квазистатическая нагрузка, линейно возрастающая от нуля до максимального значения Pmax. Затем при постоянной нагрузке активизируется алгоритм накопления усталостных повреждений.
Number of cycles (N) Рис. 2. Огибающая циклического нагружения
Для поцикловой скорости роста расслоения можно воспользоваться следующей формулой [11]:
da _у _ Ifpz
dN .¿PL dN le dN'
(2)
le _ umax u0 ^e
(3)
Изменение параметра повреждения в когезионном элементе е в зависимости от числа циклов нагружения представим следующим образом [12; 14]:
дйе дйе дГа
dN dled dN Первую производную здесь можно выразить как
dde dde du°
(4)
dld due dld
d max d
Согласно соотношению (1)
dde и un
due (ue )2 (u -u„)
max V max' \ c 0)
Из формулы (3) можно записать
due
dle
uc - u0 le
(5)
(6)
(7)
Подставляя выражения (6) и (7) в уравнение (5), получим
dde
(ue )2le
v max '
(8)
Вместо д1ей / дЫ будем использовать среднее значение, которое согласно уравнению (2) равно
dN
le da
lFPZ dN
(9)
Тогда с учётом выражений (8) и (9) соотношение (4) примет вид
dde uc u0 da )
ЗЫ «ах)2йЫ'
В качестве уравнения, описывающего поцикловую скорость роста дефекта, выберем одну из модификаций формулы Пэриса, предложенную в статье [15]:
da _ с (Gma dN \ Gc
b
1-R2
(11)
1е
где I - длина когезионного элемента е в направлении распространения расслоения; 1ей - длина усталостной трещины в элементе е; д1й/ дЫ - средняя скорость роста усталостных трещин; /РР2 - длина зоны процесса разрушения, равная расстоянию от фронта дефекта до точки, где напряжение сцепления достигает максимума [13]. Если предположить, что отношение повреждённой
1е ^ 1е
длины элемента I й к его полной длине I соответствует отношению рассеянной в этом элементе энергии к вязкости разрушения, то можно прийти к выражению
где С и Ь - постоянные материала, определяемые экспериментально по кинетической диаграмме усталостного разрушения; Я = Р:т / Р:ах - коэффициент асимметрии цикла; Р:ш и Р:ах - минимальное и максимальное значения циклической нагрузки; 0:ях - значение интенсивности освобождения энергии, соответствующее нагрузке Р:ах.
Подставляя выражение (11) в уравнение (10), можно получить формулу для вычисления приращения параметра повреждения в рассматриваемом элементе е в зависимости от данного увеличения числа циклов:
Ade (AN) _
и u„
(
(ue )2lF
v max^ F
-C
g.
X
G„
AN.
(12)
Описанная выше формулировка реализована в МКЭ-пакете ANSYS. Когезионную зону можно определять с использованием как интерфейсных, так и контактных элементов. При этом интерфейсными элементами можно соединять только элементы твёрдого тела, а контактными - ещё и оболочечные элементы. Второй подход представляется более предпочтительным, поскольку он даёт возможность рассматривать расслоения в реальных тонкостенных элементах конструкций сложной конфигурации.
Таким образом, эффективный метод моделирования дефектов данного типа заключается в следующем. Сначала задаются две совпадающие поверхности, которые будут представлять отсчётные поверхности двух соединяемых оболочек. Если в качестве таковых выбираются нижние поверхности, то нормали к ним должны быть направлены наружу. Данные поверхности покрываются одинаковыми сетками оболочечных элементов. Затем между ними в месте расположения начального расслоения определяется стандартный контакт. В оставшейся части, представляющей область возможного распространения дефекта, задаётся прикреплённый (bonded) контакт, причём для
b
установления когезионной зоны здесь необходимо с контактными элементами связать С2М-материал.
Алгоритм внесения усталостных повреждений в когезионные контактные элементы рассмотрим на примере I типа нагружения (нормальный отрыв). После выполнения очередного шага во всех узлах ] каждого когезионного элемента е, расположенного вблизи фронта расслоения и ещё не полностью разрушенного (т. е. иетах < ис), определяются относительные перемещения uej. По ним рассчитывается среднее
значение
=1X u
4 j
(13)
характеризующее раскрытие дефекта в центре элемента.
Если ие > ил (где ил - пороговое значение), то для рассматриваемого элемента е с использованием найденного на предыдущем шаге значения иетах определяется параметр повреждения йГ0и. Далее вычисляется напряжение, после чего рассчитывается значение интенсивности освобождения энергии Оетях. Затем по формуле (12) определяется приращение ДйТ.
Это позволяет найти новое значение параметра повреждения:
«С = <М +^е. (14)
Следует отметить, что в настоящей работе для задания изменений параметров повреждения когезион-ных элементов предлагается использовать дополнительные относительные перемещения. Для элемента е в соответствии с формулой (1) данное перемещение вычисляется как
Аи" = -
— d" (и — и )
new V c 0 /
которое назначается каждому из его узлов:
Аи" = Аи".
(15)
(16)
После определения перемещений для всех элементов, находящихся в зоне процесса разрушения, выполняется их осреднение в узлах:
Аи. =1X Аи"
(17)
где суммирование осуществляется по всем элементам из БР2, сходящимся в узле j; п - общее число этих элементов.
Таким образом, каждый шаг накопления усталостных повреждений здесь состоит из двух шагов нагружения. Сначала к верхней и нижней оболочкам в соответствующих узлах в дополнение к уже имеющимся перемещениям прикладываются нормальные перемещения ±Аи. / 2 и для элементов e рассчитываются новые значения и"тах, характеризующие накопленные повреждения. Затем производится удаление такой дополнительной нагрузки.
Численные примеры. Проверка работоспособности предложенного способа моделирования роста расслоений в квазистатической постановке проводилась на образцах-балках DCB и ENF. Благодаря симметрии данных образцов здесь можно ограничиться рассмотрением лишь половины конструкции шириной B/2, задавая в плоскости симметрии соответствующие граничные условия. Во всех случаях моделирование слоёв осуществлялось при помощи оболочечных элементов SHELL181. Во всех примерах нормальная Kn и тангенциальная Kt контактные жёсткости взяты равными 1105 Н/мм3.
Для более точного определения несущей способности, а также для возможности исследования закри-тического поведения вместо приращений сил задавались приращения перемещений и вычислялись соответствующие реакции, характеризующие внешнюю нагрузку.
Результаты расчётов, выполненных в программе сравнивались с решениями плоских задач, представленными в работах [1; 6], где моделирование рассматриваемых образцов осуществлялось с использованием двухмерных элементов плоского напряжённого состояния.
DCB-образец. DCB (Double Cantilever Beam) - образец в виде двухслойной консольной балки с начальным расслоением длиной а0 на свободном конце, нагруженный двумя противоположно направленными силами P, приложенными к верхнему и нижнему слоям (рис. 3). Он предназначен для экспериментального определении вязкости разрушения GIc для I типа деформирования.
Используемые в расчётах свойства материала и размеры DCB-образца представлены в табл. 1 [6]. Предполагается, что слои балки изготовлены из однонаправленного материала, армированного волокнами, которые ориентированы вдоль её оси.
ANSYS®
Ш
РЪ
2h
Рис. 3. DCB-образец
Свойства материала и размеры DCB-образца [6]
Таблица 1
Еь ГПа E2, ГПа G12, ГПа V12 GIc, Н/мм а0, МПа L, мм а0, мм 2h, мм B, мм
135,3 9,0 5,2 0,24 0,28 11,4 100 30 3 20
"
Таблица 2
Свойства материала и размеры ENF-oбpaзцa [6]
Е, ГПа V GIIc, Н/мм т0, МПа Ь, мм а0, мм 2Н, мм В, мм
70 0,33 1,45 57 100 30 3 10
Рис. 4. Зависимость «нагрузка-перемещение» для БСБ-образца
Р, 5
2Ь
1 °о А ^ 1
Рис. 5. БОТ-образец
300
0е--
О 2 4 6 8 10 12
5, ММ
Рис. 6. Зависимость «нагрузка-перемещение» для БОТ-образца
Результаты расчётов приведены на рис. 4 в виде зависимости силы P от перемещения точки её приложения 5. При этом сплошная кривая соответствует решению, полученному с помощью программы ANSYS®, а штриховая линия представляет результаты работы [6].
ENF-образец. ENF (End Notch Flexure) - трёхто-чечный изгиб двухслойной балки с начальным концевым расслоением (рис. 5). Данная конфигурация соответствует II типу деформирования и обычно используется для экспериментального определении GIIc. Свойства материала и размеры исследуемого образца приведены в табл. 2 [6]. Следует отметить, что для такого образца при задании в программе ANSYS® свойств CZM-материала параметр ß (flag for tangential slip under compressive normal contact stress) нужно положить равным 1. Результаты расчётов представлены на рис. 6.
В качестве примера оценки долговечности изделий из композиционного материала с расслоением рассмотрим усталостное разрушение DCB-образца [12]. Он изготовлен из двух слоёв углепластика и имеет следующие размеры: L = 125 мм; а0 = 47 мм; 2h = 5,4 мм; B = 25 мм. Упругие свойства однонаправленных слоёв углепластика приведены в табл. 3. Как и в работе [12], предполагается, что GIc = 0,43 Н/мм и с0 = 30 МПа, а параметры модифицированного закона Пэриса (11) равны C = 16/25 = 0,64 мм/цикл и b = 6,0.
Рассматривается отнулевой цикл нагружения (R = 0) с постоянной амплитудой, когда сила Р изменяется по гармоническому закону в пределах от нуля до максимального значения Pmax = 75 Н.
Таблица 3
Упругие свойства углепластика
= 150 ГПа Vl2 = 0,34 G12 = 4,315 ГПа
E2 = 8,819 ГПа V13 = 0,34 G,3 = 4,315 ГПа
E3 = 8,819 ГПа V23 = 0,38 G23 = 3,200 ГПа
Помимо этих данных примем также, что коэффициент жёсткости когезионного материала К составляет 1105 Н/мм3, а пороговое значение раскрытия и^ равно и0, где
30
= 3 -10-4 мм.
0 К 1105
Следуя стратегии огибающей профиля нагружения, на первом этапе прикладывалась квазистатическая нагрузка, линейно возрастающая от нуля до Р:ах. Затем на втором этапе при постоянной нагрузке Р:ах активизировался алгоритм внесения усталостных повреждений.
На рис. 7 сплошной линией представлены результаты, рассчитанные в программе ANSYS по описанной выше методике, а штриховая линия соответствует решению плоской задачи, полученному в работе [12]. Видно, что, как и в статье [12], в нашем случае долговечность рассматриваемого образца составляет около 70000 циклов.
Рис. 7. Изменение относительного перемещения точек
приложения сил как функция числа циклов нагружения
Заключение. Таким образом, представленные выше результаты позволяют сделать вывод о работоспособности и достаточной точности предложенного подхода моделирования задачи о квазистатическом и усталостном разрушении образцов из композиционных материалов с расслоением. Он реализован в CAE-системе ANSYS® на базе её стандартных возможностей, что упрощает его применение. Использование для моделирования слоёв оболочечных элементов, взаимодействующих посредством прикреплённого контакта, обладающего свойствами CZM-материала, позволяет с его помощью рассматривать расслоения в реальных тонкостенных конструкциях.
Благодарности. Работа выполнена при финансовой поддержке Министерства образования и науки РФ в рамках базовой части государственного задания. Числовые результаты, представленные в работе, были получены с использованием ресурсов суперкомпьютерного центра СГАУ.
Acknowledgements. The study was supported by the Ministry of Education and Science of the Russian Federation. The data given in the work were got with the help of the resources of the supercomputer centre SGAU.
References
1. Liu P. F., Islam M. M. A nonlinear cohesive model for mixed-mode delamination of composite laminates. Composite Structure, 2013, Is. 106, p. 47-56.
2. Krueger R. The virtual crack closure technique: history, approach and applications. Applied Mechanics Reviews, 2013, Is. 57(1-2), p. 109-143.
3. Senthil K., Arockiarajan A., Palaninathan R., Santhosh B., Usha K. M. Defects in composite structures: Its effects and prediction methods - A comprehensive review. Composite Structure, 2013, Is. 106, p. 139-149.
4. Landry B., La Plante G. Modeling delamination growth in composites under fatigue loadings of varying amplitudes. Composites Part, 2012, Is. 43(2), p. 533-541.
5. Muñoz J. J., Galvanetto U., Robinson P. On the numerical simulation of fatigue driven delamination with interface elements. International Journal of Fatigue, 2006, Is. 28(10), p. 1136-1146.
6. Xie D., Waas A. M. Discrete cohesive zone model for mixed-mode fracture using finite element analysis.
Engineering Fracture Mechanics, 2006, Is. 73, p. 17831796.
7. Chandra N., Scheider I., Ghomen K. H. Some issues in the application of cohesive zone models for metal-ceramic interfaces. International Journal of Solids and Structures, 2002, Is. 39 (11), p. 2827-2855.
8. Cornec A., Scheider I., Schwalbe K.H. On the practical application of the cohesive model. Engineering Fracture Mechanics, 2003, Is. 70, p. 1963-1987.
9. De Borst R. Numerical aspects of cohesive zone models. Engineering Fracture Mechanics, 2003, Is. 70, p. 1743-1757.
10. Yang Q. D., Cox B. N. Cohesive models for damage evaluation in laminated composites. International Journal of Fracture, 2005, Is. 133, p. 107-137.
11. Harper P. W., Hallet S. R. A fatigue degradation law cohesive interface elements - Development and application to composite materials. International Journal of Fatigue, 2010, no. 32(11), p. 1774-1787.
12. De Moura MFSF., Gonçalves JPM. Cohesive zone model for high-cycle fatigue of adhesively bonded joints under mode I loading. International Journal of Solids and Structures, 2014, Is. 51, p. 1123-1131.
13. Naghipour P., Bartch M., Vonggenreiter H. Simulation and experimental validation of mixed mode delamination in multidirectional CF/PEEK laminates under fatigue loading. International Journal of Solids and Structures, 2011, Is. 48, p. 1070-1081.
14. Turon A, Costa J, Camanho P. P., Dâvila C. G. Simulation of delamination in composites under high-cycle fatigue. Composites Part A, 2007, no. 38(11), p. 2270-2282.
15. Li C., Teng T., Wan Z., Li G., Rans C. Fatigue delamination growth for an adhesively-bonded composite joint under mode I loading. 27th ICAF Symposium. Jerusalem, June 2013.
© Чернякин С. А., Скворцов Ю. В., 2014
