О присоединённых массах панелей при колебаниях в несжимаемой среде

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

Неотъемлемой частью многих КА являются крупногабаритные панели солнечных батарей (СБ). Для синтеза динамической модели КА разрабатываются отдельные конечноэлементные модели (КЭМ) панелей СБ. Несмотря на то, что их получают на современных мощных вычислительных программных комплексах типа NASTRAN , возникает необходимость уточнения этих КЭМ, обусловленная применением в конструкции СБ композитных материалов, наличием узлов соединений и раскрытия, сотовых элементов. Для верификации КЭМ СБ проводят модальные испытания [1]. Кроме того, по результатам модальных испытаний определяются диссипативные характеристики упругих колебаний панелей СБ.

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

Проведение модальных испытаний в вакуумных камерах весьма трудоёмко и ограничено размерами этих камер. Для современных конструкций панелей СБ влияние воздушной среды может привести к уменьшению собственных частот на 5…15% и к довольно сильному увеличению коэффициентов демпфирования упругих колебаний, а влияние на собственные формы упругих колебаний практически отсутствует. С увеличением габаритных размеров панелей СБ и уменьшением их погонной массы влияние воздушной среды станет более существенным.

Инженерные оценки влияния воздушной среды на динамические характеристики панелей СБ при модальных испытаниях приведены в работе [2]. Некоторые результаты теоретического и экспериментального исследований влияния воздушной среды, силы тяжести и упругости подвески методической панели содержатся в работе [3]. Решения модельных задач о колебаниях прямоугольных пластинок в жидкости (а воздух при модальных испытаниях можно считать несжимаемой жидкостью) представлены в работе [4]. Основное внимание в этой статье уделялось определению демпфирования, обусловленного периодическим образованием и сходом вихрей с кромок пластины, но приведены также результаты определения присоединённых масс. Применялся метод граничных элементов (МГЭ) с кусочно-постоянной аппроксимацией решения на элементах.

Описание МГЭ с более точной кусочнолинейной аппроксимацией решения на элементах содержится в работе [5]. Результаты численной реализации этого метода на примере определения поля скоростей жидкости при обтекании шара и прямоугольной пластины приведены в работе [6], но присоединённые массы не представлены. Вариант МГЭ с квадратичной аппроксимацией решения, позволяющий получить также и надёжные оценки величины нелинейного вихревого демпфирования, представлен в работе [7], но только для плоских задач. Этот метод может быть применён для конструкций панелей СБ с рулонным типом раскрытия и сильным удлинением в одном направлении.

В настоящей статье рассматривается применение МГЭ к задаче о колебаниях плоской прямоугольной панели в воздушной среде. Для различных форм колебаний панели как недеформируемой конструкции исследуется зависимость присоединённой массы воздуха от отношения её геометрических размеров и относительной величины просвета между её фрагментами.

Постановка задачи

Кинематическая вязкость воздуха в обычных лабораторных условиях испытаний ν ≈ 15 мм²/с [8]. Принимая низшую круговую частоту колебаний панели ω ≈ 1 рад/с, следуя источнику [9], находим характерную толщину пограничного слоя на её поверхности (2 ν / ω )^1/2 ≈ 5,5 мм. По сравнению с характерным размером панели R ₀ ~2 м эта толщина мала.

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

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

Решения некоторых задач об определении демпфирования колебаний тонких пластин в несжимаемой среде приведены в работе [4]. Установлено, что нелинейная часть демпфирования, обусловленная вихревым движением в окрестности кромок пластин, зависит от амплитуды колебаний в степени 2/3. Для переноса полученных для тонких пластин результатов на панели СБ, кромки которых отличаются конструктивными особенностями, требуется экспериментальное обоснование.

При малых амплитудах колебаний существенное вихревое движение локализуется в тонком пограничном слое и малых областях в окрестности кромок панели. В основном объёме среды, окружающей панель, её движение потенциально. Экспериментальными исследованиями подтверждено [4, 10], что при малых амплитудах колебаний инерционные силы с достаточной точностью определяются потенциальным движением среды.

