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

Вычислительные технологии

Том 3, № 6, 1998

НЕЛИНЕЙНАЯ ЗАДАЧА О ВОЛНАХ, ВОЗНИКАЮЩИХ НА ГРАНИЦЕ РАЗДЕЛА СРЕД, ПРИ ОДНОВРЕМЕННЫХ РАЗГОННЫХ И КОЛЕБАТЕЛЬНЫХ ДВИЖЕНИЯХ КРУГОВОГО ЦИЛИНДРА *

С. И. Горлов Омский филиал Института математики им. С. Л. Соболева СО РАН, Россия e-mail: gorlov@iitam.omsk.net.ru

Method for the solution of non-linear non-stationary problems of the circular cylinder motion near the interface of two heavy fluids is developed and numerically realized. The initial boundary value problem is reduced to a system of singular integral equations concerning singularities simulating the liquid and rigid boundaries and a function describing the shape of the fluids interface. The problem of vertical and horizontal oscillations of the circular cylinder accelerating from rest under a free surface of heavy fluid is investigated in detail. The calculation results of the generated wave structures as well as of the disrtibuted and total hydrodynamic characteristics of the outline are shown. The purely oscillatory motions are also considered. For some cases the obtained results are compared with the linear theory.

1. Введение

Задача о колебательном и разгонном движении тела под границей раздела двух жидких сред является традиционной в нестационарной волновой гидродинамике. Наибольшие успехи в этой области достигнуты при помощи линейной теории [3, 4]. Результаты, полученные при решении задачи в полной нелинейной постановке, немногочисленны и содержатся в обзоре [5]. Из результатов, не вошедших в данный обзор, нужно указать следующие работы. В [14] рассмотрена плоская нелинейная нестационарная задача о колебаниях контура под свободной поверхностью. Построен итерационный процесс по времени, на каждом шаге которого решается краевая задача для потенциала скорости. Использованы методы конечных и граничных элементов. Приведены результаты по расчету сил, действующих на тело, и профилей генерируемых волн. Колебания полностью погруженного кругового цилиндра под свободной поверхностью жидкости конечной глубины рассмотрены в [12]. С помощью интегральной теоремы Коши о вычетах задача сведена к сингулярному интегральному уравнению. Интегрирование по времени проведено с использованием

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, грант №96-01-00093.

© С. И. Горлов, 1998.

рядов Тейлора с учетом членов пятого порядка. Представлены результаты численного эксперимента. При помощи разработанного метода решена также задача о поступательном движении кругового цилиндра по окружности. В [10] исследованы гармонические колебания цилиндра под свободной поверхностью. Метод использует снос граничных условий, возникающих на свободной поверхности, на горизонтальную плоскость путем разложения в ряды Тейлора. На каждом шаге по времени задача относительно потенциала скорости решена при помощи метода граничных элементов в сочетании с быстрым преобразованием Фурье. Для шага по времени использован метод Рунге — Кутта четвертого порядка. Приведены результаты расчетов волновых профилей. Показана применимость разработанного метода для решения целого класса задач. Более сложная модель отрывного обтекания цилиндра, совершающего колебания под свободной поверхностью, рассмотрена в [15]. Использован метод случайных вихрей. Представлены результаты по расчету возбуждаемых волн. Получено хорошее соответствие с результатами лабораторных измерений.

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

2. Постановка начально-краевой задачи

Рассмотрим разгонное и колебательное движения кругового цилиндра L0(t) под границей раздела сред D1 и D2 (контур L0(t) расположен в области Di). Жидкость в каждом слое Dk является идеальной, несжимаемой, весомой и однородной с плотностью pk (k = 1, 2). Ось x системы координат располагается вдоль невозмущенной в начальный момент времени границы раздела L1(t).

Центр кругового цилиндра имеет координаты

Xc(t) = -—t2S 1 + — cosM)(1 - ¿2), 0 < t < T, 2t Ш

Uot - — T ) ¿1 + — cos(wt)(1 - ¿2), t > T, 2 ) Ш

yc(t) = -h - — cos(wt)¿2, Ш

где ¿1 = 1 соответствует разгонному движению, ¿1 = 0 — его отсутствию, ¿2 = 1 — вертикальным колебаниям, ¿2 = 0 — горизонтальным, ш — частота колебаний, T — время разгона, h — ордината среднего положения центра цилиндра при колебаниях.

