CC BY
58
Вычислительные технологии
Том 17, № 1, 2012
Математическое моделирование сейсмостойкости плотины с учетом влияния водной среды
И. Д. Мулл кв. Н. И. Музаев, Б. А. Дзебоев Центр геофизических исследований ВНЦ РАН и РСО-А, Владикавказ, Россия
е-таП: dba01@yandex.ru
Поставлена и решена контактная краевая задача совместных сейсмических колебаний плотины и воды в водохранилище. Получены расчетные формулы для вычисления частоты основных форм собственных колебаний указанной системы.
Ключевые слова: приплотинный рукав, формпараметр, базовая функция, метод Ритца, интегро-дифференциальное уравнение.
В горных условиях высоконапорные плотины строятся в наиболее узких ущельях рек. Как правило, приплотиппая область водохранилища имеет непризматическую конфигурацию, где глубина воды сравнительно большая, в некоторых водоемах достигающая 300 м. На Сарезском озере (в горах Памира) глубина воды у завальной плотины превышает 500 м. При землетрясении плотина и вода в рукаве совершают колебательные движения, и частота собственных колебаний системы на некоторых объектах может достичь околорезонансного значения, что в свою очередь может привести к аварии гидротехнического сооружения с трудно предсказуемыми последствиями.
Учитывая приведенные аргументы, исследования в данном направлении представляются весьма актуальными и имеют существенное значение при проектировании, строительстве и эксплуатации гидротехнических сооружений.
Предположим, что в прямоугольной системе координат хОуг часть пространства,
Б(х, г) Б(х,г)
ограниченная условиями 0 < х < Ь, — -
< У <
0 < г < И, представляет
2 * 2
собой узкий глубокий непризматический приплотинный рукав горного водохранилища, где Ь — дайна, Б(х, г) — шири на, И — глубина воды в рукаве. В створе х = 0 находится плотина, а в створе х = Ь рукав сопрягается с основным объемом водохранилища.
На рис. 1 представлена схема рассматриваемого гидрообъекта. Будем считать, что движение воды в рукаве является безвихревым. Тогда потенциал средней по ширине рукава скорости должен удовлетворять дифференциальному уравнению [1-3]
дх2
с краевыми условиями
д2^р д2^р
1 дВ(х, г) д^р 1 дБ д^р ^
дг2 Б(х, г) дх дх Б дг дг
д^ дх
ду
х=0
х=0
ду_ 0,
ду
Ж
х=Ь
ду дг
(1)
(2) (3)
г=И
0
-1 -**
н
Расчетная схема задачи
В(х, г)
В(Х, ж)
Рис. 1. Схема приплотинного рукава водохранилища: плановая (а) и глубинная (б) конфигурации
Поперечные колебания плотины моделируются следующей краевой задачей:
12 Эг2
= 0,
дЬ2 д3 V
дЬ
г=0
г=Н
= 0,
д*3 г=0
дV(г,Ь)'
= 0,
дг
= 0,
(4)
(5)
(6)
г=Н
где ф(х, г, Ь) потенциал скорости движения воды; V(г, Ь) поперечные перемещения в теле плотины, определяемые вертикальной координатой г и временем Ь; В(х, г) ширина приплотинного рукава, зависящая от продольной х и вертикальной г координат; Н ширина плотины по направлению оси Ох; Е модуль Юнга (модуль упругости) материала плотины; р плотность в оды; р2 плотность материала плотины.
Легко заметить, что непризматическая конфигурация приплотинного рукава моделируется через функцию В(х, г). Функция В(0,г) характеризует изменение длины
д^, дV
плотины в вертикальном направлении. 1 раничное условие —— _п = —— выражает ра-
дх 1х=0 дЬ
х=0 дР I
ничное условие —— _т = 0 означает, что в створе сопряжения приплотинного рукава дЬ 1х=ь
с основным объемом водохранилища давление по глубине распределено по гидроста-
1
тистическому закону, граничное условие -—— _п = 0 — что на свободной поверхности
дЬ 1г=0
д^ 1
воды г = 0 пренебрегаются волнообразования. Граничное условие \г_н = 0 выража-
дг
рукава.
г=0
на нее изгибающий момент и поперечная сила равны нулю, граничные условия (6) что основание плотины г = Н жестко защемлено.
Выражения (1) (6) в совокупности представляют краевую задачу совместных свободных колебаний плотины и воды в рукаве. Наиболее полный анализ исследований в данном направлении представлен в монографиях [4-7].
В дифференциальных уравнениях (1) и (4) коэффициенты являются переменными, зависящими от пространственных координат х и г, В связи с этим решение поставленной краевой задачи (1)-(6) аналитическими методами связано с большими математическими трудностями, легко преодолеваемыми в частном случае, когда ширина водоема и створа плотины аппроксимируются экспоненциальной функцией вида
B{x,z) = B0es1xe
S1 = const, s2 = const, B0 = const.
(7)
При такой аппроксимации непризматической конфигурации водоема коэффициенты уравнений (1) и (4) становятся постоянными величинами, в силу чего решение краевой задачи существенно облегчается.
Приступая к решению контактной краевой задачи (1)-(6), искомые функции <р и V представим в виде
ip(x,z,t) = ф(х, z)e(s2/2)z sinut, V(z,t) = W(z) cos ut,
(8)
где ш — круговая частота собственных колебаний системы. Подставив выражения (7), (8) в (1)-(6) относительно введенных искомых функций ф и Ш, получим следующую контактную краевую задачу:
d4W п d3W - 2s2-
д2ф 82ф dx2 dz2
d2W
дф дх
+ Si— -J-ф = О,
4
dz4
dz3 дф
дх
ф
+ s2
- 12
u2W
12
dz2 a2h2 h3E -uW(z)e-S2/2)z, ф
p1ue
(s2/2)z
ф(0, z),
x=0 z=0
d2W
x=L
0,
dz2
W
, 52 / \ I n
d3W
z=0
z=H
0,
dz3
dW'
z=0
dz
(9)
(10) (11)
(12)
(13)
(14)
z=H
Решение дифференциального уравнения (9) будем искать в виде тригонометрического ряда по синусам
ф(х, z) = Фп(х) sin (j^nj^j
n= 1
где ^n — положительные корни трансцендентного уравнения
^П I S2 • n 1 о о
— eos ¡jLn Н--sm ¡jLn = ü, n = 1, 2, 3,..
H 2
или
tg Vn = -
S2H
(15)
(16)
(17)
0
2
Легко проверить, что если цп — корень характеристического уравнения (16), то выражение (15) автоматически удовлетворяет граничным условиям (12) [8].
Разложим функцию Ш(г)е-(з2/2)* в тригонометрический ряд по ортогональной системе функции в промежутке (0,Ь) [8]:
г г г
вт//! —, 8т/Л2—, ЙШ/1П—, ... (18)
те
= ^„зт^, (19)
H
n=1
H
W(z)e~{s2/2)z sin nn^dz
an = ^---. (20)
J sin2 ¡in^dz 0
Подставим разложения (15) и (19) в дифференциальное уравнение (9) и в граничные условия (11) и, приравняв коэффициенты при синусах в левых и правых частях, получим
d\¡)n , dipn
71 TF? + ~Г V'n = 0, (21)
dx2 dx \ H2 4 dф
dx
= 0. (22)
x=L
= -ШО.п, ф,
х=0
Решение дифференциального уравнения (21) с граничными условиями (22) имеет следующий вид:
Фп(х) =--^- е-^х8ЪХп(х-Ь)8т11п^-, (23)
\п сЬ АПЬ + — эЬ АПЬ н
A + ^ + п= 1,2,3... (24)
Для функции ф(х,г) получим ряд
те
(x,z) = -\- e~{si/2)xsh\n(x- L) sin¡in—.
n= i An ch AnL + — sh AnL
Подставим значения an и гф(0,г) из (20) и (25) в дифференциальное уравнение (10)
и ПОЛУЧИМ
d4W d3W 2 d2W u2W
dz4 dz3 dz2 a2h2
z H
io f) sh XnL sin r
-_ H , x W(z)e~^2>Smi,n^dz, (26)
haP2 tÍ[Kch\nL + ^sh\nL)pn J H
^=L2ñL-sin2^_ (27)
2 2^n
Выражение (26) представляет собой интегро-дифференциальное уравнение типа уравнения Фредгольма относительно искомой функции W(z).
Для определения собственной частоты основной формы колебания системы далее будет использован приближенный метод Ритца [9].
Предположим, что на плотину действует равномерно распределенная статическая нагрузка
q(z) = q0 = const. (28)
Поперечный изгиб плотины как непризматической балки описывается дифференциальным уравнением
^ # (0 о
12В5?Ге 1ё)=в"е *
с граничными условиями
(29)
d2W
dz2 W
z=0
0,
0,
d3W
dz3 dW'
z=0
z=H dz
Решение краевой задачи (29)-(31) имеет вид
0.
z=H
(30)
(31)
(32)
где
^i(z)
to to| ^ zeS2Z + s3 s2 3 -enx + s2 HeS2H 2 s2 H
- s2 Не«я _ si H2 M1PS2H 2 e . s2 (33)
Теперь рассмотрим случаи, когда на напорную грань плотины действует гидростатическое давление
(34)
Ф) = Qojj,
W(z)
(35)
^2(z)
6s2
— + — - —eS2Z
S3 04 ^5
2 s2 s2
4е«я _ 4Яе,аЯ
s2 S2
,S2 H 2
4 „3е „3 о „2
S о S о ¿/s о
3
■e
S2H
Я2
c>3
tf3 3s2'
(36)
'2 °52
Согласно методу Ритца [9] в качестве базовых можно использовать функции и Составим линейную комбинацию от базовых функций
W (z) = <£i(z)ai + P2(z)a,2,
(37)
0
z
где а1 и а2 — произвольные постоянные. Также согласно методу Ритца [9] выражение (37) следует подставить в интегро-дифференциальное уравнение (26), В результате получим
Г 12, .2
^4)(г) - + в&((г) - ^М*)
+
- 282$'+ -
,(3)
Л2 а2 12,2
а1+
а2 =
2 ЛП тз
Е 7Г (Х+ 81п 17*'
Р2 а2 Л,3 ' вп Н
п=1
Н Н
к1>п = J (р!(г)е {82/2)хътк2,п = ^ (39)
о о
г
Я
Лга сЬ Ага1/ + — вЬ ХПЬ
2
Легко проверить достоверность следующих равенств:
(38)
(40)
<14)(г) - 282<Г;(г) + в2<1/(г) = 1,
(3)
(41)
<24)(г) - 2^ (*) + (*) = г, 2 ^ , ,2 ^
(3)
12,
г —
й2а2
12а;2 1ъ2а2
п=1
г Я
<2(г)
Р1 а;2 12 ^л Д
р2 а2 И3 Д
п=1
г Н
а1 +
а2 = 0.
(42)
(43)
Согласно методу Ритца обе стороны выражения (43) последовательно умножим на базовые функции ^(г) и <2 (г) и проинтегрируем по переменной г в пределах от 0 до Н,
а1 а2
Н
12а;2 /¿2а2
Н
+
о
Н
/Р1 ,2 12 Я
<^2(г) сЬ + —— — V
Р2 а2 Л,3 «=1 в«
Н
о о п=1
а1+
а2 = 0,
Н Н
2
о о п=1
а1+
+
Н
12а;2 /¿2а2
Н
2/ 4 7 , Р1 ,2 12 ^ Д
2 и Е
Р2 а2 п=1 вп
а2 = 0,
(44)
где
H H
di,n = j e{s2/2)z^i(z) sin finjjdz, d%n = J e{s2/2)z^2(z) dz. (45)
00 Введем обозначения
H H
0 0 n=1
H H
/12 f pi 12 R
ZLpi(z)dz, m4 = — ipi(z)ip2(z) dz---йУ^т^Л».
h2 J P2 h3 Pn
0 0 n=1
H H
pi= <f2(z) dz, P2 = I Lpi(z)ip2(z) dz - — — V -^khnd2yn,
J h J p2 h n=1 Pn
0 0 n=1
H H
Рз= ZLp2(z)dz, P¿± = — / <fl(z) dz - — — У2 -k2ynd2yn. (46)
J h J p2 h n=1 Pn
0 0 n=1
В принятых обозначениях система уравнений (44) запишется в виде
w2 N f и2 .
тi---т2 I ai + I m3---т4 ) а2 = О,
22 и2 \ Í и2
(47)
Р\ ~ -^Р*) аг + - ^РА) а2 = 0.
Приравняв к нулю определитель системы (47), получим частотное уравнение (■т2Р4: - т±Р2) ^^ - {ГП\РА + т2Рз - т^Р2 - гпаР\) х
Х (^)2 + т1рз-т3Р1 = 0. (48)
, 4
. а -
^ 2
а
Решение уравнения (48) относительно круговой частоты собственных колебаний системы имеет вид
_ . miP4 + т2Р3 - тзР2 - m4Pi_^ 2(т2Р4: - т4Р2)
' (miP4 + т2Рз - тзР2 - m4Pi)2 miP3 - тзР1
1/2
±* Г 1 * ' 3— д " ' ^ ^ - ———^^ . (49)
4(ш2Р± - ш±Р2)2 т2Р4 - т4Р2 '
На рис, 2 представлены графики зависимости частоты собственных колебаний плотины от ее высоты, на рис, 3 — графики зависимости частот собственных колебаний системы от формпараметров 51 и з2 при высоте плотины Н = 150 м. Остальные входные параметры имеют те же значения, что приняты к рис, 2,
Вычисления проводились с помощью специально разработанных алгоритмов в среде МАТЬАВ в следующей последовательности.
Рис. 2. Зависимость частоты собственных колебаний плотины V = — от се высоты Н: а —
2п
первая форма колебаний, б — вторая форма колебаний; входные параметры: Ь = 5000 м, р2 = 3000 кг/м3, К = 10 м, 51 = 1п(2)/Ь, в2 = 1п(2)/Н, Е = 2.5 • 1010 Н/м2; сплошные линии — с учетом, пунктирные — без учета гидродинамического давления воды на плотине
„ a
v, Гц v, Гц
Рис. 3. Зависимость частоты собственных колебаний плотины v от формпараметров si, в2'-a — si, первая (слева) и вторая (справа) формы; б— S2, первая (слева) и вторая (справа); H = 150 м, L = 5000 м, р2 = 3000 кг/м3, h = 10 м, E = 2, 5 • 1010 Н/м2; сплошные линии — с учетом, пунктирные — без учета гидродинамического давления воды на плотине
I этап. В характеристическое уравнение (16) подставляются значения величин s2, H и вычисляется последовательность чисел (n = 200), представляющих положительные корни характеристического уравнения (16) поставленной краевой задачи. Трансцендентное уравнение (16) решается с помощью встроенной в MATLAB функции "'\fzero"'\
Вычисляютея последовательности чисел
Ai, А2, . . . , An, . . . , в1, ^2, . . . , fin, . . . , R1, R2, . . . , Rn ... (n = 200)
по формулам (24), (27) и (40) соответственно.
II этап. Вычисляются последовательности чисел
ki,n, k2,n, di,n, d2,n, mi, m2, тз, m4, Pi, P2, P3, P4 (n = 200)
по формулам (39), (45), (46) соответственно. Определенные интегралы находятся по формуле трапеции.
Вычислительные эксперименты показали, что в бесконечных суммах, содержащихся
n=
200 n
III этап. Полученные значения вышеперечисленных величин подставляются в формулу (49), В результате получаются два значения круговой частоты, меньшее из которых соответствует первой (основной) форме колебания рассматриваемой системы.
Таким образом, как показали результаты проведенных вычислительных экспериментов, собственные частоты колебаний системы существенно зависят от формпара-метров si и s2, характеризующих геометрическую конфигурацию чаши водохранилища и створа плотины. Для повышения сейсмостойкости рассматриваемого гидротехнического сооружения в проектных разработках габаритные размеры плотины и водохранилища необходимо подбирать так, чтобы частоты собственных колебаний системы резко отличались от частоты сейсмических колебаний местности.
Список литературы
[1] Музаев И.Д., Кусраев А.Г., Созанов В.Г. Математическое моделирование некоторых задач волновой гидродинамики применительно к горным водоемам // Тр. Ин-та вычисл. технологий СО РАН. Новосибирск. 1995. Т. 4, № 11. С. 164-168.
[2] Музаев И.Д., Созанов В.Г. К теории поверхностных гравитационных волн Коши — Пуассона в узких глубоких непризматических водоемах // Изв. вузов. Естеств. науки. 1995. № 3. С. 40-43.
[3] Музаев Н.И., Музаев И.Д. Постановка и решение начально-краевой задачи поверхностных гравитационных волн в водохранилище узкоканьонного типа // Там же. 2009. № 2. С. 22-24.
[4] Гольденвлат И.И., Николаенко H.A. Расчет конструкций на действие сейсмических и импульсивных нагрузок. M.: Госстройиздат, 1961. 317 с.
[5] Кульм a4 Т.П. Гидродинамика гидротехнических сооружений. M.: Изд-во АН СССР, 1963. 190 с.
[6] Уразваев М.Т. Сейсмостойкость упругих и гидроупругих систем. Ташкент: Фан, 1966. 252 с.
[7] Шульман С.П. Расчеты сейсмостойкости гидросооружений с учетом влияния водной среды. М.: Энергия, 1976. 335 с.
[8] Кошляков Н.С., Глинер Э.Б., Смирнов М.М. Уравнения в частных производных математической физики. М.: Высшая школа, 1970. 710 с.
[9] Тимошенко С.П. Колебания в инженерном деле. М.: Физматгиз, 1959. 440 с.
Поступила в редакцию 18 июня 2011 г., с доработки — 8 августа 2011 г.
