Прикладная математика и механика, 2022, T. 86, № 2, стр. 263-275

1. Введение. Одной из модельных задач, применяемых для исследования динамики тоннелей мелкого заложения под воздействием транспортной нагрузки (нагрузки от движущегося внутритуннельного транспорта), является задача о действии на расположенную в упругом полупространстве круговую цилиндрическую оболочку нагрузки, равномерно движущейся по внутренней поверхности оболочки вдоль ее образующей, параллельной свободной границе полупространства. В отличие от аналогичной задачи для упругого пространства, моделирующей тоннель глубокого заложения, данная задача является более сложной, так как возникает необходимость учитывать отражаемые границей полупространства волны. Поэтому количество публикаций, посвященных исследованию этой проблеме, немногочисленно и охватывает, в основном, последние годы, в частности [1–8]. В настоящей работе построена математическая модель динамики подкрепленного трехслойной обделкой тоннеля мелкого заложения под воздействием транспортной нагрузки. С использованием метода неполного разделения переменных и переразложения цилиндрических и плоских волн построено аналитическое решение задачи, на основе которого разработан алгоритм и программный комплекс на языке ФОРТРАН для расчета напряженно-деформированного состояния обделки тоннеля и породного массива с учетом скорости движения транспортной нагрузки, глубины заложения тоннеля и физико-механических свойств массива и материала оболочек. Рассмотрен дозвуковой случай, когда скорость движения меньше скоростей распространения упругих волн в массиве, в том числе релеевской для упругого полупространства. Приводятся и анализируются результаты численных экспериментов.

2. Постановка и аналитическое решение задачи. В качестве расчетной схемы подкрепленного трехслойной обделкой тоннеля мелкого заложения рассматривается бесконечно длинная круговая цилиндрическая трехслойная оболочка в линейно-упругом, однородном и изотропном полупространстве (массиве), отнесенному к неподвижным цилиндрической $r$, $\theta $, $z$ и декартовой $x$, $y$, $z$ системам координат, ось $z$ которых совпадает с осью оболочки и параллельна свободной от нагрузок горизонтальной границе полупространства (земной поверхности), ось $x$ – перпендикулярна к этой границе: $x \leqslant h$ (рис. 1). Внутренним слоем оболочки является толстостенная оболочка (заполнитель), а внешние слои (обшивка) представляют собой тонкостенные оболочки с радиусами срединных поверхностей ${{R}_{1}}$, ${{R}_{2}}$ и толщинами ${{h}_{{01}}}$, ${{h}_{{02}}}$. В силу малости толщин составляющих обшивку слоев допускается, что они контактируют с заполнителем и окружающим массивом вдоль своих срединных поверхностей. Контакт между слоями оболочки полагается жестким, а контакт между оболочкой и массивом полагается либо жестким, либо скользящим при двусторонней связи в радиальном направлении.

По внутренней поверхности оболочки в направлении ее оси $z$ с постоянной скоростью c движется нагрузка интенсивностью $P$, вид которой не меняется с течением времени (стационарная нагрузка). Скорость движения нагрузки принимается дозвуковой, т.е. меньше скоростей распространения волн сдвига в заполнителе и массиве (характерной для современных транспортных средств). Физико-механические свойства массива и заполнителя характеризуются соответственно следующими постоянными: ${{\nu }_{1}}$, ${{\mu }_{1}}$, ${{\rho }_{1}}$; ${{\nu }_{2}}$, ${{\mu }_{2}}$, ${{\rho }_{2}}$, где ${{\nu }_{k}}$ – коэффициент Пуассона, ${{\mu }_{k}}$ – модуль сдвига, ${{\rho }_{k}}$ – плотность ($k = 1,\;2$). В дальнейшем индекс $k = 1$ относится к массиву, а $k = 2$ – к заполнителю.