Скорость в точках граничной поверхности панели представим в виде v ( t ) f ( R ), где f ( R ) — векторная форма колебаний панели; v ( t ) — функция времени t . Для гармонических колебаний v ( t ) = v ₀cos( ω t ). Определение потенциала скоростей Φ ( t , R ) = v ( t ) ϕ ( R ) приводит к решению уравнения Лапласа с граничным условием на поверхности S панели, т. е. к решению краевой задачи Неймана [9]

∇ ² ϕ = 0; ∂ϕ / ∂ν| _S = νf , (1) где ν — орт внешней, относительно среды, нормали к поверхности панели. Если принять условие нормировки формы колебаний νf = 1 в некоторой точке панели, то v ( t ) — проекция скорости на направление нормали в этой точке.

Кинетическая энергия воздушной среды

T = — p j ( УФ )² dV = — ц v ²;

2 v             2

ц = p j ( Уф )² dV = p f ф _d V- dS ,

где ρ — плотность воздуха, а коэффициент μ, стоящий перед квадратом скорости, называется присоединённой массой. Динамическое давление на поверхности панели определяется выражением [9]

∂Φ

p ( t , R ) = – ρ       = – ρ v ^.( t ) ϕ ( R ).

∂ t

Для учёта влияния воздушной среды на форму колебаний панели СБ необходимо решать совместную задачу упругости и гидродинамики. Решения этих задач можно разделить. Провести решение задачи упругости и определить форму колебаний панели. Решить задачу гидродинамики, используя полученную форму колебаний. Провести решение задачи упругости с учётом гидродинамического давления (3) на панель, что эквивалентно добавлению присоединённой массы к элементам панели. Если форма колебаний сильно изменится, то процесс повторить.

∂ w

w = ei

∂ e

1    R ^. e

.

.

4 π    R ³

Функция w_ei определена в некоторой точке i внутри объёма среды в связанной с ней системе координат. Она удовлетворяет уравнению ∇ ² w_ei = – ∂δ _i / ∂ e , где δ -функция задана в точке i .

Применяя формулу интегрирования по частям [12] к левой части уравнения метода взвешенных невязок (4), получим

f

S

—— w dS - f V— w dV = w dS - f fv w dS , ei                      ei                         ei                    ei

∂ν         V           S ∂ν S

где вместо весовой функции w подставлена функция w_ei . Два поверхностных интеграла взаимно сокращаются. Последний поверхностный интеграл в правой части равен нулю, так как в соответственных точках на разных сторонах поверхности вектор нормали ν имеет противоположные направления.

Повторное интегрирование по частям даёт

Метод решения задачи гидродинамики

Применяя к задаче (1) метод взвешенных невязок [11], запишем

V

J w V 2 — dV = f —— - vf wdS, S ∂ν

где w — некоторая весовая функция. Как правило, для решения МГЭ задач о потенциале применяют фундаментальное решение

11 w =    ^.

4π R уравнения Лапласа ∇2w = –δ, где δ-функ-ция определена внутри объёма V вне граничной поверхности; R — радиус-вектор от этой особой точки.

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

В качестве весовой функции используем производную фундаментального решения уравнения Лапласа по направлению единичного вектора e

- f ф —e- ds + J фУ 2 w i dv = o.

S ∂ν V ^e

Используя свойства δ -функции [12], вычисляем интеграл по объёму

1 фу² w_edV =- j ф ^eiⁱ dV = j -| ^ - sdV= , (7)

V              V ∂e V ∂e         ∂e где ∂ϕi /∂e — производная потенциала в точке i по направлению вектора e. Подставляя интеграл (7) в выражение (6) и переходя к интегрированию по одной стороне граничной поверхности, приходим к интегральному уравнению

∂ϕ _i

----= I u

∂ e    S

∂ w s dS , ∂ν

где u = ϕ ₊ – ϕ _– — разность значений потенциала на разных сторонах граничной поверхности. Это уравнение позволяет определить производную потенциала по любому направлению, если известна функция u на граничной поверхности.

Выберем точку i на малом расстоянии ε от граничной поверхности и определим весовую функцию w_ni , заменяя вектор e в уравнениях (5) и (8) на орт внешней нормали n к границе. Переходя в уравнении (8) к пределу ε→ 0, получим граничное интегральное уравнение

∂ϕ     ∂w i = u    n i dS; w = ni

∂ n   S ∂ν

^{1 Rn} 4 π ^. R ^{3 ,}

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

