Analysis of delamination propagation in composite structures

全文:

Введение. В композитных многослойных элементах конструкций часто встречаются дефекты типа расслоения, при которых нарушается связь между слоями, что нередко приводит к снижению их жёсткости и прочности. В условиях нормального отрыва (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. Билинейная диаграмма деформирования Вначале, когда относительное перемещение u не превышает значения u0, материал ведёт себя линейно упруго, и его жёсткость определяется коэффициентом K = о0 / u0. При u = u0 напряжение сцепления достигает локального предела поверхностной прочности о0. Затем идёт участок разупрочнения материала, характеризуемый накоплением необратимых повреждений. В дальнейшем (после разрушения материала) о не изменяется и остаётся равным нулю. Следует отметить, что площадь под диаграммой деформирования численно равна вязкости разрушения Gc. При квазистатическом нагружении параметр повреждения d можно однозначно связать с максимальным относительным перемещением umax, достигнутым за всю предыдущую историю нагружения: d =u - u„ (1) Чтобы эта формула была справедливой не только для участка AC, но и для участка OA, минимальное значение umax следует принять равным u0. 250 Технологические процессы и материалы При усталостном нагружении формулировка CZM должна учитывать деградацию свойств когезионного материала как функцию числа циклов N. Как отмечалось ранее, для задач многоцикловой усталости следует использовать стратегию огибающей профиля нагружения. Такая огибающая для отнулевого цикла с постоянной амплитудой показана на рис. 2. Здесь сначала прикладывается квазистатическая нагрузка, линейно возрастающая от нуля до максимального значения Pmax. Затем при постоянной нагрузке активизируется алгоритм накопления усталостных повреждений. Number of cycles (N) Рис. 2. Огибающая циклического нагружения Для поцикловой скорости роста расслоения можно воспользоваться следующей формулой [11]: da _ У _ lFPZ dld dN JfZ dN le dN ’ (2) le umax u0 le (3) Изменение параметра повреждения в когезионном элементе e в зависимости от числа циклов нагружения представим следующим образом [12; 14]: dde dde dled dN dled dN Первую производную здесь можно выразить как dde dde du* (4) dld du* dl*, d max d Согласно соотношению (1) dd* uc u0 du* (u* )2 (u - u0 ) max V max Z V c 0 / Из формулы (3) можно записать du* uc - u0 l* (5) (6) (7) Подставляя выражения (6) и (7) в уравнение (5), получим dde dl* (ue )2le v max Z (8) Вместо dl*d / dN будем использовать среднее значение, которое согласно уравнению (2) равно dN le da lFPZ dN (9) Тогда с учётом выражений (8) и (9) соотношение (4) примет вид dd* dN da (umax)2 lFPZ dN (10) В качестве уравнения, описывающего поцикловую скорость роста дефекта, выберем одну из модификаций формулы Пэриса, предложенную в статье [15]: da = C Gm G„ b 1-R2 (11) где l* - длина когезионного элемента e в направлении распространения расслоения; led - длина усталостной трещины в элементе e; dld/ dN - средняя скорость роста усталостных трещин; lFPZ - длина зоны процесса разрушения, равная расстоянию от фронта дефекта до точки, где напряжение сцепления достигает максимума [13]. Если предположить, что отношение повреждённой длины элемента led к его полной длине le соответствует отношению рассеянной в этом элементе энергии к вязкости разрушения, то можно прийти к выражению где C и b - постоянные материала, определяемые экспериментально по кинетической диаграмме усталостного разрушения; R = Pmin / Pmax - коэффициент асимметрии цикла; Pmin и Pmax - минимальное и максимальное значения циклической нагрузки; Gmax - значение интенсивности освобождения энергии, соответствующее нагрузке Pmax. Подставляя выражение (11) в уравнение (10), можно получить формулу для вычисления приращения параметра повреждения в рассматриваемом элементе e в зависимости от данного увеличения числа циклов: Ad* (AN) _ и un ( (u* )2 lF, С max z F] -C G* X G„ AN. (12) Описанная выше формулировка реализована в МКЭ-пакете ANSYS. Когезионную зону можно определять с использованием как интерфейсных, так и контактных элементов. При этом интерфейсными элементами можно соединять только элементы твёрдого тела, а контактными - ещё и оболочечные элементы. Второй подход представляется более предпочтительным, поскольку он даёт возможность рассматривать расслоения в реальных тонкостенных элементах конструкций сложной конфигурации. Таким образом, эффективный метод моделирования дефектов данного типа заключается в следующем. Сначала задаются две совпадающие поверхности, которые будут представлять отсчётные поверхности двух соединяемых оболочек. Если в качестве таковых выбираются нижние поверхности, то нормали к ним должны быть направлены наружу. Данные поверхности покрываются одинаковыми сетками оболочечных элементов. Затем между ними в месте расположения начального расслоения определяется стандартный контакт. В оставшейся части, представляющей область возможного распространения дефекта, задаётся прикреплённый (bonded) контакт, причём для b 251 Вестник СибГАУ. 2014. N 4(56) установления когезионной зоны здесь необходимо с контактными элементами связать CZM-материал. Алгоритм внесения усталостных повреждений в когезионные контактные элементы рассмотрим на примере I типа нагружения (нормальный отрыв). После выполнения очередного шага во всех узлах j каждого когезионного элемента e, расположенного вблизи фронта расслоения и ещё не полностью разрушенного (т. е. uemax < uc), определяются относительные перемещения uej. По ним рассчитывается среднее значение u = 4E uj , (13) 4 j характеризующее раскрытие дефекта в центре элемента. Если ue > uth (где uth - пороговое значение), то для рассматриваемого элемента e с использованием найденного на предыдущем шаге значения uemax определяется параметр повреждения de old. Далее вычисляется напряжение, после чего рассчитывается значение интенсивности освобождения энергии Gemax. Затем по формуле (12) определяется приращение Ad\ Это позволяет найти новое значение параметра повреждения: dneeW = dOld +Ade. (14) Следует отметить, что в настоящей работе для задания изменений параметров повреждения когезионных элементов предлагается использовать дополнительные относительные перемещения. Для элемента е в соответствии с формулой (1) данное перемещение вычисляется как Aue =-UcU0-т- ue, (15) uc - dneew (uc - u0 ) которое назначается каждому из его узлов: Au1 = Aue. (16) После определения перемещений для всех элементов, находящихся в зоне процесса разрушения, выполняется их осреднение в узлах: AUj = - EAU* , (17) где суммирование осуществляется по всем элементам из FPZ, сходящимся в узле j; n - общее число этих элементов. Таким образом, каждый шаг накопления усталостных повреждений здесь состоит из двух шагов нагружения. Сначала к верхней и нижней оболочкам в соответствующих узлах в дополнение к уже имеющимся перемещениям прикладываются нормальные перемещения ±Auj / 2 и для элементов e рассчитываются новые значения uemax, характеризующие накопленные повреждения. Затем производится удаление такой дополнительной нагрузки. Численные примеры. Проверка работоспособности предложенного способа моделирования роста расслоений в квазистатической постановке проводилась на образцах-балках DCB и ENF. Благодаря симметрии данных образцов здесь можно ограничиться рассмотрением лишь половины конструкции шириной B/2, задавая в плоскости симметрии соответствующие граничные условия. Во всех случаях моделирование слоёв осуществлялось при помощи оболочечных элементов SHELL181. Во всех примерах нормальная Kn и тангенциальная Kt контактные жёсткости взяты равными 1-105 Н/мм3. Для более точного определения несущей способности, а также для возможности исследования закри-тического поведения вместо приращений сил задавались приращения перемещений и вычислялись соответствующие реакции, характеризующие внешнюю нагрузку. Результаты расчётов, выполненных в программе ANSYS®, сравнивались с решениями плоских задач, представленными в работах [1; 6], где моделирование рассматриваемых образцов осуществлялось с использованием двухмерных элементов плоского напряжённого состояния. DCB-образец. DCB (Double Cantilever Beam) - образец в виде двухслойной консольной балки с начальным расслоением длиной а0 на свободном конце, нагруженный двумя противоположно направленными силами P, приложенными к верхнему и нижнему слоям (рис. 3). Он предназначен для экспериментального определении вязкости разрушения GIc для I типа деформирования. Используемые в расчётах свойства материала и размеры DCB-образца представлены в табл. 1 [6]. Предполагается, что слои балки изготовлены из однонаправленного материала, армированного волокнами, которые ориентированы вдоль её оси. 1 - а° » 1 Р,5 2h г L Р,5 Рис. 3. DCB-образец Таблица 1 Свойства материала и размеры DCB-образца [6] Ei, ГПа 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 252 Технологические процессы и материалы Таблица 2 Свойства материала и размеры ENF-образца [6] E, ГПа V GIIc, Н/мм т0, МПа L, мм а0, мм 2h, мм B, мм 70 0,33 1,45 57 100 30 3 10 Рис. 4. Зависимость «нагрузка-перемещение» для DCB-образца Р, 8 т 1 2h У і L °с А \ J 1 Рис. 5. ENF-образец Рис. 6. Зависимость «нагрузка-перемещение» для ENF-образца 253 Вестник СибГАУ. 2014. N 4(56) Результаты расчётов приведены на рис. 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 Упругие свойства углепластика Ei = 150 ГПа ^r o' II > G12 = 4,315 ГПа E2 = 8,819 ГПа ^r o' II r<1 > G13 = 4,315 ГПа E3 = 8,819 ГПа OO o' II r*i > G23 = 3,200 ГПа Помимо этих данных примем также, что коэффициент жёсткости когезионного материала K составляет 1105 Н/мм3, а пороговое значение раскрытия uth равно u0, где W 30 о 1П-4 u0 = -- =- = 3 -10 мм. 0 K 1-105 Следуя стратегии огибающей профиля нагружения, на первом этапе прикладывалась квазистатиче-ская нагрузка, линейно возрастающая от нуля до Pmax. Затем на втором этапе при постоянной нагрузке Pmax активизировался алгоритм внесения усталостных повреждений. На рис. 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.

作者简介

Sergey Chernyakin

Samara State Aerospace University

Email: chemyakin-sa@mail.ru
assistant of Space engineering department

Yury Skvortsov

Samara State Aerospace University

Email: proch@ssau.ru
Cand. Sc., Docent of Space engineering department

参考

补充文件