Поскольку рассматривается установившийся процесс, то картина деформаций стационарна по отношению к движущейся нагрузке. Поэтому можно перейти к связанной с нагрузкой подвижной декартовой ($x,y,\eta = z - ct$) или цилиндрической (r, θ, η = = $z - ct$) системе координат.

Для описания движения массива и заполнителя используются динамические уравнения теории упругости в подвижной системе координат [4]

где ${{M}_{{pk}}} = c{\text{/}}{{c}_{{pk}}}$, ${{M}_{{sk}}} = c{\text{/}}{{c}_{{sk}}}$ – числа Маха; ${{c}_{{pk}}} = \sqrt {{{\left( {{{\lambda }_{k}} + 2{{\mu }_{k}}} \right)} \mathord{\left/ {\vphantom {{\left( {{{\lambda }_{k}} + 2{{\mu }_{k}}} \right)} {{{\rho }_{k}}}}} \right. \kern-0em} {{{\rho }_{k}}}}} $, ${{c}_{{sk}}} = \sqrt {{{{{\mu }_{k}}} \mathord{\left/ {\vphantom {{{{\mu }_{k}}} {{{\rho }_{k}}}}} \right. \kern-0em} {{{\rho }_{k}}}}} $ – скорости распространения волн расширения–сжатия и сдвига в массиве и заполнителе, ${{\lambda }_{k}} = 2{{\mu }_{k}}{{\nu }_{k}}{\text{/}}(1 - 2{{\nu }_{k}})$; ${{{\mathbf{u}}}_{k}}$ – векторы смещений точек массива и заполнителя, ${{\nabla }^{2}}$ – оператор Лапласа.

Колебания слоев обшивки описываются классическими уравнениями теории тонких оболочек в подвижной системе координат [2–4]

Здесь для наружного слоя обшивки $k = 1$, для внутреннего – $k = 2$; ${{\nu }_{{0k}}}$, ${{\mu }_{{0k}}}$, ${{\rho }_{{0k}}}$ – соответственно коэффициент Пуассона, модуль сдвига и плотность материалов слоев обшивки; ${{u}_{{0\eta k}}}$, ${{u}_{{0\theta k}}}$, ${{u}_{{0rk}}}$ – перемещения точек срединных поверхностей слоев обшивки; ${{q}_{{j{{R}_{2}}}}} = {{\left. {{{\sigma }_{{rj2}}}} \right|}_{{r = {{R}_{2}}}}}$, ${{q}_{{j1}}} = {{\left. {{{\sigma }_{{rj2}}}} \right|}_{{r = {{R}_{1}}}}}$, ${{q}_{{j{{R}_{1}}}}} = {{\left. {{{\sigma }_{{rj1}}}} \right|}_{{r = {{R}_{1}}}}}$ – составляющие реакции заполнителя и массива, $j = \eta ,\theta ,r$ (при скользящем контакте оболочки с массивом ${{q}_{{\eta {{R}_{1}}}}} = {{q}_{{\theta {{R}_{1}}}}} = 0$), ${{\sigma }_{{rj1}}}$, ${{\sigma }_{{rj2}}}$ – компоненты тензоров напряжений в массиве и заполнителе, ${{q}_{{j2}}} = {{P}_{j}}\left( {\theta ,\eta } \right)$, ${{P}_{j}}\left( {\theta ,\eta } \right)$ – составляющие интенсивности подвижной нагрузки $P\left( {\theta ,\eta } \right)$, $j = \eta ,\theta ,r$.

Поскольку граница полупространства свободна от нагрузок, то при $x = h$

При различных контактных условиях оболочки с массивом граничные условия имеют вид:

– для скользящего контакта оболочки с массивом

– для жесткого контакта оболочки с массивом

где ${{u}_{{jk}}}$ – компоненты векторов ${{{\mathbf{u}}}_{k}}$, $k = 1,2$.

Векторы ${{{\mathbf{u}}}_{k}}$ можно выразить через потенциалы Ламе [1, 4]

которые, как следует из (2.1) и (2.6), удовлетворяют уравнениям

