CC BY
30
УДК 624.131
ЛАНИС АЛЕКСЕЙ ЛЕОНИДОВИЧ, канд. техн. наук, доцент, alangeo@bk. ru
Сибирский государственный университет путей сообщения,
620023, г. Новосибирск, ул. Д. Ковальчук, 191,
ХАН ГИЛ НАМ, канд. физ.-мат. наук, ст. научный сотрудник,
khan. igd@gmail. com
Институт горного дела СО РАН,
630091, г. Новосибирск, Красный проспект, 54
МОДИФИКАЦИЯ МОДЕЛИ ГЕОСРЕДЫ ДЛЯ РЕШЕНИЯ ЗАДАЧ МЕХАНИКИ ГРУНТОВ МЕТОДОМ ДИСКРЕТНЫХ ЭЛЕМЕНТОВ
В рамках метода дискретных элементов разработана математическая модель грунта, учитывающая его состав и особенности взаимодействия между его структурными элементами на их контактах. Показана перспектива использования метода дискретных элементов для решения задач геомеханики.
Ключевые слова: инъектирование; усиление; геосреда; метод дискретных элементов; уравнения движения; напряженно-деформированное состояние; контактные силы; механика грунтов.
LANIS, ALEKSEY LEONIDOVICH, Ph.D., Assoc. Prof., alangeo@bk. ru
Siberian State Transport University,
191 D. Kovalchuk st., Novosibirsk, 630049, Russia HAN, GIL NAM, Ph.D., khan. igd@gmail. com
Institute of Mining. Siberian Branch of the Russian Academy of Sciences,
54 Krasny Prospect st., Novosibirsk, 630091, Russia
MODIFICATION OF GEOENVIRONMENT MODEL FOR SOLVING THE PROBLEM OF SOIL MECHANICS BY DISCRETE ELEMENT METHOD
Within the discrete element method the mathematical model of soil, taking into account its composition, and features of interaction between soil’s structural elements has been developed. The prospect for the use of discrete element method to solve the problems of geomechanics has been shown.
Key words: Injection; strengthening; geoenvironment; discrete element method; the equations of motion; the stress-strain state; the contact forces; soil mechanics.
Повышение экономической эффективности использования объектов транспортной инфраструктуры неразрывно связано со снижением эксплуатационных расходов на их содержание и оптимизацией расходов на капитальный ремонт и реконструкцию таких объектов. Последний тип расходов нередко является определяющим в общих расходах. Экономически оптимальными явля-
ются технологии, позволяющие выполнение работ по капитальному ремонту, реконструкции без ограничения движения транспорта. К таким технологиям можно отнести методы, основанные на инъектировании твердеющих растворов
[1], распространение которых связано с неоднородностью грунтов, принимаемой по результатам инженерно-геологического обследования. При применении указанных технологий необходимо минимизировать объемы строительномонтажных работ без ухудшения несущей способности объектов транспортной инфраструктуры. Для этого на стадии проектирования надо иметь возможность адекватно оценивать эволюцию напряженно-деформированного состояния (НДС) геосреды, протекающую при ведении работ по реконструкции объекта. В частности, для земляного полотна путей сообщения необходимо дальнейшее совершенствование расчетно-теоретической базы, позволяющей прогнозировать изменение свойств земляного массива и его устойчивость при усилении грунтов современными технологиями.
На сегодняшний день существует множество теоретических описаний отдельных аспектов поведения грунтов, основанных на механике сплошных сред. Но общей теории, способной предсказать изменение их поведения, всё ещё нет. В такой ситуации компьютерное моделирование является практически единственным инструментом, с помощью которого возможно изучение поведения грунтов во всём его многообразии. Среди известных компьютерных программ, используемых для расчета НДС геосред, большой популярностью у исследователей пользуются пакеты программ PLAXIS и ANSYS, основанные на методе конечных элементов (МКЭ). Но смоделировать с их помощью такие процессы, как большие деформации и их разрывы, возникающие, например, при потере устойчивости откосами грунтового полотна, проблематично. Связано это с тем, что МКЭ основан на представлении исследуемой среды как сплошной. В рамках сплошной среды для описания указанных процессов, помимо определяющих уравнений, необходимо еще и использование специальных теорий неупругого деформирования материала и критериев разрушения, а также дополнительных уравнений, учитывающих дилатансионные проявления, возникающие при деформировании сред со структурой, к каковым относится грунт.
Полагаем, что для решения указанных задач большие перспективы могут быть связаны с использованием численного метода частиц [2], который в последние годы широко применяется в различных областях химии и физики, однако относительно мало используется в механике. Метод дискретных элементов (МДЭ) [3], предназначенный для решения задач механики сыпучих сред, грунтов, горных пород и механики деформируемого твердого тела, является одним из вариантов метода частиц. Как принципиально дискретный метод, МДЭ не имеет недостатков континуальных моделей, проявляющихся при нарушении сплошности среды или в результате того, что среда обладает внутренней структурой. В МДЭ дилатансионные свойства грунтов, а также развитие деформаций в зонах их разрушения реализуются автоматически в результате интегрирования уравнений движения отдельных частиц. При этом МДЭ позволяет отслеживать изменения НДС среды во времени. Данный метод основан на уравнениях движения отдельных элементов, из которых со-
стоит среда. Взаимодействуют элементы посредством сил, возникающих на их контактах. Перемещения элементов возникают под действием контактных и гравитационных сил, а также статических и динамических нагрузок, приложенных к отдельным элементам. В данной работе применительно к грунтам предлагается модификация модели горных пород, приведенной в [4].
В механике грунтов под потерей их прочности подразумевается момент начала проскальзываний одной части среды относительно другой. Прочность при этом определяется сопротивлением сдвигу, которое представляет сумму контактных сопротивлений ему, отнесенную к единице площади поверхности разрушений. Общепринятым условием прочности грунтов является условие Кулона - Мора
^n = °n^Ф + C, (1)
где Tn, <5n - касательные и нормальные напряжения на площадках с нормалью n (везде далее жирными символами обозначаются векторы); ф - угол внутреннего трения; C - сцепление.
В соответствии с (1), одной из основных характеристик, определяющих прочность грунтов, является сцепление C. Тип взаимодействия, реализующийся между частицами грунта вследствие действия сцепления, может быть невосстанавливаемым и восстанавливаемым. В первом случае оно возникает в результате межмолекулярной связи между частицами и высохшим межпо-ровым веществом. Такой тип сцепления после разрушения связи не восстанавливается. Во втором случае сцепление между частицами грунта, составляющими его твердую фазу, возникает благодаря присутствию в порах жидкой фазы. Ее заметное влияние на прочность грунта связано с тем, что его твердая фаза в основном состоит из относительно мелких, вплоть до дисперсных, частиц. В результате при удалении частиц друг от друга на некоторое критическое расстояние происходит разрыв такого типа сцепления, но при повторном их сближении оно между ними вновь восстанавливается благодаря силам поверхностного натяжения. Таким образом, в грунтах существенным может оказаться как невосстанавливаемый тип сцепления Сс, так и восстанавливаемый Са (назовем их когезией и адгезией соответственно), т. е. C = C0+Ca.
С учетом последнего соотношения для определения контактных сил рассмотрим модель грунта, состоящую из неразрушаемых частиц шаровой формы. На рис. 1 приведены фрагмент реального грунта, сложенного частицами неправильной формы (рис. 1, а), и моделирующая его система из шаров, равновеликих соответствующим частицам грунта (рис. 1, б). Полагаем, что между ними при взаимодействии могут возникать силы отталкивания Fnij, трения Fxij, невосстанавливаемого Fcij и восстанавливаемого Faij сцепления.
Каждый i-й элемент определяется радиусом ri и следующими физикомеханическими параметрами: плотностью рг-, углом контактного трения фг-, модулем упругости Еи когезией Cc, адгезией Са, модулями упругости Еа и Еы виртуальных пружин когезии и адгезии, вводимых для удобства интерпретации действия сил Fcij и Faij.
а б
Рис. 1. Схема дискретной модели:
а - фрагмент грунта; б - дискретная модель
Полагаем, что
Е. << Е., Е. << Е.. (2)
с/ i5 ai i V '
Рассмотрим контакт /-го и /-го шаров с радиусами г, Г/ и координатами X X] центра масс. В общем случае, как уже отмечено, между элементами возникают следующие типы контактных сил: отталкивания, трения, невосстанав-ливаемого и восстанавливаемого сцепления (рис. 2).
а 6 е-
Рис. 2. Схема действия контактных сил
Так как силы на контактах всегда представлены противоположно направленными парами одинаковой величины, достаточно рассмотреть передачу усилия от /-го элемента на 7-й элемент. Для обеспечения правильной ориентации векторов определим для 7-го элемента нормаль контакта п7/ как единичный вектор ортогональный плоскости Пу контакта, направленный от центра 7-го элемента
пу = / Ы, (3)
где 1^1 = |х/ - х7| = ёу - расстояние между центрами шаров, а плоскость Доопределяется расстоянием 1у до нее от центра 7-го элемента
=(4+ г> - г/)' (Ч )• <4>
Для отыскания силы Fnij отталкивания введем нормальное сжатие 5. шаров как величину их взаимного перекрытия (рис. 2, а):
5р = Г + Г - ^. (5)
При этом 5. можно интерпретировать как величину сжатия виртуальной «пружины» отталкивания, действующей между элементами при их частичном наложении друг на друга. Сила Fщ отталкивания всегда направлена против
нормали п.. Модуль силы Fnij = ^п/| находится из соотношения
Г5Ь, 5 .. > 0,
Е = I 4444 (6)
п |0, 5 .. < 0,
I ’ / ’
где k, а, Ь - коэффициенты; Е . = 2Е {Е}- / (Е + Е. ) ; г. = / (г , г. ). В частности,
при k = 1/3, а = 1/2, Ь = 3/2, г. = 2гг. / (г + г. ) получим силу отталкивания между шарами по Герцу, а при k = 1, а = 1, Ь = 1, г. = 2шт ( г/. ) - линейную зависимость Епц от 5..
Сила Fay адгезии ориентирована по нормали п. (рис. 2, б, в). Она возникает между элементами в случае, когда 5. >-5., где 5. = г.-Са.. / Еа.,
Са . = Ш1П (Са., Са. ), Еа. = 2Еа ,Еа, / (Еа, + Еа. ). Считаем, что при 5. < -5]. ЗДге-
зионная связь между элементами прекращается. В данном случае 5^ можно интерпретировать как величину растяжения виртуальной пружины адгезии, возникающей между взаимодействующими элементами. Сила Faij по модулю равна
|Е -г- (5*. +5 ..), 5 .. >-5*.,
г I аУ У \ У УР У У’ ,.-,4
*ад = | * (7)
1о, 5 .. <-5...
I ’ У У
Для определения силы Fсi] когезии в начальный момент времени 70 площадку контакта Fсij стянем в точку Р. = Р., расположенную на пересечении dij и плоскости П. контакта шаров. Точку Р. жестко свяжем с . -м элементом, а точку Р.. - с .-м элементом. В процессе поступательных и вращательных движений шаров эти точки расходятся (рис. 2 а, б), в результате на -й элемент действует сила Fсу. Введем виртуальную пружину? сцепления, растянутую в направлении е. на величину в. = |е|, где е. = Р.Рц (рис. 2, б). Для .-го элемента сила Fсij действует от точки Р. к точке Р. и по модулю равна
„ .,в. <в*,
ар I п ^ * ( )
I0, в. <в.,
где s^rAj /Ecj; Ccj =min(Q,Cg); ^ =2EC1ECJ /(Ecl + ECJ). Прив,, >в*.
происходит разрушение когезионной составляющей сцепления. При повторном перекрытии элементов она не восстанавливается.
Найдем силу Fxij трения (рис. 2, в). В начальный момент времени t0 для перекрывающихся шаров стянем площадку контакта FXy в точку Р. = Р., которая расположена на пересечении d. и плоскости П. контакта шаров. Точку Р. жестко свяжем с i-м элементом, а точку Р.г - с .-м. В процессе поступательных и вращательных движений элементов эти точки разойдутся. В результате между элементами возникнет сила трения Fxij. Примем ее равной проекции вектора Ftij, направленного по е. = P.P.i, на касательный к контакту вектор т.,
который ортогонален нормали n. и лежит в плоскости Ftlj и п.. По модулю сила трения равна
tgm F ..5 .. / 5* .., 5 .. <5* .., M.., M . е ^..,
от, т. хп х v ’ хи х. ’ .i v’
Fхv■ =
. т. х. х. ’ х. х. ’ . ’ j г у'
К , 5х. >5х., M ., M. е^., (9)
0, M .. й¥ й или M„ й¥ я,
где 5хv. = е. cos а. (5. - угол наклона е. к т.); M., M.. - проекции точек P.,
P. на плоскость контакта трения; кф - коэффициент, зависящий от выбора типа трения и выбираемый из диапазона 0 * М1 (при кф = 0 получим трение
типа stick-sleep); tg9. =mln(tgфl.,tgфJ.). При выполнении условия 5х. >5х.
точки Р. и Р., ранее жестко зафиксированные на соответствующих шарах, отделяются от них и занимают такие новые положения, что удовлетворяют следующим условиям: 1) направление вектора е. сохраняется; 2) величина вектора е. уменьшается до в. = к5х. / |cos а. |; 3) середина вектора е. пересекается с осью п..
В случае, когда M . й ^. или M. й ^. , т. е. элементы теряют контакт,
сила трения приравнивается к нулю. При повторном пересечении шаров между ними вновь, как и в начальный момент времени, вводятся точки Р., Р.. контакта трения. Из соотношений (9) следует, что сила трения до достижения своего максимального значения в процессе относительного смещения элементов может при увеличении в. не только возрастать, но и с уменьшением в. падать.
Для описания перемещений элементов используются уравнения движения твердых тел
d2X m =
. (10)
d юг ^ n,
= У Nl F. + m . g.
/_(.=i v i°i
J —L = y N (F • r.),
l dt ^—‘V = 1' v V '
где . = 1,..., N N - число элементов); т. - масса; X. - радиус-вектор, проведенный от начала О глобальной системы координат ОХХ2Х3 к центру масс О.
. -го элемента; Fij - сила, действующая со стороны .-го элемента на . -й (. = 1,...,Ni, N. - число контактов . -го элемента); g - ускорение свободного падения; Ji - тензор инерции; ю. - угловая скорость; г. - радиус-вектор, проведенный из центра масс . -го элемента к точке приложения силы Fу. Начальными условиями системы (3) являются:
X (7 = 0 ) = X(0), и ,(( = 0 ) = и(0), ю. (7 = 0) = ю(0), (11)
где и. - линейная скорость центра масс . -го элемента. Как правило, в пределах области О также вводится неподвижная поверхность, ограничивающая движение частиц. Такую поверхность можно задать аналитически или составить из неподвижных дискретных элементов.
Таким образом, получена математическая модель грунта, учитывающая его структуру и специфику взаимодействия между его структурными элементами: соотношения (6) - (9) определяют силы отталкивания, восстанавливаемого сцепления, невосстанавливаемого сцепления и трения, которые могут возникать на контактах между элементами; система уравнений (10) описывает перемещение центра масс и вращение элементов; соотношения (11) задают начальные условия системы уравнений (10). Из приведенной модели следует, что при растяжении прочность грунта определяется когезией и адгезией, при сжатии - в основном отталкиванием и трением. Причем прочность на сжатие оказывается значительно большей, чем на растяжение. Это следует из условий
(2) - деформация материала при сжатии в основном определяется модулем упругости ЕI самих дискретных элементов, при растяжении - модулями упругости Ес., Еа. виртуальных пружин когезии и адгезии. Таким образом, модельная геосреда, как и реальный грунт, является разномодульной и разнопрочностной.
Согласно (1), характеристикой грунта, от которой зависит его прочность, наряду со сцеплением, является угол ф внутреннего трения. В работе [5] показано, что ф можно представить в виде суммы ф = фц + ф„, где фц -
угол контактного трения между частицами, фу - угол дилатансии. Последний существенно зависит от формы и упаковки частиц грунта, а угол фц контактного трения между частицами в разных ситуациях при различных видах нагрузки на среду также может варьироваться в широком диапазоне значений фц в зависимости от формы частиц грунта. Выше в данной работе, как чаще
всего принято и в методе дискретных элементов, рассматривались частицы шаровой формы. Связано это с минимальным объемом вычислений при поиске нормалей п к поверхности контакта шаров, а также наиболее простым представлением тензора инерции, у которого диагональные элементы равны между собой, а остальные - нули. Иногда в МДЭ используется различная полиэдральная форма элементов, суперэллипсоиды и др. Их применение существенно усложняет поиск нормалей п и тензоров инерции таких тел, соответственно, значительно удлиняется время счета задачи. Ниже приведем описание элементов неправильной формы, составленных из элементов-шаров, жестко связанных между собой. Назовем такой комбинированный элемент кластером К.
В случае, когда К составлен из одинаковых шаров, его массу М, радиус-вектор X центра масс и тензор инерции J можно вычислить аналитически [6]. В общем случае, когда кластер состоит из шаров различных размеров, параметры М, X, J найдем численно следующим образом. Пусть кластер К состоит из п произвольных шаров. Считаем, что плотность р у всех шаров одинакова. Их радиус и массу обозначим как г. и т. соответственно, х. - радиус-вектор центра .-го шара в декартовой системе координат 0х1х2х3, . = 1,., п. Пусть кластер целиком лежит внутри прямоугольного параллелепипеда О со сторонами Rk, k = 1, 2, 3. Разобьем О на сетку с кубическими ячейками О.^ с размером Ая, где .к = 1, ., пк, пк - число ячеек по осям прямоугольной системы:
пк = Тгипс (Як / А) +1. (12)
Массу М кластера найдем из соотношения
М = ^ Ат., Ат. =р(Ая)3, (13)
где . е Nm, Nm, - номера ячеек О. ^ ^, центры которых принадлежат кластеру;
АУ. - объем .-й ячейки. Координаты Хк центра масс кластера равны
X = Х. ( X. Ат.) / М, (14)
где хк. - координаты центра .-й ячейки. Компоненты тензора инерции J кластера находятся из соотношений
*^11 = Х. [(Х2. + ХЪ} )Ат. ], J12 = J21 = -Х. (Х1}Х2} Ат. ),
J22 = X . |^(Х1} + Х3. ) Ат. ], J23 = J32 = -Х . (Х2.ХУАт. ), (15)
/33 =Х . [(Х1. + Х2. )Ат. ], J13 = /31 =-Х . (X1jX3j Ат. ).
Приведем тензор инерции / к диагональному виду. Он является симметричным, поэтому для его диагонализации можно применить метод Якоби:
/* = Д/Д, и = Щ Д, и = Пп= Д, (16)
где пв - число итераций, необходимых для достижения заданной точности решения рассматриваемой задачи; /* - диагональная матрица, у которой диагональные компоненты являются собственными значениями / и моментами инерции относительно главных осей инерции кластера; Дг - матрица плоского вращения, проводимого относительно максимального по модулю недиагонального компонента матрицы /1 = Д1-1 /1 Д1 _1,и,-1 - транспонированная
матрица для и,-1. Такое плоское вращение проводится на угол а, который находится из формулы tg(2а) = 2а.. /(а^. -а ). Главные оси инерции кластера направлены вдоль собственных векторов, определяемых столбцами матрицы Д. Таким образом, в базисе }, составленном из найденных собственных векторов, тензор инерции кластера имеет диагональный вид. Базис
{е° }, жестко связан с кластером K, его начало совпадает с центром масс K.
В этом базисе координаты i-го шара кластера находятся из соотношения
x. = Uxt. (17)
Таким образом, для решения задач механики грунтов методом дискретных элементов получены математическая модель среды со структурой и возможность использования в такой модели комбинированных элементов неправильной формы в виде кластеров, составленных из жестко связанных шаров.
Библиографический список
1. Ланис, А.Л. Использование метода напорной инъекции при усилении земляного полотна железных дорог : дис. ... канд. техн. наук. - МГУПС. - М., 2009. - 152 с.
2. Хокни, Р. Численное моделирование методом частиц / Р. Хокни, Дж. Иствуд. - М. : Мир, 1987. - 640 с.
3. Cundall, P.A. A discrete numerical model for granular assemblies / P.A. Cundall,
O.D.L. Strack // Geotechnique. - 1979. - V. 29. - No. 1. - P. 47-65.
4. Хан, Г.Н. О несимметричном режиме разрушения массива горных пород в окрестности полости / Г.Н. Хан // Физическая мезомеханика. - 2008. - Т. 11. - № 1. - С. 109-114.
5. Ревуженко, А.Ф. Об учете дилатансии в основных справочных формулах механики сыпучих сред / А.Ф. Ревуженко, С.Б. Стажевский // ФТПРПИ. - 1986. - № 4. - С. 13-16.
6. Дорофеенко, С.О. Моделирование сыпучих сред методом дискретных элементов : дис. ... канд. физ.-мат. наук. - Черноголовка, 2008. - 115 с.
References
1. Lanis A.L. Ispol'zovanie metoda napomoj in#ekcii pri usilenii zemljanogo polotna zheleznyh dorog : [The use of injection pressure at strengthening of an earthen cloth of railways] dis. ... kand. tehn. nauk. - MGUPS. - Moscow, 2009. - 152 p.
2. Hokni R., Istvud Dzh.Chislennoe modelirovanie metodom chastic [Numerical modelling by a method of particles]. - Moscow : Mir, 1987. - 640 p.
3. Cundall P.A., Strack O.D.L. A discrete numerical model for granular assemblies // Geotechnique. - 1979. - V. 29, No. 1. - P. 47-65.
4. Han G.N. O nesimmetrichnom rezhime razrushenija massiva gornyh porod v okrestnosti po-losti [About the unbalanced destruction of rock mass in the near vicinity of cavity] // Fizi-cheskaja mezomehanika [ Physical Mezomechanics]. - 2008. - T. 11. - No. 1. - P. 109-114.
5. Revuzhenko,A.F., Stazhevskij S.B. Ob uchete dilatansii v osnovnyh spravochnyh formulah me-haniki sypuchih sred [On the inclusion of dilatancy in the basic reference formulas of soil mechanics] // FTPRPI. - 1986. - No. 4. - P. 13-16.
6. Dorofeenko S.O. Modelirovanie sypuchih sred metodom diskretnyh jelementov [Simulation of granular media by method of discrete element] : dis. ... kand. fiz.-mat. nauk. - Chernogolovka, 2008. - 115 p.