Координаты точек, принадлежащих контуру L0(t), определяются формулами

x(e,t)= xc (t) - R cos в, у(в, t) = yc(t) + R sin в, в е [0°, 360°].

Здесь R — радиус цилиндра, угол в отсчитывается по часовой стрелке.

Будем рассматривать задачу относительно потенциала скорости ^k(x,y, t), удовлетворяющего в областях Dk уравнению Лапласа:

xc(t) = -

A^k(x,у, t) = 0, (x,у) е Dk(t)\Lo(t), k = 1, 2.

(1)

На границе раздела сред Li(t) выполняются кинематическое и динамическое условия

Öfi

V<£i(x,y, t) ■ ni = V^(x,y, t) ■ ni = — ■ ni, (x,y) G Li(t), (2)

Pi(x,y,t) = P2(x,y,t), (x,y) G Li(t), (3)

где ni — нормаль к Li(t) в точке (x,y), ri — радиус-вектор точки (x,y), pk — гидродинамическое давление в области Dk.

В точках контура Lo(t) выполняется условие непротекания

(V^i(x,y,t) - VLo(t)) ■ no = 0, (x,y) G Lo(t), (4)

где n0 — нормаль к Lo(t), VLo(t) = (Уь0Х,Уь0У) = (¿с,Ус) — вектор скорости движения контура в выбранной системе координат.

В бесконечно удаленных точках областей Di, D2 выполняются условия отсутствия возмущений скоростей и границы раздела сред:

lim V^fc(x, y, t) = 0, k = 1, 2, (5)

lim y = 0, (x,y) G Li(t). (6)

Начальные условия для границы раздела и потенциала имеют вид

y = 0, (x,y) G Li(0),

V^(x,y, 0) = 0, (x,y) G Dk(0)\Lo(0). (7)

Гидродинамическое давление в области Dk определяется при помощи интеграла Ко-ши — Лагранжа. Запишем этот интеграл в точке границы раздела Li(t), перемещающейся со скоростью V^o(x,y,t):

d^k (x,y,t) (V^k (x,y,t))2 ^ , +\Y7 i -dt-+-2--V^o(x,y,t)V^k (x,y,t) +

. Pk (x,y,t) , , Pk^(t) , s _ т +-+ gy(x,t) =-, (x,y) G Li(t) (8)

Pk Pk

где Pkrc>(t) — давление в бесконечно удаленной точке области Dk (k = 1, 2), ö/öt — обозначение для производной, вычисляемой в подвижной системе координат, g — ускорение свободного падения.

Динамическое условие (3) с учетом (8) можно записать в виде

d^i(x,y,t) d^(x,y,t) , (V^i(x,y,t))2 (V^2(¿,У,t))2 Pi-Kl--P2-TT7--Г Pi----P2

öt öt 2 2 -PiV^o(x, У, t) V^i(x, y, t) + P2V^o(x, У, t) V^2(x, y, t)+

+ (Pi - P2)gy(x,t) = 0, (x,y) G Li(t). (9)

3. Вывод системы интегро-дифференциальных уравнений

Сведем начально-краевую задачу относительно потенциалов скоростей ^ (x, y, t) к системе интегро-дифференциальных уравнений относительно интенсивностей особенностей, моделирующих жидкие и твердые границы. С этой целью рассмотрим интенсивность вихревого слоя Yi(si, t), расположенного вдоль контура Li(t), и слоя источников q(so,t) на контуре Lo(t). Предположим, что Yi(±^,t) = 0. Возмущенное движение жидкости в областях Dk (k = 1, 2) будет описываться функцией