Здесь ${{{\mathbf{e}}}_{\eta }}$ – орт оси $\eta $, ${{M}_{{1k}}} = {{M}_{{pk}}}$, ${{M}_{{2k}}} = {{M}_{{3k}}} = {{M}_{{sk}}}$.

Через эти же потенциалы, используя (2.6) и закон Гука, можно выразить компоненты тензоров напряжений $\sigma _{{lmk}}^{{}}$ в массиве ($k = 1$) и заполнителе ($k = 2$) в цилиндрической $\left( {l,m = r,\theta ,\eta } \right)$ системе координат, а также $\sigma _{{lm1}}^{{}}$ в декартовой ($l,m = x,y,\eta $) системе координат.

Таким образом, для определения компонент напряженно-деформированного состояния (НДС) массива и заполнителя необходимо решить уравнения (2.7), используя граничные условия (2.3) и, в зависимости от условия контакта оболочки с массивом, (2.4) или (2.5).

Рассмотрим случай действия на оболочку синусоидальной по $\eta $ подвижной нагрузки с произвольной зависимостью от угловой координаты

где константа $\xi $ определяет период $T = 2\pi {\text{/}}\xi $ действующей нагрузки.

В установившемся состоянии зависимость всех величин от $\eta $ имеет вид (2.8), поэтому

Подставляя (2.9) в (2.7), получим

где ${{m}_{{jk}}} = {{\left( {1 - M_{{jk}}^{2}} \right)}^{{1/2}}}$, $\;{{m}_{{1k}}} = {{m}_{{pk}}}$, $\;{{m}_{{2k}}} = {{m}_{{3k}}} = {{m}_{{sk}}}$, $\nabla _{2}^{2}$ – двумерный оператор Лапласа.

Используя (2.9) можно получить выражения для перемещений $u_{{lk}}^{*}$ и напряжений $\sigma _{{lmk}}^{*}$ ($l,m = r,\theta ,\eta $) в массиве ($k = 1$) и заполнителе ($k = 2$), а также $u_{{l1}}^{*}$, $\sigma _{{lm1}}^{*}$ ($l,m = x,y$, η) в массиве от синусоидальной нагрузки как функции от ${{\Phi }_{{jk}}}$ (знак * означает, что данные компоненты найдены при действии на оболочку синусоидальной подвижной нагрузки).

При дозвуковой скорости движения нагрузки ${{M}_{{sk}}} < 1$, $\;{{m}_{{sk}}} > 0$, $k = 1,\;2$, и решения уравнений (2.11) можно представить в виде [3, 4]

где:

для массива

для заполнителя

Здесь ${{I}_{n}}({{k}_{j}}r)$, ${{K}_{n}}({{k}_{j}}r)$ – соответственно модифицированные функции Бесселя и функции Макдональда, ${{k}_{{j1}}} = \left| {{{m}_{{j1}}}\xi } \right|$, ${{k}_{{j2}}} = \left| {{{m}_{{j2}}}\xi } \right|$; ${{g}_{j}}\left( {\xi ,\zeta } \right)$, ${{a}_{{n1}}}, \ldots ,{{a}_{{n9}}}$ – неизвестные функции и коэффициенты, подлежащие определению, $j = 1,2,3$.

Как показано в [1, 2, 4], представление потенциалов для полупространства в форме (2.12) приводит к их следующим выражениям в декартовой системе координат:

где ${{f}_{j}} = \sqrt {{{\zeta }^{2}} + k_{{j1}}^{2}} $, ${{\Phi }_{{nj}}} = {{\left[ {{{\left( {\zeta + {{f}_{j}}} \right)} \mathord{\left/ {\vphantom {{\left( {\zeta + {{f}_{j}}} \right)} {{{k}_{{j1}}}}}} \right. \kern-0em} {{{k}_{{j1}}}}}} \right]}^{n}}$; $j = 1,2,3$.

