CC BY
70
УДК 532.5: 519.652
НАПРАВЛЕНИЯ НАУЧНЫХ ИССЛЕДОВАНИЙ КАФЕДРЫ ЮНЕСКО ПО НОВЫМ ИНФОРМАЦИОННЫМ ТЕХНОЛОГИЯМ. ЧАСТЬ 2. БЕССЕТОЧНЫЕ МЕТОДЫ
К. Е. Афанасьев, С. Н. Карабцев, Р. С. Макарчук, Т. С. Рейн
RESEARCH AREAS OF THE UNESCO CHAIR ON NEW INFORMATION TECHNOLOGIES.
PART 2. MESHLESS METHODS K. E. Afanasyev, S. N. Karabtsev, R. S. Makarchuk, T. S. Rein
В работе приводится обзор основных результатов решения нелинейных задач гидродинамики со свободными границами бессеточными методами. Особое внимание уделяется расчетам, полученным методами естественных соседей и методом сглаженных частиц. Описываются основные результаты, полученные в научной школе за период с 2004 года.
This paper provides an overview of the main results of solving the nonlinear free-boundary hydrodynamic problems with meshless methods. Special attention is paid to computations using the Natural Element Methods and Smoothed Particle Hydrodynamics. The main results obtained by the scientific school since 2004 are described.
Ключевые слова: бессеточные методы, гидродинамические нагрузки, свободные границы, несжимаемая жидкость, метод сглаженных частиц, метод естественных соседей, высокопроизводительные вычисления.
Keywords: SPH, NEM, GNEM, incompressible fluid flows, hydrodynamic loads, meshless methods, free surface, high-performance computation.
Введение
Данная статья является продолжением работы «Направления научных исследований кафедры ЮНЕСКО по новым информационным технологиям», опубликованной по результатам исследований, выполненных до 2004 года [2]. Этот период характеризуется разработкой численных методов решения задач динамики жидкости со свободными границами методами конечных элементов (МКЭ) [47; 65], граничных элементов (МГЭ) [43; 47] и его разновидностями -комплексный метод граничных элементов (КМГЭ) [34] для решения задач в пространственной постановке [4; 46].
Существенным недостатком МКЭ, МГЭ и КМГЭ является невозможность дальнейшего продолжения расчетов после изменения связности расчетной области. Например, обрушение волны можно рассчитать только до момента соприкосновения гребня волны с ее подошвой, далее проводить расчет становится невозможным в силу изменения связности расчетной области и перехлеста границ, а также возможного взаимного перемешивания слоев жидкости.
Состояние исследований в области нелинейной гидродинамики [29], в частности, в задачах со свободными границами, где происходит взаимодействие потока жидкости с твердыми телами и прибрежными сооружениями, показывает, что новые методы численного исследования и численные алгоритмы часто приходится разрабатывать для конкретных условий и особенностей задачи. Например, наиболее интересными и актуальными на сегодняшний день являются задачи исследования поверхностных волн, порожденных различными факторами, и их взаимодействия с твердыми и упругими поверхностями [41; 42; 43; 49; 55; 56; 62].
Для комплексного исследования таких нелинейных процессов потребовалась разработка численных методов, учитывающих возникающие сложности расчета. В связи с этим были начаты работы по исследо-
ванию и развитию условно-бессеточных и полностью бессеточных методов для решения двух- и трехмерных задач динамики несжимаемой жидкости со свободными границами, в которых происходят существенные деформации расчетной области, сопровождаемые разрывами основных характеристик течения. К таким задачам относятся: задачи о взаимодействии твердых тел с жидкостью, выходе волн на мелководье, разрушении плотины, взаимодействии жидкости с днищем судов, перемешивании жидкостей разной плотности, падении капли на твердые и жидкие поверхности и др. Для решения указанных сложных задач используются условно-бессеточные численные методы: методы естественных соседей - Natural Element Method (NEM) [30; 31; 32; 50; 51; 52] и General NEM (GNEM) [6; 7; 23; 37; 38; 40; 62], полностью бессеточные: методы сглаженных частиц -Smooth Particle Hydrodynamics (SPH) [5; 8; 21; 22; 23; 25; 26] и Incompressible SPH (ISPH) [5; 8; 9; 10; 19; 20]. Обзор других бессеточных методов можно найти, например, в [17; 20; 66].
В данной статье приводятся результаты решения задач с использованием методов естественных соседей и метода сглаженных частиц, делаются ссылки на оригинальные работы.
Основные направления исследований
Работы коллектива научной школы профессора К. Е. Афанасьева направлены на комплексное исследование задач о движении линейных и нелинейных волн на поверхности жидкости в бассейне со сложной геометрической конфигурацией дна, дальнейшем их распространении при выходе на мелководье и о нелинейном взаимодействии жидкости с погруженными и плавающими телами.
В обзорной статье [2] описывались результаты работ, выполненных до 2004 года. В последующий период методом граничных элементов был решен ряд
новых задач, нашедших отражение в коллективной монографии [1].
Ниже кратко описываются исследования, результаты которых приведены в [45].
Взаимодействие поверхностных волн с препятствием
Приведено описание решения ряда тестовых и практических задач гидродинамики идеальной однородной несжимаемой жидкости с поверхностными гравитационными волнами в областях со сложной формой границ.
Представлены результаты численного моделирования процессов взаимодействия поверхностных волн с частично погруженным в жидкость телом прямоугольного сечения. Постановка задачи была взята из монографии [69], где авторами проводится исследование влияния варьируемых параметров на величину заплесков, амплитуды отраженной и прошедшей волн с помощью конечно-разностных методов расчета на адаптивных сетках. Задача решается в полной нелинейной постановке методом граничных элементов. Исследуется влияние ширины и заглубления тела, а также амплитуды волны на волновые режимы.
Численное моделирование взаимодействия со-литона с подводным препятствием
Представлены результаты расчетов волнового движения идеальной однородной несжимаемой жидкости при взаимодействии уединенной волны с телом прямоугольного сечения, расположенным на дне. Исследуется влияние ширины и высоты тела, а также амплитуды волны на основные характеристики возникающего течения жидкости. Проводится анализ кинематических и динамических характеристик волны при ее движении над подводным препятствием, расположенным на горизонтальном дне. Изучается изменение динамической нагрузки при взаимодействии поверхностных волн с твердыми границами. Постановка задачи была взята из работы [68], где данная задача исследуется с помощью конечно-разностных методов расчета на адаптивных сетках.
Численное моделирование генерации поверхностных волн движением оползня
Изучаются основные результаты исследования генерации поверхностных волн движением оползня. Для изучения закономерностей волнообразования исследуется зависимость характеристик процесса от основных параметров задачи: длины и ширины оползня, глубины его заглубления и закона его движения.
Генерация волн цунами движением оползня с помощью различных математических моделей является актуальной задачей. Следует отметить работы по исследованию различных моделей движения оползней и работы по совместному лабораторному и вычислительному моделированию рассматриваемых волновых процессов. Эти исследования показали, что самые общие характеристики изучаемых волновых режимов могут быть определены с помощью простейших мо-
делей теории мелкой воды. Это касается волн, распространяющихся в сторону берега, противоположную направлению движения оползня. Однако волны, распространяющиеся в открытую зону в попутном оползню направлении, требуют для своего воспроизведения учета вертикальных процессов. Этот факт подтверждается близостью результатов, полученных с помощью полной гидродинамической модели и экспериментальных данных.
Несогласованная линейная аппроксимация в методе граничных элементов для решения пространственных задач
Продолжались исследования по трехмерному моделированию задач динамики газовых пузырей [4]. Вместе с тем большое внимание было уделено несогласованной изопараметрической линейной аппроксимации при изучении трехмерных задач динамики парогазового пузыря. Данная аппроксимация применялась в методе конечных элементов (МКЭ), её преимущества по сравнению с классической линейной аппроксимацией известны и описаны для уравнений Навье-Стокса [64]. Тем не менее, работ, в которых бы данная аппроксимация применялась к трехмерным задачам метода граничных элементов, нет. Исследованию этого вопроса посвящена работа [46], в которой основное внимание уделяется сравнению несогласованной линейной аппроксимации с классической линейной аппроксимацией. В работе рассмотрена методика вычисления регулярных и сингулярных интегралов, участвующих в формировании коэффициентов результирующей СЛАУ. Большое внимание уделяется методике вычисления скоростей. Показано преимущество несогласованной аппроксимации перед согласованной для метода граничных элементов.
ПОСТАНОВКА ЗАДАЧИ.
НЕСЖИМАЕМАЯ ЖИДКОСТЬ
Период исследований с 2004 года характеризуется использованием современных численных методов: NEM, GNEM, SPH, ISPH. Существенная особенность этих методов - возможность решения задач динамики вязкой и идеальной несжимаемой жидкости со свободными границами с сильными деформациями расчетной области.
Появление аппроксимаций, базирующихся на новых способах представления функций, позволило как усовершенствовать классические методы математического моделирования, так и построить ряд новых. Разработка новых бессеточных и условно-бессеточных методов позволила пересмотреть список актуальных научно-технических задач, доступных для моделирования. Особая эффективность применения этих методов достигается в связи с новой парадигмой программирования - параллельными вычислениями [49; 51; 62].
Выпишем полную нелинейную постановку задачи
о движении несжимаемой жидкости со свободной границей [6; 52].
Рис. 1. Схема области течения
Пусть в области D (рис. 1), ограниченной свободной поверхностью Г0 и твердыми границами Г!, задано течение несжимаемой жидкости, описываемое уравнениями движения Эйлера (1E) либо Навье-Стокса (1NS). Помимо уравнений движения должно выполняться условие несжимаемости (2E), (2NS). Также необходимо задать кинематические (3E), (3NS) и граничные условия (4E), (4NS) (таблица 1).
Таблица 1
Общая постановка нелинейных задач гидродинамики несжимаемой жидкости [48, 67]
Du = — — Vp + f , x Є D (1E) Dt p — =— — Vp + ^ Ди + f (1NS) Dt p p
V ■ u = 0, x Є D (2E) V- u = 0, x Є D (2NS)
dx / dt = u, x Є D (3E) dx / dt = u, x Є D (3NS)
{ P (x ’!^ ^7 = Patm0 (4E) [ P ( x , t ) ■ n |x є Г — = 0 { Pn (x, t)Ur0 = T(x,t) ■ n — patmn (4NS) [u(x, t) Ur— = 0 ( )
В нестационарной задаче требуется задание начальных условий для расчетных узлов х(0) = х0 и неизвестных функций во всей области течения: и(х,0) = и0(х), р(х,0) = р0(х).
Здесь х = (хь х2) - радиус-вектор расчетных узлов, и = (%, щ) - вектор скорости, р - давление, т(х, Ь) - тензор вязких напряжений, рп (х, Ь) - вектор напряжений на свободной поверхности в направлении внешней нормали, р - плотность, { = (/1, /2) =
(0, — g) - вектор плотности внешних сил, п = (%, п2) - внешняя нормаль к границе области.
МЕТОДЫ КЕМ и ОКЕМ
В основе методов естественных соседей лежит идея проекционно-сеточных методов [57] при специальном выборе базисных функций, определенных на диаграмме Вороного - это функции Сибсона [28; 30] и Лапласа [31, 44]. Базисные функции в проекционносеточных методах должны удовлетворять правилу разложения единицы. Это гарантирует точное удовлетворение граничных условий первого рода, а условия второго рода входят в правую часть результирующих уравнений естественным образом.
Алгоритм движения по времени
Для реализации численного решения пространственно-временной задачи используется идея расщепления определяющих уравнений по физическим фак-
торам, предложенная в работе [59, 70], а для задач вязкой жидкости использованная в [11].
Устойчивость решения системы уравнений Навье-Стокса методами, основанными на методе Галеркина, обеспечивается выбором конечно-элементных пространств для скорости и давления: степени интерполяционных полиномов для компонент вектора скорости и давления должны удовлетворять условию Ла-дыженской-Бабушки-Бреззи (ЛББ). В данной работе для аппроксимации уравнения неразрывности используются линейные базисные функции (функции формы расширенной интерполяции Лапласа), для уравнения движения применяются квадратичные базисные функции (функции формы Сибсона). Построение так называемого обобщенного метода естественных соседей (GNEM - General Natural Element Method), гибрида методов MFEM и NEM, приводит к удовлетворению условий ЛББ для совместной аппроксимации, что гарантирует невырожденность решения.
Суть метода расщепления заключается в разбиении описываемого физического процесса на два: конвекцию-диффузию и вклад давления. На первом этапе в уравнениях движения учитываются только конвективные члены, в результате чего выделяется фиктивная переменная u*(x, t) и записываются выражения для предиктора и корректора скорости:
u* = un + fAt + Au +1 At/I / р, (5) un+1 = u* — Vpn+1At / р . (6)
Здесь и*+1/2 = 2(+1 + и*), ЛЬ = Ьп+1 — Ьп -
шаг по времени.
На втором этапе решается уравнение Пуассона для давления:
Лрп+1 = V ■ и*р / ЛЬ. (7)
Основной алгоритм движения по времени состоит из следующих шагов:
I) определяются границы области и строятся интерполяционные функции;
II) вычисляется предиктор скорости и*(х, Ь) из системы (5);
III) решается уравнение Пуассона (7) для определения давления рп+1(х, Ь);
IV) вычисляется новое значение скорости ип+1(х, Ь) из уравнения (6) с учетом найденного на шаге III давления;
V) вычисляется новое положение узлов на (п + 1) -ом временном шаге:
х(Ь)п+1 = хп (Ь) + ип+1(х, Ь)ЛЬ ;
VI) Если не достигнут критерий останова, осуществляется переход к шагу I.
Построение межузловых связей и определение границ расчетной области
Идея методов КЕМ (метод естественных соседей) и МБЕМ (бессеточный метод конечных элементов) заключается в применении несибсоновских интерполяций [44] или Лапласа [30], построение которых осуществляется на основе координат естественных соседей. В вычислительной геометрии концепция естественных соседей связана с понятием ячейки Вороного первого или второго порядка - выпуклого многоугольника, определяемого множеством точек на плоскости (в пространстве). Дискретизация области ячейками Вороного позволяет осуществлять быстрый поиск естественных соседей для заданной точки. Нахождение естественных соседей на основе информации о смежных узлах, полученных из диаграммы Вороного, существенно снижает временные затраты, необходимые для построения интерполирующих функций. Скорость выполнения этого шага является критической для КЕМ и МБЕМ, т. к. данные методы, основываясь на методе Галеркина, требуют вычисления функций формы для каждой точки при численном интегрировании по расчетной области задачи. Соответственно, скорость работы метода естественных соседей и бессеточного метода конечных элементов зависит от эффективности алгоритма, реализующего разбиение области ячейками Вороного, что особенно актуально при численном решении задач механики жидкости со свободными границами.
Для построения диаграммы Вороного используется алгоритм «sweepline» [16; 53].
На рис. 2 а приведено разбиение области ячейками Вороного, в которой положение частиц получено из решения нестационарной задачи о движении соли-тона по дну с расположенным на нем уступом в момент опрокидывания гребня волны [50]. Как видно из рисунка, граничные ячейки имеют бесконечные ребра
на плоскости. Для решения этой проблемы область заключается в прямоугольник большего размера и ребра отсекаются по границе этого прямоугольника.
Эффективным методом определения границ расчетной области является метод “ a -shape”, который основывается на понятии а -формы [15; 16].
Задача нахождения границы по заданному множеству точек подразумевает распознавание его формы и не имеет единственного решения. a -форма парамет-ризируется действительным числом a G [0; те) и представляет собой семейство геометрических фигур в диапазоне от точки до выпуклой оболочки множества. На практике для нахождения границы используется триангуляция Делоне, а параметр a является значением критерия: если радиус описанной окружности треугольника превышает значение a , то треугольник исключается из расчетной области.
Рисунок 2 иллюстрирует процесс определения границы множества точек. На первом шаге из диаграммы Вороного (рис. 2 а) восстанавливается первоначальная триангуляция Делоне (рис. 2 б). Затем осуществляется очистка области от элементов, не удовлетворяющих критерию a (рисунки 2 в - г). Для определения геометрии границы заданной области необходимо подобрать значение параметра a , дающего “приемлемые” результаты. Универсальных методов выбора значения а нет, существуют лишь некоторые рекомендации, полученные из опыта применения алгоритма и основанные на соотношениях расстояний между самыми близкими и самыми удаленными точками области.
Рис. 2. а) диаграмма Вороного; б) первоначальная триангуляция Делоне; в) триангуляция Делоне при а = 0.25; г) триангуляция Делоне при а = 0.059
МЕТОД ISPH
Основу метода ISPH составляет формула усреднения функции по Стеклову:
f (x) = Г f (x)W(x - x, h)dx.
D (8)
Присутствующая в интеграле весовая функция W, называемая функцией ядра или ядром усреднения, может выбираться произвольно, однако должна удовлетворять некоторому обязательному набору критериев:
1. W(х — х,h) = 0, ||x — х|| > kh .
2. J W (x — x, h)dx = 1.
D
3. lim W (x — x, h) = ö(x — x).
h 0
Здесь ö - дельта-функция Дирака, kh - радиус сглаживания (носитель функции ядра W ).
Для того чтобы вычислить интеграл (8), область расчета представляется набором частиц (малых объемов сплошной среды). Функции формы строятся на основе множества ближайших соседей - таких частиц, чьи координаты для произвольной частицы i определяются как Рг = {j G Md : || xi — xj || < kh } . Отыскав такое множество для каждой частицы, можно вычислить интеграл в формуле усреднения функции по Стеклову (8):
n
f(x) = J f(x)W(x — x, h)dx и £ fW(x — xj, h)AV3,
d j=1 (9)
которая, по сути, представляет собой SPH-аппроксимацию произвольной функции f(x). Построив на основе формулы (9) аппроксимации дифференциальных операторов, входящих в уравнения (1NS)-(4NS), можно получить их дискретный аналог.
В отличие от методов естественных соседей, в методах сглаженных частиц малый объем Д Vj вычисляется не геометрически (на основе ячеек сетки), а как отношение массы к плотности, первая из которых задается в качестве начального условия и не меняется в течение всего времени проведения вычислительного эксперимента.
При моделировании течений жидкости методом ISPH часто приходится сталкиваться с проблемой образования групп (кластеров) частиц по всей области расчета. Это серьезно сказывается на качестве вычисления кинематических и динамических характеристик исследуемых процессов. Причина кроется в особенностях большинства широко используемых функций ядра, которые, строго говоря, не подходят для вычисления градиента давления, поскольку значения их первых производных стремятся к нулю по мере приближения к центру области носителя. Эта их особенность приводит к ослабеванию сил отталкивания между парами частиц по мере их сближения, обеспечивая благоприятные условия для их кластеризации. Решением этой проблемы может стать использование для аппроксимации градиента давления такой функции ядра, поведение первой производной которой на
всей области ее носителя будет монотонным, что, соответственно, обеспечит монотонное возрастание отталкивающих сил между парами частиц по мере их сближения. В связи с тем, что в настоящей работе для проведения численных расчетов используется сплайн Морриса 4-ой степени с радиусом носителя равным 5 /2^, для вычисления давления в работе используется функция ядра подобного вида:
Г(5/2-д)4 -(3/2-д)4, 0<д <3/2; (до)
Ш(х-х,к)= 60 2 1(5/2-д)4, 3/2<д <5/2;
931жк I
[0, д > 5/2
Несмотря на то, что вычисление градиента давления на основе функции ядра (10) позволяет получать на каждом временном шаге относительно равномерную расчетную «сетку» внутри области течения, тем не менее, во многих задачах возникают все же проблемы, связанные с близким расположением частиц свободной поверхности к внутренним частицам среды, за счет отсутствия нужного количества ближайших соседей у первых. Это приводит к осцилляциям давлений в приграничных областях. Для решения этой проблемы предлагается алгоритм корректировки свободной границы, заключающийся в добавлении в правую часть уравнений движения дополнительных сил отталкивания, действующих между внутренними частицами и частицами свободной поверхности. В настоящей работе для этих целей используется потенциал упругих шаров:
Г1/2К(г-г0)2, г < г0;
П(г) = \ '-о'’ <”
[0, г > г0
где К - коэффициент жесткости взаимодействующих упругих шаров (в нашем случае - частиц жидкости), г - расстояние между ними, г0 - радиус действия потенциала. Данный подход позволяет стабилизировать вычисления поля давления и получить достоверные кривые гидродинамических нагрузок.
Для интегрирования системы уравнений (ШБ)-(4МБ) по времени в методе ШРН используется схема расщепления подобная (5)-(7), с той лишь разницей, что
и*+1/2 полагается равным и" .
НЕВЯЗКАЯ ЖИДКОСТЬ. МЕТОД КЕМ
Взаимодействие уединенной волны с подводной ступенькой
Вопросам взаимодействия уединенных волн с препятствиями посвящено большое количество теоретических и экспериментальных работ в связи с важностью вопросов по определению воздействия этих волн на гидротехнические и береговые сооружения. Физические эксперименты для изучения подобных явлений оказываются сложными и дорогостоящими, а быстрота протекания реальных процессов делает численные методы практически единственным источником информации о картине течения. В задачах, имеющих практический интерес, размеры длин волн обычно бывают большими, поэтому можно пренебрегать влиянием вязкости и рассматривать задачи в поста-
новке, основанной на модели идеальной несжимаемой жидкости.
Для численного моделирования задачи о взаимодействии уединенной волны со ступенькой используется модифицированный метод естественных соседей для решения системы уравнений Эйлера, позволяющий получать кинематические и динамические характеристики исследуемой задачи. Решению данной задачи посвящено достаточное количество работ, но большая часть из них описывает численное моделирование исследуемого явления на основе потенциальной модели идеальной несжимаемой жидкости [34, 43]. Решение системы уравнений Эйлера позволило обнаружить образование вихревых структур вблизи передней границы ступеньки, что также отмечено в экспериментальной работе [29, 49].
Пусть Б - область течения идеальной несжимаемой жидкости, ограниченная свободной поверхностью Г0 и твердыми границами Г (рис. 1).
Основными определяющими параметрами задачи являются амплитуда набегающей волны А и высота подводной ступеньки й. При накате набегающей волны на ступеньку на ее переднем фронте образуется волновой сгусток. Затем амплитуда волны начинает расти, на ее поверхности формируется двойной горб, который в дальнейшем разделяется на отраженную и прошедшую волны. Тестирование алгоритма решения задачи проводилось для таких значений параметров А и й, при которых не происходит обрушения прошедшей волны. Результаты решения сравнивались с численными [68; 69] и экспериментальными [29] данными. Численное моделирование показало, что при дальнейшем движении прошедшей волны по каналу
ее форма трансформируется, волна увеличивается по амплитуде, и от нее отходит четко сформировавшаяся вторая волна, бегущая вслед за первой (рис. 3 а).
Из рисунка 3 а видно, что в диапазоне
10 < х < 11 изменения абсциссы поверхность жидкости возмущена волновой рябью. Появление волновой ряби связано с образованием вихревого течения над ступенькой. Существование вихря отмечено в экспериментальной работе, однако его влияние на амплитуды прошедших и отраженных волн до сих пор не исследовано. Для определения влияния вихря на амплитуды прошедших и отраженных волн была проведена серия расчетов. Высота ступеньки задавалась равной й = 0.5 и й = 0.7, амплитуды волн -А = 0.18; 0.3; 0.4; 0.5; 0.6. Проведенная серия расчетов показала наличие вихрей вблизи передней границы ступеньки для всех рассматриваемых значений амплитуды волны и высоты ступеньки. На рисунке 3 б приведен фрагмент картины наката уединенной волны амплитуды А = 0.6 на подводную ступеньку высотой й = 0.7 до момента обрушения гребня волны в интервале времени 0 < Ь < 9.11. После момента обрушения происходит интенсивное перемешивание жидкости, сопровождающееся появлением брызг. Картины течений на рисунке 4 представлены множеством расчетных узлов. Гребень обрушающейся волны сильно бьет в подошву, выталкивая перед собой движущуюся вперед с большой скоростью массу жидкости, которая в дальнейшем также обрушается. В последующие моменты времени свободная поверхность подвержена сильному волнению, а движение жидкости представлено множеством взаимодействующих между собой волн различных форм и размеров.
а) б)
Рис. 3. Профиль свободной границы: а) А = 0.218; й = 0.5 (режим течение без опрокидывания волны), б) А = 0.6; й = 0.7 (режим течение с опрокидыванием волны)
Рис. 4. Картины течения в различные моменты времени: А = 0.6, (I = 0.7
Рис. 5. Гидродинамическая нагрузка на правую стенку бассейна: а) й = 0.5, б) й = 0.7
Динамическая нагрузка, создаваемая волной на правую твердую стенку, представлена на рисунке 5 для высоты ступеньки ( = 0.5 и ( = 0.7 и амплитуды А = 0.4; 0.5; 0.6 . На рисунке 6 а приведена картина вихревого течения, а на рисунках 6 б и в представлено значение циркуляции вихря над ступенькой, отнесенное к глубине жидкости ¿0 для различных амплитуд набегающей волны и высоты ступеньки ( = 0.7. Максимальное значение модуля циркуляции достигается в момент отделения отраженной волны от набегающей.
Рис.6. а) вихревое течение: А = 0.5; ( = 0.7 ; б), в) циркуляция вихря над ступенькой
Результаты проведенной серии расчетов сравнивались с результатами, полученными КМГЭ на основе потенциальной модели идеальной жидкости в работе [49] по следующим параметрам: форма свободной поверхности, максимальное значение амплитуды волны Ая при накате на ступеньку, амплитуде отраженной волны Аг, максимальной амплитуде волны АЬ перед моментом обрушения.
На рисунке 7 в области изменения координаты 8 < х < 10 виден значительный изгиб свободной
границы в сторону дна для расчета методом №Ж. Подобное поведение свободной границы наблюдается для всех проведенных расчетов, параметры которых указаны выше. Причиной изгиба свободной границы является вихревое течение, образующееся над ступенькой.
Рис. 7. Профили свободной границы: кривая 1 — КМГЭ, кривая 2 — НЕМ
В результате сравнения был выявлен интересный факт, касающийся амплитуды отраженной волны: для значений амплитуд набегающей волны А > 0.4 амплитуда отраженной волны, полученной в результате численного моделирования методом КЕМ (кривая 1), превышает значение амплитуды отраженной волны для метода КМГЭ (кривая 2) на 15 - 20 % (рис. 8).
При этом значения амплитуд Ам- и А-1 отличаются незначительно (менее 5 %).
В дальнейшем условия задачи меняются: вместо уступа рассматривается расположенное на дне тело прямоугольной формы. Исследуется влияние ширины Ь и высоты препятствия ё, а также амплитуды волны на основные характеристики возникающего течения жидкости. Расчеты проводились для Ь = 2; 7, ё = 0.3; 0.5; 0.7 и амплитуд волн А = 0.2; 0.3; 0.4; 0.5; 0.6 в области, где —8 < х < 21. Решению данной задачи посвящена работа [68], в которой рассматриваются режимы движения без обрушения гребня волны при прохождении над препятствием. На рисунке 9 приведены фрагменты течения в различные моменты времени, полученные методом естественных соседей, для амплитуды А = 0.5 и размеров тела Ь = 7; й = 0.7. В области задано 14000 расчетных узлов. В зависимости от значений размеров тела и амплитуды набегающей волны обрушение можно классифицировать как скользящий или ныряющий бурун. При набегании уединенной волны на препятствие квадратного сечения волновые картины взаимодействия весьма схожи с теми, которые возникают, когда подводным препятствием служит полукруговой выступ на дне. Серия расчетов для различных амплитуд и Ь = ё = 0.3; 0.5; 0.7 показала, что в таком случае на заднем фронте волны формируется всплеск, который опрокидывается “против движения” основной волны.
Рис.8. Амплитуда отраженной волны
б)
]' 1 \~т .V. # .... ■„ V'', ' ,*■ ' /.,':у
¡1 "V/ \ *
5 10 1? А
в)
г)
Рис. 9. Картина течения в различные моменты времени: А = 0.5; Ь = 7; й = 0.7
Рис. 10. Вихри вблизи препятствия
В результате проведенной серии расчетов методом естественных соседей установлено наличие вихрей вблизи препятствия, среди которых можно выделить два, обладающих наибольшей циркуляцией (рис. 10). Один из этих вихрей располагается над передней стенкой препятствия аналогично случаю, возникающему при набегании волны на подводную ступеньку. Второй вихрь находится непосредственно за препятствием. На рисунке 11 приведены графики циркуляции этих вихрей в зависимости от высоты ё для амплитуды набегающей волны А = 0.5. Как видно из рисунков, вихрь, располагающийся за препятствием (рисунок 11 б), обладает большей интенсивностью, которая довольно продолжительное время остается на одном уровне. В результате проведенных расчетов установлено (рис. 12), что с увеличением амплитуды набегающей волны амплитуда отраженной волны для метода КЕМ (кривая 1) превышает значение амплитуды для КМГЭ (кривая 2) на 40 — 60 %.
Рис. 12. Амплитуда отраженной волны
ВЯЗКАЯ ЖИДКОСТЬ. МЕТОД ОКЕМ
Задача о разрушении плотины
Представим решение модельной задачи о разрушении плотины. Задача формулируется следующим образом (рис. 13).
Рис. 11. Циркуляция вихрей вблизи препятствия
Рис. 13. Схема расчетной области
В начальный момент времени прямоугольный столб вязкой несжимаемой жидкости находится в состоянии покоя. Под действием силы тяжести столб начинает обрушаться. В определенный момент времени жидкость ударяется о противоположную вертикальную стенку бассейна, а затем волна, образовавшаяся при откате от стенки, обрушается, нарушая связность расчетной области.
В таблице 2 приведены значения параметров, которые использовались для расчета задачи.
Таблица 2
Параметры задачи о разрушении плотины
Параметр Значение
Плотность жидкости 1000 кг/м3
Коэффициент динамической вязкости Порядок числа Рейнольдса 3.2 ' 10-2 кг/(м-с) 4 • 104
Ускорение силы тяжести 9.81 м/с2
Длина резервуара 5.84 м
Ширина столба жидкости 1.46 м
Высота жидкости в верхнем бьефе 2.92 м
Как уже было сказано выше, при расчете подобных задач классическими сеточными методами лагранже-
вой природы происходит аварийное завершение расчетов.
Численные результаты. Кинематическая картина
На рис. 14 приведены картины течений для различных моментов времени с количеством расчетных частиц равным 2500. При t = 0.077 с наступает процесс обрушения отраженной от твердой стенки волны, и дальнейшее моделирование течения жидкости класси-
ческими сеточными лагранжевыми методами становится невозможным, однако в№М успешно справляется с данной задачей (рис. 14 г). Из рисунка 14 видно, что течение жидкости в данной задаче представлено сложной многосвязной областью с сильно нелинейной геометрией свободной поверхности, что показывает численную устойчивость метода естественных соседей при решении такого рода задач.
а)
б)
в)
Рис. 14. Задача о разрушении плотины. 2500 узлов: а) і = 0 с, б) і = 0.2 с, в) і = 0.4 с, г) і = 0.8 с
На рисунке 15 для сравнения показаны фрагменты течений для 4500 частиц. Рисунки 14 и 15 демонстрируют хорошее качественное совпадение результатов.
а)
б)
в)
Рис. 15. Задача о разрушении плотины. 4500 узлов: а) і = 0 с, б) і = 0.2 с, в) і = 0.6 с, г) і = 0.8 с
Вычисление гидродинамических характеристик. Вычисление нагрузки
Актуальность вычисления нагрузок обусловлена необходимостью определения характеристик защитной стенки для предотвращения ее разрушения при воздействии на нее волны прорыва, образующейся
при квазимгновенном разрушении вертикальных стен резервуаров.
Для задачи об обрушении плотины вычислялись нагрузки (рисунок 16) на правую (кривая 1) и левую (кривая 2) вертикальные стенки бассейна для области
Рис. 16. Хронограммы гидродинамических нагрузок на правой (кривая 1) и левой (кривая 2) стенках области:а) 2500улов, в) 4500узлов
Характерным для такой расчетной области является то, что длина бассейна более чем в 3 раза превышает ширину столба жидкости. Видно, что в момент взаимодействия волны, формирующейся в результате обрушения плотины, с правой границей области нагрузка на ней достигает максимального значения. Также можно отметить, что расчетные хронограммы нагрузок на правую стенку имеют два локальных максимума. Это явление можно объяснить действием сил инерции. Кроме того, первый максимум имеет большее значение, чем второй, при этом момент максимального заплеска на правую стенку не совпадает с моментами локальных максимумов давления. Эти особенности были отмечены в экспериментах, полу-
ченных С. В. Манойлиным при накате солитонов на вертикальную стенку [58]. При дальнейшем откате волны от правой стенки нагрузка на ней начинает убывать.
Для дальнейшего анализа нагрузок на правую стенку области в задаче об обрушении столба жидкости варьируемыми параметрами задачи были выбраны высота столба жидкости Н и длина бассейна Ь. Группа численных расчетов выполнялась для следующего диапазона параметров: Ь є [1,3] метров. На рисунке 17 представлено сравнение кривых нагрузок на вертикальные стенки области для различных значений длин бассейна.
Рис. 17. Хронограммы гидродинамически а) правая стенка области
При Ь = 1 м расстояние от правого края столба жидкости до правой границы совпадает с шириной столба. Кривая с номером 1 соответствует Ь = 1 м, с номером 2 - Ь = 2 м ис номером 3 - Ь = 3 м (рис. 17).
Расчеты распространения волны, формирующейся при обрушении столба жидкости, показали, что длина бассейна, кроме случая, когда расстояние до правой стенки не превышало ширины столба жидкости, слабо влияет на характер гидродинамической нагрузки на стенку бассейна. Из приведенных рисунков видно, что при различных значениях длины бассейна нагрузки на вертикальные стенки принимают близкие значения. Это можно объяснить тем, что в случае, когда размер бассейна превышает две ширины столба жидкости, образующаяся в результате распада волна те-
б)
; нагрузок: Ь = 1м, Ь = 2 м, Ь = 3 м. б) левая стенка области
ряет больше энергии до момента наката на правую стенку, чем в случае, когда вертикальная стенка находится ближе к правой границе столба жидкости. Тогда расстояние, которое необходимо пройти волне излива до соприкосновения с правой стенкой, сокращается.
Обрушение плотины при наличии слоя жидкости на основании
Расчетная область представляет собой бассейн с ровным дном и твердыми непроницаемыми стенками, заполненный однородной вязкой несжимаемой жидкостью и разделенный в начальный момент времени тонкой непроницаемой перегородкой, создающей перепад уровня жидкости (рис. 18). В начальный момент времени перегородка начинает равномерно двигаться вверх с заданной скоростью, формирующийся
при этом столб жидкости начинает обрушаться под действием силы тяжести.
Задача является неустойчивой, если используется постоянный шаг по времени. Шаг по времени при расчете данной задачи, как и при задаче об обрушении столба жидкости, выбирался из условия (8).
О 0.2 0.4 X 0.6 0.8 1
Таблица 3
Параметры задачи о разрушении плотины при наличии слоя жидкости на дне
Параметр Значение
Скорость крышки плотины 1,5 м/с
Плотность жидкости 1000 кг/м3
Коэффициент динамической вязкости Порядок числа Рейнольдса СМ ^ 15 3к
Сила тяжести 9,81 м/с2
Высота слоя жидкости при основании Высота столба жидкости Ширина столба жидкости 1,8 ■ 10-2 м 0.15 м 0.38 м
Рис. 18. Схема расчетной области
В таблице 3 приведены значения параметров, которые использовались для расчета поставленной задачи.
Численные результаты. Кинематическая картина
В работе [62] представлено сравнение полученных картин течения задачи об обрушении плотины при наличии слоя жидкости на основании для различного числа узлов области.
На рисунке 19 приведены картины течений для различных моментов времени с количеством расчетных частиц 4030.
Рис. 19. Задача о разрушении плотины. 4030 узлов: а) і = 0.156 с, б) і = 0.219 с, в) і = 0.281 с,
г) і = 0.343 с, д) і = 0.406 с, е) і = 0.468 с
МАТЕМАТИКА |
На рисунке 20 приведено сопоставление полученных результатов с данными эксперимента из работы [11].
Рис. 20. Сравнение результатов работы авторов с экспериментальными данными в моменты времени:
а) Ь = 0.156 ^ б) Ь = 0.343 ^ в) Ь = 0.406 ^ г) Ь = 0.468 c
На рисунке 21 представлено графическое наложение результатов численного расчета задачи об обрушении плотины методом естественных соседей с экспериментальными данными. Внешней сплошной линией на рисунке нанесены оцифрованные экспериментальные данные.
а)
в)
д)
Рис. 21. Задача о разрушении плотины. Сравнение с экспериментом: а) Ь = 0.156 ^ б) Ь = 0.219 ^ в) Ь = 0.281 ^ г) Ь = 0.343 ^ д) Ь = 0.406 ^ е) Ь = 0.468 c
Вычисление гидродинамических характеристик. Вычисление гидродинамического давления На рисунке 22 представлены характерные картины распределения поля давления во всей области течения в различные моменты времени.
а)
б)
в)
г)
Рис. 22. Распределение поля давления: а) і = 0.156 с, б) і = 0.219 с, в) і = 0.281 с, г) і = 0.343 с,
д) і = 0.406 с, е) і = 0.468 с
Вычисление гидродинамических характеристик. Вычисление нагрузки
Для значений Ь = 0.018 м, Ь = 0.038 м, Ь = 0.09 м и Ь = 0.12 м вычислялись нагрузки на левую (рис. 23 а) правую (рис. 23 б) вертикальные стенки бассейна.
При Ь = 0.12 м толщина слоя практически соответствует высоте столба жидкости и выбрана эмпирически на основе проведенной серии расчетов с тем, чтобы для заданной длины бассейна не нарушалась связность области и не возникали моменты обруше-
ний. Соответственно график нагрузок для Ь = 0.12 м наиболее гладкий и имеет наименьшее число осцилляций в отличие от хронограмм нагрузок для других значений толщины слоя при основании. При значении высоты столба Ь = 0.09 м момент обрушения свободной границы возникает уже после отката образующейся в результате разрушения плотины волны излива, поэтому на хронограмме максимальное значение нагрузки лежит ниже по сравнению с режимами, для которых моменты обрушений возникают до наката волны излива на вертикальную стенку.
Рис. 23. Хронограммы гидродинамических нагрузок на вертикальные стенках области: Для значений Л = 0.018 м (кривая 1), Л = 0.038 м (кривая 2), Л = 0.09 м (кривая 3), Л = 0.12 м (кривая 4). а) правая стенка области, б) левая стенка области
Взаимодействие волнового движения с днищем плавающих тел
Подобные задачи представляют собой практический интерес, например, при моделировании посадки летательных аппаратов на воду. Определение возможных нагрузок, создаваемых плесканием жидкости на границах внешних конструкций, расположенных над свободной поверхностью жидкости, необходимо при моделировании движения этих конструкций с учетом их прочности.
Как было показано ранее, в процессе обрушения столба жидкости формируются ударные волны, амплитуда и скорость движения которых зависит от многих факторов, таких как высота столба жидкости, глубина слоя жидкости при основании, скорость движения заслонки плотины (в задаче о разрушении плотины). Подобные волны используются для моделирования процесса наката волны на горизонтальную преграду. В рассматриваемой далее задаче о разрушении плотины при наличии слоя жидкости на основании со стороны правой вертикальной стенки дополнительно добавляется горизонтальный уступ (рис. 24). Параметры жидкости выбираются такие же, как и в предыдущей задаче.
Высота слоя жидкости при основании к равняется 0.018 м. Шаг по времени переменный, выбирается из условия (8).
Приведем картины течений в различные моменты времени для случая, когда высота горизонтальной поверхности в четыре раза превышает глубину слоя. Та-
кое значение параметра Л было выбрано по причине того, что, в этом случае при обрушении столба жидкости, максимальная амплитуда волны в четыре раза превышает начальное значение глубины слоя при основании. При этом длина горизонтальной преграды была выбрана такой, чтобы момент максимальной амплитуды волны совпадал с моментом столкновения волны с уступом.
На рис. 25 представлены характерные картины течения в различные моменты времени.
Из представленной картины видно, что после удара о преграду, жидкость заполняет уступ таким образом, что удары формирующихся в процессе движения волн по горизонтальной границе возникают только до момента, когда жидкость достигнет правой стенки бассейна. Далее воздействие на преграду остается равномерным и зависит от скорости и завихренности течения.
0.3
о 0.2 0.4 Чй 0.6 0.8 ;1;
Рис. 24. Схема расчетной области
Для оценки возможной силы удара жидкости по горизонтальной поверхности была проведена серия расчетов для значений высоты преграды кд = 4Л, 6Л, 8Л, 10Л, 12Л, где Л - высота слоя жидкости при основании. На рис. 26 приведены картины течения в основные моменты взаимодействия волны обрушения с горизонтальным уступом для некоторых значений параметра Лд .
Из результатов расчетов следует, что при 1гд = 4Н и более волна обрушения не соприкасается с
горизонтальной поверхностью, а их взаимодействие
происходит только при откате волны от правой вертикальной стенки (рис. 27).
В случае, когда горизонтальная преграда находится на достаточной высоте от свободной поверхности жидкости, наблюдается плескание жидкости по преграде. Так как жидкость вязкая, то на твердой границе ставится условие прилипания, а при «стекании» жидкости с горизонтальной поверхности вниз происходит «отрыв» частиц жидкости от твердой границы. В моменты «отрыва» возникают так называемые зоны отрицательного давления (рис. 28).
а)
б)
0.2
0.1
0.6 х 0 8 1
в)
Рис. 26. Взаимодействие волны с горизонтальным уступом: а) Лд = 4 Л ; б) Лд = 6 Л ; в) Лд = 8 Л
Рис. 27. Задача о взаимодействии волны с горизонтальным уступом (Лд = 10 Л): а) Ь = 0.65 с, б)
Ь = 0.729
Рис. 28. Задача о взаимодействии волны с горизонтальным уступом (кд = 12к): а) Ь = 0.69 с, б) Ь = 0.75 в), г) фрагменты, характеризующие наличие отрицательного (отрывного) давления
Вычисление нагрузки
Для различных значений параметра кд (высота горизонтальной преграды от свободной поверхности жидкости) вычислялись нагрузки на горизонтальный уступ и правую вертикальную стенку бассейна.
а) б) в)
Рис. 29. Хронограммы гидродинамических нагрузок: кд = 4 к (кривая 1), кд = 6 к (кривая 2), кд = 10 к (кривая 3), кд = 12 Н (кривая 3). а) правая вертикальная граница, б) горизонтальная граница уступа, в) вертикальная граница уступа
На рис. 29 приведены хронограммы нагрузок на вертикальные границы области для значений 1гд = 4Н, !гд = 6Н, !гд = 8Н и !гд = 12Н соответственно на правой вертикальной стенке бассейна (рис. 29 а) горизонтальной границе уступа (рис. 29 б) и вертикальной границе уступа (рис. 29 в).
При !гд = 4Н и Нд = 6Н после наката уступ полностью заполняется жидкостью, отсюда можно наблюдать резкий скачок нагрузки, соответствующий моменту столкновения волны обрушения с уступом.
При этом при кд = 4к основной удар приходится на
вертикальную стенку уступа (рис. 29 в). С увеличением значения высоты уступа над поверхностью жидкости пики на графике нагрузки, соответствующие ударам жидкости по горизонтальной границе, появляются уже в момент отката жидкости от правой границы бассейна и принимают меньшее значение по сравнению с ударом жидкости по уступу для значений параметра 1гд = 4к и !гд = 6к.
ВЯЗКАЯ ЖИДКОСТЬ. МЕТОД КРН
Сравнение методов КРН и ОКЕМ
Численное решение задачи о разрушении плотины также было получено с использованием метода ШРН. Сравнение картин течения и поля давления для
обоих методов (ШРН и вМБМ) представлено на рисунке 30.
Рис. 30. Сравнение картин течения и поля давления: а, в, д, ж) метод ОМЕМ, б, г, е, з) метод ШРИ
Ниже также приведено сравнение хронограмм гидродинамических нагрузок (рис. 31). Рисунки 30 и 31 свидетельствуют о качественном совпадении результатов, а присутствующие в полученных результатах различия могут быть обоснованы тем, что при мо-
делировании процесса разрушения плотины методом ШРН использовалась модель турбулентности Прандт-ля. В расчетах с использованием метода вМБМ турбулентные эффекты не учитывались.
а)
б)
Рис. 31. Хронограммы гидродинамических нагрузок: а) метод ОНЕШ, б) метод №РИ
ПРОЦЕССЫ ВЗАИМОДЕЙСТВИЯ ЖИДКОСТИ С ПОГРУЖЕННЫМИ ТЕЛАМИ
Уравнения движения твердого тела
Известно, что скорость движения твердого тела
складывается из скорости vc поступательного движения его центра масс и скорости ы его вращательного движения относительно некоторого полюса (центра масс тела):
V = Vc + Ы X х*, где X* = X - Хс .
Скорость поступательного движения тела, погруженного в жидкость, может быть найдена из второго закона Ньютона:
(IV р р
М& - J рп(Г, + J Т ■ т(Г, ,
Г, г,
И-
Іі
где интегралы в правой части уравнения представляют собой силы давления и силы вязкого трения, действующие на твердое тело со стороны жидкости. Скорость вращательного движения тела может быть получена аналогичным образом из уравнения сохранения момента импульса:
1(ы = х * X(рпГ, + ^х * х(т ■ tГ, ,
(Ь г, г,
Для чего необходимо знать момент инерции тела:
I = Рь/||х *\\ІУЬ.
Значения параметров задачи приведены в таблице 4.
Таблица 4
Параметры задачи о всплытии цилиндра в жидкости
Параметр Значение
Плотность жидкости 1000 кг/м3
Плотность цилиндров 500 и 750 кг/м3
Коэффициент динамической вязкости 10-3 кг/(м-с)
Ускорение силы тяжести 9.81 м/с2
Длина бассейна 0.12 м
Глубина бассейна 0.02 м
Радиус основания цилиндра 0.005 м
Начальное положение центра масс цилиндра (0.06, 0.01) м
Здесь п и і - нормальный и касательный к поверхности тела векторы, Гь - поверхность тела.
Задача о всплытии цилиндра в жидкости
В бассейн с жидкостью погружен бесконечный абсолютно твердый круговой цилиндр меньшей по сравнению с жидкостью плотности. Жидкость и цилиндр находятся в состоянии покоя. В начальный момент времени начинается процесс всплытия цилиндра за счет действия на него архимедовой силы. Расчетная область в начальный момент времени представлена на рисунке 32.
Задача обладает свойством симметрии относительно плоскости, проходящей через ось цилиндра ортогонально покоящейся свободной поверхности жидкости в начальный момент времени. Симметричность постановки задачи означает, что моменты сил, действующих на цилиндр, равны нулю. Тем не менее, задача решается в полной двумерной постановке.
Численные результаты
На рисунке 33 приведено сравнение картин течения и поля давления для цилиндров плотности 500 кг/м3 (слева) и 750 кг/м3 (справа). В расчетах принимало участие 17376 частиц. Представленные результаты соответствуют следующим моментам времени: 0.044 с (а, б), 0.108 с (в, г), 0.232 с (д, е), 0.3 с (ж, з).
0.04
0.04 ■ г
■1 | шШш Ш и
Р: 10 40 70 100 130 160 190 0.02 I Р: 10 40 70 100 130 160 190
0.12 б) 0 0.02 0.04 0.06 0.08 0.1 0.12
0.04 I
0.02 Г
О о|_________________________________________________________________________________
в) О 0.02 0.04 0.06 0.08 0.1 0.12 г) 0 0.02 0.04 0.06 0.08 0.1 0.12
0.04
0.02
0.04
0.02
О И- . I , I ... I ... I....................I . . . I 01______________________________________________________
Д) 0 0.02 0.04 0.06 0.08 0.1 0.12 е) 0 0.02 0.04 0.06 0.08 0.1 0.12
0.04
0.02
0.04
0.02
ш и
Р: 10 40 70 100 130 160 190
ж) 0 0.02 0.04 0.06 0.08 0.1 0.12 3) 0 0.02 0.04 0.06 0.08 0.1 0.12
Рис. 33. Картины течения и поле давления в задаче о всплытии цилиндра в различные моменты времени:
а, в, д, ж) р = 500 кг/м3, б, г, е, з) р = 750 кг/м3
Ниже приведены графики зависимости от времени (рис. 34) гидродинамической силы, действующей на первый - 500 кг/м3 (а) и второй - 750 кг/м3 (б) цилиндры. Гидродинамическая сила представляет собой сумму поверхностных сил давления и вязкого трения. Максимальные значения силы соответствуют моментам наибольшего погружения цилиндров в жидкость,
минимальные - их наивысшему подъему. Осцилляции гидродинамической силы в начальный момент времени обусловлены особенностями численного счета, поскольку в начальный момент времени происходит резкий переход от первоначально заданного гидростатического распределения давления к гидродинамическому.
Рис. 34. Гидродинамическая сила, действующая на цилиндр: а) р = 500 кг/м3, б) р = 750 кг/м3
Задача о входе и погружении цилиндра в жидкость
В бассейн, заполненный жидкостью, первоначально находящейся в состоянии покоя, падает бесконечный абсолютно твердый цилиндр. Задача решается для цилиндров с круглой (рис. 35) и квадратной формами основания. Движение цилиндра с квадратным основанием происходит вдоль диагонали его основания, т.е. процесс входа на начальных этапах аналогичен входу в жидкость клина с углом килеватости 45°. Из условия задачи следует, что она обладает тем же свойством симметрии, что и предыдущая, и также решается в полной двумерной постановке.
Рис. 35. Схема расчетной области
Параметры для проведения численных расчетов приведены в таблице 5.
Таблица 5
Параметры задачи о входе и погружении цилиндра в жидкость
Параметр Значение
Плотность жидкости 1000 кг/м3
Плотность цилиндров 500 кг/м3
Коэффициент динамической вязкости 10-3 кг/(м-с)
Ускорение силы тяжести 9.81 м/с2
Длина бассейна 0.12 м
Глубина бассейна 0.02 м
Радиус (длина стороны) основания цилиндра 0.005 м (0.01 м)
Начальное положение центра масс цилиндра (0.06, 0.03) м
Численные результаты
На рисунке 36 приведено сравнение картин течения с раскраской по полю давления, полученные для кругового цилиндра (количество расчетных частиц -17885) и цилиндра с квадратным основанием (количество частиц - 17905). Фрагменты течения соответствуют моментам времени: 0.044 с (а), 0.072 с (б), 0.3 с (в).
0.04
0.02
0.04
0.02
: + . I_____I . I I 1,1 о|
0 0.02 0.04 0.06 0.08 0.1 0.12 0 0.02 0.04 0.06 0.08 0.1 0.12
а) б)
0.04
0.02
В)
0.04
0.02
0.12
Г)
0.12
Рис. 36. Картины течения и поле давления в задаче о входе и погружении цилиндра в жидкость в различные моменты времени: а, в, д) цилиндр с круглым основанием; б, г, е) цилиндр с квадратным основанием
Далее приводятся графики зависимостей от времени гидродинамической силы, действующей на цилиндры (рис. 37). В моменты времени 1 = 0.03 с (для
кругового цилиндра) и 1 = 0.024 с (для цилиндра с квадратным основанием) сила, действующая на цилиндры резко возрастает, что соответствует момен-
там удара цилиндров о поверхность заполняющей бассейн жидкости. Затем происходит некоторый спад нагрузок, после чего они возрастает вплоть до максимального погружения цилиндров в жидкость. Следующая затем стадия всплытия характеризуется постепенным уменьшением действующих на цилиндры сил, пока не будет достигнут максимальный подъем. Затем нагрузки со стороны жидкости вновь начинают возрастать по мере погружения в нее цилиндров. В течение вычислительного экспе-
римента значения нагрузки на круговой цилиндр принимают меньшие значения, чем на цилиндр с квадратным основанием, что, по-видимому, и обусловливает большую глубину повторного погружения кругового цилиндра. Несмотря на то, что задачи о всплытии и погружении цилиндров в жидкость решались в полной двумерной постановке, результаты вычислительных экспериментов демонстрируют высокую степень симметрии.
і -
I
а)
Рис. 37. Гидродинамическая си а) с круглым основанием, б)
Далее рассматриваются процессы входа и погружения в жидкость цилиндров различной массы. Параметры, при которых проводились численные расчеты, приведены в таблице 6.
Таблица 6
Параметры задачи о входе и погружении цилиндрав жидкость
Параметр Значение
Плотность жидкости 1000 кг/м3
Плотность цилиндров 250, 500 и 750 кг/м3
Коэффициент динамической вязкости 10-3 кг/(м-с)
Ускорение силы тяжести 9.81 м/с2
Длина бассейна 1 м
Глубина бассейна 0.3 м
Радиус основания цилиндра 0.055 м
t
б)
га, действующая на цилиндр: с квадратным основанием
Общее количество расчетных частиц составило 14515. Задача решается в постановке, использованной в экспериментальной работе [18]. Полученные методом сглаженных частиц результаты вычислительных экспериментов сравниваются с результатами численных расчетов методом СІР, представленными в работе [33], и данными лабораторных экспериментов [18].
Численные результаты
На рисунке 35 а, б, в приведены картины течений при входе и погружении в жидкость цилиндра плотности рЬ = 500 кг/м3 в момент времени 1 = 0,005 с, на
рисунке 35 г, д, е - в момент времени 1 = 0,030 с, на рисунке 35 ж, з, и - в момент времени 1 = 0,085 с, на рисунке 35 к, л, м - в момент времени 1 = 0,120 с. За начальный момент времени выбран момент касания цилиндром свободной поверхности жидкости в бассейне.
Результаты сравнения показывают, что метод сглаженных частиц достаточно адекватно описывает процесс погружения тела в жидкость.
а)
л)
Рис. 38. Картины течения при погружении в жидкость цилиндра плотности р = 500 кг/м : а, г, ж, к) метод ISPH; б, д, з, л) метод CIP; в, е, и, м) эксперимент
На рисунке 39 представлен график зависимости от времени глубины погружения цилиндра в жидкость. Глубина рассматривается как разность начального невозмущенного положения свободной поверхности жидкости (у = 0.3 м) в бассейне и вертикальной координаты нижней точки поверхности цилиндра. Пунктирная линия - результат, полученный методом сглаженных частиц, сплошная линия - методом СІР в работе [33], квадраты - результаты эксперимента [18].
Рис. 39. Глубина погружения цилиндра в жидкость
Для детального изучения процессов волнообразования при входе, погружении и последующем всплы-
тии цилиндров в жидкости, в бассейне были расположены «волномеры» вдоль прямых х = 0.8 м и х = 0.9 м (рис. 40). Кроме того, измерялся заплеск жидкости на правую вертикальную стенку бассейна, т.е. вдоль прямой х = 1 м.
На рисунке 41 а, б, в представлены графики зависимости от времени положения свободной поверхности жидкости в бассейне. Здесь и далее на графиках кривая 1 соответствует цилиндру плотности 250 кг/м3, кривая 2 - цилиндру плотности 500 кг/м3 и кривая 3 - цилиндру плотности 750 кг/м3. Благодаря своей малой массе, цилиндр плотности 250 кг/м3 за время своего падения не приобретает импульса, достаточного для его полного погружения в жидкость. Напротив, более тяжелые цилиндры (плотности 500 кг/м3 и 750 кг/м3) погружаются на достаточную глубину, чтобы образовавшаяся на поверхности жидкости в результате процесса погружения каверна схлопнулась, образовав вертикальную кумулятивную струю.
После удара цилиндра о поверхность жидкости в бассейне происходит формирование первых волн, расходящихся от цилиндра в направлении вертикальных стенок бассейна, которые на рисунке 41 обозначаются как тах1.
После всплытия за счет повторного погружения в жидкость цилиндров плотности 250 кг/м3 и 500 кг/м3 происходит формирование вторых волн (тах2).
Иной механизм формирования вторых волн можно наблюдать в случае цилиндра плотности 750 кг/м3: в этом случае волны формируются за счет распада кумулятивной струи, образующейся в результате схлопывание каверны над цилиндром во время его погружения, следующего за ударом о поверхность жидкости. Несмотря на различный механизм формирования вторых волн, можно отметить, что во всех случаях первый «волномер» (х = 0.8 м) отмечает большие значения их амплитуд в сравнении с первыми. При этом измеренные обоими «волномерами» амплитуды первых волн одинаковы. В случае же вторых волн можно наблюдать уменьшение их амплитуд на уровне второго «волномера» (х = 0.9 м). Связано это с тем, что на подходе к нему происходит их обрушение. Уровень заплеска на стенку бассейна во всех случаях тем выше, чем больше масса погружающегося в жидкость цилиндра. На рисунке 41 г приведены кривые гидродинамических нагрузок на правую стенку бассейна, обусловленные накатом на нее сгенерированных цилиндрами волн.
Рис. 41. Уровень свободной поверхности: а) первый волномер, б) второй волномер, в) третий волномер (правая стенка бассейна). г) гидродинамическая нагрузка на правую стенку области
Заключение
Результаты проведенных вычислительных экспериментов демонстрируют эффективность применения бессеточных методов для решения задач, в которых основной интерес представляют не пиковые нагрузки на твердые границы тела, а процессы волнообразования и гидродинамические нагрузки на конструкции, расположенные вблизи мест взаимодействия с жидкостью.
Решаемый круг задач требует для своего решения значительных вычислительных ресурсов, которые доступны при использовании СуперЭВМ. Так как программный код, реализующий используемые для решения поставленных задач современные численные методы, объемный и сложный по структуре, то реализация методов для многопроцессорных систем является не менее актуальной задачей. Использование параллельной реализации ускоряет проведение численных расчетов и дает значительно больший объем информации за счет существенной детализации параметров расчетов.
В процессе разработки параллельной программной реализации бессеточных методов была создана библиотека параллельных вычислительных программ для задач гидродинамики со свободными границами «HydroParaLib» [54; 63].
В настоящий момент на кафедре ЮНЕСКО по НИТ ведется подготовка студентов по направлениям бакалавриата и магистратуры: прикладная математика, математическое обеспечение и администрирование информационных систем, фундаментальная информатика и информационные технологии. Также существует действующая аспирантура по специальностям: 01.02.05 - механика жидкости, газа и плазмы; 05.13.18 - математическое моделирование, численные методы и комплексы программ. Под руководством заведующего кафедрой, профессора К. Е. Афанасьева защищено 13 кандидатских диссертаций, 6 из них по специальности 01.02.05 и 4 - по 05.13.18. Основное направление отмеченных научных исследований -моделирование задач со свободными поверхностями и высокопроизводительные вычисления.
Результаты работы научного коллектива неоднократно представлялись на тематических конференциях самого высокого уровня, среди которых следует отметить: “Гидродинамика больших скоростей и математическое моделирование” (Кемерово, 2006),
XV Всероссийская научно-методическая конференция «Телематика - 2008», III-я Всероссийская конференция с участием зарубежных ученых «Задачи со свободными границами: теория, эксперимент и приложения» (Бийск, 2008), X Международная летняя научная школа «Гидродинамика больших скоростей» и Международная научная конференция «Гидродинамика. Механика. Энергетические установки» (Чебоксары, 2008), 24th International Workshop on Water Waves and Floating Bodies (Зеленогорск, 2009), International conference «Mathematical and Informational Technologies» (Сербия, г. Копаоник, Черногория, г. Будва, 2009), Девятая международная конференция-семинар «Высокопроизводительные параллельные вычисления на кластерных системах» (Владимир, 2009), “Всероссийский съезд по теоретической и прикладной механике”
(Нижний Новгород, 2006 и 2011), X Всероссийская конференция «Прикладные технологии гидроакустики и гидрофизики» (Санкт-Петербург, 2010),
XVII Всероссийская научно-методическая конференция «Телематика - 2010» (Санкт-Петербург, 2010), VII Международная конференция «Лаврентьевские чтения по математике, механике и физике» (Новосибирск, 2010), X Всероссийский съезд по фундаментальным проблемам теоретической и прикладной механики (Нижний Новогород, 2011), “Задачи со свободными границами: теория, эксперимент и приложения” (Новосибирск - Бийск, 2009 и 2011), XI Всероссийская конференция «Прикладные технологии гидроакустики и гидрофизики» (Санкт-Петербург, 2012), “Информационно-вычислительные технологии и математическое моделирование” (Кемерово, 2013), Second International Workshop “Hydrodynamics of Moving Objects” (Kiev 2013) и др. Полное описание результатов исследовательской работы научной школы изложено в сборнике: “Информационные и вычислительные технологии в численных расчетах и управлении вузом” [1] и в работах [2; 3].
Основные результаты работы были использованы при выполнении следующих проектов.
Проект № 4829 «Численное моделирование течений жидкости со свободными границами современными численными методами на многопроцессорных вычислительных системах» (2005 год) по ведомственной научной программе федерального агентства по образованию «Развитие научного потенциала высшей школы».
Интеграционный проект фундаментальных исследований Объединенного ученого совета по механике и энергетике СО РАН (2006-2008 годы) по теме «Численное моделирование нестационарного взаимодействия сложных упругих конструкций с жидкостью или газом», блок 2: «Нестационарное взаимодействие нелинейных поверхностных волн с плавающими и закрепленными упругими конструкциями", пункт 1: «Развитие методов расчета гидродинамических нагрузок при резко нестационарном воздействии волн с большими деформациями области течения».
Проект № 4256 «Создание типового информационно-вычислительного портала для организации учебной и научной деятельности вуза» по ведомственной научной целевой программе федерального агентства по образованию «Развитие научного потенциала высшей школы (2006 - 2008 годы)» (2006 -2008 годы).
Проект 1.3368.2011 Государственное задание на проведение научно-исследовательских работ (2011 -2014 годы).
Литература
1. Информационные и вычислительные технологии в численных расчетах и управлении вузом: сборник трудов научной школы «Информационно-образовательные технологии» / под общ. ред. проф. К. Е. Афанасьева.
- Кемерово: ИНТ, 2010. - 604 с.
2. Афанасьев, K. E. Направления научных исследований кафедры ЮНЕСКО по новым информационным технологиям / К. Е. Афанасьев // Вестник Кемеровского государственного университета. - 2004. - № 3(19). -С. 3 - 12.
3. Afanasiev, K. E. Simulation of problems with free surfaces by a boundary element method / K. E. Afanasiev, S. V. Stukolov // Springer: Computational Science and High Performance Computing. Series: Notes on Numerical Fluid Mechanics and Multidisciplinary Design (NNFM). - 2005. - V. 88. - P. 307 - 338.
4. Afanasiev, K. E. Numerical investigation of three-dimensional bubble dynamics / K. E. Afanasiev, I. V. Grigorieva // Journal of Engineering Mathematics. Springer. - 2006. - V. 55. - № 1 - 4. - P. 65 - 80.
5. Afanasiev, К. Е. Comparative analysis of the SPH and ISPH Methods Comparative analysis of the SPH and ISPH Methods / К. Е. Afanasiev, R. S Makarchuk., A. Yu. Popov // Computational Science and high Performance Computing III. The 3rd Russian - German Advanced Research Workshop. - Novosibirsk, Russia, 2007. - Page 206 -223.
6. Afanasiev, К. Е. Numerical simulation of the dam break problem by general natural element method / К. Е. Afanasiev, T. S. Rein // Proceedings of 24th International Workshop on Water Waves and Floating Bodies. St. Petersburg, Russia. April, 19 - 22. - 2009. - Page 5 - 9.
7. Afanasiev, K. Е. Numerical simulation of hydrodynamics problems with free surfaces by natural element method / К. Е. Afanasiev, S. N. Karabcev, T. S. Rein // Zbornik radova konferecije MIT 2009. - Serbia, Beograd. - Page 13 - 23.
8. Afanasiev, К. Е. Hydrodynamic Loads Computation Using the Smoothed Particle Methods / K. Е. Afanasiev, R. S. Makarchuk, A. Popov // Hydrodynamics - Optimizing Methods and Tools. - 2011. - Page 51 - 78.
9. Afanas'ev, K. E. Calculation of hydrodynamic loads at solid boundaries of the computation domain by the ISPH method in problems with free boundaries / K. E. Afanas'ev, R. S. Makarchuk // Russian Journal of Numerical Analysis and Mathematical Modelling. - 2011. - Vol. 26. - № 5. - Р. 447 - 464.
10. Afanas'ev, K. E. Application of the natural element method for solution of problems of fluid dynamics with free boundaries / K. E. Afanas'ev, Yu. I. Shokin // Russian Journal of Numerical Analysis and Mathematical Modelling.
- 2011. - Vol. 26. - № 5. - Р. 515 - 531.
11. Chorin, A. Numerical solution of the Navier-Stokes equations / A. Chorin // Math. Comp. - 1968. - Vol. 22. -Р. 745 - 762.
12. Crespo, J. S. Effect of wet bottom on dam break evolution / J. S. Crespo // SPH European research interest community SIG. - 2007. - № 6.
13. Del Pin, F. The meshless finite element method applied to a lagrangian particle formulation of fluid flows // Instituto de Desarrollo tecnologico para la industria quimica (INTEC) universidad nacional del litoral noviembre. - 2003.
- 157 p.
14. De Berg, M. Computational Geometry. Algorithms and Applications. Second, Revised Edition / M. De Berg, M. Van Kreveld. - Berlin: Springer-Verlag, 2000.
15. Edelsbrunner, H. Three-dimensional alpha shapes / H. Edelsbrunner, E. P. Macke // ACM Trans. Graph. -1994. - Vol. 13. - № 1. - P. 43 - 72.
16. Fortune, S. J. A sweepline algorithm for Voronoi diagrams / S. J. Fortune // J. Algorithmica. - 1987. - № 2. -P. 153 - 174.
17. Idelsohn, S. A Lagrangian Meshless Finite Element Method applied to Fluid Structure Interaction Problems / S. Idelsohn, E. Onate, F. Del Pin // J. Computer and Structures, 81, 655 - 671 (2003).
18. Greenhow, M. Non-linear free surface effects: Experiments and theory / M. Greenhow, W. M. Lin // Rep.
№. 83 - 19, Dept. of Ocean Engineering, MIT, Cambridge, MA. - 1983. - Режим доступа:
http://www.dtic.mil/cgibin/GetTRDoc?AD=ADA161079 (датаобращения: 10.08.2012).
19. Koshizuka, S. A particle method for incompressible viscous flow with fluid fragmentation / S. Koshizuka,
H. Tamako, Y. Oka // Computational Fluid Dynamics Journal. - 1995. - P. 29 - 46.
20. Li, S. Meshfree and particle methods and their applications / S. Li, W. K. Liu // Appl. Mech. Rev. - 2002. -№ 55. - P. 1 - 34.
21. Liu, M. B. Smoothed particle hydrodynamics for numerical simulation of underwater explosion / M. B. Liu
[etal.] // Comput. Mech. - 2003. - 30. - P. 106 - 118.
22. Liu, G. R. Smoothed particle hydrodynamics: a meshfree particle method / G. R. Liu, M. B. Liu // World
Scientific Publishing Company. - 2003. - 472 p.
23. Lucy, L. B. A numerical approach to the testing of the fission hypothesis / L. B. Lucy// Astron. J. - 1977. -82(12). - P. 1013 - 1024.
24. Melenk, J. M. The Partition of Unity Finite Element Method: Basic Theory and Applications / J. M. Melenk,
I. Babuska // Comp. Methods Appl. Mech. Engrg. - 1996. - V. 139. - P. 289 - 314.
25. Monaghan, J. Smoothed particle hydrodynamics / J. Monaghan// Ann. Rev. Astron and Astrophysics. - 1992. -№ 30. - P. 543 - 574.
26. Monaghan, J. J. Simulation of free surface flows with SPH / J. J. Monaghan, M. C. Thompson, K. Hourigan // J. Comput. Phys. - 1994. - 110. - P. 399 - 406.
27. Onate, E. A finite point method in computational mechanics. Applications to convective transport and fluid flow / E. Onate, S. R. Idelsohn, O. C. Zienkiewicz // International Journal for Numerical Methods in Engineering. -1996. - № 39. - P. 3839 - 3866.
28. Sibson, R. A vector identity for the Dirichlet Tessellation / R. Sibson // Math. Proc. Cambridge Phil. Soc. 87(1), 1980. - P. 151 - 155.
29. Seabra-Santos, F. J. Numerical and experimental study of the transformation of a solitary wave over a shelf or isolated obstacle / F. J. Seabra-Santos, D. P. Renouard, A. M. Temperville // J. Fluid Mech. - 1987. - Vol. 176. - P. 117
- 134.
30. Sukumar, N. The natural element method in solid mechanics / N. Sukumar, B. Moran, T. Belytschko // Int. J. Num. Methods Eng. - 1998. - № 43(5). - P. 839 - 887.
31. Sukumar, N. Natural Neighbor Galerkin Methods / N. Sukumar, B. Moran, A. Yu. Semenov [et al.] // International Journal for Numerical Methods in Engineering. - 2001. - Vol. 50. - № 1. - P. 1 - 27.
32. Traversoni, L. Natural neighbor finite elements / L. Traversoni // Computational Mechanics Publications. -1994. - № 2. - P. 291 - 297.
33. Zhu, X. Application of the CIP Method to Strongly Nonlinear Wave-Body Interaction Problems / X. Zhu // Doctoral thesis for the degree of doctor ingenior. - 2006. - Режим доступа: http://ntnu.diva-portal.org/smash/get/-diva2:122822/FULLTEXT01 (дата обращения: 10.08.2012).
34. Афанасьев, К. Е. КМГЭ для решения плоских задач гидродинамики и его реализация на параллельных компьютерах: учеб. пособие / К. Е. Афанасьев, С. В. Стуколов. - Кемерово: Изд-во КемГУ, 2001. - 206 с.
35. Афанасьев, К. Е. Анализ динамических характеристик при взаимодействии уединенной волны с препятствием / К. Е. Афанасьев, Е. Н. Березин // Вычислительные технологии. - 2004. - Т. 9.- № 3. - С. 22 - 37.
36. Афанасьев, К. Е. Решение модельных задач гидродинамики методом естественных соседей/ К. Е. Афанасьев, Т. С. Рейн // Труды VII Всероссийской научно-практической конференции «Инновационные Недра Кузбасса. IT-технологии». - Кемерово: ИНТ. - 2008. - С. 286 - 291.
37. Афанасьев, К. Е. Метод естественных соседей на основе интерполяции Сибсона / К. Е. Афанасьев, С. Н. Карабцев, Т. С. Рейн [и др.] // Вестник Томского государственного университета. Приложение. - 2006. -№ 19. - С. 210 - 219.
38. Афанасьев, К. Е. Метод естественных соседей для решения задач вязкой несжимаемой жидкости / К. Е. Афанасьев, Т. С. Рейн // Вестник Новосибирского государственного университета. - 2008. - Т. 8. - Вып. 2.
- С. 30 - 38. - (Серия: Математика, механика, информатика).
39. Афанасьев, К. Е. Алгоритм поиска ближайших соседей в методе сглаженных частиц и его параллельная реализация / К. Е. Афанасьев, Р. С. Макарчук, А. Ю. Попов // Вычислительные технологии. Специальный выпуск 5. - 2008. - Т. 13. - С.9 - 14.
40. Афанасьев, К. Е. Моделирование задач гидродинамики вязкой несжимаемой жидкости со свободными границами бессеточным методом естественных соседей / К. Е. Афанасьев, Т. С. Рейн // Вычислительные технологии. - 2008. - Т. 13. - № 4. - С. 7 - 24
41. Афанасьев, К. Е. Численное и экспериментальное исследование процессов обрушения и последующего распространения нелинейных удаленных волн в прибрежной зоне / К. Е. Афанасьев, Е. Ю. Камынин, С. Н. Ка-рабцев [и др.] // Труды X Всероссийской конференции «Прикладные технологии гидроакустики и гидрофизики». - СПб.: Наука, 2010. - С. 304 - 307.
42. Афанасьев, К. Е. Численное моделирование работы опытового волнопродуктора одиночных волн / К. Е. Афанасьев, В. В. Максимов, И. С. Нуднер [и др.] // Труды XI Всероссийской конференции «Прикладные технологии гидроакустики и гидрофизики». - СПб.: Наука. - 2012. - С. 201 - 203.
43. Афанасьев, К. Е. Моделирование задач методом граничных элементов / К. Е. Афанасьев, Е. Н. Березин // LAPLAMBERTAcademicPublishingGmbH&Co.KG. - 2012. - 98 с.
44. Беликов, В. В. Несибсоновская интерполяция - новый метод интерполяции значений функции на произвольной системе точек / В. В. Беликов, В. Д. Иванов, В. К. Конторович [и др.] // Вычислительная математика и математическая физика. - 1997. - Т. 37. - № 1. - С. 11 - 17.
45. Березин, Е. Н. Численное моделирование задач идеальной несжимаемой жидкости со свободными границами методом граничных элементов: автореф. дис. ... канд. физ.-мат. наук / Е. Н. Березин. - Кемерово, 2006. -26 с.
46. Григорьева, И. В. Несогласованная линейная аппроксимация в методе граничных элементов для решения пространственных задач / И. В. Григорьева, С. А. Томилов // Гидродинамика больших скоростей и численное моделирование. - Кемерово: ИНТ, 2006. - C. 351 - 357.
47. Коннор, Дж. Метод конечных элементов в механике жидкости / Дж. Коннор, К. Бреббия. - Л.: Судостроение, 1979. - 204 с.
48. Лойцянский, Л. Г. Механика жидкости и газа / Л. Г. Лойцянский. - М.: Наука: Гл. ред. Физ.-мат. лит., 1987. - 840 с.
49. Карабцев, С. Н. Метод естественных соседей для решения задач идеальной несжимаемой жидкости со свободными границами: автореф. дис. ... канд. физ.-мат. наук / С. Н. Карабцев. - Кемерово, 2008. - 24 с.
50. Карабцев, С. Н. Эффективный алгоритм генерации конечноэлементной сетки для метода естественных соседей / С. Н. Карабцев, С. В. Стуколов // Материалы III Международной научной летней школы «Гидродинамика больших скоростей и численное моделирование». - Кемерово, ИНТ, 2006. - С. 401 - 409.
51. Карабцев, С. Н. Применение метода естественных соседей к решению задач механики жидкости со свободными поверхностями / С. Н. Карабцев, Т. С. Рейн // Материалы III Международной научной летней школы «Гидродинамика больших скоростей и численное моделирование». - Кемерово: ИНТ, 2006. - С. 393 - 401.
52. Карабцев, С. Н. Численное моделирование задачи о взаимодействии уединенной волны с подводной ступенькой методом естественных соседей / С. Н. Карабцев, С. В. Стуколов // Вестник Новосибирского государственного университета. - 2008. - Т. 8. - № 2. - С. 120 - 127. - (Серия: Математика, механика, информатика).
53. Карабцев, С. Н. Построение диаграммы Вороного и определение границ области в методе естественных соседей / С. Н. Карабцев, С. В. Стуколов // Вычислительные технологии. - 2008. - Т. 13. - № 3. - С. 65 - 81.
54. Карабцев, С. Н. Параллельная реализация метода естественных соседей для решения задач гидродинамики со свободными границами / С. Н. Карабцев, Т. С. Рейн // Вычислительные технологии. - 2008. - Т. 13, Спец. выпуск № 5: Избранные доклады VI Международной научно-практической конференции «Информационные технологии и математическое моделирование». - С. 56 - 60.
55. Коробкин, А. А. Движение вертикальной стенки, закрепленной на пружинах, под действием поверхностны волн / А. А. Коробкин, С. В. Стуколов, И. В. Стурова // ПМТФ. 2009.- Т. 50.- № 5. - С. 132 - 142.
56. Макарчук, Р. С. Матаматическое моделирование течений вязкой несжимаемой жидкости со свободными границами методом сглаженных частиц: автореф. дис. ... канд. физ.-мат. наук / Р. С. Макарчук. - Кемерово, 2012. - 24 с.
57. Марчук, Г. И. Введение в проекционно-сеточные методы / Г. И. Марчук, В. И. Агошков.- М.: Наука, 1981. - 416 с.
58. Манойлин, С. В. Некоторые экспериментально-теоретические методы определения воздействия волн цунами на гидротехнические сооружения и акватории морских портов: препринт / С. В. Манойлин.- Красноярск: ВЦ СО АН СССР, 1989. - 50 с.
59. Марчук, Г. И. Применение метода расщепления (дробных шагов) для решения задач математической физики / Г. И. Марчук, Н. Н. Яненко // Некоторые вопросы прикл. и вычисл. матем. - Новосибирск, 1966. -C. 5 - 22.
60. Овсянников, Л. В. Общие уравнения и примеры / Л. В. Овсянников // Задача о неустановившемся движении жидкости со свободной границей. - Новосибирск: Наука, 1967.
61. Патанкар, С. Численные методы решения задач теплообмена и динамики жидкости / С. Патанкар. - М.: Энергоатомиздат, 1984. - 152 с.
62. Рейн, Т. С. Численное моделирование движения вязкой несжимаемой жидкости со свободными границами обобщенным методом естественных соседей: автореф. дис. ... канд. физ.-мат. наук / Т. С. Рейн. - Кемерово, 2008. - 28 с.
63. Свидетельство об официальной регистрации программы для ЭВМ «Библиотека параллельных вычислительных программ для задач гидродинамики со свободными границами» («HydroParaLib»)» / К. Е. Афанасьев, С. В. Стуколов, Е. Н. Березин [и др.]; правообладатель Государственное образовательное учреждение высшего профессионального образования «Кемеровский государственный университет». - №2008611207, выдано свидетельство 07.03.2008 г.
64. Темам, Р. Уравнения Навье-Стокса, теория и численный анализ / Р. Темам. - М.:Мир, 1981.
65. Терентьев, А. Г. Численные методы в гидродинамике: учебное пособие / А. Г. Терентьев, К. Е. Афанасьев // Чуваш.гос. ун-т. - Чебоксары, 1987. - 80 с.
66. Франк, А. М. Дискретные модели несжимаемой жидкости / А. М.Франк. - М.: Физматлит, 2001. - 208 с.
67. Фабер, Т. Е. Гидроаэродинамика / Т. Е. Фабер. - М.: Постмаркет, 2001. - 559 с.
68. Хажоян, М. Г. Численное моделирование взаимодействия поверхностных волн с подводными препятствиями / М. Г. Хажоян, Г. С. Хакимзянов // Вычислительные технологии. - 2003. - Т. 8. - № 4. - С. 108 - 123.
69. Хакимзянов, Г. С. Численное моделирование течений жидкости с поверхностными волнами / Г. С. Ха-кимзянов [и др.]. - Новосибирск: Изд-во СО РАН, 2001. - 393 с.
70. Яненко, Н. Н. Об экономичных неявных схемах (метод дробных шагов) / Н. Н. Яненко. - Докл. АН СССР. - 1960. - Т. 134. - № 5. - С. 1034 - 1036.
Информация об авторах:
Афанасьев Константин Евгеньевич - доктор физико-математических наук, профессор, проректор по информатизации КемГУ, 8(3842) 58-32-89, afa@kemsu.ru.
Konstantin E. Afanasiev - Doctor of Physics and Mathematics, Professor, Vice-RectorforInformatization, Kemerovo State University.
Карабцев Сергей Николаевич - кандидат физико-математических наук, доцент кафедры ЮНЕСКО по НИТ КемГУ, 8 (3842)54-64-69, scarab@kemsu.ru.
Sergey N. Karabtsev- Candidate of of Physics and Mathematics, Assistant Professor at the UNESCO Chair for New ¡information Technologies, Kemerovo State University.
Макарчук Роман Сергеевич - кандидат физико-математических наук, доцент кафедры ЮНЕСКО по НИТ КемГУ, 8 (3842) 54-32-47, mak@kemsu.ni.
Roman S. Makarchuk - Candidate of Physics and Mathematics, Assistant Professor at the UNESCO Chair for New Information Technologies, Kemerovo State University.
Рейн Татьяна Сергеевна - кандидат физико-математических наук, доцент кафедры ЮНЕСКО по НИТ КемГУ, 8 (3842)54-64-69, rein@kemsu.ru.
Tatiana S. Reyn- Candidate of Physics and Mathematics, Assistant Professor at the UNESCO Chair for New Information Technologies, Kemerovo State University.
