Моделирование внутреннего и внешнего нестационарного взаимодействия корпуса летательного аппарата с жидкостью методом граничных элементов Текст научной статьи по специальности «Математика»

УДК 629.762; 532.5.031

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

© А.В. Плюснин

ОАО «ВПК машиностроения», г. Реутов Московской обл., 143966, Россия

Рассмотрены внутренние (колебания топлива в баках) и внешние (определение присоединенных масс и моментов инерции) задачи нестационарного взаимодействия корпуса ЛА и идеальной несжимаемой жидкости и их формулировки в виде граничных интегральных уравнений. Приведены формулы эффективного решения указанных задач методом граничных элементов применительно к телам вращения и примеры расчетов.

Ключевые слова: потенциальные течения, гидродинамические силы, присоединенные массы и моменты инерции, метод граничных элементов.

Имеются важные технические приложения гидромеханики, в которых математическая модель безвихревых течений идеальной несжимаемой жидкости дает удовлетворительные результаты. Такими, в частности, являются задачи определения нестационарных сил на поверхности ЛА, движущихся в воде с небольшими скоростями [1-6], а также задачи определения возмущающих сил на корпус жидкостного ЛА, вызванных малыми колебаниями топлива в баках [7, 8].

В указанных задачах в области QeК3, занятой жидкостью, поле скоростей v удовлетворяет одновременно условиям несжимаемости и потенциальности течения:

где ф — потенциал скоростей жидкости. Следовательно, в области О потенциал ф должен удовлетворять уравнению Лапласа

Дополняя уравнение (2) подходящими условиями на границе области О, получаем классические задачи теории потенциала [1-10].

В рассматриваемых задачах для вычисления сил, моментов сил, действующих на конструкцию со стороны жидкости, и инерционных характеристик жидкости (присоединенных масс и моментов инерции) требуется знать решение (потенциал ф и его нормальную производную) только на границе расчетной области. Поскольку применительно к задачам для уравнения Лапласа имеются эквивалентные формулировки в виде граничных интегральных уравнений (ГИУ) [9], привлека-

div v = 0, v = grad ф,

(1)

Дф = 0.

(2)

тельным методом решения таких задач является метод граничных элементов (МГЭ) [4, 5,11-16].

Решение задач для уравнения Лапласа методом граничных элементов. Пусть рассматривается некоторая задача для уравнения Лапласа в области Ос К3, границей которой является поверхность £. Условимся (рис. 1), что единичная нормаль п (М) к поверхности £ в точке М е £ будет внешней по отношению к области О для внутренних задач («¿О) и внутренней по отношению к области О для внешних задач («е О).

Х\,г а

Х\=2

Х2 = г сое 9 х3 = г вт 0

Рис. 1. Геометрия внешней и внутренней задач: а — поверхность ЛА; б — бак с топливом

В [9] решения граничных задач Неймана, соответствующих зада-

нию на поверхности £ нормальной производной

дф дп

ставлены в виде потенциалов простого слоя

ф<

(М ) = {—)(Р )& (Р ).

£ ГМР

= I), пред-

(3)

Решения граничных задач Дирихле, соответствующих заданию на поверхности £ самой функции ф|^ = /(°), представлены в виде потенциалов двойного слоя (для внешних задач к потенциалу двой-

£

ного слоя добавляется еще источник, центр которого не принадлежит множеству О и £ )

где М — произвольная точка области О; Р — точка интегрирования по поверхности 5; гМР — расстояние между точками М и Р. В дальнейшем точку М условимся называть точкой наблюдения, точку Р — точкой интегрирования (по поверхности 5 ).

Функции (3) и (4) определены на всем пространстве. Потенциал простого слоя — непрерывная функция, а потенциал двойного слоя при переходе через поверхность 5 изменяется скачком, причем его предельные значения и прямое значение (случай М е 5) связаны между собой соотношениями типа формул Сохоцкого [9].

Требование удовлетворения граничным условиям приводит к граничным интегральным уравнениям:

' (М)±(г^V'1 (р)"5(р)=±1 (М) <5>

ГМР у

' 5

для задач Неймана и

'" 1 Г д ( 1 ^

x ..........

v гмр у

2 л 5 дп (М) на и

2л'

для задач Дирихле [9], причем верхний знак относится к внутренним задачам, нижний знак — к внешним.

Из перечисленных четырех задач (рассматриваемых в пространстве К3) внутренняя задача Неймана либо не имеет решения, либо ее решение определено с точностью до константы. Соответствующее

однородное ГИУ (5) имеет одну собственную функцию ^)(Р), определяющую потенциал Робэна,

Ф<* )(М ) = { )(р )"5 (Р ),

5 ГМР

где ф(Л)(М)=1 при Ме Ои5. Союзным к ГИУ (5) для внутренней

задачи Неймана является ГИУ (6) для внешней задачи Дирихле. Собственная функция соответствующего однородного уравнения есть

д ( 1 ^

) = 1 (при *'") = 1 имеем [—-—- -I"5 (Р) = - 2л ). Это дает не-

5 дп (Р)( гмр У

обходимое и достаточное условие разрешимости внутренней задачи Неймана [9]:

J f( N V *) dS SJ f (N ) dS = 0. (7)

S S

Остальные три задачи разрешимы единственным образом при любой правой части. При этом ГИУ для внешней задачи Неймана и внутренней задачи Дирихле союзны между собой, не имеют собственных функций, а потому также разрешимы единственным образом при любой правой части [9]. В общем случае решение внешней задачи Неймана стремится на бесконечности к нулю обратно пропорционально расстоянию от начала координат. При выполнении условия (7), которое в данном случае не является условием разрешимости, решение стремится на бесконечности к нулю обратно пропорционально квадрату расстояния или быстрее.

Несколько другой подход к решению граничных задач для уравнения Лапласа основан на формуле Грина [9-15] (M е S):

f i ^

/ т^Ч то / т~») _

ф(М )±— ΗW )dS (P ) ' 2л J ôn(P)tMГ

_ rj—ô^(p )dS (р )=о. (g)

2лSmp ôn(P)W W w

