Метод синтеза сейсмических колебаний, соответствующих заданному семейству спектров ответа Текст научной статьи по специальности «Математика»

УДК 517.9:519.8 Вестник СПбГУ. Сер. 10, 2013, вып. 2

С. А. Дурновцева

МЕТОД СИНТЕЗА СЕЙСМИЧЕСКИХ КОЛЕБАНИЙ, СООТВЕТСТВУЮЩИХ ЗАДАННОМУ СЕМЕЙСТВУ СПЕКТРОВ ОТВЕТА

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

Проектирование АС выполняют с учетом двух уровней сейсмичности: проектного землетрясения (ПЗ) и максимального расчетного землетрясения (МРЗ) [1]. ПЗ - максимальная интенсивность сейсмических воздействий, которая с определенной вероятностью может возникнуть в течение срока службы АС, но не причинить ей каких-либо существенных повреждений и не привести к остановке функционирования. МРЗ - максимальная расчетная интенсивность сейсмических воздействий, вероятность возникновения которой хотя и очень мала, но макросейсмические повреждения могут оказаться значительными, вызвать существенные повреждения и/или частичную потерю устойчивости сооружения, но не должны привести к его обрушению или создать серьезную угрозу для окружающей среды.

Наглядным примером, демонстрирующим актуальность сейсмического проектирования, может служить авария на АС «Фукусима-1», произошедшая 11 марта 2011 г. Первый блок станции был спроектирован на пиковое ускорение ПЗ 0,18 g (1.74 м/с2) [2] и выдержал землетрясение, оцененное в 6 баллов [3] по японской шкале JMA, соответствующее 9 баллам по шкале MSK-64, т. е. приблизительный уровень ускорения доходил до 4 м/с2. Известные повреждения были вызваны потерей электроснабжения (в том числе и от резервных дизельных электростанций) вследствие цунами. Однако в момент землетрясения энергоблок был остановлен действием системы аварийной защиты, которая сработала в штатном режиме. Это означает, что станция выдержала землетрясение, как минимум на балл (по шкале MSK-64) превысившее максимальное, при котором ее проект мог гарантировать безопасный останов.

2. Необходимость формирования акселерограмм. Современные подходы к сейсмическому проектированию ядерных энергетических объектов подразумевают, что помимо балльности, характеризующей уровень сейсмической опасности зоны строительства, в качестве исходного сейсмического воздействия задаются проектные спектры ответа [4].