Воспользуемся переписанными для $\sigma _{{xx1}}^{*}$, $\sigma _{{xy1}}^{*}$, $\sigma _{{x\eta 1}}^{*}$ граничными условиями (2.3), с учетом (2.15). Выделяя коэффициенты при ${{e}^{{iy\zeta }}}$ и приравнивая, в силу произвольности y, их нулю, получим систему трех уравнений, из которой выражаем функции ${{g}_{j}}\left( {\xi ,\zeta } \right)$ через неизвестные коэффициенты ${{a}_{{n1}}}$, ${{a}_{{n2}}}$, ${{a}_{{n3}}}$:

Вид определителя $\Delta {\text{*}}$ и алгебраических дополнений $\Delta _{{jl}}^{*}$ совпадает с аналогичными определителями для неподкрепленной полости в упругом полупространстве и определен в [2, 4]. В частности, здесь $\Delta {\text{*}}$ – это определитель Рэлея, который в данном случае имеет вид

и не обращается в ноль при любых $\zeta $, если скорость движения нагрузки меньше скорости поверхностных волн Рэлея ${{c}_{R}}$, которую условимся называть рэлеевской скоростью. В противном случае в точках $\zeta = \pm \zeta * = \pm \left| \xi \right|\sqrt {M_{R}^{2} - 1} $, ${{M}_{R}} = c{\text{/}}{{c}_{R}}$ он обращается в ноль, и интегралы в формуле (2.15) становятся расходящимися.

Ограничимся случаем $c < {{c}_{R}}$. Тогда все подынтегральные функции в (2.15) непрерывны и экспоненциально стремятся к нулю на бесконечности. С учетом (2.16), потенциалы (2.15) имеют вид

Следует отметить, что рэлеевская скорость ${{c}_{R}}$ несколько ниже скорости волн сдвига в массиве.

Используя известное при $x < h$ соотношение [1, 2]

представим ${{\Phi }_{{j1}}}$ (2.12) в цилиндрической системе координат

Подставляя в последнее выражение из (2.16) ${{g}_{j}}\left( {\xi ,\zeta } \right)$, для $c < {{c}_{R}}$ получим

где ${{b}_{{nj}}} = \sum\nolimits_{l = 1}^3 {\sum\nolimits_{m = - \infty }^\infty {{{a}_{{ml}}}A_{{nj}}^{{ml}}} } $, $A_{{nj}}^{{ml}} = \int_{ - \infty }^\infty {\frac{{\Delta _{{jl}}^{*}}}{{\Delta {\kern 1pt} *}}{{\Phi }_{{ml}}}{{\Phi }_{{nj}}}{{e}^{{ - h({{f}_{l}} + {{f}_{j}})}}}} d\zeta $.

Подставляя (2.18), при $k = 1$ и (2.12), при $k = 2$ в выражения для $u_{{lk}}^{*}$, $\sigma _{{lmk}}^{*}$ ($l,m = r,\theta ,\eta $), получим новые выражения для компонент НДС массива и заполнителя в цилиндрических координатах при $c < {{c}_{R}}$, где неизвестными будут только коэффициенты ${{a}_{{n1}}}, \ldots ,{{a}_{{n9}}}$.

Подставляя (2.10) в (2.2) и разрешая полученную для n-го члена разложения систему уравнений относительно $u_{{0n\eta k}}^{{}}$, $u_{{0n\theta k}}^{{}}$, $u_{{0nrk}}^{{}}$ можно найти их выражения.

