Научная статья на тему 'Метановыделение из отбитого угля'

Метановыделение из отбитого угля Текст научной статьи по специальности «Науки о Земле и смежные экологические науки»

CC BY
273
69
i Надоели баннеры? Вы всегда можете отключить рекламу.
Ключевые слова
УГОЛЬНЫЙ ПЛАСТ / ОЧИСТНОЙ ЗАБОЙ / КОМБАЙН / МЕТАНОВЫДЕЛЕНИЕ / ВРЕМЯ ТРАНСПОРТИРОВАНИЯ

Аннотация научной статьи по наукам о Земле и смежным экологическим наукам, автор научной работы — Полевщиков Г. Я., Шинкевич М. В.

В статье приведен алгоритм расчета метановыделения из отбитого и транспортируемого угля на основе совре-менных представлений о свойствах и состояниях углеметанового пласта. Показана возможность определения пре-дельно допустимой скорости движения добычного комбайна для исключения загазирования исходящей струи возду-ха.

i Надоели баннеры? Вы всегда можете отключить рекламу.

Похожие темы научных работ по наукам о Земле и смежным экологическим наукам , автор научной работы — Полевщиков Г. Я., Шинкевич М. В.

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.
i Надоели баннеры? Вы всегда можете отключить рекламу.

Текст научной работы на тему «Метановыделение из отбитого угля»

ФИЗИЧЕСКИЕ ПРОЦЕССЫ ГОРНОГО ПРОИЗВОДСТВА УДК 622.272.6

Г.Я. Полевщиков, М.В. Шинкевич МЕТАНОВЫДЕЛЕНИЕ ИЗ ОТБИТОГО УГЛЯ

Процессы поступления метана в горные выработки и разработка соответствующих методов прогноза относится к одной из основных в рудничной аэрогазодинамике. Уровень решений ее научно-технических задач в значительной мере определяет себестоимость добычи угля, производительность и безопасность труда. С достижением скоростей подвигания забоев в десятки метров в сутки, а производительности комбайнов десятки тонн угля в минуту, динамика изменения газовой обстановки на выемочном участке выходит за пределы оперативности традиционных методов ее прогноза, что отрицательно влияет на ритмичность работы забоя. Применение высокопроизводительной техники отработки газоносных угольных пластов, вывело пласт на уровень основного, наиболее трудно управляемого, источника выделения метана.

В последние десятилетия российскими учеными активно развиваются знания о свойствах углеметанового пласта, как твердого углегазового раствора (ТУГР), газосодержание которого связано с действующими вертикальными напряжениями. Тем самым, в расчетную модель газокинетического процесса введена непосредственная связь с геомеханическими следствиями горных работ.

В работе [1] дана эмпирическая функция распада ТУГР во времени. Из ее анализа следует, что распад основной части раствора, за вычетом остаточной газоносности, происходит за 3 часа. Этот период характеризуется высоким снижением скорости газоистощения угля во времени. Дальнейшее снижение газосодержания следует, на наш взгляд, отнести к режиму диффузии под влиянием градиента концентрации газа.

Из физических особенностей твердого углегазового раствора (ТУГР) следует, что скорость га-зовыделения из отбитого угля должна быть равна скорости распада раствора. Большее выделение газа из внутренней части куска не сможет профильтроваться через соответствующий слой угля, возрастет давление газа в его пустотности (напряжения в угольном скелете), снижающее скорость распада, а при достаточном, относительно прочности скелета, давлении в нем разовьются микротрещины и скорость фильтрации придет в соответствие со скоростью распада.

В условиях длинного очистного забоя процесс распада ТУГР имеет существенные особенности.

Прежде всего, на транспортной линии выемочного участка располагаются сотни тонн угля, отбитых в различные моменты времени. С учетом этого получено выражение для определения скорости газовыделения из отбиваемого с постоянной производительностью угля

= V - 1,34С3 (1+),

м3/мин (1)

где Ук - скорость движения добычного комбайна, м/мин; тв - вынимаемая мощность пласта, м; 1ст

