Sloshing of a liquid fuel in toroidal tanks with account for capillary effect.pdf В современных ракетах и космических аппаратах имеются баки, частично заполненные жидким топливом. Явление плескания топлива в ракетах при больших эффективных силах тяготения и последствия этого явления широко известны и вполне изучены [1-3]. В настоящее время в связи с созданием орбитальных станций и разгонных блоков важную роль приобретает проблема плескания топлива в условиях микрогравитации (g = 10-6-10-4go, где go = 9.81 м/сек2), когда заметно проявляется влияние силы поверхностного натяжения. В российских и зарубежных монографиях [4-5] обобщены результаты исследования статики и динамики жидкости в условиях, близких к невесомости. Равновесная свободная поверхность жидкости в условиях микрогравитации принимает криволинейную форму, которая отличается от плоской в наземных условиях и является функцией ряда параметров, в том числе геометрии полости, угла смачивания, числа Бонда и объема заполнения. До сих пор разработаны только численные методы для определения формы свободной поверхности капиллярной жидкости в осесимметричных сосудах простой формы (цилиндр, сфера и эллипсоид) [6-7]. В статьях [8-11] выражения потенциала скоростей жидкости и поля смещения представлены по ряду характерных функций и получена задача на определение собственных частот и форм колебаний капиллярной жидкости. Предложенные в этих работах приближенно-аналитические методы подходят только для сосудов простой формы. Для решения подобных задач в статьях [12-14] были использованы модифицированный метод Галеркина и метод конечных элементов. Следует отметить, что в последнее время стали использоваться топливные баки более сложной формы - в виде коаксиального цилиндра и тороидальные, однако в них поведение жидкости с учетом капиллярного эффекта исследовано недостаточно. 152 Юй Ч. Волновые движения жидкого топлива в тороидальных сосудах Топливный бак представляет собой гидромеханическую систему газ-жидкость-твердая стенка. Вывод условий равновесия данной системы в условиях микрогравитации был подробно описан в статье [15], а в работе [16] представлено решение задачи о равновесии и колебаниях капиллярной жидкости в коаксиальноцилиндрических сосудах. В данной работе исследованы собственные колебания капиллярной жидкости в тороидальных сосудах на основе метода конечных элементов. В дальнейшем будем говорить о капиллярной или тяжелой жидкости в зависимости от того, учитывается капиллярный эффект или нет. Постановка задачи о малых колебаниях капиллярной жидкости Пусть вектор ускорения g действует параллельно продольной оси симметрии сосуда. Используем длину дуги s в качестве переменной для описания формы равновесной свободной поверхности. Введем цилиндрическую систему координат o^z (рис. 1), а также на свободной поверхности Го криволинейную систему координат osѲh таким образом, чтобы поверхность Г0 имела уравнение h = 0, а координатные линии h были направлены по внешней нормали поверхности Г0. Рис. 1. Основные обозначения параметров жидкости: zo(s) и ro(s) - функции равновесной свободной поверхности; ао - угол смачивания жидкости; h (s, Ѳ, t) - отклонение возмущенной свободной поверхности Г от Го по внешней нормали; ү - линия трехфазного контакта; Z - смачиваемая поверхность сосуда; А - область, которую занимает жидкость; ri - радиус осевой окружности тора; п - радиус окружности меридиана; n и ni - нормали свободной поверхности и смачиваемой поверхности соответственно; e - внешняя нормаль линии контакта в касательной плоскости свободной поверхности Fig. 1. General designations for parameters of a liquid: zo(s) and r0(s) are the functions of an equilibrium free surface; a0 is the contact angle; h (s, Ѳ, t) is the deviation of the perturbed free surface Г from the equilibrium one Г0 along the outward normal; ү is the contact line; Z is the wetted surface of a vessel; G is the the area occupied by the liquid; ri is the radius of a torus axial circle; n is the radius of a torus meridian circle; n and ni are the normals to the free surface and to the wetted surface, respectively; and e is the outward normal to the contact line in a tangent plane of the free surface Будем считать, что равновесная поверхность жидкости Г0 в неподвижном сосуде определена (z = z0(s), r = r0(s)), и жидкость совершает малые колебания, которые подчиняются: 153 - уравнению Лапласа д2ф 1 Эф 1 д2ф д 2ф Дф = ^ + -^ + ^^- + ^г = 0 в Q, (1) ф дг2 r dr r2 дѲ2 dz2 где ф - потенциал скоростей жидкости; - условию непротекания на смачиваемой поверхности Эф = 0 на S. (2) дп Возмущенная свободная поверхность Г может быть представлена в виде: h = h (s, Ѳ, t), и линеаризованное кинематическое условие запишется как: (3) Механика / Mechanics дп д/ Динамическое условие на свободной поверхности вытекает из линеаризованного интеграла Коши-Лагранжа: Эф + Ed. + ^ = 0 на Г„, д/ р где pd - динамическая составляющая давления и zd - проекция отклонения свободной поверхности вдоль оси z. В предположении отклонения свободной поверхности вдоль нормали zd выглядит следующим образом: zd = hn • ez = hr0s , где ez - единичный вектор оси z, а ros - первая производная функции ro(s) по длине дуги s. Выражение динамической составляющей давления pd частиц жидкости на свободной поверхности определяется из условия Лапласа для скачка давления p0 - p = 2сК, где с - коэффициент поверхностного натяжения и K - средняя кривизна свободной поверхности. Для равновесной и возмущенной свободной поверхности условие Лапласа имеет выражение: Eo -Es = -(ki + k2), Eo -E = -{(ki + k2) + [(k2 + к2)h + ДГh]}, где к1 и к2 - главные кривизны равновесной свободной поверхности, определяемые уравнениями: кі = roszoss - rosszos, к2 = zos/ro, здесь нижние индексы s и ss означают первую и вторую производные функции по длине дуги s соответственно; д2 r д 1 д2 Дг = -- + - + - -- - оператор Лапласа-Бельтрами на поверхности Гo. дs2 r дs r2 дѲ2 Из соотношения pd = p - ps получим выражение динамической составляющей давления на свободной поверхности: Pd = --[( к2 + к2) h + Дг h]. Тогда получим динамическое условие на свободной поверхности в виде: ^7 + 8hros - - [(кі2 + к2)h + Дгh] = 0 на Го • (4) д/ р[ ] Для вывода краевого условия на линии трехфазного контакта используем условие Дюпре-Юнга: nni = cos(ao), и предположение сохранения угла смачива-154 Юй Ч. Волновые движения жидкого топлива в тороидальных сосудах ния 5(ппі) = 0 в процессе плескания жидкости. Тогда получим граничное условие на линии контакта [4. С. 103]: - + %h = 0, где х д е k1 cosа0 - kf sin а0 на у, (5) где &і - кривизна твердой стенки сосуда на линии контакта. Введем радиус окружности меридиана тора Г2 как характерный размер длины и характерные значения для времени и потенциала: t* = (pr23/c)12, ф* = r22/ t*. Подставив величины x = r2x\\ t = t*t' и ф = ф*ф' в уравнения (1-5), получим безразмерную формулировку задачи, в которой опускаем верхний знак ' над буквами: * „ _ Эф Л _ Эф дИ „ Дф = 0 в Q, - = 0 на f, - = - на Г„, дп дп dt °ф + B Иг dt Bo 0 -\\_(kl + k2)И + ДгИ] = 0 на Го, 0К + %И = 0 на у де (6) где Bo = pgr22/o - число Бонда, характеризующее соотношение массовой силы и силы поверхностного натяжения. По сравнению с постановкой задачи о малых колебаниях жидкости без учета силы поверхностного натяжения капиллярная специфика в поставленной задаче (6) проявляется в динамическом условии на свободной поверхности и в граничном условии на линии трехфазного контакта. Вариационная формулировка краевой задачи Для определения конкретного решения поставленной задачи (6) следует задать начальные условия, например поле смещений h и скорость dh/dt свободной поверхности в некоторый момент времени. В данной работе рассматриваются собственные колебания, т.е. решения, зависящие от времени t по гармоническому закону. Тогда потенциал скоростей жидкости ф(г, z, Ѳ, t) и поле смещений свободной поверхности h(s, Ѳ, t) можно представить в виде: да да ф(г,z,Ѳ,t) = ^фп (г,z,Ѳ)е‘°ѵ, h(s,Ѳ,t) = ^hn (s,Ѳ)emn. n=l n=l где юя - безразмерная частота n-го тона колебаний жидкости. Исключив переменную времени, получим задачу о собственных колебаниях жидкости с учетом силы поверхностного натяжения: дф дф Дф = 0 в Q, -- = 0 на f, -- = /юиИ на Г, дп дп /Юпфп + B0hnr0s - \\(k2 + k2)К, + дгК ] = 0 на Г0, + хИп = 0 на у. Запишем вариационную формулировку задачи в виде: -j Дфп§фп^^ + J Фndf + j ^°фП + І®пК j ЪфndГ0 + j {[B0r0s Г0 (k 2 + k 2)] hn -ДГ hn + *Чфп }8МГ0 + + Xh„ у = 0. 155 Механика / Mechanics Используя интегрирование по частям dE - J ^5фп^Г0 + J Ѵфп8Ѵфп^а dn -[Дфи5Фи^= Гдф, E -J Дг hnbhnd Го = -jf^M y+J Ѵг НпЪѴг hnd Го, Го y Го получим формулировку задачи в вариационном виде: 5Пі = 0, где Пі - функционал, определяемый выражением П1 = J^n^ndQ + J {[B0r0s -(k2 + k 2)] hl +ѴГҺп ѴГhn } dГ0 + Jxhn'dy. Q Го y Пі = При использовании кинематического условия на свободной поверхности функционал Пі приобретает иной вид: d Г + + d у-^п JѴФnѴФnd Q Q J* у где In = an2. Чтобы установить связь между потенциалом скоростей жидкости фп и скоростью движения частиц жидкости на свободной поверхности дфп/дп, необходимо ввести вспомогательную задачу Неймана: Дф = 0 в Q, дфпі = 0 на E, дфпі = f на Г„, on дп где f - заданная функция на свободной поверхности жидкости. Вариационная формулировка задачи Неймана имеет вид: 8фnd Г0 = °. 0, где П2 - Аналогично можно преобразовать задачу Неймана в виде: 5П2 функционал, имеющий вид: 1 J^n^nd Q. ^ г. п = J дф^ Фndг Сосуд осесимметричный, и волны в окружном направлении обладают свойством периодичности. Тогда можно записать фп (r, z, Ѳ) и hn (s, Ѳ) в виде: ф(r,z,Ѳ) = Ф(г,z)cos(тѲ), hn (s,Ѳ) = H(s)cos(тѲ), m = 0,1,2... Отметим, что колебания, отвечающие числу волн в окружном направлении m =1, имеют антисимметричную моду колебаний, и в этом случае возникает переменное результирующее давление на сосуд в поперечном направлении. Данный случай вызывает большой интерес в инженерной практике при проектировании систем управления космических аппаратов. Важное внимание также уделяется колебаниям, отвечающим значению m = 0, в которых переменна составляющая давления вдоль оси сосуда. При заданном значении m можно преобразовать осесимметричную область жидкости из трехмерной в двухмерную: 156 Юй Ч. Волновые движения жидкого топлива в тороидальных сосудах (24) 8Ц = 0 и 8П2 = 0, rdrdz. Определение собственных частот и форм колебаний методом конечных элементов Вариационная формулировка задачи (7) пригодна для применения метода конечных элементов, с помощью которого задача с распределенными параметрами преобразуется в задачу с конечными степенями свободы. В данной работе выбираем треугольный элемент для дискретизации области, занимаемой жидкостью, в правом сечении тороидального сосуда на рис. 1. Допустим, что на контуре свободной поверхности Го имеется p элементов линии, и ;-й элемент имеет и узлов, а в плоском сечении области Q, занятой жидкостью, получены q элементов треугольника, и ;-й элемент имеет v узлов. Функции (Ф); и (5Ф/5и),- в ;-м элементе на свободной поверхности Г0 можно представить по узловым значениям данного элемента [17-18]: где Mj - функция формы j-го узла элемента линии, Ф j и (дФ/ди)у - j-е узловое значение функций Ф и 5Ф/5и в ;-м элементе линии. Таким образом, функция (Ф); в ;-м элементе в области жидкости Q имеет вид: v (ф) = Y^Nj > j=1 где Nj - функция формы j-го узла элемента треугольника, Фу - j-е узловое значение функции Ф в ;-м элементе треугольника. Поставив эти функции в функционал П2, rdrdzфft, получим функционал П2 в матричном виде: 157 Механика / Mechanics Здесь Г0i и Qi - i-й элемент на свободной поверхности и в области жидкости соответственно. В дальнейшем обозначим п/ - общее количество узлов на свободной поверхности, а т - общее количество узлов в области жидкости, Ф'=(Ф.1’ Ф/2. -> ФГ = (Ф,1, Ф, (f- \\dn)t [\\дп]п v дп )i2 {дп ,, •••, Ф,ѵ), В,. =[в'*] ив'*= j MjMkrds, 0i Г А, = [ Ajk ] и Ajk = J ц ф:=(ф,.ф„-,ф,„), (£) dNj dNk dNj dNk m NjNk dr dr dz dz r2 T f ӘФ dn ): d Ф dn rdrdz, ӘФ dn 'f J Фг =(®1, ф2> ф,, • фч) = (фГ. ф[) При составлении матриц элемента Bi и Ai применяется метод Гаусса, позволяющий повысить точность вычисления путем специального выбора узлов интегрирования [18]. Матрицы B и A составляются из матриц Bi и Ai соответственно, причем матрица B в порядке п/ х п/и матрица A может быть записана в виде: А = Аз А4 где Ai в порядке п/ х п/, A2 в порядке п/ х (щ-п/), A3 в порядке (т-п/) х п/ и A4 в порядке (п-п/) х (щ-п/). Выполнив вариацию функционала 5П2 = 0 5П2 =5ф[ (В|Ф - А^ - А2Ф2 j - 5ф[ (АзФ2 + А4Ф2) = 0, получим следующую систему уравнений: dФ А,Ф! + А2Ф2 = В-, ЛзФ2 + А4Ф2 = 0. dn После решения системы уравнений представим Фі, Ф2 через 5Ф/5п: ,dO dФ ф = А-1 -, Ф --А,-, dn dn где А5 = В-1 (А1 - Л2 А-1А3), А6 = А-1А3 А-1. Аналогично преобразуем функционал Пі в матричной форме: П' -(£I(яоВі-В2+Вз+в4)f(ф 1, ф)(А; А41ф1 где матрицы Bi, B2, B3 и B4 имеют порядок п/ х п/ и составлены из матриц Bii, B21, B3 i и B4 і соответственно, bjk - символ Кронекера, Ві, = [Bj] и Bj = J roMjMkrds, Г 0i В, = [в j] и B j - J (k 1 + k2)MjMkrds, Г 0i 158 Юй Ч. Волновые движения жидкого топлива в тороидальных сосудах Вэ,- = [в %] и в % = j ЭМ% дЫк тгЫ%Ык rds, ds ds B4i = [В % ] и В % = (xr )у § %. Исключив Фі, Ф2 через ЭФ/Эп, запишем функционал Пі в виде: ( эф Y , ч эф П'=Ь,Г 1 (К-».„ м)-, где к = в0в - в2 + В + В, A6A3 Л-1 + А6А4 Л6. м = ( Л-1)Г AjAs1 -(А-f л2 Л6 При составлении матрицы M применяется метод Guyan reduction [19] для уменьшения числа степеней свободы задачи. Данный подход позволяет получить задачу, степени свободы которой совпадают с количеством узлов на свободной поверхности. Из вариационной формулировки задачи 5Пі = 0 получим задачу на определение собственных частот и форм колебаний капиллярной жидкости: (К ->-м £="■ Из вышесказанного вытекает, что методика решения задачи о малых колебаниях капиллярной жидкости в любом осесимметричном сосуде состоит в следующем: 1) определение формы равновесной свободной поверхности капиллярной жидкости z = z0(s), r = r0(s) и получение функций r0s, k1, k2, x; 2) создание геометрической модели области жидкости и выполнение ее дискретизации треугольными элементами в программе PARTRAN; 3) обработка информации о координатах узлов и элементов в программе MATLAB; 4) получение матриц A, B, Bi, B2, B3 и B4 с использования численного интегрирования функции методом Гаусса; 5) решение задачи на определение собственных частот и форм на основе метода Guyan reduction. Обсуждение результатов Отметим, что в данной работе исследованы только колебания, отвечающие числу волн в окружном направлении m = 1. В этом случае мода колебаний жидкости антисимметрична (cosѲ), и возникает переменное результирующее давление на сосуд в поперечном направлении. Так как в литературе отсутствуют данные о колебаниях жидкости в тороидальных сосудах в условиях микрогравитации, проверка достоверности и сходимости разработанного алгоритма проведена в примере сферического сосуда. При в = 50% и Оо = 90° свободная поверхность жидкости в сферическом сосуде принимает плоскую форму, и полученная частота первого тона Х1* хорошо совпадает с численным значением, приведенным в литературе [12]. По сравнению с численными значениями из работы [13] видно, что только при 159 Механика / Mechanics в = 78% есть заметное рассогласование, а при остальных объемах заполнения сосуда жидкостью имеется хорошее совпадение (табл. 1). Можно сделать вывод, что разработанный алгоритм может быть пригоден для исследований колебаний капиллярной жидкости в топливных баках. Таблица 1 Проверка достоверности алгоритма в случае сферического сосуда, Xi* = Xi/Бо и р - относительный объем заполнения сосуда жидкостью Во = 1, ао = 5° В0 = 2, а0 = 5° В0 = 1 000, а0 = 90° в, % 25 50 78 25 50 78 50 Ал*, полученные в данной статье 1.13 1.27 1.55 1.14 1.28 1.53 1.56 Аі в литературе 1.33 1.41 2.03 1.11 1.23 1.83 1.54 _[13_ _[12_ Данный алгоритм на основе метода конечных элементов также сходится при увеличении степени дискретизации области жидкости; соответствующие результаты приведены в табл. 2. Т аблица 2 Анализ сходимости алгоритма при увеличении степени дискретизации области жидкости В0 = 1, № = 5° В0 = 2, а0 = 5° В0 = 1 000, а0 = 90° в, % 25 50 78 25 50 78 50 Шаг 1 1.148 1.304 1.811 1.165 1.320 1.784 1.561 Шаг 2 1.140 1.290 1.691 1.155 1.310 1.623 1.560 Шаг 3 1.132 1.277 1.584 1.149 1.288 1.560 1.559 Шаг 4 1.129 1.267 1.550 1.143 1.280 1.529 1.558 На рис. 2 показаны формы равновесной свободной поверхности капиллярной жидкости в зависимости от числа Бонда. Так как сосуд осесимметричный и вектор ускорения параллелен оси симметрии, равновесная свободная поверхность представляет собой осесимметричную криволинейную поверхность, и здесь показан только ее контур в правом сечении сосуда. При уменьшении числа Бонда свободная поверхность более искривлена, и ее ориентация изменяется против хода часовой стрелки в тороидальных сосудах. Возможной причиной подобного явления является то, что при уменьшении числа Бонда влияние силы поверхностного натяжения становится больше и свободная поверхность стремится принимать форму с большей кривизной. Подобное явление наблюдалось в экспериментах, проводимых в башне невесомости, при объемах заполнения сосуда жидкостью в < 20% [20]. При этом с увеличением размера сосуда появляется тенденция увеличения диапазона в, в котором свободная поверхность принимает форму, представленную на рис. 2 при В0 = 0. На самом деле в башне невесомости жидкость не успевает успокоиться до нужной степени за несколько секунд, особенно при большом ее объеме. Можно предположить, что при исключении ограничения времени эксперимента в башне невесомости форма установившейся свободной поверхности будет совпадать с численными результатами. Реальное ракетное топливо имеет угол смачивания на стенке топливных баков, близкий к нулю, и это вызывает большое искривление свободной поверхности 160 Юй Ч. Волновые движения жидкого топлива в тороидальных сосудах около линии трехфазного контакта. На рис. 3 показана схема дискретизации области жидкости, и видно, что около линии трехфазного контакта требуется большая степень дискретизации для точного моделирования поведения жидкости в условиях микрогравитации. Рис. 2. Формы равновесия свободной поверхности в разных числах Бонда Во при ri = 1.87, ао = 5°, в = 50% Fig. 2. Free surface equilibrium shapes at various Bond numbers B0 and n = 1.87, ас = 5°, в = 50% Рис. 3. Схема дискретизации области правого сечения жидкости в тороидальном сосуде Fig. 3. Discretization scheme for a right cross-section of the liquid in a toroidal vessel В табл. 3 приведены собственные частоты колебаний жидкости в условиях микрогравитации, и с повышением числа Бонда В0 значения частот увеличиваются. Для каждого объема заполнения в существуют свои характерные числа Бонда, при превышении которых собственные частоты приближаются к экспериментальным значениям, полученным в наземных условиях [21]. Таблица 3 Собственная частота первого тона Ii* колебаний жидкости при ri = 1.87, ао = 5° Во в \\ 1 2 3 4 5 6 8 10 1000 [21] 20% 0.051 0.074 0.085 0.089 - - - - 0.103 0.084 50% 0.010 0.092 0.134 0.153 0.167 0.176 - - 0.220 0.226 80% - - - 0.046 0.085 0.127 0.191 0.217 0.384 0.406 Собственные формы колебаний жидкости относительно равновесной свободной поверхности в условиях микрогравитации показаны на рис. 4 для правого сечения тороидального сосуда. Отметим, что здесь исследованы только первые четыре тона (п = 1, 2, 3, 4) колебаний капиллярной жидкости при антисимметричной моде m = 1, т.е. cosѲ. Сплошная линия показывает равновесное положение свободной поверхности, а пунктир - собственную форму. 161 Механика / Mechanics Рис. 4. Собственные формы первых четырех тонов колебаний жидкости при п = 1.87, Во = 5, ао = 5°, в = 50% Fig. 4. Natural modes of the first four liquid vibration modes at n = 1.87, B0 = 5, a0 = 5°, в = 50% С точки зрения количества узлов пересечения равновесной поверхности и собственной формы каждого тона колебаний жидкости полученные результаты похожи на экспериментальные результаты для тяжелой жидкости [21]. Заключение В настоящей работе численно определены собственные частоты и формы колебаний капиллярной жидкости в тороидальных сосудах. Чтобы учитывать силу поверхностного натяжения и граничное условие жидкости на линии трехфазного контакта, введена вариационная формулировка задачи и для ее решения разработан численный алгоритм на основе метода конечных элементов. Из результатов исследований следует, что в условиях микрогравитации сила поверхностного натяжения вызывает сильное искривление свободной поверхности и скачок давления на ней. Искривление равновесной свободной поверхности и поведение жидкости около линии трехфазного контакта имеют большое влияние на собственные частоты и формы колебаний капиллярной жидкости, особенно в тороидальных сосудах. При увеличении числа Бонда свободная поверхность становится перпендикулярной к вектору ускорения и приближается к плоской, в этих условиях можно пренебречь влиянием капиллярного эффекта, и полученные результаты сходятся с экспериментальными значениями для тяжелой жидкости.

Скачать электронную версию публикации