CC BY
41
5. Welch P.D. The Use of Fast Fourier Transform for the Estimation of Power Spectra: A Method Based on Time Averaging Over Short, Modified Periodograms // IEEE Trans. Audio Electroacoust. - 1967. - V. AU-15. - P. 70-73.
6. Percival D.B. and Walden A.T. Spectral Analysis for Physical Applications: Multitaper and Conventional Univariate Techniques. Cambridge: Cambridge University Press, 1993.
7. Бондарев В.Н., Трестер Г., Чернега В.С. Цифровая обработка сигналов: методы и средства. -Севастополь: СевГТУ, 1999. - 398 с.
Демин Анатолий Владимирович - Санкт-Петербургский государственный университет информационных
технологий, механики и оптики, доктор технических наук, профессор, dav_60@mail.ru
Войтюк Татьяна Евгеньевна - Санкт-Петербургский государственный университет информационных
технологий, механики и оптики, аспирант, tanya_4ever@mail.ru Климанов Виталий Александрович - Санкт-Петербургский государственный университет информационных
технологий, механики и оптики, кандидат технических наук, доцент, v_klimanov@mail.ru
УДК 517.958:532.5
РЕАЛИСТИЧНОЕ МОДЕЛИРОВАНИЕ СВОБОДНОЙ ВОДНОЙ ПОВЕРХНОСТИ НА АДАПТИВНЫХ СЕТКАХ ТИПА «ВОСЬМЕРИЧНОЕ ДЕРЕВО» К.Д. Никитин
Предлагается эффективная технология моделирования течения вязкой несжимаемой жидкости со свободной границей. В технологии объединены проекционный метод для решения уравнений Навье-Стокса и метод функции уровня и частиц для работы со свободной поверхностью. Расчет проводится на адаптивно сгущающихся гексаэдральных сетках типа «восьмеричное дерево».
Ключевые слова: вычислительная технология, течения со свободной поверхностью, компьютерная графика.
Введение
Моделирование течений несжимаемой жидкости привлекает большой интерес в научном сообществе. Рассматриваемые задачи включают моделирование морских волн, заливание и обтекание объектов, падение капель, образование брызг и многое другое. Моделирование подобных явлений - технологически сложная задача ввиду того, что необходимо отслеживать динамику жидкости в постоянно изменяющейся области.
Одна из важных сфер, где востребовано моделирование течений со свободной поверхностью - это компьютерная графика. Задачи компьютерной анимации ставят зачастую противоречивые требования. С одной стороны, технология должна быть эффективной с вычислительной точки зрения, поскольку необходимо проводить множество запусков с разными параметрами за разумное время. С другой стороны, детализация должна быть достаточно высокой, чтобы достичь визуальной реалистичности.
Течение вязкой несжимаемой жидкости описывается уравнениями Навье-Стокса. Возможным компромиссом между сложностью, реалистичным поведением и контролем над потоком является технология, разрабатываемая в последние десятилетия группами специалистов, занимающихся компьютерной графикой [1-4]. В ее основе лежит проекционный метод [5, 6], используемый для приближенного решения системы уравнений Навье-Стокса. Для работы со свободной поверхностью предлагается использовать два взаимодополняющих подхода: метод частиц и метод функции уровня [7].
Предлагаемая технология включает в себя задание начальной трехмерной сцены, моделирование и реалистичное отображение жидкости. Поскольку положение свободной поверхности оказывает большое влияние на динамику жидкости, а также является важным результатом моделирования, предлагается использовать адаптивные сетки, сгущающиеся к поверхности. Использование гексаэдральных разнесенных сеток, построенных по принципу восьмеричного дерева, позволяет динамически перестраивать сетку, сгущая ее к свободной поверхности в каждый момент времени. Кроме того, такие сетки имеют оптимальное число степеней свободы, что также влияет на численную эффективность предлагаемой технологии.
В данной работе последовательно предлагаются базовые уравнения, описывается алгоритм решения и численные методы, используемые для приближенного решения уравнений, приводятся результаты, иллюстрирующие работу алгоритма.
Постановка задачи и базовые уравнения
Рассматривается система уравнений Навье-Стокса, описывающая движение несжимаемой жидкости в области Q . Она состоит из уравнений момента и несжимаемости:
--V Ди + (и-V) + Vp = f,
З/ (1)
V-и = 0,
где V - кинематическая вязкость, / - время, и = (и, V, ') - поле скоростей, р - давление, а f - объемные силы, действующие на жидкость (например, сила тяжести).
Стоит отметить, что граница ЗО области О , которую занимает жидкость, состоит из двух частей: неподвижная граница (стенки сосуда) и перемещающаяся по инерции свободная поверхность. На неподвижной границе поле скоростей удовлетворяет условию Дирихле. Оно может быть как однородным (условие прилипания), так и неоднородным (например, втекающий поток). Для свободной поверхности может быть задано соотношение между давлением жидкости и силами поверхностного натяжения.
Для описания положения свободной границы вводятся функция уровня и частицы. Поверхность
определяется множеством точек х е Ш , где функция уровня ф (х) = 0. Для области, заполненной жидкостью, верно ф (х) < 0, в то время как для воздушной части выполняется ф (х) > 0. Продвижение свободной поверхности описывается уравнением переноса функции уровня:
ф+ ^ф = 0. (2) В ряде приложений требуется инициализировать функцию уровня как расстояние до поверхности со знаком и поддерживать это свойство на протяжении расчета. Помимо определения положения поверхности, функция уровня также несет в себе информацию о ее геометрии: единичная нормаль к поверхности определяется по формуле п = Vф / |Vф|, а локальная кривизна есть к = V - п.
Чтобы отслеживать мелкие элементы, которые не могут быть описаны функцией уровня, используются специальные невесомые частицы, переносимые полем скоростей. Частицы перемещаются с потоком жидкости по закону ёхр / Ж = и(хр) и несут с собой небольшой объем воздуха или жидкости, благодаря чему могут корректировать функцию уровня в случаях, когда это необходимо.
Описание алгоритма
Прежде чем начать расчет, необходимо задать первоначальное состояние задачи: границы расчетной области, начальное положение свободной поверхности, поле скоростей в объеме жидкости и на границе. Задание начального состояния и границ возможно как с помощью аналитических функций, так и с помощью поверхностных триангуляций для объектов и начального объема жидкости.
Для приближенного решения уравнений Навье-Стокса (1) с подвижной границей (2) существует несколько подходов разной степени сложности: от метода дробных шагов [3] до полностью неявных схем [8]. В данной работе используется метод дробных шагов как наименее сложный с вычислительной точки зрения. Подробное описание шагов алгоритма можно найти в [9, 10]. Временной шаг метода состоит из следующих подшагов:
1. обновление поля скоростей:
- решение уравнения моментов (конвективный перенос, действие диффузии и объемных сил);
- проекция на подпространство бездивергентных скоростей;
2. обновление положения свободной поверхности:
- продвижение функции уровня;
- продвижение частиц;
- взаимная корректировка частиц и функции уровня;
- реинициализация (восстановление функции уровня как расстояния до поверхности со знаком);
3. продвижение области;
4. перестроение сетки.
Отметим, что первый подшаг представляет собой вариант проекционного метода для задачи На-вье-Стокса с фиксированной границей. В качестве таковой выступает свободная поверхность, взятая с предыдущего момента времени. Проекционный метод применяется в два этапа. На шаге-предикторе решается уравнение моментов. Конвективный перенос осуществляется с помощью полу-Лагранжева метода, диффузия и объемные силы накладываются явным образом. Для того чтобы сделать полученное поле скоростей и° бездивергентным, вводится коррекция давления, получаемая при решении сеточного уравнения Пуассона, и с ее помощью корректируется поле скоростей:
-V-р) = -—(V-и°), и = и°-Д/р.
Д/
На втором подшаге алгоритма меняется положение свободной границы. Используя новое положение поверхности, определяется часть расчетной области, заполненная жидкостью, а также зоны, в которых должна сгущаться сетка.
Одной из самых важных частей алгоритма является шаг реинициализации. Напомним, что функция уровня отрицательна внутри жидкости и положительна вне нее. Многие приложения требуют, чтобы, помимо этого, функция уровня являлась расстоянием до свободной поверхности со знаком, поэтому имеет смысл стартовать с такой функции и периодически восстанавливать это свойство. Чтобы восстановить функцию уровня как расстояние до поверхности со знаком, сохранив при этом положение свободной границы, необходимо решить уравнение эйконала:
||Уф(x)| = 1, x eQuQffl-r , [ф (x) = 0, x eSQ, где Q.air - область, заполненная воздухом.
Сначала получается начальное приближение в поверхностных ячейках, для чего строится триангуляция свободной поверхности по методу марширующих кубов (Marching Cubes) [11]. После этого решается уравнение эйконала во всей остальной области при помощи быстро идущего метода (Fast Marching Method) [5].
Для некоторых шагов алгоритма, таких как перенос скоростей, функции уровня и частиц, а также перестроения сетки, написана параллельная версия с использованием технологии OpenMP. Для решения системы, которую дает сеточное уравнение Пуассона на шаге проекции, используется пакет PETSc.
Результаты
Описанная технология моделирования использовалась при создании всех графических примеров, приведенных в этой работе. На рис. 1 можно видеть процесс заливания стаканом воды. На рис. 1 изображена модель Армадильо (Броненосца). Модель задана при помощи триангуляции (рис. 1, а), состоящей из 346 тыс. треугольников. В качестве начального значения функции уровня (рис. 1, б) задается расстояние до поверхности модели со знаком. После этого в полученный объем жидкости ударяет струя воды (рис. 1, в, г).
Рис. 1. Модель Армадильо (а), преобразованная в жидкость (б), в которую ударяет струя воды (в, г)
В таблице приведена зависимость числа элементов (N - общее число расчетных ячеек; - число ячеек, заполненных водой) и времени работы подшагов алгоритма от шага сетки (гшг - время построения начального положения и инициализации данных; гШа1 - время одного шага алгоритма; -перенос скоростей; грГО]- - проекция на бездивергентное подпространство скоростей; г- экстрапо-
ляция скоростей с поверхности в воздух; l - перенос функции уровня; tpart - перенос частиц и коррекция функции уровня; tmesh - перестроение сетки; treinit - реинициализация). Видно, что время выполнения одного шага алгоритма пропорционально числу ячеек, заполненных водой, и обратно пропорционально квадрату шага сетки вблизи поверхности. Расчеты проводились на рабочей станции с двумя 4-ядерными процессорами Intel Xeon Х5355@2.66ГГц.
hmin 1/32 1/64 1/128 1/256 1/512
N N ), число элементов 37 409 (1768) 54 377 (7842) 121 101 (36 022) 381305 (159 211) 1 404 439 (684 521)
tinit , с 10,5 13,2 24,2 66,9 237,2
tadv , с 0,18 0,77 3,53 15,74 66,1
t с 1 pro] ■> ^ 0,20 0,60 2,10 8,66 35,9
tfillair , с 0,11 0,35 1,32 5,12 19,9
tls , с 0,03 0,06 0,23 0,98 3,9
t с 1 part ' v 0,07 0,24 1,03 4,37 18,4
t с lreimt ' v 0,12 0,26 1,00 4,43 19,9
tmesh , с 0,10 0,21 0,68 2,42 9,6
ttotal , с 0,81 2,49 9,89 41,72 173,7
Таблица. Зависимость числа элементов и времени работы подшагов алгоритма от шага сетки А,
в г
Рис. 2. Модель Армадильо, заданная в качестве неподвижного препятствия и обтекаемая потоком жидкости: итоговое изображение (а) и срезы расчетной сетки с поверхностью (б, в) и без (г)
В примере на рис. 2 используется модель Армадильо из предыдущего примера для задания неподвижных границ расчетной области. Кубическая область с заданным препятствием заполняется потоком воды. Показано итоговое изображение (рис. 2, а) и разрезы расчетной сетки с поверхностью (рис. 2,б, в) и без (рис. 2, г).
Заключение
В работе предложена эффективная вычислительная технология моделирования трехмерных течений со свободной поверхностью. Благодаря использованию метода дробных шагов и динамически сгущающихся и разгрубляющихся расчетных сеток удается проводить симуляции за приемлемое время (от нескольких часов до суток). Разработанная технология предлагает удобный инструментарий, который может быть полезен в различных приложениях, начиная от компьютерной графики и заканчивая медициной, например, моделирование течений в кровеносных сосудах.
Работа выполнена при финансовой поддержке грантов РФФИ 08-01-00159-a, 09-01-00115-а, программы Президиума РАН 21П «Фундаментальные науки - медицине» и ФЦП «Научные и научно-педагогические кадры инновационной России» 2009-2013 годы.
Автор выражает благодарность К.М. Терехову за помощь в доработке программного кода, а также Ю.В. Василевскому, А.А. Данилову, М.А. Ольшанскому и В.А. Лещинскому за ценные советы и участие в обсуждении спорных моментов и постановки задачи.
Литература
1. Carlson M. Rigid, Melting, and Flowing Fluid, PhD thesis, Georgia Institute of Technology, 2004.
2. Osher S., Fedkiw R. Level Set Methods and Dynamic Implicit Surfaces. - Springer-Verlag, 2002.
3. Osher S. and J.Sethian. Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi equations // Jour. Comp. Phys. - 1988. - V. 79. - Р. 12-49.
4. Sethian J.A. Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and Materials Science, Cambridge University Press, Cambridge, 1999.
5. Chorin A. Numerical solution of the Navier-Stokes equations // Math. Comp. - 1968. - V. 22. - Р. 745-762.
6. Яненко Н.Н. Метод дробных шагов решения многомерных задач математической физики. -Новосибирск: Наука, 1967.
7. Enright D., Fedkiw R., Ferziger J., Mitchell I. A hybrid particle level set method for improved interface capturing // J. Comp. Phys. - 2002. - V.183. - Р. 83-116.
8. Gross S., Reichelt V., Reusken A. A Finite Element Based Level Set Method for Two-Phase Incompressible Flows. Computing and Visualization in Science, 2006.
9. Никитин К.Д. Технология расчета течений со свободной границей с использованием динамических гексаэдральных сеток. Численные методы, параллельные вычисления и информационные технологии // Сборник научных трудов / Под ред. Вл.В. Воеводина и Е.Е. Тыртышникова. - М.: Изд-во МГУ, 2008. - С. 183-198.
10. Nikitin K., Vassilevski Yu. Free surface flow modelling on dynamically refined hexahedral meshes // RJNAMM. - 2008. - V. 23. - Р. 469-485.
11. Lachaud J.-O. Topologically defined iso-surfaces // DGCI. - 1996. - Р. 245-256.
Никитин Кирилл Дмитриевич - Институт вычислительной математики РАН, аспирант, nikitink@dubki.ru УДК 004.622
МЕТОДЫ МАРКИРОВАНИЯ ЦИФРОВЫХ ИЗОБРАЖЕНИЙ В ЧАСТОТНОЙ ОБЛАСТИ А.Ю. Тропченко, Ван Цзянь
Рассматриваются методы маркирования изображений цифровыми водяными знаками, например, для авторизации продукции мультимедиа. Для совмещения реализованных алгоритмов сжатия изображений на основе вейвлет-преобразования с алгоритмами защиты авторских прав наиболее подходит метод вейвлет-маркирования Ли Хуа. Именно в этом методе цифровая подпись помещается в восприимчиво значимую часть преобразованного изображения, которую рассмотренные алгоритмы сжатия стараются сохранить наиболее полно. Ключевые слова: цифровой водяной знак, маркирование, вейвлет, устойчивость.
Введение
Проблема защиты авторского права на мультимедиа-информацию привела к необходимости разработки технологии защиты авторского права и технологии защиты от копирования мультимедиа-информации. Одной из таких технологий является цифровое маркирование данных. При этом необходимо создать алгоритмы, которые позволяли бы подписывать или маркировать мультимедиа-данные без потери информативности и качества воспроизведения, однако при этом в любой момент можно было бы определить, кому принадлежат авторские права на тот или иной мультимедиа-продукт. Данные, скрыто встроенные в мультимедиа-продукты, называют цифровым водяным знаком - ЦВЗ [1, 2].
