CC BY
17
УДК 519.6:532.529.5
Гибридные методы вычислительной диагностики двухфазного потока в циркуляционном контуре
© В.Д. Сулимов, П.М. Шкапов МГТУ им. Н.Э. Баумана, Москва, 105005, Россия
Рассмотрены задачи вычислительной диагностики потока теплоносителя в замкнутом циркуляционном контуре. Разработаны математические модели акустических колебаний в двухфазном потоке. Использована косвенная диагностическая информация, которую содержат спектры колебаний потока, регистрируемые штатными системами. Сформулирована обратная задача на собственные значения, при решении которой реализован оптимизационный подход. Предполагается, что частные критерии представлены непрерывными, липшицевыми, не всюду дифференцируемыми, многоэкстремальными функциями. Поиск глобальных решений проведен с использованием новых гибридных алгоритмов, интегрирующих стохастический алгоритм сканирования пространства переменных и детерминированные методы локального поиска. Приведен численный пример модельного диагностирования фазового состава теплоносителя в циркуляционном контуре ядерной реакторной установки.
Ключевые слова: двухфазный поток, обратная задача, регуляризация, глобальная оптимизация, алгоритм Метрополиса, гибридный алгоритм.
Введение. Математическое моделирование многофазных потоков является актуальным направлением современных исследований. Следует отметить, что к настоящему времени, например для сжимаемых двухфазных потоков, отсутствуют общепринятые формулировки разрешающих уравнений как основы формирования полной математической теории начально-краевых задач, а также численных методов их решения [1]. Стандартный подход к выводу таких уравнений для двухфазного потока основан на процедуре осреднения, что ведет к системе уравнений в форме законов баланса массы, момента и энергии. В работе [2] представлена одномерная модель вязкого двухфазного потока, разработанная с использованием условия гидродинамического замыкания. Так, для модели вязкого газожидкостного потока сформулированы условия существования и единственности слабых решений. Проблема вибрации, возникающей при реализации математической модели двухфазного потока, исследована в работе [3]. Отмечено, что наличие вибрации снижает вычислительную эффективность модели и ограничивает возможности ее применения. В работе [4] дано описание математической модели двухфазного потока при наличии растворимых частиц (без перемешивания). Приведены результаты численного моделирования, соответствующие различным
значениям начальных концентраций частиц в двухфазном потоке. Проблеме существования и единственности глобальных слабых решений при моделировании двухфазных потоков посвящены работы [5, 6]. Значительный интерес представляют исследования, в которых сопоставляются результаты численного моделирования с использованием коммерческих комплексов программ, и экспериментальные данные. Результаты такого анализа для стратифицированного двухфазного потока представлены в [7].
Практическое значение в аспекте обеспечения безопасной длительной эксплуатации ядерных реакторов под давлением имеют исследования двухфазных газожидкостных потоков теплоносителя в циркуляционных контурах [8—11]. К числу актуальных направлений относится вычислительная диагностика фазового состава теплоносителя. Такая диагностика включает в себя методы и средства, предназначенные для определения характеристик исследуемых объектов по некоторой косвенной информации о них, измеряемых штатными средствами. Подход основан на формулировке и последующем решении соответствующей обратной задачи, обычно связанной с минимизацией некоторого функционала невязки. Принципиальной особенностью вычислительной диагностики является возможность использования значительных объемов информации об исследуемых объектах, для обработки и интерпретации которой применяется специализированное алгоритмическое и программное обеспечение, реализуемое на высокопроизводительных компьютерах. Выбор диагностической информации определяется, в частности, наличием штатных систем, регистрирующих полезные сигналы. Такую информацию содержат, например, спектры акустических колебаний в двухфазном потоке теплоносителя. Поэтому значительное внимание при моделировании двухфазных потоков уделяется разработке эффективных методов решения прямой задачи [12]. Естественными критериями качества математической модели диагностируемого объекта являются ее точность, вычислительная эффективность, способность корректно воспроизводить свойства объекта в требуемых пределах изменения переменных модели (переменных управления). Коррекцию моделей проводят с использованием результатов численного моделирования и соответствующих экспериментальных данных [13]. Примеры коррекции моделей двухфазных потоков приведены в [14, 15]. Следует отметить, что процедура коррекции модели связана, в свою очередь, с решением некоторой обратной задачи. При формулировке обратных задач коррекции моделей и диагностирования двухфазных потоков, в частности по спектральным данным, необходимо обеспечить корректность постановки задачи, а также учесть неполноту косвенной информации, наличие в спектрах систем кратных частот, зашумленность измеря-
емых данных и др. [16-22]. Как следствие, критериальные функции обратных задач в общем случае являются непрерывными, липшице-выми, многоэкстремальными и не всюду дифференцируемыми. Примеры использования методов глобальной оптимизации в задачах идентификации переходных процессов и диагностирования ядерных реакторов представлены в работах [23, 24].
При обеспечении безопасной эксплуатации реакторных установок значительное внимание уделяется исследованиям переходных процессов в циркуляционных контурах реакторов под давлением, в том числе контролю фазового состава теплоносителя. Появление второй фазы в потоке теплоносителя приводит, в частности, к изменению значений относительной скорости звука на участках локализации газожидкостной смеси. Это проявляется в соответствующих изменениях спектра колебаний потока, что может быть использовано в качестве косвенной информации для диагностирования фазового состава газожидкостной смеси. Критериальные функции обратной задачи определяются рассогласованием спектральных составляющих, полученных для математической модели потока, и соответствующих данных, регистрируемых штатными системами. При минимизации критериальных функций в общем случае требуются методы глобальной недифференцируемой оптимизации. Некоторые современные методы недифференцируемой оптимизации, основанные на построении сглаживающих аппроксимаций критериальных функций, представлены в работах [25—27]. Детерминированные методы решения задач глобальной оптимизации многоэкстремальных функций к настоящему времени достаточно хорошо разработаны и получили широкое распространение. Следует отметить, что эффективность детерминированных алгоритмов существенно ограничена их зависимостью от размерности задачи. В случае большого числа переменных применяют алгоритмы стохастической глобальной оптимизации. К ним относятся, например, популяционные алгоритмы [28, 29]. Вместе с тем чувствительность к выбору параметров алгоритмов этого типа, устанавливаемых пользователем или обусловленных содержанием задачи, во многом определяет скорость сходимости итерационного процесса. Этого недостатка лишен кратный алгоритм столкновения частиц M-PCA (Multi-Particle Collision Algorithm), который основан на алгоритме Метрополиса и входит в число наиболее мощных современных стохастических алгоритмов глобальной оптимизации [30]. Существенно, что применение стохастических алгоритмов глобальной оптимизации требует значительных вычислительных ресурсов. Одним из путей повышения эффективности таких алгоритмов является совершенствование процедуры локального поиска. В работе [31] представлены гибридные алгоритмы, объединяющие генетический алгоритм сканирования пространства переменных и детерминиро-
ванные методы локального поиска; отмечен также ряд недостатков описанных гибридных алгоритмов.
Целью настоящей работы является разработка новых гибридных алгоритмов глобальной недифференцируемой оптимизации, ориентированных на решение задач вычислительной диагностики гидромеханических систем.
Представлены математические модели акустических колебаний двухфазного газожидкостного потока в циркуляционном контуре, сформулирована обратная задача вычислительной диагностики фазового состава теплоносителя по спектральным данным. Предполагается, что регистрируемые данные могут быть неполными, спектры колебаний потока содержат кратные частоты, шумы отсутствуют. Описаны методы недифференцируемой оптимизации, используемые при локальном поиске в гибридных алгоритмах. Процедуры локального поиска основаны на построении сглаживающих аппроксимаций критериальных функций. Приведено описание двух новых гибридных алгоритмов глобальной недифференцируемой оптимизации, а также результаты решения модельной задачи вычислительной диагностики фазового состава теплоносителя в циркуляционном контуре ядерной реакторной установки.
Математические модели. Важной характеристикой объектов, имеющих в своем составе гидравлические системы, является спектр собственных частот колебаний рабочей жидкости в этих системах и ее отдельных контурах. Для разветвленных гидросистем, включающих различные элементы и агрегаты, объединенные в единую систему гидролиниями с потоком рабочей жидкости, методика расчета собственных частот может быть построена, например, на объединенном импедансно-матричном способе получения характеристического уравнения для расчета собственных частот [20]. В этом случае находятся передаточные матрицы В у по отдельным у-м ветвям гидросистемы как результат перемножения передаточных матриц отдельных элементов гидросистемы, представляемых в форме многополюсников, заключенных между узловыми точками на входе и выходе общего участка в каждой ветви. В месте соединения элементов рассчитываются переходные матрицы стыков и ответвлений. Потери напора на разных участках могут быть учтены введением эквивалентных сосредоточенных гидравлических сопротивлений на стыках элементов гидросистемы. Тогда
В у =П В л,
1
где В л — матрицы перехода отдельных элементов, стыков и ответвлений в каждой ветви гидросистемы, имеющие вид квадратной матрицы порядка п:
\\п шк 1•
Здесь Ьшк — элементы матрицы перехода /-го элемента у-й ветви; п = 1 для элементов с сосредоточенными параметрами, представляемых в виде двухполюсников; п = 2 для систем с распределенными параметрами (для четырехполюсников) и дополнительно с учетом тепловых процессов п = 3 (для шестиполюсников).
В узловых точках, т. е. в местах соединения отдельных ветвей гидросистемы, должны выполняться условия баланса расходов и равенства давлений в потоке рабочей жидкости.
При вычислениях следует учитывать, что скорость распространения малых возмущений в гидролиниях и элементах системы зависит как от температуры теплоносителя, так и (в большей степени) от его фазового состава. Приведенная скорость распространения малых возмущений в гидролиниях, например при пузырьковой структуре двухфазного потока, может быть рассчитана по формуле
Здесь а/ — скорость распространения малых возмущений на рассматриваемом участке гидролинии с капельной жидкостью, определяемая по формуле Жуковского с учетом податливости стенок трубопровода; при этом плотность смеси р определяется в виде
где р/, р^ — плотности капельной жидкости и парогазовой фазы соответственно; Р — объемное паросодержание в смеси. При этом следует отметить, что даже малое присутствие свободной парогазовой фазы в потоке рабочей жидкости приводит к существенному снижению скорости звука в ней, в то время как плотность смеси практически остается равной плотности капельной жидкости.
Граничные условия в концевых сечениях незамкнутых ветвей гидросистемы могут быть заданы в форме граничных импедансов, в общем случае зависящих от частоты и имеющих комплексную форму [20].
Используя известные соотношения для расчета эквивалентных передаточных матриц, параллельно и последовательно соединенных контуров гидросистемы, граничные условия и условия в узлах ветвления гидросистемы, можно получить систему уравнений, из которой, приравнивая к нулю ее главный определитель, находится характеристическое уравнение для поиска собственных частот колебаний рабочей жидкости в гидросистеме.
Р = Р/ (1 -Р) + РеР,
Для закольцованных систем типа главных циркуляционных контуров АЭС [21] с эквивалентной матрицей перехода В характеристическое уравнение для отыскания собственных частот имеет общий вид
[В - Е] = 0,
где Е — единичная матрица.
Решение последнего уравнения возможно, например, методом начального параметра с уточнением по методу Ньютона. Расчеты, выполненные таким способом, в диапазоне частот до 50 Гц хорошо
коррелируют с результатами расчетов на основе метода конечных
®
элементов в пакете прикладных программ АКБУБ .
Постановка обратной задачи. Рассматривается обратная задача вычислительной диагностики, которая в рамках выбранной математической модели описывается операторным уравнением
Ах = у, х еX, у еУ,
где X, У — гильбертовы пространства; А — компактный линейный оператор, действующий из X в У. Правая часть возмущенного операторного уравнения представляет приближенные входные данные
у , определенные из эксперимента. Предполагается, что погрешность задания входной информации 8 известна и имеет место < 8. Требуется определить устойчивые приближенные ре-
*6" *
g
шения по заданной приближенной информации * . Существенно, что во многих приложениях обратные задачи являются некорректно поставленными. Далее реализуется подход, основанный на методе регуляризации [17, 18].
Приближенное решение рассматриваемой обратной задачи связано с поиском минимума функционала Тихонова:
g
xa = arg min Ja (x), а > 0.
xeX
Здесь xa — регуляризованное решение уравнения Ax = * с параметром регуляризации а; при этом минимизируемый функционал определен в виде
8 2
Ja (x) = Ax - *
+ а || x X
где
g2
Ax - * — функционал невязки (квадрат нормы в простран-
Y
X
«1|2
x X — стабилизирующий функционал.
Реализуемый метод называют сходящимся методом регуляризации, если выполнено условие
sup
II у8-y IM
Ray8 - x ^ 0 при 8 ^ 0.
Здесь введено семейство линейных ограниченных операторов Ra : Y ^ X с параметром регуляризации а > 0, таких, что
lim RaAx = x Ух e X.
a^0
При использовании заданной априори информации х e X1 с X погрешность решения в наихудшем случае составляет [16]
где
А(5,XьR,) = sup {|| R,y5 - х 11: x eXьy eY, || y5 - y
Кроме того, для любого Rа (0) = 0 имеет место А(5,X!,Ra) > Ао(5,X!),
А 0(5, Xj) = sup HI х 11: х e X1, || Ax | |<5].
<8
Задача вычислительной диагностики системы как обратная задача на собственные значения связана с поиском вектора переменных управления, при котором первые N собственных частот (или соответствующих им собственных значений) модели совпадают с составляющими некоторого заданного ограниченного спектра или достаточно близки к ним. Для оценки уровня рассогласования сравниваемых характеристик объекта используется векторный способ описания. Так как информация о формах колебаний объекта зачастую отсутствует или является существенно неполной, далее рассматривается только рассогласование между частотными составляющими нормального и заданного спектров. Возможные подходы основаны на минимизации квадратичной функции рассогласования или минимизации максимальной из функций рассогласования спектральных составляющих. Так, для попарно сравниваемых спектральных составляющих может быть построено следующее конечное множество критериев рассогласования
f (x) =Ci (x) -Ci (x)
e X с Rn, i e J,
где С ^ (х), С (х) — собственные значения, относящиеся к исходному (текущему) и заданному спектрам; х — вектор переменных управления; X — допустимая область; п — размерность задачи; Яп —
n -мерное вещественное линейное пространство; J = {1, ... , nj. Необходимо найти такой вектор переменных управления, который приводит к наименьшим различиям между сравниваемыми спектрами, т. е. следует произвести настройку модели объекта на заданный спектр. Задача диагностирования формулируется в следующем виде: определить вектор переменных управления х е X , который минимизирует максимальное значение критерия рассогласования, т. е. требуется найти
min max {f (x)j. (P)
хеХ cRn iel
Решением сформулированной дискретной минимаксной задачи
" * / * *\т " (P) является такой вектор х =(хь ..., хп) , принадлежащий множеству допустимых значений, при котором скалярная критериальная функция f (х) = max {f1(х), ..., fN(x)J принимает минимальное значение. В случае, когда f (х*) = 0, спектр частот настраиваемой модели полностью совпадает с заданным спектром по N низшим частотам. Последнее условие вследствие неполноты экспериментальных данных и погрешностей, полученных при измерениях, не выполняется. Ниже рассматривается регуляризованная задача (P)a с многоэкстремальной не всюду дифференцируемой критериальной функцией f (х) и параметром регуляризации a > 0.
В обобщение постановок экстремальных задач вычислительной диагностики формулируется следующая задача глобальной оптимизации: требуется найти
f (х*) = min f (х), (1)
хеХ cRn
где
X ={ х е D : g (х) < 0, i е I j; (2)
D = {х е Rn : üj < ху < bj, j е J j, (3)
f (х) — целевая функция; gi (х) — функции ограничений задачи, i еI; I = {1, ..., mj — конечное множество индексов; D — область
поиска; х* — глобальное решение. Функции f (х), gi (х), i е I, задачи (1)—(3) полагаются непрерывными липшицевыми. Предполагается также, что действительная функция f :Rn ^ R является многоэкстремальной, не всюду дифференцируемой и для нее задана вычислительная процедура, позволяющая определять значения функции в точ-
ках допустимой области. Необходимо также учесть возможную высокую трудоемкость вычисления критериальных функций, что может потребовать значительных вычислительных ресурсов.
Методы сглаживающих аппроксимаций и локальная минимизация критериальных функций. Рассмотрим задачу (1)—(3), ограничившись поиском локального решения. Предварительно исследуется задача поиска минимума действительной функции
f :Яп ^Я, определенной в виде
f(x) = max {фу(x)}, i e Im = {1, ..., M}. (4)
xeX cRn
Здесь X — допустимое множество; предполагается, что все функции фу (x), i e IM, выпуклы и непрерывно дифференцируемы.
Задачи, формулируемые в минимаксной форме, относят к классу задач недифференцируемой оптимизации. Для их решения применяют специальные методы, например модифицированный метод сопряженных градиентов, метод гиперболической сглаживающей функции и др. [25, 26]. Рассматриваемый далее подход основан на построении сглаживающих аппроксимаций критериальных функций с последующим применением эффективных методов, разработанных для задач дифференцируемой оптимизации. Подход предполагает замену каждой недифференцируемой функции некоторой ее аппроксимацией, которая была бы выпуклой и дифференцируемой в области допустимых значений переменных управления. Далее в качестве процедур локального поиска гибридных алгоритмов рассматривают метод гиперболической сглаживающей функции [25, 26], а также метод построения двухпараметрических сглаживающих (p, q) -аппроксимаций критериальной функции [27].
Применительно к задаче недифференцируемой оптимизации подход с использованием гиперболического сглаживания основан на введении новой критериальной функции [26]
F(x, t) = t max{0, f (x) -t}
ieI0
с последующей аппроксимацией функции F(x, t) в виде
f (x) -1+ V(f (x) -1)2 +T2
Фх(x, t) = t
iel0 2
где т> 0 — параметр точности; / е Я. Следует отметить, что функции fi (х), / е 10, должны быть непрерывно дифференцируемыми; параметр I определяется числом ограничений.
Предложение [26]. Для любых х е Rn и t е R справедлива оценка 0 <ФХ (х, t) - F(х, t) < у.
Далее рассматриваемая задача заменяется последовательностью задач минимизации функций Фх^ (х, f (х)) и имеет место тк ^ 0 при
k ^да. В работе [26] предложен алгоритм решения указанной последовательности задач и доказана его сходимость.
Второй подход состоит в следующем. Целевую функцию (4) можно определить в эквивалентной форме
f (х) = Ф1( х) + У(Ф2( х) -Ф1( х) + (5) + У {... + У [Фм-1 (х) - Фм-2 (х) + У (Фм (х) + Фм -i(х)).]},
где
У(Фг(х)) = max {0, Фг (х)}, i е 1м. (6)
хеХ cRn
Основная идея рассматриваемого метода состоит в том, чтобы каждую функцию у (фг- (х)), i е IM , входящую в (5), заменить некоторой гладкой функцией, построить сглаженную приближенную целевую функцию, а затем применить эффективные методы гладкой минимизации. При возрастании точности аппроксимации функций (6) имеет место сходимость приближенного решения к точному.
Существенно, что уже в одномерном случае функция у(х) =
= max {0, х} в точке х = 0 дифференцируема только по направлени-
хеХ с R 1 '
ям. Возможен следующий подход: на числовой оси выделяется отрезок [p, q], содержащий точку, в которой функция у(х) имеет указанную особенность, и на этом отрезке исходная функция заменяется некоторой приближенной функцией, выпуклой и дифференцируемой в каждой точке по построению. Пусть выбраны числа p < 0 и q > 0 . Вводится двухпараметрическая аппроксимация функции у : R ^ R
0, х < p;
У(p, q, х) = Ь(p, q, х), p < х < q; х, х > q.
Здесь p, q — параметры аппроксимации, определяющие соответственно левую и правую границы отрезка [p, q], на котором задана сглаживающая функция s (p, q, х) . Приближенная функция У (p, q, х)
совпадает с исходной у(х) всюду, за исключением отрезка [p, q].
Потребуем, чтобы функция s (p, q, х) была выпуклой и по крайней
мере один раз дифференцируемой на [p, q]. При этом
s(p, q, 0) = -pn(p, q), где p, q) определяется свойствами сглаживающей функции.
Теорема 1 [27]. Пусть х* е R" и х еR" — суть точки минимума для f (х) и f (p, q, х) соответственно. Тогда 0 < f (p, q, X)-
-f(x*) <-p min {1, (M - 1)^(p, q)}.
хеХ cRn 1 П
Теперь, с использованием сглаживающих аппроксимаций, могут быть сформулированы необходимые условия Каруша — Куна — Таккера. Пусть рассматривается задача оптимизации со смешанными функциональными ограничениями
f (х) ^ min, х еD; (7)
D = { х е Rn : ^(х) = 0, О(х) < 0 }, (8)
где f : Rn ^ R — заданная функция; F : Rn ^ Rl и G : Rn ^ Rm — заданные отображения; l — число ограничений в форме равенств; m — число ограничений-неравенств.
Пусть х е D — локальное решение задачи (7), (8), причем функция f дифференцируема в точке х , отображения F и G удовлетворяют условиям гладкости. Вводится обобщенная функция Лагранжа задачи (7), (8)
L0: Rn х R х Rl x Rm ^ R; L)(х, ^0, X, ц) = hf (х) + (X, F(х)) + (ц, G(х)).
С использованием сглаживающих аппроксимаций можно определить
L0(p, q; х, X0, X, ц) = X0f (p, q; х) + (X,F(p, q; х)) + (ц,G(p, q; х)^, при этом
sL
— (p, q; х X ц) = X0f'(p, q; х) + (F'(p, q; х))г X + (G'(p, q; х))г Ц
ох
х е Rn, X0 > 0, X е Rl, ц е Rm.
Теорема 2. Пусть выбраны параметры р, д, функция / : Я" ^ Я
и отображение О : Я" ^ Ят , дифференцируемое в точке х : Я", а
отображение Р : Я" ^ Я1 дифференцируемо в некоторой окрестности этой точки, причем его производная непрерывна в точке х .
Тогда если х является локальным решением задачи (7), (8), то
найдутся число Х0 > 0 и элементы X е Я1 и р еЯ+*, не равные нулю одновременно, и такие, что
Жо / л" л" \ /л
(р, д; х, Хо, X, р) = 0,
ох
^ °(p, д; = 0.
Д о к а з а т е л ь с т в о. Доказательство получается прямой ссылкой на теорему 3 [31, с. 42] и теорему 1.
Далее рассмотрим важный практический случай задачи вычислительной диагностики — задачу минимизации (1), (3) для случая простых ограничений (на переменные управления). Требуется найти
Ш1П
х
{7 (р, д, х): ау < ху < Ъу, у е 3}, (9)
где / (р, д, х) — выпуклая функция; допустимая область X совпадает с областью поиска В. Вспомогательная задача квадратичного программирования с вектором „ е Я" формулируется в виде: найти
Щ^А „ + 2(„у)2: ау < ху < Ъ,, у е (10)
Решение задачи (10) дает „у, после чего определяются множители Каруша — Куна — Таккера и- и и+, соответствующие неравенствам ху + „у - Ъу < 0 и - ху - „у + Ъу < 0, у е 3 . Существенно, что
для минимизируемой в задаче (9) целевой функции должны выполняться условия [27]:
с¥ (р, д, х) + -—^- " и у + иу + „у = 0, у е 3;
и+ >0, и+ (аз- хз- „у ) =
> 0, и- (Х] + „у - Ъу ) = ^ ] е 3 .
Пусть требуется решить задачу (9). Выбраны числа а у, Ьу, у е J,
а также число Р, 0 < Р < 1, и параметры аппроксимации р < 0, q > 0. Алгоритм минимизации включает следующие основные шаги.
0. Выбрать точку х° , а у < х0 < Ьу, у е J.
1. Если точка хк уже построена, то вычислить вектор м>к = w (хк |.
2. Определить первое значение г = 0, 1, ..., при котором для а = (1/ 2)г будет выполнено неравенство
7(р, ^ хк +аwk)</(р, ^ хк)-ра
wk
если такое г = г0 найдено, то положить а к = 2-г°, хк+1 = хк +аkWk.
Перейти к шагу 1.
3. Критерий останова: wk = 0 .
Локальную сходимость алгоритма минимизации при использовании сглаживающих аппроксимаций критериальных функций для случая простых ограничений устанавливает приведенное далее утверждение.
Теорема 3 [27]. Пусть выбраны параметры р < 0, q > 0. Если
числа а у, Ьу, у е J, конечны и градиент функции 7 (р, q, х) удовлетворяет условию Липшица, то во всякой предельной точке последовательности хк, к = 0, 1, ... , удовлетворяются необходимые условия минимума.
Гибридные алгоритмы. Структуры алгоритмов глобальной минимизации построены на основе стохастического алгоритма М-РСА [20], объединенного с процедурами поиска локальных минимумов не всюду дифференцируемых функций. Работа современного алгоритма глобальной оптимизации М-РСА основана на использовании аналогии с физическими процессами абсорбции и рассеяния частиц при ядерных реакциях. В простейшей версии алгоритма для исследования области поиска используется одна частица. На начальном шаге выбирается пробное решение (ОЫ_Соп11§), которое затем модифицируется посредством стохастического возмущения (Рег!игЬа1;юп(.)), что позволяет найти новое решение (Ке—_СопП§). С помощью функции БкпеввО дается сравнительная оценка нового и предыдущего решений, на основании которой новое решение может быть принято или отвергнуто. Если новое решение отвергнуто, то происходит переход к функции Беайег-т§(.), реализующей схему Метрополиса. Для сканирования области, перспективной на минимум, применяются функции Рег1игЬа1юп(.) и 8ша11_Рег1игЬа1юп(.). Новое решение принимается, если оно лучше предыдущего (абсорбция); если найденное решение хуже предыдущего, то происходит переход в отдаленную область пространства поиска (рас-
сеяние), что позволяет преодолевать локальные минимумы. Эффективность описанного поиска глобального решения алгоритмом значительно повышается за счет одновременного использования большого числа частиц. Такой подход реализует алгоритм M-PCA, который непосредственно ориентирован на применение в среде параллельных вычислений. Наилучшее решение определяется с учетом данных обо всех частицах, участвующих в процессе. Единственным задаваемым параметром для алгоритма M-PCA является число итераций.
Предложены гибридные алгоритмы, интегрирующие алгоритм M-PCA, и детерминированные алгоритмы локальной минимизации. В работе [26] реализован метод гиперболической сглаживающей функции. Первый гибридный алгоритм объединяет стохастический алгоритм M-PCA сканирования пространства переменных и детерминированный градиентный локальный поиск GHS, использующий метод гиперболической сглаживающей функции. В работе [27] представлен двухпараметрический метод построения сглаживающих аппроксимаций не всюду дифференцируемых функций и предложен вариант метода линеаризации LMS со сглаживанием. Второй гибридный алгоритм объединяет алгоритм M-PCA глобального сканирования и детерминированный метод LMS локального поиска. Результирующие гибридные алгоритмы M-PCAGHS и M-PCALMS реализованы в виде прикладного программного обеспечения. Рассмотрим фрагмент псевдокода гибридного алгоритма M-PCALMS.
1. Generate an initial solution Old_Config Best_Fitness = Fitness (Old_Config) Update Blackboard
For n = 0 to # of particles For n = 0 to # of iterations Update Blackboard Perturbation ( )
If Fitness (New_Config) > Fitness (Old_Config) If Fitness (New_Config) > Best_Fitness
Best_Fitness := Fitness (New_Config) End If
Old_Config := New_Config Exploration ( ) Else
Scattering ( ) End If
End For End For
2. Exploration ( )
For n = 0 to # of iterations Small_Perturbation ( )
Local search
using Linearization Method with Smoothing Approximations Check stopping criterion: Find global solution Best Fitness Else continue
If Fitness (New_Config) > Best_Fitness Best_Fitness := Fitness (New_Config) End If
Old_Config := New_Config End For
Return
3. Scattering ( )
Pscatt = 1 - ( Fitness (New_Config)) / (Best_Fitness) If Pscatt > random(0, 1)
Old_Config := random solution Else
Exploration ( ) End If Return
В состав алгоритма M-PCALMS входят также стандартные процедуры Perturbation(.) и Small_Perturbation(.) [31].
Численный пример. Рассматривается модельная задача вычислительной диагностики фазового состава теплоносителя в замкнутом циркуляционном контуре ядерной реакторной установки [23, 32, 33]. Диагностирование проводится по косвенной информации о частотах акустических колебаний в газожидкостном потоке, полученных при измерениях. Переменными модели потока являются относительные значения скорости звука xi,% на участках контура, соответствующих: зоне нагрева теплоносителя в напорном баке системы компенсации объема (СКО) (x^; выходному объему реактора (x2); активной зоне реактора (Х3) ; проточной части главного циркуляционного насоса циркуляционной петли с СКО (Х4). При отсутствии в теплоносителе второй фазы представленный далее нормальный спектр ш j,
j = 1, 10, соответствует максимальным значениям скорости звука на выделенных участках контура.
Нормальный и аномальный спектры частот колебаний теплоносителя
j ..................... 1 2 3 4 5 6 7 8 9 10
шj, Гц............ 0,89 6,77 9,82 15,44 15,96 18,94 24,56 26,69 27,07 30,52
ш* , Гц............ 0,81 6,76 9,33 15,15 15,74 18,80 20,79 26,63 26,89 29,32
В модельной задаче аномальный спектр ш*, у = 1, 10, получен
при наличии двухфазной смеси в напорном баке СКО, в выходном объеме и активной зоне реактора, а также в проточной части главного циркуляционного насоса циркуляционной петли с СКО, при этом
х* = 74,5 %; х2 = 87,25 %; х3 = 81,5 %; х4 = 91,0 % . Критериальная функция определена с учетом десяти низших спектральных составляющих. Для решения задачи вычислительной диагностики используется гибридный алгоритм М-РСАЬМБ. После определения области переменных модели, содержащей глобальный минимум, завершающие итерации алгоритма проводятся с использованием градиентной информации для сглаживающих аппроксимаций критериальной функции. Сходимость решения иллюстрируют рис. 1 и 2.
Х\, ХЪ *3> 9,50£+01
9,00£+01
7,00£+01
Рис. 1. Изменение значений переменных управления х, с ростом числа итераций Ni
1,00£+04
Рис. 2. Зависимость критериальной функции 7 (х) и нормы вектора N от числа итераций Ni
По завершении локального поиска получено приближенное решение: х* « 74,48 % ; х2 - 87,35 %; х33 - 81,77 % ; х4 « 90,12 %. Относительная погрешность определения значений переменных модели
не превышает 1,0 % при точности настройки спектра частот порядка 0,01 Гц. Итак, по завершении настройки спектра частот математической модели газожидкостного потока на заданный аномальный спектр установлено, что имеет место появление второй фазы в потоке теплоносителя на выделенных участках циркуляционного контура.
Выводы. Разработана математическая модель акустических колебаний в двухфазном потоке теплоносителя, циркулирующем в замкнутом контуре ядерной реакторной установки. Предложен подход к решению обратных задач вычислительной диагностики фазового состава теплоносителя с использованием новых гибридных алгоритмов глобальной оптимизации. Исследование пространства переменных модели проводится стохастическим методом, реализуемым кратным алгоритмом столкновения частиц. При локальном поиске в гибридных алгоритмах градиентная информация определяется для сглаживающих аппроксимаций не всюду дифференцируемых критериальных функций. Модельное диагностирование показало возможность идентификации аномалий фазового состава теплоносителя в контуре реакторной установки с достаточной для приложений точностью.
Работа выполнена при финансовой поддержке Министерства образования и науки РФ (грант Президента РФ по поддержке научных исследований ведущих научных школ РФ, код НШ-4058.2014.8).
ЛИТЕРАТУРА
[1] Romenski E., Drikakis D., Toro E. Conservative models and numerical methods for compressible two-phase flow. J. of Scientific Computing, 2010, vol. 42, no. 1, pp. 68-95.
[2] Yao L., Zhu C. Free boundary value problem for a viscous two-phase model with mass-dependent viscosity. J. of Differential Equations, 2009, vol. 247, no. 10, pp. 2705-2739.
[3] Bonilla J., Yebra L.J., Dormido S. Chattering in dynamic mathematical two-phase flow models. Applied Mathematical Modelling, 2012, vol. 36, no. 5, pp. 2067-2081.
[4] Khatri S., Tornberg A.-K. An imbedded boundary method for soluble surfactants with interface tracking for two-phase flows. J. of Computational Physics, 2014, vol. 256, pp. 768-790.
[5] Evje S. Global weak solutions for a compressible gas-liquid model with wellformation interaction. J.of Differential Equations, 2011, vol. 251, no. 8, pp. 2352-2386.
[6] Yao L., Zhu C.J. Existence and uniqueness of global weak solution to a two-phase flow model with vacuum. Mathematische Annalen, 2011, vol. 349, no. 4, pp. 903-928.
[7] Vallée C., Hohne T., Prasser H.-M., Suhnel T. Experimental investigation and CFD simulation of horizontal stratified two-phase flow phenomena. Nuclear Engineering and Design, 2008, vol. 238, no. 3, pp. 637-646.
[8] Poullikas A. Effects of two-phase flow on the performance of nuclear reactor cooling pumps. Progress in Nuclear Energy, 2003, vol. 42, no. 1, pp. 3-10.
[9] Yang X., Schiegel J.P., Liu Y., Paranjape S., Hibiki T., Ishii M. Experimental study of interfacial area transport in air-water two phase flow in a scaled 8*8 BWR rod bundle. J. of Multiphase Flow, 2013, vol. 50, pp. 16-32.
[10] Семченков Ю.М., Мильто В.А., Шумский Б.Е. Внедрение методики контроля кипения теплоносителя в активной зоне ВВЭР-1000 в систему внут-риреакторной диагностики. Атомная энергия, 2008, т. 105, № 2, с. 79-82.
[11] Корнюшин Ю.П., Егупов Н.Д., Корнюшин П.Ю. Идентификация параметров исполнительных устройств регуляторов паровой энергетической турбины с использованием аппарата матричных операторов. Математическое моделирование и численные методы, 2015, № 2(6), с. 73-86.
[12] Kumbaro A. Simplified eigenstruture decomposition solver for the simulation of two-phase flow systems. Computers and Fluids, 2012, vol. 64, pp. 19-33.
[13] Zeidan D., Slaouti A. Validation of hyperbolic model for two-phase flow in conservative form. International J. Computational Fluid Dynamics, 2009, vol. 23, no. 9, pp. 623-641.
[14] Oberkampf W.L., Barone M.F. Measures of agreement between computation and experiment: Validation metrics. J. of Computational Physics, 2006, vol. 217, no. 1, pp. 5-36.
[15] Manes J.P., Espinoza V.H.S., Vicent S.C., Böttcher M., Stieglitz R. Validation of NEPTUNE-CFD two-phase flow models using experimental data. Science and Technology of Nuclear Installations, 2014, vol. 2014, article ID 185950, 19 p.
[16] Lippert R.A. Fixing multiple eigenvalues by a minimal perturbation. Linear Algebra and its Applications, 2010, vol. 432, pp. 1785-1817.
[17] Kirsch A. An introduction to the mathematical theory of inverse problems. 2nd edition. New York et al.: Springer, 2011, 308 p.
[18] Chu D., Lin L., Tan R.C.E., Wei Y. Condition numbers and perturbation analysis for the Tikhonov regularization of discrete ill-posed problems. Numerical Linear Algebra with Applications, 2011, vol. 18, no. 1, pp. 87-103.
[19] Гончарский А.В., Романов С.Ю. О двух подходах к решению коэффициентных обратных задач для волновых уравнений. Журнал вычислительной математики и математической физики, 2012, т. 52, № 2, с. 263-269.
[20] Гликман Б.Ф. Математические модели пневмогидравлических систем. Москва, Наука, 1986, 386 с.
[21] Кинелев В.Г., Сулимов В.Д., Шкапов П.М. Дигностирование гидросистемы на основе анализа ее частотного спектра. Изв. РАН. Энергетика, 1998, № 6, с. 112-119.
[22] Bai Z.-J., Ching W.-K. A smoothing Newton's method for the construction of a damped vibrating system from noisy test eigendata. Numerical Linear Algebra with Applications, 2009, vol. 16, no. 2, pp. 109-128.
[23] Kinelev V.G., Shkapov P.M., Sulimov V.D. Application of global optimization to VVER-1000 reactor diagnostics. Progress in Nuclear Energy, 2003, vol. 43, no. 1-4, pp. 51-56.
[24] Medeiros J.A.C., Schirru R. Identification of nuclear power plant transients using the Particle Swarm Optimization algorithm. Annals of Nuclear Energy, 2008, vol. 35, no. 4, pp. 576-582.
[25] Chen X. Smoothing methods for nonsmooth, nonconvex minimization. Mathematical Programming, 2012, vol. 134, no. 1, pp. 71-99.
[26] Bagirov A.M., Al Nuaimat A., Sultanova N. Hyperbolic smoothing function method for minimax problems. Optimization: A Journal of Mathematical Programming and Operations Research, 2013, vol. 62, no. 6, pp. 759-782.
[27] Сулимов В. Д. Локальная сглаживающая аппроксимация в гибридном алгоритме оптимизации гидромеханических систем. Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки, 2010, № 3, с. 3-14.
[28] Карпенко А.П. Современные алгоритмы поисковой оптимизации. Алгоритмы, вдохновленные природой. Москва, Изд-во МГТУ им. Н.Э. Баумана, 2014, 446 с.
[29] Wang Y., Garcia A. Interactive model-based search for global optimization. J. of Global Optimization, 2015, vol. 61, no. 3, pp. 479-495.
[30] Luz E.F.P., Becceneri J.C., de Campos Velho H.F. A new multi-particle collision algorithm for optimization in a high performance environment. J. of Computational Interdisciplinary Sciences, 2008, vol. 1, pp. 3-10.
[31] Voglis C., Parsopoulos K.E., Papageorgiou D.G., Lagans I.E., Vrahatis M.N. MEMPSODE: A global optimization software based on hybridization of population-based algorithms and local searches. Computer Physics Communications, 2012, vol. 183, no. 2, pp. 1139-1154.
[32] Измаилов А.Ф., Солодов М.В. Численные методы оптимизации. Москва, Физматлит, 2005, 304 с.
[33] Sulimov V.D., Shkapov P.M. Application of hybrid algorithms to computational diagnostic problems for hydromechanical systems. J. of Mechanics Engineering and Automation, 2012, vol. 2, no. 12, pp. 734-741.
Статья поступила в редакцию 19.08.2015
Ссылку на эту статью просим оформлять следующим образом: Сулимов В. Д., Шкапов П.М. Гибридные методы вычислительной диагностики двухфазного потока в циркуляционном контуре. Математическое моделирование и численные методы, 2015, № 3, с. 68-88.
Сулимов Валерий Дмитриевич родился в 1950 г., окончил МВТУ им. Н.Э. Баумана в 1973 г. Старший преподаватель кафедры «Теоретическая механика» им. Н.Е. Жуковского МГТУ им. Н.Э. Баумана. Автор более 60 научных работ в области математического моделирования и оптимизации динамических систем. e-mail:[email protected]
Шкапов Павел Михайлович родился в 1954 г., окончил МВТУ им. Н.Э. Баумана в 1977 г. Д-р техн. наук, заведующий кафедрой «Теоретическая механика» им. Н.Е. Жуковского МГТУ им. Н.Э. Баумана. Автор более 100 работ по математическому моделированию и оптимизации механических и гидромеханических систем. e-mail:[email protected]
Hybrid methods of computer diagnosis of two-phase flow in the circulation loop
© V.D. Sulimov, P.M. Shkapov
Bauman Moscow State Technical University, Moscow, 105005, Russia
The article considers the problems of coolant flow computational diagnostics in a closed circulation loop. The mathematical models of acoustic waves in two-phase flow are developed. Indirect diagnostic information, contained in the flow vibrational spectra recorded by regular systems is used. The inverse eigenvalue problem is formulated. Solving it the optimization approach is implemented. It is supposed that partial criteria are presented by continuous, Lipschitz, not everywhere differentiable, multi-extremal functions. Search of global solutions was performed using a new hybrid algorithms integrating stochastic algorithm of variable space viewing and deterministic methods of local search. A numerical example of model diagnosing the phase composition of the coolant in the circulation loop of nuclear reactor plant is presented.
Keywords: two-phase flow, inverse problem, regularization, global optimization, Metropolis algorithm, hybrid algorithm.
REFERENCES
[1] Romenski E., Drikakis D., Toro E. Conservative models and numerical methods for compressible two-phase flow. J. of Scientific Computing, 2010, vol. 42, no. 1, pp. 68 -95.
[2] Yao L., Zhu C. Free boundary value problem for a viscous two-phase model with mass-dependent viscosity. J. of Differential Equations, 2009, vol. 247, no. 10, pp. 2705-2739.
[3] Bonilla J., Yebra L.J., Dormido S. Chattering in dynamic mathematical two-phase flow models. Applied Mathematical Modelling, 2012, vol. 36, no. 5, pp. 2067-2081.
[4] Khatri S., Tornberg A.-K. An imbedded boundary method for soluble surfactants with interface tracking for two-phase flows. J. of Computational Physics, 2014, vol. 256, pp. 768-790.
[5] Evje S. Global weak solutions for a compressible gas-liquid model with wellformation interaction. J. of Differential Equations, 2011, vol. 251, no. 8, pp. 2352-2386.
[6] Yao L., Zhu C.J. Existence and uniqueness of global weak solution to a two-phase flow model with vacuum. Mathematische Annalen, 2011, vol. 349, no. 4, pp. 903-928.
[7] Vallée C., Höhne T., Prasser H.-M., Sühnel T. Experimental investigation and CFD simulation of horizontal stratified two-phase flow phenomena. Nuclear Engineering and Design, 2008, vol. 238, no. 3, pp. 637-646.
[8] Poullikas A. Effects of two-phase flow on the performance of nuclear reactor cooling pumps. Progress in Nuclear Energy, 2003, vol. 42, no. 1, pp. 3-10.
[9] Yang X., Schiegel J.P., Liu Y., Paranjape S., Hibiki T., Ishii M. Experimen-tal study of interfacial area transport in air-water two phase flow in a scaled 8*8 BWR rod bundle. J. of Multiphase Flow, 2013, vol. 50, pp. 16-32.
[10] Semchenkov Yu.M., Milto V.A., Shumskiy B.E. Atomnaya energiya - Nuclear Power, 2008, vol. 105, no. 2, pp. 79-82.
[11] Kornyushin Yu.P., Egupov N.D., Kornyushin P.Yu. Matematicheskoe modeliro-vanie i chislennye menody - Mathematical Modeling and Numerical Methods, 2015, no. 2(6), pp. 73-86.
[12] [Kumbaro A. Simplified eigenstruture decomposition solver for the simula-tion of two-phase flow systems. Computers and Fluids, 2012, vol. 64, pp. 19-33.
[13] Zeidan D., Slaouti A. Validation of hyperbolic model for two-phase flow in conservative form. International J. Computational Fluid Dynamics, 2009, vol. 23, no. 9, pp. 623-641.
[14] Oberkampf W.L., Barone M.F. Measures of agreement between computation and experiment: Validation metrics. J. of Computational Physics, 2006, vol. 217, no. 1, pp. 5-36.
[15] Manes J.P., Espinoza V.H.S., Vicent S.C., Böttcher M., Stieglitz R. Validation of NEPTUNE-CFD two-phase flow models using experimental data. Science and Technology of Nuclear Installations, 2014, vol. 2014, article ID 185950, 19 p.
[16] Lippert R.A. Fixing multiple eigenvalues by a minimal perturbation. Linear Algebra and its Applications, 2010, vol. 432, pp. 1785-1817.
[17] Kirsch A. An introduction to the mathematical theory of inverse problems. 2nd edition. New York et al., Springer, 2011, 308 p.
[18] [Chu D., Lin L., Tan R.C.E., Wei Y. Condition numbers and perturbation analysis for the Tikhonov regularization of discrete ill-posed problems. Numerical Linear Algebra with Applications, 2011, vol. 18, no. 1, pp. 87-103.
[19] Goncharskiy A.V., Romanov S.Yu. Zhurnal vychislitelnoy matematiki i ma-tematicheskoi fiziki - Journal of Computational Mathematics and Mathematical Physics, 2012, vol. 52, no. 2, pp. 263-269.
[20] Glikman B.F. Matematicheskie modeli pnevmogidravlicheskikh system [Mathematical Models of Pneumohydraulic Systems]. Moscow, Nauka Publ., 986, 386 p.
[21] Kinelev V.G., Sulimov V.D., Shkapov P.M. Izvestiya RAN.Energetika - Proceedings of the RAS. Power Engineering, 1998, no. 6, pp. 112-119.
[22] Bai Z.-J., Ching W.-K. A smoothing Newton's method for the construction of a damped vibrating system from noisy test eigendata. Numerical Linear Algebra with Applications, 2009, vol. 16, no. 2, pp. 109-128.
[23] Kinelev V.G., Shkapov P.M., Sulimov V.D. Application of global optimization to VVER-1000 reactor diagnostics. Progress in Nuclear Energy, 2003, vol. 43, no. 1-4, pp. 51-56.
[24] Medeiros J.A.C., Schirru R. Identification of nuclear power plant transients using the Particle Swarm Optimization algorithm. Annals of Nuclear Energy, 2008, vol. 35, no. 4, pp. 576-582.
[25] Chen X. Smoothing methods for nonsmooth, nonconvex minimization. Mathematical Programming, 2012, vol. 134, no. 1, pp. 71-99.
[26] Bagirov A.M., Al Nuaimat A, Sultanova N. Hyperbolic smoothing function method for minimax problems. Optimization: A Journal of Mathematical Programming and Operations Research, 2013, vol. 62, no. 6, pp. 759-782.
[27] Sulimov V.D. Vestnic MGTU im. N.E. Baumana. Seria Estestvennye nauki - Herald of the Bauman Moscow State Technical University. Series: Natural Sciences, 2010, no. 3, pp. 3-14.
[28] Karpenko A.P. Sovremennye algoritmy poiskovoy optimizatsii. Algoritmy, vdokh-novlennye prirodoy [Modern Algorithms of Search Optimization. Algorithms Inspired by Nature]. Moscow, BMSTU Publ., 2014, 446 p.
[29] Wang Y., Garcia A. Interactive model-based search for global optimization. J. of Global Optimization, 2015, vol. 61, no. 3, pp. 479-495.
[30] Luz E.F.P., Becceneri J.C., de Campos Velho H.F. A new multi-particle collision algorithm for optimization in a high performance environment. J. of Computational Interdisciplinary Sciences, 2008, vol. 1, pp. 3-10.
[31] Voglis C., Parsopoulos K.E., Papageorgiou D.G., Lagaris I.E., Vrahatis M.N. MEMPSODE: A global optimization software based on hybridization of population-based algorithms and local searches. Computer Physics Communications, 2012, vol. 183, no. 2, pp. 1139-1154.
[32] Izmailov A.F., Solodov M.V. Chislennye metody optimizatsii [Numerical Optimization Techniques]. Moscow, Fismatlit Publ., 2005, 304 p.
[33] Sulimov V.D., Shkapov P.M. Application of hybrid algorithms to computational diagnostic problems for hydromechanical systems. J. of Mechanics Engineering and Automation, 2012, vol. 2, no. 12, pp. 734-741.
Sulimov V.D. (b.1950) graduated from Bauman Moscow Higher Technical School in 1973. Senior lecturer of the Theoretical Mechanics Department named after N.E. Zhukovsky at Bauman Moscow State Technical University. Author of over 60 research publications in the field of mathematical modeling and dynamic system optimization. e-mail:[email protected]
Shkapov P.M. (b. 1954) graduated from Bauman Moscow Higher Technical School in 1977. Dr. Sci. (Eng.), Head of the Theoretical Mechanics Department named after N.E. Zhukovsky at Bauman Moscow State Technical University. Author of over 100 research publications in the field of mathematical modeling and optimization of mechanical and hydromechanical systems. e-mail:[email protected]