Поверхностный интеграл в уравнении (9) имеет сильную особенность типа Адамара [6]. Преобразуем этот интеграл к интегралу с более слабой особенностью типа Коши [12], считая ε ≠ 0. Используя тождество [13]

nR      n x R

V --- = - Vx----:

R ³           R ³

которое справедливо, так как для выбранной точки i вектор n — постоянный, находим

∫ ^u

S

∂ w ni

∂ν

1     ∂ Rn dS = –    ∫ u        dS =

4 π S ∂ν R ³

4 π S

n × R

равны их значениям в центре ГЭ. В таком случае на каждом ГЭ ∇ u = 0 и в выражении (11) исчезает поверхностный интеграл. Записывая уравнение (9) для производной потенциала по нормали в точках центра каждого ГЭ и учитывая выражение (11), приходим к системе линейных алгебраических уравнений

∂ϕ 1 N      n × R

—= — S u f — d R ; i = 1, 2, ..., N .(12)

∂ n    4 π j = 1 ^j R ³

l j

Интеграл в уравнении (12) вычисляется по контуру ГЭ j ; векторы n и R заданы в локальной системе координат, связанной с центром ГЭ i ; N — число ГЭ. Геометрическая форма граничной поверхности может быть произвольной. На поверхности каждого ГЭ может быть задана своя локальная криволинейная система координат.

n × R

= u ∇×      ν dS .

4 π S

_R 3

Применяя теорему Стокса

_rr ( n xR)"l          n x R

Vx u ---- v dS = ^u----d R

Л I r 3 J ^f   r ³

и учитывая тождество [13]

( n x R)      n x R

Vx I u---- = V ux----+ u Vx

R ³            R ³

n x R

R ^{3 ,}

приведём интеграл (10) к требуемому виду

∫ ^u

S

∂ w

-dS dS =--f u dv       4n l

nxR _ dR R3