- глубина вынимаемой стружки, м; Ру - плотность угля, т/м3; Гд.с.у. - газоносность отбиваемой полосы угля за вычетом газоносности угля выдаваемого за пределы выемочного участка, м3/т; (х -время транспортирования угля в пределах выемочного участка, мин; С3 -коэффициент скорости распада ТУГР.

Здесь

ГА.о.у. Го.у . Гост.техн.. , м /т, (2)

где Го,.у - газоносность пласта в снимаемой стружке, м3/т; Гост..техн. - остаточная газоносность угля выдаваемого за пределы выемочного участка, м3/т определяется по формуле [1]

Го™. = Гоу. • (С3 -1,34 • (1 + (о.- )(Сз-1)), м/т,

(3)

где (об- общее время транспортирования отбитого

угля по выемочному участку, мин.

ч, ■ /т

О 0

О 3 6 9 12 15 1,мкн

— по нормативу

— по зависимости (1)

Рис. 1. Рост относительной метанообильности отбитого угля во времени

Для проверки адекватности зависимостей (1-3) проведено сравнение с общепринятыми норма-

I, м'/мин

І II І І II

I Метонообильность фактическая

□ 1>

асчетная метанообильность I

61

205

Ь„, м

393

230 345 369

а — вентиляционные сбойки

Рис. 2. Метановыделение I из угля за время его транспортирования по конвейерному штреку протяженностью Ьв, имеющему вентиляционные сбойки с подсвежающей струей воздуха (пласт Поленов-

ский Ленинского месторождения Кузбасса)

тивными положениями [2]. На рис. 1 представлены результаты расчетов.

Из графиков видно, что, несмотря на отмеченную выше высокую начальную скорость газоис-тощения угля, зависимость (1) имеет существенно меньшие значения выделения метана, чем нормативный метод. Эта особенность потребовала ее промышленной апробации. Исследования выполнялись по данным электронных систем газового контроля, обеспечивающим измерение притока метана непосредственно из отбитого угля с квантованием в десятки секунд, на пяти выемочных участках шахт Кузбасса.

На рис. 2 представлен один из результатов сравнения расчетных и фактических значений выделения метана из угля за время его транспортировании по конвейерному штреку, от лавного перегружателя до ближайшей к очистному забою

I, м3/мин

вентиляционной сбойки с подсвежающей струей воздуха.

В период промышленной апробации установлены более значимые расхождения фактических и расчетных значений концентраций метана при отработке выбросоопасного пласта. В выявленных методом электротомографии потенциально опасных зонах из очистного забоя бурились разгрузочные скважины для предотвращения газодинамических явлений, что и привело к существенно неравномерной дегазации пласта на различных интервалах отработки выемочного столба и, соответственно, росту отклонения между расчетными и фактическими величинами выделения метана из отбитого угля.

Однако из представленных данных видно, что и в этих, очень сложных, условиях сходимость прогноза с фактическими данными существенно

Рис. 3. Фактическое выделение метана и его расчетные значения по предлагаемому и нормативному методам для условий лавы пласта 5 Чертинского месторождения Кузбасса

выше известных методов. Более того, можно сделать предварительный вывод, что результаты исследования позволяют относить ситуации превышения данных систем газового мониторинга над результатами прогноза к признакам повышенной газодинамической активности отрабатываемого интервала выемочного столба.

Сравнение расчетного метановыделения с фактическим, в лаве пласта 5 Чертинского месторождения при средней нагрузке 1000 т/сут по концентрации метана в исходящей струе воздуха из забоя приведено на рис. 4 без учета изменений коэффициента распределения воздуха.

Видно, что предлагаемый метод обеспечивает более высокую точность прогноза.

Отметим, что превышения фактических зна-

чений над расчетными в конце отработки столба обусловлены снижением коэффициента распределения воздуха и, соответственно, большим поступлениям метана на вентиляционный штрек из выработанного пространства очистного забоя.

Известно, что разработка пласта вызывает интенсивное выделение метана в очистном забое. При современных возможностях добывающих машин движение комбайна даже на небольшое расстояние с технически возможной скоростью нередко вызывает превышение предельно допустимой концентрации метана в исходящей из забоя струе воздуха.

Используя установленную зависимость, появляется возможность определять допустимую скорость движения комбайна с учетом среднесуточ-

V к, м/мин

Рис.4.. Максимально допускаемая рабочая скорость комбайна при различной среднесуточной производительности участка и коэффициенте дегазации пласта Кд = 0,3

V к, м/мин

Рис. 5. Максимально допускаемая рабочая скорость комбайна при различной среднесуточной производительности участка и коэффициенте дегазации пласта Кд = 0,5

ной нагрузки на выемочный участок.

Исходя из количества воздуха и предельно допустимой концентрации метана в исходящей струе, рассчитаем допустимое метановыделение из отбитого угля по формуле

I

о.у.д.

iGG • Im - Cд • Iпл - Сд Q

a -igg

м3/мин, (4)

где /пл - метановыделение с плоскости забоя, м3/т; Сд - предельно допустимая концентрация метана в исходящей струе участка, %; Qиcx - количество исходящего воздуха из очистного забоя, м3/мин.

Зная предельно допустимое количество метана, выделяющееся из отбиваемого и транспортируемого угля в исходящей струе воздуха, найдем допускаемую скорость движения комбайна

I,.,,

V«. =

mjpд о.у. (i -1,34с3(i + tx)C3-1)

м/мин, (5)

Результаты расчетов с различным коэффициентом предварительной дегазации пласта Кд=0,3 и Кд=0,5 представлены на рис. 4 и 5. Видим, что повышение коэффициента дегазации позволит значительно повысить допускаемую рабочую скорость комбайна. Однако такое значение коэффициента предварительной дегазации может быть достигнуто лишь с применением технологий гидроразрыва пласта.

Сравнение графиков на рис. 4 и 5 показывает, что для производительной работы очистного забоя необходимы инвестиции в разработки технологий высокоэффективной дегазации угольных пластов, которые обеспечат возврат вложений за счет повышения темпов ведения горных работ.

Таким образом, метановыделение из отбитого угля при его транспортировании соответствует газокинетическим показателям распада остаточного газосодержания твердого углегазового раствора. Установленные на этой основе зависимости обеспечивают определение основных характеристик процесса с учетом горногеологических условий и технологического режима работы забоя: глубины стружки, рабочей скорости комбайна, скорости подвигания очистного забоя, времени транспортирования угля в пределах выемочного участка. Разработанный алгоритм расчета метановыделения из отбитого угля обеспечивает определение его параметров с погрешностью не более 17 %. Появляется возможность корректировать скорость комбайна для исключения превышения максимально допустимой концентрации метана в исходящей струе воздуха.

Работа выполнена при финансовой поддержке интеграционного проекта СО РАН № 60 и грантов РФФИ №№ 10-05-90001-Бел а, 10-05-98009-рсибирьа.

СПИСОК ЛИТЕРАТУРЫ

1. Малышев, Ю.Н. Фундаментально-прикладные методы решения проблемы угольных пластов / Ю.Н. Малышев, К.Н. Трубецкой, А.Т Айруни - М.:ИАГН, 2000. - 516 с.

2. Руководство по проектированию вентиляции угольных шахт. - Макеевка-Донбасс: Ротапринт МакНИИ, 1989. - 319 с.

□ Авторы статьи:

Полевщиков Геннадий Яковлевич

- докт.техн.наук, проф., зав. лабораторией гозодинамики и геомеханики угольных месторождений Института угля СО РАН, Е-mail: Gas coal@icc.kemsc.ru.

Шинкевич Максим Валериевич

- канд. техн. наук, м. н. с лаборатории гозодинамики и геомеханики угольных месторождений Института угля СО РАН, Е-mail: Gas coal@icc.kemsc.ru.

УДК 622.235(088.8): 519.21

Д.Ю. Сирота, В.В Иванов

ЗАВИСИМОСТЬ ПОТЕНЦИАЛА ЭЛЕКТРИЧЕСКОГО ПОЛЯ ГОРНЫХ ПОРОД ОТ ПЕРИОДИЧЕСКИ ИЗМЕНЯЕМОГО ДАВЛЕНИЯ

Введение

В работах [1,2,3] Б.Г. Тарасов, В.В. Дырдин, В.В. Иванов, Д.В. Алексеев экспериментально, а затем теоретически установили взаимосвязь между градиентом среднего напряжения в горных породах Земли и напряженностью электрического поля (далее ЭП):

где Q ~ 10

3

шетки, м ; q

-> Q E = — VP . q

(i)

-29

дилатация кристаллической ре-« 1,602 • 10-19 - заряд вакансии,

e

Кл; E - напряженность ЭП, В /м; VP - градиент механических напряжений, Па /м.

Как показал Д.В. Алексеев в [3] зависимость (1) получается в результате решения системы дифференциальных уравнений (непрерывности тока, амбиполярной диффузии, Пуассона):

dn ± dt

= -div i:

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

(2)

- an1D_— _ D_n-Q_VP

i _ = _-------E- D Vn---------------------—-,(3)

kT

——

diV E =

kT

q(n + - n )

єє0

(4)

где д - заряд точечных дефектов (катионных и анионных вакансий), Кл; п ±, - соответствен-

но концентрация и дилатация для каждого типа

носителей, м-3, м3; Б ± - коэффициент самодиф-фузии дефектов положительного и отрицательно-

——

го знака, м2/с; Е - напряженность поля проводимости, В/м; Р- среднее давление в каждой точке кристалла, Па; 550 - диэлектрическая проницае-

мость кристалла, Ф/м; 1± - вектор потока точечных дефектов положительного и отрицательного знаков.

В уравнениях (2 - 4) предполагается, что напряжение, а следовательно, и его градиент не изменяются во времени.

Данная статья посвящена обобщению зависимости (1) на случай периодически измененяемого напряжения.

§1. Основные закономерности

В случае периодически изменяемого напряжения, его величина будет определяться следующим образом: Р сов(f • г), где f - частота изменения

давления. Если частота f =0 , то периодически изменяемое напряжение переходит в стационарный режим.

Запишем систему (2,3) с учетом дефектов структуры образца по 8сИойку и Френкелю относительно новых переменных N = п + п и

5Ы = п + - п - .

Для случая дефектов структуры по 8сИойку:

'(Ж); =- А2Ж + БЧ 2(Ж) + вБЧ2 N

-РАЪV2 Р сова • г) (5)

=-рА2Ж + БЧ2 N + вБЧ 2(Ж)

- А3Ч2Р соб( f • г)

для дефектов структуры по Френкелю:

где

(SN )t = -A2SN + DV 2(SN) + eDV2 N + A3V2P cos( f • t)

(N )t =-eA2SN + DV2 N + eDV 2(SN)

+ eA3V2P cos( f • t)

A2 = 2/(рбє0), A3 = 2Q/(p q2),

(6)

1/ p = q2 nD /(kT), в = [ D + - D - ]/[ D + + D ], D = (D + + D—)/2 .

Применим к системам (5, 6) преобразование Fourier по пространственным координатам x, y, z .

Учитывая, что оператор V2 преобразуется в

2 2 2,2,2 — со , где со = о1 + со2 + со3 - параметр

Fourier-преобразования, получим следующие системы дифференциальных уравнений:

Для дефектов структуры по Schottky:

'(5N)'t = —(4 + Do2)5N — (3Dm2 N +

+ PA° P c°s(/ • t) (7)

(N )t =—P(A2 + Da2)5N — Do2 N +

+ A32 P cos( f • t) для дефектов структуры по Френкелю:

'(5N)[ = —(A2 + Do2)5N —

— pDrn 2N — A3o2P cos(f • t)

(N)t =—e(A2 + Da2)5N — ()

— Do2N — pA3m2Pcos(f • t)

Преобразуем систему (7) к обыкновенному дифференциальному уравнению второго порядка.

(SN)ntt + и • (SN)[ + г/ • SN = —f вA3a>2 Psin(f •t)

(9)

где и = A2 + 2Da>2 ,

n = (1 — в2) Do2( A2 + Do2).

Отметим, что неоднородное дифференциальное уравнение (9) при частотеf =0 преобразуется к однородному уравнению, рассмотренному в работе [4].

Запишем общее решение уравнения (9) в виде суммы решения однородного уравнения и частного решения, соответствующего правой части уравнения (9).

SN(t, о) = c1 exp(k1t) + c2 exp(k2t) + (SN)* ,(10)

где

- A2 - 2Da2 _ [A22 +

+ 4Є2 Da2( A2 + a2 D)]

Л

С1, С2 - константы интегрирования.

Выражение для (SN) в случае дефектов структуры по 8сИойку будет иметь вид:

/слгч* eA3 uf22 n ч

(SN) = -----^2 P cos(f • t) —

[n —f ]2 + U f

[n — f2] f eA^Q [n — f2]2 + u2 f

(11)

P sin(f •t)

(n — f 2)(в2 — 1) Do2 — uf2

+

[n — f2]2 +u2 f2 (n — f2) + и(в2 — 1) Do2

о2 A3 P cos( f • t) +

[n — f2]2 +u2 f2

f A3rn2 P sin( f • t)

(13)

§ 2. Некоторые частные случаи

Дальнейшее решение задачи заключается в обратном преобразовании Fourier выражений (11, 13) с дальнейшей подстановкой полученных выражений в (4). Однако, рассмотрение этой задачи в общем случае весьма затруднительно. Поэтому рассмотрим несколько частных случаев, в которых можно достаточно просто установить требуемую взаимосвязь между градиентом механических напряжений и напряженностью ЭП.

Рассмотрим следующие случаи: 1) D+ * D и

2) D+ >> D , D + << D— . В зависимости от соотношений между коэффициентами самодиффу-

зии D ± изменяется значение параметра в. А

именно, если D+ * D , то в * 0; если

D+ >> D~, то в * 1; если D+ << D—, то

в * — 1.

Пусть в * 0. В этом случае выражение (11) обратится в нуль, что будет означать независимость величины потенциала ЭП от периодического давления в случае дефектов структуры по Schottky.

Выражение (13) преобразуется к виду:

(SN) =

* _ (f2 — n)Do2 —uf2

+

[n — f2]2 +u2f2 (n — f 2) — uD°2 . o2P f t

nrf+uf fA°<‘sin( f •t)

о2 A3 P cos( f • t) +

Рассмотрим случай дефектов по Френкелю и систему уравнений (8).

Обыкновенное дифференциальное уравнение второго порядка в этом случае имеет вид:

(Ж х; + и (Ж х + пж =

= (в2 -\)Б®4А3Рсов(f • г) + f А3®2Рвт(f •г)

(12)

где

/и = А2 + 2Б®2 ,

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

П = (1 -в2) Б®2( А2 + Б®2).

Отличие (12) от (9) только лишь в правой части, а следовательно, в частном решении (Ж) .

Выражение для (SN) в случае дефектов структуры по Френкелю будет иметь вид:

(Ж )* =

(14)

при этом

п = (1 - р2)Б®2(А2 + Б®2) * П®2(А2 + В®2) Сделаем некоторые порядковые оценки величин, входящих в выражение (14). Величину Б оценим из выражения [5]:

В = В, ехр(-и0 / квТ), (15)

где В0 * 2,0 -10-4, м2/с; и0 - энергия активации, Дж; кв * 1,3805 -10 23 - постоянная

Больцмана, Дж/°К; Т=293 К - температура .

Энергия активации по данным работы [5] для ЫБ, ЫС1, ЫВг варьируется в пределах и0 * (0,496 ~ 1,5) • 10-19 Дж. Тогда величина

D * [9,46 • 10—12 ~ 15,74 • 10—23], м2/с.

Так как ss0 * 8,854 -10—12 Ф/м, то при

р = 10 ~ 10 Ом^м

2

A2 =-

Р^0

2,2589 • (10 ~ 109).

Следовательно, вкладом в общую сумму А2 + Б®2 произведения Б®2 по сравнению с А2 , при условии значимости только первого десятка гармоник, можно пренебречь. Тогда выражения для г/ и ц преобразуются к виду: П ~ Б®2А2, /и ~ А2. Преобразуем (14) с учетом полученных упрощений:

Жгч * (f2 - °®2 А2)Б®2 - А2 f2 2,с

(Ж) = гп 2.—^2,2—®2А3Рсов(f • г) + [Б®2А2 - f2]2 + А22f2

(Da2 A2 — f2) — A2 Dm1 [Dm2A2 — f2]2 + A2f2 J

2 тл „2 л \тл „2 л -t2

f A302 P sin( f • t) =

= (f2 — DoA2)Doz — A, f = [ Do2 A2 — f2]2 + 4f2 °

со2 A3 P cos( f • t) +

+

— f 2

[Drn2A2 — f2]2 + Alf2J 3

f A3a>2P sin( f • t).

Дальнейшие упрощения можно произвести, сделав предположения о частоте приложения внешнего давления на образец.

Если внешнее приложение нагрузки является высокочастным, например, с периодом Т=1 с, то f2 * 39,47.

Произведение же Б®2А2 варьируется в пределах 3,557 • 10-19 ~ 2,138, поэтому им можно пренебречь.

Получим следующее:

(SN)* = f 2'I?2 , 2 a2A3Pcos(ft) +

+ -

f4 + A2f

f4 + A2f

f A3®2 P sin(ft) = -

-A

f2 + Л22

(SN) = + Л—f p cos(ft) + p sin(ft)

A2 — f A2 — f

Производя обратное преобразование Fourier и учитывая (4), получим выражение для потенциала в случае дефектов по Schottky:

x a2Л3Pcos(ft) +

-1

f2 + A

f Л3ю2P sin(ft)

p = +

—[Л2 cos( ft) + f sin(ft)]

єєп

Обратное преобразование Fourier дает:

і 2 ,2-P(X, У, z)

(І8)

A22 - f 2

SN (x, y, z) =

+ -

f2 + Л2 І

■ A3 cos(ft)V2P(x, y, z) +

f 2 A 2 fA3sin(ft)V P(x, У, z)

f + A2

Подставляя полученное выражение в (4), получим уравнение для потенциала p :

V2p = —— A cos(ft) + f sin(ft)] x

Далее, в случае дефектов по Френкелю выражение (13) для (8М) после некоторых упрощений и преобразований примет вид:

* А А

(5Ы) = —2^—®2Р со^(/і) . Следовательно, А2 /

после подстановки в (4), получим следующее выражение для потенциала:

єєп

p =

q

-cos(f )P(x, y, z) (І9)

V 2 P(x, y, z)

f 2 + A22

Простейшее решение последнего уравнения имеет вид:

р = -є x A cos(f •t) + f sin( f •t)]

A

x

єє

(І7)

f2 + Л2

x P( x, y, z)

880 А2 - У2

Выводы

1) В случае равенства коэффициентов само-диффузии положительно и отрицательно заряженных дефектов структуры получим, что в случае дефектов по 8сИойку потенциал ЭП не зависит от периодически изменяемого давления; в случае же дефектов по Френкелю указанная зависимость присутствует.

2) В случае повышенной самодиффузии того или иного знака зависимость между потенциалом и периодически изменяемом давлении присутствует как в случае дефектов по 8сИойку, так и в случае дефектов по Френкелю.