Для определения коэффициентов ${{a}_{{n1}}}, \ldots ,{{a}_{{n9}}}$ воспользуемся, в зависимости от условия сопряжения оболочки с массивом, переписанными для $u_{{lk}}^{*}$ ($l = r,\theta ,\eta $) и $\sigma _{{r\eta 1}}^{*}$, $\sigma _{{r\theta 1}}^{*}$ граничными условиями (2.4) или (2.5). Подставляя в граничные условия соответствующие выражения и приравнивая коэффициенты рядов при ${{e}^{{in\theta }}}$, получим бесконечную систему ($n = 0, \pm 1, \pm 2, \ldots $) линейных алгебраических уравнений, для решения которой можно использовать метод редукции или более удобный для решения поставленной задачи метод последовательных отражений [2], позволяющий при каждом последовательном отражении решать систему линейных уравнений блочно-диагонального вида с матрицами размером 9 × 9 и определителями ${{\Delta }_{n}}\left( {\xi ,c} \right)$ вдоль главной диагонали.

Зная решение задачи для синусоидальной нагрузки, реакцию оболочки и окружающей ее среды на движущуюся с постоянной скоростью апериодическую (локальную) нагрузку вида $P\left( {\theta ,\xi } \right) = p\left( \theta \right)p\left( \eta \right)$ (характерного для транспортных средств) можно найти при помощи суперпозиции, используя представление нагрузки и компонент НДС массива и заполнителя в виде интегралов Фурье

Здесь $p{\kern 1pt} *{\kern 1pt} \left( \xi \right) = \int_{ - \infty }^\infty {p\left( \eta \right){{e}^{{ - i\xi \eta }}}d\eta } $.

Для вычислений перемещений и напряжений (2.19) можно использовать любой численный метод интегрирования, если определители ${{\Delta }_{n}}\left( {\xi ,c} \right)$ (n = 0, ±1, ±2, …) отличны от нуля, т.е. когда скорость движения нагрузки c меньше ее критических скоростей ${{c}_{{(n)*}}}$. Значения ${{c}_{{(n)*}}}$ определяются из дисперсионных уравнений ${{\Delta }_{n}}\left( {\xi ,c} \right) = 0$ [3] и могут оказаться меньше рэлеевской скорости. Окончательное решение будет зависеть от конкретного вида движущейся нагрузки.

Следует заметить, что исключая из постановки задачи граничные условия (2.3) и исключая из (2.12) $\Phi _{{j1}}^{{(2)}}$, получим решение аналогичной задачи для упругого пространства.

3. Численные эксперименты. Рассмотрим подкрепленный трехслойной обделкой тоннель глубиной заложения h = 6 м в породном массиве со следующими характеристиками: ${{\nu }_{1}}$ = 0.294, ${{\mu }_{1}} = \mu $ = 1.094 × 10⁸ Па, ${{\rho }_{1}}$ = 1.5 × 10³ кг/м³. Расчетные параметры для обделки: обшивка – стальные (ν₀₁ = ν₀₂ = 0.3, μ₀₁ = μ₀₂ = 8.08 × 10¹⁰ Па, ρ₀₁ = ρ₀₂ = = 7.8 × 10³ кг/м³) тонкостенные оболочки одинаковой толщины h₀₁ = h₀₂ = 0.02 м с радиусами срединных поверхностей R₁ = 3.0 м и R₂ = 2.5 м; заполнитель – толстостенная бетонная (ν₂ = 0.2, μ₂ = 1.21 × 10¹⁰ Па, ρ₂ = 2.5 × 10³ кг/м³) оболочка. Контакт между породным массивом и обделкой, а также между ее слоями полагаем жестким.

Исследуем напряженно-деформированное состояние рассматриваемого тоннеля при действии изображенной на рис. 2 цилиндрической нормальной нагрузки от внутритуннельного транспорта, движущейся по его лотку со скоростью c = 100 м/с. Нагрузка равномерно распределена вдоль оси η в интервале $\left| \eta \right| \leqslant {{l}_{0}}$ = 0.2 м.

Рассмотрим два случая нагружения: симметричное и несимметричное. В первом случае будем считать, что интенсивность нагрузки постоянна по всей поверхности ее приложения, то есть ${{P}_{r}} = q$.