V( t) = [ Yi(si,t) dsi + f g(so,t) dso ( )

V (Z,t) = 2n^ z - Z(si) +2W z - Z(so) , ( )

Li(t) Lo(t)

которая удовлетворяет условию отсутствия возмущений в бесконечно удаленных точках (5):

lim V(z,t) = 0. (11)

Предельные значения функции V(z,t) при подходе к z(si) G Li(t) сверху и снизу, а также к z(so) G Lo (t) с внешней стороны контура:

V+(z(si),t) = Vi(z(si),t) - 1 Yi(si,t)e-iÖ1 (si'4), (12)

V-(z(si),t)= V i(z(si),t) + 1 Yi(si,t)e-iöl(si'i), (13)

V + (z(so),t) = Vo(z(so),t) - 2q(so,t)e-iöo(so'i), (14)

где Vj(z(sj),t) определяется из (10) для z(sj) G Lj(t), 6j(sj, t) — угол между касательной в точке z(sj) G Lj(t) и осью x (j = 0,1). При z(sj) G Lj(t) (j = 0,1) несобственные интегралы, входящие в выражение (10) для комплексной скорости, следует понимать в смысле главного значения по Коши.

Из (12), (13) следует, что нормальная компонента скорости при переходе через границу раздела сред Li(t) непрерывна.

Примем следующее предположение относительно скорости V^o(x,y,t) ((x,y) G Li(t)):

dx--* dy = V i(z(Si),t)- (15)

С учетом (10), (12)-(15) граничные условия (2), (4), (9) будут иметь вид

dz (si)

dt

= Vi(z(si),t), z(si) G Li(t), (16)

dG(si,t) /Vi(z(si),t)Vi(z(si),t) + T ( ) Yi(si,t)'. ( L (t) (17) P*-ö--+ g Tmz(si)---- , z(si) G Li(t), (17)

dt 2 ' * v ^ 8

G(si,t)= / fY^ + 2p,Vis(ai,t^ dai, p, = ^,

J V 2 у Pi + Pi

—<x

,*) = Ив (Ц(;(в,),*)е*', ) е (*), 7 = 0,1, 1т ((Уо(;фо),*) - Цьо(¿))е^о'4)) , г(во) е По(*) = ^х - ¿У^у• (18)

Таким образом, начально-краевая задача (1)-(7) свелась к определению функций 71(^1, ¿) и д(зо,£), а также границы раздела из системы сингулярных интегро-дифференци-

альных уравнений (16)-(18) с учетом (10), (11) и начальных условий:

1т;(в1) = 0, ;(в1) е ¿1(0), 71(31, 0) = фо, 0) = 0.

(19)

Гидродинамическое давление в точках контура ; (зо) определяется при помощи интеграла Коши — Лагранжа, записанного в подвижной системе координат, связанной с кон-

туром:

- /= -Р1

д_

т

sо

И)« (ао, ¿) ^ао -

- Ив (V¿0(¿)Уо(;Ы,*)) +

V о (;ЫЖо (;Ы,*)

(20)

Здесь /(¿) — некоторая функция, зависящая только от времени.

Суммарные гидродинамические нагрузки Дх и можно получить интегрированием давления р(зо,£) вдоль контура:

(р(во,*) - /(¿))в-г0о(«о'4)

о.

(21)

2

о

2

4. Решение системы интегро-дифференциальных уравнений

Полученная система интегро-дифференциальных соотношений (16)—(18) нелинейна. Эта нелинейность обусловлена двумя факторами: интенсивности особенностей 71(в1, ¿) и д(зо, ¿) входят в граничные условия нелинейным образом и неизвестна граница раздела Ь1(^). Подобное обстоятельство вносит определенные трудности при решении системы (16)—(18).

Для решения системы интегро-дифференциальных уравнений (16)—(18) будем использовать метод коллокаций. Для этого на каждом шаге по времени ¿п (п = 0,1,...) будем рассматривать границу раздела ¿П на конечном интервале (верхний индекс служит для обозначения значения функции на п-м шаге по времени). Разобьем контуры ¿П на интервалы 1, в" »] (г = 1, ... , I) и ¿П на [во1, ] (7 = 1, ... , 7). Параметризация контура длиной дуги также зависит от времени. Рассмотрим на полученных интервалах точки кол-локации ;га(зП*) е Ь" К,* е [в?^, в"»]) и ;га(5о ,,) е ¿П (4, е [во,,-1, во,]). Будем требовать выполнения динамического и кинематического условий на границе раздела (16), (17) в точках (вП**) (г = 1, ... , I), условия непротекания на контуре £о(£) (18) в точках ,) (7 = 1, ... , 7). Решение полученной системы строится с использованием двух итерационных процессов, один из которых связан с интегрированием по времени уравнений (16) и (17) при помощи метода Рунге — Кутта — Фельберга пятого порядка точности [6]. При этом на каждом шаге по времени ¿п (п = 1, 2, ...) получаем значение функции Сга(зП*) и точки границы раздела ;га(зП*) е Вторая итерационная процедура связана с методом

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

+ = , Л*м) € ¿П, г =1, ...,/, (22)

= 1т ((К)) - ^)е^0 )) , ) € ^ = 1, ..., ^ (23)

