Mathematical Models and Evaluation Software for Stress-Strain State of the Earth’s Lithosphere

Толық мәтін

Введение

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

Преимущество метода математического моделирования перед другими методами исследования геодеформационных процессов является значительным [1–2], особенно если этот метод применяется в рамках комплексного подхода при оценке риска чрезвычайных ситуаций геодинамического характера [1]. Важную роль играют экономическая эффективность, оперативность и быстроокупаемость этого метода. Кроме того, метод математического моделирования позволяет достаточно легко менять объект исследования.

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

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

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

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

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

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

1) вертикальные и горизонтальные смещения в геологической среде;

2) скорости вертикальных и горизонтальных движений в геологической среде;

3) напряжения и деформации в геологической среде, возникающие вследствие осуществления вертикальных и горизонтальных движений;

4) относительная плотность потенциальной энергии деформируемых пород геологической среды;

5) плотность геологической среды;

6) упругие модули геологической среды;

7) вязкость геологической среды;

8) ускорение свободного падения.

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

Обзор литературы

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

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

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

Детерминированные модели по своему пространственному разрешению бывают либо двумерные, либо трехмерные. Двумерные модели позволяют решать ограниченный круг практических задач, обычно связанных с оценкой геодинамического риска вдоль какого-либо профильного направления, т. е. в вертикальных разрезах земной коры. В настоящее время подобные модели практически не используются, и основное внимание в оценках геодинамического риска уделяется трехмерным моделям. Они более эффективны, поскольку позволяют охватывать уже не отдельные слои, а целые блоки земной коры. Трехмерные модели, как правило, содержат граничные условия следующего вида [3; 4]:

$\left\{\begin{array}{l}\rho g{u}_z(x,y,h)-{\sigma }_z(x,y,h)=\mathrm{0,}\hfill \\ {\sigma }_z(x,y\mathrm{,0})+({\rho }_0-\rho )g{u}_z(x,y\mathrm{,0})=\rho h_{\max }P(x,y),\hfill \\ {\tau }_{xz}(x,y,h)=\mathrm{0,}\hfill \\ {\tau }_{xz}(x,y\mathrm{,0})=\mathrm{0,}\hfill \end{array}\right.$

где τ_xz – касательные напряжения в плоскости XZ; σ_z – вертикальная составляющая нормальных напряжений; u_z – вертикальная составляющая полного вектора смещений в геологической среде; P(x, y) – величина распределенной нагрузки в точке с координатами (x, y); ρ − средняя плотность вещества геологической среды; ρ₀ − средняя плотность вещества ниже условной модельной границы приложения внешней распределенной нагрузки; g – ускорение силы тяжести; h_max – наибольшая величина отклонения модельной нижней границы от нулевого уровня; h – средняя толщина рассматриваемого модельного упругого пространства.

Аналитические выражения для компонент тензора напряжений и составляющих вектора смещений в геологической среде в обобщенном виде в таких моделях обычно записываются в следующем виде [3; 4]:

${\sigma }_{\alpha }(x,y,z)=\frac{1}{4}{\sigma }_{\alpha }^{(00)}+\frac{1}{2}{\displaystyle \sum _{m=1}^s{\sigma }_{\alpha }^{(m0)}}+\frac{1}{2}{\displaystyle \sum _{n=1}^l{\sigma }_{\alpha }^{(0n)}}+{\displaystyle \sum _{m=1}^s{\displaystyle \sum _{n=1}^l{\sigma }_{\alpha }^{(mn)}}},$

$u_{\beta }(x,y,z)=\frac{1}{4}{u}_{\beta }^{(00)}+\frac{1}{2}{\displaystyle \sum _{m=1}^s{u}_{\beta }^{(m0)}}+\frac{1}{2}{\displaystyle \sum _{n=1}^l{u}_{\beta }^{(0n)}}+{\displaystyle \sum _{m=1}^s{\displaystyle \sum _{n=1}^l{u}_{\beta }^{(mn)}}},$  (2)

где $\alpha \in \left\{x,y,z,xy,xz,yz\right\}$ , $\beta \in \left\{x,y,z\right\}$ .

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

Для количественной оценки сейсмической опасности на региональном масштабном уровне существует компьютерная программа GeoTens [4; 5].

Компьютерная программа GeoVector позволяет численно реализовать вертикально-горизонтальную математическую модель оценки геодинамической устойчивости среды распределенной природно-технической системы [4; 5].

В программе GeoVector имеется возможность учета сброса напряжений (рассеяния сейсмодеформационной энергии) посредством задания условия в окне «Сброс напряжений».

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

Недостаток вероятностных моделей заключается в том, что их использование, как и использование детерминированных моделей, на неоднородных в геологическом и тектоническом отношении территориях не всегда позволяет адекватно выполнить оценки геодинамического риска. Это относится и к территориям, где имеется существенная неопределенность в отношении информации о физико-механических характеристиках геологической среды и (или) распределенных геофизических полях-источниках геодинамических аномалий. Все это послужило веским доводом в пользу того, что были предприняты весьма удачные попытки разработки моделей, основанных на теории нечетких множеств [3; 6].

Как известно, при использовании этой теории в практических целях особенно важно грамотно сформировать базу правил системы нечеткого вывода, выбрать входные и выходные лингвистические переменные [7–10]. Поэтому в работах, посвященных построению моделей, основанных на теории нечетких множеств, в качестве входных переменных используются характеристики геологической среды, полученные как на основании анализа макросейсмических и других геофизических данных, так и информация о напряжениях, смещениях и деформациях, рассчитанных в результате применения детерминированных и вероятностных математических моделей. После этого формируется набор выходных переменных, который затем подвергается процедуре фаззификации [6].

В ряде других работ представлены модели нечеткой оценки состояния геологической среды, основанные на прямом нечетком выводе и алгоритме Мамдани, а также модели, использующие алгоритм Сугено [2].

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

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

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

Описанные программные продукты разработаны в среде программирования Borland Delphi 7.0.

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

Материалы и методы

Основным фактором, влияющим на формирование и развитие геодинамических неустойчивостей (т. е. землетрясений), как показывают результаты исследований [11–16], является гравитационное поле планеты. При математическом моделировании напряжений в статье используются данные аномального гравитационного поля в изостатической редукции.

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

В качестве такого параметра авторы предлагают использовать данные о скоростях движений на поверхности земной коры (рис. 1).

 

 

 

Рис. 1. Карта скоростей горизонтальных движений на поверхности земной коры (мм/год)

Fig. 1. The map of velocities of horizontal movements on the crustal surface (in mm/year)

 

Представленная на рис. 1 карта является авторской разработкой ученых Калифорнийского Технологического института, созданной по данным измерений, выполненных группировкой из 30 спутников, входящих в глобальную систему позиционирования¹. Эти данные успешно применяются для навигации и точного измерения геодезических координат различных объектов.

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

На основании метода Крайгинга, реализованного в программном комплексе Surfer 10.0, построено поле распределения векторов скоростей горизонтальных и вертикальных движений на поверхности земной коры.

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

Математически модель напряженно-деформированного состояния литосферы Земли представляет собой следующие системы уравнений:

1) $\left\{\begin{array}{l}{\sigma }_z(x,y,0)={\rho }_H(x,y)\Delta g_u(x,y)\chi H(x,y),\hfill \\ u_z(x,y,h)=v_z(x,y,h)t,\hfill \\ u_x(x,y,h)=v_x(x,y,h)t,\hfill \\ {\tau }_{xz}(x,y,h)=\mathrm{0,}\hfill \end{array}\right.$(3)

2) $\left\{\begin{array}{l}{\sigma }_z(x,y,0)={\rho }_H(x,y)\Delta g_u(x,y)\chi H(x,y),\hfill \\ u_z(x,y,h)=v_z(x,y,h)t,\hfill \\ u_y(x,y,h)=v_y(x,y,h)t,\hfill \\ {\tau }_{yz}(x,y,h)=\mathrm{0,}\hfill \end{array}\right.$(4)

где ${\sigma }_z(x,y,0)$  – вертикальная компонента нормальных напряжений на границе «литосфера – литосферная мантия»; $u_z(x,y,h)$  – вертикальная составляющая вектора смещений на поверхности земной коры; $u_x(x,y,h)$  – x-компонента горизонтальной составляющей вектора смещений на поверхности земной коры; $u_y(x,y,h)$  – y-компонента горизонтальной составляющей вектора смещений на поверхности земной коры; ${\tau }_{yz}(x,y,h)$  – компонента сдвиговых напряжений в плоскости XZ на поверхности земной коры; ${\tau }_{yz}(x,y,h)$  – компонента сдвиговых напряжений в плоскости YZ на поверхности земной коры; $\Delta g_u(x,y)$  – величина аномалии гравитационного поля в изостатической редукции в точке с координатами (x, y); H(x, y) – толщина литосферы в точке с координатами (x, y); $\chi$  – безразмерный коэффициент, равный 10⁶; h – толщина модельного слоя геологической среды (литосферы); t – временной интервал (1 год); ρ_H (x, y) – плотность геологической среды в точке с координатами (x, y), вычисляемая по соотношению ${\rho }_H(x,y)={\displaystyle \underset{0}{\overset{H}{\int }}\rho (x,y,z)dz}$ .

Для расчета компонент тензора напряжений и составляющих вектора смещений в литосфере достаточно воспользоваться одним из соотношений, поскольку выражения для компонент сдвиговых напряжений и компонент горизонтальной составляющей вектора смещений формально будут одинаковыми, различаясь только координатами x и y. Будем строить выражения для компонент тензора напряжений и составляющих вектора смещений, основываясь на системе уравнений (1) и следующих соотношениях (для k-ой гармоники), общепринятых в литературе, посвященной оценке геодинамического риска [5; 6; 15]:

 

$\left\{\begin{array}{l}{\sigma }_z(x,y,z)=k^2\left\{\left[\frac{\mu }{\lambda +\mu }C-k(B+Dz)\right]shkz+\left.\left[\frac{\mu }{\lambda +\mu }D-k(A+Cz)\right]chkz\right\}\cos k_{x}x\cos k_{y}y,\right.\hfill \\ u_z(x,y,z)=\frac{k}{2\mu }\left\{\left[\frac{2\mu }{\lambda +\mu }D-k(A+Cz)\right]shkz\left.+\left[\frac{2\mu }{\lambda +\mu }C-k(B+Dz\right]chkz\right\}\cos k_{x}x\cos k_{y}y,\right.\hfill \\ u_x(x,y,z)=\frac{k_x}{2\mu }\left\{\left[k(B+Dz)+C\right]shkz\left.+\left[k(A+Cz)+D\right]chkz\right\}\text{sin}{k}_x\cos k_{y}y,\right.\hfill \\ {\tau }_{xz}(x,y,z)=k{k}_x\left\{\left[k(A+Cz)+\frac{\lambda }{\lambda +\mu }D\right]shkz\left.+\left[k(B+Dz)+\frac{\lambda }{\lambda +\mu }C\right]chkz\right\}\text{sin}{k}_{x}x\cos k_{y}y,\right.\hfill \end{array}\right.$(5)

где A, B, C, D – неизвестные коэффициенты, определяемые из условий математической модели (3).

Подставляя уравнения (5) в соответствующие соотношения системы (3) и выполнив необходимые преобразования, получим следующую систему уравнений:

 

 $\left\{\begin{array}{l}{k}^2\left(\frac{\mu }{\lambda +\mu }D-kA\right)={\rho }_H(x,y)\Delta g_u(x,y)\chi H(x,y),\hfill \\ \frac{k}{2\mu }\left\{\left[\frac{2\mu }{\lambda +\mu }D-k(A+Ch)\right]shkh\left.+\left[\frac{2\mu }{\lambda +\mu }C-k(B+Dh)\right]chkh\right\}=v_z(x,y,h)t,\right.\hfill \\ \frac{k_x}{2\mu }\left\{\left[k(B+Dh)+C\right]shkh+\left.\left[k(A+Ch)+D\right]chkh\right\}=v_x(x,y,h)t,\right.\hfill \\ \left[k(A+Ch)+\frac{\lambda }{\lambda +\mu }D\right]shkh+\left[k(B+Dh)+\frac{\lambda }{\lambda +\mu }C\right]chkh=0.\hfill \end{array}\right.$(6)

Выражая из первого уравнения системы (6) неизвестный коэффициент A, получим:

$\left\{\begin{array}{l}A=\frac{\mu }{k(\lambda +\mu )}D-\frac{{\rho }_H(x,y)\Delta g_u(x,y)\chi H(x,y)}{k^3},\hfill \\ \left[\frac{2\mu }{\lambda +\mu }D-kA-Ckh\right]shkh+\left[\frac{2\mu }{\lambda +\mu }C-kB-Dkh\right]chkh=\frac{2\mu v_z(x,y,h)t}{k},\hfill \\ (kB+Dkh+C)shkh+(kA+Ckh+D)chkh=\frac{2\mu v_x(x,y,h)t}{k_x},\hfill \\ \left(kA+Ckh+\frac{\lambda }{\lambda +\mu }D\right)shkh+\left(kB+Dkh+\frac{\lambda }{\lambda +\mu }C\right)chkh=0.\hfill \end{array}\right.$(7)

Подставляя выражение для коэффициента A в остальные уравнения системы (7) и вводя для удобства расчетов обозначение $P_g(x,y)={\rho }_H(x,y)\Delta g_u(x,y)\chi H(x,y)$  , получим следующую систему уравнений:

 

$\left\{\begin{array}{l}\left(\frac{\mu }{\lambda +\mu }D+\frac{P_g(x,y)}{k^2}-Ckh\right)shkh+\left(\frac{2\mu }{\lambda +\mu }C-kB-Dkh\right)chkh=\frac{2\mu v_z(x,y,h)t}{k},\hfill \\ (kB+Dkh+C)shkh+\left[D\left(\frac{\mu }{\lambda +\mu }+1\right)-\frac{P_g(x,y)}{k^2}+Ckh\right]chkh=\frac{2\mu v_x(x,y,h)t}{k_x},\hfill \\ \left(D-\frac{P_g(x,y)}{k^2}+Ckh\right)shkh+\left(kB+Dkh+\frac{\lambda }{\lambda +\mu }C\right)chkh=0.\hfill \end{array}\right.$(8)

Выразим из третьего уравнения системы (8) коэффициент B, а именно:

$B=\frac{P_g(x,y)}{k^3}th(kh)-[kh+th(kh)]D-\left(\frac{\lambda }{\lambda +\mu }+khth(kh)\right)C.$(9)

Подставив выражение для коэффициента B в первое и второе уравнения системы (9) и выполнив необходимые преобразования, получим систему, состоящую из двух уравнений:

 

$\left\{\begin{array}{l}C+Dth(kh)=\frac{2\mu (\lambda +\mu )v_z(x,y,h)}{k(\lambda +2\mu )},\\ \left(\frac{\lambda +2\mu }{\lambda +\mu }ch(kh)-th(kh)\cdot sh(kh)\right)D+(kh\cdot ch(kh)+\\ +\frac{\mu }{\lambda +\mu }sh(kh)\left.-kh\cdot th(kh)\cdot sh(kh)\right)(kh\cdot ch(kh)+\\ +\frac{\mu }{\lambda +\mu }sh(kh)-kh\cdot th(kh)\cdot sh(kh))C=\\ =\frac{2\mu v_x(x,y,h)}{k_x}+\frac{P_g(x,y)}{k^2}\times [ch(kh)-th(kh)\cdot sh(kh)].\end{array}\right.$(10)

Получим теперь из первого уравнения системы (10) соотношение для коэффициента C:

$C=\frac{2\mu (\lambda +\mu )v_z(x,y,h)}{k(\lambda +2\mu )}-D\cdot th(kh).$(11)

С учетом соотношения (11) второе уравнение системы (10) теперь будет содержать только один неизвестный коэффициент D. В таком случае приведем теперь окончательное выражение для нахождения коэффициента D:

$\begin{array}{l}D=\left(\frac{2\mu v_x(x,y,h)}{k_x}\right.+\frac{P_g(x,y)}{k^2}[ch(kh)-th(kh)\cdot sh(kh)]\left.-\frac{2\mu (\lambda +\mu )v_z(x,y,h)}{k(\lambda +2\mu )}\right)\times \\ \times \left(kh\cdot ch(kh)+\frac{\mu }{\lambda +\mu }sh(kh)-\right.\left(-kh\cdot th(kh)\cdot sh(kh)\right)/\frac{\lambda +2\mu }{\lambda +\mu }\times [ch(kh)-th(kh)\cdot sh(kh)]+kh\cdot sh(kh)\left(\frac{sh(kh)}{ch(kh)}\cdot th(kh)-1\right).\end{array}$(12)

Таким образом, используя соотношения (9), (11), (12) и первое уравнение из системы (7), можно рассчитать компоненты ${\sigma }_z(x,y,z)$  , ${\tau }_{xz}(x,y,z)$ тензора напряжений и составляющие $u_z(x,y,z)$  и $u_x(x,y,z)$  вектора смещений в литосфере, представленные в системе (3).

Приведем также выражения для $u_y(x,y,z)$  – составляющей вектора смещений и остальных компонент тензора напряжений [4; 5; 15; 16]:

$u_y(x,y,z)=\frac{k_y}{2\mu }\left\{\left[k(B+Dz)+C\right]shkz\right.\left.+\left[k(A+Cz)+D\right]chkz\right\}\cos k_{x}x\sin k_{y}y,(13.1)$

$\begin{array}{l}{\sigma }_x(x,y,z)=\left\{\left[k{k}_x^2(B+Dz)+C\left(k_x^2+\frac{\lambda k^2}{\lambda +\mu }\right)\right]shkz+\right.\left.\left[k{k}_x^2(A+Cz)+D\left(k_x^2+\frac{\lambda k^2}{\lambda +\mu }\right)\right]chkz\right\}\times \cos k_{x}x\cos k_{y}y,\text{(13}\text{.2)}\\ \end{array}$${\sigma }_y(x,y,z)=\left\{\left[k{k}_y^2(B+Dz)+C\left(k_y^2+\frac{\lambda k^2}{\lambda +\mu }\right)\right]shkz+\right.\left.\left[k{k}_y^2(A+Cz)+D\left(k_y^2+\frac{\lambda k^2}{\lambda +\mu }\right)\right]chkz\right\}\times \cos k_{x}x\cos k_{y}y,\text{(13}\text{.3)}$

$\begin{array}{r}{\tau }_{yz}(x,y,z)=k{k}_y\left\{\left[k(A+Cz)+\frac{\lambda }{\lambda +\mu }D\right]shkz+\left.\left[k(B+Dz)+\frac{\lambda }{\lambda +\mu }C\right]chkz\right\}\times \cos k_{x}x\sin k_{y}y,(13.4)\right.\\ \end{array}$$\begin{array}{r}{\tau }_{xy}(x,y,z)=-k_x{k}_y\left\{\left[k(B+Dz)+C\right]shkz+\left.\left[k(A+Cz)+D\right]chkz\right\}\times \sin k_{x}x\sin k_{y}y.(13.5)\right.\\ \end{array}$

Используя соотношения (13.1)–(13.5) и выражения системы (5), мы можем восстановить во всем объеме литосферы поля напряжений и смещений, соответствующие упругой модели, учитывающей пространственные вариации упругих модулей, поскольку величины λ и µ входят в выражения для компонент тензора напряжений и составляющих вектора смещений.

Однако нас интересуeт не теоретическое пространственное распределение напряжений, смещений и деформаций в литосфере Земли, а распределение, соответствующее реальной геологической среде, состояние которой в значительной степени зависит от пространственного распределения температуры в литосфере и, как следствие, от пространственной вариативности вязкости геологической среды.

Для оценки пространственной вариативности вязкости геологической среды в зависимости от пространственного распределения температуры в литосфере Земли была использована методика определения величины вязкости, основанная на достаточно хорошо известном соотношении, связывающим вязкость, температуру, упругие модули и энергию активации² [1]:

 $\eta (z)=A_k\exp \left[\frac{U_a}{RT}\exp \left(\frac{2.5P}{\lambda +2/3\mu }\right)\right],$(14)

где R – газовая постоянная, равная 8,31 $\frac{ДЖ}{моль\cdot °К}$  ; T – температура; P – давление; U_a – энергия активации; A_k – некоторый коэффициент пропорциональности.

Численные значения всех величин, входящих в данное соотношение, были либо определены нами ранее, либо заимствованы из работ других авторов³; ⁴ [1–4; 17].

Согласно упруго-вязкой модели Максвелла, общая деформация складывается из упругой и вязкой составляющих⁵; ⁶ [17; 18]:

$\gamma =\frac{P}{E}+\frac{P}{\eta }\cdot \tau$ , (15)

где $\gamma$  – деформация; P – напряжение; $E=\frac{(3\lambda +2\mu )\mu }{\lambda +\mu }$  – упругий модуль (модуль Юнга); $\eta$  – вязкость; $\tau$  – время.

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

Результаты исследования

На основании соотношений (5), (13.1)–(13.5) и с учетом выражения (14) произведен расчет напряженно-деформированного состояния литосферы Земли на различных глубинных уровнях по упругой и упруго-вязким моделям.

Численные расчеты напряжений, смещений и деформаций в литосфере Земли были выполнены с помощью компьютерной программы GeoTens (рис. 2).

 

 

Рис. 2. Рабочее окно программы GeoTens

 

 

 

Fig. 2. Working window of the GeoTens program

 

 

После запуска программа запрашивает имя файла изостатических аномалий, которое нужно обязательно ввести, поскольку именно в этом файле в формате системы GEOGRID размещена информация о внешней вертикальной нагрузке, действующей на среду РПТС. Также необходимо ввести глубины (в км), на которых расположены условная граница возмущений и расчетная плоскость. После ввода с клавиатуры длин (в км) долготного и широтного профилей необходимо нажать на кнопку «Прочитать». При этом происходит считывание данных из файла изостатических аномалий, и программа готова к выполнению расчетов. Расчеты производятся после нажатия на кнопку «Рассчитать».

После этого необходимо записать полученные результаты во внешние файлы (по каждой компоненте тензора напряжений и составляющей вектора смещений в геосреде записывается отдельный файл), нажав на кнопку «Записать». Далее можно продолжить расчеты по какой-либо другой территории (нажав на кнопку «Продолжить»), либо выйти из программы (кнопка «Выход»).

Полученные результаты исследований с помощью программного комплекса Surfer в виде изолиний выводятся на карты и представлены на рис. 3; 4.

 

 

 

Рис. 3. Эквипотенциальное распределение сдвиговых упругих деформаций в литосфере
на глубине 10 км (сечение изолиний 2·10‒8), осредненных за 1 год. Кружками красного цвета
изображено распределение эпицентров сейсмических событий на глубине 10 км
за период времени 1914–2014 гг. (карта построена авторами)

Fig. 3. The equipotential distribution of the shear elastic deformations in the lithosphere at a depth
of 10 km (the section of the isolines is 2·10‒8) averaged over 1 year. Mugs of red color depict
the distribution of epicenters of seismic events at a depth of 10 km for the period of 1914–2014
(the map was built by the authors)

 

 

 

 

Рис. 4. Эквипотенциальное распределение сдвиговых упруго-вязких деформаций
в литосфере на глубине 10 км (сечение изолиний 2·10‒8), осредненных за 1 год
(карта построена авторами)

Fig. 4. The equipotential distribution of the shear elastic-viscous deformations in the lithosphere
at a depth of 10 km (the section of the isolines is 2·10‒8) averaged over 1 year
(the map was built by the authors)

 

 

Обсуждение и заключение

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

 

¹           GPS Time Series. URL: https://sideshow.jpl.nasa.gov/post/series.html

²           Минаев В. А., Фаддеев А. О. Оценки геоэкологических рисков. Моделирование безопасности туристско-рекреационных территорий. М. : Финансы и статистика, изд. дом ИНФРА-М, 2009. 370 с.

³           Математическое моделирование геодинамических рисков: оценки и перспективы / В. А. Минаев [и др.] ; под ред. В. А. Минаева, А. О. Фаддеева, К. М. Бондаря. Хабаровск : ДВЮИ МВД России, 2015. 211 с.

⁴           Минаев В. А., Фаддеев А. О., Кузьменко Н. А. Моделирование и оценка геодинамических рисков. М. : «РТСофт»–«Космоскоп», 2017. 256 с.

⁵           Математическое моделирование геодинамических рисков: оценки и перспективы.

⁶           Минаев В. А., Фаддеев А. О., Кузьменко Н. А. Моделирование и оценка геодинамических рисков.

 

Авторлар туралы

Alexander Faddeev

The Academy of the FPS of Russia

Email: fao1@mail.ru
ORCID iD: 0000-0002-7259-1693
ResearcherId: I-3739-2018

Professor, Chair of Mathematics and Information Technologies,D.Sc. (Engineering), Associate Professor

Ресей, 1 Sennaya St., Ryazan 390000

Svetlana Pavlova

The Academy of the FPS of Russia

Хат алмасуға жауапты Автор.
Email: pavlovhome23@rambler.ru
ORCID iD: 0000-0001-8634-9163
ResearcherId: I-2969-2018

Associate Professor, Chair of Mathematics and Information Technologies, Ph.D. (Engineering)

Ресей, 1 Sennaya St., Ryazan 390000

Tatiana Nevdakh

The Academy of the FPS of Russia

Email: ntmterrano@mail.ru
ORCID iD: 0000-0002-4799-0748
ResearcherId: I-3611-2018

Lecturer, Chair of Mathematics and Information Technologies

Ресей, 1 Sennaya St., Ryazan 390000

Әдебиет тізімі

Қосымша файлдар