Во втором случае интенсивность равномерно распределенной по угловой координате нагрузки, расположенной слева от вертикальной диаметральной плоскости тоннеля, в два раза больше интенсивности такой же нагрузки, действующей справа от этой плоскости, то есть при 90° ≤ θ ≤ 150° ${{P}_{r}} = q$; при 210° ≤ θ ≤ 270° ${{P}_{r}} = 2q$.

Параметр нагрузки $q$ (Па) подбираем таким образом, чтобы общая нагрузка по всей длине $2{{l}_{0}}$ (м) участка нагружения ${{P}_{r}} = q$ равнялась эквивалентной сосредоточенной нормальной кольцевой нагрузке интенсивностью P°° (Н/м), то есть $q = P{\text{/}}2{{l}_{0}}$.

Численные исследования соответствующих данному случаю дисперсионных уравнений показали, что в дозвуковом интервале скоростей они корней не имеют.

Введем обозначения: $u_{r}^{^\circ }$ = u_rμ/P° (м), (м), $\sigma _{{r\theta }}^{^\circ }$ = σ_rθ/P°, $\sigma _{{\theta \theta }}^{^\circ }$ = σ_θθ/P°, $\sigma _{{\eta \eta }}^{^\circ }$ = σ_ηη/P°, $u_{x}^{^\circ }$ = u_xμ/P° (м), (м), $u_{y}^{^\circ }$ = u_yμ/P° (м), $\sigma _{{yy}}^{^\circ }$ = σ_yy/P°, где P° = P°°/м (Па).

В табл. 1, 2 помещены данные расчетов НДС поперечного сечения (η = 0) тоннеля при действии симметричной и несимметричной движущейся нагрузки.

Согласно данным таблиц 1, 2 на рис. 3 построены эпюры нормальных напряжений $\sigma _{{\theta \theta }}^{^\circ }$, $\sigma _{{\eta \eta }}^{^\circ }$ и радиальных перемещений $u_{r}^{^\circ }$ на контуре r = R₂ контактирующего с внутренней оболочкой обшивки бетонного слоя при η = 0.

Из анализа результатов расчетов следует, что в случае несимметричной нагрузки, т.е. при увеличении интенсивности левой половины симметричной нагрузки в два раза, симметричный характер распределения перемещений и напряжений по контурам сечения тоннеля нарушается. При этом, в заполнителе (бетонном слое) экстремальные радиальные перемещения ${{u}_{r}}$ при θ = –120° возрастают в 3.4 раза, а экстремальные напряжения ${{\sigma }_{{r\theta }}}$, ${{\sigma }_{{\theta \theta }}}$ и ${{\sigma }_{{\eta \eta }}}$ увеличиваются соответственно в 2.0 (при θ = –140°), 1.5 (при θ = 180°) и 2.0 (при θ = –120°) раза. На контактирующей с обделкой поверхности массива при любом нагружении тоннеля экстремальные напряжения значительно ниже, чем в бетонном слое обделки.

На земной поверхности также нарушается симметрия в распределении напряжений и перемещений. Результаты расчетов НДС земной поверхности в координатной плоскости $xy$ ($\eta = 0$) при действии на тоннель симметричной и несимметричной движущейся нагрузки представлены в табл. 3.

На рис. 4 в координатной плоскости $xy$ ($\eta = 0$) показаны изменения компонент НДС земной поверхности при действии на тоннель симметричной и несимметричной подвижных нагрузок. Обозначения кривых: симметричное нагружение (2.1); несимметричное нагружение (2.2). Как следует из анализа результатов расчетов, во втором случае нагружения (несимметричное нагружение) максимальный прогиб ${{u}_{x}}$ земной поверхности в 1.5 раза больше, чем в первом случае (при симметричном нагружении), а ее наибольшее горизонтальное смещение ${{u}_{y}}$ возрастает в 15 раз. Экстремальные нормальные напряжения ${{\sigma }_{y}}$ и ${{\sigma }_{{\eta \eta }}}$ увеличиваются соответственно в 1.7 и 1.5 раза.

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

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