CC BY
39
УДК 622.273.2:622.831
А.П.ГОСПОДАРИКОВ, д-р техн. наук, профессор, kafmat_@mail.ru М.А.ЗАЦЕПИН, канд. физ.-мат. наук, ассистент, michael_zatsepin@mail.ru А.В.МЕЛЕШКО, аспирант, meleshkoandrey@rambler.ru Санкт-Петербургский государственный горный университет
A.P.GOSPODARIKOV, Dr. in eng. sc., professor, kafmat_@mail.ru M.A.ZATSEPIN, PhD in phys. & math., assistant lecturer, michael_zatsepin@mail.ru A.V.MELESHKO, post-graduate student, meleshkoandrey@rambler.ru Saint Petersburg State Mining University
ОБ ОДНОМ АЛГОРИТМЕ ЧИСЛЕННОГО РЕШЕНИЯ НЕЛИНЕЙНЫХ КРАЕВЫХ ЗАДАЧ ГЕОМЕХАНИКИ
Изложены алгоритмы прогноза напряженно-деформированного состояния неоднородного слоистого физически нелинейного породного массива. Моделирование базируется на применении комплекса вычислительных методов: вариационного, дискретного продолжения по числовому параметру, квазилинеаризации нелинейных краевых задач, конечных разностей, матричной прогонки и общего итерационного процесса.
Ключевые слова: прогноз, геомеханика, массив горных пород, напряженно-деформированное состояние, вычислительные методы.
THE ALGORITHM OF NUMERICAL SOLUTION THE NON-LINEAR BOUNDARY VALUE PROBLEMS OF GEOMECHANICS
This article is devoted to the algorithm of the forecast of stress and strain state nonho-mogeneous layered physical nonlinear rock massif. The modeling based on applying the whole package of computational methods: variational method, discrete continuation on numeric parameter method, quasi-linearization of nonlinear value boundary problem, finite difference method, Thomas algorithm and iterative process.
Key words: forecast, geomechanics, rock massif, stress and strain state, computational methods.
Разработка пологих месторождений полезных ископаемых приводит исследователей к необходимости решения важных прикладных задач геомеханики: изучение вопросов проявления горного давления и напряженно-деформированного состояния (НДС) массива горных пород (МГП) в окрестности выработок различного назначения и произвольной конфигурации. Специалистам необходимо иметь информацию о протекании в этих массивах разнообразных физико-механических процессов, которые зависят не только от специфики горно-геологических условий отработки различных ме-
сторождений полезных ископаемых, но и от изменчивости этих условий в пределах шахтных полей.
Очевидно, что изучение параметров механических процессов в МГП не может быть однозначно определено только на основе использования данных натурных экспериментов, либо данных лабораторных исследований или результатов аналитических расчетов. По-видимому, постановка такой задачи и нецелесообразна, так как даже ведущим представителям горной науки, специализирующимся в области горной геомеханики, не удалось до сих пор установить
306 _
ISSN 0135-3500. Записки Горного института. Т. 196
общие физико-механические законы, адекватно отражающие реальное поведение МГП. Поэтому по-прежнему остается актуальной задача исследования механических процессов в породных массивах различными методами (от экспериментальных натурных измерений до вычислительных процедур).
Важные прикладные задачи горной геомеханики решались многими исследователями различными аналитическими и численными методами. В частности, хорошо известны результаты численно-аналитических расчетов оседания многослойного массива над выработанным пространством и опорного давления впереди очистного забоя на основе вариационного метода В.З.Власова в рамках линейного закона Гука (А.А.Борисов, Н.Н.Кайдалов, В.Г.Лабазин и др.), а также примененного в работах А.П.Госпо-дарикова нелинейного закона деформирования горных пород [2]. Если же массив горных пород ослаблен выработками непрямоугольного сечения, то использовать в явном виде уравнения равновесия, полученные В.З.Власовым, не представляется возможным. В этом случае необходимо получить основные соотношения для определения НДС МГП в полярной системе координат, базируясь на идеях технической теории В.З.Власова.
Рассмотрим массив горных пород мощностью h, ограниченный двумя горизонтальными плоскостями и боковой цилиндрической поверхностью. Для определения НДС упругого МГП применим метод перемещений, приняв за основные неизвестные перемещения и = и( х, у, z), V = у( х, у, z), ^ = w(х, у, z) - проекции вектора перемещений на координатные оси Ох, Оу, Oz произвольной точки М(х, у, z) массива. Искомые перемещения ищем в виде следующих конечных разложений:
и = ц(х у )ф1(г) + ц (x, у )ф 2 (z) +... +
П
+ ц(x, у )ф п (= Ец (x, у )фг (
i=1
V = V1 (^ у (Ю + V2 (x, у Н2 (Ю + . +
т
+ Vт (x, уНт (= ^; (x, у ; ((1)
1=1
w = wl(х у )Х1(+ ^(x, у )х 2(+ . +
р
+wp(x, у к р(= Е wk(x, у к к(z),
k=1
где ф Фп (4 (z),. ■ ■> (z),. ■ ■> 1р (г) -
заданные определенным образом безразмерные координатные функции, характеризующие распределение горизонтальных и вертикальных перемещений породного массива; координата г отсчитывается от дневной поверхности (г = 0); функции ц. (х, у), V1 (х, у), wк (х, у) являются в данном случае искомыми.
Был использован обобщенный закон Гука:
= 18 к1 , (2)
где коэффициенты С.] к1 в уравнениях (1)
суть упругие постоянные данного материала и зависят как от декартовых координат, так и от направления.
Следовательно, формулы (1) относятся к случаю анизотропного и неоднородного тела. Модели анизотропных линейно-деформируемых тел широко используются в геомеханике, например, для описания поведения сланцеватых тонкослоистых пород. Среди них наибольшее распространение получила транверсально-изотропная среда, характерной чертой которой является постоянство свойств в различных направлениях в плоскости изотропии и различие в свойствах в направлении, перпендикулярном к этой плоскости.
Относительные деформации (в рамках геометрически линейных задач) выражаются при помощи зависимостей Коши:
1
Е» = 2
(
ды.
\
1
ды
—- +
дх. дх ,
V 1 ' J
., ] = 1,2,3. (3)
Исходная система дифференциальных уравнений равновесия, описывающая НДС МГП, имеет вид
с дах дг ху , , дф.
I ^ ф. ^ +1 —^ Ф. йг -1 т +
дх
ду
+ тхгФ. |0 +1ХФ.^ = М = 1, 2,1
п;
дг da y
f—— w dz + f—— wdz-fr
1 rbr 1 1 Fh, 1 1 yz
h dx h dy h
dW i
dz
dz +
i h
+ г yz w j| 0 + J YWjdz = 0, j = 1, 2,
, m;
(4)
fdixz h dx
h dy
dy
Zkdz -Ja z—^dz + h dz
+ °гXк|° ХХ= ° к = 1 2 —> Р.
h
Уравнения (4) называются вариационными уравнениями равновесия в форме В.З.Власова. Они получены путем умножения уравнений равновесия пространственной теории упругости соответственно на координатные функции фг., уj, х к с последующим интегрированием по координате 2. Внеинтегральные члены определяются значениями напряжений на граничных плоскостях 2 = 0, 2 = И. В частности, подстилающий пласт, на котором лежит слоистый породный массив, может рассматриваться как упругое основание с заданным коэффициентом жесткости к .
При этом на нижней граничной плоскости выполняются условия
a z XkL =h = - kw( x, y, z)y,(z)
i = 1, 2,
z=h
а на свободно зависающей части (к = 0) породного массива они принимают вид
a z У
k I z=h
= 0, i = 1, 2, ..., n.
В общем случае на боковой цилиндрической поверхности рассматриваемого участка породного массива формулируются условия для напряжений или перемещений (возможен случай смешанных граничных условий).
Введение компонент вектора перемещений в виде разложений (1) приводит к необходимости оценки решения и, V, w на основе вариационного метода. В известных работах Л.В.Канторовича и В.И.Крылова, С.Г.Михлина и других авторов на многочисленных примерах показано, что даже при выборе одной-двух координатных функций погрешность решения соответствует необходимой точности. Несмотря на то, то для
308 _
вариационного метода в форме В.З.Власова общие погрешности не получены, на практике подтверждается факт: даже при небольшом числе выбранных координатных функций, учитывающих особенности процесса деформирования породного массива, погрешность метода, как правило, не превышает 1° %, что вполне допустимо для инженерных приложений в задачах геомеханики.
Отметим, что координатные функции (при условии их полноты) в разложениях (1), их вид и количество должны быть выбраны с учетом физической природы явления в исследуемой задаче (в частности, для многих задач, связанных с определением напряженного состояния слоистого массива, эти функции можно принять кусочно-линейными).
Для случая изгиба породного массива по цилиндрической поверхности трехмерная задача становится двумерной: горизонтальные и вертикальные перемещения и( х, у) и v(х, у) МГП (аналог формул (1)) примут вид
u( x y) = Z u(х)ф;(y);
i=1
V = Zv j(x )w 1(y),
j=1
(5)
где координатные функции ф1 (у),..., фп (у); у1 (у),..., ут (у) выбраны определенным образом по всей толщине породного массива, т.е. при 0 < у< И.
В этом случае система (4) после проведения всех необходимых выкладок превращается в систему п + т обыкновенных линейных или нелинейных (физически нелинейные задачи) дифференциальных уравнений второго порядка относительно искомых функций, зависящих только от одного переменного. Во многих интересных для практики случаях (в рамках линейного закона Гука) эта система интегрируется в замкнутой форме, и решение задачи удается получить в виде, пригодном для инженерных расчетов.
Разрешающие системы уравнений для случаев как линейного закона Гука, так и
n
нелинейного приведены в работе [2]. В частности, для линейного закона Гука получим симметричную систему п + т обыкновенных линейных дифференциальных уравнений второго порядка с постоянными коэффициентами. Наряду с методами численного интегрирования системы, одним из эффективных методов решения данной системы (в случае задания небольшого числа координатных функций) является метод А.Н.Крылова, позволяющий привести эту систему к одному эквивалентному ей дифференциальному уравнению порядка 2(п + т).
Сконструированные разрешающие системы уравнений [1], дополненные надлежащими граничными условиями, образуют в общем случае нелинейную краевую задачу на отыскание некоторой системы функций переменной х: на промежутке [0, /] изменения независимой переменной х отыскивается вектор-функция V (х) неизвестных, которая всюду внутри (0, I) удовлетворяет системе 2(т + п) обыкновенных дифференциальных уравнений первого порядка:
— = g (V, х), 0 < х < l. dx
(6)
Граничные условия в рассматриваемой двухточечной нелинейной краевой задаче связывают значения V0 = lim V (х) и
V = lim V (х) некоторой системой трансцендентных уравнений общего вида
D (V),V) = 0. (7)
В уравнениях (6), (7) задаваемые вектор-функции g, D считаются непрерывными вместе с первыми производными по всем указанным аргументам.
Решение нелинейной краевой задачи (6),(7) осуществляется с помощью введения числового параметра [3]. Для этого в нелинейное операторное уравнение вида T(v) = 0, символизирующее краевую задачу (6), (7), вводится числовой параметр X:
T (v, X) = 0, X 0 < А,<Л,
причем выполняются следующие равенства:
T (V0, X 0) = 0, T(v, Л) - T (v),
где элемент v0 = v(X0) является решением параметрического уравнения при X = X0 .
Решение v(X 0) строится «по точкам» с помощью итерационного процесса дискретного продолжения по параметру X, начиная с известной точки (v0,X0). В данной работе используется одношаговый итерационный процесс, у которого при каждом шаге k = 1, 2, ... в качестве начальной итерации принимается решение vk_ = v(Xк_), решение рассматриваемой задачи на предыдущем шаге имеет вид
v0 = vk _1, vk" = Ä(vF, Xk), n = 1,2,...;
\vk = lim v"n, к = 1,2,....
Конструирование итерационного процесса осуществляется на основе метода линеаризации [1], в котором оператор A(v, X) строится в неявной форме на основе записи эквивалентного операторного уравнения вида
[Tj(v, X)]v + T2(v, X) = 0.
(8)
Для дальнейшего построения итерационного процесса воспользуемся методом Ньютона - Канторовича (больше известного в западной литературе как метод Ньютона -Рафсона) [1], по которому, вводя индексы итерации п, п - 1 в уравнение (8) при каждом п = 1, 2, ... приходим к соотношению
т>П-1, х к К +
+ тК-1, хк) - [т>П-1, хк ж-1 = 0. (9)
Очевидно, что последовательность уП определяется путем решения при каждом соответствующего линейного операторного уравнения (линейной краевой задачи), а с учетом необходимого для численной реализации перехода к соответствующей дискретной задаче - из решения системы алгебраических уравнений.
Таким образом, постановка нелинейной краевой задачи по определению напряжен_ 309
но-деформированного состояния массива горных пород соответствует некоторому процессу квазистатического нагружения и автоматически включает задание некоторой программы нагружения.
Следовательно, имеем задачу на определение V (х, X):
^ = ¿(у,х,X), 0 < х < I, dx
Б(У0,У1,х,X) = 0; X0 <Х<Л. (10)
Применяя метод линеаризации к краевой задаче типа (10), переходим к уравнениям, тождественно равным уравнениям (10) и имеющим вид
^ = О, (V, X, х )У + ^ (У0, х, X);
dx
ад,VI,X)у + Д (У0,У1,X)у +
+ б (У0,У, X) I = 0,
а итерационный процесс метода линеаризации запишем в виде
dv
dx
= G4-lVn +
D0"-'v0 + Dr'v, + Dn1 = 0.
(11)
Здесь индекс к продолжения по числовому параметру X опущен.
Таким образом, система уравнений (11) является реализацией метода итерационного продолжения по числовому параметру X для получения решения нелинейной краевой задачи.
Подлежащая численному решению линейная краевая задача метода линеаризации (11) заменяется соответствующей линейной разностной задачей. Последняя (система линейных алгебраических уравнений) возникает при использовании универсальной двухточечной разностной схемы, для решения которой применяется метод матричной прогонки [2].
ЛИТЕРАТУРА
1. Беллман Р. Квазилинеаризация и нелинейные краевые задачи / Р.Беллман, Р.Калаба. М., 1968. 184 с.
2. Господариков А.П. Метод решения нелинейных задач механики горных пород при подземной разработке пластовых месторождений / Санкт-Петербургский горный ин-т. СПб, 1999. 129 с.
3. Григолюк Э.И. Проблемы нелинейного деформирования / Э.И.Григолюк, В.И.Шалашилин. М., 1988. 232 с.
REFERENCES
1. Bellman R., Kalaba R. Quasi-linearization and boundary-value problem. Moscow, 1968. 184 p.
2. Gospodarikov A.P. The method of the computation of the nonlinear problems of rock massifs / Saint Petersburg State Mining Institute. Saint Petersburg, 1999. 129 p.
3. Grigoluk E.I., Shalashilin V.I. Problems of nonlinear deformation. Moscow, 1988. 232 p.