- [ V u x

4 n J _

n x R

_R 3

ν dS .

В методе граничных элементов поверхность представляется треугольными или четырехугольными граничными элементами (ГЭ). Если применяется вариант МГЭ, в котором выполняется межэлементная непрерывность функции u , то контурный интеграл в формуле (11) исчезает, так как на краях панели разрыв потенциала равен нулю, а внутренние стороны ГЭ обходятся дважды в противоположных направлениях.

Наиболее проста численная реализация варианта МГЭ, в котором функция разрыва потенциала u и нормальная производная потенциала на каждом ГЭ

Результаты расчётов и анализа

Приведём результаты определения присоединённых масс для различных форм колебаний плоской прямоугольной панели, показанной на рис. 1. Оси x и y декартовой системы координат Oxyz находятся в плоскости панели и являются осями симметрии панели. Панель может состоять из двух фрагментов, разделённых зазором 2 Δ . Длина каждого фрагмента вдоль оси x равна a , вдоль оси y — 2 b .

Рис. 1. Схема плоской прямоугольной панели с просветом 2 Δ между одинаковыми фрагментами панели: a — длина фрагмента; b — полуширина фрагмента

Панель представлялась равномерной сеткой прямоугольных ГЭ. Применялся вариант МГЭ с кусочно-постоянной аппроксимацией решения на ГЭ, основанный на уравнениях (12). Интегралы по контурам прямоугольных ГЭ в этих уравнениях определены аналитически. Расчёты выполнены для форм колебаний панели, симметричных

и антисимметричных относительно осей координат. Граничные элементы и неизвестные переменные u_j в их центре задавались на четверти панели. На остальной части панели ГЭ и переменные u_j определялись с учётом свойств симметрии, что в четыре раза снизило размерность системы линейных уравнений (12). Расчёты проведены по программе, составленной на языке FORTRAN .

Численное решение системы уравнений (12) с увеличением количества ГЭ сходится медленно, что объясняется сингулярной особенностью потенциала скоростей вблизи краёв панели. Для повышения точности вычислений присоединённой массы μ , определяемой по формуле (2), применялось нелинейное преобразование Шенкса [14]

Ц =

и „ц - ц/. n +2 n n +1

ц + ц - 2и.

n +2       n n +1

где μ _n , μ _{n +1}, μ _{n +2} — значения присоединённых масс, полученные на трёх последовательных геометрически подобных сетках ГЭ. Начальная сетка ГЭ задавалась, а каждая следующая сетка строилась путём деления прямоугольных ГЭ предыдущей сетки на четыре одинаковых прямоугольника. Это преобразование может использоваться также для уточнения значений разрыва потенциала u_j в точках j на панели. На четверти панели использовалось до 2 048 граничных элементов.

В табл. 1 представлены результаты определения присоединённых масс для поступательных колебаний панели, не имеющей зазора между фрагментами, в зависимости от её удлинения. Форма колебаний панели f ( x , y ) = 1. Приведены безразмерные коэффициенты присоединённых масс C_T = μ /(2 πρ b ² a ). Присоединённые массы вычисляются по формулам

μ = C_T ( a/b )2 πρ b ² a , a ≥ b ;

(13) μ = C_T ( b/a ) 2 πρ a ² b , a ≤ b .

Если a ≤ b , то коэффициенты C_T определяются из табл. 1 заменой a ↔ b . Как видно, коэффициент присоединённой массы существенно зависит от длины панели. Если a → ∞, то C_T → 1, так как присоединённая масса на единице длины бесконечно длинной панели равна πρ b ² [15].

Таблица 1

Зависимость присоединённой массы от длины панели. Форма колебаний панели f(x, y) = 1

a/b	1,0	1,5	2,0	3,0	4,0	5,0
C T	0,579	0,689	0,756	0,832	0,872	0,896

В табл. 2 представлены результаты определения присоединённых масс для колебаний вокруг оси x панели, не имеющей зазора между фрагментами, в зависимости от удлинения панели. Форма колебаний панели принималась в виде f ( x , y ) = y/b , поэтому вычислялись присоединённые массы при нормировке на максимум формы колебаний. Приведены безразмерные коэффициенты присоединённых масс C_R = μ /( πρ ab ²/4). Присоединённая масса при колебаниях панели вокруг оси x вычисляется по формуле

μ = C_R ( a/b ) πρ ab ²/4. (14)

Как известно [15], присоединённый момент инерции на единице длины бесконечно длинной панели равен πρ b ⁴/8. Для определения присоединённого момента инерции следовало задавать форму колебаний в виде f ( x , y ) = y . Так как принята форма колебаний, которая в b раз меньше, то соответствующая ей присоединённая масса равна πρ b ²/8. Если a → ∞, то C_R → 1.

Таблица 2

Зависимость присоединённой массы от длины панели. Форма колебаний панели f(x, y) = y/b

a/b	0,125	0,25	0,50	1,00	2,00	4,00
C R	0,274	0,455	0,659	0,816	0,906	0,952

Для колебаний вокруг оси x панели удлинённой формы оценку присоединённой массы можно получить, применяя простую зависимость b

μ ≈ πρ a ² ∫ ( y/b )² dy .            (15)

–b

Применение этой зависимости для панели с отношением сторон a/b = 0,125, при котором расстояние от оси вращения до краёв панели в четыре раза больше ширины панели, даёт величину присоединённой массы 16πρa3/3. По формуле (14), в которой из табл. 2 CR = 0,274, получаем величину присоединённой массы 4,384πρa3. Погрешность приближённой оценки составляет 21,7% в сторону завышения. Для отношения сторон a/b = 0,25 погрешность приближённой оценки увеличивается до 46,5%.

Присоединённая масса для колебаний панели вокруг оси y , если форма колебаний задана в виде f ( x, y ) = x / a , определяется по формуле

μ = C_R ( b / a ) πρ ba ²/4,         (16)

получаемой путём замены a ↔ b в уравнении (14). Используя зависимости (13), (14), (16), применением принципа линейной суперпозиции можно определить присоединённую массу панели для любой плоской формы колебаний.

Пусть ϕ ₁( R ), ϕ ₂( R ), ϕ ₃( R ) — потенциалы скоростей для форм колебаний панели f ₁( x , y ) = 1, f ₂( x, y ) = y / b , f ₃( x , y ) = x / a , соответственно. Тогда для потенциала

ϕ ( R ) = α ₁ ϕ ₁( R ) + α ₂ ϕ ₂( R ) + α ₃ ϕ ₃( R )

на поверхности панели выполняется граничное условие

f ( x, y ) = α ₁ + α ₂ y/b + α ₃ x/a .     (17)

Из свойств симметрии следуют условия ортогональности

∫ ∇ϕ _i ∇ϕ _jdV = 0, i ≠ j ( i , j = 1, 2, 3). V

Учитывая эти условия ортогональности и граничное условие для потенциала ϕ ( R ), получим следующее выражение присоединённой массы панели

μ = ρ ∫(∇ϕ)2dV = α12μ1 + α22μ2 + α32μ3, (18) V где μ1, μ2, μ3 определяются по формулам (13), (14), (16), соответственно.

Приведём примеры применения принципа линейной суперпозиции. Рассмотрим колебания панели вокруг её свободного края, параллельного оси x . Нормированная форма колебаний панели определяется выражением (17), в котором следует положить α ₁ = 0,5; α ₂ = 0,5; α ₃ = 0. Подставляя эти значения в выражение (18) и учитывая зависимости (13), (14), получим

μ = 0,25 μ ₁ + 0,25 μ ₂ =

= C_T ( b / a ) πρ a ² b /2 + C_R ( a / b ) πρ ab ²/16, если a ≤ b , так как использована вторая формула из выражений (13). Для панели с отношением сторон b / a = 4 из табл. 1 и 2 находим C_T (4) = 0,872 и C_R (0,25) = 0,455, и величина присоединённой массы μ = 2,199 πρ a ³.

Выполняя аналогично зависимости (15) приближённую оценку, получим μ = 8 πρ a ³/3 с погрешностью 21% в сторону завышения.

Рассмотрим колебания панели квадратной формы вокруг её диагонали. Нормированная форма колебаний панели определяется выражением (17), в котором следует положить α ₁ = 0; α ₂ = α ₃ = 0,5. Для панели такой формы μ ₂ = μ ₃, и из формулы (18) следует μ = 0,5 μ ₂, где μ ₂ определяется выражением (16), в котором полагается b = a и C_R (1) = 0,816 (табл. 2). Присоединённая масса колебаний панели вокруг диагонали оказывается в два раза меньше присоединённой массы панели при колебаниях вокруг оси x или y , если используется одинаковая нормировка на максимум формы колебаний.

Выше рассматривались колебания панели без просвета между фрагментами, наличие которого влияет на величину присоединённой массы. Ограничимся исследованием влияния просвета на присоединённую массу панели с фрагментами квадратной формы, совершающей поступательные колебания по нормали к её плоскости. Форма колебаний панели f ( x , y ) = 1. Присоединённую массу при условии a ≥ b аналогично зависимости (13) представим в виде

μ = C_T ( a / b , Δ / a )2 πρ b ² a , a ≥ b .

В рассматриваемом случае b = a /2. В табл. 3 приведены значения безразмерного коэффициента присоединённых масс C_T (2, Δ / a ) = 2 μ /( πρ a ³).

Таблица 3

Зависимость присоединённой массы панели с фрагментами квадратной формы от величины просвета между фрагментами. Форма колебаний панели f ( x , y ) = 1

Δ/ a	0,0125	0,025	0,050	0,125	0,250
CT (2, Δ/ a )	0,640	0,627	0,613	0,597	0,587

Заметим, что C_T (2, Δ / a ) → C_T (2) = 0,756, если Δ / a → 0. Если Δ / a → ∞, то фрагменты не влияют друг на друга, поэтому C_T (2, ∞) → C_T (1) = 0,579, так как они имеют квадратную форму. Для относительно малого расстояния между краями фрагментов 2 Δ = 0,025 a их взаимное влияние приводит к увеличению присоединённой массы на 10,5%, и оно уменьшается с увеличением расстояния.

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

Модальные испытания панели СБ могут проводиться возле силовой стены. Приведённые в табл. 3 данные позволяют оценить влияние силовой стены на величину присоединённой массы панели. Можно принимать, что колебания каждого из фрагментов происходят возле силовой стены в плоскости x = 0, так как движение среды симметрично относительно этой плоскости. Если край фрагмента панели находится на расстоянии Δ = 0,25 a от стены, то её влияние приводит к увеличению присоединённой массы фрагмента всего на 1,4%.

Заключение

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