РЕАЛИСТИЧНОЕ МОДЕЛИРОВАНИЕ СВОБОДНОЙ ВОДНОЙ ПОВЕРХНОСТИ НА АДАПТИВНЫХ СЕТКАХ ТИПА «ВОСЬМЕРИЧНОЕ ДЕРЕВО»
РЕАЛИСТИЧНОЕ МОДЕЛИРОВАНИЕ СВОБОДНОЙ ВОДНОЙ ПОВЕРХНОСТИ ...
УДК 517.958:532.5
РЕАЛИСТИЧНОЕ МОДЕЛИРОВАНИЕ СВОБОДНОЙ ВОДНОЙ ПОВЕРХНОСТИ НА АДАПТИВНЫХ СЕТКАХ ТИПА «ВОСЬМЕРИЧНОЕ ДЕРЕВО»
К.Д. Никитин Предлагается эффективная технология моделирования течения вязкой несжимаемой жидкости со свободной границей. В технологии объединены проекционный метод для решения уравнений Навье–Стокса и метод функции уровня и частиц для работы со свободной поверхностью. Расчет проводится на адаптивно сгущающихся гексаэдральных сетках типа «восьмеричное дерево». Ключевые слова: вычислительная технология, течения со свободной поверхностью, компьютерная графика.
Введение Моделирование течений несжимаемой жидкости привлекает большой интерес в научном сообществе. Рассматриваемые задачи включают моделирование морских волн, заливание и обтекание объектов, падение капель, образование брызг и многое другое. Моделирование подобных явлений – технологически сложная задача ввиду того, что необходимо отслеживать динамику жидкости в постоянно изменяющейся области. Одна из важных сфер, где востребовано моделирование течений со свободной поверхностью – это компьютерная графика. Задачи компьютерной анимации ставят зачастую противоречивые требования. С одной стороны, технология должна быть эффективной с вычислительной точки зрения, поскольку необходимо проводить множество запусков с разными параметрами за разумное время. С другой стороны, детализация должна быть достаточно высокой, чтобы достичь визуальной реалистичности. Течение вязкой несжимаемой жидкости описывается уравнениями Навье–Стокса. Возможным компромиссом между сложностью, реалистичным поведением и контролем над потоком является технология, разрабатываемая в последние десятилетия группами специалистов, занимающихся компьютерной графикой [1–4]. В ее основе лежит проекционный метод [5, 6], используемый для приближенного решения системы уравнений Навье–Стокса. Для работы со свободной поверхностью предлагается использовать два взаимодополняющих подхода: метод частиц и метод функции уровня [7]. Предлагаемая технология включает в себя задание начальной трехмерной сцены, моделирование и реалистичное отображение жидкости. Поскольку положение свободной поверхности оказывает большое влияние на динамику жидкости, а также является важным результатом моделирования, предлагается использовать адаптивные сетки, сгущающиеся к поверхности. Использование гексаэдральных разнесенных сеток, построенных по принципу восьмеричного дерева, позволяет динамически перестраивать сетку, сгущая ее к свободной поверхности в каждый момент времени. Кроме того, такие сетки имеют оптимальное число степеней свободы, что также влияет на численную эффективность предлагаемой технологии. В данной работе последовательно предлагаются базовые уравнения, описывается алгоритм решения и численные методы, используемые для приближенного решения уравнений, приводятся результаты, иллюстрирующие работу алгоритма.
Постановка задачи и базовые уравнения Рассматривается система уравнений Навье–Стокса, описывающая движение несжимаемой жидкости в области . Она состоит из уравнений момента и несжимаемости:
60 Научно-технический вестник Санкт-Петербургского государственного университета
информационных технологий, механики и оптики, 2010, № 6 (70)
К.Д. Никитин
u t
u
+
u
u
+ p
=
f
,
(1)
u = 0,
где – кинематическая вязкость, t – время, u = (u, ν, w) – поле скоростей, p – давление, а f – объемные
силы, действующие на жидкость (например, сила тяжести). Стоит отметить, что граница области , которую занимает жидкость, состоит из двух частей:
неподвижная граница (стенки сосуда) и перемещающаяся по инерции свободная поверхность. На неподвижной границе поле скоростей удовлетворяет условию Дирихле. Оно может быть как однородным (условие прилипания), так и неоднородным (например, втекающий поток). Для свободной поверхности может быть задано соотношение между давлением жидкости и силами поверхностного натяжения.
Для описания положения свободной границы вводятся функция уровня и частицы. Поверхность
определяется множеством точек x 3 , где функция уровня (x) = 0 . Для области, заполненной жид-
костью, верно (x) < 0 , в то время как для воздушной части выполняется (x) > 0 . Продвижение сво-
бодной поверхности описывается уравнением переноса функции уровня:
t +u = 0 .
(2)
В ряде приложений требуется инициализировать функцию уровня как расстояние до поверхности
со знаком и поддерживать это свойство на протяжении расчета. Помимо определения положения по-
верхности, функция уровня также несет в себе информацию о ее геометрии: единичная нормаль к по-
верхности определяется по формуле n = / , а локальная кривизна есть k = n .
Чтобы отслеживать мелкие элементы, которые не могут быть описаны функцией уровня, используются специальные невесомые частицы, переносимые полем скоростей. Частицы перемещаются с пото-
ком жидкости по закону dx p / dt = u x p и несут с собой небольшой объем воздуха или жидкости, благо-
даря чему могут корректировать функцию уровня в случаях, когда это необходимо.
Описание алгоритма
Прежде чем начать расчет, необходимо задать первоначальное состояние задачи: границы расчет-
ной области, начальное положение свободной поверхности, поле скоростей в объеме жидкости и на гра-
нице. Задание начального состояния и границ возможно как с помощью аналитических функций, так и с
помощью поверхностных триангуляций для объектов и начального объема жидкости.
Для приближенного решения уравнений Навье–Стокса (1) с подвижной границей (2) существует
несколько подходов разной степени сложности: от метода дробных шагов [3] до полностью неявных
схем [8]. В данной работе используется метод дробных шагов как наименее сложный с вычислительной
точки зрения. Подробное описание шагов алгоритма можно найти в [9, 10]. Временной шаг метода со-
стоит из следующих подшагов:
1. обновление поля скоростей:
решение уравнения моментов (конвективный перенос, действие диффузии и объемных сил);
проекция на подпространство бездивергентных скоростей; 2. обновление положения свободной поверхности:
продвижение функции уровня;
продвижение частиц;
взаимная корректировка частиц и функции уровня;
реинициализация (восстановление функции уровня как расстояния до поверхности со знаком); 3. продвижение области;
4. перестроение сетки.
Отметим, что первый подшаг представляет собой вариант проекционного метода для задачи На-
вье–Стокса с фиксированной границей. В качестве таковой выступает свободная поверхность, взятая с
предыдущего момента времени. Проекционный метод применяется в два этапа. На шаге-предикторе решается уравнение моментов. Конвективный перенос осуществляется с помощью полу-Лагранжева мето-
да, диффузия и объемные силы накладываются явным образом. Для того чтобы сделать полученное поле скоростей u бездивергентным, вводится коррекция давления, получаемая при решении сеточного урав-
нения Пуассона, и с ее помощью корректируется поле скоростей:
(p)
1 t
(
u)
,
u u t p .
На втором подшаге алгоритма меняется положение свободной границы. Используя новое положе-
ние поверхности, определяется часть расчетной области, заполненная жидкостью, а также зоны, в кото-
рых должна сгущаться сетка.
Научно-технический вестник Санкт-Петербургского государственного университета информационных технологий, механики и оптики, 2010, № 6 (70)
61
РЕАЛИСТИЧНОЕ МОДЕЛИРОВАНИЕ СВОБОДНОЙ ВОДНОЙ ПОВЕРХНОСТИ ...
Одной из самых важных частей алгоритма является шаг реинициализации. Напомним, что функция уровня отрицательна внутри жидкости и положительна вне нее. Многие приложения требуют, чтобы, помимо этого, функция уровня являлась расстоянием до свободной поверхности со знаком, поэтому имеет смысл стартовать с такой функции и периодически восстанавливать это свойство. Чтобы восстановить функцию уровня как расстояние до поверхности со знаком, сохранив при этом положение свободной границы, необходимо решить уравнение эйконала:
(x) 1, (x) 0,
x air , x,
где air – область, заполненная воздухом.
Сначала получается начальное приближение в поверхностных ячейках, для чего строится триангуляция свободной поверхности по методу марширующих кубов (Marching Cubes) [11]. После этого решается уравнение эйконала во всей остальной области при помощи быстро идущего метода (Fast Marching Method) [5].
Для некоторых шагов алгоритма, таких как перенос скоростей, функции уровня и частиц, а также перестроения сетки, написана параллельная версия с использованием технологии OpenMP. Для решения системы, которую дает сеточное уравнение Пуассона на шаге проекции, используется пакет PETSc.
Результаты
Описанная технология моделирования использовалась при создании всех графических примеров, приведенных в этой работе. На рис. 1 можно видеть процесс заливания стаканом воды. На рис. 1 изображена модель Армадильо (Броненосца). Модель задана при помощи триангуляции (рис. 1, а), состоящей из 346 тыс. треугольников. В качестве начального значения функции уровня (рис. 1, б) задается расстояние до поверхности модели со знаком. После этого в полученный объем жидкости ударяет струя воды (рис. 1, в, г).
аб
вг Рис. 1. Модель Армадильо (а), преобразованная в жидкость (б), в которую ударяет струя воды (в, г)
В таблице приведена зависимость числа элементов ( N – общее число расчетных ячеек; Nw – число ячеек, заполненных водой) и времени работы подшагов алгоритма от шага сетки ( tinit – время построения начального положения и инициализации данных; ttotal – время одного шага алгоритма; tadv – перенос скоростей; t proj – проекция на бездивергентное подпространство скоростей; t fillair – экстрапо-
62 Научно-технический вестник Санкт-Петербургского государственного университета
информационных технологий, механики и оптики, 2010, № 6 (70)
К.Д. Никитин
ляция скоростей с поверхности в воздух; tls – перенос функции уровня; t part – перенос частиц и коррек-
ция функции уровня; tmesh – перестроение сетки; treinit – реинициализация). Видно, что время выполнения одного шага алгоритма пропорционально числу ячеек, заполненных водой, и обратно пропорционально квадрату шага сетки вблизи поверхности. Расчеты проводились на рабочей станции с двумя 4ядерными процессорами Intel Xeon X5355@2.66ГГц.
hmin
N Nw ,
число элементов tinit , с
tadv , с t proj , с
t fillair , с
tls , с t part , с
treinit , с tmesh , с
ttotal , с
1/32
37 409 (1768) 10,5
0,18 0,20
0,11
0,03 0,07
0,12 0,10
0,81
1/64
54 377 (7842) 13,2
0,77 0,60
0,35
0,06 0,24
0,26 0,21
2,49
1/128
121 101 (36 022)
24,2
3,53 2,10
1,32
0,23 1,03
1,00 0,68
9,89
1/256
381 305 (159 211)
66,9
15,74 8,66
5,12
0,98 4,37
4,43 2,42
41,72
1/512
1 404 439 (684 521)
237,2
66,1 35,9
19,9
3,9 18,4
19,9 9,6
173,7
Таблица. Зависимость числа элементов и времени работы подшагов алгоритма от шага сетки hmin
аб
вг
Рис. 2. Модель Армадильо, заданная в качестве неподвижного препятствия и обтекаемая потоком жидкости: итоговое изображение (а) и срезы расчетной сетки с поверхностью (б, в) и без (г)
В примере на рис. 2 используется модель Армадильо из предыдущего примера для задания неподвижных границ расчетной области. Кубическая область с заданным препятствием заполняется потоком воды. Показано итоговое изображение (рис. 2, а) и разрезы расчетной сетки с поверхностью (рис. 2,б, в) и без (рис. 2, г).
Научно-технический вестник Санкт-Петербургского государственного университета информационных технологий, механики и оптики, 2010, № 6 (70)
63
РЕАЛИСТИЧНОЕ МОДЕЛИРОВАНИЕ СВОБОДНОЙ ВОДНОЙ ПОВЕРХНОСТИ ...
Заключение
В работе предложена эффективная вычислительная технология моделирования трехмерных течений со свободной поверхностью. Благодаря использованию метода дробных шагов и динамически сгущающихся и разгрубляющихся расчетных сеток удается проводить симуляции за приемлемое время (от нескольких часов до суток). Разработанная технология предлагает удобный инструментарий, который может быть полезен в различных приложениях, начиная от компьютерной графики и заканчивая медициной, например, моделирование течений в кровеносных сосудах.
Работа выполнена при финансовой поддержке грантов РФФИ 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
64 Научно-технический вестник Санкт-Петербургского государственного университета
информационных технологий, механики и оптики, 2010, № 6 (70)
УДК 517.958:532.5
РЕАЛИСТИЧНОЕ МОДЕЛИРОВАНИЕ СВОБОДНОЙ ВОДНОЙ ПОВЕРХНОСТИ НА АДАПТИВНЫХ СЕТКАХ ТИПА «ВОСЬМЕРИЧНОЕ ДЕРЕВО»
К.Д. Никитин Предлагается эффективная технология моделирования течения вязкой несжимаемой жидкости со свободной границей. В технологии объединены проекционный метод для решения уравнений Навье–Стокса и метод функции уровня и частиц для работы со свободной поверхностью. Расчет проводится на адаптивно сгущающихся гексаэдральных сетках типа «восьмеричное дерево». Ключевые слова: вычислительная технология, течения со свободной поверхностью, компьютерная графика.
Введение Моделирование течений несжимаемой жидкости привлекает большой интерес в научном сообществе. Рассматриваемые задачи включают моделирование морских волн, заливание и обтекание объектов, падение капель, образование брызг и многое другое. Моделирование подобных явлений – технологически сложная задача ввиду того, что необходимо отслеживать динамику жидкости в постоянно изменяющейся области. Одна из важных сфер, где востребовано моделирование течений со свободной поверхностью – это компьютерная графика. Задачи компьютерной анимации ставят зачастую противоречивые требования. С одной стороны, технология должна быть эффективной с вычислительной точки зрения, поскольку необходимо проводить множество запусков с разными параметрами за разумное время. С другой стороны, детализация должна быть достаточно высокой, чтобы достичь визуальной реалистичности. Течение вязкой несжимаемой жидкости описывается уравнениями Навье–Стокса. Возможным компромиссом между сложностью, реалистичным поведением и контролем над потоком является технология, разрабатываемая в последние десятилетия группами специалистов, занимающихся компьютерной графикой [1–4]. В ее основе лежит проекционный метод [5, 6], используемый для приближенного решения системы уравнений Навье–Стокса. Для работы со свободной поверхностью предлагается использовать два взаимодополняющих подхода: метод частиц и метод функции уровня [7]. Предлагаемая технология включает в себя задание начальной трехмерной сцены, моделирование и реалистичное отображение жидкости. Поскольку положение свободной поверхности оказывает большое влияние на динамику жидкости, а также является важным результатом моделирования, предлагается использовать адаптивные сетки, сгущающиеся к поверхности. Использование гексаэдральных разнесенных сеток, построенных по принципу восьмеричного дерева, позволяет динамически перестраивать сетку, сгущая ее к свободной поверхности в каждый момент времени. Кроме того, такие сетки имеют оптимальное число степеней свободы, что также влияет на численную эффективность предлагаемой технологии. В данной работе последовательно предлагаются базовые уравнения, описывается алгоритм решения и численные методы, используемые для приближенного решения уравнений, приводятся результаты, иллюстрирующие работу алгоритма.
Постановка задачи и базовые уравнения Рассматривается система уравнений Навье–Стокса, описывающая движение несжимаемой жидкости в области . Она состоит из уравнений момента и несжимаемости:
60 Научно-технический вестник Санкт-Петербургского государственного университета
информационных технологий, механики и оптики, 2010, № 6 (70)
К.Д. Никитин
u t
u
+
u
u
+ p
=
f
,
(1)
u = 0,
где – кинематическая вязкость, t – время, u = (u, ν, w) – поле скоростей, p – давление, а f – объемные
силы, действующие на жидкость (например, сила тяжести). Стоит отметить, что граница области , которую занимает жидкость, состоит из двух частей:
неподвижная граница (стенки сосуда) и перемещающаяся по инерции свободная поверхность. На неподвижной границе поле скоростей удовлетворяет условию Дирихле. Оно может быть как однородным (условие прилипания), так и неоднородным (например, втекающий поток). Для свободной поверхности может быть задано соотношение между давлением жидкости и силами поверхностного натяжения.
Для описания положения свободной границы вводятся функция уровня и частицы. Поверхность
определяется множеством точек x 3 , где функция уровня (x) = 0 . Для области, заполненной жид-
костью, верно (x) < 0 , в то время как для воздушной части выполняется (x) > 0 . Продвижение сво-
бодной поверхности описывается уравнением переноса функции уровня:
t +u = 0 .
(2)
В ряде приложений требуется инициализировать функцию уровня как расстояние до поверхности
со знаком и поддерживать это свойство на протяжении расчета. Помимо определения положения по-
верхности, функция уровня также несет в себе информацию о ее геометрии: единичная нормаль к по-
верхности определяется по формуле n = / , а локальная кривизна есть k = n .
Чтобы отслеживать мелкие элементы, которые не могут быть описаны функцией уровня, используются специальные невесомые частицы, переносимые полем скоростей. Частицы перемещаются с пото-
ком жидкости по закону dx p / dt = u x p и несут с собой небольшой объем воздуха или жидкости, благо-
даря чему могут корректировать функцию уровня в случаях, когда это необходимо.
Описание алгоритма
Прежде чем начать расчет, необходимо задать первоначальное состояние задачи: границы расчет-
ной области, начальное положение свободной поверхности, поле скоростей в объеме жидкости и на гра-
нице. Задание начального состояния и границ возможно как с помощью аналитических функций, так и с
помощью поверхностных триангуляций для объектов и начального объема жидкости.
Для приближенного решения уравнений Навье–Стокса (1) с подвижной границей (2) существует
несколько подходов разной степени сложности: от метода дробных шагов [3] до полностью неявных
схем [8]. В данной работе используется метод дробных шагов как наименее сложный с вычислительной
точки зрения. Подробное описание шагов алгоритма можно найти в [9, 10]. Временной шаг метода со-
стоит из следующих подшагов:
1. обновление поля скоростей:
решение уравнения моментов (конвективный перенос, действие диффузии и объемных сил);
проекция на подпространство бездивергентных скоростей; 2. обновление положения свободной поверхности:
продвижение функции уровня;
продвижение частиц;
взаимная корректировка частиц и функции уровня;
реинициализация (восстановление функции уровня как расстояния до поверхности со знаком); 3. продвижение области;
4. перестроение сетки.
Отметим, что первый подшаг представляет собой вариант проекционного метода для задачи На-
вье–Стокса с фиксированной границей. В качестве таковой выступает свободная поверхность, взятая с
предыдущего момента времени. Проекционный метод применяется в два этапа. На шаге-предикторе решается уравнение моментов. Конвективный перенос осуществляется с помощью полу-Лагранжева мето-
да, диффузия и объемные силы накладываются явным образом. Для того чтобы сделать полученное поле скоростей u бездивергентным, вводится коррекция давления, получаемая при решении сеточного урав-
нения Пуассона, и с ее помощью корректируется поле скоростей:
(p)
1 t
(
u)
,
u u t p .
На втором подшаге алгоритма меняется положение свободной границы. Используя новое положе-
ние поверхности, определяется часть расчетной области, заполненная жидкостью, а также зоны, в кото-
рых должна сгущаться сетка.
Научно-технический вестник Санкт-Петербургского государственного университета информационных технологий, механики и оптики, 2010, № 6 (70)
61
РЕАЛИСТИЧНОЕ МОДЕЛИРОВАНИЕ СВОБОДНОЙ ВОДНОЙ ПОВЕРХНОСТИ ...
Одной из самых важных частей алгоритма является шаг реинициализации. Напомним, что функция уровня отрицательна внутри жидкости и положительна вне нее. Многие приложения требуют, чтобы, помимо этого, функция уровня являлась расстоянием до свободной поверхности со знаком, поэтому имеет смысл стартовать с такой функции и периодически восстанавливать это свойство. Чтобы восстановить функцию уровня как расстояние до поверхности со знаком, сохранив при этом положение свободной границы, необходимо решить уравнение эйконала:
(x) 1, (x) 0,
x air , x,
где air – область, заполненная воздухом.
Сначала получается начальное приближение в поверхностных ячейках, для чего строится триангуляция свободной поверхности по методу марширующих кубов (Marching Cubes) [11]. После этого решается уравнение эйконала во всей остальной области при помощи быстро идущего метода (Fast Marching Method) [5].
Для некоторых шагов алгоритма, таких как перенос скоростей, функции уровня и частиц, а также перестроения сетки, написана параллельная версия с использованием технологии OpenMP. Для решения системы, которую дает сеточное уравнение Пуассона на шаге проекции, используется пакет PETSc.
Результаты
Описанная технология моделирования использовалась при создании всех графических примеров, приведенных в этой работе. На рис. 1 можно видеть процесс заливания стаканом воды. На рис. 1 изображена модель Армадильо (Броненосца). Модель задана при помощи триангуляции (рис. 1, а), состоящей из 346 тыс. треугольников. В качестве начального значения функции уровня (рис. 1, б) задается расстояние до поверхности модели со знаком. После этого в полученный объем жидкости ударяет струя воды (рис. 1, в, г).
аб
вг Рис. 1. Модель Армадильо (а), преобразованная в жидкость (б), в которую ударяет струя воды (в, г)
В таблице приведена зависимость числа элементов ( N – общее число расчетных ячеек; Nw – число ячеек, заполненных водой) и времени работы подшагов алгоритма от шага сетки ( tinit – время построения начального положения и инициализации данных; ttotal – время одного шага алгоритма; tadv – перенос скоростей; t proj – проекция на бездивергентное подпространство скоростей; t fillair – экстрапо-
62 Научно-технический вестник Санкт-Петербургского государственного университета
информационных технологий, механики и оптики, 2010, № 6 (70)
К.Д. Никитин
ляция скоростей с поверхности в воздух; tls – перенос функции уровня; t part – перенос частиц и коррек-
ция функции уровня; tmesh – перестроение сетки; treinit – реинициализация). Видно, что время выполнения одного шага алгоритма пропорционально числу ячеек, заполненных водой, и обратно пропорционально квадрату шага сетки вблизи поверхности. Расчеты проводились на рабочей станции с двумя 4ядерными процессорами Intel Xeon X5355@2.66ГГц.
hmin
N Nw ,
число элементов tinit , с
tadv , с t proj , с
t fillair , с
tls , с t part , с
treinit , с tmesh , с
ttotal , с
1/32
37 409 (1768) 10,5
0,18 0,20
0,11
0,03 0,07
0,12 0,10
0,81
1/64
54 377 (7842) 13,2
0,77 0,60
0,35
0,06 0,24
0,26 0,21
2,49
1/128
121 101 (36 022)
24,2
3,53 2,10
1,32
0,23 1,03
1,00 0,68
9,89
1/256
381 305 (159 211)
66,9
15,74 8,66
5,12
0,98 4,37
4,43 2,42
41,72
1/512
1 404 439 (684 521)
237,2
66,1 35,9
19,9
3,9 18,4
19,9 9,6
173,7
Таблица. Зависимость числа элементов и времени работы подшагов алгоритма от шага сетки hmin
аб
вг
Рис. 2. Модель Армадильо, заданная в качестве неподвижного препятствия и обтекаемая потоком жидкости: итоговое изображение (а) и срезы расчетной сетки с поверхностью (б, в) и без (г)
В примере на рис. 2 используется модель Армадильо из предыдущего примера для задания неподвижных границ расчетной области. Кубическая область с заданным препятствием заполняется потоком воды. Показано итоговое изображение (рис. 2, а) и разрезы расчетной сетки с поверхностью (рис. 2,б, в) и без (рис. 2, г).
Научно-технический вестник Санкт-Петербургского государственного университета информационных технологий, механики и оптики, 2010, № 6 (70)
63
РЕАЛИСТИЧНОЕ МОДЕЛИРОВАНИЕ СВОБОДНОЙ ВОДНОЙ ПОВЕРХНОСТИ ...
Заключение
В работе предложена эффективная вычислительная технология моделирования трехмерных течений со свободной поверхностью. Благодаря использованию метода дробных шагов и динамически сгущающихся и разгрубляющихся расчетных сеток удается проводить симуляции за приемлемое время (от нескольких часов до суток). Разработанная технология предлагает удобный инструментарий, который может быть полезен в различных приложениях, начиная от компьютерной графики и заканчивая медициной, например, моделирование течений в кровеносных сосудах.
Работа выполнена при финансовой поддержке грантов РФФИ 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
64 Научно-технический вестник Санкт-Петербургского государственного университета
информационных технологий, механики и оптики, 2010, № 6 (70)