Для задач Неймана имеем отсюда ГИУ

ф(М)±— Г-^Шф(р)^ (р) = П ' 2л^ дп (Р^ М Г

=^ ^ I№ (Р ). (9)

2Л ^ ГМР

В формулах (8), (9) верхний знак соответствует внутренним задачам, нижний — внешним задачам. Уравнение (9) для внешней задачи аналогично ГИУ (6) для внутренней задачи Дирихле, а потому всегда единственным образом разрешимо. Уравнение (9) для внутренней задачи аналогично ГИУ (6) для внешней задачи Дирихле. Поскольку в силу свойств потенциала Робэна

}(р И (Р ){—I](Р (Р )=

^ игРР

=| I ^ >(Р (Р ) -±- }(Р (Р ) =

5 8ГРР

= | I(* }(Р )Ф( *}(Р ) (Р ) = { /(*)(Р ) (Р ),

его условие разрешимости совпадает с условием (7).

Для того чтобы решить ГИУ с достаточной точностью, необходимо обеспечить подходящую аппроксимацию поверхности 5 объединением граничных элементов (ГЭ), аппроксимацию искомых и известных функций на ГЭ и численное интегрирование по ГЭ потенциалов простого и двойного слоя.

Многие ЛА и топливные баки жидкостных ЛА являются телами вращения. Рассмотрим контур С, получаемый пересечением поверхности 5 произвольным меридиональным сечением. Интегральное уравнение каждой из рассматриваемых здесь задач, заданное на поверхности 5, может быть преобразовано в интегральное уравнение, заданное на контуре С. Хотя расчетные формулы усложняются, имеется выигрыш в размерности задачи. В частности, более наглядными получаются результаты расчетов. Ограничимся этим случаем.

Рассмотрим прямоугольную Ох1 х2х3 и цилиндрическую Oг0z системы координат (см. рис. 1). Контур С выберем в меридиональном сечении 0 = 0. Ориентацию контура С выберем так, чтобы нормаль п к поверхности 5 при положительном направлении обхода контура оставалась справа. Начальную точку контура С совмещаем с началом координат О. Контур С будем задавать параметрическими уравнениями

г=Гс ( 5); z = zc ( 5), (10)

где 5 — длина дуги, отсчитываемая от точки О. Для компонент нормали имеем тогда соотношения

п

= zC( 5); Щ =-гС( 5). (11)

Решения ф рассматриваемых задач, как будет показано далее, имеют специальный вид, так что

Ф (г, ^ 0) - у(0) (г, ^ 0) = Ф(0) (г, z) 0 (0);

(12)

^(г, ^ 0)^«(г, ^ 0) = Ф«(г, z)0(0), дп

где 0(0) — известные функции простого вида. Поэтому достаточно

располагать точки наблюдения М на контуре С. Их положение будет определяться значением длины дуги 5, где 5 е[0; Ьс ], ЬС — полная длина контура. Точки Р интегрирования по поверхности 5 будем задавать угловой координатой $е[0; 2л] соответствующего меридионального сечения и длиной дуги а вдоль контура, являющегося следом поверхности 5 в этом сечении. Следовательно, цилиндрическими координатами этих точек будут значения р = Гс (а), £ = zc (а) и О.

Элемент поверхности в точке Р можно определить как dS = = Гс (a)d9dа. Тогда потенциалы простого и двойного слоя в ГИУ можно записать следующим образом:

где

|К(^5\Ы, Р)М/(°'1)(Р^(Р) =

S

2л

{ф(°Д)(а)гс (а)а{ К(й?'5а, 9)0(9)9 =

с °

= | 5)( а)Ф(°д)(а)а, (13)

Ф(°'1)(а) = Ф(°'1)(гс (а), ¿с (а)); К«(М,Р)=

ГМР

,(d „Ч д Г 1 ^

К(а '(М, Р ) =

дп (Р )

ГМР у

К^''>(, а, 9)=К(*5)((.), ¿с (.), гс (а), ¿с (а), КV,Р);

2л

5)( а) = гс (а) | К(,5)(, а, 9)0(9)9.

°

Основное интегральное соотношение (9) в результате преобразования к контурному интегрированию принимает следующий вид:

0 (°) Ф(0) (5)+— Г К ^) ( а) ф(0) (а)d а + 2л с

+ — Га)ф(1)(а)dа = °. (14)

2л

с

^, ^

Явные формулы для выражений 0(0) и К( ,будут приведены ниже.

К уравнениям вида (14) МГЭ применяется следующим образом. Контур с разбиваем точками деления а° =°<а1 <...<аN_1 <аN = Ьс на N гладких дуг, которые и будут ГЭ. Точки наблюдения М - выберем в центрах ГЭ: а-+у2 = °,5(а- + а-+1), где °< ]<N _ 1. Выписываем соотношения типа (14) для каждой точки наблюдения:

Л N-1

®(0)Ф(0)(оу+1,2 )± ±£1<?Ф(°>К+1/2 ) +

2- к = 0

. N-1

+ 2-1^к^К+1/2 ) = 0, (15)

2- к = 0

где

N-1

|Кк)(ау+1/2, а)Ф(0'1)(а)^а. £ I(кк)Ф(0'1)(а)

С к = 0

ак+1

I(кк )= ¡К (,к )(а ,+1/2, а) а.

(16)

Соотношения (15) образуют систему N линейных алгебраических уравнений (с плотной несимметричной матрицей), решаемую методом исключения Гаусса [17, 18]. Коэффициенты (16) вычисляем по квадратурным формулам, например по формулам Гаусса наивысшего порядка точности [12, 19].

Расчет присоединенных масс и моментов инерции ЛА. Пусть тело конечных размеров движется по произвольному закону в безграничной идеальной несжимаемой жидкости. Свяжем с телом прямоугольную систему координат ОХ1Х2Х3, а ее начало выберем в качестве полюса (рис. 1, а). Скорость полюса \0 и угловую скорость © вращения точек тела вокруг полюса разложим по ортам е, связанных осей: У0 (X)е,, © = ш, (X) е,. Тогда потенциал скоростей безотрывного течения, вызванного движением тела, имеет вид [1, 2, 4, 20, 21]

3

ф(хь Х2, Хз, X)=Х[ф} Х2, Хз) (X) + ф]+3 Х2, Х3)©] (X)],

] = 1

