The stability of monolithic lining of the vertical mine workings with the initial porosity of the material and non-elastic compression work skeleton

Full Text

Введение. При добыче полезных ископаемых должно быть пройдено большое количество вертикальных и горизонтальных выработок, которые являются долговременными и дорогостоящими инженерными сооружениями. Состояние горных выработок в зависимости от их назначения должно удовлетворять различным требованиям, основным из которых является обеспечение безопасных условий труда и сохранности оборудования. В связи с этим возникают требования к проведению укрепительных работ горных выработок, то есть созданию крепежных конструкций - крепей. Одним из главных условий безопасной и надежной работы горной выработки является устойчивость ее крепи. В большинстве случаев процесс потери устойчивости происходит при неупругих деформациях, поэтому при моделировании этого явления используются модели, учитывающие одновременно упругие и пластические свойства материалов. Необходимо также отметить, что большинство материалов имеет пористую структуру [1], поэтому при решении задач устойчивости наряду со сложной реологией необходимо также учитывать пористые свойства материалов. Решение такого класса задач состоит из двух этапов. На первом этапе определяется основное напряженно-деформированного состояние (НДС) крепи. На втором - решается задача об устойчивости основного состояния в рамках динамического подхода, который, в свою очередь, сводится к определению величины критического давления, равномерно распределенного по внутреннему контуру крепи. В отличие от [2-5], в настоящей работе на основе точных трехмерных уравнений [6] исследуется вопрос об устойчивости состояния равновесия монолитной крепи вертикальной горной выработки большой протяженности. При этом материал крепи моделируется пористой средой, сжатый скелет которой обладает одновременно упругими и пластическими свойствами. Основное НДС крепи вертикальной горной выработки. Деформирование пористого материала с начальным раствором пор ε0 можно разделить на два взаимосвязанных этапа [7]. Первый - упругое деформирование сжимаемой пористой среды, второй - неупругое деформирование сжатого скелета с упрочняющимися упругопластическими свойствами. Связь между напряжениями и деформациями на первом этапе берется в виде закона Гука для сжимаемого тела e e σjβ = λ1 εαα gjβ + 2µ1 εβj gjβ , σjβ , e εβj e при - εαα < ε0 , (1) где - смешанные компоненты метрического тензора, тензора напряжений и тензора упругих деформаций соответственно, λ1 , µ1 - параметры Ламе сжимаемого тела. 458 Устойчивость монолитной крепи вертикальной горной выработки . . . На втором этапе упругие деформации сжатого скелета подчиняются закону Гука для несжимаемого тела [8] e e Sjβ = 2 (µ0 + µ1 ) εβj -2µ0 εβj 0 e 2 + µ1 ε0 gjβ , при - εαα = ε0 , 3 (2) e где Sjβ - компоненты тензора девиатора напряжений; εβj 0 - компоненты тензора упругих деформаций, вычисленные на момент полного сжатия пор, то e есть при выполнении условия εαα = -ε0 ; µ0 + µ1 - модуль сдвига несжимаемого тела. В зоне пластического деформирования сжатого скелета будем использовать модель несжимаемого упрочняющегося упругопластического тела [9] с поверхностью нагружения p p F = Sjβ - c εβj Sβj - c εjβ - k 2 , (3) p где εβj - компоненты тензора пластических деформаций, c - коэффициент упрочнения, k - предел текучести материала. Полная деформация в пластической зоне складывается из упругой и пластической составляющих εβj = e εβj p + εβj , (4) причем пластическая и упругая составляющие объемной деформации соответственно удовлетворяют условиям несжимаемости p εnn = 0, e εnn = -ε0 . (5) В (1)-(5) и далее индексы «e» и «p» вверху величин обозначают их принадлежность соответственно к упругой и пластической зонам деформирования сжатого скелета. Далее рассмотрим задачу определения НДС цилиндрического тела (рис. 1), являющегося крепью вертикальной выработки. Обозначим через b и a соответственно внешний и внутренний радиусы крепи. Действие массива горных пород на крепь заменим сжимающей нагрузкой интенсивностью qb , равномерно распределенной по внешней поверхности. Сжимающая нагрузка интенсивностью qa , равномерно распределенная по внутренней поверхности, моделирует собой давление жидкости или газа на крепь. Так как выработка имеет большую протяженность, то, согласно [5], при определении НДС не будем учитывать эффекты, связанные с тем, что выработка имеет конечную глубину. НДС монолитной крепи вертикальной горной выработки в рамках плоского деформированного состояния в цилиндрической системе координат (r, θ, z) моделируется следующими соотношениями геометрически линейной теории: - уравнение равновесия σr - σθ dσr + = 0; (6) dr r 459 Г о ц е в Д. В., Б у н т о в А. Е. Рис. 1. Монолитная крепь вертикального шахтного ствола под действием радиального сжатия [Figure 1. Monolithic lined circular vertical shaft under the radial contraction action] - соотношения Коши u du , εθ = , dr r где u - радиальная составляющая вектора перемещений; - граничные условия в напряжениях εr = σr r=b = -qb , σr r=a = -qa (qa > 0, qb > 0). (7) (8) Связь между напряжениями и деформациями при упругом деформировании пористой среды возьмем в виде соотношений (1), которые при принятых допущениях перепишутся в виде σr = (λ1 + 2µ1 ) εr + λ1 εθ , σθ = λ1 εr + (λ1 + 2µ1 ) εθ , σz = λ1 (εr + εθ ) . (9) Упругие деформации сжатого скелета связаны с напряжениями соотношениями (2), которые в нашем случае примут вид 2 Sr = 2(µ0 + µ1 )εr - 2µ0 εr0 + µ1 ε0 , 3 2 Sθ = 2(µ0 + µ1 )εθ - 2µ0 εθ0 + µ1 ε0 , 3 2 S z = µ 1 ε0 . 3 (10) В (10) и далее индекс «0» внизу компонент деформаций, напряжений и перемещений обозначает, что они вычислены на момент полного сжатия пор. Функция нагружения (3), соотношения для полных деформаций в пластической зоне сжатого скелета (4), условия несжимаемости (5) в случае плоского деформированного состояния для рассматриваемой задачи перепишутся 460 Устойчивость монолитной крепи вертикальной горной выработки . . . соответственно следующим образом: (Sr - cεpr )2 + (Sθ - cεpθ )2 + Sz2 = 2k 2 , εer + εr = εer + εpr , εeθ + εez = -ε0 , εθ = εpr + εpθ , εpθ + εpz εeθ + (11) (12) = 0. (13) На границе γ раздела сред упругого и пластического деформирования сжатого скелета должны выполняться условия непрерывности перемещений и напряжений: [u]|r=γ = 0, [σr ]|r=γ = 0, [σθ ]|r=γ = 0. (14) В (14) и далее квадратные скобки обозначают разность значений выражений, соответствующих упругой и пластической областям на границе r = γ. Соотношения (6)-(14) представляют собой математические модели, описывающие НДС монолитной крепи вертикального шахтного ствола на этапах упругого деформирования пористой среды и неупругого деформирования сжатого скелета. НДС монолитной крепи на первом этапе, то есть при наличии несхлопнутых пор, согласно (6)-(9) определяется соотношениями (qb - qa a2 )r (qb - qa )a2 + , 2(λ1 + 1)(a2 - 1) 2(a2 - 1)r qb - qa a2 (qb - qa )a2 εr = - , 2(λ1 + 1)(a2 - 1) 2(a2 - 1)r2 qb - qa a2 (qb - qa )a2 εθ = + , 2 2(λ1 + 1)(a - 1) 2(a2 - 1)r2 qa a2 (r2 - 1) qb (a2 - r2 ) σr = 2 + 2 , r (1 - a2 ) r (1 - a2 ) qa a2 (r2 + 1) qb (r2 + a2 ) σθ = 2 - 2 , r (1 - a2 ) r (1 - a2 ) qa λ1 a2 qb λ 1 σz = - . 2 (λ1 + 1)(1 - a ) (λ1 + 1)(1 - a2 ) u= (15) В (15) и далее все соотношения записаны в безразмерном виде, при этом все величины имеющие размерность напряжений отнесены к величине µ1 , а имеющие размерность длины - к радиусу b. Из (15) следует, что объемная деформация при упругом сжатии пор определяется в виде qb - qa a 2 εαα = εr + εθ = . (λ1 + 1)(a2 - 1) Следовательно, достижение величины начального раствора пор нулевого значения (иначе - достижение объемной деформацией величины ε0 ) при упругом деформировании материала происходит одновременно во всей крепи под действием нагрузок, удовлетворяющих условию qb = ε0 (λ1 + 1)(1 - a2 ) + qa f (ε0 )a2 , 461 Г о ц е в Д. В., Б у н т о в А. Е. 1, если ε0 = 0, 0, если ε0 = 0. При этом НДС (15) на момент полного закрытия пор перепишется в виде где f (ε0 ) = ε0 r (qa f (ε0 ) - ε0 (λ1 + 1))a2 + , 2 2r ε0 (qa f (ε0 ) - ε0 (λ1 + 1))a2 , εr0 = - - 2 2r2 ε0 (qa f (ε0 ) - ε0 (λ1 + 1))a2 εθ0 = - + , 2 2r2 (qa f (ε0 ) - ε0 (λ1 + 1))a2 σr0 = -ε0 (λ1 + 1) - , r2 (qa f (ε0 ) - ε0 (λ1 + 1))a2 σθ0 = -ε0 (λ1 + 1) + , r2 σz0 = -λ1 ε0 . u0 = - (16) Таким образом, если qb < ε0 (λ1 + 1)(1 - a2 ) + qa f (ε0 )a2 , то полного закрытия пор не происходит и материал ведет себя как сжимаемая упругая среда с параметрами λ1 , µ1 = 1, ε0 . НДС монолитной крепи на втором этапе деформирования согласно (6)- (14), (16) определяется следующим образом: - в упругой области (γ < r < 1): σr = χγ 2 σθ = χγ 2 ε20 1 1 - 2 - qb , 3 r 2 ε 1 k 2 - 0 1 + 2 - qb ; 3 r k2 - (17) - в пластической области (a < r < γ): εpr = -εpθ = 1 χ c + 2µ σr = -qa + χ σθ = -qa + χ ε20 γ2 1- 2 , 3 r γ2 γ2 2µ r2 r - 2 + 1 - 2 + 2 ln , 2 a r c + 2µ a a γ2 γ2 4µ 1 r2 γ2 r + + 3 - - + ln 2 2 2 2 a r c + 2µ 2 a r a k2 - ε20 3 ε2 k2 - 0 3 k2 - (18) . Перемещения и полные деформации в упругой и пластической областях определяются соотношениями u= 462 D ε0 r - , r 2 εr = - D ε0 - , r2 2 εθ = D ε0 - . r2 2 (19) Устойчивость монолитной крепи вертикальной горной выработки . . . Здесь в (17)-(19) введены следующие обозначения: χ = (qa - qb ), µ = 1 + µ0 , ε2 1 χγ 2 k 2 - 0 + µ0 qa f (ε0 ) - ε0 (λ1 + 1) a2 . 2µ 3 Радиус γ раздела зон упругого и пластического деформирования определяется из решения уравнения D= qb - qa + χ k2 - ε20 3 γ2 2µ γ γ2 2 + - γ + 2 ln 1 - a2 c + 2µ a2 a = 0. Результаты численного эксперимента, проводимого в рамках полученных решений, описывающих НДС крепи вертикальной выработки, представлены на рис. 2, 2. На рис. 2, a кривая 1 соответствует k = 0.01, кривая 2 - k = 0.015, кривая 3 - k = 0.02. На рис. 2, b кривая 1 соответствует µ = 1, кривая 2 - µ = 2, кривая 3 - µ = 3. На рис. 3, a и b кривые 1 соответствуют k = 8.5·10-3 , кривые 2 - k = 0.01, кривые 3 - k = 0.012. На рис. 3, c и d кривые 1 соответствуют ε0 = 10-4 , кривые 2 - ε0 = 7 · 10-3 , кривые 3 - ε0 = 0.01. Безразмерные значения других физико-механических и геометрических параметров, если не оговорено особо, в расчетах брались следующими: a = = 0.5, b = 1, qa = 0.001, qb = 0.012, c = 0.005, λ1 = 3, µ1 = 1, k = 0.01, ε0 = 0.001, µ = 2. Знание основного НДС крепи позволяет перейти к решению задачи устойчивости этого состояния. Пространственная устойчивость монолитной крепи вертикальной горной выработки. Исследование устойчивости основного состояния (17)-(19) монолитной цилиндрической крепи со сжатым скелетом при принятии обобщенной концепции продолжающегося нагружения [6] сводится к решению системы дифференциальных уравнений в вариациях при соответствующих граничных условиях [5]. Уравнения равновесия для областей пластического и упругого деформирования материала крепи со сжатыми порами имеют вид ∇i (σji + σα0i ∇α uj ) = 0. (20) Здесь и далее ∇ - символ ковариантного дифференцирования, кружок вверху соответствует компонентам основного невозмущенного состояния, определенного соотношениями (17)-(19). Граничные условия на внутренней и внешней поверхностях крепежной конструкции запишем в виде 0 Ni σji + σαi ∇α uj = 0. (21) Условия непрерывности на упругопластической границе γ имеют вид 0 Ni σji + σαi ∇α uj = 0, [uj ] = 0. (22) 463 Г о ц е в Д. В., Б у н т о в А. Е. Рис. 2. Зависимость радиуса упругопластической границы γ от начального раствора пор ε0 [Figure 2. The dependence of the elastic-plastic border radius γ on the value of ε0 ] Рис. 3. Зависимость напряжений σr и σθ от радиальной координаты r [Figure 3. The dependence of the stresses σr and σθ on the radial coordinate value r] 464 Устойчивость монолитной крепи вертикальной горной выработки . . . Связь между амплитудными значениями напряжений и перемещений для сжатого скелета, обладающего упругопластическими свойствами и свойством дальнейшей несжимаемости в пластической и упругой областях, представима в форме [5] σjβ = (xβα g αα ∇α uα + p)gjβ + (1 - gjβ )g ββ µ(∇β uj + ∇j uβ ). (23) Здесь в (23) отсутствует суммирование по индексам j, β и производится по индексу α. Величины xβα в пластической области определяются соотношением xβα = 2µgβα - υfαα fββ , (24) где υ= 4µ2 , k 2 (2µ + c) 0 0p fij = sij -c εij , в упругой области - соотношением (24), где надо положить υ = 0. Условие несжимаемости для материала крепи со сжатым скелетом представимо в форме ∇α uα = 0. (25) Уравнения (20)-(23) с учетом условия несжимаемости (25) в областях пластического и упругого деформирования материала со сжатым скелетом представляют собой взаимосвязанную замкнутую систему уравнений для исследования устойчивости основного состояния крепи вертикальной цилиндрической выработки, когда имеется граница раздела областей упругого и пластического поведения материала при нагружении. Система уравнений (20), (23), (25) представляет собой систему дифференциальных уравнений в частных производных относительно амплитудных значений векторов перемещений u, ν, w и гидростатического давления для пластической и упругой зон крепи. Нетривиальное решение этой задачи соответствует потере устойчивости основного состояния. Для нахождения собственных значений перемещения и гидростатические давления в каждой из зон упругого и пластического деформирования аппроксимируем двойными тригонометрическими рядами ∞ ∞ n ∞ m ∞ u= Anm (r) cos(mθ) cos(nz), ν= Bnm (r) sin(mθ) cos(nz), n m ∞ ∞ w= (26) Cnm (r) cos(mθ) sin(nz), n ∞ m ∞ n m p= Dnm (r) cos(mθ) cos(nz), где n, m - параметры волнообразования. Так как система уравнений (20)-(23), (25) является линейной и однородной, ее можно записать для каждого члена с одинаковыми значениями n, m. Для упрощения записи в дальнейшем индексы n и m будем опускать. 465 Г о ц е в Д. В., Б у н т о в А. Е. Подставляя u, ν, w, p, определяемые равенствами (26), в краевую задачу (20)-(22) и учитывая (23), (25), после ряда преобразований получим краевую задачу в терминах функций A(r), B(r), D(r). Уравнения равновесия имеют следующий вид: 0 0 0 0 ξ1 (r)A(r) + ξ2 (r)A (r) + ξ3 (r)A (r) + ξ4 (r)B(r)+ 0 + ξ5 (r)B (r) + rD (r) = 0, 0 0 0 0 ξ6 (r)A(r) + ξ7 (r)A (r) + ξ8 (r)B(r) + ξ9 (r)B (r)+ (27) 0 + ξ10 (r)B (r) + mD(r) = 0, 0 0 0 0 ξ11 (r)A(r) + ξ12 (r)A (r) + ξ13 (r)A (r) + ξ14 (r)rA (r)+ 0 0 0 + ξ11 (r)mB(r) + ξ15 (r)B (r) + ξ14 (r)mB (r) + D(r)r3 n2 = 0, где 0 ξ1 (r) = 0 0 1 2 2 0 0 r ρω - (σθ +µ)(1 + m2 ) + υ ψ2 -r2 (σz +µ)n2 , r 0 0 0 0 0 0 ξ2 (r) = σr +r σr,r -υ ψ1 +µ, ξ3 (r) = r(-υ ψ1 +µ + σr ), 0 0 0 0 m 0 ξ4 (r) = (-2 σθ -2µ + υ ψ2 ), ξ5 (r) = -mυ ψ3 , r 0 0 0 0 m 0 ξ6 (r) = - (2 σθ -υ ψ2 +4µ), ξ7 (r) = m(υ ψ3 -2µ), r 0 0 1 2 2 0 0 0 ξ8 (r) = r ρω - σθ -r2 (µ + σz )n2 - (3µ - υ ψ2 + σθ )m2 - µ , r 0 0 0 0 ξ9 (r) = σr +r σr,r +µ, 0 0 ξ10 (r) = r(σr +µ), 0 0 0 0 ξ11 (r) = r2 ρω 2 - n2 (2µ + σz ) - m2 (µ + σθ ) + µn2 r2 + µ + σr -r σr,r , 0 0 0 0 0 ξ12 (r) = r r2 ρω 2 - n2 (2µ + σz ) - m2 (µ + σθ ) + µn2 r2 - µ - σr +r σr,r , 0 0 0 0 ξ13 (r) = r2 2(µ + σr ) + r σr,r , 0 0 ξ14 (r) = r2 (µ + σr ), 0 0 ξ15 (r) = -rm(µ + σr -r σr,r ), 0 ψ1 = µ 1 ε0 +χ 3 1 k 2 - µ21 ε20 3 0 ψ3 = 0 2 , ψ2 = 2µ1 ε0 3 2 µ 1 ε0 -χ 3 1 k 2 - µ21 ε20 3 2 , - k2 . При этом соотношения (27) будут справедливы для пластической области (a < r < γ) деформирования крепи, если в них всем переменным величинам приписать вверху индекс «p», и для упругой (γ < r < 1), если в этих соотношениях положить υ = 0 и всем переменным величинам приписать вверху индекс «e». 466 Устойчивость монолитной крепи вертикальной горной выработки . . . Граничные условия на внутренней и внешней поверхностях крепи соответственно запишутся в виде 0 0 (Ap (a) + B p (a)m) φ1 + φ2 A p (a) + Dp (a) = 0, 0 (28) mµAp (a) + µB p (a) + φ3 B p (a) = 0, 0 φ4 Ap (a) + (A p (a)a + A p (a)a2 - B p (a)m + B 0 A e (1) φ6 +De (1) = 0, 0 mµAe (1) + µB e (1) - φ5 B e (1) 0 Ae (1) φ7 +(A e (1) 0 p (a)am) φ 5 = 0; (29) = 0, + A e (1) - B e (1)m + B 0 e (1)m) φ 5 = 0, где 0 0 φ1 0 υ ψ3 =- , a 0 0 φ2 = -qa - υ ψ1 +2µ, 0 0 φ 4 = n 2 a 2 µ - µ - qa , φ3 = -a(µ - qa ), 0 φ 5 = µ - qa , 0 φ 7 = n 2 µ - µ - qa . φ6 = 2µ - qb , Условия непрерывности на границе γ зон упругого и пластического деформирования примут вид 0 0 0 φ8 [A ] + [D] + φ9 (Ap + mB p ) - A p υ ψ1 = 0, 0 (30) mµ[A] + µ[B] + φ10 [B ] = 0, 0 φ11 [A] - γ[A ] - γ 2 [A ] + m[B] - γm[B ] = 0, где 0 φ8 0 = 0 0 σre (γ) 0 + 2µ, φ10 = -γ σre (γ) + µ , 0 φ9 υ ψ3 , =- γ 0 φ11 = 1 - µn2 γ 2 0 σre (γ) . +µ Таким образом, математическая модель отказа монолитной крепи вертикальной горной выработки представляется в виде бесконечной системы обыкновенных дифференциальных уравнений (27) с краевыми условиями (28), (29) и условиями сопряжения (30). Для решения этой системы уравнений будем использовать метод конечных разностей [1], суть которого заключается в замене производных разностными отношениями, что, в свою очередь, приводит к замене дифференциальных уравнений разностными, решение которых в любой конечной области сводится к решению конечной системы алгебраических уравнений. Для непрерывной на отрезке функции y = y(r) это осуществляется следующим образом: отрезок [a, b] разбивается точками i = 0, 1, . . . , M на S равных промежутков длины τ = (b - a)/M . Производные функций на 467 Г о ц е в Д. В., Б у н т о в А. Е. данном интервале заменяются во внутренних точках i = 2, 3, . . . , M - 2 центральными разностями второго порядка аппроксимации: 1 (yi+1 - yi-1 ), 2τ 1 yi = 2 (yi+1 - 2yi + yi-1 ), τ 1 yi = 3 (yi+2 - 2yi+1 + 2yi-1 - yi-2 ). 2τ yi = (31) В приграничных точках r1 и rM -1 производные будем аппроксимировать следующими конечными разностями: 1 (yi+1 - yi-1 ), 2τ 1 yi = 2 (yi+1 - 2yi + yi-1 ), τ 1 yi = ± 3 (-yi±3 + 6yi±2 - 12yi±1 + 10yi - 3yi∓1 ), 2τ yi = (32) здесь верхние знаки соответствуют i = 1, нижние i = M - 1. В крайних точках r0 и rM производные будем аппроксимировать следующими конечными разностями: 1 yi = ± (-3yi + 4yi±1 - yi±2 ), 2τ 1 yi = 2 (2yi - 5yi±1 + 4yi±2 - yi±3 ), τ 1 yi = ± 3 (-5yi + 18yi±1 - 24yi±2 + 14yi±3 - 3yi±4 ), 2τ (33) где верхние знаки соответствуют i = 0, нижние - i = M . Обозначим через M p и M e количество точек разбиения в областях пластического и упругого деформирования материала крепи соответственно. Заменяя производные функций A(r), B(r) и D(r) в (27) для пластической области V p в точках j = 2, 3, 4, . . . , M p - 1 и для упругой области V e в точках j = M p + 3, M p + 4, . . . , M p + M e через конечные разности (31), получим 0 0 0 0 0 0 (2ξ3 -τ ξ2 )Aj-1 + (2τ 2 ξ1 -4 ξ3 )Aj + Aj+1 (τ ξ2 +2 ξ3 )- 0 0 0 0 0 -τ ξ5 Bj-1 + 2τ 2 ξ4 Bj + τ ξ5 Bj+1 - rτ ξ2 Dj-1 + rτ ξ2 Dj+1 = 0, 0 0 0 0 0 -τ ξ7 Aj-1 + 2τ 2 ξ6 Aj + τ ξ7 Aj+1 + (2 ξ10 -τ ξ9 )Bj-1 + 0 0 0 0 +2(τ 2 ξ8 -2 ξ10 )Bj + (τ ξ9 +2 ξ10 )Bj+1 + mDj = 0, 0 - rξ14 Aj-2 0 0 0 0 0 + (2τ ξ13 -τ 2 ξ12 + 2rξ14 )Aj-1 + (2τ 3 ξ11 -4τ ξ13 )Aj + 0 0 0 0 +(τ 2 ξ12 +2τ ξ13 -2 rξ14 )Aj+1 + rξ14 Aj+2 + 0 0 0 0 +(2τ mξ14 -τ 2 ξ15 )Bj-1 + (2τ 3 m ξ11 -4τ mξ14 )Bj + 0 0 +(τ 2 ξ15 +2τ mξ14 )Bj+1 + 2τ 3 r3 n2 Dj = 0. 468 (34) Устойчивость монолитной крепи вертикальной горной выработки . . . В точках j = M p , M p + M e + 1, согласно тому, что в приграничных точках производные будем аппроксимировать схемой (32) (нижние знаки), для уравнений равновесия (27) запишем 0 0 0 0 0 0 (2ξ3 -τ ξ2 )Aj-1 + (2τ 2 ξ1 -4 ξ3 )Aj + (τ ξ2 +2 ξ3 )Aj+1 - 0 0 0 0 0 -τ ξ5 Bj-1 + 2τ 2 ξ4 Bj + τ ξ5 Bj+1 - rτ ξ2 Dj-1 + rτ ξ2 Dj+1 = 0, 0 0 0 0 0 -τ ξ7 Aj-1 + 2τ 2 ξ6 Aj + τ ξ7 Aj+1 + (2 ξ10 -τ ξ9 )Bj-1 + 0 0 0 0 +2(τ 2 ξ8 -2 ξ10 )Bj + (τ ξ9 +2 ξ10 )Bj+1 + mDj = 0, (35) 0 r ξ14 Aj-3 0 0 0 0 - 6r ξ14 Aj-2 + (2τ ξ13 -τ 2 ξ12 +12r ξ14 )Aj-1 + 0 0 0 0 0 0 +2(τ 3 ξ11 -2τ ξ13 -5r ξ14 )Aj + (τ 2 ξ12 +2τ ξ13 +3r ξ14 )Aj+1 + 0 0 0 0 +τ (2m ξ14 -τ ξ15 )Bj-1 + 2τ m(τ 2 ξ11 -2 ξ14 )Bj + 0 0 +τ (τ ξ15 +2m ξ14 )Bj+1 + 2τ 3 r3 n2 Dj = 0. В точках j = 1, M p + 2, согласно тому, что в приграничных точках производные будем аппроксимировать схемой (32) (верхние знаки), для уравнений равновесия (27) запишем 0 0 0 0 0 0 Aj-1 (2ξ3 -τ ξ2 ) + Aj (2τ 2 ξ1 -4 ξ3 ) + Aj+1 (τ ξ2 +2 ξ3 )- 0 0 0 0 0 -τ ξ5 Bj-1 + 2τ 2 ξ4 Bj + τ ξ5 Bj+1 - rτ ξ2 Dj-1 + rτ ξ2 Dj+1 = 0, 0 0 0 0 0 -τ ξ7 Aj-1 + 2τ 2 ξ6 Aj + τ ξ7 Aj+1 + (2 ξ10 -τ ξ9 )Bj-1 + 0 0 0 0 +2(τ 2 ξ8 -2 ξ10 )Bj + (τ ξ9 +2 ξ10 )Bj+1 + mDj = 0, (2τ 0 0 0 0 0 0 ξ13 -τ 2 ξ12 -3r ξ14 )Aj-1 + 2(τ 3 ξ11 -2τ ξ13 +5r ξ14 )Aj + 0 0 0 0 0 +(τ 2 ξ12 +2τ ξ13 -12r ξ14 )Aj+1 + 6r ξ14 Aj+2 - r ξ14 Aj+3 + 0 0 0 0 +τ (2m ξ14 -τ ξ15 )Bj-1 + 2τ m(τ 2 ξ11 -2 ξ14 )Bj + 0 0 +τ (τ ξ15 +2m ξ14 )Bj+1 + 2τ 3 r3 n2 Dj = (36) 0. Шаги разностной сетки в каждой из зон пластического и упругого деформирования следующие: τp = γ-a , Mp + 1 τe = 1-γ . Me + 1 Отметим, что если соотношения (34)-(36) записываютcя для пластической области, то в них надо приписать всем неизвестным и переменным величинам вверху индекс «p», для упругой области - индекс «e». Производные функций A(r), B(r) и D(r) в условиях (28)-(30) аппроксимируются конечными разностями типа (33). При этом условиям (28) соответ469 Г о ц е в Д. В., Б у н т о в А. Е. ствует разностный вид 0 0 0 0 0 (2τ p φ1 -3 φ2 )Ap0 + 4 φ2 Ap1 - φ2 Ap2 + 2τ p m φ1 B0p + 2τ p D0p = 0, 0 0 0 2τ p mµAp0 + (2τ p µ - 3 φ3 )B0p + 4 φ3 B1p - φ3 B2p = 0, (37) 0 0 2τ 2 (φ4 / φ5 ) - a(3τ p - 4a) Ap0 + 2a(2τ p - 5a)Ap1 + a(8a - τ p )Ap2 - -2a2 Ap3 - τ p m(3a + 2τ p )B0p + 4τ p amB1p - τ p amB2p = 0, условиям (29) - вид 0 e = 0, -AeM -2 + 4AeM -1 - 3AeM - (2τ e / φ6 )DM 0 0 0 e e e e 2τ e mµAeM - φ5 BM -2 + 4 φ5 BM -1 + (2τ µ - 3 φ5 )BM = 0, 0 (38) 0 -2AeM -3 + (8 + τ e )AeM -2 + 4 + 3τ e + 2(τ e )2 (φ7 / φ5 ) AeM - e e e -2(5 + 2τ e )AeM -1 + τ e mBM -2 - 4τ mBM -1 + e = 0, +τ e m(3 - 2τ e )BM где M = M p + M e + 2; условиям (30) - вид 0 0 0 0 0 0 τ e (φ8 +υ ψ1 )(-Apj-2 + 4Apj-1 ) + 2τ e τ p φ9 -3τ e (φ8 +υ ψ1 ) - 3τ p φ8 Aj + 0 0 +τ p φ8 (4Aej+1 - Aej+2 ) + 2τ e τ p m φ9 Bj = 0, p p e e ) = 0, τ e (-Bj-2 + 4Bj-1 ) - 3(τ e + τ p )Bj + τ p (4Bj+1 - Bj+2 2γ(τ e )2 Apj-3 - + 8γ)Apj-2 + 2(τ e )2 (2τ p + 5γ)Apj-1 + +(4γ((τ p )2 - (τ e )2 ) - 3τ e τ p (τ e + τ p ))Aj + +2(τ p )2 (2τ e - 5γ)Aej+1 + (τ p )2 (8γ - τ e )Aej+2 - 2γ(τ p )2 Aej+3 (39) (τ e )2 (τ p = 0. Здесь в (39) j = M p + 1. В результате получаем систему линейных однородных алгебраических уравнений, которая в матричной форме может быть представлена в виде [Xij ] × [Yj ] = 0. (40) Здесь [Yj ] - вектор неизвестных, который имеет следующие компоненты: [Yj ] = (Ap0 , Ap1 , . . . , ApM p , ApM p +1 , AeM p +1 , AeM p +2 , . . . , AeM p +M e +1 , AeM p +M e +2 p p e e e e B0p , B1p , . . . , BM p , BM p +1 , BM p +1 , BM p +2 , . . . , BM p +M e +1 , BM p +M e +2 , p p p p e e e e D0 , D1 , . . . , DM p , DM p +1 , DM p +1 , DM p +2 , . . . , DM p +M e +1 , DM p +M e +2 ). Таким образом, определение критического значения внутреннего давления интенсивностью qa , соответствующего потере устойчивости монолитной крепи вертикальной горной выработки, сводится к разрешимости матричного уравнения, что, в свою очередь, соответствует равенству нулю определителя 470 Устойчивость монолитной крепи вертикальной горной выработки . . . конечной разностной системы уравнений (34)-(36), (37)-(39) при ограничении числа членов в рядах (26). При вычислении определителя наряду с нахождением основного напряженно-деформированного состояния крепи для каждой из областей упругого и пластического деформирования необходимо учитывать уравнение, определяющее положение упругопластической границы γ. Минимизация должна производиться по шагу разностной сетки, параметрам волнообразования по контуру m и образующей n, параметрам материала и конструкции λj . Таким образом, получаем задачу многомерной оптимизации величины qa в зависимости от m, n при условии равенства нулю определителя полученной алгебраической системы (40). Результаты вычислительного эксперимента представлены на рис. 4, 5. На рис. 4 показана зависимость критической величины внутреннего давления интенсивностью qa от внешней нагрузки интенсивности qb , моделирующей собой давление массива горных пород на крепь при различных значениях внутреннего радиуса a монолитной крепи. При этом кривая 1 соответствует a = 0.6, кривая 2 - a = 0.5, кривая 3 - a = 0.4. На рис. 5 представлена зависимость критического давления интенсивностью qa от относительного внутреннего радиуса a монолитной крепи при различных значениях коэффициента упрочнения c материала крепи с полностью сжатой матрицей. При этом кривая 1 соответствует c = 0.01, кривая 2 - c = 0.04, кривая 3 - c = 0.08. В расчетах принималось qb = 0.03. Всем расчетам, результаты которых отражены на рис. 4, 5, соответствуют значения параметров волнообразования m = n = 3. Безразмерные значения других физико-механических и геометрических параметров, если не оговорено особо, брались следующими: k = 0.007, λ1 = 2, µ1 = 1, ε0 = 8 · 10-4 . Анализ численных расчетов позволяет сделать следующие выводы: - потеря устойчивости монолитной крепи вертикальной горной выработки происходит по осесимметричной форме, которой соответствуют значения параметров волнообразования m = n = 3; Рис. 4. Зависимость критической величины внутреннего давления от внешней нагрузки [Figure 4. The dependence of the critical value of internal pressure on the external load value] Рис. 5. Зависимость критической величины внутреннего давления от внутреннего радиуса монолитной крепи [Figure 5. The dependence of the critical value of internal pressure on the inner radius of the monolithic lining] 471 Г о ц е в Д. В., Б у н т о в А. Е. - область устойчивости монолитной крепи вертикальной горной выработки существенно зависит как от физико-механических, так и от геометрических параметров конструкции, при этом как при увеличении относительной толщины крепи, так и с ростом относительного коэффициента упрочнения c материала крепи с полностью сжатой матрицей область устойчивости расширяется; - с ростом давления горного массива интенсивностью qb на внешнюю поверхность крепи значения внутреннего критического давления, соответствующие потере устойчивости, увеличиваются. Декларация о финансовых и других взаимоотношениях. Исследование не имело спонсорской поддержки. Все авторы принимали участие в разработке концепции статьи и в написании рукописи. Авторы несут полную ответственность за предоставление окончательной рукописи в печать. Окончательная версия рукописи была одобрена всеми авторами. Авторы не получали гонорар за статью.

About the authors

Dmitry V Gotsev

Russian Air Force Military Educational and Scientific Center of the “N. E. Zhukovskiy and Yu. A. Gagarin Air Force Academy”

Email: rbgotsev@mail.ru
(Dr. Phys. & Math. Sci.; rbgotsev@mail.ru), Professor, Dept. of Mathematics 54 a, Staryh Bolshevikov str., Voronezh, 394064, Russian Federation

Alexey E Buntov

Russian Air Force Military Educational and Scientific Center of the “N. E. Zhukovskiy and Yu. A. Gagarin Air Force Academy”

Email: alexey.buntov@mail.ru
(alexey.buntov@mail.ru; Corresponding Author), Postgraduate Student in a Military Academy, Dept. of the Airport Engineering Support 54 a, Staryh Bolshevikov str., Voronezh, 394064, Russian Federation

References

Supplementary files