Спектром ответа (СО) называют совокупность значений (график, таблица) абсолютных максимальных ответных ускорений одномассовой колебательной системы (сейсмо-

Дурновцева Светлана Александровна — аспирант, Санкт-Петербургский государственный университет; e-mail: [email protected].

© С. А. Дурновцева, 2013

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

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

Кроме того, часто возникает необходимость провести расчет конструкции с уровнем демпфирования когда для нее задан СО с уровнем демпфирования £2 (£1 = (2). В таких случаях СО пересчитывают для уровня £1 и синтезируют акселерограмму. Но определенный по такой акселерограмме СО с исходным уровнем демпфирования £2 может значительно отличаться от заданного СО, т. е. не обеспечивается соответствие акселерограммы и ее СО.

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

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

тх + Ьх + кх = -та(Ь), (1)

в котором т - масса; Ь - коэффициент затухания; к - жесткость; х, х и х - ускорение, скорость и перемещение массы; а(1) - переносное ускорение.

После деления (1) на т получаем

х + 2£х + и2х = -а(Ь), (2)

где ш = у/к/т - круговая частота осциллятора без демпфирования; £ = Ь/(2тси) = - безразмерный коэффициент демпфирования.

С математической точки зрения СО представляет собой функцию

Яах(и) = тах \хш (¿)|.

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

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

На рис. 1 показано семейство исходных обобщенных СО (жирные линии) для уровней демпфирования 0.5, 2, 5, 10 и 20%. На рис. 2 приведена акселерограмма, синтезированная для исходного СО с уровнем демпфирования 5%. Тонкими линиями показаны СО, рассчитанные для этой синтезированной акселерограммы. Как видно на рис. 1, спектры синтезированной акселерограммы значительно отличаются от исходных данных (стрелками обозначены отклонения исходных спектров от рассчитанных по акселерограмме, указана относительная погрешность). Таким образом, методы синтеза акселерограмм, работающие только с одним СО, не обеспечивают достаточной точности задания сейсмической нагрузки.

10 1 10" 101 102 ю

Рис. 1. Семейство спектров, соответствующих акселерограмме, синтезированной по СО с уровнем демпфирования 5%

"0 2 4 6 8 10 12 14 16 18 20 *

Рис. 2. Акселерограмма для СО с уровнем демпфирования 5%, входящего в семейство на рис. 1

5. Решение поставленной задачи. Поиск функции во временной области крайне сложен, так как представление сигнала во временной области не дает информации о наличии в нем тех или иных частотных составляющих. Для нахождения сейсмического возмущения, соответствующего семейству СО, используется преобразование Фурье [7], которое позволяет получить спектральное представление сигнала (зависимость амплитуды от частоты).

Предлагаемый подход к синтезу акселерограмм носит итеративный характер и состоит из следующих этапов:

• выбор начального приближения;

• формирование вектора параметров оптимизации;

• вычисление целевой функции;

• минимизация целевой функции;

• нахождение решения.

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

Так, согласно [8, 9], шаг дискретизации (сИ) акселерограмм во времени должен составлять не более 0.01 с. Параметры огибающей колебаний (сглаженной функции закономерного изменения пиковых амплитуд во времени) определены в [10], ее длительности временных участков в зависимости от магнитуды землетрясения указаны в таблице.

Параметры огибающей функции

Магнитуда Время нарастания Продолжительность Время затухания

амплитуды, 1г, с максимальных амплитуд, с колебаний, t¿, с

7.0-7.5 2 13 9

6.5-7.0 1,5 10 7

6.0-6.5 1 7 5

5.5-6.0 1 6 4

5.0-5.5 1 5 4

В работе использована огибающая функция ¡епу вида

{вт(^), если 0 < Ь < ¿г, 1, если < Ь < + ,

вт(|21), если и +гт < t < и +гт -И<г-

При известных шаге дискретизации и огибающей функции количество точек будущей акселерограммы определяется как пр = (Ьг + Ьт + Ь^/А. Поскольку алгоритмы быстрого преобразования Фурье (БПФ) [11] работают с максимальной скоростью при количестве точек, равном степени двойки, то положим число точек акселерограммы N = 2к так, чтобы 2к-1 < пр ^ 2к. При этом формально длительность акселерограммы (а значит, и огибающей функции) увеличится на N — пр)А секунд. Однако будем считать, что (Ь) = 0 при + Ьт + <Ь ^ + Ьт + + N — пр)А. Такое допущение позволяет сохранить форму огибающей функции, удовлетворяющую [10], и получить оптимальное количество параметров для БПФ.

Очевидно, что сейсмическое воздействие по своей природе является суперпозицией гармоник. Пусть число составляющих гармоник равно Ц- + 1, тогда можно сформировать начальное приближение к решению, положив амплитуды и фазы гармоник случайными числами в диапазонах от 0 до 1 и от 0 до 2п соответственно. Составив суперпозицию с положенными амплитудами и фазами и умножив ее на огибающую функцию, найдем максимальное по модулю значение нулевого приближения к акселерограмме т&Х1 ао(Ь). Окончательно скорректируем начальное приближение к решению, умножив амплитуды на отношение ZPA/ т&Х1 ао(Ь), где ZPA - так называемое ускорение нулевого периода (ускорение СО в асимптотической или твердотельной области спектра, которая обычно лежит в диапазоне частот более 33 Гц).

5.2. Формирование вектора параметров оптимизации. Применив преобразование Фурье к начальному приближению, получим

п

а(г) = £ Хк в*"»*,

к=1

где Хк - к-тая комплексная амплитуда; ки - к-тая круговая частота гармонического колебания.

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

Ак = ±у/Ъе(Хк)+1т(Хк),

срк = arctan к = I, N/2.

Тогда вектор x параметров оптимизации будет состоять из следующих компонентов: x = [Ai ,A2, ...,An/2,¥>1, P2, ■■■, ^N/2]N xi-

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

/ \ V2

f (a(t)) = [J2J2(RSlarget^) - RSLrrentH)) ■

^ j и '

Здесь RSjarget - j-й исходный (целевой) спектр ответа, RS3current - текущее приближение к j-му целевому спектру ответа. Нормативные документы [6, 8-10] говорят о том, что СО должны быть рассчитаны с достаточно малым шагом по частоте (для 5%-го спектра количество частот в диапазоне от 0-1 до 50 Гц должно быть не менее 100) и логарифмически расположены- На практике исходные спектры зачастую заданы небольшим количеством точек, поэтому перед вычислением целевой функции СО интерполируются по новому частотному диапазону, найденному по правилу ^i+i = , где Zj -коэффициент демпфирования j-го СО.

Для вычисления СО уравнение (2) численно интегрируется [12] для всего диапазона частот и всех уровней демпфирования. Функция правой части a(t) восстанавливается обратным преобразованием Фурье по вектору параметров на текущей итерации. Дополнительно текущее приближение к акселерограмме a(t) «огибается» - умножается на огибающую функцию и уже в таком виде передается в качестве аргумента в целевую функцию.

Очевидно, что наличие высоких частот в спектре Фурье приводит к недопустимому росту СО в высокочастотной области. Для исключения подобной ситуации перед восстановлением a(t) коэффициенты Xk, соответствующие высоким частотам, полагаются равными нулю. Рассмотрим подробнее процедуру «обрезки» высоких частот.

Согласно теореме Найквиста-Котельникова [7], аналоговый сигнал, имеющий ограниченный по ширине спектр, может быть восстановлен однозначно и без потерь по своим дискретным отсчетам, взятым с частотой, строго большей удвоенной верхней частоты спектра. Приведенная теорема дает возможность вычислить частоты, по которым исходный сигнал раскладывается в спектр Фурье. Пусть dt - шаг дискретизации, N - количество отсчетов, тогда частота Найквиста равна VNyq =

Учитывая особенности алгоритма БПФ, которые приводят к смещению отрицательной части спектра вправо, можем построить спектр Фурье для N +1 частоты от —VNyq до +VNyq с шагом 5v = -щ^. При этом Xq будет соответствовать частоте v = 0, диапазон номеров от 1 до N/2 занимают положительные частоты, оставшиеся значения соответствуют отрицательным частотам.

Задав частоту «обрезки» vcut, найдем в диапазоне отрицательных частот такой номер к, чтобы выполнялось Vk ^ vcut ^ Vk+i. Тогда для исключения нежелательного роста СО правее vcut достаточно положить Xi = 0 при i = N/2 — ncut...N/2 + ncut, где

ncut = N — к.

Применение описанного приема для иси = 30 Гц демонстрируют рис. 3-5.

а

■2-'-'-'-'-'-'-'-'-'-1 ,

0123456789 10 *

Рис. 3. Удаление высоких частот из акселерограммы

Рис. 4. Фурье-спектры акселерограмм на рис. 3

Сплошная линия — спектр исходной акселерограммы, пунктирная — спектр акселерограммы после удаления высоких частот.

Рис. 5. СО, рассчитанные по акселерограммам на рис. 3

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

5.4. Многомерная минимизация целевой функции. Для минимизации используется модифицированный алгоритм Хука-Дживса [13]. Поиск оптимального

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

С учетом «обрезки» высоких частот количество параметров оптимизации сокращается до N — 2пси1. Вычисляем значение целевой функции в начальной точке х = [А1 ,А2, ...,Ак-2пи ,Р1,Р2, ■■■, -2пи](м-2пс„4)х1. Задаем вектор шагов по компонентам вектора х:

На первом этапе каждая компонента вектора x поочередно изменяется прибавлением и вычитанием соответствующего шага. Если значение целевой функции после изменения вектора параметров в г-й компоненте улучшилось, то заменяем эту компоненту на измененную: Xi := xi ± stepi.

В классическом алгоритме в случае, если шаги по всем координатам вектора x не привели к улучшению целевой функции, то исследование повторяется вокруг той же точки x, но с вдвое уменьшенной длиной шага. Модифицируем алгоритм Хука-Дживса следующим образом: делим пополам шаг stepi в случае, если из исходной точки удалось сделать лишь небольшое количество успешных шагов (менее 20% от общего числа параметров) или ^ < 0.05. Шаг называем успешным, если значение целевой функции после такого шага меньше, чем в исходной точке. Отношение ^ j^1 показывает, насколько быстро убывает целевая функция, в случае падения этой величины ниже 0.05 полагаем длину шага неэффективной для дальнейшей оптимизации.

Во втором этапе алгоритма Хука-Дживса - «поиске по образцу» - используется

информация, полученная в процессе «исследующего поиска». Пусть x1 - базисная точка

до выполнения первого этапа алгоритма, x2 получена в результате «исследующего поис-

2 1

ка», тогда разумно двигаться дальше в направлении x2 — x1 , поскольку поиск в этом направлении уже привел к улучшению целевой функции.

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

5.5. Нахождение решения. Искомая акселерограмма вычисляется как точка минимума целевой функции

6. Полученные результаты. Пример расчета по разработанному методу приведен на рис. 6, 7.

На всем диапазоне частот, для которого приведены СО на рис. 6, были вычислены максимальные отклонения СО, рассчитанных по синтезированной акселерограмме (маркированные линии) от исходных СО (жирные линии). Приведем оценки этих отклонений:

если г = 1, N/2 — ncut,

если г = N/2 — ncut, N — 2п,

cut.

(3)

a(t) = argmin(/(a(t))).

Демпфирование, % Максимальная относительная погрешность, %

0.05 2

18.1 35.8 28.2 11.9 6.1

5

10 20

12 г

Рис. 6. Исходное семейство СО, приведенное на рис. 1 (жирная линия) и семейство СО, рассчитанное по акселерограмме, синтезированной разработанным методом (маркированная линия)

а

для семейства СО, приведенного на рис. 1

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

Разработанная модификация существующего оптимизационного алгоритма позволяет ускорить процесс синтеза акселерограмм.

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

Уменьшение относительных погрешностей по сравнению с существующими методами (см. рис. 1) демонстрирует эффективность предложенного алгоритма синтеза.

Литература

1. НП-031-01. Нормы проектирования сейсмостойких атомных станций. Введ. 01.01.2002. М.: Госатомнадзор России, 2001. 27 с.

2. Ellingwood B. An Investigation of the Miyagi-ken-oki, Japan, earthquake of June 12, 1978: NBS special publication. Tokyo: NBS, 1980. 225 p.

3. Earthquake hazards division, disaster prevention research institute, Kyoto university // URL: http://www.eqh.dpri.kyoto-u.ac.jp/ masumi/ecastweb/110311/index.htm.

4. Бирбраер А. Н. Расчет конструкций на сейсмостойкость. СПб.: Наука, 1998. 354 с.

5. Клаф Р., Пензиен Дж. Динамика сооружений / пер. с англ. Л. Ш. Килимника, А. В. Швецовой; под ред. Г. А. Казиной. М.: Стройиздат, 1979. 320 с. (Clough Ray W., Penzien J. Dynamics of structures.)

6. РБ-06-98. Определение исходных сейсмических колебаний грунта для проектных основ. Введ. 01.07.1999. М.: Госатомнадзор России, 2000. 76 с.

7. Марпл-младший С. Л. Цифровой спектральный анализ и его приложения / пер. с англ. О. И. Хабарова, Г. А. Сидоровой; под ред. И. С. Рыжака. М.: Мир, 1990. 265 с. (Marple S. Lawrence, jr. Digital spectral analysis.)

8. МР 1.5.2.05.999.0025-2011. Расчет и проектирование сейсмостойких атомных станций. Введ. 17.10.2011. СПб.: ОАО «Концерн Росэнергоатом», 2011. 140 с.

9. ASCE/SEI 43-05. Seismic design criteria for structures, systems, and components in nuclear facilities. Reston: American society of civil engineers, 2005. 81 p.

10. ASCE-4-98. Seismic analysis of safety-related nuclear structures and commentary. Reston: American society of civil engineers, 1998. 118 p.

11. Cooley J. W., Tukey J. W. An algorithm for the machine computation of the complex Fourier series // Mathematics of Computation. 1965. Vol. 19. P. 297-301.

12. Gupta Ajaya K. Response spectrum method in seismic analysis and design of structures. Boston: Blackwell Scientific Publications, 1990. 170 p.

13. Банди Б. Методы оптимизации. Вводный курс / пер. с англ. В. А. Волынского. М.: Мир, 1988. 127 с. (Bunday Brian D. Basic optimisation methods.)

Статья рекомендована к печати проф. Е. И. Веремеем. Статья поступила в редакцию 20 декабря 2012 г.