и определяет зависимость между потенциалом ЭП и давлением в случае дефектов структуры по Френкелю.

Пусть теперь в = ± 1. Тогда Г) = 0 и учитывая, что [Л ~ А2, (11) преобразуется к виду:

СПИСОК ЛИТЕРАТУРЫ

1. Тарасов, Б. Г. Применение метода электрометрии для контроля за состоянием горных выработок в условиях рудника «Октябрьский» [Текст] / Б. Г. Тарасов и др. - В сб. «Вопросы рудничной аэрологии». -Кемерово, КузПИ, 1976, вып 4. - с. 250 - 257.

2. Тарасов, Б. Г. Геотектонические процессы и аномалии квазистационарного электрического поля в земной коре [Текст] / Б. Г. Тарасов, В. В. Дырдин, В. В. Иванов // ДАН СССР 1990. Т. 312. №5. - с. 1092 - 1095.

3. Алексеев, Д. В. Баротоки в твердых телах с диффузионным механизмом проводимости [Текст] / Д. В. Алексеев // ФТТ 1991, т. 33, №10 -с. 1456 - 1476

4. Тарасов, Б. Г. Физический контроль массивов горных пород [Текст] / Б. Г. Тарасов, В. В. Дырдин, В. В. Иванов, А. Н. Фокин. - М.: Недра, 1994. - 238 с.

5. Китель, Ч. Введение в физику твердого тела [Текст] / Ч. Китель. - М.:Наука, 1978. - 789 с.

х

2

3

X

X

□ Авторы статьи:

Сирота Иванов

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

Дмитрий Юрьевич Вадим Васильевич

- ст. преп. каф. математики КузГТУ - докт. техн. наук,. проф. каф.

e-mail: sirotadm@kuzbass.net теоретической и горной механики

КузГТУ Е-mail: raen@kuzstu.ru

i Надоели баннеры? Вы всегда можете отключить рекламу.