Чисельне моделювання гiдродинамiчної взаємодiї тiл, що рухаються в рiдинi

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

Повний опис

Збережено в:

Бібліографічні деталі
Дата:	2006
Автори:	Горбань, В.О., Маcюк, С.В.
Формат:	Стаття
Мова:	Ukrainian
Опубліковано:	Інститут гідромеханіки НАН України 2006
Онлайн доступ:	http://dspace.nbuv.gov.ua/handle/123456789/4761
Теги:	Додати тег Немає тегів, Будьте першим, хто поставить тег для цього запису!
Назва журналу:	Digital Library of Periodicals of National Academy of Sciences of Ukraine
Цитувати:	Чисельне моделювання гiдродинамiчної взаємодiї тiл, що рухаються в рiдинi / В.О. Горбань, С.В. Маcюк // Прикладна гідромеханіка. — 2006. — Т. 8, № 3. — С. 27-49. — Бібліогр.: 23 назв. — укр.

Репозитарії

Digital Library of Periodicals of National Academy of Sciences of Ukraine

id	irk-123456789-4761
record_format	dspace
spelling	irk-123456789-47612009-12-23T12:01:05Z Чисельне моделювання гiдродинамiчної взаємодiї тiл, що рухаються в рiдинi Горбань, В.О. Маcюк, С.В. В работе развивается метод граничных интегральных элементов для решения задач определения параметров потенциальных полей, генерированных при движении нескольких тел в безграничной жидкости или одного тела в ограниченной области. Построен численный алгоритм, использующий аппроксимацию поверхности трехмерного тела системой треугольных элементов. Проведена апробация алгоритма, сравнение с известными результатами и экспериментальными данными. Получены зависимости коэффициентов присоединенных масс для трехосного эллипсоида, движущегося в трапецевидном канале, от параметров эллипсоида и канала. Проведены расчеты гидродинамических сил и моментов, возникающих при движении двух эллипсоидов с различными скоростями. У роботi розвивається метод граничних iнтегральних елементiв для розв'язання задач визначення параметрiв потенцiйних полiв, генерованих при руховi декiлькох тiл у безмежнiй рiдинi або одного тiла в обмеженiй областi. Побудовано чисельний алгоритм, що використовує апроксимацiю поверхнi тривимiрного тiла системою трикутних елементiв. Проведена апробацiя алгоритму, порiвняння з вiдомими розв'язками та експериментальними даними. Одержанi залежностi коефiцiєнтiв приєднаних мас для тривiсного елiпсоїда, що рухається в трапецевидному каналi, вiд параметрiв елiпсоїда та каналу. Проведенi розрахунки гiдродинамiчних сил i моментiв, що виникають при руховi двох елiпсоїдiв з рiзними швидкостями. A method of boundary integral elements is developed to obtain parameters of the potential fields generated when either two bodies move in an unbounded region or a body moves in the region with a complex boundary. The numerical algorithm that uses approximation of 3-D body surface with a system of triangular elements is constructed. Approbation of the algorithm and comparison of the present results with known numerical and experimental data are carried out. Added mass coefficients of triaxial ellipsoid that moves in the trapezoidal channel are obtained against different parameters of both the ellipsoid and the channel. Hydrodynamics forces and moments acting on the system of two ellipsoids moving with different velocities are calculated. 2006 Article Чисельне моделювання гiдродинамiчної взаємодiї тiл, що рухаються в рiдинi / В.О. Горбань, С.В. Маcюк // Прикладна гідромеханіка. — 2006. — Т. 8, № 3. — С. 27-49. — Бібліогр.: 23 назв. — укр. 1561-9087 http://dspace.nbuv.gov.ua/handle/123456789/4761 uk Інститут гідромеханіки НАН України
institution	Digital Library of Periodicals of National Academy of Sciences of Ukraine
collection	DSpace DC
language	Ukrainian
description	В работе развивается метод граничных интегральных элементов для решения задач определения параметров потенциальных полей, генерированных при движении нескольких тел в безграничной жидкости или одного тела в ограниченной области. Построен численный алгоритм, использующий аппроксимацию поверхности трехмерного тела системой треугольных элементов. Проведена апробация алгоритма, сравнение с известными результатами и экспериментальными данными. Получены зависимости коэффициентов присоединенных масс для трехосного эллипсоида, движущегося в трапецевидном канале, от параметров эллипсоида и канала. Проведены расчеты гидродинамических сил и моментов, возникающих при движении двух эллипсоидов с различными скоростями.
format	Article
author	Горбань, В.О. Маcюк, С.В.
spellingShingle	Горбань, В.О. Маcюк, С.В. Чисельне моделювання гiдродинамiчної взаємодiї тiл, що рухаються в рiдинi
author_facet	Горбань, В.О. Маcюк, С.В.
author_sort	Горбань, В.О.
title	Чисельне моделювання гiдродинамiчної взаємодiї тiл, що рухаються в рiдинi
title_short	Чисельне моделювання гiдродинамiчної взаємодiї тiл, що рухаються в рiдинi
title_full	Чисельне моделювання гiдродинамiчної взаємодiї тiл, що рухаються в рiдинi
title_fullStr	Чисельне моделювання гiдродинамiчної взаємодiї тiл, що рухаються в рiдинi
title_full_unstemmed	Чисельне моделювання гiдродинамiчної взаємодiї тiл, що рухаються в рiдинi
title_sort	чисельне моделювання гiдродинамiчної взаємодiї тiл, що рухаються в рiдинi
publisher	Інститут гідромеханіки НАН України
publishDate	2006
url	http://dspace.nbuv.gov.ua/handle/123456789/4761
citation_txt	Чисельне моделювання гiдродинамiчної взаємодiї тiл, що рухаються в рiдинi / В.О. Горбань, С.В. Маcюк // Прикладна гідромеханіка. — 2006. — Т. 8, № 3. — С. 27-49. — Бібліогр.: 23 назв. — укр.
work_keys_str_mv	AT gorbanʹvo čiselʹnemodelûvannâgidrodinamičnoívzaêmodiítilŝoruhaûtʹsâvridini AT macûksv čiselʹnemodelûvannâgidrodinamičnoívzaêmodiítilŝoruhaûtʹsâvridini
first_indexed	2025-07-02T07:58:20Z
last_indexed	2025-07-02T07:58:20Z
_version_	1836521201336320000
fulltext	ISSN 1561 -9087 Прикладна гiдромеханiка. 2006. Том 8, N 3. С. 27 – 49 УДК 539.59 ЧИСЕЛЬНЕ МОДЕЛЮВАННЯ ГIДРОДИНАМIЧНОЇ ВЗАЄМОДIЇ ТIЛ, ЩО РУХАЮТЬСЯ В РIДИНI В. О. Г О РБ А Н Ь, С. В. М АС ЮК Iнститут гiдромеханiки НАН України, Київ Отримано 04.03.2006 У роботi розвивається метод граничних iнтегральних елементiв для розв’язання задач визначення параметрiв по- тенцiйних полiв, генерованих при руховi декiлькох тiл у безмежнiй рiдинi або одного тiла в обмеженiй областi. Побудовано чисельний алгоритм, що використовує апроксимацiю поверхнi тривимiрного тiла системою трикутних елементiв. Проведена апробацiя алгоритму, порiвняння з вiдомими розв’язками та експериментальними даними. Одержанi залежностi коефiцiєнтiв приєднаних мас для тривiсного елiпсоїда, що рухається в трапецевидному кана- лi, вiд параметрiв елiпсоїда та каналу. Проведенi розрахунки гiдродинамiчних сил i моментiв, що виникають при руховi двох елiпсоїдiв з рiзними швидкостями. В работе развивается метод граничных интегральных элементов для решения задач определения параметров потен- циальных полей, генерированных при движении нескольких тел в безграничной жидкости или одного тела в огра- ниченной области. Построен численный алгоритм, использующий аппроксимацию поверхности трехмерного тела системой треугольных элементов. Проведена апробация алгоритма, сравнение с известными результатами и экспе- риментальными данными. Получены зависимости коэффициентов присоединенных масс для трехосного эллипсоида, движущегося в трапецевидном канале, от параметров эллипсоида и канала. Проведены расчеты гидродинамических сил и моментов, возникающих при движении двух эллипсоидов с различными скоростями. A method of boundary integral elements is developed to obtain parameters of the potential fields generated when either two bodies move in an unbounded region or a body moves in the region with a complex boundary. The numerical algorithm that uses approximation of 3-D body surface with a system of triangular elements is constructed. Approbation of the algorithm and comparison of the present results with known numerical and experimental data are carried out. Added mass coefficients of triaxial ellipsoid that moves in the trapezoidal channel are obtained against different parameters of both the ellipsoid and the channel. Hydrodynamics forces and moments acting on the system of two ellipsoids moving with different velocities are calculated. ВСТУП Вивчення гiдродинамiчної взаємодiї тiл, що ру- хаються в рiдинi, а також взаємодiї тiл i гра- ниць течiї має важливе практичне значення для вирiшення багатьох задач гiдромеханiки судна. До такого класу задач належать задачi про рух судна в обмеженому фарватерi (зокрема, на мiлкiй водi, в каналi, бiля берега), про спiльний рух i гiдро- динамiчну взаємодiю декiлькох суден, про вплив рухомого судна на берeговi споруди тощо. Акту- альнiсть цих задач зростає з пiдвищенням iнтен- сивностi судноплавства. Зокрема, з початку ми- нулого столiття в морських акваторiях та рiчках зафiксовано велику кiлькiсть зiткнень суден при обгонах i зустрiчах. Слiд зазначити, що при ру- ховi суден у каналах виникає явище притягання до стiнок каналу, викликане несиметричнiстю по- ля тиску з правого i лiвого бортiв судна при його змiщеннi вiд площини симетрiї та наближеннi до одного з вiдкосiв каналу. При руховi суден на мiл- кiй водi виникає явище притягання судна до дна, спостерiгається збiльшення середньої просадки су- дна (тобто його занурення зростає вiдносно стану спокою). Тому окремi частини судна заглиблюю- ться значно бiльше, нiж їхня просадка. Часто це зумовлює просiдання судна на грунт. Останнiм часом все бiльшу роль вiдiграють за- дачi впливу судноплавства в рiчках на навколи- шнє середовище, зокрема, на екологiчну ситуацiю в прибережнiй зонi, на процеси перенесення шкi- дливих домiшок, на змiну рiвня води бiля бере- га тощо. Важливими для оптимiзацiї режимiв ру- ху суден є також задачi взаємодiї судна з нерiв- ностями дна, наприклад, з локальним пiдйомом дна (зменшенням глибини) або локальними змiна- ми перерiзу каналу (звуження, розширення, змi- на форми). Як свiдчить досвiд, локальнi неодно- рiдностi (береговi споруди, пришвартованi судна та iншi) приводять до трансформацiї гiдродинамi- чних полiв, викликаних рухом судна, зумовлюють виникнення додаткових гiдродинамiчних сил як на береговi споруди, так i на судно. Це пояснює не- абиякий iнтерес iнженерiв–суднобудiвникiв та гi- дротехнiкiв до результатiв вiдповiдних гiдродина- мiчних дослiджень. Гiдродинамiчна взаємодiя тiл, що рухаються в рiдинi, розглядалaся у роботах Г. Г. Павленка, А. А. Костюкова, Ю. М. Мастушкiна, Ньюмена та iнших дослiдникiв [6-8, 10]. Для аналiзу авто- ри використовували експериментальнi методи або c© В. О. Горбань, С. В. Маcюк, 2006 27 ISSN 1561 -9087 Прикладна гiдромеханiка. 2006. Том 8, N 3. С. 27 – 49 спрощенi напiвемпiричнi пiдходи [16-18, 21]. Новий напрямок дослiджень пов’язаний з розвитком ме- тодiв чисельного моделювання [3, 12-16, 23]. В бiльшостi робiт стосовно взаємодiї суден роз- глядалися два судна, що рухаються паралельно з постiйними швидкостями V1 i V2, зокрема випадки V1 = V2 (моделювання операцiй перевантаження чи дозаправки пiд час руху суден у морi), V1 = −V2 (при зустрiчному розходженнi), V1 ≈ V2 (при обгонi суден i при маневруваннi) i V1 = 0, V2 6= 0 (при проходженнi одного судна повз пришварто- ване судно). Зокрема, в роботi Коена i Бека [12] викладено результати дослiджень чотирьох рiзних ситуацiй: задачi про поворот судна на певний кут, про рух судна бiля вертикальної стiнки, про рух судна в каналi з двома вертикальними стiнками та про рух судна повз iнше нерухоме судно. Всi експе- рименти проводились з врахуванням дна, змiню- валися глибина води, вiдстань вiд судна до стiнки або вiдстань мiж двома суднами. В експериментах вимiрювалися гiдродинамiчнi сили, що дiють на моделi суден. Викладенi також чисельнi результа- ти, одержанi з застосуванням наближення “тонко- го” тiла. Проведене в роботi порiвняння чисельних результатiв з експериментом свiдчить про їх повну якiсну кореляцiю, хоча кiлькiсно результати були близькими лише в окремих випадках. Експеримен- ти виявили також деякi цiкавi ефекти. Наприклад, показали, що боковi гiдродинамiчнi сили i момен- ти значно зростають при зменшеннi глибини води, а також при наближеннi судна до однiєї зi стiнок каналу. Для гiдродинамiчної взаємодiї двох тiл ха- рактерним є те, що сили, якi дiють на рухоме су- дно, значно меншi вiд сил, якi дiють на нерухоме судно. Це свiдчить, зокрема, про те, що при руховi судна в каналах чи в iнших обмежених акваторi- ях бiльшу увагу слiд придiляти гiдродинамiчним навантаженням на береговi споруди та пришвар- тованi або заякоренi судна. В роботi [15] на основi теорiї “тонкого” тiла одер- жанi оцiнки для додаткових гiдродинамiчних бо- кових сил i моментiв, що дiють на судно, при про- ходженнi повз рiзнi нерiвномiрностi в каналi (зву- ження каналу, виступи) i показано, що цi сили i моменти можуть бути значними. Гiдродинамiчнi явища помiтно вiдрiзняються залежно вiд швидкостi судна i його розмiрiв, а та- кож залежно вiд числа Фруда Fr. При малих чи- слах Фруда Fr << 1 хвилi практично не утворю- ються, вiдбувається денiвеляцiя (локальна змiна рiвня) поверхнi, зумовлена перерозподiлом тиску на вiльнiй поверхнi при обтiканнi тiла, виникають додатковi пiдсмоктуючi сили до дна i до стiнки каналу. При бiльших числах Фруда значну роль вiдiграє хвилеутворення, яке зумовлює рiзке зрос- тання сили опору руховi судна. Розв’язок задачi гi- дродинамiчної взаємодiї суден, як правило, прово- дився без врахування хвилеутворення. Вiдомi ли- ше поодинокi розрахунки гiдродинамiчних полiв з врахуванням хвиль. Це пов’язано з iстотними тру- днощами як методологiчного, так i обчислюваль- ного характеру. В статтi обмежимось аналiзом задач без враху- вання процесiв хвилеутворення. Методологiя дослiджень. Для аналiзу гiдро- динамiчної взаємодiї тiл та оцiнки сил i моментiв застосовувались рiзнi пiдходи i методи. В багатьох роботах використовували експериментальнi коефi- цiєнти для сил та моментiв, а також напiвемпiри- чнi формули для просадки судна. Часто викори- стовувався наближений метод Блоха–Гiневського [4, 6, 8]. Потенцiал швидкостей системи тiл, що ру- хаються в iдеальнiй безмежнiй рiдинi, визначався наближено: Φ = n ∑ k=1 6 ∑ i=1 vkiϕki. Тут vki – проекцiї на осi нерухомої системи коорди- нат вектора лiнiйної швидкостi ~vk = (vk1, vk2, vk3) i вектора кутової швидкостi ~ωk = (ωk4, ωk5, ωk6), ϕki – потенцiал швидкостей, що визначається по- ступальним рухом k-го тiла вздовж однiєї з осей координат з одиничною швидкiстю i = 1, 2, 3 або обертанням його навколо однiєї з осей з одиничною кутовою швидкiстю i = 4, 5, 6 за умови, що всi iншi тiла є нерухомими. Для визначення ϕki використо- вуються спiввiдношення: ϕ0 ki = ϕki + 6 ∑ r=1 n ∑ j 6=k ∆v (k) jr · ϕjr, де ϕ0 ki – потенцiал, що визначає рух k-го тiла в безмежнiй рiдинi при вiдсутностi iнших тiл; ∆v (k) jr – iндуктивна швидкiсть, викликана рухом j-го тi- ла в центрi k-го тiла. Як показано Ю.М. Масту- шкiним [8], такий пiдхiд дає непоганi результати, а отриманi залежностi якiсно близькi до експери- ментальних. Останнiм часом широкого розповсюдження на- були чисельнi методи, зокрема, панельний метод (метод граничних елементiв – МГЕ). Основа цього методу – теорема Грiна, в якiй потенцiал швидко- стi в кожнiй точцi рiдини визначається розподiлом гiдродинамiчних особливостей по граничнiй по- верхнi [1, 2]. В загальному випадку задача визна- чення потенцiалу ϕ приводиться до iнтегрального рiвняння вiдносно невiдомої iнтенсивностi джерел 28 В. О. Горбань, С. В. Маcюк ISSN 1561 -9087 Прикладна гiдромеханiка. 2006. Том 8, N 3. С. 27 – 49 (або диполiв), розподiлених по граничнiй поверх- нi. Такий пiдхiд мав головним чином теоретичний iнтерес, доки Хесс i Смiт [14] не розвинули панель- ний метод i не показали його ефективнiсть для тривимiрних тiл, що рухаються в необмеженiй рi- динi. Враховуючи, що областi вихрового слiду бi- ля поверхнi судна та пограничного шару мають незначнi розмiри, при аналiзi гiдродинамiчної вза- ємодiї тiл, як правило, вони не враховуються. Тодi задача приводиться до розв’язання тривимiрного рiвняння Лапласа для потенцiалу швидкостi. Основнi кроки панельного методу полягають у наступному: 1. Гiдродинамiчний потенцiал представлено роз- подiлом джерел невiдомої iнтенсивностi по грани- чнiй поверхнi (непрямий метод граничних елемен- тiв). 2. Поверхня тiла апроксимується великим чи- слом елементарних панелей (у Смiта i Хесса це плоскi чотирикутники). 3. Iнтенсивностi джерел i диполiв приймаються постiйними на кожнiй панелi i є невiдомими. 4. У випадку задачi Неймана нормальна похiдна потенцiалу швидкостей в центрi (геометричному) кожної панелi прирiвнюється до нормальної швид- костi в цiй точцi. У випадку задачi Дiрiхле без- посередньо в цих точках визначається потенцiал швидкостей. Рiвняння для невiдомих iнтенсивно- стей джерел (або потенцiалiв) розв’язуються стан- дартними методами лiнiйної алгебри. 5. Сили та моменти визначають, застосовуючи рiвняння Бернуллi або рiвняння Лагранжа 2-го ро- ду. Дослiдниками було докладено багато зусиль, щоб застосувати описаний метод до задач вiльної поверхнi, якi включають у себе рух плаваючих та занурених тiл. Використання панельних методiв для вiльної поверхнi пов’язане з труднощами че- рез складнiсть вiдповiдної функцiї Грiна та обме- женi можливостi обчислювальної технiки. Числовi апроксимацiї функцiї Грiна для деяких випадкiв руху тiла бiля вiльної поверхнi розвивались Нью- меном [19, 20]. Незважаючи на велику кiлькiсть натурних та експериментальних дослiджень, слiд зазначити, що структура течiї при обтiканнi берегових струк- тур та кораблiв у каналах вивчена недостатньо. Зокрема, це стосується нерiвномiрних режимiв ру- ху та маневрування судна, особливостей роботи рушiйно–стернового комплексу, впливу ширини та посадки судна, якi в моделi “тонкого” тiла опису- ються наближено. У данiй роботi детально розглянуто непрямий метод граничних елементiв, з використанням три- кутних панелей в зручнiй для алгоритмiзацiї фор- мi. Проведена апробацiя цього методу для задачi просторового руху сфери та елiпсоїда у необмеже- нiй рiдинi. За допомогою вказаного методу розв’я- зана задача про рух тривiсного елiпсоїда в каналi. Проведенi дослiдження структури течiї, яка ви- никає при обтiканнi елiпсоїда при малих числах Фруда та денiвеляцiї вiльної поверхнi. Проведено також аналiз гiдродинамiчної взаємодiї двох три- вiсних елiпсоїдiв. 1. ПОСТАНОВКА ЗАДАЧI Розглянемо рух тiла в iдеальнiй нестисливiй рi- динi. Позначимо область, у якiй рухається тiло, че- рез Ω, а її граничну поверхню – через S. У випадку стацiонарної безвихрової течiї цей рух можна опи- сати за допомогою потенцiалу ϕ(~x), де ~x = (x, y, z) ∈ Ω. Для нестисливостi рiдини потенцiал ϕ(~x) за- довiльняє рiвнянню Лапласа: ∆ϕ(~x) = 0. (1) Використаємо ортогональну декартову систему координат 0xyz, зв’язану з рухомим судном: го- ризонтальна площина 0xy спiвпадає з поверхнею ватерлiнiї судна, тобто з незбуреною поверхнею во- ди. Будемо розглядати випадок малих чисел Фру- да Fr << 1 (Fr < 0.1); вплив хвилеутворення на сили гiдродинамичної взаємодiї судна з границя- ми (наприклад, з стiнкою каналу) чи з iншими су- днами не будемо розглядати. Якщо поверхня води незбурена, для аналiзу полiв, пов’язаних з судном, можна скористатись дзеркальним вiдображенням зануреної частини судна вiдносно площини 0xy. Гi- дродинамiчнi поля, генерованi одержаним таким чином тiлом у безграничнiй рiдинi, еквiвалентнi аналогiчним полям, що викликанi рухом судна, яке плаває на поверхнi води. Проводиться також перехiд вiд форми судна зi складними обрисами до еквiвалентного елiпсоїда. Для вибору розмiрiв цього елiпсоїда можна використати, наприклад, пiдхiд, коли довжина судна, його ширина та пpо- садка дорiвнюють вiдповiдно 2a, 2b i c (a, b, c – пiв- осi елiпсоїдa). Iнший можливий алгоритм: шири- на судна дорiвнює 2b, пpосадка спiвпадає з c, а величина a визначається iз умови рiвностi об’ємiв пiдводної частини судна i половини елiпсоїда. За- значимо, що використання моделi еквiвалентного елiпсоїдa для чисельного алгоритму, який побудо- вано далi, не є принциповим обмеженням. Нижче розглядається гiдродинамiчна взаємодiя декiлькох тiл, що рухаються в нерухомiй рiдинi, або взає- В. О. Горбань, С. В. Маcюк 29 ISSN 1561 -9087 Прикладна гiдромеханiка. 2006. Том 8, N 3. С. 27 – 49 модiя тiл з гладкими непроникними границями. Чисельнi результати приведенi для випадку три- вiсних елiпсоїдiв. Закон руху тiла ~V (~x0, t) вважа- ється заданим: ∂ϕ ∂n ∣ ∣ ∣ ∣ S = ~V (~x0, t) · ~n(~x0), (2) ~x0 ∈ S, ~n(~x0) = [nx, ny, nz] – одинична зовнiшня нормаль до поверхнi S, ~x0 = (x0, y0, z0) ∈ S. Сингулярнi розв’язки рiвняння Лапласа. Фундаментальнi розв’язки задачi дають значення потенцiалу ϕ в довiльнiй точцi спостереження ~x, викликанi одиничним джерелом, яке знаходиться в точцi ~ξ = (ξ, η, ζ). У випадку необмеженої рiди- ни, коли ~a – довiльний вектор, маємо: ϕ(~x) = G(~x, ~ξ), G(~x, ~ξ) = 1 4πr(~x, ~ξ) , (3) r(~x, ~ξ) = √ (x− ξ)2 + (y − η)2 + (z − ζ)2, ∂ϕ(~x) ∂a = ∂G ∂a = Ga(~x, ~ξ) = ~V (~x)~a(~x), ~V (~x) = 1 4πr3 ( ξ − x, η − y, ζ − z ) , Ga(~x, ~ξ) = (ξ − x)ax + (η − y)ay + (ζ − z)az 4πr3 . Для побудови розв’язку рiвняння Лапласа при ру- ховi тiла бiля вiльної поверхнi рiдини використо- вується обмеження, що тiло рухається з постiйною швидкiстю v0 в напрямку осi x. Якщо система ко- ординат зв’язана з тiлом, то в такiй рухомiй систе- мi координат потенцiал швидкостi ϕ є функцiєю лише координат x, y, z точок простору i не зале- жить вiд часу. Введемо також нерухому в просторi систему ко- ординат x1, y1, z1, яка в момент часу t = 0 збiгає- ться з рухомою системою x, y, z, тодi: x1 = x+ v0t, y1 = y, z1 = z; ϕ1(x1, y1, z1, t) = ϕ(x1 − v0t, y1, z1), звiдки ∂ϕ1 ∂t = −v0 ∂ϕ ∂x , ∂2ϕ1 ∂t2 = v2 0 ∂2ϕ ∂x2 . Якщо розглядати замiсть руху тiла з швидкiстю v0 обтiкання його потоком зi швидкiстю −v0 на не- скiнченностi, то потенцiал швидкостi ϕ∗ такої течiї можна записати у виглядi ϕ∗(x, y, z) = −v0x+ ϕ(x, y, z), де ϕ – потенцiал швидкостей, викликаний прису- тнiстю тiла в потоцi. Запишемо кiнематичну частину граничної умо- ви на вiльнiй поверхнi рiдини для потенцiала швидкостей ϕ∗. Нехай рiвняння вiльної поверхнi буде z = f(x, y). Тодi: −v0 ∂ϕ ∂x ∂f ∂x + ∂ϕ ∂y ∂f ∂y − ∂ϕ ∂z = 0, при z = f(x, y). Динамiчна умова на вiльнiй по- верхнi випливає з iнтеграла Ейлера-Бернулi: z − v0 g ∂ϕ ∂x + + 1 2g [( ∂ϕ ∂x )2 + ( ∂ϕ ∂y )2 + ( ∂ϕ ∂z )2] = 0. Тут g – прискорення вiльного падiння. В загально- му випадку гранична умова є нелiнiйною i викону- ється на невiдомiй хвильовiй поверхнi z = f(x, y). У випадку лiнiйних хвиль маємо: ∂2ϕ ∂x2 + g v2 0 ∂ϕ ∂z = 0, z = 0 при обтiканнi тiла потоком рiдини або: ∂2ϕ1 ∂t2 + g ∂ϕ1 ∂z1 = 0, z1 = 0 в нерухомiй системi координат 0x1y1z1. На нескiнченностi виконується умова затухання хвиль, що генеруються тiлом. Для такого випадку розв’язок побудовано Хавелоком i Сретенським [6, 11]: G(~x, ~ξ) = 1 4π [ 1 r(~x, ~ξ) − 1 r′(~x, ~ξ) +G1(~x, ~ξ) ] , r′(~x, ~ξ) = √ (x− ξ)2 + (y − η)2 + (z + ζ)2, G1(~x, ~ξ) = ν π × × π ∫ −π ∞ ∫ 0 eλ(ζ+z+i(x−ξ) cos θ+i(y−η) sin θ) ν − iµ cos θ − λ cos2 θ dλdθ, де ν = g/v2 0 , µ→ 0. У випадку малих чисел Фруда розв’язок набуває вигляду [6]: G(~x, ~ξ) = 1 4π [ 1 r(~x, ~ξ) + 1 r′(~x, ~ξ) ] . (4) Чисельнi алгоритми для обчислення функцiї G1(~x, ~ξ) побудованi у роботах Ньюмена [19, 20]. 30 В. О. Горбань, С. В. Маcюк ISSN 1561 -9087 Прикладна гiдромеханiка. 2006. Том 8, N 3. С. 27 – 49 У випадку обмеженої глибини води (рух судна на мiлководдi), коли на днi виконується умова не- протiкання, функцiю Грiна G(~x, ~ξ) можна побу- дувати, застосовуючи метод дзеркальних вiдобра- жень: G(~x, ~ξ) = +∞ ∑ k=−∞ 1 4π [ 1 rk(~x, ~ξ, h) + 1 r′k(~x, ~ξ, h) ] , (5) rk(~x, ~ξ, h) = √ (x− ξ)2 + (y − η)2 + (z − ζ + 2kh)2, r′k(~x, ~ξ, h) = √ (x− ξ)2 + (y − η)2 + (z + ζ + 2kh)2, Ga(~x, ~ξ) = ∂G ∂a = = +∞ ∑ k=−∞ [ (ξ − x)ax + (η − y)ay 4πr3k(~x, ~ξ, h) + + (ζ − z − 2kh)az 4πr3k(~x, ~ξ, h) + (ξ − x)ax 4πr′3k(~x, ~ξ, h) + + (η − y)ay + (ζ + z + 2kh)az 4πr′3k(~x, ~ξ, h) ] . Тут h− глибина води. 2. ЧИСЕЛЬНА МОДЕЛЬ Чисельна модель передбачає застосування мето- ду граничних елементiв, описаного нижче. Поверхня тiла S моделюється розподiлом дже- рел iнтенсивностi q(~ξ). Реакцiя рiдини в певнiй точцi спостереження ~x0 на розподiлення джерел q(~ξ), тобто значення в цiй точцi потенцiалу ϕ(~x0) та значення швидкостi в довiльному напрямку ~a(~x0), визначається через вiдповiднi iнтеграли по поверхнi S: ϕ(~x) = ∫ S G(~x, ~ξ)q(~ξ)dS +C. Враховуючи умови на нескiнченностi: r → ∞, ϕ→ 0, C = 0. Аналогiчно для ∂ϕ/∂a маємо: ∂ϕ(~x) ∂a = ∫ S Ga(~x, ~ξ)q(~ξ)dS. Якщо точка ~x0 розташована на поверхнi S, отри- маємо: ϕ(~x0) = ∫ S G(~x0, ~ξ)q(~ξ)dS, ∂ϕ(~x0) ∂n = ∗ ∫ S Gn(~x0, ~ξ)q(~ξ)dS. (6) Рис. 1. Апроксимацiя поверхнi граничними елементами Тут ∗ ∫ S – позначає сингулярний iнтеграл з особливi- стю вGn(~x0, ~ξ) при ~x0 → ~ξ i визначається як голов- не значення iнтегралу Кошi з додатковим членом, зумовленим особливiстю. У коректно поставленiй задачi одна з фун- кцiй ϕ(~x0) або ∂ϕ(~x0)/∂n повинна бути вiдомою в кожнiй точцi границi S. Рiвняння (6) можна розглядати як систему двох iнтегральних рiв- нянь вiдносно єдиної невiдомої функцiї q(~ξ). Пiсля знаходження q(~ξ) значення ϕ(~x) та швидкостi ~V (~x) в довiльнiй точцi ~x ∈ Ω розрахувати нескладно. Дискретизацiя iнтегралiв. Для розв’язання рiвнянь (6) будемо використовувати чисельнi ал- горитми, а саме дискретизацiю поверхневих iнте- гралiв плоскими трикутними граничними елемен- тами, на кожному з яких, наприклад j-елементi, iнтенсивнiсть поверхневих джерел q(~ξj) постiйна. Якщо поверхню S апроксимувати N - граничними трикутниками, як показано на рис. 1, то можна записати дискретнi аналоги рiвнянь (6). Для потенцiала ϕ(~xi 0) i нормальної складової швидкостi Vn(~xi 0) на i-му граничному елементi будемо мати: ϕ(~xi 0) = N ∑ j=1 q(~ξj) ∫ ∆Sj G(~xi 0, ~ξj)dS, (7) Vn(~xi 0) = N ∑ j=1 q(~ξj) ∗ ∫ ∆Sj Gn(xi 0, ~ξj)dS, де xi 0 – координати середньої точки i-го гранично- го елемента, ∆Sj – площа j-го граничного елемен- та. Формування матриць системи. Вирази (7) визначають потенцiал i нормальну швидкiсть у се- реднiй точцi (центрi мас) граничного елемента з номером i, викликанi дiєю всiх джерел з iнтенсив- ностями q(~ξ). Перепишемо рiвняння (7) у бiльш В. О. Горбань, С. В. Маcюк 31 ISSN 1561 -9087 Прикладна гiдромеханiка. 2006. Том 8, N 3. С. 27 – 49 зручному виглядi: ϕi = ( ∫ ∆Sj GijdS ) ~qj , V i n = ( ∗ ∫ ∆Sj Gij n dS ) ~qj, (8) де ~qj – вектор стовпчик розмiрностi N , а члени в дужках – вектор-рядки вiдповiдних розмiрностей. Кожний елемент цих вектор-рядкiв виходить у ре- зультатi iнтегрування. Таким чином, першим еле- ментом першого вектор-рядка в спiввiдношеннi (8) буде: Gi1 = ∫ ∆S1 G ( ~xi 0, ~ξ1 ) dS. Обчислення цих промiжних iнтегралiв буде де- тально описано нижче, а зараз зазначимо, що оста- точна форма спiввiдношень (8) еквiвалентна на- ступнiй: ϕi = (Gij)~qj, V i n = (Gij n )~qj . Якщо аналогiчнi операцiї провести для всiх гра- ничних елементiв (i = 1, 2, ...N), то одержану в ре- зультатi повну систему рiвнянь для ~ϕ i ~Vn можна записати у виглядi: ~ϕ = GS~q, ~Vn = GS n~q, (9) де очевидно ~ϕ, ~Vn, ~q – N -вимiрнi вектори грани- чних значень, а GS i GS n – матрицi розмiром N×N . Знаходження дотичних швидкостей на граничнiй поверхнi. Єдина невiдома в рiвнян- нi (9) ~q може бути знайдена звичайними метода- ми матричної алгебри шляхом обертання матрицi розмiром N ×N . Далi, пiдставляючи ~q почергово в кожне з рiвнянь (9), ми можемо визначити зна- чення потенцiалу та нормальної швидкостi на всiй поверхнi S. Також за допомогою другого з рiвнянь (9) може бути вирахована величина швидкостi у довiльному напрямку ~k(~xi 0). Для знаходження дотичних швидкостей на гра- ничнiй поверхнi скористаємось рiвнянням V i τ = ( ~Gij τ ) ~qj , де V i τ – невiдома дотична швидкiсть у напрям- ку ~τ (~x0) = (τx, τy, τz) на довiльному елементi; ~qj − вiдомий вектор розподiлу щiльностi фiктив- них джерел на поверхнi S, знайдений з (9); Gij τ = ∫ ∆Sj Gτ (~xi 0, ~ξj)dS. Знаходження потенцiалу та швидкостi у внутрiшнiх точках областi. Значення потенцi- алу у внутрiшнiх точках ~xi областi Ω вираховую- ться за вiдомими значеннями ~q за допомогою рiв- нянь (9). Компоненти швидкостi в тих самих то- чках ~xi в певному напрямку ~k(~xi) також можуть бути знайденi з рiвнянь (9). Визначення матриць Gij n та Gij є ключовим мо- ментом панельного методу. Видiлимо три типи елементiв матриць Gij n та Gij. 1) ДiагональнiGjj n таGjj – точка ~x0, в якiй обчи- слюється швидкiсть вiд j-го дискретного елемен- та, розташована на цьому елементi. 2) Елементи матрицi близькi за величиною до дiагональних: точка ~x0 – розташована на сусiднiх з j-тим дискретних елементах (випадок ближнього поля), до таких елементiв вiднесемо також iншi дискретнi елементи, розташованi на малiй вiдстанi вiд j-го елементa. Такi ситуацiї зустрiчаються для тiл складної форми з локальними звуженнями або при дискретнiй апроксимацiї поверхонь декiлькох тiл, близько розташованих одне до одного. 3) Всi iншi елементи матриць Gij n та Gij – так званий випадок дальнього поля, коли точка ~x0, в якiй обчислюється швидкiсть, розташована доста- тньо далеко вiд j-го елемента поверхнi. Алгоритми обчислення iнтегралiв Gij n та Gij у цих випадках iстотно вiдрiзняються, зокрема, для дiагональних елементiв Gjj n маємо сингулярнi iнте- грали. Точнiсть визначення матриць Gij n i Gij та необхiдний для цього об’єм обчислень визначають ефективнiсть чисельного методу. Тому дослiдни- ки придiляють велику увагу процедурi приведен- ня двовимiрних iнтегралiв до одновимiрних. Для деяких з цих iнтегралiв вдається побудувати ана- лiтичнi вирази. Враховуючи важливiсть цього ал- горитму, перетворення iнтегралiв Gij n , Gij τ та Gij докладно описано нижче. Перетворення координат. Для знаходження матриць коефiцiєнтiв у рiвняннях (9) потрiбно ви- раховувати iнтеграли по плоских елементах. Оче- видно, що для зручностi обчислення цих iнтегра- лiв доцiльно ввести локальну систему координат, зв’язану з елементом, i застосувати геометричне перетворення координат. Розглянемо для початку перетворення коорди- нат у двовимiрному випадку. Нехай ми маємо три- кутний елемент A1A2A3 (рис. 2) та P (x, y). Потрi- бно виконати таке перетворення: P (x, y) → P (τ∗1 , τ ∗ 2 ) → P (τ1, τ2), щоб: A3(x3, y3) → A3(0, 0) → A3(0, 0), 32 В. О. Горбань, С. В. Маcюк ISSN 1561 -9087 Прикладна гiдромеханiка. 2006. Том 8, N 3. С. 27 – 49 A1(x1, y1) → A1(τ ∗ 1 , 0) → A1(1, 0), A2(x2, y2) → A2(0, τ ∗ 2 ) → A2(0, 1), a б с Рис. 2. Перетворення координат: а – на площинi; б – в просторовому випадку (тетраедр); с – в просторовому випадку для трикутника з довiльною орiєнтацiєю тобто потрiбно виконати таке перетворення систе- ми координат 0xy у неортогональну (косокутну) систему координат 0τ1τ2, щоб її осi були спрямова- нi по сторонах трикутника, а довжини цих сторiн у новiй системi координат були рiвними одиницi (рис. 2, а) x− x3 = τ∗1 cosα1 + τ∗2 cosα2 = = τ∗1 x1 − x3 l1 + τ∗2 x2 − x3 l2 y − y3 = τ∗1 sinα1 + τ∗2 sinα2 = = τ∗1 y1 − y3 l1 + τ∗2 y2 − y3 l2 . Якщо ввести позначення τi = τ∗i /li, i = 1, 2, то матимемо: x = τ1(x1 − x3) + τ2(x2 − x3) + x3, y = τ1(y1 − y3) + τ2(y2 − y3) + y3. Узагальнюючи отримане перетворення, для тривимiрного випадку одержимо (рис. 2, б): x = τ1(x1 − x4) + τ2(x2 − x4) + τ3(x3 − x4) + x4, y = τ1(y1 − y4) + τ2(y2 − y4) + τ3(y3 − y4) + y4, z = τ1(z1 − z4) + τ2(z2 − z4) + τ3(z3 − z4) + z4. Осi нової системи координат 0τ1τ2τ3 спрямованi по ребрах тетраедра A1A2A3A4. Розглянемо випадок, коли плоский (трикутний) елемент знаходиться в просторi (рис. 2, с). Можна показати, що формули для перетворення координат мають вигляд: x = τ1(x1 − x3) + τ2(x2 − x3) + x3, y = τ1(y1 − y3) + τ2(y2 − y3) + y3, z = τ1(z1 − z3) + τ2(z2 − z3) + z3, або у матричному виглядi:   x y z   =   x1 x2 x3 y1 y2 y3 z1 z2 z3     τ1 τ2 1 − τ1 − τ2   . (10) В результатi осi нової системи координат τ1 i τ2 будуть спрямованi по сторонах трикутника A3A1 та A3A2, а вiсь τ3 буде перпендикулярною до пло- щини трикутника A1A2A3 Обчислення iнтегралiв. Розглянемо iнтегра- ли, якi потрiбно обчислити для знаходження ма- триць коефiцiєнтiв у рiвняннях (9): 1) ∫ ∆Sj G(~xi 0, ~ξj)dS, В. О. Горбань, С. В. Маcюк 33 ISSN 1561 -9087 Прикладна гiдромеханiка. 2006. Том 8, N 3. С. 27 – 49 2) ∗ ∫ ∆Sj Gn(~xi 0, ~ξj)dS, 3) ∫ ∆Sj Gτ(~xi 0, ~ξj)dS. (11) Випадок несингулярних iнтегралiв. Роз- глянемо спочатку iнтеграли 1) i 2) у несингуляр- ному випадку, тобто коли ~xi 0 6= ~ξj. Застосуємо до змiнної ~ξj = (ξ, η, ζ) координатнi перетворення (10). Таким чином, маємо: ξ = (ξ1 − ξ3)τ1 + (ξ2 − ξ3)τ2 + ξ3, η = (η1 − η3)τ1 + (η2 − η3)τ2 + η3, ζ = (ζ1 − ζ3)τ1 + (ζ2 − ζ3)τ2 + ζ3, (ξk, ηk, ζk) – координати т. Ak (рис. 2, с), ∫ ∆Sj G(~xi 0, ~ξj)dS = = ∫ ∆Sj dS 4π √ (x0 − ξ)2 + (y0 − η)2 + (z0 − ζ)2 = = √ ẼG̃− F̃ 2 4π 1 ∫ 0 1−τ1 ∫ 0 dτ1dτ2 r(τ1, τ2) , (12) r(τ1, τ2) = {[ (ξ1 − ξ3)τ1 + (ξ2 − ξ3)τ2 + ξ3 − x0 ]2 + + [ (η1 − η3)τ1 + (η2 − η3)τ2 + η3 − y0 ]2 + + [ (ζ1 − ζ3)τ1 + (ζ2 − ζ3)τ2 + ζ3 − z0 ]2}1/2 , Ẽ, G̃, F̃ – коефiцiєнти квадратичної форми для j-го дискретного елементу поверхнi ∆Sj : Ẽ = ( ∂ξ ∂τ1 )2 + ( ∂η ∂τ1 )2 + ( ∂ζ ∂τ1 )2 = = (ξ1 − ξ3) 2 + (η1 − η3) 2 + (ζ1 − ζ3) 2, G̃ = ( ∂ξ ∂τ2 )2 + ( ∂η ∂τ2 )2 + ( ∂ζ ∂τ2 )2 = = (ξ2 − ξ3) 2 + (η2 − η3) 2 + (ζ2 − ζ3) 2, F̃ = ∂ξ ∂τ1 ∂ξ ∂τ2 + ∂η ∂τ1 ∂η ∂τ2 + ∂ζ ∂τ1 ∂ζ ∂τ2 = (ξ1−ξ3)(ξ2−ξ3)+ +(η1 − η3)(η2 − η3) + (ζ1 − ζ3)(ζ2 − ζ3). (13) Аналогiчно: ∫ ∆Sj Gn(~xi, ~ξj)dS = = √ ẼG̃− F̃ 2 4π 1 ∫ 0 1−τ1 ∫ 0 ~r(τ1, τ2) · ~n(~x0) r3(τ1, τ2) dτ1dτ2. (14) Випадок дального поля. У випадку дальньо- го поля, тобто при умовi, що: r(~xi 0, ~ξj c) = = √ (x0 − ξc)2 + (y0 − ηc)2 + (z0 − ζc)2 ≥ 4 √ 2∆Sj , де ~ξj c = (ξc, ηc, ζc) – геометричний центр j-го еле- мента, для виразiв (12) i (14) можна застосувати чисельне iнтегрування по трикутнику [1, ст. 479]. Вираз для iнтегралу 3) в (11) отримується з (14) замiною нормалi ~n0 = (nx, ny, nz) на дотичний на- прямок ~τ0 = (τx, τy, τz). Випадок ближнього поля. У випадку бли- жнього поля, тобто при умовi ~r(~xi 0, ~ξj c) < 4 √ 2∆Sj вирази (12) та (14) потрiбно iнтегрувати з доста- тньо високою точнiстю. Спробуємо проiнтегрува- ти (14) по змiннiй τ2, для цього введемо наступнi позначення: a1 = ξ1 − ξ2; b1 = η1 − η3; c1 = ζ1 − ζ3; a2 = ξ2 − ξ3; b2 = η2 − η3; c2 = ζ2 − ζ3; a3 = ξ3 − x0; b3 = η3 − y0; c3 = ζ3 − z0; Ẽ1 = a1a3 + b1b3 + c1c3; G̃1 = a2a3 + b2b3 + c2c3; F̃1 = a2 3 + b23 + c23; A = a1nx + b1ny + c1nz; B = a2nx + b2ny + c2nz; C = a3nx + b3ny + c3nz; (15) тодi з спiввiдношень (13) маємо: Ẽ = a2 1 + b21 + c21; G̃ = a2 2 + b22 + c22; F̃ = a1a2 + b1b2 + c1c2. Тодi з спiввiдношень (14) матимемо: ∫ ∆Sj Gn(~xi, ~ξj)dS = 34 В. О. Горбань, С. В. Маcюк ISSN 1561 -9087 Прикладна гiдромеханiка. 2006. Том 8, N 3. С. 27 – 49 = √ ẼG̃− F̃ 2 4π 1 ∫ 0 1−τ1 ∫ 0 ~r∗(τ1, τ2)~n(~x0) (r∗(τ1, τ2))3 dτ1dτ2, (16) r∗(τ1, τ2) = [ (a1τ1 + a2τ2 + a3) 2+ +(b1τ1 + b2τ2 + b3) 2 + (c1τ1 + c2τ2 + c3) 2 ]1/2 . Нехай: c̃ = √ ẼG̃− F̃ 2 2π , c(τ1) = Ẽτ2 1 + 2Ẽ1τ1 + F̃1, r∗∗(τ1, τ2) = √ G̃τ2 2 + 2(F̃ τ1 + G̃1)τ2 + c(τ1). Розкриваючи дужки в спiввiдношеннi (16) i пере- творюючи його з використанням (15), матимемо: Gij n = c̃ 2 1 ∫ 0 1−τ1 ∫ 0 (Bτ2 + (Aτ1 + C))dτ1dτ2 √ (r∗∗(τ1, τ2))3 , (17) позначимо через ∆ – дискримiнант пiдкореневого виразу, який є квадратним тричленом вiдносно τ2. Таким чином: ∆ = G̃(Ẽτ2 1 + 2Ẽ1τ1 + F̃1) − (F̃ τ1 + G̃1) 2 = = (G̃Ẽ − F̃ 2)τ2 1 + 2(Ẽ1G̃− F̃ G̃1)τ1 + G̃F̃1 − G̃2 1. I. У випадку, якщо ∆ 6= 0, iнтегруючи (17) по τ2, отримаємо: Gij n = −c̃ 1 ∫ 0 B ( (F̃ τ1 + G̃1)τ2 + c(τ1) ) ∆ · r∗∗(τ1, τ2) ∣ ∣ ∣ ∣ τ2=1−τ1 τ2=0 dτ1+ +c̃ 1 ∫ 0 (Aτ1 + C) ( G̃τ2 + F̃ τ1 + G̃1 ) ∆ · r∗∗(τ1, τ2) ∣ ∣ ∣ ∣ τ2=1−τ1 τ2=0 dτ1 = = −c̃ 1 ∫ 0 B ( (F̃ τ1 + G̃1)(1 − τ1) + c(τ1) ) ∆ · r∗∗(τ1, 1− τ1) dτ1+ +c̃ 1 ∫ 0 (Aτ1 + C)(G̃(1 − τ1) + F̃ τ1 + G̃1) ∆ · r∗∗(τ1, 1− τ1) dτ1− −c̃ 1 ∫ 0 [ A(τ1 + C)(F̃ τ1 + G̃1) − Bc(τ1) ] ∆ · √ c(τ1) dτ1. (18) II. У випадку, коли ∆ = 0 маємо: c̃∗ = − √ ẼG̃− F̃ 2 8πG̃ · √ G̃ , F̃ ∗ = F̃ /G̃, G̃∗ 1 = G̃1/G̃, Gij n = −c̃∗ 1 ∫ 0 [ Aτ1 + C (τ2 + τ1F̃ ∗ + G̃∗ 1) · ∣ ∣τ2 + τ1F̃ ∗ + G̃∗ 1 ∣ ∣ − − B(2τ2 + τ1F̃ ∗ + G̃∗ 1) ( τ2 + τ1F̃ ∗ + G̃∗ 1 )∣ ∣τ2 + τ1F̃ ∗ + G̃∗ 1 ∣ ∣ ] ∣ ∣ ∣ ∣ τ2=1−τ1 τ2=0 dτ1 = = c̃∗ 1 ∫ 0 [ Aτ1 +C +B(τ1F̃ ∗ + G̃∗ 1) (τ1F̃ ∗ + G̃∗ 1) · ∣ ∣τ1F̃ ∗ + G̃∗ 1 ∣ ∣ − − Aτ1 + C +B ( 2(1 − τ1) + τ1F̃ ∗ + G̃∗ 1 ) ( 1 − τ1 + τ1F̃ ∗ + G̃∗ 1 ) · ∣ ∣1 − τ1 + τ1F̃ ∗ + G̃∗ 1 ∣ ∣ ] dτ1. (19) Отриманi iнтеграли (18) та (19) у випадку вiдсу- тностi особливостi, тобто при умовi, що ~ξ 6= ~x0 мо- жна проiнтегрувати з заданою точнiстю методом Гауса або Симпсона. Проводячи аналогiчнi мiркування, перетворимо вираз (12) до вигляду: Gij = c̃ 2 1 ∫ 0 1−τ1 ∫ 0 dτ1dτ2 √ G̃τ2 2 + 2(F̃ τ1 + G̃1)τ2 + c(τ1) . I. У випадку ∆ 6= 0: Gij = c̃ 2 √ G̃ 1 ∫ 0 ln \|Q\| ∣ ∣ ∣ ∣ τ2=1−τ1 τ2=0 dτ1 = = c̃ 2 √ G̃ 1 ∫ 0 ( ln ∣ ∣ ∣ ∣ Q∗ √ G̃c(τ1) + F̃ τ1 + G̃1 ∣ ∣ ∣ ∣ ) dτ1, (20) Q = 2 √ G̃(G̃τ2 2 + 2(F̃ τ1 + G̃1)τ2 + c(τ1))+ +2G̃τ2 + 2(F̃ τ1 + G̃1), Q∗ = G̃(1 − τ1) + F̃ τ1 + G̃1+ + √ G̃(G̃(1 − τ1)2 + 2(F̃ τ1 + G̃1)(1 − τ1) + c(τ1)). II. У випадку ∆ = 0: Gij = c̃ 2 √ G̃ 1 ∫ 0 ln ∣ ∣2G̃τ2 + 2(F̃ τ1 + G̃1) ∣ ∣× ×sign(2G̃τ2 + 2(F̃ τ1 + G̃1)) ∣ ∣ ∣ ∣ τ2=1−τ1 τ2=0 dτ1 = = c̃ 2 √ G̃ 1 ∫ 0 ln ∣ ∣2G̃(1 − τ1) + 2(F̃ τ1 + G̃1) ∣ ∣× В. О. Горбань, С. В. Маcюк 35 ISSN 1561 -9087 Прикладна гiдромеханiка. 2006. Том 8, N 3. С. 27 – 49 ×sign ( G̃(1 − τ1) + F̃ τ1 + G̃1 ) · dτ1 = = c̃ 2 √ G̃ 1 ∫ 0 ln ∣ ∣2(F̃ τ1 + G̃1) ∣ ∣sign ( F̃ τ1 + G̃1 ) dτ1. (21) Iнтеграли (20) та (21), так само як i (18) та (19) у випадку вiдсутностi особливостi, тобто при умовi, що ~xi 0 6∈ ∆Sj , можна проiнтегрувати, застосовую- чи чисельнi методи Гауса або Симпсона. Iнтеграл Gij τ = ∫ ∆Sj Fτ(~xi 0, ~ξj)dS знаходиться з формул (18) та (19) замiною нор- малi ~n0 = (nx, ny, nz) на дотичний напрямок ~τ0 = (τx, τy, τz). Випадок сингулярних iнтегралiв. У ви- падку, коли пiдинтегральна функцiя в iнтегра- лах (11) необмежено зростає в серединi областi iнтегрування, застосування простих схем чисель- ного iнтегрування неможливе. Можна дещо змi- стити точку спостереження в середину областi Ω, так що пiдинтегральна функцiя буде обмеженою i її можна буде проiнтегрувати за допомогою ква- дратурної формули бiльш високого порядку то- чностi, але такi процедури застосовувати не реко- мендується. Такi iнтеграли з особливiстю пiдинте- гральної функцiї краще всього обчислювати ана- лiтично. Розглянемо iнтеграл ∫ ∆Sj G(~xi 0, ~ξj)dS у випадку, коли ~xi 0 ∈ ∆Sj , тобто коли пiдинтегральна функ- цiя необмежено зростає при i = j. Для його обчи- слення введемо полярну локальну систему коор- динат 0Rθ, зв’язану з j-тим елементом, (рис. 3). Розiб’ємо область iнтегрування на трикутники 0A1A2, 0A2A3 та 0A3A1. Тодi загальний iнтеграл – це сума iнтегралiв по кожному з цих трикутни- кiв. Нехай θk та Rk – змiннi iнтегрування для k-го трикутника (k = 1, 2, 3). Таким чином маємо: Gjj = 1 4π 3 ∑ k=1 π−αk−βk ∫ 0 RL k ∫ 0 RkdRkdθk Rk = = 1 4π 3 ∑ k=1 π−αk−βk ∫ 0 dθk RL k ∫ 0 dRk, (22) де L – контур, що обмежує плоский елемент j, тоб- то: L = L1 +L2 +L3, де L1 = A1A2, L2 = A2A3, L3 = A1A3; RL k = \|~x0 − ~ξLk \|, ~ξLk ∈ Lk. (23) Рис. 3. Схема до обчислення сингулярного iнтегралу Gjj Тодi: Gjj = 1 4π 3 ∑ k=1 π−αk−βk ∫ 0 Rkdθk ∣ ∣ ∣ ∣ Rk=RL k Rk=0 = = 1 4π 3 ∑ k=1 π−αk−βk ∫ 0 RL kdθk. Зазначимо, що: RL k = hk cosϕk = = hk cos [ π/2 − (αk + θk) ] = hk sin(αk + θk) , (24) де hk – перпендикуляр, опущений з точки 0 на сторону Lk, a ϕk – кут, утворений цим перпеди- куляром з радiусом iнтегрування Rk. Враховуючи вираз (24), маємо: Gij = 1 4π 3 ∑ k=1 hk π−αk−βk ∫ 0 dθk sin(αk + θk) = = 1 4π 3 ∑ k=1 hk ln tg ( αk + θk 2 )∣ ∣ ∣ ∣ θk=π−αk−βk θk=0 = = 1 4π 3 ∑ k=1 hk ln 1 tg αk 2 tg βk 2 . (25) Розглянемо тепер iнтеграл ∫ ∆Sj Gn(~xi 0, ~ξj)dS у випадку, коли ~xi 0 ∈ ∆Sj , коли пiдинтегральна функцiя необмежено зростає при i = j. Для його обчислення, як i в попередньому випадку, введемо 36 В. О. Горбань, С. В. Маcюк ISSN 1561 -9087 Прикладна гiдромеханiка. 2006. Том 8, N 3. С. 27 – 49 Рис. 4. Схема до обчислення сингулярного iнтегралу Gjj n полярну локальну систему координат, зв’язану з елементом j, та розiб’ємо область iнтегрування на три трикутники 0A1A2, 0A2A3 та 0A3A1 (дивись рис. 4). Проведемо наступнi мiркування. Нехай точка ~x0 знаходиться на деякiй нескiнченно малiй вiдстанi δ → 0 вiд площини елемента ∆Sj в серединi обла- стi Ω. Тодi матимемо: Gjj n = = 1 4π 3 ∑ k=1 π−αk−βk ∫ 0 RL k ∫ 0 [(~ξ − ~x0)~n(~x0)]RkdRkdθk ( √ R2 k + δ2)3 = = 1 4π 3 ∑ k=1 π−αk−βk ∫ 0 dθk RL k ∫ 0 RkδdRk ( √ R2 k + δ2 )3 = = − 1 4π 3 ∑ k=1 π−αk−βk ∫ 0 δ √ R2 k + δ2 ∣ ∣ ∣ ∣ Rk=RL k Rk=0 dθk = = − 1 4π 3 ∑ k=1 ( π−αk−βk ∫ 0 δdθk √ (RL k )2 + δ2 − π + αk + βk ) . Враховуючи те, що δ → 0, одержимо: Gjj n = 1 4π 3 ∑ k=1 (π − αk − βk) = = 1 4π [ 3π − 3 ∑ k=1 (αk + βk) ] = 1 4π (3π − π) = 1 2 . (26) Рис. 5. Схема до обчислення сингулярного iнтегралу Gjj τ Тепер розглянемо iнтеграл ∫ ∆Sj Gτ(~xi 0, ~ξj)dS у ви- падку ~xi 0 ∈ ∆Sj , тобто коли iснує особливiсть при i = j. Для його обчислення введемо локальну по- лярну систему координат, зв’язану з j-м елемен- том, та розiб’ємо область iнтегрування на три три- кутники 0A1A2, 0A2A3 та 0A3A1, (рис. 5). Нехай θk та Rk – вiдповiднi змiннi iнтегрування для k-го трикутника (k = 1, 2, 3); β∗ k – кут, утворений ра- дiусом RL k та вектором (~ξ − ~x0), а α∗ k – кут мiж радiусом RL k та вектором ~τ (~x0), γk – кут мiж ве- кторами (~ξ − ~x0) та ~τ (~x0). Нехай ~τ (x0) – одиничний вектор, що лежить у площинi елемента j i напрямлений, наприклад, па- ралельно ~OA1. Проведемо наступнi мiркування. Нехай точка ~x0 знаходиться на деякiй нескiнчен- но малiй вiдстанi δ → 0 вiд площини елемента ∆Sj назовнi областi Ω. Тодi матимемо: Gjj τ = 1 4π 3 ∑ k=1 π−αk−βk ∫ 0 RL k ∫ 0 [ (~ξ − ~x0)~τ (~x0) ] RkdRkdθk (√ R2 k + δ2 )3 = = − 1 4π 3 ∑ i=1 π−αk−βk ∫ 0 RL k ∫ 0 ( √ R2 k + δ2Rk cos γkdRkdθk (√ R2 k + δ2)3 ) . Враховуючи те, що: cos γk = cosα∗ k · cos β∗ k = = − cos θ∗k cos β∗ k = − cos θ∗k Rk√ Rk + δ2 , В. О. Горбань, С. В. Маcюк 37 ISSN 1561 -9087 Прикладна гiдромеханiка. 2006. Том 8, N 3. С. 27 – 49 де θ∗1 = θ1; θ∗2 = θ2 + π − α1 − β1; θ∗3 = θ3 + 2π − α1 − β1 − α2 − β2; будемо мати: Gjj τ = 1 4π 3 ∑ k=1 π−αk−βk ∫ 0 RL k ∫ 0 R2 k cos θ∗kdRkdθk (√ R2 k + δ2 )3 = = 1 4π 3 ∑ k=1 π−αk−βk ∫ 0 [ − Rk (√ R2 k + δ2 )+ + ln ( Rk + √ R2 k + δ2 ) ]∣ ∣ ∣ ∣ RL k Rk=0 cos θ∗kdθk = = 1 4π 3 ∑ k=1 π−αk−βk ∫ 0 − ln δ cos θ∗kdθk− − π−αk−βk ∫ 0 RL k cos θ∗k √ (RL k )2 + δ2 dθk+ + π−αk−βk ∫ 0 ln ( RL k + √ (RL k )2 + δ2 ) cos θ∗kdθk = = 1 4π 3 ∑ k=1 [ − π−αk−βk ∫ 0 RL k cos θ∗kdθk √ (RL k )2 + δ2 + + π−αk−βk ∫ 0 ln ( RL k + √ (RL k )2 + δ2 ) cosθ∗kdθk ] . Враховуючи те, що δ → 0, а також формулу (24), отримаємо: Gjj τ = = − 1 4π 3 ∑ k=1 π−αk−βk ∫ 0 ln ( sin(αk + θk) hk ) cos θ∗kdθk. (27) Вираз (27) вже не мiстить особливостей, тому мо- же бути проiнтегрований методом Симпсона або Гауса з заданою точнiстю. Iдея Хесса i Смiта щодо обчислення iнте- гралiв. Для обчислення iнтегралiв типу (11) Хесс i Смiт [14] запропонували описаний нижче пiд- хiд. Спочатку зробимо перетворення системи ко- ординат так, щоб початок нової системи коорди- нат A4 спiвпадав з точкою перетину перпендику- ляра, проведеного з т. ~x0 до площини трикутника Рис. 6. Схема до обчислення iнтегралiв за методом Хеса та Смiта A1A2A3 див. (рис. 6), вiсь τ3 нової системи коорди- нат була спрямована по нормалi до площини три- кутника A1A2A3, а осi τ1 та τ2 лежали будь-яким чином у площинi цього трикутника. Знайдемо координати точки A4 в старiй системi координат 0xyz. Нехай ~n(nξ, nη, nζ) – нормаль до площини A1A2A3, тодi рiвняння площини A1A2A3 буде мати вигляд: (x− ξ3)nξ + (y − η3)nη + (z − ζ3)nζ = 0, (28) а рiвняння перпендикуляра до площини A1A2A3, що проходить через точку ~x0, запишеться у вигля- дi: (x− x0) nξ + (y − y0) nη + (z − z0) nζ = −p. (29) З рiвнянь (28) i (29) знаходимо координати т. A4 : ξ4 = x0 − nξ p, η4 = y0 − nη p, ζ4 = z0 − nζ p, (30) де p = x0nξ + y0nη + z0nζ − (ξ3nξ + η3nη + ζ3nζ). Виконаємо перетворення системи координат на- ступним чином: O(0, 0, 0) → A4(ξ4, η4, ζ4), ~k(0, 0, 1) → ~n(nξ, nη, nζ), ~i(1, 0, 0) → ~A3A1 \|A3A1\| ( ξ1 − ξ3 \|A3A1\| , η1 − η3 \|A3A1\| , ζ1 − ζ3 \|A3A1\| ) , ~j(0, 1, 0) → ~n × ~A3A1 \|A3A1\| , 38 В. О. Горбань, С. В. Маcюк ISSN 1561 -9087 Прикладна гiдромеханiка. 2006. Том 8, N 3. С. 27 – 49 де \|A3A1\| = √ (ξ1 − ξ3)2 + (η1 − η3)2 + (ζ1 − ζ3)2. Тобто початок нової системи повинен знаходитись у точцi A4, вiсь τ3 спрямована по нормалi до пло- щини трикутника A1A2A3, вiсь τ1 – паралельна сторонi трикутника A1A3, а вiсь τ2 спрямована так, що система координат τ1τ2τ3 – прямокутна правостороння. Введемо наступнi позначення для косинусiв кутiв мiж новими та старими осями ко- ординат: l1 = ξ1 − ξ3 \|A3A1\| ; m1 = η1 − η3 \|A3A1\| ; n1 = ζ1 − ζ3 \|A3A1\| ; l2 = m3n1 − n3m1; m2 = n3l1 − l3n1; n2 = l3m1 −m3l1; l3 = nξ; m3 = nη; n3 = nζ ; тодi перетворення координат запишеться насту- пним чином: τ1 = l1x+m1y + n1z − ξ4; τ2 = l2x+m2y + n2z − η4; τ3 = l3x+m3y + n3z − ζ4. Тепер перейдемо безпосередньо до обчислення iн- тегралiв. Зазначимо, що ∫ ∆A1A2A3 f(R, θ)dR dθ = = S12 ∫ ∆A4A1A2 f(R, θ)dR dθ+ +S23 ∫ ∆A4A2A3 f(R, θ)dR dθ+ +S31 ∫ ∆A4A3A1 f(R, θ)dR dθ, (31) де Sij = 1, якщо центр нової системи коорди- нат A4 i будь-яка точка в серединi трикутника A1A2A3, наприклад, точка перетину його медiан M(τ1M , τ2M , 0), лежить по одну сторону прямої, проведеної через точки Ai та Aj, i Sij = −1, якщо – по рiзнi сторони. З курсу аналiтичної геометрiї вiдомо, що Sij = −sign(aijτ1M + bijτ2M − cij); cij > 0; i, j = 1, 2, 3, i 6= j, де aij, bij, cij− коефiцiєнти в нормальному рiвнян- нi прямої, що проходить через AiAj. Враховуючи викладки, приведенi в попередньому пунктi, фор- мули (22)–(27), а також формулу (31), можемо за- писати: Gij = 1 4π 3 ∑ k=1 S∗ k βk ∫ αk RL k ∫ 0 RkdRk dθk √ R2 k +H2 = = 1 4π 3 ∑ k=1 S∗ k βk ∫ αk √ R2 k +H2 ∣ ∣ Rk=RL k Rk=0 dθk = = 1 4π 3 ∑ k=1 S∗ k βk ∫ αk √ (RL k )2 +H2 dθk, де S∗ 1 = S12, S ∗ 2 = S23, S ∗ 3 = S31, H – вiдстань мiж точками ~x0 i A4. Враховуючи те, що: RL k = ∣ ∣ ∣ ∣ ∣ c∗k a∗k cos θk + b∗k sin θk ∣ ∣ ∣ ∣ ∣ , (32) де a∗1 = a12, a∗2 = a23, a∗3 = a31; b∗1 = b12, b∗2 = b23, b∗3 = b31; c∗1 = c12, c∗2 = c23, c∗3 = c31, остаточно будемо мати: Gij = 1 4π 3 ∑ k=1 S∗ k× × βk ∫ αk √ ( c∗k a∗k cos θk + b∗k sin θ )2 +H2 dθk. Аналогiчним чином iнтегруємо функцiю Gij n : Gij n = 1 4π 3 ∑ k=1 S∗ k× × βk ∫ αk RL k ∫ 0 ( ~ξ − ~x0)~n(~x0 ) RkdRkdθk (√ R2 k +H2 )3 = = 1 4π 3 ∑ k=1 S∗ k βk ∫ αk [ − Hnz √ R2 k +H2 − Rknx cos θk √ R2 k +H2 + +ln ( Rk + √ R2 k +H2 ) nx cos θk − Rkny sin θk √ R2 k +H2 + +ln ( Rk + √ R2 k +H2 ) ny sin θk ] dθk ∣ ∣ ∣ ∣ Rk=RL k Rk=0 = В. О. Горбань, С. В. Маcюк 39 ISSN 1561 -9087 Прикладна гiдромеханiка. 2006. Том 8, N 3. С. 27 – 49 = 1 4π 3 ∑ k=1 S∗ k βk ∫ αk [( 1 − H √ (RL k )2 +H2 ) nz+ + ( ln ( \|RL k + √ (RL k )2 +H2 ) − RL k √ (RL k )2 +H2 ) × × ( nx cos θk + ny sin θk ) ] dθk, (33) де RL k – визначається формулою (23). Iнтеграли (32) та (33) не мають особливостей. Для їх визна- чення використовується метод Симпсона або Гау- са. 3. ВИЗНАЧЕННЯ СИЛ ГIДРОДИНАМIЧНОЇ ВЗАЄМОДIЇ ТIЛ При маневруваннi декiлькох суден гiдродинамi- чна взаємодiя мiж ними може викликати значнi додатковi навантаження (сили i моменти). Цi гi- дродинамiчнi навантаження особливо небезпечнi, коли вiдстань мiж суднами невелика, а також при руховi в каналi або на мiлководдi. Найбiльших значень додатковi навантаження набувають при обгонi одного судна iншим, коли швидкoстi V1 i V2 близькi, або коли одне з суден пришвартоване (V1 = 0, V2 6= 0). В данiй статтi ми розглянемо саме останнiй випадок, хоча методологiя, яку ми буде- мо застосовувати, є загальною для всiх взаємних рухiв суден. Щоб спростити викладки, введемо узагальненi координати xi, i = 1, · · · , 12, де x1, x2, x3 – змiщен- ня геометричного центру першого тiла вiдносно початкового положення, x4, x5, x6 – кути поворо- ту першого тiла вiдносно його центру; x7, · · · , x12 – вiдповiдно змiщення i кути повороту другого тi- ла, узагальненi швидкостi ui = ẋi. Аналогiчним чином введемо потенцiали ϕi та ni = ∂ϕ/∂n, i = 1, · · · , 12. Функцiї ϕ1, ϕ2, · · · , ϕ12 вiдповiдають рi- зним випадкам руху тiл. (Наприклад, ϕ1− випад- ку руху тiл, коли u1 = 1, u2 = u3 = · · · = u12 = 0). Кожний потенцiал ϕi задовiльняє рiвнянню Ла- пласа (1). Узагальненi сили Fi можна знайти як iнтеграл вiд тиску по зануренiй поверхнi тiла SB : Fi = ∫ ∫ SB p ni dS = = −ρ ∫ ∫ SB ( ∂ ϕ ∂ t + 1 2 ∇ϕ · ∇ϕ ) ni dS. (34) Тут ∂/∂t – похiдна по часу в нерухомiй системi координат. Якщо перейти до системи координат, зв’язаної з тiлом, то рiвняння (34) запишеться у виглядi: Fi = −ρ d dt ∫ ∫ SB ϕnidS+ +ρ ∫ ∫ SB ( (~v0 + [~ω0 × ~r])∇ϕ− 1 2 ∇ϕ∇ϕ ) nidS, (35) де ~v0 та ~ω0 вiдповiдно поступальна та обертальна швидкостi рухомої системи координат. Другий пiдхiд до визначення сил i моментiв по- в’язаний з застосуванням рiвнянь Лагранжа 2-го роду [3, 4, 9, 16]. Якщо границi тiл та каналу не- проникнi, для кiнетичної енергiї рiдини маємо: T = λjkujuk 2 , λjk = −ρ ∫ SB ϕj ∂ϕk ∂n dS (j, k = 1, 2, . . . , 12), (36) де λjk – тензор приєднаних мас, елементи якого за- лежать вiд форми границь. Тут i далi по iндексах, що повторюються, потрiбно виконувати пiдсумо- вування. Узагальненi сили обчислюються iз рiвняння Ла- гранжа: Fi = − d dt ∂T ∂ui + ∂T ∂xi , ∂T ∂xi = 1 2 ujuk ∂λjk ∂xi , ∂T ∂ui = 1 2 λjk(ukσij + ujσik) = mijuj , σij та σik – символи Кронекера, { σij = 0, i 6= j, σij = 1, i = j, d dt ∂T ∂ui = λij u̇j + ujuk ∂λij ∂xk . Пiсля ряду тотожних перетворень одержимо: Fi = −λij u̇j − ujuk (∂λij ∂xk − 1 2 ∂λjk ∂xi ) . (37) У випадку, коли одне тiло рухається з постiйною швидкiстю Vi, а iнше залишається нерухомим, ма- ємо: Fi = −V 2 1 (∂λi1 ∂x1 − 1 2 ∂λ11 ∂xi ) = 40 В. О. Горбань, С. В. Маcюк ISSN 1561 -9087 Прикладна гiдромеханiка. 2006. Том 8, N 3. С. 27 – 49 a б Рис. 7. Порiвняння гiдродинамiчних потенцiалiв ϕ̄, обчислених аналiтично та за допомогою МГЕ на лiнiї перетину площини 0xz та поверхнi сфери (а) та елiпсоїда (б): 1 – аналiтичний розв’язок; 2 – розвязок МГЕ = −V 2 1 (∂λ1i ∂x1 − 1 2 ∂λ11 ∂xi ) . Коли два тiла рухаються з постiйними швидко- стями V у необмеженiй рiдинi або вздовж осi ка- налу постiйного поперечного перерiзу, кiнетична енергiя рiдини буде залежати лише вiд їхнього вза- ємного розташування: Fi = −1 2 V 2 ∂ ∂xi (λ11 + 2λ17 + λ77). Пiдхiд, що грунтується на використаннi систе- ми рiвнянь Лагранжа 2-го роду, часто є простi- шим i ефективнiшим. (Саме вiн застосовувався в розрахунках, результатi яких наведенi нижче). Якщо швидкiсть визначена експериментально, з напiвемпiричних мiркувань чи вирахувана з доста- тньою точнiстю, для визначення сил краще скори- статись рiвнянням (35). a б Рис. 8. Порiвняння дотичних швидкостей Vτ , обчислених аналiтично та за допомогою МГЕ на лiнiї перетину площини 0xz та поверхнi сфери (а) та елiпсоїда (б): 1 - аналiтичний розв’язок; 2 - розв’язок МГЕ 4. РЕЗУЛЬТАТИ ЧИСЕЛЬНОГО МОДЕЛЮВАННЯ Апробацiя алгоритму. Чисельний алгоритм, описаний в роздiлi 2, було запрограмовано та апро- бовано для наступних тестових задач. 1. Сфера рухається в необмеженiй рiдинi в на- прямку осi x системи координат 0xy. Для введе- ння безрозмiрних величин використовується радi- ус сфери R та її швидкiсть V0. Кiлькiсть плоских елементiв на сферi n змiнювалась (для приведених результатiв розрахункiв n = 3120). 2. Елiпсоїд рухається в необмеженiй рiдинi у на- прямку осi x. Пiвосi елiпсоїда: b/a = 1/6, c/a = = 1/12. Система координат 0xyz вибрана таким чином, що велика пiввiсь елiпсоїда a збiгається з вiссю x, пiввiсь b – з вiссю y, а пiввiсь c – з вiссю z. Кiлькiсть плоских елементiв на елiпсоїдi n = 3120. Для дискретiзацiї поверхнi тiла використовува- лась система кутових координат (ϕ ′ , ψ ′ ): коорди- В. О. Горбань, С. В. Маcюк 41 ISSN 1561 -9087 Прикладна гiдромеханiка. 2006. Том 8, N 3. С. 27 – 49 нати точки на поверхнi вираховувались за фор- мулами: x = cosψ ′ , y = b sinψ ′ cosϕ ′ та z = = c sinψ ′ sinϕ ′ . Причому кути ϕ ′ i ψ ′ змiню- вались рiвномiрно з одинаковим кроком π/30. Рис. 9. Порiвняння дотичних швидкостей з врахуванням знаку на лiнiї перетину 0xz та поверхнi елiпсоїда поблизу дна, обчислених за допомогою рiзних пiдходiв. 1 – алгоритм МГЕ з дискретизацiєю поверхнi дна; 2 – використання спецiальної функцiї Грiна i побудованої з застосуванням вiдображень a б Рис. 10. Поперечний перерiз симетричного (а) та несиметричного (б) каналу Одержанi значення гiдродинамiчних потенцiа- лiв ϕ̄ = ϕ/V0R для сфери i ϕ̄ = ϕ/V0a для елi- псоїда та безрозмiрних дотичних швидкостей (на поверхнi тiла) V̄τ = Vτ/V0 викликаних рухом цього тiла, порiвнювалися з вiдповiдними аналiтичними розв’язками, рис. 7–8. Точки, позначенi на малюн- ках, вiдповiдають геометричним центрам дискре- тних елементiв поверхнi (панелей), в яких розра- ховувались значення величин. Значення гiдродинамiчних потенцiалiв для сфе- ри та елiпсоїда наведенi на рис. 7. Вони майже точ- но збiгаються з аналiтичними (вiдносна похибка у випадку елiпсоїда складає близько 3%). Стосов- но дотичних швидкостей (рис. 8), розбiжнiсть мiж чисельними та аналiтичними розв’язками бiльша, (на елiпсоїдi похибка досягає 10%). Це пояснює- ться особливостями розбиття поверхнi на дискре- тнi елементи. Рух елiпсоїда поблизу дна. Для розрахунку поля швидкостей при руховi елiпсоїда поблизу дна при малих числах Фруда можна використовувати два еквiвалентнi пiдходи [16]. 1. Розбити дно та нижню половину елiпсоїда на плоскi панелi, ввести на кожнiй панелi фiктивнi джерела iнтенсивностi q, виконати граничнi умови на днi та на елiпсоїдi i включити одержанi рiвнян- ня в систему (9). При необхiдностi, таким чином можна врахувати i деформацiї вiльної поверхнi. a б Рис. 11. Форма денiвеляцiї вiльної поверхнi при руховi елiпсоїда в симетричному (а) та несиметричному (б) каналi 2. Iнший пiдхiд передбачає використання функ- цiї Грiна джерела, яка задовiльняє граничну умо- ву непротiкання на плоскому днi. Така функцiя 42 В. О. Горбань, С. В. Маcюк ISSN 1561 -9087 Прикладна гiдромеханiка. 2006. Том 8, N 3. С. 27 – 49 a б Рис. 12. Лiнiї течiї на поверхнi елiпсоїда при руховi в необмеженiй рiдинi в проекцiї: а – на вертикальну площину 0yz; б – на горизонтальну площину 0xy a б Рис. 13. Лiнiї течiї на поверхнi елiпсоїда в проекцiї на вертикальну площину 0yz (а) та горизонтальну площину 0xy (б) при руховi в симетричному каналi а б Рис. 14. Лiнiї течiї на поверхнi елiпсоїда в проекцiї на вертикальну (а) i горизонтальну площину 0xy (б) в несиметричному каналi поблизу стiнки а б Рис. 15. Лiнiї течiї на днi для випадку обтiкання елiпсоїда: а – в симетричному каналi; б – в несиметричному каналi поблизу стiнки В. О. Горбань, С. В. Маcюк 43 ISSN 1561 -9087 Прикладна гiдромеханiка. 2006. Том 8, N 3. С. 27 – 49 будується методом дзеркальних вiдображень i має вигляд нескiнченного ряду (5). Як показує досвiд проведення розрахункiв, щоб досягнути задовiль- ної точностi, можна обмежитись 41-м членом ряду (5), тобто k = −20, . . . , 20. У цьому випадку на па- нелi слiд розбивати лише нижню половину елiпсо- їда. Порiвняння дотичних швидкостей Vτ на лiнiї перетину площини 0xz та поверхнi елiпсоїда при- веденi на рис. 9. При цьому вибирались наступнi параметри: a = 1, b/a = 1/6, c/a = 1/12, h/c = 1, 1, розмiр дiлянки дна 6a×4a, h− глибина води. Осi системи координат спрямованi таким чином: велика пiввiсь елiпсоїда a збiгається з вiссю 0x, пiввiсь b – з вiссю 0y, а пiввiсь c – з вiссю 0z. Кiль- кiсть елементiв (панелей), розташованих на поло- винi елiпсоїда – 1740, на днi – 6912. З рис. 9 видно, що результати обох пiдходiв повнiстю спiвпада- ють. Розбиття дна на панелi дає змогу розв’язува- ти задачi з нерiвним дном, наприклад, при руховi елiпсоїда у каналi з вiдкосами. Слiд зауважити, що поздовжня дотична до поверхнi елiпсоїда швид- кiсть (вздовж стiнгера) бiля дна має двi особливi точки (Vτ = 0): на носику елiпсоїда та на певнiй невеликiй вiдстанi вiд носика, в нашому випадку це точка x/a = 0.8. Виникнення другої особливої точки пов’язане з безпосередньою близькiстю дна. При вiддаленнi вiд дна на певну вiдстань, що зале- жить вiд параметрiв елiпсоїда, ця особлива точка зникає. Гiдродинамiчнi поля та структура течiї при руховi елiпсоїда в симетричному i не- симетричному каналах. Розглянемо рух тривi- сного елiпсоїда в каналi при малих числах Фруда (рис. 10). Розрахунки проводились для руху тривiсного елiпсоїда. Розглядались канали з двома типами поперечних перерiзiв – симетричним (pис. 10, а) i несиметричним (pис. 10, б), вздовж яких рухає- ться елiпсоїд 1 з пiвосями b/a = 1/6, c/a = 1/12, (рис. 10), 1 – перерiз елiпсоїда, 2 – вiдкiс каналу; 3 – дно каналу; 4 – вiльна поверхня; 5 – вертикальна стiнка каналу. Прямокутна система координат ви- брана таким чином, що напрямок осi 0x збiгається з напрямком руху елiпсоїда, вiсь 0z спрямована вгору перпендикулярно до вiльної поверхнi, а вiсь 0y спрямована так, щоб утворилась правосторон- ня система координат. Задачi про рух елiпсоїда в симетричному (рис.10, а) та несиметричному (рис. 10, б) каналах розв’язувались чисельно методом граничних еле- ментiв з використанням описаного алгоритму. Роз- мiри каналiв вибирались наступним чином: h/c = 1, 1; ∆h/c = 1/12, d/b = 6, ∆d/b = 3, α = π/30. В обох задачах на плоскi елементи розбивалась нижня частина елiпсоїда, дно i вiдкоси каналу 2. Довжина дiлянки каналу, на якiй проводилась дискретизацiя, L = 6a. У випадку симетричної за- дачi для моделювання вiльної поверхнi була ви- користана функцiя Грiна, отримана на основi вi- дображення каналу та нижньої частини елiпсоїда вiдносно площини 0xy. У випадку несиметричного каналу додатково проводилось вiдображення вiд- носно вертикальної стiнки 5. В результатi розрахункiв були отриманi поля швидкостей при обтiканнi елiпсоїда, розташовано- го в симетричному та несиметричному каналах. На основi результатiв розрахункiв потенцiала та швидкостi, застосовуючи рiвняння Бернуллi, не- важко побудувати залежностi для розподiлу тис- ку, зокрема на поверхнi води (вiльнiй поверхнi). Змiна тиску бiля елiпсоїда зумовлює вiдповiдний пiдйом (чи падiння) рiвня води hw – денiвеля- цiю вiльної поверхнi. Одержанi форми денiвеляцiї вiльних поверхонь, рис. 11 (hw вiдхилення рiвня води вiд площини 0xy). Картини лiнiй течiї для рiзних випадкiв руху елiпсоїда показанi на рис. 12–15. Порiвняння картини течiї бiля елiпсоїда в без- межнiй рiдинi (рис. 12) та в каналах (рис. 13, 14) показує, що в каналi вiдбувається iнтенсивне роз- тiкання води вiд поздовжньої площини симетрiї елiпсоїда. Тiло обтiкається бiльшою мiрою з бокiв. Асиметрiя каналу зумовлює несиметрiю течiї нав- коло елiпсоїда: спостерiгається перетiкання рiдини через зазор мiж дном та елiпсоїдом (рис. 14). Розрахунок приєднаних мас. Основними складовими при розрахунку сил гiдродинамiчної взаємодiї є приєднанi маси тiл та їх похiднi. При- єднанi маси в загальному випадку визначаються з рiвняння (36). Враховуючи дискретизацiю поверх- нi, iнтеграл у правiй частинi рiвняння (36) вира- ховується наступним чином: λjk = NB ∑ i=1 ϕi j ∂ϕi k ∂n ∆Si, де NB – кiлькiсть панелей на тiлi; ∆Si – площа i-тої панелi; ϕi j – значення потенцiалу ϕj на i-тiй панелi (визначається з рiвнянь (9)). Порiвняння коефiцiєнтiв приєднаних мас розра- хованих чисельно для тривiсного елiпсоїда, що ру- хається в необмеженiй рiдинi, з точними значен- нями [5] та вiдноснi похибки наведенi в табл. 1. З даних таблицi видно, що найбiльша вiдносна по- хибка досягає 13, 5% при обчисленнi λ55. Це по- в’язано зi способом розбиття поверхнi елiпсоїда на дискретнi елементи. Цю похибку можна зменши- ти, розбивши поверхню iншим способом. Для всiх 44 В. О. Горбань, С. В. Маcюк ISSN 1561 -9087 Прикладна гiдромеханiка. 2006. Том 8, N 3. С. 27 – 49 Табл 1. Приєднанi маси елiпсоїда b/a = 1/6, c/a = 1/12, що рухається в необмеженiй рiдинi кiлькiсть час розрахунку λ11/m δλ11,% λ22/m δλ22,% λ33/m δλ33,% площадок n t точне значення – 0.046 – 0.842 – 3.342 – n = 180 t = 5c 0.047 1.4 0.863 2.5 3.436 2.8 n = 760 t = 60c 0.047 1.9 0.862 2.3 3.422 2.4 n = 1740 t = 360c 0.047 1.5 0.858 1.9 3.395 1.6 Кiлькiсть час розрахунку λ44/Jxx δλ44,% λ55/Jyy δλ55,% λ66/Jzz δλ66,% площадок n t точне значення – 0.443 – 1.698 – 0.408 – n = 180 t = 5c 0.435 1.8 1.469 13.5 0.404 0.9 n = 760 t = 60c 0.448 1.1 1.485 12.6 0.414 1.5 n = 1740 t = 360c 0.448 1.3 1.498 11.8 0.414 1.5 Табл 2. Приєднанi маси елiпсоїда b/a = 1/6, c/a = 1/12 залежно вiд глибини для випадку симетричного каналу h/c λ11/m λ22/m λ33/m λ44/Jxx λ55/Jyy λ66/Jzz 1.1 0.266 3.843 8.842 0.814 3.430 0.920 1.3 0.216 2.483 6.848 0.629 2.878 0.721 1.5 0.182 1.935 5.866 0.553 2.600 0.623 1.7 0.157 1.635 5.277 0.513 2.334 0.560 1.9 0.139 1.450 4.885 0.488 2.192 0.521 3.1 0.086 1.039 3.944 0.458 1.900 0.446 5.1 0.062 0.912 3.595 0.450 1.774 0.421 7.1 0.054 0.881 3.495 0.449 1.743 0.416 9.1 0.051 0.870 3.455 0.449 1.733 0.415 11.1 0.049 0.865 3.435 0.448 1.728 0.414 iнших випадкiв похибка у розрахунку приєднаних мас не перевищує 3%. В табл. 2 наведенi приєднанi маси тривiсного елiпсоїда, що рухається в симетричному каналi (рис. 10, а), залежно вiд глибини каналу. Всi при- єднанi маси вiднесенi до маси витiсненої рiдини m = 2 3πρabc або до вiдповiдного моменту iнерцiї: Jxx = 2 15 πρabc(b2 + c2) Jyy = 2 15 πρabc(a2 + c2) Jzz = 2 15 πρabc(a2 + b2) На рис. 16 показанi залежностi приєднаних мас вiд глибини каналу. Видно, що при збiльшеннi гли- бини каналу приєднанi маси рiзко зменшуються, а при досягненнi певної глибини майже не змiнюю- ться (асимптотично наближаються до певного зна- чення). Результати чисельного моделювання гi- дродинамiчної взаємодiї двох тривiсних елi- псоїдiв у безграничнiй рiдинi. Для оцiнки то- чностi побудованих чисельних алгоритмiв розра- хунку сил та моментiв, що виникають внаслiдок гiдродинамiчної взаємодiї тiл, проведено порiвнян- ня результатiв розрахункiв з даними аеродинамi- чного експерименту роботи [22], в якiй представле- нi результати випробувань в аеродинамiчнiй трубi двох елiпсоїдiв обертання зi спiввiдношенням осей 6:1. Сили i моменти вимiрювались за допомогою високочутливої тензосистеми з шiстьма ступеня- ми свободи. Обидвi моделi встановлювались в тру- бi паралельно з рiзними вiдстанями ∆η i ∆ξ (рис. 17). В експериментах тiла обдувалися при рiзних поздовжнiх змiщеннях ∆ξ. Така постановка задачi В. О. Горбань, С. В. Маcюк 45 ISSN 1561 -9087 Прикладна гiдромеханiка. 2006. Том 8, N 3. С. 27 – 49 а б Рис. 16. Залежностi коефiцiєнтiв приєднаних мас елiпсоїда вiд глибини h/c у випадку руху елiпсоїда в симетричному каналi Рис. 17. Схема руху двох елiпсоїдiв Рис. 18. Залежнiсть коефiцiєнта бокової гiдродинамiчної сили Cy при руховi двох елiпсоїдiв обертання (∆η/2a = 0.288) однакових розмiрiв в безграничнiй рiдинi вiд безрозмiрного поздовжнього змiщення ∆ξ/2a: 1– результати розрахункiв; 2 – експериментальнi данi Рис. 19. Залежнiсть коефiцiєнтiв бокових сил, що виникають у випадку, коли один з елiпсоїдiв нерухомий, а iнший з однаковою швидкiстю рухається повз нього в безграничнiй рiдинi (∆η/b = 4) вiд безрозмiрного поздовжнього змiщення ∆ξ/2a: 1 – коефiцiєнт сили, що дiє на рухомий елiпсоїд; 2 – коефiцiєнт сили, що дiє на нерухомий елiпсоїд вiдповiдає руховi двох тiл з однаковими швидко- стями. Проведенi в [22] розрахунки з використан- ням моделi “тонкого” тiла показали їх якiсну вiд- повiднiсть експериментальним даним. Для таких же умов були проведенi розрахунки з використанням побудованого вище чисельного ал- горитму МГЕ, який дозволяє виконати умови не- протiкання набагато точнiше, нiж теорiя “тонкого” тiла. На рис. 18 наведенi графiки порiвняння ко- ефiцiєнтiв бокової сили Cy = Fy/2ρV 2ab залежно вiд поздовжнього змiщення ∆ξ. З рисунка видно, що експериментальнi данi добре узгоджуються з 46 В. О. Горбань, С. В. Маcюк ISSN 1561 -9087 Прикладна гiдромеханiка. 2006. Том 8, N 3. С. 27 – 49 а б Рис. 20. Вплив вiдстанi мiж тiлами (1−∆η/b = 3; 2−∆η/b = 4; 3−∆η/b = 5; 4−∆η/b = 8) у випадку “тонких” тiл (b/a = 0.1; c/a = 0.05) на гiдродинамiчнi коефiцiєнти бокової сили Cy (a) та моменту Cm (б) на нерухоме тiло при проходженнi повз нього iншого тiла, таких же розмiрiв y безграничнiй рiдинi результатами чисельного моделювання. Незначне перевищення експериментальних даних для коефi- цiєнта бокової сили над розрахунковими значення- ми може бути пов’язане з впливом пограничного шару та вiдповiдним зменшенням вiдстанi мiж тi- лами. Цiкавим є випадок, коли тiла рухаються з рiзни- ми швидкостями або коли одне тiло нерухоме, а iнше рухається повз нього. Як показують дослiди, проведенi в роботi [12], за умови однакових роз- мiрiв обох тiл сила, що дiє на нерухоме тiло, на- багато бiльша, нiж сила, яка дiє на рухоме тiло. Цей неочевидний факт пiдтверджується наближе- ними розрахунками (рис. 19). Тому сили, якi дiють а б Рис. 21. Вплив вiдстанi мiж тiлами (1−∆η/b = 3; 2−∆η/b = 4; 3− η/b = 5; 4−∆η/b = 8) у випадку “товстих” тiл (b/a = 0.25; c/a = 0.0125) на гiдродинамiчнi коефiцiєнти бокової сили Cy (a) та моменту Cm (б) на нерухоме тiло при проходженнi повз нього iншого тiла, таких же розмiрiв в безграничнiй рiдинi на нерухоме тiло, викликають бiльший iнтерес, на- приклад, у випадку руху одного судна поблизу за- якореного (або пришвартованого) судна чи гiдро- технiчної споруди. Розглядаються лише сили, якi дiють на нерухоме тiло. На рис. 20 наведенi результати розрахункiв ко- ефiцiєнтiв бокових гiдродинамiчних сил Cy = = Fy/2ρV 2ab та моментiв Cm = Mz/4ρV 2a2b, що дiють на нерухомi “тонкi” тiла при рiзних вiд- станях мiж рухомим та нерухомим тiлами η. Тi- лами в даному випадку служать тривiснi елiпсо- їди однакових розмiрiв зi спiввiдношеннями осей b/a = 0.1, c/a = 0.05. З цих графiкiв видно, що В. О. Горбань, С. В. Маcюк 47 ISSN 1561 -9087 Прикладна гiдромеханiка. 2006. Том 8, N 3. С. 27 – 49 а б Рис. 22. Гiдродинамiчнi коефiцiєнти бокової сили Cy (а) i моменту Cm (б), що дiють на нерухомий елiпсоїд (b1/a1 = 0.25; c1/a1 = 0.0125) при проходженнi повз нього iншого елiпсоїда з параметрами b2/a2 = 0.1; c2/a2 = 0.05; a2/a1 = 0.05 на рiзних вiдстанях: 1−∆η/b1 = 3; 2−∆η/b1 = 4; 3−∆η/b1 = 5; 4−∆η/b1 = 8 при зближеннi тiл сили взаємодiї мiж ними рiз- ко зростають. На рис. 21 наведенi аналогiчнi роз- рахунки для “товстого” елiпсоїда зi спiввiдношен- нями осей b/a = 0.25, c/a = 0.0125. Порiвнюю- чи вiдповiднi коефiцiєнти сил та моментiв, можна помiтити, що залежностi для “тонких” i “товсти- х” тiл дещо вiдрiзняються. У випадку “товсти- х” елiпсоїдiв гiдродинамiчнi коефiцiєнти в декiль- ка разiв бiльшi. На рис. 22 наведенi результати чисельних розрахункiв бокових сил та моментiв, що дiють на нерухомий елiпсоїд зi спiввiдношен- ням осей b1/a1 = 0.25, c1/a1 = 0.0125 пiд час проходження повз нього елiпсоїда iнших розмiрiв b2 = a2/10, c2 = a2/20. До характерних особливостей отриманих зале- жностей слiд вiднести iстотне змiщення максиму- му функцiї Cm вздовж осi ∆ξ при змiнi довжини рухомого елiпсоїда. ВИСНОВКИ У роботi побудована чисельна модель гiдроди- намiчної взаємодiї тiл (зокрема суден), що рухаю- ться в рiдинi. Модель грунтується на застосуван- нi методу граничних елементiв. Поверхня тiла та непроникнi дiлянки границi апроксимуються си- стемою трикутних елементiв. На кожному дискре- тному елементi поверхнi розподiленi джерела по- стiйної iнтенсивностi. Значну увагу придiлено по- будовi матрицi, яка описує взаємний вплив три- кутних дискретних елементiв. У роботi проведена апробацiя чисельної моделi. Виконанi розрахунки для потенцiалу ϕ i швидкостi рiдини на поверх- нi рухомого елiпсоїда в необмеженiй рiдинi та по- близу плоскої стiнки. Одержанi данi, якi опису- ють денiвеляцiю поверхнi води при руховi судна (що моделюється еквiвалентним тривiсним елiпсо- їдом) у каналi з трапецевидним перерiзом. Побудо- ванi залежностi коефiцiєнтiв приєднаних мас елi- псоїда вiд його розташування в каналi, параметрiв каналу та елiпсоїда. Розрахованi сили гiдродина- мiчної взаємодiї двох елiпсоїдiв, якi рухаються в безмежнiй рiдинi. Розглянутi випадки однакових та рiзних (за розмiрами) елiпсоїдiв, коли вони ру- хаються в одному напрямку, та коли один елiпсо- їд рухається, а другий нерухомий. Побудована чи- сельна модель та отриманi залежностi дозволяють з бiльшою точнiстю проводити аналiз спецiальних режимiв руху суден, зокрема в рiчках, каналах, поблизу гiдротехнiчних споруд тощо. 1. Бенержи П., Баттерфилд Р. Метод граничных эле- ментов в прикладных науках.– М.: Мир, 1984.– 494 с. 2. Бреббия К., Теллес Ж., Вроубел Л. Методы гра- ничных элементов.– М.: Мир, 1987.– 524 с. 3. Горбань И. Н. Исследование вихревых течений в ограниченных областях. Диссерт. на соиск. уч. ст. к.ф.-м. н.– К.: ИГМ НАНУ, 1993.– 236 с. 4. Горбань В. А. О динамике систем тел в идеальной жидкости.– К.: Гидромеханика, 1983.– 236 с. 5. Короткин А.И. Присоединенные мacсы судна.– Л.: Судостроение, 1986.– 312 с. 6. Костюков А.А. Взаимодействие тел, движущихся в жидкости.– Л.: Судостроение, 1972.– 310 с. 7. Костюков А.А. Теория корабельных волн и волно- вого сопротивления.– Л.: Судпромгиз, 1987.– 123 с. 48 В. О. Горбань, С. В. Маcюк ISSN 1561 -9087 Прикладна гiдромеханiка. 2006. Том 8, N 3. С. 27 – 49 8. Macтушкин Ю.М. Гидродинамическое взаимодей- ствие судов при встречах и обгонах.– Л.: Судостро- ение, 1987.– 311 с. 9. Милн-Томсон Л.М. Теоретическая гидромехани- ка.– М.: Мир, 1964.– 555 с. 10. Павленко Г. Е. Избранные труды.– К.: Наук. дум- ка, 1978.– 495 с. 11. Сретенский Л. Н. Теория волновых движений жидкости.– М.: Наука, 1977.– 815 с. 12. Cohen S., Beck R. Experimental and theoretical hydrodynamic forces on a mathematical model in confined waters // Journal of Ship Research.– 1983.– 27, N 2.– P. . 13. Q.Y. Gui, J.M. Chuang and C.C. Hsiung A discussi- on on irregularities which occur in computation of ship manoeuvring in a restricted waterway using the numerical conformal mapping method // Int. Shipbuild.Progr.– 1992.– 39, N 417.– P. 55-67. 14. Hess J. L., Smith A. M. O. Calculations of nonlifti- ng potential flow about arbitrary three-dimensional bodies // Journal of Ship Research.– 1964.– 8, N 2.– P. 22–44. 15. C.C. Hsiung and Qianyi Gui Computing interacti- on forces and moments on a ship in restricted waterways // Int. Shipbuild. Progr.– 1988.– 35, N 403.– P. 219-254. 16. Korsmeyer F.T., Lee C.-H., Newman J.N. Computation of Ship Interaction forces in restricted waters // Journal of Ship Research.– 1993 .– 37, N 4.– P. 298–306. 17. P.Krishnankutty and K.S. Varyani Force on the mooring lines of a ship due the hydrodynamic interaction effects of a passing ship // Int. Shipbuild. Progr.– 2004.– 51, N 1.– P. 33-57. 18. Kyulevcheliev S., Georgiev S. Experimental observations of ship wavemaking at trans - and supecritical speeds // Euro-Conference HADMAR 2001,Varna,Bulgaria.– 2001.– N .– P. . 19. Newman N. J. Evaluation of the wave-resistance Green function: part 1 - the double integral // Journal of Ship Research.– 1987.– 31, N 2.– P. 79– 90. 20. Newman N. J. Evaluation of the wave-resistance Green function: part 2 – the single integral on the centerplane // Journal of Ship Research.– 1987.– 31, N 3.– P. 145–150. 21. Vantorre M., Laforce E., Verzhbitskaya E. Model tests based formulations of ship-ship interaction forces for simulation purposes // IMSF 28th Annual General Meeting,Genova.– 2001.– N .– P. . 22. Weihs D., Ringel M., Victor M. Aerodynamic interactions between adjacent slender Bodies // AI- AA Journal.– 2006.– 44, N 3.– P. 481-484. 23. Zhi Guo and Allen T. Cnwang Oblique Impact of Two Cylinders in a Uniform Flow // Journal of Ship Research.– 1991.– 35, N 3.– P. 219-229. В. О. Горбань, С. В. Маcюк 49

Чисельне моделювання гiдродинамiчної взаємодiї тiл, що рухаються в рiдинi

Репозитарії

Схожі ресурси