Дискретизация (22), (23) основана на использовании метода панелей высокого порядка [9]. Граница раздела ¿П на г-м интервале [з^^ь^] (г = 1, ..., I) и контур ¿П на ]-м интервале [з0а-1, 50^] (^ = 1, ..., /) аппроксимируются параболой, а функции и 5га(з0) на этих же интервалах — линейной функцией. После решения системы линейных алгебраических уравнений относительно значений и 5га(з0) на концах интервалов, из (10) определяются значения комплексной скорости V (г) в точках коллокации ^) € ¿П, а из (20), (21) — распределенные и суммарные гидродинамические характеристики контура.

5. Результаты численного эксперимента

На основании разработанного метода проведен численый эксперимент по решению задачи о вертикальных и горизонтальных колебаниях кругового цилиндра, совершающего разгон, под свободной поверхностью тяжелой жидкости (р* = 1). Рассмотрен также случай чисто колебательных движений.

Во всех рассмотренных случаях расчетная область рассматривалась в интервале |ж/Д| < 25. Число узлов на свободной поверхности и контуре выбиралось равным 500 и 80 соответственно. С целью исключения отражения волн от границ вычислительной области в интервале 20 < |ж/Д| < 25 введен демпфирующий слой по методике, предложенной в [7]. Шаг интегрирования менялся динамически в процессе счета от Дт = 0, 05 до 0,005 (т = — безразмерное время). Устранение коротковолновой неустойчивости, возника-

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

в выражение (21) для суммарных гидродинамических нагрузок, вычислялись при помощи кубических сплайнов [2]. Разгон контура приводит к обрушению генерируемых волн. Для корректного моделирования этого процесса с момента появления неоднозначности на свободной поверхности и до момента полного обрушения применялось переразбиение, основанное на использовании параметрических сплайнов по технологии, описанной в [13].

Процесс численного решения задачи контролировался с помощью интегрального закона сохранения энергии [7]

~ +те

У ^о(Зь£)71 (51, г) ^1 +

Е

Р1 + Р2 4

+те +те

2

+ 2р* / (51,г) ^ + д у (в1,г)жв1 (зьг) ^

У*

-те

Рис. 1. Вертикальные (I) и горизонтальные (II) колебания кругового цилиндра при Ег = 0.5,

а = 2, Н/К = 2.

а — генерируемые поверхностные волны, б — суммарная гидродинамическая нагрузка Су (I) и Сх, Су (II), в — распределенная гидродинамическая нагрузка Ср для в = 0° (1), 90° (2), 270° (3), 270° (4) (— — настоящий метод, — — линейное дипольное приближение [5]).

где (х(в1, £), у(в1, £)) е ¿1 (£), £) и ^о(в1,£) — полусумма функций тока и потенциалов

при подходе к границе раздела ¿1(£) сверху и снизу. Использование перечисленных выше методик и процедур позволило проводить расчеты с изменением энергии, не превышающим 1 %.

Результаты численного решения задачи о вертикальных и горизонтальных колебаниях кругового цилиндра под свободной поверхностью жидкости представлены на рис. 1. Рассчитывались форма свободной поверхности, распределение давления по контуру СР = (р - /(¿))/(р1 У)2), суммарные гидродинамические характеристики контура (Сх, Су) = 2(Ях, Яу)/(р!Уо2Я). В качестве системы безразмерных параметров были приняты следующие величины: Ег=Уо/\/дД, а = ^Я/Ц, Л./Д.

На рис. 1, II, а показаны также волновые профили, полученные при решении линейной задачи о горизонтальных колебаниях кругового цилиндра в дипольном приближении. Безразмерное возвышение свободной поверхности определяется выражением [5]

п(£,т ) =

кйку Бт ат б1п о

соб ат

а

- £

соб (т — т) йт,

(24)