где потенциалы фj и фj+3 удовлетворяют краевым задачам (] = 1, 2, 3):

Аф, = 0|0;

5ф j

—=п •е з

дп

аф j+3=01 □;

дф j+3

дп

(17)

= (г • п )• е,

к

5

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

остальные три потенциала описывают течение жидкости, вызванное вращением тела с единичной угловой скоростью вокруг оси Oxj. Поскольку JndS = J grad(1)dS = 0, Jr• ndS = - J rotrdS = 0, единич-

S R3\Q S R3\Q

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

Основные величины, характеризующие возмущенное движение жидкости в целом как механической системы, в частности, главный вектор и главный момент действующих на ЛА гидродинамических сил, кинетическая энергия жидкости, выражаются через параметры движения ЛА и коэффициенты

hj = -Рж [фг^ = hji, (18)

S an

называемые коэффициентами присоединенных масс [1, 2, 20].

Для тела вращения (см. рис. 1) орты цилиндрической системы координат выражаются через векторы ei следующим образом:

e r = cos 9 e2 + sin 9 e3; e9 = -sin 9 e2 + cos 9 e3;

e z = e1.

Кроме того, справедливы следующие выражения для радиуса-вектора r произвольной точки пространства и для вектора n нормали к поверхности тела:

r = r er + z ez; n = nrer + nzez; r • n = (znr - rnz ) ee.

Используя данные формулы, получим следующие выражения для правых частей краевых задач (17):

аф1=n = n • e1 =(nre r + nze z )e z = nz =-r¡C(s); on

аФ2,3 / ч fcos 91 , /4ícos 91

__=n2,3 = n• e2,3 =(nrer + nzez)e2,3 = nr1 . Af = zc(s) . A;

On [ sin 9J [ sin 9J

°Ф4

on

=(r • n )• e1 =(znr - rnz )ee • e z = 0;

дф5,

дп

- = (r • n)• e2,з =(znr -rnz)e0 • e2,3 =

= {znr - rnz )

-ГН ** ми

Из вида правых частей и из выражения оператора Лапласа в цилиндрических координатах

Аф j =

д2фj +1 diL+± д!ф> i д2фj

дг2 r дг r2 д02 dz2

(19)

вытекают следующие выражения для единичных потенциалов:

>1 (r, 0, z) = Oi (r, z); ф2 (r, 0, z) = Ф2 (r, z)cos 0; фз (r, 0, z) = Ф2 (r, z) sin 0;

ф4 (r, 0, z) = 0; ф5(r, 0, z) = -Ф6(r, z)sin0;

ф6 (r, 0, z) = Ф6 (r, z) cos 0.

Подстановка полученных выражений в формулы (18) дает следующее обобщенное выражение:

2 л Lc

X¡j =A,ji =-Рж j (sin 0)a (cos 0)Pj d0 J fij (s)ds, 0 0

где показатели степени a i и рi принимают значения либо нуль, либо единица. Отсюда легко видеть, что лишь следующие коэффициенты присоединенных масс будут отличны от нуля: Lc

^11 = 2лРж J rC(s )Ф1 (s )ds;

0

Lc

^22 = ^33=- лРж J zC(s)Ф 2 (s)ds;

0 Lc

^26 =-^35 =-Лрж J [zCzC(s) + rCrC (s)] Ф2 (s) d

Lc 0

= -лР ж J zC(s)Ф6 (s)ds;

0 Lc

^66 = ^55 = - лРж J [zCzC (s) + rCrC (s)] Ф6 (s)ds

S =

(20)

J (zC2 irC2)Ф6 (s)ds.

2

Итак, для вычисления коэффициентов (20) надо решить краевые задачи для потенциалов фь ф2 и ф6. В формулах (12)-(14) для потенциала ф1 имеем 0(0) = 1, а для потенциалов ф2 и ф6 имеем 0(0) = cos 0. Интегрирование по угловой координате О в формуле

(13) сводится к вычислению эллиптических интегралов [19, 22], для которых известны простые аппроксимации [23]. Соответственно этому имеем:

2л

K (s)

ф1

1

(а,+1/2, а)=rc (а)/ — 0|ф1 (9)d9 =

0 rMP

= 2 л rc (а) КППС (rC (а j+1/2 ), zC (а j+1/2 ), rC (а), zC (а)) ;

1(+./2> а)=rc(а)^(M)0L(9)d9 =

"2 л rC (а)КПДС (rC (j+1/2 ), zC (аj+1/2 ), rC (а) zc (а) rc (а) zc (а));

K (s)

2л

ф2, фб

(а j+1/2' с) = rc (а)/ 0|ф2, фб (9)d 9 =

0 MP

= 2 л rc (а)ЛППС (rC (аj+1/2 )' zC (аj+1/2 )' rC (а)' zC (а));

K(d) (а j+1/2, а) = rc (а)? —^-i—Id (9)d 9=

V j+ ' ^ dn (PrMP J |ф2,^ '

ф2, фб

= 2 л rc (а)ЛПДС (rc (аj+1/2 ) zc (аj+1/2 ), rc (а) zc (а), rc (а), zc (а)),

где [4, 15, 19, 22]

Kппс (rM, zm, r, z) = 2rA(rM, zM,r,z)к(к2);

KПДС (rM, zM, r, z, r\ z') =

= (rM, ZM, r, z, r', z')(k2)-z'[k(к2)-Е(к2)}a(rM, zm, r, z); Л ппс (rM, zm , r, z) = 2 rA ((, Zm , r, z) 2D (к2 )-к(к2) ;

ЛПДС (rM, zM , r, z, r', z') = = 2 [Г1 (rM, ZM, r, Z r', z')— Г 2 (rM, ZM, r, z, r', Z')]A (rM, ZM, r, z);

Г1 (м, 2м,г, ^ г , 2 --2-Е(к );

(м - г) + (2м - 2)

г2 ((, 2м, г, г, г' , г' ) =

' 2 В (гм, 2м, г, 2, г \ 2 ')- гмг2 ' 0 к ).

(гм + г )2 +(м - 2 )2 1

А (гм, 2м, г, 2)=-

В (гм, 2м, г, 2, г' , 2 ') =

к2 =

-\/(гм + г)2 + (м - 2)2

' ' ) = 2г[г '(гм -г)-г'(2м -2)] .

(гм - г) +(м - 2) 4 гмг

2 2 (гм + г) + (2м - 2)

-/ 2 7 1

к(к2 )=/ ^Ф ;

о ->у1 - к2б1п2 ф

-/ 2 --

Е(к2)= | -к^т2фйф;

о

0(к2) в1и2ФаФ к(к2)-Е(к2) 1 ' /о Ф- к281и2 ф к2 ■

Для вычисления значений Ф12 6 (о j+1/2 ) составим СЛАУ

Ф 1,2,6 (+1/2)-Т1- Е ■ Ф1,2,6 ((+1/2)

2 — к = 0 Ф1,2,6

N-1

1 N-1 (\

=__^ V г )

= ^к

к = о

Фи

- гс

г

2с

2с2с+ Сс

(Ок+1/2 ),

где 0 < ■ < N -1 и коэффициенты

т(а, л

Чк

Ф1,2,6

I К

(а,.

Ф1,:

(ч-^ О)О

вычисляем по приведенным выше выражениям (с использованием стандартных квадратурных формул и аппроксимаций для эллиптических интегралов).

к

В качестве простого примера найдем потенциалы обтекания и коэффициенты присоединенных масс тела вращения, состоящего из передней полусферической части радиуса Я и цилиндрической части длины Ь (см. рис. 1, а). Параметрические уравнения такой поверхности задаются соотношениями, представленными в табл. 1. Результаты расчета функций Ф1 2 6 в зависимости от длины дуги приведены на рис. 2

(распределения единичных потенциалов ф1,2,6 при 0 = 0). Расчетные

значения коэффициентов присоединенных масс приведены в табл. 2 в сравнении с вычислениями по методу плоских сечений (МПС) [4, 20].

Таблица 1

Параметрические уравнения контура ЛА (Ьс = (1 + 0,5 п) Я + Ь)

s rc (s) zc ( s ) rC (s) zC( s )

0 < s < R s 0 1 0

R < s < R+L R s - R 0 1

R + L < s < LC R ( s-R - L A R cos I-I I R J R . (s-R -LA R sin I-I I R J . ( s - R - L A - sin I-I I R J ( s - R - L A cos I-I I R J

Расстояние вдоль контура тела

Рис. 2. Изменения единичных потенциалов вдоль образующей тела вращения (от хвоста к носку, Я = 1, Ь = 1, 100 ГЭ):

1 — потенциал ф1; 2 — ф2; 3 — ф6

Таблица 2

Пример расчета коэффициентов присоединенных масс тела вращения

Размеры тела R=1, L=1 R=1, L=9

Метод расчета МГЭ(100 ГЭ) МПС МГЭ (200 ГЭ) МПС

^ и/ Р ж 2,924 - 3,326 -

^ 22 / Р ж зз/ Р ж 2,826 5,236 26,515 30,37

^ 26 / Р ж 35 / Р ж 2,190 4,886 127,94 149,7

^ бб/ Р ж 55/ Р ж 1,999 5,835 776,88 975,2

Расчет гидродинамических коэффициентов, характеризующих силовые и инерционные свойства жидкости при поперечных колебаниях в топливных баках ЛА. В задаче о малых поперечных колебаниях топлива в баках жидкостных ЛА (см. рис. 1, б) частью общей границы S жидкого объема О служит твердая стенка топливного бака (обозначим ее 5^), а остальная часть границы — это свободная поверхность жидкости (обозначим ее £п). Основной математической моделью в инженерных расчетах является модель идеальной несжимаемой жидкости.

В полете бак совершает (благодаря работе системы управления) малое возмущенное движение относительно некоторого заданного (программного) движения. Строго говоря, на каждую жидкую частицу действует сила земного притяжения, на частицы зеркала топлива действует давление наддува бака, на частицы, контактирующие со стенками бака, действует давление их взаимодействия. Но характерное время возмущений движения бака на несколько порядков больше времени распространения волн давления в объеме топлива. Стенки бака не препятствуют движению жидкости в поле тяготения (условно говоря, бак свободно падает вместе с топливом). Поэтому при отсутствии возмущений, т. е. при программном движении, достаточно рассматривать жидкость в предельном состоянии покоя относительно бака и считать, что ее свободная поверхность, так называемое зеркало топлива, перпендикулярна оси бака. Этому состоянию соответствует плотность массовых сил w = -vXl g е1, где g — ускорение силы тяжести; Vх — продольная перегрузка ЛА. Поперечное возмущенное движение бака приводит к колебаниям топлива в баке относительно этого положения покоя, т. е. в том же поле массовых сил.

За исключением переходных процессов, вызванных внезапными возмущениями, и явлений резонанса возмущенное движение топлива в баке сопровождается малыми перемещениями и частиц относительно положения равновесия. Тогда движение жидкости удобнее описывать потенциалом перемещений Охоцимского х [7, 8], т. е. вместо соотношений (1), (2) использовать соотношения

Шу и = 0; и = grad х; Ах = 0, (21)

причем V «Эи/ д1 и ф«Эх/Э.

С программным движением бака будем связывать семейство неподвижных («абсолютных») систем координат О *х*х2 х3. Для связанной системы координат сохраняем прежнее обозначение Ох1 х2 х3. При отсутствии возмущений система Ох1 х2х3 совпадает с актуальной (соответствующей данному моменту времени) абсолютной системой.

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

и=иО + их г,

где и^ =и2 (г)е2 + и3 (г)е3 — смещение полюса О относительно его

невозмущенного положения О*; и = и(г)е2 + и(г)ез — вектор поворота бака на малый угол; г — радиус-вектор точек бака в связанных осях. Тогда условие непротекания на смоченной поверхности бака принимает вид:

Яу

—=и х п=и2п2 + и3п3 + и2 (х3п1 - х1п3) + и3 (х1п2 - х2п1). дп

На поверхности £п, аналогично классическому подходу, принятому в теории малых гравитационных волн на поверхности жидкости [25], ставятся динамическое условие, соответствующее известному давлению в газовой подушке над зеркалом топлива, и кинематическое условие, определяющее возмущенное положение зеркала. Применительно к рассматриваемому случаю основное соотношение, получаемое в результате линеаризации указанных граничных условий и исключения высоты ^ = ду/дх1= 0 подъема зеркала топлива над

невозмущенным уровнем, имеет вид [7, 8, 25]

д2у ду

—г+v хЕ—

дг2 х дх1

= 0.

Для малых возмущенных движений бака и топлива абсолютная и относительная производные по времени [1, 2, 26] отличаются на величины второго порядка малости, поэтому их различием будем пренебрегать.

Таким образом, требуется найти решение следующей внутренней краевой задачи

АУ=0о;

ду/дп = и2п2 + изпз + «2 (хзп1 - х1пз ) + из (х1п2 - х2п1 )) ; (22)

д 2у/дг 2 + v х. Е д%1дх1

= 0.

В отличие от задачи обтекания ЛА, которая привела к краевым задачам (17) для единичных потенциалов, зависящих только от внешней геометрии тела, решение задачи (22) зависит не только от

5П

внутренней геометрии бака, но и от текущего уровня топлива H, снижающегося по мере выработки, а также от «предыстории» колебаний зеркала топлива.

При фиксированном уровне топлива соответствующая задаче (22) однородная задача (отличие состоит в условии 8%/8п|„ = 0) является

задачей на собственные колебания. Имеется счетное множество собственных частот колебаний ю2 и соответствующих им собственных

функций уп (хъ x2, Х3), которые находят, решая задачу на собственные значения:

Ау=0q ;

8v/8n = 0 Si; (23)

dy¡dxi -v'xlg-1ш2у s = 0.

sII

Функции уп нормируют условием 8у n/ dx1 ^„м) =1 в некоторой фиксированной точке

Q (норм)е S n SIb и вместе с функцией у0 Q=1

они образуют полную ортогональную систему в области Q [7, 8].

Общее решение однородной задачи (22) и соответствующая ему форма свободной поверхности топлива имеют вид:

Х(одн )(Х1, Х2, Хз, t) = ^( An1 cos Юп/ + A(2) sin Юп/)уп X2, X3 );

n = 1

С (2, X3, t ) = Эх(одн V^X1 =¿( аП1 cos Юп/ + 42) sin Юп/)^п (Х2, Х3 ),

Sl1 п = 1

где ^п (Х2, Х3)=8уп/8X1| = у~^-1ш2 уп (H, Х2, Х3), — собственные формы колебаний зеркала. Функции ¥п удовлетворяют условию ¥ пи(норм) = 1 и вместе с функцией ¥0 L = 1 образуют полную ортогональную систему на поверхности Sn [7, 8].

Разложение произвольной функции f, определенной на поверхности SII, в обобщенный ряд Фурье по системе ¥п имеет вид

ю С ¥ f = V^ ^п х п

N2 '

п=1 п

где Сп = j f¥пdS, N2n = ¡¥2ndS [7, 8].

Sii

Sii

Положим Б = Б, и Б,, и рассмотрим краевые задачи (17) как внутренние задачи для потенциала перемещений (заведомо разрешимые с точностью до постоянной). Для индексов у = 2,3 соответствующими

«единичными» потенциалами являются %2 3 = х2 3 и %у+3 = Ту (решения задачи Жуковского) [7, 8]. Функция

%

(неодн)=

£ \Uj()^+0(ОТ]+£^п() (24)

У=2,3

п = 1

является гармонической и удовлетворяет неоднородному граничному условию задачи (22). Подставим ее в граничное условие этой задачи на поверхности Б11, где х1 = Н; п = е1; д/дп = д/дХ^ При этом учтем, что

дТ у

дх1

= дТ j

дп

= (г • п)-е ,-|„ =

(еу-г)-п|

-Х1П3 + Х3П1 (у =2)"1=| Х3 (у = 2)1

Х1П2 -Х2П1 (у = 3) [-Х2 (у = 3)/'

и воспользуемся разложениями

м Л( • )ш о Х. = У п ^п . РжХу - —2 ; п = 1 —п

м л(у )щ о Т = У Оп п

у ¿а N2 : п = 1 —п

где

Ж )=0ж|

Б,

Ху¥пёБ; 4о2 = Ож | Т• ^пйБ.

Б,,

(25)

Имеем:

• (неодн)

£ Ну(t)xj +(3У(0^1 + £

'П у=2,3

-ьг-2 V- Гл(;и. ^) + .(• )

У Х

Б 2 "п

- 1 Ю

п = 1 шп

п $„ (t) =

= £Ьж—-2 £ Ки(0+••• (t)

п = 1 | у = 2,3

у х1 £

2п

(t )к п;

_д%_

дХ1

/ \ дТ •

= £° • ( ) аТ

у = 2,3 °Х1

+ £^ п*п (t ) =

Б„ п = 1

О (t) Х3 - 03 (t) Х2 + £ ^ ^п () =

п=1

и

11

и

м

м

-- Е кМ 2 [^2 (г) -42)из (г)]+^ (г )К;

п = 1

Е ^ и + (/)

У = 2,3

, Р Ж У яМ,

I(У

+ У я

"^3)»2 «)-^П2)из (/)

юп

2 - ( +юП^п ) = 0.

Отсюда получаем обыкновенные дифференциальные уравнения для функций лп (г):

Дп (+®2Лп)=- Е и (г(г)

У = 2,3 Ь

"У Я

(г )-42 >»з (г)'

(26)

где

Д п

РжУ

Юп

Таким образом, если функции лп (г) удовлетворяют уравнениям (26), то функция х( неодн) является решением задачи (22).

Величины дп,

"У1

и "¿п называют соответственно обобщенными массами и коэффициентами гидродинамических силы и момента, соответствующих п-й форме колебаний топлива в баке. Через эти коэффициенты и присоединенные моменты инерции

т(У)=Р ГТ ^

¿П =ржJ1У

дп

сК

(27)

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

Таким образом, основные характеристики воздействия на бак со стороны колеблющегося топлива находят из решения (для каждого уровня заполнения бака) задачи на собственные значения (2з) и внутренних задач Неймана (17) для потенциалов У2,з. Первостепенное значение имеет первый тон антисимметричных колебаний, при которых зеркало топлива сохраняет плоскую форму [8]. Остальные формы колебаний часто не учитывают.

Рассмотрим бак в форме тела вращения. Тогда потенциалы Т2,3, соответствующие потенциалам ф56 в формулах (19), можно выбрать в виде

, ч , ч fsin0] Т2,з(r, 0, z) = + T(r, z) L. (28)

[cos0J

Строго говоря, следует полагать Тj = Т j + Cj, где Tj есть правая часть формулы (28), а Cj — произвольные константы. Нетрудно видеть, что появятся слагаемые Cj и j (t) в решении (24) и Cj ü j (t) во втором граничном условии задачи (22). Через функции Tj определяются

присоединенные моменты инерции (27) и коэффициенты (вторая из формул (25)). Вектор перемещений u и, в частности, форма свобод-

ной границы g =

С = Эх(неодн )/дх1

от произвольных констант не зависят.

Из выполнения условий разрешимости для внутренних задач (17) следует независимость присоединенных моментов инерции от Су:

г~ дТу п~ чдТ, ,дТ, , дТу

Г Т = Г (Т 1 + С,-)ёБ - С,- ёБ = Г Т у—J-dБ. Б у дп Б дп уБ дп Б у дп

Ввиду равенств

Г ¥пёБ = Г^Б = Г^шёБ = Г^Б - Г^шёБ =

Б Б дХ1 Б дп Б дп Б дп

-A^dS = ÍAyndS = 0;

i дп ¿

\ Т у ¥ п ёБ = \ Т у ¥ п ёБ - Су | ¥ пёБ = \ Т у ¥ п ёБ

Б11 Б11 БИ Б11

определение коэффициентов ^¿1 и, следовательно, функций 8п () уравнениями (26) также не содержит произвола.

Наконец, появление слагаемых вида Су0у (t) в условии на свободной границе просто подчеркивает тот факт, что для внутренних задач потенциал скоростей/перемещений остается определенным (из вторых по счету уравнений (1) и (21)) лишь с точностью до произвольной функции времени. При постановке внешних задач (в М3) этот произвол обычно устраняется условием стремления потенциала к нулю на бесконечности.

Sii

С помощью МГЭ функция Т( r, z) находится на контуре C меридионального сечения области Q аналогично функции Ф6 (r, z), только в формулах (14) и (15) надо брать верхние знаки. В результате имеем СЛАУ (о < j < N -1)

1 N-1 . . i N-1

+1/2 ) + j j+1/2 ) = j (zczC+ rcrC )) .

k=о J ф, k=о

Заметим, что неоднозначность решения исходной граничной задачи и соответствующего ГИУ была устранена выбором решения в виде (28). Присоединенный момент инерции вычисляем, используя аналогию между выражениями (18) и (27). Тогда по последней из формул (2о) имеем

Lc

jc (zc2+c2) )&

2 о

При решении МГЭ задачи (23) на собственные значения ограничимся расчетом антисимметричных тонов. В этом случае собственные функции уп имеют вид

, ч , ч f sin 0] уп (r, 0, z) = y(r, z) J, (29)

[cos 0J

где sin 0 соответствует колебаниям в плоскости Ox1x3, а cos 0 соответствует колебаниям в плоскости Ox1 x2. Нас интересует нахождение функций у на контуре c. Рассмотрение второго случая оставляет точку наблюдения M в сечении 0 = 0 и не требует вывода новых формул, связанных с интегрированием по угловой координате.

Соответственно внутренней задаче, выберем в ГИУ (14) верхние

-л

знаки, а также положим 0(0) = cos 0. На поверхности Sj —п=0, а

дп

на поверхности -Уп. Следовательно, согласно обо-

дп Vxxg

значениям (12), имеем Ф(0) (а)=у(а) на всем контуре c, Ф(1)(а)=0

2

на части cj = c о Sj и Ф(1)(а) = Юп у(а) на остальной его части

У X1g

cjj = c о Sjj. Тогда y(s) + _L jK(d) (s,o)y(a)do--— j K(s) (s,o)y(a)do = 0.

/ TT " rr\ / 7T \ f ГГ " rr\

c Ф6 2 X gc

фб

В результате дискретизации МГЭ получаем однородную СЛАУ

1 N _1 2 N _1 ij )+ ¿I j г (о j+1/2)-т-^— I j r( j+1/2 )о, (30)

фб 2Л к=N

2л к = 0

фб

где N — число ГЭ на части С контура С. Задача определения решения СЛАУ (30) есть обобщенная задача на собственные значения и решается стандартными численными методами линейной алгебры [16, 17].

Рассмотрим в качестве примера расчет гидродинамических характеристик колебаний топлива в баке со сферическим днищем радиусом Я, боковой цилиндрической поверхностью того же радиуса и уровнем топлива Н (см. рис. 1, б). Параметрические уравнения контура С задаются соотношениями, представленными в табл. 3. Результаты расчета функций у (для трех тонов антисимметричных колебаний) и Т приведены на рис. 3. Расчетные значения гидродинамических коэффициентов даны в табл. 4. Отметим, что эти коэффициенты мало отличаются от точных значений для цилиндрического бака [8], рассчитанных при тех же условиях Я = 1, Н = 2.

Таблица 3

Параметрические уравнения контура С (Ьс = (0,5 л Я + Н)

s re (s) Ze (s ) rC( s )

%R 0 < s <- 2 s R sin— R R ('- cos Rj s cos— R

—< s <í*-0 R + H 2 12 J R s-(! - 1JR 0

^ y-1j R + H < s < LC Le - s H -1

0,6 0,4 0,2 0 -0,2 -0,4 -0,6 -0,8 -1,0 -1,2

0 0,4 0,8 1,2 1,6 2,0 2,4 2,8 3,2 3,6 Расстояние вдоль контура расчетной тела

Рис. 3. Изменения функций у и Т вдоль контура меридионального сечения области, занятой топливом (Я = 1; Н = 2; 180 ГЭ): 1 — у, первый тон; 2 — у, второй тон; 3 — у, третий тон; 4 — точка нормировки;

5 — Т

Таблица 4

Коэффициенты, характеризующие колебания топлива по трем антисимметричным тонам, и присоединенный момент инерции для бака осесимметричной формы (R = 1, H = 2, 180 ГЭ)

Тон (n) Гидродинамические коэффициенты

Ц n ^ n ^On < g-1®n J

1 0,6063 0,9284 1,437 1,834 7,141

2 0,2879 0,1114 0,2019 5,335

3 0,1853 0,04379 -0,08234 8,548

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

ЛИТЕРАТУРА

[1] Седов Л.И. Механика сплошной среды. Т. 2. Санкт-Петербург, Лань, 2004, 560 с.

[2] Лойцянский Л.Г. Механика жидкости и газа. Москва, Наука, 1987, 840 с.

[3] Batchelor G.K. An Introduction to Fluid Dynamics. Cambridge, University Press, 2000, 615 p.

[4] Дегтярь В.Г., Пегов В.И. Гидродинамика подводного старта ракет. Москва, Машиностроение; Машиностроение-Полет, 2009, 448 с.

[5] Якимов Ю.Л. Сборник статей. Москва, Изд-во Московского университета, 2013, 536 с.

[6] Короткин А.И. Присоединенные массы судостроительных конструкций: справочник. Санкт-Петербург, Мор Вест, 2007, 448 с.

[7] Рабинович Б.И. Введение в динамику ракет-носителей космических аппаратов. Москва, Машиностроение, 1975, 416 с.

[8] Микишев Г.Н. Экспериментальные методы в динамике космических аппаратов. Москва, Машиностроение, 1978, 248 с.

[9] Соболев С.Л. Уравнения математической физики. Москва, Наука, 1966, 444 с.

[10] Димитриенко Ю.И. Механика сплошной среды. В 4 т. Т. 1. Тензорный анализ. Москва, Изд-во МГТУ им. Н.Э. Баумана, 2011, 463 с.

[11] Cruse T.A., Rizzo F.J., eds. Boundary integral equation method: computational applications in applied mechanics! Ed. ву New York, ASME, 1975, 390 p.

[12] Baneijee P.K., Butterfield R. Boundary Element Methods in Engineering Science. London, McGraw-Hill Book Company (UK) Limited, 1981, 494 p.

[13] Brebbia C.A., Telles J.C.F., Wrobel L.C. Boundary Element Techniques. Berlin, Springer-Verlag, 1984, 523 p.

[14] Cebeci T., Shao J.P., Kafyeke F., Laurendeau E. Computational Fluid Dynamics for Engineers. Long Beach: Horizons Publishing Inc., 2005, 396 p.

[15] Власова Е.А., Зарубин В.С., Кувыркин Г.Н. Приближенные методы математической физики. Москва, Изд-во МГТУ им. Н.Э. Баумана, 2001, 700 с.

[16] Зарубин В.С., Кувыркин Г.Н. Особенности математического моделирования технических устройств. Математическое моделирование и численные методы, 2014, № 1, с. 5-18.

[17] Плюснин А.В. Расчет нестационарной гидродинамической нагрузки на крышку пускового устройства при подводном выбросе ЛА избыточным давлением газа. Инженерный журнал: наука и инновации, 2013, №7(19). URL: http://engjournal.ru/ catalog/mathmodel/aero/847.html

[18] Demmel J.W. Applied Numerical Linear Algebra. Philadelphia, SIAM, 1997, 430 p.

[19] Калиткин Н.Н. Численные методы. Москва, Наука, 1978, 512 с.

[20] Риман И.С., Крепс Р.Л. Присоединенные массы тел различной формы. Труды ЦАГИ, 1947, № 635, 47 с.

[21] Шашин В.М. Гидромеханика. Москва, Высшая школа, 1990, 360 с.

[22] Фихтенгольц Г.М. Курс дифференциального и интегрального исчисления. Т. 2. Москва, Наука, 1970, 800 с.

[23] Abramowicz M., Stegun I.A. Handbook of mathematical function. New York, Dover, 1965, 832 p.

[24] Плюснин А.В., Доденко И.А. Исследование точности метода распределенной присоединенной массы при расчете нестационарной поперечной нагрузки на деформируемый корпус ЛА при подводном выбросе. Инженерный журнал: наука и инновации, 2013, №7(19). URL: http://engjournal.ru/ articles/841/841.pdf

[25] Сретенский Л.Н. Теория волновых движений жидкости. Москва, Наука, 1977, 816 с.

[26] Димитриенко Ю.И. Механика сплошной среды. В 4 т. Т. 2. Универсальные законы механики и электродинамики сплошной среды. Т. 2. Москва, Изд-во МГТУ им. Н.Э. Баумана, 2011, 560 с.

Статья поступила в редакцию 02.09.2014

Ссылку на эту статью просим оформлять следующим образом:

Плюснин А.В. Моделирование внутреннего и внешнего нестационарного взаимодействия корпуса летательного аппарата с жидкостью методом граничных элементов. Математическое моделирование и численные методы, 2014, № 2, с. 77-100.

Плюснин Андрей Владимирович родился в 1959 г., окончил Днепропетровский государственный университет в 1982 г., аспирантуру МГУ им. М.В. Ломоносова — в 1987 г. Канд. физ-мат. наук, зам. начальника отдела ОАО «ВПК «НПО машиностроение». Автор более 20 научных публикаций. e-mail: andrey.plusnin@gmail.com

Boundary-element-method modelling of inside and outside non-stationary interaction of aircraft body and liquid

© A.V. Plyusnin

Joint Stock Company Military-Industrial Corporation "NPO Mashinostroenia", Reutov town, 143966, Russia

The article considers inner and outer problems of non-stationary interaction between aircraft body and incompressible ideal fluid and statement of the problems in the form of

boundary integral equations. By inner problems we mean vibration of fuel in tanks and by outer problems we mean determination of additional masses and moments of inertia. We provide formula of efficient solutions for these problems by the boundary element method as applied to bodies of revolution and examples of calculations.

Keywords: potential flow, hydrodynamic forces, additional masses and moments of inertia, boundary element method.

REFERENCES

[1] Sedov L.I. Mekhanika sploshnoy sredy [Continuum Mechanics]. Vol. 2. Saint-Petersburg, Lan' Publ., 2004, 560 p.

[2] Loitsyanskiy L. Mekhanika zhidkosty i gaza [Mechanics of Liquid and Gas]. Moscow, Nauka Publ., 1987, 840 p.

[3] Batchelor G.K. An Introduction to Fluid Dynamics. Cambridge, University Press, 2000, 615 p.

[4] Degtyar' V.G., Pegov V.I. Gidrodinamika podvodnogo starta raket [Hydrodynamics of Submarine-Launched Missiles]. Moscow, Mashinostroenie Publ.; Mashinostroenie-Polet Publ., 2009, 448 p.

[5] Yakimov Yu.L. Sbornik Statey [Miscellany]. Moscow, Moscow State University Publ., 2013, 536 p.

[6] Korotkin A.I. Prisoedinennye massy sudostroitel'nykh konstruktsiy: spravochnik [Added Masses of Ship-Building Structures: a Handbook]. Saint-Petersburg, Mor Vest Publ., 2007, 448 p.

[7] Rabinovich B.I. Vvedenie v dinamiku raket-nositeley kosmicheskikh apparatov [Introduction to the Dynamics of Launch Vehicles and Spacecraft]. Moscow, Mashinostroenie Publ., 1975, 416 p.

[8] Mikishev G.N. Eksperimental'nye metody v dinamike kosmicheskikh apparatov [Experimental Techniques in the Dynamics of Space Vehicles]. Moscow, Mashinostroenie Publ., 1978, 248 p.

[9] Sobolev S.L. Uravneniya matematicheskoy fiziki [Equations of Mathematical Physics]. Moscow, Nauka Publ., 1966, 444 p.

[10] Dimitrienko Yu.I. Mekhanika sploshnoy sredy [The Mechanics of a Continuous Medium].Vol. 1. Tenzorny analiz [Tensor Analysis]. Moscow, BMSTU Publ., 2011, 463 p.

[11] Cruse T.A., Rizzo F.J., eds. Boundary integral equation method: computational applications in applied mechanics. New York, ASME, 1975, 390 p.

[12] Baneijee P.K., Butterfield R. Boundary Element Methods in Engineering Science. London, McGraw-Hill Book Company (UK) Limited, 1981, 494 p.

[13] Brebbia C.A., Telles J.C.F., Wrobel L.C. Boundary Element Techniques. Berlin, Springer-Verlag, 1984, 523 p.

[14] Cebeci T., Shao J.P., Kafyeke F., Laurendeau E. Computational Fluid Dynamics for Engineers. Long Beach, Horizons Publishing Inc., 2005, 396 p.

[15] Vlasova E.A., Zarubin V.S., Kuvyrkin G.N. Priblizhennye Metody Matematicheskoy Phiziki [Approximate Methods of Mathematical Physics]. Moscow, BMSTU Publ., 2001, 700 p.

[16] Zarubin V.S., Kuvyrkin G.N. Matematicheskoe modelirovanie i chislennye metody — Mathematical Modelling and Numerical Methods, 2014, no. 1, pp. 5-18.

[17] Plyusnin A.V. Inzhenerny zhurnal: nauka i innovatsii — Engineering Journal: Science and Innovation, 2013, no. 7(19). Available at: http://engjournal.ru/ catalog/mathmodel/aero/847.html

[18] Demmel J.W. Applied Numerical Linear Algebra. Philadelphia, SIAM, 1997, 430 p.

[19] Kalitkin N.N. Chislennye Metody [Numerical Methods]. Moscow, Nauka Publ., 1978, 512 p.

[20] Rieman I.S., Kreps R.L. Trudy TsAGI — TsAGI Science Journal, 1947, no. 635, pp. 47.

[21] Shashin V.M. Gidromekhanika [Hydromechanics]. Moscow, Vysshaya Shkola Publ., 1990, 360 p.

[22] Fikhtengol'ts G.M. Kurs differentsial'nogo i integral'nogo ischisleniya [The Course of Differential and Integral Calculus]. Vol. 2. Moscow, Nauka Publ., 1970, 800 p.

[23] Abramowicz M., Stegun I.A. Handbook of mathematical function. New York, Dover, 1965, 832 p.

[24] Plyusnin A.B., Dodenko I.A. Inzhenerny zhurnal: nauka i innovatsii — Engineering Journal: Science and Innovation, 2013, no. 7(19). Available at: http://engjournal.ru/articles/841/841 .pdf

[25] Sretensky L.N. Teoriya volnovykh dvizheniy zhidkosti [The Theory of Wave Motions of a Fluid]. Moscow, NaukaPubl., 1977, 816 p.

[26] Dimitrienko Yu.I. Mekhanika sploshnoy sredy. Tom 2. Universalnye zakony mekhaniki i elektrodinamiki sploshnoy sredy. [The Mechanics of a Continuous Medium. Vol. 2. Universal Laws of Mechanics and Electrodynamics of Continuous Media]. Moscow, BMSTU Publ., 2011, 560 p.

Plyusnin A.V. (b.1959) graduated from the Dnepropetrovsk State University, speciality "Hydroaerodynamics" and took a post-graduate course at the Wave and Gas Dynamics Department of Lomonosov Moscow State University. Ph.D., assoc. professor of the Computational Mathematics and Mathematical Physics Department at Bauman Moscow State Technical University, Deputy Chief of Section at JSC "MIC "NPO Mashi-nostroenia". Author of several articles on unsteady problems in hydrogasdynamics and theory of elasticity. e-mail: andrey.plusnin@gmail.com