CC BY
14УДК 620.1.08
МОДЕЛИРОВАНИЕ ТРЕЩИНООБРАЗОВАНИЯ В ЖЕЛЕЗОБЕТОННЫХ КОНСТРУКЦИЯХ НА ОСНОВЕ ДЕФОРМАЦИОННОЙ РАСЧЕТНОЙ МОДЕЛИ
канд. техн. наук, доц. Д.О. ГЛУХОВ, д-р техн. наук, проф. Д.Н. ЛАЗОВСКИЙ, канд. техн. наук, доц. Е.Д. ЛАЗОВСКИЙ, Т.М. ГЛУХОВА (Полоцкий государственный университет)
Предлагается метод моделирования трещинообразования в железобетонных конструкциях, основанный на интеграции нелинейной деформационной модели железобетонного элемента по наклонному сечению, и, как частный случай, сечению нормальному к продольной оси элемента, и метода конечных элементов. В рамках такой интеграции появляется возможность учесть не только изменение жестко-стей конечных элементов под действием нагрузки и соответствующее перераспределение внутренних усилий, но и перераспределение напряжений по сечениям, изменение поля углов поворота главных напряжений в процессе распространения нормальных и наклонных трещин. В работе предложен метод комбинированных стратегий, позволяющий эффективно решать систему уравнений равновесия и совместности деформаций расчетной модели.
Ключевые слова: метод конечных элементов, пластические деформации, наклонные трещины, процесс трещинообразования, перераспределение усилий, изгибная и продольная жесткости.
Введение. Метод конечных элементов (МКЭ) устойчиво вошел в практику расчетов строительных конструкций. В настоящее время существует огромное количество программных продуктов конечно-элементного анализа, разработаны высокоточные стержневые, пластинчатые, объемные конечные элементы, конечные элементы, моделирующие пластические шарниры, элементы с меняющимся поперечным сечением и т.д. МКЭ имеет насыщенную историю, которая началась с работ Дж.У. Рэлея и В. Ритца и получила развитие в работах И.Г. Бубнова, Б.Г. Галеркина, Л.В. Контровича, О. Зенкевича, Р. Галлагера, И.Чанга, А.В. Перельмутера, В.И. Сливкера и многих других.
Для моделирования напряженно-деформированного состояния строительных конструкций, а также моделирования важного для железобетонных конструкций процесса трещинообразования существует несколько основных подходов. Прежде всего, это упрощенные инженерные методы поиска критических сечений [1]. Недостатком данного подхода является невозможность учета перераспределения усилий в процессах трещинообразования, необходимость пересмотра упрощенных формул при изменении схем армирования, работа с ограниченным рядом простых симметричных сечений. Преимуществом является простота расчетов и достаточная точность.
Также все большую популярность приобретают модели на основе сеток и объемных конечных элементов [2-4]. Преимущества данного подхода очевидны. Появляется возможность моделирования произвольной формы поперечного сечения и армирования, в том числе предварительно напряженной железобетонной конструкции, верификации и эффектной визуализации результатов расчета. Основными недостатками данного подхода являются высокая вычислительная сложность и ресурсоемкость, сложность перехода от конструкции сплошного поперечного сечения к конструкции, состоящей из отдельных блоков, разделенных трещинами, необходимость существенного сгущения сетки в области априорного существования трещин.
Основная часть. Предлагаемый метод основан на использовании МКЭ для стержневых систем (в том числе и для учета пластических свойств материалов) и нелинейных расчетных моделей сечений железобетонного элемента в сложном напряженно-деформированном состоянии при одновременном воздействие продольной силы, изгибающего момента и перерезывающей силы. Данный метод позволяет получить параметры напряженно-деформированного состояния железобетонного элемента: распределение сдвиговых деформаций, главных сжимающих и растягивающих деформаций и напряжений, поле углов поворота направлений главных сжимающих/растягивающих напряжений, визуализировать процесс перераспределения усилий в конструкции, связанный с изменением соответствующих жесткостей конечного элемента, а также влияние процесса перераспределения усилий на процесс трещинообразования.
Нелинейные расчетные модели железобетонного элемента по сечению, наклонному к его продольной оси, и алгоритмы решения систем интегральных сингулярных уравнений исследованы нами в течение последних 10 лет. Они основаны на модифицированной теории полей сжатия, разработанной канадскими учеными Майклом Коллинзом, Фрэнком Веккьо и Денисом Митчелом. MCFT (Modified Compression Field Theory) вошла в строительные нормы и правила Канады CSA-A23.3 1994 (CSA1994) [5-8].
Система уравнений исследуемой нами модели сложно-деформированного состояния сечения железобетонного элемента представлена в виде
My =-ÍMs(х' У))(У - Уо) dxdy+N (y- - Уо)
C
Mx =-lfCT(S( X' У))( X - Xo ) dxdy + N ( X- - Xo )
C
, (1)
N = ff a(s( x, у )) dxdy
C
V = ff t (x, y) dxdy
C
где o(s( x, y)) - нормальные напряжения в бетоне или арматуре, s( x, y) - продольные деформации бетона или арматуры (проекция на ось х суперпозиции главных сжимающих и главных растягивающих деформаций, полученная с учетом поворота элементарных площадок под действием перерезывающей силы);
xa и уа - расстояние от центра изгиба сечения, нейтральной точки, до осей x и y соответственно; хае и Уое - расстояние от места приложения продольного усилия до осей x и y соответственно; N - продольная сила от действия внешней нагрузки; Mx, My - изгибающие моменты, действующие в плоскостях xz, yz; t( x, y) - касательные напряжения в бетоне.
Уравнения совместности деформаций выглядят следующим образом:
ax = ac2 COS2 a+ ac1sin2 a + Px asx, a = ac2 sin2 a + acl cos2 a +
y c2 c1 ^^^ ^ ry^syi
t = (-ac2 + ac1)sin a- cos a;
= sc2 cos2 a + scl sin2 a,
(2)
S = ec2 sin2 a + scl cos2 a, (3)
У xy ¿ \ •
— = (-sc2 + scl )sin a - cos a,
где acx, acy, ac1, ac2, Txy - напряжения в бетоне;
asx, a - напряжения в продольной и поперечной арматуре;
е*, еу, Уху - относительные деформации элементарного участка;
6 - угол, определяющий положение наклонных трещин и главных напряжений.
Для определения распределения касательных напряжений, обеспечивающих равновесие по перерезывающей силе, нами использовался метод двух сечений [7]. Особенности данной математической модели и алгоритм решения системы уравнений равновесия и совместности деформаций подробно описаны в диссертации Е.Д. Лазовского [9] и в работе Т.М. Глуховой [10].
Алгоритм определения напряженно деформированного состояния элементов под воздействием среза с изгибом и продольным сжатием строится нами по следующей схеме:
1) задается начальное распределение касательных деформаций;
2) по методу двух сечений [7] определяется распределение касательные напряжений, обеспечивающее равновесие по перерезывающей силе;
3) для заданного изгибающего момента и усилия продольного сжатия/растяжения находится распределение продольных деформаций;
4) по полученным распределениям продольных деформаций разработанным нами методом комбинированных стратегий находится такое распределение деформаций, при котором главные сжимающие напряжения, определенные диаграммным методом, равны значениям, определенным по кругам Мора;
5) имея распределения всех составляющих тензора деформаций, определяются внутренние усилия, возникающие в сечении, нормальном к продольной оси элемента;
6) если определенные внутренние усилия не удовлетворяют уравнениям равновесия, то соответствующим образом корректируется искомое распределение деформаций и выполняется переход к пункту 2, иначе - успешное завершение итерационного процесса.
Пункт 4 алгоритма является ключевым для решения системы уравнений. Для понимания задачи, разрешаемой в пункте 4, приведем математическое описание зависимости главных сжимающих напряжений, определенных по диаграммам и кругу Мора, с заданным распределение продольных, главных и касательных деформаций.
Угол наклона осей главных напряжений по отношению к продольной оси определяется по формуле
-b,sT, + tan q =-
b • s • f
4ASwEs (sy - Sx, ) • (A^Es + ' is f ) + b2s2T,2
1+ J200 • eu
2 • ASWES (Sl, -Sx, ) Q -bi"i 1 >/4ASwES (S1,i - Sx,i ) • (ASwS1,iES + b • Sy • s • Еа ) + b,.2S\¡2
tan q =-5-, если s
2 • ASwES (S1,i - Sx,i )
Относительные деформации по направлению главных сжимающих напряжений
s2 ; = sx i (1 + tan2 q ) — S ¿ • tan2 q .
Относительные деформации по направлению оси y
S +S •tan2 q
S = —--
y1 + tan2 q
Напряжения в поперечной арматуре
а,, = s . • e .
Sy, y, i S
Главные растягивающие напряжения в бетоне
ctci,Í = ei,í • Ec, если ei,,- <s,rt ; fct
стс1 i =-1 , если e г < e .
i+7 200 -e 1,1
Касательные напряжения, способные передаться через наклонную трещину,
0,18 -f
t„- =
24 • w '
0,3 + -
«max + 16
Ограничиваем значение главных растягивающих напряжений в бетоне
= , если ^ ^ •tan 0 ;
= Tci ■tan 0 , если > Tci ■tan 0 .
Главные сжимающие напряжения в бетоне по измененной диаграмме деформирования бетона с учетом двухосного напряженного состояния
(аc2,i У =а
с 2 max,/
-0,002 -0,002
Главные сжимающие напряжения в бетоне из кругов Мора для напряжений
)'' = • (tan0 + -^).
' tan 9 '
В модели присутствует несколько точек разрыва поверхности решения для главных сжимающих напряжений. Во-первых, наличие области в окрестности нуля, для которой значение является невычислимым (точка разрыва второго рода). Во-вторых, характеристическим значением продольной деформации на растяжение является предельная деформация на растяжение scrt, определяющая скачкообразное изменение напряжения в результате трещинообразования (точка разрыва первого рода), выражающееся в замене математического определения угла наклона направления главных напряжений.
Проиллюстрируем данную ситуацию на семействе графиков (рисунок 1). Наглядно видно, что для определения границ области, в которой локализовано решение, необходимо применять различные стратегии поиска границ и решения как такового.
а - точки совпадения значении главных сжимающих напряжении, определенных диаграммным методом и геометрическим способом по кругу Мора; Ь - характер зависимости при больших значениях деформации, соответствующих главным сжимающим напряжениям; с - скачкообразное изменение напряжения при достижении деформацией значения ; й- аналогичный скачек для напряжении, определенных по кругу Мора
Рисунок 1. - Зависимость главных сжимающих напряжений от деформаций при заданных распределениях продольных и сдвиговых деформации для значении изгибающего момента до момента образования трещины
Для «больших» изгибающих моментов характер зависимости меняется. Точки решения «уходят» в область, значительно превышающую точку разрыва, что видно из рисунка 2.
а)
/ / В 10(19 С 1счг19
/ / У В 1ое4 С1ов4
/ с" ДО В1ов12
/ / ' / С1ов12 В1ов16
/ в 1ое19 С 1.1819
7У
| 1
Ь) 30
— В 1ор1
- С1о81
0.01 0.02 0.03 0.04 0.05
0.0» 10" 5.0*10"' 1.0» КГ* 13* Ю-4 2.0«! О""1 2.5*10"« 3.0» 10^
с)
/в\оц4 -/- С1ов4 / В 1ое8
- В1о612 - С1оВ12 В1о816
/ С (0816
/ / / / - С10819 - В1ов1 - С 1ов 1
/
с1)
В 1ое4 С1О84
1 --■—.
0.006 0.0065 0.007 Г) Й075 0.008
0.002 0.004 0.006 0СИ1К 0.01
Рисунок 2. - Зависимость главных сжимающих напряжений от деформаций при заданных распределениях продольных и сдвиговых деформаций для «больших» значений изгибающего момента
С учетом выявленных особенностей математической модели разработан алгоритм поиска решения для задачи определения распределения деформаций, соответствующих главным сжимающим напряжениям, комбинирующий различные стратегии поиска границ области решения и решения как такового:
1) методом последовательных приближений ищем нижнюю границу области определения функции главных сжимающих напряжений, определяемых из круга Мора;
2) если найденная нижняя граница меньше есг(, то априорная верхняя граница принимается равной
есг(, иначе принимается как максимально возможная сжимающая деформация (0,06).
3) методом дихотомии корректируем верхнюю границу области решения (область отрицательных значений) главных сжимающих напряжений, определяемых диаграммным методом;
4) методом секущих в пределах определенной области решения ищем точку совпадения искомых напряжений.
Данный алгоритм реализован нами на языке программирования С++. Сравнительные тесты с канадским программным комплексом Response-2000 [11] университета Торонто, показал, что предложенный нами алгоритм обладает более высокой скоростью сходимости и точностью расчета, дает более стабильные результаты расчета (отсутствие «выколотых» точек на графиках распределений) на всей области определения.
С целью реализации предложенной конечно-элементной модели в рамках исследования разработана библиотека функций по формированию, расчету и визуализации результатов расчета пространственных стержневых конструкций, построенных на сжато -изогнутых стержневых конечных элементах. Для плоского случая с 6 степенями свободы матрица жесткости и вектор реакций имеет следующий вид:
" ЕЛ Т 0 0 ЕЛ Т 0 0 2
0 12Ы l3 6Ы l2 0 12EJ 13 6EJ 12 ЧуЬ 2 щ
[к ], = 0 ел 6Ы 12 0 4£/ l 0 0 ЕЛ 6Ы 12 0 2EJ 1 0 , Ш, = ЧуЬ 12 дхЬ V 6, щ
l 1 2 V,
0 12Ш l3 6Ы l2 0 12EJ 13 6EJ 12 ЧуЬ 2 к.
0 6Ы 12 2Ы l 0 6Ы 12 4EJ 1 _ ЧуЬ2 . 12
где Е1 - изгибная жесткость сечения;
ЕА - продольная жесткость сечения;
Ь - длина конечного элемента;
и, V, 6 - горизонтальное, вертикальное перемещения и угол поворота узлов конечного элемента.
При формировании глобальной матрицы жесткости, локальные матрицы умножаются на матрицы трансформаций с целью приведения их к глобальной системе координат и соответствующим степеням свободы глобальной системы. Далее локальные матрицы объединяются по правилу суперпозиции по глобальным степеням свободы. К глобальному вектору {/} добавляются сосредоточенные узловые нагрузки и граничные условия. Полученную систему линейных алгебраических уравнений решают прямыми или итерационными методами.
Традиционно изгибная жесткость железобетонного элемента определяется как произведение изгибающего момента на радиус кривизны нейтральной оси элемента. Однако для усиленного под нагрузкой элемента мы оперируем понятием кривизны и момента до и после усиления. В таком случае изгибная и продольная жесткость определяется интегрированием:
Шу =-[[E(х,у)(У -У0)2dxdy;
C
EA = -Ц E ( x, у ) dxdy.
C
Основное уравнение метода конечных элементов для вычисления перемещений имеет вид
{8}=[ * Г },
где {8} - вектор перемещений;
- вектор внешней нагрузки;
[ К ] - глобальная матрица жесткости системы.
Вектор внутренних усилий отдельно взятого конечного элемента определяют как произведение локальной матрицы жесткости и вектора перемещений узлов конечного элемента.
Построенная конечно элементная модель подвергается адаптивной дискретизации: дальнейшему автоматическому разбиению с шагом, соответствующим среднему расстоянию между трещинами. Каждому конечному элементу назначается диаграмма зависимости изгибной жесткости от изгибающего момента и продольной силы в виде кусочно-линейной 3D-интерполяции по выделенным характеристическим точкам зависимости, полученной в модуле нелинейного моделирования.
В процессе моделирования нагружения организуется цикл последовательного увеличения нагрузок с фиксацией минимальной жесткости каждого сечения. Результаты моделирования выводятся в виде анимации расчетной схемы, анимации изменения эпюр моментов, перерезывающей силы и прогиба. На эпюре моментов цветом иллюстрируется изменение жесткости конечных элементов. Процесс анимации останавливается при достижении предельного момента, определенного в материале, в каком-либо конечном элементе. Отличия в эпюрах моментов, определенных по линейной и нелинейной модели, могут достигать 30%. Результат моделирования 3-пролетной неразрезной балки показан на рисунке 3.
а
б
Рисунок 3. - Результат расчета по линейной (а) и нелинейной (б) моделям. Отличие в определении максимального момента составило 16,98%
Кроме этого мы видим, что замена областей трещинообразования единственным пластическим шарниром, устанавливаемым в середине пролета или над опорами, не отражает реального процесса перераспределения напряжений и, соответственно, усилий. Данное обстоятельство имеет и экспериментальное подтверждение в целом ряде исследований, выполненных в ИИЦ Полоцкого государственного университета.
Определив для каждого конечного элемента внутренние усилия, по нелинейной деформационной модели мы имеем возможность вычислить распределение деформаций и напряжений по сечению железобетонного элемента и построить поле соответствующих параметров в конструкции в целом. Пример вычисленных полей приводится на рисунке 4.
Рисунок 4. - Поля угла наклона а направления главных растягивающих напряжений, нормальных (проекция на ось х) напряжений Ох, касательных напряжений т
Анализ выполненных исследований показал, что в большинстве литературных источников [9] ширина раскрытия наклонной трещины V (мм) определяется из относительного удлинения поперечной арматуры е^у, мм/мм, на определенном участке I, мм:
w _ е I.
(4)
Длину участка l, на котором определяется деформация поперечной арматуры, разные научные школы предлагают принимать по-разному. По ЕКБ-ФИБ - это расстояние между наклонными трещинами, мм:
^ __1_
m0 sin 0 cos 0 ,
где ^, Бт - шаг наклонных трещин соответственно в продольном, поперечном направлениях, мм:
S_ _ 2| C + + 0,25kj ^ ; Smy _2íCy +S] + 0,25k ^
Ps
где Сх, С - расстояние от центра тяжести сечения до центра тяжести соответственно продольной,
поперечной арматуры, мм;
6 - угол наклона главных напряжений по отношению к продольной оси элемента, рад; ^ - расстояние между стержнями продольной арматуры, мм; ^ - шаг поперечной арматуры, мм;
к - коэффициент, учитывающий форму поверхности арматурных стержней (профильная/гладкая); 0х, 0 - диаметр соответственно продольной, поперечной арматуры, мм;
- степень продольного армирования, мм2/мм2; р^ - степень поперечного армирования, мм2/мм2;
А^ - площадь сечения продольной арматуры, мм2; к - высота сечения элемента, мм; Ь - ширина сечения элемента, мм; А - площадь сечения поперечной арматуры, мм2;
р = —;
Г SX 77' hb
Psw =
Asw
bs
Более поздние публикации ЕКБ-ФИБ рекомендуют вместо относительного удлинения поперечной арматуры е принимать относительные деформации по направлению главных растягивающих напряжений в бетоне ес1, мм/мм. Таким образом, выражение (4) приобретает вид w = ес1 .
Однако обращаем внимание на то, что в области трещинообразования определить угол наклона главных напряжений по отношению к продольной оси элемента 9 не представляется возможным, поскольку напряжения в бетоне отсутствуют в связи с тем, что деформации в бетоне растянутой зоны превысили предельную деформацию на растяжение. Поэтому мы фиксируем угол поворота, предшествовавший моменту образования трещины.
С учетом изменения жесткостей конечных элементов под нагрузкой и перераспределения усилий, фиксации углов поворота главных напряжений, а следовательно, направлений образующихся трещин, нами реализован алгоритм пластического расчета трещинообразования в железобетонной конструкции.
Сравнение картин трещинообразования, полученных пластическим и упругим методами расчета, приводится на рисунке 5.
Рисунок 5.
б
■ Сравнение расчетных картин трещинообразования пластическим (а) и упругим (б) методами. Углы наклона главных напряжений показаны контур-линиями
Заключение.
1. В результате выполненного исследования построена система конечно-элементного анализа пространственных статически неопределимых строительных конструкций с учетом физической нелинейности, использующей реальные графики зависимости жесткости железобетонного элемента от вектора внутренних усилий.
2. Эксперименты показали, что моделирование физической нелинейности путем искусственного точечного позиционирования специальных конечных элементов, имитирующих пластический шарнир, некорректно, поскольку для оценки реального напряженно-деформированного состояния существенное снижение изгибной и продольной жесткости происходит по всей длине нагруженных конечных элементов и процесс трещинообразования идет на достаточно протяженных участках конструкции.
3. Исследована зависимость точности моделирования от степени дискретизации модели и количества шагов последовательного наращивания нагрузок.
4. Построены алгоритмы аппроксимации и визуализации поля деформаций и поля напряжений в железобетонном элементе на основе расчета серии сечений (нормальных и наклонных к продольной оси элемента).
ЛИТЕРАТУРА
1. Herbrand, Martin (2017). Shear strength models for reinforced and prestressed concrete members. DOI: 10.18154/RWTH-2017-06170.
2. Xenos, D., Modelling the failure of reinforced concrete with nonlocal and crack band approaches using the damage-plasticity model CDPM2 [Electronic resource] / D. Xenos, P. Grassl // Finite Elements in Analysis and Design. - 2016. -Vol. 117-118. - P. 11-20. - DOI: https://doi.org/10.1016/j.finel.2016.04.002.
3. Uncertainty of the smeared crack model applied to RC beams [Electronic resource] / A. Rimkus [et al.] // Engineering Fracture Mechanics (IF 3.426). - 2020. - Vol. 233. - DOI: https://doi.org/10.1016/j.engfracmech.2020.107088.
4. Sasmal, S. Nonlinear FE simulations of structural behavior parameters of reinforced concrete beam with epoxy-bonded FRP [Electronic resource] / S. Sasmal, S. Kalidoss // J. of the Mechanical Behavior of Materials. - 2015. - Vol. 24, № 1-2. - P. 35-46. - DOI: https://doi.org/10.1515/jmbm-2015-0004.
5. Collins, M.P. Pre-stressed Concrete Structures / M.P. Collins, D. Mitchell. - N.J. : Prentice-Hall Inc. : Englewood Cuffs, 1991. - 766 p.
6. General Shear Design Method / M.P. Collins [et al.] // ACI Struct. Journ. - 1996. - V. 93, № 1. - P. 36-45.
7. Vecchio, F.J. Analysis based on the Modified Compression Field Theory / F.J. Vecchio // IABSE Colloq. On Structural Concrete. - 1991. - V. 62. - P. 321-326.
8. Vecchio, F.J. The Modified Compression Field Theory for Reinforced Concrete Elements Subjected to Shear / F.J. Vecchio, M.P. Collins // ACI Journal. - 1986. - Vol. 83, № 2, - P. 219-231.
9. Лазовский, Е.Д. Усиление изгибаемых железобетонных элементов в зоне среза : автореф. дис. ... канд. техн. наук : 05.23.01 / Е.Д. Лазовский ; Бел. нац. техн. ун-т. - Минск, 2014. - 24 с.
10. Глухова, Т.М. Комбинирование стратегий поиска решений в задачах расчета систем нелинейных сингулярных уравнений /Вестн. Полоц. гос. ун-та. Сер. С, Фундам. науки. - 2013. - № 4. - С. 22-27.
11. Bentz, E.C. Sectional Analysis of Reinforced Concrete Members : PhD Thesis / E.C. Bentz ; Department of Civil Engineering, University of Toronto. - 2000. - 310 p.
Поступила 13.05.2021
MODELING OF CRACKING IN REINFORCED CONCRETE STRUCTURES BASED ON DEFORMATION CALCULATION MODEL
D. GLUKHOV, D. LAZOUSKI, Y. LAZOUSKI, T. GLUKHOVA
The article proposes a method for modeling crack formation in reinforced concrete structures, based on the integration of a nonlinear deformation model of a reinforced concrete element and the finite element method. Within the framework of such integration, it becomes possible not only to take into account the change in the stiffness of finite elements under the action of the load and the corresponding redistribution of internal forces, but also to take into account the redistribution of stresses over the sections, the change in the field ofrotation of the principal stresses in the process of crack propagation. The paper proposes a method of combined strategies that makes it possible to effectively solve the system of equations of equilibrium and compatibility of deformations of the design model.
Keywords: finite element method, plastic deformations, inclined cracks, crack formation process, redistribution of forces, bending and longitudinal stiffness.
