CC BY
24
Вычислительные технологии Том 20, № 4, 2015
Математическое моделирование сейсмических колебаний системы, состоящей из водохранилища, плотины, фундаментного блока и подфундаментных слоев грунта
И. Д. МУЗАЕВ1'2'*, Н.И. МУЗАЕВ1
1 Геофизический институт Владикавказского НЦ РАН, Россия
2Финансовый университет при Правительстве РФ. Владикавказский филиал, Россия *Контактный e-mail: muzaevid@mail.ru
Разработана математическая модель сейсмических колебаний системы, состоящей из водохранилища, плотины, фундаментного блока и подфундаментных слоев грунтового массива. Модель представляет собой контактную краевую задачу для дифференциального уравнения потенциального движения воды в водохранилище, дифференциального уравнения колебания консольной балки (балочная модель плотины), а также для дифференциальных уравнений поперечных (сдвиговых) колебаний двух слоев подфундаментного массива грунта. Эти уравнения взаимосвязаны через граничные условия на контактных поверхностях. Краевая задача решена аналитически. Полученные расчетные зависимости и результаты выполненных по ним вычислительных экспериментов позволяют на стадии изыскания и проектирования рассматриваемого гидротехнического сооружения подбирать место строительства и габаритные размеры сооружения, которые минимизируют динамическое воздействие падающей на систему сейсмической волны.
Ключевые слова: сейсмические колебания системы, подфундаментный слой грунтового массива, контактная краевая задача, сдвиговые волны, нейтральная ось, потенциальное движение воды.
В прямоугольной системе координат Охух сейсмические колебания системы, состоящей из узкоглубокого непризматического водохранилища, плотины, фундаментного блока и подфундаментного двухслойного массива грунта моделируются следующей контактной краевой задачей математической физики (рис. 1) [1-8]:
д2ш д2ш 1 дВдш 1 дВдш
—- +--- +----- +----- = 0 (1)
дх2 дг2 В дх дх В дг дг
0 < х < то, В (х,г) В(х,г)
пРи ^--< У < ,
,0 < г < Ь, д2 ( д2У\ д2У
д.V2 \Е1~сЫ + = В-0,")Р(0,Z,t), 0 < * < Ц -2)
ИВТ СО РАН, 2015
а
Рис. 1. Схематический чертеж гидротехнического сооружения: а — продольное сечение; б — конфигурация в плане; в — поперечное сечение
— а
д2щ 2 д2щ д 2
дг2
— а
2 __
1 дг2
2 д2 щ
дг2
= 0 при - Н < < 0, = 0 при 0 < гх < ж,
Р(0,z,í) = -ро V, ;, 0 < г < Ь,
дг дУ
-—-, (р ограничена при х ^ ж,
х=0
дер
V(г, ¿) = щ (гх ,г)
г=Ь
= 0,
х=0
дУ (г, г)
др дг
= 0,
^=ь
г!=-Н
дг
0 д2у д3у 0 0
= 0, —— = -г-г- = 0 при г = 0,
г=Ь
дг2 дг3
т-
д 2 щ
дг2
- Ох во
дщ
X! = -Н
дгл
=-Н
- ^ (Е1—1
дг \ дг2 )
= 0,
^=ь
дщ дщ ь1=&2^~ , = Щ при = 0,
П2,
дщ дщ
ограничены гх ^ ж.
(3)
(4)
(5)
(6)
(7)
(8)
(9) (10) (11)
дЪ ' дг
Легко заметить, что краевая задача (1)—(11) записана в двух параллельных координатных системах. В первой системе начало берется на гребне плотины, в другой начало
выбрано в плоскости раздела слоев. Координаты г и связаны между собой зависимостью
Такой подход значительно облегчает ход решения поставленной контактной краевой задачи.
В выражениях (1)—(12) приняты следующие обозначения и предположения: Ь — время; г и — вертикальные координаты, отсчитываемые от точек О и 01 соответственно; — потенциал скорости безвихревого движения воды в непризматическом водохранилище; и1(г1,Ь) и и2(х1,1) — поперечные сдвиговые перемещения в верхнем ограниченном и нижнем полуограниченном слоях грунта; С1,р1,С2,р2 — модули упругости и плотности в верхнем и нижнем слоях соответственно; V(г,Ь) — перемещение центральной оси сооружения при его поперечном колебании; т — масса фундаментного монолитного массива; Е1 — жесткость высотного сооружения; Е — модуль упругости материала; I — момент инерции поперечного сечения сооружения относительно нейтральной оси; в — площадь поперечного сечения сооружения; во — площадь планового сечения фундаментного массива; Н — мощность (толщина) верхнего слоя грунта; р — плотность материала сооружения; а\ и а2 — скорости распространения поперечной сдвиговой волны в верхнем и нижнем слоях грунта
Выражение (1) представляет собой дифференциальное уравнение двумерного безвихревого движения воды в узком и глубоком полуограниченном водохранилище непризматической конфигурации с глубиной Ь и шириной В(х,х), зависящее от продольной х и вертикальной г координат. Выражение (2) представляет собой дифференциальное уравнение поперечных упругих колебаний балки-плотины с учетом действия на ее напорную грань гидродинамического давления воды. Выражения (3) и (4) — дифференциальные уравнения поперечных сдвиговых колебаний слоев грунта под фундаментом. Первое из граничных условий (6) выражает равенство скоростей горизонтальных перемещений сооружений и воды у напорной грани сооружения, второе — ограниченность потенциала скоростей воды в бесконечности. Первое из граничных условий (7) выражает пренебрежение волнообразованием на свободной поверхности воды в водохранилище, второе — непроницаемость дна водохранилища. Первое из граничных условий (8) выражает равенства горизонтальных перемещений подошвы плотины, фундаментного блока и верхней грани верхнего слоя грунта, второе — жесткую заделку подошвы плотины с фундаментным блоком, третье — равенство нулю изгибающего момента и поперечной силы на гребне плотины. Условие (9) представляется вторым законом Ньютона для фундаментного монолитного блока, на который кроме силы инерции действуют касательная сила со стороны массива грунта и поперечная сила со стороны подошвы плотины. Граничные условия (10) выражают равенство перемещений и касательных напряжений на поверхности контакта слоев грунта с разными физико-механическими характеристиками.
Непризматическую конфигурацию приплотинного участка водохранилища аппроксимируем через экспоненциальную функцию следующего вида:
г = ¿1 + Ь + Н.
(12)
(13)
В(х,г) = Вое31Хе-
(14)
Постоянная ^ характеризует непризматическую конфигурацию водохранилища в плане, — сужение водохранилища по глубине, В0 — длина плотины на отметке уреза воды. В дальнейшем 51 и будем называть формпараметрами каньона на при-плотинном участке водохранилища.
На первом этапе моделирования решается гидродинамическая часть поставленной контактной краевой задачи.
При гармоническом колебании всей системы можно применить следующие подстановки:
ф, г, г) = ф(х, г)е<ш*е42*, V(г, ¿) = У/(г)еш\ (15)
С учетом (14) и (15) краевая задача относительно потенциала скорости сводится к виду
дф дх
д2ф + д2ф + ^ дф s2 ^ = 0 (16)
дх2 dz2 дх 4 '
s2
= iuW(z)e-z, ф ^ 0 при х ^ ж, (17)
х=0
ф
~ 0 (t +1 *)
= 0. (18)
х=Ь
Неизвестную функцию гф будем искать в виде следующего тригонометрического ряда по переменной г [6]:
те
ф(х,г) = ^21фп (x)sin^n— , (19)
п=1
где — положительные корни трансцендентного уравнения
— cos + — sin = 0
— 2
или
2 . .
tg =--п = 1, 2, 3,... (20)
2 —
Легко проверить, что если — корень характеристического уравнения (20), то
выражение (19) автоматически удовлетворяет граничным условиям (18).
s2
Разложим формально неизвестную функцию W(z)е~в тригонометрический ряд по ортогональной системе функции в промежутке (0; —) [6]:
— , sin ^2 — , ..., sin^n — ,
(21)
s2 X—^ 6
W (z) е-= 2_^an sin —, (22)
ra=1
J W(z)e 2 * sin ^n^-dz
an = --z-. (23)
sin2 ^n^rdz
L
Подставим разложения (19) и (22) в дифференциальное уравнение (16) и в граничные условия (17) и, приравнивая коэффициенты при синусах в левых и правых частях, получим
(Р фп + ¿фп (ДI + 32\,
+ 811ъ \Т2 + = 0, (24)
¿ф„
dx
= гшап, фп ^ 0 при х ^ ж. (25)
х=0
Решение дифференциального уравнения (24) с граничными условиями (25) имеет вид
фп(х) = - г-Ш е_Кх, (26)
= £ + (27)
2 V L2 4 4 Для потенциала скорости получается следующая формула:
<fl(x, z, t) = ( — ше^г V ^е-Кх sin pj- ) ешК (28)
V t=iXn L)
Гидродинамическое давление на напорной грани плотины вычисляется по формуле
dip(x, z, t)
Р(z, t) = -Ро-
<х
= рош2е* - V ^ sin ^е^. (29)
, Лп L
х=0 п=1
На втором этапе решения поставленной краевой задачи (1)—(12) рассматривается дифференциальное уравнение (2), моделирующее поперечные колебания подвижной консольной балки (балочная модель плотины). В этом уравнении момент инерции I и площадь вычисляются по формулам
к3
I = В(0, г) —, в = В(0, г)к, (30)
где В(0, г) — ширина каньона в створе плотины; к — толщина плотины.
В(0, г) = Вое _32". (31)
Подставив выражения (15), (23), (30) и (31) в уравнение (2), получим
¿3Ш 2(12Ш 12 р 2ттг - 2 + ^ - =
* ь
24Рош2_е *, ^ | 81П ^ 1 ч*
Ё W(Z)e-f* sin Un jdz I , (32)
Eh3L ¿—s \ Pn^nJ " L
n=1 * о
Pn = 1 - ^. (33)
2 ¡Лп
Решения волновых дифференциальных уравнений (3) и (4) представим в виде сумм прямых и обратных волн:
I) = Ле^) +в1еК'"^), (34)
и2 (г) = Л ^2) + в ^2), (35)
где А\, В\, А2, В — амплитуды падающих и отраженных поперечных (сдвиговых) волн в верхнем и нижнем слоях соответственно.
Для неизвестной функции Ш(г) получено обыкновенное интегродифференциальное уравнение типа уравнения Фредгольма (32). Путь аналитического решения интегро-дифференциального уравнения (32) является весьма кропотливым и связан с математическими трудностями. Они сравнительно легко преодолеваются, если принять следующее упрощающее предположение.
В правой части интегродифференциального уравнения (33) искомая функция содержится под интегралом. Сама функция V(г, ¿) = Ш(г)ешЬ выражает процесс сейсмического колебания плотины. В первом приближении будем считать, что она не зависит от координаты г. Механически это допущение означает, что плотина совершает колебательное движение как монолитное тело. Такое допущение принимается в нескольких аналогичных работах, например в [1, 2, 7]:
V(х,г) « V(ь,г) = иц(-н,г). (36)
При таком предположении интегродифференциальное уравнение (33) вырождается в обыкновенное дифференциальное уравнение с заранее неизвестной правой частью:
№ № 12 р 2
X
= -те** Ue""% + Вгегш^ V -L-In sin^, (37)
V / —* л„ р„ ь
п=1
£2^ . z
In = J е 2 z sin ^n—dz. (38)
о
Общее решение дифференциального уравнения (37) имеет следующий вид:
s2
W(z) = (C1 sh Xiz + C2 ch Xiz + C3 sin \2z + c4 cos \2z)e^z+
+W*(z) (Aie% + Biегш, (39)
где W* (z) — его частное решение без множителя (^А1е ш ai + В1 еш a.
Постоянные C1, C2, C3 и C4 определяются в процессе удовлетворения граничным условиям. Всего имеются семь граничных условий. Два — на гребне плотины
-=/1?4 *=/У1Ч■ (40)
L
Три условия — на месте стыковки сооружения, фундаментного блока и верхнего слоя грунта. Два условия — на месте контакта верхнего и нижнего слоев. Произвольных постоянных — восемь: С\, С2, С3, С4, А, В\, А2 и В2.
Однако в дальнейшем заданной величиной будем считать амплитуду падающей на нижний слой сейсмической волны, т.е. А2. Тогда количество неизвестных постоянных и количество линейных алгебраических уравнений будет 7. Следовательно, задача замыкается. После определения постоянных их значения подставляются в выражение (39). Затем из выражения V(z, t) = W(z)егш* выделяется действительная часть и с ее помощью определяется амплитуда колебаний гребня плотины. Она получается прямо пропорциональной амплитуде падающей на нижний слой сейсмической волны, т. е. величине А2.
С целью сокращения математических выкладок рассмотрим упрощенный вариант решаемой контактной краевой задачи. Принимаются предположения о том, что припло-тинная область каньона имеет призматическую конфигурацию, т.е. si = s2 = 0. При этом для амплитуды колебания гребня плотины получаются следующие зависимости:
у = 2 GjS— , (41)
ai у/К2у2 + у2
ух = [(тш2 + El-fp'''(0) cos AL)(1 + ch XL cos AL) + EIA3R] cos —+
ü\
+G\s—(1 + ch AL cos AL) sin —, (42)
a\ a\
. и
у* = [(тш2 + EI^p'''(0) cos AL)(1 + ch AL cos AL) + EIA3R] sin--
ü\
-Gvs — (1 + ch AL cos AL) cos —, (43)
ai a\
z = (1 + ф)(еЪ AL + cos AL) - T* (sh AL + sin AL), (44)
ф = '''(0) sin AL + (L), T* = '''(0) cos AL, (45)
A3 A3
R =(1 + ф)^ -T* sh AL sin AL, R = sh AL cos AL + ch AL sin AL, (46)
7 = 24P0— A = 4 12p0 — К = ^ (47)
7 LEh3 , A V Eh2 , К P2ü2 , ( )
2ж
ш = t = 2l{V,
«'"(«) = - £ ^. *= £ «у- (48)
Легко заметить, что при т = 0 имеет следующие значения:
(2п - 1)тт Д = ап =-—-, п = 1, 2, 3,...
Результаты вычислительных экспериментов показывают, что интенсивность сейсмического колебания плотины существенно зависит от габаритных размеров сооружения, физико-механических характеристик и мощностей подфундаментных слоев грунта, а также от амплитуды и частоты (периода) падающей на систему сейсмической волны.
Реализацию полученных расчетных формул на ЭВМ и автоматизацию вычислительных процедур целесообразно выполнять в такой последовательности.
— В базу данных вводят идентификаторы входных параметров системы, числовые значения которых устанавливаются сейсмическим микрорайонированием местности, а также в соответствии со строительными нормами и правилами на сейсмоактивных территориях.
— По зависимостям (47), (48) последовательно вычисляют 7, Л, к, а также ряды @'''(0) и [3(Ь). Эти ряды быстро сходятся, и их можно обрывать на сотом члене.
— По зависимостям (42)-(46) вычисляют величины ф, Т*, К1 и Д.
— По зависимости (41) определяют значения у, изменяя в цикле высоту сооружения Ь от 1 до 300 м с шагом 1 м. Круговой частоте ш присваивают числовые значения, характерные для сильных землетрясений. В результате получится график зависимости величины у от высоты сооружения, как это представлено на рис. 2. Величина у равна кратности увеличения амплитуд перемещений, скорости и ускорений гребня плотины по сравнению с теми же амплитудами для падающей на систему сейсмической волны.
График, представленный на рис. 2, получен при следующих числовых значениях параметров системы: частота сейсмической волны V =1 Гц (ш = 2ж), мощность верхнего слоя Н = 30 м, плотность материала сооружения р = 2500 кг/м3, плотности верхнего и нижнего слоев грунта р1 = 2000 кг/м3 и р2 = 2200 кг/м3, скорости распространения поперечной волны в слоях а1 = 600 м/с и а2 = 1000 м/с, модуль упругости тела плотины Е = 2 • 1010 н/м2, толщина плотины К =10 м. Фундаментный блок в данном примере считаем частью плотины.
7 16 14 12 10
8
6
4 2
0 50 100 150 200 250 300
Высота, м
Рис. 2. График зависимости относительных амплитуд колебания гребня плотины от ее высоты
1.0 Гц
Без учета воды
_ 1 - С учетом воды
Согласно рис. 2 при изменении высоты плотины в промежутке 65 < L < 75 м сооружение попадает в околорезонансное состояние и амплитуда его колебания увеличивается в 10-16 раз по сравнению с амплитудами падающей сейсмической волны, тогда как при высоте 110 < L < 170 м сооружение выходит из околорезонансного состояния и ее амплитуды увеличиваются только в 3.5-4 раза.
Как ни парадоксально, но при вышеприведенных значениях входных параметров системы увеличение высоты сооружения в пределах указанного промежутка минимизирует сейсмическое воздействие местности на сооружение.
Полученные в статье формулы могут быть приложены для решения весьма актуальной проблемы сейсмологии — автоматизированного подбора места строительства, а также габаритных размеров сооружения, при которых минимизируется динамическое воздействие падающей на систему сейсмической волны. Эта проблема решается путем дальнейшего исследования на условно-безусловный минимум полученной функции нескольких переменных
у = y(L,h,H,Gi,G2,ai,a2), где L и h — габаритные размеры сооружения. Остальные переменные характеризуют геологическое строение подфундаментного массива грунта. Круговой частоте ш как внешнему параметру в цикле придаются числовые значения, характерные для сильных землетрясений.
Таким образом, полученные в статье расчетные зависимости и результаты выполненных вычислительных экспериментов позволяют на стадии проектирования гидротехнических сооружений обоснованно выбирать место строительства и их габаритные размеры, при которых минимизируются колебательные воздействия в случае возникновения сейсмических событий.
Список литературы / References
[1] Лейбензон Л.С. О натуральных периодах колебания плотины, подпирающей реку: Труды. Т. 1. М.: Изд-во АН СССР, 1951. С. 157-161.
Lybenzon, L.S. On the natural periods of fluctuation of the dam which is propping up the river. Proceedings. V. 1. Moscow: Izdatelstvo AN SSSR, 1951. P. 157-161. (in Russ.)
[2] Кульмач П.П. Гидродинамика гидротехнических сооружений. М.: Изд-во АН СССР, 1963. 190 с.
Kulmach, P.P. Hydrodynamics of hydraulic engineering constructions. Moscow: Izdatelstvo AN SSSR, 1963. 190 p. (in Russ.)
[3] Westergaardt, Y.M. Water pressures on dams during earthquakes // Proceedings of the American Society of Civil Engineers. 1932. Vol. 58, No. 2. P. 274-279.
[4] Музаев И.Д., Музаев Н.И., Дзебоев Б.А. Математическое моделирование сейсмостойкости плотины с учетом влияния водной среды // Вычисл. технологии. 2012. Т. 17, № 1. С. 90-99.
Muzaev, I.D., Muzaev, N.I., Dzeboyev, B.A. Mathematical modeling of seismic stability of a dam, taking into account the influence of the aquatic environment // Computational Technologies. 2012. T. 17, No. 1. P. 90-99. (in Russ.)
[5] Музаев И.Д., Харебов К.С., Музаев Н.И. Постановка и решение многослойной краевой задачи сейсмических колебаний высотного сооружения // Сейсмостойкое строительство. Безопасность сооружений. 2014. № 1. С. 22-25.
Muzaev, I.D., Harebov, K.S., Muzaev, N.I. The multilayer boundary-value problem of the high-rise building seismic vibrations setting and the solution // Earthquake Engineering. Safety of Structures. 2014. No. 1. P. 22-25. (in Russ.)
[6] Кошляков Н.С., Глинер Э.Б., Смирнов М.М. Уравнения в частных производных математической физики. М.: Высшая школа, 1970. 710 с.
Koshlyakov, N.S., Gliner, E.B., Smirnov, M.M. The partial derivatives equations in mathematical physics. Moscow: Vysshaya Shkola. 1970. 710 p.(in Russ.)
[7] Шульман С.Г. Расчеты сейсмостойкости гидросооружений с учетом влияния водной среды. М.: Энергия, 1976. 335 с.
Shulman, S.G. Calculations of seismic stability for hydroconstructions taking into account influence of the water environment. Moscow: Energiya, 1976. 335 p. (in Russ.)
[8] Хачиян Э.Е. Задача усиления или ослабления эффекта сейсмического воздействия на поверхности земли // Вест. ЦНИИСК им. В.А. Кучеренко. Исследования по теории сооружений. 2009. № 1(XXVI). С. 67-80.
Hachiyan, E.E. The problem of strengthening or weakening of the effect of seismic impact on the ground // Vestnik TsNIISK im. V.A. Kucherenko. Issledovaniya po teorii sooruzheniy. 2009. No. 1(XXVI). P. 67-80. (in Russ.)
Поступила в 'редакцию 27 мая 2015 г., с доработки — 8 июля 2015 г.
Mathematical modeling of seismic fluctuations in a system consisting of reservoir, dam, base block and subbase layers of the earth
MUZAEV, ILLARION D.1'2'*, MUZAEV, NUGZAR I.1
1Geophysical Institute of the Vladikavkaz Scientific Center RAS, Vladikavkaz, 362002, Russia
2Financial University of the Government of Russian Federation. Vladikavkaz Branch, Vladikavkaz, 362002, Russia
* Corresponding author: Muzaev, Illarion D., e-mail: muzaevid@mail.ru
In this article we develop a mathematical model for seismic fluctuations of a system consisting of schematized reservoirs, dams, base block and subbase layers of the soil massif. The model is a contact boundary value problem for a differential equation of the potential movement of water in the reservoir, the differential equation fluctuations cantilever beam (beam model of the dam), as well as differential equations of transverse (shear) oscillations of two layers subbase solid ground. These equations are interconnected through the boundary conditions at the contact surfaces. The boundary value problem is solved analytically. Initially we solve the hydrodynamic part of the initial-boundary value problem. It consists of a differential equation of elliptic type with variable coefficients. The variability owes to factors due to the non prismatic configuration of the dam part of the reservoir, which is approximated by an exponential function of two spatial variables. With this approximation, an exact analytic solution of this differential equation with variable coefficients was found.
© ICT SB RAS, 2015
In the process of solving the boundary problem we obtain an integral-differential equation of Fredholm type, which was solved using approximations.
Numerical experiments have shown that the inclusion of the hydrodynamic water pressure at the upstream side of the dam significantly increases the amplitude of the displacement velocity and acceleration for oscillations of the dam. This effect must be considered when designing the construction.
Results of the computational experiments allowed us to research and design the hydraulic structures and to properly choose the place of construction and dimensions of the structures that minimize the dynamic effects of the incident on the system of seismic waves.
Keywords: seismic fluctuations of system, subbase layer of the soil massif, contact regional task, shift waves, neutral axis, potential movement of water.
Received 27 May 2015 Received in revised form 8 July 2015