где (£, п) = (х,у)/Д. Вычисление двойного интеграла в (24) проводилось на основе куба-турных формул, описанных в [1]. Следует отметить, что волновые профили, полученные на основе линейной и нелинейной теорий, хорошо согласуются между собой. Подобный вывод был сделан также в [5] на основании сравнения дипольного приближения (24) и численных расчетов, представленных в [8].

т

к

о

Рис. 2. Вертикальные (I) и горизонтальные (II) колебания кругового цилиндра, совершающего

разгон, при Рг = 0.5, а = 2, К/К = 2, И0/У0 = 1. а — генерируемые поверхностные волны, б — суммарная гидродинамическая нагрузка Сх и Су, в — распределенная гидродинамическая нагрузка Ср для в = 0° (1), 90° (2), 180° (3), 270° (4).

На рис. 2 приведены расчеты волновых профилей и гидродинамических нагрузок контура, совершающего одновременные разгонные и колебательные движения. Разгон осуществлялся до безразмерного момента времени тх = ТУ0/Л. Наблюдается интересная картина обрушения волн, вызванных разгоном тела. Настоящий метод позволяет рассчитывать этот процесс до момента близкого взаимодействия образовавшейся струи со свободной поверхностью.

6. Заключение

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

Список литературы

[1] Бахвалов Н.С., Жидков Н.П., Кобельков Г. М. Численные методы. Наука, М., 1987.

[2] РУСАКОВ С. В. Методы сплайн-функций в вычислительной гидродинамике. Изд-во Пермского гос. ун-та, Пермь, 1987.

[3] СТЕПАНЯНЦ Ю.А., СтУРОВА И. В., ТЕОДОРОВИЧ Э. В. Линейная теория генерации поверхностных и внутренних волн. В "Итоги науки и техники. Сер. МЖГ", ВИНИТИ, М., 21, 1987, 93-179.

[4] СтУРОВА И. В. Линейная теория генерации поверхностных и внутренних волн локальными возмущениями: Дис. ... докт. физ.-мат. наук. Новосибирск, 1994.

[5] СТУРОВА И. В. Численные расчеты в задачах генерации плоских поверхностных волн. Препринт №5, ВЦ СО АН СССР, Красноярск, 1990.

[6] Форсайт Дж., Малькольм М., МОУЛЕР К. Машинные методы математических вычислений. Мир, М., 1980.

[7] Baker G. R., Meiron D. I., Orszag S. A. Application of a generalized vortex method to nonlinear free-surface flows. Proc. of the 3rd Int. Conf. on Numerical Ship Hydrodinamics, Paris, 1981, 179-191.

[8] HAUSSLING H. J., Coleman R. M. Finite-difference computations using boundary-fitted coordinates for free-surface potential flows generated by submerged bodies. Proc. of the 2nd Int. Conf. on Numerical Ship Hydrodinamics, Berkeley, 1977, 221-233.

[9] HESS J. L. Higher-order numerical solution of the integral equation for the two-dimensional Neumann problem. Comput. Meth. Appl. Mech. and Engng., 2, No. 1, 1973, 1-15.

[10] Kim Y. J., Kim D. J., Hwang J. H. Calculation of nonlinear free-surface flows using two-dimensional numerical wave tank. Proc. of the 4th Int. Offshore and Polar Eng. Conf., Osaka, 3, 1994, 25-31.

[11] Longuet-Higgins M.S., Cokelet E. D. The deformation of steep surface waves on water. I. A numerical method of computations. Proc. Roy. Soc., London, A350, 1976, 1-26.

[12] Teles DA Silva A. F., Peregrine D. H. Nonlinear perturbations on a free surface induced by a submerged body: a boundary integral approach. Engng. Anal., 7, No. 4, 1990, 214-222.

[13] TSAI W. T., Yue D. K. P. Interactions between a free surface and a vortex sheet shed in the wake of a surface-piercing plate. J. Fluid Mech., 257, 1993, 691-721.

[14] Wu G. X., EATOCK Taylor R. Time stepping solutions of the two-dimensional nonlinear wave radiation problem. Ocean Engng., 22, No. 8, 1995, 785-798.

[15] Yeung R. W., Vaidhyanathan M. Hyghly separated flows near a free surface. Int. Conf. Hydrodyn, Wuxi, 1994, 118-128.

Поступила в редакцию 12 февраля 1998 г., в переработанном виде 13 мая 1998 г.