CC BY
19
ИРКУТСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ ПУТЕЙ СООБЩЕНИЯ ЙН|
Быкова Н.М., Каргапольцев С.К., Пыхалов А.А., Милов А.Е.
УДК 693.546.4
НЕКОТОРЫЕ ПРИНЦИПЫ РАЗРАБОТКИ КОНЕЧНО - ЭЛЕМЕНТНЫХ МОДЕЛЕЙ ТОННЕЛЕЙ В СТРУКТУРНО НЕОДНОРОДНЫХ ГОРНЫХ МАССИВАХ
СНиП 32-04-97 «Тоннели железнодорожные и автодорожные» [1] регламентирует необходимость соответствия расчетных моделей тоннельных обделок условиям работы сооружений, технологии их возведения, характеру взаимодействия элементов конструкций между собой и окружающим грунтом с учетом неблагоприятных сочетаний нагрузок и воздействий во время строительства и эксплуатации сооружения. Качество и надежность проектных решений при этом зависят от степени приближения построения расчетной модели к реальным условиям работы конструкций.
Работа подземного сооружения во время эксплуатации во многом определяется правильностью учета при проектировании внешних воздействий, статики и динамики взаимодействия конструкций, использования приближенных к реальным физических моделей работы материалов. В сухих однородных скальных грунтах транспортные тоннели могут нормально эксплуатироваться сотни лет. Совсем иная картина наблюдается в тоннелях, построенных в условиях изменения горного, гидростатического давления, температурных, сейсмических и геодеформационных воздействий. Под геодеформационными воздействиями понимаются смещения горных блоков по разломам. Сами по себе изменения внешних условий дополнительно влекут изменения физического состояния, как грунтов, окружающих тоннель, так и бетонов обделки.
Для обеспечения безопасности и долговечности эксплуатации транспортных тоннелей, развития эффективных методов содержания и ремонта конструкций подземных сооружений необходимо иметь достоверную оценку их работы и состояния. Контроль напряженно-деформированного состояния (НДС) подземных сооружений, прогнозирование их
поведения в зависимости от изменения каких-либо условий можно осуществлять с помощью виртуальной физико - математической модели, как глобальной, так и фрагментарной. Глобальная модель охватывает работу объекта в целом и может быть выполнена без детализации отдельных конструкций. Из глобальной модели вырезаются любые сложные участки с заданием граничных условий с целью анализа НДС более детальных конструктивных моделей, например, совместной работы конструкций обделки и верхнего строения пути. Математическое моделирование должно сопровождаться натурными методами сбора информации.
На первом этапе поставлена цель построения и анализа глобальной физико-математической конечно-элементной модели. В качестве объекта был выбран Северо-Муйский железнодорожный тоннель, отличающийся уникальностью своей работы в геодинамичес-ки-активном районе Байкальской рифтовой зоны.
Для достижения цели в работе решалось две задачи. Первая из них состояла в получении основных зависимостей МКЭ и построение конечно-элементной (КЭ) глобальной модели рассматриваемого объекта, куда входят: ее геометрические характеристики; процедура дискретизации на конечные элементы; применяемые материалы и их свойства; граничные условия выделения симметричной части тоннеля и сегмента земной коры; величины и особенности используемых типов силового воздействия. Моделирование горного массива осуществляется относительно данных геодезических съемок инженерно-геологических разрезов.
Вторая задача состояла в проведении анализа полученной КЭ модели тоннеля, пути и
горного массива, состоявшего из решения ряда физических задач, в качестве которых приведены расчеты напряженно - деформированного состояния (НДС) при воздействии гравитационной нагрузки; определение температурного поля в условиях стационарной теплопроводности, для максимальных условий зимнего периода эксплуатации; НДС при воздействии полученного поля температур.
Результатом представленных типов анализа является оценка НДС конструкции, представляющую собой полную протяженность обделки тоннеля пути и горного массива.
Основные зависимости МКЭ [2] построены для объемного НДС, с использованием конечного элемента в виде изопараметричес-кого гексаэдра первого порядка аппроксимации, с дополнительным подключением алгоритма несвязных функций формы [3]. Необходимость этого алгоритма обусловлена повышением уровня точности решения, получаемого без увеличения размерности глобальной системы алгебраических уравнений, что является решающим фактором для анализа сложных объектов моделирования.
Однако сама процедура решения задачи может оказаться эффективнее, если интерполировать перемещения полиномами более высокой степени, чем координаты. Удобным способом улучшения изопараметрического элемента является введение несовместных функций формы.
В этом случае поле перемещений деформируемой системы представляется в виде:
(1)
где: и1 — и8, у1 — у8 , - - перемещения, связанные с 8 узлами конечного элемента сетки; и9 — и11, - у11, - ш11 -дополнительные степени свободы, определяющие более высокую степень интерполяции перемещений.
Вектор-столбец функций формы в этом случае имеет вид:
- —л Т
и 1 >1 ^
и = {ТУ }< и 2 II > > 2
и11 > 11 .
{#} м2 ... мп},
(2)
вышенные деформации сдвига, возникающие при моделировании изгиба толстых пластин и оболочек объемными изопараметрическими конечными элементами. В рассматриваемом объекте моделирования к ним можно отнести обделки тоннеля и другие конструктивные элементы. Явление завышения деформации сдвига возникает из-за неспособности узловых функций форм конечного элемента отразить деформированное состояние, вызванное изгибающим воздействием в моделируемом объекте.
На рис. 1а показаны деформации бесконечно малого объема, вырезанного, например, из обделки тоннеля, в сравнении с деформациями конечного элемента, того же объема (рис. 1б), находящихся под действием одного и того же изгибающего момента. Из анализа рисунков видно, что конечный элемент с линейной аппроксимацией не способен отразить деформации изгиба и заменяет их деформациями сдвига. Завышением деформаций сдвига "отбирается" часть энергии деформации, что приводит, в частности, к равновесию работ деформирования твердых тел при меньших величинах перемещений узлов.
Тем самым КЭ модель ужесточается по сравнению с реальной конструкцией.
Выражение для вычисления матрицы жесткости конечного элемента [2], с использованием несовместных функций формы [3] за-
писывается в виде:
И - М \р №.
V(е)
Или, в квадратурной форме записи
[к ](-I[В,,С., V)]т[в(т.)]>
г =1 }=1 т =1
(4)
(5)
[ В (§, , Л ,, Ст , Гт )] [ } (5, , Л ,, Ст )],
а)
где: Ng, N11 - несовместные функции формы, задаваемые в виде:
= 1 2; = 1 -Л2; = 1 -<;2. ( 3) Использование несовместных функций формы при составлении матрицы жесткости конечного элемента позволяет уменьшить за-
б)
Рис. 1. Явление завышения деформаций сдвига при моделировании изгиба толстых пластин.
ИРКУТСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ ПУТЕЙ СООБЩЕНИЯ
где: Т.т - значение температуры в текущей
точке интегрирования:
Т.. =
ут
{Ж£ . , Л. ,5т )}
Т
(6)
дЖ„ дЖ,„ дЖ,
дг дг дг
дЖ 9 дЖ 10 дЖ 11
д0 д0 д0
дЖ 9 дЖ 10 дЖ 11
дг
дг
дг
]
0 0
0 "2Л 0
0 0
"2?
(8)
Рассматривая
Выражение для матрицы дифференцирования, с использованием несовместных функций формы, приобретает вид:
[в] = [в в2 ... вп]. (7)
Для вычисления производных от несовместных функций формы N3, N по глобальным координатам х, у, г используется прием аналогичный узловым составляющим, где имеем:
где {5} е - содержит узловые значения перемещений, 5 - содержит переменные, не связанные с узлами КЭ,
[К](е) =[КККТК] . (10)
Из второй строки в выражении (9) находим
5 = "К 1К Т {5}(е). (11)
После подстановки этого выражения в первую строку выражения (7), получаем известное соотношение:
[К ](е) {5} ={^}(е), (12)
где [К](е) - редуцированная матрица жесткости конечного элемента порядка 24, определяемая выражением вида:
[К](е) = К " КК " К Т. (13)
Нередуцированная матрица жесткости элемента может быть записана в виде
[КГ =
К11 К12 ...К111
К К К
К 21 К 22 ...К 2 11
К К К К 11 1 К 112...К 11 11
(14)
где каждая из подматриц имеет размерность [3x3] и строится следующим образом
8
к„ = 1
¿=1 .=1
т=1
Кц к12 к1з
1 Г£ 1 Г£ 1 Г£
К 21 К 22 К 23
К" К™ К™
К 31 К 32 К 33
det [ ^., Л., С т )], (15)
матрицу жесткости элемента, необходимо заметить, что использование выражения (4) приводит к матрице [К] )
33-го порядка, а степени свободы и9,и10,..жп не связаны с узлами конечного элемента. Для снижения размерности матрицы до 24-х, дополнительные степени свободы исключаются путем статической конденсации, которая эквивалентна минимизации полной потенциальной энергии конечного элемента, как деформируемой системы, относительно несвязных перемещений и9,и10,..жп. Рассмотрим уравнение равновесия:
[К](е) г_|(е) !=*И(е) к (9)
к™ к™ к15
к11 к12 к13
к15 к15 к15
21 22 23
к15 к15 к15
31 32 33
(16)
Бг
г 'Л . ■ Ст )Г [ о (. )]в (¡^ , Л ., Ст )].
После подстановки в выражение (13) матрицы дифференцирования, содержащей производные от узловых и несовместных функций формы, а также выражения для матрицы упругости объемного НДС, элементы подматрицы Кге записываются в явном виде:
к? =(б;б; + в2в 2 )о +(б[б2 + в: в 2 )о2 + +(вг3 в 3 + в4 в 4 )ю3-, =(в1 о + в2 о1)в 3 + в3 в о к;; =(в; + в 2 )в 2 + в;в ;в 3;
К21 =(в;в 2 + в в )в; + в ;в; в ъ; (17)
к 22 = (в ;в; + в ;в 5 )в; + в ;в в; к 23 = в\в\в 2 + в ;в;в;; К 31 =(в; + в 2 )в; в 2 + в;в ;в;;
К32 = в в в 2 + в ;3 в;в з;
к ;з =(в;в; + в ;в; )в; + в 4 в 4в.
Поскольку матрица жесткости конечного элемента [К](е) симметрична, естественно, что работать и хранить ее в памяти вычислитель-
Рис. 2. Поперечное сечение тоннеля.
ной машины выгодно в виде ее верхней (нижней) части матрицы, включая главную диагональ.
Таким образом, столбцы верхнего треугольника матрицы с диагональными элементами можно последовательно поместить в век-
тор-столбец {K}
( e)
{к}( е 1 =
K K K K
K
10
K 5 K 9
K
.300
K
K K
K
277
. (18)
Для доступа к г-му, j-му элементу используется следующая FORTRAN конструкция:
Value = Ke(j*(j + 1)/2+1-i), где: Value - i-ый, j-ый элемент матрицы [K](e).
Характеристика тоннеля и пути.
Анализируемый тоннель, а также железнодорожный путь и горный массив являются протяженным объектом. Геометрическая характеристика поперечного сечения обделки тоннеля, используемого для моделирования, представлена на рис. 2. Обделка тоннеля, пути и горного массива моделируются в полную длину 15,3 км, с чередованием по длине горных блоков различной направленности, заполненных дробленным перетертым до пес-чанно-глинистого состояния грунтом №1 и крепкими гранитами — грунт №2, 3. Модель допускает задание каждому блоку своих физико - мехимических характеристик.
Исходная информация, необходимая для моделирования геометрических характеристики горного массива, а также расположение
Рис. 4. Общий вид конечно-элементной (КЭ) модели горного массива, обделки тоннеля и пути.
в нем слоев грунта, представлена на рис. 3 и построена в системе AutoCAD.
Для выявления общих закономерностей изменения НДС тоннельной обделки используется упругая модель Мизеса - Губера. Применение представленной модели допустимо, поскольку рассматриваемые конструкции обделки тоннеля и пути преимущественно находятся в условиях сжатия. Напряженно-деформированное состояние в горном массиве используется, в данном случае, в качестве гравитационной нагрузки на представленные конструкции.
Разбиение геометрической модели тоннеля и горного массива на конечные элементы осуществлялось с необходимым сгущением конечно-элементной (КЭ) сетки наиболее на-
Рис. 5. Фрагмент модели с плавным переходом от КЭ сетки обделки тоннеля к КЭ сетке горного массива.
груженных мест деформируемого объекта. Общий вид КЭ модели горного массива, обделки тоннеля и пути представлен на рис. 4.
Дискретизация горного массива в КЭ модели была выполнена более крупными конечными элементами. Сетка конечных элементов подстраивалась под форму и размеры слоев грунта, которые в модели (рис. 4, 5) представлены цветом, соответствующим их свойствам. Для моделирования связи между КЭ сеткой горного массива и КЭ сеткой обделки тоннеля используется плавный переход, выверенный относительно градиента изменения поля напряжений.
Дискретизация обделки тоннеля, а также, области, непосредственно прилегающие к ней, имеют максимальную плотность разбиения. По длине обделки тоннеля в сетке используется один конечный элемент на три погонных метра длины тоннеля. Другие размеры конечных элементов выдержаны относительно требований МКЭ по пропорции их сторон и углов между ребрами. Основной подбор необходимого сгущения КЭ сетки проводился в области обделки тоннеля, его сопряжения с путевым бетоном и рельсом. На рис. 6 представлено моделируемое сечение и сопряжение путевого бетона и рельса.
Дискретизация рельс осуществляется с применением балочных конечных элементов. Соединение балочных КЭ с узлами КЭ сетки путевого бетона осуществлялось посредством конечных элементов типа жесткая вставка. Соединение осуществляется между центром масс сечения балочного КЭ и узлом КЭ сетки
МЕХАНИКА. ТРАНСПОРТ. МАШИНОСТРОЕНИЕ
О
Рис. 6. - Фрагменты сечения тоннеля.
путевого бетона. Моделирование шпал в представленной КЭ сетке проводится с учетом усреднения свойств путевого бетона и бетона шпал, моделируемых с применением конечных элементов объемного напряженно-деформируемого состояния. Позднее, при КЭ -моделировании отдельных участков с решением задач анализа НДС верхнего строения пути, планируется большая конструктивная детализация, чем в глобальной модели.
Всего в КЭ модели представленной конструкции горного массива, обделки тоннеля и пути использовано более 2 000 000 узлов КЭ сетки.
Граничными условиями выделяются в плоскостях симметрии степени свободы с взаимным нулевым перемещением на гранях симметричных частей конструкции горного массива, обделки тоннеля и пути. Это - защемление нормальных степеней свободы в горизонтальном направлении оси Ъ -поперек тоннеля, по оси X - вдоль тоннеля и вертикальном направлении по оси У - нижней плоскости горного массива. Представленные типы граничных условий используется для всех узлов КЭ -сетки рассматриваемой модели.
Свойства материалов в модели для первого приближения приняты изотропными. Значения свойств материалов пересчитаны в соответствии с заданными геометрическими единицами модели в миллиметрах.
Рельсы приняты марки Р-65 со следующими свойствами: Модуль Юнга Е =210000 МПа,
плотность р = 7,8*10-9т/мм3, коэффициент Пуассона л = 0,3, коэффициент теплопроводности X = 12*103 Вт/мм*К°, коэффициент линейного температурного расширения а = 30*10-61/С°.
Материал бетона обделки - В30, со следующими свойствами: модуль Юнга Е = 0,324*105 МПа, плотность р = 2,5*10-9 т/мм3, Коэффициент Пуассона л = 0,2, коэффициент теплопроводности Х=1.5*103 Вт/мм*К°, коэффициент линейного температурного расширения а = 12*10-61/С°.
Материал путевого бетона В15 со следующими свойствами: Модуль Юнга Е = 0,27*105 МПа, плотность р = 2,5*10-9 т/мм3, Коэффициент Пуассона л = 0,2, коэффициент теплопроводности Х=1.5*103 Вт/мм*К°, коэффициент линейного температурного расширения а=12*10-6 1/С°.
Материалы горного массива были условно разделены на три группы согласно коэффициенту крепости.
Материал "крепкого" скального грунта (синий цвет на рис. 4, 5) со свойствами: Модуль Юнга Е = 0,46*105 МПа, плотность р = 2,68*10-9 т/мм3, Коэффициент Пуассона Л = 0,2, коэффициент теплопроводности Х = 3*103 Вт/мм*К°, коэффициент линейного температурного расширения а = 8.1*10-61/С°.
Материал "слабого" грунта (красный цвет на рис. 4,5) со свойствами: Модуль Юнга Е = 0,13*105МПа, плотность р = 2,71*10-9т/мм3, Коэффициент Пуассона л = 0,3, коэффициент теплопроводности Х = 3*103 Вт/мм*К°, коэффициент линейного температурного расширения а =8.1*10-6 1/С°.
Материал "сыпучего" очень слабого грунта (темнозеленый цвет на рис. 4,5) со свойствами: Модуль Юнга Е = 0,011*105 МПа, плотность р = 2,65*10-9 т/мм3, Коэффициент Пуассона л = 0,3, коэффициент теплопроводности Х = 3*103 Вт/мм*К°, коэффициент линейного температурного расширения а = 8.1*10-61/С°.
Нагрузки включают в себя гравитационные силы горного массива, обделки тоннеля, путевого бетона и рельса пути. Эти силы определяется в КЭ модели относительно ускорения свободного падения, плотности используемого материала, а также относительно геометрических размеров моделируемых ко-
Рис. 7. Геологический разрез и перемещения в масштабе (мм) при воздействии гравитационных сил а)-геологический разрез; б) суммарные перемещения в модели горной перемычки; б) то же в тоннеле.
нструкций. Дополнительно рассматриваются различные варианты температурных и геодеформационных воздействий.
Расчет модели выполняется 17-20 часов. На рис. 7 показаны инженерно-геологический разрез, суммарные перемещения в мм всей модели и в точках верха и низа поперечного сечения тоннеля. Данный расчет выполнялся на приложение всей гравитационной нагрузки и служит иллюстрацией работы модели, а также на качественном уровне показывает неоднородный характер работы тоннеля по длине. Для определения количественных показателей НДС необходимы дополнительные приемы, учитывающие специфику передачи гор-
ного давления на обделку при различных строительных технологиях.
БИБЛИОГРАФИЯ
1. СНиП 32-04-97. Тоннели железнодорожные и автодорожные /Госстрой России. -М.: ГП ЦПП, 1997. - 19 с.
2. Зенкевич О.С. Метод конечных элементов в технике. - М.: Мир. -1975. -542 с.
3. Пыхалов А. А. Контактная задача статического и динамического анализа сборных роторов турбомашин/Автореферат диссертации на соискание ученой степени доктора технических наук. -М.:МАИ. - 2006. -32 с.
